修正Cholesky分解

[←]　[↑]　[→]

---

修正Cholesky分解のアルゴリズム A=R^†DR 変数を節約するため対角rii=diとする。さらに作業領域として下三角rji=sij (i<j) と使ってもよい。
【外積形式】 /* Aの上三角をRに入れる */ for i=1 to n   for j=1 to n      if(i<=j) rij = aij      if(i>j )  rij = 0 /* Rの計算 */ for k=1 to n   for i=k+1 to n      for j=i to n          rij = rij - rki^*rkj/rkk   for j=k+1 to n     rkj = rkj/rkk	【内積形式、行ごと】 /* Aの上三角をRに入れる */ for i=1 to n   for j=1 to n      if(i==j) rij = aij      if(i<j)  sij = aij      if(i>j)  rij = 0 /* Rの計算 */ for k=1 to n     for i=1 to k-1         rkk = rkk - rik^*sik     for j=k+1 to n         for i=1 to k-1             skj = skj - rik^*sij         rkj = skj/rkk	【内積形式、列ごと】 /* Aの上三角をRに入れる */ for i=1 to n   for j=1 to n      if(i==j) rij = aij      if(i<j)  sij = aij      if(i>j)  rij = 0 /* Rの計算 */ for k=1 to n     for i=1 to k-1         for j=1 to i-1             sik = sik - rji^*sjk         rik = sik/rii         rkk = rkk - rik^*sik
連立1次方程式Ax=b を解くには /* R^†y=bをyについて解く(前進消去) */ for i=1 to n    tmp = 0    for j=1 to i-1       tmp = tmp + rji^* yj    yi = bi - tmp /* Rx=D^-1yをxについて解く(後退代入) */ for i=n to 1    tmp = 0    for j=i+1 to n       tmp = tmp + rij xj    xi = yi / rii - tmp

正定値のエルミート行列Aを三角行列に分解するアルゴリズムを導く。

Aの対角要素a_ii はエルミート行列だから実数であり

さらに正定値だからa_ii>0であることに注意。

◇

まずn行n列のAを1行1列目とそれ以外に分けて次のように書く

対角要素a₁₁の下を消去すると

さらにa₁₁の右側も消去できて、Aは3つの行列に分解できる。

ここで

このn-1行n-1列の行列も正定値のエルミート行列になっている。なぜなら

の右辺は明らかに正定値のエルミート行列だからである。

したがって同様に分解を続ければ、対角行列Dと対角要素が1である下三角R^†と上三角Rに分解される。これを修正コレスキー分解という。

A=R^†DR

具体的には、次の手順になる。ただし変数を節約するためr_ii=d_i とする。

・Aの上三角を作業領域Rに入れる。 ・ステップk (k=1 to n)において 　k+1行目から下のr を求める 　　rij = rij - rki^*rkj/rkk　(i=k+1 to n, j=i to n) 　 (LU分解のj=k+1 to nに比べて半分で済む) 　k行目のrを求める 　　rkj = rkj/rkk (i=k+1 to n)

【外積形式】 /* Aの上三角をRに入れる */ for i=1 to n     for j=1 to i-1         rij = 0     for j=i to n         rij = aij /* Rの計算 */ for k=1 to n     for i=k+1 to n         for j=i to n             rij = rij - rki^*rkj/rkk     for j=k+1 to n         rkj = rkj/rkk

◇◇

上記のアルゴリズムは外積aa^†を使うので外積形式とよばれるが、内積形式のアルゴリズムもある。

アルゴリズムを導く方針は、A=R^†DRと分解できたとして要素ごとに式を書き下し、これをd,rについて解く。

簡単のためn=3の場合を書き下すと

d,rについて解く手順として、行ごとか列ごとかの2とおりがある。

◇◇◇

行ごとに解くと

1行目　d₁=a₁₁ , r₁₂=a₁₂/d₁ , r₁₃=a₁₃/d₁

2行目　d₂=a₂₂ - d₁ r₁₂^*r₁₂ , r₂₃=(a₂₃ - d₁ r₁₂^*r₁₃)/d₂

3行目　d₃=a₃₃ - d₁ r ₁₃^*r₁₃ - d₂ r₂₃^*r₂₃

この手順を一般化すると次のようになる。ただし変数を節約するためr_ii=d_i とする。

