0

I've a piece of code to find the inverse of a matrix using Gaussian elimination. To check if the solution is correct, I match it with the numpy linalg.inv solution. They appear to be same, but it returns false when I check. This is what I mean

e = inverse(zeze)

[[ 0.11111111  0.37037037 -0.2962963 ]
 [ 0.88888889 -0.03703704 -0.37037037]
 [-0.22222222 -0.07407407  0.25925926]]

And this is the other one

d = np.linalg.inv(zeze)
print(d)

[[ 0.11111111  0.37037037 -0.2962963 ]
 [ 0.88888889 -0.03703704 -0.37037037]
 [-0.22222222 -0.07407407  0.25925926]]

print(d== e)

[[False False False]
 [False False False]
 [False False False]]

The matrix zeze is

zeze = np.array([[1,2,4],[4,1,6],[2,2,9]])

And the code snippet for the inverse function is

def copier(a):
    row = len(a)
    col = len(a[0])

    copi = np.zeros((row,col))
    for i in range(row):
        for j in range(col):
            copi[i][j] = a[i][j]
    return copi
def inverse(a):
    if len(a) == len(a[0]):
        detm = np.linalg.det(a)
        if detm == 0:
            print("If the determinant of a matrix equals zero, it's impossible to find it's inverse.")
            print("And guess what Sherlock? ")
            print("This one has a determinant zero")
        elif detm != 0:
            
            mat_copy = copier(a)
            id_mat = np.identity(len(a))
            id_mat_copy = copier(id_mat)
            rows_mat = list(range(len(a)))

            for diag in range(len(a)):
                scaler = 1.0/mat_copy[diag][diag]
                for i in range(len(a)):
                    mat_copy[diag][i] *= scaler
                    id_mat_copy[diag][i] *= scaler
                for j in rows_mat[0:diag] + rows_mat[diag+1:]:
                    crScaler = mat_copy[j][diag]
                    for k in range(len(a)):
                        mat_copy[j][k] = mat_copy[j][k] - crScaler * mat_copy[diag][k]
                        id_mat_copy[j][k] = id_mat_copy[j][k] - crScaler * id_mat_copy[diag][k]
        print(id_mat_copy)
    else:
        print(len(a))
        print(len(a[0]))
        print("First of all, you're stupid enough to try to invert a non-square matrix")
        print("Secondly, if you were not sure, why didn't you check the invertibility first?")
        print("Overconfident kids these days")
Patrick Artner
  • 50,409
  • 9
  • 43
  • 69
Anweshan Goswami
  • 75
  • 8
  • 3
    Most probably because you are comparing floating point values. .. 0.111111111111111112 != 0.111111111111111111 - and only 0.11111111 are shown. – Patrick Artner Oct 04 '21 at 05:21
  • 3
    The obvious thing to do here would be to print the difference between the two arrays, to see the exact magnitude of the difference. It's probably going to be something on the order of 1e-15; expecting exact equality of floats is unreasonable. This is exactly why numpy has an `.allclose()` method. – jasonharper Oct 04 '21 at 05:23
  • Try to compare the first elements manually and see the outcome. If they are not, try truncating them using *round()* or *floor()*. I'm no expert, but I can't help noticing that you are using floating point numbers that can sometimes be bitches. https://stackoverflow.com/questions/783897/how-to-truncate-float-values – Alex Oct 04 '21 at 05:26

1 Answers1

2

You are comparing floating point values. Only some of the actual digits are shown - and 0.11111111111111111111 is not the same as 0.11111111111111111112

Use numpy.allclose to compare them.

import numpy as np

a = np.array([[0.11111111, 0.37037037, -0.2962963],
              [ 0.88888889, -0.03703704, -0.37037037],
              [ -0.22222222, -0.07407407, 0.2592592600000]])

b = np.array([[0.11111111, 0.37037037, -0.2962963],
              [ 0.88888889, -0.03703704, -0.37037037],
              [ -0.22222222, -0.07407407, 0.2592592600001]])

print(a == b)

print(np.allclose(a,b))  # accepts delta for comparison, see documentation

Output:

# a == b
[[ True  True  True]
 [ True  True  True]
 [ True  True False]]

# allclose
True
Patrick Artner
  • 50,409
  • 9
  • 43
  • 69