셀 번호만 섞었더니 LU-SGS가 2.8배 느려졌다 — RCM 재정렬과 시작 셀이 정하는 스윕 방향

LU-SGS에서 셀 번호는 정보가 한 스윕에 얼마나 멀리 가는지를 정한다. RCM은 대역폭을 24로 줄이지만, 반복 횟수 9회와 21회를 가르는 것은 시작 셀과 흐름 방향의 관계다.

셀 번호를 섞으면 무엇이 바뀔까#

같은 격자다. 셀 번호만 무작위로 섞었다. 행렬은 행과 열의 순서만 바뀐다. 수렴한 해는 그대로다. 그런데 LU-SGS(Lower-Upper Symmetric Gauss–Seidel, 하삼각·상삼각 두 번의 스윕으로 음해 시스템을 근사하는 기법) 반복은 9회에서 25회가 됐다. 차이는 "정보가 한 스윕에 얼마나 멀리 가는가"에 있다. RCM(Reverse Cuthill–McKee, 너비 우선 탐색으로 행렬 대역폭을 줄이는 재번호 기법)이 고치는 것과 못 고치는 것을, 576셀 출력으로 가른다.

2017년 노트의 질문: 어느 이웃이 Lower인가#

2017년 말 비정렬 3D 코드에 LU-SGS를 넣으면서 남긴 메모가 있다. 요지는 이렇다. LU-SGS는 Lower 부분과 Upper 부분을 나눠 계산한다. 그런데 비정렬 격자에서는 주변 셀 중 어느 것이 Lower이고 어느 것이 Upper인지 알기 어렵다. 그래서 RCM 오픈소스를 찾아 초기화 단계에서 셀을 재정렬했다.

지금 보면 이 질문의 앞부분은 반만 맞다. Lower와 Upper는 격자가 정하지 않는다. 셀에 붙인 번호가 정한다. 번호가 나보다 작은 이웃이 Lower, 큰 이웃이 Upper다. 어떤 번호든 붙어 있기만 하면 LU-SGS는 돈다. 진짜 질문은 "어느 번호가 좋은가"다.

메모에는 "이 RCM 코드는 사면체 격자에만 적용된다는데, 다른 격자에도 될 것 같다"는 줄도 있다. 이건 맞다. RCM은 셀 모양을 보지 않는다. 면을 공유하는 이웃 목록, 즉 그래프만 본다.

random labels와 RCM을 바꾸면, 20×20 셀의 1차 풍상 대류–확산에서 색이 퍼지는 길이 달라진다.

RCM start:

flow angle = 30° diffusion D = 0.010 (cell Pe = 100) cells/frame = 16

sweeping| iter = 0 | residual = 1.0e+0natural: 9 it, bw 20random: 21 it, bw 392RCM SW: 9 it, bw 20RCM SE: 18 it, bw 20RCM NW: 18 it, bw 20RCM NE: 9 it, bw 20

무작위 번호는 전진 스윕 한 번에도 색이 군데군데만 남는다. RCM SW는 흐름 화살표를 따라 한 스윕에 거의 다 찬다. RCM start를 SE로 바꾸면 bw는 그대로인데 반복 횟수만 길어진다.

LU-SGS 두 번의 스윕을 식으로 쓰면#

정상 상태 선형 시스템 $A\phi = b$ 를 의사시간(pseudo-time) 전진으로 푼다고 하자. 행렬을 대각과 비대각으로 나눈다.

A = D - L - U, \qquad (D + \tfrac{1}{\Delta\tau}I - L)\,(D + \tfrac{1}{\Delta\tau}I)^{-1}\,(D + \tfrac{1}{\Delta\tau}I - U)\,\Delta\phi = r

$D$ 는 대각 $a_P$ , $L$ 은 번호가 작은 이웃의 계수 $a_{nb}$ , $U$ 는 번호가 큰 이웃의 계수, $r = b - A\phi$ 는 잔차, $\Delta\tau$ 는 의사시간 간격이다.

