115

第 10 章 連立一次方程式の解法 3 — Krylov 部
分空間法
共役勾配法は 1952 年に，M.R.Hestenes と E.Stiefelによって発表された。
それは，実に華やかなデビューだった。新聞やラジオは「科学史上に残る大発明」として，ビッグ・ニュースの形で報道した。当時，筆者は高校生であったから，詳しい内容は全く分からなかったが連立方程式を解くための画期的な方法が発明された。これまで何百年もの間，n 元の連立方程式を解くには n3 に比例する演算が必要なもの，とされてきたが，新しい方法によれば，これを n2 に比例する演算回数で解くことができる・・・。という記事を (中略) はっきり記憶している。数値計算法の一解法が，これほど大々的に報道されたことは，前にも後にも例がない。
戸川隼人「共役勾配法」（教育出版）

傾斜法は反復法の一種であるが，前章の反復法とは原理的に著しく異なる点がある。そこで新たに１章を起こした。
近年よく取り上げられる共役勾配法とその変種は，普通，傾斜法の一種として導出される。ここでは傾斜法の概観を述べた後，その代表的なアルゴリズムである最急降下法と共役勾配法について述べる。
内容

1. 傾斜法 (Gradient Method) —最適化問題のための傾斜法を連立一次方程式に適用できる条件と，その原理を述べる。

2. 最急降下法 (Steepest Desent Method) —最急降下法のアルゴリズムを述べる。

3. 共役勾配法 (Conjugate-Gradient Method) —共役勾配法のアルゴリズムを述べる。

第 10 章 連立一次方程式の解法 3 — Krylov 部分空間法

116

10.1 傾斜法
行列 A ∈ Mn (R) が正定値対称，即ち任意の Rn ベクトル x , 0 に対し(x, Ax) = (Ax, x) ≥ 0であるとき，連立一次方程式

Ax = b(10.1)

は 2 次形式の極値問題

f (y) ===1(y − x, A(y − x))21(y − x, Ay − b)211(y, Ay) − (y, b) + (x, b)
22(10.2)の解でもある。このとき，極値問題 (10.2) を解くための次のアルゴリズムを，傾斜法という。
アルゴリズム 16 (傾斜法 (Gradient Method))

1. 初期値 x0 ∈ Rn を決める。
2. r0 := b − Ax0 ，p0 := r0 とする。
3. k = 0, 1, 2, · · · に対して以下を計算する。
(rk , pk )(pk , Apk )(b) xk+1 := xk + αk pk(a) αk :=(c) rk+1 := rk − αk Apk (又は := b − Axk+1 )(d) 収束判定(e) pk+1 を決める
ここで rk ∈ Rn を残差ベクトル (residual vector)，pk ∈ Rn を修正ベクトル (modification vector) という。
残差ベクトル rk は

rk = b − Axkという量であるが，この形だと反復一回毎に Axk を計算せねばならず，不経済である。従って普通はrk=b − Axk=b − A(xk + αk pk )=(b − Axk ) − αk Apk=

rk − αk Apkという形で再帰的に計算される。
この残差ベクトルの幾何学的意味は次の補題で示される。

10.1. 傾斜法117

補題 10.1.1

 ∂f
 ∂x (y)  ∂ f1
∂x2 (y)  = −(b − Ay) = −r
∇ f (y) = ..
.

∂f∂xn (y)(証)

f (y)

1
1(y, Ay) − (y, b) + (x, b)
22∑
1∑
1ai j yi y j −yi bi + (x, b)
2 i, j
2i

==から

∂f(y)
∂yl

==

1∑
1∑
al j y j +a jl y j − bl
2 j
2 j
∑
al j y j − bljより∇ f (y) = −(b − Ax) = −rが示された。

(証明終)

この補題より，xk に対する残差 rk は，点 xk における最も勾配のきつい，谷 ( f (y) = 0) 方向への傾斜を示していることがわかる。
また αk はつぎのような性質を持つ。
補題 10.1.2∂f(xk+1 ) = 0∂αk(証)f (xk+1 )

==
=
=

から

