前処理付き共役勾配法

係数行列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 736件 [詳細] filels_pcg.eq15.gif 746件 [詳細] filels_pcg.eq2.gif 802件 [詳細] filels_pcg.eq13.gif 825件 [詳細] filels_pcg.eq50.gif 709件 [詳細] filels_pcg.eq25.gif 704件 [詳細] filels_pcg.eq31.gif 674件 [詳細] filels_pcg.eq6.gif 604件 [詳細] filels_pcg.eq39.gif 610件 [詳細] filels_pcg.eq16.gif 668件 [詳細] filels_pcg.eq57.gif 595件 [詳細] filels_pcg.eq46.gif 689件 [詳細] filels_pcg.eq37.gif 736件 [詳細] filels_pcg.eq38.gif 743件 [詳細] filels_pcg.eq24.gif 776件 [詳細] filels_pcg.eq44.gif 660件 [詳細] filels_pcg.eq66.gif 632件 [詳細] filels_pcg.eq62.gif 615件 [詳細] filels_pcg.eq55.gif 643件 [詳細] filels_pcg.eq63.gif 613件 [詳細] filels_pcg.eq1.gif 787件 [詳細] filels_pcg.eq60.gif 588件 [詳細] filels_pcg.eq35.gif 617件 [詳細] filels_pcg.eq33.gif 656件 [詳細] filels_pcg.eq34.gif 641件 [詳細] filels_pcg.eq43.gif 631件 [詳細] filels_pcg.eq4.gif 681件 [詳細] filels_pcg.eq48.gif 638件 [詳細] filels_pcg.eq22.gif 601件 [詳細] filels_pcg.eq5.gif 632件 [詳細] filels_pcg.eq42.gif 594件 [詳細] filels_pcg.eq14.gif 686件 [詳細] filels_pcg.eq41.gif 623件 [詳細] filels_pcg.eq40.gif 556件 [詳細] filels_pcg.eq12.gif 643件 [詳細] filels_pcg.eq51.gif 644件 [詳細] filels_pcg.eq53.gif 561件 [詳細] filels_pcg.eq32.gif 569件 [詳細] filels_pcg.eq17.gif 663件 [詳細] filels_pcg.eq52.gif 619件 [詳細] filels_pcg.eq8.gif 604件 [詳細] filels_pcg.eq67.gif 504件 [詳細] filels_pcg.eq27.gif 682件 [詳細] filels_pcg.eq18.gif 628件 [詳細] filels_pcg.eq10.gif 663件 [詳細] filels_pcg.eq21.gif 695件 [詳細] filels_pcg.eq59.gif 554件 [詳細] filels_pcg.eq11.gif 637件 [詳細] filels_pcg.eq54.gif 643件 [詳細] filels_pcg.eq65.gif 595件 [詳細] filels_pcg.eq49.gif 644件 [詳細] filels_pcg.eq9.gif 529件 [詳細] filels_pcg.eq58.gif 607件 [詳細] filels_pcg.eq47.gif 541件 [詳細] filels_pcg.eq29.gif 680件 [詳細] filels_pcg.eq45.gif 591件 [詳細] filels_pcg.eq68.gif 589件 [詳細] filels_pcg.eq26.gif 663件 [詳細] filels_pcg.eq28.gif 654件 [詳細] filels_pcg.eq61.gif 562件 [詳細] filels_pcg.eq7.gif 593件 [詳細] filels_pcg.eq20.gif 505件 [詳細] filels_pcg.eq64.gif 556件 [詳細] filels_pcg.eq36.gif 598件 [詳細] filels_pcg.eq30.gif 682件 [詳細] filels_pcg.eq56.gif 634件 [詳細] filels_pcg.eq23.gif 709件 [詳細] filels_pcg.eq3.gif 606件 [詳細]

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