前処理付き共役勾配法

前処理付き共役勾配法†

係数行列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法で得られる係数行列，変数，右辺項などを以下のように表す．

共役勾配法の各手順について変換した式を求める．

• 残差ベクトルについて
  
  よって，

• 方向ベクトルについて
  
  とおくと，より，
  
  よって，

• 修正係数について
  
  なので，
  
  よって，

• 修正係数について
  
  と同様に変換すると以下が得られる．

• について
  
  よって，

• 残差ベクトルの更新値について
  
  よって，

• 方向ベクトルの更新値について
  
  よって，

不完全コレスキー分解付き共役勾配法の計算手順†

不完全コレスキー分解付き共役勾配法の計算手順を以下に示す(としている)．

1. 初期近似解を適当に設定
2. 不完全コレスキー分解によりL,Dを計算
3. 初期近似解に対する残差を計算，修正方向ベクトルを);で初期化
   
4. 以下の計算を収束するまで繰り返す(k=0,1,2,...)
5. を計算
   1. 修正係数を計算
      
   2. k+1ステップの近似値を算出
      
   3. k+1の近似値に対する残差を計算
      
   4. ならば収束したとして計算を終了
   5. を残差に掛ける
      
   6. 方向ベクトルの修正係数を計算
      
   7. k+1ステップの方向ベクトルを計算
      

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

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

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

不完全コレスキー分解付きの共役勾配法の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;
    }
}