이 근사 분해는 두 번의 스윕으로 풀린다. 전진 스윕은 번호 순서대로 지나가며 작은 번호 이웃의 새 값만 쓴다. 후진 스윕은 역순으로 지나가며 큰 번호 이웃의 값만 더한다.

\Delta\phi^{*}_c = \frac{r_c + \sum_{k<c} a_{ck}\,\Delta\phi^{*}_k}{d}, \qquad \Delta\phi_c = \Delta\phi^{*}_c + \frac{\sum_{k>c} a_{ck}\,\Delta\phi_k}{d}

$k<c$ 는 "순번이 $c$ 보다 앞선 이웃", $d = a_P + 1/\Delta\tau$ 다.

여기서 갈린다. 셀 $c$ 의 상류 이웃이 $c$ 보다 앞 순번이면, 그 이웃이 이번 스윕에 받은 수정값이 바로 $c$ 로 넘어온다. 뒤 순번이면 다음 반복까지 기다려야 한다. 순수 풍상 대류라면, 흐름 방향으로 번호가 매겨진 경우 $L$ 이 행렬 전체를 담는다. 전진 스윕 한 번이 정확한 해가 된다.

Cuthill–McKee는 너비 우선 탐색이다#

Cuthill과 McKee가 1969년에 낸 알고리즘은 단순하다. 시작 셀 하나를 고른다. 너비 우선 탐색(BFS, 가까운 이웃부터 층층이 방문)으로 번호를 붙인다. 같은 층에서는 이웃 수가 적은 셀부터 붙인다. 이렇게 하면 이웃끼리 번호 차이가 한 층 너비를 넘지 않는다.

이 번호 차이의 최댓값이 대역폭이다.

\beta = \max_{(i,j)\,\in\,\text{faces}} \lvert p(i) - p(j) \rvert

$p(i)$ 는 셀 $i$ 의 새 번호, 최댓값은 면을 공유하는 모든 셀 쌍에 대해 잡는다.

1971년 George는 이 순서를 뒤집기만 해도 직접 분해의 채움(fill-in)이 줄어든다는 것을 보였다. 그게 RCM의 R이다. 시작 셀은 보통 그래프의 "가장자리"에서 고른다. BFS를 반복해 가장 먼 셀로 옮겨가는 의사 주변 셀(pseudo-peripheral node) 탐색이 표준이다.

코드: 576셀, 번호 네 가지#

24×24 격자에 두 문제를 건다. 하나는 30° 방향 흐름에 확산이 작은 대류 지배 문제다(셀 Péclet 수 100). 다른 하나는 순수 확산이다. 격자 생성기가 번호를 아무렇게나 붙인 상황을 흉내 내려고 라벨을 섞는다. 그다음 RCM으로 다시 매긴다. 대류 문제에서는 RCM 시작 셀을 네 모서리로 바꿔 본다.

import math, random
from collections import deque
 
def build_grid(n, ux, uy, diff):
    # 셀 (i,j) -> 이웃 리스트와 계수 a_nb, 대각 a_P (1차 풍상 + 확산)
    nbr, coef = [], []
    for j in range(n):
        for i in range(n):
            lst, cf = [], []
            for di, dj, nx_, ny_ in ((1,0,1,0),(-1,0,-1,0),(0,1,0,1),(0,-1,0,-1)):
                ii, jj = i+di, j+dj
                if 0 <= ii < n and 0 <= jj < n:
                    lst.append(jj*n+ii)
                    cf.append(diff + max(-(ux*nx_+uy*ny_), 0.0))
            nbr.append(lst); coef.append(cf)
    ap = 4*diff + abs(ux) + abs(uy)
    return nbr, coef, ap
 
def shuffle_labels(nbr, coef, rng):
    m = len(nbr)
    perm = list(range(m)); rng.shuffle(perm)   # old -> new
    inv = [0]*m
    for o, nw in enumerate(perm): inv[nw] = o
    nb2 = [[int(perm[k]) for k in nbr[inv[c]]] for c in range(m)]
    cf2 = [coef[inv[c]] for c in range(m)]
    return nb2, cf2, perm
 
