조건수와 안정성

문제의 민감도와 알고리즘의 품질

(가)에서는 계산이 완벽해도 소용이 없습니다. 입력 $x$ 자체가 $10^{-8}$만큼 불확실하고, $\ln$이 1 근처에서 입력의 상대 변화를 1만 배로 키우기 때문입니다. 이것은 문제의 성질입니다. (나)에서는 문제 $x \mapsto \ln(1 + x)$가 전혀 예민하지 않습니다. 입력이 $10^{-8}$ 바뀌면 답도 약 $10^{-8}$ 바뀝니다. 오차는 계산 순서 탓입니다. $1 + x$를 먼저 반올림하는 순간 $x$의 아래 자리가 $u \approx 10^{-16}$만큼 잘려 나가는데, 그 잘림은 $x$에 견주면 $10^{-16}/10^{-10} = 10^{-6}$ 정도입니다. 이것은 알고리즘의 성질입니다.

이 구분은 실무에서 질문을 바꿔 줍니다. 결과가 이상하면 먼저 "조건수가 얼마인가"를 묻습니다. 조건수가 $10^{12}$이면 배정밀도의 16자리 중 12자리는 원래 믿을 수 없는 것이니, 알고리즘을 고칠 것이 아니라 문제를 바꿔야 합니다(측정을 더 정확히, 모델을 다르게, 변수를 바꿔서). 조건수가 작은데 결과가 나쁘면 그때는 알고리즘을 바꿉니다. MathBook의 조건수 맛보기와 LinearBook의 행렬 조건수를 여기서 일반적인 틀로 다시 세웁니다.

상대 조건수: 입력 오차가 몇 배로 퍼지는가

입력이 $x$에서 $x + \Delta x$로 바뀌면 출력은 $f(x + \Delta x) - f(x) \approx f'(x)\,\Delta x$만큼 바뀝니다(1차 테일러 근사, CalculusBook의 미분과 오차 전파). 양변을 $f(x)$로 나누고 오른쪽에 $x/x$를 끼워 넣으면 상대 변화끼리의 관계가 나옵니다.

표에서 몇 가지 규칙이 보입니다. 곱셈형 연산(거듭제곱, 역수)은 언제나 조건이 좋습니다. 조건이 나빠지는 곳은 거의 언제나 출력이 0에 가까워지는 곳($\ln 1 = 0$, $\sin\pi = 0$, $\cos(\pi/2) = 0$)이거나 출력이 무한대로 가는 곳($\tan$의 극점)입니다. 0 근처에서 상대 오차는 작은 절대 변화에도 크게 바뀌기 때문입니다. 지수 함수는 특별합니다. $e^x$의 조건수 $|x|$는 $x = 700$이면 700이라, 입력의 $u$가 출력에서는 $700u \approx 8\times10^{-14}$가 됩니다. 공학의 아레니우스 식 $e^{-E/RT}$처럼 지수가 큰 모델이 온도 측정 오차에 예민한 이유입니다.

함수별 조건수: 빼기, 여러 변수, 방정식의 근

입력이 여러 개인 문제도 같은 생각으로 다룹니다. 각 입력이 상대적으로 $\delta_i$만큼 흔들릴 때 출력의 상대 변화는 1차까지 $\sum_i (x_i\,\partial f/\partial x_i)/f \cdot \delta_i$입니다.

뺄셈의 조건수 $(|x| + |y|)/|x - y|$는 1장 정리 1.9의 증폭 인자 그대로입니다. 1장에서는 "뺄셈이 오차를 키운다"고 했지만, 이제 더 정확히 말할 수 있습니다. 뺄셈 연산은 정확하거나 거의 정확합니다. 키우는 것은 $x \approx y$에서 조건이 나쁜 뺄셈이라는 문제이고, 그 문제에 들어가는 입력이 이미 오차를 품고 있을 때 피해가 드러납니다. 예를 들어 $1.0001 - 1$의 조건수는 $2.0001/0.0001 \approx 2.0\times10^4$입니다.

