共役勾配法

数理最適化アルゴリズム

与えられた線形システムに関連する二次関数を最小化する、最適ステップサイズ（緑）と共役ベクトル（赤）を用いた勾配降下法の収束の比較。共役勾配法は、正確な演算を仮定すると、最大nステップで収束します。ここで、nはシステムの行列のサイズです（ここではn  = 2）。

数学において、共役勾配法は、特定の線形方程式系、特に半正定値行列を持つ線形方程式系の数値解を求めるアルゴリズムである。共役勾配法は、反復アルゴリズムとして実装されることが多く、直接実装やコレスキー分解などの直接法では扱えないほど大規模な疎行列に適用できる。大規模な疎行列は、偏微分方程式や最適化問題を数値的に解く際によく用いられる。

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

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

共役勾配法で扱う問題の説明

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

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

ベクトル に対して、既知の行列は対称行列（すなわち）、正定値行列（すなわち内の零でないすべてのベクトルに対して）、実数であり、 も既知である。この系の唯一の解を と表記する。 ${\displaystyle \mathbf {x} }$ ${\displaystyle n\times n}$ ${\displaystyle \mathbf {A} }$ ${\displaystyle \mathbf {A} ^{\mathsf {T}}=\mathbf {A} }$ ${\displaystyle \mathbf {x} ^{\mathsf {T}}\mathbf {Ax} >0}$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbb {R} ^{n}}$ ${\displaystyle \mathbf {b} }$ ${\displaystyle \mathbf {x} _{*}}$

直接的な方法としての導出

共役勾配法は、最適化のための共役方向法の特化や、固有値問題のためのアーノルディ／ランチョス反復法の変形など、様々な観点から導出可能です。アプローチは異なりますが、これらの導出には共通のテーマ、すなわち残差の直交性と探索方向の共役性の証明があります。この2つの性質は、この手法のよく知られた簡潔な定式化を開発する上で極めて重要です。

2つの非零ベクトルとが（に関して）共役であるとは、 ${\displaystyle \mathbf {u} }$ ${\displaystyle \mathbf {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} }}}.}$

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

反復的な方法として

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

の初期推定値を と表記する（一般性を失うことなく と仮定できる。そうでない場合は、代わりに を考える）。 から始めて解を探索し、各反復において、解に近づいているかどうか（これは未知である）を示す指標が必要となる。この指標は、解が次の二次関数の唯一の最小値でもあるという事実から得られる。 ${\displaystyle \mathbf {x} _{*}}$ ${\displaystyle \mathbf {x} _{0}}$ ${\displaystyle \mathbf {x} _{0}=\mathbf {0} }$ ${\displaystyle \mathbf {Az} =\mathbf {b} -\mathbf {Ax} _{0}}$ ${\displaystyle \mathbf {x} _{0}}$ ${\displaystyle \mathbf {x} _{*}}$ ${\displaystyle \mathbf {x} _{*}}$

${\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 \mathbb {R} ^{n}\,.}$

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

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

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

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

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

番目のステップにおける残差を次のようにします。 ${\displaystyle \mathbf {r} _{k}}$ ${\displaystyle 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 {r} _{k}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{i}}{\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}}}$

結果として得られるアルゴリズム

上記のアルゴリズムは、共役勾配法を最も簡単に説明しています。 一見すると、このアルゴリズムは、以前のすべての探索方向と留数ベクトルの保存、および多くの行列とベクトルの乗算を必要とするため、計算コストが高くなる可能性があります。 しかし、アルゴリズムを詳しく分析すると、 は に直交していること、つまりの場合 であることがわかります。また、はに -直交していること、つまり の場合 です。 これは、アルゴリズムが進むにつれて とみなすことができ、は同じクリロフ部分空間に広がります。ここで、 は標準の内積 に関する直交基底を形成し、は によって誘導される内積 に関する直交基底を形成します。 したがって、は のクリロフ部分空間への射影と見なすことができます。 ${\displaystyle \mathbf {r} _{i}}$ ${\displaystyle \mathbf {r} _{j}}$ ${\displaystyle \mathbf {r} _{i}^{\mathsf {T}}\mathbf {r} _{j}=0}$ ${\displaystyle i\neq j}$ ${\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]を解くアルゴリズムは以下で詳述されます。ここでは実対称正定値行列です。入力ベクトルは近似初期解または です。これは、上で説明した手順の正確な表現とは異なります。 ${\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 \mathbf {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}}}$

