共役勾配法シミュレーター 戻る
数値線形代数シミュレーター

共役勾配法シミュレーター — CG 法と最急降下法の収束比較

対称正定値行列 A に対する Ax=b を共役勾配法で反復解。残差の対数履歴を最急降下法と比較し、なぜ CG が n 反復で厳密収束するのかを学べます。

パラメータ設定
反復回数 k_max
A[0][0] (対角)
A[1][1] (対角)
b[0] (右辺)

A は対称三対角 4×4 行列。固定要素: A[1][2]=A[2][1]=1, A[2][3]=A[3][2]=1, A[0][1]=A[1][0]=1, A[2][2]=2, A[3][3]=1。b の残り: b[1]=2, b[2]=3, b[3]=4。初期推定 x_0=0。

A =
b =
計算結果
解 x_1
残差 ‖r‖_2
厳密収束反復 k*
条件数 cond(A)
残差ノルムの収束履歴 log10(‖r_k‖_2)

青実線=CG 法/赤実線=最急降下法/横軸=反復 k/縦軸=log10(‖r_k‖_2)

理論・主要公式

CG 法は、残差 $r_k = b - A x_k$ と探索方向 $p_k$ を A-直交(共役)に保ちながら更新する反復法です。

初期化: $r_0 = b - A x_0$, $p_0 = r_0$。反復ステップ:

$$\alpha_k = \frac{r_k^\top r_k}{p_k^\top A p_k}, \quad x_{k+1} = x_k + \alpha_k p_k$$ $$r_{k+1} = r_k - \alpha_k A p_k, \quad \beta_k = \frac{r_{k+1}^\top r_{k+1}}{r_k^\top r_k}$$ $$p_{k+1} = r_{k+1} + \beta_k p_k$$

理想算術では n 次元の対称正定値系は高々 n 反復で厳密解に到達します。誤差の収束率は条件数 $\kappa = \lambda_\max/\lambda_\min$ により $\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^k$ で抑えられます。

共役勾配法シミュレーターとは

🙋
FEM の参考書を読んでたら、剛性方程式 Kd=f を CG 法で解くって何度も出てくるんですけど、CG 法って何ですか?
🎓
ざっくり言うと、対称正定値の行列 A に対して Ax=b を「反復」で解く方法だよ。各ステップで残差 $r_k = b - A x_k$ と、過去のステップと A-直交(共役)になるような探索方向 $p_k$ を作って、その方向に最適な歩幅 $\alpha_k$ だけ進む。シミュレーターで「反復回数」を1から10まで動かして、青いカーブが下がっていくのを見てみて。残差ノルムが対数スケールでガクッと落ちる。
🙋
隣に赤い線もありますね。これは何ですか?
🎓
それが最急降下法だ。毎回「今の残差方向」だけを見て動く素朴な方法で、CG 法と直接比べると差が一目瞭然。条件数の大きい行列(A[0][0]を10にしてA[1][1]を1にしてみて)だと、最急降下法は谷の中をジグザグして全然下がらない。CG 法は「同じ方向を二度探さない」性質があるから、4次元の問題ならたった4反復で厳密解にピタッと着くんだ。
🙋
本当だ!4反復目で青いマーカーが急に落ちて、残差が 1e-15 とかになってますね。これって偶然ですか?
🎓
偶然じゃない、理論なんだ。$n$ 次元の対称正定値系では、CG 法は理想算術なら高々 $n$ 反復で厳密解に到達する——これは Krylov 部分空間 $\mathrm{span}\{r_0, Ar_0, A^2r_0, \ldots\}$ が n 次元になったら、その中に必ず真の解が含まれるから。実機の浮動小数点演算では丸め誤差で少し増えるけど、それでも最急降下法より段違いに速い。
🙋
条件数 cond(A) のカードが11くらいになってますね。これって何の数字ですか?
🎓
条件数は $\kappa = \lambda_\max / \lambda_\min$、最大固有値と最小固有値の比だ。これが大きいと等高線が細長い谷になって最急降下法はジグザグする。CG 法の誤差収束率は $\left(\frac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\right)^k$ で抑えられる——平方根が効くから、条件数100の問題でも、CG は最急降下法の10倍くらい速く下がる計算になるね。実務では「前処理付き CG (PCG)」で条件数を実効的に下げてさらに加速するのが標準だよ。

よくある質問

