4.2.QR法の概要

QR法は、以下に示すn次正方行列\( A \)の三角化を繰り返すことで、最終的に対角要素が固有値となる上三角行列に収束する性質を利用します。 \( A \)の三角化には、ユニタリ行列\( Q \)による相似変換を用います。

<QR法のステップ>

  1. \( A^{(0)} = A \)とします。
  2. \( A^{(0)} = Q_0 R_0 \)から\( Q_0 \)、\( R_0 \)を算出します(QR分解)。
    \( \rightarrow R_0 = Q_0^t A^{(0)} = Q_0^t A \)
  3. \( R_0 Q_0 = A^{(1)} \ \)を計算します。
    \( \rightarrow A^{(1)} = Q_0^t A Q_0 \)
  4. \( A^{(1)} = Q_1 R_1 \)から\( Q_1 \)、\( R_1 \)を算出します(QR分解)。
    \( \rightarrow R_1 = Q_1^t A^{(1)} = Q_1^t Q_0^t A Q_0 \)
  5. \( R_1 Q_1 = A^{(2)} \ \)を計算します。
    \( \rightarrow A^{(2)} = Q_1^t Q_0^t A Q_0 Q_1 \)
  6. ・・・
  7. \( A^{(k-1)} = Q_{k-1} R_{k-1} \ \)から\( Q_{k-1} \)、\( R_{k-1} \ \)を算出します(QR分解)。
    \( \begin{eqnarray} & \rightarrow R_{k-1} & = Q_{k-1}^t A^{(k-1)} \\ & & = Q_{k-1}^t \cdots Q_0^t A Q_0 \cdots Q_{k-2} \end{eqnarray} \)
  8. \( R_{k-1} Q_{k-1} = A^{(k)} \ \)を計算します。
    \( \rightarrow A^{(k)} = Q_{k-1}^t \cdots Q_0^t A Q_0 \cdots Q_{k-1} \)
  9. ・・・

これを\( k \rightarrow \infty \)とすることで、\( A \rightarrow \)上三角行列に収束します。
なお、ユニタリ行列\( Q \)は次の性質を持ちます。
  • \( Q^t = Q^{-1} \)(逆行列と転置行列は等しい)
  • \( Q \)の各列ベクトルは大きさ1で、それぞれが直交する
この計算は、\( k \rightarrow \infty \)で\( A^{(k)} \)の左下半分が0に収束し上三角化され、対角要素が\( A \)の固有値に収束します。 実際の数値計算(コンピュータでの計算)では∞の繰り返しはできないため、\( A^{(k)} \)の変化が収束判定値以下になったところを近似解として採用します。

ここからは、QR法の中身について話を進めます。
QR法は次の手順で行います。

(1)行列\( A \)を、ハウスホルダー行列\( H \)の相似変換を用いてヘッセンベルグ化します。

なぜヘッセンベルグ化が必要か?というと、行列\( A \)が大規模なときQR分解の時間短縮に有効とされるためです。

\( B = H^t A H \)(\( H \):ハウスホルダー行列)
ヘッセンベルグ型: \( \begin{pmatrix} & b_{11} & b_{12} & b_{13} & \cdots & b_{1n} & \\ & b_{21} & b_{22} & b_{23} & \cdots & b_{2n} & \\ & 0 & b_{32} & b_{33} & \cdots & b_{3n} & \\ & \vdots & \ddots & \ddots & \ddots & \vdots & \\ & 0 & \cdots & 0 & b_{n-1n} & b_{nn} & \end{pmatrix} \)

(2)ヘッセンベルグ化した行列\( B \)を、上記<QR法のステップ>で上三角化します。

このとき、\( B^{(1)}, B^{(2)}, \cdots \)もまたヘッセンベルグ行列になります(証明略)。 ヘッセンベルグ行列をQR分解する際は、ギブンス行列\( G \)を用いて、以下の手順で行います (k回目を例にとります)。

