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://jamesmccaffrey.wordpress.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()

You must be logged in to post a comment.