行列 A が対称(A^T = A)かつ正定値(任意の非零ベクトル x について x^T A x > 0)であることが必要です。実務では FEM の剛性行列、ラプラシアン、共分散行列など、自然に SPD になる問題が多くあります。非対称な系には Bi-CGSTAB や GMRES、対称だが不定な系には MINRES など、CG の派生・関連手法が使われます。
CG 法は行列 A をベクトルに掛ける操作(matVec)だけで進められるため、疎行列の構造を保ったまま解けます。直接法は分解後に「フィルイン」が発生してメモリが爆発するため、数百万自由度の FEM では現実的でないことが多く、PCG(前処理付き CG)が標準です。一方、複数の右辺 b に対して同じ A を解く場合は LU 分解の使い回しが効くため、直接法の方が有利な場面もあります。
M は「A に近く、かつ M^{-1} v が安く計算できる」もの。Jacobi(対角スケーリング)は最安だが効果は限定的、不完全コレスキー分解 (IC(0), ICT) は構造解析で広く使われます。流体や複雑な物理では代数的多重格子 (AMG) が条件数によらず O(n) で収束する強力な前処理子として選ばれます。Trilinos, PETSc などのライブラリには多数の前処理子が用意されています。
最も一般的なのは相対残差 ‖r_k‖_2 / ‖b‖_2 < tol で、tol は 1e-6 〜 1e-10 を取ります。ただし条件数が大きい場合、残差が小さくても真の誤差 ‖x_k - x_*‖ は大きいことがあるため、A-ノルム誤差や有効桁数の見積もりも併用します。最大反復回数(例: 2n, あるいは 500〜2000)でも停止させ、収束しない場合は前処理子の見直しが必要です。

実世界での応用

有限要素法 (FEM) による構造解析: 自動車・航空機・建築の応力解析では、節点変位を未知数とする大規模な剛性方程式 Kd=f を解きます。K は対称正定値かつ疎な行列で、数百万自由度の問題が日常的に現れます。前処理付き CG (IC-PCG, AMG-PCG) は商用ソルバー (Abaqus, Ansys, Nastran) や OSS (FreeFEM, FEniCS) で標準的に使われる主力解法です。

偏微分方程式の離散化: 熱伝導、ポアソン方程式、楕円型 PDE を有限差分や有限体積法で離散化すると、対称正定値の系が現れます。マルチグリッド法と CG の組み合わせは、地球流体や気象シミュレーションで条件数に依らない O(n) スケーラビリティを実現しています。

機械学習・最適化: 二次形式の最小化 min (1/2)x^T A x - b^T x は Ax=b の解と等価なので、CG 法は本質的に最適化アルゴリズムです。ニュートン法の各ステップでヘッシアン方程式を解く「ヘッシアンフリー最適化」では、行列を陽に持たずに CG で解く工夫が深層学習や強化学習で使われています。

画像処理・トモグラフィー: CT 再構成や画像復元では、観測モデルから生じる大規模な正規方程式 (A^T A) x = A^T b を解きます。A^T A は対称正定値(または半正定値+正則化)なので、共役勾配法系の解法 (CGLS, LSQR) が高速・低メモリの定番手法です。

よくある誤解と注意点

最も多い誤解は、「CG 法は n 反復で必ず厳密解になる」と無条件に信じることです。これは「理想算術」の話で、浮動小数点演算では丸め誤差で共役性 $p_i^\top A p_j = 0$ が徐々に崩れ、n 反復経過後も残差が残ることがあります。実機では「相対残差が許容値以下になるまで反復」という基準で停止し、必要に応じて「再起動 (restart)」や再直交化を行います。シミュレーターでは小さな 4 次元問題なので、丸め誤差の影響はほぼ見えませんが、数百万次元の実問題では条件数次第で挙動が大きく変わります。

次に多いのが、「条件数が大きいと CG はダメ」と諦めてしまうことです。確かに条件数 $\kappa$ が大きいほど収束は遅くなりますが、解は「前処理子 M」で大きく変えられます。実務では生 CG はほとんど使わず、不完全コレスキー (IC)、対称 Gauss-Seidel (SSOR)、代数的多重格子 (AMG) などで実効的な条件数を $\sqrt{\kappa}$ や $O(1)$ にまで下げます。「CG が遅い」のではなく「前処理子の選択が悪い」場合が大半です。シミュレーターで A[0][0]=10、A[1][1]=1 と差を付けると、最急降下法と CG の差が広がるのが見えますが、これに良い前処理子を入れると CG はさらに 1-2 反復で収束します。

最後に、非対称・不定の系に CG を適用してしまう失敗。CG は A が対称正定値であることを前提に設計されています。非対称行列 (例: 流体の対流項を含む系) に CG を素朴に当てると、$p_k^\top A p_k$ が負になったり、分母がゼロになって発散します。非対称には Bi-CGSTAB, GMRES, IDR(s)、対称不定には MINRES, SYMMLQ といった専用の Krylov 部分空間法を使うべきです。CG はあくまで SPD 専用の「最適化された解法」だと理解することが、安全な実装の第一歩です。