前処理付き共役勾配法

係数行列Aの条件数が大きい場合,つまり,ls_pcg.eq1.gifが平べったい楕円体になっていると, 誤差が大きくなり,収束が遅くなる. 収束を早めたい場合は,条件数を小さくするように前処理としてAを変形すればよい. つまり,

ls_pcg.eq2.gif

ここで,ls_pcg.eq3.gifls_pcg.eq4.gifの正則行列である. このときにls_pcg.eq5.gifの条件数が1に近いと高速な収束することが期待できる. このように前処理を施すことで収束を早める方法が前処理付きクリロフ部分空間法(Preconditioned Krylov subspace method)である. 共役勾配法に適用した場合は,前処理付き共役勾配法(Preconditioned Conjugate Gradiate method : PCG法)と呼ぶ. ちなみに,固有値の最大値と最小値の比が収束性に影響するのは共役勾配法だけでなく, 他のクリロフ部分空間法でも同様であるので,上記の前処理をそれらに適用することもできる. ここではよく用いられる不完全コレスキー分解付き共役勾配法(Incomplete Cholesky Conjugate Gradient method : ICCG法)について述べる.

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

不完全コレスキー分解により係数行列は以下のように分解される.

ls_pcg.eq6.gif

対角行列ls_pcg.eq7.gifに対して,ls_pcg.eq8.gifを定義する. Dは対角行列であるのでls_pcg.eq9.gifである. Aをこれを用いて表すと,

ls_pcg.eq10.gif

ここで,

ls_pcg.eq11.gif

と置いた. Uを使ってls_pcg.eq12.gifの式を書き換える.

ls_pcg.eq13.gif

ls_pcg.eq14.gifよりls_pcg.eq15.gifなので,

ls_pcg.eq16.gif

両辺にls_pcg.eq17.gifを掛ける.

ls_pcg.eq18.gif

この式から,

ls_pcg.eq19.gif

と計算するのが不完全コレスキー分解付共役勾配法(Incomplete Cholesky Conjugate Gradient method : ICCG法)である (ちなみに,コレスキー分解できればよいので修正コレスキー分解でも同じ). 前節の式に当てはめた場合,ls_pcg.eq20.gifである.

変形した係数行列ls_pcg.eq21.gifについて考える. 係数行列を変形すると,

ls_pcg.eq22.gif

ls_pcg.eq23.gifより,ls_pcg.eq24.gifなので,

ls_pcg.eq25.gif

と単位行列となる.単位行列の最小・最大固有値はともに1なので, もし正確に分解されているならば,係数行列の条件数は1になっている. 修正コレスキー分解ならばこの処理で常に条件数1が得られるが,計算速度の問題から近似的な分解ではあるが, 不完全コレスキー分解がよく用いられる.

次に共役勾配法の各手順を変換して,ICCG法に対応させる. ICCG法で得られる係数行列,変数,右辺項などを以下のように表す.

ls_pcg.eq26.gif

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

  • 残差ベクトルls_pcg.eq27.gifについて
    ls_pcg.eq28.gif
    よって,
    ls_pcg.eq29.gif
  • 方向ベクトルls_pcg.eq30.gifについて
    ls_pcg.eq31.gif
    とおくと,ls_pcg.eq32.gifより,
    ls_pcg.eq33.gif
    よって,
    ls_pcg.eq34.gif
  • 修正係数ls_pcg.eq35.gifについて
    ls_pcg.eq36.gif
    ls_pcg.eq37.gifなので,
    ls_pcg.eq38.gif
    よって,
    ls_pcg.eq39.gif
  • 修正係数ls_pcg.eq40.gifについて
    ls_pcg.eq41.gif
    ls_pcg.eq35.gifと同様に変換すると以下が得られる.
    ls_pcg.eq42.gif
  • ls_pcg.eq43.gifについて
    ls_pcg.eq44.gif
    よって,
    ls_pcg.eq45.gif
  • 残差ベクトルの更新値ls_pcg.eq46.gifについて
    ls_pcg.eq47.gif
    よって,
    ls_pcg.eq48.gif
  • 方向ベクトルの更新値ls_pcg.eq49.gifについて
    ls_pcg.eq50.gif
    よって,
    ls_pcg.eq51.gif

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

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

  1. 初期近似解ls_pcg.eq53.gifを適当に設定
  2. 不完全コレスキー分解によりL,Dを計算
  3. 初期近似解に対する残差を計算,修正方向ベクトルを70%&ref(ls_pcg.eq30.gif,nolink);ls_pcg.eq54.gifで初期化
    ls_pcg.eq55.gif
  4. 以下の計算を収束するまで繰り返す(k=0,1,2,...)
  5. ls_pcg.eq56.gifを計算
    1. 修正係数ls_pcg.eq35.gifを計算
      ls_pcg.eq57.gif
    2. k+1ステップの近似値を算出
      ls_pcg.eq45.gif
    3. k+1の近似値に対する残差を計算
      ls_pcg.eq58.gif
    4. ls_pcg.eq59.gifならば収束したとして計算を終了
    5. ls_pcg.eq60.gifを残差に掛ける
      ls_pcg.eq61.gif
    6. 方向ベクトルls_pcg.eq62.gifの修正係数ls_pcg.eq40.gifを計算
      ls_pcg.eq63.gif
    7. k+1ステップの方向ベクトルls_pcg.eq49.gifを計算
      ls_pcg.eq64.gif

