べき乗法・逆反復法の概要とQR法・逆反復法を用いた固有値・固有ベクトルの数値計算

QR法と逆反復法を組み合わせることで固有値・固有ベクトルの近似計算を行うことができます。QR法は別記事で取り扱ったので、当記事ではべき乗法や逆反復法に関して取りまとめ、QR法で得られた固有値を元に逆反復法を用いて固有ベクトルの近似計算を行います。

・参考
数値計算講義ノート$6$ 固有値問題の解法

逆反復法の理解

べき乗法

べき乗法は$n$次正方行列$A$の最大固有値と最大固有値に対応する固有ベクトルを計算する手法である。べき乗法の仕組みを理解するにあたっては、$Ax_{i} = \lambda_{i} x_{i}$の式と任意の大きさ$1$のベクトル$y$が大きさ$1$の固有ベクトル$x_i$を用いて下記のように表されることに基づいて考えればよい。
$$
\large
\begin{align}
y &= c_1 x_1 + \cdots c_n x_n \\
&= \sum_{i=1}^{n} c_i x_i
\end{align}
$$

上記に対して$Ay$を計算すると以下のように計算できる。
$$
\large
\begin{align}
Ay &= A(c_1 x_1 + \cdots c_n x_n) \\
&= c_1 A x_1 + \cdots + c_n A x_n \\
&= c_1 \lambda_1 x_1 + \cdots + c_n \lambda_n x_n
\end{align}
$$

上記に基づいて、$A^{k}y$は下記のように計算できる。
$$
\large
\begin{align}
A^{k}y &= A^{k-1} (c_1 \lambda_1 x_1 + \cdots + c_n \lambda_n x_n) \\
&= A^{k-2} (c_1 \lambda_1^2 x_1 + \cdots + c_n \lambda_n^2 x_n) \\
&= \cdots \\
&= A (c_1 \lambda_1^{k-1} x_1 + \cdots + c_n \lambda_n^{k-1} x_n) \\
&= c_1 \lambda_1^k x_1 + \cdots + c_n \lambda_n^k x_n \quad (1)
\end{align}
$$

ここで$|\lambda_1|>|\lambda_2|>\cdots|\lambda_{n}|$が成立すると考える場合、ベクトル$A^{k}y$は最大固有値$\lambda_1$に対応する$x_1$の方向に収束すると考えられる。このことは$\displaystyle \lim_{k \to \infty} \frac{A^{k}y}{|A^{k}y|}$を下記のように考えることで示せる。
$$
\large
\begin{align}
\lim_{k \to \infty} \frac{A^{k}y}{|A^{k}y|} &= \lim_{k \to \infty} \frac{c_1 \lambda_1^k x_1 + \cdots + c_n \lambda_n^k x_n}{(c_1^2 \lambda_1^{2k} + \cdots + c_n^2 \lambda_n^{2k})^{1/2}} \\
&= \lim_{k \to \infty} \frac{c_1 x_1 + \cdots + c_n (\lambda_n/\lambda_1)^k x_n}{(c_1^2 + \cdots + c_n^2 (\lambda_n/\lambda_1)^{2k})^{1/2}} \\
&= \frac{c_1 x_1}{c_1} = x_1 \quad (2)
\end{align}
$$

上記の計算にあたっては「$|x_i|=1$かつ$x_i$が固有ベクトルである $\iff$ $x_i$が正規直交基底である」ことに基づいて計算を行なった。ここまでの内容に基づいて下記のような繰り返し演算で最大固有値$\lambda_{1}$に対応する固有ベクトル$x_1$を求めることができる。
$$
\large
\begin{align}
y_0 &= \frac{u}{|u|} \\
y_{k+1} &= \frac{Ay_k}{|Ay_k|} \\
\lim_{k \to \infty} y_{k} &= x_1
\end{align}
$$

上記では$\displaystyle y_{k+1} = \frac{Ay_k}{|Ay_k|}$のように繰り返しごとに$|Ay_k|$で割り算を行うが、ベクトルの方向が変わらないことから$(2)$式と同じ結果が得られる。

逆反復法

固有値・固有ベクトルに関する$Ax = \lambda x$の式は下記のように変形を行うことができる。
$$
\large
\begin{align}
Ax &= \lambda x \\
A^{-1}Ax &= \lambda A^{-1} x \\
A^{-1} x &= \frac{1}{\lambda} x
\end{align}
$$

