ヤコビ法と並列計算

ヤコビ法（Jacobi反復法）は、巨大かつ疎な連立一次方程式 $A\mathbf{x}=\mathbf{b}$ を、行列の対角成分を用いた反復更新で解く定常反復法である。差分化・有限要素化で現れるポアソン型・拡散型・弾性型の線形方程式を、単純な構造で扱える点が特徴である。

参考ドキュメント

• SolverIterative-25-29（ヤコビ法・GS法の収束条件と対角優位性, 東京大学） https://nkl.cc.u-tokyo.ac.jp/15n/SolverIterative-25-29.pdf
• R. Barrett et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods（反復法の定番資料, PDF） https://www.netlib.org/templates/templates.pdf
• MPI programming by the example of Jacobi computations（MPIと領域分割の例, PDF） https://www.uio.no/studier/emner/matnat/ifi/IN3200/v22/teaching-material/in3200-jacobi-2022.pdf

1. 取り扱う問題と記号

$n$ 元連立一次方程式

$$A\mathbf{x}=\mathbf{b},A\in {\mathbb{R}}^{n\times n},\mathbf{x},\mathbf{b}\in {\mathbb{R}}^n$$

を考える。$A$ は疎であり、$n$ が大きくても $A$ の全要素を保持しない（スパース形式で保持する）状況を想定する。

行列分割として

$$A=D+L+U$$

を用いる。ここで $D$ は対角成分のみ、$L$ は下三角（対角成分を除く）、$U$ は上三角（対角成分を除く）である。

2. ヤコビ法の反復式

ヤコビ法は、$D$ のみを左辺に残して反復更新する方法である。

$$(D+L+U)\mathbf{x}=\mathbf{b}\Rightarrow D\mathbf{x}=\mathbf{b}-(L+U)\mathbf{x}$$

よって、反復 $k=0,1,2,\dots$ に対して

$${\mathbf{x}}^{(k+1)}=D^{-1}(\mathbf{b}-(L+U){\mathbf{x}}^{(k)})$$

となる。

成分表示は

$$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j\ne i}{a}_{ij}{x}_j^{(k)})(i=1,\dots ,n)$$

である。更新に必要なのは「前の反復 ${\mathbf{x}}^{(k)}$」のみであり、各 $i$ の更新は互いに独立である。

3. 誤差伝播と収束条件

真の解を ${\mathbf{x}}^{\ast }$、誤差を ${\mathbf{e}}^{(k)}={\mathbf{x}}^{(k)}-{\mathbf{x}}^{\ast }$ とする。ヤコビ法の誤差は

$${\mathbf{e}}^{(k+1)}=B{\mathbf{e}}^{(k)},B=-D^{-1}(L+U)$$

に従う。したがって収束条件は

$$\rho (B)<1$$

（$\rho$ はスペクトル半径）である。

よく使われる十分条件として「対角優位」がある。例えば行ごとに

$$|a_{ii}|>\sum _{j\ne i}|a_{ij}|$$

が成り立つと、ヤコビ法は収束しやすい（ただし十分条件であり必要条件ではない）。

4. 緩和ヤコビ法

ヤコビ法は単純である一方、収束が遅い場合が多い。そこで緩和係数 $\omega$ を導入した緩和ヤコビ法（Weighted / Relaxed Jacobi）が用いられる。

$${\mathbf{x}}^{(k+1)}=(1-\omega ){\mathbf{x}}^{(k)}+\omega D^{-1}(\mathbf{b}-(L+U){\mathbf{x}}^{(k)})$$

$\omega =1$ で基本のヤコビ法に一致する。$\omega$ を適切に選ぶことで、特に多重格子法の平滑化（高周波誤差の減衰）に有用な更新となることが多い。

5. 反復停止条件

計算では誤差 ${\mathbf{e}}^{(k)}$ を直接評価できないため、残差

$${\mathbf{r}}^{(k)}=\mathbf{b}-A{\mathbf{x}}^{(k)}$$

で停止判定する。典型例は

$$\frac{\parallel {\mathbf{r}}^{(k)}{\parallel }_2}{\parallel \mathbf{b}{\parallel }_2}<\epsilon \text{または}\parallel {\mathbf{r}}^{(k)}{\parallel }_2<\epsilon$$

である。疎行列では ${\mathbf{r}}^{(k)}$ の計算は「疎行列ベクトル積（SpMV）」として実装されることが多い。

6. ガウス・ザイデル法との違い

ヤコビ法とガウス・ザイデル法の反復式を並べると差が見えやすい。

• ヤコビ法（全成分を前反復から計算）

