反復法と疎行列 — Jacobi・CG・GMRES・前処理子・スパース直接法

概要

反復法は連立方程式 $A \mathbf{x} = \mathbf{b}$ を近似解 $\mathbf{x}_k$ の列で逐次更新する手法で、 直接法 (LU・Cholesky) と異なり $A$ を陽に分解しない。 $A$ が大規模疎行列で、$n$ が数百万〜数億の場合は直接法のメモリが破綻するため、反復法が唯一の現実的選択肢になる。

反復法のクラスは大きく次の通り:

• 古典反復: Jacobi / Gauss-Seidel / SOR — 行ごとの更新規則、収束は遅いが前処理子・平滑化として現役
• Krylov 部分空間法: CG・MINRES・GMRES・BiCGStab・IDR・QMR — 行列ベクトル積のみで $\mathcal{K}_k(A, \mathbf{r}_0)$ 内で最適解を構成
• 多重格子 (multigrid): 階層的粗化と平滑化、PDE 離散化で漸近最適 $O(n)$

本記事では各反復法、前処理子、疎行列格納、スパース直接法、停止条件と再開戦略を解説する。 API は 反復法ソルバー、 疎行列直接法、 疎行列反復法 を参照。

古典反復: Jacobi・Gauss-Seidel・SOR

$A$ を $A = D + L + U$（対角・厳密下三角・厳密上三角）と分解し、$\mathbf{x}_{k+1}$ を $\mathbf{x}_k$ から逐次更新する。

Jacobi 法

$$\mathbf{x}_{k+1} = D^{-1} (\mathbf{b} - (L + U) \mathbf{x}_k)$$

各成分は他の成分のみに依存するため並列化が容易。 収束条件はスペクトル半径 $\rho(D^{-1}(L+U)) < 1$ で、対角優位行列で十分。 収束率は $\rho$ で支配され、離散 Poisson では $\rho = 1 - O(h^2)$（$h$ は格子幅）と非常に遅い。

Gauss-Seidel 法

$$(D + L) \mathbf{x}_{k+1} = \mathbf{b} - U \mathbf{x}_k$$

更新済みの $\mathbf{x}_{k+1}$ の成分を逐次使うため、Jacobi より速いが本質的に逐次的。 対角優位・対称正定値で収束保証。スペクトル半径は Jacobi のおよそ二乗。

SOR (Successive Over-Relaxation)

緩和パラメータ $\omega$ を導入:

$$(D + \omega L) \mathbf{x}_{k+1} = \omega \mathbf{b} + ((1 - \omega) D - \omega U) \mathbf{x}_k$$

$\omega = 1$ が Gauss-Seidel、$\omega \in (1, 2)$ で過緩和。 離散ラプラシアンに対する最適 $\omega$ は

$$\omega_{\mathrm{opt}} = \frac{2}{1 + \sqrt{1 - \rho_J^2}}, \quad \rho_J = \rho(D^{-1}(L+U))$$

で、$\omega_{\mathrm{opt}}$ を用いれば収束率は Jacobi の平方根オーダーまで改善する。 実用上は最適 $\omega$ の解析的計算が困難な問題が多く、Krylov 法に置き換えられているが、AMG の平滑化や前処理子としては現役。

対称正定値の Krylov 法: CG / PCG

共役勾配法 (Conjugate Gradient, CG) は対称正定値 $A$ に対して、 二次形式 $\varphi(\mathbf{x}) = \tfrac{1}{2} \mathbf{x}^\top A \mathbf{x} - \mathbf{b}^\top \mathbf{x}$ の最小化として $A\mathbf{x} = \mathbf{b}$ を解く。

CG の基本

各反復で残差 $\mathbf{r}_k = \mathbf{b} - A \mathbf{x}_k$ から $A$-共役な探索方向 $\mathbf{p}_k$ を構築し、$\varphi$ を最小化するステップ幅 $\alpha_k$ で進む。 重要な性質として、$k$ 回反復後の解 $\mathbf{x}_k$ は Krylov 部分空間 $\mathcal{K}_k(A, \mathbf{r}_0)$ 内で $A$-ノルム残差を最小化する最良近似であり、 厳密算術では高々 $n$ ステップで厳密解に至る (有限終了)。 実用上は丸め誤差で完全終了せず、$k \ll n$ で十分な精度に達することが目的。

