マルチグリッド法の一般的な考え方について

1 はじめに大規模連立一次方程式に対する有効な解法の一つとしてマルチグリッド法が知られている. また,マルチグリッド法は近似的に用いることで, CG 法や Bi-CGSTAB 法などの Krylov 部分空間法の前処理としても用いられる.
マルチグリッド法は, 偏微分方程式の離散化の情報を元に設計される幾何的マルチグリッド法と代数的に行う代数的マルチグリッド法に分けられる. 以下では, 幾何的マルチグリッド法および代数的マルチグリッド法の概略について紹介する.
マルチグリッド法に関する詳細については [4, 3.7 節] や [5], また Krylov 部分空間法やその前処理に関しては [1–4] 等を参照されたい.

2 幾何的マルチグリッド法
境界 Γ において Dirichlet 条件を課した, 単位正方領域 Ω = [0, 1] × [0, 1] 内部の Laplace 方程式∂2u ∂2u+ 2 = 0 ((x, y) ∈ (0, 1) × (0, 1))∂x2∂yu(x, y) = g(x, y) ((x, y) ∈ Γ)を刻み幅 h = 1/N の格子を用いた 5 点中心差分法

4ui,j − ui+1,j − ui−1,j − ui,j+1 − ui,j−1 = 0(1 ≤ i, j ≤ N − 1)

で離散化して得られる方程式

Ah uh = bh(1)を例に, 幾何的マルチグリッド法の概略を示す.

2.1 2 段グリッド法
偏微分方程式を離散化して得られる方程式 (1) に対して Jacobi 法や Gauss-Seidel 法などの定常反復法を適用すると, 誤差の持つ空間的高周波成分が速く減衰すること (平滑化作用) が知られてい1る. このため空間的低周波成分で構成される粗い格子に対応する低次元の方程式

A2h u2h = b2h(2)を用いて, 方程式 (1) の近似解の精度を効率的に改善できることが期待される. この方法を 2 段グリッド法と呼ぶ.
細かい格子上の近似解 uh から粗い格子上の近似解 u2h を作ることを制限 (restriction) と呼び,一般には横長の長方行列 Rh を用い u2h = Rh uh とする. 制限行列 Rh の設定の最も単純な方法は
hu2hi,j = u2i,2j と設定することであるが, 通常は近傍の 9 点を用い,u2hi,j =

1 h
1 ∑ h
1 ∑ h
1 ∑ hu2i,2j +u2i+p,2j +u2i,2j+q +u4
8 p=±18 q=±1
16 q=±1 2i+p,2j+qとなるように設定する.
逆に, 粗い格子上の近似解 u2h から細かい格子上の近似解 uh を作ることを補完 (interpolation)または延長 (prolongation) と呼び, 一般には縦長の長方行列 P h を用い uh = P h u2h とする. 補完
行列 P h の設定の標準的な方法としては,uh2i,2j = u2hi,j ,1
uh2i+1,2j = (u2h+ u2hi+1,j ),2 i,j1
uh2i,2j+1 = (u2h+ u2hi,j+1 ),2 i,j
12h2h+ u2h
uh2i+1,2j+1 = (u2hi+1,j + ui,j+1 + ui+1,j+1 )4 i,jとすることである.

2 段グリッド法のアルゴリズムを Algorithm 1 に記す.
Algorithm 1 2 段グリッド法
1: 方程式 Ah uh = bh に対し定常反復法を数反復適用する.
2: 残差 r h := bh − Ah uh の粗い格子への制限 r 2h = Rh r h を計算する.
3: 粗い格子の方程式 A2h v 2h = r 2h を解く.

e h = uh + P h v 2h と修正する.
4: v 2h の補完 P v 2h を用い, 細かい格子上の近似解を ue h を初期近似解として, 定常反復法を数反復適用する.
5: 方程式 Ah uh = bh に対し u

2.2 マルチグリッド法
マルチグリッド法は上記の 2 段グリッド法を再帰的に行う方法である. マルチグリッド法を連立一次方程式 (1) に対する解法として用いる場合, 粗い方程式 (2) が十分小さくなるまで再帰的に 2段グリッド法を適用し, 十分小さい方程式を LU 分解等の直接法で解く. 再帰の具体的な方法については V サイクルや W サイクルなどいくつかの方法が提案されている.

