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

中島 研吾

東京大学情報基盤センター

同 大学院情報理工学系研究科数理情報学専攻 数値解析 (科目番号 500080)

• べき乗法

行列の固有値問題

0

x

x

Ax





,



を満足する



と _{x を求める} – ： 固有値（_{eigenvalue）} – x ： 固有ベクトル（eigenvector）

一般固有値問題（_{General Eigenvalue Problem）}

0

x

Mx

Ax





,



ここでは標準固有値問題を扱う 固有値 • 固有振動数 • 行列の性質に影響：スペクトル半径，条件数

固有値問題の例（1/3）

2

x

1 x k m k _m k 運動方程式 2 2 1 2 1 2 1 1 ) ( ) ( x m kx x x k x m x x k kx           x x _         m k m k m k m k dt d / 2 / / / 2 2        2 1 x x x

固有値問題の例（2/3）

x x _          m k m k m k m k dt d / 2 / / / 2 2 t j

e

m

k

m

k

m

k

m

k



a

x

A





















/

2

/

/

/

2

振動的な解を仮定

a

Aa





2 ω（固有円振動数）

固有値問題の例（3/3）

固有振動数(_{Natural Frequency)} （構造物などの）力学システムには，固有振動数が存在する. 固有振動数あるいは，それに近い周波数で 力学システムを加振すると，システムは共振を起す. 共振したシステムは，非常に大きな変位，ひずみ，応力を 生じて，システムが崩壊，破損する！ 共振を避けたり，抑制したりする設計が必要 （耐震設計・免振設計など）

固有値問題の計算（1/3）



















2

1

1

1

A

の固有値・固有ベクトルを求めよ.

0

x

x

Ax





,





A





I





x



0

,

x



0





0

det







A



I

特性方程式 0 1 3 2 1 1 1 ) det(  2              I A 特性方程式=0 2 5 3    2 5 3 , 2 5 3 2 1      

固有値問題の計算（2/3）

x

Ax





より 2 2 1 1 2 1

2

x

x

x

x

x

x















この連立方程式は、必ず不定 したがって，_x1,x2のどちらか一方を定数をおく. たとえば x₁=c₁とおけば x₂=（1-λ）c₁ 固有ベクトル：                                      1 5 1 1 7 . 2 1 2 5 1 1 1 1 2 1 1 1 c c c c x x    

固有値問題の計算例（3/3）

一般のn元の正方行列Aの固有値，固有ベクトルは，前述した ような方法で求めることができる

0

)

det(

A

 I





特性方程式は固有値λについてのn次の代数方程式（非線形） 大規模な次元（_>106_{)を有する行列の固有値問題も扱える} 方法が開発されている：実に様々な解法がある 実用上重要なのは（絶対値）最大・最小固有値 重根があると特別な扱い必要 - 本講義では基本的に重根は無しとする

• べき乗法

べき乗法（Power Method）

