休憩室

科学・技術・芸術

数値解析ノート: 連立1次方程式の定常反復解法(2). Gauss-Seidel法

はじめに

本記事では,$\mathbf{A}\mathbf{x} = \mathbf{b}$ 型の連立一次方程式を解くための反復法であるGauss-Seidel法の導出と,その性質について述べる.

Gauss-Seidel法は,Jacobi法と並んで定常反復法という大きな枠組みに分類される方法である.ここで,定常反復法というのは,方程式 $\mathbf{A}\mathbf{x} = \mathbf{b}$ を満たす解の候補としてある初期ベクトル $\mathbf{x} = \mathbf{x}^{(0)}$ を考え,ある反復行列 $\mathbf{B}$ および反復ベクトル $\mathbf{c}$ を用いて $$\mathbf{x}^{(k + 1)} = \mathbf{B} \mathbf{x}^{(k)} + \mathbf{c}; \;k = 0, 1, 2, \ldots$$ という漸化式にしたがって解ベクトル $\mathbf{x}$ の更新を行い,真の解に収束させていく方法であった.「定常反復法」というネーミングは,更新に用いられる $(\mathbf{B}, \mathbf{c})$ が変化しないことを指している1. 定常反復法に関する一般的については,前回の記事で取り扱ったのでそちらも適当に参照して欲しい.

数値解析ノート: 連立1次方程式の定常反復解法(1). Jacobiの反復法 - 休憩室

Jacobi法は,反復行列と反復ベクトルをそれぞれ $$\mathbf{B} = - \mathbf{D}^{-1} (\mathbf{L} + \mathbf{U}), \;\mathbf{c} = \mathbf{D}^{-1} \mathbf{b}$$ ととった場合に相当するが,Gauss-Seidel法ではこれが $$\mathbf{B} = - (\mathbf{D} + \mathbf{L})^{-1} \mathbf{U}, \;\mathbf{c} = (\mathbf{D} + \mathbf{L})^{-1}\mathbf{b}$$ に変更される.


定常反復法としてのGauss-Seidel法

Jacobi法は,更新したいベクトル $\mathbf{x}^{(k)}$ の各成分 ${x^{(k)}_{j}}_{j = 1, 2, \ldots, N}$ が,$\mathbf{x}^{(k)}$ から $\mathbf{x}^{(k + 1)}$ への更新の途中で互いに干渉しないという性質のおかげで,素直な並列化ができた.この辺に関しては,実装と数値実験の結果を示した次の記事を参照して欲しい.

数値解析ノート1. Jacobiの反復法 - Qiita

一方で,Gauss-Seidel法はこの並列化を見かけ上犠牲にする代わりに,ほんの少し収束性を良くする仕組みを備えている2.また,収束性のみならず実装が自然になるというメリットもある.この辺りの実用上の事情は次の記事に記しておいた. qiita.com


Gauss-Seidel法のアルゴリズム

Gauss-Seidel法は,成分表示で見ると

$$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};\; k = 0, 1, 2, \ldots$$

という更新規則に基づくアルゴリズムである.一方,Jacobi法の更新規則は

$$x^{(k + 1)}_{i} = (b_{i} - \sum_{j < i} a_{i, j} x^{(k)}_{j} - \sum_{j > i} a_{i, j} x^{(k)}_{j}) / a_{i, i};\; k = 0, 1, 2, \ldots$$

で与えられており,これらはよく似ている.違いは更新したばかりのベクトルの成分を即座に用いるかどうかである.

このことを更に一般化すると,何らかの順列 ${p_{j}}_{j = 1, 2, \ldots, N}$ が与えられており,この順序に従って

$$x^{(k + 1)}_{p_{i}} = (b_{p_{i}} - \sum_{p_{j} < p_{i}} a_{p_{i}, p_{j}} x^{(k)}_{p_{j}} - \sum_{p_{j} > p_{i}} a_{p_{i}, p_{j}} x^{(k)}_{p_{j}}) / a_{p_{i}, p_{i}};\; k = 0, 1, 2, \ldots$$

という更新を施すのもまたGauss−Seidel法の資格を持つ.例えば本来の(素朴な)Gauss-Seidel法を逆順にした

$$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};\; k = 0, 1, 2, \ldots$$

という更新規則はその例である.

ちなみに,Gauss-Seidel法は逐次的過緩和(SOR)法という比較的有名なアルゴリズムを導出する上で重要である.


  1. 一方で,共役勾配法(CG法)などに代表される非定常反復法では,これらの反復行列・反復ベクトルにもう少し複雑な情報をもたせることで収束性を良くする.一方,そのような場合には収束性に関する数学的な考察が難しくなる.

  2. この,見かけ上犠牲というのが微妙なステートメントになってしまったが,要するにJacobi法のようなstraightforwardな並列化はできないよ,という意味でしかない.実際,順序付けに関する工夫というのを行えば並列性が回復できるらしいが,その順序付けに関しては筆者も未だによく理解していないので,とりあえず大阪電通大の伊藤祥司氏の学位論文へのリンクを貼っておく: 大規模線形方程式に対するスーパーコンピュータ向きの前処理法 - 国立国会図書館デジタルコレクション(直でダウンロードはできないが,章毎にGoogle検索でPDFが引っかかり,その中の一つに赤黒順序付けという最も簡単な順序付けに関する解説がある)