수치적으로 안정한 계산 패턴 모음

교차 엔트로피 손실을 직접 짰습니다. 식이 $-\log p_y$ 이고 $p$ 는 softmax이니 두 줄이면 끝납니다. 그런데 몇백 스텝 뒤에 손실이 inf가 되고, F.cross_entropy로 바꾸니 똑같은 데이터에서 멀쩡히 돕니다. 코드는 수식 그대로였습니다.

p = np.exp(z) / np.exp(z).sum()
loss = -np.log(p[y])

로짓이 [1000, 3, 5]일 때 이 코드는 inf를 내놓습니다. [-800, -900, -1000]이면 nan입니다. 안정한 구현은 각각 995.0과 200.0을 정확히 돌려줍니다. 수학적으로 같은 식이 계산으로는 다른 식이라는 것이고, 프레임워크는 그 차이를 알고 다른 쪽을 씁니다.

지난 글에서 그 차이가 어디서 나오는지 — 반올림, 파괴적 상쇄, 흡수 — 를 셋으로 갈랐습니다. 이 글은 그 셋을 막는 처방 일곱 개를 목록으로 세우고, 각각이 무엇을 막는지 하나씩 재 봅니다.

최댓값 빼기

softmax는 로짓 전체를 같은 상수만큼 밀어도 값이 안 변합니다.

$$\frac{e^{z_i - c}}{\sum_j e^{z_j - c}} = \frac{e^{-c}e^{z_i}}{e^{-c}\sum_j e^{z_j}} = \frac{e^{z_i}}{\sum_j e^{z_j}}$$

분자와 분모에 똑같이 $e^{-c}$ 가 곱해져 약분됩니다. 이 자유를 써서 $c = \max_j z_j$ 로 두면 모든 지수가 0 이하가 되어 $e^{z_i-c} \le 1$ 이고, 오버플로가 원천적으로 사라집니다.

아래쪽도 같이 고쳐집니다. 로짓이 전부 $-1000$ 근처면 그대로 계산했을 때 모든 지수가 0으로 언더플로해서 $0/0$ 이 되는데, 최댓값을 빼면 적어도 하나는 $e^0 = 1$ 이라 분모가 절대 0이 아닙니다.

같은 트릭이 log-sum-exp에도 그대로 붙습니다.

$$\log\sum_j e^{z_j} = c + \log\sum_j e^{z_j - c}, \qquad c = \max_j z_j$$

log-sum-exp를 다룬 글에서 유도한 그대로입니다. 그리고 교차 엔트로피는 softmax를 거치지 않고 이것 하나로 바로 씁니다.

$$-\log p_y = -\left(z_y - \log\sum_j e^{z_j}\right)$$

나눗셈도 로그도 사라졌습니다. 서두의 inf와 nan이 이 한 줄로 없어집니다. F.cross_entropy가 softmax를 따로 계산하지 않고 log_softmax와 nll_loss를 묶어 쓰는 이유가 이것입니다.

막는 고장: 오버플로, 언더플로, $0/0$.

로그 공간에서 계산하기

확률을 여러 개 곱하면 값이 지수적으로 작아집니다. 0.7을 2,000번 곱하면

>>> np.prod(np.full(2000, 0.7))
1.57065220561795e-310

fp64의 최소 정규수가 $2.2\times10^{-308}$ 이라 이미 비정규 영역에 들어갔고, 항이 조금만 더 많으면 그냥 0입니다. 0이 되고 나면 되돌릴 방법이 없습니다.

로그를 취하면 곱셈이 덧셈이 됩니다.

$$\log\prod_i p_i = \sum_i \log p_i$$

같은 계산이 $-713.35$ 로 나옵니다. fp64가 $\pm1.8\times10^{308}$ 까지 담으니 항이 수백만 개여도 여유가 있습니다. 로그 공간은 곱셈을 덧셈으로 바꾸는 것이 아니라, 지수적으로 퍼지는 값의 범위를 선형으로 눌러 담는 것입니다.

곱셈은 덧셈이 되고 나눗셈은 뺄셈이 되는데, 덧셈만 곤란합니다. $\log(a+b)$ 를 $\log a$ 와 $\log b$ 로 구해야 하는데 그것이 정확히 log-sum-exp입니다. 그래서 로그 공간으로 옮기면 최댓값 빼기가 반드시 따라옵니다.

