Jacobi 1万回でも消えない誤差を、V-サイクル10回で消す ― マルチグリッドの本当の秘密

ω=2/3 なら最高周波の減衰は 1/3 で、Jacobi 80 回のあと残る 1e-2 は粗い格子へ渡す低周波です。

Jacobi 1万回が終わらない理由#

エンジニア B は 1024×1024 格子で圧力 Poisson を点 Jacobi で解いていました。200 秒たっても残差は同じ桁で揺れます。CFL は小さく、格子は直交です。コードは正しい。それでも Jacobi が決して消さない誤差があります。

同じ誤差でも格子が変われば周波数が変わります。V-サイクルはその一行から出ます。

Pre/Post を 1 と 5 の間で動かすと、同じ 64 セル・同じゼロ右辺で V-サイクルの残差だけが急になり、緩和を足した利得はすぐ折れます。

Pre-smooth iterations: 2Post-smooth iterations: 2

Both solvers attack the same homogeneous Poisson problem with random initial guess. Jacobi (orange) does 4 sweeps per step; V-cycle (cyan) does 2 pre + 2 post sweeps per level on a 6-level hierarchy.

上のオレンジは Jacobi です。振動は収まり、大きな山は残ります。シアンの対数残差はほぼ直線で降り、オレンジは低周波の平坦線に捕まります。

誤差の周波数 ― 早く死ぬものと永遠に生きるもの

線形系 $A u = f$ を解いているとします。現在の近似 $u^{(k)}$ の誤差は

e^{(k)} = u_{\text{exact}} - u^{(k)}

一本のベクトルです。残差は $r^{(k)} = f - A u^{(k)} = A e^{(k)}$ です。

重み付き Jacobi の誤差伝播では、Fourier モード $\sin(j \pi x)$ は一回の掃引で

e^{(k+1)}_j = \left(1 - \omega + \omega \cos(j \pi h)\right) e^{(k)}_j

だけ減衰します。 $\omega$ は緩和係数、 $j$ はモード番号、 $h$ は格子間隔です。減衰係数の絶対値は次の表です。

モード	$j h$	減衰係数 ( $\omega = 2/3$ )
最も振動的	1.0	$\approx 1/3$
中間	0.5	$\approx 0.67$
最も緩やか	$1/N$	$\approx 1.0 - O(h^2)$

短波長は一回で 1/3 に縮みます。長波長はほとんど動きません。Jacobi が消せないのは、格子全体に広がる滑らかな低周波です。

粗い格子で見るとその低周波は高周波になる

ここでマルチグリッドの判断が分かれます。 $N$ セルで表現できる最長波長は $\lambda \sim N h$ です。同じ波長を $N/2$ セルへ移すと、新しい格子から見てセル数は半分です。短い波長に見えます。正確には、細かい格子の最低モードは粗い格子の中間モードです。中間モードはそこでの Jacobi がよく消します。

もう一段粗くすれば、波長はさらに短くなります。 $\log_2 N$ 段まで降りれば、すべてのモードはどこかで高周波になります。各段で Jacobi を数回ずつ回せば、理屈の上ではすべての誤差が消えます。低周波を反復法側に残して見比べる記録が Jacobi から BiCGSTAB までの収束です。

V-サイクル ― 降りて再び上る道#

この判断を手順にしたのが V-サイクルです。一回の中身は次です。

前平滑 (pre-smoothing) — 細かい格子で Jacobi を $\nu_1$ 回。
残差 — $r_h = f_h - A_h u_h$ 。
制限 (restriction) — $r_{2h} = R\, r_h$ 。残差を粗い格子へ移します。
再帰 — 粗い格子で $A_{2h} e_{2h} = r_{2h}$ を V-サイクルで解きます。十分小さければ直接解きます。
延長 (prolongation) — $e_h = P\, e_{2h}$ 。補正を細かい格子へ上げます。
補正 — $u_h \leftarrow u_h + e_h$ 。
後平滑 (post-smoothing) — 細かい格子で Jacobi を $\nu_2$ 回。

降りるときは残差だけ、上るときは補正だけを持ちます。道は V の字です。一サイクルの費用は $O(N)$ です。細かい格子を一回掃く費用に、幾何級数で小さくなる項が足されるからです。一サイクルで残差はおよそ一桁落ちます。

制限と延長 ― 二つの格子の為替レート

二つの格子を結ぶ演算子は対でなければなりません。1D の vertex-centered 格子で最も単純な対は次です。

