共役勾配法

数学において、共役勾配法は、行列が半正定値である特定の線形方程式系の数値解法のアルゴリズムです。共役勾配法は、多くの場合、反復アルゴリズムとして実装され、直接実装やコレスキー分解などの他の直接的な方法で処理するには大きすぎるスパースシステムに適用できます。大規模なスパースシステムは、偏微分方程式や最適化問題を 数値的に解くときによく発生します。

共役勾配法は、エネルギー最小化などの制約のない最適化問題を解くためにも使用できます。この法則は、マグナス・ヘステネスとエドゥアルド・シュティーフェル^{[1] [2]が}Z4でプログラムし、^[3]広範囲に研究したことで広く知られています。^{[4] [5]}

二重共役勾配法は、非対称行列への一般化を提供します。さまざまな非線形共役勾配法は、非線形最適化問題の最小値を求めます。

線形方程式の連立方程式を解きたいとします。

${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$

ベクトル に対して、既知の行列は対称行列(つまり、A ^T = A )、正定値行列(つまり、R ⁿのすべての非ゼロベクトルに対してx ^T Ax > 0 )、実数であり、既知です。このシステムの唯一の解を と表します。 ${\displaystyle \mathbf {x} }$${\displaystyle n\times n}$${\displaystyle \mathbf {A} }$${\displaystyle \mathbf {x} }$${\displaystyle \mathbf {b} }$${\displaystyle \mathbf {x} _{*}}$

共役勾配法は、最適化のための共役方向法の特殊化や、固有値問題のためのArnoldi / Lanczos反復法のバリエーションなど、いくつかの異なる観点から導出できます。アプローチは異なりますが、これらの導出には、残差の直交性と探索方向の共役性の証明という共通のトピックがあります。これら 2 つの特性は、この方法のよく知られた簡潔な定式化を開発する上で非常に重要です。

2つの非ゼロベクトルuとvが（に関して）共役であるとは、 ${\displaystyle \mathbf {A} }$

${\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =0.}$

は対称かつ正定値なので、左辺は内積を定義する。${\displaystyle \mathbf {A} }$

${\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} } :=\langle \mathbf {A} \mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} ^{\mathsf {T}}\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle .}$

2つのベクトルが共役であるのは、この内積に関して直交する場合に限ります。共役であることは対称的な関係です。つまり、がと共役であれば、 はと共役です。 ${\displaystyle \mathbf {u} }$${\displaystyle \mathbf {v} }$${\displaystyle \mathbf {v} }$${\displaystyle \mathbf {u} }$

${\displaystyle P=\{\mathbf {p} _{1},\dots ,\mathbf {p} _{n}\}}$

は、に関して相互に共役なベクトルの集合、つまりすべての に対して です。すると はの基底を形成し、この基底で の解を次のように表すことができます。${\displaystyle n}$${\displaystyle \mathbf {A} }$${\displaystyle \mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0}$${\displaystyle i\neq j}$${\displaystyle P}$${\displaystyle \mathbb {R} ^{n}}$${\displaystyle \mathbf {x} _{*}}$${\displaystyle \mathbf {Ax} =\mathbf {b} }$

${\displaystyle \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}\Rightarrow \mathbf {A} \mathbf { x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}.}$

この問題にベクトルを左 掛けすると次のようになる。${\displaystyle \mathbf {Ax} =\mathbf {b} }$${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}}$

${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {b} =\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}=\sum _{i=1}^{n}\alpha _{i}\left\langle \ mathbf {p} _{k},\mathbf {p} _{i}\right\rangle _{\mathbf {A} }=\alpha _{k}\left\langle \mathbf {p} _{k},\mathbf {p} _{k}\right\rangle _{\mathbf {A} }}$

など

${\displaystyle \alpha _{k}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\langle \mathbf {p} _{k},\mathbf { p} _{k}\rangle _{\mathbf {A} }}}.}$

これにより、方程式Ax = bを解くための次の方法^[4]が得られます。共役方向のシーケンスを見つけて、係数を計算します。 ${\displaystyle n}$${\displaystyle \alpha_{k}}$

