Linear Algebra
SICP の問題 2.37 を基に図解でわかる線形代数に出てるナニをプログラムしてみる。
とりあえず
2.37 な実装が以下。
(define (accumulate op initial sequence)
(if (null? sequence)
initial
(op (car sequence)
(accumulate op initial (cdr sequence)))))
(define (accumulate-n op init sequence)
(if (null? (car sequence))
'()
(cons (accumulate op init (map car sequence))
(accumulate-n op init (map cdr sequence)))))
(define (dot-product v w)
(accumulate + 0 (map * v w)))
(define (m-*-v m v)
(map (lambda (row-m) (dot-product row-m v)) m))
(define (transpose mat)
(accumulate-n cons '() mat))
(define (m-*-m m n)
(let ((cols (transpose n)))
(map (lambda (row-m)
(map (lambda (row-cols)
(dot-product row-m row-cols))
cols))
m)))試験が以下。
(use gauche.test)
(add-load-path ".")
(load "2.37")
(test-start "dot-product")
(test-section "dot-product")
(test* "(dot-product '(0 9) '(2 0))"
0
(dot-product '(0 9) '(2 0)))
(test* "(dot-product '(2 0) '(2 0))"
4
(dot-product '(2 0) '(2 0)))
(test* "(dot-product '(1 1) '(3 1))"
4
(dot-product '(1 1) '(3 1)))
(test-section "m-*-v")
(test* "(m-*-v '((1 2 3) (4 5 6)) '(1 2 3))"
'(14 32)
(m-*-v '((1 2 3) (4 5 6)) '(1 2 3)))
(test* "(m-*-v '((1 2) (0 -1)) '(3 4))"
'(11 -4)
(m-*-v '((1 2) (0 -1)) '(3 4)))
(test-section "transpose")
(test* "(transpose '((1 2 3) (4 5 6) (7 8 9)))"
'((1 4 7) (2 5 8) (3 6 9))
(transpose '((1 2 3) (4 5 6) (7 8 9))))
(test-section "m-*-m")
(test* "(m-*-m '((1 3) (2 4)) '((4 2) (3 1)))"
'((13 5) (20 8))
(m-*-m '((1 3) (2 4)) '((4 2) (3 1))))
(test* "(m-*-m '((4 2) (3 1)) '((1 3) (2 4)))"
'((8 20) (5 13))
(m-*-m '((4 2) (3 1)) '((1 3) (2 4))))
(test* "(m-*-m '((1 2) (0 -1)) '((3) (4)))"
'((11) (-4))
(m-*-m '((1 2) (0 -1)) '((3) (4))))
(test* "(m-*-m '((2 1) (-1 3)) '((0 1) (1 -2)))"
'((1 0) (3 -7))
(m-*-m '((2 1) (-1 3)) '((0 1) (1 -2))))
(test-end)
まず trace
試験から書くか。
(test-section "trA")
(test* "(trA '((1 2 3) (4 5 6) (7 8 9)))"
15
(trA '((1 2 3) (4 5 6) (7 8 9))))これ、どーすりゃ良いのやら。R5RS 見たら list-ref って手続きがあるな。以下で試験パスしてるみたい。基本的に正方行列でないと error になるはずですが、これはこれでヨシって事に。
(define (trA l)
(let trA-inner ((idx 0) (rslt 0) (l l))
(if (null? l)
rslt
(trA-inner (+ 1 idx) (+ rslt (list-ref (car l) idx)) (cdr l))))
)
2 次正方行列の determinant
ええと試験は以下?
(test-section "detA-2")
(test* "(detA-2 '((1 2) (3 4)))"
-2
(detA-2 '((1 2) (3 4))))
(test* "(detA-2 '((2 3) (2 3)))"
0
(detA-2 '((2 3) (2 3))))実装は簡単ですな。でもこんなので良いのだろうか。
(define (detA-2 l)
(- (* (car (car l)) (cadr (cadr l)))
(* (car (cadr l)) (cadr (car l)))))
Inverse
逆行列。その前に行列のスカラ倍な手続きが必要ですな。なんて名前にするべきか。
(test-section "scalar-*-m")
(test* "(scalar-*-m 2 '((1 2) (3 4)))"
'((2 4) (6 8))
(scalar-*-m 2 '((1 2) (3 4))))実装が以下で一応試験パス。
(define (scalar-*-m s l)
(define (row-calc row)
(if (null? row)
'()
(cons (* s (car row)) (row-calc (cdr row)))))
(let s-*-m-inner ((l l) (rslt '()))
(if (null? l)
rslt
(s-*-m-inner (cdr l)
(append rslt (list (row-calc (car l))))))))なんか微妙だなぁ。こうしてみた。
(define (scalar-*-m s l)
(define (row-calc row)
(if (null? row)
'()
(cons (* s (car row)) (row-calc (cdr row)))))
(define (s-*-m-inner l rslt)
(if (null? l)
rslt
(s-*-m-inner (cdr l)
(append rslt (list (row-calc (car l)))))))
(s-*-m-inner l '()))そーゆー問題ではないな。row-calc が末尾呼び出しでないあたりが気持ち悪い。これ、map 使えば良いのか。ちょい再検討。まず、scalar-*-v を作るか。
(test-section "scalar-*-v")
(test* "(scalar-*-v 2 '(1 2))"
'(2 4)
(scalar-*-v 2 '(1 2)))map 使えば簡単。
(define (scalar-*-v s v) (map (lambda (x) (* s x)) v))
これを使う scalar-*-m を作れば良い。ここでも map が使えるはず。
(define (scalar-*-m s m) (map (lambda (x) (scalar-*-v s x)) m))
試験パス。いやはや。
で
これで逆行列の公式を適用する準備ができた訳ですな。detA が 0 じゃなければ云々、とテキストには書いてありますな。まず試験を準備。
(test-section "inverse")
(test* "(inverse '((2 3) (2 3)))"
#f
(inverse '((2 3) (2 3))))
(test* "(inverse '((2 5) (1 2)))"
'((-2 5) (1 -2))
(inverse '((2 5) (1 2))))
(test* "(inverse '((-2 6) (2 -4)))"
'((1 3/2) (1/2 1/2))
(inverse '((-2 6) (2 -4))))detA が 0 なら #f を戻す方向。実装は以下?
(define (inverse m)
(let ((detA (detA-2 m)))
(if (= 0 detA)
#f
(scalar-*-m (/ 1 detA) m))))試験パスしない。定義を良く見たらスカラ倍するマトリクスにも変更があるな。これはどうしたものやら。2 次の正方行列限定で無理矢理以下?
(define (inverse-2 m)
(let ((detA (detA-2 m))
(lst (list (list (cadr (cadr m)) (- (cadr (car m))))
(list (- (car (cadr m))) (car (car m))))))
(if (= 0 detA)
#f
(scalar-*-m (/ 1 detA) lst))))これは酷い。試験にはパスしてますが。とりあえずこの時点でエントリ投入。