これは最も一般的に用いられるアルゴリズムです。 の同じ式は、フレッチャー・リーブスの非線形共役勾配法でも用いられます。 ${\displaystyle \beta _{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}}$

アルファとベータの計算

アルゴリズムでは、が に直交するように選択される。分母は次のように簡略化される。 ${\displaystyle \alpha _{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}}}}$

であるからである。 はと共役となるように選ばれる。最初はは ${\displaystyle \mathbf {r} _{k+1}=\mathbf {p} _{k+1}-\mathbf {\beta } _{k}\mathbf {p} _{k}}$ ${\displaystyle \beta _{k}}$ ${\displaystyle \mathbf {p} _{k+1}}$ ${\displaystyle \mathbf {p} _{k}}$ ${\displaystyle \beta _{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}),}$

の分子は次のように書き直される。 ${\displaystyle \beta _{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}}$

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

サンプルコードJulia（プログラミング言語）

線形代数を使用する 「」 x = 共役勾配(A, b, x0 = ゼロ(b); atol = 長さ(b) * eps(ノルム(b))共役勾配法を使用して、解を `A * x = b` に返します。`A` は正定値行列またはその他の線形演算子である必要があります。`x0` は解の初期推定値です (デフォルトはゼロ ベクトルです)。`atol`は残差`b - A * x`の大きさの絶対許容値である。収束のため（デフォルトはマシンイプシロン）。近似解ベクトル `x` を返します。「」関数共役勾配(  A , b :: AbstractVector , x0 :: AbstractVector = zero ( b ); atol = length ( b ) * eps ( norm ( b ))     ） x = copy ( x0 ) # ソリューションを初期化する    r = b - A * x0 # 初期残差        p =コピー( r ) # 初期探索方向    r²old = r ' * r # 残差の2乗ノルム      k = 0   while r²old > atol ^ 2 # 収束するまで反復する     Ap = A * p # 検索方向      α = r²old / ( p ' * Ap ) # ステップサイズ        @。x += α * p # 解を更新する       # 残差を更新: if ( k + 1 ) % 16 == 0 # 16回の反復ごとに残差を最初から再計算する         r .= b .- A * x # 数値誤差の蓄積を避けるため        それ以外 @. r -= α * Ap # 行列ベクトル積を1つ節約する更新式を使用する       終わり r²new = r ' * r     @. p = r + ( r²new / r²old ) * p # 検索方向を更新           r²old = r²new # 残差ノルムの二乗を更新    k += 1   終わり xを返す 終わり

サンプルコードMATLAB

関数 x =共役勾配( A, b, x0, tol )  % 共役勾配法を使用して `A * x = b` の解を返します。% 注意: A は対称かつ正定値である必要があります。 ナルギン< 4の場合    tol = eps ;   終わり r = b - A * x0 ;       p = r ;   rsold = r ' * r ;     x = x0 ;   sqrt ( rsold ) > tolであるとき    Ap = A * p ;     アルファ= rsold / ( p ' * Ap );       x = x +アルファ* p ;       r = r -アルファ* Ap ;       rsnew = r ' * r ;     p = r + ( rsnew / rsold ) * p ;         rsold = rsnew ;   終わり終わり

数値例

次式で表される線形システム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 _{0として使用します。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 ₁を計算することができます

${\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 ₁を次の式で 計算します。

${\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 ₁を決定するために使用されるスカラーβ ₀を計算することです。

${\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 ₁を次の関係式を使って 計算することができます。

${\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 ₁を使用してスカラーα ₁を計算します。

${\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 ₁を求めたのと同じ方法を使用して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はシステムの次数）で正確な解に到達していたはずです。

有限終了性

正確な演算では、必要な反復回数は行列の次数以下です。この動作は、共役勾配法の有限停止性として知られています。これは、正確な演算を用いると、線形システムの厳密解に有限回数（システムの次元に等しい回数まで）で到達できるという手法の能力を指します。この特性は、各反復において、この手法がそれ以前のすべての残差と直交する残差ベクトルを生成するという事実に由来します。これらの残差は互いに直交する集合を形成します。

n次元空間では、 n個以上の線形独立かつ互いに直交するベクトルを構成することは、そのうちの1つが零ベクトルでない限り不可能である。したがって、零残差が現れると、この方法は解に到達したことになり、必ず終了する。これにより、共役勾配法は最大nステップで収束することが保証される。

これを説明するために、次のシステムを考えてみましょう。

${\displaystyle A={\begin{bmatrix}3&-2\\-2&4\end{bmatrix}},\quad \mathbf {b} ={\begin{bmatrix}1\\1\end{bmatrix}}}$

初期推定値 から始めます。は対称正定値であり、系は2次元であるため、共役勾配法は2ステップ以内で厳密解を求めるはずです。次のMATLABコードはこの動作を示しています。 ${\displaystyle \mathbf {x} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}}$ ${\displaystyle A}$

A = [ 3 , - 2 ; - 2 , 4 ]; x_true = [ 1 ; 1 ]; b = A * x_true ;            x = [ 1 ; 2 ]; % 初期推測値r = b - A * x ; p = r ;            k = 1 : 2 の場合Ap = A * p ; alpha = ( r ' * r ) / ( p ' * Ap ); x = x + alpha * p ; r_new = r - alpha * Ap ; beta = ( r_new ' * r_new ) / ( r ' * r ); p = r_new + beta * p ; r = r_new ;終了                                                  disp ( '正確な解:' ); disp ( x );

出力は、この手法が2回の反復後に到達し、理論予測と一致していることを示しています。この例は、共役勾配法が理想的な条件下で直接法としてどのように動作するかを示しています。 ${\displaystyle {\begin{bmatrix}1\\1\end{bmatrix}}}$

スパースシステムへの応用

有限停止性は、科学技術の応用において頻繁に生じる大規模な疎システムを解く際にも実用的な意味を持ちます。例えば、一様格子上で有限差分を用いて2次元ラプラス方程式を離散化すると、疎線形システム が得られます。ここで、 は対称かつ正定値です。 ${\displaystyle \nabla ^{2}u=0}$ ${\displaystyle A\mathbf {x} =\mathbf {b} }$ ${\displaystyle A}$

内部グリッドを用いると、係数行列は5点ステンシルパターンを持つ系が得られます。各行には、中心点とその近傍に対応する最大5つの非零要素が含まれます。例えば、このようなグリッドから生成される行列は次のようになります。 ${\displaystyle 5\times 5}$ ${\displaystyle 25\times 25}$ ${\displaystyle A}$ ${\displaystyle A}$

${\displaystyle A={\begin{bmatrix}4&-1&0&\cdots &-1&0&\cdots \\-1&4&-1&\cdots &0&0&\cdots \\0&-1&4&-1&0&0&\cdots \\\vdots &\vdots &\ddots &\ddots &\ddots &\vdots \\-1&0&\cdots &-1&4&-1&\cdots \\0&0&\cdots &0&-1&4&\cdots \\\vdots &\vdots &\cdots &\cdots &\cdots &\ddots \end{bmatrix}}}$

システムの次元は25ですが、共役勾配法は理論上、正確な演算の下では最大25回の反復で収束することが保証されています。実際には、行列のスペクトル特性により、はるかに少ないステップで収束することがよくあります。この効率性により、共役勾配法は、熱伝導、流体力学、静電気学などの偏微分方程式から生じる大規模システムを解く際に特に魅力的です。

収束特性

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

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

しかし、大きな対称行列で固有値が対数的に分布している場合には、興味深いケースが生じます。たとえば、 がランダム直交行列で が対角行列で固有値が から までの範囲にあり、対数的に分布しているとします。CGM には有限停止特性があり、理論的には最大 ステップで正確な解に到達するはずですが、この手法では収束が停滞する場合があります。このようなシナリオでは、さらに多くの反復処理 (たとえば、行列サイズの 10 倍) を行った後でも、誤差はわずかにしか減少しない場合があります (たとえば、 から)。さらに、反復誤差が大きく振動する場合があり、停止条件としては信頼できません。この収束不良は条件数だけ (たとえば) によって説明されるのではなく、固有値の分布自体によって説明されます。固有値がより均等に間隔をあけられているかランダムに分布している場合、このような収束の問題は通常は発生せず、CGM のパフォーマンスは だけでなく、固有値の分布方法にも依存することが強調されます。^[9] ${\displaystyle A=QDQ^{T}}$ ${\displaystyle Q}$ ${\displaystyle D}$ ${\displaystyle \lambda _{n}=1}$ ${\displaystyle \lambda _{1}=10^{6}}$ ${\displaystyle n}$ ${\displaystyle 10^{-5}}$ ${\displaystyle \kappa _{2}(A)=10^{6}}$ ${\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] [10]} ${\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\,.}$

この極限は、 としてスケールするJacobiやGauss–Seidelの反復法に比べて、より速い収束速度を示します。 ${\displaystyle \approx 1-{\frac {2}{\kappa (\mathbf {A} )}}}$

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

実践的な収束

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

前処理共役勾配法

ほとんどの場合、共役勾配法の高速収束を保証するためには前処理が必要です。が対称正定値で、かつ条件数が前処理共役勾配法よりも優れている場合、前処理共役勾配法を使用することができます。これは以下の形式をとります。^[11] ${\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}

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

${\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} }$

前処理行列は対称正定値かつ固定でなければなりません。つまり、反復ごとに変化してはいけません。前処理行列に関するこれらの仮定のいずれかが破られると、前処理付き共役勾配法の挙動は予測不可能になる可能性があります。 ${\displaystyle \mathbf {M} }$

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

前処理子の実際の使用

重要なのは、この処理で使用するために行列を明示的に逆行列化する必要はないということです。逆行列化は、共役勾配法自体を解くよりも多くの時間と計算リソースを必要とするからです。例えば、不完全コレスキー分解から得られる前処理行列を使用するとします。結果の行列は下三角行列 であり、前処理行列は次のようになります。 ${\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} }$

およびは既知であり、 は上三角関数であるため、 を解くことは、後方代入を使用することで簡単かつ計算コストが低くなります。 ${\displaystyle \mathbf {a} }$ ${\displaystyle \mathbf {L} ^{\mathsf {T}}}$ ${\displaystyle \mathbf {L} ^{\mathsf {T}}}$ ${\displaystyle \mathbf {z} }$

この方法を使用すると、または を明示的に反転する必要はまったくなく、 が得られます。 ${\displaystyle \mathbf {M} }$ ${\displaystyle \mathbf {L} }$ ${\displaystyle \mathbf {z} }$

柔軟な前処理共役勾配法

数値的に難しい応用では、複雑な前処理が用いられますが、これは反復ごとに変化する前処理変数につながる可能性があります。前処理変数が反復ごとに対称正定値であっても、変化する可能性があるという事実は上記の議論を無効にし、実際のテストでは上記のアルゴリズムの収束を著しく遅くします。Polak -Ribièreの式 を用いると、

${\displaystyle \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\left(\mathbf {z} _{k+1}-\mathbf {z} _{k}\right)}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}}}$