ここで,ls_pcg.eq65.gifls_pcg.eq66.gifの計算は, ls_pcg.eq60.gifを直接計算するのではなく,

ls_pcg.eq67.gif

を前進代入,後退代入で解く.

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

不完全コレスキー分解付きの共役勾配法の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
/*!
 * 不完全コレスキー分解による前処理付共役勾配法により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関数は前進代入,後退代入でls_pcg.eq68.gifを計算する関数である.

  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
/*!
 * (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;
    }
}

添付ファイル: filels_pcg.eq19.gif 895件 [詳細] filels_pcg.eq15.gif 905件 [詳細] filels_pcg.eq2.gif 961件 [詳細] filels_pcg.eq13.gif 973件 [詳細] filels_pcg.eq50.gif 861件 [詳細] filels_pcg.eq25.gif 814件 [詳細] filels_pcg.eq31.gif 775件 [詳細] filels_pcg.eq6.gif 687件 [詳細] filels_pcg.eq39.gif 704件 [詳細] filels_pcg.eq16.gif 767件 [詳細] filels_pcg.eq57.gif 677件 [詳細] filels_pcg.eq46.gif 788件 [詳細] filels_pcg.eq37.gif 845件 [詳細] filels_pcg.eq38.gif 851件 [詳細] filels_pcg.eq24.gif 879件 [詳細] filels_pcg.eq44.gif 767件 [詳細] filels_pcg.eq66.gif 716件 [詳細] filels_pcg.eq62.gif 696件 [詳細] filels_pcg.eq55.gif 748件 [詳細] filels_pcg.eq63.gif 694件 [詳細] filels_pcg.eq1.gif 903件 [詳細] filels_pcg.eq60.gif 675件 [詳細] filels_pcg.eq35.gif 702件 [詳細] filels_pcg.eq33.gif 749件 [詳細] filels_pcg.eq34.gif 719件 [詳細] filels_pcg.eq43.gif 712件 [詳細] filels_pcg.eq4.gif 765件 [詳細] filels_pcg.eq48.gif 716件 [詳細] filels_pcg.eq22.gif 745件 [詳細] filels_pcg.eq5.gif 712件 [詳細] filels_pcg.eq42.gif 683件 [詳細] filels_pcg.eq14.gif 770件 [詳細] filels_pcg.eq41.gif 710件 [詳細] filels_pcg.eq40.gif 640件 [詳細] filels_pcg.eq12.gif 729件 [詳細] filels_pcg.eq51.gif 724件 [詳細] filels_pcg.eq53.gif 644件 [詳細] filels_pcg.eq32.gif 644件 [詳細] filels_pcg.eq17.gif 749件 [詳細] filels_pcg.eq52.gif 704件 [詳細] filels_pcg.eq8.gif 668件 [詳細] filels_pcg.eq67.gif 573件 [詳細] filels_pcg.eq27.gif 765件 [詳細] filels_pcg.eq18.gif 724件 [詳細] filels_pcg.eq10.gif 744件 [詳細] filels_pcg.eq21.gif 778件 [詳細] filels_pcg.eq59.gif 618件 [詳細] filels_pcg.eq11.gif 722件 [詳細] filels_pcg.eq54.gif 741件 [詳細] filels_pcg.eq65.gif 652件 [詳細] filels_pcg.eq49.gif 723件 [詳細] filels_pcg.eq9.gif 595件 [詳細] filels_pcg.eq58.gif 672件 [詳細] filels_pcg.eq47.gif 626件 [詳細] filels_pcg.eq29.gif 761件 [詳細] filels_pcg.eq45.gif 677件 [詳細] filels_pcg.eq68.gif 655件 [詳細] filels_pcg.eq26.gif 756件 [詳細] filels_pcg.eq28.gif 780件 [詳細] filels_pcg.eq61.gif 619件 [詳細] filels_pcg.eq7.gif 657件 [詳細] filels_pcg.eq20.gif 588件 [詳細] filels_pcg.eq64.gif 621件 [詳細] filels_pcg.eq36.gif 679件 [詳細] filels_pcg.eq30.gif 774件 [詳細] filels_pcg.eq56.gif 734件 [詳細] filels_pcg.eq23.gif 797件 [詳細] filels_pcg.eq3.gif 704件 [詳細]

トップ   編集 凍結 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2012-07-03 (火) 11:04:04