・Aの上三角を作業領域Rに入れる。 ・ステップk (k=1 to n)において k行目のrkk を求める k行目のrkj を求める ここでriiで割ってから再びかけるのは無駄なので 中間変数sij=riirij を使うとよい

【内積形式,行ごと】 /* Aの上三角をRに入れる */ for i=1 to n     for j=1 to i-1         rij = 0     rii = aii     for j=i+1 to n         sij = aij /* Rの計算 */ for k=1 to n     for i=1 to k-1         rkk = rkk - rik^*sik     for j=k+1 to n         for i=1 to k-1             skj = skj - rik^*sij         rkj = skj/rkk

◇◇◇◇

列ごとに解くと

1列目　d₁=a₁₁

2列目　r₁₂=a₁₂/d₁ ,d₂=a₂₂ - d₁ r₁₂^*r₁₂

3列目　r₁₃=a₁₃/d₁ , r₂₃=(a₂₃ - d₁ r₁₂^*r₁₃)/d₂ , d₃=a₃₃ - d₁ r ₁₃^*r₁₃ - d₂ r₂₃^*r₂₃

この手順を一般化すると次のようになる。ただし変数を節約するためr_ii=d_i とする。

・Aの上三角を作業領域Rに入れる。 ・ステップk (k=1 to n)において k列目のrik を求める k列目のrkk を求める ここでrii で割ってから再びかけるのは無駄なので 中間変数sij=riirij を使うとよい

【内積形式,列ごと】 /* Aの上三角をRに入れる */ for i=1 to n     for j=1 to i-1         rij = 0     rii = aii     for j=i+1 to n         sij = aij /* Rの計算 */ for k=1 to n     for i=1 to k-1         for j=1 to i-1             sik = sik - rji^*sjk         rik = sik/rii          rkk = rkk - rik^*sik

◇◇◇◇◇

上記のアルゴリズムを検証するために使ったCコードは次のとおり。変数s_ijは下三角r_jiを使っている。

/* choleskey.cpp */

#include <iostream>

#include <complex>

using namespace std;

static const int N = 3;     // 行列のサイズ

typedef complex<float> cpx_t;      // 複素数型

/* 行列を表示するテンプレート */

template <int M,int N>

void print_matrix(cpx_t a[M][N]) {

    for(int i=0; i<M; i++) {

        cout << "    ";

        for(int j=0; j<N; j++) { cout << a[i][j] << "\t"; }

        cout << endl;

    }

}

/* ベクトルを表示するテンプレート */

template <int M>

void print_vector(cpx_t a[M]) {

    for(int i=0; i<M; i++) { cout << "    " << a[i] << endl; }

}

/***************************************************************

 * サブルーチン - 修正コレスキー分解,外積形式 (R)=(U_conj)(D)(U)

***************************************************************/

void cholesky_outer(cpx_t r[N][N]) {

    for(int k=0; k<N; k++) {

        for(int i=k+1; i<N; i++) {

            for(int j=i; j<N; j++) {

                r[i][j] -= conj(r[k][i]) * r[k][j] / r[k][k];

            }

        }

        for(int i=k+1; i<N; i++) {

            r[k][i] = r[k][i] / r[k][k];

        }

    }

}

/***************************************************************

 * サブルーチン - 修正コレスキー分解,内積形式,行ごと (R)=(U_conj)(D)(U)

***************************************************************/

void cholesky_inner_kj(cpx_t r[N][N]) {

    for(int k=0; k<N; k++) {

        for(int i=0; i<k; i++) {

            r[k][k] -= conj(r[i][k]) * r[k][i];

        }

        for(int j=k+1; j<N; j++) {

            for(int i=0; i<k; i++) {

                r[j][k] -= conj(r[i][k]) * r[j][i];

            }

            r[k][j] = r[j][k] / r[k][k];

        }

    }

}

/***************************************************************

 * サブルーチン - 修正コレスキー分解,内積形式,列ごと (R)=(U_conj)(D)(U)

 ***************************************************************/

void cholesky_inner_ik(cpx_t r[N][N]) {

    for(int k=1; k<N; k++) { // k=0のときは何もしないのでk=1から始める

        for(int i=0; i<k; i++) {

            for(int j=0; j<i; j++) {

                r[k][i] -= conj(r[j][i]) * r[k][j];

            }

            r[i][k] = r[k][i] / r[i][i];

            r[k][k] -= conj(r[i][k]) * r[k][i];

        }

    }

}

