数値計算,洲之内著,石渡改訂,サイエンス社 を読む (2)
9.1節「対称行列の固有値」を読み、ヤコビ法を実装した。
ヤコビ法
実対称行列を対角化するアルゴリズム。
ギブンスの回転行列Pを利用して実対称行列Aの非対角成分を0に近づけていく。
この操作を繰り返し、行列Aを対角化する。（詳しい説明は省略）
ヤコビ法をschemeで実装し、以下の例題9.1でプログラムの動作をテストした。
固有値を近似計算できている。ちゃんと動作しているっぽい。
（対角化された行列が微妙に対称行列になってない。非対角成分のわずかな誤差。おそらく行列の積の計算時に発生した誤差が原因）
例題9.1
解答:
先に結果を書き、後でソースコードを書く。
固有値は対角成分
((3.005096954913813 9.346597978251512e-7 9.898932550527318e-5) (9.346597978216164e-7 1.9999802019410715 3.469446951953614e-18) (9.898932550526628e-5 0.0 0.994922843145116))
対応する固有ベクトルを並べた行列
((0.998704179891266 -0.010888886258671288 -0.04971310917413542) (0.010431980378144533 0.9999010148892374 -0.009441091508041592) (0.04981099128810623 0.008910251372375937 0.9987189156951907))
(define (my-map-item2 op a b) (define (my-map-item2-p op a b tmp) (if (and (null? a) (null? b)) tmp (my-map-item2-p op (cdr a) (cdr b) (append tmp (list (op (car a) (car b))))))) (my-map-item2-p op a b '() )) (define (l-accum2 op1 op2 ini a b) (define (l-accum2-p op1 op2 ini a b) (if (and (null? a) (null? b)) ini (l-accum2 op1 op2 (op1 ini (op2 (car a) (car b))) (cdr a) (cdr b)))) (l-accum2-p op1 op2 ini a b)) (define (do-n-times op ini n) (define (do-n-times-p op ini n) (if (= n 0) ini (do-n-times-p op (op ini) (- n 1)))) (do-n-times-p op ini n)) (define (zero-line n) (do-n-times (lambda (x) (append x '(0))) '() n)) (define (zero-mat m n) (do-n-times (lambda (x) (append x (list (zero-line n)))) '() m)) (define (do-n-times-with-filter-counter op ini n filter? op-c counter) (define (do-n-times-with-filter-counter-p op ini n filter? op-c counter) (if (= n 0) ini (do-n-times-with-filter-counter-p op (op ini filter? counter) (- n 1) filter? op-c (op-c counter)))) (do-n-times-with-filter-counter-p op ini n filter? op-c counter)) (define (eye-line n loc1) (do-n-times-with-filter-counter (lambda (ini filter? counter) (if (filter? counter) (append ini '(1)) (append ini '(0)))) '() n (lambda (counter) (= counter 1)) (lambda (counter) (- counter 1)) loc1)) (define (eye-mat n) (do-n-times-with-filter-counter (lambda (ini filter counter) (append ini (list (filter counter)))) '() n (lambda (counter) (eye-line n counter)) (lambda (counter) (+ counter 1)) 1)) ;;和 (define (add-mat a1 a2) (my-map-item2 (lambda (s t) (my-map-item2 + s t)) a1 a2)) ;;スカラー倍 (define (scalar-mat a AA) (define (scalar-each a a1) (map (lambda (x) (* a x)) a1)) (map (lambda (x) (scalar-each a x)) AA)) ;;引き算 (define (sub-mat A1 A2) (add-mat A1 (scalar-mat -1 A2))) (define (sum-of-mul-each a1 a2) (l-accum2 + * 0 a1 a2)) ;;転置行列 (define (transposed-mat a) (apply map list a)) ;;行列の積 (define (mul-mat A1 A2) (define (mul-oneline-mat oneline A2) (map (lambda (x) (sum-of-mul-each oneline x)) (transposed-mat A2))) (define (mul-mat-pre A1 A2 tmp) (if (null? A1) tmp (mul-mat-pre (cdr A1) A2 (append tmp (list (mul-oneline-mat (car A1) A2)))))) (mul-mat-pre A1 A2 '() )) ;;行列のサイズ ;;(行数. 列数) (define (size-mat A) (define (N-mat A j) (if (null? (car A)) j (N-mat (list (cdr (car A))) (+ j 1)))) (define (M-mat A i) (if (null? A) i (M-mat (cdr A) (+ i 1)))) (cons (M-mat A 0) (N-mat A 0))) ;;行列の左右連結 (define (connect-mat-side A B) (define (connect-mat-side-p A B tmp) (if (null? A) tmp (connect-mat-side-p (cdr A) (cdr B) (append tmp (list (append (car A) (car B))))))) (if (= (car (size-mat A)) (car (size-mat B))) (connect-mat-side-p A B '()) (error "CONNECT-MAT-SIDE: (connect-mat-side A B) WRONG SIZE MATRIXES !"))) ;;行列から行列を抽出 (define (extract-mat A i j k l) (define (extract-from-list x1 x2 list-a) (define (extract-from-list-p x1 x2 list-a tmp count) (cond ((> count x2) tmp) ((and (>= count x1) (<= count x2)) (extract-from-list-p x1 x2 (cdr list-a) (append tmp (list (car list-a))) (+ count 1))) ((< count x1) (extract-from-list-p x1 x2 (cdr list-a) tmp (+ count 1))))) (if (and (> x1 0) (>= x2 x1) (>= (length list-a) x2)) (extract-from-list-p x1 x2 list-a '() 1) (error "EXTRACT-FROM-LIST: (extract-from-list x1 x2 list-a) WRONG X1 X2!"))) (define (extract-mat-p D k l tmp) (if (null? D) tmp (extract-mat-p (cdr D) k l (append tmp (list (extract-from-list k l (car D))))))) (extract-mat-p (extract-from-list i j A) k l '())) (define (cdr*n x n) (define (cdr*n-p x n) (if (= n 0) x (cdr*n-p (cdr x) (- n 1)))) (cdr*n-p x n)) ;;nondiag-biggest-abs ;;入力:対称行列A ;;非対角成分の絶対値が最大となる成分をa_{ij}としたとき ;;出力:(a_{ij} i . j) ;;ただし上三角だけ調べる (define (nondiag-biggest-abs A) ;;内部定義 (define (lbai l r) (define (lbai-p l tmp r c) (if (null? l) tmp (lbai-p (cdr l) (if (> (abs (car l)) (abs (car tmp))) (cons (car l) (cons r c)) tmp) r (+ c 1)))) (lbai-p l (cons 0 0) r (+ r 1))) (define (nondiag-biggest-abs-p A tmp row) (if (= (car (size-mat A)) 1) tmp (nondiag-biggest-abs-p (cdr A) (if (> (abs (car (lbai (cdr*n (car A) row) row))) (abs (car tmp))) (lbai (cdr*n (car A) row) row) tmp) (+ row 1)))) ; (if (equal? A (transposed-mat A)) (nondiag-biggest-abs-p A (cons 0 "all-nondiag") 1)) ; (error "nondiag-biggest-abs: input matrix should be a symmetric matrix"))) ;;ヤコビ法 ;;入力: A:実対称行列, tolerance:非対角成分の許容誤差（正の実数） ;;出力:(対角行列 . 対応する固有ベクトルを並べた行列) (define (jacobi-method A tolerance) ;;内部定義 (define (norm-plane a b) (sqrt (+ (* a a) (* b b)))) ;;ギブンスの回転行列(準備) ;;入力:行列A ;;出力: '(k m cos(psi) sin(psi)) ;;ただしk,m は、いずれも回転する行と列 (define (Pdat-from-mat A) (let ((nba (nondiag-biggest-abs A))) (let ((k (cadr nba)) (m (cddr nba)) (akm (car nba))) (let ((akk (caar (extract-mat A k k k k))) (amm (caar (extract-mat A m m m m)))) (let ((s (/ (- akk amm) 2.0))) (if (>= s 0) (let ((c-psi (sqrt (/ (+ 1 (/ s (norm-plane s akm))) 2)))) (let ((s-psi (/ akm (* 2 (norm-plane s akm) c-psi)))) (list k m c-psi s-psi))) (let ((c-psi (sqrt (/ (- 1 (/ s (norm-plane s akm))) 2)))) (let ((s-psi (/ (* -1.0 akm) (* 2 (norm-plane s akm) c-psi)))) (list k m c-psi s-psi))))))))) ;;入力: '(k m cos(psi) sin(psi)), 行列A ;;出力: ギブンスの回転行列 (define (make-P-from-Pdat-A Pdat A) (let ((k (car Pdat)) (m (car (cdr*n Pdat 1))) (c-psi (car (cdr*n Pdat 2))) (s-psi (car (cdr*n Pdat 3))) (n (car (size-mat A)))) (let ((m1 (- k 1)) (m2 (- m k 1)) (m3 (- n m))) (let ((els (lambda (lines counter-ini) (do-n-times-with-filter-counter (lambda (ini filter counter) (append ini (list (filter counter)))) '() lines (lambda (counter) (eye-line n counter)) (lambda (counter) (+ counter 1)) counter-ini)))) (let ((em1 (els m1 1)) (em2 (els m2 (+ k 1))) (em3 (els m3 (+ m 1)))) (let ((linek (map + (map (lambda (x) (* c-psi x)) (eye-line n k)) (map (lambda (x) (* -1 s-psi x)) (eye-line n m)))) (linem (map + (map (lambda (x) (* s-psi x)) (eye-line n k)) (map (lambda (x) (* c-psi x)) (eye-line n m))))) (append em1 (list linek) em2 (list linem) em3))))))) (define (jacobi-method-p A tolerance tmpA tmpP) (if (< (abs (car (nondiag-biggest-abs tmpA))) tolerance) (cons tmpA tmpP) ; (begin (display tmpA) (newline) ;;デバッグ用 A行列(tmpA)を各ステップで表示 (jacobi-method-p A tolerance (mul-mat (mul-mat (transposed-mat (make-P-from-Pdat-A (Pdat-from-mat tmpA) tmpA)) tmpA) (make-P-from-Pdat-A (Pdat-from-mat tmpA) tmpA)) (mul-mat tmpP (make-P-from-Pdat-A (Pdat-from-mat tmpA) tmpA)))));) ;;displayデバック時 ;;ここから本体 ヤコビ法 (if (equal? (transposed-mat A) A) (jacobi-method-p A tolerance A (eye-mat (car (size-mat A)))) (error "jacobi-method: (jacobi-method A tolerance) A should be a symetric matrix!!")))
次のように実行した。
(define A '((3 0.01 0.1) (0.01 2 0.01) (0.1 0.01 1))) ;;固有値が並んだ対角行列 (car (jacobi-method A 0.0001)) ;;固有ベクトルが並んだ行列 (cdr (jacobi-method A 0.0001))
