gijyutsu-keisan.com

3.ヤコビ法

ヤコビ法はn次実対称行列Aのすべての固有値と固有ベクトルを求める方法です。対称行列限定であることに注意が必要です。

3.1.ヤコビ法の概要

n次実対称行列Aに対して、非対角成分を0にする直交変換を繰り返し行うことで、Aを対角化させます (1章(1)の性質によります)。 このとき得られたAの対角成分がAの固有値となり、それぞれに対応する固有ベクトルは変換行列の積の各列ベクトルが対応します。

ヤコビ法で用いる直交変換Pは次の特徴を持つものを用います。
n=2を例にとり、A、Pを次のように設定します。
\[ A = \begin{pmatrix} a & b \\ b & c \end{pmatrix} , \ P = \begin{pmatrix} \cos{ \theta } & -\sin{ \theta } \ \\ \sin{ \theta } & \cos{ \theta } \end{pmatrix} \]
このとき\( P^t = P^{-1} \)(転置行列と逆行列は同じ)が成り立つので、Pは直交行列であることがわかります。 そこで、直交変換である\( P^t A P \)の非対角成分が0になるよう\( \theta \)を決定します。
\[ \begin{align} P^t A P & = \begin{pmatrix} \cos{ \theta } & \sin{ \theta } \\ -\sin{ \theta } & \cos{ \theta } \ \end{pmatrix} \begin{pmatrix} a & b \\ b & c \end{pmatrix} \begin{pmatrix} \cos{ \theta } & -\sin{ \theta } \ \\ \sin{ \theta } & \cos{ \theta } \end{pmatrix} \\ & = \begin{pmatrix} a' & b' \\ b' & c' \end{pmatrix} \end{align} \]
非対角成分b'を0にしたいので、
\[ \begin{align} b' & = ( c - a ) \sin{ \theta } \cos{ \theta } + b ( \cos^2{ \ \theta } - \sin^2{ \ \theta }) \\ & = \frac{ c - a }{ 2 } \ \sin{ 2 \theta } + b \cos{ 2 \theta } \\ & = 0 \end{align} \\ \]
\[ \begin{eqnarray} \therefore \begin{cases} & \tan{ 2 \theta } = \displaystyle \frac{ 2b }{ a - c } & ( a \neq c ) \\ & \theta = \displaystyle \frac{ \pi }{ 4 } & ( a = c ) \end{cases} \end{eqnarray} \]
この性質を利用して、n次に適用を拡大します。

準備として単位行列Eを列ベクトルで表します。
\[ E = ( \boldsymbol{ e_1 }, \ldots, \boldsymbol{ e_n } ) \]
直交行列Pを以下で定義します。

  1. \( P = E \)とおく
  2. \( P \)の\( s \)列目と\( t \)列目の\( s,t \)行成分を以下で置き換え\( ( s < t ) \)
    \[ p_{ss} = \cos{ \theta } , \ p_{st} = \sin{ \theta } , \ p_{is} = 0 ( i \neq s,t ) \\ p_{ts} = -\sin{ \theta } , \ p_{tt} = \cos{ \theta } , \ p_{it} = 0 ( i \neq s,t ) \]
    このとき、
    \[ \begin{eqnarray} \theta = \begin{cases} & \displaystyle \frac{ 1 }{ 2 } \tan^{-1}{ \frac{ 2 a_{st} }{ a_{ss} - a_{tt} } } & ( a \neq c ) \\ & \displaystyle \frac{ \pi }{ 4 } & ( a = c ) \end{cases} \end{eqnarray} \]

この定義によって\( P^t = P^{ -1 } \)を満たすことは、計算すれば確認できます。

次にsとtの決め方について見ていきます。
Aの上三角成分の中から絶対値最大となる成分\( |a_{st} | \)の行番号を\( s \)、列番号を\( t \)とします。 この\( s, t \)によって決まった直交行列\( P \)を\( P_0 \)とおいて、次の計算を行います。
\[ A^{(1)} = P_0^t A P_0 = ( a_{ij}^{(1)} \ ) \]
このとき、\( a_{st}^{(1)} = a_{ts}^{(1)} = 0 \)となることは成分計算することで確かめられます。

次に、\( A^{(1)} \)の上三角非対角成分の中から絶対値最大となる成分\( |a_{st} | \)の行番号を\( s \)、列番号を\( t \)とします。
このとき、\( \displaystyle max_{ ( s < t ) \ }⁡|a_{st} | < \epsilon \)(\( \epsilon \):収束判定値)であれば\( A^{(1)} \)は対角化されたとみなせ、\( A^{(1)} \)の対角要素\( a_{ii} \)が固有値\( \lambda_i, \ P_0 \)の列ベクトルが固有値に対応する固有ベクトルとなり、計算は完了です。\( \displaystyle max_{ ( s < t ) \ }⁡|a_{st} | \geq \epsilon \)なら次の計算を行います。
\[ A^{(2)} = P_1^t A^{(1)} P_1 = P_2^t P_1^t A P_1 P_2 \]
この操作を\( \displaystyle max_{ ( s < t ) \ }⁡|a_{st} | < \epsilon \)になるまで繰り返し計算を行います。
\[ A^{(k)} = P_{k-1}^t A^{(k-1)} P_{k-1} = P_{k-1}^t \ldots P_1^t A P_1 \ldots P_{k-1} \]
このとき\( A^{(k)} \)の対角成分が\( A \)の固有値\( a_{ii}^{(k)} , \ X = (P_0 \ldots P_{k-1} ) \)の各列ベクトルが\( A \)の固有値に対応する固有ベクトルになります。

3.2.ヤコビ法の計算フロー

ヤコビ法の計算フローは次の通りです。

ヤコビ法の計算フロー
図3.2-1 ヤコビ法の計算フロー

参考文献