연립방정식의 직접 해법

가우스 소거의 연산 수

열 우선 순서로 쓴 LU 분해의 핵심은 세 겹의 반복문입니다. $k$번째 단계에서 피벗 아래의 $n - k$개 행마다 승수를 하나 구하고(나눗셈), 남은 $(n - k)\times(n - k)$ 부분 행렬을 갱신합니다(곱셈과 뺄셈).

for k = 1 … n−1
  for i = k+1 … n
    a[i,k] ← a[i,k] / a[k,k]              // 승수 ℓᵢₖ를 그 자리에 저장
    for j = k+1 … n
      a[i,j] ← a[i,j] − a[i,k]·a[k,j]    // 슈어 보수 갱신: 곱 1 + 뺄셈 1
// 끝나면 a의 아래 삼각 = L(대각 1 생략), 위 삼각 = U

그래서 문제 상자의 답은 $10^3$배, 약 100초입니다. 메모리는 $n^2$이라 100배(배정밀도 $10^4\times10^4$ 행렬은 800 MB)입니다. 실제 라이브러리(LAPACK의 dgetrf)는 같은 연산을 블록 단위로 묶어 행렬–행렬 곱(BLAS 3단계)으로 바꿉니다. 연산 수는 같지만 캐시에 올린 데이터를 여러 번 재사용하므로, 위의 단순한 세 겹 반복보다 수십 배 빠르게 돕니다.

표의 넷째 줄이 "역행렬을 계산하지 말라"는 격언의 근거입니다. $A^{-1}\mathbf b$를 구하려고 역행렬을 만들면 분해만 할 때의 세 배를 쓰고, 정확도도 대개 더 나쁩니다. 우변이 $m$개 있으면 분해는 한 번만 하고 대입을 $m$번 하면 됩니다($\tfrac23n^3 + 2mn^2$).

피벗팅: 작은 피벗과 성장 인자

부분 피벗팅은 승수를 $|\ell_{ik}| \le 1$로 묶습니다. 그런데 행 하나에 큰 수를 곱해 두면 이 장치가 속습니다. $10x + 10^5y = 10^5$, $x + y = 2$(첫 식에 $10^5$를 곱한 것과 같은 계)에서 부분 피벗팅은 10을 피벗으로 골라 3자리 계산에서 다시 $x = 0$을 냅니다. 행마다 가장 큰 원소 $s_i = \max_j|a_{ij}|$로 나눈 상대 크기 $|a_{ik}|/s_i$로 피벗을 고르면($10/10^5$ 대 $1/1$) 둘째 행이 선택되어 $x = y = 1.00$이 나옵니다. 이것이 척도 부분 피벗팅(scaled partial pivoting)입니다. 실무에서는 행렬을 미리 평형화(행·열 척도 조정)하고 부분 피벗팅을 씁니다.

상한에 닿는 예가 윌킨슨의 행렬입니다. 대각이 1, 대각 아래가 모두 $-1$, 마지막 열이 모두 1인 행렬은 부분 피벗팅이 행을 바꾸지 않고(첫 열의 최댓값이 이미 대각에 있음), 마지막 열이 단계마다 두 배가 되어 $u_{nn} = 2^{n-1}$입니다. 조건수는 $n$ 정도로 작으니 문제는 전혀 예민하지 않은데, 알고리즘이 답을 망칩니다.

표는 참해를 무작위로 정하고 $\mathbf b = A\mathbf x$로 만든 실험입니다. 그런데 이런 행렬은 실제 응용에서 거의 나오지 않습니다. 무작위 50×50 행렬 200개에 부분 피벗팅을 쓰면 성장 인자의 중앙값은 7.5, 최댓값은 15.8이었습니다. 이 경험적 사실 때문에 LAPACK의 일반 풀이기는 부분 피벗팅을 기본으로 쓰고, 완전 피벗팅은 거의 쓰지 않습니다.

PA = LU와 여러 우변

부분 피벗팅을 쓴 소거는 행 교환을 기록하는 치환 행렬 $P$와 함께 $PA = LU$를 만듭니다(LinearBook 4장의 정리 4.3). 구현에서는 $P$를 행렬로 저장하지 않고 정수 배열 piv로 두며, 승수 $\ell_{ik}$를 소거로 0이 된 칸에 덮어써서 $L$과 $U$를 원래 배열 하나에 담습니다. 행을 바꿀 때는 이미 저장한 승수도 함께 바뀐다는 것이 손계산의 함정입니다.

