Finding inverse of a matrix using Gauss-Jordan Elimination in Python

Question

So I am trying to find inverse of a matrix (using Python lists) by Gauss-Jordan Elimination. But I am facing this peculiar problem. In the code below, I apply my code to the given matrix and it reduces to the identity matrix as intended.

M = [[0, 2, 1], [4, 0, 1], [-1, 2, 0]]
P = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
n = len(P)
def inverse(a):
    for k in range(n):
        if abs(a[k][k]) < 1.0e-12:
            for i in range(k+1, n):
                if abs(a[i][k]) > abs(a[k][k]):
                    for j in range(k, n):
                        a[k][j], a[i][j] = a[i][j], a[k][j]
                    break
        pivot = a[k][k]
        for j in range(k, n):
            a[k][j] /= pivot
        for i in range(n):
            if i == k or a[i][k] == 0: continue
            factor = a[i][k]
            for j in range(k, n):
                a[i][j] -= factor * a[k][j]
    return a

inverse(M)

The output is

[[1.0, 0.0, 0.0], [0, 1.0, 0.0], [0.0, 0.0, 1.0]]

But when I apply the same code, after adding lines of code for my identity matrix (which is part of the augmented matrix with the given matrix), It is not giving me the correct inverse when it should (as I am applying same operation on it as I'm applying on the given matrix).

M = [[0, 2, 1], [4, 0, 1], [-1, 2, 0]]
P = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
n = len(P)
def inverse(a, b):
    for k in range(n):
        if abs(a[k][k]) < 1.0e-12:
            for i in range(k+1, n):
                if abs(a[i][k]) > abs(a[k][k]):
                    for j in range(k, n):
                        a[k][j], a[i][j] = a[i][j], a[k][j]
                        b[k][j], b[i][j] = b[i][j], b[k][j]
                    b[k], b[i] = b[i], b[k]
                    break
        pivot = a[k][k]
        for j in range(k, n):
            a[k][j] /= pivot
            b[k][j] /= pivot
        for i in range(n):
            if i == k or a[i][k] == 0: continue
            factor = a[i][k]
            for j in range(k, n):
                a[i][j] -= factor * a[k][j]
                b[i][j] -= factor * b[k][j]
    return a, b

inverse(M, P)

The output is not the inverse matrix, but something else (though the last column has correct entries).

([[1.0, 0.0, 0.0], [0, 1.0, 0.0], [0.0, 0.0, 1.0]],
 [[0.0, 0.25, 0.3333333333333333],
  [1, 0.0, 0.6666666666666666],
  [0.0, 0.25, -1.3333333333333333]])

I tried using print statements while debugging, and I think the line where I divide the rows by the pivot, has some issue. It works fine in the original matrix but doesn't work with the identity matrix. Also, note that only the last column of the identity matrix gets converted to the correct entries of the inverse matrix.

The correct inverse matrix for reference is

I = [[-1/3, 1/3, 1/3], [-1/6, 1/6, 2/3], [4/3, -1/3, -4/3]]

Any help would be much appreciated. Thanks in advance!

Answer 1

This is my implementation for gaussian elimination. Forgive my bad variable assignment:

import numpy as np
def matrix_inverter(matrix):
    if len(matrix[0])==len(matrix):
        order= len(matrix)
        unity= np.zeros([order,order])
        for l in range(order):
            unity[l,l]= 1.0
        for k in range(order):
            pivot= matrix[k,k]
            for j in range(order):
                if j==k:
                    continue
                else:
                    sec= matrix[j,k]
                    mult= sec/pivot
                    for i in range(order):
                        matrix[j,i] = matrix[j,i] - mult*matrix[k,i]
                        unity[j,i] =  unity[j,i] - mult*unity[k,i]
            matrix[k]= matrix[k]/pivot
            unity[k]= unity[k]/pivot
        return unity
    else:
        return 'invalid order'

Answer 2

I could not figure out the problem in my code though I found a workaround and I'll be instead posting that. Hope that might help someone else.

Instead of doing the same operations on the two matrices of which I formed the augmented matrix, I decided to deal it as a n x 2n matrix. Here is the fully working code:

def inverse(a):
    n = len(a) #defining the range through which loops will run
    #constructing the n X 2n augmented matrix
    P = [[0.0 for i in range(len(a))] for j in range(len(a))]
    for i in range(3):
        for j in range(3):
            P[j][j] = 1.0
    for i in range(len(a)):
        a[i].extend(P[i])
    #main loop for gaussian elimination begins here
    for k in range(n):
        if abs(a[k][k]) < 1.0e-12:
            for i in range(k+1, n):
                if abs(a[i][k]) > abs(a[k][k]):
                    for j in range(k, 2*n):
                        a[k][j], a[i][j] = a[i][j], a[k][j] #swapping of rows
                    break
        pivot = a[k][k] #defining the pivot
        if pivot == 0: #checking if matrix is invertible
            print("This matrix is not invertible.")
            return
        else:
            for j in range(k, 2*n): #index of columns of the pivot row
                a[k][j] /= pivot
            for i in range(n): #index the subtracted rows
                if i == k or a[i][k] == 0: continue
                factor = a[i][k]
                for j in range(k, 2*n): #index the columns for subtraction
                    a[i][j] -= factor * a[k][j]
    for i in range(len(a)): #displaying the matrix
        for j in range(n, len(a[0])):
            print(a[i][j], end = " ")
        print()