부동소수점은 어디서 새는가: 반올림, 상쇄, 더하는 순서

같은 스크립트를 두 번 돌렸는데 학습 손실이 소수 넷째 자리에서 갈립니다. 시드도 같고 데이터도 같고 순서도 같은데 GPU 개수만 둘에서 넷으로 바꿨습니다. 「부동소수점 오차겠지」로 넘기고 나면, 몇 주 뒤 fp16으로 추론을 내리다가 어텐션 출력이 통째로 nan이 되는 것을 보고 다시 같은 자리로 돌아옵니다.

두 현상은 같은 뿌리에서 나옵니다. 컴퓨터가 담는 실수는 실수가 아니라 유한한 개수의 눈금이고, 그 눈금이 어디에 얼마 간격으로 서 있는지를 알면 두 가지 모두 예측 가능한 사건이 됩니다.

초급 47번 「자릿수와 오차」에서 유효숫자와 상대오차를 손계산으로 다뤘는데, 밑이 2인 진법을 세우지 않아 「0.1 + 0.2가 왜 0.3이 아닌가」와 머신 엡실론을 이 글로 넘겼습니다. 여기서 받습니다.

0.1은 애초에 0.1이 아니다

>>> f"{0.1:.20f}"
'0.10000000000000000555'
>>> f"{0.1 + 0.2:.20f}"
'0.30000000000000004441'
>>> f"{0.3:.20f}"
'0.29999999999999998890'

덧셈이 틀린 것이 아닙니다. 입력이 이미 0.1과 0.2가 아니었습니다.

이유는 진법입니다. 십진법에서 $1/3$ 이 $0.333\ldots$ 로 안 끝나듯, 이진법에서는 $1/10$ 이 안 끝납니다. 분모가 밑의 소인수만으로 이뤄져야 유한 소수가 되는데, 10의 소인수는 2와 5이고 이진법의 밑은 2뿐이라 5를 처리할 방법이 없습니다. 그래서

$$0.1 = 0.0001100110011\overline{0011}\ldots_{(2)}$$

이고, 유한한 비트에 담으려면 어딘가에서 잘라야 합니다. 파이썬이 실제로 담고 있는 값은 정확히

$$\frac{3602879701896397}{36028797018963968} = 0.1000000000000000055511151231257827\ldots$$

입니다. 0.2도 마찬가지로 조금 크고, 둘을 더한 값이 0.3에 가장 가까운 눈금보다 한 칸 위에 떨어집니다.

여기서 나오는 첫 번째 규칙입니다. 부동소수점 값을 == 로 비교하지 않습니다. 두 계산 경로가 수학적으로 같아도 눈금 하나만큼 갈릴 수 있습니다.

부호 · 지수 · 가수, 그리고 머신 엡실론

IEEE 754는 실수 하나를 세 조각으로 담습니다.

$$x = (-1)^{s} \times 1.f \times 2^{e}$$

$s$ 는 부호 1비트, $e$ 는 지수, $1.f$ 는 가수(mantissa)입니다. 맨 앞의 1은 항상 1이라 저장하지 않고, 소수점 아래 $p$ 비트만 담습니다. 이 $p$ 가 정밀도를 정합니다.

지수 비트는 범위를, 가수 비트는 정밀도를 정합니다. 지수 비트를 하나 늘리면 표현 가능한 범위가 제곱으로 넓어지고, 가수 비트를 하나 늘리면 눈금 간격이 절반이 됩니다.

가수가 $p$ 비트라는 말은 $[1, 2)$ 구간을 $2^p$ 등분한다는 뜻입니다. 그 구간의 눈금 간격이 $2^{-p}$ 이고, 이 값을 머신 엡실론이라고 부릅니다 — 1에 더했을 때 결과가 1과 달라지는 가장 작은 수입니다. fp32에서 확인하면 이렇습니다.

>>> np.float32(1.0) + np.float32(2.0 ** -23) == np.float32(1.0)
False
>>> np.float32(1.0) + np.float32(2.0 ** -24) == np.float32(1.0)
True

$2^{-23} \approx 1.19\times10^{-7}$ 이 fp32의 머신 엡실론입니다. 이보다 작은 것을 1에 더하면 아무 일도 안 일어납니다.

