Jacobi 1만 번이 못 푸는 것을 V-사이클 10번이 푼다

ω=2/3이면 최고주파 감쇠는 1/3이고, 80번 Jacobi 뒤에 남는 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(weighted Jacobi)의 오차 전파를 분석하면, 오차의 푸리에 모드 $\sin(k \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$ 셀의 성긴 격자로 옮기면 같은 파장이 격자 입장에서는 절반의 셀만 차지한다. 더 짧은 파장이 된다. 더 정확히, fine 격자의 가장 낮은 모드는 coarse 격자에서 중간 모드로 보인다. 중간 모드는 그곳의 Jacobi가 잘 죽인다.

격자를 한 단계 더 성기게 만들면 한 단계 더 짧은 파장이 된다. $\log_2 N$ 단계까지 내려가면 모든 모드는 어딘가에서 고주파가 된다. 모든 층에서 Jacobi를 몇 번씩만 돌리면 모든 오차가 죽는다. 이론상 그렇다. 저주파를 반복법 쪽에 남겨 두고 견주는 기록이 Jacobi부터 BiCGSTAB까지의 수렴이다.

V-사이클 — 내려갔다 올라오는 길#

이 판단을 알고리듬으로 만든 것이 V-사이클이다. 한 V-사이클은 다음과 같다.

사전 완화 (pre-smoothing) — fine 격자에서 Jacobi $\nu_1$ 회.
잔차 계산 — $r_h = f_h - A_h u_h$ .
제한 (restriction) — $r_{2h} = R\, r_h$ . 잔차를 성긴 격자로 옮긴다.
재귀 호출 — coarse 격자에서 보정 $A_{2h} e_{2h} = r_{2h}$ 를 V-사이클로 푼다 (성긴 격자가 충분히 작으면 직접 풀이).
연장 (prolongation) — $e_h = P\, e_{2h}$ . 보정을 fine 격자로 올린다.
보정 — $u_h \leftarrow u_h + e_h$ .
사후 완화 (post-smoothing) — fine 격자에서 Jacobi $\nu_2$ 회.

내려가는 도중엔 잔차만 들고 가고, 올라오는 도중엔 보정만 들고 온다. V자처럼 한 번 깊이 들어갔다가 다시 올라온다. 한 사이클의 전체 비용은 $O(N)$ 이다. fine 격자 한 번 돌리는 비용에 기하급수로 작아지는 항을 더한 것이기 때문이다. 그리고 한 사이클이 잔차를 한 자릿수씩 떨어뜨린다.

Restriction과 Prolongation — 두 격자 사이의 환율#

두 격자 사이를 옮기는 두 연산자는 서로 짝을 이뤄야 한다. 1D vertex-centered 격자에서 가장 단순한 짝은 다음이다.

Full weighting 제한: coarse 점 $i$ 의 잔차는 fine 점 $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}

선형 연장: coarse 격자 점 위에서는 그대로, 사이 점에서는 두 이웃의 평균.

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

이 둘은 transpose 관계 ( $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:
    """한 격자 수준의 상태: 해, 우변, 격자간격."""
    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용 가중 잡코비. 양 끝은 디리클레 0."""
    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-사이클은 fine 격자 기준 4번(2 pre + 2 post) 사용한다. 비용은 비슷한데 잔차 떨어지는 속도는 비교가 안 된다. 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-사이클은 fine에서 pre 2회와 post 2회, 잔차는 보통 $10^{-10}$ 아래다. $10^{-2}$ 와 $10^{-10}$ 사이는 여덟 자리다.

가장 성긴 층은 $n \le 4$ 에서 Jacobi 50회다. 64셀은 64, 32, 16, 8, 4로 네 번 제한된다. full weighting은 $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-사이클이 마법은 아니다. 실무에서 자주 깨지는 두 지점.

거친 계수. 이방성 또는 점프성 계수의 Poisson에서는 Jacobi가 적당한 smoother가 되지 못한다. $-\nabla \cdot (k(x) \nabla u) = f$ 에서 $k$ 가 자릿수로 점프하면, 표준 V-사이클은 잔차가 어느 자릿수에서 멈춘다. 해결책은 line-Jacobi, 정렬 격자, 또는 Algebraic Multigrid(AMG)다. 격자가 아니라 행렬에서 직접 계층을 만든다.

비등각 격자. 비정렬 격자나 AMR(Adaptive Mesh Refinement)에서는 fine에서 coarse로 가는 대응이 자명하지 않다. agglomeration multigrid는 fine 셀들을 묶어 coarse 셀을 정의하고, restriction은 부피 가중 평균, prolongation은 cell-centered 보간을 쓴다. 구현은 복잡하지만 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에서 읽은 숫자다.

Jacobi 1만 번이 못 푸는 것을 V-사이클 10번이 푼다 — 멀티그리드의 진짜 비밀