How to find the inverse of a matrix? (step by step solution)

Question

How to find the inverse of a matrix in Python by hand?

Reasons you might want to do this: I need to learn how to inverse a matrix for a school project. I'm interested in calculating the inverse of a matrix using a programming language.

Answer 1

I recently wrote a code to find the inverse of a matrix in Python. It gives a step by step explanation as you run the code. It also determines whether the inverse exists. I hope you enjoy it!

This code is for educational purposes. This might not be the most efficient way.

# Import packages
from numpy import *
from random import *

# Defining functions (skip this cell for readability)

def Gauss_explain():
    global A, A_inv
    for j in range(i+1,len(A)):
        if (A[j][i]):
            factor=A[j][i]/A[i][i]
            for k in range(len(A)):
                A[j][k]-=factor*A[i][k]
                A_inv[j][k]-=factor*A_inv[i][k]
            print()
            if (factor>0): print('--- Gauss elimination: row',j+1,'-',factor,'* row',i+1,'---\n')
            else: print('--- Gauss elimination: row',j+1,'+',-factor,'* row',i+1,'---\n')
            print(A)
            print()
            print(A_inv)
        
def Switch_explain():
    global A, A_inv
    det=0
    for j in range(i+1,len(A)):
        if (A[j][i]):
            temp=A[i].copy()
            A[i]=A[j]
            A[j]=temp
            temp=A_inv[i].copy()
            A_inv[i]=A_inv[j]
            A_inv[j]=temp
            det=1
            print()
            print('--- Switch rows',i+1,'and',j+1,'---\n')
            print(A)
            print()
            print(A_inv)
            break
    return det

def Gauss_up_explain():
    global A, A_inv
    if (A[j][i]):
        factor=A[j][i]/A[i][i]
        A[j][i]-=factor*A[i][i]
        for k in range(len(A)):
            A_inv[j][k]-=factor*A_inv[i][k]
        print()
        if (factor>0): print('--- Gauss elimination (up): row',j+1,'-',factor,'* row',i+1,'---\n')
        else: print('--- Gauss elimination (up): row',j+1,'+',-factor,'* row',i+1,'---\n')
        print(A)
        print()
        print(A_inv)
    
def Divide_explain():
    global A, A_inv
    if (A[i][i]!=1):
        factor=A[i][i]
        A[i][i]=1
        A_inv[i]/=factor
        print()
        print('--- Divide row',i+1,'by',factor,'---\n')
        print(A)
        print()
        print(A_inv)

# Choose a random seed to generate a matrix
seed(28)

# Interesting examples:
# Seed 4 is not invertable (det=0)
# Seed 28 uses the switch rows operator (zero value on diagonal)

# Generate a matrix A with random input
A_len=3
A=zeros((A_len,A_len))
for i in range(A_len):
    for j in range(A_len):
        A[i][j]=int(random()*11)
A_0=A.copy()
        
print('Matrix A:')
print(A)

# Generate the identity tensor (starting point for the inverse)
A_inv=eye(len(A))
print()
print('Inverse A (starting point):')
print(A_inv)

# Start solving for the inverse
# All operations are performed on the matrix A, as well as the identity tensor to find its inverse. 
# While det=1 the inverse exists, if det=0 the operation is aborted. (note: this function does not find the determinant)
det=1

for i in range(len(A)):
    
    # If the number on the diagonal is nonzero, apply Gauss elimination to triangualize the matrix. 
    if (A[i][i]): 
        Gauss_explain()
        
    # If the number on the diagonal is zero, check for nonzero terms below in the same column.
    # In that case switch the rows and perform the Gauss elimination nonetheless. 
    elif (Switch_explain()):
        Gauss_explain()
        
    # If all numbers below the diagonal are also zero, the determinant = 0, and the inverse doesn't exist.
    # The operation is aborted.
    else:
        det=0
        print()
        print('--- Det = 0, not invertable. ---')    
        break

if (det):
    # Now we know the inverse exists
    # We apply again Gauss elimination, this time to diagonalize the matrix.
    for i in range(len(A)-1,0,-1):
        for j in range(i-1,-1,-1):
            Gauss_up_explain()
    
    # Divide the rows of the matrix, so that we find the Identity tensor.
    # This results also in the inverse for the second matrix. 
    for i in range(len(A)):
        Divide_explain()

# Check if the matrix times its inverse gives the identity tensor
A_0@A_inv

The result might be:

Matrix A:
[[1. 1. 6.]
 [1. 1. 5.]
 [4. 2. 4.]]

Inverse A (starting point):
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

--- Gauss elimination: row 2 - 1.0 * row 1 ---

[[ 1.  1.  6.]
 [ 0.  0. -1.]
 [ 4.  2.  4.]]

[[ 1.  0.  0.]
 [-1.  1.  0.]
 [ 0.  0.  1.]]

--- Gauss elimination: row 3 - 4.0 * row 1 ---

[[  1.   1.   6.]
 [  0.   0.  -1.]
 [  0.  -2. -20.]]

[[ 1.  0.  0.]
 [-1.  1.  0.]
 [-4.  0.  1.]]

--- Switch rows 2 and 3 ---

[[  1.   1.   6.]
 [  0.  -2. -20.]
 [  0.   0.  -1.]]

[[ 1.  0.  0.]
 [-4.  0.  1.]
 [-1.  1.  0.]]

--- Gauss elimination (up): row 2 - 20.0 * row 3 ---

[[ 1.  1.  6.]
 [ 0. -2.  0.]
 [ 0.  0. -1.]]

[[  1.   0.   0.]
 [ 16. -20.   1.]
 [ -1.   1.   0.]]

--- Gauss elimination (up): row 1 + 6.0 * row 3 ---

[[ 1.  1.  0.]
 [ 0. -2.  0.]
 [ 0.  0. -1.]]

[[ -5.   6.   0.]
 [ 16. -20.   1.]
 [ -1.   1.   0.]]

--- Gauss elimination (up): row 1 + 0.5 * row 2 ---

[[ 1.  0.  0.]
 [ 0. -2.  0.]
 [ 0.  0. -1.]]

[[  3.   -4.    0.5]
 [ 16.  -20.    1. ]
 [ -1.    1.    0. ]]

--- Divide row 2 by -2.0 ---

[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0. -1.]]

[[ 3.  -4.   0.5]
 [-8.  10.  -0.5]
 [-1.   1.   0. ]]

--- Divide row 3 by -1.0 ---

[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

[[ 3.  -4.   0.5]
 [-8.  10.  -0.5]
 [ 1.  -1.  -0. ]]