QRアルゴリズム

数値線形代数において、QRアルゴリズムまたはQR反復法は固有値アルゴリズム、つまり行列の固有値と固有ベクトルを計算する手順です。QRアルゴリズムは、1950年代後半にJohn GF FrancisとVera N. Kublanovskayaによって独立して開発されました。^{[1] [2] [3]} 基本的な考え方は、 QR分解を実行し、行列を直交行列と上三角行列の積として表し、逆の順序で因子を乗算して反復することです。

正式には、A を固有値を計算したい実数行列とし、A ₀  := Aとします。k番目のステップ ( k = 0から開始) では、 QR 分解 A _k = Q _k  R _kを計算します。ここで、Q _kは直交行列(つまり、Q ^T = Q ⁻¹ )、R _kは上三角行列です。次に、A _{k +1} = R _k  Q _kを形成します。 すべてのA _{k は}類似しており、したがって同じ固有値を持つことに注意してください。このアルゴリズムは、直交相似変換によって進行するため、数値的に安定しています。 $${\displaystyle A_{k+1}=R_{k}Q_{k}=Q_{k}^{-1}Q_{k}R_{k}Q_{k}=Q_{k}^{-1} A_{k}Q_{k}=Q_{k}^{\mathsf {T}}A_{k}Q_{k},}$$

ある条件の下では、^[4]行列A _kは三角行列、つまりAのシュア形式に収束する。三角行列の固有値は対角線上に並べられ、固有値問題が解かれる。収束のテストでは正確なゼロを要求することは非現実的であるが^{[引用が必要]}、ゲルシュゴリンの円定理は誤差の上限を与える。

上記の粗い形式では、反復処理のコストが比較的高くなります。この問題は、まず行列Aを上ヘッセンベルグ形式(ハウスホルダー削減に基づく手法を使用した算術演算が必要) にすることで軽減できます。これには、直交相似変換の有限シーケンスが使用され、両側 QR 分解に似ています。^{[5] [6]} (QR 分解では、ハウスホルダー反射子は左側でのみ乗算されますが、ヘッセンベルグの場合は、左右両方で乗算されます。) 上ヘッセンベルグ行列の QR 分解を決定するには、算術演算が必要です。さらに、ヘッセンベルグ形式は既にほぼ上三角であるため (各対角線の下に 0 以外の要素が 1 つだけある)、これを開始点として使用すると、QR アルゴリズムの収束に必要なステップ数が削減されます。 ${\textstyle {\tfrac {10}{3}}n^{3}+{\mathcal {O}}(n^{2})}$${\textstyle 6n^{2}+{\mathcal {O}}(n)}$

元の行列が対称であれば、上ヘッセンベルク行列も対称で三角行列となり、すべてのA _kも対称となる。この場合、ヘッセンベルク形式に到達するには、ハウスホルダー還元に基づく手法を使用した算術演算が必要となる。^{[5] [6]}対称三角行列の QR 分解を決定するには、演算が必要となる。^[7]${\textstyle {\tfrac {4}{3}}n^{3}+{\mathcal {O}}(n^{2})}$${\displaystyle {\mathcal {O}}(n)}$

ヘッセンベルク行列が何らかの に対して要素を持つ場合、つまり、対角線のすぐ下の要素の 1 つが実際にゼロである場合、その行列はブロックに分解され、そのブロックの固有値問題は個別に解くことができます。固有値は、最初の行と列の部分行列の固有値か、残りの行と列の部分行列の固有値のいずれかです。QR 反復ステップの目的は、これらの要素の 1 つを縮小して、実質的に対角線に沿った小さなブロックが行列の大部分から分離されるようにすることです。実固有値の場合、それは通常、右下隅のブロックです (この場合、要素はその固有値を保持します)。一方、共役複素固有値のペアの場合は、それは右下隅のブロックです。 ${\displaystyle A}$${\displaystyle a_{k,k-1}=0}$${\displaystyle k}$${\displaystyle k-1}$${\displaystyle a_{k,k-1}}$${\displaystyle 1\times 1}$${\displaystyle a_{nn}}$${\displaystyle 2\times 2}$

収束率は固有値間の分離に依存するため、実用的なアルゴリズムでは、明示的または暗黙的なシフトを使用して分離を増やし、収束を加速します。一般的な対称 QR アルゴリズムは、1 回または 2 回の反復で各固有値を分離し (その後、行列のサイズを縮小します)、効率的かつ堅牢になります。^{[説明が必要]}

