べき乗法による近似固有値の計算

今回の授業では、行列の固有値と固有ベクトルの計算方法を検討します。

以下の式を満たす$\lambda$と$x$を$A$の固有値と固有ベクトルという。

$$ Ax=\lambda x $$

$n \times n$行列 $A$ の固有値が以下のように並んでいる場合、べき乗法によって $A$ の固有値を計算することができます。

$$ |\lambda_1| \le |\lambda_2| \le \cdots \le |\lambda_{n-1}| < |\lambda_{n}| $$

行列の固有値$\lambda_i$に対応する単位固有ベクトルを$\phi_i$とします。以下$\lambda_n$と$\phi_n$の反復計算法を説明します。

1. 反復法による近似ベクトルを求める

行列Aは固有値n個を持っているの場合、行列の固有ベクトくは$R^n$の基底をなす。よって、与えられる$R^n$のベクトル$x$について、以下のように分解することができます。

$$ x=c_1 \phi_1 + c_2 \phi_2 + \cdots + c_n \phi_n $$

この時、行列$A$をベクトル$x$にかけると、$Ax$は以下の形となります。

$$ Ax = c_1 \lambda_1 \phi_1 + c_2 \lambda_2 \phi_2 + \cdots + c_n \lambda_n \phi_n $$

$Ax$の各成分の変化を見ると、$\phi_i$ の成分が$\lambda_i$倍となっていることがわかります。$\lambda_n$は絶対値の一番大きい固有値なので、 $\phi_n$の成分の変化が大きいです。よって、$x\to Ax$の計算を繰り返すと、$\phi_n$の成分が他の成分と比べて大きくなり、主要項となっています。 固有ベクトルの長さを単位化して、即ち、 $x\to x/\|x\|$,以下の反復法が得られます。

$$ \tilde{x}=Ax^{\{k\}}, x^{\{k+1\}} = \frac{\tilde{x}}{\|\tilde{x}\|} $$

最初のベクトル $x$ の $\phi_{n}$ 成分の係数$c_n$が非ゼロである場合、反復回数$k$を大きくなると、以下の結果が得られる。

$$\lim_{ k\to \infty} x^{\{k\}} = \phi_n \text{ OR } \lim_{k=\infty} x^{\{k\}} = -\phi_n $$

べき乗法のアルゴリズム

近似固有ベクトルの計算

$$ x^{\{k+1\}}= \frac{Ax^{\{k\}}}{\| Ax^{\{k\}} \|_2} $$

実際の計算では、以下のように$x$の更新を繰り返しします。

$$ \frac{Ax}{\| Ax \|_2} \to x $$

2. 近似固有ベクトルによる近似固有値を求める

厳密な固有ベクトル $x$と固有値$\lambda$は$Ax-\lambda x=0$を満たします。 よって、近似固有ベクトル $\tilde{x}$ に対して、以下の最小化問題を解いて$\lambda$の近似値を算出することができます。

$$ \min_{\lambda} \| A\tilde{x} - \lambda \tilde{x} \| $$

上記の最小化問題の目標関数$\| A\tilde{x} - \lambda \tilde{x} \|$ の平方は$\lambda$の2次多項式となります。即ち、

$$ \| A\tilde{x} - \lambda \tilde{x} \|^2 = \left( A\tilde{x} - \lambda \tilde{x} \right)^T (A\tilde{x} - \lambda \tilde{x} ) = (A^TAx)^2 - \lambda (Ax + A^T x) + \lambda^2 \tilde{x}^T\tilde{x} $$

この最小化問題が容易に容易に解けます。最小化する$\lambda$は以下の式で得られます。

$$ \lambda= \frac{1}{2} \frac{x^T(A+A^T)x}{x^Tx} $$

特に、行列 $A$ が対称であれば、

$$ \lambda = \frac{x^TAx}{x^Tx} $$

演習1

次の行列 $A$ の最大固有値と固有ベクトルを計算しなさい。 初期ベクトル$x^{\{0\}}$は乱数で作成してください。

$$ A= \left( \begin{array}{ccccc} 2 & -1 & 0 & 0 & 0\\ -1 & 2 & -1 & 0 & 0\\ 0 & -1 & 2 & -1 & 0\\ 0 & 0 & -1 & 2 & -1\\ 0 & 0 & 0 & -1 & 2 \end{array} \right) $$

(1) 収束の確認

(2) 近似固有値の誤差 $Err_1(k)$ と近似固有ベクトルの誤差 $Err_2(k)$ について、以下の評価式があります。

$$ Err_1(k) \approx c_1\alpha^{k},\quad Err_2(k) \approx c_2\beta^{k} $$

数値計算の結果によって、上記の誤差評価式の$\alpha$と$\beta$を求めてください。

(3) 異なる初期ベクトルを使うとき、$\alpha$と$\beta$が変わるかどうかについて検討してください。

補足1

MATLABが提供する行列Aのすべての固有値と固有ベクトルの計算方法

 [eig_vec, eig_value] = eig(A)

ここで、eig_vecの列ベクトルは固有ベクトル、eig_valueの対角成分が固有値である。具体的に、eig_vecのk列目のベクトルとeig_valueの(k,k)対角成分に対して,$A*eig_vec(:,k)=eig_value(k,k)*eig_vec(:,k)$が成り立つします。

補足2

近似固有ベクトルと近似固有値の誤差を計算します。特に、真の固有ベクトル(精度の良い固有ベクトル)と近似固有ベクトルの差を計算するとき、ベクトルの向き方向を確認するのは必要です。

演習2

行列Aの固有値$\lambda=1$に対応する固有ベクトル$v_1$を求めてください。(eig関数または手計算で求めること) ベクトル$v_1$をべき乗法の初期ベクトルとして使う場合、べき乗法でのベクトルが最大固有値の固有ベクトルに収束するかどうかについて確認してください。収束の理由を検討してください。