막는 고장: 언더플로, 오버플로.

log1p와 expm1

$x$ 가 아주 작을 때 $\log(1+x)$ 를 그대로 계산하면 두 번 손해를 봅니다. 먼저 $1+x$ 를 만들면서 $x$ 의 아랫자리가 1의 눈금에 흡수되고, 그다음 $\log$ 가 1 근처에서 $0$ 을 돌려줍니다.

$x$	log(1 + x)	log1p(x)	참값
10⁻⁸	1.000000 × 10⁻⁸	1.000000 × 10⁻⁸	1.0 × 10⁻⁸
10⁻¹²	1.000089 × 10⁻¹²	1.000000 × 10⁻¹²	1.0 × 10⁻¹²
10⁻¹⁶	0	1.000000 × 10⁻¹⁶	1.0 × 10⁻¹⁶

log1p는 $1+x$ 를 만들지 않고 급수 $\log(1+x) = x - x^2/2 + x^3/3 - \cdots$ 를 직접 씁니다. 중간에 1이 등장하지 않으니 흡수될 것이 없습니다. expm1(x)는 $e^x - 1$ 의 짝으로, $e^x$ 가 1에 가까울 때 1을 빼면서 일어나는 상쇄를 피합니다.

쓰는 자리가 생각보다 많습니다. 로지스틱 손실 $\log(1+e^{-z})$ 가 그렇고, 작은 확률의 로그 오즈, 이자·성장률 계산이 그렇습니다.

막는 고장: 흡수, 파괴적 상쇄.

분산의 두 공식

분산 공식은 두 가지로 씁니다.

$$s^2 = \frac{1}{n-1}\sum_i (x_i - \bar{x})^2 \qquad\text{또는}\qquad s^2 = \frac{1}{n-1}\left(\sum_i x_i^2 - \frac{(\sum_i x_i)^2}{n}\right)$$

오른쪽이 매력적입니다. $\sum x_i$ 와 $\sum x_i^2$ 만 들고 다니면 되니 데이터를 한 번만 훑고 스트리밍으로도 됩니다. 왼쪽은 평균을 먼저 구해야 해서 두 번 훑습니다.

그런데 오른쪽은 파괴적 상쇄를 정면으로 밟습니다. $\sum x_i^2$ 와 $(\sum x_i)^2/n$ 은 둘 다 크고 거의 같은데, 답인 분산은 그 차이라서 작습니다. 데이터 $1,2,3,4,5$ 에 상수를 더해 가며 fp64로 재 보면 이렇습니다.

데이터	한 번 훑기	두 번 훑기	Welford	참값
1~5	2.5	2.5	2.5	2.5
+10⁴	2.5	2.5	2.5	2.5
+10⁸	2.0	2.5	2.5	2.5
+10⁹	0.0	2.5	2.5	2.5

값 자체는 안 변하고 평균만 옮겼는데 한 번 훑기 공식은 $10^9$ 근처에서 0을 내놓습니다. fp64인데도 그렇습니다. 데이터가 평균에서 조금씩 떨어져 있는 형태 — 온도, 센서 값, 타임스탬프, 큰 상수가 실린 특성 — 가 정확히 이 모양입니다.

한 번 훑기와 정확성을 둘 다 가지려면 Welford 알고리즘을 씁니다. 평균과 제곱합을 매 항마다 갱신하되, 큰 값끼리 빼는 자리를 만들지 않습니다.

def welford(xs):
    n = 0
    mean = 0.0
    M2 = 0.0
    for x in xs:
        n += 1
        d = x - mean          # 평균에서의 편차 — 크기가 작다
        mean += d / n
        M2 += d * (x - mean)  # 갱신 전 편차 × 갱신 후 편차
    return M2 / (n - 1)

핵심은 d = x - mean입니다. 큰 수 둘을 빼는 것은 여기서도 마찬가지지만, 그 결과를 제곱해서 쌓기 전에 빼기 때문에 상쇄가 한 번만 일어나고 확대되지 않습니다. 위 표의 마지막 열이 그 결과이고, $10^9$ 에서도 2.5입니다. numpy와 pandas의 분산은 두 번 훑기를, 스트리밍 통계 라이브러리는 Welford를 씁니다.

