수치 선형대수

부동소수점: 컴퓨터의 실수는 격자다

컴퓨터의 실수는 연속된 직선이 아니라 유한한 점들의 격자입니다. IEEE 754 배정밀도(double)는 수를 $\pm(1.f_1f_2\cdots f_{52})_2 \times 2^{e}$로 저장합니다. 가수부 52비트, 지수 11비트, 부호 1비트, 모두 64비트입니다. 지수가 같은 구간 $[2^e, 2^{e+1})$ 안에서 점들은 같은 간격 $2^{e-52}$으로 놓이고, 구간이 두 배 커질 때마다 간격도 두 배가 됩니다. 그래서 절대 간격은 수의 크기에 따라 달라지지만 상대 간격은 거의 일정합니다.

연산 하나는 이렇게 정확한데 왜 계산이 틀릴까요? 범인은 대부분 상쇄 오차(catastrophic cancellation)입니다. 거의 같은 두 수 $a \approx b$를 빼면, 뺄셈 자체는 정확해도 a와 b에 이미 들어 있던 반올림 오차가 결과의 크기에 비해 엄청나게 커집니다. 예를 들어 $x = 10^{-8}$일 때 $\cos x = 1 - 5\times10^{-17}$인데 double은 이 값을 1과 구별하지 못하므로 $1 - \cos x$는 0이 됩니다. 정답은 $5\times10^{-17}$인데 말이죠.

해결책은 대수적으로 같지만 뺄셈이 없는 식으로 바꾸는 것입니다. $1 - \cos x = 2\sin^2(x/2)$, $\sqrt{1+x} - 1 = x/(\sqrt{1+x}+1)$, 그리고 $e^x - 1$은 라이브러리의 expm1이 정확히 계산합니다. 아래 시뮬레이터에서 x를 0으로 보내며 두 식을 비교해 보세요.

더 자세한 부동소수점 이야기(0.1의 2진 표현, 반올림 오차가 쌓여 생긴 실제 사고)는 MathBook "컴퓨터는 어떻게 틀리는가"에 있습니다. 이 장의 관심은 행렬입니다. 연산 하나의 오차가 u라면, n×n 행렬을 다루는 알고리즘 전체의 오차는 얼마일까요?

조건수: 문제 자체가 얼마나 예민한가

벡터의 크기는 1장의 노름으로 잽니다. 행렬의 크기는 "벡터를 최대 몇 배까지 늘이는가"로 재는 것이 자연스럽습니다.

등호가 성립하는 경우도 있습니다. $\mathbf{x}$가 A가 가장 많이 늘이는 방향이고 $\delta\mathbf{b}$가 가장 적게 늘이는 방향일 때입니다. 그러니 κ는 과장된 상한이 아니라 실제로 일어날 수 있는 최악의 증폭률입니다. 기하로 보면, 우변 오차가 반지름 r인 원 위에 있을 때 해의 오차는 $A^{-1}$이 그 원을 보낸 타원 위에 있고, 타원의 긴 축이 $r/\sigma_n$, 짧은 축이 $r/\sigma_1$입니다.

조건수가 큰 행렬의 고전적인 예는 힐베르트 행렬(Hilbert matrix) $H_{ij} = 1/(i+j-1)$입니다. 구간 [0, 1]에서 다항식 $1, x, x^2, \dots$로 함수를 최소제곱 근사할 때 나오는 행렬이고($H_{ij} = \int_0^1 x^{i-1}x^{j-1}\,dx$), 모든 원소가 깔끔한 분수인데도 n이 커질수록 조건수가 대략 $e^{3.5n}$꼴로 폭발합니다. 아래에서 참해가 $\mathbf{x} = (1, 1, \dots, 1)$이 되도록 $\mathbf{b} = H\mathbf{1}$을 만들고 부분 피벗팅 소거로 풀어 봅니다.

안정성: 알고리즘이 문제를 망치지 않는가

답 $\hat{\mathbf{x}}$의 오차를 재는 방법은 두 가지입니다. 전방 오차(forward error)는 답이 참값에서 얼마나 떨어졌나 $\|\hat{\mathbf{x}} - \mathbf{x}\|/\|\mathbf{x}\|$입니다. 후방 오차(backward error)는 거꾸로 묻습니다. "$\hat{\mathbf{x}}$가 정확한 답이 되려면 문제를 얼마나 바꿔야 하나?" 즉 $(A + \delta A)\hat{\mathbf{x}} = \mathbf{b}$를 만족하는 가장 작은 $\|\delta A\|/\|A\|$입니다.

