CODEWITHSUNDEEP
pythonnumpymatrix
Yoana G
Yoana G

Reputation: 600

Matrix multiplied by inverse does not yield Identity

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)

Upvotes: 6

Views: 3871

Answers (5)

Ahmad Khan
Ahmad Khan

Reputation: 2693

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

Upvotes: 10

Ethan Koch
Ethan Koch

Reputation: 277

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

Upvotes: 0

Ethan Koch
Ethan Koch

Reputation: 277

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.

Upvotes: 3

Matthieu Brucher
Matthieu Brucher

Reputation: 22023

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.

Upvotes: 2

Mad Physicist
Mad Physicist

Reputation: 114290

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.

Upvotes: 4

Related Questions