막는 고장: 파괴적 상쇄.

정규화의 ε

LayerNorm이나 BatchNorm의 서명에는 eps=1e-5 같은 인자가 있습니다.

$$\hat{x}_i = \frac{x_i - \mu}{\sqrt{\sigma^2 + \varepsilon}}$$

분모가 0이 되는 것을 막는 값이라고 설명되는데, 그것은 절반입니다. 실제로 두 가지를 막습니다.

입력	$\varepsilon = 0$	$\varepsilon = 10^{-5}$
[3, 3, 3, 3] (완전히 같음)	[nan, nan, nan, nan]	[0, 0, 0, 0]
[1.0000001, 1, 1, 1] (거의 같음)	[2, 0, 0, 0]	[0, 0, 0, 0]

첫 줄이 알려진 쪽입니다 — 분산이 정확히 0이라 $0/0$ 이 됩니다.

둘째 줄이 진짜 문제입니다. nan이 아니라 그럴듯한 숫자가 나옵니다. 입력의 차이는 $10^{-7}$ 밖에 안 되는데 출력은 2가 됐습니다. 분산이 아주 작으면 $\sqrt{\sigma^2}$ 도 아주 작고, 그 작은 수로 나누면 입력에 실려 있던 반올림 오차가 그대로 확대됩니다. $10^{-7}$ 짜리 잡음이 2가 되는 증폭입니다. nan은 눈에 띄어 잡히기라도 하는데, 이쪽은 그냥 잘못된 값으로 흘러 다음 층으로 갑니다.

$\varepsilon$ 은 분모에 하한을 두어 이 증폭에 천장을 씌웁니다. 분산이 $\varepsilon$ 보다 훨씬 크면 아무 일도 안 하고, 작아지면 나눗셈이 $\sqrt{\varepsilon}$ 으로 나누는 것에 수렴합니다. 작은 상수를 더해 두는 것이 아니라, 나눗셈의 이득에 상한을 거는 장치입니다.

$\sqrt{\sigma^2+\varepsilon}$ 이지 $\sqrt{\sigma^2}+\varepsilon$ 이 아니라는 점도 중요합니다. 제곱근 안에 넣어야 $\sigma^2$ 와 같은 단위에서 비교되고, $\varepsilon$ 을 「분산의 하한」으로 읽을 수 있습니다.

막는 고장: $0/0$, 그리고 작은 분모에 의한 오차 증폭.

누산은 fp32에, 저장은 bf16에

지난 글에서 확인한 대로, 정밀도가 필요한 곳은 값을 담는 자리가 아니라 값이 쌓이는 자리입니다. 한 번 곱하고 버릴 값은 bf16으로 충분하고, 수천 개를 더하거나 작은 갱신을 오래 쌓는 자리는 fp32여야 합니다.

가장 단적인 예가 가중치 갱신입니다. 크기 1인 가중치에 학습률 $10^{-4}$ 짜리 갱신을 천 번 더해 봤습니다.

w32 = np.float32(1.0)
w16 = bf16(np.float32(1.0))
for _ in range(1000):
    w32 = np.float32(w32 - lr * g)
    w16 = bf16(np.float32(w16 - lr * g))
# w32 → 0.89998,  w16 → 1.0  (참값 0.9)

bf16 가중치는 천 스텝 동안 한 번도 안 움직였습니다. 머신 엡실론이 $7.8\times10^{-3}$ 이라 $10^{-4}$ 짜리 갱신이 매번 통째로 흡수됩니다. 학습이 진행되는 것처럼 보이면서 실제로는 아무 일도 안 일어나는 고장이고, 손실 곡선만 봐서는 안 보입니다.

그래서 혼합정밀도는 「무엇을 쓰나」가 아니라 「무엇을 어디에 두나」입니다.

마스터 가중치를 fp32로 따로 들고 다니면서 순전파용 bf16 사본을 매 스텝 만드는 것이 표준 구성입니다. 메모리를 더 쓰는 것처럼 보이지만, 활성값이 차지하는 양이 가중치보다 훨씬 크기 때문에 전체로는 여전히 이득입니다.

막는 고장: 흡수, 오버플로.

역행렬을 만들지 않는다