지수가 1 오르면 같은 개수의 눈금이 두 배 넓은 구간을 나눠 맡으므로 간격도 두 배가 됩니다.

여기서 두 번째 규칙이 나옵니다. 부동소수점의 절대오차는 값에 비례해 커지고, 일정한 것은 상대오차뿐입니다. 반올림은 가장 가까운 눈금으로 가므로 오차가 간격의 절반을 넘지 않고, 따라서

$$\frac{|\hat{x} - x|}{|x|} \le \frac{\varepsilon}{2} = 2^{-(p+1)}$$

입니다. $\hat{x}$ 는 저장된 값, $x$ 는 참값입니다. 이 부등식이 이 글에 나오는 모든 사고의 출발점입니다.

형식	지수 비트	가수 비트 $p$	머신 엡실론 $2^{-p}$	상대오차 상한	표현 범위	유효 십진자리
fp64	11	52	2.22 × 10⁻¹⁶	1.11 × 10⁻¹⁶	±1.8 × 10³⁰⁸	약 16
fp32	8	23	1.19 × 10⁻⁷	5.96 × 10⁻⁸	±3.4 × 10³⁸	약 7
bf16	8	7	7.81 × 10⁻³	3.91 × 10⁻³	±3.4 × 10³⁸	약 2
fp16	5	10	9.77 × 10⁻⁴	4.88 × 10⁻⁴	±6.6 × 10⁴	약 3

같은 16비트, 다른 고장

bf16과 fp16은 둘 다 16비트인데 나누는 방식이 정반대입니다. bf16은 fp32의 지수 8비트를 그대로 두고 가수만 잘랐고, fp16은 지수를 5비트로 줄이는 대신 가수를 10비트 남겼습니다. 그래서 서로 다른 지점에서 터집니다.

값	fp16	bf16
1 × 10⁻⁸	0 (언더플로)	1.0012 × 10⁻⁸
1.00391	1.00391	1.0 (가수 부족)
30,000	30,000	29,952
70,000	inf (오버플로)	70,144
1 × 10¹⁰	inf	9.9992 × 10⁹

fp16은 범위에서 터집니다. 최대값이 65,504라 그 위는 전부 inf입니다. 어텐션이 실제로 이 벽을 만납니다 — softmax를 계산하면서 $e^{z}$ 를 구하는데, $z = 12$ 만 돼도 $e^{12} \approx 162{,}755$ 로 fp16의 상한을 넘습니다. bf16에서는 162,816으로 무사합니다. 로짓 자체도 시퀀스가 길어지고 스케일이 커지면 $-155{,}219$ 같은 값이 나오는데, fp16에서는 이 값이 이미 -inf입니다.

bf16은 정밀도에서 터집니다. 가수가 7비트뿐이라 $1 + 2^{-8}$ 은 그냥 1이 됩니다. 유효 십진자리가 두 자리 남짓이라, 큰 누적합에 작은 갱신을 더하는 일이 통째로 무시됩니다. 학습률 $10^{-4}$ 로 크기 1짜리 가중치를 갱신하면 bf16에서는 아무 일도 안 일어납니다.

그래서 실무의 정답은 「무엇을 쓰나」가 아니라 「무엇을 어디에 두나」입니다. 행렬 곱은 bf16으로 하고(범위가 넓어 오버플로 걱정이 없음) 누산과 가중치 갱신은 fp32로 하는 것이 혼합정밀도의 기본형입니다. 어디에 무엇을 두는지는 다음 글 「수치적으로 안정한 계산 패턴 모음」이 목록으로 맡습니다.

파괴적 상쇄

반올림 오차는 상대오차로 보면 $10^{-7}$ 수준이라 대개 무시할 만합니다. 그런데 그 오차가 결과의 맨 앞자리로 승격되는 연산이 하나 있습니다.

비슷한 크기의 두 수를 뺄 때입니다.

$a$ 와 $b$ 는 각각 일곱 자리를 믿을 수 있었습니다. 그런데 앞의 여섯 자리가 서로 같아서 빼는 순간 0으로 사라지고, 남는 것은 마지막 한 자리뿐입니다. 원래 그 자리 아래에 있던 반올림 오차가 이제 결과의 두 번째 자리가 됩니다. 이것을 파괴적 상쇄(catastrophic cancellation)라고 합니다.

