Jacobi法

ここまで連立一次方程式を解く手法としてガウスの消去法をベースとしたものを見てきましたが、ここからいくつかは反復法と呼ばれるものを見ていきたいと思います。LU分解などでは行列の値を直接操作することによって、解を求めていました。このような手法での計算量はO(N^3)であり、行列のサイズが大きくなると解くのが困難となってきます。また、係数行列が疎行列(行列の内部の要素のほとんどがゼロであるような行列)であるような場合には、直接法ではフィルイン(ゼロだった要素がゼロでなくなる、つまり疎行列であることの利点を活かしづらくなる)が発生してしまいます。反復法というのはある適当な初期解を設定した後、繰り返し計算を行うことによって、求めるべき解に近づけていくというものです。正しい解を得られるか、十分に解との差が小さくなるまで処理を行います。反復法には大きく分けて、定常反復法と非定常反復法の二つがあります。今回書くJacobi法は定常反復法の方になります。反復法として広く知られている共役勾配法は非定常反復法にカテゴライズされます。では早速Jacobi法について、見ていきたいと思います。

Jacobi法の流れ

Jacobi(ヤコビ)法の大まかなイメージとしては、繰り返し処理を行う中で、ひとつ前の繰り返しで得られた近似解をもとに新たな近似解を生成するというものです。
数式的に表すと、Ax=bを解く上で、まず係数行列Aを
A=L + D + U
という形で分解をします。なお、Lは下三角行列(対角成分もゼロ)、DはAの対角成分を並べた対角行列、Uは上三角行列(対角成分はゼロ)となります。そして、k回目の反復で得られた近似解x_kからk+1回目の近似解を次のように定義するのがヤコビ法です。
x^{k+1} = D^{-1}(-(L + U)x^k + b)
ここでDの逆行列についてですが、対角行列ですので、単純にDの各要素の逆数を取ったものがDの逆行列となります。数式のイメージとしてはxのi行目の値の更新を行う際には、ひとつ前の近似解のi行目以外のデータをもとにi行目の方程式を満たすように解を求めてあげるというようなものになります。繰り返し処理を行っていくと最終的に近似解が固まるようになります。つまり、]
x^* =  D^{-1}(-(L + U)x^* + b)
となります。
この式の両辺にDを掛けて項を移動すればAx^*=bとなり、x^*がこの方程式の解であることがわかります。
このように流れとしては非常に単純でわかりやすいものとなっています。

Jacobi法のプログラム

Jacobi法のプログラムが下のものになります。中身の処理自体に関しては上に書いた流れのとおりです。追加して書くべきこととしては、最大の反復回数を設定してあげることと、どこで反復を打ち切るのかについてです。最大の反復回数に関しては対象とする行列のサイズに併せて設定してあげればいいと思われます。どこで反復を打ち切るのか、つまり、近似解と解との誤差が十分に小さいと判断するのはどうやって行うかについてですが、ここでは新たに得られた近似解とその前の反復で得られた解について、差の絶対値をとったものを採用し、それが十分に小さいか(例えば、EPS=1.0e-5、EPSは精度を表すイプシロンの略)を判定し、小さければ反復終了、大きければ反復を継続というプログラムになっています。数式で表すと
|x^{k}-x^{k-1}| < \mathrm{EPS}
というような感じになります。

void solve_jacobi(real *dmat, real *x, real *y, int M, int N) {
  
  int i, j, k, it;
  real r; //residual
  real *x_old;
  x_old = (real *)malloc(sizeof(real) * N);

  /* Set initial solution */
  for (i = 0; i < N; i++) {
    x[i] = 0;
    x_old[i] = 0;
  }
  
  /* Iteration */
  for (it = 0; it < MAXIT; it++) {
    
    for (i = 0; i < M; i++) {
      x[i] = y[i];
      for (j = 0; j < N; j++) {
	if (i != j) {
	  x[i] -= dmat[i * M + j] * x_old[j];
	}
      }
      x[i] /= dmat[i * M + i];
    }

    r = 0;
    for (i = 0; i < N; i++) {
      r += fabs(x[i] - x_old[i]);
    }
    
    if (r < EPS) {
      break;
    }

    cout << it << ": ";
    print_vector(x, M);
    for (i = 0; i < N; i++) {
      x_old[i] = x[i];
    }

  }

}

ヤコビ法の利点としては各行の処理を独立して行えるので並列化を行いやすいというところですね。ただ収束までの反復数が多いというのが難点です。
この次は順番的に反復数の削減を可能にしたガウス・ザイデル法についてかな。