プログラミング演習（P3）：固有値問題、冪乗法

1. n x n の正方行列 A について、
   
   A x =　a x 
   
   となるような 1 x n のベクトル x と、値 a を求める問題を固有値問題という。ここでは、行列 Ａ　の i 行 j　列成分 A_ij について、
   
   A_ij = A_ji　, A_ij は実数
   
   となる実対称行列に対する冪乗法を学ぶ。
   
   
2. 行列 A の固有値に同じ値のものがないとして、それらの絶対値の大きい順に
   
   |a₁| > |a₂| > ... > |a_n|
   
   と並べて、対応する固有ベクトル x₁, x₂, ..., x_n が一次独立であるとする。あるベクトル y　を x₁, x₂, ..., x_n の線形結合で
   
   y = c₁x₁ + c₂x₂ + ... +  c_nx_n
   
   とすると、これに行列 A を１回かけたものは Ax_j=a_jx_j, (j=1,...,n)より
   
   y^₍₁₎ = A y = c₁a₁x₁ + c₂a₂x₂ +  ... + c_na_nx_n
   
   さらに Ｎ 回かけたものは
   

   y^(N) = A y^(N-1) = A^N y
          
           = c₁a₁^Nx₁ + c₂a₂^Nx₂ +  ... + c_na_n^Nx_n

           = a₁^N [c₁x₁ + c₂(a₂/a₁)^Nx₂ + ... + c_n(a_n/a₁)^Nx_n]
   
   となる。i≠1 の場合、 |a_i/a₁| < 1 なので、 y^(N) (N->∞) ->  a₁^N c₁x₁ に収束し、
   
   [y^{(N) T} y^(N-1)] / [y^{(N-1) T} y ^(N-1)] -> a₁ ... (*)
   
   の最大固有値に収束する。 ここで T は転置を表す。[y^{(N) T} y^(N-1)]  は ベクトル y^(N)  とベクトル y^(N-1) の内積である。
   
   
3. アルゴリズムは
   


1. ベクトル　y を入力し、初期値として、y⁽⁰⁾=y/sqrt[y^Ty] として  y^(0)T y⁽⁰⁾ =1 となるように規格化する。 k=0 とする。
   
2. z^(k) = A y^(k) をつくる。（行列ベクトル積）
3. b^(k) = y^(k)Tz^(k) をつくる。（内積 ... 上の式(*)の左辺に対応）
4. y^(k+1) = z^(k)/ sqrt[z^(k)T z^(k)] として y を更新する。（規格化）
   
5. 十分小さい e に対して 3 の内積について | (b^(k)-b^(k-1))/b^(k) |　< e となるまで繰り返す。（収束判定）
6. そのときの b^(k) と y^(k+1) が固有値と固有ベクトルとなる。
   

 （例題P３−１）以下のプログラムは 4 x 4 の行列を matrix.d から、4行のベクトルを vector.d から読み込んで、冪乗法のアルゴリズムに従って、最大固有値と固有ベクトルの収束性をみるプログラムである。作成して実行してみましょう。

        ( 5 2 3 4 )
        ( 2 5 4 3 )
A =   ( 3 4 5 2 )
        ( 4 3 2 5 )

 (実行結果)
gcc pex3-1.c -o pex3-1 -lm
./pex3-1

b = 11.733333,   x =  0.485516 0.514076 0.485516 0.514076
b = 13.990212,   x =  0.497955 0.502037 0.497955 0.502037
b = 13.999800,   x =  0.499708 0.500291 0.499708 0.500291
b = 13.999996,   x =  0.499958 0.500042 0.499958 0.500042
b = 14.000000,   x =  0.499994 0.500006 0.499994 0.500006
b = 14.000000,   x =  0.499999 0.500001 0.499999 0.500001
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000
b = 14.000000,   x =  0.500000 0.500000 0.500000 0.500000

（課題P３−１）

上のプログラムを作成し、matrix.d と vector.d を読み込んで、実行しなさい。初期値を変えても固有値と固有ベクトルが求まるか確かめなさい

（課題P３−２）
得られた固有値と固有ベクトルが、もとの固有値問題を満足するか確かめるよう、プログラムを拡張しなさい。すなわち行列 A とベクトル x の積 A x と ベクトル x に固有値をかけたもの a x が一致するか確かめる。


（課題P３−３）
上の例題プログラムでは、行列ベクトルの積をkmax=20 回繰り返す。b の値を一回前とのものと比較して、十分違いが小さくなったら繰り返しをやめて（アルゴリズムの5）収束した固有値と固有ベクトルを出力するプログラムを作 成し、いくつかの行列と初期値について試してみなさい。収束の仕方を調べなさい。 (ヒント) プログラム中の /* 内積 */ の部分の b の値を新しい変数に代入しておく。そして、新しく計算された b の値と比較する。


（余裕がある人へ）
（課題P３−４）
n x n の実対称行列を入力すると、その絶対値の最大固有値と対応する固有ベクトルを出力する関数を作成せよ。

(さらに余裕がある人へ。成績評価にはあまり影響しない課題)

上の例題では、行列 A の最大固有値 a₁ と対応する固有ベクトル x₁ を求めた。ここでは、２番目に大きい固有値と対応する固有ベクトルを求める方法について考える。新しく行列

B = A - a₁ x₁ x₁^T

とすると、 B の固有値は、 0, a₂, a₃, ..., a_nとなり、対応する固有ベクトルは x₁, x₂, x₃, ... , x_n となる。これを用いて２番目に大きい固有値 a₂ と対応する固有ベクトル x₂ を冪乗法により求めることができる。

（課題P３−５）
例題プログラムに
B = A - a₁ x₁ x₁^T
を作って２番目に大きい固有値と固有ベクトルを求めるプログラムを作成し、実行しなさい。

（課題P３−６）
4 x 4 の行列の全ての固有値と固有ベクトルを求めるプログラムを作成し、 matrix.d について実行しなさい。