上記より$A^{-1}$の固有値$1/\lambda$に対する固有ベクトルが$x$であると考えられる。この式にべき乗法を適用すると$1/\lambda_{i}$が最大、すなわち$\lambda_{i}$が最小になる固有値に対応する固有ベクトルが得られる。

ここでQR法などによって$A$の固有値の近似値$\hat{\lambda}_{i}$が得られているとき、$A-\hat{\lambda}_{i}I$に関して下記が成立する。
$$
\large
\begin{align}
(A-\hat{\lambda}_{i}I) x_i &= Ax_i – \hat{\lambda}_{i} x_{i} \\
&= (\lambda_{i} – \hat{\lambda}_{i}) x_{i} \\
(A-\hat{\lambda}_{i}I)^{-1} x_{i} &= \frac{1}{\lambda_{i} – \hat{\lambda}_{i}} x_{i}
\end{align}
$$

上記にべき乗法を適用すると、$\hat{\lambda}_{i}$に対応する固有ベクトルを得ることができる。上記は漸化式を用いて下記のように表記できる。
$$
\large
\begin{align}
y_0 &= \frac{u}{|u|} \\
y_{k+1} &= \frac{(A-\hat{\lambda}_{i}I)^{-1}y_k}{|(A-\hat{\lambda}_{i}I)^{-1}y_k|} \\
\lim_{k \to \infty} y_{k} &= x_i
\end{align}
$$

上記の逆行列の計算はLU分解を元に行うことができる。詳細は下記で取り扱った。

また、LU分解を用いた「逆行列」の計算は「連立方程式の解」の計算に基づいて計算を行う。よって、下記のように連立方程式の解を解くことによって$x_i$を計算することもできる。
$$
\large
\begin{align}
y_0 &= \frac{u}{|u|} \\
(A-\hat{\lambda}_{i}I)z_{k+1} &= y_{k}, \quad y_{k+1} = \frac{z_{k}}{|z_{k+1}|} \\
\lim_{k \to \infty} y_{k} &= x_i
\end{align}
$$

上記の計算は$y_{k}$が既知のとき、連立方程式の$(A-\hat{\lambda}_{i}I)z_{k+1} = y_{k}$を$z_{k+1}$に関して解くにあたって、LU分解などを用いることを意味する。詳しくは「LU分解を用いた連立方程式の解法」で取りまとめた。

固有値・固有ベクトルの計算

べき乗法を用いた最大固有値と対応する固有ベクトルの計算

行列$\displaystyle A = \left(\begin{array}{cc} 2 & 1 \\ 1 & 2 \end{array} \right)$に対してべき乗法を適用すると下記のように計算を行うことができる。

import numpy as np

A = np.array([[2., 1.], [1., 2.]])
u = np.array([2., 1.])

y_k = u
for i in range(10):
    y_k = np.dot(A,y_k)/np.linalg.norm(np.dot(A,y_k))
    print(y_k)

・実行結果

[ 0.78086881  0.62469505]
[ 0.73279349  0.6804511 ]
[ 0.71578195  0.69832385]
[ 0.71001067  0.70419091]
[ 0.70807608  0.70613615]
[ 0.70743003  0.70678338]
[ 0.70721455  0.706999  ]
[ 0.70714271  0.70707086]
[ 0.70711876  0.70709481]
[ 0.70711077  0.70710279]

上記はベクトル$y_{k}$の初期値を$\displaystyle y_{0} = \left(\begin{array}{c} 2 \\ 1 \end{array} \right)$とおき、繰り返し計算を行なったが、$y_{0}$が他の値でも収束することを以下、確認する。

import numpy as np

A = np.array([[2., 1.], [1., 2.]])
u = np.array([[2., 1.], [1., 2.], [1., 0.], [2., 2.], [-1., 2]])

for i in range(u.shape[0]):
    y_k = u[i,:]
    for j in range(10):
        y_k = np.dot(A,y_k)/np.linalg.norm(np.dot(A,y_k))
    print("・v_0: {}".format(u[i,:]))
    print("eigen value: {:.1f}".format(np.dot(y_k,np.dot(A,y_k))))
    print("eigen vector: {}".format(y_k))

