LOBPCG

局所最適ブロック前処理共役勾配法（LOBPCG ）は、対称一般化固有値問題の最大（または最小）固有値とそれに対応する固有ベクトルを見つけるための行列を使用しない方法です。

${\displaystyle Ax=\lambda Bx,}$

与えられた複素エルミート行列または実対称行列のペアに対して、行列は正定値であると仮定されます。 ${\displaystyle (A,B)}$${\displaystyle B}$

1948 年にKantorovich は、スカラー積におけるレイリー商のスケール勾配の方向を用いた最急降下法によって対称行列の最小の固有値 を計算することを提案しました。このとき、ステップ サイズは、ベクトル と の線形スパンでレイリー商を最小化することによって計算されます(つまり、局所的に最適な方法)。Samokish ^{[1] は}、残差ベクトルに前処理子を適用して前処理された方向を生成し、固有ベクトルに近づくにつれて漸近収束率の境界を導出することを提案しました。D'yakonov ^[2]は、スペクトル的に等価な前処理を提案し、非漸近収束率の境界を導出しました。固有値問題に対するブロック局所最適多段階最急降下法は、^[3]で説明されました。現在の近似値、現在の残差、および前回の近似値によって張られる部分空間上のレイリー商の局所最小化、およびそのブロックバージョンは、^[4]で紹介されました。前処理バージョンは、 ^[5]および^{[6]で分析されました。}${\displaystyle \lambda_{1}}$${\displaystyle A}$${\displaystyle r=Ax-\lambda (x)x}$ ${\displaystyle \lambda (x)=(x,Ax)/(x,x)}$ ${\displaystyle (x,y)=x'y}$${\displaystyle x}$${\displaystyle w}$ ${\displaystyle T}$${\displaystyle r}$${\displaystyle w=Tr}$${\displaystyle x}$

出典: ^[7]

• 行列フリー、つまり係数行列を明示的に保存する必要はなく、行列ベクトル積を評価することによって行列にアクセスできます。
• 因数分解不要、つまり一般化固有値問題であっても行列分解を必要としません。
• 反復あたりのコストとメモリ使用量は、対称行列の単一の極値固有値対を計算するLanczos 法と同等です。
• 線形収束は理論的に保証され、実際に観察されます。
• ランチョス法とは対照的に、可変および非対称の前処理と固定および正定値の前処理を含む直接の前処理により、収束が加速されます。
• 事前処理により、効率的な領域分割とマルチグリッド技術を簡単に組み込むことができます。
• ウォームスタートを実行し、各反復ごとに固有ベクトルの近似値を計算します。
• ランチョス法に比べて数値的に安定しており、低精度のコンピュータ演算でも動作可能です。
• 実装が簡単で、すでに多くのバージョンが登場しています。
• ブロッキングにより、 BLAS 3などの非常に効率的な行列間演算を利用できるようになります。
• ブロック サイズは、収束速度と、各反復における直交化およびレイリー リッツ法のコンピューター コストとのバランスをとるように調整できます。

この方法は、一般化レイリー商の反復最大化（または最小化）を実行する。

${\displaystyle \rho (x):=\rho (A,B;x):={\frac {x^{T}Ax}{x^{T}Bx}},}$

その結果、最大（または最小）の固有値対が見つかる。${\displaystyle Ax=\lambda Bx.}$

一般化レイリー商の最も急な上昇の方向、つまり勾配は、ベクトル

${\displaystyle r:=Ax-\rho (x)Bx,}$

固有ベクトル残差と呼ばれる。前処理 が利用可能な場合は、残差に適用され、ベクトル ${\displaystyle T}$

${\displaystyle w:=Tr,}$

を前処理残差と呼ぶ。前処理をしない場合は、 と設定する。反復法 では、${\displaystyle T:=I}$${\displaystyle w:=r}$

${\displaystyle x^{i+1}:=x^{i}+\alpha^{i}T(Ax^{i}-\rho (x^{i})Bx^{i}),}$

つまり、

${\displaystyle x^{i+1}:=x^{i}+\alpha^{i}w^{i},\,}$

${\displaystyle w^{i}:=Tr^{i},\,}$

${\displaystyle r^{i}:=Ax^{i}-\rho (x^{i})Bx^{i},}$

