反復法による連立一次方程式の求解(CG・GMRES)

なぜ直接法では足りないのか

有限要素法や有限差分法で偏微分方程式を離散化すると、未知数が数万から数十億に及ぶ連立一次方程式 Ax = b が現れます。行列 A は疎（ほとんどの成分がゼロ）ですが、これを直接法（ガウス消去法・LU分解）で解こうとすると フィルイン（fill-in） が問題になります。消去の過程で元々ゼロだった成分が非ゼロに変わり、疎性が失われていく現象です。

次数付け（reordering、nested dissection等）を工夫すればフィルインはある程度抑えられますが、次元が上がるほど効果は薄れます。目安として、nested dissectionを使った疎コレスキー分解の計算量は2次元格子で未知数数のおよそ1.5乗、3次元格子ではおよそ2乗に達し、メモリ使用量も次元が上がるほど急増します。3次元問題では数千万自由度を超えると単一計算では手に負えなくなります。

これに対し反復法は、行列 A を明示的に分解せず、行列ベクトル積 Av の計算だけを使って解に近づいていきます。疎行列なら Av は非ゼロ要素数に比例するコストで済み、行列を分解形式で保持する必要もありません。大規模疎行列系で反復法が主役になる理由はここにあります。

共役勾配法(CG) ── 対称正定値行列に特化した最適性

CG法は A が対称かつ正定値（symmetric positive definite、SPD）であることを前提に、Ax = b を二次形式 f(x) = (1/2)xᵀAx − bᵀx の最小化問題に読み替えて解きます。この二次形式の勾配がゼロになる点が方程式の解と一致するため、最適化として解けるわけです。

単純な最急降下法は勾配方向にしか進まず、細長い谷（条件数の大きい行列）でジグザグに振動して収束が遅くなります。CG法はこれを回避するため、過去の探索方向すべてに対して A 内積の意味で直交（共役）な新しい方向を毎回選びます。これにより、理論上は丸め誤差がなければ未知数の次元数以内（n回）で厳密解に到達することが保証されます。

CG法の1反復で行う主な計算(概念)
  1. 残差 r = b − Ax を計算
  2. 直前までの探索方向と共役になるよう新方向 p を更新
  3. ステップ幅 α を、残差を最小化するよう解析的に決定
  4. x ← x + α·p, r ← r − α·A p で更新
  → 必要なのは Av 型の行列ベクトル積とベクトル内積のみ

実務上重要なのは、有限回収束の保証は理想条件下の話であり、実際の収束速度を左右するのは 条件数 κ(A)（最大固有値と最小固有値の比） だということです。CG法の誤差減少率は概ね次のように評価されます。

CG法の収束レート(SPD行列、誤差はAノルムで評価)
  反復k回後の誤差 ≦ 2 * ((√κ − 1)/(√κ + 1))^k * 初期誤差

  → κ が大きい(悪条件)ほど収束が遅い
  → κ = 1(単位行列相当)なら1回で収束

条件数が大きい（悪条件）行列では、理論上の有限回収束を待つ余裕はなく、実用上は数十〜数百反復で許容誤差まで落とすことが目標になります。

GMRES ── 非対称行列のための一般化

多くの実問題（移流拡散方程式、非対称な離散化スキーム、非対称な物理系）では A が対称ではありません。CG法の共役方向という考え方はそのままでは使えないため、GMRES（Generalized Minimal Residual） が使われます。

GMRESは、初期残差から生成される クリロフ部分空間（Krylov subspace） {b, Ab, A²b, …, A^(k-1)b} の中で、残差ノルム ||b − Ax|| を最小にする近似解を毎反復探します。内部ではアーノルディ法（Arnoldi process）で正規直交基底を作りながら空間を拡張していきます。

GMRESの反復イメージ
  k反復目までに構築したクリロフ空間 K_k = span{b, Ab, ..., A^(k-1)b}
  x_k = argmin_{x ∈ x0 + K_k} ||b − Ax||₂
  → 反復が進むほど探索空間が広がり残差は単調非増加

CG法との決定的な違いは、GMRESは過去の基底ベクトルすべてを保持し、直交化する必要がある点です。反復 k 回目には k 本のベクトルを記憶し、直交化コストは k に比例して増えていきます。これはメモリと計算コストが反復回数とともに増大することを意味し、実用上は一定反復数（例えば30〜50回）でクリロフ空間を打ち切って再起動する リスタート版GMRES(m) が使われます。ただしリスタートは過去の探索方向の情報を捨てるため、収束が停滞する場合があります。

前処理(preconditioning)が実務上ほぼ必須である理由

CG法もGMRESも、収束速度は行列 A の固有値分布（実質的には条件数）に支配されます。離散化の格子を細かくするほど、多くの偏微分方程式由来の行列は条件数が悪化します（典型的にはラプラシアン型で格子幅の2乗に反比例して条件数が増大）。素の A に反復法を適用すると、格子を細かくするたびに反復回数が増え続け、実用にならなくなります。

前処理は、Ax = b を 同じ解を持つがはるかに条件の良い系 M⁻¹Ax = M⁻¹b に変換する手法です。M は A を近似しつつ、M⁻¹v の適用が安価に計算できる行列として選びます。

前処理の基本的な考え方
  元の系:      Ax = b            (κ(A) が大きく収束が遅い)
  前処理後:    M⁻¹A x = M⁻¹b     (κ(M⁻¹A) ≈ 1 に近づける)

  M の条件:
    ・A を良く近似する(M⁻¹A が単位行列に近い)
    ・M⁻¹v の計算が Av 同様に安価

代表的な前処理には、不完全LU分解（ILU、フィルインを制限したLU分解）、不完全コレスキー分解（IC、SPD行列向け）、ヤコビ前処理（対角成分のみを使う最も単純な形）、そして多重格子法（multigrid、異なる解像度の格子を往復して低周波誤差と高周波誤差を別々に消す手法）があります。多重格子法は理想的には反復回数が問題サイズに依存しない（格子を細かくしても反復回数がほぼ一定の）性質を持ち、大規模系での前処理として特に重視されます。

まとめ

• 大規模疎行列系では直接法のフィルインがメモリ・計算量を破綻させるため、行列ベクトル積のみで進む反復法が必須になる。
• CG法は対称正定値行列専用で、過去の探索方向と共役な方向を選ぶことで理論上n回以内の収束を保証し、メモリも反復あたり一定。
• GMRESは非対称行列にも使えるが、クリロフ部分空間の基底を保持・直交化するためメモリとコストが反復数に比例して増え、実務ではリスタート版が使われる。
• 両手法とも収束速度は行列の条件数に支配され、格子を細かくするほど条件数が悪化するため、ILUや多重格子法などの前処理なしでは大規模計算に実用耐性がない。