Understanding and Mitigating Numerical Sources of Nondeterminism in LLM Inference

In plain words

A language model's answer can change when nothing about the question changes. If you run the same model with the same prompt, the same seed and greedy decoding (always pick the most likely next word), the answer can still differ depending on how many GPUs you use, which GPU model you use, or how many questions the server processes at once. The cause is arithmetic: computers round numbers, and the order in which they add a long list of rounded numbers changes the final total slightly. A tiny difference in one early step can push a long reasoning chain onto a different path, so the final score moves too.

The authors show that this is common, and that it hits long reasoning answers hardest. Using a number format with more precision (32-bit instead of 16-bit) removes most of the problem, but doubles the memory for the weights. Their proposed LayerCast keeps the weights in 16-bit to save memory, and converts each weight back to 32-bit only at the moment it is used, so the arithmetic itself is done at the higher precision.

TL;DR

• Under BF16 with greedy decoding, the reported numbers depend on GPU count, GPU type and batch size, both in accuracy and in response length §3.2.
• FP32 makes greedy runs nearly identical across these settings §3.2. FP16 is in between §3.2.
• LayerCast stores linear-layer weights in BF16 and runs the computation in FP32, aiming for FP32-level stability at lower memory use §4.
• The numbers are author-reported on four 7–8B models and five benchmarks. The code has a broken import and some undocumented differences from the paper (see Pitfalls).

Problem

Benchmark scores for LLMs are assumed to be reproducible. Two common evaluation styles are greedy decoding with one run, and random sampling averaged over several runs §1 §2.1. The authors argue that greedy decoding is not deterministic across hardware and system configurations in practice, because of rounding in floating-point arithmetic §1 §2.2. Random sampling is also affected: the spread across system configurations is added to the spread from sampling itself §3.3. The problem is worst for reasoning models, whose long outputs let early rounding differences accumulate §1 §3.2. Common fixes do not cover this. Fixed seeds do not ensure determinism under greedy decoding §2.1, and batch-invariant kernels only address changes in the batch dimension, not changes in tensor-parallel size or GPU type §2.2.

Method

The paper has two parts: an analysis of how number format and runtime configuration change outputs §2–§3, and an inference pipeline, LayerCast §4.

Number formats. BF16, FP16 and FP32 differ in mantissa bits. The authors state FP16 has 10 and BF16 has 7 bits, with FP32 at 23 §3.2 fig. 2. A small example shows how the same value rounds differently in each format tab. 1. Floating-point addition is not associative, so the order of a sum changes its rounding §2.2 tab. 2. In a GPU reduction, the order depends on thread scheduling §2.2.

Sources of order change in serving. The authors list continuous batching, the choice of matmul kernel (for example split-K or not), kernel block sizes, collective operations such as AllReduce, and tensor-parallel layout as sources of changed summation order §2.2.

Analysis protocol. For each model and benchmark, the authors run 12 configurations: 2 GPU types (A100, L40S) × 2 GPU counts (2, 4) × 3 batch sizes (8, 16, 32) §3.1. Greedy metrics are computed over these 12 runs:

• Std@Acc: sample standard deviation of accuracy across the 12 runs §3.1.
• Avg_Std@Output_Length: per problem, standard deviation of output length across the 12 runs, averaged over the dataset §3.1.
• Div_Index: the first token position where runs that agree so far produce a different token, averaged over problems §3.1 app. C.
• Avg_Std@top1_prob: for positions before divergence, the standard deviation of the top-1 token probability across runs, averaged over positions and problems §3.2.

For random sampling, the paper reports Pass@1 (mean accuracy over runs) and the spread across six setups that vary batch size and GPU count. The spread is not taken over repeated sampling runs of a single setup §3.3.

Ablations. The authors vary batch size, GPU count and GPU type one at a time on DeepSeek-R1-Distill-Qwen-7B and report Avg_Std@top1_prob §3.4 app. D.

LayerCast. The procedure in words §4:

1. Load all parameters in FP32.
2. Cast every linear layer's weights and biases to BF16 for storage.
3. Before each linear layer runs, upcast that layer's weight back to FP32, one layer at a time, and compute the matmul in FP32, so that all computation happens in FP32 and only the stored weights are 16-bit §4 fig. 7.

The authors describe the memory saving as mainly for the KV cache in long-context use §4. The release is a vLLM patch §1.

# Load every parameter in FP32  [§4]
# [code: evaluation/eval_layercast.py#L159-L163]
model = load_model(checkpoint, dtype=float32)