$A\mathbf{x} = \mathbf{b}$ 를 풀 때 수식은 $\mathbf{x} = A^{-1}\mathbf{b}$ 지만, 코드에서 inv(A) @ b를 쓰지 않고 solve(A, b)를 씁니다. 이유는 속도만이 아닙니다.

먼저 얼마나 나빠질 수 있는지를 재는 값이 필요합니다. 조건수는 입력의 상대오차가 출력에서 최대 몇 배로 커지는지를 나타내는 수이고, 행렬의 가장 큰 특잇값과 가장 작은 특잇값의 비입니다.

$$\kappa(A) = \frac{\sigma_{\max}}{\sigma_{\min}}$$

특잇값 분해를 다룬 글에서 본 대로, $A$ 를 곱하는 것은 단위원을 타원으로 만드는 일입니다. 조건수는 그 타원이 얼마나 납작한지입니다.

조건수가 $10^k$ 이면 답에서 십진 $k$ 자리를 잃습니다. fp64가 가진 것이 16자리이므로, 조건수 $10^{14}$ 짜리 문제를 fp64로 풀면 두 자리만 남습니다.

힐베르트 행렬로 두 방법을 비교했습니다. 답이 전부 1이 되도록 $\mathbf{b} = A\mathbf{1}$ 로 만들어 두고 풀었습니다.

크기	조건수	inv(A) @ b 최대오차	solve(A, b) 최대오차
6 × 6	1.5 × 10⁷	1.98 × 10⁻⁹	6.65 × 10⁻¹¹
10 × 10	1.6 × 10¹³	8.13 × 10⁻³	4.42 × 10⁻⁵
14 × 14	3.2 × 10¹⁷	1.17 × 10³	7.49

solve가 매번 30배에서 180배 정확합니다. 이유는 inv가 계산을 두 번 하기 때문입니다 — 먼저 $A^{-1}$ 을 구하면서 반올림 오차를 쌓고, 그다음 곱하면서 또 쌓습니다. solve는 LU 분해를 한 번 해서 바로 답으로 갑니다. 중간 결과인 $A^{-1}$ 을 만들지 않으므로 거기서 생길 오차도 없습니다.

같은 원리의 더 큰 사례가 최소제곱입니다. 정규방정식 $X^\top X\mathbf{w} = X^\top\mathbf{y}$ 를 만들어 푸는 것이 교과서 유도인데, $X^\top X$ 를 만드는 순간 조건수가 제곱됩니다.

$$\kappa(X^\top X) = \kappa(X)^2$$

특잇값이 제곱되기 때문입니다. $\sigma_{\max}/\sigma_{\min}$ 이 $\sigma_{\max}^2/\sigma_{\min}^2$ 이 됩니다.

거의 같은 열 두 개가 있는 $50\times8$ 설계행렬로 재 보니 $\kappa(X) = 1.78\times10^7$ 인데 $\kappa(X^\top X) = 3.31\times10^{14}$ 이었고, 계수 오차가 $4.5\times10^{-3}$ 대 $6.9\times10^{-10}$ 으로 650만 배 갈렸습니다. numpy.linalg.lstsq가 정규방정식을 만들지 않고 $X$ 를 그대로 QR 분해하는 이유입니다.

막는 고장: 조건수에 의한 오차 확대.

목록으로

일곱 개를 한자리에 놓으면 이렇습니다.

패턴	막는 고장	프레임워크에서
최댓값 빼기	오버플로 · 언더플로 · $0/0$	log_softmax, logsumexp, cross_entropy
로그 공간 계산	언더플로	log_prob, HMM·CRF의 전방 알고리즘
log1p · expm1	흡수 · 상쇄	softplus, logsigmoid
두 번 훑기 또는 Welford	파괴적 상쇄	var, std, 러닝 통계
정규화의 $\varepsilon$	$0/0$ · 오차 증폭	LayerNorm(eps=1e-5), RMSNorm
fp32 누산	흡수	텐서 코어, GradScaler, 마스터 가중치
solve · QR	조건수 확대	linalg.solve, lstsq, cholesky_solve

세 갈래로 다시 묶으면 더 짧아집니다.

