Matrix Decomposition Using Python

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()
This entry was posted in Machine Learning. Bookmark the permalink.