LU分解では係数行列をls_cholesky.eq1.gifの形に分解したが, Aが正定値対称行列であるならば,

ls_cholesky.eq2.gif

と分解できる下三角行列Lが存在する.この式はLU分解でls_cholesky.eq3.gifと置いたものに相当する. このような分解がコレスキー分解(Cholesky decomposition)である.

上式を成分で表す. ls_cholesky.eq4.gifを,

ls_cholesky.eq5.gif

とすると,

ls_cholesky.eq6.gif

となる.ls_cholesky.eq7.gifの場合の要素を分けて書くと,

ls_cholesky.eq8.gif

よって,以下の式をls_cholesky.eq9.gifと計算していくことでLを求めることができる.

ls_cholesky.eq10.gif

たとえば,

ls_cholesky.eq11.gif

という手順で計算していく. 対称行列に限定されるものの,LU分解と異なり,Lだけを求めればよいので計算手順は半分で済む.

下三角行列Lが求まったら,LU分解のときと同様に前進代入,後退代入の手順で連立1次方程式ls_cholesky.eq12.gifを解くことができる.

ls_cholesky.eq13.gif

コレスキー分解では計算に平方根が含まれており,平方根の中が正でなければならない. また,0でもls_cholesky.eq14.gifの計算でゼロ割が発生してしまうので,

ls_cholesky.eq15.gif

でなければならない. この条件は対称正定値行列のときには満たされることが証明されているが, 正定値でない対称行列では問題がある. また,計算コストがかかる平方根の計算もできれば避けたい. そのため,コレスキー分解はあまり使われず, 以下で述べる修正コレスキー分解や不完全コレスキー分解がよく用いられる.

修正コレスキー分解

コレスキー分解の平方根の問題を解決するために改良を加えられたのが 修正コレスキー分解(Modified Cholosky decomposition)である.

まず,コレスキー分解の下三角行列ls_cholesky.eq4.gifを対角成分が1の下三角行列と対角行列の積に分解する.

ls_cholesky.eq16.gif

この式をコレスキー分解に当てはめると,

ls_cholesky.eq17.gif

ここで,ls_cholesky.eq18.gifls_cholesky.eq19.gifとおくと,

ls_cholesky.eq20.gif

と分解できる.これが修正コレスキー分解である.

この式を成分で表すと以下となる.

ls_cholesky.eq21.gif

コレスキー分解の時と同様にls_cholesky.eq7.gifの要素を分けて書くと,

ls_cholesky.eq22.gif

ここで,ls_cholesky.eq23.gifなので,

ls_cholesky.eq24.gif

よって,以下の式をls_cholesky.eq25.gifと計算していくことでL,Dを求めることができる.

ls_cholesky.eq26.gif

ただし,ls_cholesky.eq27.gifである.

L,Dが求まったら,LU分解のときと同様に前進代入,後退代入の手順で連立1次方程式ls_cholesky.eq12.gifを解くことができる.

ls_cholesky.eq28.gif

計算式をみてもわかるように,L,Dに分解することで対角成分の2乗項ls_cholesky.eq29.gifがなくなり, 最終的な計算手順での平方根の計算が必要なくなる. ここではls_cholesky.eq23.gifとなるように分解したが,ls_cholesky.eq30.gifとしてもよい (つまり,ls_cholesky.eq31.gif). この場合は,ls_cholesky.eq32.gifに関する式が以下のようになる.

ls_cholesky.eq33.gif

よって,以下の式をls_cholesky.eq25.gifと計算していくことでL,Dを求めることができる.

ls_cholesky.eq34.gif

ここで,ls_cholesky.eq35.gifである. 実際には上記の2式は1つにまとめることができる.そのため,式がシンプルになり, 条件分けがない分,実装コードも短くなるので私はこちらの方が好みである.

修正コレスキー分解の実装

修正コレスキー分解をC++で実装した例を以下に示す. それぞれ,ls_cholesky.eq23.gifの場合,ls_cholesky.eq36.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
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
/*!
 * 修正コレスキー分解(modified Cholesky decomposition)
 *  - 対称行列A(n×n)を下三角行列(L:Lower triangular matrix)と対角行列の積(LDL^T)に分解する
 *  - l_ii = 1とした場合
 *  - L: i > jの要素が非ゼロで対角成分は1
 * @param[in] A n×nの対称行列
 * @param[out] L 対角成分が1の下三角行列
 * @param[out] d 対角行列(対角成分のみ)
 * @param[in] n 行列の大きさ
 * @return 1:成功,0:失敗
 */
