Matrix Determinant And Inverse
May 23, 2017
We begin with a function that, given a matrix A and a cell at location Aij in the matrix, returns a newly-allocated matrix from which row i and column j have been deleted:
(define (sub-matrix a i j)
(let ((r (matrix-rows a)) (c (matrix-cols a)))
(let ((m (make-matrix (- r 1) (- c 1))) (new-i -1))
(for (old-i c)
(when (not (= old-i i))
(set! new-i (+ new-i 1))
(let ((new-j -1))
(for (old-j r)
(when (not (= old-j j))
(set! new-j (+ new-j 1))
(matrix-set! m new-i new-j
(matrix-ref a old-i old-j)))))))
m)))
That function doesn’t require A to be square, but the remaining functions do. Here we calculate the determinant by taking the first row of the matrix; the calculation is recursive, with a 2 × 2 matrix being the base case:
(define (matrix-determinant a) ; assume a is square
(let ((n (matrix-rows a)))
(if (= n 2)
(- (* (matrix-ref a 0 0) (matrix-ref a 1 1))
(* (matrix-ref a 1 0) (matrix-ref a 0 1)))
(let loop ((j 0) (k 1) (d 0))
(if (= j n) d
(loop (+ j 1) (* k -1)
(+ d (* k (matrix-ref a 0 j)
(matrix-determinant
(sub-matrix a 0 j))))))))))
The matrix of cofactors replaces each cell in a matrix with the determinant of the sub-matrix that eliminates the cell:
(define (matrix-cofactors a) ; assume a is square
(let* ((n (matrix-rows a)) (cof (make-matrix n n)))
(if (= n 2)
(for (i n)
(for (j n)
(matrix-set! cof i j
(* (expt -1 (+ i j))
(matrix-ref a (- 1 i) (- 1 j))))))
(for (i n)
(for (j n)
(matrix-set! cof i j
(* (expt -1 (+ i j))
(matrix-determinant (sub-matrix a i j)))))))
cof))
The adjugate is the transpose of the cofactor matrix:
(define (matrix-adjugate a) (matrix-transpose (matrix-cofactors a)))
Finally, the inverse is the scalar product of the inverse of the determinant times the adjugate:
(define (matrix-inverse a)
(matrix-scalar-multiply
(/ (matrix-determinant a))
(matrix-adjugate a)))
Who makes up all these words? Here’s an example showing that the product of a matrix and its inverse is the identity matrix:
> (define a '#( #(6 24 1) #(13 16 10) #(20 17 15))) > (matrix-multiply a (matrix-inverse a)) #(#(1 0 0) #(0 1 0) #(0 0 1))
We used the matrix datatype from the Standard Prelude and basic arithmetic on matrices from a previous exercise. You can run the program at http://ideone.com/nZGfo3.
Programming Praxis 
No worse than addend, summand, subtrahend, minuend, multiplicand.
In Python, the determinant part.
def determinant(matrix): if len(matrix) == 2: return (matrix[0][0]*matrix[1][1]-matrix[0][1]*matrix[1][0]) else: negator = 1 result = 0 for col in range(len(matrix)): sub_matrix = [[] for i in range(len(matrix)-1)] for sub_row in range(1,len(matrix)): for sub_col in range(len(matrix)): if sub_col != col: sub_matrix[sub_row-1].append(matrix[sub_row][sub_col]) result += negator * matrix[0][col] * determinant(sub_matrix) negator *= -1 return result print(-240 == determinant([[1,-4,9],[-6,7,3],[1,2,3]])) print(441 == determinant([[6,24,1],[13,16,10],[20,17,15]]))I’ve found these two links with explanations and code. It’s quite a hassle for larger dimensions….
Inverse of a matrix: http://www.geeksforgeeks.org/adjoint-inverse-matrix/
Determinant of a matrix: http://www.geeksforgeeks.org/determinant-of-a-matrix/
Though, I haven’t tried their code yet.
Here in Python an implementation of an LU decomposition. Once the matrix is decomposed, the determinant and inverse are straightforward.
from copy import deepcopy from functools import reduce from operator import mul def trans(A): return [list(a) for a in zip(*A)] def dot(a, b): return sum(ai*bi for ai, bi in zip(a, b)) def matvec(A, b): return [dot(a, b) for a in A] def matmul(A, B): return trans(matvec(A, b) for b in zip(*B)) def zeros(m, n): return [[0]*n for i in range(m)] def lu_decompose(A, tol): """ returns matrix LU with upper and lower matrix vector P: P[:n] contains the ordering of rows P[n] = n + numbeer of swaps in P[:n] """ LU = deepcopy(A) # the original matrix is not changed n = len(LU) P = list(range(n+1)) for i in range(n): imax = max(range(i, n), key=lambda m: abs(LU[m][i])) if abs(LU[i][imax]) < tol: raise ValueError("Matrix is degenerate") if imax != i: P[i], P[imax] = P[imax], P[i] LU[i], LU[imax] = LU[imax], LU[i] P[n] += 1 for j in range(i+1, n): LU[j][i] /= LU[i][i] for k in range(i+1, n): LU[j][k] -= LU[j][i] * LU[i][k]; return LU, P def lu_solve(LU, P, b): """solve Ax=b A is te original matrix""" n = len(LU) x = [b[p] for p in P[:-1]] for i in range(n): x[i] = b[P[i]] for k in range(i): x[i] -= LU[i][k] * x[k] for i in range(n-1, -1, -1): for k in range(i + 1, n): x[i] -= LU[i][k] * x[k] x[i] /= LU[i][i] return x def lu_invert(LU, P): """returns the inverse of the original matrix""" n = len(LU) IA = zeros(n, n) for j in range(n): for i in range(n): IA[i][j] = 1.0 if P[i] == j else 0.0 for k in range(i): IA[i][j] -= LU[i][k] * IA[k][j] for i in range(n-1, -1, -1): for k in range(i+1, n): IA[i][j] -= LU[i][k] * IA[k][j] IA[i][j] /= LU[i][i] return IA def lu_determinant(LU, P): """returns the determinant of the original matrix""" n = len(LU) det = reduce(mul, (LU[i][i] for i in range(n))) return det if (P[n]-n) % 2 == 0 else -det def print_matrix(A): print() for row in A: for ai in row: print("{:8.4f}".format(ai), end=" ") print() b = [[6, 24, 1], [13, 16, 10], [20, 17, 15]] A, P = lu_decompose(b, 1e-8) ib = lu_invert(A, P) print_matrix(ib) print_matrix(matmul(b, ib)) print(lu_determinant(A, P))Klong (for determinant)
steve@steve-Satellite-L555D:~$ rlwrap kg Klong 20161212 :" Beginning of code" det::{:[2=#x; cal1(x); cal2(x)]} cal1::{((x:@[0 0])*(x:@[1 1]))-((x:@[0 1])*(x:@[1 0]))} mults::{[a b]; a::1; b::[]; {b::b,(x*a); a::-a}'*x; b} makerow::{[a pos rn row t]; a::[]; t::x; rn::y; pos::z; row::t@rn; ((1+&#row):=0,pos){:[x=1; a::a,y; a]}'row; a} maketbl::{[pos tbl]; tbl::x; pos::y; {makerow(tbl; x; pos)}'(1+!(#tbl)-1)} maketbls::{[t]; t::x; {maketbl(t; x)}'!#t} cal2::{[a t]; a::0; t::x; mults(t){a::a+(x*det(y))}'maketbls(t); a} :" End of code" t2::[[1 2] [3 4]] t3::[[1 2 3] [2 3 4] [3 4 5]] t4::[[1 2 3 4] [2 3 4 5] [3 4 55 666] [7 88 999 1234]] :" Should be -2" det(t2) -2 :" Should be 0" det(t3) 0 :" Should be 498600" det(t4) 498600 :" Should be 18 per Wikipedia" det([[-2 2 -3] [-1 1 3] [2 0 -1]]) 18 :" Per @Rutger should be -240" det([[1 -4 9] [-6 7 3] [1 2 3]]) -240 :" Per @Rutger shuld be 441" det([[6 24 1] [13 16 10] [20 17 15]]) 441