가우스 소거법

1. 개요

연립방정식? 그냥 위에서 아래로 하나씩 소거하면 되는 거 아니었어? 응, 맞아. 근데 그게 O(n^3)이라서 문제야.

가우스 소거법(Gaussian elimination)은 선형 연립방정식 $A\mathbf{x} = \mathbf{b}$를 풀기 위해 계수행렬을 기본 행 연산(elementary row operation)으로 상삼각행렬(upper triangular matrix) 꼴로 만든 뒤, 후진대입(back substitution)으로 해를 구하는 대표적인 직접법(direct method)이다. 중학교 때 배운 “가감법”을 행렬로 체계화한 것이라고 보면 되는데, 이 단순해 보이는 알고리즘이 수치해석에서 희소행렬 이전 시대를 지배한 만능키였다.

이름은 카를 프리드리히 가우스(Carl Friedrich Gauss)에게서 왔지만, 소거 자체는 훨씬 오래됐다. 한나라 때 편찬된 중국 수학서 구장산술(九章算術) 제8장 “방정(方程)“이 이미 같은 절차를 표로 기술하고 있다.¹ 가우스가 남긴 것은 방법이 아니라 표기법이었다 — 1810년 소행성 팔라스의 궤도를 최소제곱으로 맞추는 계산에서 그가 쓴 대괄호 기호가 19세기 직업 계산수(human computer)들 사이에 퍼지면서, 르장드르가 “보통의(ordinary)” 소거라 부르던 것이 “가우스” 소거법이 되었다.² 어쨌든 오늘날 거의 모든 유한요소법·유한체적법 솔버의 밑바닥 어딘가에는 이 소거 과정이 LU 분해의 탈을 쓰고 돌아가고 있다.

2. 알고리즘

$n \times n$ 행렬에 대한 가우스 소거법은 두 단계로 나뉜다. 첫째는 전진 소거(forward elimination), 둘째는 후진 대입(back substitution)이다.

전진 소거에서는 $k$번째 열 아래쪽 성분들을 0으로 만든다. $k$번째 행을 기준(pivot row)으로 삼아 $i > k$인 각 행에서 배수를 빼주는데, 그 배수(multiplier)는 다음과 같다.

$$m_{ik} = \frac{a_{ik}}{a_{kk}}, \qquad a_{ij} \leftarrow a_{ij} - m_{ik}\,a_{kj}$$

이 과정을 $k = 1$부터 $n-1$까지 반복하면 행렬이 상삼각행렬 $U$가 된다. 우변 벡터 $\mathbf{b}$도 같이 갱신해준다.

전진 소거가 끝나면 마지막 방정식은 미지수가 하나뿐이므로 바로 풀 수 있다. 이걸 위로 거슬러 올라가며 대입하는 것이 후진 대입이다.

$$x_i = \frac{1}{u_{ii}}\left( b_i - \sum_{j=i+1}^{n} u_{ij}\,x_j \right), \qquad i = n, n-1, \dots, 1$$

여기서 눈치 빠른 사람은 알아챘을 것이다. 저 배수 $m_{ik}$ 값들을 버리지 않고 하삼각 자리에 차곡차곡 저장하면 그게 바로 LU 분해다. 가우스 소거법과 LU 분해는 사실상 같은 계산의 두 얼굴이고, 그래서 실제 라이브러리에는 “가우스 소거법”이라는 이름의 함수가 없다 — LAPACK에서 이 일을 하는 루틴 이름은 dgetrf(분해)와 dgetrs(삼각 풀이)다.

3. 부분 피벗팅

순진하게 위 공식을 그대로 돌리면 언젠가 큰코다친다. 피벗 $a_{kk}$가 0이면 0으로 나누는 참사가 벌어지고, 0은 아니어도 아주 작으면 배수 $m_{ik}$가 폭발해서 반올림 오차(round-off error)가 걷잡을 수 없이 증폭된다.³ 그래서 실무에서는 반드시 부분 피벗팅(partial pivoting)을 쓴다.

부분 피벗팅은 $k$번째 소거 단계에서 $k$열의 대각 아래 성분 중 절댓값이 가장 큰 놈을 찾아 그 행을 피벗 행으로 끌어올리는 것이다.

$$p = \arg\max_{i \ge k} |a_{ik}|$$