def bfs_levels(nbr, s):
    dist = [-1]*len(nbr); dist[s] = 0; q = deque([s]); last = s
    while q:
        c = q.popleft(); last = c
        for k in nbr[c]:
            if dist[k] < 0:
                dist[k] = dist[c]+1; q.append(k)
    return dist, last
 
def peripheral_start(nbr):
    s, ecc = 0, -1
    while True:
        dist, far = bfs_levels(nbr, s)
        if dist[far] <= ecc: return s
        ecc, s = dist[far], far
 
def cuthill_mckee(nbr, start):
    seen = [False]*len(nbr); seen[start] = True
    order, q = [], deque([start])
    while q:
        c = q.popleft(); order.append(c)
        for k in sorted(nbr[c], key=lambda k: len(nbr[k])):
            if not seen[k]:
                seen[k] = True; q.append(k)
    return order
 
def bandwidth(nbr, order):
    pos = {c: p for p, c in enumerate(order)}
    return max(abs(pos[c]-pos[k]) for c in range(len(nbr)) for k in nbr[c])
 
def lusgs_solve(nbr, coef, ap, order, dtau, tol=1e-8, itmax=4000):
    m = len(nbr); x = [0.0]*m; b = [1.0]*m
    d = ap + 1.0/dtau
    rank = [0]*m
    for p, c in enumerate(order): rank[c] = p
    bnorm = m**0.5
    for it in range(1, itmax+1):
        r = [b[c] - ap*x[c] + sum(a*x[k] for a, k in zip(coef[c], nbr[c])) for c in range(m)]
        if sum(v*v for v in r)**0.5/bnorm < tol: return it-1
        dx = [0.0]*m
        for c in order:                      # forward: 낮은 rank 이웃만 (L)
            s = r[c] + sum(a*dx[k] for a, k in zip(coef[c], nbr[c]) if rank[k] < rank[c])
            dx[c] = s/d
        for c in reversed(order):            # backward: 높은 rank 이웃만 (U)
            s = sum(a*dx[k] for a, k in zip(coef[c], nbr[c]) if rank[k] > rank[c])
            dx[c] += s/d
        for c in range(m): x[c] += dx[c]
    return itmax
 
rng = random.Random(7)
n = 24
th = math.radians(30)
for name, ux, uy, df, dtau in (("convection", math.cos(th), math.sin(th), 0.01, 1e3),
                               ("diffusion", 0.0, 0.0, 1.0, 0.5)):
    nb0, cf0, ap = build_grid(n, ux, uy, df)
    nb, cf, perm = shuffle_labels(nb0, cf0, rng)
    m = len(nb)
    s = peripheral_start(nb)
    cm = cuthill_mckee(nb, s)
    rcm = cm[::-1]
    # 자연 순서는 셔플 전 격자를 그대로 쓴다
    orders = {
        "natural (i,j)": (nb0, cf0, list(range(m))),
        "random labels": (nb, cf, list(range(m))),
        "Cuthill-McKee": (nb, cf, cm),
        "reverse CM": (nb, cf, rcm),
    }
    print(f"[{name}] cells={m}")
    for k, (g, c, o) in orders.items():
        it = lusgs_solve(g, c, ap, o, dtau)
        print(f"  {k:15s} bandwidth={bandwidth(g, o):5d}  LU-SGS iters={it:5d}")
    if name == "convection":
        for corner, old in (("SW", 0), ("SE", n-1), ("NW", n*n-n), ("NE", n*n-1)):
            o = cuthill_mckee(nb, perm[old])[::-1]
            print(f"  RCM start {corner}    bandwidth={bandwidth(nb, o):5d}  LU-SGS iters={lusgs_solve(nb, cf, ap, o, dtau):5d}")