아래 지도는 입력이 둘인 문제의 조건수를 평면 전체에 색으로 칠합니다. 조건이 나쁜 곳은 대개 곡선 하나를 따라 몰려 있습니다. 뺄셈은 대각선 $x = y$, 이차방정식의 근은 판별식이 0인 포물선 위입니다.

방정식의 근은 얼마나 예민한가

3장에서 푸는 문제 $f(x) = 0$도 하나의 문제입니다. 입력은 함수 $f$(실제로는 그 계수나 매개변수), 출력은 근 $r$입니다. $f$를 $f + \varepsilon g$로 조금 바꾸면 근이 $r + \Delta r$로 움직이고, $0 = f(r + \Delta r) + \varepsilon g(r + \Delta r) \approx f'(r)\Delta r + \varepsilon g(r)$에서 $\Delta r$을 구할 수 있습니다.

분모 $f'(r)$이 핵심입니다. 근에서 그래프가 납작하면, 곧 $f'(r)$이 작으면 $f$를 위아래로 조금만 흔들어도 그래프가 0을 지나는 자리가 크게 미끄러집니다. 3장 멈춤 조건에서 본 "불확실 구간"이 바로 이것입니다. $f'(r) = 0$인 중근에서는 정리 2.5가 아예 쓸모없고, 더 나쁜 일이 일어납니다.

윌킨슨 다항식: 서로 멀리 떨어진 단순근도 예민할 수 있다

결과는 충격적입니다. 아래 표는 이 장의 시뮬레이터와 같은 방법(곱 꼴로 정확히 계산한 $w$에 교란 항을 더하고 에를리히–애버스 반복으로 20개 근을 동시에 구함)으로 얻은 값입니다. 근 10부터 19까지가 복소수 쌍 다섯 개로 바뀌고, 그중 하나는 실수축에서 2.81이나 떨어집니다. 이 값은 윌킨슨이 발표한 $16.730737 \pm 2.812625i$와 일치합니다.

정리 2.5가 이유를 말해 줍니다. 근 16에서 $\partial r/\partial a_{19} = -r^{19}/w'(r) = -16^{19}/(15!\cdot4!) \approx -2.4\times10^9$입니다. 분자 $16^{19} \approx 7.6\times10^{22}$이 분모 $|w'(16)| = 15!\,4! \approx 3.1\times10^{13}$보다 압도적으로 큽니다. 상대 조건수 $3.2\times10^{10}$에 교란 $5.7\times10^{-10}$을 곱하면 18로, 1차 근사가 의미 없을 만큼 큰 변화입니다. 근이 서로 충돌해 복소수로 튀어나간 것입니다.

더 놀라운 사실이 있습니다. 계수 21개 가운데 $x^3$부터 $x^7$까지의 계수 다섯 개는 $2^{53}$보다 커서 배정밀도에 정확히 담기지 않습니다. 계수를 double로 저장만 해도(상대 교란 $10^{-16}$ 이하) 근 13의 오차가 $6.2\times10^{-4}$가 됩니다. "모든 계수" 열의 조건수 $5.4\times10^{13}$과 $u$를 곱한 $6\times10^{-3}$이 그 상한입니다. 그래서 고유값을 특성 다항식의 근으로 구하지 않습니다(14장). 고유값 문제 자체는 조건이 좋을 수 있는데, 다항식 계수로 옮기는 순간 훨씬 예민한 문제로 바꾸게 됩니다. 반대로 동반 행렬의 고유값으로 다항식의 근을 구하는 방법(4장)은 계수에 대해 후방 안정해서, 오차가 이 조건수가 허락하는 만큼만 생깁니다.

전방 오차와 후방 오차, 후방 안정

