7

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.

ShanZhengYang
  • 16,511
  • 49
  • 132
  • 234
  • I don't believe Python uses standard IEEE floating points, I thought they used arbitrary precision decimals for things that size. Then again, I could be mistaken on what that implies. – Rob Foley Jul 02 '15 at 15:44
  • 2
    What are the determinant (`numpy.linalg.det(m)`) and the condition number (`numpy.linalg.cond(m)`) of your matrix? – DSM Jul 02 '15 at 15:48
  • @DSM `print np.linalg.cond(matrix)` outputs `2.72245023716e+20`. For the inverse, `print np.linalg.cond(new_matrix)` outputs `1.74731032174e+19`. – ShanZhengYang Jul 02 '15 at 15:56
  • @DSM For the determinants, there's a problem with numpy: `print np.linalg.det(matrix)` outputs `0.0` and for the inverse, `print np.linalg.det(new_matrix)` outputs `inf`. What would be the correct `numpy.linalg` operations to give the correct values? – ShanZhengYang Jul 02 '15 at 15:59
  • 2
    If the determinant of the matrix is zero, then it doesn't have an inverse. Are you really sure the zero determinant on A is incorrect? I did some test with matrices of random numbers 400x400 and it works, with minor errors due to floating point operations. –  Jul 02 '15 at 17:38
  • 2
    @Sauruxum numpy.linalg.det() often gives the incorrect answer. The operation does slogdet() first and then takes the exponential. See the code. Also, the matrix DOES have an inverse. There's no reason why not. – ShanZhengYang Jul 02 '15 at 17:43
  • 1
    @RobFoley no, what gave you that idea? Python and numpy both use standard IEEE floats. – ali_m Jul 02 '15 at 19:08
  • I was thinking of Python's integers. Though Python does have the decimal.Decimal class. – Rob Foley Jul 02 '15 at 19:20
  • 1
    Do you really *need* to invert the matrix? If you're trying to solve a linear system of equations it's [faster and more accurate](http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/) to use `np.linag.solve` to compute the solution using factorization than to invert the matrix directly. – ali_m Jul 04 '15 at 21:51
  • Generally speaking, consider using a pseudo-inverse ( np.pinv() ), when unsure about the Matrix properties and whether it is invertible or not. np.pinv() uses Singular value decomposition to calculate inverse of A which is mathematically close to A^-1. It is especially useful when ill-conditioned matrix or non square matrix. some cases the assertion np.allclose will fail, read [here](https://stackoverflow.com/questions/71719948/when-will-numpy-linalg-inv-and-numpy-linalg-pinv-give-very-different-values) for why. – Ritwick Jha Aug 01 '23 at 10:52

2 Answers2

6

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.

rth
  • 10,680
  • 7
  • 53
  • 77
  • 2
    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` – ali_m Jul 02 '15 at 19:11
  • 1
    Sorry, that should be `np.finfo(matrix.dtype).eps` – ali_m Jul 05 '15 at 11:54
1

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

Community
  • 1
  • 1
gg349
  • 21,996
  • 5
  • 54
  • 64
  • 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? – ShanZhengYang Jul 02 '15 at 16:31
  • 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. – gg349 Jul 02 '15 at 16:46
  • See my edit above. I'm pretty sure `numpy.linalg.inv()` results in numerical errors with matrices this large. – ShanZhengYang Jul 02 '15 at 17:04
  • 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) – gg349 Jul 02 '15 at 18:03
  • What techniques do you recommend to invert large matrices to a high-level of precision? – ShanZhengYang Jul 03 '15 at 17:45
  • 3
    See what a preconditioner is [here](https://en.wikipedia.org/wiki/Preconditioner). I recommend to not invert the matrix, depending on the problem you're trying to solve.. – gg349 Jul 03 '15 at 21:53
  • 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? – ShanZhengYang Jul 04 '15 at 15:08
  • @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? – ShanZhengYang Jul 05 '15 at 21:37
  • Consider posting another question explaining better theproblem – gg349 Jul 06 '15 at 05:20
  • 1
    @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]))`). – ali_m Jul 08 '15 at 18:13
  • @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. – rth Jul 08 '15 at 18:58
  • 1
    At least in the current dev version of numpy, `inv` also uses `gesv`. See [here](https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src#L1560-L1580) - the `lapack_func` used in `inv` is `gesv` - it just passes an identity matrix as the "B" parameter. – ali_m Jul 08 '15 at 19:30