これは前処理付き最急上昇（または最急降下）として知られており、スカラーは ステップサイズと呼ばれます。最適なステップサイズはレイリー商を最大化することで決定できます。つまり、 ${\displaystyle \alpha^{i}}$

${\displaystyle x^{i+1}:=\arg \max _{y\in span\{x^{i},w^{i}\}}\rho (y)}$

（または最小化する場合）、その場合、その方法は局所最適と呼ばれます。 ${\displaystyle \arg \min}$

局所的に最適な前処理付き最急上昇（または最急降下）の収束を劇的に加速するには、2 項の再帰関係に 1 つのベクトルを追加して3 項にします。

${\displaystyle x^{i+1}:=\arg \max _{y\in span\{x^{i},w^{i},x^{i-1}\}}\rho (y)}$

（最小化する場合に使用する）。3次元部分空間におけるレイリー商の最大化/最小化は、レイリー・リッツ法によって数値的に実行できます。リチャードソン外挿など、ベクトルを追加しても大幅な加速は得られませんが^[8]、計算コストが増加するため、一般的には推奨されません。 ${\displaystyle \arg \min}$

反復が収束するにつれて、ベクトルと は にほぼ線形従属するようになり、精度の低下を招き、丸め誤差がある場合にはレイリー・リッツ法は数値的に不安定になります。精度の低下は、3次元部分空間 の基底で、ベクトル をからさらに離れたベクトル に置き換え、部分空間を変更せずに直交化やその他の追加操作を回避することで回避できます。^{[8]さらに、3次元部分空間の基底を直交化することは、}条件の悪い固有値問題の安定性と達成可能な精度を向上させるために必要になる場合があります。 ${\displaystyle x^{i}}$${\displaystyle x^{i-1}}$${\displaystyle x^{i-1}}$${\displaystyle p^{i}}$${\displaystyle x^{i}}$${\displaystyle span\{x^{i},w^{i},x^{i-1}\}}$

これはLOBPCG法の単一ベクトル版であ​​り、対称固有値問題の場合に前処理付き 共役勾配線形ソルバーを一般化できる方法の1つです。 ^[8]とが自明な場合でも、両方の近似値は同じクリロフ部分空間に属しますが、結果の近似値はランチョスアルゴリズムで得られる近似値とは異なります。 ${\displaystyle T=I}$${\displaystyle B=I}$${\displaystyle i>3}$

LOBPCG の単一ベクトル バージョンは極めてシンプルで効率が高いため、組み込みASICまたはFPGAでのグラフ分割によるスペクトルクラスタリングベースのリアルタイム異常検出から、エクサスケールTOP500スーパーコンピューターでの記録的な計算複雑度の物理現象のモデリングまで、厳しいハードウェア制限下での固有値関連のアプリケーションに最適です。

後続の固有ペアは、直交デフレーションを補足した単一ベクトル LOBPCG を介して 1 つずつ計算することも、ブロックとして同時に計算することもできます。前者のアプローチでは、すでに計算された近似固有ベクトルの不正確さが、その後に計算される固有ベクトルの精度に加算的に影響し、新しい計算ごとにエラーが増加します。LOBPCG のブロック バージョンでは、複数の近似固有ベクトルをブロック内で局所的に最適な方法で反復します。^{[8] により、}単一ベクトル LOBPCG では収束が遅い、ほぼ複数の固有値に対応するものも含め、高速で正確で堅牢な固有ベクトルの計算が可能になります。ブロック サイズを調整して、数値安定性、収束速度、直交化とレイリー リッツ法の計算コストのバランスをとることができます。

LOBPCG のブロック アプローチでは、単一ベクトルとをブロック ベクトル、つまり行列とに置き換えます。ここで、たとえば、 のすべての列は、固有ベクトルの 1 つを近似します。すべての列は同時に反復され、近似固有ベクトルの次の行列は、行列およびのすべての列によって張られる部分空間上でレイリー–リッツ法によって決定されます。 の各列は、 のすべての列について前処理された残差として単純に計算されます。行列は、の列と の列によって張られる部分空間が同じになるように決定されます。 ${\displaystyle x^{i},\,w^{i},}$${\displaystyle p^{i}}$${\displaystyle X^{i},\,W^{i},}$${\displaystyle P^{i}}$${\displaystyle X^{i}}$${\displaystyle X^{i+1}}$${\displaystyle X^{i},\,W^{i},}$${\displaystyle P^{i}}$${\displaystyle W^{i}}$${\displaystyle X^{i}.}$${\displaystyle P^{i}}$${\displaystyle [X^{i},\,P^{i}]}$${\displaystyle [X^{i},\,X^{i-1}]}$