범위가 모자란 것은 좌표를 옮겨 고칩니다. 최댓값 빼기와 로그 공간이 그렇고, 둘 다 「값 자체를 바꾸지 않는 변환을 골라 안전한 구간으로 옮긴다」는 같은 일입니다.

상쇄가 있는 것은 식을 바꿔 고칩니다. log1p, Welford, solve가 그렇고, 정밀도를 올려서는 안 고쳐집니다.

쌓이는 것은 누산기를 키워 고칩니다. fp32 누산이 그렇고, 정규화의 $\varepsilon$ 은 나눗셈의 이득에 천장을 씌운다는 점에서 조금 다른 갈래입니다.

그래서 직접 짜면 안 되나

서두의 교차 엔트로피로 돌아갑니다. 코드가 틀린 것이 아니라, 수식을 그대로 옮기면 좌표를 옮기는 단계가 빠집니다. 수학에는 $e^{1000}$ 이라는 값이 문제없이 존재하고, 계산에는 없습니다.

그래서 직접 짜지 말라는 뜻은 아닙니다. 짤 때 확인할 목록이 있다는 뜻입니다.

• 지수를 계산하는가 → 지수가 양수가 되지 않도록 상수를 뺄 수 있는가
• 확률을 여러 개 곱하는가 → 로그 공간으로 옮길 수 있는가
• 비슷한 크기의 두 수를 빼는가 → 그 뺄셈을 없앤 동치 식이 있는가
• 작은 값으로 나누는가 → 분모에 하한이 있는가
• 많은 항을 더하는가 → 누산기가 입력보다 넓은가
• 행렬의 역을 만드는가 → 그것 없이 푸는 함수가 있는가

여섯 줄이고, 대부분의 자리에서는 「없다」가 답이라 그냥 지나갑니다. 걸리는 자리가 정확히 프레임워크가 몰래 다르게 계산하고 있던 자리입니다.

정리

• 최댓값 빼기는 softmax가 로짓 전체를 밀어도 값이 안 변한다는 성질을 쓴다. 오버플로와 $0/0$ 을 동시에 막고, 교차 엔트로피는 아예 $-(z_y - \operatorname{logsumexp}(z))$ 로 써서 나눗셈과 로그를 없앤다.
• 로그 공간은 0.7을 2,000번 곱해 $1.6\times10^{-310}$ 이 되는 것을 $-713.35$ 로 바꾼다. 대신 덧셈이 log-sum-exp가 되므로 최댓값 빼기가 따라온다.
• log1p·expm1은 1 근처에서 흡수와 상쇄를 피한다. $x = 10^{-16}$ 에서 log(1+x)는 0이고 log1p(x)는 정확하다.
• 분산의 한 번 훑기 공식은 파괴적 상쇄를 정면으로 밟는다. 데이터에 $10^9$ 을 더하기만 해도 fp64에서 0이 나온다. 스트리밍이 필요하면 Welford를 쓴다.
• 정규화의 $\varepsilon$ 이 막는 것은 $0/0$ 만이 아니다. 분산이 아주 작을 때 $10^{-7}$ 짜리 잡음이 2로 증폭되는 것을 막는다. nan보다 이쪽이 위험하다 — 그럴듯한 숫자라서 안 잡힌다.
• 정밀도는 쌓이는 곳에 둔다. bf16 가중치에 $10^{-4}$ 갱신을 천 번 더해도 값이 1.0에서 안 움직인다. 마스터 가중치와 누산기는 fp32다.
• 조건수 $\sigma_{\max}/\sigma_{\min}$ 은 단위원이 얼마나 납작한 타원이 되는지다. $10^k$ 이면 십진 $k$ 자리를 잃는다.
• inv(A) @ b 대신 solve(A, b). 힐베르트 행렬에서 30~180배 정확했다. 중간 결과 $A^{-1}$ 을 안 만들면 거기서 생길 오차도 없다.
• 정규방정식은 조건수를 제곱한다. $\kappa(X)=1.78\times10^7$ 이 $\kappa(X^\top X)=3.31\times10^{14}$ 이 되고, 계수 오차가 650만 배 갈렸다. lstsq는 $X$ 를 그대로 QR 분해한다.
• 셋으로 묶으면 — 범위가 모자라면 좌표를 옮기고, 상쇄가 있으면 식을 바꾸고, 쌓이는 것은 누산기를 키운다.

---

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