![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
共役勾配法†n元連立1次方程式  の係数行列Aが正値対称行列であるならば，
その解  ここで，(*,*)はベクトルの内積を表す．
まず，  となる．最小値を持つ点では関数の偏微分値(傾き)が0となる．
 Aは対称行列であるので，  この式は 共役勾配法の計算手順†共役勾配法の計算手順を以下に示す．
共役勾配法の実装†共役勾配法のC++による実装例を以下に示す． /*!
* 共役勾配法によりA・x=bを解く
* @param[in] A n×n正値対称行列
* @param[in] b 右辺ベクトル
* @param[out] x 結果ベクトル
* @param[in] n 行列の大きさ
* @param[inout] max_iter 最大反復数(反復終了後,実際の反復数を返す)
* @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す)
* @return 1:成功,0:失敗
*/
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);
// 第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){
// y = AP の計算
for(int i = 0; i < n; ++i){
y[i] = dot(A[i], p, n);
}
// alpha = r*r/(P*AP)の計算
alpha = rr0/dot(p, y, n);
// 解x、残差rの更新
for(int i = 0; i < n; ++i){
x[i] += alpha*p[i];
r[i] -= alpha*y[i];
}
// (r*r)_(k+1)の計算
rr1 = dot(r, r, n);
// 収束判定 (||r||<=eps)
e = sqrt(rr1);
if(e < eps){
k++;
break;
}
// βの計算とPの更新
beta = rr1/rr0;
for(int i = 0; i < n; ++i){
p[i] = r[i]+beta*p[i];
}
// (r*r)_(k+1)を次のステップのために確保しておく
rr0 = rr1;
}
max_iter = k+1;
eps = e;
return 1;
}
ここで，std::vectorの内積を求めるdot関数を定義して用いている． 条件数†ここでは前処理付き共役勾配法について説明する前に， 共役勾配法の収束性と係数行列の条件数について述べる． 共役勾配法は正定値対称行列にのみ適用できる．
なぜ対称行列出なければならないのかというのは方程式  は2次式であり，個々の式は  と変形する(要は解の公式を導出する過程の式)．
このとき 一方，  であるということである．ここで， 行列が正定値であること = その行列の固有値がすべて正 であるので，行列の固有値がすべて正ならば最小値を持つことを証明する． n元連立1次方程式  となる．ここで，  よって，  である．方程式  ここで，  と変形できるので，  となる．さらに，  Dは対角行列なので  右辺第1項は2次形式である．そして，Dは対角行列で，その対角成分，つまり，固有値が これでAが正定値でなければならないことは証明できた．
次に，  の形になる．この式は  と書け，2つのノルムの積となる． 正定値実対称行列の場合，最大固有値を最小固有値で割った比が条件数になる(すべての固有値は正であることに注意)． 条件数の最小値は1であり，条件数が1に近いほど誤差が小さく，アルゴリズムの収束性がよくなるため， 連立1次方程式の精度や収束性を評価する上でとても重要な指標となる． 前処理付き共役勾配法†前処理付き共役勾配法†係数行列Aの条件数が大きい場合，つまり，  ここで， 不完全コレスキー分解付き共役勾配法†不完全コレスキー分解により係数行列は以下のように分解される．  対角行列  ここで，  と置いた．
Uを使って 
 両辺に  この式から，  と計算するのが不完全コレスキー分解付共役勾配法(Incomplete Cholesky Conjugate Gradient method : ICCG法)である
(ちなみに，コレスキー分解できればよいので修正コレスキー分解でも同じ)．
前節の式に当てはめた場合， 変形した係数行列 
 と単位行列となる．単位行列の最小・最大固有値はともに1なので， もし正確に分解されているならば，係数行列の条件数は1になっている． 修正コレスキー分解ならばこの処理で常に条件数1が得られるが，計算速度の問題から近似的な分解ではあるが， 不完全コレスキー分解がよく用いられる． 次に共役勾配法の各手順を変換して，ICCG法に対応させる． ICCG法で得られる係数行列，変数，右辺項などを以下のように表す．  共役勾配法の各手順について変換した式を求める．