レイリー・リッツ法の結果は、行列およびのすべての列によって張られる部分空間によって決定されます。ここで、部分空間の基底は理論的には任意です。ただし、不正確なコンピュータ演算では、基底ベクトルの一部が近似的に線形従属している場合、レイリー・リッツ法は数値的に不安定になります。数値的不安定性は通常、たとえば反復ブロック内の固有ベクトルの一部が、特定のコンピュータ精度で達成可能な精度にすでに達している場合に発生し、特に単精度などの低精度で顕著になります。 ${\displaystyle X^{i},\,W^{i},}$${\displaystyle P^{i}}$

LOBPCG の複数の異なる実装の技術は、 部分空間の適切な基底を選択することで、最小限の計算コストでレイリー・リッツ法の数値安定性を確保することです。グラム・シュミット過程などによって基底ベクトルを直交させる、おそらく最も安定したアプローチは、最も計算コストが高いものでもあります。 たとえば、LOBPCG の実装^{[9] [10]は、}正規行列の不安定だが効率的なコレスキー分解を利用しています。これは、部分空間全体ではなく、個々の行列 と に対してのみ実行されます。 コンピュータのメモリ量は絶えず増加しているため、今日では一般的なブロック サイズを の範囲に収めることができます。この範囲では、直交化とレイリー・リッツ法に費やされる計算時間の割合が支配的になり始めます。 ${\displaystyle W^{i}}$${\displaystyle P^{i}}$${\displaystyle 10^{3}-10^{4}}$

サブスペースを反復する固有値問題のブロック法では、一般的に、一部の反復固有ベクトルが他のものよりも速く収束するため、すでに収束した固有ベクトルをロックする、つまり反復ループから削除して、不要な計算を排除し、数値安定性を向上させる必要があります。 固有ベクトルを単純に削除すると、反復中のベクトルにその重複が形成される可能性があります。 対称固有値問題の固有ベクトルがペアワイズ直交であるという事実は、すべての反復ベクトルをロックされたベクトルに対して直交に保つことを示唆しています。

ロックは、計算コストを最小限に抑えながら数値精度と安定性を維持するさまざまな方法で実装できます。たとえば、LOBPCG の実装^{[9] [10]は、 [8] [11]}に従って、ロックされた固有ベクトルがコード入力として機能し、変更されないハード ロック、つまり制限によるデフレーションを、ロックされたベクトルが残差を計算するという通常最もコストのかかる反復ステップには参加しないが、レイリー リッツ法には完全に参加し、レイリー リッツ法によって変更できるソフト ロックから分離します。

LOBPCG は、行列のすべての列をレイリー・リッツ法に含め、その結果、最大で- 行 - 列の固有値問題を解き、各反復で最大で- 行 - 列のドット積を計算する必要があります。ここで、 はブロック サイズ、つまり列の数を表します。ブロック サイズが大きい場合、これにより計算コストと I/O コストが支配的になり、複数の計算デバイスが同時に実行される並列化が制限されます。 ${\displaystyle X^{i},\,W^{i},}$${\displaystyle P^{i}}$${\displaystyle 3k}$${\displaystyle 3k}$${\displaystyle 9k^{2}}$ ${\displaystyle k}$${\displaystyle k}$