>>> a, b = np.float32(1.2345678), np.float32(1.2345670)
>>> float(a - b)
7.152557373046875e-07        # 참값은 8.0e-07

상대오차가 11%입니다. 뺄셈 자체는 완벽하게 정확했습니다 — 두 저장값의 차이를 정확히 계산했습니다. 새어 나간 것은 입력이 이미 갖고 있던 오차이고, 상쇄가 그것을 확대했을 뿐입니다.

교과서적인 예가 이차방정식의 근입니다. $ax^2+bx+c=0$ 의 근을

$$x = \frac{-b \pm \sqrt{b^2-4ac}}{2a}$$

로 구하는데, $b^2 \gg 4ac$ 이면 $\sqrt{b^2-4ac} \approx |b|$ 라서 두 근 중 하나에서 $-b$ 와 $\sqrt{\cdot}$ 가 거의 같은 크기로 서로 뺍니다. $a=1$, $b=10^8$, $c=1$ 로 fp64에서 계산하면

	작은 근	큰 근
공식 그대로	−7.4506 × 10⁻⁹	−1.0 × 10⁸
참값	−1.0 × 10⁻⁸	−1.0 × 10⁸

작은 근이 25%나 틀렸습니다. fp64인데도 그렇습니다.

고치는 방법은 상쇄가 일어나는 뺄셈을 아예 없애는 것입니다. 근과 계수의 관계 $x_1x_2 = c/a$ 를 쓰면 됩니다.

def roots_stable(a, b, c):
    d = math.sqrt(b * b - 4 * a * c)
    q = -0.5 * (b + math.copysign(d, b))   # 부호를 맞춰 '더하기'로 만든다
    return q / a, c / q                    # 작은 근은 나눗셈으로 얻는다

copysign으로 $b$ 와 같은 부호를 붙였기 때문에 $b + \operatorname{sign}(b)\sqrt{\cdot}$ 는 같은 부호끼리의 덧셈이라 상쇄가 없습니다. 그렇게 안전하게 구한 큰 근에서 작은 근을 나눗셈으로 끌어냅니다. 이 방법으로 계산하면 작은 근이 정확히 $-1.0\times10^{-8}$ 입니다.

여기서 나오는 규칙입니다. 상쇄를 만나면 정밀도를 올리지 말고 식을 바꿉니다. fp32를 fp64로 바꾸는 것은 자릿수를 몇 개 더 벌어 줄 뿐 구조를 안 고칩니다. 위 예가 fp64에서도 25% 틀린 것이 그 증거입니다.

더하는 순서가 결과를 바꾼다

실수 덧셈에는 결합법칙이 있습니다. $(a+b)+c = a+(b+c)$ 입니다. 부동소수점 덧셈에는 없습니다. 각 단계마다 눈금으로 반올림이 들어가는데 그 눈금이 중간합의 크기에 따라 달라지기 때문입니다.

가장 뚜렷한 예입니다. fp32에서 $10^8$ 에 1을 백만 번 더해 봅니다.

f32 = np.float32
big, n = f32(1e8), 1_000_000

s = big
for _ in range(n):
    s = f32(s + f32(1.0))
print(float(s))            # 100000000.0  — 아무것도 안 더해졌다

s = f32(0.0)
for _ in range(n):
    s = f32(s + f32(1.0))
print(float(s + big))      # 101000000.0  — 참값

이유는 눈금 간격입니다. fp32에서 $10^8$ 근처의 눈금 간격은 8입니다. 거기에 1을 더하면 가장 가까운 눈금은 여전히 자기 자신이라 반올림이 원래 값을 돌려줍니다. 백만 번을 반복해도 마찬가지입니다. 반면 작은 것부터 더하면 중간합이 $10^6$ 근처에 있는 동안 간격이 0.0625라 1이 제대로 쌓입니다.

이 현상에 이름이 있습니다 — 더할 값이 중간합의 눈금 간격보다 작아 통째로 버려지는 것을 흡수(absorption)라고 합니다.

