2022年度プログラミング演習A・B

第11回：行列の固有値の計算1

1. 行列の固有値と固有ベクトル

今回と次回の授業では、与えられた行列の固有値と固有ベクトルを近似的に計算する方法について扱います。

1.1. 数学的な定義

$n$ 次正方行列 $A$ に対して、

$$ Ax=\lambda x $$

を満たすスカラー $\lambda$ と零ベクトルでない $n$ 次元ベクトル $x$ を、それぞれ $A$ の固有値（eigenvalue）、$A$ の $\lambda$ に対応する固有ベクトル（eigenvector）と言います。

$\lambda$ が $A$ の固有値であることは

$$ \det \left(A-\lambda I_n\right)=0 $$

（ただし、$\det$ は行列式、$I_n$ は $n$ 次単位行列）と同値であり、これは $\lambda$ についての $n$ 次方程式であることから、$A$ は重複を許して $n$ 個の固有値を持つことが分かります。

また、$x$ が $A$ の $\lambda$ に対応する固有ベクトルであるならば、その $0$ でないスカラー倍もまた $A$ の $\lambda$ に対応する固有ベクトルとなります。よって、$1$ 個の固有値に対応する固有ベクトルは無数に存在しますが、重要なのはその方向です。この後では、特に大きさが $1$ であるような固有ベクトル（単位固有ベクトル）を考えることになります。

1.2. MATLABで用意された計算方法

MATLABでは、行列の固有値と固有ベクトルを計算するための関数eigが用意されています。以下がその使用例です。

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を返す

for k=1:size(A,1)
    A*X(:,k)-D(k,k)*X(:,k)  %Ax-λxの値を表示して、零ベクトルになるか確認する
end

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 =

   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

ans =

  -1.8041e-16
  -3.3307e-16
   9.7145e-16
  -1.6653e-16
  -3.7470e-16

ans =

   0.0000e+00
   2.2204e-16
   1.1562e-16
  -2.2204e-16
  -4.4409e-16

ans =

   0.0000e+00
   0.0000e+00
  -2.2204e-16
  -3.5942e-16
  -2.2204e-16

ans =

  -2.2204e-16
   6.6613e-16
   2.7073e-16
   2.2204e-16
   4.4409e-16

ans =

   0.0000e+00
  -2.2204e-16
   4.4409e-16
   0.0000e+00
   2.2204e-16

2. べき乗法

べき乗法（power iteration, power method）は、行列の絶対値最大固有値および対応する固有ベクトルを反復計算によって近似的に求めるアルゴリズムです。

2.1. べき乗法の考え方

$A$ を対角化可能な $n$ 次正方行列とし（特に、$A$ が実対称行列であれば対角化可能）、$A$ が重複を許して

$$ |\lambda_1|>|\lambda_2|\geq |\lambda_3|\geq \cdots \geq |\lambda_n| $$

を満たす $n$ 個の固有値 $\lambda_1,\ldots,\lambda_n$ を持つと仮定します。

さらに、$u_1$ を $A$ の $\lambda_1$ に対応する単位固有ベクトルとし、与えられた初期値ベクトル $x^{(0)}$ が

$$ u_1{}^T x^{(0)}\neq 0 $$

を満たすと仮定します。

このとき、

$$ \lambda^{(k)}=x^{(k)}{}^T A x^{(k)},\quad x^{(k+1)}=\frac{Ax^{(k)}}{\left\|Ax^{(k)}\right\|} $$

によって数列 $\left\{\lambda^{(k)}\right\}_{k=0}^{\infty}$ とベクトル列 $\left\{x^{(k)}\right\}_{k=0}^{\infty}$ を定めると、

$$ \lim_{k\to\infty}\lambda^{(k)}=\lambda_1 $$

が成り立ち、$x^{(k)}$ は $k\to\infty$ のとき $u_1$ と同じ方向に近付きます。

※細かく言うと、$u_1{}^T x^{(0)}>0$ かつ $\lambda_1>0$ の場合には $\lim_{k\to\infty}x^{(k)}=u_1$、$u_1{}^T x^{(0)}<0$ かつ $\lambda_1>0$ の場合には $\lim_{k\to\infty}x^{(k)}=-u_1$、$\lambda_1<0$ の場合には $x^{(k)}$ は $u_1$ と $-u_1$ に交互に近付きます。

（証明）

$A$ が対角化可能であることから、ベクトルの集合 $\{u_1,\ldots,u_n\}$ が一次独立、すなわち基底となるように、固有値 $\lambda_2,\ldots,\lambda_n$ にそれぞれ対応する固有ベクトル $u_2,\ldots,u_n$ を取ることができます。

このとき、初期値ベクトル $x^{(0)}$ を

$$ x^{(0)}=c_1u_1+c_2u_2+\cdots+c_nu_n $$

の形で表すことができ、仮定より $c_1=u_1{}^T x^{(0)}\neq 0$ です。

$x^{(k)}$ の定義より