촐레스키 분해: 대칭 양정치 행렬

대칭이지만 양정치가 아닌 행렬(구조물의 좌굴 해석, 제약 조건이 있는 최적화의 KKT 행렬)에는 대칭성을 살리면서 1×1과 2×2 블록 피벗을 쓰는 $PAP^{\mathsf T} = LDL^{\mathsf T}$ 분해(번치–카우프만 피벗팅)를 씁니다. LinearBook 4장의 $LDL^{\mathsf T}$에 피벗팅을 더한 것입니다.

띠 행렬과 토머스 알고리즘

삼중대각계 $a_ix_{i-1} + b_ix_i + c_ix_{i+1} = d_i$($a_1 = c_n = 0$)에 피벗팅 없는 소거를 쓰면 $L$과 $U$도 두 줄짜리 띠로 남습니다. 행마다 위쪽 행을 한 번 빼서 정규화한 계수 $c_i'$, $d_i'$를 만들고(전진 소거), 아래에서 위로 대입합니다.

function thomas(a, b, c, d)        // 길이 n, a[1] = c[n] = 0
  c'[1] ← c[1] / b[1];   d'[1] ← d[1] / b[1]
  for i = 2 … n                     // 전진 소거: 행마다 상수 번
    m ← b[i] − a[i]·c'[i−1]           // 이번 행의 피벗
    c'[i] ← c[i] / m
    d'[i] ← (d[i] − a[i]·d'[i−1]) / m
  x[n] ← d'[n]
  for i = n−1 … 1                   // 후진 대입
    x[i] ← d'[i] − c'[i]·x[i+1]
  return x                          // 총 약 8n 연산, 메모리 O(n)

반띠폭이 $w$인 띠 행렬로 일반화하면 연산 $O(nw^2)$, 저장 $O(nw)$입니다(LinearBook 정리 4.5). 부분 피벗팅을 쓰면 $U$의 위쪽 띠가 $2w$까지 넓어질 수 있으므로 LAPACK의 띠 풀이기(dgbsv)는 그만큼 여분의 자리를 마련합니다.

노름, 조건수, 그리고 오차 한계

정리 12.3과 합치면 실용적인 경험 법칙이 나옵니다. 후방 안정한 풀이기는 $\|\delta A\|/\|A\| \approx u \approx 10^{-16}$인 문제를 정확히 푸는 셈이므로, 상대 오차 ≈ κ·u, 즉 대략 $\log_{10}\kappa$자리를 잃는다. κ = 10⁶이면 16자리 중 약 10자리가 남습니다. LAPACK은 κ를 직접 계산하지 않고($A^{-1}$이 필요하므로) LU 분해를 재사용해 $O(n^2)$에 $\|A^{-1}\|_1$을 추정합니다(dgecon).

힐베르트 행렬: 깔끔한 분수, 사라지는 자릿수

$H_{ij} = 1/(i + j - 1)$인 힐베르트 행렬(Hilbert matrix)은 구간 $[0, 1]$에서 단항식 $1, x, \dots, x^{n-1}$으로 최소제곱 근사를 할 때 나오는 그람 행렬입니다(7장). 단항식들이 $[0, 1]$에서 서로 거의 같은 모양이라(높은 차수일수록 $x^k$과 $x^{k+1}$은 구별하기 어렵다) 열벡터들이 거의 일차 종속이 됩니다. 역행렬은 정수 원소를 가지는 공식이 알려져 있어 조건수를 정확히 계산할 수 있습니다.

표가 정리 12.8을 그대로 보여 줍니다. 잔차는 n과 무관하게 반올림 수준(후방 안정)인데, 오차는 κ·u를 따라 커져 n = 12에서 한 자리도 남지 않습니다. 단, 표의 $\mathbf b = H\mathbf 1$ 자체도 반올림된 값이고 $H$의 원소($1/3$ 등)도 정확히 저장되지 않으므로, 이 오차의 상당 부분은 풀이기 이전에 이미 입력 단계에서 생긴 것입니다. 어떤 알고리즘도 되돌릴 수 없는 몫입니다.