2

一方, マルチグリッド法を前処理として用いる場合は, 2 段グリッド法を数段階だけ再帰的に行い, この時, 最も粗い方程式に対しては定常反復法を数反復適用することで近似解を得る.

3 代数的マルチグリッド法
上述のように, 幾何的マルチグリッド法は, 偏微分方程式の離散化から得られる方程式 (1) を対象とした手法であり, 粗い方程式 (2) の係数行列 A2h , 制限行列 Rh , 補完行列 P h は対象とする偏微分方程式の離散化に基づき決定される. 一方, 一般の線形方程式に対しても同様の手法を構成することが出来る. この手法を代数的マルチグリッド (AMG) 法と呼ぶ.
一般の n 次の連立一次方程式

Au = bに対する代数的マルチグリッド法は, 幾何的マルチグリッド法と同様に, Algorithm 2 に示す 2 段グリッド法を再帰適用することで得られる.

Algorithm 2 代数的マルチグリッド法における 2 段グリッド法
1: 方程式 Au = b に対し定常反復法を数反復適用する.

b = Rr を計算する.
2: 残差 r := b − Au の粗い格子への制限 r
bv = rb を解く.
3: 粗い格子の方程式 Abb の補完 P v
b を用い, 細かい格子上の近似解を u
e = u + Pvb と修正する.
4: v
e を初期近似解として, 定常反復法を数反復適用する.
5: 方程式 Au = b に対し u以下では, 代数的マルチグリッド法における 2 段グリッド法で現れる粗い格子の方程式

bv = rbAb(3)b および制限行列 R, 補完行列 P の設定法の概略について示す.の係数行列 A
幾何的マルチグリッド法は, 定常反復法を適用すると誤差の空間的高周波成分が速く減衰するという平滑化作用に着目し, 空間的低周波成分で構成される粗い格子を用いて近似解を改善する. 一
方代数的マルチグリッド法では, この性質を利用し, 定常反復法で減衰する成分を空間的高周波成分, 減衰しない成分を空間的低周波成分 (代数的に滑らかな成分) と定義する. つまり, 定常反復法の反復行列を S とすると, Se ≈ e を満たすベクトル e を代数的に滑らかな成分と定義する. ここで, 定常反復法として (減速)Jacobi 法を用い, また係数行列 A の対角要素の値が要素毎に同程度であると仮定すると,Ae ≈ 0(4)を満たす.
また, 行列 A の第 i 行における比較的大きな非対角要素の列番号を, ある閾値 0 < θ ≤ 1 を用い

Si := {j|j ̸= i, |aij | ≥ θ max |aik |}k̸=iと定義する.

この時, 行列 A が対称な M 行列*1 である場合, 式 (4) より, 各 j*1 全ての非対角要素が零以下: aij ≤ 0, i ̸= j, かつ全ての固有値の実部が正となる行列.

3∈ Si に対して,

