One of the fundamental algorithms in machine learning and numerical programming is matrix decomposition. The idea is to break a source matrix M into two matrices A and B so that A * B = M. Matrix decomposition doesn’t have any direct use, but decomposition is the key component of matrix inversion and several other essential algorithms.

One rainy weekend in the Pacific Northwest where I live, (yes, the “rainy” in the PNW is redundant) I decided to implement a matrix decomposition function using Python. There are several decomposition algorithm. I used Crout’s algorithm.

My demo program started out with matrix M:

[3.0  7.0  2.0  5.0] 
[1.0  8.0  4.0  2.0] 
[2.0  1.0  9.0  3.0] 
[5.0  4.0  7.0  1.0]])

The resulting LUM (“lower-upper”) decomposition and auxiliary permutation array are:

[ 5.0000  4.0000  7.0000  1.0000]
[ 0.2000  7.2000  2.6000  1.8000]
[ 0.4000 -0.0833  6.4167  2.7500]
[ 0.6000  0.6389 -0.6017  4.9048]

[3 1 2 0]

The lower part of LUM (elements below the diagonal with dummy 1s on the diagonal) is the first part of the decomposition:

[ 1.0000  0.0000  0.0000  0.0000]
[ 0.2000  1.0000  0.0000  0.0000]
[ 0.4000 -0.0833  1.0000  0.0000]
[ 0.6000  0.6389 -0.6017  1.0000]

The upper part of LUM is:

[ 5.0000  4.0000  7.0000  1.0000]
[ 0.0000  7.2000  2.6000  1.8000]
[ 0.0000  0.0000  6.4167  2.7500]
[ 0.0000  0.0000  0.0000  4.9048]

When you multiply lower times upper, the result is almost the original M, except some of the rows are permuted:

[ 5.0000  4.0000  7.0000  1.0000]
[ 1.0000  8.0000  4.0000  2.0000]
[ 2.0000  1.0000  9.0000  3.0000]
[ 3.0000  7.0000  2.0000  5.0000]

To recover the original matrix you can rearrange rows according to the information in the permutation array [3 1 2 0].

When I implemented Crout’s decomposition, I ran into a nasty syntax issue when swapping two rows of a matrix. It took me an hour of debugging. The details could fill an entire blog post so instead I’ll just note that swapping two rows of a matrix requires care.

Decomposition in art. Left: By artist Cane Dojcilovic. Center: By artist Alberto Seveso. Right: By artist Antonio Saraiva.

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

# crout.py
# crout's decompossition

import numpy as np

def mat_decompose(m):
  # Crout's LU decomposition for matrix determinant and inverse
  # stores combined lower & upper in lum[][]
  # stores row permuations into perm[]
  # toggle is +1 (even) or -1 number of row permutations
  # lower gets dummy 1.0s on diagonal (0.0s above)
  # upper gets lum values on diagonal (0.0s below)

  toggle = +1  # even
  n = len(m)
  lum = np.copy(m)
  perm = np.arange(n)

  for j in range(0,n-1):  # process 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 get_lower(lum):
  n = len(lum)
  result = np.zeros((n,n))
  for i in range(n):
    for j in range(n):
      if i == j:
        result[i][j] = 1.0
      elif i "gt" j:
        result[i][j] = lum[i][j]
  return result

def get_upper(lum):
  n = len(lum)
  result = np.zeros((n,n))
  for i in range(n):
    for j in range(n):
      if i "lte" j:
        result[i][j] = lum[i][j]
  return result

def unwind(lu, perm):
  result = np.copy(lu)
  for i in range(len(perm)):
    j = perm[i]
    result[i] = lu[j]
  return result
    

def main():
  print("\nBegin Crout's matrix decomposition demo ")
  np.set_printoptions(formatter={'float': '{: 0.4f}'.format})
 
  m = np.array([[3.0, 7.0, 2.0, 5.0],
                [1.0, 8.0, 4.0, 2.0],
                [2.0, 1.0, 9.0, 3.0],
                [5.0, 4.0, 7.0, 1.0]])

  print("\nm = ")
  print(m)

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

  print("\nlum = ")
  print(lum)
  print("\nperm = ")
  print(perm)
  
  lower = get_lower(lum)
  print("\nlower = ")
  print(lower)

  upper = get_upper(lum)
  print("\nupper = ")
  print(upper)

  lu = np.matmul(lower, upper)
  print("\nlower * upper = ")
  print(lu)

  original = unwind(lu, perm)
  print("\noriginal = ")
  print(original)

  print("\nEnd demo ")

if __name__ == "__main__":
  main()