べき乗法 - nibosiiwasi’s blog

今日は、べき乗法（power iteration）のプログラムをschemeで書く。

べき乗法は、行列の、絶対値が最大となる固有値と、
それに対応する固有ベクトルを求める。

（べき乗法では、固有値を全て求めることはできない。
まずはプログラムを書いて、実際に動作を確認してみる。）

いつものようにschemeで書いた。

$A= \begin{bmatrix} 1 & 2 \\ 1 & 1 \end{bmatrix}$
の絶対値最大の固有値および、それに対応する固有ベクトルをべき乗法で求めよ。

(解答)

(define (disp-mat a) (for-each (lambda (x) (display x) (newline)) a))

(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 (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 (power-iter A ini-v counter-limit)
  ;;ここから内部定義
  (define (power-iter-p A v1 v2 counter counter-limit)
    (if (>= counter counter-limit)
      (let ((l1 (car (transposed-mat v1)))
            (l2 (car (transposed-mat v2))))
        (let ((max-eiglist (map (lambda (x1 x2)
                                  (/ x2 x1))
                                l1
                                l2)))
          (begin
            (display "\n絶対値最大の固有値の近似値\n")
            (display max-eiglist)
            (display "\n\nそれに対応する固有ベクトル\n")
            (disp-mat (transposed-mat (list
                                       (map (lambda (x y)
                                              (/ x (expt y counter)))
                                            l2
                                            max-eiglist)))))))
      (power-iter-p A v2 (mul-mat A v2) (+ counter 1) counter-limit)))
  ;;ここまで内部定義
  (power-iter-p A '() ini-v 0 counter-limit))

本の回答と同じく、５回繰り返してみる。

(power-iter '((1 2) (1 1)) (transposed-mat '((1 1))) 5)

結果、一応、本の回答と同じ答えが出ることを確認した。

（「絶対値最大の固有値の近似値」をリストで表示する。
本来同じ値に収束すべき物であるが、並べて表示することにした。）

絶対値最大の固有値の近似値
(99/41 70/29)

それに対応する固有ベクトル
(115856201/96059601)
(20511149/24010000)

次は正確数でなく、不正確数（小数）で計算してみる。

(power-iter (scalar-mat 1.0 '((1 2) (1 1)))
            (transposed-mat '((1 1)))
            5)

結果

絶対値最大の固有値の近似値
(2.4146341463414633 2.413793103448276)

それに対応する固有ベクトル
(1.2060866357335795)
(0.8542752603082052)

さらに回数を増やして計算してみる。

(power-iter (scalar-mat 1.0 '((1 2) (1 1)))
            (transposed-mat '((1 1)))
            150)

結果。ほとんど収束した模様。

絶対値最大の固有値の近似値
(2.414213562373095 2.414213562373095)

それに対応する固有ベクトル
(1.2071067811865552)
(0.8535533905932793)

注１：初期値のベクトルは、非ゼロベクトルで、
絶対値最大の固有値に対応する固有ベクトルの成分を持っている必要がある。

注２：行列が、絶対値最大の固有値に、
絶対値が近い固有値を持つ場合、べき乗法の収束は非常に遅くなる。

べき乗法は、固有値を1個しか計算できないという意味で、
あまり使えないように見えるかもしれない。

しかし、べき乗法はかなり重要なアルゴリズムだと思う。

理由は２つある。
１つ目は、べき乗法のアルゴリズムは極めて単純でプログラミングが簡単。
２つ目は、絶対値最大の固有値、固有ベクトルが意味を持つ局面がある。
例えば、（今更ではあるが）GoogleのPageRankでは、べき乗法が重要な役割を果たす。

ここからは愚痴です

べき乗法について私が調べた中では、英語のWikipediaのページ
Power iteration - Wikipedia
が、かなり正確で詳しい。今日のブログを書く際、参考にさせてもらった。

現在の日本語のWikipediaでは、
べき乗法 - Wikipedia
不十分に感じる。理由は、Jordan標準形について一切触れていないから。

べき乗法については
"洲之内著,石渡改訂:数値計算,サイエンス社"
も同じ程度の説明しかない。
つまり行列が対角化可能であると（仮定？）して証明している。

（しかし線形代数の本に書かれていることだが、
一般的に、正方行列は対角化可能であるとは限らない。
正方行列はJordan標準形に相似変換が可能である。）

将来的には私も線形代数を完璧にマスターして、
日本語のWikipediaのページの修正に参加したい気分だ。
線形代数なら私でもたぶんなんとかなるでしょ。