前処理付き共役勾配法

前処理付き共役勾配法†

係数行列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,...)
を計算
1. 修正係数を計算
2. k+1ステップの近似値を算出
3. k+1の近似値に対する残差を計算
4. ならば収束したとして計算を終了
5. を残差に掛ける
6. 方向ベクトルの修正係数を計算
7. k+1ステップの方向ベクトルを計算

ここで，との計算は，を直接計算するのではなく，

を前進代入，後退代入で解く．

↑

不完全コレスキー分解付き共役勾配法の実装†

不完全コレスキー分解付きの共役勾配法のC++による実装例を以下に示す．


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関数は前進代入，後退代入でを計算する関数である．


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;
    }
}

Menu