[Pythonによる科学・技術計算]べき乗法による行列の固有値問題の数値解法，数値線形代数

#はじめに

理・工系では量子力学や振動子系の問題で固有値問題を解くことが多いです。

べき乗法は，標準固有値問題

$$A\mathbf{u} = \lambda \mathbf{u}$$
の絶対値最大の固有値(優越固有値)を求めるアルゴリズムの一つです。**
その辺の入門書を開くと最初に目に入ることが多いと思います。

固有値問題の初歩的な数値解法を学ぶうえでたいへん教育的な内容となっています。
本稿ではこの方法について考えていきます。

この方法はとても分かりやすいです。

計算手続き・ポイントは，

1. 適当な初期推定固有ベクトル$u_0$を設定する。
2. $u_0$ にひたすらAを掛ける ($AAAAAA...\mathbf{u_0} = A^k\mathbf{u_0} \ 　(k-> ∞)$。
3. すると，方向ベクトルが一定ベクトル$\mathbf{u'}$に収束する。
4. 多くの場合，その$\mathbf{u'}$が絶対値最大の固有値に対応する固有ベクトルになっている！
5. $u'$が固有ベクトルなら，それに対応する固有値$\lambda$は$\mathbf{u^T}A\mathbf{u}/(|u|^2)$であたえられる。

ということで，Aを掛けて方向ベクトルの収束性をチェックすることが主に行うことです。
それから，「べき乗法」の名前の由来は，上記2でわかるとおり，係数行列Aの累乗(べき)を評価する方法のためです。

＃なお，上記の4の数学的根拠が気になる方は本稿の補遺を参考にしてください。

また，べき乗法はAの絶対値最大の固有値を求める方法ですが，
**いくつか工夫することによって絶対値最小の固有値や，与えられた複素数$z$に最も近い固有値・2番目の固有値...なども求めることが可能です。**おいおい記事を投稿する予定です。

---

#内容

べき乗法を適用した簡単な固有値問題を解き，
ただしく絶対値最大の固有値とそれに対応する固有ベクトルが得られることを確認します。

$A\mathbf{u} = \lambda \mathbf{u}$の$A$として，

A=
\left[
\begin{matrix}
1 & 2 \\
3 & 4 
\end{matrix}
\right]

を考えます。

絶対値最大の固有値の厳密解は$\frac{5+\sqrt{33}}{2} = 5.372281323...$です。

---

#コード
収束判定条件として，
繰り返しステップ$k$と$k+1$における$\lambda$の変化の絶対値
$\frac{|\lambda_{k+1}-\lambda_k|}{| \lambda_k|} \lt \epsilon=0.0001$となるまで繰り返し計算を行います。

"""
行列の固有値問題: べき乗法:
"""

import numpy as np

A=np.array([[1,2],[3,4]])

x0 = np.array([1,0]); x1 = np.array([0,1])
u = 1.0*x0+2.0*x1 #初期固有ベクトルです。適当です。

rel_eps = 0.0001 # 固有値の収束条件
#クリロフ列生成

rel_delta_u=100.0
while rel_delta_u >= rel_eps :  # メインループ
    uu = u/np.linalg.norm(u) #正規化(ノルムを1にする)
    print("u=",uu)
    
    u = np.dot(A,uu.T)

    eigen_value=np.dot(uu,u)/(np.dot(uu,uu.T))

    print("eigen_value=",eigen_value)
    
    delta_u_vec = uu-u/np.linalg.norm(u)
    abs_delta_u_value= np.linalg.norm(delta_u_vec)
    rel_delta_u=abs_delta_u_value/np.abs(eigen_value) # 繰り返しステップに対する固有値の相対変化量
    print("rel_delta_u_vec = ",rel_delta_u)

#結果

u= [ 0.41612395  0.9093079 ]
eigen_value= 5.37244655582
rel_delta_u_vec =  3.29180183204e-05

厳密解 5.372281323...と非常に良く一致していることがわかります。

---

#参考文献
本稿を執筆するうえで参考になった書籍は以下の通りです。[1]は平易な記述で分かりやすいです。[2]はnumpyやscipyを利用した固有値問題の解法についてまとめたものです。

[1] ギルバート ストラング，『世界標準MIT教科書 ストラング:線形代数イントロダクション 』,近代科学社,2015.

[2] [[Pythonによる科学・技術計算] numpy・scipyを用いた(一般化)固有値問題の解法，ライブラリ利用]
(http://qiita.com/sci_Haru/items/034c6f74d415c1c10d0b)

---

#補遺
べき乗法は，標準固有値問題

$$A\mathbf{u} = \lambda \mathbf{u}$$

において，

初期ベクトル$u_0$から出発して

$$ u_i = Au_{i-1}, (i=1,2,...)$$
により絶対値が最大の固有値・固有ベクトルを求める方法です。
(この$u_0, Au_0, A^2u_0, A^3u_0$,...はクリロフ列と呼ばれています)

さて，固有値方程式
$$Au_i = \lambda_i u_i$$

において，固有値に縮退がない，つまり，
$$\lambda_i \neq \lambda_j \ (i \neq j)$$

と仮定します。
また，固有値の大きさの順番については

$$|\lambda_1| > |\lambda_2|>|\lambda_3|...$$

と仮定します。$|\lambda_1|$が絶対値最大の固有値で，べき乗法で求めたいものです。

さて，いま適当なベクトル$u_0$を考えます。$u_0$が以下のように展開できるように係数$c_i$が決まっているとします。

$$ u_0 = c_1 \mathbf{u_1} + c_2 \mathbf{u_2}+c3 \mathbf{u_3}+...$$

これに$A$のべき$A^k$を作用させると，

$$A^k u_0 = c_1 A^k \mathbf{u_1} + c_2 A^k \mathbf{u_2} + c_3 A^k \mathbf{u_3}+... $$
$$ = c_1 \lambda_1^k \mathbf{u_1} + c_2 \lambda_2^k \mathbf{u_2} + c_3 \lambda_3^k\mathbf{u_3} +...$$
$$ =\lambda_1^k (c_1 \mathbf{u_1} + c_2 \lambda_2^k/\lambda_1^k \mathbf{u_2} + c_3 \lambda_3^k/\lambda_1^k\mathbf{u_3}) +... $$

となります。

$|\lambda_1|$がもっとも大きい固有値なので，大きなべき数kに対して$\lambda_2^k/\lambda_1^k$ や$\lambda_3^k/\lambda_1^k$ ...はゼロに収束するはずです。つまり

$$ (k-> ∞ )A^k u_0 => \lambda_1^k (c_1 \mathbf{u_1} ) $$

となると期待できます。

こうして初期試行ベクトル$u_0$に繰り返し$A$を作用させることにより，絶対値最大の固有値に対応する固有ベクトルに平行なベクトルが得らます。

固有ベクトルは常に定数倍(任意複素数値倍)の不定性があります(固有方程式の両辺を複素数倍しても固有値は変わらないということです)ので，目的に応じたものが選ばれます。多くの場合，ノルムを1にするものが選ばれます。

絶対値最大の固有値$\lambda$は，十分大きな$k$に対して，固有値方程式から

$$\lambda = \frac{u_k^T A u_k}{u_k^T u_k} =\frac{u_k^T u_{k+1}}{u_k^T u_k} $$

として求めることができます。これはレイリー(Rayleigh)商と呼ばれるものです。

以上のようにして，固有値問題の絶対値最大の固有値およびそれに対応する固有ベクトルを計算すること可能です。

##注意

1. べき乗法を使う場合，$A$の固有値にダブりがないという条件が満たされている必要があります。
   固有値が縮退または近い値の場合，解が収束しないかまたは収束までに非常に多くの繰り返しステップが必要になります。
2. 実数行列が複素数固有値$z$を持つ場合，その共役複素数$z*$も固有値となります。そのときは$|z|=|z*|$となるため，固有値が同じになってしまいます。それらが絶対値最大の場合，上記1の状況になってしまうため，解が収束しません。