LOBPCG のオリジナルの論文^[8]では、LOBPCG II と呼ばれる修正法について説明しており、これは、各所望の固有値対に対して LOBPCG 法の単一ベクトル版を実行し、3 x 3 の投影固有値問題をレイリー・リッツ法で解くという問題に対処するものです。すべての固有値対に対するグローバルなレイリー・リッツ法は、各反復で実行されますが、行列 の列に対してのみ実行されるため、反復ごとにから への必要なドット積の数が から に削減され、グローバル投影固有値問題のサイズが-by-から-by-に削減されます。参考文献^[12]では、さらに進んで、LOBPCG アルゴリズムを各近似固有ベクトルに個別に適用します。つまり、固定回数の反復で、各所望の固有値対に対して LOBPCG 法のブロック化されていないバージョンを実行します。これらの実行におけるレイリー・リッツ法では、3 x 3 の投影固有値問題のセットを解くだけで済みます。すべての必要な固有ペアに対するグローバル レイリー リッツ手順は、ブロックされていない LOBPCG 反復の固定数の最後にのみ定期的に適用されます。 ${\displaystyle k}$${\displaystyle k}$${\displaystyle X^{i}}$${\displaystyle k^{2}+3k}$${\displaystyle 9k^{2}}$${\displaystyle k}$${\displaystyle k}$${\displaystyle 3k}$${\displaystyle 3k}$

このような変更は、元の LOBPCG に比べて堅牢性が低くなる可能性があります。単一ベクトル LOBPCG の個別に実行されるブランチは、連続した反復パスに従わず、代わりに反転して、同じ固有ベクトルへの重複した近似値を作成する可能性があります。単一ベクトル LOBPCG はクラスター化された固有値には適していない可能性がありますが、個別の小ブロック LOBPCG 実行では、固有値のクラスターの数とサイズが事前に不明な場合があるため、反復プロセス中にブロック サイズを自動的に決定する必要があります。

LOBPCG は、その構造上、包括的な収束理論を持つブロック最急勾配降下法よりも遅くない速度でレイリー商を最小化することが保証されています^{[8] 。すべての}固有ベクトルは、勾配が消失するレイリー商の停留点です。したがって、勾配降下法は、任意の固有ベクトルの近傍で速度が低下する可能性がありますが、線形収束率で固有ベクトルに収束するか、この固有ベクトルが鞍点である場合は、反復レイリー商が対応する固有値を下回り、次の固有値に線形収束し始める可能性が高くなることが保証されています。線形収束率の最悪値は決定されており^{[8] 、固有値と行列}スペクトルの残りの部分との相対的なギャップと、前処理子 （存在する場合）の品質に依存します。

一般行列の場合、固有ベクトルを予測して常に適切に機能する初期近似値を生成する方法は明らかにありません。LOBPCG による反復解法は、初期固有ベクトル近似値の影響を受けやすい可能性があります。たとえば、中間の固有ペアを通過するにつれて収束に時間がかかり、速度が低下します。さらに、理論上は、ミスの確率がゼロであっても、必ずしも最小の固有ペアへの収束を保証することはできません。LOBPCGでは、平均がゼロの良質なランダム ガウス関数が初期近似値を生成するためのデフォルトとして一般的に使用されています。初期近似値を固定するには、乱数ジェネレータに固定シードを選択できます。

Lanczos 法とは対照的に、LOBPCG は実際には 漸近的な超線形収束を示すことはほとんどありません。

LOBPCG は、例えば、共分散行列D ^T D を明示的に計算することなく、つまり行列を使わない方法で、平均がゼロのデータ行列Dに対してPCA の反復計算を行うなど、いくつかの最大特異値とそれに対応する特異ベクトル (部分 SVD) を計算するために簡単に適応できます。主な計算は、共分散行列D ^T DとブロックベクトルXの積D ^T (DX)の関数の評価で、これは目的の特異ベクトルを反復的に近似します。PCA は共分散行列の最大固有値を必要としますが、LOBPCG は通常、最小固有値を計算するために実装されます。簡単な回避策は、関数を否定し、D ^T (DX ) を-D ^T ( DX)に置き換えて固有値の順序を逆にすることです。これは、LOBPCG が固有値問題の行列が正定値であるかどうかを気にしないためです。^[9]

PCAとSVDのためのLOBPCGはSciPyのリビジョン1.4.0以降に実装されている^[13]

LOBPCGの発明者であるAndrew Knyazevは、 PETSc、hypre 、およびParallel Hierarchical Adaptive MultiLevel method (PHAML)へのインターフェースを備えたBlock Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) ^{[14] [15]}と呼ばれる参照実装を公開しました。^[16]その他の実装としては、GNU Octave、^[17] MATLAB（分散配列やタイル配列を含む）、^[9] Java、^[18] Anasazi（Trilinos）、^[19] SLEPc、^{[20] [21]} SciPy 、^[10] Julia、^[22] MAGMA、^[23] Pytorch、^[24] Rust、^[25] OpenMPおよびOpenACC、^[26] CuPy（CUDAによって高速化されたNumPy互換の配列ライブラリ）、^[27] Google JAX、^[28] NVIDIA AMGX などがあります。 ^{[29] LOBPCGは}TensorFlowに実装されていますが、^[30]含まれていません。

