べき乗法のアルゴリズムを実現するコードを書き、
$$ A= \begin{pmatrix} 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{pmatrix} $$の絶対値最大固有値および対応する固有ベクトルの近似値を求めてください。ただし、
$$ x^{(0)}= \begin{pmatrix} 1\\ 0\\ 0\\ 0\\ 0 \end{pmatrix} ,\quad EPS=10^{-2},\quad K=100 $$とします。
A = [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];
EPS = 1E-2;
K = 100;
x = [1;0;0;0;0];
r = x'*A*x;
k = 0;
while norm(A*x-r*x) >= EPS && k <= K
x = A*x/norm(A*x); %xの更新
r = x'*A*x; %rの更新
k = k+1;
end
if k > K
printf("%d回の反復では、収束しなかった。\n",K)
else
printf("%d回の反復で、絶対値最大固有値の近似値%.5f\n",k,r)
printf("および対応する固有ベクトルの近似値\n")
printf("%.5f\n",x)
printf("を求めることができた。")
end
23回の反復で、絶対値最大固有値の近似値3.73196 および対応する固有ベクトルの近似値 0.29437 -0.50568 0.57731 -0.49426 0.28295 を求めることができた。
演習1と同じ行列 $A$ の固有値 $1$ に対応する固有ベクトルをeig関数または手計算で求め、求めた固有ベクトルを初期値ベクトル $x^{(0)}$ としてべき乗法を実行してください（$EPS$ と $K$ は演習1と同じにすること）。その結果になった理由について考察してください。
A = [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];
[X,D] = eig(A) %Aの固有値からなる対角行列Dと、対応する固有ベクトル（大きさ1）を列に持つ行列Xを返す
X =
2.8868e-01 -5.0000e-01 5.7735e-01 5.0000e-01 -2.8868e-01
5.0000e-01 -5.0000e-01 1.1489e-17 -5.0000e-01 5.0000e-01
5.7735e-01 2.2665e-16 -5.7735e-01 -1.5971e-16 -5.7735e-01
5.0000e-01 5.0000e-01 1.7971e-16 5.0000e-01 5.0000e-01
2.8868e-01 5.0000e-01 5.7735e-01 -5.0000e-01 -2.8868e-01
D =
Diagonal Matrix
0.26795 0 0 0 0
0 1.00000 0 0 0
0 0 2.00000 0 0
0 0 0 3.00000 0
0 0 0 0 3.73205
EPS = 1E-2;
K = 100;
x = X(:,2); %上の出力より、Xの2列目が固有値1に対応する固有ベクトル（計算誤差を含む）
%プログラムによって、Xの列から固有値1に対応する固有ベクトルを見つけることも可能
r = x'*A*x;
k = 0;
while norm(A*x-r*x) >= EPS && k <= K
x = A*x/norm(A*x);
r = x'*A*x;
k = k+1;
end
if k > K
printf("%d回の反復では、収束しなかった。\n",K)
else
printf("%d回の反復で、絶対値最大固有値の近似値%.5f\n",k,r)
printf("および対応する固有ベクトルの近似値\n")
printf("%.5f\n",x)
printf("を求めることができた。")
end
0回の反復で、絶対値最大固有値の近似値1.00000 および対応する固有ベクトルの近似値 -0.50000 -0.50000 0.00000 0.50000 0.50000 を求めることができた。
固有値 $1$ に対応する正規化された固有ベクトルを初期値ベクトル $x^{(0)}$ にしたとき、べき乗法では反復が行われなかった。
このような結果になるのは、$x^{(0)}$ から $\lambda^{(0)}=1$ が直ちに計算され、固有値と固有ベクトルの定義より
$$ \left\| Ax^{(0)}-\lambda^{(0)}x^{(0)}\right\| = 0 < EPS, $$すなわち反復の終了条件が最初の時点で満たされるためである（実際には計算誤差が含まれるため、「$=$」は「$\fallingdotseq$」）。
仮にこの反復の終了条件を考えない場合にも、「実対称行列において異なる固有値に対応する固有ベクトルは直交する」ことから
$$ u_1{}^T x^{(0)}=0 $$（ただし、$u_1$ は絶対値最大固有値に対応する固有ベクトル）であるため、べき乗法の仮定が満たされずうまく収束しないと予想される。
演習2において、次のように固有値 $1$ に対応する正規化されていない固有ベクトルを初期値ベクトル $x^{(0)}$ にする場合には、べき乗法で反復が $1$ 回だけ行われる。
A = [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];
EPS = 1E-2;
K = 100;
x = [-1;-1;0;1;1]; %固有値1に対応する正規化されていない固有ベクトル
r = x'*A*x;
k = 0;
while norm(A*x-r*x) >= EPS && k <= K
x = A*x/norm(A*x);
r = x'*A*x;
k = k+1;
end
if k > K
printf("%d回の反復では、収束しなかった。\n",K)
else
printf("%d回の反復で、絶対値最大固有値の近似値%.5f\n",k,r)
printf("および対応する固有ベクトルの近似値\n")
printf("%.5f\n",x)
printf("を求めることができた。")
end
1回の反復で、絶対値最大固有値の近似値1.00000 および対応する固有ベクトルの近似値 -0.50000 -0.50000 0.00000 0.50000 0.50000 を求めることができた。
これは、
$$ \lambda^{(0)}=x^{(0)}{}^T A x^{(0)}=1\cdot\left\|x^{(0)}\right\|^2\neq 1 $$より最初の時点では反復の終了条件が満たされないが、
$$ x^{(1)}=\frac{Ax^{(0)}}{\left\|Ax^{(0)}\right\|} $$が正規化されていることから $\lambda^{(1)}=1$ となり、$1$ 回目の反復の後に終了条件が満たされるためである。
そもそも、固有値をレイリー商による
$$ \lambda^{(k)}=\frac{x^{(k)}{}^T A x^{(k)}}{x^{(k)}{}^T x^{(k)}} $$ではなく
$$ \lambda^{(k)}=x^{(k)}{}^T A x^{(k)} $$で計算できるのは $x^{(k)}$ が正規化されていることが前提となっているため、厳密に言えばxの初期値に対して正規化x = x/norm(x)を行う方が良い。
ただし、もともと影響があるのは $\lambda^{(0)}$ のみであるため、そうしなくても基本的にべき乗法の実行に問題は生じない。
演習2の考察で言及したように、反復の終了条件
$$ \left\| Ax^{(k)}-\lambda^{(k)}x^{(k)}\right\| < EPS $$を考えずに単純に $K$ 回反復した後の結果を表示してみると、次のようになる。
A = [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];
K = 100;
x = [-0.5;-0.5;0;0.5;0.5]; %固有値1に対応する固有ベクトルの厳密値
%[X,D] = eig(A);
%x = X(:,2); %eig関数によって求めた固有値1に対応する固有ベクトル（計算誤差を含む）
r = x'*A*x;
k = 0;
while k < K
x = A*x/norm(A*x);
r = x'*A*x;
k = k+1;
end
printf("%d回の反復で、絶対値最大固有値の近似値%.5f\n",k,r)
printf("および対応する固有ベクトルの近似値\n")
printf("%.5f\n",x)
printf("を求めることができた。")
100回の反復で、絶対値最大固有値の近似値1.00000 および対応する固有ベクトルの近似値 -0.50000 -0.50000 0.00000 0.50000 0.50000 を求めることができた。
初期値ベクトルとして固有値 $1$ に対応する固有ベクトルの厳密値を用いる場合には、$100$ 回反復しても結果の固有値および固有ベクトルの値に変化はない。
一方で、初期値ベクトルとしてeig関数によって求めた固有ベクトルを用いる場合には、計算誤差の影響で $28$ 回目付近から固有値・固有ベクトルが大きく変化し、$100$ 回反復した結果としては絶対値最大固有値および対応する固有ベクトルに収束する。
上で扱ったべき乗法のアルゴリズムは、行列の積を何度も計算するためあまり計算効率が良くありません。効率が良くなるようにアルゴリズムを改良して、そのコードを書いてください（考える問題は演習1と同じにすること）。反復の終了条件を変更してもよいです。
演習1のコードではA*xという計算を何度も行っていることが計算効率の観点で良くないため、A*xの値を保持する変数を用意してコードを改良する。
A = [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];
EPS = 1E-2;
K = 100;
x = [1;0;0;0;0];
y = A*x; %A*xの値を保持する変数yを用意し、以後の計算ではなるべくyを使用する
r = x'*y;
k = 0;
while norm(y-r*x) >= EPS && k <= K
x = y/norm(y); %xの更新
y = A*x; %xの更新直後にyの更新が必要
r = x'*y; %rの更新
k = k+1;
end
if k > K
printf("%d回の反復では、収束しなかった。\n",K)
else
printf("%d回の反復で、絶対値最大固有値の近似値%.5f\n",k,r)
printf("および対応する固有ベクトルの近似値\n")
printf("%.5f\n",x)
printf("を求めることができた。")
end
23回の反復で、絶対値最大固有値の近似値3.73196 および対応する固有ベクトルの近似値 0.29437 -0.50568 0.57731 -0.49426 0.28295 を求めることができた。
反復回数を $k$ とするとき、演習1のコードではA*xの計算を $(4k+1)$ 回行っていたが、改良したこのコードでは $(k+1)$ 回に抑えることができている。