Eigenvalue algorithm
数学 において 、 べき乗反復法（ べき乗法 とも呼ばれる ）は 固有値アルゴリズム です。 対角化可能な 行列 が与えられると、このアルゴリズムは、（絶対値で）最大の 固有値 である 数と、 （対応する 固有ベクトル である ）非零ベクトル （つまり）を生成します 。このアルゴリズムは、 フォン ・ミーゼス 反復法 としても知られています。 [1] A {\displaystyle A} λ {\displaystyle \lambda } A {\displaystyle A} v {\displaystyle v} λ {\displaystyle \lambda } A v = λ v {\displaystyle Av=\lambda v}
べき乗反復法は非常に単純なアルゴリズムですが、収束が遅い場合があります。このアルゴリズムで最も時間のかかる演算は、行列とベクトルの乗算であるため、適切に実装すれば、非常に大きな 疎行列 に効果的です 。収束速度は のようになります。 ここで は反復回数です（後のセクションを参照）。言い換えれば、収束は スペクトルギャップ を底とする指数関数的です。 A {\displaystyle A} ( λ 2 / λ 1 ) k {\displaystyle (\lambda _{2}/\lambda _{1})^{k}} k {\displaystyle k}
法 2x2行列上のべき乗反復法アルゴリズムを視覚化するアニメーション。行列は2つの固有ベクトルで表されます。誤差は次のように計算されます | | approximation − largest eigenvector | | {\displaystyle ||{\text{approximation}}-{\text{largest eigenvector}}||} べき乗反復アルゴリズムはベクトル から開始します。ベクトル は、支配的な固有ベクトルの近似値、またはランダムベクトルである可能性があります。この方法は、 再帰関係によって記述されます。 b 0 {\displaystyle b_{0}}
b k + 1 = A b k ‖ A b k ‖ {\displaystyle b_{k+1}={\frac {Ab_{k}}{\|Ab_{k}\|}}} したがって、反復ごとに、ベクトルは 行列 で乗算され 、正規化されます。 b k {\displaystyle b_{k}} A {\displaystyle A}
が 他の固有値よりも絶対値が大きい固有値を持ち、開始ベクトルが 支配的な固有値に関連付けられた固有ベクトルの方向に非ゼロの成分を持つと仮定すると、部分列は 支配的な固有値に関連付けられた固有ベクトルに収束します。 A {\displaystyle A} b 0 {\displaystyle b_{0}} ( b k ) {\displaystyle \left(b_{k}\right)}
上記の2つの仮定がない場合、列は 必ずしも収束しません。この列では、 ( b k ) {\displaystyle \left(b_{k}\right)}
b k = e i ϕ k v 1 + r k {\displaystyle b_{k}=e^{i\phi _{k}}v_{1}+r_{k}} 、 ここで、 は支配的な固有値に関連付けられた固有ベクトルであり、 です 。項の存在は、 でない限り収束しない ことを意味します。上記の2つの仮定の下では、 で定義される 列は v 1 {\displaystyle v_{1}} ‖ r k ‖ → 0 {\displaystyle \|r_{k}\|\rightarrow 0} e i ϕ k {\displaystyle e^{i\phi _{k}}} ( b k ) {\displaystyle \left(b_{k}\right)} e i ϕ k = 1 {\displaystyle e^{i\phi _{k}}=1} ( μ k ) {\displaystyle \left(\mu _{k}\right)}
μ k = b k ∗ A b k b k ∗ b k {\displaystyle \mu _{k}={\frac {b_{k}^{*}Ab_{k}}{b_{k}^{*}b_{k}}}} 支配的な固有値（ レイリー商 ）に収束します。 [ 説明が必要 ]
これは以下のアルゴリズム（PythonとNumPyで示されています）で計算できます。
#!/usr/bin/env python3 import numpy as np def power_iteration ( A : np . ndarray , num_iterations : int ) -> np . ndarray : # 理想的にはランダムベクトルを選択します # ベクトルが # 固有ベクトルに直交する可能性を減らすため b_k = np . random . rand ( A . shape [ 1 ] ) for _ in range ( num_iterations ): # 行列とベクトルの積Abを計算します b_k1 = np . dot ( A , b_k ) # ノルムを計算します b_k1_norm = np.linalg.norm ( b_k1 ) # ベクトル b_k = b_k1 / b_k1_norm b_k 返す power_iteration ( np . array ([[ 0.5 , 0.5 ], [ 0.2 , 0.8 ]]), 10 ) ベクトルは関連する固有ベクトルに収束します。理想的には、 関連する固有値を得るために レイリー商を 使用する必要があります b k {\displaystyle b_{k}}
このアルゴリズムは、 Google PageRank を 計算 するために使用されます
この方法は、レイリー商を計算することでスペクトル半径 （正方行列の場合、最大の大きさを持つ固有値） を計算するためにも使用できます。
ρ ( A ) = max { | λ 1 | , … , | λ n | } = b k ⊤ A b k b k ⊤ b k . {\displaystyle \rho (A)=\max \left\{|\lambda _{1}|,\dotsc ,|\lambda _{n}|\right\}={\frac {b_{k}^{\top }Ab_{k}}{b_{k}^{\top }b_{k}}}.}
分析 をジョルダン標準形 に分解するとします 。 ここ で、 の最初の列は の固有ベクトルで、これは 支配的固有値 に対応します 。 一般に の支配的固有値は 一意なので、 の最初のジョルダンブロックは の 最大の固有値の大きさである 行列 です 。 開始ベクトルは の列の線形結合として表すことができます 。 A {\displaystyle A} A = V J V − 1 {\displaystyle A=VJV^{-1}} V {\displaystyle V} A {\displaystyle A} λ 1 {\displaystyle \lambda _{1}} A {\displaystyle A} J {\displaystyle J} 1 × 1 {\displaystyle 1\times 1} [ λ 1 ] , {\displaystyle [\lambda _{1}],} λ 1 {\displaystyle \lambda _{1}} A {\displaystyle A} b 0 {\displaystyle b_{0}} V {\displaystyle V}
b 0 = c 1 v 1 + c 2 v 2 + ⋯ + c n v n . {\displaystyle b_{0}=c_{1}v_{1}+c_{2}v_{2}+\cdots +c_{n}v_{n}.} 仮定により、 は支配的固有ベクトルの方向に非ゼロの成分を持つため、 となります 。 b 0 {\displaystyle b_{0}} c 1 ≠ 0 {\displaystyle c_{1}\neq 0}
の計算上有用な 漸化式は 次 のように書き直すことができます。 b k + 1 {\displaystyle b_{k+1}}
b k + 1 = A b k ‖ A b k ‖ = A k + 1 b 0 ‖ A k + 1 b 0 ‖ , {\displaystyle b_{k+1}={\frac {Ab_{k}}{\|Ab_{k}\|}}={\frac {A^{k+1}b_{0}}{\|A^{k+1}b_{0}\|}},} ここで、式 は 次の解析により適しています。 A k + 1 b 0 ‖ A k + 1 b 0 ‖ {\displaystyle {\frac {A^{k+1}b_{0}}{\|A^{k+1}b_{0}\|}}}
b k = A k b 0 ‖ A k b 0 ‖ = ( V J V − 1 ) k b 0 ‖ ( V J V − 1 ) k b 0 ‖ = V J k V − 1 b 0 ‖ V J k V − 1 b 0 ‖ = V J k V − 1 ( c 1 v 1 + c 2 v 2 + ⋯ + c n v n ) ‖ V J k V − 1 ( c 1 v 1 + c 2 v 2 + ⋯ + c n v n ) ‖ = V J k ( c 1 e 1 + c 2 e 2 + ⋯ + c n e n ) ‖ V J k ( c 1 e 1 + c 2 e 2 + ⋯ + c n e n ) ‖ = ( λ 1 | λ 1 | ) k c 1 | c 1 | v 1 + 1 c 1 V ( 1 λ 1 J ) k ( c 2 e 2 + ⋯ + c n e n ) ‖ v 1 + 1 c 1 V ( 1 λ 1 J ) k ( c 2 e 2 + ⋯ + c n e n ) ‖ {\displaystyle {\begin{aligned}b_{k}&={\frac {A^{k}b_{0}}{\|A^{k}b_{0}\|}}\\&={\frac {\left(VJV^{-1}\right)^{k}b_{0}}{\|\left(VJV^{-1}\right)^{k}b_{0}\|}}\\&={\frac {VJ^{k}V^{-1}b_{0}}{\|VJ^{k}V^{-1}b_{0}\|}}\\&={\frac {VJ^{k}V^{-1}\left(c_{1}v_{1}+c_{2}v_{2}+\cdots +c_{n}v_{n}\right)}{\|VJ^{k}V^{-1}\left(c_{1}v_{1}+c_{2}v_{2}+\cdots +c_{n}v_{n}\right)\|}}\\&={\frac {VJ^{k}\left(c_{1}e_{1}+c_{2}e_{2}+\cdots +c_{n}e_{n}\right)}{\|VJ^{k}\left(c_{1}e_{1}+c_{2}e_{2}+\cdots +c_{n}e_{n}\right)\|}}\\&=\left({\frac {\lambda _{1}}{|\lambda _{1}|}}\right)^{k}{\frac {c_{1}}{|c_{1}|}}{\frac {v_{1}+{\frac {1}{c_{1}}}V\left({\frac {1}{\lambda _{1}}}J\right)^{k}\left(c_{2}e_{2}+\cdots +c_{n}e_{n}\right)}{\left\|v_{1}+{\frac {1}{c_{1}}}V\left({\frac {1}{\lambda _{1}}}J\right)^{k}\left(c_{2}e_{2}+\cdots +c_{n}e_{n}\right)\right\|}}\end{aligned}}} 上記の式は次のように簡略化されます。 k → ∞ {\displaystyle k\to \infty }
( 1 λ 1 J ) k = [ [ 1 ] ( 1 λ 1 J 2 ) k ⋱ ( 1 λ 1 J m ) k ] → [ 1 0 ⋱ 0 ] as k → ∞ . {\displaystyle \left({\frac {1}{\lambda _{1}}}J\right)^{k}={\begin{bmatrix}[1]&&&&\\&\left({\frac {1}{\lambda _{1}}}J_{2}\right)^{k}&&&\\&&\ddots &\\&&&\left({\frac {1}{\lambda _{1}}}J_{m}\right)^{k}\\\end{bmatrix}}\rightarrow {\begin{bmatrix}1&&&&\\&0&&&\\&&\ddots &\\&&&0\\\end{bmatrix}}\quad {\text{as}}\quad k\to \infty .} この極限は、 の固有値の大きさが1未満である という事実から導かれます。したがって、 1 λ 1 J i {\displaystyle {\frac {1}{\lambda _{1}}}J_{i}}
( 1 λ 1 J i ) k → 0 as k → ∞ . {\displaystyle \left({\frac {1}{\lambda _{1}}}J_{i}\right)^{k}\to 0\quad {\text{as}}\quad k\to \infty .} したがって、
1 c 1 V ( 1 λ 1 J ) k ( c 2 e 2 + ⋯ + c n e n ) → 0 as k → ∞ {\displaystyle {\frac {1}{c_{1}}}V\left({\frac {1}{\lambda _{1}}}J\right)^{k}\left(c_{2}e_{2}+\cdots +c_{n}e_{n}\right)\to 0\quad {\text{as}}\quad k\to \infty } この事実を用いて、 が大きい 場合 との関係を強調する形で を書くことができます。 b k {\displaystyle b_{k}} v 1 {\displaystyle v_{1}} k {\displaystyle k}
b k = ( λ 1 | λ 1 | ) k c 1 | c 1 | v 1 + 1 c 1 V ( 1 λ 1 J ) k ( c 2 e 2 + ⋯ + c n e n ) ‖ v 1 + 1 c 1 V ( 1 λ 1 J ) k ( c 2 e 2 + ⋯ + c n e n ) ‖ = e i ϕ k c 1 | c 1 | v 1 ‖ v 1 ‖ + r k {\displaystyle {\begin{aligned}b_{k}&=\left({\frac {\lambda _{1}}{|\lambda _{1}|}}\right)^{k}{\frac {c_{1}}{|c_{1}|}}{\frac {v_{1}+{\frac {1}{c_{1}}}V\left({\frac {1}{\lambda _{1}}}J\right)^{k}\left(c_{2}e_{2}+\cdots +c_{n}e_{n}\right)}{\left\|v_{1}+{\frac {1}{c_{1}}}V\left({\frac {1}{\lambda _{1}}}J\right)^{k}\left(c_{2}e_{2}+\cdots +c_{n}e_{n}\right)\right\|}}\\[6pt]&=e^{i\phi _{k}}{\frac {c_{1}}{|c_{1}|}}{\frac {v_{1}}{\|v_{1}\|}}+r_{k}\end{aligned}}} ここで 、 および として e i ϕ k = ( λ 1 / | λ 1 | ) k {\displaystyle e^{i\phi _{k}}=\left(\lambda _{1}/|\lambda _{1}|\right)^{k}} ‖ r k ‖ → 0 {\displaystyle \|r_{k}\|\to 0} k → ∞ {\displaystyle k\to \infty }
数列は 有界であるため、収束する部分数列を含みます。支配的な固有値に対応する固有ベクトルはスカラーまでしか一意ではないため、数列は 収束しない可能性がありますが、 が大きいに対して はほぼ の固有ベクトルになります 。 ( b k ) {\displaystyle \left(b_{k}\right)} ( b k ) {\displaystyle \left(b_{k}\right)} b k {\displaystyle b_{k}} A {\displaystyle A} k {\displaystyle k}
あるいは、が 対角化 可能であれば 、次の証明で同じ結果が得られます。 A {\displaystyle A}
を の固有値（重複度付き） と し 、 を対応する固有ベクトルとします。 が 支配的な固有値であると仮定すると、 に対してとなります 。 λ 1 , λ 2 , … , λ m {\displaystyle \lambda _{1},\lambda _{2},\ldots ,\lambda _{m}} m {\displaystyle m} A {\displaystyle A} v 1 , v 2 , … , v m {\displaystyle v_{1},v_{2},\ldots ,v_{m}} λ 1 {\displaystyle \lambda _{1}} | λ 1 | > | λ j | {\displaystyle |\lambda _{1}|>|\lambda _{j}|} j > 1 {\displaystyle j>1}
初期ベクトルは次のよう に書けます。 b 0 {\displaystyle b_{0}}
b 0 = c 1 v 1 + c 2 v 2 + ⋯ + c m v m . {\displaystyle b_{0}=c_{1}v_{1}+c_{2}v_{2}+\cdots +c_{m}v_{m}.} をランダムに（一様確率で）選択する と、 確率1 で収束します 。ここで、 b 0 {\displaystyle b_{0}} c 1 ≠ 0 {\displaystyle c_{1}\neq 0}
A k b 0 = c 1 A k v 1 + c 2 A k v 2 + ⋯ + c m A k v m = c 1 λ 1 k v 1 + c 2 λ 2 k v 2 + ⋯ + c m λ m k v m = c 1 λ 1 k ( v 1 + c 2 c 1 ( λ 2 λ 1 ) k v 2 + ⋯ + c m c 1 ( λ m λ 1 ) k v m ) → c 1 λ 1 k v 1 | λ j λ 1 | < 1 for j > 1 {\displaystyle {\begin{aligned}A^{k}b_{0}&=c_{1}A^{k}v_{1}+c_{2}A^{k}v_{2}+\cdots +c_{m}A^{k}v_{m}\\&=c_{1}\lambda _{1}^{k}v_{1}+c_{2}\lambda _{2}^{k}v_{2}+\cdots +c_{m}\lambda _{m}^{k}v_{m}\\&=c_{1}\lambda _{1}^{k}\left(v_{1}+{\frac {c_{2}}{c_{1}}}\left({\frac {\lambda _{2}}{\lambda _{1}}}\right)^{k}v_{2}+\cdots +{\frac {c_{m}}{c_{1}}}\left({\frac {\lambda _{m}}{\lambda _{1}}}\right)^{k}v_{m}\right)\\&\to c_{1}\lambda _{1}^{k}v_{1}&&\left|{\frac {\lambda _{j}}{\lambda _{1}}}\right|<1{\text{ for }}j>1\end{aligned}}} 一方、
b k = A k b 0 ‖ A k b 0 ‖ . {\displaystyle b_{k}={\frac {A^{k}b_{0}}{\|A^{k}b_{0}\|}}.} したがって、 は固有ベクトル （の倍数）に収束します 。収束は 幾何収束 で、比は b k {\displaystyle b_{k}} v 1 {\displaystyle v_{1}}
| λ 2 λ 1 | , {\displaystyle \left|{\frac {\lambda _{2}}{\lambda _{1}}}\right|,} です。ここで、 は 2番目の支配的な固有値を表します。したがって、支配的な固有値と大きさが近い固有値がある場合、この方法はゆっくりと収束します。 λ 2 {\displaystyle \lambda _{2}}
応用 べき乗反復法は行列の 1 つの固有値のみを近似しますが、特定の 計算問題 には依然として有用です。たとえば、 Google は 検索エンジンでドキュメントの PageRank を 計算するために使用しています [2]。 また、 Twitter は フォローすべき人物の推薦をユーザーに表示するために使用しています べき乗反復法は 、Web 行列などの 疎行列 に特に適しています。また、係数行列を明示的に格納する必要がなく 、代わりに行列ベクトル積を評価する関数にアクセスできる 行列フリー法としても使用できます。条件 が整えられた 非対称行列の場合、べき乗反復法はより複雑な Arnoldi 反復法 よりも優れたパフォーマンスを発揮します 。対称行列の場合、反復あたりのコストが小さいことを犠牲にすることなく収束速度を簡単に上げることができるため、べき乗反復法はほとんど使用されません。たとえば、 Lanczos 反復法 や LOBPCG を 参照してください。 A {\displaystyle A} A x {\displaystyle Ax}
より高度な固有値アルゴリズムのいくつかは、べき乗反復法のバリエーションとして理解できます。例えば、 逆反復 法は行列にべき乗反復法を適用します 。他のアルゴリズムは、ベクトルによって生成された部分空間全体を調べます 。この部分空間は クリロフ部分空間として知られています。これは アーノルディ反復法 または ランチョス反復法 によって計算できます 。グラム反復法 [4] は、最大固有値対を計算するための超線形かつ決定論的な方法です。 A − 1 {\displaystyle A^{-1}} b k {\displaystyle b_{k}}
参照
参考文献 ^ Richard von Mises and H. Pollaczek-Geiringer, Praktische Verfahren der Gleichungsauflösung , ZAMM - Zeitschrift für Angewandte Mathematik und Mechanik 9, 152-164 (1929). ^ Ipsen, Ilse 、およびRebecca M. Wills (2005年5月5～8日). 「第7回IMACS国際科学計算反復法シンポジウム」 (PDF) . フィールズ研究所、カナダ、トロント. {{cite news }}: CS1 maint: multiple names: authors list (link )^ Delattre, B.; Barthélemy, Q.; Araujo, A.; Allauzen, A. (2023)、「グラム反復法による畳み込み層のリプシッツ定数の効率的な境界」、 第40回国際機械学習会議論文集 ： 7513– 7532
![{\displaystyle [\lambda_{1}],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0f3c83a8bae265cd5523c330d3b5f9901fee6ab0)
![{\displaystyle \left({\frac {1}{\lambda _{1}}}J\right)^{k}={\begin{bmatrix}[1]&&&&\\&\left({\frac {1}{\lambda _{1}}}J_{2}\right)^{k}&&&\\&&\ddots &\\&&&\left({\frac {1}{\lambda _{1}}}J_{m}\right)^{k}\\\end{bmatrix}}\rightarrow {\begin{bmatrix}1&&&&\\&0&&\\&&\ddots &\\&&&0\\\end{bmatrix}}\quad {\text{as}}\quad k\to \infty .}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cec7d61aee7e07ab39497f6cf782f074c415a7fd)
![{\displaystyle {\begin{aligned}b_{k}&=\left({\frac {\lambda _{1}}{|\lambda _{1}|}}\right)^{k}{\frac {c_{1}}{|c_{1}|}}{\frac {v_{1}+{\frac {1}{c_{1}}}V\left({\frac {1}{\lambda _{1}}}J\right)^{k}\left(c_{2}e_{2}+\cdots +c_{n}e_{n}\right)}{\left\|v_{1}+{\frac {1}{c_{1}}}V\left({\frac {1}{\lambda _{1}}}J\right)^{k}\left(c_{2}e_{2}+\cdots +c_{n}e_{n}\right)\right\|}}\\[6pt]&=e^{i\phi _{k}}{\frac {c_{1}}{|c_{1}|}}{\frac {v_{1}}{\|v_{1}\|}}+r_{k}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/024832f7d5c7731b36d56df90fde0723dd480484)