基本的な QR アルゴリズムは、 A が正定値対称行列である場合に視覚化できます。その場合、 A は2 次元の楕円または高次元の楕円体として表すことができます。アルゴリズムへの入力と 1 回の反復の関係は、図 1 のように表すことができます (クリックするとアニメーションが表示されます)。LR アルゴリズムが QR アルゴリズムと並べて表されていることに注意してください。

1 回の反復で、楕円は x 軸に向かって傾くか「倒れ」ます。楕円の長い半軸がx 軸に平行な場合、QR の 1 回の反復では何も行われません。アルゴリズムが「何もしない」別の状況は、長い半軸が x 軸ではなく y 軸に平行な場合です。その場合、楕円はどちらの方向にも倒れることなく不安定にバランスをとっていると考えられます。どちらの状況でも、行列は対角です。アルゴリズムの反復が「何もしない」状況を固定点と呼びます。アルゴリズムが採用する戦略は、固定点に向かう反復です。1 つの固定点は安定していますが、もう 1 つは不安定であることに注意してください。楕円が不安定な固定点からわずかに傾いている場合、QR の 1 回の反復で、楕円は固定点に向かうのではなく、固定点から離れる方向に傾きます。最終的には、アルゴリズムは別の固定点に収束しますが、それには長い時間がかかります。

対称行列の固有ベクトルを1つでも見つけることは計算不可能である（計算可能解析の定義に従った正確な実数演算では）ことを指摘しておく価値がある。 ^[8]この困難は、行列の固有値の重複度が分からない場合に常に存在する。一方、固有値を見つける場合には同じ問題は存在しない。行列の固有値は常に計算可能である。

ここで、これらの困難が基本的な QR アルゴリズムでどのように現れるかについて説明します。これは、図 2 に示されています。楕円は正定値対称行列を表していることを思い出してください。入力行列の 2 つの固有値が互いに近づくと、入力楕円は円に変わります。円は単位行列の倍数に対応します。近似円は、固有値が行列の対角要素にほぼ等しい単位行列の近似倍数に対応します。したがって、その場合、固有値を近似的に見つける問題は簡単であることが示されています。ただし、楕円の半軸に何が起こるかに注意してください。QR (または LR) の反復により、入力楕円が円に近づくにつれて、半軸の傾きがますます小さくなります。半軸が x 軸と y 軸に平行である場合にのみ、固有ベクトルを知ることができます。入力楕円がより円形になるほど、ほぼ平行状態を達成するために必要な反復回数は無制限に増加します。

任意の対称行列の固有分解を計算することは不可能かもしれませんが、行列を任意の小さな量だけ摂動させて、結果の行列の固有分解を計算することは常に可能です。行列がほぼ円として描かれている場合、その行列を完全な円として描かれた行列に置き換えることができます。その場合、行列は単位行列の倍数であり、その固有分解は即時に行われます。ただし、結果の固有基底は元の固有基底からかなり離れる可能性があることに注意してください。

楕円が円形に近づくにつれて速度が低下することには逆のことが言えます。楕円が引き伸ばされて円形でなくなると、楕円の回転が速くなります。楕円が表す行列を に置き換えると、このような引き伸ばしが誘発されます。ここで、はの最小の固有値にほぼ相当します。この場合、楕円の 2 つの半軸の比は に近づきます。高次元では、このようにシフトすると、楕円体の最小の半軸の長さが他の半軸に比べて小さくなり、最小の固有値への収束が速くなりますが、他の固有値への収束は速まりません。最小の固有値が完全に決定されると、これは役に立たなくなるため、行列を に縮小する必要があります。これは、最後の行と列を削除することを意味します。 ${\displaystyle M}$${\displaystyle M-\lambda I}$${\displaystyle \lambda}$${\displaystyle M}$${\displaystyle \infty}$

不安定な固定点の問題にも対処する必要があります。シフト ヒューリスティックは、この問題にも対処するために設計されることがよくあります。実際のシフトは不連続でランダムであることが多いです。ウィルキンソンのシフトは、私たちが視覚化しているような対称行列に適していますが、特に不連続です。

