# [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.
>>>
>>>     ===
>>>
>>>     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
```