用Racket计算行列式

最基础的是拉普拉斯展开。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
; 计算行列式
(define (determinant matrix)
(define (remove-row-col m row col)
(for/list ([r (in-range (length m))]
#:when (not (= r row)))
(for/list ([c (in-range (length (m 0)))]
#:when (not (= c col)))
(list-ref (list-ref m r) c))))

(cond
[(= (length matrix) 1) (list-ref (list-ref matrix 0) 0)] ; 1x1 矩阵
[else
(for/sum ([i (in-range (length matrix))])
(let ([sign (if (even? i) 1 -1)])
(* sign
(list-ref (list-ref matrix 0) i)
(determinant (remove-row-col matrix 0 i)))))]))

(define example-matrix [[1 2 3] [0 1 4] [5 6 0]])
(determinant example-matrix) ; 应该返回 -1

但是这么写效率很差,可以使用高斯消元法。高斯消元法通过将矩阵转换为行最简形(或阶梯形),然后计算对角线元素的乘积来计算行列式。这种方法的复杂度远低于直接使用拉普拉斯展开。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

(define (swap-rows matrix row1 row2)
(let ([temp (list-ref matrix row1)])
(list-set matrix row1 (list-ref matrix row2))
(list-set matrix row2 temp)))

(define (scale-row matrix row factor)
(list-set matrix row (map (λ (x) (* x factor)) (list-ref matrix row))))

(define (add-scaled-row matrix src-row dest-row factor)
(list-set matrix dest-row (map (λ (x y) (+ x (* y factor)))
(list-ref matrix dest-row)
(list-ref matrix src-row))))

(define (gaussian-elimination matrix)
(let ([n (length matrix)])
(for/fold ([det 1]
[m matrix])
([i (in-range n)])
(match-let* ([(pivot-row (for/or ([r (in-range i n)])
(when (not (zero? (list-ref (list-ref m r) i)))
r)))
(pivot (list-ref (list-ref m pivot-row) i))]
(unless pivot-row (return (values 0 m)))
(when (not (= pivot-row i))
(set! det (- det))
(set! m (swap-rows m i pivot-row)))
(set! m (scale-row m i (/ 1 pivot)))
(for ([j (in-range (+ i 1) n)])
(set! m (add-scaled-row m i j (- (list-ref (list-ref m j) i)))))
(values (* det pivot) m)))))

(define (determinant matrix)
(let-values ([(det _) (gaussian-elimination matrix)])
det))

(define example-matrix [[1 2 3] [0 1 4] [5 6 0]])
(determinant example-matrix) ; 应该返回 -1