はじめに
先日，$\mathbf{A}\mathbf{x}=\mathbf{b}$ という形の連立一次方程式の反復解法で，特に定常反復法というクラスに属する方法のうち最も単純なJacobi法と呼ばれるものを扱った:
数値解析ノート1. Jacobiの定常反復法
このJacobi法は，素直な拡張によってGauss-Seidel法に変形することができる．実際の導出については次の記事をご覧いただきたい．
数値解析ノート: 連立1次方程式の定常反復解法(2). Gauss-Seidel法
本記事では，Gauss-Seidel法の実装とその実験結果，そしてJacobi法との比較を記す．
Gauss-Seidel法の実装と実験結果
ひとつ注意しておくと，Gauss-Seidel法はJacobi法の実装の一つの単純化なのに，収束性が良いのである．よって，これを定常反復法の出発点と捉えても良い．
まず，Jacobi法をC++で実装した時のコア部分は以下の通りである．
for (size_t k = 0; k < ITERS; k++)
{
double y[N] = {0.0};
for (size_t i = 0; i < N; i++)
{
y[i] = x[i];
}
for (size_t i = 0; i < N; i++)
{
for (size_t j = 0; j < i; j++)
{
x[i] -= A[i * N + j] * y[j]; /* use old vector x^(k - 1) */
}
for (size_t j = i + 1; j < N; j++)
{
x[i] -= A[i * N + j] * y[j]; /* use old vector x^(k - 1) */
}
x[i] /= A[i * N + i];
}
/* calculate residue */
double Ax[N] = {0.0};
double res = residue_normal_axb(A, x, b, Ax, N);
}
一方，Gauss-Seidel法を同じくC++で実装すると，更新前のベクトル $\mathbf{x}^{(k)}$ を格納するための配列 y を用意しなくても良い分，いくぶんシンプルなコードになる．
for (size_t k = 0; k < ITERS; k++)
{
for (size_t i = 0; i < N; i++)
{
x[i] = b[i];
for (size_t j = 0; j < i; j++)
{
x[i] -= A[i * N + j] * x[j]; /* use already refreshed vector x^(k) */
}
for (size_t j = i + 1; j < N; j++)
{
x[i] -= A[i * N + j] * x[j]; /* use old vector x^(k - 1) */
}
x[i] /= A[i * N + i];
}
/* calculate residue */
double Ax[N] = {0.0};
double res = residue_normal_axb(A, x, b, Ax, N);
}
Gauss-Seidel法では具体的に何をやっているのかというと
Jacobi法では単に
$$x^{(k + 1)}_{i} = (b_{i} - \sum_{j\neq i}a_{i, j}x^{(k)}_{j}) / a_{i, i}$$
という規則に従ってベクトル $\mathbf{x}$ を更新していたのを，要素ごとの更新順序を $j$ の昇順に固定し
$$x^{(k + 1)}_{i} = (b_{i} - \sum_{j< i}a\_{i, j}x\^{(k + 1)}\_{j} - \sum\_{j> i}a_{i, j}x^{(k)}_{j}) / a_{i, i}$$
とすることで更新直後のベクトルの成分を即座に使えるようにする
という処方である．言ってみれば一つの貪欲法である．
そこで，前回記事の慣習通りにソースコードをコンパイル・実行すると，誤差 $|\mathbf{A}\mathbf{x}_{k} - \mathbf{b}|^{2}$ の挙動は次のようになる．
比較のためにJacobi法の場合も再掲しておく．
確かに，Gauss-Seidel方のほうが収束性は若干良い1．
最後に
今回の数値実験も，次のGithubリポジトリに含まれているので，その気になれば追試できると思う．
https://github.com/tarotene/iterative-methods