반복 해법

큰 희소 행렬: 왜 분해하지 않는가

격자점을 행 순서로 번호 매기면 $n = N^2$차 행렬 $A$의 각 행에는 대각 원소 4와 이웃 네 개의 $-1$만 있습니다. 0이 아닌 원소 수(nnz)는 약 $5n$입니다. 이런 행렬을 희소 행렬(sparse matrix)이라 합니다. 문제는 분해입니다. 위아래 이웃이 $N$칸 떨어져 있어 띠 폭이 $N$이고, 촐레스키 인수는 띠 안을 거의 다 채웁니다(채움(fill-in), 12장).

반복법은 $A$를 바꾸지 않습니다. 필요한 것은 "벡터 $\mathbf v$를 주면 $A\mathbf v$를 돌려주는" 함수뿐이고, 한 번의 비용은 nnz에 비례합니다. 그래서 반복 횟수를 $k$라 하면 총비용은 $k\cdot\mathrm{nnz}$이고, 이 장의 모든 아이디어는 결국 $k$를 줄이는 기술입니다. 같은 문제에서 야코비는 수백만 번, 최적 SOR과 CG는 수천 번, 다중 격자는 열 번 안팎의 V-사이클이 필요하다는 것을 이 장 끝에서 보게 됩니다.

야코비와 가우스–자이델

두 가지 선택이 있습니다. 한 바퀴를 도는 동안 이전 값만 쓰는 것이 야코비 방법(Jacobi method), 방금 계산한 $x_1^{(k+1)}, \dots, x_{i-1}^{(k+1)}$을 즉시 쓰는 것이 가우스–자이델 방법(Gauss–Seidel method)입니다.

function gauss_seidel(A, b, x, tol, maxit)
  for k = 1 .. maxit
    for i = 1 .. n                          // 희소 행렬이면 i행의 0 아닌 원소만
      s = b[i]
      for j in nonzeros(row i), j != i
        s = s - A[i][j] * x[j]               // j < i 는 이미 새 값
      x[i] = s / A[i][i]                     // 야코비는 여기서 xnew[i]에 쓴다
    if norm(b - A x) <= tol * norm(b) then return x
  return x                                  // 수렴 실패: 호출자가 확인

2 × 2에서 기하가 분명해집니다. 두 방정식은 평면의 두 직선이고 해는 교점입니다. 가우스–자이델의 첫 갱신은 $x_2$를 고정한 채 $x_1$을 바꿔 첫째 직선 위로 수평 이동하고, 둘째 갱신은 둘째 직선 위로 수직 이동합니다. 계단을 타고 교점으로 내려가는 셈입니다. 야코비는 두 좌표를 동시에 바꾸므로 두 직선 어느 쪽에도 놓이지 않은 점으로 대각선 이동합니다.

분할과 반복 행렬, 스펙트럼 반지름

세 방법을 한 틀에 넣습니다. $A$를 쉽게 풀리는 부분 $M$과 나머지로 분할해 $A = M - N$이라 쓰면 $A\mathbf x = \mathbf b$는 $M\mathbf x = N\mathbf x + \mathbf b$이고, 이를 반복식으로 읽습니다.

오른쪽 꼴이 핵심입니다. 모든 정상 반복법은 "잔차 $\mathbf r_k = \mathbf b - A\mathbf x_k$를 $M^{-1}$로 근사 풀이해 더한다"는 한 줄이고, $M$이 $A$에 가까울수록 한 번에 많이 나아갑니다. 오차 $\mathbf e_k = \mathbf x_k - \mathbf x$는 $\mathbf b$와 무관하게 $\mathbf e_{k+1} = G\mathbf e_k$, 곧 $\mathbf e_k = G^k\mathbf e_0$를 따릅니다.

필요충분조건은 고유값을 알아야 쓸 수 있습니다. 실무에서는 행렬을 보고 바로 판정할 수 있는 충분조건이 더 유용합니다.

SOR과 최적 이완 계수