반복 개선: 싼 분해로 비싼 정확도를

근사해 $\hat{\mathbf x}$의 오차 $\mathbf e = \mathbf x - \hat{\mathbf x}$는 정확히 $A\mathbf e = \mathbf r$을 만족합니다. 이 보정 방정식을 이미 있는 (부정확한) LU로 풀어 더하면, 보정값의 상대 오차가 처음 해와 같은 수준이어도 보정 자체가 작으므로 해는 한 단계 더 정확해집니다. 핵심은 잔차 $\mathbf r = \mathbf b - A\hat{\mathbf x}$를 높은 정밀도로 계산하는 것입니다. 거의 같은 두 수의 뺄셈이라 여기서 정밀도를 아끼면 자릿수 상쇄(1장)로 정보가 사라집니다.

function mixed_refine(A, b)
  (L, U, P) ← lu(single(A))              // O(n³): 단정밀도, 싸다
  x ← solve_single(L, U, P, b)            // O(n²)
  repeat
    r ← b − A·x                           // O(n²): 배정밀도 (원래 A로)
    d ← solve_single(L, U, P, r)          // 같은 분해 재사용
    x ← x + d                             // 배정밀도로 더하기
  until ‖d‖ ≤ 작은 수·‖x‖
  return x

표는 특이값이 $1$부터 $1/\kappa$까지 등비로 퍼진 행렬로 만든 실험입니다(상대 오차, ∞-노름). κ = 10⁴에서는 단 두 번 만에 오차가 $6.8\times10^{-5}$에서 배정밀도 바닥 $\kappa u \approx 10^{-12}$ 근처로 떨어집니다. 단정밀도의 $u_f \approx 6\times10^{-8}$이므로 κ = 10⁸은 $\kappa u_f \approx 6$으로 정리의 조건 밖입니다. 이 실험에서는 반복마다 약 한 자리씩 느리게 줄었지만, 다른 행렬에서는 정체하거나 발산할 수 있습니다. 최근의 혼합 정밀도 풀이기는 보정 방정식을 GMRES(13장)로 풀어 이 한계를 넓힙니다.

희소 행렬과 채움, 재배열

희소 행렬은 0 아닌 원소만 저장합니다. 가장 흔한 형식 CSR(compressed sparse row)은 값 배열, 열 번호 배열, 각 행의 시작 위치 배열 세 개로 행렬을 담아 저장량이 원소 수 $\mathrm{nnz}$에 비례합니다. 문제는 분해입니다. 소거 중 $a_{ij} \leftarrow a_{ij} - \ell_{ik}u_{kj}$에서 $a_{ij} = 0$이어도 $\ell_{ik}$와 $u_{kj}$가 둘 다 0이 아니면 새 원소가 생깁니다.

대표적인 재배열은 셋입니다. 역 커트힐–맥키(reverse Cuthill–McKee, RCM)는 너비 우선 탐색으로 번호를 매겨 띠폭을 줄입니다. 최소 차수(minimum degree)는 매 단계 소거 그래프에서 차수가 가장 작은 정점을 고릅니다(실무에서는 근사 최소 차수 AMD). 중첩 분할(nested dissection)은 격자를 분리자(separator) 한 줄로 반씩 나누고, 분리자를 맨 마지막에 번호 매기는 일을 재귀적으로 반복합니다. 2차원 $k\times k$ 격자에서 중첩 분할의 채움은 $O(n\log n)$, 연산은 $O(n^{3/2})$($n = k^2$)로 점근적으로 최적이라는 것이 알려져 있습니다.

표는 소거 그래프로 기호적 채움을 센 값입니다(숫자 상쇄는 무시). 문제 상자의 답은 자연 순서 8,019칸, 최소 차수 3,729칸입니다. 이 크기에서는 탐욕적 최소 차수가 가장 적지만, 격자가 커질수록 중첩 분할의 점근적 우위가 드러나고 병렬화에도 유리해 대규모 FEM 풀이기(MUMPS, PARDISO 등)는 두 방법을 섞어 씁니다. 채움을 아무리 줄여도 3차원 문제에서는 $L$이 너무 커지는 경우가 많아, 그때는 반복 해법(13장)으로 넘어갑니다.

이 도구가 쓰이는 곳