フルウェイティング制限: 粗い点 $i$ の残差は、細かい点 $2i-1, 2i, 2i+1$ の加重平均です。

r_{2h, i} = \frac{1}{4} r_{h, 2i-1} + \frac{1}{2} r_{h, 2i} + \frac{1}{4} r_{h, 2i+1}

線形延長: 粗い点の真上はそのまま、間の点は両隣の平均です。

e_{h, 2i} = e_{2h, i}, \qquad e_{h, 2i+1} = \frac{1}{2}\left(e_{2h, i} + e_{2h, i+1}\right)

二つは転置の関係です( $P = c R^T$ 、 $c$ は格子比の定数)。V-サイクルの収束の保証にこの対が入っています。非対称な対は発散することがあります。

Python で追う 2 格子 V-サイクル#

1D Poisson $-u''(x) = f(x)$ 、 $u(0)=u(1)=0$ を V-サイクルで解きます。64 セルから始めます。

import numpy as np
 
class GridLevel:
    """1 つの格子レベルの状態：解、右辺、格子間隔。"""
    def __init__(self, n):
        self.n = n              # セル数
        self.h = 1.0 / n
        self.u = np.zeros(n + 1)
        self.f = np.zeros(n + 1)
 
def jacobi_relax(u, f, h, iters):
    """1D Poisson 用の重み付き Jacobi。両端は Dirichlet ゼロ。"""
    n = len(u) - 1
    nxt = np.zeros_like(u)
    for _ in range(iters):
        nxt[1:n] = 0.5 * (u[0:n-1] + u[2:n+1] + h * h * f[1:n])
        u[1:n] = nxt[1:n]
    return u
 
def residual_norm(u, f, h):
    n = len(u) - 1
    r = np.zeros_like(u)
    inv = 1.0 / (h * h)
    r[1:n] = f[1:n] - (-u[0:n-1] + 2 * u[1:n] - u[2:n+1]) * inv
    return float(np.sqrt(np.mean(r[1:n] ** 2)))
 
def restrict_full_weighting(rh):
    n = len(rh) - 1
    n2 = n // 2
    r2 = np.zeros(n2 + 1)
    r2[1:n2] = 0.25 * rh[1:n-1:2] + 0.5 * rh[2:n:2] + 0.25 * rh[3:n+1:2]
    return r2
 
def prolong_linear(e2h):
    n2 = len(e2h) - 1
    n = 2 * n2
    eh = np.zeros(n + 1)
    eh[2:n:2] = e2h[1:n2]
    eh[1:n+1:2] = 0.5 * (e2h[0:n2] + e2h[1:n2+1])
    return eh
 
def v_cycle(u, f, h, n_pre=2, n_post=2):
    n = len(u) - 1
    if n <= 4:
        return jacobi_relax(u, f, h, 50)
    jacobi_relax(u, f, h, n_pre)
    inv = 1.0 / (h * h)
    r = np.zeros_like(u)
    r[1:n] = f[1:n] - (-u[0:n-1] + 2 * u[1:n] - u[2:n+1]) * inv
    r2 = restrict_full_weighting(r)
    e2 = np.zeros_like(r2)
    v_cycle(e2, r2, 2 * h, n_pre, n_post)
    e = prolong_linear(e2)
    u += e
    jacobi_relax(u, f, h, n_post)
    return u
 
# --- 実行 ---
n = 64
h = 1.0 / n
x = np.linspace(0, 1, n + 1)
f = (np.pi ** 2) * np.sin(np.pi * x)   # 厳密解 u = sin(πx)
 
u_jac = np.zeros(n + 1)
u_mg = np.zeros(n + 1)
 
for k in range(20):
    jacobi_relax(u_jac, f, h, 4)
    v_cycle(u_mg, f, h, n_pre=2, n_post=2)
    print(f"step {k:3d} | jacobi r={residual_norm(u_jac, f, h):.3e} "
          f"v-cycle r={residual_norm(u_mg, f, h):.3e}")

step 一回で Jacobi は 4 回、V-サイクルは細かい格子基準で 4 回(pre 2 + post 2)です。費用は近く、残差の落ち方は比べ物になりません。20 ステップ後、V-サイクルの残差は通常 $10^{-10}$ 未満、Jacobi は $10^{-2}$ 付近で止まります。