축차 과이완(successive over-relaxation, SOR)은 가우스–자이델 값 $\tilde x_i$와 이전 값 $x_i^{(k)}$의 가중 평균을 씁니다.

놀랍게도 많은 실용 행렬에서 $\rho(G_\omega)$를 $\omega$의 함수로 정확히 쓸 수 있습니다. 조건은 행렬이 일관된 순서(consistently ordered)를 갖는 것인데, 삼중대각 행렬과 5점 격자 행렬(자연 순서나 적흑 순서)이 여기에 속합니다.

2차원 푸아송 문제(한 변 $N$점, $h = 1/(N+1)$)에서는 $\rho_J = \cos\pi h$이므로 $\sqrt{1 - \rho_J^2} = \sin\pi h$이고

1에서 떨어진 거리가 $h^2$에서 $h$로 바뀝니다. 오차를 일정 배 줄이는 반복 수가 $O(N^2)$에서 $O(N)$으로 줄어든다는 뜻입니다.

이차 형식과 최급 강하

$\phi$의 등고선은 $A$의 고유벡터를 축으로 하는 타원이고, 긴 축과 짧은 축의 길이 비는 $\sqrt\kappa$($\kappa = \lambda_{\max}/\lambda_{\min}$)입니다. 잔차 $\mathbf r$은 가장 가파른 내리막 방향입니다. 그 방향의 직선 $\mathbf x + \alpha\mathbf r$ 위에서 $\phi$를 최소로 하는 $\alpha$는 $\frac{d}{d\alpha}\phi = -\mathbf r^{\mathsf T}\mathbf r + \alpha\,\mathbf r^{\mathsf T}A\mathbf r = 0$에서 나옵니다.

function steepest_descent(A, b, x, tol)
  r = b - A x
  while norm(r) > tol * norm(b)
    q = A r                                // 반복마다 행렬–벡터 곱 한 번
    alpha = dot(r, r) / dot(r, q)          // 직선 위의 정확한 최솟점
    x = x + alpha * r
    r = r - alpha * q                      // b - A x 를 다시 계산하지 않는다
  return x

켤레 기울기법

지그재그의 원인은 새 방향으로 최소화하면서 이전 방향에서 얻은 최소성을 망치는 데 있습니다. 해결책은 방향들을 $A$에 대해 서로 "수직"으로 고르는 것입니다. $\mathbf p_i^{\mathsf T}A\mathbf p_j = 0$($i \ne j$)이면 $\mathbf p_i, \mathbf p_j$를 A-켤레(A-conjugate)라 합니다. $A = R^{\mathsf T}R$로 좌표를 $\mathbf y = R\mathbf x$로 바꾸면 타원 등고선이 원이 되고, A-켤레 방향은 그 좌표에서 평범한 직교 방향입니다. 원 위에서 서로 수직인 방향으로 한 번씩 최소화하면 $n$번 만에 중심에 닿습니다.

A-켤레 방향을 미리 구하는 것은 비쌉니다(그람–슈미트에 모든 이전 방향이 필요). 켤레 기울기법(conjugate gradient, CG)의 기적은 새 잔차를 바로 직전 방향 하나에 대해서만 A-켤레로 만들어도 모든 이전 방향과 저절로 A-켤레가 된다는 것입니다.

function cg(A, b, x, tol)
  r = b - A x;  p = r;  rr = dot(r, r)
  while sqrt(rr) > tol * norm(b)
    q = A p
    alpha = rr / dot(p, q)                 // p 방향의 정확한 최소화
    x = x + alpha * p
    r = r - alpha * q
    rr_new = dot(r, r)
    beta = rr_new / rr                     // 새 방향을 p와 A-켤레로
    p = r + beta * p
    rr = rr_new
  return x                               // 반복당: 곱 1번, 내적 2번, 벡터 갱신 3번