이렇게 하면 모든 배수가 $|m_{ik}| \le 1$을 만족해 소거 중에 값이 튀는 폭이 제한된다. 행 교환 정보는 치환행렬(permutation matrix) $P$로 기록되며, 그래서 정확한 표현은 $PA = LU$가 된다. 완전 피벗팅(complete pivoting)은 열까지 교환해 이론적으로 더 안정하지만, 매 단계 남은 부분행렬 전체를 뒤져야 해서 비교 횟수가 $O(n^3)$으로 늘고 BLAS 3 수준의 블록 갱신과도 궁합이 나쁘다. 그래서 실전에서는 거의 안 쓰인다. 부분 피벗팅으로도 대부분의 실제 문제는 충분히 안정적으로 풀린다는 게 검증 및 확인을 거친 국룰이다.

4. 계산복잡도

가우스 소거법의 악명 높은 지점이 바로 여기다. 단계 $k$에서는 배수 $n-k$개를 만든 뒤, 남은 부분행렬 $(n-k) \times (n-k)$와 우변의 $n-k$개 성분을 각각 “곱하고 빼는” 한 쌍으로 갱신한다. 곱셈·나눗셈과 덧셈·뺄셈을 각각 1 flop으로 세는 관례를 따르면

$$\sum_{k=1}^{n-1} 2(n-k)(n-k+1) = \frac{2}{3}(n-1)n(n+1) \approx \frac{2}{3}n^3$$

이다. 배수를 만드는 나눗셈은 다 합쳐도 $O(n^2)$이라 선두항에 묻힌다. 즉 지배적인 항이 $O(n^3)$이고, 관례적으로 인용되는 숫자가 이 $\tfrac{2}{3}n^3$ flops다. 후진 대입은 $2n^2$ flops뿐이라 소거 앞에서는 반올림 오차 수준으로 무시된다. 이 세제곱이 왜 무서운가 하면, 미지수가 10배 늘면 계산량은 1000배 늘기 때문이다.

미지수 $n$	대략 flops ($\tfrac{2}{3}n^3$)	체감
$10^3$	약 $6.7 \times 10^8$	눈 깜짝
$10^4$	약 $6.7 \times 10^{11}$	커피 한 잔
$10^6$	약 $6.7 \times 10^{17}$	인생 끝날 때까지

메모리도 같이 죽는다. $n = 10^6$짜리 밀집 행렬은 배정밀도로 $8 \times 10^{12}$바이트, 즉 8 TB다. 연산이 아니라 저장부터 이미 불가능하다.

그래서 전산유체역학이나 구조해석처럼 미지수가 수백만 개인 실제 문제에서는 밀집(dense) 가우스 소거법을 통째로 돌리는 짓은 절대 하지 않는다. 대신 계수행렬이 희소행렬이라는 점을 이용해 희소 직접 솔버로 가거나, 아예 직접법을 버리고 반복법과 크리로프 부분공간법으로 갈아탄다. 밀집 소거법은 어디까지나 수백~수천 규모의 부분 문제나 전처리기 내부, 그리고 희소 솔버의 초노드(supernode) 블록 안에서 조용히 쓰인다.

5. 후방 안정성과 성장 인자

가우스 소거법의 오차 거동은 성장 인자(growth factor)로 분석한다. 소거 도중 등장한 모든 중간 성분의 최대 절댓값을 원래 행렬의 최대 절댓값으로 나눈 값이다.

$$\rho_n = \frac{\max_{i,j,k} |a_{ij}^{(k)}|}{\max_{i,j} |a_{ij}|}$$

부동소수점으로 계산한 해 $\hat{\mathbf{x}}$는 원래 문제가 아니라 살짝 흔들린 문제의 정확한 해임이 보장된다. 즉 $(A + \Delta A)\hat{\mathbf{x}} = \mathbf{b}$를 만족하는 $\Delta A$가 존재하고, 그 크기는

$$\frac{\|\Delta A\|_\infty}{\|A\|_\infty} \lesssim c_n \, \rho_n \, \varepsilon_{\text{mach}}$$

로 묶인다($c_n$은 $n$의 낮은 차수 다항식, 배정밀도에서 $\varepsilon_{\text{mach}} = 2^{-53} \approx 1.1 \times 10^{-16}$). 이런 성질을 후방 안정(backward stable)이라 부른다. 핵심은 후방 안정성이 곧 정확한 답을 뜻하지는 않는다는 것이다. 해의 상대오차는 여기에 조건수가 한 번 더 곱해져

$$\frac{\|\hat{\mathbf{x}} - \mathbf{x}\|}{\|\mathbf{x}\|} \lesssim \kappa(A) \, c_n \, \rho_n \, \varepsilon_{\text{mach}}$$

