Microsoft PowerPoint - Eigen.ppt [互換モード]

Transcription

1 固有値解析 中島研吾 東京大学情報基盤センター同大学院情報理工学系研究科数理情報学専攻数値解析 ( 科目番号 58)

2 行列の固有値問題 べき乗法 対称行列の固有値計算法 Eige

3 Eige A 行列の固有値問題 標準固有値問題 (Stdrd Eigevle Problem を満足する と を求める : 固有値 (eigevle) : 固有ベクトル (eigevetor) 一般固有値問題 (Geerl Eigevle Problem) A M ここでは標準固有値問題を扱う固有値 固有振動数 行列の性質に影響 : スペクトル半径 条件数

4 固有値問題の例 (/) 4 Eige m m 運動方程式 ) ( ) ( m m m m m m dt d / / / /

5 固有値問題の例 (/) 5 Eige m m m m dt d / / / / t j e m m m m A / / / / 振動的な解を仮定 A ω( 固有円振動数 )

6 Eige 6 固有値問題の例 (/) 固有振動数 (Ntrl Freqe) ( 構造物などの ) 力学システムには 固有振動数が存在する. 固有振動数あるいは それに近い周波数で力学システムを加振すると システムは共振を起す. 共振したシステムは 非常に大きな変位 ひずみ 応力を生じて システムが崩壊 破損する! 共振を避けたり 抑制したりする設計が必要 ( 耐震設計 免振設計など )

7 固有値問題の計算 (/) 7 Eige A の固有値 固有ベクトルを求めよ. A I A det I A 特性方程式 ) det( I A 特性方程式 = 5 5 5

8 固有値問題の計算 (/) 8 Eige A より この連立方程式は 必ず不定したがって のどちらか一方を定数をおく. たとえば = とおけば =(-λ) 固有ベクトル :

9 Eige 9 固有値問題の計算例 (/) 一般の 元の正方行列 A の固有値 固有ベクトルは 前述したような方法で求めることができる 特性方程式は固有値 λ についての 次の代数方程式 ( 非線形 ) det( A I) 大規模な次元 (> 6 ) を有する行列の固有値問題も扱える方法が開発されている : 実に様々な解法がある 実用上重要なのは ( 絶対値 ) 最大 最小固有値重根があると特別な扱い必要 - 本講義では基本的に重根は無しとする

10 行列の固有値問題 べき乗法 対称行列の固有値計算法 Eige

11 Eige べき乗法 (Power Method) 絶対値最大の実固有値とそれに対応する固有ベクトルを求める方法 適当な初期ベクトル () から始めて () () ( ) A A () () A ( ) A をどんどん乗じていく但し 単に乗じていくだけでは 発散したり 原点に収束したりしてしまうので 常に () の大きさを一定 ( 例えば =) に保つ必要がある. () は絶対値最大の固有値に対応する固有ベクトルに収束していく

12 Eige べき乗法のアルゴリズム Ste : () = である初期ベクトル () を選び = とする Ste : 以下のように (+) を更新する : ( ) A ( ) Ste : =+ として Ste を繰り返す ( ) ( ) ( ) ( ) ( ) () :A の絶対値最大の実固有値に収束 :A の絶対値最大の実固有値に対応する固有ベクトルに収束

13 Eige べき乗法が最大固有値に収束する理由 (/) () A () () A A () ) ( ) ( 固有値 ( 絶対値の大きさ順 ) それに対応する固有ベクトル ( 一次独立と仮定 )

14 Eige 4 べき乗法が最大固有値に収束する理由 (/) if ) ( lim i i i ) ( : if べき乗法によって求められるベクトル () の 方向 が最大固有値 に対応する固有ベクトル のそれに収束していく

15 Eige 5 べき乗法が最大固有値に収束する理由 (/) ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( Α ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (

16 Eige べき乗法の収束 i i i lim i / が より充分小さいことが収束に影響 特に以下の成立が高速な収束に必要 6

17 Eige A べき乗法の例 (/) の絶対値最大の固有値およびその固有ベクトルをべき乗法により求めよ. 回目 A 7

18 Eige べき乗法の例 (/) 回目 A

19 Eige べき乗法の例 (/) 回目 A 前述した厳密解

20 Eige 逆べき乗法 絶対値 最小 の実固有値とそれに対応する固有ベクトルを求める方法 A A A A として A にべき乗法を適用する LU A として LU 分解を求めておくと効率が良い

21 Eige べき乗法の加速手法 : 原点移動 (Shift) / の値を小さくすることにより収束を加速する A A B B : B I I A I where : ostt : 行列 Bの固有値 (: 行列 Aの固有値 ) 行列 Bの固有ベクトル (Aの固有ベクトルに一致) 適当な定数 を選択することにより行列 Bの絶対値最大 / 番目に大きな固有値の比を小さくできれば 行列 Bにべき乗法を適用した方が良い 行列 B の固有値行列 A

22 Eige 原点移動の効果 下記の条件において A の絶対値最大の固有値およびその固有ベクトルをべき乗法 原点移動付きべき乗法により求めよ. A (). 4 原点移動無し 原点移動有り.E+.E+.5E E+.6585E+.684E E+ 5.68E E+

23 Eige べき乗法 原点移動付きべき乗法の例 べき乗法 do iter= Y()= A()*X() + A()*X() Y()= A()*X() + A()*X() EIGEN= X()*Y() + X()*Y() DL= dsqrt(y()**+y()**) X()= Y()/DL X()= Y()/DL eddo 原点移動付きべき乗法 X()=.d; X()=.d A()= A() - SHIFT A()= A() - SHIFT do iter= Y()= A()*X() + A()*X() Y()= A()*X() + A()*X() EIGEN= X()*Y() + X()*Y() + SHIFT DL= dsqrt(y()**+y()**) X()= Y()/DL X()= Y()/DL eddo

24 行列の固有値問題 べき乗法 対称行列の固有値計算法 Eige 4

25 Eige 5 対称行列の固有値計算法 実対称行列の固有値 実数 弾性振動問題などで工学的に重要な実対称行列の固有値計算法として代表的な手法について紹介する : ハウスホルダ変換 (Hoseholder) による三重対角化 ( tridigoliztio) 二分法 (Bi-Setio) による固有値計算 逆反復法による固有ベクトル計算

26 Eige 6 相似変換 (Similr Trsformtio) N Nの正方行列 A Bに対して以下を満たすような正則行列 Pが存在するとする : B= P - A P このときAとBは相似 (similr) であると呼び BはAを相似変換した行列であると言う AとBが相似であればそれらの固有値は一致する 任意の固有値に対するBの固有ベクトルを とすると A の固有ベクトルは P となる

27 Eige 7 Hoseholder 変換 : 三重対角化 (/6) N 次のベクトル に対して以下の行列 Q を定義するとき 行列 Q による相似変換をハウスホルダー変換 (Hoseholder) と呼ぶ : Q I T T 変換行列 Q は対称かつ直交 : Q T T T T I I T T I T Q Q T Q QQ T T I I I T T T T I 4

28 Hoseholder 変換 : 三重対角化 (/6) 8 Eige 以下に示す対称行列 A を Q によって三重対角化する : ~ A A

29 Hoseholder 変換 : 三重対角化 (/6) 9 Eige N 次のベクトル を以下のように置く : s s s s

30 Hoseholder 変換 : 三重対角化 (4/6) Eige 変換行列 Q を以下のように置く : T I Q s Q

31 Hoseholder 変換 : 三重対角化 (5/6) Eige s s Q AQ AQ Q B i i s sig s は以下のようにとられる 桁落ちを防ぐため と s の符号は同じになるようにする :

32 Hoseholder 変換 : 三重対角化 (6/6) Eige s s Q AQ AQ Q B この操作を (-) 回繰り返すことによって行列 A は三重対角行列に変換可能される A ~ 新たな A とする

33 Hoseholder 変換 : 非対称行列の場合 Eige 三重対角行列ではなく 下記に示すような上ヘッセンベルク行列 (Hesseberg) となる * * * * * * * * * * * * * * * * * * ~ A

34 Eige 4 スツルム列 (Strm Chi/Seqee) 実区間 [b] において 実係数を持つ多項式 f() が与えられた場合 以下の4 条件を満たす実係数多項式の列 f() f () f () f ()... f l () は実区間 [b] においてスツルム列をなすという 但し f ()=f() 実区間 [b] 内の全ての点 に対して 隣り合うつの多項式 f () f + () は同時にとならない 実区間 [b] 内のある点 で f ( )= ならば f - ( ) f + ( )< 列の最後の式 f l () は実区間 [b] において一定の符号を持つ 4ある点 で f( )= ならば f ( ) f ( )>である

35 Eige 5 スツルムの定理 (Strm s Theorem) 多項式の列 f() f () f () f ()... f l () が実区間 [b] においてスツルム列をなし f() f(b) とする を固定して関数列 f() f () f () f ()... f l () を左から右に見ていったときの符号の変化の回数を N() とする 多項式 f() の実区間 [b] に存在する零点 ( 解 ) の個数 は以下の式で与えられる ( 証明略 ): = N() - N(b)

36 二分法 (/4) 6 Eige 三重対角行列に対して行列を考え その第 主小行列を と置く : A ~ A I ~ これを最後の行に関して展開すると以下の漸化式を得る : = について成立するように下記のように仮定しておく :

37 Eige 7 二分法 (/4) = のとき以下の 次多項式の根が A ~ の固有値 Aの固有値 : I ~ A 上記多項式の以下の列はスツルム列を構成する ( 証明付録 ) 対称行列の固有値は全て実数であり 以下を仮定すると : 実区間 [b] に存在する零点 ( 固有値 ) の個数 は : = N() - N(b) = なら実区間 [b] に固有値が 個存在 より大きい固有値の個数は N() 証明略 スツルムの定理より導かれる

38 Eige 8 二分法 (/4) 二分法では スツルムの定理を用いて行列の特性方程式の根の存在範囲を狭めて行くことで固有値の近似解を得る ある適当な実定数 [b] に関して もし ( 番目に大きい固有値 ) が区間 [b] の間に存在するのであれば 以下が成立 : N N b b 区間 [b] を半分に狭めるために 点の中点 を考える もし が区間 [] に存在するならば 下記が成立する : N N そうでなければ は区間 [b] に存在する の存在する区間を改めて [b] と設定し以上を繰り返す 正の微少量 に対して -b <ならば = (+b)/として終了

39 Eige 9 二分法 (4/4) [b] の初期値は前述のゲルシュゴリンの定理 ( 次頁 ) を使用して以下のように設定することができる : r r m i i i i b r 予め b を固定して絞りこめば最大固有値を最初に求められる (+) 番目に大きい固有値は を上限値として繰り返し適用することで計算できる 逆に を固定して絞りこめば最小固有値を最初に求めることができ 番目に小さい固有値を下限として (+) 番目に小さい固有値を求められる

40 Eige 4 ゲルシュゴリンの定理 (Gershgori) 中心が ii 半径 r i ij i j の円で囲まれた複素平面内の領域を S i このとき 行列 A( ij ) の全ての固有値 は和集合存在する すなわち以下を満たす行番号 iが存在 : ii ( 証明 ) i j ij S i i の内部に を A= を満たすAの固有ベクトルとする の絶対値最大の成分を i とするとき A= の第 i 行を書き下すと以下を得る ii i j これから直ちに結論を得る ij j i

41 Eige 4 逆反復法による固有ベクトル計算 Iverse Itertio 二分法によって求めた固有値を とすると適当な初期ベクトル () について以下の方程式を解いていくと : i i I A i のとき (i) は固有値 の固有ベクトルに収束していくことが期待される

42 計算例 (/) 4 Eige A

43 Eige 4 計算例 (/) =.7E+ {5.57E- 5.87E E-.678E-.578E-.7E-} =.988E+ {5.87E-.578E- -.7E E E E-} = 7.747E- {4.565E- -.7E E E-.678E- 5.87E-} 4 = 4.46E- {.678E E E- 5.87E-.7E E-} 5 =.89E- {.578E E-.678E-.7E E E-} 6 =.65E- {.7E E- 5.87E E E E-}

44 44 本講義のまとめ スーパーコンピューティングへの招待 連立一次方程式の解法 ( 直接法 反復法 ) 偏微分方程式の数値解法 固有値解法 C 言語によるプログラミング ( 入門編 ) 基礎的な事項 ( 様々な原理 ) の説明 証明 数学的な背景をしっかりと理解した上で自分でプログラムを作って動かして見ることが重要 色々なことにチャレンジしてほしい 計算機を使いこなせること ( 数学的背景を理解した上でプログラムを作れること ) は チャレンジ可能性の幅を大きく広げることになる

45 もし ()= - ()= が成立すると 下記漸化式より - ()= となる : 従って 全ての j について j ()= となってしまうため 下記よりこの仮定はあり得ない : スツルム列を構成することの証明 (/) 45 Eige 実区間内の全ての点 に対して 隣り合う つの多項式 () + () は同時に とならない 実区間内のある点 で ( )= ならば - ( ) + ( )< if

46 スツルム列を構成することの証明 (/) 46 Eige 列の最後の式 () は実区間において一定の符号を持つ 4 ある点 で ( )= ならば ( ) - ( )> である これは下記より明らか : (*)

47 スツルム列を構成することの証明 (/) 47 Eige ここで下記のように q を定義すると (*) は (*) のように表される : (*) q q q ところで 以下が成立する : q したがって (*) より以下が成立する : q q q