不完全コレスキー分解付き共役勾配法の計算手順†不完全コレスキー分解付き共役勾配法の計算手順を以下に示す(
ここで，  を前進代入，後退代入で解く． 不完全コレスキー分解付き共役勾配法の実装†不完全コレスキー分解付きの共役勾配法のC++による実装例を以下に示す． /*!
* 不完全コレスキー分解による前処理付共役勾配法によりA・x=bを解く
* @param[in] A n×n正値対称行列
* @param[in] b 右辺ベクトル
* @param[out] x 結果ベクトル
* @param[in] n 行列の大きさ
* @param[inout] max_iter 最大反復数(反復終了後,実際の反復数を返す)
* @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す)
* @return 1:成功,0:失敗
*/
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);
// 第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_0 = (LDL^T)^-1 r_0 の計算
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){
// y = AP の計算
for(int i = 0; i < n; ++i){
y[i] = dot(A[i], p, n);
}
// alpha = r*r/(P*AP)の計算
alpha = rr0/dot(p, y, n);
// 解x、残差rの更新
for(int i = 0; i < n; ++i){
x[i] += alpha*p[i];
r[i] -= alpha*y[i];
}
// (r*r)_(k+1)の計算
ICRes(L, d, r, r2, n);
rr1 = dot(r, r2, n);
// 収束判定 (||r||<=eps)
e = sqrt(rr1);
if(e < eps){
k++;
break;
}
// βの計算とPの更新
beta = rr1/rr0;
for(int i = 0; i < n; ++i){
p[i] = r2[i]+beta*p[i];
}
// (r*r)_(k+1)を次のステップのために確保しておく
rr0 = rr1;
}
max_iter = k;
eps = e;
return 1;
}
ここで，dot関数はstd::vectorの内積を求める関数， IncompleteCholeskyDecomp2関数は三角分解のところで説明した不完全コレスキー分解のコードである． また，ICRes関数は前進代入，後退代入で /*!
* (LDL^T)^-1 r の計算
* @param[in] L,d IC分解で得られた下三角行列と対角行列(対角成分のみのベクトル)
* @param[in] r 残差ベクトル
* @param[in] u (LDL^T)^-1 rを計算した結果
*/
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;
}
}
参考文献†
|
を求めることは，
以下の方程式の最小値探索と同等である．
を直接解く代わりに，
方程式
の最小値探索により解を求める方法を
共役勾配法(conjugate gradient method:CG法)と呼ぶ．
で偏微分すると，右辺第1項はi=kとj=kの要素のみが残るので，
である．よって，
を適当に設定
);
を計算
を計算
ならば収束したとして計算を終了
の修正係数
を計算
を計算
の最小化に関しての記述のところで説明した．一方でなぜ正定値(positive definite)でなければならないのかを考える．
方程式
の形で表すことができる．
この2次式が最小値を持つ条件を考える．
まず，2次式
にかかる係数
が
ならば2次式のグラフが上に開いた形になり，最小値を持つ．
逆に
だと2次式は最大値を持ち，
で最小値も最大値も持たない．
大事なのは係数
の場合はどうだろうか．
まず，Aが正定値であるということは，
の任意の実ベクトル
について，
を行列Aの2次形式と呼ぶ．
の係数行列Aの
モード行列をR，固有値を対角成分にした対角行列をDとする．
モード行列とは行列の固有ベクトルを並べた行列で，
の直交行列である．
係数行列は
で表されるので，連立1次方程式の式は，
とおくと，
を
を使って表すと，
とおくと，
である．よって，
の係数aになる．
上で示したように，2次式が最小値を持つ条件は
に関する式が実際にどのようになっているかを考える．
2次元(
)のとき，
を軸とする楕円(
)であり，
固有値(
)の比は，軸の長さ
の平方根の比と等しい．
つまり，固有値の最大値と最小値の差が大きいと平べったい楕円になってしまう．
ここで，行列Aの条件数(condition number)というものを考える．
条件数は，
が平べったい楕円体になっていると，
誤差が大きくなり，収束が遅くなる．
収束を早めたい場合は，条件数を小さくするように前処理としてAを変形すればよい．
つまり，
は
の正則行列である．
このときに
の条件数が1に近いと高速な収束することが期待できる．
このように前処理を施すことで収束を早める方法が前処理付きクリロフ部分空間法(Preconditioned Krylov subspace method)である．
共役勾配法に適用した場合は，前処理付き共役勾配法(Preconditioned Conjugate Gradiate method : PCG法)と呼ぶ．
ちなみに，固有値の最大値と最小値の比が収束性に影響するのは共役勾配法だけでなく，
他のクリロフ部分空間法でも同様であるので，上記の前処理をそれらに適用することもできる．
ここではよく用いられる不完全コレスキー分解付き共役勾配法(Incomplete Cholesky Conjugate Gradient method : ICCG法)について述べる．
に対して，
を定義する．
Dは対角行列であるので
である．
Aをこれを用いて表すと，
の式を書き換える．
より
なので，
を掛ける．
である．
について考える．
係数行列を変形すると，
より，
なので，
について
について
より，
について
なので，
について
について
について
について
としている)．
を適当に設定
);
を計算
ならば収束したとして計算を終了
を残差に掛ける
の修正係数
と
の計算は，
を計算する関数である．