CODEWITHSUNDEEP
pythonnumpymatrix
lamb_da_calculus
lamb_da_calculus

Reputation: 295

Why does numpy say this matrix is singular *sometimes*?

Part of my code inverts a matrix (really an ndarray) using numpy.linalg.inv. However, this frequently errors out as follows:

numpy.linalg.linalg.LinAlgError: Singular matrix

That would be fine if the matrix was actually singular. But that doesn't seem to be the case.

For example, I'm printing the matrix before trying to invert it. So right before the error it prints this:

[[ 0.76400334  0.22660491]
[ 0.22660491  0.06721147]]

... and then returns the above singularity error when it tries to invert that matrix. But from what I can tell this matrix is invertible. Numpy even seems to agree when asked later.

>>> numpy.linalg.inv([[0.76400334, 0.22660491], [0.22660491,    0.06721147]])
array([[  2.88436275e+07,  -9.72469076e+07],
   [ -9.72469076e+07,   3.27870046e+08]])

Here's the exact code snippet:

print np.dot(np.transpose(X), X)
print np.linalg.inv(np.dot(np.transpose(X),X))

The first line prints the matrix above; the second line fails.

So what distinguishes the two actions above? Why does the stand-alone code work even though it errors out in my script?

EDIT: Per Colonel Beauvel's request, if I do

try:
    print np.dot(np.transpose(X), X)
    z = np.linalg.inv(np.dot(np.transpose(X), X))
except:
    z = "whoops"
print z

it outputs

[[ 0.01328185  0.1092696 ]
[ 0.1092696   0.89895982]]
whoops

but trying this on its own I get

>>> numpy.linalg.inv([[0.01328185, 0.1092696], [0.1092696, 0.89895982]])
array([[  2.24677775e+08,  -2.73098420e+07],
   [ -2.73098420e+07,   3.31954382e+06]])

Upvotes: 4

Views: 3559

Answers (2)

user2379410
user2379410

Reputation:

It's a matter of printing precision. The IEEE 754 doubles, that you're most likely using, have about 16 decimal digits of precision and you need to write out 17 to preserve the binary value.

Here's a small example. First create a singlular matrix:

In [1]: import numpy as np

In [2]: np.random.seed(0)

In [3]: a, b, c = np.random.rand(3)

In [4]: d = b*c / a

In [5]: X = np.array([[a, b],[c, d]])

Print and try to invert it:

In [6]: X
Out[6]: 
array([[ 0.5488135 ,  0.71518937],
       [ 0.60276338,  0.78549444]])

In [7]: np.linalg.inv(X)
LinAlgError: Singular matrix

Try to invert the printed matrix:

In [8]: Y = np.array([[ 0.5488135 ,  0.71518937],
   ...:               [ 0.60276338,  0.78549444]])

In [9]: np.linalg.inv(Y)
Out[9]: 
array([[-85805775.2940297 ,  78125795.99532071],
       [ 65844615.19517545, -59951242.76033063]])

Succes!

Increase printing precision and try again:

In [10]: np.set_printoptions(precision=17)

In [11]: X
Out[11]: 
array([[ 0.54881350392732475,  0.71518936637241948],
       [ 0.60276337607164387,  0.78549444195576024]])

In [12]: Z = np.array([[ 0.54881350392732475,  0.71518936637241948],
    ...:               [ 0.60276337607164387,  0.78549444195576024]])

In [13]: np.linalg.inv(Z)
LinAlgError: Singular matrix

Upvotes: 4

Colonel Beauvel
Colonel Beauvel

Reputation: 31161

I just compute the determinant:

In [130]: m = np.array([[ 0.76400334, 0.22660491],[ 0.22660491,0.06721147]])

In [131]: np.linalg.det(m)
Out[131]: 2.3302017068132921e-09

# which is in fact for a 2D matrix 0.76400334*0.06721147 - 0.22660491*0.22660491

Which is already quit close to 0.

If a matrix m can be inverted, mathematically you can compute the adjoint and divide by the determinant to get the inverted matrix. Numerically if the determinant is too small, this can entail the kind of error you have ...

Upvotes: 1

Related Questions