コレスキー分解

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

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

上式を成分で表す． を，

とすると，

となる．の場合の要素を分けて書くと，

よって，以下の式をと計算していくことでLを求めることができる．

たとえば，

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

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

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

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

修正コレスキー分解†

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

まず，コレスキー分解の下三角行列を対角成分が1の下三角行列と対角行列の積に分解する．

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

ここで，，とおくと，

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

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

コレスキー分解の時と同様にの要素を分けて書くと，

ここで，なので，

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

ただし，である．

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

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

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

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

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

修正コレスキー分解をC++で実装した例を以下に示す． それぞれ，の場合，の場合の実装である．

/*!
 * 修正コレスキー分解(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;
}

/*!
 * 修正コレスキー分解(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が非常に大きい場合が多い． たとえば，数値流体力学において有限差分法を用い， 計算空間をに分割した場合， 線形システムとして解かなければならない圧力のポアソン方程式の行列のサイズはにもなる． この分解を修正コレスキー分解などで正確に計算しようとすると非常に時間がかかる．

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

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

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

• の時は，
• 右辺に対応する要素が0ならば，も0

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

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

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

/*!
 * 不完全コレスキー分解(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;
}

以下はとした場合．

/*!
 * 不完全コレスキー分解(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;
}