Numerical Precision and Stability in Practice

1. What is it?#

The practical discipline of keeping an inference system numerically correct while running it in the lowest precision you can get away with, with fused kernels, on multiple GPUs, at varying batch sizes.

The core uncomfortable fact: your inference system is nondeterministic and slightly different from the reference implementation, and that is normal. The engineering task is bounding the difference, not eliminating it.

2. Why does it exist?#

Because floating point addition is not associative:

a, b, c := 0.1, 0.2, 0.3
fmt.Println((a+b)+c == a+(b+c)) // false

And because everything in a GPU inference stack changes the order of additions: batch size (different GEMM kernel → different reduction order), tensor parallel degree (different partial sums), fusion (different intermediate rounding), library version.

The result: the same prompt can produce different text on the same hardware with a different batch composition. Users notice. You need to be able to explain it.

3. Simple analogy#

Rounding a long bill. Add the items in a different order and round to cents at each step, and you can end up a cent apart. On a small bill it doesn’t matter. But if a cent difference flips a decision — “does this exceed the threshold?” — the outcome changes categorically.

In LLMs, sampling is the threshold. A 1e-4 difference in a logit rarely matters, until two candidate tokens are nearly tied, and then it changes the word — and every word after it.

4. Tiny example#

// order.go — the same dot product, summed in two different orders.
package main

import (
	"fmt"
	"math/rand"
)

const d = 4096

// One accumulator, left to right — how a simple kernel sums when the row is computed alone.
func dotSequential(x, w []float32) (s float32) {
	for i := range x {
		s += x[i] * w[i]
	}
	return s
}

// Eight accumulators combined at the end — how a wide, batched kernel sums.
func dotLanes(x, w []float32) float32 {
	var lane [8]float32
	for i := range x {
		lane[i%8] += x[i] * w[i]
	}
	return ((lane[0] + lane[1]) + (lane[2] + lane[3])) + ((lane[4] + lane[5]) + (lane[6] + lane[7]))
}

func main() {
	rng := rand.New(rand.NewSource(0))
	x := make([]float32, d)
	for i := range x {
		x[i] = float32(rng.NormFloat64())
	}
	alone, batched := make([]float32, d), make([]float32, d) // d "logits"
	var maxDiff float32
	identical := true
	for row := 0; row < d; row++ {
		w := make([]float32, d)
		for i := range w {
			w[i] = float32(rng.NormFloat64())
		}
		alone[row], batched[row] = dotSequential(x, w), dotLanes(x, w)
		diff := alone[row] - batched[row]
		maxDiff = max(maxDiff, diff, -diff)
		identical = identical && diff == 0
	}
	fmt.Println("max abs diff:", maxDiff)
	fmt.Println("identical   :", identical)

	argmax := func(v []float32) (best int) {
		for i := range v {
			if v[i] > v[best] {
				best = i
			}
		}
		return best
	}
	fmt.Println("argmax same :", argmax(alone) == argmax(batched))
}

Output: max diff around 5e-4 in FP32, and identical: false. Nothing about the inputs changed — only the order of the additions.

On a GPU this happens without you asking for it. Compute one row alone, then the same row inside a batch of 32: in FP16 the results differ by ~1e-2. The same input produced a different result because other requests were in the batch. A different batch size selects a different cuBLAS kernel with a different reduction order.

Now the consequence:

fmt.Println("argmax same :", argmax(alone) == argmax(batched))
// usually true, but not always — and when it's false, the whole generation diverges

5. Technical explanation#

The sources of nondeterminism, ranked#

1. Batch composition        different GEMM kernel/reduction order        ~1e-3 in FP16
2. Tensor parallelism       partial sums combined in nondeterministic order ~1e-3
3. Atomics in reductions    split-K GEMM, some attention kernels          ~1e-4
4. Fusion                   different intermediate rounding               ~1e-4
5. Library/driver version   different kernel selection                    varies
6. Sampling RNG             seed and generator state                      total
7. Prefix cache hits        recomputed vs cached KV differ slightly       ~1e-4

Items 1-5 are all “the arithmetic is done in a different valid order.” Item 6 is intentional. Item 7 surprises people: a request served from a prefix cache gets slightly different KV than one that recomputed it, so the same prompt can generate differently depending on cache state.

The stability checklist for LLM inference#

✓ Softmax with max subtraction, in FP32
✓ Attention score accumulation in FP32
✓ RMSNorm/LayerNorm reduction in FP32
✓ Residual stream in BF16 or FP32, not FP16 (magnitudes grow with depth)
✓ Logits upcast to FP32 before sampling
✓ Masks as -inf, not a large finite negative
✓ eps values matching the reference implementation
✓ No accumulation of many values in FP16/BF16
✓ RoPE sin/cos tables computed in FP32

Each line is a real bug someone has shipped.

Why FP16 fails at depth#

With pre-norm transformers the residual stream is never normalized, so its magnitude grows roughly with sqrt(L):

