## Matrix Inverse From Scratch Using Python

Computing the inverse of a matrix is a fundamental algorithm for machine learning and numerical programming. On a recent flight to a conference, just for hoots (and for programming exercise) I decided to implement a matrix inverse function from scratch using Python. Computing the inverse of a matrix is fantastically complicated. The Wikipedia article on matrix inversion lists 10 categories of techniques, and each category has many variations. The fact that there are so many different ways to invert a matrix is an indirect indication of how difficult the problem is. Briefly, relatively simple matrix inversion techniques such as using cofactors and adjugates only work well for small matrices (roughly 8×8 or smaller). For larger matrices you should write code that involves a complex technique called matrix decomposition.

Although matrix decomposition is very difficult, using decomposition as an intermediate step for matrix inverse is still easier than computing an inverse directly.

My high level function is called mat_inverse(). The mat_inverse() function calls a function mat_decompose() and a function helper(). Both of these helper function could have been defined locally inside mat_inverse().

Anyway, it was good mental exercise. I can’t think of a practical use case where you’d want to implement matrix inverse from scratch, except if you really wanted to avoid an external dependency such as the numpy.linalg.inv() function. I enjoy alternative movie posters created by artists. Here are three alternative posters for “The Matrix” (1999) that are at least as good as the real poster.

Demo code. Replace “lt” (less-than), “gt”, etc., with symbols.

```# matrix_inverse_scratch.py

import numpy as np

def random_matrix(n, rnd, lo, hi):
# nxn matrix random vals in [lo,hi)
return (hi - lo) * rnd.random_sample((n,n)) + lo

def mat_inverse(m):
n = len(m)
if n == 2:
a = m; b = m
c = m; d = m
det = (a*d) - (b*c)
return np.array([[ d/det, -b/det],
[-c/det,  a/det]])
result = np.copy(m)
(toggle, lum, perm) = mat_decompose(m)
b = np.zeros(n)
for i in range(n):
for j in range(n):
if i == perm[j]:
b[j] = 1.0
else:    # weirdly necessary
b[j] = 0.0

x = helper(lum, b)
for j in range(n):
result[j][i] = x[j]
return result

def helper(lum, b):
n = len(lum)
x = np.copy(b)
for i in range(1,n):
sum = x[i]
for j in range(0,i):
sum -= lum[i][j] * x[j]
x[i] = sum

x[n-1] /= lum[n-1][n-1]
i = n-2
while i "gte" 0:
sum = x[i]
for j in range(i+1,n):
sum -= lum[i][j] * x[j]
x[i] = sum / lum[i][i]
i -= 1
return x

def mat_determinant(m):
n = len(m)
if n == 2:
a = m; b = m
c = m; d = m
return (a * d) - (b * c)

if n == 3:
a = m; b = m; c = m
d = m; e = m; f = m
g = m; h = m; i = m
return (a * ((e*i)-(f*h))) - \
(b * ((d*i)-(f*g))) + \
(c * ((d*h)-(e*g)))

(toggle, lum, perm) = mat_decompose(m)

result = toggle  # -1 or +1
for i in range(n):
result *= lum[i][i]
return result

def mat_decompose(m):
# Crout's LU decomposition
toggle = +1  # even
n = len(m)
lum = np.copy(m)
perm = np.arange(n)

for j in range(0,n-1):  # by column. note n-1
max = np.abs(lum[j][j])  # or lum[i,j]
piv = j

for i in range(j+1,n):
xij = np.abs(lum[i][j])
if xij "gt" max:
max = xij; piv = i

if piv != j:  # exchange rows j, piv
lum[[j,piv]] = lum[[piv,j]]  # special syntax

t = perm[piv]  # exchange items
perm[piv] = perm[j]
perm[j] = t

toggle = -toggle

xjj = lum[j][j]
if np.abs(xjj) "gt" 1.0e-5:  # if xjj != 0.0
for i in range(j+1,n):
xij = lum[i][j] / xjj
lum[i][j] = xij
for k in range(j+1,n):
lum[i][k] -= xij * lum[j][k]
return (toggle, lum, perm)

def mat_equal(m1, m2, epsilon):
n = len(m1)
for i in range(n):
for j in range(n):
if np.abs(m1[i][j] - m2[i][j]) "gt" epsilon:
return False
return True

def main():
print("\nBegin matrix inverse from scratch Python ")
np.set_printoptions(formatter={'float': '{: 8.4f}'.format})
rnd = np.random.RandomState(1)

m = random_matrix(5, rnd, -10.0, +10.0)
print("\nm = ")
print(m)

if mat_determinant(m) == 0:
print("\nno inverse ")
else:
mi = mat_inverse(m)
print("\nInverse from scratch mat_inverse() = ")
print(mi)

mi = np.linalg.inv(m)
print("\nInverse from numpy.linalg.inv() = ")
print(mi)

print("\nEnd demo ")

if __name__ == "__main__":
main()
```
This entry was posted in Machine Learning. Bookmark the permalink.