Disagreement in obtaining inverse of Sympy matrix in two ways

Advertisements

I am attempting to invert a simple numerical matrix in Sympy. I start with a matrix:

In [56]: A_MATRIX_M
Out[56]: 
⎡0.174683794941032   0.174688696013867   0.174688696013957   0.174683794941032   0.174688696013867   0.174688696013957 ⎤
⎢                                                                                                                      ⎥
⎢0.152985639387222   0.153053070078189   0.153053070078123   0.152985677401143   0.153053070078189   0.153053070078123 ⎥
⎢                                                                                                                      ⎥
⎢0.122289719307085    0.12255894161074    0.12255894161075   0.122289780080333    0.12255894161074    0.12255894161075 ⎥
⎢                                                                                                                      ⎥
⎢0.0888623464387251  0.0894599879866164  0.0894599879864983  0.0888623905998649  0.0894599879866165  0.0894599879864984⎥
⎢                                                                                                                      ⎥
⎢0.0584375462953993  0.0593271595778474  0.0593271595779078  0.0584375753365935  0.0593271595778473  0.0593271595779078⎥
⎢                                                                                                                      ⎥
⎣0.034711924284652   0.0356325950288525  0.035632595028856   0.0347119415351321  0.0356325950288525  0.035632595028856 ⎦

I can get the inverse with .inv():

In [59]: A_M_INV.inv()
Out[59]: 
⎡-4.20862006694537e-38  -3.64678015602041e-38  -4.89393679687261e-39  -5.26008046298975e-38  -1.29781892788192e-38  1.12107780599699e-37 ⎤
⎢                                                                                                                                        ⎥
⎢-1.19250990841802e-38  -1.17838767512258e-38  -5.49930182641769e-39  -1.79347381150519e-38  -3.76673961453884e-39  3.67594291864767e-38 ⎥
⎢                                                                                                                                        ⎥
⎢1.27521806859415e-38   1.38241733868129e-38   1.33016305680983e-39   1.89098712400354e-38   7.63333345652249e-39   -4.07574230117266e-38⎥
⎢                                                                                                                                        ⎥
⎢-3.71566091342595e-39  -3.69178932808298e-39  -1.37147644288023e-39  -4.68481431180264e-39  -3.97158856857654e-39  8.39451101073365e-39 ⎥
⎢                                                                                                                                        ⎥
⎢1.47794212524182e-38   1.21859956435038e-38   7.46216287435121e-40   2.08223967916098e-38   8.25097707225834e-39   -4.0465253879286e-38 ⎥
⎢                                                                                                                                        ⎥
⎣-2.5158985894826e-39   -3.62021834376681e-39  -3.99023376900161e-40   -5.881653010039e-39   -3.44892577031322e-39  6.50635833262352e-39 ⎦

However, I should also be able to get the inverse by dividing the adjugate matrix by the matrix determinant. However, it gives a completely different answer:

In [57]: A_M_INV = A_MATRIX_M.adjugate() / A_MATRIX_M.det()

In [58]: A_M_INV
Out[58]: 
⎡-1.03014934740997e+38  9.76285991918083e+37   -2.61712726811554e+38  -2.91755290755826e+38  2.02091254231468e+37   8.60970411767247e+37 ⎤
⎢                                                                                                                                        ⎥
⎢-1.0350460657459e+38   1.09196345494244e+38   1.65391718422458e+38   -4.62921116834264e+37  -4.10946591995717e+38  -2.93534695043802e+38⎥
⎢                                                                                                                                        ⎥
⎢3.44664755653471e+37   -1.86912642859017e+38  -8.76351586189679e+37  -2.26077970623935e+38  -1.67328014514703e+37  1.00789976403353e+38 ⎥
⎢                                                                                                                                        ⎥
⎢1.82916593489888e+37   -8.10152924803704e+37  -1.33638262519626e+38  2.88888937621382e+38   1.19814518627368e+38   -3.22157029803083e+38⎥
⎢                                                                                                                                        ⎥
⎢6.14049097203001e+37   2.94210971355518e+37   5.88396356342516e+37   -2.59752889211905e+38   1.036688778158e+38    1.24211877586664e+38 ⎥
⎢                                                                                                                                        ⎥
⎣-4.62263666898771e+37  2.94055576923073e+37   -1.04165393746291e+38  -2.89791317977202e+37  -5.86035000752247e+37  -1.95539341633241e+38⎦

Here’s an example where someone else uses this method to invert the matrix:

Why is sympy.Matrix.inv slow?

Why are they different? Attempting to use this inverse to solve a simple equation shows that neither seems to be the correct inverse. I’m not sure what could be going wrong; is it because of the large exponents or precision needed to describe the matrix?

Thank you in advance for your help!

EDIT: with @Jared’s help below, I did a little more digging and found this question: Inverse matrix results different in MATLAB and Python now that I know my matrix is ill-conditioned/essentially singular, which precludes it from having a well-defined inverse.

>Solution :

Here are my comments posted as an answer:

First of all, notice that your "inverse" is just the zeros matrix. I created your matrix and calculated the determinant, which turns out to be approximately 0.

import sympy as smp

A = smp.Matrix([[0.174683794941032, 0.174688696013867, 0.174688696013957, 0.174683794941032, 0.174688696013867, 0.174688696013957],
                [0.152985639387222, 0.153053070078189, 0.153053070078123, 0.152985677401143, 0.153053070078189, 0.153053070078123],
                [0.122289719307085, 0.12255894161074, 0.12255894161075, 0.122289780080333, 0.12255894161074, 0.12255894161075],
                [0.0888623464387251, 0.0894599879866164, 0.0894599879864983, 0.0888623905998649, 0.0894599879866165, 0.0894599879864984],
                [0.0584375462953993, 0.0593271595778474, 0.0593271595779078, 0.0584375753365935, 0.0593271595778473, 0.0593271595779078],
                [0.034711924284652, 0.0356325950288525, 0.035632595028856, 0.0347119415351321, 0.0356325950288525, 0.035632595028856]])

print(A.det())  # -9.86673942982700e-57

When you compute the inverse the second way, you get large numbers because you’re dividing by a small number. Your matrix is ill-conditioned and essentially singular (i.e. it doesn’t have an inverse).

Remember, in this case you’re working with floating point values, so the determinant value, which sympy computes to be approximately 0, is 0. I’m not sure why sympy doesn’t consider this matrix singular and returns an inverse (I don’t know what check it performs under the hood). Sympy seems to think it is invertible (sympy says it is full rank when checking A.rank()), but the results seem to say otherwise.

We can also check the eigenvalues of the matrix. (The eigenvalues are the keys to the dictionary returned by Matrix.eigenvals().)

print(A.eigenvals()) # {0.633575259554063: 1, 4.07013965394438e-9: 1, -6.35328826990267e-17: 1, -2.95260889975432e-15: 1, -3.59523027736480e-14: 1, 0.000542688212375667: 1}

Notice that 3 (arguably 4) of them are essentially 0. A matrix must have non-zero eigenvalues to be invertible.

Here is a list of conditions for a matrix to be invertible: https://mathworld.wolfram.com/InvertibleMatrixTheorem.html.

Leave a ReplyCancel reply