두 정리를 합치면 실무의 어림셈이 나옵니다. 후방 안정한 알고리즘의 후방 오차는 $u \approx 10^{-16}$ 정도이므로 정리 16.2에 의해

거듭제곱법: 곱하기만 해서 고유벡터 찾기

아무 벡터 $\mathbf{x}_0$에 A를 계속 곱해 봅시다. A가 고유벡터 $\mathbf{v}_1,\dots,\mathbf{v}_n$으로 대각화되고 $|\lambda_1| \gt |\lambda_2| \ge \cdots$라 하면, $\mathbf{x}_0 = c_1\mathbf{v}_1 + \cdots + c_n\mathbf{v}_n$에서

수렴 속도가 $|\lambda_2/\lambda_1|$로 정해진다는 것이 약점입니다. 두 고유값이 가까우면 한없이 느립니다. 해결책은 A를 바꿔서 비율을 좋게 만드는 것입니다. $(A - \sigma I)^{-1}$의 고유값은 $1/(\lambda_i - \sigma)$이므로, 이동량 σ를 원하는 고유값 근처에 두면 그 고유값이 압도적으로 커집니다. 이것이 이동 역반복(shifted inverse iteration)이고, 매 단계 곱셈 대신 $(A - \sigma I)\mathbf{y} = \mathbf{x}$를 풉니다(LU는 한 번만 분해하면 됩니다). 수렴비는 $|\lambda_j - \sigma| / |\lambda_k - \sigma|$(가장 가까운 것 ÷ 두 번째로 가까운 것)입니다.

한 걸음 더 나가면, 매 단계 σ를 현재의 레일리 몫으로 바꿀 수 있습니다. 이것이 레일리 몫 반복(Rayleigh quotient iteration)이고, 대칭 행렬에서는 대개 3차 수렴합니다. 맞는 자릿수가 매 단계 세 배씩 늘어 서너 번이면 기계 정밀도에 닿습니다. 대신 매번 새 행렬을 분해해야 하고, 어느 고유값으로 갈지는 시작점에 달려 있습니다.

거듭제곱법은 단순하지만 아주 큰 행렬에서 여전히 쓰입니다. 곱셈만 하면 되므로 행렬을 저장하지 않고 "곱하는 방법"만 알아도 됩니다. 웹 페이지의 링크 구조에서 중요도를 매긴 초기 페이지랭크가 대표적인 예로, 전이 행렬의 가장 큰 고유값(=1)의 고유벡터를 거듭제곱법으로 구했습니다(12장).

QR 알고리즘: 모든 고유값을 한꺼번에

아이디어는 놀랄 만큼 짧습니다. A를 QR 분해하고, 두 인수를 거꾸로 곱합니다.

실제 라이브러리의 QR 알고리즘은 세 가지를 더합니다. ① 먼저 직교 닮음 변환으로 A를 헤센베르크 꼴(대각 바로 아래 한 줄 밑은 모두 0)로 만들어 한 단계의 비용을 $O(n^3)$에서 $O(n^2)$로 줄입니다. ② 역반복처럼 이동 $A_k - \sigma_k I = Q_kR_k,\; A_{k+1} = R_kQ_k + \sigma_kI$를 써서 수렴을 2차~3차로 빠르게 합니다. ③ 대각 아래 원소가 충분히 작아지면 그 고유값을 떼어 내고(감축, deflation) 더 작은 행렬로 계속합니다. 그 결과 n×n 행렬의 고유값 전부를 대략 $10n^3$ flop 정도에 구합니다. 복소 고유값 쌍은 실수 연산으로 2×2 블록을 남기는 "실수 슈르 꼴"로 처리합니다.

반복법: 야코비, 가우스–자이델, SOR

A를 쉽게 풀리는 부분 M과 나머지로 쪼갭니다: $A = M - N$. 그러면 $A\mathbf{x} = \mathbf{b}$는 $M\mathbf{x} = N\mathbf{x} + \mathbf{b}$이고, 이것을 반복식으로 바꿉니다.

오차 $\mathbf{e}_k = \mathbf{x}_k - \mathbf{x}$는 $\mathbf{e}_{k+1} = G\mathbf{e}_k$, $G = M^{-1}N$을 따릅니다. 그러므로 수렴 여부는 반복 행렬 G의 고유값 크기가 결정합니다. G의 고유값 절댓값의 최댓값을 스펙트럼 반지름(spectral radius) $\rho(G)$라 합니다.