共役ベクトルを注意深く選択すれば、解の近似値を得るためにすべての共役ベクトルが必要なくなるかもしれません。したがって、共役勾配法を反復法として考えます。これにより、 nが非常に大きく、直接法では時間がかかりすぎるような システムも近似的に解くことができます。${\displaystyle \mathbf {p} _{k}}$${\displaystyle \mathbf {x} _{*}}$

x _∗の初期推定値をx ₀で表します（一般性を失うことなくx ₀ = 0と仮定できますが、そうでない場合は代わりにシステムAz = b − Ax ₀を検討します）。x _{0から始めて解を探しますが、各反復で、解}x _∗ （私たちには未知）に近づいているかどうかを示すメトリックが必要です。このメトリックは、解x _∗が次の2次関数の唯一の最小値でもあるという事実から来ています。

${\displaystyle f(\mathbf {x} )={\tfrac {1}{2}}\mathbf {x} ^{\mathsf {T}}\mathbf {A} \mathbf {x} -\mathbf {x} ^{\mathsf {T}}\mathbf {b} ,\qquad \mathbf {x} \in \mathbf {R} ^{n}\,.}$

2次導関数のヘッセ行列が対称正定値である ことから、一意の最小化子の存在は明らかである。

${\displaystyle \mathbf {H} (f(\mathbf {x} ))=\mathbf {A} \,,}$

そして、最小化器（D f ( x )=0を使用）が初期問題を解くことは、その一次導関数から導かれる。

${\displaystyle \nabla f(\mathbf {x} )=\mathbf {A} \mathbf {x} -\mathbf {b} \,.}$

これは、最初の基底ベクトルp _{0 を}x = x ₀におけるfの勾配の負数とすることを示唆しています。 fの勾配はAx − bに等しくなります。初期推定値x ₀から始めて、これはp ₀ = b − Ax ₀とすることを意味します。基底の他のベクトルは勾配と共役になるため、共役勾配法と呼ばれます。 p _{0 は}、アルゴリズムのこの初期ステップによって提供される 残差でもあることに注意してください。

r _{k を}k番目のステップでの残差とします。

${\displaystyle \mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}.}$

上で観察されたように、は におけるの負の勾配であるため、勾配降下法ではr _{k の}方向に移動する必要があります。ただし、ここでは方向が互いに共役でなければならないことを主張します。これを実施する実用的な方法は、次の検索方向が現在の残差と以前のすべての検索方向から構築されるように要求することです。共役制約は正規直交型の制約であるため、アルゴリズムはグラム・シュミット正規直交化の例として見ることができます。これにより、次の式が得られます。 ${\displaystyle \mathbf {r} _{k}}$${\displaystyle f}$${\displaystyle \mathbf {x} _{k}}$${\displaystyle \mathbf {p} _{k}}$

${\displaystyle \mathbf {p} _{k}=\mathbf {r} _{k}-\sum _{i<k}{\frac {\mathbf {p} _{i}^{\mathsf {T }}\mathbf {A} \mathbf {r} _{k}}{\mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}}}\mathbf {p} _{i}}$

（共役制約が収束に与える影響については、記事の冒頭の図を参照）。この方向に従うと、次の最適位置は次のように与えられる。

${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}$

と

${\displaystyle \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{ \mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {p} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A } \mathbf {p} _{k}}},}$

ここで最後の等式は の定義から導かれる。 の式は、 x _{k +1の式を}fに代入し、 について最小化することで得られる。${\displaystyle \mathbf {r} _{k}}$${\displaystyle \alpha_{k}}$${\displaystyle \alpha_{k}}$