収束率

条件数 $\kappa = \kappa_2(A)$ に対して

$$\frac{\|\mathbf{x}_k - \mathbf{x}^*\|_A}{\|\mathbf{x}_0 - \mathbf{x}^*\|_A} \leq 2 \left(\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}\right)^k$$

悪条件 ($\kappa \gg 1$) では収束が遅いため、ほぼ常に前処理子 $M \approx A$ を併用する PCG が標準。 前処理付き残差 $\mathbf{z}_k = M^{-1} \mathbf{r}_k$ で探索方向を構築すれば、有効条件数が $\kappa(M^{-1} A)$ に下がる。

用途

有限要素法・有限差分法による楕円型 PDE（電位、温度、構造静解析）、グラフラプラシアン、最小二乗の正規方程式、機械学習の二次最適化部分問題。 メモリは $O(n)$（3 つの $n$ ベクトルのみ）、反復毎のコストは SpMV と数本の内積で $O(\mathrm{nnz})$。

対称不定値: MINRES と SYMMLQ

$A$ が対称だが正定値でないとき、CG の根拠（凸二次形式の最小化）は崩れる。 Paige-Saunders による MINRES と SYMMLQ は、対称三項漸化 (Lanczos) を用いつつ正定値性を要求しない。

MINRES

残差の 2 ノルム $\|\mathbf{r}_k\|_2$ を Krylov 部分空間内で最小化する。 Lanczos 三項漸化で構築した三重対角行列を直接 QR 分解して $\mathbf{x}_k$ を更新するので、メモリは CG とほぼ同じ $O(n)$ で済む。 対称不定値の鞍点系（KKT 行列、Stokes 方程式）の標準解法。

SYMMLQ

誤差 $\|\mathbf{x}_k - \mathbf{x}^*\|_2$ を最小化する変種で、$A$ が特異・ほぼ特異な場合に MINRES より安定。 ただし停止判定がやや難しい。

非対称: GMRES・BiCGStab・IDR(s)・QMR

非対称 $A$ では対称三項漸化が成り立たず、より複雑な反復が必要になる。

GMRES (Generalized Minimal Residual)

Arnoldi 反復で Krylov 部分空間の正規直交基底 $Q_k$ と上 Hessenberg $H_k$ を構築し、 $\|\mathbf{b} - A \mathbf{x}_k\|_2$ を $\mathcal{K}_k$ 内で最小化する。 単調収束保証があるが、メモリは反復回数 $k$ に比例し $O(nk)$、各反復の MGS コストも $O(nk)$ で増大する。 実用的には再開 GMRES (GMRES(m))で $m = 20$–$50$ で再開する。 再開ごとに収束履歴がリセットされ、最悪収束が停滞するが、前処理が効いていれば収束は十分速い。

BiCGStab

非対称行列に対して短い漸化（CG 風の 3-term recurrence）で動作する。 Bi-Lanczos の不安定性を補正するために、各 BiCG ステップの後に GMRES(1) ステップを追加して残差を安定化する。 メモリ $O(n)$、反復あたり SpMV 2 回。安定性は GMRES に劣り、収束が振動・停止することもある。

IDR(s) (Induced Dimension Reduction)

Sonneveld-van Gijzen (2008) による比較的新しい方法。$s$ 個のシャドウベクトルを用いて Krylov 部分空間を効率的に縮約し、 BiCGStab より滑らかな収束と少ない SpMV 回数を達成する。 $s = 4$–$8$ が典型値。GMRES に対してメモリ優位、BiCGStab に対して安定性優位。 (現行 sangi には未実装、将来追加予定)

QMR (Quasi-Minimal Residual)

Bi-Lanczos に基づき、転置 $A^\top$ も必要だが疑似的に残差を最小化する。 BiCGStab が出る前は非対称の標準だった。Free et al. (TFQMR) や Freund (TFQMR) の派生で転置を回避する変種もある。 (現行 sangi には未実装、将来追加予定)

前処理子