フレッチャー・リーブス式の代わりに

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

この場合、収束性が大幅に改善される可能性がある。^[14]この前処理付き共役勾配法のバージョンは、可変の前処理を許容するため、 ^[15] 柔軟であると言える。また、この柔軟版は、前処理が対称正定値（SPD）でなくても堅牢であることが示されている^[16]。

柔軟なバージョンの実装には、追加のベクトルを保存する必要があります。固定SPD前処理の場合、β _kの両式は正確な算術演算、つまり丸め誤差なしで等価です。 ${\displaystyle \mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k}=0,}$

ポラック・リビエールの式を用いた方法の収束挙動が優れている理由の数学的説明は、この場合、この方法が局所最適であり、特に、局所最適な最急降下法よりも収束が遅くならないということである。^[17]

局所最適最急降下法との比較

元の共役勾配法と前処理付き共役勾配法の両方において、 直線探索法、最急降下法を用いて局所的に最適にするには、を設定するだけで済みます。この置換により、ベクトルpは常にベクトルzと同じになるため、ベクトルpを保存する必要がなくなります。したがって、これらの最急降下法の各反復処理は、共役勾配法に比べて若干コストが低くなります。ただし、後者は、（高度に）変数の多い前処理やSPD以外の前処理を使用しない限り、収束が速くなります（上記参照）。 ${\displaystyle \beta _{k}:=0}$