int ModifiedCholeskyDecomp(const vector< vector<double> > &A, vector< vector<double> > &L, vector<double> &d, int n)
{
    if(n <= 0) return 0;
 
    d[0] = A[0][0];
    L[0][0] = 1.0;
 
    for(int i = 1; i < n; ++i){
        // i < k の場合
        for(int j = 0; j < i; ++j){
            double lld = A[i][j];
            for(int k = 0; k < j; ++k){
                lld -= L[i][k]*L[j][k]*d[k];
            }
            L[i][j] = (1.0/d[j])*lld;
        }
 
        // i == k の場合
        double ld = A[i][i];
        for(int k = 0; k < i; ++k){
            ld -= L[i][k]*L[i][k]*d[k];
        }
        d[i] = ld;
        L[i][i] = 1.0;
    }
 
    return 1;
}
  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
/*!
 * 修正コレスキー分解(modified Cholesky decomposition)
 *  - 対称行列A(n×n)を下三角行列(L:Lower triangular matrix)と対角行列の積(LDL^T)に分解する
 *  - l_ii * d_i = 1とした場合
 *  - L: i > jの要素が非ゼロで対角成分は1
 * @param[in] A n×nの対称行列
 * @param[out] L 対角成分が1の下三角行列
 * @param[out] d 対角行列(対角成分のみ)
 * @param[in] n 行列の大きさ
 * @return 1:成功,0:失敗
 */
int ModifiedCholeskyDecomp2(const vector< vector<double> > &A, vector< vector<double> > &L, vector<double> &d, int n)
{
    if(n <= 0) return 0;
 
    L[0][0] = A[0][0];
    d[0] = 1.0/L[0][0];
 
    for(int i = 1; i < n; ++i){
        for(int j = 0; j <= i; ++j){
            double lld = A[i][j];
            for(int k = 0; k < j; ++k){
                lld -= L[i][k]*L[j][k]*d[k];
            }
            L[i][j] = lld;
        }
        d[i] = 1.0/L[i][i];
    }
 
    return 1;
}

不完全コレスキー分解

実際の問題を線形システムでモデル化した場合,nが非常に大きい場合が多い. たとえば,数値流体力学において有限差分法を用い, 計算空間をls_cholesky.eq37.gifに分割した場合, 線形システムとして解かなければならない圧力のポアソン方程式の行列のサイズはls_cholesky.eq38.gifにもなる. この分解を修正コレスキー分解などで正確に計算しようとすると非常に時間がかかる.

一方,こういった問題にみられる行列の性質として,疎行列というものがある. 疎行列は対角成分付近にのみ値が有り,非対角成分のほとんどが0の行列である. 例としてあげた圧力のポアソン方程式の行列も対称な疎行列となる. このような疎行列に対し,その要素が0ならば,その位置のLの要素を0にするといった近似を置くことで, 高速に,かつ,疎行列の性質を保ったまま変換するのが, 不完全コレスキー分解(Incomplete Cholosky decomposition : IC法)である.

たとえば,係数行列Aの0要素の位置のLの要素を0とする不完全コレスキー分解の場合, 修正コレスキー分解の計算式(ls_cholesky.eq36.gifの場合),

ls_cholesky.eq39.gif

に対して,以下の条件付で計算を行う.

  • ls_cholesky.eq40.gifの時は,ls_cholesky.eq41.gif
  • 右辺に対応する要素ls_cholesky.eq42.gifが0ならば,ls_cholesky.eq43.gifも0

なお,この方法はあくまで近似値を求めているので,正確な分解ではないことに注意.

不完全コレスキー分解の実装

以下は不完全コレスキー分解のコード例である.

  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
/*!
 * 不完全コレスキー分解(incomplete Cholesky decomposition)
 *  - 対称行列A(n×n)を下三角行列(L:Lower triangular matrix)と対角行列の積(LDL^T)に分解する
 *  - l_ii = 1とした場合
 *  - L: i > jの要素が非ゼロで対角成分は1
 *  - 行列Aの値が0である要素に対応する部分を飛ばす
 * @param[in] A n×nの対称行列
 * @param[out] L 対角成分が1の下三角行列
 * @param[out] d 対角行列(対角成分のみ)
 * @param[in] n 行列の大きさ
 * @return 1:成功,0:失敗
 */
