6

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)

Mad Physicist
  • 107,652
  • 25
  • 181
  • 264
Yoana G
  • 580
  • 5
  • 8

5 Answers5

9

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]])
Ahmad Khan
  • 2,655
  • 19
  • 25
3

When you print out I, it looks like this:

array([[ 1.00000000e+00, -5.55111512e-17, -1.66533454e-16, 0.00000000e+00],
       [ 6.38378239e-16,  1.00000000e+00, -5.55111512e-17, 0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00, 0.00000000e+00],
       [-5.55111512e-17, -1.11022302e-16, -1.24900090e-16, 1.00000000e+00]])

However, the 1.00 entries are not exact. When you print 1 - I, you can see this:

array([[-2.22044605e-16,  1.00000000e+00,  1.00000000e+00, 1.00000000e+00],
       [ 1.00000000e+00,  2.22044605e-16,  1.00000000e+00, 1.00000000e+00],
       [ 1.00000000e+00,  1.00000000e+00,  0.00000000e+00, 1.00000000e+00],
       [ 1.00000000e+00,  1.00000000e+00,  1.00000000e+00, 0.00000000e+00]])

The positive diagonal entries represent values in I that are slightly less than one. When you do integer truncation (which is what astype(int) does), you will set those elements to zero. Instead, round the values to the nearest integer instead of truncating:

np.around(I).astype(int)

However, you will not always have integer inputs like this, in which case rounding will be misleading. Numpy provides the allclose funcion for comparing values within a tolerance:

np.allclose(I, np.identity(I.shape[0]), rtol=0, atol=1e-15)

You can also do an element-wise check using isclose:

np.isclose(I, np.identity(I.shape[0]), rtol=0, atol=1e-15)

I set the relative tolerance to zero since it is multiplied by the elements of the second matrix, making it useless in this situation.

Mad Physicist
  • 107,652
  • 25
  • 181
  • 264
2

The issue is that the astype function does not round, it simply truncates. So, the reason you do not see the identity matrix is that the other values that should be 1 were somewhere around 0.99999. You can use this:

import numpy as np
a = np.random.random((4,4))
b = np.linalg.inv(a)
# 0.00001 is the tolerance about which you which to consider values to be == 1
c = np.array(a@b + 0.00001, dtype=int)
print(c)

If you want to simply round (high tolerance == 0.5) use this instead:

import numpy as np
a = np.random.random((4,4))
b = np.linalg.inv(a)
c = np.array(np.round(a@b), dtype=int)
print(c)

Additionally, it is likely best practice to use the full np.linalg.inv function.

Ethan Koch
  • 267
  • 1
  • 6
  • It's not the best solution, as a result that is 50% off would consider the result as good when it's not. – Matthieu Brucher Oct 29 '18 at 14:16
  • I have updated it -- the first solution includes a tolerance threshold – Ethan Koch Oct 29 '18 at 15:02
  • This is still far from correct and doesn't properly check that the result is close to the identity matrix. The proper way to do it is @MadPhysicist answer. – Matthieu Brucher Oct 29 '18 at 15:04
  • "correct" is subjective -- and my solution is significantly more readable than the @MadPhysicist answer, which also matters. – Ethan Koch Oct 29 '18 at 15:07
  • the author wanted to know how to get the identity, not check if it was the identity -- that is a different question – Ethan Koch Oct 29 '18 at 15:08
  • Yes, that's what he wanted, but used the wrong way to check for it: "when trying to multiply the original matrix by inverse the output is not the identity as supposed by the inverse matrix definition". Let's give good practices and the proper way of doing things. – Matthieu Brucher Oct 29 '18 at 15:10
  • understood -- while improve in the future. – Ethan Koch Oct 29 '18 at 15:13
1

The question is, is the matrix close to np.eye(4) or not.

This is how you should check:

I = M@Mi
EPS = 1e-8
r = np.all(np.abs(I - np.eye(4)) < EPS)

r will indicate if the two matrices (I and the identity) are close up to 1e-8.

Mad Physicist
  • 107,652
  • 25
  • 181
  • 264
Matthieu Brucher
  • 21,634
  • 7
  • 38
  • 62
0

What is your goal here?

It seems like you just want to know how to get the identity matrix and compare if the results of your matrix multiplication are close to it.

If so, this is what you should do:

import numpy as np

matrix1 = np.random.random((4,4))
matrix1_inv = np.linalg.inv(matrix1)

# get what may be identity matrix
ident_pred = matrix1 @ matrix1_inv

# get what is the identity matrix
ident_true = np.identity(matrix1.shape[0])

# check that they are the same
print(np.allclose(ident_pred, ident_true))
Ethan Koch
  • 267
  • 1
  • 6