가장 단순한 예는 사칙연산 자체입니다. 표준 모형 $\mathrm{fl}(x - y) = (x - y)(1 + \delta)$는 $= x(1 + \delta) - y(1 + \delta)$로 읽을 수 있습니다. 입력을 각각 상대적으로 $u$ 이하만큼 바꾼 문제의 정확한 답이므로 뺄셈은 후방 안정합니다. 그런데도 $x \approx y$이면 전방 오차가 큽니다. 문제(뺄셈)의 조건수가 크기 때문입니다. "뺄셈은 불안정하다"는 말은 틀린 말이고, "$x \approx y$인 뺄셈은 조건이 나쁘다"가 맞는 말입니다.

$x = 2.01$에서 계산값은 부호조차 틀렸지만 상한 안에 있습니다. 호너 방법은 할 일을 다 했고, 전개된 계수로 근 근처의 값을 묻는 문제가 $10^{23}$만큼 조건이 나쁜 것입니다. 3장 예제 3.7의 $(x - 1)^7$ 잡음과 같은 현상입니다. 처방은 알고리즘이 아니라 표현을 바꾸는 것입니다. 곱 꼴 $(x - 2)^9$이나 근 근처의 테일러 계수로 저장하면 조건수가 9가 됩니다.

아래 시뮬레이터는 여러 계산에서 전방 오차를 재고, 정리 2.9를 거꾸로 써서 후방 오차를 $\eta \approx (\text{전방 오차})/\kappa$로 추정합니다. 안정한 알고리즘은 $\eta$가 $u$ 근처에 붙어 있고, 불안정한 알고리즘은 $\eta$가 $u$보다 훨씬 위로 올라갑니다.

안정한 알고리즘과 불안정한 알고리즘

조건이 좋은 문제를 불안정하게 계산하는 알고리즘은 대개 조건이 나쁜 중간 단계를 거칩니다. $\ln(1 + x)$를 "$w = 1 + x$, 그다음 $\ln w$"로 계산하면 둘째 단계의 조건수가 $1/|\ln w| \approx 1/|x|$입니다. 첫째 단계에서 생긴 $u$의 오차가 둘째 단계에서 $1/|x|$배로 커집니다. 문제 전체의 조건수는 1인데 알고리즘이 스스로 예민한 하위 문제를 만든 것입니다.

점화식: 같은 식, 반대 방향, 정반대의 결과

원인을 추적해 봅시다. 계산값 $\hat I_{n-1} = I_{n-1} + e_{n-1}$을 점화식에 넣으면 $\hat I_n = 1 - n\hat I_{n-1} = I_n - n e_{n-1}$이므로 오차는 $e_n = -n e_{n-1}$, 곧

