Is numpy.linalg.inv() giving the correct matrix inverse? EDIT: Why does inv() gives numerical errors?

Question

I have a matrix shaped (4000, 4000) and I would like to take the inverse. (My intuition of the inverting matrices breaks down with such large matrices.)

The beginning matrix has values of the magnitude e-10, with the following values: print matrix gives an output

[[  2.19885119e-10   2.16462810e-10   2.13062782e-10 ...,  -2.16462810e-10
   -2.19885119e-10  -2.16462810e-10]
 [  2.16462810e-10   2.19885119e-10   2.16462810e-10 ...,  -2.13062782e-10
   -2.16462810e-10  -2.19885119e-10]
 [  2.13062782e-10   2.16462810e-10   2.19885119e-10 ...,  -2.16462810e-10
   -2.13062782e-10  -2.16462810e-10]
 ..., 
 [ -2.16462810e-10  -2.13062782e-10  -2.16462810e-10 ...,   2.19885119e-10
    2.16462810e-10   2.13062782e-10]
 [ -2.19885119e-10  -2.16462810e-10  -2.13062782e-10 ...,   2.16462810e-10
    2.19885119e-10   2.16462810e-10]
 [ -2.16462810e-10  -2.19885119e-10  -2.16462810e-10 ...,   2.13062782e-10
    2.16462810e-10   2.19885119e-10]]

I then use NumPy's numpy.linalg.inv() to invert the matrix.

import numpy as np
new_matrix = np.linalg.inv(matrix)
print new_matrix

This is the output I get in return:

[[  1.95176541e+25   9.66643852e+23  -1.22660930e+25 ...,  -1.96621184e+25
   -9.41413909e+24   1.33500310e+25]
 [  2.01500967e+25   1.08946558e+24  -1.25813014e+25 ...,  -2.07717912e+25
   -9.86804459e+24   1.42950556e+25]
 [  3.55575106e+25   2.11333704e+24  -2.25333936e+25 ...,  -3.68616202e+25
   -1.72651875e+25   2.51239524e+25]
 ..., 
 [  3.07255588e+25   1.61759838e+24  -1.95678425e+25 ...,  -3.15440712e+25
   -1.47472306e+25   2.13570651e+25]
 [ -7.24380790e+24  -8.63730581e+23   4.90519245e+24 ...,   8.30663797e+24
    3.70858694e+24  -5.32291734e+24]
 [ -1.95760004e+25  -1.12341031e+24   1.23820305e+25 ...,   2.01608416e+25
    9.40221886e+24  -1.37605863e+25]]

That's a huge difference! How could that be? A matrix of magnitude e-10 is inverted to a matrix of magnitude e+25?

Is this mathematically correct, or are the IEEE floating point values breaking down?

If this is mathematically correct, could someone explain to me the mathematical intuition behind this?

EDIT:

Following the comments below, I decided to test.

np.dot(matrix, new_matrix) should give the identity matrix, A * A^T = Identity.

This is my output:

[[  0.   -3.  -16.  ...,  16.    8.   12. ]
 [-24.   -1.5  -8.  ...,  32.   -4.   36. ]
 [ 40.    1.  -64.  ...,  24.   20.   24. ]
 ..., 
 [ 32.   -0.5  48.  ..., -16.  -20.   16. ]
 [ 40.    7.   16.  ..., -48.  -36.  -28. ]
 [ 16.    3.   12.  ..., -80.   16.    0. ]]

Why does numpy.linalg.inv() result in numerical errors?

np.allclose( np.dot(matrix, new_matrix), np.identity(4000) )

gives False.

Answer 1

Your matrix is ill-conditionned, since

np.linalg.cond(matrix) > np.finfo(matrix.dtype).eps

According to this answer you could consider using Singular Value Decomposition to inverse such matrices.

Answer 2

For the determinant of the 2 matrices, you have that

det(A) * det(A^{-1}) = 1

so that if det(A) is big, then det(A^{-1}) is small. For the norm of the 2 matrices, (if you pick a sub-multiplicative norm), you have:

1  =  |A*A^{-1}| >= |A| |A^-1|

where || is a reasonable choice of a norm that is sub-multiplicative. Here you have the intuition of what you are observing numerically: if the >= sign is actually a ~=, you recover the same observation that is strictly true for the determinant.

The same reasoning applies if you consider the product

A * A^{-1} = 1

for a matrix A with all positive elements. For the elements on the diagonal of 1 at the RHS, you need very small numbers from A^{-1} if the elements of A are very big.

PS: Notice however that this does not PROVE that this trend always holds. This just provides mathematical intuition of why you observe this scaling.

EDIT, in reply to the comments:

Originally the question was "If this is mathematically correct, could someone explain to me the mathematical intuition behind this?". And indeed, it is mathematically correct and sound that given a matrix with small numbers, the inverse will have large numbers. Above I explain why this is the case.

To answer the other question that came up in the OP's edit, which is why inv() results in numerical errors: inverting matrices is a HARD problem. That is why every time we can, we avoid inverting them. For example, for the problem

A x = b

we do not calculate the inverse of A, but use other algorithms (in practise, you'd call scipy.linalg.solve in python for example).