読者です 読者をやめる 読者になる 読者になる

数値計算,洲之内著,石渡改訂,サイエンス社 を読む (2)

9.1節「対称行列の固有値」を読み、ヤコビ法を実装した。

ヤコビ法
実対称行列を対角化するアルゴリズム
ギブンスの回転行列Pを利用して実対称行列Aの非対角成分を0に近づけていく。

A := P^T A P
この操作を繰り返し、行列Aを対角化する。(詳しい説明は省略)

ヤコビ法をschemeで実装し、以下の例題9.1でプログラムの動作をテストした。
固有値を近似計算できている。ちゃんと動作しているっぽい。
(対角化された行列が微妙に対称行列になってない。非対角成分のわずかな誤差。おそらく行列の積の計算時に発生した誤差が原因)

例題9.1

\displaystyle
A=
\begin{bmatrix}
 3 & 0.01 & 0.1 \\
 0.01 & 2 & 0.01 \\
 0.1 & 0.01 & 1
\end{bmatrix}
固有値および固有ベクトルをヤコビ法により求めよ.

解答:
先に結果を書き、後でソースコードを書く。
固有値は対角成分

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


ソースコード scheme

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