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

べき乗法は、行列の、絶対値が最大となる固有値と、

それに対応する固有ベクトルを求める。

（べき乗法では、固有値を全て求めることはできない。

まずはプログラムを書いて、実際に動作を確認してみる。）

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

"洲之内著,石渡改訂:数値計算,サイエンス社"

例題9.2



の絶対値最大の固有値および、それに対応する固有ベクトルをべき乗法で求めよ。

(解答)

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

絶対値最大の固有値の近似値

" ) ( display max-eiglist ) ( display "



それに対応する固有ベクトル

" ) ( 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では、べき乗法が重要な役割を果たす。





