${\displaystyle {\begin{aligned}f(\mathbf {x} _{k+1})&=f(\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k})=:g(\alpha _{k})\\g'(\alpha _{k})&{\overset {!}{=}}0\quad \Rightarrow \quad \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathsf {T}}(\mathbf {b} -\mathbf {Ax} _{k})}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}\,.\end{aligned}}}$

上記のアルゴリズムは、共役勾配法を最も簡単に説明しています。 一見すると、このアルゴリズムは、以前のすべての検索方向と留数ベクトルの保存、および多くの行列とベクトルの乗算を必要とするため、計算コストが高くなる可能性があります。 ただし、アルゴリズムを詳しく分析すると、 はに直交していること、つまり、 i ≠ j の場合であることがわかります。 また、 はに -直交していること、つまり、の場合であることがわかります。 これは、アルゴリズムが進むにつれて、と が同じクリロフ部分空間に広がると見なすことができます。ここで、 は標準の内積に関して直交基底を形成し、は によって誘導される内積に関して直交基底を形成します。 したがって、 はのクリロフ部分空間への射影と見なすことができます。 ${\displaystyle \mathbf {r} _{i}}$${\displaystyle \mathbf {r} _{j}}$${\displaystyle \mathbf {r} _{i}^{\mathsf {T}}\mathbf {r} _{j}=0}$${\displaystyle \mathbf {p} _{i}}$${\displaystyle \mathbf {A} }$${\displaystyle \mathbf {p} _{j}}$${\displaystyle \mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0}$${\displaystyle i\neq j}$${\displaystyle \mathbf {p} _{i}}$${\displaystyle \mathbf {r} _{i}}$${\displaystyle \mathbf {r} _{i}}$${\displaystyle \mathbf {p} _{i}}$${\displaystyle \mathbf {A} }$${\displaystyle \mathbf {x} _{k}}$${\displaystyle \mathbf {x} }$

つまり、CG法が から始まる場合、^[6]を解くアルゴリズムは以下で詳しく説明されています。ここで、は実数、対称、正定値行列です。入力ベクトルは、近似初期解または0になります。これは、上で説明した正確な手順の異なる定式化です。 ${\displaystyle \mathbf {x} _{0}=0}$$${\displaystyle x_{k}=\mathrm {argmin} _{y\in \mathbb {R} ^{n}}{\left\{(x-y)^{\top }A(x-y):y\in \operatorname {span} \left\{b,Ab,\ldots ,A^{k-1}b\right\}\right\}}}$$${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$${\displaystyle \mathbf {A} }$${\displaystyle \mathbf {x} _{0}}$

${\displaystyle {\begin{aligned}&\mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}\\&{\hbox{if }}\mathbf {r} _{0}{\text{ is sufficiently small, then return }}\mathbf {x} _{0}{\text{ as the result}}\\&\mathbf {p} _{0}:=\mathbf {r} _{0}\\&k:=0\\&{\text{repeat}}\\&\qquad \alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}\\&\qquad \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}\\&\qquad \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}\\&\qquad {\hbox{if }}\mathbf {r} _{k+1}{\text{ is sufficiently small, then exit loop}}\\&\qquad \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}}\\&\qquad \mathbf {p} _{k+1}:=\mathbf {r} _{k+1}+\beta _{k}\mathbf {p} _{k}\\&\qquad k:=k+1\\&{\text{end repeat}}\\&{\text{return }}\mathbf {x} _{k+1}{\text{ as the result}}\end{aligned}}}$

これは最も一般的に使用されるアルゴリズムです。β _kの同じ式は、フレッチャー・リーブスの非線形共役勾配法でも使用されます。

は、 に勾配降下法を適用して計算されることに注意してください。 を設定すると、 は同様にから勾配降下法によって計算されます。つまり、 は共役勾配反復の再開の単純な実装として使用できます。^{[4]再開により収束が遅くなる可能性がありますが、共役勾配法が}丸め誤差 などにより誤動作した場合には、安定性が向上する可能性があります。 ${\displaystyle \mathbf {x} _{1}}$${\displaystyle \mathbf {x} _{0}}$${\displaystyle \beta _{k}=0}$${\displaystyle \mathbf {x} _{k+1}}$${\displaystyle \mathbf {x} _{k}}$

式と はどちらも正確な算術で成り立つため、式と は数学的に等価です。前者は、 を評価するためにベクトルがすでに計算されているため、による余分な乗算を避けるためにアルゴリズムで使用されます。後者は、丸め誤差の蓄積を伴う再帰によって暗黙的な計算を明示的な計算に置き換えるため、より正確である可能性があり、そのため、時々評価することを推奨します。^[7]${\displaystyle \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}$${\displaystyle \mathbf {r} _{k}:=\mathbf {b} -\mathbf {Ax} _{k}}$${\displaystyle \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}}$${\displaystyle \mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+1}}$${\displaystyle \mathbf {A} }$${\displaystyle \mathbf {Ap} _{k}}$${\displaystyle \alpha _{k}}$${\displaystyle \mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+1}}$

