固有値解析 中島研吾 東京大学情報基盤センター同大学院情報理工学系研究科数理情報学専攻数値解析 ( 科目番号 -58)
行列の固有値問題 べき乗法 対称行列の固有値計算法 : ヤコビ法
A 行列の固有値問題 標準固有値問題 (Stndrd vlue Prolem を満足する と を求める : 固有値 (eigenvlue) : 固有ベクトル (eigenvector) 一般固有値問題 (Generl vlue Prolem) A M ここでは標準固有値問題を扱う固有値 固有振動数 行列の性質に影響 : スペクトル半径 条件数
固有値問題の計算 (/) 4 A の固有値 固有ベクトルを求めよ. A I A det I A 特性方程式 ) det( I A 特性方程式 = 5 5 5
固有値問題の計算 (/) 5 A より この連立方程式は 必ず不定したがって のどちらか一方を定数をおく. たとえば =c とおけば =(-λ)c 固有ベクトル :.6 5.7 5 c c c c
6 固有値問題の計算例 (/) 一般の n 元の正方行列 A の固有値 固有ベクトルは 前述したような方法で求めることができる 特性方程式は固有値 λ についての n 次の代数方程式 ( 非線形 ) det( A I) 大規模な次元 (> 6 ) を有する行列の固有値問題も扱える方法が開発されている : 実に様々な解法がある 実用上重要なのは ( 絶対値 ) 最大 最小固有値重根があると特別な扱い必要 - 本講義では基本的に重根は無しとする
7 べき乗法 逆べき乗法 標準固有値問題の解法 小規模問題 : ヤコビ法 中規模問題 : ハウスホルダー法 ( 対称 )QR 法 大規模問題 : 逆反復法 同時反復法 ( べき乗法の拡張 )
行列の固有値問題 べき乗法 対称行列の固有値計算法 : ヤコビ法 8
9 べき乗法 (Power Method) 絶対値最大の実固有値とそれに対応する固有ベクトルを求める方法 適当な初期ベクトル () から始めて () () ( ) A A () () A ( ) A をどんどん乗じていく但し 単に乗じていくだけでは 発散したり 原点に収束したりしてしまうので 常に () の大きさを一定 ( 例えば =) に保つ必要がある. () は絶対値最大の固有値に対応する固有ベクトルに収束していく
べき乗法のアルゴリズム Step : () = である初期ベクトル () を選び = とする Step : 以下のように (+) を更新する : y ( ) A ( ) Step : =+ として Step を繰り返す ( ) y ( ) ( ) y y ( ) ( ) () :A の絶対値最大の実固有値に収束 :A の絶対値最大の実固有値に対応する固有ベクトルに収束
べき乗法が最大固有値に収束する理由 (/) n c n c c c y () n n c n c c c Ay () () n n n c c c c y A Ay () ) ( ) ( n n 固有値 ( 絶対値の大きさ順 ) それに対応する固有ベクトル ( 一次独立と仮定 )
べき乗法が最大固有値に収束する理由 (/) c if n n n c c c c c c c y ) ( lim i i n i ) ( : y c if べき乗法によって求められるベクトル () の 方向 が最大固有値 に対応する固有ベクトル のそれに収束していく
べき乗法が最大固有値に収束する理由 (/) ) ( ) ( ) ( ) ( ) ( c c y y y y ) ( ) ( ) ( c y Α y c c c c y y y y y ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( y y
べき乗法の収束 i i i n lim i / が より充分小さいことが収束に影響 特に以下の成立が高速な収束に必要 4
A べき乗法の例 (/) の絶対値最大の固有値およびその固有ベクトルをべき乗法により求めよ. 回目 y A y y 5
べき乗法の例 (/) 回目 y A y y. 5 5 6
べき乗法の例 (/) 回目 y A 5 8 y 5 8 5 8 5 8 89 y 5 8. 65 89 4 前述した厳密解. 684 5 7
8 逆べき乗法 (Inverse Power Method) 絶対値 最小 の実固有値とそれに対応する固有ベクトルを求める方法 A A A ' A ' として ' A ' にべき乗法を適用する LU A として LU 分解を求めておくと効率が良い LU 分解による前進後退代入 :A - を乗じるのと同じこと
べき乗法の加速手法 : 原点移動 (Shift) / の値を小さくすることにより収束を加速する A A B B : p B pi pi A pi where p : constnt p : 行列 Bの固有値 (: 行列 Aの固有値 ) 行列 Bの固有ベクトル (Aの固有ベクトルに一致) 適当な定数 pを選択することにより行列 Bの絶対値最大 / 番目に大きな固有値の比を小さくできれば 行列 Bにべき乗法を適用した方が良い p p 行列 B の固有値行列 A 9
Solver-Direct Scil によるプログラム (/) http://nl.cc.u-toyo.c.jp/6n/ クリックして.zip をダウンロード
Solver-Direct Scil によるプログラム (/) http://nl.cc.u-toyo.c.jp/5n/.zip を解凍 直下に というディレクトリ Scilプログラムのソースファイル (*.sce*.sci) README.tt( 下記 : 必ず読むこと ) eigen-p.sce: Power Method with/without shift eigen-p.sce: Power Method 66 with/without shift eigen-p.sce: Inverse Power Method 66 with/without shift lu.sci : LU fctoriztion fs.sci: Forwrd/Bcwrd Sustitution LU 分解前進後退代入 プログラム名 内容
原点移動の効果 :eigen-p.sce 下記の条件において A の絶対値最大の固有値およびその固有ベクトルをべき乗法 原点移動付きべき乗法により求めよ. A () p. 4 原点移動無し 原点移動有り.E+.E+.5E+.67647E+.6585E+.684E+ 4.67978E+ 5.68E+ 6.684E+ SHIFT? -->.4
べき乗法の例 ( シフト無し ) eigen-p.sce n=; A=[-;-]; X=[;]; printf("shift?? n"); Shift= scnf("%f"); printf(" n"); printf("shift %e n" Shift); printf(" n"); printf("#originl n"); for iter= : Y()= A()*X() + A()*X(); Y()= A()*X() + A()*X(); = X()*Y() + X()*Y(); DL= sqrt(y()*y()+y()*y()); X()= Y()/DL; X()= Y()/DL; printf ("Iters %d %e %e %e n"iter X() X()); end y A ( ) ( )
べき乗法の例 ( シフト無し ) eigen-p.sce n=; A=[-;-]; X=[;]; printf("shift?? n"); Shift= scnf("%f"); printf(" n"); printf("shift %e n" Shift); printf(" n"); printf("#originl n"); for iter= : Y()= A()*X() + A()*X(); Y()= A()*X() + A()*X(); = X()*Y() + X()*Y(); DL= sqrt(y()*y()+y()*y()); X()= Y()/DL; X()= Y()/DL; printf ("Iters %d %e %e %e n"iter X() X()); end ( ) y ( ) 4
べき乗法の例 ( シフト無し ) eigen-p.sce n=; A=[-;-]; X=[;]; printf("shift?? n"); Shift= scnf("%f"); printf(" n"); printf("shift %e n" Shift); printf(" n"); printf("#originl n"); for iter= : Y()= A()*X() + A()*X(); Y()= A()*X() + A()*X(); = X()*Y() + X()*Y(); DL= sqrt(y()*y()+y()*y()); X()= Y()/DL; X()= Y()/DL; printf ("Iters %d %e %e %e n"iter X() X()); end ( ) y y ( ) ( ) 5
べき乗法の例 ( シフト有り ) eigen-p.sce X()=.; X()=.; A()= A() - Shift; A()= A() - Shift; for iter= : Y()= A()*X() + A()*X(); Y()= A()*X() + A()*X(); EIGEN= X()*Y() + X()*Y() + Shift; DL= sqrt(y()^+y()^); X()= Y()/DL; X()= Y()/DL; printf ("Iters %d %e %e %e %e n"iter EIGEN X() X() Shift); end A A B B : p B pi pi : A pi where p p p p : constnt 行列 Bの固有値 (: 行列 Aの固有値 ) 行列 Bの固有ベクトル (Aの固有ベクトルに一致) 6
計算例 6 5 4 A =.7E+ {5.57E- 5.87E- 4.565E-.678E-.578E-.7E-} =.988E+ {5.87E-.578E- -.7E- -4.565E- -5.57E- -.678E-} = 7.747E- {4.565E- -.7E- -5.57E- -.578E-.678E- 5.87E-} 4 = 4.46E- {.678E- -4.565E- -.578E- 5.87E-.7E- -5.57E-} 5 =.89E- {.578E- -5.57E-.678E-.7E- -5.87E- 4.565E-} 6 =.65E- {.7E- -.678E- 5.87E- -5.57E- 4.565E- -.578E-} 5 5 4 4 4 4 7
べき乗法の例 ( シフト無し ) eigen-p.sce n=6; X=[;;;;;]; Y=[;;;;;]; EIG=.; for iter= : for i=:n Y(i)=.; end for i=:n for j=:n Y(i)= Y(i) + A(ij)*X(j); end end RDL=./sqrt(DL); for i= :n X(i)= Y(i)*RDL; end RESID= sqrt((-eig)*.. (-EIG)/(EIG*EIG)); if RESID < EPS then re; end; EIG= ; end =.; DL =.; for i= :n = + X(i)*Y(i); DL = DL + Y(i)*Y(i); end RESID 8
原点移動付きべき乗法の例 eigen-p.sce n=6; X=[;;;;;]; Y=[;;;;;]; for i=:n A(ii)= A(ii) - Shift; end EIG=.; for iter= : for i=:n Y(i)=.; end for i=:n for j=:n Y(i)= Y(i) + A(ij)*X(j); end end = Shift; DL =.; for i= :n = + X(i)*Y(i); DL = DL + Y(i)*Y(i); end RDL=./sqrt(DL); for i= :n X(i)= Y(i)*RDL; end RESID= sqrt((-eig)*.. (-EIG)/(EIG*EIG)); if RESID < EPS then re; end; EIG= ; end 9
逆べき乗法の例 ( シフト無し ) eigen-p.sce eec( C: eigen lu.sci ); eec('c: eigen fs.sci'); n=6; X=[;;;;;]; Y=[;;;;;]; 外部関数 [A]=lu(An); LU 分解 ( 関数呼び出し ) EIG=.e-; for iter= : for i=:n Y(i)=X(i); end [Y]=fs (AYn); REIG=.; DL =.; for i= :n REIG= REIG + X(i)*Y(i); DL = DL + Y(i)*Y(i); end =./REIG; 前進後退代入 RDL=./sqrt(DL); for i=:n X(i)= Y(i)*RDL; end RESID= sqrt((-eig)*.. (-EIG)/(EIG*EIG)); if RESID < EPS then re; end; EIG= ; end A A ' ' A A ' ' A にべき乗法適用
逆べき乗法の例 ( シフト無し ) eigen-p.sce eec( C: eigen lu.sci ); eec('c: eigen fs.sci'); n=6; X=[;;;;;]; Y=[;;;;;]; 外部関数 [A]=lu(An); LU 分解 ( 関数呼び出し ) EIG=.e-; for iter= : for i=:n Y(i)=X(i); end [Y]=fs (AYn); REIG=.; DL =.; for i= :n REIG= REIG + X(i)*Y(i); DL = DL + Y(i)*Y(i); end =./REIG; 前進後退代入 RDL=./sqrt(DL); for i=:n X(i)= Y(i)*RDL; end RESID= sqrt((-eig)*.. (-EIG)/(EIG*EIG)); if RESID < EPS then re; end; EIG= ; end y A ' A 逆行列をかけるかわりに LU 分解の前進後退代入を実施する
逆べき乗法の例 ( シフト無し ) eigen-p.sce eec( C: eigen lu.sci ); eec('c: eigen fs.sci'); n=6; X=[;;;;;]; Y=[;;;;;]; 外部関数 [A]=lu(An); LU 分解 ( 関数呼び出し ) EIG=.e-; for iter= : for i=:n Y(i)=X(i); end [Y]=fs (AYn); REIG=.; DL =.; for i= :n REIG= REIG + X(i)*Y(i); DL = DL + Y(i)*Y(i); end =./REIG; 前進後退代入 RDL=./sqrt(DL); for i=:n X(i)= Y(i)*RDL; end RESID= sqrt((-eig)*.. (-EIG)/(EIG*EIG)); if RESID < EPS then re; end; EIG= ; end [A]= lu(an) 出力変数 関数入力変数 [y...yn]= foo(...m)
行列の固有値問題 べき乗法 対称行列の固有値計算法 : ヤコビ法
4 対称行列の固有値計算法 実対称行列の固有値 実数 弾性振動問題等 相似変換 N N の正方行列 A B に対して以下を満たすような正則行列 P が存在するとする : B= P - A P このときAとBは相似 (similr) であると呼び BはAを相似変換した行列であると言う AとBが相似であればそれらの固有値は一致する 任意の固有値に対するBの固有ベクトルを とすると Aの固有ベクトルは P となる ( 証明略 )
ヤコビ法 (Jcoi)(/5) 実対称行列 Aに対して二次元の回転に相当する相似変換を繰り返しながら対角行列へ近づけて行く方法 結果は合成されたn 次元の回転となる 繰り返しの段階で得られる行列において指定した行列成分がになるような相似変換 B= P - A P p 列目 q 列目 P =P = P = P =- それ以外は単位行列と同じ P - = P T 直交行列 B= P - A P 直交変換 B T = (P - A P) T = P T A (P - ) T = P - A P = B 対称 P 5 p 行目 q 行目
ヤコビ法 (Jcoi)(/5) 6 q p q q q p p p ij ij q p q p q p j i 4 tn tn if 繰り返すことによって行列全体の非対角成分が最終的に全て に収束することが期待される
7 q p B= P - A P A = P - A q p q p q q q p p p q p q q p p ' ' ' ' A :p 行 P: 列 ( = 他は ) A :q 行 P: 列 ( = 他は ) A :q A :p
8 q p B= P - A P A = P - A q p q p q q q p p p q p q q p p ' ' ' ' A :p 行 P: 列 ( = 他は ) A :q 行 P: 列 ( = 他は ) A :q A :p
' ' ' ' p q ' ' 9 B= P - A P A = P - A A :p A :q A :p 行 P: 列 ( = 他は ) A :q 行 P: 列 ( = 他は ) ' ' p q p q p q p p ' ' p q q q p q p q
4 B= P - A P A = P - A ' ' ' ' (A :p P:p 列 ) (A :p )
B= P - A P A = P - A ' ' ' ' ' ' (A :p ) 4 ' ' (A :p P:p 列 )
B= P - A P A = P - A ' ' ' ' (A :q P:q 列 ) (A :q ) 4
B= P - A P A = P - A ' ' ' ' ' ' ' ' (A :q P:q 列 ) (A :q ) 4
44 B= P - A P A = P - A ' ' ' ' (A :p P:q 列 ) (A :p )
B= P - A P A = P - A ' ' ' ' 45 ' ' (A :p ) ' ' (A :p P:q 列 )
46 B= P - A P A = P - A ' ' ' ' (A :q P:p 列 ) (A :q )
B= P - A P A = P - A ' ' ' ' 47 ' ' (A :q ) ' ' (A :q P:p 列 )
48 tr B tr n 実行列の要素の平方和 T T T T A A A A ij A ij A ji ij ji ij P AP i tr: trce 対角和 P AP n n i j n n i j T T T T T T T T B B tr P A P P AP tr P A AP tr P A AP tr T T A AP P tr tr T AB trba n n i j T T T T B B tr A AP P tr A A ij ij n n i j n n i j 直交変換によって行列成分の平方和 ( 二乗和 ) は一定に保たれる
ヤコビ法 (Jcoi)(/5) 49 i j p q ij ij ij ij p q p q p q 直交変換 ( 相似変換 ) によって対角成分の二乗和が増加 行列全体の成分の二乗和は直交変換 ( 相似変換 ) によって一定に保たれている 非対角成分の二乗和が減少している 繰り返すことによって に近づく
ヤコビ法 (Jcoi)(4/5) 5 以下のように相似変換 行列に番号をつける (A=A ) A T これまでの結果から : Λ m lim m T P m P P P AT A P m m P m m... ここでは下記の対角行列である 行列の固有値は相似変換により不変に保たれるため の対角成分が行列 Aの固有値に他ならない Λ n
ヤコビ法 (Jcoi)(5/5) 5 前ページの相似変換は下記のように表される : TΛ AT 行列 T の第 列からなる列ベクトルを t とすると 下記が成立 : At t 列ベクトル t は固有値 に対応する固有ベクトルとなる
A 回目の反復で になった () () 成分の絶対値が再び大きくなっている ### Originl Mtri.E+.E+.E+.E+.E+.E+.E+.E+.E+.E+.E+.E+.E+.E+.E+.E+ の絶対値最大成分 () -7.88E- 回転角度 ### ITERATIONS 5.4776E+ -.54E+ -4.44E-6.89E+ -5.97E- -.E-6 4.54E+.99E+.7E+.89E+.99E+.E+.E+ -5.97E-.7E+.E+.E+ の絶対値最大成分 () -.9474E- ### ITERATIONS 4.489E+ -.495E+ -7.457E-.E-6-8.75E- -7.457E- 4.54E+.79E+.7E+.E+.79E+.954E+ 7.94E- -8.75E-.7E+ 7.94E-.E+ 4 4.4464E- 5
A 回 になった成分の絶対値が再び大きくなっているが全体としては に近づいている ### ITERATIONS.7559E+ -.495E+ -.48E+.E-6-4.69E- -.48E+ 5.57E+.96E+ -.97E-6.E+.96E+.954E+ -.8E- -4.69E-.E+ -.8E- -.94E- 5.86E- ### ITERATIONS 4.684E+ -.495E+ -8.7E- 5.8E- -4.69E- -8.7E- 6.85E+ -.E-6-6.54E- 5.8E- -4.44E-6.675E+ -9.848E- -4.69E- -6.54E- -9.848E- -.94E- 9.585E- ### ITERATIONS 5.688E+ -.576E+.E+ 5.776E- -4.7E-.665E-6 6.9E+ -5.5E- -.88E- 5.776E- -5.5E-.675E+ -9.848E- -4.7E- -.88E- -9.848E- -.94E- 4-4.7655E-8 5
A ### ITERATIONS 7.94E- -.77E+ -.6E-.46E-6 -.77E- -.6E- 6.94E+.99E-6-4.9E-7 -.E- 4.45E-8.745E+ -7.456E- -.77E- -5.94E- -7.458E-.867E-.6797E- 54 ### ITERATIONS 8.4596E- -.77E+ -.46E-6.46E-6 -.77E- -6.46E-7 6.94E+.99E-6-4.9E-7 -.E- 4.45E-8.745E+ -7.456E- -.77E-.E- -7.458E-.867E- 4 6.9487E- ### ITERATIONS 9.54565E- -.77E+ -.46E-6.46E-6.E-7-6.767E-7 6.94E+.99E-6-4.9E-7 -.785E- 4.45E-8.745E+ -7.456E-.E+.E- -7.458E-.867E- 4 絶対値最大成分 <.e- のため終了
55 A ### ITERATIONS 9.54565E- -.77E+ -.46E-6.46E-6.E-7-6.767E-7 6.94E+.99E-6-4.9E-7 -.785E- 4.45E-8.745E+ -7.456E-.E+.E- -7.458E-.867E- 4 ### vlues/vectors -.77E+ 7.7E- -5.84E- -4.77E-.9E- 6.94E+ 5.8E- 4.984E- 5.45E-.56E-.745E+.994E- -6.E- 7.8E- -.766E- 4.867E- -.999E- -.7E-.E- 8.677E- 固有値 ~ に対応する固有ベクトル に対応する固有ベクトル に対応する固有ベクトル に対応する固有ベクトル
レポート課題 56 提出期限 : 月 日 ( 月 ):( レポートボックス設置 ) 下記の対称行列の固有値をヤコビ法によって求めよ 算出手順の概要 プログラムリスト ( プログラムを作成する場合 ) p.5-55 に相当する算出経過を添付すること 更に固有ベクトルを求めた場合は加点する 固有値のみ :7 点 ( 合格 ) + 固有ベクトル : 点 (T を保存しておく ) 4 4 4 4 5 5 4 5 6 A
プログラム例 Fortrn(/4) suroutine JACOBI (AEVNOUTMAX) implicit REAL*8 (A-HO-Z) dimension A(NN) E(N) V(NN) integer P Q NM= N do iter= MAX 57 P= Q= do i= NM do j= i+ N if (ds(a(ij)).gt.ds(a(pq))) then P= i Q= j endif enddo enddo 最大の A(PQ) を探索 (P Q)
do iter= MAX ( ) プログラム例 Fortrn(/4) if (ds(a(pq)).lt..d-) eit if (ds(a(pp)-a(qq)).lt..d-4) then T= dtn(.d) else R=.d*A(PQ)/(A(PP)-A(QQ)) T= dtn(r)*.5d endif T : if tn 4 tn. 58 最大の A(PQ) が - 未満だったら終了
do iter= MAX プログラム例 Fortrn(/4) 59 ( ) S= d(t) C= d(t) do = N Ap= A(P) Aq= A(Q) A(P)= Ap*C + Aq*S A(Q)= -Ap*S + Aq*C enddo do = N Ap= A(P) Ap= A(Q) A(P)= Ap*C + Aq*S A(Q)= -Ap*S + Aq*C enddo
do iter= MAX プログラム例 Fortrn(/4) 6 ( ) S= d(t) C= d(t) do = N Ap= A(P) Aq= A(Q) A(P)= Ap*C + Aq*S A(Q)= -Ap*S + Aq*C enddo B= P - A P A = P - A p q
do iter= MAX プログラム例 Fortrn(/4) 6 ( ) S= d(t) C= d(t) B= P - A P A = P - A B = A P do = N Ap= A(P) Ap= A(Q) A(P)= Ap*C + Aq*S A(Q)= -Ap*S + Aq*C enddo ' ' ' p ' ' ' q
do iter= MAX プログラム例 Fortrn(4/4) 6 ( ) VAL=.d do i= N do j= N if (i.ne.j) VAL= VAL + A(ij)** enddo enddo [A] の非対角成分の二乗和の平方根を計算 VAL= dsqrt(val) write (*'(i8 pe6.6)') iter VAL enddo return end