ソフトウェアパッケージscikit-learnとMegaman ^[31]はLOBPCGを使用して、スペクトルクラスタリング^[32]とラプラシアン固有マップを介したマニフォールド学習^[33]を大規模なデータセットに拡張します。NVIDIAはCUDA 8で導入されたnvGRAPHライブラリにLOBPCGを実装しました^{[34 ]。Sphynx [35]}は、分散メモリと共有メモリのハイブリッドに対応した並列グラフパーティショナーで、分散メモリ設定のGPUで動作する最初のグラフパーティショニングツールであり、グラフパーティショニングにスペクトルクラスタリングを使用し、AnasaziパッケージのLOBPCGを使用してグラフのラプラシアン行列上の固有ベクトルを計算します。

LOBPCGは、ABINIT ^[36]（CUDA版を含む）とOctopus ^[37]に実装されています。これは、日本の地球シミュレータスーパーコンピュータ上で、ゴードン・ベル賞の最終候補者によって数十億サイズの行列に使用されています。 ^{[38] [39]}超伝導の背後にあるメカニズムを理解するための強相関電子系のハバードモデルでは、LOBPCGを使用して、京コンピュータ^[40]とマルチGPUシステム上でハミルトニアンの基底状態を計算しています。 ^[41]

LOBPCGには、平面波基底を用いたコーン・シャム方程式と密度汎関数理論（DFT） 用のMATLAB ^[42]とJulia ^{[43] [44]}バージョンがあります。最近の実装には、TTPY、 ^[45] Platypus‐QM、^[46] MFDn、^[47] ACE-Molecule、^{[48] LACONIC [49]}などがあります。

BLOPEXのLOBPCGは、制約によるマルチレベルバランス領域分割（BDDC）ソルバーライブラリBDDCMLの前処理設定に使用され、OpenFTL（オープン有限要素テンプレートライブラリ）および破砕多孔質媒体内の地下水の流れ、溶質および熱輸送のFlow123dシミュレータに組み込まれています。LOBPCGはLS-DYNAに^[50]実装されており、間接的にANSYSに実装されています。^[51]

LOBPCG は、PYFEMax および高性能マルチフィジックス有限要素ソフトウェア Netgen/NGSolve のコア固有値ソルバーの 1 つです。hypre の LOBPCG は、有限要素法MFEM用のオープンソースの軽量でスケーラブルなC++ライブラリに組み込まれており、 BLAST、XBraid、VisIt 、xSDK、 SciDACの FASTMath 研究所、およびエクサスケール コンピューティングプロジェクトの共同設計 Center for Efficient Exascale Discretizations (CEED)など、多くのプロジェクトで使用されています。

反復LOBPCGベースの近似ローパスフィルタはノイズ除去に使用できます。 ^[52]を参照。例えば、全変動ノイズ除去を加速します。

スペクトルクラスタリングによる画像セグメンテーションは、ピクセル間のアフィニティ行列を使用して低次元の埋め込みを実行し、続いて低次元空間で固有ベクトルの成分をクラスタリングします（たとえば、バイラテラルフィルタにグラフラプラシアンを使用）。マルチグリッド前処理によるLOBPCGによるスペクトルグラフ分割による画像セグメンテーションは、[ 53 ]で最初に提案され、[ ⁵⁴ ]^{と[ 55]}で実際にテストされました。後者のアプローチは、後にPython scikit-learn ^[56]で実装され、グラフラプラシアンの固有値問題を解くために代数マルチグリッド前処理を備えたSciPyのLOBPCGを使用しています。

• MATLABでの LOBPCG
• Octaveの LOBPCG
• SciPyでの LOBPCG
• Google CodeのJavaでの LOBPCG
• GitHubのブロック局所最適前処理固有値ソルバー (BLOPEX) の LOBPCGとGoogle Codeにアーカイブされています