# Store each projection weight and bias in BF16  [§4]
# [code: evaluation/patch_layercast.py#L48-L54]
# (the code matches modules whose name contains "proj")
for module in model.projection_modules():
    module.weight = to_bfloat16(module.weight)
    if module.bias is not None:
        module.bias = to_bfloat16(module.bias)

# Before a linear layer runs, make an FP32 copy of its weight  [§4]
# [code: evaluation/patch_layercast.py#L92-L102]
def linear_forward(layer, inputs):
    # Inputs are already FP32  [code: evaluation/patch_layercast.py#L96-L97]
    weight = to_float32(layer.weight)     # temporary copy for this call
    bias = to_float32(layer.bias)         # None if the layer has no bias
    return linear(inputs, weight, bias)   # output stays FP32

Equations: the paper gives the non-associativity of a reduction operator in display form, without equation numbers §2.2. There are no other equations.

Hyperparameters

Results

Author-reported. Setup: 4 models, 5 benchmarks, 12 runtime configurations per model and benchmark, vLLM backend, greedy decoding unless noted §3.1. The FP32 and BF16 rows are the paper's baselines; there is no separate external baseline.

Greedy decoding: accuracy standard deviation (Std@Acc, %) across 12 runs tab. 3

Greedy decoding: standard deviation of output length (tokens) across 12 runs tab. 4

Random sampling: standard deviation of Pass@1 (%) across 6 configurations tab. 5

Headline numbers from the text

LayerCast vs BF16 vs FP32, accuracy std across 12 configurations, DeepSeek-R1-Distill-Qwen-7B app. F. Values are fractions here, not percents.

Other reported results

• Qwen3-32B, AIME'24, 9 FP32 and 12 BF16 configurations: BF16 Std@Acc 5.02%, FP32 1.11% (text) app. C.1.
• DeepSeek-R1-Distill-Llama-8B, GPQA Diamond, 12 configurations: BF16 Std@Acc 3.03%, FP32 1.96% (text) app. C.1.
• HF Transformers backend, AIME'24, 12 configurations: BF16 Std@Acc 2.23% (Qwen2.5-7B-Instruct) and 1.51% (Llama-3.1-8B-Instruct); FP32 0% for both app. C.2.
• LayerCast on DeepSeek-R1-Distill-Llama-8B, AIME'24: 36.67% in all 12 configurations, identical to FP32; BF16 ranges from 23.33% to 40.00% app. F.

Rounding example tab. 1: the value 1.00012 is stored as about 1.00012004 in FP32 (error ≈ +4.4e-8) and as 1.0 in FP16 and BF16 (error ≈ −0.00012).

Paper vs code

Code: https://github.com/nanomaoli/llm_reproducibility at code_commit 09ed20f0 (HEAD on the clone, dated 2026-05-29). The paper links this repository in its abstract. The repository's latest README also describes a later method (TBIK, arXiv:2511.17826) and later kernels. I did not check whether the evaluation files at this commit are the ones used for the paper's numbers.

• The LayerCast script cannot run as shipped. evaluation/eval_layercast.py imports patch_vllm eval_layercast.py, but the repo has no patch_vllm module. The LayerCast code is in evaluation/patch_layercast.py patch_layercast.py. The README tells the user to run eval_layercast.py README.md.
• Baseline and LayerCast differ in attention backend. The patch sets VLLM_ATTENTION_BACKEND to XFORMERS patch_layercast.py patch_layercast.py. The baseline script eval_main.py has no patch import (a grep of the repo finds none) and makes its own vLLM call eval_main.py, so it uses vLLM's default backend. The paper does not mention an attention backend §3.1. Comparisons between BF16 or FP32 baselines and LayerCast therefore change the backend too. The effect is not measured here.
• KV cache is FP32 in the code. The patch sets the cache dtype to FP32 patch_layercast.py and sizes cache blocks as FP32 patch_layercast.py. The paper says LayerCast helps mainly the KV cache for long contexts §4. The code does not reduce KV-cache memory. The 34% memory saving §4 therefore cannot come from the KV cache in this code, and the source of the figure is not stated.
• Not every linear layer is cast to BF16. The paper says all linear weights and biases are cast §4. The conversion matches module names that contain proj patch_layercast.py. The LM head is converted by calling the same function on the head module itself patch_layercast.py. Its own module name has no proj, so as I read it, the head weight stays FP32 (inferred from code, not run).
• Logits are scaled by 1.2. Both patched constructors create the logits processor with scale=1.2 patch_layercast.py patch_layercast.py. The paper does not mention this. A positive scale does not change which token is the greedy argmax, but it does change the logprobs saved by the script and any top-1 probability statistics computed from them. I did not check this against vLLM's LogitsProcessor.
• Eager mode. The LayerCast run sets enforce_eager=True eval_layercast.py. The paper does not say whether the baselines used the same setting. The code for eval_main.py was not compared on this point.
• Model coverage of the patch. The patch replaces forward methods for the Qwen2 and Llama classes only patch_layercast.py patch_layercast.py. The paper's LayerCast results cover DeepSeek-R1-Distill-Qwen-7B and DeepSeek-R1-Distill-Llama-8B app. F. Qwen3-32B is not covered by this code as written, and LayerCast results are not given for it.
• Sampling settings match. Temperature 0.7 and top-p 0.95 for sampling eval_passk.py match §3.1. Greedy uses temperature 0 eval_main.py and matches §3.1.
• vLLM version is inconsistent. The README says the experiments used vLLM 0.8.2 README.md, but its install line pins vLLM 0.11.0 README.md. The patch imports vllm.worker.cache_engine and vllm.attention patch_layercast.py patch_layercast.py. I did not check whether these modules exist in either version.
• Seed. The paper does not state a seed value §1. The code default is 42 eval_main.py.
• Prompt templates were not compared line by line with the code.
• "Few lines of code change" §1 does not match the code: the LayerCast module is about 310 lines and replaces cache allocation, rotary embedding, linear application and attention forward patch_layercast.py.

