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

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

以下の式を満たす$\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^{\{0\}}$:初期ベクトル　（$x^{\{0\}}$の$\phi_n$方向における成分が$0$ではないことを仮定します。）
  • 反復計算
$$ x^{\{k+1\}}= \frac{Ax^{\{k\}}}{\| Ax^{\{k\}} \|_2} $$

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

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

（２）近似固有値の計算

  • 近似固有ベクトル $\tilde{x}$ に対して、以下の計算式で $\|A\tilde{x} - \lambda \tilde{x}\|$ を最小化する $\lambda$ を算出できます。
$$ \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$が変わるかどうかについて検討してください。

In [6]:
A=diag(ones(5,1)*2) - diag(ones(4,1),1) - diag(ones(4,1),-1) 
x =rand(5,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 =

   0.23549
   0.46998
   0.17980
   0.30208
   0.78517

補足1

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

 [eig_vec, eig_value] = eig(A)

ここで、eig_vecの列ベクトルは固有ベクトル、eig_valueの対角成分が固有値である。

In [7]:
[V,D]=eig(A)

V =

   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

補足2

厳密な固有値と固有ベクトルは未知ですが、MATLABのeig関数またはべき乗法を使って、近似固有値と近似固有ベクトルの算出ができます。

以下は反復計算を１００回を行って、良い近似固有値と固有ベクトルを算出しています。

In [8]:
%初期値ベクトル
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

solution_x =

   0.28868
  -0.50000
   0.57735
  -0.50000
   0.28868

solution_lambda =  3.7321

補足3

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

  • ベクトル$v_1$とベクトル$v_2$の内積が負であれば、二つのベクトルが成す角は$\pi/2$以上なります。　（異なる方向と考えられる）
  • ベクトル$v_1$とベクトル$v_2$の内積が正であれば、二つのベクトルが成す角は$\pi/2$以下なります。　（同じ方向と考えられる）
In [10]:
%初期値ベクトル
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-')
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -40 -35 -30 -25 -20 -15 -10 -5 0 0 20 40 60 80 100 gnuplot_plot_1a gnuplot_plot_2a

演習２

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

In [ ]:
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)