[convection] cells=576
  natural (i,j)   bandwidth=   24  LU-SGS iters=    9
  random labels   bandwidth=  553  LU-SGS iters=   25
  Cuthill-McKee   bandwidth=   24  LU-SGS iters=    9
  reverse CM      bandwidth=   24  LU-SGS iters=    9
  RCM start SW    bandwidth=   24  LU-SGS iters=    9
  RCM start SE    bandwidth=   24  LU-SGS iters=   21
  RCM start NW    bandwidth=   24  LU-SGS iters=   21
  RCM start NE    bandwidth=   24  LU-SGS iters=    9
[diffusion] cells=576
  natural (i,j)   bandwidth=   24  LU-SGS iters= 1549
  random labels   bandwidth=  561  LU-SGS iters= 1747
  Cuthill-McKee   bandwidth=   24  LU-SGS iters= 1549
  reverse CM      bandwidth=   24  LU-SGS iters= 1549

같은 576셀, 같은 행렬, 같은 수렴 기준( $10^{-8}$ )이다. 바뀐 것은 셀을 지나가는 순서뿐이다.

25회와 9회를 가른 것은 상류 이웃의 순번이다#

흐름이 30°면 각 셀의 상류 이웃은 서쪽과 남쪽이다. 자연 순서( $i$ 먼저, $j$ 다음)에서는 두 이웃 모두 앞 순번이다. 그래서 전진 스윕 한 번이 남서쪽에서 북동쪽까지 정보를 끝까지 나른다. 남은 9회는 작은 확산 항과 의사시간 항을 걷어내는 몫이다.

무작위 번호에서는 상류 이웃이 앞 순번일 확률이 반반이다. 정보가 한 스윕에 몇 셀 가다가 끊긴다. 24셀을 가로지르려면 여러 반복이 필요하다. 그 결과가 25회, 약 2.8배다.

대역폭도 553에서 24로 줄었다. 이 숫자는 반복 횟수와는 다른 이야기다. 실제 코드에서 이웃 셀 데이터가 메모리에서 얼마나 멀리 흩어져 있는지를 뜻한다. 대역폭이 작으면 캐시 적중이 좋아져 반복 한 번이 빨라진다. RCM은 반복 횟수와 반복 비용을 동시에 줄인 셈이다. 다만 Python 실험은 캐시 효과를 재지 못하므로 여기서는 반복 횟수만 비교한다.

같은 대역폭 24, 반복은 9회와 21회#

네 모서리 실험이 더 흥미롭다. 어느 모서리에서 시작해도 대역폭은 24다. 그런데 SW와 NE는 9회, SE와 NW는 21회다.

RCM의 BFS 층은 시작 모서리를 중심으로 한 대각선 띠다. SW에서 시작하면 띠가 흐름에 수직으로 놓이고, 순번이 흐름을 따라 증가한다. NE에서 시작해도 순번이 흐름을 거슬러 증가할 뿐이다. 이때는 후진 스윕이 흐름 방향이 된다. SGS는 전진과 후진을 모두 하므로 둘은 대칭이다. 같은 이유로 Cuthill–McKee와 reverse CM의 반복 횟수도 같다.

SE에서 시작하면 순번이 북서쪽으로 증가한다. 흐름의 $x$ 성분과는 반대, $y$ 성분과는 같은 방향이다. 전진 스윕도, 후진 스윕도 상류 이웃 하나씩만 제때 받는다. 그래서 21회다. 대역폭은 행렬의 모양만 본다. 흐름 방향은 모른다.

25를 9로 나누면 2.78이다. 제목의 2.8배는 그 반올림이다. 같은 576셀, 같은 $10^{-8}$ 이다. 대역폭은 553에서 24로 내려간다. 553은 셀 수 576 옆에 붙는다. 이웃 번호가 행렬 거의 전체를 가로지른다는 뜻이다. 24는 자연 순서가 이미 닿는 $\beta$ 의 바닥이다. 24×24에서 세로 이웃은 $j$ 가 하나 늘 때 번호가 $n=24$ 만큼 떨어진다. RCM은 그 바닥까지 줄였고, 그 아래는 없다.

