[plt-scheme] Matrix Determinants
So, when I get bored doing math assignments, I write short programs to
double-check my answers. It's somewhat of a bad habit actually, as it
often ends up taking long enough that I don't have enough time to
finish the assignment...
Anyways, during my last assignment I came up with the following for
taking matrix determinants and multiplying matrices. Matrices are
lists of rows, rows are lists of numbers, all rows of a matrix are the
same length. I wasn't really worrying about numerical error, since the
matrices I gave it consisted of exact numbers. I'm pretty happy with
how simple matrix-mult is ((apply map list m) is neat!), but taking
the determinant of a matrix proved messy. I add multiples of the first
row to the rest to zero out the rest of the first column, so I only
need to do ~n recursive calls. It gets messy when I need to deal with
the case that the first row has a zero as the first element, swapping
list elements is tricky!
Can anyone come up with a simple matrix determinant implementation?
(other than the obvious O(n!) one of course)
Henk
(module math mzscheme
(require (lib "list.ss" "srfi" "1"))
(provide det matrix-mult)
(define (find-index pred l)
(let loop ((index 0) (l l))
(cond ((null? l) #f)
((pred (car l)) index)
(else (loop (add1 index) (cdr l))))))
(define (swap l i1 i2)
(let loop ((l l) (i1 (min i1 i2)) (i2 (max i1 i2)))
(if (= 0 i1)
(cons (list-ref l i2)
(replace-index (cdr l) (sub1 i2) (car l)))
(cons (car l)
(loop (cdr l) (sub1 i1) (sub1 i2))))))
(define (replace-index l i v)
(if (= i 0)
(cons v (cdr l))
(cons (car l) (replace-index (cdr l) (add1 i) v))))
; O(n^3) where n is the height of m
(define (det m)
(let ((index (find-index (lambda (r) (not (= 0 (car r)))) m)))
(if (not index)
0
(let ((swap-multiplier (if (= index 0) 1 -1))
(m-swapped (if (= index 0) m (swap m 0 index))))
(let ((submatrix (gaussian-submatrix m-swapped)))
(if (null? submatrix)
(caar m-swapped)
(* swap-multiplier
(caar m-swapped)
(det submatrix))))))))
; Zeros entries under the top-left entry, by adding multiples of the first
; row. Returns the submatrix ommitting the first row and column. The top left
; entry must be nonzero. I couldn't really think of a good name, "gaussian
; submatrix" doesn't actually mean anything as far as I know.
; O(w*h) where w and h are the width and height of m
(define (gaussian-submatrix m)
(let ((first-row (car m)))
(map (lambda (row)
(let ((mult (/ (car row)
(car first-row))))
(map (lambda (above current)
(- current (* mult above)))
(cdr first-row)
(cdr row))))
(cdr m))))
(define (matrix-mult x y)
(let ((y-cols (apply map list y)))
(map (lambda (x-row)
(map (lambda (y-col) (fold + 0 (map * x-row y-col)))
y-cols))
x)))
) ; end module
Examples:
> (det '((1 2 3) (4 5 6) (7 8 9)))
0
> (det '((1)))
1
> (det '((10 20) (30 41)))
-190