Pitfalls

• Inconsistent numbers in the paper.
  • Text says only 2.2% of FP32 MATH500 examples diverge for DeepSeek-R1-Distill-Qwen-7B §3.2. The table for the same model and benchmark gives 5.80% app. C.
  • In the random-sampling appendix, the BF16 MATH500 row for Qwen2.5-7B-Instruct repeats the DeepSeek-R1-Distill-Qwen-7B values (90.65, 91.15, …) and its stdev 0.3158 app. E. The main table gives 0.4663 for that cell tab. 5. Treat the Qwen2.5 BF16 MATH500 random-sampling figure as unresolved.
• "Larger batch sizes consistently lead to lower Avg_Std@top1_prob" app. D is not borne out by the appendix table for every row. For example, DeepSeek-R1-Distill-Qwen-7B BF16 LCB-Easy is 56.00 (BS=8), 58.99 (BS=16), 51.39 (BS=32), and Llama-3.1-8B-Instruct BF16 AIME'24 is 53.86, 58.55, 47.23 app. D. Use the claim as a trend, not a rule.
• Units differ between tables. Most tables report percentages. The LayerCast table reports fractions (0.4333 = 43.33%) app. F.
• Small-sample statistics. Std@Acc uses 12 configurations, and the random-sampling stdev uses 6 §3.1 §3.3. The AIME'24 set has 30 problems, so one answer is 3.33 percentage points §3.3. The paper accepts this noise for the n=64 FP32 exception §3.3.
• Table numbers here come from LaTeX source order. I counted tables and figures in the order the source defines them. Check against the rendered PDF before citing a number to a reader.
• Zero-like values. Entries such as 5.8e-17 and 1.45e-17 are floating-point noise for zero tab. 3. Read them as zero.
• LayerCast cost claims lack data. "Doubles memory and inference time" for FP32 versus BF16 and the 34% saving are stated in §4 without a table or method. Divergence "below 3.4%" appears only in §4 text and a figure.
• Greedy outputs need the same code path. The paper's claim is about the effect of hardware and batch configuration under BF16. The code shows that the LayerCast run also changes the attention backend and eager mode (see Paper vs code). Do not attribute a LayerCast-vs-BF16 difference to the number format alone without controlling these.
• Accuracy is not the only effect. Response length varies by thousands of tokens under BF16 for reasoning models tab. 4. Efficient-reasoning evaluations that compare output lengths are affected §3.2.
• Mixed-precision memory. LayerCast's per-call FP32 weight copy is temporary. The paper does not measure its peak memory §4.
• Prompts and datasets come from SkyThought-evals with stated templates app. A. Changing the template changes the output set, so reproducing the numbers needs the same templates.

Lineage

• Builds on: none stated.
• Related work the paper names: batch-invariant kernels for FlexAttention, RMSNorm and MatMul (cited to a blog post, not an arXiv paper) §2.2; PyTorch deterministic-algorithm flags §2.1; FP32 inference practice §5.
• Compared against: the same models run in BF16, FP16 and FP32 under the same configurations §3.2 tab. 3.
• Follow-up named in the code README: TBIK, arXiv:2511.17826, for tensor-parallel invariance README.md. Not reviewed here, and no card exists yet.