固有値計算
4．2．QR法の概要
QR法は、以下に示すn次正方行列\( A \)の三角化を繰り返すことで、最終的に対角要素が固有値となる上三角行列に収束する性質を利用します。 \( A \)の三角化には、ユニタリ行列\( Q \)による相似変換を用います。＜QR法のステップ＞
- \( A^{(0)} = A \)とします。
- \( A^{(0)} = Q_0 R_0 \)から\( Q_0 \)、\( R_0 \)を算出します（QR分解）。
\( \rightarrow R_0 = Q_0^t A^{(0)} = Q_0^t A \) - \( R_0 Q_0 = A^{(1)} \ \)を計算します。
\( \rightarrow A^{(1)} = Q_0^t A Q_0 \) - \( 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 \) - \( R_1 Q_1 = A^{(2)} \ \)を計算します。
\( \rightarrow A^{(2)} = Q_1^t Q_0^t A Q_0 Q_1 \) - ・・・
- \( 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} \) - \( 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} \) - ・・・
- \( Q^t = Q^{-1} \)（逆行列と転置行列は等しい）
- \( Q \)の各列ベクトルは大きさ1で、それぞれが直交する
ここからは、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分解＞
- 2行1列目の要素を0にする：\( G_1^t \ B^{(k-1)} = B_1 \)
- 3行2列目の要素を0にする：\( G_2^t \ B_1 = G_2^t \ G_1^t B^{(k -1)} = B_2 \)
- ・・・
- 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}
\]
以上がQR法による固有値算出の概要になります。
4．3．QR法の計算フロー
QR法の計算フローは次の通りです。
図4.3-1 QR法の計算フロー