가 된다. $\kappa(A)$가 $10^{16}$ 근처면 배정밀도의 유효숫자가 전부 날아간다는 뜻이고, 이건 알고리즘 탓이 아니라 문제 자체가 삐뚤어진(ill-conditioned) 것이다. 소거법을 아무리 정성껏 돌려도 답이 없고, 전처리기로 문제를 다듬거나 정식화 자체를 다시 봐야 한다.

남는 변수는 $\rho_n$인데, 부분 피벗팅에서는 이론 상한이 $\rho_n \le 2^{n-1}$이고 이 값을 실제로 달성하는 윌킨슨(Wilkinson) 행렬이 존재한다. $n = 60$이면 $2^{59} \approx 5.8 \times 10^{17}$, 즉 배정밀도가 통째로 증발하는 크기다. 완전 피벗팅은 상한이 $n$의 다항식보다 살짝 빠르게 자라는 정도로 훨씬 얌전하지만, “상한이 딱 $n$일 것”이라던 윌킨슨의 추측은 틀렸다 — 1991년 굴드(N. Gould)가 $13 \times 13$ 행렬에서 $\rho = 13.0205$를 보고했고, 이듬해 에델만(A. Edelman)이 정확 연산으로 반례임을 확정했다.⁴

6. 실전에서 정말 터지는가

“부분 피벗팅은 상한이 지수인데 왜 다들 마음 놓고 쓰는가”는 수치해석의 오래된 질문이었다. 무작위 행렬에서는 $\rho_n$이 $n^{2/3}$ 수준에 머물러 통계적으로 안전하다는 것이 정설이지만, 실무 문제에서 터진 사례가 없다는 말은 사실이 아니다. 라이트(Wright, 1993)는 상미분방정식의 2점 경계값 문제를 표준 방식으로 이산화하면 행마다 비영 성분이 3개뿐인 성긴 띠행렬인데도 성장 인자가 차원에 대해 지수적으로 커지는 예를 만들었다. 윌킨슨 행렬처럼 인위적으로 꾸민 것이 아니라, 평범한 경계값 문제 풀이에서 자연히 나오는 행렬이라는 점이 요점이다. 볼테라 적분방정식 이산화에서도 비슷한 보고가 있다.

실무적 결론은 이렇다. 부분 피벗팅은 기본값으로 두되, 믿지 말고 재라. 대부분의 라이브러리는 잔차 $\|\mathbf{b} - A\hat{\mathbf{x}}\|$나 조건수 추정치($\kappa$ 추정 루틴)를 같이 뱉어주고, 그 두 숫자가 위 부등식의 $\rho_n$과 $\kappa(A)$ 자리를 사후에 채워준다. 잔차는 작은데 답이 이상하면 조건수를, 잔차마저 크면 피벗팅·성장 인자를 의심하는 순서다. 필요하면 반복 정련(iterative refinement)으로 한두 자리를 되찾을 수도 있다.

7. 관련 문서

• LU 분해
• 촐레스키 분해
• 반복법
• 희소행렬 · 희소 직접 솔버
• 전처리기
• 크리로프 부분공간법
• 조건수
• 부동소수점 연산
• LAPACK
• 수치해석
• 고유값 문제
• 유한요소법 · 강성행렬

8. Footnotes

1. 구장산술 제8장 “방정”의 계산법이 사실상 가우스 소거법과 동일하다. 편찬 시기는 한나라 때로, 그르차(Grcar)는 통용되는 주석본의 연도를 179년으로 잡는다. 어느 쪽이든 가우스보다 1600년 이상 앞선 셈이다. 수학사에서 이름은 마지막에 만진 사람이 다 가져가는 법이다. ↩

2. 뉴턴은 출판하기 싫어했고, 오일러는 권하지 않았고, 르장드르는 “보통의” 방법이라 불렀고, 가우스 본인도 “흔한(common)” 방법이라고 썼다. 그런데 이름은 가우스가 가져갔다. 그르차의 논문 제목이 「어쩌다 보통의 소거가 가우스 소거법이 되었나」인 이유다. ↩

3. 작은 피벗으로 나눈다는 건 큰 수를 만들었다가 다시 빼는 과정을 반복한다는 뜻인데, 부동소수점 세계에서 “비슷한 크기 큰 수끼리의 뺄셈”은 유효숫자가 증발하는 재앙적 소거(catastrophic cancellation)의 지름길이다. ↩

4. 굴드의 반례는 부동소수점 탐색으로 찾은 것이라 “정확 연산에서도 반례인가”가 한 해 동안 쟁점이었다. 참고로 완전 피벗팅 상한 자체는 윌킨슨의 고전적 결과가 반세기를 버텼고, 2020년대에 들어서야 개선된 상한이 나왔다. ↩