처방은 점화식을 거꾸로 쓰는 것입니다. $I_{n-1} = (1 - I_n)/n$이면 오차는 $e_{n-1} = -e_n/n$, 단계마다 $n$으로 나뉩니다. 시작값을 아무렇게나, 예를 들어 $I_N = 0$으로 놓아도 오차 $|e_N| = I_N \lt 1/(N + 1)$이 $N!/n!$로 나뉘어 곧 사라집니다. 베셀 함수처럼 점화식으로 계산하는 특수 함수 라이브러리가 쓰는 밀러 알고리즘(Miller's algorithm)이 이 생각입니다.

뒤로 가는 계산은 시작값이 완전히 틀렸는데도($I_{20} = 0$, 상대 오차 100%) 열 단계 만에 12자리가 맞습니다. 앞으로 가는 계산은 시작값이 16자리 정확한데도 18단계 만에 부호가 틀립니다. 같은 수학 공식이 계산 방향에 따라 안정하기도, 불안정하기도 합니다.

급수: x < 0에서 eˣ의 테일러 급수

항의 크기를 보면 답이 나옵니다. $x = -20$이면 $|x|^k/k!$는 $k = 19, 20$에서 최대 $4.31\times10^7$이 되고, 부호가 번갈아 바뀌는 이 큰 항들이 서로 지워져 $2\times10^{-9}$이 남아야 합니다. 그런데 $4.3\times10^7$ 크기의 수를 더할 때마다 반올림 오차가 $u\times4.3\times10^7 \approx 5\times10^{-9}$씩 생깁니다. 답보다 큰 오차입니다. 1장 정리 1.10의 합산 오차 상한 $\propto u\sum|t_k|$가 정확히 이 상황입니다.

고치는 법은 간단합니다. $x \lt 0$이면 $e^{x} = 1/e^{|x|}$로 계산합니다. $e^{|x|}$의 급수는 항이 모두 양수라 상쇄가 없고(합의 조건수 1), 역수는 조건수 1입니다.

다시 쓰기로 정확도 구하기: log1p, expm1, 근의 공식

정리 2.10으로 진단합니다. 첫 단계 $w = 1 + 0.05/n$에서 $0.05/n$의 아래 자리가 잘려 나가고($10^{-15}$이면 $u$와 비슷한 크기), 다음 단계 $w^n$은 $w$에 대해 조건수가 $n$입니다. 처방은 $w^n = \exp\bigl(n\ln(1 + r/n)\bigr)$으로 쓰고, $\ln(1 + z)$를 $1 + z$를 만들지 않고 계산하는 함수 log1p(z)를 쓰는 것입니다.

(표의 상대 오차는 같은 $n$에서의 정확한 $(1 + 0.05/n)^n$에 대한 것이다.) 수치 라이브러리가 log1p와 expm1($e^x - 1$)을 따로 제공하는 이유입니다. 두 함수는 작은 $x$에서 조건이 좋은 문제($\kappa \approx 1$)인데, 이름 그대로의 식으로 계산하면 상쇄가 끼어드는 경우입니다.

라이브러리가 없다면 어떻게 할까요? 카한이 알려 준 놀라운 요령이 있습니다. $w = \mathrm{fl}(1 + x)$는 틀렸지만, 틀린 $w$에 대해 $\ln w$와 $w - 1$은 모두 정확히 계산됩니다($w - 1$은 스테르베즈 보조정리로 정확). 그러면 매끄러운 함수 $\ln w/(w - 1)$은 $w$ 대신 $1 + x$에서의 값과 거의 같으므로, 여기에 $x$를 곱하면 오차가 상쇄됩니다.

// 카한의 log1p: ln(1 + x)를 상쇄 없이 (작은 |x|)
w ← fl(1 + x)
if w = 1 then return x          // |x| < u: ln(1 + x) ≈ x
return ln(w) · x / (w − 1)        // ln(w)/(w − 1)은 w의 오차에 둔감

// 같은 생각의 expm1: e^x − 1
v ← exp(x)
if v = 1 then return x
if v − 1 = −1 then return −1     // x가 매우 작은 음수
return (v − 1) · x / ln(v)

앞의 시뮬레이터에서 '다시 쓴 식'이 이 두 공식입니다. 이 장과 1장에서 본 다시 쓰기를 모으면 몇 가지 원칙이 보입니다.

행렬 조건수 맛보기

스칼라 조건수를 벡터로 옮기려면 크기를 재는 도구, 곧 노름이 필요합니다. 벡터 노름 $\|\mathbf x\|$과 그에 맞는 행렬 노름 $\|A\| = \max_{\mathbf x \ne 0}\|A\mathbf x\|/\|\mathbf x\|$(12장에서 자세히)을 쓰면 다음이 성립합니다.

위 문제에서 $\|\delta\mathbf b\|_2/\|\mathbf b\|_2 = 10^{-4}/2.83 = 3.5\times10^{-5}$, $\|\delta\mathbf x\|_2/\|\mathbf x\|_2 = 1$이므로 증폭은 $2.8\times10^4$배로, 상한 $\kappa_2 = 4.0\times10^4$ 안에 있습니다. 조건이 나쁜 행렬의 대표는 힐베르트 행렬 $H_{ij} = 1/(i + j - 1)$입니다. 다항식 최소제곱 맞춤에서 자연스럽게 나오는 행렬입니다(7장). 아래 표는 $H\mathbf x = H\mathbf 1$을 부분 피벗팅 가우스 소거로 풀어 참해 $\mathbf 1$과 비교한 결과입니다($\kappa_2$는 $H$와 정확한 정수 행렬 $H^{-1}$의 최대 고유값의 곱).

이 표는 이 장 전체의 요약입니다. 잔차는 $n$과 상관없이 $10^{-16}$ 수준입니다. 부분 피벗팅 가우스 소거가 (실제로는) 후방 안정하다는 뜻입니다. 해의 오차는 $\kappa_2 u$를 거의 그대로 따라갑니다(정리 2.9). $n = 12$이면 배정밀도로 한 자리도 믿을 수 없는데, 알고리즘 탓이 아니라 문제 탓입니다. 그래서 잔차가 작다고 해가 정확한 것은 아닙니다. 반복 해법의 멈춤 조건(13장), 반데르몽드 행렬과 보간(5장), 정규 방정식의 조건수 제곱 문제(7장)가 모두 이 사실에서 출발합니다. 노름의 정의, 오차 한계의 증명, 반복 개선은 12장에서 다룹니다.

반올림 오차 모형과 혼합 정밀도

지금까지 "알고리즘의 후방 오차가 $u$ 정도"라고 말했습니다. 이것을 실제로 증명하는 도구가 1장의 표준 반올림 오차 모형입니다. 부동소수점 수 $x, y$와 $\circ \in \{+, -, \times, \div\}$에 대해

이 결과들은 정밀도를 바꿔 쓸 때 진가를 발휘합니다. 그래픽 처리 장치와 기계 학습 가속기는 16비트 형식을 즐겨 씁니다. 반정밀도(binary16)의 $u = 2^{-11} \approx 4.9\times10^{-4}$, bfloat16의 $u = 2^{-8} \approx 3.9\times10^{-3}$입니다(1장 표). 정리 2.14의 상한 $\gamma_n$은 반정밀도에서 $n = 1000$이면 $0.95$, $n \ge 2048$이면 아예 무의미합니다. 실제로도 그렇습니다. 덧셈 결과가 2048을 넘으면 반정밀도의 간격이 2라 1보다 작은 항은 더해지지 않습니다.

혼합 정밀도(mixed precision)는 이 상황의 표준 처방입니다. 데이터는 낮은 정밀도로 저장하고(메모리와 대역폭이 절반 이하), 누적만 높은 정밀도로 합니다. 16비트 두 수의 곱은 가수가 22비트 이하라 단정밀도에 정확히 담기므로, 곱은 정확하고 덧셈만 단정밀도 $u = 6\times10^{-8}$으로 반올림됩니다. 오차 상한의 $u$가 누산기의 $u$로 바뀌는 것입니다. 행렬 곱 전용 하드웨어(텐서 코어 등)가 "16비트 입력, 32비트 누적"으로 설계된 이유입니다.

bfloat16 누적은 합이 256에 이르자 멈췄습니다. 256 근처에서 bfloat16의 간격은 2이고, 1보다 작은 항은 반 간격보다 작아 반올림으로 사라집니다. 1장의 조화급수가 멈춘 것과 같은 현상입니다. 누산기를 32비트로 바꾸면 같은 입력으로 오차가 백만 배 줄어듭니다. 이 생각을 선형계로 넓힌 것이 혼합 정밀도 반복 개선입니다. 낮은 정밀도로 LU 분해를 하고, 잔차 $\mathbf r = \mathbf b - A\hat{\mathbf x}$만 높은 정밀도로 계산해 보정 $A\,\delta\mathbf x = \mathbf r$을 되풀이하면, $\kappa(A)$와 낮은 정밀도의 $u$의 곱이 1보다 충분히 작은 한 높은 정밀도 수준의 해를 얻습니다(12장).

이 도구가 쓰이는 곳

"문제의 조건수 × 알고리즘의 후방 오차"라는 틀은 이 책의 모든 장에서 오차를 읽는 기본 문법입니다.