二重積分器の最適フィードバック制御器としての共役勾配法

共役勾配法は最適制御理論を用いて導くこともできる。^[18]このアプローチでは、共役勾配法は二重積分器システムに対する最適フィードバック制御器として導出される。量と量は可変フィードバックゲインである。^[18] $${\displaystyle u=k(x,v):=-\gamma _{a}\nabla f(x)-\gamma _{b}v}$$ $${\displaystyle {\dot {x}}=v,\quad {\dot {v}}=u}$$ ${\displaystyle \gamma _{a}}$ ${\displaystyle \gamma _{b}}$

正規方程式の共役勾配

共役勾配法は、任意のn行m列の行列に適用することができます。これは、正規方程式 A ^T Aと右辺ベクトルA ^T bに適用することで可能です。これは、 A ^T Aが任意のAに対して対称な半正定値行列であるため です。結果は、正規方程式（CGNまたはCGNR）上の共役勾配となります。

A ^T Ax = A ^T b

反復法であるため、メモリ上でA ^T Aを明示的に形成する必要はなく 、行列-ベクトル乗算と転置行列-ベクトル乗算を実行するだけで済みます。したがって、CGNR法は、Aが疎行列である場合に特に有用です。これらの演算は通常非常に効率的だからです。しかし、正規方程式を形成することの欠点は、条件数κ( A ^T A ) が κ ² ( A )に等しいため、CGNR法の収束速度が遅くなり、近似解の品質が丸め誤差の影響を受けやすくなる可能性があることです。適切な前処理を見つけることは、CGNR法を使用する上でしばしば重要な部分となります。