최소성에서 수렴 속도가 나옵니다. $\mathbf x_k - \mathbf x^* = P_k(A)(\mathbf x_0 - \mathbf x^*)$이고 $P_k$는 $P_k(0) = 1$인 $k$차 다항식 중 무엇이든 될 수 있으므로, 고유벡터로 전개하면

최급 강하의 $(\kappa - 1)/(\kappa + 1) \approx 1 - 2/\kappa$가 CG에서는 $1 - 2/\sqrt\kappa$로 바뀝니다. $\kappa = 10^4$이면 오차를 10배 줄이는 데 최급 강하는 약 11500회, CG는 약 115회입니다. 푸아송 문제의 $\kappa \approx 4/(\pi h)^2$를 넣으면 CG의 반복 수는 $O(N)$, 최적 SOR과 같은 차수이면서 $\omega$를 고를 필요가 없습니다.

2차원에서는 CG가 언제나 두 걸음이니 차이가 극적이지만, 실제로 관심 있는 것은 $n$이 클 때 $n$보다 훨씬 적은 반복으로 충분히 정확해지는가입니다. 정리 13.9의 두 한계는 같은 $\kappa$라도 고유값이 어떻게 흩어져 있느냐에 따라 CG가 아주 다르게 행동한다는 것을 말해 줍니다.

전처리: 문제를 풀기 쉬운 꼴로 바꾸기

$A$와 비슷하면서 $M\mathbf z = \mathbf r$이 싸게 풀리는 SPD 행렬 $M$을 전처리기(preconditioner)라 합니다. $M^{-1}A\mathbf x = M^{-1}\mathbf b$는 같은 해를 갖고, $M \approx A$이면 $M^{-1}A \approx I$라 고유값이 1 근처에 모입니다. 대칭을 유지하려면 $M = LL^{\mathsf T}$로 쓰고 $L^{-1}AL^{-\mathsf T}$에 CG를 적용하는데, 식을 정리하면 $L$이 사라지고 매 반복 $M\mathbf z = \mathbf r$ 하나만 더 푸는 전처리 CG(PCG)가 됩니다.

function pcg(A, b, x, M, tol)
  r = b - A x;  z = solve(M, r);  p = z;  rz = dot(r, z)
  while norm(r) > tol * norm(b)
    q = A p
    alpha = rz / dot(p, q)
    x = x + alpha * p;  r = r - alpha * q
    z = solve(M, r)                        // 전처리: M z = r (싸게!)
    rz_new = dot(r, z);  beta = rz_new / rz
    p = z + beta * p;  rz = rz_new
  return x

$M = A$면 한 번에 끝나지만 $M$을 푸는 것이 원래 문제이고, $M = I$면 공짜지만 효과가 없습니다. 쓸모 있는 전처리기는 그 사이의 절충입니다.

표가 전처리의 두 얼굴을 보여 줍니다. 균일한 푸아송 문제에서 야코비 전처리는 대각이 모두 4라 아무 일도 하지 않고, IC(0)는 반복을 절반 아래로 줄이지만 $N$에 비례해 느는 성질은 그대로입니다. 계수가 1000배 뛰는 문제에서는 야코비만으로도 반복이 7–11분의 1이 됩니다. 반복 수가 $N$과 무관해지는 전처리기는 다음 절의 다중 격자에서 나옵니다.

크릴로프 부분공간과 GMRES

CG가 한 일을 다시 보면, $\mathbf x_k$를 크릴로프 공간 $\mathbf x_0 + \mathcal K_k(A, \mathbf r_0)$에서 골랐습니다. 이 공간은 행렬–벡터 곱만으로 만들어지므로 비대칭이어도 쓸 수 있습니다. 바뀌는 것은 "가장 좋은 점"의 기준입니다. A-노름은 SPD에서만 노름이므로, 대신 잔차의 2-노름을 최소화합니다. 이것이 GMRES(generalized minimal residual)입니다.

