LU分解

前回はガウスの消去法について書きましたが、今回はそこから一歩進めた方法であるLU分解について記します。
簡単に説明するとある方程式Ax=bの係数行列Aについて、下三角行列Lと上三角行列Uを用いて

という形で分解してあげるのがLU分解となります。あとは方程式のほうに戻ってあげるとLUx=bとなりますので、まずUx=cとでも置いて、Lc=bを解きます。cが求まった後にUx=cを解くと、無事にxを求めることが出来ます。Lは下三角行列なので、Lc=bを解く操作を前進代入、Uは上三角行列なので、Ux=cを解く操作を後進代入といいます。このように分解をすると何が嬉しいのかについてですが、ガウスの消去法では方程式の右辺ベクトルを含めて求解を進めていきます。そのため、行列が同じでも右辺ベクトルが異なるような方程式が現れた場合にはもう一度最初から解く必要が出てきます。このような行列が同じで右辺ベクトルだけが異なるというようなシチュエーションは多いらしく、ガウスの消去法では非常にコストが大きいという問題があります。LU分解を用いて予め行列Aを分解しておくことによって、右辺ベクトルが変わった際に前進代入と後進代入をさくっと行うだけで解を求めることができるようになります。そのため、方程式を解く上で非常に有用な手段であるといえます。また、この分解操作をどのように行うかはこの後書いていく予定ですが、この分解操作をうまくやれば行列行列積を使って行うことが出来るようです。行列行列積はメモリアクセス量に対して演算が多く、キャッシュで再利用させるなどうまく処理してあげることによって高い性能を出すことが可能な演算です。なので、単純に入れ替えだけを行う操作に比べて高い性能を出すことがLU分解では出来るようです。

LU分解のアルゴリズム

LU分解のアルゴリズムとして、right-looking、left-looking、Crout法の3つがあります。ここではそれぞれどのようなアルゴリズム化を見ていきたいと思います。併せてプログラムもここで記します。

Right-looking

大まかな流れとしては
①対象となる行を決める
②対象となる行の対角成分を見る
③対角成分をもとに、その成分と同じ列の要素を更新(なお対象となるのはその成分より下の行の要素)
④対象となる行の要素を用いて、対角成分より右下の要素を更新していく
⑤この操作をすべての行に対して行う

対象となる対角成分の右部分を参照しながら右下部分の更新を行うことからright-lookingと呼ばれています。なお、Lの対角成分が1となることを仮定しています。

プログラムで表すと下のような形になります。ちなみに行列はMxNとします。

void lu_decompose_right(real *dmat, int M, int N) {
  
  int i, j, k;
  real temp, vij;

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

Left-looking

Right-lookingでは右側を見て右側を更新していたのに対して、left-lookingでは左側を見て更新を行っていきます。これもLの対角成分が1になることを仮定しています。
こちらもおおまかな流れを書くと
①対象となる行を決める
②その行の対角成分を見る
③その対角成分と同じ列の上の部分と左側の部分を元にして上の部分を更新する
④最後に対角成分と同じ列の下の部分を対角成分を元にして更新
⑤これを繰り返す

④はright-lookingと同じですが、①～③の流れは完全に別物になっています。

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

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

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

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

今回は行方向に連続になるように行列のメモリ配置を取っていましたが、列方向に連続になるように取ってあげれば、更新の部分がメモリ連続になって良さそうだなあとか思ったり。

Crout法

Crout法では右だけとか左だけとではなく、様々なところを参照しながら更新を行っていきます。また、Crout法ではUの対角成分を1とします。
流れとしては
①対象となる行を決める
②その行の対角成分を見つける
③その対角成分の左下と右上を元に、対角成分と同じ行と列の更新を行う
④これをすべての行について繰り返す

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

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

  int i, j, k;
  real vii;  

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

Crout法では更新する際に2箇所の参照領域が必要となります。そのため、行列全体を分割して並列に処理を行うような場合には、参照領域を求めてプロセス間で通信が発生してしまいます。そのため、並列計算にはあまり向いてはいません。ただ、ノード内でメモリが共有されているような状態でしたら問題はありません。

ここまでで3つのアルゴリズムの紹介が終わりました。どうやら並列計算とかを念頭に置くと、right-lookingが最も良いらしいですね。あと、最初に行列行列積がうんぬんと言いましたが、どうやらブロック形式のアルゴリズムで用いられるらしいです。今回は力尽きたのでここまでとしますが、やる気が出た時にでもブロック形式は勉強したいなあというところです。