ei ≈ ej が成り立つ. (これを「Si の方向に滑らか」と表現する.)変数の番号の集合を V := {1, 2, . . . , n}, 粗い格子に属する番号の集合を C ∈ V , 粗い格子に含まれない集合を F := V \ C と置く. この時, 代数的に滑らかな成分が Si 方向に滑らかである点に着目し, n × |C| の補完行列 P を{b)i =(P vvbi∑

(i ∈ C)(i ∈ F )bjj∈Ci πij v(5)のように設定する. ただし, Ci = C ∩ Si であり, πij は適当な重みである.

b を選んで新しい近似解 ue = u+Pvb
近似解 u の誤差 e が代数的に滑らかな時, 適当なベクトル v
e = e + Pv
b を 0 にできるためには e ∈ Ran(P ) が成り立つ必要がある. このため, 任意のの誤差 e
代数的に滑らかな成分が (少なくとも近似的に)Ran(P ) に属するように, 重み πij を設定する.
式 (4) は, 集合 Fi := F ∩ Si , Wi := Ni \ Si を用い,

∑aii ei ≈ −∑aij ej −j∈Ci∑aik ek −k∈Fiail ell∈Wiと書くことができる. ここで, ail , l ∈ Wi は比較的小さい値であることに注意すると,(aii +∑)∑ei ≈ −ailj∈Cil∈Wi∑

aij ej −aik ekk∈Fiと近似できる. また, 各 k ∈ Fi に対し,ek ≈∑akj ej /j∈Ci∑akj(6)j∈Cjと近似すると,ej ≈∑πij ej ,πij = −j∈Ciaii +[1∑l∈Wi aikaij +∑k∈Fi(a a∑ ik kjm∈Ci akm)]と表現できる. つまり, この重み πij を用いて式 (5) により補完行列 P を定義することで,e ∈ Ran(P ) が近似的に成り立つ.
集合 C の設定に際し, 式 (6) において∑j∈Cj akj ≈ 0 とならないように,

条件 1 i, k ∈ F に対し, k ∈ Si ならば Ci ∩ Ck ̸= ∅との条件を課す. また, 計算量の観点から |C| が小さくなるよう,
条件 2 C は「任意の i, j ∈ C に対して j ̸∈ Si 」を満たす集合の中で極大との条件を課す. しかし, 一般には条件 1 と条件 2 を同時に満たす集合 C が存在するとは限らないため, 条件 1 を満たす集合の中で条件 2 をできるだけ満たすよう集合 C を設定する.

bはまた, 一般的に, 粗い方程式 (3) の係数行列 Ab = RAP,A4表1集合定義VCFNiSiCiFiWi{1, 2, . . . , n}C∈VV \C{j|j ̸= i, aij ̸= 0}{j|j ̸= i, |aij | ≥ θ maxk̸=i |aik |}C ∩ SiF ∩ Si (= Si \ Ci )N i \ Si

各種集合の概略

意味
全ての変数番号の集合
粗い格子に含まれる変数番号の集合
粗い格子に含まれない変数番号の集合
第 i 行中の非零要素の列番号 (̸= i) の集合
第 i 行中の比較的大きな非対角要素の列番号の集合
第 i 行中で比較的大きく, かつ粗い格子に含まれる列番号の集合
第 i 行中で比較的大きく, かつ粗い格子に含まれない列番号の集合
第 i 行中で比較的小さい非零要素の列番号の集合,

制限行列 R は

R = PTe = u + Pv
b の残差が Ran(P ) と直交するようにするたと設定する. この設定法は, 新しい近似解 ue の誤差 ee = e + P v
b が最小となることとめであり, これは行列 A が正定値対称の時には, 近似解 u等価である.
各種集合の概略を表 1 に示す. また, 具体的な集合 C, F の設定法の一例を Algorithm 3 に示す.

Algorithm 3 代数的マルチグリッド法の集合 C, F, Ci , Fi , Wi の設定法の一例
1: 集合 Si , SiT , Ni の設定

Si := {j||aij | ≥ θ maxk̸=i |aik |}SiT := {j|i ∈ Sj }Ni := {j|j ̸= i, ai,j ̸= 0}
2: 集合 C, F の設定

C = ∅, F = ∅, V = {1, . . . , n}, λi = |SiT |Pick an i ∈ V with maximal λi , and set C = C ∪ {i}, V = V − {i}For all j ∈ SiT ∩ VSet F = F ∪ {j}, V = V − {j}For all l ∈ SjT ∩ V, λl = λl + 1
End forFor all j ∈ Si ∩ V, λj = λj − 1If V = ∅ stop, otherwise, go to (⋆)
3: 集合 Ci , Fi , Wi の設定

Ci := C ∩ Si
Fi := F ∩ Si
Wi := Ni \ Si5(⋆)

参考文献[1] Richard Barrett, Michael W. Berry, Tony F. Chan, James Demmel, June Donato, JackDongarra, Victor Eijkhout, Roldan Pozo, Charles Romine and Henk A. van der Vorst,Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods,
SIAM, Philadelphia, PA, 1993.
[2] 藤野清次, 張紹良, 反復法の数理, 朝倉書店, 東京, 1996.
[3] Yousef Saad, Iterative methods for sparse linear systems. 2nd edition, SIAM, Philadelphia,PA, 2003.
[4] 杉原正顯, 室田一雄, 線形計算の数理, 岩波書店, 東京, 2009.
[5] Ulrich Trottenberg and Cornelis W. Oosterlee and Anton Schüller, Multigrid, AcademicPress, San Diego, CA, 2001.

6