이를 위해 먼저 크릴로프 공간의 정규직교 기저를 만듭니다. $\mathbf r_0, A\mathbf r_0, A^2\mathbf r_0, \dots$를 그대로 쓰면 거듭제곱법처럼 모두 지배 고유벡터 쪽으로 쏠려 거의 평행해지므로, 새 벡터를 만들 때마다 그람–슈미트로 직교화합니다. 이것이 아놀디 과정(Arnoldi iteration)입니다.

function gmres(A, b, x0, m, tol)          // m: 재시작 길이
  r = b - A x0;  beta = norm(r);  v[1] = r / beta
  for j = 1 .. m                           // 아놀디: A V_m = V_(m+1) H
    w = A v[j]
    for i = 1 .. j
      H[i][j] = dot(w, v[i]);  w = w - H[i][j] * v[i]
    H[j+1][j] = norm(w);  v[j+1] = w / H[j+1][j]
    // 작은 최소제곱 문제: min_y ‖beta e1 - H y‖, 기븐스 회전으로 갱신
    if |잔차 추정| <= tol * norm(b) then break
  x = x0 + V y
  if 수렴하지 않았으면 x0 = x 로 다시 시작   // GMRES(m)

비용이 계속 커지므로 실무에서는 $m$번마다 지금까지의 해를 새 출발점으로 삼아 다시 시작합니다. 이 재시작 GMRES(m)은 메모리를 $m$개 벡터로 묶는 대가로 크릴로프 공간에 쌓인 정보를 버리므로 더 많은 반복이 필요하고, 때로는 정체합니다. 대칭 행렬에서는 아놀디의 $H$가 삼중대각이 되어 직교화가 직전 두 벡터로 충분하며(란초시 과정), 이것이 CG가 짧은 점화식으로 끝나는 이유입니다. 비대칭에서 짧은 점화식을 원하면 BiCGSTAB 같은 방법을 쓰지만 최소성은 포기합니다.

다중 격자: 매끄러운 오차는 거친 격자에서

1차원 모델 $-u'' = f$($N$점, $h = 1/(N+1)$)에서 오차를 사인 모드 $\sin(k\pi x)$로 전개해 봅시다. 가중 야코비 $\mathbf x \leftarrow \mathbf x + \omega D^{-1}\mathbf r$는 각 모드를 독립적으로 줄이며, 그 배율은

그래서 완화 반복은 해법으로는 느리지만 평활화(smoothing)로는 훌륭합니다. 몇 번만 돌리면 고주파 오차가 사라지고 매끄러운 오차가 남습니다. 매끄러운 오차는 격자 간격을 두 배로 한 거친 격자에서도 잘 표현되며, 거기서는 같은 모드가 상대적으로 두 배 높은 주파수라 평활화가 다시 효과적입니다. 이 관찰을 재귀로 엮은 것이 다중 격자법(multigrid)입니다.

function vcycle(h, u, f)
  if 격자가 가장 거칠면 return solve_exact(A_h, f)
  u = smooth(A_h, u, f, nu)                // 가우스–자이델 2번: 고주파 제거
  r = f - A_h u
  rc = restrict(r)                         // 2h 격자로 (가중 평균)
  ec = vcycle(2h, 0, rc)                   // 오차 방정식 A_2h e = rc 를 재귀로
  u = u + prolong(ec)                      // h 격자로 (선형 보간)
  u = smooth(A_h, u, f, nu)
  return u

마지막으로 판의 온도 문제에서 이 장의 방법을 모두 경주시켜 봅니다. 해만 보면 차이가 잘 보이지 않지만 오차를 보면 각 방법의 성격이 드러납니다. 가우스–자이델의 오차는 매끄러운 언덕으로 남아 천천히 가라앉고, 다중 격자의 오차는 사이클마다 전체가 한 자릿수씩 내려갑니다.

이 도구가 쓰이는 곳

큰 희소 선형계는 편미분방정식, 그래프, 영상, 최적화의 내부 어디에나 있습니다. 어떤 반복법을 고를지는 행렬의 성질이 정합니다.