現代の計算実践では、QR アルゴリズムは、多重シフトの使用を導入しやすくする暗黙のバージョンで実行されます。^[4]行列は、最初に明示バージョンと同様に上ヘッセンベルグ形式になります。次に、各ステップで、 の最初の列が小さなサイズのハウスホルダー相似変換によって の最初の列に変換されます^{[明確化が必要]} (または)。ここで、次数 の は、シフト戦略を定義する多項式です (多くの場合、ここでと は、の末尾の主部分行列の 2 つの固有値、いわゆる暗黙の二重シフトです)。次に、サイズのハウスホルダー変換を連続して実行して、作業行列を上ヘッセンベルグ形式に戻します。この操作は、アルゴリズムのステップに沿った行列の非ゼロ要素の特異な形状のため、バルジ追跡として知られています。最初のバージョンと同様に、 の下対角要素の 1 つが十分に小さくなるとすぐにデフレーションが実行されます。 ${\displaystyle A_{0}=QAQ^{\mathsf {T}}}$${\displaystyle A_{k}}$${\displaystyle p(A_{k})}$ ${\displaystyle p(A_{k})e_{1}}$${\displaystyle p(A_{k})}$${\displaystyle r}$${\displaystyle p(x)=(x-\lambda )(x-{\bar {\lambda }})}$${\displaystyle \lambda}$${\displaystyle {\bar {\lambda }}}$${\displaystyle 2\times 2}$${\displaystyle A_{k}}$${\displaystyle r+1}$${\displaystyle A_{k}}$${\displaystyle A_{k}}$

この手順の現代の暗黙的バージョンではQR分解が明示的に実行されないため、Watkins ^[9]などの一部の著者は、その名前をFrancisアルゴリズムに変更することを提案しました。GolubとVan Loanは、Francis QRステップという用語を使用しています。

QR アルゴリズムは、基本的な「べき乗」固有値アルゴリズムのより洗練されたバリエーションと見ることができます。べき乗アルゴリズムは、Aに単一のベクトルを繰り返し掛け、各反復の後に正規化することを思い出してください。ベクトルは、最大の固有値の固有ベクトルに収束します。代わりに、QR アルゴリズムは、QR 分解を使用して再正規化 (および直交化) し、ベクトルの完全な基底で動作します。対称行列Aの場合、収束すると、AQ = QΛになります。ここで、ΛはAが収束した固有値の対角行列であり、Qはそこに到達するのに必要なすべての直交相似変換の合成です。したがって、 Qの列は固有ベクトルです。

QR アルゴリズムの前には、QR 分解の代わりにLU 分解を使用するLR アルゴリズムがありました。QR アルゴリズムの方が安定しているため、LR アルゴリズムは現在ではほとんど使用されていません。ただし、これは QR アルゴリズムの開発における重要なステップを表しています。

LRアルゴリズムは、当時ETHチューリッヒのエドゥアルト・シュティーフェルの研究助手として働いていたハインツ・ルティスハウザーによって1950年代初頭に開発されました。シュティーフェルは、ルティスハウザーに、モーメントのシーケンスy ₀ ^T A ^k x ₀、k = 0、1、...（ここでx ₀とy _{0は任意のベクトル）を使用して}Aの固有値を見つけることを提案しました。ルティスハウザーは、このタスクのためにアレクサンダー・エイトケンのアルゴリズムを採用し、それを商差アルゴリズムまたはqdアルゴリズムに発展させました。計算を適切な形に整理した後、彼はqdアルゴリズムが実際には三角行列に適用された反復A _k = L _k U _k（LU分解）、A _{k +1} = U _k L _kであることを発見しました。そこからLRアルゴリズムが続きます。^[10]

QR アルゴリズムの 1 つの変種であるGolub -Kahan-Reinschアルゴリズムは、一般行列を二重対角行列に縮小することから始まります。^[11]特異値の計算のためのこのQR アルゴリズムの変種は、Golub と Kahan (1965) によって最初に説明されました。LAPACKサブルーチンDBDSQR は、特異値が非常に小さい場合をカバーするためにいくつかの変更を加えたこの反復法を実装しています (Demmel と Kahan 1990)。ハウスホルダー反射と、適切な場合はQR 分解を使用する最初のステップと合わせて、これは特異値分解を計算するための DGESVD ルーチンを形成します。QR アルゴリズムは、対応する収束結果を持つ無限次元でも実装できます。^{[12] [13]}

• デメル、ジェームズ、カハン、ウィリアム(1990)。「二重対角行列の正確な特異値」。SIAM Journal on Scientific and Statistical Computing。11 ( 5): 873–912。CiteSeerX 10.1.1.48.3740。doi : 10.1137/ 0911052。 
• Golub, Gene H. ; Kahan, William (1965). 「行列の特異値と擬似逆行列の計算」. Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis . 2 (2): 205–224. Bibcode :1965SJNA....2..205G. doi :10.1137/0702016. JSTOR  2949777.

• PlanetMathにおける固有値問題。
• ピーター・J・オルバーによる直交基底と QR アルゴリズムの動作に関するメモ
• QRメソッド用モジュール
• C++ ライブラリ