1(xk+1 , Axk+1 ) − (xk+1 , b) + (x, b)21(xk + αk pk , Axk + αk Apk ) − (xk + αk pk , b) + (x, b)2
1
1
1
1(xk , Axk ) + (xk , αk Apk ) + (αk pk , Axk ) + (αk pk , αk Apk ) − (xk , b) − (αk pk , b) + (x, b)
22
2
2
1 2
1(xk , Axk ) + αk (xk , Apk ) + αk (pk , Apk ) − (xk , b) − αk (pk , b) + (x, b)
22

第 10 章 連立一次方程式の解法 3 — Krylov 部分空間法

118∂f(xk+1 )∂αk=(xk , Apk ) + αk (pk , Apk ) − (pk , b)

=
=
=
=
=

(xk , Apk ) − (pk , b) + (rk , pk )(pk , Axk ) − (pk , b) + (rk , pK )(pk , −(b − Axk )) + (rk , pk )−(rk , pk ) + (rk , pk )0
(証明終)

この補題から，αk が，pk という方向において最も極値に近い位置を決めていることがわかる。
傾斜法は，前述のアルゴリズムの修正ベクトル pk の決め方によりさまざまな変種を作ることができる。

10.2 最急降下法 (Steepest Decent)
最急降下法は前節の傾斜法において，修正ベクトル pk をpk := rk(10.3)とする。アルゴリズムは次のようになる。
アルゴリズム 17 (最急降下法)

1. 初期値 x0 ∈ Rn を決める。
2. r0 := b − Ax0 とする。
3. k = 0, 1, 2, · · · に対して以下を計算する。
∥rk ∥22(rk , Ark )(b) xk+1 := xk + αk rk(a) αk :=(c) rk+1 := rk − αk Ark (又は := b − Axk+1 )(d) ∥xk+1 − xk ∥ < ε∥xk ∥ ならば終了。そうでなければ (a) へ。
このアルゴリズムは基本的には，前節の反復法と数値的な性質は類似点が多く，停止則も同じものが使える。

10.3 共役勾配法
共役勾配法 (Conjugate-Gradient Method, CG 法) のアルゴリズムは以下の通りである。
アルゴリズム 18 (共役勾配法 (Conjugate-Gradient Method))

1. 初期値 x0 ∈ Rn を決める。

10.3. 共役勾配法

119

2. r0 := b − Ax0 ，p0 := r0 とする。
3. k = 0, 1, 2, · · · に対して以下を計算する。
(rk , pk )(pk , Apk )(b) xk+1 := xk + αk pk(a) αk :=(c) rk+1 := rk − αk Apk (又は := b − Axk+1 )∥r∥2(d) βk := ∥rk+1∥2 2k 2(e) 収束判定(f) pk+1 := rk+1 + βk pk
このアルゴリズムの特徴は(1) ベクトル計算が多い(2) 反復法に似てはいるが，有限回の操作で収束する(3) 数値的性質がよく把握されていない点にある。(1) は現在流行の並列計算機にとっては都合が良い。1CPU あれば n 回の計算が必要なところを，n CPU ならば 1 回の計算で済む1 。前章の反復法よりも並列化しやすいのである。この点を踏まえて，共役勾配法系の解法をベクトル反復法という名で分類している本もある。

(2) は次の補題から示される。
補題 10.3.1
共役勾配法において，0 ≤ i < j ≤ l ≤ n のとき(1) (r j , pi ) = 0(2) (r j , p j ) = ∥r j ∥22(3) (pi , Ap j ) = 0(4) (ri , r j ) = 0が成立する。
(証)
ここで 0 ≤ i < j ≤ k とする。k ≤ l とし，k に対する帰納法で証明する。まず命題 (Ak ) を(1) (r j , pi ) = 0(2) (r j , p j ) = ∥r j ∥22(3) (pi , Ap j ) = 0(4) (ri , r j ) = 0

1 CPU 間の同期をとらなればいけないから，そう単純ではないが。

第 10 章 連立一次方程式の解法 3 — Krylov 部分空間法

120とする。
(A0 ) は (2) が r0 = p0 から自明。それ以外は命題の条件を満たしていない。
(Ak ) を仮定すると (Ak+1 ) は次のようにして示される。

