N1 · 부동소수점과 오차를 설명하는 언어#
1. 변동이 작은 관측값의 분산이 왜 0으로 나오는가#
같은 단위로 기록한 세 관측값이 \(x_1=10^8-1\), \(x_2=10^8\), \(x_3=10^8+1\)이라고 합시다. 예를 들어 큰 기준 잔액 주위에서 하루마다 조금씩 달라지는 계정 금액이며, 입력의 측정오차는 우선 제외합니다. 표본평균은 \(\bar x=10^8\), 불편 표본분산의 정의는
평균은 금액 단위, 분산은 그 단위의 제곱입니다. 실수의 대수법칙으로 분자를 전개하면 \(\sum x_i^2-3\bar x^2\)와 같습니다. 그러나 컴퓨터에서 제곱을 먼저 더하면 \(3\times10^{16}\) 크기의 두 수를 뺀 뒤 2만 남겨야 합니다. binary64에서는 이 크기에서 인접한 표현 가능 수의 간격이 이미 2보다 클 수 있습니다. 정확한 정보가 제곱·합 단계에서 없어지면 마지막 뺄셈으로 복원할 수 없습니다.
import numpy as np
from fractions import Fraction
x = np.array([1e8-1, 1e8, 1e8+1], dtype=np.float64)
centered = np.sum((x-x.mean())**2)/(len(x)-1)
raw = (np.sum(x*x)-len(x)*x.mean()**2)/(len(x)-1)
exact = sum((Fraction(int(v))-10**8)**2 for v in x)/2
print("정확한 값, 중심화, 원시 제곱합:", exact, centered, raw)
assert exact == 1 and centered == 1 and raw == 0
정확한 값, 중심화, 원시 제곱합: 1 1.0 0.0
이 예에서 모든 입력 정수는 정확히 저장됩니다. 따라서 문제는 입력을 십진수에서 이진수로 바꾸는 단계가 아니라 계산 과정에 있습니다. 입력오차와 계산오차를 분리해야 개선 방법을 정할 수 있습니다.
그림 88 관측값은 \(c-1,c,c+1\)이다. 가로축 \(c\)는 로그 눈금이며 세로축은 계산한 표본분산이다. 표시한 입력 범위에서는 세 관측 정수가 정확히 저장된다. 점들을 잇는 선은 서로 다른 계산 실험의 순서를 표시한다.#
2. 저장 가능한 수와 한 번의 반올림#
binary64의 정상 유한수는 부호를 제외하면 \((1.b_1b_2\cdots b_{52})_2\,2^e\), \(-1022\le e\le1023\) 형태입니다. 유효 이진 자릿수는 숨은 첫 자리까지 53개입니다. 정상수 아래에는 간격 \(2^{-1074}\)의 부분정상수가 있으며, 0도 표현됩니다. 범위를 넘는 결과에는 무한대 등이 쓰입니다.
1 바로 다음의 수는 \(1+2^{-52}\)입니다. 여기서 간격 \(\epsilon=2^{-52}\)와 최근접 반올림의 단위 반올림오차 \(u=2^{-53}\)를 구별합니다. 최근접에서 정확한 값이 중간이면 마지막 유효 비트가 짝수인 쪽을 고릅니다. 아래의 상대오차 모형은 결과가 정상수 범위에 있고 오버플로·유해한 언더플로가 없다는 가정 아래 사용합니다.
나눗셈에는 \(b\ne0\)도 필요합니다. 정확한 결과가 0이면 저장 결과도 0이며 \(\delta=0\)으로 둘 수 있습니다. 부분정상수에서는 상대오차 대신 절대 반올림오차 \(2^{-1075}\) 이하라는 표현을 써야 합니다. 이 장에서는 특별히 말하지 않으면 정상 범위, 최근접 반올림, 각 연산 뒤의 binary64 저장을 가정합니다. 표준 모형의 범위는 LAPACK의 부동소수점 설명과도 대조할 수 있습니다. 그 문헌의 \(\epsilon\)은 여기의 \(u\)를 뜻하므로 기호보다 정의를 확인해야 합니다. 연산을 합치는 FMA나 합산 순서를 바꾸는 병렬 라이브러리는 별도 알고리즘입니다.
\(0.1\)은 유한 이진분수로 표현되지 않습니다. 실제 저장값을 유리수로 읽으면
세 저장값으로 실행한 (0.1 + 0.2) - 0.3의 결과는 \(2^{-54}\)입니다. 이는 십진 실수 식이 0이라는 사실과 모순되지 않습니다. 컴퓨터는 반올림된 입력에 반올림되는 연산을 적용했습니다.
덧셈의 결합법칙도 더 이상 보장되지 않습니다. \(a=2^{53}\), \(b=1\), \(c=-2^{53}\)이면 \(a+b\)는 중간점 반올림으로 \(a\)가 되어 \((a+b)+c=0\)입니다. 반면 \(b+c=-(2^{53}-1)\)은 정확히 표현되어 \(a+(b+c)=1\)입니다. 같은 항을 어떤 순서로 더했는지까지 알고리즘의 일부입니다.
3. 큰 두 수를 빼는 것과 큰 오차를 만드는 것#
정확한 양수 \(a,b\)를 빼서 \(d=a-b\ne0\)를 구한다고 합시다. 입력에 \(|\Delta a|\le\eta|a|\), \(|\Delta b|\le\eta|b|\)가 이미 있으면
가까운 두 수의 차이는 입력의 상대오차에 민감한 문제입니다. 마지막 뺄셈 자체가 큰 반올림오차를 낸다는 뜻은 아닙니다. 예컨대 저장된 가까운 두 수의 뺄셈은 정확할 수도 있지만, 그 전에 각 수에 들어간 오차는 차이에 비해 커집니다. 원시 제곱합의 분산은 큰 중간 결과의 작은 상대오차를 최종 분산의 큰 상대오차로 바꾸었습니다.
중심화는 같은 분산을 계산하는 다른 식입니다. 자료를 한 번씩 읽어야 한다면 Welford 갱신을 쓸 수 있습니다. \(m_{k-1}\)이 앞 \(k-1\)개 평균이고 \(S_{k-1}\)이 평균 주위 제곱합일 때
로 갱신합니다. 정확한 산술에서 왜 같은 양인지 마지막 절에서 증명합니다. 이 식도 오차가 전혀 없는 것은 아니며 입력 범위·정밀도·순서의 영향을 받습니다. 큰 원시 제곱합 두 개를 직접 빼는 계산을 피한다는 개선입니다.
4. 여러 연산의 오차를 합하는 방법#
각 오차가 \(u\) 이하라도 천 번의 연산오차를 그냥 \(u\)라고 쓸 수는 없습니다. \(nu<1\)에서
를 정의합니다. \(n\)개의 \((1+\delta_i)\) 또는 그 역수를 곱한 결과는 \(1+\theta_n\), \(|\theta_n|\le\gamma_n\)로 표현할 수 있습니다. 작은 \(nu\)에서는 \(\gamma_n\approx nu\)이지만 증명에는 분모까지 있는 값을 씁니다.
순서대로 더한 \(\widehat s=\operatorname{fl}(x_1+\cdots+x_n)\)에 대해
상대오차는 여기에 \(\sum|x_i|/|\sum x_i|\)를 곱해야 합니다. 양·음이 상쇄되는 합에서는 이 비가 큽니다. 쌍을 묶어 균형 이진트리로 더하면 각 입력에서 최종 결과까지 연산 횟수는 \(\lceil\log_2n\rceil\) 이하이므로 상한의 \(\gamma_{n-1}\)을 \(\gamma_{\lceil\log_2n\rceil}\)로 바꿀 수 있습니다.
내적은 곱셈도 포함합니다. 곱을 따로 반올림하고 차례로 더하면
행렬곱에서는 각 내적에 적용하여 \(|\widehat C-AB|\le\gamma_n|A||B|\)를 얻습니다. 여기의 절댓값과 부등식은 성분별입니다. 이것만으로 모든 출력 성분에 공통인 작은 \(\Delta A\) 하나가 있어 \(\widehat C=(A+\Delta A)B\)라고 결론낼 수는 없습니다.
보상합산은 잃은 낮은 자리 정보를 다음 덧셈에 반영합니다. 아래 코드는 널리 쓰이는 Kahan 갱신입니다. 이 장에서 완전히 증명하는 것은 순차·쌍별 합의 위 상한이며, 보상합산의 일반 전진오차 정리는 별도 조망입니다. 특정 실험에서 개선된 결과와 모든 입력에 대한 보장을 구별합니다.
def sequential(values):
total = 0.
for value in values:
total = total + float(value)
return total
def pairwise(values):
values = list(values)
if len(values) <= 1:
return float(values[0]) if values else 0.
mid = len(values)//2
return pairwise(values[:mid]) + pairwise(values[mid:])
def kahan(values):
total = compensation = 0.
for value in values:
corrected = float(value) - compensation
updated = total + corrected
compensation = (updated-total) - corrected
total = updated
return total
values = [1.] + [2.**-53]*4096
exact = 1. + 2.**-41
answers = [sequential(values), pairwise(values), kahan(values)]
print("순차·쌍별·보상 절대오차:", [abs(v-exact) for v in answers])
assert answers[0] == 1. and answers[2] == exact
순차·쌍별·보상 절대오차: [4.547473508864641e-13, 0.0, 0.0]
그림 89 입력은 \(1\) 뒤에 \(2^{-53}\)을 \(n\)개 붙인 수열이다. 양축은 로그 눈금이다. 정확히 0인 오차는 표시 하한에 빈 표식으로 따로 나타내며 양의 오차로 해석하지 않는다.#
5. 잔차가 작다는 말에 무엇이 빠져 있는가#
두 부문의 균형을 \(Ax=b\)로 계산했다고 합시다. \(A\)는 정해진 계수, \(b\)는 관측한 순수요, \(x\)는 금액 단위의 상태입니다. 계산 결과 \(\widehat x\)를 대입한 잔차는 \(r=b-A\widehat x\)입니다. 정확한 해와의 전진오차는 \(\widehat x-x=-A^{-1}r\)이므로 잔차가 작아도 \(A^{-1}\)이 크면 해의 오차가 커질 수 있습니다.
여기서 \(r=(0,\varepsilon)^T\)이고 상대 전진오차는 \(\|\widehat x-x\|_2/\|x\|_2=1\)입니다. 잔차는 \(\varepsilon\to0\)에서 0으로 가지만 해의 상대오차는 그대로입니다. 두 식이 거의 같아지면서 두 성분을 따로 결정하기 어려워진 탓입니다.
유도노름에서 \(\kappa(A)=\|A\|\|A^{-1}\|\)라 쓰면
이것은 문제의 민감도와 계산 결과의 잔차를 결합한 상한입니다. \(A\)를 고정하고 \(b\)만 상대적으로 바꿀 때 해당 \(b\)에서의 정확한 국소 조건수는 \(\|A^{-1}\|\|b\|/\|x\|\)이며, 모든 \(b\) 중 최악이 \(\kappa(A)\)입니다. 조건수는 문제와 입력·출력 노름의 성질이고, 안정성은 알고리즘의 성질입니다.
그림 90 \(A_\varepsilon\)의 정확한 해는 \((1,1)\), 검사할 근삿값은 \((2,0)\)으로 고정했다. 곡선은 이 두 벡터의 정확한 식으로 계산했으며 특정 선형계 풀이기의 성능 실험이 아니다.#
6. 작은 입력 변화의 정확한 해로 설명하기#
후진오차는 \(\widehat x\)가 정확한 해가 되도록 입력을 최소한 얼마나 바꾸어야 하는지 묻습니다. 실수 또는 복소수의 2-노름을 사용하고 \(\widehat x\ne0\)라 합시다. \(b\)는 고정하고 \(A\)만 바꾸면
어떤 허용 변화도 \(\Delta A\widehat x=r\)를 만족하므로 이보다 작은 노름은 불가능합니다. 따라서 이 rank-one 변화는 최소값을 실제로 달성합니다. \(A\)와 \(b\)를 모두 상대적으로 같은 비율까지 바꿀 수 있게 하면 최소 후진오차는
대칭성·희소성·비음수 같은 추가 구조는 여기서 요구하지 않았습니다. 구조를 보존하는 후진오차는 다른 최적화 문제입니다. 잔차 자체를 부동소수점으로 계산하면 \(\widehat r\)에도 오차가 있어, 대략적인 크기는 \(\gamma_{n+1}(|b|+|A||\widehat x|)\)로 제한할 수 있습니다. 이 규모 이하의 계산 잔차를 정확한 0으로 해석해서는 안 됩니다.
문제 자체의 민감도는 단위를 바꾸어도 수치가 달라질 수 있습니다. 관측 시점 \(1999,2000,2001\)에 선형 추세 \(y=\alpha+\beta t\)를 맞추는 설계행렬 \(X=[\mathbf1,t]\) 대신 \(\tau=t-2000\)으로 써서 \(y=a+b\tau\)를 맞추면 \(a=\alpha+2000\beta\), \(b=\beta\)입니다. 적합 가능한 관측벡터의 공간은 같지만 열의 크기·각도와 계수의 해석은 달라집니다. 중심화한 행렬의 열은 직교하고 제곱길이는 3과 2입니다. 따라서 조건수는 \(\sqrt{3/2}\)입니다. 로그를 취하거나 차분하는 것은 이와 달리 관측모형 자체를 바꿀 수 있습니다.
7. 연습과 전체 풀이#
1. \(u\), 1 다음 수와의 간격, 가장 작은 양의 부분정상수를 구별하세요.
풀이. binary64 최근접에서 \(u=2^{-53}\)은 정상수 상대 반올림 상한입니다. 1 다음 수와의 간격은 \(2^{-52}=2u\)입니다. 가장 작은 양의 부분정상수는 \(2^{-1074}\)이며 상대오차 \(u\) 모형이 그 범위 전체에 적용되지는 않습니다.
2. \((2^{53}+1)-2^{53}\)을 정확한 산술과 binary64로 계산하세요.
풀이. 정확한 값은 1입니다. \(2^{53}\) 위의 표현 간격은 2이고 \(2^{53}+1\)은 중간점입니다. 짝수 유효 끝비트를 갖는 \(2^{53}\)으로 반올림되어 저장 계산은 0입니다. 괄호를 \(2^{53}+(1-2^{53})\)로 바꾸면 안쪽 정수는 정확하게 저장되어 1이 됩니다.
3. \(A=\operatorname{diag}(1,10^{-6})\), \(b=(1,0)^T\), \(\Delta b=(0,10^{-6})^T\)를 사용해 민감도를 확인하세요.
풀이. 원래 해는 \((1,0)^T\), 바뀐 해는 \((1,1)^T\)입니다. 상대 입력변화는 \(10^{-6}\), 상대 해 변화는 1입니다. \(\kappa_2(A)=10^6\)이므로 조건수 상한에 정확히 도달합니다.
4. 위 \(A_\varepsilon\), \(\widehat x=(2,0)^T\)에서 \(A\)만 바꾸는 최소 후진오차를 직접 구성하세요.
풀이. \(r=(0,\varepsilon)^T\)이고 \(\|\widehat x\|^2=4\)이므로 \(\Delta A=\begin{pmatrix}0&0\\\varepsilon/2&0\end{pmatrix}\)입니다. 둘째 행에 \(\widehat x\)를 곱하면 \(2+\varepsilon\)가 되어 \(b\)와 일치합니다. \(\|\Delta A\|_2=\varepsilon/2=\|r\|/\|\widehat x\|\)여서 최소입니다.
5. \(X=[\mathbf1,t]\)와 중심화한 \(Z=[\mathbf1,t-2000]\)의 관계를 쓰세요.
풀이. \(X=Z\begin{pmatrix}1&2000\\0&1\end{pmatrix}\)입니다. 오른쪽 행렬은 가역이므로 열공간은 같습니다. \(X(\alpha,\beta)^T=Z(\alpha+2000\beta,\beta)^T\)여서 적합벡터는 같고 계수 좌표만 다릅니다. 작은 조건수는 새 좌표의 상대오차를 말하며 원래 절편으로 돌아가는 변환의 민감도까지 없애지는 않습니다.
6. \(x=(1,-1+\epsilon)\)의 합에서 입력 상대오차의 조건수를 구하세요. \(0<\epsilon<1\)입니다.
풀이. 정확한 합은 \(\epsilon\)이고 절댓값 합은 \(2-\epsilon\)입니다. 성분별 상대오차를 각각 \(\eta\) 이하로 허용하면 합의 최악 상대오차는 \(\eta(2-\epsilon)/\epsilon\)입니다. 첫 성분 변화는 \(+\eta\), 둘째 변화도 \(+\eta(1-\epsilon)\)로 택하면 등호를 달성합니다. 이 문제의 민감도는 합산 알고리즘을 바꾼다고 사라지지 않습니다.
8. 지금까지의 내용을 수학의 언어로 정리해 봅시다#
이제 반올림의 국소 모형에서 합산오차를 유도하고, 잔차와 해의 오차를 연결합니다. 모든 상대오차에는 분모가 0이 아니라는 조건을, 모든 \(\gamma_k\)에는 \(ku<1\)을 붙입니다.
보조정리 1. 반올림과 곱의 누적#
정상 범위의 최근접 이진 \(p\)자리 반올림에서 \(\operatorname{fl}(z)=z(1+\delta)\), \(|\delta|\le2^{-p}\)입니다. 또한 \(|\delta_i|\le u\), \(e_i\in\{-1,1\}\)이면 \(\prod_{i=1}^n(1+\delta_i)^{e_i}=1+\theta_n\), \(|\theta_n|\le\gamma_n\)입니다.
증명. \(2^e\le|z|<2^{e+1}\)인 구간의 간격은 \(2^{e-p+1}\)입니다. 최근접 오차는 그 절반 \(2^{e-p}\) 이하이며 \(|z|\ge2^e\)로 나누면 \(2^{-p}\)입니다. 경계에서는 작은 쪽 간격이 더 작아 같은 상한을 만족합니다.
곱의 각 인자는 \(1-u\) 이상, \((1-u)^{-1}\) 이하입니다. \(e_i=-1\)일 때도 \((1+u)^{-1}\ge1-u\)이고 \(e_i=1\)일 때도 \(1+u\le(1-u)^{-1}\)이기 때문입니다. Bernoulli 부등식 \((1-u)^n\ge1-nu\)는 귀납으로 얻습니다. 따라서 전체 곱은 \(1-nu\) 이상, \((1-nu)^{-1}\) 이하입니다. 1을 빼면 아래쪽 오차는 \(nu\le\gamma_n\), 위쪽은 \(nu/(1-nu)=\gamma_n\)입니다. ∎
정리 2. 순차합·쌍별합·내적의 오차#
4절의 세 상한이 성립합니다.
증명. 순차합에서 첫 덧셈의 결과는 \((x_1+x_2)(1+\delta_2)\), 다음은 \(((x_1+x_2)(1+\delta_2)+x_3)(1+\delta_3)\)입니다. 계속 분배하면 각 \(x_i\)에 붙는 인자는 자기 투입 시점부터 마지막까지의 반올림 인자의 곱입니다. 개수는 최대 \(n-1\)이므로 결과는 \(\sum_i x_i(1+\theta_i)\), \(|\theta_i|\le\gamma_{n-1}\)입니다. 정확한 합을 빼고 삼각부등식을 적용합니다.
균형 쌍별합도 같은 전개를 트리의 잎에서 뿌리까지 합니다. 잎마다 인자의 수가 최대 \(\lceil\log_2n\rceil\)이므로 같은 증명에서 그 값으로 바뀝니다. 내적에서는 각 입력항 \(x_iy_i\)에 곱셈의 반올림 인자가 하나 더 붙어 최대 \(n\)개가 됩니다. 성분별 행렬곱 상한은 행과 열의 각 내적에 이 결과를 적용한 것입니다. ∎
정리 3. Welford 갱신의 정확한 항등식#
\(m_k=k^{-1}\sum_{i\le k}x_i\), \(S_k=\sum_{i\le k}(x_i-m_k)^2\)이면 3절의 갱신이 정확한 산술에서 성립합니다.
증명. 평균식에서 \(m_k-m_{k-1}=(x_k-m_{k-1})/k=d/k\)입니다. 앞 \(k-1\)개 항을 새 평균으로 바꾸면
교차항은 \(\sum_{i<k}(x_i-m_{k-1})=0\)이어서 사라집니다. 새 항은 \((x_k-m_k)^2=d^2(k-1)^2/k^2\)입니다. 두 증가량의 합은 \(d^2(k-1)/k=d(x_k-m_k)\)이므로 결론입니다. \(k=1\)에서는 \(m_1=x_1\), \(S_1=0\)으로 시작합니다. ∎
정의와 정리 4. 조건수와 선형계의 전진오차#
미분 가능한 \(f\)의 상대 노름 조건수는 \(x\ne0\), \(f(x)\ne0\)에서
고정된 가역 \(A\)의 \(f(b)=A^{-1}b\)에서는 앞서 제시한 조건수를 얻고, \(b\) 전체에서의 최댓값은 \(\kappa(A)\)입니다. 또한 \((A+E)(x+\Delta x)=b+f\), \(Ax=b\), \(\|E\|\le\eta\|A\|\), \(\|f\|\le\eta\|b\|\), \(\eta\kappa(A)<1\)이면
증명. 미분의 정의에서 나머지는 \(o(\|h\|)\)이며 작은 공 전체에서 그 비가 0으로 갑니다. 상한은 유도노름으로 얻고, 유한차원 단위구면에서 \(Df(x)\)의 노름을 달성하는 방향으로 \(h\)를 택하면 같은 하한을 얻습니다. 선형 \(f\)의 미분은 \(A^{-1}\)입니다. \(\|b\|/\|A^{-1}b\|\le\|A\|\)이며 \(\|Av\|=\|A\|\|v\|\)인 \(v\)에 \(b=Av\)를 넣으면 등호이므로 최악은 \(\kappa(A)\)입니다.
섭동식에서 \((I+A^{-1}E)\Delta x=A^{-1}f-A^{-1}Ex\)입니다. \(\|A^{-1}E\|\le\eta\kappa(A)<1\)이므로 H8의 Neumann 급수로 역행렬 노름은 \(1/(1-\eta\kappa(A))\) 이하입니다. 따라서 분자는 \(\|A^{-1}\|(\eta\|b\|+\eta\|A\|\|x\|)\le2\eta\kappa(A)\|x\|\)이고 결론을 얻습니다. 잔차 상한은 \(\widehat x-x=-A^{-1}r\)와 \(\|b\|\le\|A\|\|x\|\)를 결합한 특수한 경우입니다. ∎
정리 5. 2-노름 후진오차와 특이행렬까지의 거리#
\(\widehat x\ne0\)이고 \(a=\|A\|_2>0\), \(c=\|b\|_2\)라 합시다. \(\|E\|_2\le\eta a\), \(\|f\|_2\le\eta c\), \((A+E)\widehat x=b+f\)가 가능한 최소 \(\eta\)는 6절의 식입니다. 가역행렬에서 \(\min_{A+E\text{ 특이}}\|E\|_2=\sigma_{\min}(A)\)이며 따라서 상대거리는 \(1/\kappa_2(A)\)입니다.
증명. \(r=E\widehat x-f\)에서 \(\|r\|\le\eta(a\|\widehat x\|+c)\)이므로 제시한 값이 하한입니다. \(s=\|\widehat x\|\), \(d=as+c\)라 하고
를 택합니다. \(E\widehat x-f=r\)이고 \(\|E\|=a\|r\|/d\), \(\|f\|=c\|r\|/d\)이므로 하한을 달성합니다. \(r=0\)일 때에도 영행렬·영벡터로 성립합니다.
\(A+E\)가 특이면 핵에 속한 단위벡터 \(v\)가 있어 \(Av=-Ev\)입니다. 그러면 \(\|E\|\ge\|Av\|\ge\sigma_{\min}(A)\)입니다. 최소 특이벡터 쌍 \(Av=\sigma_{\min}u\)에 \(E=-\sigma_{\min}uv^*\)를 택하면 \((A+E)v=0\)이고 \(\|E\|=\sigma_{\min}\)입니다. H4의 SVD에서 \(\|A\|_2=\sigma_{\max}\), \(\|A^{-1}\|_2=1/\sigma_{\min}\)이므로 상대거리 식을 얻습니다. 같은 rank-one 변화는 Frobenius 노름도 \(\sigma_{\min}\)여서 그 절대거리 역시 같습니다. ∎
입력의 민감도와 알고리즘의 반올림오차를 구별했습니다. 다음 장에서는 소거법의 각 연산을 이 관점으로 다시 읽고, 피벗과 행렬 구조가 정확도와 비용에 미치는 영향을 봅니다. N2로 이어 읽기.