<QR分解>

  1. 2行1列目の要素を0にする:\( G_1^t \ B^{(k-1)} = B_1 \)
  2. 3行2列目の要素を0にする:\( G_2^t \ B_1 = G_2^t \ G_1^t B^{(k -1)} = B_2 \)
  3. ・・・
  4. n行n-1列目の要素を0にする:\( G_{n-1}^t \ B_{n-2} \ \)\( = G_{n-1}^t \cdots \ G_1^t \ B^{(k -1)} \ \)\(= R_k \)

このときギブンス行列\( G_i \)はユニタリ行列であることから、\( G_1 \cdots G_{n-1} = Q_k \ \)とでき(ユニタリ行列)、
\[ G_{n-1}^t \cdots \ G_1^t \ B^{(k -1)} = Q_k^t B^{(k -1)} = R_k \]
となります。
QR分解にギブンス行列を用いるのは、\( i \)行\( i+1 \)列成分を0にするための計算を、\( i, i+1 \)の二行分のみ行えばよい、つまり計算量を減らせる利点があるためです。
例えば2行1列目を0にする場合、ギブンス変換では1,2行目のみ計算することになります(計算部分を赤字にしています)。
\[ \begin{eqnarray} & G_1^t B & = \begin{pmatrix} & \cos{\theta} & -\sin{\theta} & 0 & \cdots & \\ & \sin{\theta} & \cos{\theta} & 0 & \cdots & \\ & 0 & 0 & 1 & & \\ & \vdots & & & \ddots & \end{pmatrix} \begin{pmatrix} & b_{11} & b_{12} & b_{13} & \cdots & \\ & b_{21} & b_{22} & b_{23} & \cdots & \\ & 0 & b_{32} & b_{33} & & \\ & \vdots & & \ & \ddots & \end{pmatrix} \\ & & = \begin{pmatrix} & \txc{red}{ b_{11} \cos{\theta} - b_{21} \sin{\theta} } \quad & \txc{red}{ b_{11} \cos{\theta} - b_{21} \sin{\theta} } & & \txc{red}{ \cdots } & \\ & \txc{red}{ b_{11} \sin{\theta} + b_{21} \cos{\theta} } \quad & \txc{red}{ b_{11} \sin{\theta} + b_{21} \cos{\theta} } & & \txc{red}{ \cdots } & \\ & 0 & b_{32} & b_{33} & & \\ & \vdots & & & \ddots & \end{pmatrix} \end{eqnarray} \]
このとき、
\[ \cos \theta = \frac{ b_{11} }{ \sqrt{ b_{11}^2 + b_{21}^2 } } , \quad \sin \theta = \frac{ -b_{21} }{ \sqrt{ b_{11}^2 + b_{21}^2 } } \]
とすれば、目的は達せられます。
\[ \begin{eqnarray} & G_1^t B & = \begin{pmatrix} & \txc{red}{ b_{11}' } \ & \txc{red}{ b_{12}' } \ & \txc{red}{ b_{13}' } \ & \txc{red}{ \cdots } & \\ & \txc{red}{ 0 } \ & \txc{red}{ b_{22}' } \ & \txc{red}{ b_{23}' } \ & \txc{red}{ \cdots } & \\ & 0 \ & b_{32} \ & b_{33} & & \\ & \vdots \ & \ & \ & \ddots & \end{pmatrix} \end{eqnarray} \]


ハウスホルダー行列、ギブンス行列ともにユニタリ行列なので、ヘッセンベルグ行列\( B \) と\( A \)の固有値は一致します。
以上がQR法による固有値算出の概要になります。

4.3.QR法の計算フロー

QR法の計算フローは次の通りです。

QR法の計算フロー
図4.3-1 QR法の計算フロー

参考文献

固有値計算は、CAE、CFD、統計解析、機械学習など実用的にも重要な役割を担っています。