数値解析 第14回 (2) 累乗法
累乗法のアルゴリズム
上述の通り、手計算どおりの計算を実現するためにはたくさんの処理が必要です。
そこで、シンプルな反復法で固有値・固有ベクトルを計算する方法を紹介しましょう。
Th.5 対角化可能な $n$ 次行列 $A$ の固有値が条件
$\dps{|\,\lambda_1\,| \gt \max(|\,\lambda_2\,|,\cdots, |\,\lambda_n\,|)}$
$\cdots\cdots$ $(\ast)$
を満たすとする。( $\lambda_1$ を「最大固有値」と呼ぶ。)
このとき、次のアルゴリズムによって $(\lambda_1, \vvv_1)$ を求めることができる。
Algorithm 6
- $\vvv=$ ランダムな単位ベクトル
- $\www=$ ( $A\vvv$ 方向の単位ベクトル )
- if ( $\www$ ≒ $\vvv$ or $\www$ ≒ $-\vvv$ ) then 4° へ else (新$\vvv)=\www$ として 2° へ
- $\lambda=(\vvv,A\vvv)/(\vvv,\vvv)$ として $(\lambda,\vvv)$ を出力して終了
ただし $(\xxx,\yyy)$ は標準内積 $\dps{(\xxx,\yyy)=\sum_{i=1}^n\,x_i\,y_i}$ を表す。
実行例
Ex.7 5次行列
$\dps{
A=\mat{rrrrr}{\
-6 & -4 & -7 \ \\[-0.2em]
-4 & -3 & -3 \ \\[-0.2em]
-7 & -3 & -8 \ \\[-0.2em]
}
}$ の場合
実行例を折りたたむ
v = [ 0.4984323604, 0.80980835590, 0.3094763461]
v = [-0.6447177816, -0.41093312268, -0.6445719128]
v = [ 0.6310132733, 0.36167194503, 0.6863058014]
v = [-0.6301415521, -0.35584999826, -0.6901394084]
v = [ 0.6300168634, 0.35524585187, 0.6905643609]
v = [-0.6300047155, -0.35518066806, -0.6906089715]
v = [ 0.6300033728, 0.35517374127, 0.6906137588]
v = [-0.6300032313, -0.35517300111, -0.6906142685]
v = [ 0.6300032162, 0.35517292218, 0.6906143229]
v = [-0.6300032145, -0.35517291376, -0.6906143287]
v = [ 0.6300032144, 0.35517291286, 0.6906143293]
v = [-0.6300032144, -0.35517291276, -0.6906143294]
v = [ 0.6300032144, 0.35517291275, 0.6906143294]
λ1 = -15.928508004
v1 = [ 0.6300032144, 0.35517291275, 0.6906143294]
ただし縦ベクトルを横に表示しています。
Th.5 の証明
$A$ は対角化可能と仮定していますので、
固有ベクトルたち $\vvv_1$, $\cdots$, $\vvv_n$ は一次独立であり、
$\vvv$ の初期値 $\vvv_0$ は
$\dps{\vvv_0=\sum_{i=1}^n c_i\,\vvv_i=c_1\,\vvv_1+\cdots+c_n\,\vvv_n}$
と書くことができます。すると
$\dps{A\vvv_j=\lambda_j\,\vvv_j}$
ゆえ、$k=1,2,\cdots$ に対して
\begin{align}
A^k\vvv_0
&=c_1A^k\vvv_1+\cdots+c_nA^k\vvv_n \\
&=c_1\lambda_1^k\vvv_1+\cdots+c_n\lambda_n^k\vvv_n \\
\end{align}
が成り立ちます。仮定 $(\ast)$ より
$\dps{\lim_{k\rightarrow\infty}\left(\frac{\lambda_j}{\lambda_1}\right)^k=0}$
( $\forall\,j \gt 1$ )
が成り立ちますので、
$A^k\vvv_0$ の方向は $k\rightarrow\infty$ のとき $\vvv_1$ 方向へ収束してゆきます。
Alg.6 の $\vvv$ は $A^k\vvv_0$ 方向の単位ベクトルですので、
( $\pm 1$ 倍を除いて ) $\vvv_1$ 方向の単位ベクトルに収束します。
あとは
$A\vvv=\lambda\vvv$
から
$(\vvv,A\vvv)/(\vvv,\vvv)=(\vvv,\lambda\vvv)/(\vvv,\vvv)=\lambda(\vvv,\vvv)/(\vvv,\vvv)=\lambda$
が得られます。