表の減衰は $\omega = 2/3$ で閉じます。最も速いモード $jh = 1$ は $\cos(\pi) = -1$ なので、係数は $1 - 2\omega = -1/3$ です。絶対値は表の $1/3$ と一致します。そのモードだけを十回緩和すると $(1/3)^{10} \approx 1.7 \times 10^{-5}$ まで死にます。中間モードは表の $\approx 0.67$ です。0.67 を 20 回掛けると $0.67^{20} \approx 3.3 \times 10^{-4}$ です。最も遅いモードは $1 - O(h^2)$ です。 $N = 64$ 、 $h = 1/64$ なら $h^2 = 1/4096 \approx 2.4 \times 10^{-4}$ です。一回で削れる量がその桁です。

出力と噛み合います。step 一回の Jacobi は 4 回です。20 step なら 80 回です。残差は $10^{-2}$ 付近に残ると書いてあります。高周波は $1/3$ の減衰で、その前に死にます。残った $10^{-2}$ は表の低周波です。同じ 20 step の V-サイクルは細かい格子で pre 2 回と post 2 回、残差は通常 $10^{-10}$ 未満です。 $10^{-2}$ と $10^{-10}$ の間は八桁です。

最も粗い層は $n \le 4$ で Jacobi を 50 回です。64 セルは 64, 32, 16, 8, 4 と四回制限されます。フルウェイティングは $1/4$ , $1/2$ , $1/4$ です。粗い層の費用の和は $N/2 + N/4 + \cdots$ で、 $N$ を超えません。V-サイクル 10 回はその和を十回まわす費用です。Jacobi の 1 万回とは、回数の桁が先に違います。

最も遅いモードに同じ式を入れます。 $N = 64$ 、 $jh = 1/N$ 、 $\omega = 2/3$ なら減衰は約 $0.9992$ です。一覧の 80 回を掛けても約 $0.94$ です。一桁削るのに約 2900 回かかります。 $N = 1024$ に同じ式を入れると減衰は約 $0.999997$ で、一桁に約 $7.3 \times 10^{5}$ 回です。1 万回掛けてもそのモードは約 $0.97$ に残ります。1024×1024 で 200 秒、残差が同じ桁にいた場所と同じです。この式は 1D の重み付き Jacobi です。エンジニアの格子は 2D でした。一サイクルで一桁落ちる V-サイクルなら、その遅いモードを掃引が動かさない仕事を、およそ十回で進めます。

$0.9992$ は測った残差ではありません。表の式を最も遅いモードに入れた値です。80 回掛けても動くのはおよそ 6% です。20 ステップ後に $10^{-2}$ 付近と書かれた Jacobi の残差は、減衰 $1/3$ の速いモードが先に死んだあと、1 に近い係数だけが残ったときの姿です。 $N = 1024$ では同じ係数が 1 万回のあと約 $0.97$ を残すので、桁は動きません。V-サイクルが 1 回でおよそ一桁落とすという見積もりと、この 1 に近い係数は別の時間スケールです。

それでも常に速いとは限らない ― 二つの罠

V-サイクルは魔法ではありません。実務でよく折れる場所が二つあります。

粗い係数。 異方性、あるいは係数が跳ぶ Poisson では Jacobi は適切な smoother になりません。 $-\nabla \cdot (k(x) \nabla u) = f$ で $k$ が桁で跳ぶと、標準 V-サイクルは残差をどこかの桁に止めます。対処は line-Jacobi、整列格子、あるいは代数的マルチグリッド(AMG)です。格子ではなく行列から階層を作ります。

非整合な格子。 非構造格子や AMR では fine から coarse への対応が自明ではありません。agglomeration multigrid は fine セルを束ねて coarse セルを定義し、制限は体積加重平均、延長はセル中心の補間を使います。実装は重いですが、OpenFOAM の GAMG がこの方式です。

収束率が格子独立になるには、smoother と演算子と格子比の三つが揃う必要があります。一つでもずれると、反復回数は格子サイズとともに少しずつ増えます。そのとき W-サイクルか F-サイクルを見ます。

標準 V-サイクルを下ろす場所#

$k$ が桁で跳ぶなら、標準 Jacobi は smoother ではありません。残差はどこかの桁に捕まります。line-Jacobi か AMG へ乗り換えます。引き伸ばした境界層でのその乗り換えを書いたのがラインレット前処理です。

非構造や AMR では、添字を半分にして coarse を作れません。OpenFOAM の GAMG はセルを束ねます。反復回数が $N$ に沿って増えるなら W-サイクルか F-サイクルです。64 セル、 $n \le 4$ での Jacobi 50 回、20 ステップの $10^{-10}$ 対 $10^{-2}$ は、整列した 1D Poisson で読んだ数字です。