ハウスホルダー行列について
ハウスホルダー行列のおさらいから話を始めよう。
ハウスホルダー行列とはエルミート性とユニタリー性を兼ね備えた行列だった。
主な用途として、たとえばこのハウスホルダー行列により、ある正方行列が対称であれば三重対角行列に,非対称であればヘッセンベルグ行列にそれぞれ変置換することができる。
・三重対角行列：例
\begin{eqnarray}
X = \left(
\begin{array}{cccc}
x_{ 11 } & x_{ 12 } & 0 & \ldots & 0 & 0 \\
x_{ 21 } & x_{ 22 }& x_{23} & \ldots & 0 & 0 \\
0 & x_{ 32 }& x_{33} & \ldots & 0 & 0 \\
\vdots & \vdots & \vdots& \ddots & \vdots & \vdots \\
0 & 0 & 0 & \ldots & x_{ n-1,n-1 } & x_{n-1,n} \\
0 & 0 & 0 & \ldots & x_{ n,n-1 } & x_{n,n}
\end{array}
\right)
\end{eqnarray}
・ヘッセンベルグ行列：例
\begin{eqnarray}
Y = \left(
\begin{array}{cccc}
y_{ 11 } & y_{ 12 } & y_{13} & \ldots & y_{1,n-1} & y_{1,n} \\
y_{ 21 } & y_{ 22 } & y_{23} & \ldots & y_{2,n-1} & y_{2,n}\\
0 & y_{ 32 }& y_{33} & \ldots & y_{3,n-1} & y_{3,n} \\
\vdots & \vdots & \vdots& \ddots & \vdots & \vdots \\
0 & 0 & 0 & \ldots & y_{ n-1,n-1 } & y_{n-1,n} \\
0 & 0 & 0 & \ldots & y_{ n,n-1 } & y_{n,n}
\end{array}
\right)
\end{eqnarray}
またハウスホルダー行列のユニタリー性を利用して、ある行列$A$をユニタリー行列$Q$と上三角行列$R$に分解することができる。これをQR分解という。
A = QR
さらには、このQR分解を繰り返すことで、非特異行列$A$の固有値の数値計算をすることができる。これをQR法という。
このQR法についての概要が本記事の趣旨になる。
以後、行列$A$をフルランクの$\mathbb{C^{n×n}}$正規行列とし、
行列$A$の固有値について
|λ_{1}|>|λ_{2}|>…>|λ_{n}|>0
が成り立っているものとする。
QR法による固有値計算(1)
まず$A$をQR分解しよう。これを$A_{0}$とする。
A=A_{0} = Q_{0}R_{0} ...(_{1-1})\\
$Q_{0}$,$R_{0}$については以下の通りとなる。
Q_{0} = A_{0}R^{-1}_{0}...(_{1-1-a}) \\
R_{0} = Q^{H}_{0}A_{0} ...(_{1-1-b})
次に$R_{0}$,$Q_{0}$により次のような行列を作る。これを$A_{1}$としよう。
さらにこの$A_{1}$を$Q_{1}$,$R_{1}$にQR分解する。
A_{1} = R_{0}Q_{0}=Q_{1}R_{1}..._{(1-2)}\\
$Q_{1}$,$R_{1}$については以下の通りであり、
Q_{1} = A_{1}R^{-1}_{1}...(_{1-2-a}) \\
R_{1} = Q^{H}_{1}A_{1} ...(_{1-2-b})
また式(1-1-a),(1-1-b),(1-2)より次式が得られる
A_{1} = Q^{H}_{0}A_{0}Q_{0} ..._{(1-2-c)}\\
A_{1} = R_{0}A_{0}R^{-1}_{0}..._{(1-2-d)} \\
先ほどと同様に$R_{1}$,$Q_{1}$より行列$A_{2}$を作り、これを$Q_{2}$,$R_{2}$にQR分解しよう。
A_{2} = R_{1}Q_{1}=Q_{2}R_{2}..._{(1-3)}\\
$Q_{2}$,$R_{2}$については以下の通りであり、
Q_{2} = A_{2}R^{-1}_{2}...(_{1-3-a}) \\
R_{2} = Q^{H}_{2}A_{2} ...(_{1-3-b})
また、式(1-3)について$\tilde{Q_{2}}=Q_{0}Q_{1}$ $\tilde{R_{2}}=R_{1}R_{0}$とした上で、式(1-2-a),(1-2-b),(1-2-c),(1-2-d)を踏まえれば、
A_{2} = Q^{H}_{1}A_{1}Q_{1}=Q^{H}_{1}Q^{H}_{0}A_{0}Q_{0}Q_{1}=\tilde{Q^{H}_{2}}A_{0}\tilde{Q_{2}}..._{(1-3-c)}\\
A_{2} = R_{1}A_{1}R^{-1}_{1}=R_{1}R_{0}A_{0}R^{-1}_{0}R^{-1}_{1}=\tilde{R_{2}}A_{0}\tilde{R^{-1}_{2}}..._{(1-3-d)}\\
がそれぞれ得られる。
同様の手順により、
$\tilde{Q_{N}}=Q_{0}Q_{1}...Q_{N-1}$
$\tilde{R_{N}}=R_{N-1}...R_{1}R_{0}$
と置くことで、
\begin{align}
A_{N} &= R_{N-1}Q_{N-1} = Q_{N}R_{N} ..._{(1-4-a)}\\
A_{N} &= R_{N-1}Q_{N-1}=Q^{H}_{N-1}A_{N-1}Q_{N-1}=...=\tilde{Q^{H}_{N}}A_{0}\tilde{Q_{N}}
..._{(1-4-b)}\\
A_{N} &= R_{N-1}Q_{N-1}=R_{N-1}A_{N-1}R^{-1}_{N-1}=...=\tilde{R_{N}}A_{0}\tilde{R^{-1}_{N}}
..._{(1-4-c)}\\
\end{align}
がそれぞれ得られる。
さて、この$Q_{t}$について、ユニタリー行列$X_{t}$により、次の関係が成り立つとしよう。
※なお、$X_{0}=I$とする。
Q_{t} = X^{H}_{t}X_{t+1} ..._{(1-5)}\\
(t=0,1,…N-1)
式(1-5)を踏まえて$\tilde{Q_{N}}=Q_{0}Q_{1}...Q_{N-1}$を計算すれば、
\begin{align}
\tilde{Q_{N}}&=Q_{0}Q_{1}...Q_{N-1}\\
&=(X^{H}_{0}X_{1})(X^{H}_{1}X_{2})...(X^{H}_{N-1}X_{N})\\
&=X^{H}_{0}X_{N}..._{(1-6)}\\
\end{align}
が得られる。
また$A_{0}$および$X_{0}$より$\tilde{A}=X_{0}A_{0}X^H_{0}$ となる行列$\tilde{A}$を利用することで式(1-4-b),式(1-6)より次式が得られる。
\tilde{A}^{N}=(X_{0}A_{0}X^H_{0})^{N}=X_{0}A^{N}_{0}X^{H}_{0}\\
\Leftrightarrow A^{N
}_{0}=X^{H}_{N}\tilde{A}^{N}X_{N}..._{(1-7)}\\
\begin{align}
A_{N} &=\tilde{Q^{H}_{N}}A_{0}\tilde{Q_{N}}=X^{H}_{N}X_{0}A_{0}X^{H}_{0}X_{N}\\
&=X^{H}_{N}\tilde{A}X_{N}
..._{(1-8)}\\
\end{align}
ここで$\tilde{Q_{t}}\tilde{R_{t}}_{(t=0,1,…N)}$について以下のような関係が成り立っていることを確認しておこう。
\begin{align}
\tilde{Q}_{1}\tilde{R}_{1}&=Q_{0}R_{0}=A_{0}\\
\tilde{Q}_{2}\tilde{R}_{2}&=Q_{0}Q_{1}R_{1}R_{0}=Q_{0}A_{1}R_{0}=(Q_{0}R_{0})^2=A^2_{0}\\
\tilde{Q}_{3}\tilde{R}_{3}&=Q_{0}Q_{1}Q_{2}R_{2}R_{1}R_{0}=Q_{0}Q_{1}A_{2}R_{1}R_{0}\\&=Q_{0}(Q_{1}R_{1})^2R_{0}=Q_{0}A^2_{1}R_{0}=(Q_{0}R_{0})^3=A^3_{0}\\
...\\
\tilde{Q}_{N}\tilde{R}_{N}&=Q_{0}...Q_{N-1}R_{N-1}...R_{0}=(Q_{0}R_{0})^N=A^{N}_{0}\end{align}
$\tilde{Q_{N}}\tilde{R_{N}}=A^{N}_{0}$および式(1-6),(1-7)より
\tilde{Q_{N}}\tilde{R_{N}}=A^{N}_{0}=X^{H}_{0}\tilde{A^{N}}X_{0}\\
\Leftrightarrow X_{N}\tilde{R_{N}}=\tilde{A^N}X_{0}..._{(1-9)}\\
を得ることができる。
QR法による固有値計算(2)
$\tilde{A}$が正規行列であれば対角化できる。
つまり$Z$を非特異行列,$Λ$を$A$の固有値を対角成分にもつ対角行列とすれば、以下の式(2-1)が得られる。
\tilde{A}=ZΛZ^{-1}..._{(2-1)}\\
さてこの$Z$を次のように$QR$分解しよう。
Z=Q_{z}R_{z}..._{(2-2)}\\
式(2-1),(2-2)を踏まえれば、
Q^{H}_{z}\tilde{A}Q_{z}=R_{z}ΛR^{-1}_{z}
..._{(2-3)}\\
\begin{align}
\tilde{A}^{N}&=(ZΛZ^{-1})^{N}=ZΛ^{N}Z^{-1}\\
&=Q_{z}R_{z}Λ^{N}R^{-1}_{z}Q^{H}_{z}..._{(2-4)}\\
\end{align}
となる。
$L$を対角成分がすべて1からなる下三角行列,$U$を上三角行列として、$Z^{H}X_{0}$を式(2-5)のように$LU分解$する。
また$L$,$Λ^{N}$,$R_{z}$からなる$Ψ_{N}$を式(2-6)のように定める。
Z^{H}X_{0}=LU\\
\Leftrightarrow X_{0}=ZLU=Q_{z}R_{z}LU..._{(2-5)}\\
Ψ_{N} =R_{z}Λ^{N}LΛ^{-N}R^{-1}_{z}..._{(2-6)}
式(2-4),(2-5),(2-6)に基づいて$\tilde{A^{N}}X_{0}$を考えよう。
すると、
\begin{align}
\tilde{A}^{N}X_{0}&=(Q_{z}R_{z}Λ^{N}R^{-1}_{z}Q^{H}_{z})(Q_{z}R_{z}LU)\\
&=Q_{z}R_{z}Λ^{N}LU\\&=Q_{z}(R_{z}Λ^{N}LΛ^{-N}R^{-1}_{z} )(R_{z}Λ^{N}U)\\
&=Q_{z}Ψ_{N}(R_{z}Λ^{N}U)..._{(2-7)}\\
\end{align}
が得られる。
ここで$Ψ_{N}$に着目しよう。
これを$\lim_{N \to \infty} \hat{Q_{N}}=I_{n}$となるようにQR分解する。
Ψ_{N}=R_{z}(Λ^{N}LΛ^{-N})R^{-1}_{z}=\hat{Q_{N}}\hat{R_{N}}..._{(2-8)}
※筆者補足
\begin{eqnarray}
Λ^{N}LΛ^{-N} = \left(
\begin{array}{cccc}
1 & 0 & 0 & \ldots & 0 & 0 \\
l_{21}(λ_{2}/λ_{1})^{N} & 1 & 0 & \ldots & 0 & 0 \\
l_{31}(λ_{3}/λ_{1})^{N} & l_{32}(λ_{3}/λ_{2})^{N} & 1 & \ldots & 0 & 0 \\
\vdots & \vdots & \vdots& \ddots & \vdots & \vdots \\
l_{n-1,1}(λ_{n-1}/λ_{1})^{N} & l_{n-1,2}(λ_{n-1}/λ_{2})^{N} & l_{n-1,3}(λ_{n-1}/λ_{3})^{N} & \ldots & 1 & 0 \\
l_{n,1}(λ_{n}/λ_{1})^{N} & l_{n,2}(λ_{n}/λ_{2})^{N} & l_{n,3}(λ_{n}/λ_{3})^{N} & \ldots & l_{n,n-1}(λ_{n}/λ_{n-1})^{N} & 1 \\
\end{array}
\right)
\end{eqnarray}
であり、
\lim_{N \to \infty} l_{j,k}(λ_{j}/λ_{k})^{N}=0\\
_{(1 \leqq k \leqq n-1 ), (2 \leqq j \leqq n), (k<j)}
より、$\lim_{N \to \infty} Λ^{N}LΛ^{-N}= I_{n}$となる。
このため、
\lim_{N \to \infty} Ψ_{N}=\lim_{N \to \infty} R_{z}(Λ^{N}LΛ^{-N})R^{-1}_{z}\\
=R_{z}R^{-1}_{z}=I_{n}\\
\lim_{N \to \infty} Ψ_{N}=\lim_{N \to \infty} \hat{Q_{N}}\hat{R_{N}}=I_{n}\\
となる。
さて$\tilde{A}^{N}X_{0}=X_{N}\tilde{R_{N}}..._{(1-10)}$ および式(2-7),(2-8)より
\begin{align}
Q_{z}\hat{Q_{N}}\hat{R_{N}}R_{z}Λ^{N}U&=X_{N}\tilde{R_{N}}\\
\Leftrightarrow (X^{H}_{N}Q_{z}\hat{Q_{N}})(\hat{R_{N}}R_{z}Λ^{N}U\tilde{R^{-1}_{N}})&=I_{N}\\
\Leftrightarrow \breve{Q}\breve{R}&=I_{n}..._{(2-9)}
\end{align}
となる。
※$\breve{Q}$をユニタリ行列,$\breve{R}$を上三角行列として以下のように置いた。
\begin{align}
\breve{Q}&=X^{H}_{N}Q_{z}\hat{Q_{N}}\\
\breve{R}&=\hat{R_{N}}R_{z}Λ^{N}U\tilde{R^{-1}_{N}}
\end{align}
対角成分の絶対値が全て1からなる対角行列のことをユニタリー対角行列という。これを$D_{N}$としよう。
$\breve{Q}\breve{R}=I_{n}$となる場合、$\breve{Q^{H}}=\breve{R}=D_{N}$となることから、式(2-9)について、
\begin{align}
(X^{H}_{N}Q_{z}\hat{Q_{N}})D_{N}&=I_{n}\\
\Leftrightarrow Q_{z}\hat{Q_{N}}&=X_{N}D^{H}_{N}..._{(2-10)}
\end{align}
が成り立つ。そこで、式(2-9)および$\lim_{N \to \infty} \hat{Q_{N}}=I_{n}$を踏まえれば
\lim_{N \to \infty} X_{N}D^{H}_{N}=\lim_{N \to \infty} Q_{z}\hat{Q_{N}}=Q_{z}..._{(2-11)}
が得られる。
ここで式(1-8)を思い出してほしい。
A_{N} =X^{H}_{N}\tilde{A}X_{N}..._{(1-8)}\\
この式(1-8)について、左から$D_{N}$右から$D^{H}_{N}$をそれぞれかけることで、
D_{N}A_{N}D^{H}_{N}=D_{N}X^{H}_{N}\tilde{A}X_{N}D^{H}_{N}
..._{(2-12)}\\
が得ることができる。
すると式(2-3),(2-10),(2-11)より、
\lim_{N \to \infty} D_{N}A_{N}D^{H}_{N}=\lim_{N \to \infty}
D_{N}X^{H}_{N}\tilde{A}X_{N}D^{H}_{N}\\
=Q^{H}_{z}\tilde{A}Q_{z}=R_{z}ΛR^{-1}_{z}..._{(2-13)}
が成立する。
※$R_{z}ΛR^{-1}_{z}$とは以下のような上三角行列である。これを$S$と置こう。
\begin{eqnarray}
S = R_{z}ΛR^{H}_{z} = \left(
\begin{array}{cccc}
λ_{ 1 } & s_{ 12 } & s_{13} & \ldots & s_{1,n-1} & s_{1,n} \\
0 & λ_{ 2 }& s_{23} & \ldots & s_{2,n-1} & s_{2,n} \\
0 & 0 & λ_{3} & \ldots & s_{3,n-1} & s_{3,n}\\
\vdots & \vdots & \vdots& \ddots & \vdots & \vdots \\
0 & 0 & 0 & \ldots & λ_{ n-1 } & s_{n-1,n} \\
0 & 0 & 0 & \ldots & 0 & λ_{n}
\end{array}
\right)
\end{eqnarray}
式(2-13)より$\tilde{A}=Q_{z}SQ^{H}_{z}$となるが、これは$\tilde{A}$のシューア分解に他ならない.
つまるところ、QR分解を繰り返すことで、行列$A$は$Q_{z}SQ^{H}_{z}$に収束する。
この性質を利用して$A$の固有値の数値解を求めているのがQR法になるわけだ。
以上がQR法の概要となる
参考文献
・新井仁之「線形代数 基礎と応用」(日本評論社)
・柳井春夫,竹内啓「射影行列・一般逆行列・特異値分解<新装版>」 (東京大学出版会)
・森正武,杉原正顕,室田一雄「岩波講座 応用数学 方法2 線形計算」 (岩波書店)