固有値計算
2．べき乗法
べき乗法は、n次正方行列Aの絶対値最大となる固有値と、それに対応する固有ベクトルを求める方法です。 Aは対称でも非対称でも構いません。2．1．べき乗法の概要
べき乗法は、以下の前提条件をもとに、反復計算によって絶対値最大の固有値と対応する固有ベクトルを求めます （今回の場合、算出する固有値は\( \lambda_1 \)になります）。＜前提条件＞
- Aの固有値\( \lambda_i (i=1,2,\ldots,n) \)に重複はない
- 固有値\( \lambda_i \)の添字\( i \)の順番は\( |\lambda_1 | > |\lambda_2 | \geq \ldots \geq |\lambda_{n-1} | \)
- 各固有値\( \lambda_i \)に対応する固有ベクトル\( \boldsymbol{x_i} \)は一次独立
\[
\boldsymbol{ u }
=
c_1 \boldsymbol{x_1} + c_2 \boldsymbol{x_2} + \ldots + c_{n} \boldsymbol{x_{n}}
\tag{2.1-1}
\]
このベクトル\( \boldsymbol{ u } \)に対して次の反復計算を行います（理由はあとでわかります）。ただし、この計算において\( \boldsymbol{ u^{(k)} } \ \)に大きさ1（\( |\boldsymbol{ u^{(k)} } | = 1 \)）の制約を与えます。
\[
\begin{align}
& \boldsymbol{ u^{(1)} } = A \boldsymbol{ u }
\\
& \boldsymbol{ u^{(2)} } = A \boldsymbol{ u^{(1)} } = A^2 \boldsymbol{ u }
\\
& \ldots
\\
& \boldsymbol{ u^{(k)} } = \ldots = A^k \boldsymbol{ u }
\end{align}
\tag{2.1-2}
\]
ところで、固有値と固有ベクトルは\( A \boldsymbol{ x_i } = \lambda_i \boldsymbol{ x_i } \)を満たすので、\( \boldsymbol{ u^{(k)} } \ \)は次のようにも表現できます。
\[
\begin{align}
\boldsymbol{ u^{(k)} }
& =
c_1 \lambda_1^{(k)} \boldsymbol{x_1} + c_2 \lambda_2^{(k)} \boldsymbol{x_2}
+ \ldots + c_{n} \lambda_{n}^{(k)} \boldsymbol{x_{n}}
\\
& =
\lambda_1^{(k)} \
\left \{
c_1 \boldsymbol{x_1} + c_2 \left ( \frac{ \lambda_2 }{ \lambda_1 } \right )^k \boldsymbol{x_2} +
\ldots + c_{n} \left ( \frac{ \lambda_{n} }{ \lambda_1 } \right )^k \boldsymbol{x_{n}}
\right \}
\end{align}
\tag{2.1-3}
\]
ここで反復回数を\( \infty (k \rightarrow \infty) \)にすると、（2.1-3）式は前提条件（2）より\( ￥\lambda_i / \lambda_1 \rightarrow 0 \)となって
\[
\boldsymbol{ u^{(k)} }
\rightarrow
c_1 \lambda_1^k
\boldsymbol{ x_1 }
\tag{2.1-4}
\]
に収束します。さらに\( \boldsymbol{ x_i }, \boldsymbol{ x_j } \)は線形独立、つまり\( \boldsymbol{ x_i } \cdot \boldsymbol{ x_j} = 0 ( i \neq j ) \)であることを考慮に入れると、次の2つの内積は
\[
\begin{align}
& \boldsymbol{ u^{(k+1)} } \cdot \boldsymbol{ u^{(k)} }
=
\lambda_1^{ 2k + 1 }
\left \{
c_1 | \boldsymbol{x_1} |^2 + c_2 \left ( \frac{ \lambda_2 }{ \lambda_1 } \right )^{ 2k+1 } | \boldsymbol{x_2} | ^2 +
\ldots +
c_{n} \left ( \frac{ \lambda_{n} }{ \lambda_1 } \right )^{ 2k+1 } | \boldsymbol{x_{n}} | ^2
\right \}
\\
& \boldsymbol{ u^{(k+1)} } \cdot \boldsymbol{ u^{(k)} }
=
\lambda_1^{ 2k }
\left \{
c_1 | \boldsymbol{x_1} |^2 + c_2 \left ( \frac{ \lambda_2 }{ \lambda_1 } \right )^{ 2k } | \boldsymbol{x_2} | ^2 +
\ldots +
c_{n} \left ( \frac{ \lambda_{n} }{ \lambda_1 } \right )^{ 2k } | \boldsymbol{x_{n}} | ^2
\right \}
\end{align}
\]
となります。これらの内積の比は、\( k \rightarrow \infty \)とすると\( \lambda_1 \)に収束することがわかります。
\[
\frac{ \boldsymbol{ u^{(k+1)} } \cdot \boldsymbol{ u^{(k)} } }
{ \boldsymbol{ u^{(k)} } \cdot \boldsymbol{ u^{(k)} } }
\ \rightarrow \
\lambda_1
\tag{2.1-5}
\]
\( |\boldsymbol{ u^{(k)} } | = 1 \)の制約条件から、\( \boldsymbol{ u^{(k+1)} } \cdot \boldsymbol{ u^{(k)} } \rightarrow \lambda_1\)でとなることがわかります。
さらに、\( \lambda_1 \)に対応する固有ベクトル\( \boldsymbol{ x_1 } \)は（2.1-4）式から求まります。
\[
\frac{ 1 }{ c_1 \lambda_1^k } \boldsymbol{ u^{(k)} }
\rightarrow
\boldsymbol{x_1}
\]
このように、反復計算を\( \infty (k \rightarrow \infty) \)繰り返すことで、Aの絶対値最大となる固有値と、それに対応する固有ベクトルの真値が得られることになります。
実際の計算では\( k \)を∞にすることは不可能なので、
\[
| \boldsymbol{ u^{(k+1)} }
\cdot
\boldsymbol{ u^{(k)} }
-
\boldsymbol{ u^{(k)} }
\cdot
\boldsymbol{ u^{(k-1)} } \ |
<
\epsilon
\]
（\( \epsilon \)は収束判定値）
となったところで、\( A \)の絶対値最大固有値（近似）：\( \lambda_1 = \boldsymbol{ u^{(k+1)} } \cdot \boldsymbol{ u^{(k)} } \ \)
\( \lambda_1 \)に対応した近似固有ベクトル（近似）：\( \boldsymbol{ x_1 } = \boldsymbol{ u^{(k+1)} } / | \boldsymbol{ u^{(k+1)} } \ | \ \)
とします。また、収束しない場合は∞ループに陥るため、最大反復回数を設定することで強制的にプログラムを終了させる必要があります （設定しなければフリーズしてしまいます）。
最後に、固有ベクトルは任意定数\( c_1 \)の不定性が残ります。 そのため、固有ベクトル\( \boldsymbol{ x_1 } \)は大きさ1のベクトルとして表現するのが一般的です （ちなみに\( 1 / \lambda_1^k \ \)もすべてのの成分に掛かっているので、\( \boldsymbol{ x_1 } \parallel \boldsymbol{ u^{(∞)} } \)となります）。
2．2．べき乗法の計算フロー
べき乗法の計算フローは次の通りです。解変数は固有値：\( \lambda_1 \)、固有ベクトル：\( \boldsymbol{ x_1 } \)になります。 また、\( \epsilon \)は収束判定値（例えば1E-6などの微小値）、\( N_{max} \)は最大反復回数です。
図2.2-1 べき乗法の計算フロー