残差のノルムは、通常、停止基準として使用されます。明示的残差のノルムは、正確な演算と、収束が自然に停滞する丸め誤差が存在する場合の両方で、保証されたレベルの精度を提供します。対照的に、暗黙的残差は、丸め誤差のレベルをはるかに下回る振幅が小さくなり続けることが知られているため、収束の停滞を判断するために使用することはできません。 ${\displaystyle \mathbf {r} _{k+1}:=\mathbf {b} -\mathbf {Ax} _{k+1}}$${\displaystyle \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}}$

このアルゴリズムでは、α _{k は}に直交するように選択される。分母は次のように簡略化される。 ${\displaystyle \mathbf {r} _{k+1}}$${\displaystyle \mathbf {r} _{k}}$

${\displaystyle \alpha _{k}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}}$

であるからである。β k_はと共役となるように選ばれる。まず、β _kは ${\displaystyle \mathbf {r} _{k+1}=\mathbf {p} _{k+1}-\mathbf {\beta } _{k}\mathbf {p} _{k}}$${\displaystyle \mathbf {p} _{k+1}}$${\displaystyle \mathbf {p} _{k}}$

${\displaystyle \beta _{k}=-{\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}}}}$

使用して

${\displaystyle \mathbf {r} _{k+1}=\mathbf {r} _{k}-\alpha _{k}\mathbf {A} \mathbf {p} _{k}}$

そして同様に

${\displaystyle \mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}(\mathbf {r} _{k}-\mathbf {r} _{k+1}),}$

β _kの分子は次のように書き直される。

${\displaystyle \mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})=-{\frac {1}{\alpha _{k}}}\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {r} _{k+1}}$

とは設計上直交しているため、分母は次のように書き直される。 ${\displaystyle \mathbf {r} _{k+1}}$${\displaystyle \mathbf {r} _{k}}$

${\displaystyle \mathbf {p} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}=(\mathbf {r} _{k}+\beta _{k-1}\mathbf {p} _{k-1})^{\mathsf {T}}\mathbf {A} \mathbf {p} _{k}={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}(\mathbf {r} _{k}-\mathbf {r} _{k+1})={\frac {1}{\alpha _{k}}}\mathbf {r} _{k}^{\mathsf {T}}\mathbf {r} _{k}}$

探索方向p _kが共役であり、残差が直交していることを前提とします。これにより、α _{k を}キャンセルした後のアルゴリズムのβ が得られます。

「」
    共役勾配!(A, b, x)

共役勾配法を使用して、解を `A * x = b` に戻します。
「」
関数共役勾配! ( 
    A :: AbstractMatrix 、b :: AbstractVector 、x :: AbstractVector ; tol = eps ( eltype ( b ))   
）
    # 残差ベクトルを初期化する
    残差= b - A * x      
    # 検索方向ベクトルを初期化する
    search_direction =コピー(残余)  
    # 初期の二乗残差ノルムを計算する
	ノルム( x ) = sqrt (合計( x .^ 2 ))  
    old_resid_norm =ノルム(残差)  

    # 収束するまで繰り返す
    old_resid_norm > tolの場合   
        A_検索方向= A *検索方向    
        ステップサイズ=古い残差ノルム^ 2 / (検索方向' * A 検索方向)      
        # ソリューションの更新
        @. x = x +ステップサイズ*検索方向       
        # 残差を更新
        @.残差=残差-ステップサイズ* A_search_direction       
        new_resid_norm =ノルム(残差)  
        
        # 検索方向ベクトルを更新
        @ .search_direction =残余+     
            (新しい残留ノルム/古い残留ノルム) ^ 2 *検索方向    
        # 次の反復のために残差ノルムの二乗を更新します
        古い残留ノルム=新しい残留ノルム  
    終わり
    xを返す 
終わり

次式で表される 線形システムAx = bを考える。

${\displaystyle \mathbf {A} \mathbf {x} ={\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}},}$

共役勾配法の2つのステップを初期推定から実行します。

${\displaystyle \mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}}$

