Numpy matrix inversion gives wrong results for order greater than 18

user_14425061 0
Viewed 114

I am facing a problem in inverting a matrix with NumPy. Strangely, it gives correct results only up till order 18. As soon as the order becomes greater than 18, it gives erroneous results.

import numpy as np
from decimal import Decimal
import numpy.matlib

I_1=np.matlib.eye(ngrid,ngrid,k=0,dtype=Decimal)
I_2=np.matlib.eye(ngrid,ngrid,k=1,dtype=Decimal)
I_3=np.matlib.eye(ngrid,ngrid,k=2,dtype=Decimal)

B=I_1 + 10.*I_2 + I_3
B=np.divide(B,12.)

B_inv=np.linalg.inv(B)
print B_inv

C=B.dot(B_inv)
print C

The last line was included to check whether it gave the correct result.

numpy python
2 Answers

There is technically nothing wrong with your code. However, you are running into a very fundamental issue in numerical analysis. There are two effects at play here:

  1. Floating point errors
  2. The condition number of a problem

I won’t explain them here, because virtually any introductory book on numerical analysis (e.g. Kendall Atkinson) does a better job than I could do here. I’ll only leave you with this code snippet:

import numpy as np

NGRID = 8


def HilbertMatrix(N, dtype=np.float64):
    arr = np.zeros((N,N), dtype=dtype)
    for i in range(N):
        for j in range(N):
            arr[i,j] = 1 / (i + j + 1)
    return arr


H = HilbertMatrix(NGRID)
J = np.linalg.inv(H)

C = np.dot(H, J)

print(np.allclose(C, np.eye(NGRID)))
print(np.linalg.norm(C))

Try playing around with NGRID and see how far away C will be from the identity matrix.

user_106396940

This is a a numerical analysis issue.

You're holding decimal objects in your numpy array for some reason, and this gets converted into the Object dtype. When you try to perform matrix inversion, the matrix gets converted into float (I'm guessing with low precision).

The condition number of a matrix helps us understand how hard it is to invert the matrix (see wikipedia for a more formal definition). If you check B using np.linalg.cond (make sure to convert to a numpy dtype first) you will see the condition number is around 1e18 (which is very very high).

Another way you can think of it if you know some linear algebra - the determinant of a matrix allows us to determine whether the matrix is singular or not. Det(B) = 0 means the matrix is singular. In your case you created a matrix for which calculating the determinant is very easy (upper triangular).

The matrix is upper triangular, for which the determinant is simply the product of the diagonal, which is (1/12)**(n_grid). This is why as ngrid grows, the matrix is nearly singular, and for the floating point precision you have (which I'm guessing gets converted to 16 under the hood), the matrix is essentially singular / non-invertible.

If you fix your code, which I took the liberty of doing, the implicit np dtype is float64, and the inversion works as desired.

I took the liberty of slightly modifying your code to make it work (ngrid wasn't defined) (changed ngrid to 19 to check).

import numpy as np
ngrid = 19
I_1 = np.eye(ngrid, ngrid, k=0)
I_2 = np.eye(ngrid, ngrid, k=1)
I_3 = np.eye(ngrid, ngrid, k=2)
B = I_1 + 10.*I_2 + I_3
B = np.divide(B, 12.)
print(np.linalg.cond(B))
B_inv = np.linalg.inv(B)
C = B.dot(B_inv)
assert np.diag(C).sum() == ngrid
user_94050220
Related