前処理子 $M$ は $M^{-1} A \approx I$ となる近似で、$A\mathbf{x} = \mathbf{b}$ を $M^{-1} A \mathbf{x} = M^{-1} \mathbf{b}$ に変換して有効条件数を下げる。 反復ごとに $M^{-1} \mathbf{r}$ を適用するため、$M^{-1}$ は安価に計算できなければならない（$O(\mathrm{nnz})$ オーダー）。

Jacobi 前処理

$M = D = \mathrm{diag}(A)$。最も単純だが、対角優位な問題でしか効果が薄い。並列性は最大。

SSOR (Symmetric SOR)

SOR の対称版で、対称正定値系に使える前処理。$M^{-1} = (D/\omega + U)^{-1} (D/\omega) (D/\omega + L)^{-1}$ の形。 分解はゼロ fill-in で前処理セットアップが安価。

不完全 LU / 不完全 Cholesky (ILU / IC)

$A$ の非零パターンに合わせた近似分解 $\tilde{L}\tilde{U} \approx A$。fill-in レベル $p$ と閾値 $\tau$ で疎性と精度のトレードオフを制御する。 symmetric positive definite には IC、それ以外は ILU。 セットアップは $O(\mathrm{nnz})$、適用は前進・後退代入で $O(\mathrm{nnz})$。 Manteuffel シフト $\alpha I$ で安定化することも多い。LU・Cholesky の記事も参照。

AMG (代数的多重格子)

幾何構造に依存せず、行列の代数構造のみから階層的粗化と平滑化を構築する。 典型的に Gauss-Seidel・Jacobi を平滑化、強連結成分の粗化で粗格子を構成し、粗格子で直接解または再帰 AMG を呼ぶ。 楕円型 PDE で漸近最適 $O(n)$ の前処理を提供し、CG / GMRES と組み合わせると $n$ が $10^8$ オーダーでも実用速度が出る。 セットアップが重く、コードも複雑なため、HYPRE や Trilinos など外部ライブラリの呼び出しが一般的。 (現行 sangi には未実装、将来追加予定)

疎行列格納形式

非零要素のみを格納することで、$n$ が大きく非零密度 $\mathrm{nnz}/n^2$ が小さい行列のメモリと計算コストを劇的に削減する。

CSR (Compressed Sparse Row)

• values[k]: 非零値を行順に並べた長さ $\mathrm{nnz}$ の配列
• col_idx[k]: 各非零値の列番号、長さ $\mathrm{nnz}$
• row_ptr[i]: $i$ 行目の最初の非零値の values 内インデックス、長さ $n+1$

SpMV $\mathbf{y} = A \mathbf{x}$ は行ごとに values[row_ptr[i] : row_ptr[i+1]] と x[col_idx[...]] の内積で計算でき、行優先メモリアクセスが効率的。 sangi の既定形式。

CSC (Compressed Sparse Column)

CSR の列版。SpMV $A^\top \mathbf{x}$ や、LU 分解で列単位のピボット処理に向く。 同じ行列を CSR と CSC で二重に保持することもある。

COO (Coordinate)

$(i, j, v)$ のタプルリスト。読み込み・行列構築には便利だが、SpMV に不向きで、構築後に CSR/CSC に変換するのが一般的。

BSR (Block Sparse Row)

$b \times b$ ブロック単位で非零ブロックを CSR 形式で並べる。 構造化メッシュの PDE 離散化のように非零が小ブロックを成す場合、BLAS-3 を活用してスループットを上げられる。

SpMV のメモリバウンド性

SpMV は典型的にメモリ帯域がボトルネック。$\mathbf{x}$ のランダムアクセス（col_idx 経由）のキャッシュミスを減らす再順序付け (Cuthill-McKee, RCM, AMD) や、Roofline モデルでの最適化がしばしば検討される。

スパース直接法

反復法では収束が極端に遅いか保証できない問題（強い不定値、悪条件、複数右辺）では、直接法を疎行列向けに改良したスパース直接法が有力になる。

記号フェーズと数値フェーズ

スパース直接法は通常 2 段階で進む:

1. 記号分解 (symbolic factorization): $L, U$ の非零パターンを決定する。fill-in 最小化のための行・列順序付け (AMD, COLAMD, METIS) と elimination tree の構築。
2. 数値分解 (numerical factorization): 決定されたパターンに数値を流し込む。