무작위 데이터에서는 어떨까요. $[0,1)$ 균등 난수 1,000만 개를 여섯 가지 방법으로 fp32에서 더하고 fp64 참값과 비교했습니다.

방법	결과	상대오차
작은 것부터 (오름차순 정렬)	5,057,480.5	1.15 × 10⁻²
뒤에서부터	5,000,322.0	1.08 × 10⁻⁴
앞에서부터	4,999,722.5	1.14 × 10⁻⁵
numpy 기본 (쌍대 합산)	4,999,779.5	2.41 × 10⁻⁸
카한 합산	4,999,779.5	2.41 × 10⁻⁸
fp32 입력 → fp64 누산	4,999,779.62	2.24 × 10⁻¹⁵

참값은 4,999,779.6205입니다. 같은 수 1,000만 개를 순서만 바꿔 더했는데 오차가 백만 배 갈립니다.

세 가지를 짚습니다.

「작은 것부터 더하면 낫다」는 통념이 여기서는 가장 나빴습니다. 시드를 세 번 바꿔도 1.151.16%로 재현됐고, 방향도 항상 과대평가였습니다. 원인은 정렬이 반올림을 한쪽으로 몰았기 때문입니다. 마지막 백만 개는 값이 0.91.0 사이인데 그때 중간합은 500만 근처라 눈금 간격이 0.5입니다. 0.91.0을 0.5 눈금에 맞추면 나머지가 항상 0.40.5 구간에 떨어지고, 그 구간은 언제나 위쪽 눈금이 더 가깝습니다. 반올림이 매번 같은 방향으로 일어나 오차가 상쇄되지 않고 쌓입니다. 반올림이 오차를 지우는 것은 남는 자리가 고르게 흩어질 때뿐이고, 정렬은 그 고름을 깨는 조작입니다.

numpy의 기본 합산이 소박한 반복문보다 천 배 정확합니다. np.sum은 배열을 반으로 쪼개 각각 더한 뒤 합치는 쌍대 합산(pairwise summation)을 씁니다. 중간합이 언제나 원소 개수에 비례하는 크기에 머물러 흡수가 덜 일어나고, 오차가 $O(n)$ 이 아니라 $O(\log n)$ 으로 쌓입니다. 손으로 for 루프를 써서 더하는 것이 라이브러리보다 나쁜 이유가 이것입니다.

카한 합산은 버려진 부분을 다음 항에 되돌립니다.

def kahan(a):
    s = f32(0.0)
    c = f32(0.0)                    # 지난번에 흡수돼 사라진 몫
    for v in a:
        y = f32(v - c)              # 잃었던 만큼 미리 보태 둔다
        t = f32(s + y)              # 여기서 다시 일부가 흡수된다
        c = f32(f32(t - s) - y)     # 실제로 들어간 양과의 차이가 이번에 잃은 몫
        s = t
    return float(s)

핵심은 세 번째 줄입니다. t - s는 「이번에 실제로 반영된 양」이고, 거기서 y(넣으려던 양)를 빼면 흡수돼 사라진 몫이 나옵니다. 그것을 c에 담아 두었다가 다음 항에서 되돌립니다. 이 트릭 덕분에 오차가 $n$ 과 무관한 상수로 묶이지만, 덧셈이 네 번으로 늘고 순차적이라 벡터화가 안 됩니다.

그리고 가장 싼 방법은 맨 아래 줄입니다. 입력은 fp32로 두고 누산기만 fp64로 두면 오차가 $10^{-15}$ 로 떨어집니다. np.sum(x, dtype=np.float64) 한 줄이고, GPU의 텐서 코어가 bf16 입력을 fp32 누산기로 처리하는 것도 정확히 이 발상입니다. 정밀도는 저장하는 곳이 아니라 쌓이는 곳에 필요합니다.

그래서 GPU 개수를 바꾸면 왜 손실이 달라지나

서두의 두 현상으로 돌아갑니다.

