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