## 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))

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.

Pages: 1 2

### 5 Responses to “Matrix Determinant And Inverse”

1. John Cowan said

No worse than addend, summand, subtrahend, minuend, multiplicand.

2. Rutger said

In Python, the determinant part.

```def determinant(matrix):
if len(matrix) == 2:
return (matrix*matrix-matrix*matrix)
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[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]]))
```
3. Milbrae said

I’ve found these two links with explanations and code. It’s quite a hassle for larger dimensions….

Determinant of a matrix: http://www.geeksforgeeks.org/determinant-of-a-matrix/

Though, I haven’t tried their code yet.

4. Paul said

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 [*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))
```
```inverse
0.1587  -0.7778   0.5079
0.0113   0.1587  -0.1066
-0.2245   0.8571  -0.4898
b*invb
1.0000  -0.0000   0.0000
0.0000   1.0000   0.0000
0.0000   0.0000   1.0000
determinant 440.99999999999994
```
5. Steve said

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

```