(1) i < j ≤ k のときは示されているから，まず i = k, j = k + 1 のときを示せばよい。
(rk+1 , pk )===(rk , pk ) − αk (Apk , pk )(rk , pk )(rk , pk ) −(Apk , pk )(pk , Apk )0

さらに (Ak )(1)，(Ak )(3) から(rk+1 , p j ) = (rk , p j ) − αk (Apk , p j ) = 0がわかる。
(2) 上の命題を使って(rk+1 , pk+1 )==(rk+1 , rk+1 ) + βk (rk+1 , pk )(rk+1 , rk+1 )が示された。

(3)(pk+1 , Ap j )

=
=
=
=

(4) (Ak+1 )(1)(3) から(rk+1 , Ap j ) − βk (pk , Ap j )1(rk+1 , r j − r j+1 ) − βk (pk , Ap j )αj1(rk+1 , p j − β j−1 p j−1 − p j+1 + β j p j ) − βk (Apk , pk )αj{0 for j < k(Ak )(3) と (Ak+1 )(1) から。
0 for j = k αk , βk を展開し，(Ak+1 )(1)(2) から。
(r j , rk+1 ) = (r j , rk ) − αk (r j , Apk ) = 0。

よって補題は証明された。

(証明終)

この補題の (4) と第 1 章の基本定理からl≤nでなければならない。従って，共役勾配法は n 回以下の反復で真値が得られる。が，実際の数値計算においてはこの条件が満足されず，前章で示した反復法と同様，収束判定に基づいて反復を制御する必要がある。

10.4 共役勾配法の数値的性質
共役勾配法が有限回の代数的操作によって終了するアルゴリズムであることは既に証明した。
しかし，有限桁の浮動小数点数を用いて計算を行うと，大抵の場合，丸め誤差のため真値を得ることができない。そればかりか，限界精度を得ようとすると n 回以上の反復回数を必要とする。
このような，理論と現実との食い違いは割合早くから知られていた。[43] に Ginsburg が寄稿した内容は，もっと以前に同人によって知らされてきたし，[39] にもその記述がある。
共役勾配法の開発チームの一人である T. Ginsburg は [43] に次の言葉を残した。

10.4. 共役勾配法の数値的性質

121“· · · the CG-method has to be treated as if it were an infinite iterative process, proper termination is not trivial.”
実際いくつかの問題を解いてみると，限界精度に到達するまでの回数はまちまちである。Ginsburg
はその目安として，3n〜5n 回まわせばほぼ事足りるとしている。5n 回でも十分でないものは解を得ることはできないであろうとも言っている。
確かに一般的傾向としてなら間違いではない。しかし，反例はいくらでもある。あくまで一つの目安として考えるべきだろう。
例題 10.4.1

 1 1
 2
 1A =  3 .
 ..
 1n121323123..
.

1..
.

2n3n············



]
 [min(i, j) =max(i, j).. . 
11n2n3nを用いてAT Ax = AT bという方程式を解く。次数は n = 200 する。ここで真の解は

x = [1 1 · · · 1]Tとした。
これを解くと反復回数最大相対誤差1791.368 × 10−31879

1.544 × 10−6となった。つまり 10n 回反復してやっと精度限界までたどり着いたことになる。
この例でもわかるように，問題毎に必要な反復回数は全く異なり，明確な基準を設けることは不可能である。しかも，それは問題の悪条件性にはあまり関係しない [39]。
また，共役勾配法においては停止則の問題は非常に重要である。それは反復法でよく用いられる∥xk+1 − xk ∥ ≤ ε∥xk+1 ∥ になったら停止する。という，近似値を直接監視する手法が使えないことである。それは次のような例がめずらしくないことでわかる。
例題 10.4.2

 nn − 1 ··· n − 1n···A =  .
.
..
 ..
12···

 
 1 1  
 1 2  = [n − |i − j|], x =  ..  .. .  .  1n第 10 章 連立一次方程式の解法 3 — Krylov 部分空間法

122とする連立一次方程式

AT Ax = AT bを解く。ここで次元数は n = 20 とし，IEEE754 倍精度計算を用いた。

ε = 10−4 としたときには，まだ 1〜2 桁しか出ていないのに止まってしまう。
x4

限界精度は

x5

12.98681.10133145741063380E + 00
414303215180E + 01.98681.10133

276664000060E + 00
417987338890E + 01

345.10062.99568.99391

701544218720E + 01
510096151920E + 00
548805348430E + 00.10062.99568.99391

699858318230E + 01
519670928390E + 00
615683767450E + 00

67.99894.10037653557090620E + 00
572541968110E + 01.99894.10037755009384550E + 00
582440072840E + 01

8910.10039.10001.99655

092320427600E + 01
612977915360E + 01
740790353490E + 00.10039.10001.99655

099941670590E + 01
618390650230E + 01
783355162320E + 00

1112.99655.10001740790350030E + 00
612977914350E + 01.99655.10001783355158880E + 00
618390649220E + 01

131415.10039.10037.99894092320427310E + 01
572541967950E + 01
653557082920E + 00.10039.10037.99894

099941670300E + 01
582440072680E + 01
755009376880E + 00

1617.99391.99568548805344020E + 00
510096148570E + 00.99391.99568615683763070E + 00
519670925070E + 00

181920.10062.10133.98681701544218340E + 01
414303216630E + 01
145741066940E + 00.10062.10133.98681

699858317860E + 01
417987340350E + 01
276664003650E + 00

10.4. 共役勾配法の数値的性質

123x21.10000000000178860E + 01

234.99999999999774080E + 00.99999999998195110E + 00.99999999998440340E + 00

56.10000000000041760E + 01.10000000000193850E + 01

789.10000000000128360E + 01.99999999999306180E + 00.99999999998139700E + 00

1011.99999999998912740E + 00.10000000000090180E + 01

121314.10000000000187760E + 01.10000000000080810E + 01.99999999998862040E + 00

1516.99999999998185020E + 00.99999999999493800E + 00

171819.10000000000136250E + 01.10000000000169100E + 01
.10000000000019750E + 01

20.99999999998447330E + 00

で 12 桁は出る。
では ε を小さくすればよいかというとそうではない。あまり小さくすると，少し条件の悪い問題
では止まらないことがある。つまり，この停止則は共役勾配法では下手に使用しないことが望ましい。現在主に用いられているものとしては∥rk ∥ ≤ ε∥r0 ∥となったとき停止するというものである。
最後に，CG 法の計算精度に対する鋭敏性を示す例題を挙げる。
例題 10.4.3
係数行列 A と解 x を 100 99 · · · 99 99 · · ·A =  .
..
 ..
.

1
1 ···



 0 1 
 1 1  , x = ..  .. 
 . . 991とする。この問題に対して CG 法を，IEEE754 単精度 (float), 倍精度 (double)，128bit, 256bit, 512bit,1024bit でそれぞれ計算し，残差 ∥rk ∥2 = ∥b − Axk ∥2 をプロットすると図 10.1 のようになる。

第 10 章 連立一次方程式の解法 3 — Krylov 部分空間法

124||b-Axi||21IEEE754 float
1e-05
IEEE754 double
1e-10

gmp 128bits

1e-15

gmp256bits
1e-20

gmp 1024bits
02040

反復回数60gmp 512bits80

図 10.1: CG 法における丸め誤差の影響

通常は IEEE754 倍精度での数値実験結果のみを使用することが多いが，こうして多倍長計算も行ってみると，如何に CG 法が丸め誤差の影響に敏感であるかが明確となる。特に 512bit 計算でも僅かだが残差のノルムが上昇しており，完全に単調収束させるにはそれ以上の精度が必要であることがわかる。

演習問題
1. CG 法における初期残差 r0 を元に，係数行列を乗じた Ar0 , A2 r0 , ..., Ak r0 によって張られる線型部分空間 Kk を Krylov 部分空間と呼ぶ。この時pi ∈ Kk (i = 0, 1, 2, ..., k)であることを示せ2 。

参考文献
少し古いが，CG 法が考案された歴史的背景なども伺える邦文書として
共役勾配法
戸川隼人
教育出版

1977 年がある。
2 これが，CG 法を Krylov 部分空間法と呼んでいる理由である。