いくつかのアルゴリズムが提案されています（例：CGLS、LSQR）。LSQRアルゴリズムは、Aが悪条件、つまりAの条件数が大きい場合に、最も数値安定性が高いと言われています。

複素エルミート行列の共役勾配法

共役勾配法は、簡単な修正を加えることで、複素数値行列Aとベクトルbが与えられた場合に、複素数値ベクトルxに関する線形方程式系を解くことに拡張できます。ここで、Aはエルミート行列（すなわちA' = A）で正定値行列であり、記号'は共役転置を表します。簡単な修正は、実転置を共役転置ですべて置き換えるだけです。 ${\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }$

メリットとデメリット

共役勾配法の利点と欠点は、ネミロフスキーとベンタルの講義ノートにまとめられている。^{[19] ：Sec.7.3}

病的な例

この例は^[20]からの引用です。 とし、 を定義します。は可逆なので、 には一意の解が存在します。 これを共役勾配降下法で解くと、収束性が悪くなります。つまり、CGプロセス中、誤差は指数関数的に増加し、一意の解が見つかると突然ゼロになります。 ${\textstyle t\in (0,1)}$ $${\displaystyle W={\begin{bmatrix}t&{\sqrt {t}}&&&&\\{\sqrt {t}}&1+t&{\sqrt {t}}&&&\\&{\sqrt {t}}&1+t&{\sqrt {t}}&&\\&&{\sqrt {t}}&\ddots &\ddots &\\&&&\ddots &&\\&&&&&{\sqrt {t}}\\&&&&{\sqrt {t}}&1+t\end{bmatrix}},\quad b={\begin{bmatrix}1\\0\\\vdots \\0\end{bmatrix}}}$$ ${\displaystyle W}$ ${\textstyle Wx=b}$ $${\displaystyle \|b-Wx_{k}\|^{2}=(1/t)^{k},\quad \|b-Wx_{n}\|^{2}=0}$$

参照

• 二共役勾配法（BiCG）
• 共役勾配二乗法（CGS）
• 共役残差法
• ガウス確率伝播
• 反復法：線形システム
• クリロフ部分空間
• 非線形共役勾配法
• 前処理
• 疎行列とベクトルの乗算

参考文献

