直接法としての導出 共役勾配法は、最適化のための共役方向法の特殊化や、固有値 問題に対するアーノルディ ／ランチョス 反復法の変形など、いくつかの異なる観点から導出できる。これらの導出方法はアプローチこそ異なるものの、共通のテーマは、残差の直交性と探索方向の共役性を証明することである。これら2つの性質は、この方法のよく知られた簡潔な定式化を導き出す上で極めて重要である。
2 つの非ゼロベクトルはu {\displaystyle \mathbf {u} } v {\displaystyle \mathbf {v} } A {\displaystyle \mathbf {A} }
u T A v = 0. {\displaystyle \mathbf {u} ^{\mathsf {T}}\mathbf {A} \mathbf {v} =0.} 以来A {\displaystyle \mathbf {A} } 内積を定義する。
u T A v = ⟨ u 、 v ⟩ A := ⟨ A u 、 v ⟩ = ⟨ u 、 A T v ⟩ = ⟨ u 、 A v ⟩ 。 {\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 つのベクトルが共役であるのは、それらがこの内積に関して直交する場合に限る。共役であることは対称的な関係である。u {\displaystyle \mathbf {u} } v {\displaystyle \mathbf {v} } v {\displaystyle \mathbf {v} } u {\displaystyle \mathbf {u} }
P = { p 1 、 … 、 p n } {\displaystyle P=\{\mathbf {p} _{1},\dots ,\mathbf {p} _{n}\}} は、n {\displaystyle n} A {\displaystyle \mathbf {A} } p 私 T A p j = 0 {\displaystyle \mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0} 私 ≠ j {\displaystyle i\neq j} P {\displaystyle P} 基礎 を形成するR n {\displaystyle \mathbb {R} ^{n}} x * {\displaystyle \mathbf {x} _{*}} A x = b {\displaystyle \mathbf {Ax} =\mathbf {b} }
x * = ∑ 私 = 1 n α 私 p 私 ⇒ A x * = ∑ 私 = 1 n α 私 A p 私 。 {\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}.} 問題を左から乗算するA x = b {\displaystyle \mathbf {Ax} =\mathbf {b} } p k T {\displaystyle \mathbf {p} _{k}^{\mathsf {T}}}
p k T b = p k T A x * = ∑ 私 = 1 n α 私 p k T A p 私 = ∑ 私 = 1 n α 私 ⟨ p k 、 p 私 ⟩ A = α k ⟨ p k 、 p k ⟩ A {\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} }} など
α k = ⟨ p k 、 b ⟩ ⟨ p k 、 p k ⟩ A 。 {\displaystyle \alpha _{k}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} }}}.} これにより、方程式を解くための以下の方法[ 4 ]が得られる。 A x = b {\displaystyle \mathbf {Ax} =\mathbf {b} } n {\displaystyle n} α k {\displaystyle \alpha _{k}}
反復法として 共役ベクトルを選択するとp k {\displaystyle \mathbf {p} _{k}} x * {\displaystyle \mathbf {x} _{*}} n {\displaystyle n}
初期推定値を と表記します。x * {\displaystyle \mathbf {x} _{*}} x 0 {\displaystyle \mathbf {x} _{0}} x 0 = 0 {\displaystyle \mathbf {x} _{0}=\mathbf {0} } A z = b − A x 0 {\displaystyle \mathbf {Az} =\mathbf {b} -\mathbf {Ax} _{0}} x 0 {\displaystyle \mathbf {x} _{0}} x * {\displaystyle \mathbf {x} _{*}} x * {\displaystyle \mathbf {x} _{*}} 二次関数の唯一の最小値でもある。
f （ x ） = 1 2 x T A x − x T b 、 x ∈ R n 。 {\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階微分行列 が対称正定値であることから明らかである。
H （ f （ x ） ） = A 、 {\displaystyle \mathbf {H} (f(\mathbf {x} ))=\mathbf {A} \,,} ミニマイザー（使用）D f （ x ） = 0 {\displaystyle Df(\mathbf {x} )=0}
∇ f （ x ） = A x − b 。 {\displaystyle \nabla f(\mathbf {x} )=\mathbf {A} \mathbf {x} -\mathbf {b} \,.} これは、最初の基底ベクトルを取ることを示唆しています。p 0 {\displaystyle \mathbf {p} _{0}} f {\displaystyle f} x = x 0 {\displaystyle \mathbf {x} =\mathbf {x} _{0}} f {\displaystyle f} A x − b {\displaystyle \mathbf {Ax} -\mathbf {b} } x 0 {\displaystyle \mathbf {x} _{0}} p 0 = b − A x 0 {\displaystyle \mathbf {p} _{0}=\mathbf {b} -\mathbf {Ax} _{0}} 共役勾配法 と呼ばれます。p 0 {\displaystyle \mathbf {p} _{0}} 残差 でもある。
させてr k {\displaystyle \mathbf {r} _{k}} は k {\displaystyle k}
r k = b − A x k 。 {\displaystyle \mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}.} 上記のように、r k {\displaystyle \mathbf {r} _{k}} f {\displaystyle f} x k {\displaystyle \mathbf {x} _{k}} 勾配降下 法では方向r k p k {\displaystyle \mathbf {p} _{k}} グラム・シュミット正規直交化 の例として見なすことができます。これにより、次の式が得られます。
p k = r k − ∑ 私 < k r k T A p 私 p 私 T A p 私 p 私 {\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}} （共役制約が収束に及ぼす影響については、記事冒頭の図を参照してください。）この方向に従うと、次の最適位置は次のように与えられます。
x k + 1 = x k + α k p k {\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}} と
α k = p k T （ b − A x k ） p k T A p k = p k T r k p k T A 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}}},} ここで最後の等式は、の定義から導かれる。r k {\displaystyle \mathbf {r} _{k}} α k {\displaystyle \alpha _{k}} x k +1f に代入し、それを最小化すれば、を導出できる。α k {\displaystyle \alpha _{k}}
f （ x k + 1 ） = f （ x k + α k p k ） =: g （ α k ） g ′ （ α k ） = ! 0 ⇒ α k = p k T （ b − A x k ） p k T A p 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}}}
結果として得られるアルゴリズム 上記のアルゴリズムは、共役勾配法の最も分かりやすい説明です。一見すると、このアルゴリズムは、以前のすべての探索方向と残差ベクトル、および多くの行列ベクトル乗算を保存する必要があり、計算コストが高くなる可能性があります。しかし、アルゴリズムのより詳細な分析[ 6 ] ： r 私 {\displaystyle \mathbf {r} _{i}} r j {\displaystyle \mathbf {r} _{j}} r 私 T r j = 0 {\displaystyle \mathbf {r} _{i}^{\mathsf {T}}\mathbf {r} _{j}=0} 私 ≠ j {\displaystyle i\neq j} p 私 {\displaystyle \mathbf {p} _{i}} A {\displaystyle \mathbf {A} } p j {\displaystyle \mathbf {p} _{j}} p 私 T A p j = 0 {\displaystyle \mathbf {p} _{i}^{\mathsf {T}}\mathbf {A} \mathbf {p} _{j}=0} 私 ≠ j {\displaystyle i\neq j} p 私 {\displaystyle \mathbf {p} _{i}} r 私 {\displaystyle \mathbf {r} _{i}} クリロフ部分空間 を張る、r 私 {\displaystyle \mathbf {r} _{i}} p 私 {\displaystyle \mathbf {p} _{i}} A {\displaystyle \mathbf {A} } x k {\displaystyle \mathbf {x} _{k}} x {\displaystyle \mathbf {x} }
つまり、CG法がx 0 = 0 {\displaystyle \mathbf {x} _{0}=0} [ 7 ] x k = 1 r g m 私 n y ∈ R n { （ x * − y ） ⊤ A （ x * − y ） : y ∈ スパン { b 、 A b 、 … 、 A k − 1 b } } {\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\}}} x * {\displaystyle x_{*}} A x = b {\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} }
解決のためのアルゴリズムは以下に詳述する。A x = b {\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} } A {\displaystyle \mathbf {A} } x 0 {\displaystyle \mathbf {x} _{0}} 0 {\displaystyle \mathbf {0} }
r 0 := b − A x 0 もし r 0 が十分に小さい場合は、 x 0 その結果 p 0 := r 0 k := 0 繰り返す α k := r k T r k p k T A p k x k + 1 := x k + α k p k r k + 1 := r k − α k A p k もし r k + 1 が十分に小さい場合、ループを終了します。 β k := r k + 1 T r k + 1 r k T r k p k + 1 := r k + 1 + β k p k k := k + 1 繰り返し終了 戻る x k + 1 その結果 {\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 {\displaystyle \beta _{k}} 非線形共役勾配法 でも使用されます。
線形代数 を使用する """ x = conjugate_gradient(A, b, x0 = zero(b); atol=length(b)*eps(norm(b)) 共役勾配法を用いて、`A * x = b` の解を返してください。 `A`は正定値行列またはその他の線形演算子でなければなりません。 `x0`は解の初期推定値です（デフォルトはゼロベクトルです）。 `atol`は残差`b - A * x`の大きさに対する絶対許容値です。 収束のため（デフォルトはマシンイプシロン）。 近似解ベクトル`x`を返します。 """ 関数 conjugate_gradient ( A 、 b :: AbstractVector 、 x0 :: AbstractVector = zero ( b ); atol = length ( b ) * eps ( norm ( b )) ） x = copy ( x0 ) # 解を初期化する r = b - A * x0 # 初期残差 p = copy ( r ) # 初期探索方向 r²old = r ' * r # 残差の二乗ノルム 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 を返す 終わり function = conjugate_gradient ( A, b, x0, tol ) 共役勾配法を用いて、`A * x = b` の解を返します。 % 注意: A は対称かつ正定値である必要があります。 ナルギン < 4 の場合 tol = eps ; 終わり r = b - A * x0 ; p = r ; rsold = r ' * r ; x = x0 ; while sqrt ( rsold ) > tol Ap = A * p ; alpha = rsold / ( p ' * Ap ); x = x + alpha * p ; r = r - alpha * Ap ; rsnew = r ' * r ; p = r + ( rsnew / rsold ) * p ; rsold = rsnew ; 終わり 終わり
数値例 線形システムAx = b を考えます。
A x = [ 4 1 1 3 ] [ x 1 x 2 ] = [ 1 2 ] 、 {\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つのステップを実行します。
x 0 = [ 2 1 ] {\displaystyle \mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}} システムの近似解を見つけるため。
解決 参考までに、正確な解は次のとおりです。
x = [ 1 11 7 11 ] ≈ [ 0.0909 0.6364 ] {\displaystyle \mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\\\{\frac {7}{11}}\end{bmatrix}}\approx {\begin{bmatrix}0.0909\\\\0.6364\end{bmatrix}}} 最初のステップは、x 0 に関連付けられた残差ベクトルr 0 を計算することです。この残差は、式r 0 = b - Ax 0 から計算され、この場合、次のようになります。
r 0 = [ 1 2 ] − [ 4 1 1 3 ] [ 2 1 ] = [ − 8 − 3 ] = p 0 。 {\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 の 選択方法は、以降の反復で変更されます。
次に、関係式を用いてスカラーα₀を計算します 。
α 0 = r 0 T r 0 p 0 T A p 0 = [ − 8 − 3 ] [ − 8 − 3 ] [ − 8 − 3 ] [ 4 1 1 3 ] [ − 8 − 3 ] = 73 331 ≈ 0.2205 {\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 を計算できます。
x 1 = x 0 + α 0 p 0 = [ 2 1 ] + 73 331 [ − 8 − 3 ] ≈ [ 0.2356 0.3384 ] 。 {\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 1 が得られます。次に、式を使用して次の残差ベクトルr 1を計算します。
r 1 = r 0 − α 0 A p 0 = [ − 8 − 3 ] − 73 331 [ 4 1 1 3 ] [ − 8 − 3 ] ≈ [ − 0.2810 0.7492 ] 。 {\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
β 0 = r 1 T r 1 r 0 T r 0 ≈ [ − 0.2810 0.7492 ] [ − 0.2810 0.7492 ] [ − 8 − 3 ] [ − 8 − 3 ] = 0.0088。 {\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₁ を計算できます。
p 1 = r 1 + β 0 p 0 ≈ [ − 0.2810 0.7492 ] + 0.0088 [ − 8 − 3 ] = [ − 0.3511 0.7229 ] 。 {\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}}.} α 0 p 1 を使用してスカラーα 1
α 1 = r 1 T r 1 p 1 T A p 1 ≈ [ − 0.2810 0.7492 ] [ − 0.2810 0.7492 ] [ − 0.3511 0.7229 ] [ 4 1 1 3 ] [ − 0.3511 0.7229 ] = 0.4122。 {\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 2 を求めます。
x 2 = x 1 + α 1 p 1 ≈ [ 0.2356 0.3384 ] + 0.4122 [ − 0.3511 0.7229 ] = [ 0.0909 0.6364 ] 。 {\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 1 およびx 0 よりもシステムの解に対する「より良い」近似値となります。この例で限定精度ではなく厳密な算術演算を使用した場合、理論的にはn = 2 回（n はシステムの次数）の反復後に厳密解に到達したはずです。
有限終端特性 厳密な計算を用いる場合、必要な反復回数は行列の次数以下になります。この性質は、共役勾配法の有限終了性 として知られています。これは、厳密な計算を用いる場合、線形システムの厳密解に有限のステップ数（最大でもシステムの次元数）で到達できることを意味します。この性質は、各反復において、この方法がそれまでのすべての残差と直交する残差ベクトルを生成するという事実から生じます。これらの残差は互いに直交する集合を形成します。
n 次元空間では、 n 個を超える線形独立かつ互いに直交するベクトルを構成することは、そのうちの1つが零ベクトルでない限り不可能です。したがって、残差がゼロになった時点で、この方法は解に到達したとみなされ、終了しなければなりません。これにより、共役勾配法は最大でもn ステップで収束することが保証されます。
これを実証するために、次のシステムを考えてみましょう。
A = [ 3 − 2 − 2 4 ] 、 b = [ 1 1 ] {\displaystyle A={\begin{bmatrix}3&-2\\-2&4\end{bmatrix}},\quad \mathbf {b} ={\begin{bmatrix}1\\1\end{bmatrix}}}
まず最初の推測から始めますx 0 = [ 1 2 ] {\displaystyle \mathbf {x} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}} A {\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 ; end disp ( '正確な解:' ); disp ( x ); 出力は、メソッドが到達したことを確認します[ 1 1 ] {\displaystyle {\begin{bmatrix}1\\1\end{bmatrix}}}
柔軟な前処理付き共役勾配法 数値的に困難なアプリケーションでは、高度な前処理器が使用され、反復ごとに変化する可変前処理につながる可能性があります。前処理器がすべての反復で対称正定値であっても、変化する可能性があるという事実は上記の議論を無効にし、実際のテストでは、上記のアルゴリズムの収束を著しく遅くします。Polak –Ribièreの 公式を使用すると、
β k := r k + 1 T （ z k + 1 − z k ） r k T z k {\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}}}} フレッチャー・リーブスの 公式の代わりに
β k := r k + 1 T z k + 1 r k T z k {\displaystyle \beta _{k}:={\frac {\mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k+1}}{\mathbf {r} _{k}^{\mathsf {T}}\mathbf {z} _{k}}}} この場合、収束が劇的に改善される可能性があります。[ 15 ] 。 また、 前 処理行列が対称正定値（SPD）でない場合でも、柔軟であることが示されています[ 17 ]
柔軟なバージョンの実装には、追加のベクトルを保存する必要があります。固定SPD前処理器の場合、r k + 1 T z k = 0 、 {\displaystyle \mathbf {r} _{k+1}^{\mathsf {T}}\mathbf {z} _{k}=0,} したがって、 β k の2つの式は、丸め誤差 なしで、正確な算術において等価です。
Polak–Ribière の 公式を用いた方法の収束挙動が優れていることの数学的な説明は、この場合、その方法が局所的に最適で あり、特に、局所的に最適な最急降下法よりも収束が遅くならないということである。[ 18 ]
局所的に最適な最急降下法との比較 元の共役勾配法と前処理付き共役勾配法の両方で、設定する必要があるのは β k := 0 {\displaystyle \beta _{k}:=0} 線探索 法や最急降下 法を用いる。この置換により、ベクトルp は常にベクトル z p 最急降下 法の各反復は、共役勾配法に比べて若干コストが低くなる。ただし、後者は、（非常に）可変な、または非 SPD な前処理器 を使用しない限り、より速く収束する（上記参照）。
正規方程式に対する共役勾配法 共役勾配法は、任意のn × m行列に適用できます。これは、 正規方程式 A T A と右辺ベクトルA T b に適用することで実現できます。A T A は任意の A に対して対称正定値半行列であるため です 。 結果 として 得られる のは、正規方程式に対する共役勾配 ( CGN またはCGNR ) です。
A T Ax = A T b 反復法であるため、 A T A をメモリ上に明示的に生成する必要はなく 、行列とベクトルの乗算および転置行列とベクトルの乗算を実行するだけで済みます。したがって、これらの演算は通常非常に効率的であるため、 Aが 疎行列で ある場合に CGNR は特に有用です。ただし、正規方程式を生成する際の欠点は、条件数 κ( A T A ) が κ 2 ( A )と等しくなるため、CGNR の収束速度が遅くなり、近似解の精度が丸め誤差に敏感になる可能性があることです。適切な前処理行列 を見つけることは、CGNR 法を使用する上で重要な部分となることがよくあります。
いくつかのアルゴリズムが提案されている（例：CGLS、LSQR）。LSQRアルゴリズムは、 Aが 条件数が 大きい、つまり条件が悪い 場合に、最も優れた数値安定性を持つと言われている。
複素エルミート行列に対する共役勾配法 共役勾配法にわずかな修正を加えることで、複素数値行列 A とベクトル b が与えられた場合に、線形方程式系を解くことができる。A x = b {\displaystyle \mathbf {A} \mathbf {x} =\mathbf {b} } エルミート行列 (つまり A' = A) であり、正定値行列 であり、記号 ' は共役転置 を表します。簡単な変更は、実転置を共役 転置 に置き換えるだけです。
メリットとデメリット 共役勾配法の利点と欠点については、NemirovskyとBenTalの講義ノート[ 20 ] の
参考文献 ↑ Hestenes, Magnus R. ; Stiefel, Eduard (1952 年 12 月). "線形システムを解くための共役勾配法" (PDF) . Journal of Research of the National Bureau of Standards . 49 (6): 409. doi : 10.6028/jres.049.044 .↑ Straeter, TA (1971). On the Extension of the Davidon–Broyden Class of Rank One, Quasi-Newton Minimization Methods to an Infinite Dimensional Hilbert Space with Applications to Optimal Control Problems (PhD thesis). North Carolina State University. hdl : 2060/19710026200 – via NASA Technical Reports Server. ↑ シュパイザー、アンブロス (2004)。 「Konrad Zuse und die ERMETH: Ein weltweiter Architektur-Vergleich」 [ Konrad Zuse と ERMETH: 建築の世界的比較 ] 。 『Hellige』、ハンス・ディーター編著。 情報収集。 Visionen、Paradigmen、Leitmotive (ドイツ語)。ベルリン：シュプリンガー。 p. 185.ISBN 3-540-00217-0 。1 2 3 4 ポリャク、ボリス (1987)。 最適化入門 。 1 2 3 Greenbaum, Anne (1997). Iterative Methods for Solving Linear Systems . doi : 10.1137/1.9781611970937 . ISBN 978-0-89871-396-1 。↑ Botev, Zdravko I.; Kroese, Dirk P.; Taimre, Thomas (2025). Data Science and Machine Learning: Mathematical and Statistical Methods (2nd ed.). Boca Raton ; London: CRC Press. pp. 558–559 . ISBN 978-1-032-48868-4 。↑ Paquette, Elliot; Trogdon, Thomas (2023 年 3 月). "標本共分散行列に対する共役勾配法と MINRES アルゴリズムの普遍性" . Communications on Pure and Applied Mathematics . 76 (5): 1085– 1136. arXiv : 2007.00640 . doi : 10.1002/cpa.22081 . ISSN 0010-3640 . ↑ Shewchuk, Jonathan R (1994). An Introduction to the Conjugate Gradient Method Without the Agonizing Pain (PDF) . ↑ Saad, Yousef (2003). Iterative methods for sparse linear systems (2nd ed.). Philadelphia, Pa.: Society for Industrial and Applied Mathematics. pp . 195. ISBN 978-0-89871-534-7 。↑ Holmes, M. (2023). Introduction to Scientific Computing and Data Analysis, 2nd Ed . Springer. ISBN 978-3-031-22429-4 。↑ Hackbusch, W. (2016-06-21). Iterative solution of large sparse systems of equations (2nd ed.). Switzerland: Springer. ISBN 978-3-319-28483-5 OCLC 952572240 ↑ Barrett, Richard; Berry, Michael; Chan, Tony F.; Demmel, James; Donato, June; Dongarra, Jack; Eijkhout, Victor; Pozo, Roldan; Romine, Charles; van der Vorst, Henk. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods (PDF) (2nd ed.). Philadelphia, PA: SIAM. p. 13 . 2020-03-31 に取得 。 ↑ Golub, Gene H.; Van Loan, Charles F. (2013). 行列計算 （第4 版）. ジョンズ・ホプキンス大学出版局. セクション 11.5.2. ISBN 978-1-4214-0794-4 。↑ Concus, P.; Golub, GH; Meurant, G. (1985). "共役勾配法のためのブロック前処理" . SIAM Journal on Scientific and Statistical Computing . 6 (1): 220– 252. doi : 10.1137/0906018 . ↑ 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 . ↑ Notay, Yvan (2000). "Flexible Conjugate Gradients". SIAM Journal on Scientific Computing . 22 (4): 1444– 1460. CiteSeerX 10.1.1.35.7473 . doi : 10.1137/S1064827599362314 . ↑ 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 . ↑ 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 . 1 2 Ross, IM 、「加速最適化のための最適制御理論」、 arXiv : 1902.09004 、2019年。↑ Nemirovsky と Ben-Tal (2023). "最適化 III: 凸最適化" (PDF) 。 ↑ Pennington, Fabian Pedregosa, Courtney Paquette, Tom Trogdon, Jeffrey. "ランダム行列理論と機械学習チュートリアル" . random-matrix-learning.github.io . 2023-12-05 に取得. {{cite web}}: CS1 maint: 複数の名前: 著者リスト (リンク)
さらに読む アトキンソン、ケンドール A. (1988). 「第 8.9 節」.数値解析入門 版). ジョン・ワイリー・アンド・サンズ. ISBN 978-0-471-50023-0 。 アヴリエル、モルデカイ（2003）。非線形計画法：分析と手法 。ドーバー出版。ISBN 978-0-486-43227-4 。 Golub, Gene H.; Van Loan, Charles F. (2013). 「第 11 章」.行列計算 (第 4 版). ジョンズ・ホプキンス大学出版局. ISBN 978-1-4214-0794-4 。 Saad, Yousef (2003-04-01). 「第6章」 .スパース線形システムのための反復法 版）. SIAM. ISBN 978-0-89871-534-7 。 Gérard Meurant: 「共役勾配法におけるサイレントエラーの検出と修正」、Numerical Algorithms、第92巻（2023年）、869-891ページ。url= https://doi.org/10.1007/s11075-022-01380-1 Meurant, Gerard; Tichy, Petr (2024).共役勾配法における誤差ノルム推定 . SIAM. ISBN 978-1-61197-785-1 。 