・実行結果

・v_0: [ 2.  1.]
eigen value: 3.0
eigen vector: [ 0.70711077  0.70710279]
・v_0: [ 1.  2.]
eigen value: 3.0
eigen vector: [ 0.70710279  0.70711077]
・v_0: [ 1.  0.]
eigen value: 3.0
eigen vector: [ 0.70711876  0.70709481]
・v_0: [ 2.  2.]
eigen value: 3.0
eigen vector: [ 0.70710678  0.70710678]
・v_0: [-1.  2.]
eigen value: 3.0
eigen vector: [ 0.70707086  0.70714271]

大きさ$1$の固有ベクトル$y_k$に対応する固有値が$\lambda = y_{k}^{\mathrm{T}}Ay_{k}$で求められることを用いて、上記では固有値の計算も合わせて行なった。$\lambda = y_{k}^{\mathrm{T}}A y_{k}$の式は$Ay_{k} = \lambda y_{k}$の両辺に左から$y_{k}^{\mathrm{T}}$をかけることによって下記のように導出できる。
$$
\large
\begin{align}
Ay_{k} &= \lambda y_{k} \\
y_{k}^{\mathrm{T}} A y_{k} &= \lambda y_{k}^{\mathrm{T}} y_{k} \\
&= \lambda
\end{align}
$$

QR法を用いた固有値の計算

詳細は上記でまとめたが、QR法を用いて固有値は下記のように計算できる。

import numpy as np

def QR_gram(X):
    U, Q, R = np.zeros([X.shape[0],X.shape[0]]), np.zeros([X.shape[0],X.shape[0]]), np.zeros([X.shape[0],X.shape[1]])
    U[:,0] = X[:,0]
    Q[:,0] = U[:,0]/np.linalg.norm(U[:,0])
    R[0,0] = np.dot(X[:,0],U[:,0]/np.linalg.norm(U[:,0]))
    for i in range(1,X.shape[0]):
       U[:,i] = X[:,i]
       for j in range(i):
           U[:,i] += -(np.dot(X[:,i],U[:,j])/np.dot(U[:,j],U[:,j]))*U[:,j]
           R[j,i] = np.dot(X[:,i],U[:,j]/np.linalg.norm(U[:,j]))
       Q[:,i] = U[:,i]/np.linalg.norm(U[:,i])
       R[i,i] = np.dot(X[:,i],U[:,i]/np.linalg.norm(U[:,i]))
    return Q, R

A = np.array([[2., 1.], [1., 2.]])
Q, R = QR_gram(A)

A_k = A
Q_k, Q_k_inv = Q, Q.T
for i in range(10):
    Q_k, _ = QR_gram(A_k)
    Q_k_inv = Q_k.T
    A_k = np.dot(Q_k_inv,np.dot(A_k,Q_k))
    
Lambda_A = A_k
print("lambda_1: {:.1f}".format(Lambda_A[0,0]))
print("lambda_2: {:.1f}".format(Lambda_A[1,1]))

・実行結果

> print("lambda_1: {:.1f}".format(Lambda_A[0,0]))
lambda_1: 3.0
> print("lambda_2: {:.1f}".format(Lambda_A[1,1]))
lambda_2: 1.0

逆反復法を用いた固有ベクトルの計算

前項の「QR法」の計算結果を元に、下記のように逆反復法の計算を行うことができる。

lamb = np.array([Lambda_A[0,0], Lambda_A[1,1]])
n = lamb.shape[0]
u = np.array([1., 0.])

for i in range(n):
    A_ = A - lamb[i]*np.eye(n)
    inv_A_ = np.linalg.inv(A_)
    y_k = u/np.linalg.norm(u)
    for j in range(10):
        y_k = np.dot(inv_A_,y_k)/np.linalg.norm(np.dot(inv_A_,y_k))
    print("・eigen value: {:.1f}".format(lamb[i]))
    print("eigen vector: {}".format(y_k))

・実行結果

・eigen value: 3.0
eigen vector: [ 0.70710678  0.70710678]
・eigen value: 1.0
eigen vector: [ 0.70710678 -0.70710678]

上記の計算にあたっては簡易化にあたってnp.linalg.invを用いたが、逆行列の計算に関しては下記で取りまとめたLU分解などを用いることで行うことができる。