모서리 네 개의 $\beta$ 는 전부 24다. 반복만 SW 9, NE 9, SE 21, NW 21이다. 21/9는 2.33이라 무작위 25보다는 낫고, 흐름을 따라간 9보다는 두 배가 넘는다. SE와 NW가 같은 21인 것은 스윕이 두 번, 서로 역방향이라서다. 한쪽 대각선에서 9가 나오면 반대 대각선도 9다. CM과 reverse CM이 둘 다 9인 것과 같은 대칭이다.

확산 문제에서는 13%만 차이 난다#

순수 확산에서는 숫자가 다르다. 자연 순서와 RCM은 1549회, 무작위 번호는 1747회다. 차이는 13%다.

확산에는 상류가 없다. 정보는 모든 방향으로 같은 속도로 퍼진다. Gauss–Seidel이 느린 이유는 번호가 아니라 긴 파장이다. 그 파장은 순서로 안 줄고, Krylov 쪽에 넘길 일이다. 번호가 반복 횟수를 크게 가르는 쪽은 대류다. 압축성 Euler처럼 쌍곡형이 센 문제가 그 경우다.

1549와 1747의 차이는 198회다. 198/1549는 0.128, 곧 13%다. 대역폭이 561까지 벌어져도 반복은 그 정도다. 대류에서 553이 25회를 만든 것과 대비된다. $\beta$ 를 줄이는 일과 $k \lt c$ 에 상류를 넣는 일은 같은 재정렬 안에서도 따로 논다. 확산에는 상류가 없어서 두 번째가 거의 비어 있다.

그런데 왜 'Reverse'인가#

실험에서 CM과 RCM의 반복 횟수는 같았다. 뒤집는 이유가 없어 보인다. 효과가 나는 곳은 다른 곳이다.

George가 보인 것은 직접 분해의 채움이다. 대역 행렬을 Cholesky나 LU로 분해하면 각 행의 첫 비영 원소부터 대각까지가 채워진다. 이 "프로파일(profile)" 합이 RCM에서 거의 항상 CM 이하로 나온다. 그래서 ILU 전처리나 직접 해법을 쓸 계획이면 RCM이 기본값이다. LU-SGS만 쓴다면 CM이든 RCM이든 차이가 없다. 대칭인 두 스윕이 순서의 방향을 지운다.

1D에서는 왼쪽에서 오른쪽으로 번호를 매기면 끝이다. 번호 문제는 2D·3D 비정렬에서 생긴다. 격자 생성기는 셀을 만들어진 순서대로 붙인다. 점 삽입이나 프론트 전진을 따른다. 해석 코드에서 보면 무작위 번호에 가깝다.

대역폭 24인데도 21회에서 멈추면 시작 셀을 버린다#

재정렬은 초기화에서 한 번이다. 계산 중에 번호를 다시 매길 이유는 없다. 면 연결로 인접 그래프를 만들고, 셀 배열과 면의 좌우 인덱스와 경계 목록을 새 번호로 옮긴다. 사면체든 프리즘이든 그래프만 보면 된다. 비용은 BFS 한 번, 셀 수에 비례한다.

대류가 세면 시작 셀이 대역폭보다 먼저다. 이번 출력에서 네 모서리의 $\beta$ 는 모두 24인데 SW·NE는 9회, SE·NW는 21회다. 의사 주변 셀은 그래프에서 가장 먼 셀일 뿐, 유입을 모른다. 9회는 시작이 흐름 축 모서리에 걸린 경우다. 걸린다는 보장은 없다. 재정렬 전후로 해가 기계 정밀도에서 달라지면 인덱스 변환이 빠진 것이다. 번호는 해를 바꾸면 안 된다. 1D에서 왼쪽부터 매기면 끝나던 분해는 Burgers 예제의 LU-SGS에 있다. 그 번호 문제는 여기서 시작된다.