[勉強ノート] 「線形計算の数理」7.3 同時反復法
今日の記事では「線形計算の数理」の7.3で紹介されている固有値問題のアルゴリズムの同時反復法についてのアルゴリズムの説明と固有ベクトル、固有値に収束することの証明を説明します。
個人的には初見で全然収束の証明の意味がわからず、しかもネットで調べてもわかりやすい証明を見つけることができなかった部分なので、同様にわからない方の参考になれば幸いです。
証明が長いのでPython実装に関しては別の記事にします。
同時反復法とは
同時反復法は複数の固有値、固有ベクトルを同時に求める手法です。固有値問題の基本的な手法であるべき乗法は、固有値の絶対値が最大となる固有値とそれに対応する固有ベクトルを求める手法でした。一方、同時反復法はこれを拡張し、固有値の絶対値が大きい順に任意の数の固有値とそれに対応する固有ベクトルを求める手法になります。ちなみに、固有値問題のアルゴリズムとしてよく目にするQR法は同時反復法の計算を少し変えたもので、基本的な考え方は同じものになっています。
同時反復法の証明ではべき乗法と似たような流れがでてきます。べき乗法については以下の記事で説明していますので、よければ参考にしてください。
アルゴリズムの概要
まずは同時反復法のアルゴリズムについて説明します。同時反復法は入力として以下のものを取ります。
- \( A\): 固有値、固有ベクトルを計算する複素\( n \)次正方行列
- \( X_0 \): 独立な単位列ベクトルを持つ任意の\( n \times m \)行列 \(X_0^H X_0 = I_m \)
- ただし、\( m \)は求める固有値、固有ベクトルの数で \( (m \le n) \)
これらを入力として以下のように計算していきます。
- \( A \) と\( X_0 \)を入力に取る
- 以下を\( X_{k} \) が収束するまで繰り返す
- \( Y_{k + 1} := AX_k \)
- \( Y_{k + 1} = X_{k+1}R_k\) にQR分解する
- \( A_k := X_k^H A X_k\) を計算する
- \( A_k \) の対角成分を絶対値が大きい\( m \)個の固有値、\( X_k\) の列ベクトルが順に対応する固有ベクトルとして出力する
QR分解をご存じない方に簡単に説明するとQR分解は行列\( M \) を以下のように \(Q \) と \( R \) に分解する操作です。
$$ \begin{align*}
M = QR \tag{7.10.1}
\end{align*} $$
ここで\( Q\) の列ベクトルは互いに独立な単位ベクトルで、Qは\( n \times m \) 行列です。( ただし、\( m \le n \) ) つまり \( Q^H Q = I_m \) が成り立ちます。また、\( R \) は\( m \) 次上三角行列です。
こちらのアルゴリズム、私は初見でなんでうまくいくのかわからなかったので、次から頑張って説明します。
同時反復法でなぜ固有ベクトルが求められるのか？
ここから同時反復法によって固有ベクトルが求められる証明をやっていきます。まず、固有値を絶対値が大きい順に番号を振り、以下のようにします。
$$ \begin{align*}
|\lambda_1| \geq |\lambda_2| \geq … \geq |\lambda_n|\tag{7.11}
\end{align*} $$
また、固有ベクトルを単位ベクトルとし、それぞれの固有値に対応する固有ベクトルを固有値の絶対値が大きい順に \( \boldsymbol{z}_{1}, \boldsymbol{z}_{2}, …, \boldsymbol{z}_{n} \) とします。そして、\(A\) が対角化可能(つまり\(n\)個の1次独立な固有ベクトルを持つ)という仮定がおけるとします。つまり、\(n\)個の固有ベクトルを用いると、\(n\)次元の任意のベクトルを固有値ベクトルを使って表すことができます。
また、ここからある行列の列ベクトルに注目することが多くなるので、ある行列Xの\( j \) 列目のベクトルを \( X_{:, j}\) で表します。例えば、固有ベクトルを列ベクトルとして並べた行列\( Z = [\boldsymbol{z}_1,…, \boldsymbol{z}_m] \)の\(j\)列目は\(Z_{:,j}=\boldsymbol{z}_j \) となります。
さて、まずはアルゴリズムを繰り返すと\( X_k \) がどのような行列になるかを追ってみます。アルゴリズムより、以下のような漸化式が成り立つことがわかります。
$$ \begin{align*}
AX_k = X_{k+1}R_k \tag{7.10.2}
\end{align*} $$
この漸化式を使って\(A^kX_0\)を\(k=1\)から順番に見ていくと以下のようになることがわかります。
$$ \begin{align*}
AX_0 =& X_1R_0 \\
A^2X_0 =& A(AX_0) = AX_1R_0 = X_2R_1R_0 \\
…&\\
A^kX_0 =& A(X_{k-1}) = AX_{k-1}R_{k-2}…R_0 = X_{k}R_{k-1}…R_0 \\
=& X_{k} \Gamma_k \tag{7.10}
\end{align*} $$
ここで\( \Gamma_k = R_{k-1}…R_0\) と置きました。\( X_{k}\) と\(\Gamma_k\) がどういう形の行列か順にみていくと、\( X_{k}\) はQR分解した結果、列ベクトル同士が互いに直交している行列、つまり、\( X_{k}^H X_{k} = I_m\)となります。また、 \( R_k \) はQR分解により上三角行列なので、この\( R_k \) を順番にかけていった \(\Gamma_k\)も上三角行列になります。（上三角行列同士の積は上三角行列になります。）
このため、同時反復法のアルゴリズムは\( A^kX_0 \) に対して直接QR分解を適用して \(X_{k}\) と\( \Gamma_k\) に分解したのと同じ状態になっています。このため、この後は\( A^kX_0 \) を直接QR分解したとして考えます。
では、\( A^kX_0 \) をQR分解したときに \(X_{k}\)がどういう行列になっているかを見てみます。
QR分解のアルゴリズムとしてはいくつかありますが、今回はグラム・シュミットの正規直交化法を利用したＱR分解をするとします。
グラム・シュミットの正規直交化法はWikipediaがわかりやすい漸化式になっていたので、これを参考に今回の変数表記でQR分解の漸化式を示します。
https://ja.wikipedia.org/wiki/%E3%82%B0%E3%83%A9%E3%83%A0%E3%83%BB%E3%82%B7%E3%83%A5%E3%83%9F%E3%83%83%E3%83%88%E3%81%AE%E6%AD%A3%E8%A6%8F%E7%9B%B4%E4%BA%A4%E5%8C%96%E6%B3%95
$$ \begin{align*}
U_{k,:,1} =& A^kX_{0,:,1} \tag{7.10.3} \\
X_{k,:,1} =& \frac{U_{k,:,1}}{||U_{k,:,1}||} \tag{7.10.4}\\
U_{k,:,2} =& A^kX_{0,:,2} – (X_{k,:,1}^HA^kX_{0,:,2})X_{k,:,1} \tag{7.10.5}\\
X_{k,:,2} =& \frac{U_{k,:,2}}{||U_{k,:,2}||} \tag{7.10.6}\\
…&\\
U_{k,:,m} =& A^kX_{0,:,m} – \sum_{j=1}^{m} (X_{k,:,j}^HA^kX_{0,:,m})X_{k,:,j} \tag{7.10.7}\\
X_{k,:,m} =& \frac{U_{k,:,m}}{||U_{k,:,m}||} \tag{7.10.8}\\
\end{align*} $$
\( \Gamma_k\)に関してはここからの証明で利用しないため、上の漸化式から省略しています。
\(X_{k,:,j}\)を順に見ていく前に、\( X_0 \) について見ていきます。固有ベクトルを利用すると、入力の\( X_0 \) を以下のように固有ベクトルを列ベクトルとする行列\(Z = [\boldsymbol{z}_1,…, \boldsymbol{z}_m]\)と複素数を要素に持つ \( C \) に分解できます。
$$ \begin{align*}
X_0 = ZC \tag{7.10.9}
\end{align*} $$
\(X_0\) の一つの列ベクトルに注目するとわかりやすいかと思うので、\(X_0\)の\( j \) 番目の列ベクトルを固有ベクトルで展開した場合の式を以下に示します。
$$ \begin{align*}
X_{0,:,j} =& c_{1, j} \boldsymbol{z}_1 + … + c_{n, j} \boldsymbol{z}_n \\
=& \sum_{i=1}^{m} c_{i, j} \boldsymbol{z}_i \tag{7.10.10}\\
\end{align*} $$
ここで\( c_{i,j} \) は \( C \) の\( {i, j }\) 番目の要素です。これを利用すると\( A^kX_{0,:j} \) は以下のようになります。
$$ \begin{align*}
A^k X_{0,:,j} =& \sum_{i=1}^n c_{i, j} A^k \boldsymbol{z}_i \\
=& \sum_{i=1}^n c_{i, j} \lambda_i^k \boldsymbol{z}_i
\tag{7.10.11}
\end{align*} $$
これらを利用して、ここから\( X_{k,:,1}\) から順に\( k \to \infty \) のときにどういうベクトルに収束していくのかを見ていきます。ここからはべき乗法と一部似た変形をしていきます。まずは、\( X_{k,:,1}\)のベクトルがアルゴリズムに従うとどういう形になるかを見ていきます。グラム・シュミットの正規直交化法の漸化式より以下のようになります。
$$ \begin{align*}
X_{k,:,1} =& \frac{U_{k,:,1}}{||U_{k,:,1}||} \\
=& \frac{A^kX_{0,:,1}}{||A^kX_{0,:,1}||} \tag{7.10.12}
\end{align*} $$
この式はべき乗法のときに出てきた式と同じ式になっています。このため、同じように\( |\frac{\lambda_j}{\lambda_1} |< 1\) と固有ベクトルが単位ベクトルであることを利用すると、以下のように\( X_{k,:,1} \) は\( k \to \infty \) のときに\(\boldsymbol{z}_1 \)に収束します。
$$ \begin{align*}
\lim_{k \to \infty} X_{k,:,1} =& \lim_{k \to \infty} \frac{A^kX_{0,:,1}}{||A^kX_{0,:,1}||} \\
=& \lim_{k \to \infty} \frac{\sum_{i=1}^n c_{i, 1} \lambda_i^k \boldsymbol{z}_i }{||\sum_{i=1}^n c_{i, 1} \lambda_i^k \boldsymbol{z}_i ||} \\
=& \lim_{k \to \infty} \frac{c_{1,1} \lambda_1^k \boldsymbol{z}_1 + \sum_{i=2}^n c_{i, 1} \lambda_i^k \boldsymbol{z}_i }{||c_{1, 1} \lambda_1^k \boldsymbol{z}_1 +\sum_{i=2}^n c_{i, 1} \lambda_i^k \boldsymbol{z}_i ||} \\
=& \lim_{k \to \infty} \frac{\lambda_1^k \left[ c_{1, 1} \boldsymbol{z}_1 + \sum_{i=2}^n c_{i, 1} \left( \frac{\lambda_i}{\lambda_1}\right)^k \boldsymbol{z}_i \right]}{\lambda_1^k||\left[ c_{1, 1} \boldsymbol{z}_1 +\sum_{i=2}^n c_{i, 1} \left( \frac{\lambda_i}{\lambda_1}\right)^k \boldsymbol{z}_i \right] ||} \\
=& \lim_{k \to \infty} \frac{\left[ c_{1, 1} \boldsymbol{z}_1 + \sum_{i=2}^n c_{i, 1} \left( \frac{\lambda_i}{\lambda_1}\right)^k \boldsymbol{z}_i \right]}{||\left[ c_{1, 1} \boldsymbol{z}_1 +\sum_{i=2}^n c_{i, 1} \left( \frac{\lambda_i}{\lambda_1}\right)^k \boldsymbol{z}_i \right] ||} \\
=& \lim_{k \to \infty} \frac{\left[ c_{1, 1} \boldsymbol{z}_1 + \sum_{i=2}^n c_{i, 1} \left( \frac{\lambda_i}{\lambda_1}\right)^k \boldsymbol{z}_i \right]}{||\left[ c_{1, 1} \boldsymbol{z}_1 +\sum_{i=2}^n c_{i, 1} \left( \frac{\lambda_i}{\lambda_1}\right)^k \boldsymbol{z}_i \right] ||} \\
=&\frac{c_{1, 1} \boldsymbol{z}_1}{||c_{1, 1} \boldsymbol{z}_1||} \\
=&\frac{\boldsymbol{z}_1}{||\boldsymbol{z}_1||} \\
=&\boldsymbol{z}_1
\tag{7.10.13}
\end{align*}$$
ここまではある程度本にある証明ですが、ここから先は本に説明が書いてないので、間違っている可能性がありますが、ご了承ください。
次に\( X_{k,:,2} \)が\( k \to \infty \)のときにどのような値に収束していくか見ていきます。
まず、\( X_{k,:,2} \) のベクトルがグラム・シュミットの正規直交化法に従うとどのような形なのかを詳しく見ていきます。
$$ \begin{align*}
X_{k,:,2} =& \frac{U_{k,:,2}}{||U_{k,:,2}||} \\
=& \frac{A^kX_{0,:,2} – (X_{k,:,1}^HA^kX_{0,:,2})X_{k,:,1}}{||A^kX_{0,:,2} – (X_{k,:,1}^HA^kX_{0,:,2})X_{k,:,1}||}
\tag{7.10.14}
\end{align*}$$
ここで、グラム・シュミットの正規直交化法では\(X_{k,:,1} \) と \( X_{k,:,2} \) が直交になるように変形されています。つまり、内積が0になるため以下のような関係が成り立ちます。
$$ \begin{align*}
X_{k,:,1}^HX_{k,:,2} =& \frac{U_{k,:,1}^HU_{k,:,2}}{||U_{k,:,1}||||U_{k,:,2}||} = 0
\tag{7.10.15}
\end{align*}$$
また\( k \to \infty \)のとき、\(X_{k,:,1}\) が\(\boldsymbol{z}_1\)に収束します。\(X_{k,:,1} \) と \( X_{k,:,2} \) が直交関係にあることと、\(X_{k,:,1}\) が\(\boldsymbol{z}_1\)に収束することから、\( k \to \infty \)のとき\(U_{k,:,2}\) の中にある\(\boldsymbol{z}_1\)にかかっている係数は0に収束するということがわかります。
この結果を利用しつつ、\( k \to \infty \)のときに\( U_{k,:,2}\) がどのようなベクトルに収束するか説明しやすいように以下のように変形します。
$$ \begin{align*}
U_{k,:,2} =& A^kX_{0,:,2} – (X_{k,:,1}^HA^kX_{0,:,2})X_{k,:,1} \\
=& \sum_{i=1}^n c_{i, 2} \lambda_i^k \boldsymbol{z}_i – (X_{k,:,1}^HA^kX_{0,:,2})X_{k,:,1} \\
=& c_{1, 2} \lambda_1^k \boldsymbol{z}_i1+ c_{2, 2} \lambda_2^k \boldsymbol{z}_2 + \sum_{i=3}^n c_{i, 2} \lambda_i^k \boldsymbol{z}_i \\
& – (X_{k,:,1}^HA^kX_{0,:,2})X_{k,:,1} \\
=& \lambda_2^k \left[ c_{2, 2} \boldsymbol{z}_2 + \sum_{i=3}^n c_{i, 2} \left( \frac{\lambda_i}{\lambda_2}\right)^k \boldsymbol{z}_i \right. \\
& \left. + \frac{1}{\lambda_2^k} \left( c_{1, 2} \lambda_1^k \boldsymbol{z}_1 – (X_{k,:,1}^HA^kX_{0,:,2})X_{k,:,1} \right) \right]
\tag{7.10.16}
\end{align*}$$
式変形した結果、\( k \to \infty \)のときに括弧の中の各項がどうなるかを順番に見ていきます。まず、1つ目の\( c_{2, 2} \boldsymbol{z}_2 \) は \( k \) に依存してないのでそのままです。二つ目のsumの項に関しては\( i=3 \) からになっているため、\( |\frac{\lambda_i}{\lambda_2} |< 1 \) が成り立ちます。このため、0に収束します。最後に3項目ですが、\( (X_{k,:,1}^HA^kX_{0,:,2}) \) は内積の計算に相当するため、スカラになります。また、\( k \to \infty \)のとき\( X_{k,:,1} \) は\( \boldsymbol{z}_1 \) に収束します。このため、3項目は\( k \to \infty \)のときに\( \boldsymbol{z}_1 \) に関するベクトルだけとなります。ここで、\(U_{k,:,2}\) の中にある\(\boldsymbol{z}_1\)にかかっている係数は0に収束するということを思い出すと3項目の部分が0に収束することがわかります。この結果、\( k \to \infty \)のときに1つ目の項だけが残ることがわかります。
ここで、\( k \to \infty \)のときに0に収束する部分をまとめて\( \boldsymbol{e}_{k,2} \) と置きます。そしてこの結果から、\( X_{k,:,2} \) は以下のようになります。
$$ \begin{align*}
\lim_{k \to \infty} X_{k,:,2} =& \frac{U_{k,:,2}}{||U_{k,:,2}||} \\
=& \lim_{k \to \infty} \frac{\lambda_2^k \left[ c_{2, 2} \boldsymbol{z}_2 + \boldsymbol{e}_{k,2} \right] }{||\lambda_2^k \left[ c_{2, 2} \boldsymbol{z}_2 + \boldsymbol{e}_{k,2} \right]||} \\
=& \lim_{k \to \infty} \frac{c_{2, 2} \boldsymbol{z}_2 + \boldsymbol{e}_{k,2} }{||c_{2, 2} \boldsymbol{z}_2 + \boldsymbol{e}_{k,2}||} \\
=& \frac{c_{2, 2} \boldsymbol{z}_2 }{||c_{2, 2} \boldsymbol{z}_2||} \\
=& \frac{\boldsymbol{z}_2 }{||\boldsymbol{z}_2||} \\
=& \boldsymbol{z}_2
\tag{7.10.17}
\end{align*}$$
これで\( X_{k,:,2} \) についても示すことができました。\( j \ge 2 \) のとき \( X_{k,:,j} \) は同じような証明になります。つまり、\( U_{k,:,j} \) を以下のように変形すると
$$ \begin{align*}
U_{k,:,j} =& A^kX_{0,:,j} – \sum_{j^{\prime}=1}^{m} (X_{k,:,j^{\prime}}^HA^kX_{0,:,j})X_{k,:,j^{\prime}} \\
=& \lambda_j^k \left[ c_{j, j} \boldsymbol{z}_j + \sum_{i=j+1}^n c_{i, j} \left( \frac{\lambda_i}{\lambda_j}\right)^k \boldsymbol{z}_i \right. \\
& \left. + \frac{1}{\lambda_j^k} \sum_{j^{\prime}=1}^{j-1} \left( c_{j^{\prime}, j} \lambda_{j^{\prime}}^k \boldsymbol{z}_{j^{\prime}} – (X_{k,:,j^{\prime}}^HA^kX_{0,:,j})X_{k,:,j^{\prime}} \right) \right]
\tag{7.10.18}
\end{align*}$$
括弧の中の2項目と3項目は先ほどと同じ流れによって0に収束します。このため、\( k \to \infty \)のときに0に収束する部分をまとめて\( \boldsymbol{e}_{k,j} \) と置きます。そして、この結果から、\( X_{k,:,j} \) は以下のようになります。
$$ \begin{align*}
\lim_{k \to \infty} X_{k,:,j} =& \frac{U_{k,:,j}}{||U_{k,:,j}||} \\
=& \lim_{k \to \infty} \frac{\lambda_j^k \left[ c_{j, j} \boldsymbol{z}_j + \boldsymbol{e}_{k,j} \right] }{||\lambda_j^k \left[ c_{j, j} \boldsymbol{z}_j + \boldsymbol{e}_{k,j} \right]||} \\
=& \lim_{k \to \infty} \frac{c_{j, j} \boldsymbol{z}_j + \boldsymbol{e}_{k,j} }{||c_{j, j} \boldsymbol{z}_2 + \boldsymbol{e}_{k,j}||} \\
=& \frac{c_{j, j} \boldsymbol{z}_j }{||c_{j, j} \boldsymbol{z}_j||} \\
=& \frac{\boldsymbol{z}_j}{||\boldsymbol{z}_j||} \\
=& \boldsymbol{z}_j
\tag{7.10.19}
\end{align*}$$
以上の結果から \( \lim_{k \to \infty} X_k = Z \) 、つまり、\( X_k \) の各列が絶対値が大きい固有値の固有ベクトルに収束することが示されました。
同時反復法でなぜ固有値が求められるのか？
次に同時反復法でなぜ固有値が求められるか？を説明します。同時反復法では\(X_k^H A X_k\)の行列の対角成分が絶対値の大きい\( m \) 個の固有値に収束します。本にこの証明が書かれていますが、私がぱっと分からなかったところを補足しながら説明します。ここからは本と同じように\(m=n \) として説明します。
まず、固有ベクトルを先ほどと同じように\( Z = [\boldsymbol{z}_1, … ,\boldsymbol{z}_n]\) と置き、固有値を対角成分とする行列を\( \Lambda = \text{diag}(\lambda_1,…,\lambda_n ) \) と置きます。こようにすると以下のような関係になります。
$$ \begin{align*}
A = Z \Lambda Z^{-1} \tag{7.14.1}
\end{align*} $$
また、\( Z \) をQR分解して以下のような行列に分解できるとします。
$$ \begin{align*}
Z = QR \tag{7.14.2}
\end{align*} $$
そして、固有ベクトルの収束の証明の時と同じように\(X_0 \) を以下のように複素数を要素に持つ\( C\)と\(Z\) を使って以下のように表せられるとします。
$$ \begin{align*}
X_0 = ZC \tag{7.14.3}
\end{align*} $$
また、\(C \) はLU分解によって以下のように分解できるとします。
$$ \begin{align*}
C = LU \tag{7.14.4}
\end{align*} $$
LU分解は\( L \) (下三角行列。ただし、対角成分は1)と\( U \) (上三角行列)に分ける操作です。
これらを使うと\(A^k X_0 \) は以下のように変形できます。
$$ \begin{align*}
A^kX_0 =& Z \Lambda^k Z^{-1} X_0 \\
=& Z \Lambda^k Z^{-1} ZC \\
=& Z \Lambda^k C \\
=& Z \Lambda^k LU \\
=& Z \Lambda^k L(R\Lambda^k)^{-1}(R\Lambda^k)U \\
=& Z \Lambda^k L\Lambda^{-k}R^{-1}R\Lambda^kU \\
=& QR \Lambda^k L\Lambda^{-k}R^{-1}R\Lambda^kU \\
=& Q \cdot R (\Lambda^k L\Lambda^{-k})R^{-1} \cdot R\Lambda^kU
\tag{7.14}
\end{align*} $$
ここで\( R (\Lambda^k L\Lambda^{-k})R^{-1}\) を部分を以下のようにQR分解します。
$$ \begin{align*}
R (\Lambda^k L\Lambda^{-k})R^{-1} = G_k P_k \tag{7.15}
\end{align*} $$
ここで\( G_k^H G_k = I\) であり、\( P_k \) は上三角行列です。これを式(7.14) に代入します。
$$ \begin{align*}
A^kX_0 =& Q \cdot R (\Lambda^k L\Lambda^{-k})R^{-1} \cdot R\Lambda^kU \\
=& Q \cdot G_k P_k \cdot R\Lambda^kU
\tag{7.15.1}
\end{align*} $$
上の式(7.15.1)と固有ベクトルの収束の時に使った式(7.10)を利用すると以下のような関係が成り立ちます。
$$ \begin{align*}
A^kX_0 =& Q \cdot G_k P_k \cdot R\Lambda^kU \\
=& X_k \Gamma_k
\tag{7.15.2}
\end{align*} $$
この式の\( Q G_k P_k R\Lambda^kU
= X_k \Gamma_k \) の両辺を以下のように変形します。
$$ \begin{align*}
Q G_k P_k R\Lambda^kU =& X_k \Gamma_k \\
P_k R\Lambda^kU \Gamma_k^{-1}=& (QG_k)^HX_k
\tag{7.15.3}
\end{align*} $$
ここで\(Q, G_k\)はQR分解で計算されるユニタリ行列部分であり、ユニタリ行列同士の積はユニタリ行列です。また、\( QG_k\) はユニタリ行列なので、逆行列は\( (QG_k)^H\)になっています。
この式(7.15.3)の式変形の結果の左辺と右辺のそれぞれに注目します。まず、左辺の\( P_k, R, U , \Gamma_k \)は上三角行列です。また、上三角行列の逆行列もまた上三角行列であるため\( \Gamma_k^{-1} \) も上三角行列です。また、\( \Lambda^k \) は対角行列です。そして、上三角行列と上三角行列、上三角行列と対角行列はすべて上三角行列になります。このため、左辺は上三角行列になっています。
一方、右辺を見ていきます。\( Q, G_k, X_k \) はすべてQR分解の結果のユニタリ行列になっています。ユニタリ行列とユニタリ行列の積はユニタリ行列であるため右辺はユニタリ行列です。
つまり、式(7.15.3)の行列は左辺が上三角行列であり、右辺はユニタリ行列であるため、この二つの性質を併せ持つユニタリ対角行列（すなわち対角成分以外は0で対角要素の絶対値が1の対角行列）となっています。このユニタリ対角行列を\(D_k\)と置きます。これを使うと式(7.15.3)は以下のように成り立ちます。
$$ \begin{align*}
(QG_k)^HX_k =& P_k R\Lambda^kU \Gamma_k^{-1} \\
(QG_k)^HX_k =& D_k \\
X_k = QG_kD_k
\tag{7.15.4}
\end{align*} $$
ここで\( Q, G_k \) はQR分解によるユニタリ行列であるため、この積の\(QG_k\) もユニタリ行列です。そして、\(QG_k\) の逆行列は\((QG_k)^H\)です。
次に式(7.15)に出てきた\(\Lambda^k L\Lambda^{-k}\)について考えます。\(L\)はLU分解の下三角行列であり、対角成分は1になっています。また、\( \Lambda \) は対角行列であり、対角行列の逆行列は元の対角行列の対角成分の逆数です。このため、\(\Lambda^k L\Lambda^{-k}\)の対角成分は1になっています。また、対角行列と下三角行列の積は下三角行列です。このため、\(\Lambda^k L\Lambda^{-k}\)は下三角行列です。ただ、対角行列以外 (\( i > j \))の要素は\(L\) の要素を\(l_{i,j}\)とすると以下のような式になっています。
$$ \begin{align*}
l_{i,j} \left( \frac{\lambda_i}{\lambda_j}\right)^k
\tag{7.15.5}
\end{align*} $$
ここで\( |\frac{\lambda_i}{\lambda_j}| < 1\) であるため、\( k \to \infty\) のとき0に収束します。このため、\( \Lambda^k L\Lambda^{-k} \) は\( k \to \infty\) のとき以下のようになります。
$$ \begin{align*}
\lim_{k \to \infty} \Lambda^k L\Lambda^{-k} = I_n
\tag{7.15.6}
\end{align*} $$
このため、\( k \to \infty \) のとき式(7.15)は以下のようになります。
$$ \begin{align*}
\lim_{k \to \infty} R (\Lambda^k L\Lambda^{-k})R^{-1} =& R I_n R^{-1} \\
=& R R^{-1} \\
=& I_n \\
\tag{7.15.7}
\end{align*} $$
この結果から\( k \to \infty \) のとき式(7.15)の左辺が単位行列になるため、QR分解よる行列\( G_k, P_k \) それぞれ以下のように単位行列と収束します。
$$ \begin{align*}
\lim_{k \to \infty} G_k =& I_n \\
\lim_{k \to \infty} P_k =& I_n \\
\tag{7.15.8}
\end{align*} $$
この結果を利用する前に、式(7.15.4)の両辺を式変形します。
$$ \begin{align*}
X_k =& QG_kD_k \\
X_kD_k^H =& QG_kD_kD_k^H \\
X_kD_k^H =& QG_k \\
\tag{7.15.9}
\end{align*} $$
ここで、\(D_k\) はユニタリ対角行列であるため、\(D_k\)の逆行列は\(D_k^H\)です。この式(7.15.8)は\( k \to \infty \) のとき以下のようになります。
$$ \begin{align*}
\lim_{k \to \infty} X_kD_k^H =& \lim_{k \to \infty} QG_k \\
=& Q
\tag{7.16}
\end{align*} $$
また、後ほど使う式を順番に示していきます。まず、\(A_k \) を以下のように置きます。
$$ \begin{align*}
A_k = X_k^H A X_k \tag{7.17.1}
\end{align*} $$
ここで、\(A_k\) の各要素を\( a_{i,j}^{(k)} \) とします。つまり、最終的に\(k \to \infty\)のとき\( a_{i,i}^{(k)} \)が固有値に収束していればよいことになります。
次に式(7.14.2)を利用して式(7.14.1)を変形して以下の関係を示しておきます。
$$ \begin{align*}
A = Z \Lambda Z^{-1} \\
A = QR \Lambda (QR)^{-1} \\
A = QR \Lambda R^{-1} Q^H\\
Q^H A Q= R \Lambda R^{-1}
\tag{7.17.2}
\end{align*} $$
ここで\( Q \) はユニタリ行列であるため、\( Q \)の逆行列は\( Q^H \)です。
また、\(D_k\)がユニタリ対角行列なため、以下の式が成り立ちます。
$$ \begin{align*}
A_k =& D_k A_k D_k^H
\tag{7.17.3}
\end{align*} $$
この式(7.17.2)、(7.17.3)と先ほどの式(7.16)を利用すると\(k \to \infty\)のとき\(A_k\)は以下のようになります。
$$ \begin{align*}
\lim_{k \to \infty} A_k =& \lim_{k \to \infty} D_k A_k D_k^H \\
=& \lim_{k \to \infty} D_k X_k^H A X_k D_k^H \\
=& \lim_{k \to \infty} (X_k D_k^H)^H A (X_k D_k^H) \\
=& Q^H A Q \\
=& R \Lambda R^{-1}
\tag{7.16}
\end{align*} $$
ここで、\( R \Lambda R^{-1} \) の対角成分について注目します。\( R \) は上三角行列です。この\(R\)の対角要素を\( r_{i,i}\) とします。また、上三角行列の逆行列は上三角行列であり、対角行列の逆行列の対角成分は、元の対角行列の対角成分の逆数になるので\( 1/r_{i,i} \)です。一方、\( \Lambda \) は対角行列です。このため\( R \Lambda R^{-1} \) の対角成分は以下のようになります。
$$ \begin{align*}
r_{i,i} \lambda_i \frac{1}{r_{i,i}} = \lambda_i
\tag{7.16.1}
\end{align*} $$
このため、\( k \to \infty \)のとき\( A_k\)の対角成分\(a_{i,i}^{(k)}\)は以下のようになります。
$$ \begin{align*}
\lim_{k \to \infty} a_{i,i}^{(k)} =& r_{i,i} \lambda_i \frac{1}{r_{i,i}} \\
=& \lambda_i
\tag{7.16.1}
\end{align*} $$
これで\( A_k\) の対角成分が固有値に収束することが示すことができました。
終わりに
今回の記事では同時反復法のアルゴリズムの説明を行いました。正直、本を読んでも全然わからず、しかもネットで検索しても私がわかるレベルまで詳しく書かれた証明が出てこなかったため、かなり苦労しました。同様にわからず悩んでいる方の参考に少しでもなれば幸いです。
次はQR法の軽い証明と実装を紹介する記事をかければと思っています。