손실이 소수 넷째 자리에서 갈린 것은 GPU 개수가 부분합의 묶음을 바꿨기 때문입니다. GPU 두 대는 배치를 둘로 쪼개 각각 더한 뒤 합치고, 네 대는 넷으로 쪼갭니다. 수학적으로는 같은 합이지만 부동소수점에서는 결합법칙이 없으므로 다른 값이 나옵니다. 위 표에서 같은 데이터가 묶는 방법에 따라 4,999,722부터 5,057,480까지 갈렸던 것과 같은 일입니다. 이것은 버그가 아니고, 없앨 수도 없습니다. 재현성이 꼭 필요하면 누산 순서를 고정하는 결정론적 커널을 켜야 하고 그만큼 느려집니다.

fp16 어텐션이 nan이 된 것은 범위 문제입니다. 지수 5비트짜리 형식의 상한 65,504를 $e^{z}$ 가 넘었습니다. 이쪽은 고칠 수 있고 고치는 방법이 이미 표준입니다 — softmax를 계산하기 전에 최댓값을 빼면 지수가 절대 양수가 되지 않아 오버플로가 원천적으로 사라집니다.

정리하면 이 글의 세 고장은 성격이 다릅니다. 반올림은 피할 수 없고 상대오차 상한으로 관리하는 것입니다. 상쇄는 식을 바꿔 피하는 것입니다. 흡수와 순서 의존은 누산기를 어디에 두느냐로 관리하는 것입니다. 셋을 구별하는 것이 「정밀도를 올려 보자」와 「구조를 고치자」를 가르는 자리입니다.

정리

• 0.1은 이진법에서 유한 소수가 아니다. 분모 10의 소인수에 5가 있어서다. 0.1 + 0.2 != 0.3은 덧셈이 아니라 입력이 이미 틀려서 생긴다. 부동소수점을 ==로 비교하지 않는다.
• 실수는 $(-1)^s \times 1.f \times 2^e$ 로 저장된다. 지수 비트가 범위를, 가수 비트가 정밀도를 정한다.
• 가수 $p$ 비트는 $[1,2)$ 를 $2^p$ 등분한다. 그 간격이 머신 엡실론 $2^{-p}$ 이고, 반올림 상대오차는 그 절반을 넘지 않는다. fp32는 $2^{-23} \approx 1.19\times10^{-7}$ 이다.
• 지수가 1 오르면 눈금 간격이 두 배가 된다. 절대오차는 값에 비례해 커지고 상대오차만 일정하다. fp32에서 1 근처 간격은 $1.2\times10^{-7}$ 인데 $10^8$ 근처에서는 8이다.
• bf16은 fp32의 지수를 그대로 두고 가수를 잘랐고, fp16은 반대다. fp16은 65,504에서 오버플로로 터지고($e^{12}$ 가 이미 넘는다), bf16은 $1+2^{-8}$ 이 그냥 1이 되는 정밀도로 터진다.
• 파괴적 상쇄는 비슷한 두 수를 뺄 때 겹치던 앞자리가 사라지고 반올림 오차가 맨 앞자리로 올라오는 것이다. fp32에서 두 수를 빼 상대오차 11%가 나왔다.
• 상쇄를 만나면 정밀도를 올리지 말고 식을 바꾼다. 이차방정식 근의 공식은 fp64에서도 작은 근이 25% 틀렸고, copysign으로 덧셈으로 바꾸면 정확해진다.
• 부동소수점 덧셈에는 결합법칙이 없다. fp32에서 $10^8$ 에 1을 백만 번 더하면 눈금 간격이 8이라 하나도 안 쌓인다(흡수).
• 1,000만 개를 순서만 바꿔 더하니 오차가 $1.15\times10^{-2}$ 부터 $2.24\times10^{-15}$ 까지 갈렸다. 오름차순 정렬이 가장 나빴다 — 반올림이 매번 같은 방향으로 몰려 상쇄되지 않는다.
• np.sum의 쌍대 합산은 소박한 반복문보다 천 배 정확하고, 카한 합산은 흡수된 몫을 다음 항에 되돌린다. 가장 싼 방법은 누산기만 fp64로 두는 것이다 — 한 줄로 $10^{-15}$ 다.
• GPU 개수를 바꾸면 손실이 달라지는 것은 부분합의 묶음이 달라져서다. 버그가 아니고 없앨 수도 없다.

---

읽어주셔서 감사합니다. 😊