python - Numpy inverts a non-invertible matrix

Question

My friend and I executed this lines of code in Python 2 and Python 3:

import numpy as np  
mat = np.array([[1,0,0],[-1,3,3],[1,2,2]]) 
np.linalg.inv(mat)

Which returns:

array([[  1.00000000e+00,   0.00000000e+00,   0.00000000e+00],
   [  1.50119988e+16,   6.00479950e+15,  -9.00719925e+15],
   [ -1.50119988e+16,  -6.00479950e+15,   9.00719925e+15]])

Which is odd given:

np.linalg.matrix_rank(mat)

returns 2, thus indicating that the matrix is not invertible.

I understand from this thread that is probably due to the way numpy and python handle floating point numbers, although my matrix consists of whole numbers.

Is there a particular reason why mat breaks numpy's inverse implementation?

Answer 1

As DYZ pointed out the matrix is not invertible because it's rank is 2 not 3.

The reason you are getting such results is because numpy is using LU decomposition to calculate the inverse. This algorithm can yield results even in cases when your matrix is singular. Read linked wikipedia article if you are interested in details.

Note that produced 'inverse' is out of whack. So if you try to use it to solve systems of linear equations it will most likely give you bunch of NaNs and Infs.

I guess numpy does not check the quality of results which is common for high performance libraries. You can do such check yourself very cheaply by multiplying your original matrix by the supposed inverse and checking if numbers on diagonal are close to 1 and other numbers are zeroes. They will not necessarily be exactly equal to zero or one due to nature of floating point computation

As DSM pointed out the condition number of your matrix is really high.

>> cond(A)
ans =   2.4956e+16

So you are loosing 16 digits of precision due to such ill-conditioned matrix. On top of error caused by floating point imprecision.

By the way as others pointed out above your example doesn't work in Numpy 1.12.0

>>> import numpy as np
>>> np.version.version
'1.12.0'

>>> import numpy as np
>>> mat = np.array([[1,0,0],[-1,3,3],[1,2,2]])
>>> np.linalg.inv(mat)
Traceback (most recent call last):
  File "/Users/vlad/.pyenv/versions/CourseraDL/lib/python3.4/site-packages/numpy/linalg/linalg.py", line 90, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
numpy.linalg.linalg.LinAlgError: Singular matrix
>>>