Find inverse matrix using gauss jordan elimination python

To inverse square matrix of order n using Gauss Jordan Elimination, we first augment input matrix of size n x n by Identity Matrix of size n x n.

After augmentation, row operation is carried out according to Gauss Jordan Elimination to transform first n x n part of n x 2n augmented matrix to identity matrix.

Matrix Inverse Using Gauss Jordan Python Program


# Importing NumPy Library
import numpy as np
import sys

# Reading order of matrix
n = int(input('Enter order of matrix: '))

# Making numpy array of n x 2n size and initializing 
# to zero for storing augmented matrix
a = np.zeros((n,2*n))

# Reading matrix coefficients
print('Enter Matrix Coefficients:')
for i in range(n):
    for j in range(n):
        a[i][j] = float(input( 'a['+str(i)+']['+ str(j)+']='))

# Augmenting Identity Matrix of Order n
for i in range(n):        
    for j in range(n):
        if i == j:
            a[i][j+n] = 1

# Applying Guass Jordan Elimination
for i in range(n):
    if a[i][i] == 0.0:
        sys.exit('Divide by zero detected!')
        
    for j in range(n):
        if i != j:
            ratio = a[j][i]/a[i][i]

            for k in range(2*n):
                a[j][k] = a[j][k] - ratio * a[i][k]

# Row operation to make principal diagonal element to 1
for i in range(n):
    divisor = a[i][i]
    for j in range(2*n):
        a[i][j] = a[i][j]/divisor

# Displaying Inverse Matrix
print('\nINVERSE MATRIX IS:')
for i in range(n):
    for j in range(n, 2*n):
        print(a[i][j], end='\t')
    print()

Output

Enter order of matrix: 3
Enter Matrix Coefficients:
a[0][0]=1
a[0][1]=1
a[0][2]=3
a[1][0]=1
a[1][1]=3
a[1][2]=-3
a[2][0]=-2
a[2][1]=-4
a[2][2]=-4

INVERSE MATRIX IS:
 3.0	 1.0	 1.5	
-1.25	-0.25	-0.75	
-0.25	-0.25	-0.25	

Recommended Readings

  1. Matrix Inverse Using Gauss Jordan Method Algorithm
  2. Matrix Inverse Using Gauss Jordan Method Pseudocode
  3. Matrix Inverse Using Gauss Jordan C Program
  4. Matrix Inverse Using Gauss Jordan C++ Program

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!

How do you find the inverse of a matrix in python?

Python provides a very easy method to calculate the inverse of a matrix. The function numpy. linalg. inv() which is available in the python NumPy module is used to compute the inverse of a matrix.