共役勾配法†
n元連立1次方程式
の係数行列Aが正値対称行列であるならば，
その解
を求めることは，
以下の方程式の最小値探索と同等である．
ここで，(*,*)はベクトルの内積を表す．
を直接解く代わりに，
方程式
の最小値探索により解を求める方法を
共役勾配法(conjugate gradient method:CG法)と呼ぶ．
の最小値探索が元のn元連立1次方程式を解くことと同等になることの証明を以下に示す．
まず，を成分で表すと，
となる．最小値を持つ点では関数の偏微分値(傾き)が0となる．
で偏微分すると，右辺第1項はi=kとj=kの要素のみが残るので，
Aは対称行列であるので，
である．よって，
この式はを成分で表示したものである．
よって，の解はを最小にすると同じになる．
共役勾配法の計算手順†
共役勾配法の計算手順を以下に示す．
- 初期近似解
を適当に設定
- 初期近似解に対する残差を計算，修正方向ベクトルを
);で初期化
- 以下の計算を収束するまで繰り返す(k=0,1,2,...)
を計算- 修正係数
を計算
- k+1ステップの近似値を算出
- k+1の近似値に対する残差を計算
ならば収束したとして計算を終了- 方向ベクトル
の修正係数
を計算
- k+1ステップの方向ベクトル
を計算
共役勾配法の実装†
共役勾配法のC++による実装例を以下に示す．
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
| |
int CGSolver(const vector< vector<double> > &A, const vector<double> &b, vector<double> &x, int n, int &max_iter, double &eps)
{
if(n <= 0) return 0;
vector<double> r(n), p(n), y(n);
x.assign(n, 0.0);
for(int i = 0; i < n; ++i){
double ax = 0.0;
for(int j = 0; j < n; ++j){
ax += A[i][j]*x[j];
}
r[i] = b[i]-ax;
p[i] = r[i];
}
double rr0 = dot(r, r, n), rr1;
double alpha, beta;
double e = 0.0;
int k;
for(k = 0; k < max_iter; ++k){
for(int i = 0; i < n; ++i){
y[i] = dot(A[i], p, n);
}
alpha = rr0/dot(p, y, n);
for(int i = 0; i < n; ++i){
x[i] += alpha*p[i];
r[i] -= alpha*y[i];
}
rr1 = dot(r, r, n);
e = sqrt(rr1);
if(e < eps){
k++;
break;
}
beta = rr1/rr0;
for(int i = 0; i < n; ++i){
p[i] = r[i]+beta*p[i];
}
rr0 = rr1;
}
max_iter = k+1;
eps = e;
return 1;
}
|
ここで，std::vectorの内積を求めるdot関数を定義して用いている．
条件数†
ここでは前処理付き共役勾配法について説明する前に，
共役勾配法の収束性と係数行列の条件数について述べる．
共役勾配法は正定値対称行列にのみ適用できる．
なぜ対称行列出なければならないのかというのは方程式
の最小化に関しての記述のところで説明した．一方でなぜ正定値(positive definite)でなければならないのかを考える．
方程式
は2次式であり，個々の式は
の形で表すことができる．
この2次式が最小値を持つ条件を考える．
まず，2次式を，
と変形する(要は解の公式を導出する過程の式)．
このとき
にかかる係数
が
ならば2次式のグラフが上に開いた形になり，最小値を持つ．
逆に
だと2次式は最大値を持ち，
で最小値も最大値も持たない．
大事なのは係数の符号である．
一方，
の場合はどうだろうか．
まず，Aが正定値であるということは，
の任意の実ベクトル
について，
であるということである．ここで，
を行列Aの2次形式と呼ぶ．
行列が正定値であること = その行列の固有値がすべて正
であるので，行列の固有値がすべて正ならば最小値を持つことを証明する．
n元連立1次方程式
の係数行列Aの
モード行列をR，固有値を対角成分にした対角行列をDとする．
モード行列とは行列の固有ベクトルを並べた行列で，
の直交行列である．
係数行列は
で表されるので，連立1次方程式の式は，
となる．ここで，
とおくと，
よって，
である．方程式
を
を使って表すと，
ここで，
と変形できるので，
となる．さらに，
とおくと，
Dは対角行列なので
である．よって，
右辺第1項は2次形式である．そして，Dは対角行列で，その対角成分，つまり，固有値が
の係数aになる．
上で示したように，2次式が最小値を持つ条件はであるので，
すべての固有値が正ならばは最小値を持つ．
つまり，言い換えるとAが正定値ならばは最小値を持つ．
これでAが正定値でなければならないことは証明できた．
次に，
に関する式が実際にどのようになっているかを考える．
2次元(
)のとき，
の形になる．この式は
を軸とする楕円(
)であり，
固有値(
)の比は，軸の長さ
の平方根の比と等しい．
つまり，固有値の最大値と最小値の差が大きいと平べったい楕円になってしまう．
ここで，行列Aの条件数(condition number)というものを考える．
条件数は，
と書け，2つのノルムの積となる．
正定値実対称行列の場合，最大固有値を最小固有値で割った比が条件数になる(すべての固有値は正であることに注意)．
条件数の最小値は1であり，条件数が1に近いほど誤差が小さく，アルゴリズムの収束性がよくなるため，
連立1次方程式の精度や収束性を評価する上でとても重要な指標となる．
前処理付き共役勾配法†
前処理付き共役勾配法†
係数行列Aの条件数が大きい場合，つまり，
が平べったい楕円体になっていると，
誤差が大きくなり，収束が遅くなる．
収束を早めたい場合は，条件数を小さくするように前処理としてAを変形すればよい．
つまり，
ここで，
は
の正則行列である．
このときに
の条件数が1に近いと高速な収束することが期待できる．
このように前処理を施すことで収束を早める方法が前処理付きクリロフ部分空間法(Preconditioned Krylov subspace method)である．
共役勾配法に適用した場合は，前処理付き共役勾配法(Preconditioned Conjugate Gradiate method : PCG法)と呼ぶ．
ちなみに，固有値の最大値と最小値の比が収束性に影響するのは共役勾配法だけでなく，
他のクリロフ部分空間法でも同様であるので，上記の前処理をそれらに適用することもできる．
ここではよく用いられる不完全コレスキー分解付き共役勾配法(Incomplete Cholesky Conjugate Gradient method : ICCG法)について述べる．
不完全コレスキー分解付き共役勾配法†
不完全コレスキー分解により係数行列は以下のように分解される．
対角行列
に対して，
を定義する．
Dは対角行列であるので
である．
Aをこれを用いて表すと，
ここで，
と置いた．
Uを使って
の式を書き換える．
より
なので，
両辺に
を掛ける．
この式から，
と計算するのが不完全コレスキー分解付共役勾配法(Incomplete Cholesky Conjugate Gradient method : ICCG法)である
(ちなみに，コレスキー分解できればよいので修正コレスキー分解でも同じ)．
前節の式に当てはめた場合，
である．
変形した係数行列
について考える．
係数行列を変形すると，
より，
なので，
と単位行列となる．単位行列の最小・最大固有値はともに1なので，
もし正確に分解されているならば，係数行列の条件数は1になっている．
修正コレスキー分解ならばこの処理で常に条件数1が得られるが，計算速度の問題から近似的な分解ではあるが，
不完全コレスキー分解がよく用いられる．
次に共役勾配法の各手順を変換して，ICCG法に対応させる．
ICCG法で得られる係数行列，変数，右辺項などを以下のように表す．
共役勾配法の各手順について変換した式を求める．
- 残差ベクトル
について
よって，
- 方向ベクトル
について
とおくと，
より，
よって，
- 修正係数
について
なので，
よって，
- 修正係数
について
と同様に変換すると以下が得られる．
について
よって，
- 残差ベクトルの更新値
について
よって，
- 方向ベクトルの更新値
について
よって，
不完全コレスキー分解付き共役勾配法の計算手順†
不完全コレスキー分解付き共役勾配法の計算手順を以下に示す(
としている)．
- 初期近似解
を適当に設定
- 不完全コレスキー分解によりL,Dを計算
- 初期近似解に対する残差を計算，修正方向ベクトルを
);で初期化
- 以下の計算を収束するまで繰り返す(k=0,1,2,...)
を計算
- 修正係数を計算
- k+1ステップの近似値を算出
- k+1の近似値に対する残差を計算
ならば収束したとして計算を終了
を残差に掛ける
- 方向ベクトル
の修正係数を計算
- k+1ステップの方向ベクトルを計算
ここで，
と
の計算は，
を直接計算するのではなく，
を前進代入，後退代入で解く．
不完全コレスキー分解付き共役勾配法の実装†
不完全コレスキー分解付きの共役勾配法のC++による実装例を以下に示す．
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
| |
int ICCGSolver(const vector< vector<double> > &A, const vector<double> &b, vector<double> &x, int n, int &max_iter, double &eps)
{
if(n <= 0) return 0;
vector<double> r(n), p(n), y(n), r2(n);
x.assign(n, 0.0);
vector<double> d(n);
vector< vector<double> > L(n, vector<double>(n, 0.0));
IncompleteCholeskyDecomp2(A, L, d, n);
for(int i = 0; i < n; ++i){
double ax = 0.0;
for(int j = 0; j < n; ++j){
ax += A[i][j]*x[j];
}
r[i] = b[i]-ax;
}
ICRes(L, d, r, p, n);
float rr0 = dot(r, p, n), rr1;
float alpha, beta;
double e = 0.0;
int k;
for(k = 0; k < max_iter; ++k){
for(int i = 0; i < n; ++i){
y[i] = dot(A[i], p, n);
}
alpha = rr0/dot(p, y, n);
for(int i = 0; i < n; ++i){
x[i] += alpha*p[i];
r[i] -= alpha*y[i];
}
ICRes(L, d, r, r2, n);
rr1 = dot(r, r2, n);
e = sqrt(rr1);
if(e < eps){
k++;
break;
}
beta = rr1/rr0;
for(int i = 0; i < n; ++i){
p[i] = r2[i]+beta*p[i];
}
rr0 = rr1;
}
max_iter = k;
eps = e;
return 1;
}
|
ここで，dot関数はstd::vectorの内積を求める関数，
IncompleteCholeskyDecomp2関数は三角分解のところで説明した不完全コレスキー分解のコードである．
また，ICRes関数は前進代入，後退代入で
を計算する関数である．
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
| |
inline void ICRes(const vector< vector<double> > &L, const vector<double> &d, const vector<double> &r, vector<double> &u, int n)
{
vector<double> y(n);
for(int i = 0; i < n; ++i){
double rly = r[i];
for(int j = 0; j < i; ++j){
rly -= L[i][j]*y[j];
}
y[i] = rly/L[i][i];
}
for(int i = n-1; i >= 0; --i){
double lu = 0.0;
for(int j = i+1; j < n; ++j){
lu += L[j][i]*u[j];
}
u[i] = y[i]-d[i]*lu;
}
}
|
参考文献†
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