絶対値最大の実固有値と それに対応する固有ベクトルを求める方法 適当な初期ベクトル _x(0) _{から始めて} ) ( ) 1 ( ) 1 ( ) 2 ( ) 0 ( ) 1 ( k k _Ax x Ax x Ax x      Aをどんどん乗じていく 但し，単に乗じていくだけでは、 発散したり，原点に収束したり してしまうので，常に _x(k)_{の大きさを} 一定（例えば=1）に保つ必要がある. x(k)_{は絶対値最大の固有値に対応する固有ベクトルに収束していく}

べき乗法のアルゴリズム

• Step 0: ||x(0)_|| 2=1 である初期ベクトル x(0) を選び，k=0 とする • Step 1: 以下のように x(k+1) _{を更新する：} • Step 2: k=k+1としてStep 1を繰り返す





2 ) ( ) ( ) 1 ( ) ( ) ( ) ( ) (

_,

_,

k k k k k k k

y

y

x

y

x

Ax

y















：_{A の絶対値最大の実固有値に収束} x(k) _{：A の絶対値最大の実固有値に対応する} 固有ベクトルに収束

13

べき乗法が最大固有値に収束する理由（1/3）

n n c c c c x x x x y(0)  _{1 1}  _{2 2}  _{3 3} _ n n nc c c c x x x x Ay x(1)  (0)  _{1 1 1}  _{2 2 2}  _{3 3 3} _  n n k n k k k k k k _{Ay A y c x c x c x c x} x( )  ( 1)  (0)  ₁ '_{1 1}  ₂ '_{2 2}  ₃ '_{3 3} _  ' n    ₁  ₂  ₃ _  n x x x x₁, ₂, ₃,_, 固有値（絶対値の大きさ順） それに対応する固有ベクトル （一次独立と仮定）

べき乗法が最大固有値に収束する理由（2/3）

0 '₁  c if                                 n _n k n k k k k c c c c c c c x x x x y 1 1 3 1 3 1 3 2 1 2 1 2 1 1 1 ) ( ' ' ' ' ' ' '        _



2,3, ,



lim 0 1 1 1             k i k i _{i n}      1 1 1 ) ( _' : y c x k if   k   k  べき乗法によって求められるベクトル_{の「方向」が最大固有値 } x(k) に対応する 固有ベクトル _x1 のそれに収束していく

15

べき乗法が最大固有値に収束する理由（3/3）

2 ) 1 ( 1 1 1 1 2 ) 1 ( ) 1 ( ) ( 1 1 1 1 ) 1 ( _{' '} 1       _{    } k k k k k k k _{c c} y x y y x x y   2 ) 1 ( 1 1 1 ) ( ) ( 1 '  _    k k _k k c y x x Α y     





   





        1 2 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 ' , 1 ' 1 ' , 1 ' , , _                                                                  k k k k k k k k k k k k c c c c y x y x y x y x x x y x







 

,





,



1 , , , _{( ) ( )} 1 ) ( ) ( ) ( ) ( ) ( ) (    k k k k k k k k x x y x x x y x  

べき乗法の収束



2,3, ,



lim 0 1 1 1             k i k i _{i n}      |



_i/



₁|が1より充分小さいことが収束に影響，特に 以下の成立が高速な収束に必要 1 1 2   

17

べき乗法の例（1/3）



















2

1

1

1

A

の絶対値最大の固有値およびその固有ベクトルを べき乗法により求めよ.  0

_

 

₁

_,

₀

_,

_y

 0

_

_Ax

 0

_

 

₁

_,

_

₁

x

   

 

₂

 

1

,

1

1

1

,

1

1

2 0 1

_

_

_

_

y

x

1回目    



x

0

,

y

0







   

1

0

,

1

,



1





1

   

 

 

₁₃

 

2

,

3

1

3

,

2

2

1

13

2

3

,

2

2

1

1

2 1 2

_

_

_

_

_

_

y

x

2回目    









 

2

.

500

2

5

3

,

2

2

1

,

1

1

2

1

,

1 1

_

_

_











_

_



y

x

 

 

   

 

₂

_,

₃

2

1

,

1

,

1

2

1

_{1 1} 1

_

_

_y

_

_Ax

_

_

x

19

べき乗法の例（3/3）

3回目    

 

 

₈₉

 

5

,

8

1

8

,

5

13

1

89

13

8

,

5

13

1

1

2 2 3

_

_

_

_

_

_

y

x

   









 

2

.

6153

13

34

8

,

5

13

1

,

3

2

13

1

,

2 2

_

_

_











_

_



y

x

 

 

   

 

₅

_,

₈

13

1

,

3

,

2

13

1

_{2 2} 2

_

_

_y

_

_Ax

_

_

x

前述した厳密解

2

.

618034

2

5

3

1









逆べき乗法

絶対値「最小」の実固有値と それに対応する固有ベクトルを求める方法 1 ' 1 ' 1 1

,

   



















A

A

x

x

A

x

Ax

として

A

'

x





'

x

にべき乗法を適用する

A

LU



として_{LU分解を求めておくと効率が良い}

21

べき乗法の加速手法：原点移動（Shift）

|



₂/



₁|の値を小さくすることにより収束を加速する









x Ix x Bx x x I B Ax I A B x Ax p p p p where p               , : constant 行列Bの固有値（：行列Aの固有値） 行列_{Bの固有ベクトル（Aの固有ベクトルに一致）}





: : x p   適当な定数_{pを選択することにより行列Bの絶対値最大／2番目に大き} な固有値の比を小さくできれば，行列_{Bにべき乗法を適用した方が良い} 1 2 1 2        p p 行列Bの固有値 行列A

原点移動の効果

 

1

,

0

,

0

.

40

,

2

1

1

1

₍₀₎























x

p

A

下記の条件において_{Aの絶対値最大の固有値およびその固} 有ベクトルをべき乗法，原点移動付きべき乗法により求めよ. 原点移動無し 原点移動有り 1 1.000000E+00 1.000000E+00 2 2.500000E+00 2.617647E+00 3 2.615385E+00 2.618034E+00 4 2.617978E+00 5 2.618033E+00

23

べき乗法・原点移動付きべき乗法の例

べき乗法

do iter= 1, 10

Y(1)= A(1,1)*X(1) + A(1,2)*X(2) Y(2)= A(2,1)*X(1) + A(2,2)*X(2)

EIGEN= X(1)*Y(1) + X(2)*Y(2)

DL= dsqrt(Y(1)**2+Y(2)**2) X(1)= Y(1)/DL X(2)= Y(2)/DL enddo 原点移動付きべき乗法 X(1)= 1.d0; X(2)= 0.d0 A(1,1)= A(1,1) - SHIFT A(2,2)= A(2,2) - SHIFT

do iter= 1, 10

Y(1)= A(1,1)*X(1) + A(1,2)*X(2) Y(2)= A(2,1)*X(1) + A(2,2)*X(2)

EIGEN= X(1)*Y(1) + X(2)*Y(2) + SHIFT

DL= dsqrt(Y(1)**2+Y(2)**2) X(1)= Y(1)/DL

X(2)= Y(2)/DL enddo

• べき乗法

対称行列の固有値計算法

• 実対称行列の固有値⇒実数 • 弾性振動問題などで工学的に重要な実対称行列の固有 値計算法として代表的な手法について紹介する： – ハウスホルダ変換（Householder）による三重対角化（ tridiagonalization） – 二分法（Bi-Section）による固有値計算 – 逆反復法による固有ベクトル計算

• N×Nの正方行列A, Bに対して以下を満たすような正則 行列Pが存在するとする： B= P-1_{A P} • このときAとBは相似（similar）であると呼び，BはAを相 似変換した行列であると言う。 • AとBが相似であればそれらの固有値は一致する • 任意の固有値に対するBの固有ベクトルを x とすると，A の固有ベクトルは _{Px となる}

Householder変換：三重対角化（1/6）

N次のベクトルx,yに対して以下の行列Qを定義するとき，行列Q による相似変換をハウスホルダー変換（_{Householder）と呼ぶ：} 1 , 2        T _uuT y x y x u uu I Q 変換行列_{Qは対称かつ直交：}





 







 

uu u I u uu uu I uu I uu I QQ Q Q Q uu I uu I uu I Q                 T T T T T T T T T T T T T T 4 2 2 2 2 2 2 2

以下に示す対称行列_{AをQによって三重対角化する：}                                               n n n n nn nk n n nk kk k k n k n k a a a a a a a a a a a a a a a a             1 1 1 3 3 2 2 2 1 1 1 2 1 2 1 2 2 22 21 1 1 21 11 0 0 0 0 0 0 0 0 0 0 0 0 ~                             A A

Householder変換：三重対角化（3/6）

N次のベクトル x,y,u を以下のように置く ：                                                                                            1 31 21 3 2 1 1 31 21 1 11 1 31 21 11 0 1 , 0 0 0 , n n n n a a s a u u u u a a s a s a a a a a      x y u y x y x





2 1 2 31 2 21 s a an a       y _ x

変換行列 _Q1 を以下のように置く ：                                 2 2 2 2 2 2 2 2 1 2 1 2 2 0 2 2 1 2 0 2 2 2 1 0 0 0 0 1 2 n k n n n k k k n k T u u u u u u u u u u u u u u u                 uu I Q                                                           0 0 2 2 1 2 0 2 2 2 1 0 0 0 0 1 1 11 31 21 11 2 2 2 2 2 2 1                 s a a a a a u u u u u u u u u u n k k k n k x Q y

Householder変換：三重対角化（5/6）

                          nn nk n nk kk k n k a a a a a a a a a s s a ' ' ' 0 ' ' ' 0 ' ' ' 0 0 2 2 2 2 22 1 1 11 1 1 1 1 1                 AQ Q AQ Q B

 

      _{    }  



  n i i n n a a a a a a s 2 2 1 2 1 2 1 , 1 2 21 2 21 21 1 sign  s₁は以下のようにとられる。桁落ちを防ぐため， _a21 と_s1の符号 は同じになるようにする：

                          nn nk n nk kk k n k a a a a a a a a a s s a ' ' ' 0 ' ' ' 0 ' ' ' 0 0 2 2 2 2 22 1 1 11 1 1 1 1 1                 AQ Q AQ Q B この操作を（_{n-2）回繰り返すことによって行列Aは三重} 対角行列 に変換可能される_A~ 新たな_Aとする

Householder変換：非対称行列の場合

三重対角行列ではなく，下記に示すような上ヘッセンベルク行 列（_{Hessenberg）となる}                      * * 0 0 0 * * 0 0 * * * * 0 * * * * * * * * * * ~             A

実区間_{[a,b]において，実係数を持つ多項式f(x)が与えられた場} 合，以下の4条件を満たす実係数多項式の列 f(x), f₁(x), f₂(x), f₃(x), ..., f_l(x) は実区間_{[a,b]においてスツルム列をなすという。但し f0(x)=f(x)} ① 実区間_{[a,b]内の全ての点xに対して，隣り合う2つの多項式} f_k(x), f_k+1(x)は同時に0とならない ② 実区間_{[a,b]内のある点x0}で _{fk(x0)=0 ならば, fk-1(x0) fk+1(x0)<0} ③ 列の最後の式_{fl(x)は実区間[a,b]において一定の符号を持つ} ④ ある点_x0で _{f(x0)=0 ならば f’(x0) f1(x0)>0である}

スツルムの定理（Sturm’s Theorem）

• 多項式の列 f(x), f₁(x), f₂(x), f₃(x), ..., f_l(x) が実区間[a,b]におい てスツルム列をなし，_{f(a) f(b)≠0 とする} • xを固定して関数列 f(x), f₁(x), f₂(x), f₃(x), ..., f_l(x) を左から右に 見ていったときの符号の変化の回数を _{N(x) とする} • 多項式 f(x) の実区間[a,b]に存在する零点（解）の個数 n₀ は 以下の式で与えられる（証明略）： n₀ = N(a) - N(b)

• 三重対角行列 に対して行列 を考え，その第_k主小行 列を _pk



と置く：

 

k k k k k p                                   1 1 1 3 3 2 2 2 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0             A~



I  A~ • これを最後の行に関して展開すると以下の漸化式を得る：

  

  

  

  2 ₂

 

 1 1       _{k k k k} k p p p • k=2について成立するように下記のように仮定しておく：

二分法（2/4）

• k=n のとき以下のn次多項式の根が の固有値⇒Aの固有値：A~

 

  I  A~ n p • 上記多項式の以下の列はスツルム列を構成する（証明付録）

 

 , p ₁

 

 , p ₂

 

 , , p₁

   

 , p₀  p_{n n n } • 対称行列の固有値は全て実数であり，以下を仮定すると： • 実区間[a,b]に存在する零点（固有値）の個数 n₀ は： n₀ = N(a) - N(b)， n₀ =1なら実区間[a,b]に固有値が1個存在 •



より大きい固有値の個数は_N(



₎ – 証明略，スツルムの定理より導かれる n     ₁  ₂  ₃  _   _ 

• 二分法では，スツルムの定理を用いて行列の特性方程式の 根の存在範囲を狭めて行くことで固有値の近似解を得る。 • ある適当な実定数[a,b]に関して，もし



_k （_{k番目に大きい固有} 値）が区間_{[a,b]の間に存在するのであれば，以下が成立：}

 

a k N

 

b N k  ,  • 区間[a,b]を半分に狭めるために2点の中点 を考える。 • もし



_k が区間_{[a,c]に存在するならば，下記が成立する：} そうでなければ，



_k は区間_{[c,b]に存在する。} •



_kの存在する区間を改めて_{[a,b]と設定し以上を繰り返す。} • 正の微少量



に対して _|a-b|<



ならば



_{= (a+b)/2として終了。} 2 b a c  

 

a k N

 

c N k  , 

二分法（4/4）

• [a,b]の初期値は前述のゲルシュゴリンの定理（次頁）を使用し て以下のように設定することができる：





r b r a r _{i i i n} n i         _   , 0 , max _{1 0} 1      • 予めbを固定して絞りこめば最大固有値を最初に求められる – （k+1）番目に大きい固有値は _k を上限値として繰り返し適用すること で計算できる • 逆にaを固定して絞りこめば最小固有値を最初に求めることが でき，_{k番目に小さい固有値を下限として（k+1）番目に小さい} 固有値を求められる

中心が_aii，半径



の円で囲まれた複素平面内の領域を_Si   j i ij i a r このとき，行列_A（aij）の全ての固有値



_k は和集合 の内部に 存在する。すなわち以下を満たす行番号_{iが存在：}



n i i S 1 



   j i ij k ii a a  （証明） xを Ax=



_kx を満たすAの固有ベクトルとする。 xの絶対値最大の 成分を_xiとするとき， _Ax=



_{kx の第i行を書き下すと以下を得る。}



         j i _i j ij k ii x x a a 

逆反復法による固有ベクトル計算

Inverse Iteration • 二分法によって求めた固有値を



_k とすると適当な初期ベクト ル_x(0)_{について以下の方程式を解いていくと：} • k⇒∞のとき， x(i)_は固有値



kの固有ベクトルに収束していくこと が期待される。





   





 , 2 , 1 , 0 1 _{ }  i i _i kI A x x 

                                              294 . 0416 . 0 0 0 0 0416 . 398 . 117 . 0 0 0 0 117 . 641 . 318 . 0 0 0 0 318 . 47 . 1 25 . 1 0 0 2 0 25 . 1 2 . 12 42 . 7 0 0 0 0 42 . 7 00 . 6 1 1 1 1 1 1 1 2 2 2 2 2 1 2 3 3 3 3 1 2 3 4 4 4 1 2 3 4 5 5 1 2 3 4 5 6 A

計算例（2/2）



₁= 1.721E+01

{5.507E-01 5.187E-01 4.565E-01 3.678E-01 2.578E-01 1.327E-01}



₂= 1.988E+00

{5.187E-01 2.578E-01 -1.327E-01 -4.565E-01 -5.507E-01 -3.678E-01}



₃= 7.747E-01

{4.565E-01 -1.327E-01 -5.507E-01 -2.578E-01 3.678E-01 5.187E-01}



₄= 4.462E-01

{3.678E-01 -4.565E-01 -2.578E-01 5.187E-01 1.327E-01 -5.507E-01}



₅= 3.189E-01

{2.578E-01 -5.507E-01 3.678E-01 1.327E-01 -5.187E-01 4.565E-01}



₆= 2.652E-01

• スーパーコンピューティングへの招待 • 連立一次方程式の解法（直接法，反復法） • 偏微分方程式の数値解法 • 固有値解法 • C言語によるプログラミング（入門編） – 基礎的な事項（様々な原理）の説明，証明 • 数学的な背景をしっかりと理解した上で自分でプログラムを 作って動かして見ることが重要 • 色々なことにチャレンジしてほしい – 計算機を使いこなせること（数学的背景を理解した上でプログラム

• もし p_k()=p_k-1()=0 が成立すると，下記漸化式より， p_k-2()=0となる： • 従って，全ての j についてp_j()=0となってしまうため，下記よりこの仮 定はあり得ない：

スツルム列を構成することの証明（1/3）

① 実区間内の全ての点



に対して，隣り合う_{2つの多項式 pk(}



_), p_k+1(



)は同時に0とならない

  

  

  

  2 ₂

 

 1 1       _{k k k k} k p p p

 

1 0   p ② 実区間内のある点



₀で _pk(



_{0)=0 ならば, pk-1(}



_{0) pk+1(}



_0)<0

  

  

 

 

 

 

   

0 0 if 0 1 0 1 0 1 2 0 1 0 1 2 1 1                             k k k k k k k k k k k p p p p p p p p

 

 , p ₁

 

 , p ₂

 

 , , p₁

   

 , p₀  p_{n n n }

③ 列の最後の式 _p0(



_{)は実区間において一定の符号を持つ}

  

  

 

 





 

  



 

 

   

 

 

 

 

 

 





 

                   ' 1 1 ' ' 2 2 1 ' 1 1 ' 2 2 1 1                      k k k k k k k k k k k k k k k k p p p p p p p p p p p p ④ ある点



₀で _{pn (}



_{0)=0 ならば pn’(}



_{0) pn-1(}



_{0)>0である}

 

1 0   p • これは下記より明らか：

 

 , p ₁

 

 , p ₂

 

 , , p₁

   

 , p₀  p_{n n n }

スツルム列を構成することの証明（3/3）

• ここで下記のように q_k を定義すると（*1）は（*2）のように表される： （*2）

   

 

 

 

p

 

q

 

k n q p p p p q k k k k k k k k k , , 3 , 2 , 1 2 1 2 1 ' 1 1 '                    • ところで，以下が成立する：

 

   

   

 

 

 

1 0 1 ' 1 1 ' 1 ' 0 1 0 ' 1 1                    p p p p p p p q k  • したがって（*2）より以下が成立する：

 

 

   

 

 

 

   

0

 

0 0 , , 3 , 2 , 0 0 0 1 0 ' 0 ' 1 1 '                      n n n n n n n n n k p p p q p p p p q n k q  

 

 , p ₁

 

 , p ₂

 

 , , p₁

   

 , p₀  p_{n n n }