システムのおおよその解を見つけるためです。

参考までに、正確な解は

${\displaystyle \mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\\\{\frac {7}{11}}\end{bmatrix}}\approx {\begin{bmatrix}0.0909\\\\0.6364\end{bmatrix}}}$

最初のステップは、 x ₀に関連付けられた残差ベクトルr _{0 を計算することです。この残差は、} r ₀ = b - Ax ₀という式から計算され、この場合は次の式に等しくなります。

${\displaystyle \mathbf {r} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}-{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}2\\1\end{bmatrix}}={\begin{bmatrix}-8\\-3\end{bmatrix}}=\mathbf {p} _{0}.}$

これは最初の反復なので、残差ベクトルr _{0 を}初期検索方向p ₀として使用します。p _kの選択方法は、以降の反復で変更されます。

ここで、スカラーα ₀を次の関係式を使って 計算する。

${\displaystyle \alpha _{0}={\frac {\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}{\mathbf {p} _{0}^{\mathsf {T}}\mathbf {Ap} _{0}}}={\frac {{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}={\frac {73}{331}}\approx 0.2205}$

式を使って x _1を計算することができます

${\displaystyle \mathbf {x} _{1}=\mathbf {x} _{0}+\alpha _{0}\mathbf {p} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}+{\frac {73}{331}}{\begin{bmatrix}-8\\-3\end{bmatrix}}\approx {\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}.}$

この結果で最初の反復が完了し、システムx ₁の「改善された」近似解が得られます。次に、次の残差ベクトルr _{1 を}次の式を使用して 計算します。

${\displaystyle \mathbf {r} _{1}=\mathbf {r} _{0}-\alpha _{0}\mathbf {A} \mathbf {p} _{0}={\begin{bmatrix}-8\\-3\end{bmatrix}}-{\frac {73}{331}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}\approx {\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.}$

プロセスの次のステップは、最終的に次の検索方向p _{1 を}決定するために使用されるスカラーβ _{0 を}計算することです。

${\displaystyle \beta _{0}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}=0.0088.}$

さて、このスカラーβ ₀を使って、次の探索方向p _{1 を}次の関係式を使って 計算することができます。

${\displaystyle \mathbf {p} _{1}=\mathbf {r} _{1}+\beta _{0}\mathbf {p} _{0}\approx {\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}+0.0088{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}.}$

ここで、 α ₀に使用したのと同じ方法を使用して、新たに取得したp _{1 を}使用してスカラーα ₁を計算します。

${\displaystyle \alpha _{1}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {p} _{1}^{\mathsf {T}}\mathbf {Ap} _{1}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-0.3511&0.7229\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}}}=0.4122.}$

最後に、 x _{1 を}見つけるのと同じ方法を使用してx ₂を見つけます。

${\displaystyle \mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}\approx {\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}+0.4122{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}={\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.}$

結果x _2は、 x ₁やx ₀よりもシステムの解の「より良い」近似値です。この例で限定精度ではなく正確な演算を使用した場合、理論的にはn = 2 回の反復 ( nはシステムの次数) の後で正確な解に到達します。

共役勾配法は理論的には直接法とみなすことができます。丸め誤差がない場合、行列のサイズを超えない有限回数の反復後に正確な解が得られるからです。実際には、共役勾配法は小さな摂動に対しても不安定であるため、正確な解は決して得られません。たとえば、クリロフ部分空間を生成する際の退化の性質により、ほとんどの方向は実際には共役ではありません。

反復法として、共役勾配法は、正確な解への近似値を単調に（エネルギーノルムで）改善し、比較的少ない（問題のサイズに比べて）反復回数で必要な許容値に達する可能性があります。改善は通常線形であり、その速度はシステム行列の条件数によって決まります。条件数が大きいほど、改善は遅くなります。^[8]${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \kappa (A)}$${\displaystyle A}$${\displaystyle \kappa (A)}$

が大きい場合、 通常は、元のシステムをより小さくなるように置き換える前処理が使用されます(以下を参照)。 ${\displaystyle \kappa (A)}$${\displaystyle \mathbf {Ax} -\mathbf {b} =0}$${\displaystyle \mathbf {M} ^{-1}(\mathbf {Ax} -\mathbf {b} )=0}$${\displaystyle \kappa (\mathbf {M} ^{-1}\mathbf {A} )}$${\displaystyle \kappa (\mathbf {A} )}$