$$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j<i}{a}_{ij}{x}_j^{(k)}-\sum _{j>i}{a}_{ij}{x}_j^{(k)})$$

• ガウス・ザイデル法（更新済み成分を即座に利用）

$$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum _{j<i}{a}_{ij}{x}_j^{(k+1)}-\sum _{j>i}{a}_{ij}{x}_j^{(k)})$$

ガウス・ザイデル法は収束が速い場合が多い一方、$x_1^{(k+1)}\to x_2^{(k+1)}\to \dots$ の依存が生じるため、素朴には並列化が難しい。ヤコビ法はその逆で、収束は遅めでも並列化しやすい。

7. 並列化

7.1 共有メモリ並列（スレッド並列）

ヤコビ法は、反復 $k$ における更新が

• 入力：${\mathbf{x}}^{(k)}$
• 出力：${\mathbf{x}}^{(k+1)}$ という二配列モデルで完結する。

したがって、$i$ ループをスレッドで分割しやすい。計算の中心は

• 各行の非零要素 $a_{ij}$ と $x_j^{(k)}$ の積和（SpMVまたは行単位の積和）
• 最後に $x_i^{(k+1)}$ を書き込む である。

注意点として、ヤコビ法はしばしば「メモリ帯域律速」になりやすい。$A$（疎行列）と ${\mathbf{x}}^{(k)}$ の参照が不規則になり、キャッシュ効率が性能を左右する。

7.2 分散メモリ並列

空間離散化に由来する $A$ は、格子・要素が空間的に局在した近傍結合を持つことが多い。そこで領域分割（domain decomposition）により、未知数をプロセスごとに分担する。

例として2次元5点差分（ポアソン型）を想起すると、各プロセスは自領域内部の更新はローカルに完結するが、境界付近の更新には隣接領域の値が必要である。これをハロー（ゴーストセル）領域に受け取り、反復ごとに近傍通信する。

概念的には各反復で

1. 隣接プロセスと境界データを交換（ハロー更新）
2. 自領域の全点でヤコビ更新（完全にローカル）
3. 必要なら残差ノルムを全体集約（グローバル還元） という流れになる。

7.3 グローバル還元

残差ノルム $\parallel {\mathbf{r}}^{(k)}{\parallel }_2$ は

$$\parallel \mathbf{r}{\parallel }_2^2=\sum _{i=1}^n{r}_i^2$$

であり、分散メモリでは部分和を足し上げる全体通信（例：Allreduce）が必要になる。反復ごとにこれを行うと、通信遅延が反復全体を支配し得る。

対策としては

• ノルム評価の頻度を下げる（毎反復ではなく数反復ごと）
• 近似的な停止判定（ローカル基準を併用）
• 非同期反復（次節） などが検討対象となる。

7.4 非同期ヤコビ

大規模並列では「同期」のコストが顕在化する。非同期ヤコビは、厳密に同じ ${\mathbf{x}}^{(k)}$ を全プロセスで揃えてから更新するのではなく、

• 受け取れた最新の境界値を使って更新を進める
• 同期点を減らす という発想である。

ただし、非同期化は収束解析が難しく、通信遅延・近似誤差・停止判定の設計と一体で扱う必要がある。

7.5 ブロックヤコビ

ヤコビ法は反復法であると同時に「前処理」の基本形でもある。前処理行列を

$$P=D$$

とみなせば、$P^{-1}A$ の条件を改善して Krylov 法（CG/GMRES等）を加速する入口になる。

分散並列では、対角だけでなく各プロセスのローカルブロックを $P$ とするブロックヤコビ前処理も自然である。これは「各領域内はローカルに解き、領域間は反復でつなぐ」という並列向きの設計である。

8. 位置づけ

ヤコビ法単体は、難しい問題では収束が遅くなりやすい。その一方で、並列性の高さから以下の用途で重要である。

• 多重格子法の平滑化（緩和ヤコビ）
• Krylov反復法の前処理（ヤコビ前処理、ブロックヤコビ）
• 近似解で十分な場面の初期反復（粗い解の獲得）

9. 性質の比較

手法	反復更新の依存	収束傾向	並列化のしやすさ	典型的な役割
ヤコビ	前反復のみ参照（独立）	遅めになりやすい	高い（データ並列向き）	平滑化・前処理・初期反復
緩和ヤコビ	前反復のみ参照（独立）	$\omega$で調整可能	高い	多重格子の平滑化
ガウス・ザイデル	更新済み成分に依存	速い場合が多い	低い（工夫が必要）	小規模・逐次向き、または色分割等で並列化
SOR/SSOR	GSに緩和を追加	速い場合が多い	低〜中	収束加速、ただし設計が難しい