[勉強ノート] 「線形計算の数理」 7.2.1 べき乗法 基本形

今日の記事では「線形計算の数理」の7.2.1で紹介されているべき乗法について解説とPythonによる実装を紹介します。基本的には自分があとで見返したときに納得できるように説明を書いたので、かなり基本的なところから書いてあると思います。

実数に対応したPythonによる実装はこちらに置いてあります。

https://github.com/shu65/theoretical-numerical-linear-algebra/blob/main/%E7%B7%9A%E5%BD%A2%E8%A8%88%E7%AE%97%E3%81%AE%E6%95%B0%E7%90%86_7_2_1_%E3%81%B9%E3%81%8D%E4%B9%97%E6%B3%95_%E5%9F%BA%E6%9C%AC%E5%BD%A2.ipynb

固有値問題とは

まずはべき乗法の前に固有値問題について触れます。

複素$n$次正方行列 $A$ に対して以下のような関係を満たす 複素数の$\lambda$ と 複素ベクトル$z$ の組を全部、またはいくつか求める問題を固有値問題と言います。

$$\begin{align*} A \boldsymbol{z} = \lambda \boldsymbol{z} \tag{7.1} \end{align*}$$

このとき、 $\lambda$を固有値、$z$を固有ベクトルと言います。

この固有値問題を解くアルゴリズムの一つが今回紹介するべき乗法です。

べき乗法とは

べき乗法は固有値問題を解くアルゴリズムで、固有値のうち、絶対値最大となる固有値とその固有ベクトル を求めるアルゴリズムになります。

アルゴリズムの概要

具体的なPythonコードは後ほどしめしますが、基本的には以下のようなアルゴリズムになります。

1. $A$ と任意の単位ベクトル $\boldsymbol{x}_0$を入力に取る
2. 以下を$\boldsymbol{x}_{k}$ が収束するまで繰り返す
   1. $\boldsymbol{y}_{k} := A \boldsymbol{x}_{k-1}$
   2. $\boldsymbol{x}_{k} := \boldsymbol{y}_{k} / || \boldsymbol{y}_{k} ||$
3. 収束した$\boldsymbol{x}_{k}$ を$\boldsymbol{x}$ とし、$\boldsymbol{x}$ と $\boldsymbol{y}$ ($=A \boldsymbol{x}$) の要素比の絶対値 $| y_i / x_i |$ が最大となる $y_i / x_i$ (ただし、 $x_i \ne 0$ ) を求める。
4. 3で求めた$y_i / x_i$を絶対値最大となる固有値、その固有ベクトルを$\boldsymbol{x}$ として出力する

絶対値最大の固有値の固有ベクトルに収束する理由

先ほど紹介したアルゴリズムでは$\boldsymbol{x}_{k}$が絶対値最大となる固有値の固有ベクトルに収束することを利用していました。そもそもなぜ $\boldsymbol{x}_{k}$が固有ベクトルに収束していくのか？を示します。

まず、複素$n$次正方行列$A$の固有値 $\lambda_i$ を以下のように絶対値が大きい順に番号付けします。

$$\begin{align*} |\lambda_1| \geq |\lambda_2| \geq … \geq |\lambda_n|\tag{7.6} \end{align*}$$

また、固有ベクトルを単位ベクトルとし、それぞれの固有値に対応する固有ベクトルを固有値の絶対値が大きい順に $\boldsymbol{z}_{1}, \boldsymbol{z}_{2}, …, \boldsymbol{z}_{n}$ とします。ここで$A$ が対角化可能(つまり$n$個の1次独立な固有ベクトルを持つ)という仮定がおけるとします。つまり、$n$個の固有ベクトルを用いると$n$次元の任意のベクトルを固有値ベクトルを使って表すことができます。

これを利用して任意の単位ベクトル $\boldsymbol{x}_0$を以下のように固有ベクトルを使って表します。

$$\begin{align*} \boldsymbol{x}_0 = \sum_{i=1}^n c_i \boldsymbol{z}_i \tag{7.7.1} \end{align*}$$

ここで$c_i$ は複素数です。この(7.7.1)で$A^k$を左からかけ、(7.1)を利用すると以下のように変形できます。

$$\begin{align*} A^k \boldsymbol{x}_0 =& \sum_{i=1}^n c_i A^k \boldsymbol{z}_i \\ =& \sum_{i=1}^n c_i \lambda_i^k \boldsymbol{z}_i \\ =& \lambda_1^k \left[ c_1 \boldsymbol{z}_1 + \sum_{i=2}^n c_i \left( \frac{\lambda_i}{\lambda_1} \right)^k \boldsymbol{z}_i \right] \tag{7.7} \end{align*}$$

このとき$\lambda_1$ は絶対値最大の固有値であるため、 $|\frac{\lambda_i}{\lambda_1}| < 1$ が成り立ちます。このため、sumの項に関しては$k \to \infty$ のとき0に収束することがわかります。

この式を両辺を $|| A^k \boldsymbol{x}_0 ||$ で割り、$k \to \infty$ のときにどうなるかを見てみます。固有ベクトルが単位ベクトルということを利用すると以下のように固有ベクトル $\boldsymbol{z}_1$だけにすることができます。

$$\begin{align*} \lim_{k \to \infty} \frac{A^k \boldsymbol{x}_0}{|| A^k \boldsymbol{x}_0 ||} =& \lim_{k \to \infty} \frac{\lambda_1^k \left[ c_1 \boldsymbol{z}_1 + \sum_{i=2}^n c_i \left( \frac{\lambda_i}{\lambda_1} \right)^k \boldsymbol{z}_i \right]}{|| \lambda_1^k \left[ c_1 \boldsymbol{z}_1 + \sum_{i=2}^n c_i \left( \frac{\lambda_i}{\lambda_1} \right)^k \boldsymbol{z}_i \right] ||} \\ =& \lim_{k \to \infty} \frac{\lambda_1^k \left[ c_1 \boldsymbol{z}_1 + \sum_{i=2}^n c_i \left( \frac{\lambda_i}{\lambda_1} \right)^k \boldsymbol{z}_i \right]}{\lambda_1^k || \left[ c_1 \boldsymbol{z}_1 + \sum_{i=2}^n c_i \left( \frac{\lambda_i}{\lambda_1} \right)^k \boldsymbol{z}_i \right] ||} \\ =& \lim_{k \to \infty} \frac{\left[ c_1 \boldsymbol{z}_1 + \sum_{i=2}^n c_i \left( \frac{\lambda_i}{\lambda_1} \right)^k \boldsymbol{z}_i \right]}{|| \left[ c_1 \boldsymbol{z}_1 + \sum_{i=2}^n c_i \left( \frac{\lambda_i}{\lambda_1} \right)^k \boldsymbol{z}_i \right] ||} \\ =& \frac{c_1 \boldsymbol{z}_1}{||c_1 \boldsymbol{z}_1||} \\ =& \frac{c_1 \boldsymbol{z}_1}{c_1 ||\boldsymbol{z}_1||} \\ =& \boldsymbol{z}_1 \tag{7.7.2} \end{align*}$$

この式の左辺の $\frac{A^k \boldsymbol{x}_0}{|| A^k \boldsymbol{x}_0 ||}$ は$\boldsymbol{x}_k$ の定義から $\frac{A^k \boldsymbol{x}_0}{|| A^k \boldsymbol{x}_0 ||} = \boldsymbol{x}_k$です。このため、$\boldsymbol{x}_k$ が収束すると $\boldsymbol{z}_1$ の近似ベクトルになります。

固有ベクトルから固有値を求める

次に固有値の求め方を説明します。先ほどの説明で $\boldsymbol{x}_k$ が収束すると固有値ベクトル $\boldsymbol{z}_1$ に収束することがわかりました。このため、収束した$\boldsymbol{x}_k$を $\boldsymbol{x}$ とし、$\boldsymbol{y} = A\boldsymbol{x}$ と置きます。この時、$\boldsymbol{x}$ は固有ベクトルの近似なので以下の関係が成り立ちます。

$$\begin{align*} \boldsymbol{y} = A \boldsymbol{x} \approx \lambda_1 \boldsymbol{x} \tag{7.7.3} \end{align*}$$

このため、$\boldsymbol{x}$ と$\boldsymbol{y}$の要素をそれぞれ $x_i$、$y_i$とするとき以下の式が成り立ちます。

$$\begin{align*} y_i \approx \lambda_1 x_i \tag{7.7.4} \end{align*}$$

式 (7.7.4)から$x_i \ne 0$ のときは以下の関係が成り立ちます。

$$\begin{align*} \lambda_1　\approx \frac{y_i}{x_i} \tag{7.7.5} \end{align*}$$

本では詳しく理由は書かれていないですが、$|\frac{y_i}{x_i}|$ が最大となるような $i$ を選ぶとよいとあります。

また$A$ がエルミート行列の時は以下の式(レイリー商)のほうが良い近似らしいです。

$$\begin{align*} \lambda_1　\approx \frac{\boldsymbol{x}^H A \boldsymbol{x}}{\boldsymbol{x}^H \boldsymbol{x}} \tag{7.9} \end{align*}$$

他の方のblog記事を読むとレイリー商を利用して固有値を出しているケースが多い印象でした。

Pythonによる実装

アルゴリズムの概要で示したものをPythonで具体的に書きます。今回は複素数ではなく実数のみ対応しています。また、簡単な行列演算はnumpyを使用しています。

def power_method(a, eps=1e-7, max_iterations=100):
  x0 = np.ones((a.shape[1],))/np.sqrt(a.shape[1])
  x_k = x0
  for i in range(max_iterations):
    y_k = a @ x_k
    x_kp1 = y_k / np.linalg.norm(y_k)
    relative_residual_norm = np.linalg.norm(x_kp1 - x_k)
    if relative_residual_norm < eps:
      break
    x_k = x_kp1 
  x = x_k
  y = a @ x
  eigenvalue = None
  for x_i, y_i in zip(x, y):
    if np.abs(x_i) >= eps:
      tmp = y_i / x_i
      if (eigenvalue is None) or (np.abs(eigenvalue) < np.abs(tmp)):
        eigenvalue = tmp
  eigenvector = x
  return eigenvalue, eigenvector

このコードを実行すると以下のようになります。

# 入力のAの行列
a_matrix = np.array(
    [[1, 0],
     [0, 2]]
)
eigenvalue, eigenvector = power_method(a_matrix)
print("eigenvalues:", eigenvalue)
print("eigenvectors:")
print(eigenvector)

# 出力
eigenvalues: 2.0
eigenvectors:
[1.1920929e-07 1.0000000e+00]

numpyの固有値と固有ベクトルを求める np.linalg.eig()でも近い値が得られているので想定通りの動作をしていると考えています。

終わりに

最近固有値問題について勉強していて、同僚に良い本として「線形計算の数理」を教えてもらったので読んでいました。この際、アルゴリズムの背景にある数学の証明が面白いと感じたのと、他の方のべき乗法の記事を読んでもわからないところがあったので、自分があとで見返したときに納得できるように説明を書きました。

実は「線形計算の数理」で紹介されている固有値問題のアルゴリズムは大部分を勉強がてら実装してあるので、随時今回のような記事にしていければと思っています。