Layer 1:   |x| ~ 1
Layer 40:  |x| ~ 20-60
Layer 80:  |x| ~ 50-200

Still well inside FP16’s range (65504) — but intermediate values inside the FFN, especially the SwiGLU product, can be 10-100x the residual magnitude. At 80 layers with an unlucky input, FP16 overflows. BF16 has the same range as FP32 and simply doesn’t have this failure mode. That’s why BF16 won.

Validating a precision change#

LEVEL 1 — numerical
  max/mean absolute and relative difference in logits vs FP32 reference
  KL divergence between output distributions
  Target: KL < 0.01, top-1 agreement > 99%

LEVEL 2 — task
  perplexity on held-out data (necessary, far from sufficient)
  standard benchmarks (MMLU, GSM8K, HumanEval) on YOUR domain

LEVEL 3 — behavioral
  long-form generation quality (500+ tokens) — where compounding shows
  instruction following
  refusal/safety behavior (quantization can shift this!)
  structured output validity (JSON schema adherence)

LEVEL 4 — production
  A/B test with real traffic and human or LLM-judge evaluation

Level 1 alone is not enough, and it’s where most teams stop. A model can have tiny logit differences and meaningfully worse multi-step reasoning, because errors compound across 500 autoregressive steps.

Achieving determinism (when you must)#

1. Fix the batch size (or use batch-invariant kernels)
2. Fix the tensor parallel degree
3. Disable atomics / split-K in kernels
4. Fix all library versions
5. Seed the RNG per request
6. Disable prefix caching
7. Set torch.use_deterministic_algorithms(True)

Cost: substantial — often 20-50% throughput. Determinism is expensive; decide whether you actually need it. Usually you need reproducibility for debugging (achievable with a deterministic debug mode) rather than production determinism.

6. Under the hood#

Where the differences physically arise: a GEMM with M=1 might use a GEMV kernel that accumulates serially in one thread. With M=32 it uses a tiled kernel where 8 warps each compute partial sums that are then combined via shared memory. Different tree shape, different rounding.

Split-K goes further: partial results from different thread blocks are combined with atomic adds in an order determined by scheduling — genuinely nondeterministic run to run.

7. Performance implications#

• FP32 softmax/reductions cost almost nothing (they’re small relative to the tensor) and buy correctness. Always do it.
• Determinism costs 20-50%. Don’t pay it unless required.
• Precision reduction is the highest-leverage optimization available — 2-4x — which is why the validation discipline matters so much.

8. Production implications#

• Document the exact numerical configuration in your model registry: dtype, quantization method, TP degree, kernel backend, library versions.
• Never change precision without Level 3 validation.
• Expect and explain nondeterminism to stakeholders before someone files it as a bug. Have a written answer ready.
• Provide a deterministic debug mode (batch 1, TP 1, seeded, no prefix cache) for reproducing user-reported issues.
• Monitor for NaN/inf in production. A cheap periodic check on logits catches numerical degradation before users do.
• Run quality evaluation in CI on the exact serving configuration, not on the reference implementation.

9. Common mistakes#

Using FP16 for a deep model. Range failures at depth.

Validating with perplexity only. Insensitive to what users care about.

Testing only short generations. Misses compounding.

Assuming the quantized model that scores well on MMLU is fine for your task. Domain shift.

Promising determinism you can’t deliver.

Accumulating in low precision in a custom kernel.

Forgetting that prefix caching changes results.

10. Hands-on exercise#

A. Reproduce nondeterminism. Run the section 4 example. Then run the same prompt through a real inference server 20 times at different concurrency levels. How often does the output differ?

B. Build a validation harness. Write a script that, given two model configurations, reports: max logit diff, KL divergence, top-1 agreement, perplexity on a held-out set, and side-by-side generations for 20 prompts. Keep this — you’ll use it for every quantization decision.

C. FP16 depth failure. Take a deep model in FP16 and find an input that produces inf or NaN. Instrument to find which layer. Repeat in BF16.

D. Compounding. Generate 1000 tokens with FP16 and BF16 from the same seed and prompt. Plot the token index at which they first diverge, over 50 prompts. Build a histogram.

E. Determinism cost. Enable full determinism and measure the throughput loss.

11. Interview questions#

1. Why can the same prompt produce different outputs on the same hardware?
2. Rank the sources of nondeterminism in an LLM serving stack.
3. Why is BF16 preferred over FP16 for deep transformers?
4. How would you validate a quantization change before deploying it?
5. Why is perplexity insufficient for validating a precision change?
6. What does it cost to make inference deterministic, and when is it worth it?
7. Why does prefix caching affect output reproducibility?

12. Further reading#

• [FUNDAMENTAL] Goldberg, “What Every Computer Scientist Should Know About Floating-Point Arithmetic”
• [ESTABLISHED] Micikevicius et al., “Mixed Precision Training”
• [REFERENCE] PyTorch reproducibility documentation
• Next: Section V — LLM Inference