[ncl-talk] inverse_matrix problem
Walter Hannah
whannah at atmos.colostate.edu
Tue Nov 11 07:07:17 MST 2014
Maria,
Thanks for this insight. I guess I'll just have to dust off my linear
algebra book.
Walter
On 11/10/14 6:35 PM, Maria Gehne - NOAA Affiliate wrote:
> 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 <mailto: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
>> <mailto:whannah at atmos.colostate.edu>>
>> Date: Tuesday, 11 November 2014 9:52 AM
>> To: Dennis Shea <shea at ucar.edu <mailto:shea at ucar.edu>>
>> Cc: "ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>"
>> <ncl-talk at ucar.edu <mailto: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
>>> <mailto: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/20141111/52dcc9dc/attachment.html
More information about the ncl-talk
mailing list