今回の演習では、行列の固有値と固有ベクトルの計算方法を検討します。
以下の式を満たす$\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$とします。
（１）近似固有ベクトルの計算
$$ x^{\{k+1\}}= \frac{Ax^{\{k\}}}{\| Ax^{\{k\}} \|_2} $$
実際の計算では、以下のように$x$の更新を繰り返しします。
$$ \frac{Ax}{\| Ax \|_2} \to x $$
（２）近似固有値の計算
$$ \lambda= \frac{1}{2} \frac{x^T(A+A^T)x}{x^Tx} $$
特に、行列 $A$ が対称であれば、
$$ \lambda = \frac{x^TAx}{x^Tx} $$
次の行列 $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$が変わるかどうかについて検討してください。
A=diag(ones(5,1)*2) - diag(ones(4,1),1) - diag(ones(4,1),-1)
x =rand(5,1)
MATLABが提供する行列Aのすべての固有値と固有ベクトルの計算方法
[eig_vec, eig_value] = eig(A)ここで、eig_vecの列ベクトルは固有ベクトル、eig_valueの対角成分が固有値である。
[V,D]=eig(A)
厳密な固有値と固有ベクトルは未知ですが、MATLABのeig関数またはべき乗法を使って、近似固有値と近似固有ベクトルの算出ができます。
以下は反復計算を１００回を行って、良い近似固有値と固有ベクトルを算出しています。
%初期値ベクトル
x =rand(5,1);
%100回の反復計算
for k=1:100
newx=A*x;
x=/norm(newx);
end
lambda = (x'*A*x)/(x'*x);
solution_x=x
solution_lambda = lambda
近似固有ベクトルと近似固有値の誤差を計算します。特に、真の固有ベクトル（精度の良い固有ベクトル）と近似固有ベクトルの差を計算するとき、ベクトルの向き方向を確認するのは必要です。
%初期値ベクトル
x =rand(5,1);
x_error_list=[];
lambda_error_list=[];
%100回の反復計算
for k=1:100
newx=A*x;
x=/norm(newx);
%% -------- 固有ベクトルの誤差 －－－－－－－－
direction_confirm = x'*solution_x;%ベクトルの内積によって、方向を確認します。
if direction_confirm >0 %同じ方向
x_error = norm(x - solution_x);
else %逆方向の場合
x_error = norm(x + solution_x);
end
x_error_list(k) = x_error;
%% -------- 固有値の誤差 －－－－－－－－
lambda = (x'*A*x)/(x'*x);
lambda_error = abs(lambda - solution_lambda); %誤差の絶対値を取る
lambda_error_list(k) = lambda_error;
end
hold on
plot(1:100, log(x_error_list),'r-')
plot(1:100, log(lambda_error_list),'b-')
行列Aの固有値$\lambda=1$に対応する固有ベクトル$v1$を求めてください。（eig関数または手計算で求めること） ベクトル$v1$をべき乗法の初期ベクトルとして使う場合、べき乗法でのベクトルが最大固有値の固有ベクトルに収束するかどうかについて確認してください。その収束のオーダーと理由を検討してください。
A=diag(ones(5,1)*2) - diag(ones(4,1),1) - diag(ones(4,1),-1)
x =rand(5,1)
[V,D]=eig(A)
v1=V(:,2)