I was sitting at home one Sunday morning. I like to write code almost every day, to entertain myself, and to keep in practice. I realized it had been many months since I implemented a function to compute the inverse of a matrix using Gauss-Jordan elimination to reduced row echelon form.

The Wikipedia article on matrix inverse lists 14 different algorithms, and each algorithm has multiple variations, and each variation has dozens of implementation alternatives. Six examples include 1.) LUP decomposition algorithm (Crout version and others), 2.) QR decomposition algorithm (Householder version and others), 3.) SVD decomposition algorithm (Jacobi version and others), 4.) Cholesky decomposition (Banachiewicz version and others), 5.) Cayley-Hamilton algorithm (Faddeev–LeVerrier version and others), 6.) Newton Iteration. The fact that there are so many different techniques is a direct reflection of how tricky it is to compute a matrix inverse.

The Gauss-Jordan technique is arguably the one of the simplest ways to compute a matrix inverse but it is rarely used in practice because, compared to other techniques, Gauss-Jordan is susceptible to numerical instability, especially when the source matrix is nearly singular (meaning has no inverse).

Anyway, after a few minutes of coding, I had a demo up and running. The output is:

Begin matrix inverse using Gauss-Jordan elimination
  to reduced row echelon form

Source matrix:
[[3. 5. 9.]
 [1. 1. 5.]
 [2. 4. 8.]]

Computing inverse
Done

Inverse from Gauss-Jordan:
[[ 1.5000  0.5000 -2.0000]
 [-0.2500 -0.7500  0.7500]
 [-0.2500  0.2500  0.2500]]

Minv * M:
[[ 1.0000 -0.0000 -0.0000]
 [ 0.0000  1.0000  0.0000]
 [ 0.0000  0.0000  1.0000]]

M * Minv:
[[ 1.0000 -0.0000  0.0000]
 [-0.0000  1.0000  0.0000]
 [-0.0000 -0.0000  1.0000]]

End demo

Notice that some of the values are “-0.0000” which means a negative number close to 0, not some sort of weird value.

The Gauss-Jordan technique starts by creating a matrix that is the source matrix that is augmented by the Identity matrix to the right. For the demo matrix, the augmented matrix is:

 3  5  9   1  0  0
 1  1  5   0  1  0
 2  4  8   0  0  1

Then the algorithm iterates through each row and performs three possible row operations:

* exchange two rows
* multiply a row by a constant
* multiply a row by a constant times another row

The version I implemented doesn’t exchange rows, but some versions do. The goal is to transform the augmented matrix so that the left side of the augmented matrix becomes the identity matrix. For the demo, the result is:

 1  0  0    1.50   0.50  -2.00
 0  1  0   -0.25  -0.75   0.75  
 0  0  1   -0.25   0.25   0.25

When finished, the right side of the transformed augmented matrix is the inverse of the source matrix. Clever.

OK, good fun for me. But as mentioned, the Gauss-Jordan algorithm to compute a matrix inverse is not used very often, so this was essentially mental exercise for me.

I don’t know much about art or fashion. Wait, that’s not entirely correct. More accurate is: I know almost nothing about art or fashion. I did a little bit of Internet searching for “inverse black and white fashion photography” and fell down a rabbit hole of art information that I really didn’t understand.

I have no idea if these images are good or not. I prefer dealing with math and computer science where it’s easier to evaluate things.

Demo code:

# matrix_inverse_gauss_jordan.py

import numpy as np

# -----------------------------------------------------------

def mat_inverse_gauss_jordan(A):
  n = len(A)  # assumed square
  B = np.concatenate((A, np.eye(n)), axis=1)  # augmented

  for i in range(n):
    # make the diagonal element 1
    pv = B[i, i]
    if np.abs(pv) "lt" 1.0e-12:  # replace with symbol
      return None  # singular matrix
    B[i, :] = B[i, :] / pv

    # eliminate other rows
    for j in range(n):
      if i != j:
        factor = B[j, i]  # note indices
        B[j, :] = B[j, :] - factor * B[i, :]

  # extract the inverse matrix from right side of augmented
  result = B[:, n:]
  return result

# -----------------------------------------------------------

print("\nBegin matrix inverse using Gauss-Jordan " + 
  "elimination to reduced row echelon form ")

M = np.array([[3, 5, 9],
              [1, 1, 5],
              [2, 4, 8]], dtype=np.float64)
print("\nSource matrix: ")
print(M)

print("\nComputing inverse ")
Minv = mat_inverse_gauss_jordan(M)
print("Done ")

np.set_printoptions(precision=4, suppress=True,
  floatmode='fixed')
print("\nInverse from Gauss-Jordan: ")
print(Minv)

Minv_M = np.matmul(Minv, M)
print("\nMinv * M: ")
print(Minv_M)

M_Minv = np.matmul(M, Minv)
print("\nM * Minv: ")
print(M_Minv)

print("\nEnd demo ")