Gauss-Seidel法 - Slope Blog

前回はヤコビ法について書きましたが、今回はそこから一歩推し進めたGauss-Seidel(ガウス・ザイデル)法について書きます。ヤコビ法では収束までの反復数が多いというのが問題でした。これを解決するために考えだされたのがガウス・ザイデル法です。簡単に違いを述べると、ヤコビ法は一回前の反復で得られた近似解を使って新しいより良い近似解を作ろうというものでしたが、ガウス・ザイデル法は更新が終わった近似解の値、つまり新しい近似解と古い近似解の両方を使って更新を行っていくというものです。
Ax=bを解く上で、k反復目の近似解を $x^k$ とします。例えばk+1反復目の近似解のi行目を求めることを考えます。その際、1～i-1行目の近似解に関しては既に更新が終わっているのでk+1反復目のものを使い、i+1行目からN行目までの近似解に関してはまだ更新が終わってないのでk反復目までのものを用います。最も新しい近似解を随時使っていくことによって収束を早めることが可能となります。
プログラムは下のような形になります。

void solve_gauss_seidel(real *dmat, real *x, real *y, int M, int N) {
  
  int i, j, k, it;
  int M, N;
  real old_x, r; //residual

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

    if (r < EPS) {
      break;
    }

    cout << it << ": ";
    print_vector(x, M);

  }

}

前回のヤコビ法の実装では古い近似解の値を残しておきましたが、ガウス・ザイデル法では随時更新していくことから残して置く必要は無いため、消えています。手元にある適当なテストデータで実験してみたところ、収束までの反復数も削減されており、ヤコビ法と比較して収束のスピードが上がっていることを確認しました。ただ、近似解の行間に依存関係があるため、並列化が行いづらいという問題があります。

SOR法

がっつり書くつもりは無いですが、更に収束のスピードを上げる手法としてSOR(Successive Over-Relaxation)法(逐次緩和法)というものがあります。これはガウス・ザイデル法で得られた近似解をそのままその反復での解にするのではなく、その値とひとつ前の近似解を組み合わせたものを新たな近似解とする、というものです。k+1反復目でガウス・ザイデル法によって得られた解を $x^*_{k+1}$ とし、k反復目での近似解を $x_k$ とすると、新たな近似解をパラメータ $\omega$ を用いて
$x_{k+1}=(1-\omega)x_{k} + \omega x^*_{k+1}$
としたものが、SOR法です。
パラメータ $\omega$ は $0<\omega<2$ の範囲で設定しなければ収束はしないそうです。ただうまく設定してあげればガウス・ザイデル法よりも収束が高速化されます。ちなみに $\omega=1$ だと単なるガウス・ザイデル法になります。ここらへんのパラメータ調整的な話はきりが無い気がするので、SOR法についてこのぐらいで終わりにしたいと思います。