1. ヘステネス、マグナス・R. ;スティフェル、エドゥアルド(1952年12月). 「線形システムの解法における共役勾配法」(PDF) .米国国立標準局研究ジャーナル. 49 (6): 409. doi : 10.6028/jres.049.044 .
2. Straeter, TA (1971). 「ランク1のDavidon–Broydenクラスの準ニュートン最小化法の無限次元ヒルベルト空間への拡張と最適制御問題への応用について」（博士論文）ノースカロライナ州立大学. hdl : 2060/19710026200 – NASA技術報告サーバー経由.
3. シュパイザー、アンブロス(2004)。 「Konrad Zuse und die ERMETH: Ein weltweiter Architektur-Vergleich」 [Konrad Zuse と ERMETH: 建築の世界的比較]。 『Hellige』、ハンス・ディーター編著。情報収集。 Visionen、Paradigmen、Leitmotive (ドイツ語)。ベルリン：シュプリンガー。 p. 185.ISBN​ 3-540-00217-0。
4. Polyak, Boris (1987). 最適化入門.
5. Greenbaum, Anne (1997).線形システムの反復解法. doi :10.1137/1.9781611970937. ISBN 978-0-89871-396-1。
6. パケット, エリオット; トログドン, トーマス (2023年3月). 「標本共分散行列における共役勾配法とMINRESアルゴリズムの普遍性」. Communications on Pure and Applied Mathematics . 76 (5): 1085–1136 . arXiv : 2007.00640 . doi :10.1002/cpa.22081. ISSN  0010-3640.
7. Shewchuk, Jonathan R (1994). 苦痛を伴わない共役勾配法入門(PDF) .
8. Saad, Yousef (2003). 疎線形システムのための反復法（第2版）. ペンシルベニア州フィラデルフィア: Society for Industrial and Applied Mathematics. pp. 195. ISBN 978-0-89871-534-7。
9. Holmes, M. (2023). 『科学計算とデータ分析入門』第2版. Springer. ISBN 978-3-031-22429-4。
10. Hackbusch, W. (2016-06-21).大規模疎方程式系の反復解法（第2版）. スイス: Springer. ISBN 978-3-319-28483-5. OCLC  952572240.
11. Barrett, Richard; Berry, Michael; Chan, Tony F.; Demmel, James; Donato, June; Dongarra, Jack; Eijkhout, Victor; Pozo, Roldan; Romine, Charles; van der Vorst, Henk. 『線形システムの解法テンプレート：反復法のビルディングブロック』（PDF） （第2版）フィラデルフィア、ペンシルベニア州：SIAM​​。13ページ。 2020年3月31日閲覧。
12. Golub, Gene H.; Van Loan, Charles F. (2013). Matrix Computations (第4版). Johns Hopkins University Press. sec. 11.5.2. ISBN 978-1-4214-0794-4。
13. Concus, P.; Golub, GH; Meurant, G. (1985). 「共役勾配法のためのブロック前処理」. SIAM Journal on Scientific and Statistical Computing . 6 (1): 220– 252. doi :10.1137/0906018.
14. Golub, Gene H.; Ye, Qiang (1999). 「Inexact Preconditioned Conjugate Gradient Method with Inner-Outer Iteration」. SIAM Journal on Scientific Computing . 21 (4): 1305. CiteSeerX 10.1.1.56.1755 . doi :10.1137/S1064827597323415. 
15. Notay, Yvan (2000). 「柔軟な共役勾配」. SIAM Journal on Scientific Computing . 22 (4): 1444– 1460. CiteSeerX 10.1.1.35.7473 . doi :10.1137/S1064827599362314. 
16. Bouwmeester, Henricus; Dougherty, Andrew; Knyazev, Andrew V. (2015). 「共役勾配法と最急降下法1における非対称前処理」. Procedia Computer Science . 51 : 276–285 . arXiv : 1212.6680 . doi : 10.1016/j.procs.2015.05.241 . S2CID  51978658.
17. Knyazev, Andrew V.; Lashuk, Ilya (2008). 「変数前処理を用いた最急降下法と共役勾配法」. SIAM Journal on Matrix Analysis and Applications . 29 (4): 1267. arXiv : math/0605767 . doi :10.1137/060675290. S2CID  17614913.
18. Ross, IM、「加速最適化のための最適制御理論」、arXiv : 1902.09004、2019年。
19. NemirovskyとBen-Tal (2023).「最適化III：凸最適化」（PDF） .
20. ペニントン、ファビアン・ペドレゴサ、コートニー・パケット、トム・トログドン、ジェフリー。「ランダム行列理論と機械学習チュートリアル」。random -matrix -learning.github.io 。 2023年12月5日閲覧。{{cite web}}: CS1 maint: 複数の名前: 著者リスト (リンク)

さらに読む

• アトキンソン, ケンデル A. (1988). 「セクション8.9」.数値解析入門（第2版）. John Wiley and Sons. ISBN 978-0-471-50023-0。
• アヴリエル、モーデカイ（2003年）『非線形計画法：分析と手法』ドーバー出版、ISBN 978-0-486-43227-4。
• ゴルブ, ジーン・H.; ヴァン・ローン, チャールズ・F. (2013). 「第11章」.行列計算（第4版）. ジョンズ・ホプキンス大学出版局. ISBN 978-1-4214-0794-4。
• サード、ユセフ (2003-04-01). 「第6章」 .疎線形システムのための反復法（第2版）. SIAM. ISBN 978-0-89871-534-7。
• ジェラール・ムラン：「共役勾配アルゴリズムにおけるサイレントエラーの検出と修正」、数値アルゴリズム、第92巻（2023年）、pp.869-891。url=https://doi.org/10.1007/s11075-022-01380-1
• Meurant, Gerard; Tichy, Petr (2024).共役勾配法における誤差ノルム推定. SIAM. ISBN 978-1-61197-785-1。

外部リンク

• 「共役勾配法」数学百科事典、EMS Press、2001 [1994]

「https://en.wikipedia.org/w/index.php?title=共役勾配法&oldid=1318675934」より取得