/***************************************************************

 * サブルーチン - 前進消去 (U_conj)y=bをyについて解く

***************************************************************/

void forward_elimination(cpx_t r[N][N],cpx_t b[N],cpx_t y[N]) {

    cpx_t tmp;

    for(int i=0; i<N; i++) {

        tmp = 0;

        for(int j=0; j<i; j++) {

            tmp += conj(r[j][i]) * y[j];

        }

        y[i] = b[i] - tmp;

    }

}

/***************************************************************

 * サブルーチン - 後退代入 (D)(U)x=yをxについて解く

***************************************************************/

void back_substitution(cpx_t r[N][N],cpx_t y[N],cpx_t x[N]) {

    cpx_t tmp;

    for(int i=N-1; i>=0; i--) {

        tmp = 0;

        for(int j=i+1; j<N; j++) {

            tmp += r[i][j] * x[j];

        }

        x[i] = y[i] / r[i][i] - tmp;

    }

}

/***************************************************************

 * メインルーチン - (A)x = bをxについて解く

***************************************************************/

int main(int argc, char *argv[]) {

    cpx_t r[N][N];

    cpx_t x[N];

    cpx_t y[N];

    cpx_t Ax[N];

/* 行列(A)とbに値を設定 */

    cpx_t a[N][N] =

       { {  2   ,  0+1i,  1    }

       , {  0-1i,  3   ,  0-1i }

       , {  1   ,  0+1i,  4    } };

    cpx_t b[N] =

       { 1+1i

       , 2-1i

       , 3+1i };

    cout << "-- A=" << endl;

    print_matrix<N,N>(a);

    cout << "-- b=" << endl;

    print_vector<N>(b);

/* 行列(R)に(A)を代入 */

    for(int i=0; i<N; i++) {

        for(int j=0; j<N; j++) {

#ifdef CHOLESKY_OUTER

// 外積形式のとき上三角に入れて,結果も上三角に出る,下三角は不使用

            r[i][j] = (i<=j)? a[i][j] : 999;

#else

// 内積形式のとき下三角に入れて,結果は上三角に出る

            r[j][i] = (i<=j)? a[i][j] : 999;

#endif

        }

    }

    cout << "-- initial R=" << endl;

    print_matrix<N,N>(r);

/* 修正コレスキー分解 (R)=(U_conj)(D)(U) */

#ifdef CHOLESKY_OUTER

    cholesky_outer(r);      // 外積形式

#elif defined CHOLESKY_INNER_KJ

    cholesky_inner_kj(r);   // 内積形式,行ごと

#else

    cholesky_inner_ik(r);   // 内積形式,列ごと

#endif

/* 前進消去 (U_conj)y=bをyについて解く */

    forward_elimination(r,b,y);

/* 後退代入 (D)(U)x=yをxについて解く */

    back_substitution(r,y,x);

/* 分解後の行列(R)とxを表示 */

    cout << "-- R=" << endl;

    print_matrix<N,N>(r);

    cout << "-- x=" << endl;

    print_vector<N>(x);

/* 検算 (A)x=b */

    for(int i=0; i<N; i++) {

        Ax[i] = 0;

        for(int j=0; j<N; j++) {

            Ax[i] += a[i][j] * x[j];

        }

    }

    cout << "-- Ax=" << endl;

    print_vector<N>(Ax);

}

CHOLESKY_INNER_KJの場合

-- A=

    (2,0)     (0,1)  (1,0) 

    (0,-1)    (3,0)  (0,-1)

    (1,0)     (0,1)  (4,0) 

-- p=

    (1,1)

    (2,-1)

    (3,1)

-- initial R=

    (2,0)     (999,0)       (999,0)      

    (0,1)     (3,0)  (999,0)      

    (1,0)     (0,-1) (4,0) 

-- R=

    (2,0)     (0,0.5)       (0.5,0)      

    (0,1)     (2.5,0)       (0,-0.2)     

    (1,0)     (0,-0.5)      (3.4,0)      

-- x=

    (0.117647,0.176471)

    (0.588235,-0.0588235)

    (0.705882,0.0588235)

-- Ax=

    (1,1)

    (2,-1)

    (3,1)

---

[←]　[↑]　[→]