위의 온도 문제(2차원 라플라스 방정식, 한 변 n칸, $h = 1/(n+1)$)에서는 스펙트럼 반지름을 정확히 계산할 수 있습니다. 야코비는 $\rho_J = \cos(\pi h)$, 가우스–자이델은 $\rho_{GS} = \rho_J^2$(두 배 빠름), 최적 $\omega^* = 2/(1+\sin\pi h)$를 쓴 SOR은 $\rho = \omega^* - 1 \approx 1 - 2\pi h$입니다. 격자가 촘촘해질수록 $\rho_J \approx 1 - \tfrac{\pi^2h^2}{2}$는 1에 붙어 버리므로 야코비는 $O(n^2)$단계, 최적 SOR은 $O(n)$단계가 필요합니다.

켤레 기울기법: 가장 빠른 내리막을 똑똑하게

f의 등고선은 A의 고유벡터를 축으로 하는 타원이고, 축 길이 비가 $\sqrt{\kappa}$입니다. 가장 단순한 전략은 매번 가장 가파른 방향 $-\nabla f = \mathbf{r}_k$(잔차!)로 그 직선 위의 최솟점까지 가는 최급강하법(steepest descent)입니다. 그런데 타원이 길쭉하면 골짜기 벽 사이를 지그재그로 오가며 매 단계 오차가 $(\kappa-1)/(\kappa+1)$배밖에 줄지 않습니다.

1952년 헤스테네스와 슈티펠이 발표한 켤레 기울기법(conjugate gradient, CG)은 방향을 더 똑똑하게 고릅니다. 새 방향을 이전 방향들과 A-켤레($\mathbf{p}_i^{\mathsf T}A\mathbf{p}_j = 0$)가 되도록 잡으면, 한 방향에서 최소화한 것을 다음 방향이 망치지 않습니다. 타원을 원으로 펴는 좌표에서 보면 서로 수직인 방향으로 한 번씩만 가는 것과 같습니다.

최급강하법의 $(\kappa-1)/(\kappa+1)$과 비교하면 CG는 κ 대신 $\sqrt{\kappa}$가 들어갑니다. κ = 10⁴이면 최급강하법은 오차를 10배 줄이는 데 약 1만 단계, CG는 약 100단계가 듭니다. 실제 계산에서는 반올림 때문에 정확히 n단계에 끝나지 않지만, 행렬이 클 때 CG는 n보다 훨씬 적은 단계에서 충분히 정확해지므로 "직접법"이 아니라 반복법으로 씁니다.

CG의 수렴은 κ에 달려 있으니, 같은 해를 갖되 조건수가 작은 문제로 바꿔 풀면 됩니다. 이것이 전처리(preconditioning)입니다. A와 비슷하면서 풀기 쉬운 SPD 행렬 M을 골라 $M^{-1}A\mathbf{x} = M^{-1}\mathbf{b}$를 풉니다(대칭성을 지키도록 변형한 "전처리 CG"는 매 단계 $M\mathbf{z} = \mathbf{r}$ 하나를 더 풉니다). M = A면 κ = 1이지만 M을 푸는 것이 원래 문제만큼 어렵고, M = I면 쉽지만 효과가 없습니다. 그 사이의 절충으로 대각만 쓰는 야코비 전처리, 채움을 버린 불완전 콜레스키, 격자 문제에서 거의 최적인 다중 격자 같은 방법이 쓰입니다. 큰 문제를 빨리 푸는 기술의 대부분은 좋은 전처리기를 설계하는 기술입니다.

희소 행렬과 계산 규모

대부분의 원소가 0인 행렬을 희소 행렬(sparse matrix)이라 하고, 0이 아닌 원소 수를 nnz라 씁니다. 가장 널리 쓰는 저장 형식은 CSR(compressed sparse row)입니다. 0이 아닌 값을 행 순서대로 늘어놓은 values, 각 값의 열 번호 col_idx, 그리고 각 행이 values의 어디서 시작하는지 적은 row_ptr(길이 n+1) 세 배열로 행렬을 표현합니다.

희소 행렬에도 직접법을 쓸 수 있습니다. 다만 LU 분해가 0이었던 칸을 채우는 채움 때문에 미지수의 번호를 잘 매기는 것(재순서화)이 성능을 좌우합니다. 2차원 격자 문제는 희소 직접법이 여전히 강하고, 3차원 문제나 아주 큰 문제는 채움이 감당할 수 없을 만큼 커져서 전처리한 반복법이 주로 쓰입니다. 대략적인 비교는 다음과 같습니다.

이 도구가 쓰이는 곳

수치 선형대수는 계산 과학과 공학 시뮬레이션의 밑바닥입니다. 이 책의 앞 장들이 "무엇을 계산하는가"였다면 이 장은 "어떻게 믿을 만하게, 충분히 빨리 계산하는가"입니다.