多項式のサブセットを次のように定義する。

${\displaystyle \Pi _{k}^{*}:=\left\lbrace \ p\in \Pi _{k}\ :\ p(0)=1\ \right\rbrace \,,}$

ここで、最大次数の多項式の集合です。 ${\displaystyle \Pi _{k}}$${\displaystyle k}$

を厳密解の反復近似とし、誤差を と定義する。ここで、収束率は次のように近似できる^{[4] [9]}${\displaystyle \left(\mathbf {x} _{k}\right)_{k}}$${\displaystyle \mathbf {x} _{*}}$${\displaystyle \mathbf {e} _{k}:=\mathbf {x} _{k}-\mathbf {x} _{*}}$

${\displaystyle {\begin{aligned}\left\|\mathbf {e} _{k}\right\|_{\mathbf {A} }&=\min _{p\in \Pi _{k}^{*}}\left\|p(\mathbf {A} )\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq \min _{p\in \Pi _{k}^{*}}\,\max _{\lambda \in \sigma (\mathbf {A} )}|p(\lambda )|\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\left({\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\exp \left({\frac {-2k}{\sqrt {\kappa (\mathbf {A} )}}}\right)\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\,,\end{aligned}}}$

ここで はスペクトルを表し、 は条件数を表します。 ${\displaystyle \sigma (\mathbf {A} )}$${\displaystyle \kappa (\mathbf {A} )}$

これは、任意の に対して、反復が誤差をまで減らすのに十分であることを示しています。 ${\displaystyle k={\tfrac {1}{2}}{\sqrt {\kappa (\mathbf {A} )}}\log \left(\left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\varepsilon ^{-1}\right)}$${\displaystyle 2\varepsilon }$${\displaystyle \varepsilon >0}$

注意してください、重要な限界は${\displaystyle \kappa (\mathbf {A} )}$${\displaystyle \infty }$

${\displaystyle {\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\approx 1-{\frac {2}{\sqrt {\kappa (\mathbf {A} )}}}\quad {\text{for}}\quad \kappa (\mathbf {A} )\gg 1\,.}$

この極限は、 としてスケールするヤコビ法やガウス・ザイデル法の反復法と比較して、より速い収束率を示します。 ${\displaystyle \approx 1-{\frac {2}{\kappa (\mathbf {A} )}}}$

収束定理では丸め誤差は想定されていないが、アン・グリーンバウム^[5]によって理論的に説明されているように、収束限界は実際には一般的に有効である。

ランダムに初期化された場合、誤差は最初はより小さい有効条件数を反映するクリロフ部分空間内で除去されるため、最初の段階の反復は多くの場合最も高速になります。収束の2番目の段階は、通常、 による理論的な収束境界によって明確に定義されますが、行列のスペクトルの分布と誤差のスペクトル分布によっては、超線形になることがあります。^[5]最後の段階では、達成可能な最小の精度に達し、収束が停止するか、または方法が発散し始めることもあります。大規模な行列の倍精度浮動小数点形式の一般的な科学計算アプリケーションでは、共役勾配法は、最初の段階または2番目の段階で反復を終了する許容値を持つ停止基準を使用します。 ${\textstyle {\sqrt {\kappa (\mathbf {A} )}}}$${\displaystyle A}$

ほとんどの場合、共役勾配法の高速収束を確実にするために前処理が必要です。が対称正定値で、よりも条件数が良い場合は、前処理付き共役勾配法を使用できます。 これは次の形式になります。^[10]${\displaystyle \mathbf {M} ^{-1}}$${\displaystyle \mathbf {M} ^{-1}\mathbf {A} }$${\displaystyle \mathbf {A} }$

${\displaystyle \mathbf {r} _{0}:=\mathbf {b} -\mathbf {Ax} _{0}}$

${\displaystyle {\textrm {Solve:}}\mathbf {M} \mathbf {z} _{0}:=\mathbf {r} _{0}}$

${\displaystyle \mathbf {p} _{0}:=\mathbf {z} _{0}}$

${\displaystyle k:=0\,}$

繰り返す

${\displaystyle \alpha _{k}:={\frac {\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}{\mathbf {p} _{k}^{\mathsf {T}}\mathbf {Ap} _{k}}}}$

${\displaystyle \mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}}$

${\displaystyle \mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}\mathbf {Ap} _{k}}$

r _{k +1}が十分に小さい場合はループを終了し、

${\displaystyle \mathrm {Solve} \ \mathbf {M} \mathbf {z} _{k+1}:=\mathbf {r} _{k+1}}$

${\displaystyle \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}}}$

