共役勾配法 （きょうやくこうばいほう、英 : conjugate gradient method 対称 正定値 行列を係数とする連立一次方程式 を解くためのアルゴリズム である[1] [2] [3] [4] 反復法 として利用され[1] [2] [3] [4] コレスキー分解 のような直接法では大きすぎて取り扱えない、大規模な疎行列 を解くために利用される。そのような問題は偏微分方程式 などを数値的に解く際に常に現れる[1] [5] [6] [7]
線型方程式の二次形式を最小化するための、最適なステップサイズによる最急降下法 （緑）の収束と共役勾配法（赤）の収束の比較。共役勾配法は、厳密にはn 次の係数行列に対して高々n ステップで収束する（ここではn =2）。 共役勾配法は、エネルギー最小化などの最適化問題 を解くために用いることもできる[8] [9] [10]
双共役勾配法 （英語版 ） [11] 非線形共役勾配法 が提案されている[12] [13] [14]
対称正定値行列A を係数とするn 元連立一次方程式
の解をx * とする。
非零ベクトルu 、v が
u
T
A
v
=
0
{\displaystyle \mathbf {u} ^{\mathrm {T} }\mathbf {A} \mathbf {v} =\mathbf {0} }
を満たすとき、u 、v はA に関して共役であるという[2] [3] [4] A は対称正定値なので、左辺から内積
⟨
u
,
v
⟩
A
:=
⟨
A
T
u
,
v
⟩
=
⟨
A
u
,
v
⟩
=
⟨
u
,
A
v
⟩
=
u
T
A
v
{\displaystyle \langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }:=\langle \mathbf {A} ^{\mathrm {T} }\mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {A} \mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle =\mathbf {u} ^{\mathrm {T} }\mathbf {A} \mathbf {v} }
を定義することができる。この内積に関して2つのベクトルが直交するなら、それらのベクトルは互いに共役である。
この関係は対称で、u がv に対して共役なら、v もu に対して共役である（この場合の「共役」は複素共役 と無関係であることに注意）。
{p k n 個の互いに共役なベクトル列とする。p k 基底 R n Ax = b の解x *
x
∗
=
∑
i
=
1
n
α
i
p
i
{\displaystyle \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}}
と書ける。ただし係数は
A
x
∗
=
∑
i
=
1
n
α
i
A
p
i
=
b
.
{\displaystyle \mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}=\mathbf {b} .}
p
k
T
A
x
∗
=
∑
i
=
1
n
α
i
p
k
T
A
p
i
=
p
k
T
b
.
{\displaystyle \mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {x} _{*}=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{i}=\mathbf {p} _{k}^{\mathrm {T} }\mathbf {b} .}
α
k
=
p
k
T
b
p
k
T
A
p
k
=
⟨
p
k
,
b
⟩
⟨
p
k
,
p
k
⟩
A
=
⟨
p
k
,
b
⟩
‖
p
k
‖
A
2
.
{\displaystyle \alpha _{k}={\frac {\mathbf {p} _{k}^{\mathrm {T} }\mathbf {b} }{\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{k}}}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\,\,\,\langle \mathbf {p} _{k},\mathbf {p} _{k}\rangle _{\mathbf {A} }}}={\frac {\langle \mathbf {p} _{k},\mathbf {b} \rangle }{\,\,\,\|\mathbf {p} _{k}\|_{\mathbf {A} }^{2}}}.}
で与えられる。
この結果は、上で定義した内積を考えるのが最も分かりやすいと思われる。
以上から、Ax = b を解くための方法が得られる。すなわち、まずn 個の共役な方向を見つけ、それから係数αk
共役なベクトル列p k x * の良い近似を得られる可能性がある。そこで、共役勾配法を反復法として利用することを考える[2] [3] [4] n が非常に大きく、直接法では解くのに時間がかかりすぎるような問題にも適用することができる。
x * の初期値をx 0 = 0 とする。x * が二次形式
f
(
x
)
=
1
2
x
T
A
x
−
b
T
x
,
x
∈
R
n
.
{\displaystyle f(\mathbf {x} )={\frac {1}{2}}\mathbf {x} ^{\mathrm {T} }\mathbf {A} \mathbf {x} -\mathbf {b} ^{\mathrm {T} }\mathbf {x} ,\quad \mathbf {x} \in \mathbf {R} ^{n}.}
を最小化する一意な解であることに注意し、最初の基底ベクトルp 1 をx = x 0 でのf の勾配Ax 0 − b =− b共役勾配法 と呼ぶ[2] [3] [4]
r k k ステップ目での残差
r
k
=
b
−
A
x
k
{\displaystyle \mathbf {r} _{k}=\mathbf {b} -\mathbf {Ax} _{k}}
とする。r k x = x k f の負の勾配であることに注意されたい。最急降下法 はr k p k r k
p
k
+
1
=
r
k
+
1
−
p
k
T
A
r
k
+
1
p
k
T
A
p
k
p
k
{\displaystyle \mathbf {p} _{k+1}=\mathbf {r} _{k+1}-{\frac {\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {r} _{k+1}}{\mathbf {p} _{k}^{\mathrm {T} }\mathbf {A} \mathbf {p} _{k}}}\mathbf {p} _{k}}
のように表すことができる（記事冒頭の図を参照）。
Gnu Octave で書くと以下のようになる。
function [x] = conjgrad(A,b,x0)
r = b - A*x0;
w = -r;
z = A*w;
a = (r'*w)/(w'*z);
x = x0 + a*w;
B = 0;
for i = 1:size(A)(1);
r = r - a*z;
if( norm(r) < 1e-10 )
break;
endif
B = (r'*z)/(w'*z);
w = -r + B*w;
z = A*w;
a = (r'*w)/(w'*z);
x = x + a*w;
end
end
前処理行列 とは、A と同値なP -1 A (P T )-1 の条件数 がA より小さく、Ax =b よりP -1 A (P T )-1 x ′ =P -1 b ′ の方が容易に解けるような正定値行列 P .P T を指す[4] ヤコビ法 、ガウス・ザイデル法 、対称SOR法 などが用いられる[15] [16]
最も単純な前処理行列は、A の対角要素のみからなる対角行列である。これはヤコビ前処理または対角スケーリングとして知られている。対角行列は逆行列の計算が容易かつメモリも消費しない点で、入門用として優れた方法である。より洗練された方法では、κ (A )の減少による収束の高速化とP -1 の計算に要する時間とのトレードオフを考えることになる。
任意の実行列A に対してA T A は対称（半）正定値となるので、係数行列をA T A 、右辺をA T b とする正規方程式 を解くことにより、共役勾配法を任意のn ×m 行列に対して適用することができる（CGNR法[17]
反復法としては、A T A を明示的に保持する必要がなく、行列ベクトル積、転置行列ベクトル積を計算すればよいので、A が疎行列である場合にはCGNR法は特に有効である。ただし、条件数κ(A T A )がκ(A 2 )に等しいことから収束は遅くなる傾向があり、前処理行列 を使用するCGLS (Conjugate Gradient Least Squares[18] A が悪条件である場合に最も数値的に安定な解法である[19] [20]
皆本晃弥. (2005). UNIX & Informatioin Science-5 C 言語による数値計算入門.
田端正久; 偏微分方程式の数値解析, 2010. 岩波書店.
登坂宣好, & 大西和榮. (2003). 偏微分方程式の数値シミュレーション. 東京大学出版会.
Zworski, M. (2002). Numerical linear algebra and solvability of partial differential equations. Communications in mathematical physics, 229(2), 293-307.
Gill, P. E., Murray, W., & Wright, M. H. (1991). Numerical linear algebra and optimization (Vol. 1, p. 74). Redwood City, CA: Addison-Wesley.
Gilbert, J. C., & Nocedal, J. (1992). Global convergence properties of conjugate gradient methods for optimization. SIAM Journal on optimization, 2(1), 21-42.
Steihaug, T. (1983). The conjugate gradient method and trust regions in large scale optimization. SIAM Journal on Numerical Analysis, 20(3), 626-637.
Dai, Y. H. (2010). Nonlinear conjugate gradient methods. Wiley Encyclopedia of Operations Research and Management Science.
Hager, W. W., & Zhang, H. (2006). A survey of nonlinear conjugate gradient methods. Pacific journal of Optimization, 2(1), 35-58.
Dai, Y., Han, J., Liu, G., Sun, D., Yin, H., & Yuan, Y. X. (2000). Convergence properties of nonlinear conjugate gradient methods. SIAM Journal on Optimization, 10(2), 345-358.
Eisenstat, S. C. (1981). Efficient implementation of a class of preconditioned conjugate gradient methods. SIAM Journal on Scientific and Statistical Computing, 2(1), 1-4.
Kaasschieter, E. F. (1988). Preconditioned conjugate gradients for solving singular systems. Journal of Computational and Applied Mathematics, 24(1-2), 265-275.
Bjorck, A. (1996). Numerical methods for least squares problems (Vol. 51). SIAM.
Paige, C. and Saunders, M. "LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares." ACM Trans. Math. Soft. 8, 43-71, 1982.
Paige, C. C., & Saunders, M. A. (1982). Algorithm 583: LSQR: Sparse linear equations and least squares problems. ACM Transactions on Mathematical Software (TOMS), 8(2), 195-209.
