Cholesky分解

今回は特定の条件を満たした時に使えるLU分解として、Cholesky(コレスキー)分解について書きます。特定の条件というのは、方程式Ax=bについて、係数行列Aが正定値対称行列である、というものです。対称行列というのは転置を行っても同じ行列になる、つまりとなるものです。正定値行列というのは行列と同じサイズの任意のベクトルxを持ってきて常に、となる行列Aを言います。このような行列はすべての固有値が正であるなど特殊な性質を持っています。このような正定値対称行列に対して例えばLU分解を行ってあげると、となります。つまり、Lが求まれば自動的にUも求められるということになります。LU分解ではLとUのそれぞれを計算して求めていましたが、このような特殊な性質を用いてあげることによって、計算量を半分近くまで減らすことが可能となります。これがコレスキー分解の目標となります。

コレスキー分解の流れとプログラム

実質的にやっていることはLU分解とほとんど変わらず、行列の左上からひたすら値の更新を行っていくというものになります。
①対象となる行を選択
②その行の対角成分に着目
③分解後の対角成分を求める
④対角成分と同じ列の要素を更新する(下三角行列Lにおいて、着目している対角成分の下側を更新する<=>上三角行列Uにおいて、着目している対角成分の右側を更新する)
⑤①から④を繰り返す

③や④の操作に関しては、対称行列であるという特徴を活かすと、うまく計算量を削減することが可能になります。

プログラムは下のようになります。

void cholesky_decompose(real *dmat, int M, int N) {

  int i, j, k;
  real vii, vji;

  for (i = 0; i < M; i++) {
    vii = dmat[i * N + i];
    for (j = 0; j < i; j++) {
      vii -= dmat[i * N + j] * dmat[i * N + j];
    }
    dmat[i * N + i] = sqrt(vii);

    for (j = i + 1; j < M; j++) {
      vji = dmat[j * N + i];
      for (k = 0; k < i - 1; k++) {
        vji -= dmat[j * N + k] * dmat[i * N + k];
      }
      dmat[j * N + i] = vji / dmat[i * N + i];
    }
  }

}

計算式にルート計算が入るので、計算のコストが大きいというのはあります。また、正定値性を持つ行列なら問題は無いのですが、ルートの中が正でなければ実数空間での処理は出来ないのでその点は問題です。
通常のコレスキー分解を改良することで、ルートを使わず、かつ対象範囲を対称行列(正定値性は問わない)まで拡張した修正コレスキー分解というのも存在して、こちらのほうがよく使われているらしいです。ざっくり調べた感じでは、修正コレスキー分解では対角成分だけ抜き出して考えるというもののようです。つまり、対角行列Dを持ってきて、

と分解するというものです。従来のコレスキー分解同様Lに関してはLU分解の半分の計算量ですみ、対角成分の計算に関してはルート計算が発生しないという強みがあります。

ここまではLU分解とCholesky分解という二つの代表的な連立一次方程式の解法を見てきました。この二つはどちらも直接法と呼ばれる方法に分類されていて、行列を直接操作することによって、解を厳密に求めていくというものになります。これに対して、反復法と呼ばれる異なるアプローチで解を求めにいく方法があります。今後はそちらの方を見ていく予定です。