Matrix multiplied by inverse does not yield Identity

Question

I am trying to find the inverse matrix of a given matrix using the np.linalg.inv() function. inv yields a matrix which looks alright, but then when trying to multiply the original matrix by inverse the output is not the identity as supposed by the inverse matrix definition.

from numpy.linalg import inv

M = np.random.random((4, 4))
Mi = inv(M)
I = M @ Mi # using matrix multiplication operator
I.astype(int) #as the output looks like 2.77555756e-17

>>> array([[0, 0, 0, 0],
           [0, 0, 0, 0],
           [0, 0, 0, 0],
           [0, 0, 0, 1]])

Which is clearly not the identity (gives a slightly different answer when running multiple times)

Answer 1

Try rounding it first before converting to int.

np.around(I).astype(int)

Creating a random matrix:

>>> M = np.random.random((4,4))
>>> M
array([[0.51351957, 0.57864882, 0.0489495 , 0.85066216],
       [0.60052988, 0.93844708, 0.74651889, 0.17237584],
       [0.26191596, 0.46451226, 0.46514401, 0.81917544],
       [0.19247662, 0.82801899, 0.83839146, 0.08531949]])

Taking inverse:

>>> from numpy.linalg import inv
>>> Mi = inv(M)
>>> Mi
array([[-1.3515514 ,  3.53647196,  1.0391335 , -3.64654487],
       [ 2.76188122, -2.23981308, -2.74634579,  3.35680468],
       [-2.44320291,  1.47102487,  2.36135635, -1.28451339],
       [ 0.2533113 , -0.69591469,  1.10498293, -0.00818495]])

Now, multiplying M and Mi should produce identity.

>>> M @ Mi
array([[ 1.00000000e+00, -4.44089210e-16, -1.11022302e-16, -6.93889390e-18],
       [-4.16333634e-17,  1.00000000e+00, -8.32667268e-17, -8.60856525e-17],
       [ 5.55111512e-17, -2.22044605e-16,  1.00000000e+00, -1.57859836e-16],
       [ 6.24500451e-17, -8.32667268e-17, -2.35922393e-16, 1.00000000e+00]])

But this is clearly not identity. But if you look closely, the diagonal values are pretty close to 1, and all other values are really small numbers (almost zeros), like -16 or -17 in the exponent.

This error is, because float values are never the exact values, they always have some error in them. Have a look at the article 15. Floating Point Arithmetic: Issues and Limitations and Is floating point math broken?.

Now, if we just convert it to int, the chances are, that it will still not be an identity. Because a value which is really close to 1, it can actually be a little smaller than 1, resulting in 0 when casted to int.

>>> (M @ Mi).astype(int)
array([[1, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 1, 0],
       [0, 0, 0, 0]])

But if you round it first before converting to int, you'll get an identity.

>>> np.around(M @ Mi).astype(int)
array([[1, 0, 0, 0],
       [0, 1, 0, 0],
       [0, 0, 1, 0],
       [0, 0, 0, 1]])