int IncompleteCholeskyDecomp(const vector< vector<double> > &A, vector< vector<double> > &L, vector<double> &d, int n)
{
    if(n <= 0) return 0;
 
    d[0] = A[0][0];
    L[0][0] = 1.0;
 
    for(int i = 1; i < n; ++i){
        // i < k の場合
        for(int j = 0; j < i; ++j){
            if(fabs(A[i][j]) < 1.0e-10) continue;
 
            double lld = A[i][j];
            for(int k = 0; k < j; ++k){
                lld -= L[i][k]*L[j][k]*d[k];
            }
            L[i][j] = (1.0/d[j])*lld;
        }
 
        // i == k の場合
        double ld = A[i][i];
        for(int k = 0; k < i; ++k){
            ld -= L[i][k]*L[i][k]*d[k];
        }
        d[i] = ld;
        L[i][i] = 1.0;
    }
 
    return 1;
}

以下はls_cholesky.eq36.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
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
/*!
 * 不完全コレスキー分解(incomplete Cholesky decomposition)
 *  - 対称行列A(n×n)を下三角行列(L:Lower triangular matrix)と対角行列の積(LDL^T)に分解する
 *  - l_ii * d_i = 1とした場合
 *  - L: i > jの要素が非ゼロで対角成分は1
 *  - 行列Aの値が0である要素に対応する部分を飛ばす
 * @param[in] A n×nの対称行列
 * @param[out] L 対角成分が1の下三角行列
 * @param[out] d 対角行列(対角成分のみ)
 * @param[in] n 行列の大きさ
 * @return 1:成功,0:失敗
 */
int IncompleteCholeskyDecomp2(const vector< vector<double> > &A, vector< vector<double> > &L, vector<double> &d, int n)
{
    if(n <= 0) return 0;
 
    L[0][0] = A[0][0];
    d[0] = 1.0/L[0][0];
 
    for(int i = 1; i < n; ++i){
        for(int j = 0; j <= i; ++j){
            if(fabs(A[i][j]) < 1.0e-10) continue;
 
            double lld = A[i][j];
            for(int k = 0; k < j; ++k){
                lld -= L[i][k]*L[j][k]*d[k];
            }
            L[i][j] = lld;
        }
 
        d[i] = 1.0/L[i][i];
    }
 
    return 1;
}

添付ファイル: filels_cholesky.eq40.gif 433件 [詳細] filels_cholesky.eq5.gif 327件 [詳細] filels_cholesky.eq13.gif 344件 [詳細] filels_cholesky.eq21.gif 337件 [詳細] filels_cholesky.eq31.gif 284件 [詳細] filels_cholesky.eq36.gif 274件 [詳細] filels_cholesky.eq7.gif 341件 [詳細] filels_cholesky.eq30.gif 336件 [詳細] filels_cholesky.eq17.gif 317件 [詳細] filels_cholesky.eq24.gif 345件 [詳細] filels_cholesky.eq12.gif 359件 [詳細] filels_cholesky.eq19.gif 410件 [詳細] filels_cholesky.eq2.gif 338件 [詳細] filels_cholesky.eq14.gif 390件 [詳細] filels_cholesky.eq33.gif 329件 [詳細] filels_cholesky.eq43.gif 401件 [詳細] filels_cholesky.eq10.gif 358件 [詳細] filels_cholesky.eq35.gif 327件 [詳細] filels_cholesky.eq39.gif 383件 [詳細] filels_cholesky.eq20.gif 410件 [詳細] filels_cholesky.eq8.gif 316件 [詳細] filels_cholesky.eq6.gif 338件 [詳細] filels_cholesky.eq25.gif 355件 [詳細] filels_cholesky.eq29.gif 356件 [詳細] filels_cholesky.eq34.gif 322件 [詳細] filels_cholesky.eq37.gif 334件 [詳細] filels_cholesky.eq11.gif 955件 [詳細] filels_cholesky.eq16.gif 412件 [詳細] filels_cholesky.eq28.gif 401件 [詳細] filels_cholesky.eq1.gif 352件 [詳細] filels_cholesky.eq41.gif 375件 [詳細] filels_cholesky.eq22.gif 397件 [詳細] filels_cholesky.eq23.gif 328件 [詳細] filels_cholesky.eq4.gif 347件 [詳細] filels_cholesky.eq9.gif 321件 [詳細] filels_cholesky.eq27.gif 380件 [詳細] filels_cholesky.eq26.gif 390件 [詳細] filels_cholesky.eq18.gif 456件 [詳細] filels_cholesky.eq3.gif 297件 [詳細] filels_cholesky.eq42.gif 451件 [詳細] filels_cholesky.eq32.gif 430件 [詳細] filels_cholesky.eq38.gif 418件 [詳細] filels_cholesky.eq15.gif 395件 [詳細]

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