数値解析ノート2. Gauss-Seidelの定常反復法

### はじめに

先日，$\mathbf{A}\mathbf{x}=\mathbf{b}$ という形の連立一次方程式の反復解法で，特に**定常反復法**というクラスに属する方法のうち最も単純なJacobi法と呼ばれるものを扱った:
[数値解析ノート1. Jacobiの定常反復法](https://qiita.com/tarotene/items/4099eb9eae66695eb5f4)

このJacobi法は，素直な拡張によってGauss-Seidel法に変形することができる．実際の導出については次の記事をご覧いただきたい．
[数値解析ノート: 連立1次方程式の定常反復解法(2). Gauss-Seidel法](https://trtn.hatenablog.com/entry/2020/03/16/191918)

本記事では，Gauss-Seidel法の実装とその実験結果，そしてJacobi法との比較を記す．

---

### Gauss-Seidel法の実装と実験結果

ひとつ注意しておくと，Gauss-Seidel法はJacobi法の実装の一つの単純化なのに，収束性が良いのである．よって，これを定常反復法の出発点と捉えても良い．

まず，Jacobi法をC++で実装した時のコア部分は以下の通りである．

```cpp
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` を用意しなくても良い分，いくぶんシンプルなコードになる．

```cpp
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}$$
とすることで更新直後のベクトルの成分を即座に使えるようにする

という処方である．言ってみれば一つの**貪欲法**である．

そこで，[前回記事](https://qiita.com/tarotene/items/4099eb9eae66695eb5f4)の慣習通りにソースコードをコンパイル・実行すると，誤差 $\|\mathbf{A}\mathbf{x}_{k} - \mathbf{b}\|^{2}$ の挙動は次のようになる．
![002-Gauss-Seidel-residue-Poisson.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/93579/250c5660-6110-33a2-781d-e3ab2fd4a4a4.png)

比較のためにJacobi法の場合も再掲しておく．
![001-Jacobi-residue-Poisson.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/93579/9cbcc51c-844f-afff-3223-deefc3aeda82.png)

確かに，Gauss-Seidel方のほうが収束性は若干良い[^comp]．
[^comp]: この計算結果がインチキなのではないかというツッコミがあるかも知れないので，念の為に述べておくと，全く同じ行列 $\mathbf{A}$ と全く同じ統計性に従って生成されたランダムベクトル $\mathbf{b}$ に対して，Jacobi方とGauss-Seidel法のパフォーマンスの比較を行っている．行列とベクトルを生成する関数はヘッダ `iter.hpp` にまとめてあるので，興味のある人は[リポジトリ](https://github.com/tarotene/iterative-methods)に見に行って下さい．

### 最後に

今回の数値実験も，次のGithubリポジトリに含まれているので，その気になれば追試できると思う．
https://github.com/tarotene/iterative-methods