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.

Advertisements

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[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]]))
    
  3. Milbrae said

    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.

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: