プログラミング

HOME
プログラミング
[勉強ノート] 「線形計算の数理」7.4 QR法

3月 29, 2022 / 最終更新日時 : 3月 29, 2022 Shuji Suzuki (shu) プログラミング

[勉強ノート] 「線形計算の数理」7.4 QR法

今日の記事は前回に引き続き「線形計算の数理」の固有値問題パートの内容です。今回は7.4 「QR法」についてアルゴリズムの概要となぜそれがうまくいくのか？、Pythonによる実装の3つについて説明していきます。

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_4_QR%E6%B3%95.ipynb

1 QR法
2 QR法のアルゴリズム
3 QR法のアルゴリズムでなぜ固有値が求められるのか？
4 Pythonによる実装
5 おまけ：固有ベクトルを求めるには？
6 終わりに

QR法

QR法は複数の固有値、固有ベクトルを同時に求める手法で、前回紹介した同時反復法の計算を少し改良した手法です。同時反復法についてはこちらの記事で詳しく扱っていますので、参考にしてみてください。

[勉強ノート] 「線形計算の数理」7.3 同時反復法

アルゴリズムとしては基本的にはQR分解して出力される行列から次の行列を作り、またQR分解を行うということを収束するまで繰り返すということしていきます。

次にアルゴリズムについて詳しく説明していきます。

QR法のアルゴリズム

まずはQR法のアルゴリズムについて説明します。QR法は入力として以下のものを取ります。

$ A$: 固有値、固有ベクトルを計算する複素$ n $次正方行列

これらを入力として以下のように計算していきます。

$ A $ を入力にとり、$ A_0 $ とする
以下を$ A_k $ の対角成分が収束するまで繰り返す
1. $A_k = Q_k R_k$ にQR分解する
2. $A_{k+1} = R_kQ_k$ とする
$ A_k $ の対角成分を固有値として出力する

QR分解をご存じない方に簡単に説明するとQR分解は行列$ M $ を以下のように $Q $ と $ R $ に分解する操作です。
$$ \begin{align*}
M = QR
\end{align*} $$

ここで$ Q$ の列ベクトルは互いに独立な単位ベクトルで、Qは$ n \times m $ 行列です。( ただし、$ m \le n $ ) つまり $ Q^H Q = I_m $ が成り立ちます。また、$ R $ は$ m $ 次上三角行列です。今回は$ m= n $ となっています。

本のほうにちゃんとした言及がないのでここでは固有値だけを求めるアルゴリズムを示していますが、後ほどおまけのところに書いたように上のアルゴリズムを改良すると固有ベクトルも求めることができます。

QR法のアルゴリズムでなぜ固有値が求められるのか？

次にQR法がなぜ固有値を出力できるのか？について説明していきます。

基本的には同時反復法の$A_k$からQR法の以下の漸化式に式変形できることができることを示します。

$$ \begin{align*}
A_k = Q_k R_k \tag{7.22.1}\\
A_{k+1} = R_k Q_k \tag{7.22.2} \\
\end{align*} $$

これを示すことができれば、$A_k$の対角成分が固有値に収束する証明は同時反復法のときと同じになります。
まずは同時反復法の$A_k$は以下の通りです。

$$ \begin{align*}
A_k = X_k^H A X_k \tag{7.21.1}
\end{align*} $$

ここで、QR法の場合は$X_0 = I$ です。また同時反復法では$X_k$と$X_{k+1}$の間には以下の関係が成り立っています。

$$ \begin{align*}
A X_k= X_{k+1}R_k \tag{7.19}
\end{align*} $$

同時反復法の証明のときに示した通り、$ X_k $ と$ X_{k+1} $はユニタリ行列です。この式に対して左から$X_k^H$をかけると
$$ \begin{align*}
X_k^H A X_k= X_k^H X_{k+1}R_k \tag{7.19.1}
\end{align*} $$

ここで$Q_k$ を以下のように置きます。

$$ \begin{align*}
Q_k = X_k^H X_{k+1} \tag{7.20}
\end{align*} $$

この式(7.20)と式(7.21.1)を使うと以下のような関係が導けます。

$$ \begin{align*}
A_k =& X_k^H A X_k \\
=& X_k^H X_{k+1}R_k \\
=& Q_k R_k\tag{7.20.1}
\end{align*} $$

この式(7.20.1)は示したかった式の一つの(7.22.1)と同じです。

次に$ X_k $ と$ X_{k+1} $はユニタリ行列であるため$Q_k $はユニタリ行列であるということを利用して、式(7.20)を式変形します。

$$ \begin{align*}
X_k^H X_{k+1} = Q_k \\
X_{k+1} = X_k Q_k
\tag{7.20.2}
\end{align*} $$

これを利用すると式(7.21.1)から$A_{k+1}$と$A_k$の間に以下のような関係があることがわかります。

$$ \begin{align*}
A_{k+1} =& X_{k+1}^H A X_{k+1} \\
=& (X_k Q_k)^H A (X_k Q_k) \\
=& Q_k^H X_k^H A X_k Q_k \\
=& Q_k^H A_k Q_k \\
\tag{7.21}
\end{align*} $$

また、式(7.22.1)から以下の式が導けます。

$$ \begin{align*}
Q_k R_k =& A_k \\
R_k =& Q_k^H A_k \tag{7.22.3}
\end{align*} $$

そして式(7.21)と式(7.22.3)を利用すると

$$ \begin{align*}
A_{k+1} =& Q_k^H A_k Q_k \\
=& R_k Q_k \\
\tag{7.22.4}
\end{align*} $$

式(7.22.4)は示したかった式(7.22.2)と同じ式です。

これで同時反復法の$A_k$からQR法の以下の漸化式に式変形できることがわかりました。このため、QR法がなぜ固有値を求められるのか？という証明は同時反復法の証明と同じになります。

Pythonによる実装

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

import numpy as np


def qr_method(a, eps=1e-7, max_iterations=100):
  a_k = a
  q_km1 = np.identity(a_k.shape[0])
  for i in range(max_iterations):
    q_k, r_k = np.linalg.qr(a_k)
    relative_residual_norm = np.linalg.norm(q_km1 - q_k)
    if relative_residual_norm < eps:
      break
    a_k = r_k @ q_k
    q_km1 = q_k
  eigenvalues = np.diag(a_k)
  return eigenvalues

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

#入力Aの行列
import numpy as np


a_matrix = np.array([
      [1., 0., 0.,],
      [0., 2., 1.,],
      [0., 1., 1.],
  ])
eigenvalues = qr_method(a_matrix)
print("eigenvalues:", eigenvalues)

実行した際の出力は以下のようになります。

#出力
eigenvalues: [1.         2.61803399 0.38196601]

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

おまけ：固有ベクトルを求めるには？

本では$A_{k+1}$についてしか書かれておらず、QR法では固有ベクトルを求められないというネットの記事も見つかるので、私は最初、固有ベクトルが求められないのかと思っていました。

ただ、同時反復法では$X_k$が固有ベクトルに収束することが示されているので、$X_k$を求めることができれば固有ベクトルを求めることできるはずです。そう思って式を見直していると式(7.20.2)で以下のような関係が導けています。

$$ \begin{align*}
X_{k+1} = X_k Q_k
\tag{7.20.2}
\end{align*} $$

ここで$X_0$はQR法では単位ベクトルとなっています。つまり以下のようにアルゴリズムを変更すれば$X_k$を求め、$X_k$が固有ベクトルに収束することを利用して、固有ベクトルも出力することができます。

$ A $ を入力にとり、$ A_0 $ とする
$X_0 = I $ とする
以下を$ A_k $ の対角成分が収束するまで繰り返す
1. $A_k = Q_k R_k$ にQR分解する
2. $A_{k+1} = R_kQ_k$ とする
3. $X_{k+1} = X_k Q_k$ とする
$ A_k $ の対角成分を固有値、$X_k$を固有値に対応する固有ベクトルとして出力する

上記のアルゴリズムでは$X_k$に関係する部分が追加されています。

実際にPythonで実装すると以下のようになります。

import numpy as np


def qr_method_with_eigenvectors(a, eps=1e-7, max_iterations=100):
  a_k = a
  n = a_k.shape[0]
  eigenvectors = np.identity(n)
  for i in range(max_iterations):
    q_k, r_k = np.linalg.qr(a_k)
    eigenvectors = np.dot(q_k, eigenvectors)
    a_kp1 = r_k @ q_k
    relative_residual_norm = np.linalg.norm(np.diag(a_kp1) - np.diag(a_k))
    if relative_residual_norm < eps:
      break
    a_k = a_kp1
  eigenvalues = np.diag(a_kp1)
  return eigenvalues, eigenvectors

これを以下のように実行すると

import numpy as np

a_matrix = np.array([
      [1., 0., 0.,],
      [0., 2., 1.,],
      [0., 1., 1.],
  ])

eigenvalues, eigenvectors = qr_method_with_eigenvectors(a_matrix)
print("eigenvalues:", eigenvalues)
print("eigenvectors:")
print(eigenvectors)

出力は以下のようになります。

#出力
eigenvalues: [1.         2.61803399 0.38196601]
eigenvectors:
[[ 1.          0.          0.        ]
 [ 0.          0.85065394  0.52572604]
 [ 0.         -0.52572604  0.85065394]]

なので、正しい固有ベクトルが求められている気がします。

QR法と言った際は固有ベクトルは求めない手法を指すのか、単純に本で言及してないだけなのかはわかりませんが、少し改良すれば固有ベクトルは求められることがわかりました。