同じ非零構造で右辺だけ変わる問題（時間刻みの異なる ODE、ニュートン法の反復）では記号フェーズを再利用でき、数値フェーズが繰り返しコストになる。

simplicial と supernodal/multifrontal

simplicial 法は要素単位で更新（典型的に CHOLMOD の simplicial モード）、 supernodal/multifrontal はパターンが類似した連続行をスーパーノードとしてまとめ、BLAS-3 で更新する。 前者は実装が単純、後者は速いが複雑。CHOLMOD・SuperLU・MUMPS・PARDISO が代表的なライブラリで、それぞれ supernodal/multifrontal を採用。

順序付け (ordering)

fill-in は順序付けに大きく依存する。 AMD (Approximate Minimum Degree)、Nested Dissection（METIS）、COLAMD（非対称用）が標準。 最適順序付け問題は NP-hard で、ヒューリスティクスが使われる。

sangi のスパース直接法

sangi は 疎行列直接法 API で simplicial Cholesky と疎 LU を提供する。 中規模 ($n \lesssim 10^5$) の疎行列を直接解くのに使い、それを超える規模では反復法（CG/GMRES + 前処理子）が現実的選択肢。

停止条件・stagnation・再開戦略

停止条件

典型的な停止条件は次の通り:

• 絶対残差: $\|\mathbf{r}_k\|_2 < \tau_a$
• 相対残差: $\|\mathbf{r}_k\|_2 / \|\mathbf{b}\|_2 < \tau_r$ （標準）
• 後退誤差ベース: $\|\mathbf{r}_k\|_2 / (\|A\|_2 \|\mathbf{x}_k\|_2 + \|\mathbf{b}\|_2) < \tau$
• 最大反復数: $k < k_{\max}$

後退誤差ベースは Higham (2002) が推奨する条件で、悪条件問題で過剰反復を防ぐ。 相対残差のみでは前進誤差の保証にならず、$\kappa(A) \cdot \|\mathbf{r}_k\| / \|\mathbf{b}\|$ が前進相対誤差の上界。

Stagnation 検知

$\|\mathbf{r}_k\|_2$ が数反復にわたり減少しない場合、停滞 (stagnation) と判定して以下のいずれかを行う:

• 前処理子の更新: より強い前処理子に切り替える（IC(0) → IC(1) など）
• 再開: GMRES なら再開、CG なら現在解を初期値として再起動
• 反復精化 (iterative refinement): 直接法フォールバック + 残差再計算
• 諦める: 反復法が問題に合っていない可能性、直接法に切り替え

再開戦略

GMRES(m) の再開はメモリ削減のためだが、収束履歴がリセットされるため最適な $m$ の選択が重要。 $m$ が小さすぎると stagnation、大きすぎるとメモリと反復コストが増す。 実用上 $m = 20$–$50$ から始めて、収束が遅ければ $m$ を倍々で増やす。

適応的再開（収束履歴を見て $m$ を動的に変える）や、固有値情報を引き継ぐ deflated GMRES などの高度な変種もある。

sangi の選択

参考文献

• Saad, Y. (2003). Iterative Methods for Sparse Linear Systems, 2nd ed., SIAM.
• Hestenes, M. R. and Stiefel, E. (1952). "Methods of conjugate gradients for solving linear systems". J. Res. Nat. Bur. Standards, 49, 409–436.
• Saad, Y. and Schultz, M. H. (1986). "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems". SIAM J. Sci. Stat. Comput., 7(3), 856–869.
• van der Vorst, H. A. (1992). "BiCGSTAB: A fast and smoothly converging variant of Bi-CG". SIAM J. Sci. Stat. Comput., 13(2), 631–644.
• Sonneveld, P. and van Gijzen, M. B. (2008). "IDR(s): A family of simple and fast algorithms for solving large nonsymmetric systems of linear equations". SIAM J. Sci. Comput., 31(2), 1035–1062.
• Davis, T. A. (2006). Direct Methods for Sparse Linear Systems, SIAM.
• Briggs, W. L., Henson, V. E., and McCormick, S. F. (2000). A Multigrid Tutorial, 2nd ed., SIAM.