LU分解†

これまでn元連立1次方程式

ls_lu.eq1.gif

を直接法，反復法などで解くことだけを考えてきた． それでは実際の問題ではこのような線形システムはどのように扱われるだろうか． たとえば，制御の分野ではシステムの状態変化を捉えるために用いている． 係数行列Aがシステム内部を表し，右辺の定数ベクトルbで外乱などの状態を表す． システム内部が変わらず，外の状態が様々に変化したときの状態を知りたいとき， 係数行列Aは変化せず，bのみが変わるだけである． このように係数行列が固定で右辺の定数ベクトルだけが変化するということは， 物理学などの他の分野でも多くある． このとき，係数行列Aを解きやすい形に分解しておけば，計算量を大幅に減らすことができる． ここではそのような分解の一つであるLU分解について述べる．

n元連立1次方程式の係数行列を考える．

ls_lu.eq2.gif

これを以下の下三角行列(lower triangular matrix) L と 上三角行列 (upper triangular matrix) U に分解する．

ls_lu.eq3.gif
ls_lu.eq4.gif

ここで，ls_lu.eq5.gifである． このように分解することをLU分解と呼ぶ．

LU分解した行列で連立1次方程式を解くことを考える． まず，ls_lu.eq5.gifより，

ls_lu.eq6.gif

ls_lu.eq7.gifとすると，

ls_lu.eq8.gif

となる．Lは下三角行列であるので，1行目から順番に代入していくことで， 容易にls_lu.eq9.gifを算出することができる．この処理を前進代入(forward substitution)と呼ぶ． ls_lu.eq9.gifが求まったらそれを2番目の式に代入する． Uは上三角行列であるので，ガウスの消去法と同じく後退代入(backward substitution) していくことで解ls_lu.eq10.gifを求めることができる．

LU分解の手順†

係数行列AをLU分解するためには，以下の式を1行目から順番に適用していく．

ls_lu.eq11.gif

メモリを節約するために，LU分解した結果を以下のように1つの行列に格納する．

ls_lu.eq12.gif
↑

LU分解の実装†

LU分解をC++で実装した例を以下に示す．

/*!
 * LU分解(ピボット交換なし)
 *  - 行列A(n×n)を下三角行列(L:Lower triangular matrix)と上三角行列(U:Upper triangular matrix)に分解する
 *  - L: i >= j,  U: i < j の要素が非ゼロでUの対角成分は1
 *  - LとUを一つの行列にまとめた形で結果を返す
 * @param[inout] A n×nの係数行列．LU分解した結果を格納する．
 * @param[in] n 行列の大きさ
 * @return 1:成功,0:失敗
 */
int LUDecomp(vector< vector<double> > &A, int n)
{
    if(n <= 0) return 0;

    for(int i = 0; i < n; ++i){
        // l_ijの計算(i >= j)
        for(int j = 0; j <= i; ++j){
            double lu = A[i][j];
            for(int k = 0; k < j; ++k){
                lu -= A[i][k]*A[k][j];    // l_ik * u_kj
            }
            A[i][j] = lu;
        }

        // u_ijの計算(i < j)
        for(int j = i+1; j < n; ++j){
            double lu = A[i][j];
            for(int k = 0; k < i; ++k){
                lu -= A[i][k]*A[k][j];    // l_ik * u_kj
            }
            A[i][j] = lu/A[i][i];
        }
    }

    return 1;
}
↑

LU分解で連立1次方程式を解く手順†

LU分解で連立1次方程式を解く方法は上でも述べたがここではより具体的な手順について説明する． LU分解で連立1次方程式を解く手順は以下である．

  1. 下三角行列Lと上三角行列Uを求める．
  2. 前進代入によりls_lu.eq9.gifを求める．
  3. 後退代入によりls_lu.eq10.gifを求める． 1はすでに述べているので2,3について説明する．
  • 前進代入 下三角行列Lと定数ベクトルbからyを求める．
    ls_lu.eq13.gif
    Lは下三角行列であるので，1行目は単純にls_lu.eq14.gifでls_lu.eq15.gifが算出される． 得られたls_lu.eq15.gifを2行目の式に代入して同様に計算するとls_lu.eq16.gifが得られる．これをn行目まで繰り返す．式で書くと以下となる．
    ls_lu.eq17.gif
    ここで，ls_lu.eq18.gifである．
  • 後退代入 前進代入で得られたyと上三角行列Uからxを求める．
    ls_lu.eq19.gif
    基本的にはガウスの消去法の後退代入と同じである． ただし，各行は対角項で正規化されている(対角成分がすべて1)． 前進代入とは逆にn行目から順番にls_lu.eq20.gifを求めていく．
    ls_lu.eq21.gif
    ここで，ls_lu.eq18.gifである．
↑

LU分解を用いた連立1次方程式の解法の実装†

LU分解で連立1次方程式を解くコードをC++で実装した例を以下に示す．

/*!
 * LU分解した行列a(n×n)から前進代入・後退代入によりA・x=bを解く
 * @param[in] A LU分解された行列
 * @param[in] b 右辺ベクトル
 * @param[out] x 結果ベクトル
 * @param[in] n 行列の大きさ
 * @return 1:成功,0:失敗
 */
int LUSolver(const vector< vector<double> > &A, const vector<double> &b, vector<double> &x, int n)
{
    if(n <= 0) return 0;

    // 前進代入(forward substitution)
    //  LY=bからYを計算
    for(int i = 0; i < n; ++i){
        double bly = b[i];
        for(int j = 0; j < i; ++j){
            bly -= A[i][j]*x[j];
        }
        x[i] = bly/A[i][i];
    }

    // 後退代入(back substitution)
    //  UX=YからXを計算
    for(int i = n-1; i >= 0; --i){
        double yux = x[i];
        for(int j = i+1; j < n; ++j){
            yux -= A[i][j]*x[j];
        }
        x[i] = yux;
    }

    return 1;
}
↑

コレスキー分解†

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.gif，ls_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の場合の実装である．

/*!
 * 修正コレスキー分解(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が非常に大きい場合が多い． たとえば，数値流体力学において有限差分法を用い， 計算空間を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

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

↑

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

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

/*!
 * 不完全コレスキー分解(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とした場合．

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

参考文献†

  • 佐藤次男, 中村理一郎, “よくわかる数値計算 アルゴリズムと誤差解析の実際”, 日刊工業新聞社, 2001.
  • 川上一郎, “数値計算の基礎”, http://www7.ocn.ne.jp/~kawa1/
トップ   編集 凍結 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2024-03-08 (金) 18:06:11