[ncl-talk] inverse_matrix problem
Maria Gehne - NOAA Affiliate
maria.gehne at noaa.gov
Mon Nov 10 16:35:36 MST 2014
The matrix has three rows columns equal to zero, which makes it singular.
Several of the eigenvalues are zero and that's probably why inverse_matrix
doesn't work and returns all missing values. You will have the same problem
if you just replace all zeros with the same non-zero number. In that case
you have several eigenvalues that are the same and the solution to the
linear system that Will suggested does not have a unique solution.
I attached a matlab code that computes the three examples Dennis sent and I
added one more that just replaces the zeros with 7 and you still have a
singular matrix. So non-invertable.
I don't know what you are trying to do, but you could just delete those
three column/rows that have only zeros, since they don't provide any
additional information about the matrix, and invert the rest of the matrix.
So, shrinking you matrix and inverting only on the subspace where the
matrix is non-singular.
Maria
On Mon, Nov 10, 2014 at 4:25 PM, Walter Hannah <whannah at atmos.colostate.edu>
wrote:
> Will,
>
> I tried your suggestion, but all I got back was missing values.
>
> After Dennis's suggestion, I tried replacing the zeros with tiny random
> numbers, using this code:
>
> R = random_normal( 0.,1e-12, (/nlev*2,nlev*2/))
> M = where(M.eq.0,R,M)
> iM = inverse_matrix(M)
>
> Interestingly, this worked perfectly, and I got what I expected.
> I also tried using this approach to trying Will's suggestion, and that
> also worked with the random numbers!
>
> I can't figure out a mathematical explanation, but this seems like a
> pretty nice fix. Might be good to add this to the online documentation.
>
> Thanks,
> Walter
>
>
> On 11/10/14 6:17 PM, Will Hobbs wrote:
>
> Walter
>
> Instead of using the inverse matrix() function, try using
> solve_linsys(A, B) instead, for A = your matrix and B = I (the identity
> matrix). Since A * inv(A) = I, the output will be the inverse.
>
> Mathematically this is no different to inverse_matrix(), but the code is
> a little different. In particular, if the matrix is singular (i.e. the
> determinant = 0, and it cannot be inverted) then it will return an error
> code.
> https://www.ncl.ucar.edu/Document/Functions/Built-in/solve_linsys.shtml
>
> This might help you figure out if the matrix itself is the problem (I
> suspect, as Dennis has suggested, that may be the issue).
>
> Will
>
>
> From: Walter Hannah <whannah at atmos.colostate.edu>
> Date: Tuesday, 11 November 2014 9:52 AM
> To: Dennis Shea <shea at ucar.edu>
> Cc: "ncl-talk at ucar.edu" <ncl-talk at ucar.edu>
> Subject: Re: [ncl-talk] inverse_matrix problem
>
> Dennis,
>
> Wow! Thanks for going above and beyond! This helps a lot. I'll try to look
> for a way to get a modified matrix that doesn't have all those zeros.
>
> Walter
>
> On 11/10/14 5:50 PM, Dennis Shea wrote:
>
> You have rows of 0.0 ... I think that is causing issues. I'm sure some
> 'math-whiz' could answer but I attached a simple test code that reads your
> data (thx for providing). It calculates the inverse and the identity
> matrix. Note that NCL is using a direct call to the documented LAPACK
> routines.
>
> ===
> TEST1: duplicates your findings
>
> TEST 2: create a bogus 20x20 array of random numbers. No row consists of
> all constants. Identity matrix is successfully calculated.
>
> TEST 3: wherever the original array had a row of all 0.0, substitute
> random numbers. The identity matrix is successfully calculated.
>
> On Mon, Nov 10, 2014 at 1:57 PM, Walter Hannah <
> whannah at atmos.colostate.edu> wrote:
>
>> Hi,
>>
>> I'm trying to use the "inverse_matrix" function. The matrix doesn't have
>> any missing values, but I still get this error message:
>>
>> Variable: M
>> Type: float
>> Total Size: 1600 bytes
>> 400 values
>> Number of Dimensions: 2
>> Dimensions and sizes: [pertlev | 20] x [tendlev | 20]
>> Coordinates:
>> pertlev: [50..1950]
>> tendlev: [50..1950]
>> Number Of Attributes: 2
>> missing_value : 9.96921e+36
>> _FillValue : 9.96921e+36
>> (0)
>> ----------------------------------------------------------------------------
>>
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 -0.31 -0.06 0.15 -0.27 -0.56 -1.21 -1.38 -1.36 1.49
>> 0.00 0.00 4.65 0.37 0.46 0.47 0.52 0.55 0.77 2.49
>> 0.00 0.26 -0.35 -0.29 -0.65 -1.13 -2.24 -2.45 -2.47 2.91
>> 0.00 0.00 4.07 2.45 0.95 0.87 0.90 1.05 1.55 4.81
>> 0.00 -0.11 1.43 -0.91 -1.69 -1.52 -2.63 -2.81 -2.69 3.17
>> 0.00 0.00 -2.35 2.61 1.49 1.01 1.00 1.18 1.78 5.32
>> 0.00 -0.02 0.55 1.08 -1.74 -2.04 -2.59 -2.71 -2.39 2.63
>> 0.00 0.00 -0.75 -1.10 1.35 1.02 0.89 1.00 1.57 4.87
>> 0.00 -0.06 0.16 0.52 0.72 -1.86 -2.88 -2.74 -2.29 2.31
>> 0.00 0.00 0.31 -1.30 -0.26 0.87 0.75 0.82 1.27 4.50
>> 0.00 -0.13 -0.07 0.20 0.54 0.27 -2.74 -3.10 -2.27 2.04
>> 0.00 0.00 0.08 -0.49 -0.47 0.14 0.72 0.69 1.13 4.24
>> 0.00 -0.15 -0.06 0.20 -0.17 0.02 -0.77 -3.87 -3.13 2.39
>> 0.00 0.00 0.69 -0.08 0.09 0.03 0.36 0.90 1.29 4.59
>> 0.00 -0.14 0.05 0.09 -0.17 -0.29 -0.14 -0.66 -5.36 2.00
>> 0.00 0.00 0.39 0.17 0.14 0.07 -0.01 0.22 1.00 4.05
>> 0.00 0.02 -0.07 -0.01 0.00 0.02 0.43 0.94 0.77 -2.90
>> 0.00 0.00 -0.09 0.09 -0.04 -0.07 -0.17 -0.26 -0.49 0.54
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
>> 0.00 -0.01 -0.04 0.12 -0.02 -0.01 -0.03 -0.05 -0.05 0.04
>> 0.00 0.00 -2.17 -0.25 -0.14 0.01 0.02 0.02 0.02 0.12
>> 0.00 0.00 -0.19 0.11 0.28 0.06 0.11 0.06 0.07 -0.13
>> 0.00 0.00 2.75 -1.97 -0.04 -0.12 -0.05 -0.05 -0.08 -0.10
>> 0.00 -0.02 -0.14 -0.08 0.31 0.21 0.13 0.08 0.06 -0.08
>> 0.00 0.00 0.48 1.11 -1.09 -0.13 -0.10 -0.06 -0.09 -0.12
>> 0.00 -0.01 -0.03 -0.11 -0.07 0.32 0.27 0.12 0.18 -0.17
>> 0.00 0.00 0.02 0.63 0.52 -0.94 -0.02 -0.15 -0.08 -0.25
>> 0.00 0.02 0.08 -0.06 -0.24 0.10 0.43 0.26 0.16 -0.17
>> 0.00 0.00 0.07 0.31 0.35 0.40 -1.06 -0.06 -0.20 -0.26
>> 0.00 0.06 0.03 -0.10 0.07 0.11 1.03 1.38 0.53 -1.10
>> 0.00 0.00 -0.44 0.14 -0.03 0.08 0.19 -1.61 -0.76 -1.16
>> 0.00 0.11 -0.03 -0.08 0.07 0.16 0.46 1.28 1.65 -1.21
>> 0.00 0.00 -0.25 0.07 -0.06 0.00 0.14 0.41 -3.04 -0.66
>> 0.00 0.12 -0.30 -0.15 0.80 1.58 2.86 3.62 5.02 -3.15
>> 0.00 0.00 -1.23 -0.91 -0.80 -0.94 -0.85 -0.66 0.21 -10.64
>>
>> (0)
>> ----------------------------------------------------------------------------
>> (0) # missing values: 0
>> (0) # values = 0.0 : 111
>> warning:inverse_matrix: info = 1; missing values not allowed
>>
>>
>> I thought it might have something to do with those zero values, but if I
>> replace them with non-zero values I get a similar, but not identical,
>> result:
>>
>> ...
>> (0)
>> ----------------------------------------------------------------------------
>> (0) # missing values: 0
>> (0) # values = 0.0 : 0
>> warning:inverse_matrix: info = 20; missing values not allowed
>>
>>
>>
>> Thanks,
>> Walter
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20141110/935393f5/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: tst_inv_matrix.m
Type: application/octet-stream
Size: 1360 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20141110/935393f5/attachment.obj
More information about the ncl-talk
mailing list