${\displaystyle \mathbf {p} _{k+1}:=\mathbf {z} _{k+1}+\beta _{k}\mathbf {p} _{k}}$

${\displaystyle k:=k+1\,}$

繰り返し終了

結果はx _{k +1}

上記の定式化は、前処理されたシステムに通常の共役勾配法を適用することと同等である^[11]

${\displaystyle \mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf {T}}\mathbf {\hat {x}} =\mathbf {E} ^{-1}\mathbf {b} }$

どこ

${\displaystyle \mathbf {EE} ^{\mathsf {T}}=\mathbf {M} ,\qquad \mathbf {\hat {x}} =\mathbf {E} ^{\mathsf {T}}\mathbf {x} .}$

システムの対称性（および正定値性）を保つには、前処理のコレスキー分解を使用する必要があります。ただし、この分解を計算する必要はなく、 を知っていれば十分です。が と同じスペクトルを持つことを示すことができます。 ${\displaystyle \mathbf {M} ^{-1}}$${\displaystyle \mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\mathsf {T}}}$${\displaystyle \mathbf {M} ^{-1}\mathbf {A} }$

前処理行列M は対称正定値で固定されている必要があります。つまり、反復ごとに変更することはできません。前処理行列に関するこれらの仮定のいずれかが違反すると、前処理付き共役勾配法の動作が予測不可能になる可能性があります。

よく使われる前処理の例としては不完全コレスキー分解がある。^[12]

プロセスで使用するために行列を明示的に反転することは望ましくないことに留意することが重要です。反転には、共役勾配アルゴリズム自体を解くよりも多くの時間/計算リソースがかかるからです。例として、不完全コレスキー分解から得られる前処理行列を使用するとします。結果の行列は下三角行列であり、前処理行列は次のようになります。 ${\displaystyle \mathbf {M} }$${\displaystyle \mathbf {M} ^{-1}}$${\displaystyle \mathbf {M} }$${\displaystyle \mathbf {L} }$

${\displaystyle \mathbf {M} =\mathbf {LL} ^{\mathsf {T}}}$

次に解決しなければならないのは、次の点です。

${\displaystyle \mathbf {Mz} =\mathbf {r} }$

${\displaystyle \mathbf {z} =\mathbf {M} ^{-1}\mathbf {r} }$

しかし：

${\displaystyle \mathbf {M} ^{-1}=(\mathbf {L} ^{-1})^{\mathsf {T}}\mathbf {L} ^{-1}}$

それから：

${\displaystyle \mathbf {z} =(\mathbf {L} ^{-1})^{\mathsf {T}}\mathbf {L} ^{-1}\mathbf {r} }$

中間ベクトルを取りましょう: ${\displaystyle \mathbf {a} }$

${\displaystyle \mathbf {a} =\mathbf {L} ^{-1}\mathbf {r} }$

${\displaystyle \mathbf {r} =\mathbf {L} \mathbf {a} }$

およびは既知であり、 は下三角なので、 を解くのは前方置換を使用することで簡単かつ計算コストが低くなります。次に、元の方程式に を代入します。${\displaystyle \mathbf {r} }$${\displaystyle \mathbf {L} }$${\displaystyle \mathbf {L} }$${\displaystyle \mathbf {a} }$${\displaystyle \mathbf {a} }$

${\displaystyle \mathbf {z} =(\mathbf {L} ^{-1})^{\mathsf {T}}\mathbf {a} }$

${\displaystyle \mathbf {a} =\mathbf {L} ^{\mathsf {T}}\mathbf {z} }$