A few days ago I was exploring the ideas behind implementing matrix inversion from scratch using Python. There are dozens of matrix inversion algorithms but the one I usually use involves decomposing the source matrix into two other matrices. And there are several matrix decomposition algorithms; Crout’s algorithm is my usual choice. See https://jamesmccaffreyblog.com/2021/12/24/matrix-decomposition-using-python/.

One of the interesting side effects of matrix decomposition is that you can get the determinant of a matrix by multiplying the elements on the combined Lower and Upper decomposition matrix. So I decided to code up a determinant function.

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

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

  (toggle, lum, perm) = mat_decompose(m)
  print("\nThe combined LU decomposition matrix is: ")
  print(lum)
  print("The decomposition toggle is: %d " % toggle)

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

The function checks if the source matrix is 2×2 because decomposition requires 3×3 or greater. The function checks for 3×3 as a matter of efficiency because decomposition is very complex.

My demo generated a 4×4 matrix with random values between -10.0 and +10.0:

[[ -1.6596   4.4065  -9.9977  -3.9533]
 [ -7.0649  -8.1532  -6.2748  -3.0888]
 [ -2.0647   0.7763  -1.6161   3.7044]
 [ -5.9110   7.5623  -9.4522   3.4094]]

The emo computed the determinant using my from-scratch function and the library numpy.linalg.det() function and both were the same: -1553.31819174. My function displayed the computed LU decomposition matrix and toggle:

LU = 
[[ -7.0649  -8.1532  -6.2748  -3.0888]
 [  0.8367  14.3839  -4.2023   5.9936]
 [  0.2349   0.4395  -6.6768  -5.8620]
 [  0.2922   0.2196  -0.1708   2.2893]]
toggle = -1

The determinant was found by multiplying -1 * -7.0649 * 14.3839 * -6.6768 * 2.2893.

Good fun and a good daily Python programming exercise.

Two random images from an Internet search for “determinant”. The photo on the left makes sense — a video game under development. But the photo on the right — a flash mob stealing coats in San Francisco — seems to have nothing to do with “determinant”. Maybe some sort of correlation between race and crime, as the determinants of criminal behavior?

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

# matrix_det_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_determinant(m):
  n = len(m)
  if n == 2:
    a = m[0][0]; b = m[0][1]
    c = m[1][0]; d = m[1][1]
    return (a * d) - (b * c)

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

  (toggle, lum, perm) = mat_decompose(m)
  print("\nThe combined LU decomposition matrix is: ")
  print(lum)
  print("The decomposition toggle is: %d " % toggle)

  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 main():
  print("\nBegin matrix determinant from scratch Python ")
  np.set_printoptions(formatter={'float': '{: 8.4f}'.format})
  rnd = np.random.RandomState(1)

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

  d1 = mat_determinant(m)
  d2 = np.linalg.det(m)
  print("\nDeterminant from scratch   = %0.8f " % d1)
  print("Determinant from np.linalg = %0.8f " % d2)

  print("\nEnd demo ")

if __name__ == "__main__":
  main()