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.

Solution 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.

Solution 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).

Comments

[ShanZhengYang]

If this is correct, my test should be to take the product, correct? The original matrix is matrix. The inverse is new_matrix. Thus np.dot(matrix, new_matrix) should give us the identity matrix, correct?

[gg349]

it should, unless you are hitting numerical errors. Try for example from numpy.random import rand; from scipy.linalg import inv; A = rand(4,4); print(dot(A,inv(A))) and you'll get a unitary matrix.

[ShanZhengYang]

See my edit above. I'm pretty sure numpy.linalg.inv() results in numerical errors with matrices this large.

[gg349]

you shouldn't be that sure of numerical errors for matrices of this size: using a random matrix generated with A=rand(400,400) gives a result res = dot(A,inv(A)) with off-diagonal elements smaller than 1e-11 in absolute value (It depends of course on the matrix itself, not just on the size)

[ali_m]

It should depend on the epsilon for the dtype of the array (i.e. np.finfo(matrix).eps), which may not be the same as sys.float_info.epsilon

[ShanZhengYang]

What techniques do you recommend to invert large matrices to a high-level of precision?

[gg349]

See what a preconditioner is here. I recommend to not invert the matrix, depending on the problem you're trying to solve..

[ShanZhengYang]

I'll try to reformulate the problem. "in practise, you'd call scipy.linalg.solve in python for example" How would you use scipy.linalg.solve to invert a matrix?

[ali_m]

Sorry, that should be np.finfo(matrix.dtype).eps

[ShanZhengYang]

@gg349 I am obviously a bit out of my depth here. I need to invert the matrix in order to solve a problem---I'm not sure how to reformulate this equation. Do you have any resources for computing linear algebra problems for scientific computing?

[gg349]

Consider posting another question explaining better theproblem

[ali_m]

@rth "You can't use scipy.linalg.solve to invert a matrix" - Sure you can! In fact, that's how np.linalg.inv works (it does the equivalent of np.linalg.solve(A, np.eye(A.shape[0]))).

[rth]

@ali_m hmm, yes I did not think about this too much. They use different BLAS calls (gesv for solve vs getrf followed by getri for inv), but the result is indeed the same (and with only ~10% computational overhead for the solve call with 1000x1000 matrices.) and both are based on LU factorisation, so you are right.

[ali_m]

At least in the current dev version of numpy, inv also uses gesv. See here - the lapack_func used in inv is gesv - it just passes an identity matrix as the "B" parameter.