$$ x^{(k)}=\frac{Ax^{(k-1)}}{\left\|Ax^{(k-1)}\right\|}=\frac{A\frac{Ax^{(k-2)}}{\left\|Ax^{(k-2)}\right\|}}{\left\|A\frac{Ax^{(k-2)}}{\left\|Ax^{(k-2)}\right\|}\right\|}=\frac{A^2x^{(k-2)}}{\left\|A^2x^{(k-2)}\right\|}=\cdots=\frac{A^kx^{(0)}}{\left\|A^kx^{(0)}\right\|} $$

であり、各 $i$ について $Au_i=\lambda_iu_i$ であるから

$$ \begin{aligned} A^kx^{(0)} & =c_1\lambda_1^ku_1+c_2\lambda_2^ku_2+\cdots+c_n\lambda_n^ku_n\\ & =c_1\lambda_1^k\left(u_1+\frac{c_2}{c_1}\left(\frac{\lambda_2}{\lambda_1}\right)^ku_2+\cdots+\frac{c_n}{c_1}\left(\frac{\lambda_n}{\lambda_1}\right)^ku_n\right) \end{aligned} $$

が成り立ちます。

仮定より

$$ \left|\frac{\lambda_i}{\lambda_1}\right|<1,\quad i=2,\ldots,n $$

を用いれば、$x^{(k)}$ が $k\to\infty$ のとき $u_1$ と同じ方向に近付くことが分かります。

さらに、簡単のために $\lim_{k\to\infty}x^{(k)}=u_1$ とすれば、

$$ \lim_{k\to\infty}\lambda^{(k)}=\lim_{k\to\infty}\left(x^{(k)}{}^T A x^{(k)}\right)=u_1{}^T A u_1=\lambda_1\|u_1\|^2=\lambda_1 $$

が得られます。

（解釈）

上ではベクトルの正規化（大きさを $1$ にする処理）のために式が分かりづらくなっていますが、べき乗法の本質としては、

• 初期値ベクトル $x^{(0)}$ に行列 $A$ を左から何度も掛けると、絶対値最大固有値に対応する固有ベクトルの方向に近付いていく

• レイリー商 $R_A(x)=\frac{x^TAx}{x^Tx}$ の「固有ベクトルに対して固有値を返す」という性質を利用して、固有ベクトルの近似値から固有値の近似値を計算できる

ということが重要です。

2.2. 反復の終了条件

べき乗法における反復計算の終了条件は様々に考えられますが、ここでは $Ax^{(k)}$ が $2$ ノルムの意味で $\lambda^{(k)}x^{(k)}$ とほぼ等しくなったとき、すなわち、微小正数 $EPS$ に対して

$$ \left\| Ax^{(k)}-\lambda^{(k)}x^{(k)}\right\| < EPS $$

を満足したときに、「収束した」とみなして反復を終了することにします。

さらに、数列 $\left\{\lambda^{(k)}\right\}_{k=0}^{\infty}$ およびベクトル列 $\left\{x^{(k)}\right\}_{k=0}^{\infty}$ が収束しない場合もあるため、反復回数が事前に決めた上限 $K$ を超えたときには「収束しなかった」とみなして反復を終了することにします。

2.3. プログラムとしての実装

べき乗法のアルゴリズムは、次のようにまとめられます。

---

• 入力：行列 $A$、初期値ベクトル $x^{(0)}$、微小正数 $EPS$、反復回数の上限 $K$

• 出力：反復回数 $k$ と近似固有値 $r$ と近似固有ベクトル $x$、または「反復回数の上限を超えた」旨のメッセージ

• Step 1：$A$、$EPS$、$K$ を初期化する。$x=x^{(0)}$ とする。

• Step 2：$r=x^TAx$ とする。$k=0$ とする。

• Step 3：$\left\| Ax-rx\right\|\geq EPS$ かつ $k\leq K$ の間、次のStep 3-1～3-3を繰り返す。

  • Step 3-1：$x$ を更新する。

  • Step 3-2：$r$ を更新する。

  • Step 3-3：$k$ に $1$ を加算する。

• Step 4：もし $k>K$ なら、「反復回数の上限を超えた」旨のメッセージを出力する。そうでなければ、反復回数 $k$ と近似固有値 $r$ と近似固有ベクトル $x$ を出力する。

---

赤字の箇所は敢えて曖昧に書いていますので、どうしたらよいか考えてください。

演習1

べき乗法のアルゴリズムを実現するコードを書き、

$$ 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 $$

とします。

%演習1のコード

%行列の転置は「'」を使う
%2ノルムは「norm(ベクトル)」を使う

演習2

演習1と同じ行列 $A$ の固有値 $1$ に対応する固有ベクトルをeig関数または手計算で求め、求めた固有ベクトルを初期値ベクトル $x^{(0)}$ としてべき乗法を実行してください（$EPS$ と $K$ は演習1と同じにすること）。その結果になった理由について考察してください。

%演習2のコード

演習3（オプション）

上で扱ったべき乗法のアルゴリズムは、行列の積を何度も計算するためあまり計算効率が良くありません。効率が良くなるようにアルゴリズムを改良して、そのコードを書いてください（考える問題は演習1と同じにすること）。反復の終了条件を変更してもよいです。

%演習3のコード

第11回レポート課題

演習1～演習3に取り組んでください。