Numerical Reproducibility in ML Systems

A reference implementation's answer is its arithmetic, not its formula. Seven numerical and ML packages were reimplemented in a new language from the references' own source and compared against scikit-learn, scipy, pandas, networkx, PySR and PyTorch along two execution paths, interpreted and compiled, against the same recorded fixtures, so that the check runs in CI with no reference runtime present. Where bit-exact parity was achievable, every remaining difference was traced to the operation that produced it; where it was not, the comparison used explicit numerical or behavioural criteria and documented the exceptions. The central case was verified again today against the checked-in scikit-learn fixture and the compiled kernel that wrote it: of 214 TF-IDF weights, a plain sum of squares reproduces 137, a fused multiply-add accumulated in the order the corpus first met each word reproduces all 214, and the wheel's l2-normalisation kernel compiles that sum to a single fused instruction.

Key findings

• The documented formula rarely specifies the exact arithmetic the reference executes. The matching answer depended on fused multiply-add contraction introduced by the reference's compiler, on the summation order a fitted row inherited from the corpus, on compensated and pairwise sums, on a variance kernel's restart, on stopping rules and call counts, and on tie-breaking at a split.
• Three standards of parity, stated per package: bit-exact for the frames, the retrieval, the graphs and the optimisers' points and counts; a tolerance for the boosting (1e-15), the GPT (1e-9 relative) and PySR's evaluation (1e-12); and behavioural for PySR's stochastic search, compared on what it recovers.
• Both the contraction and the order are real, and separable: against 214 checked-in weights, a plain sum misses 77 whichever order it uses; a fused sum in sorted vocabulary order still misses 47; a fused sum in first-seen order misses none. The mechanism is in the object code: the reference wheel's float64 normalisation kernel holds one fused multiply-add and no separate multiply or add.
• Holding a library to a reference is also a compiler test. The seven holdings found 38 unique defects underneath them, 30 in the compiler and associated tooling and 8 in the data-frame module, each fixed with a regression program that reproduces it on both paths.
• Bit-exactness is a diagnostic, not always the target: one golden-section pick lands a step off scipy's where the reference's own last bit is the artefact, and it is explained rather than hidden.
• Every holding passes at one commit: 45 checks across seven packages in 20 seconds and the 38 frame cases in 2 seconds, re-run on 2026-10-11 with no Python in the loop. A binary three days stale failed three of them: reference version, fixtures, compiler revision and execution binary must form one reproducible unit.

1. The claim

An ML system reproduces a reference when its numbers are the reference's numbers, and in practice that is tested bit for bit: a model ported between frameworks, a metric computed by two systems, a feature recomputed in a store and in a notebook. When such a test fails in the last digit, the usual verdict is floating-point noise and a tolerance is added. The seven holdings below refused the tolerance and asked, for every last-digit difference, what the reference actually computed. In every discrepancy investigated here the answer was not noise. It was an arithmetic choice the reference's documentation does not state, and often one its authors did not make, because the compiler made it for them.

2. The method

• The reference implementation is the specification. Its source defines the intended operations; its compiled artifact determines the arithmetic actually executed. Read the kernel itself, the Cython or the C, never its docstring or a summary of it, and where the last digit is at stake, read the object code it compiled to. One summary cost an hour and a wrong port.
• Record once, compare always. A script beside each package runs the reference and writes its answers to tracked fixtures, floats as the shortest text that reads back to the same double. The holding compares every number along both execution paths, interpreted and compiled, which share a lowering and a runtime and are therefore a differential check rather than two independent implementations, and it runs in CI with no reference runtime installed.
• Hold the counts too. Iterations, function calls, fitted-row order, tie nodes. An optimiser that reaches the same point by a different route is not the reference.
• Make the contraction explicit. A sum that must match a compiled reference is written with a fused multiply-add in the language, so the choice is visible in the source rather than left to a compiler flag.
• Explain to the bit or fix, or state the standard. Where bit-exact parity is achievable a difference is traced to the operation that produced it and either the port is corrected or the README states why the reference's digit is the artefact. Where it is not, because the reference draws random numbers or accumulates in an order a port cannot see, the package states a numerical tolerance or a behavioural criterion and documents the exceptions. A tolerance is never the first resort.

3. The seven holdings and their standards of parity

4. One finding end to end: the weights that needed the compiler's arithmetic

scikit-learn's TF-IDF rows are l2-normalised: each weight is divided by the square root of the row's sum of squares. The documentation says exactly that. The kernel that computes the norm is a Cython loop, sum += data[j] * data[j] over the row's stored values, compiled to C and then to machine code. In the wheel that wrote the fixture, scikit-learn 1.9.1 for Python 3.12 on arm64 macOS, disassembling the extension shows that the float64 specialisation of the l2 row-normalisation kernel carries one fused multiply-add instruction and no separate multiply or add: the compiler contracted the written multiply and add into one operation that rounds once where the source rounds twice. The float32 specialisation of the same kernel keeps a multiply and an add. The port, summing plainly, reproduced two thirds of the weights and missed the rest in the last digit.

The accumulation loop, as the disassembler prints it; the whole function and the command that produced it are in a tracked file:

 movi   d0, #0                 ; sum = 0
ldr    d1, [x12]              ; d1 = data[j]
fmadd  d0, d1, d1, d0         ; sum = d1 * d1 + sum, rounded once
add    x13, x13, #0x1         ; j += 1
add    x12, x12, x26          ; next element
cmp    x10, x13
b.ne   0x6246c                ; back to the load
fcmp   d0, #0.0               ; a row of zeros is left alone 

The order of the sum was the second half. In the tested scikit-learn 1.9.1 build and corpus, a fitted row's stored values are not in vocabulary order: the vectorizer assigns each word an index as it first meets it in the corpus, builds the rows in that order, and afterwards renumbers the indices to the sorted vocabulary without reordering the data, so the norm of a fitted row is accumulated in the order the corpus first met its words, while a row transformed later, a query, is built fresh in vocabulary order. That is an observation about this implementation and fixture, reproduced by reconstruction, not a stable guarantee of the vectorizer; neither fact is in the documentation, and both are in the source of the tested version.

The check that was re-run today takes the checked-in fixture written by scikit-learn 1.9.1 and reconstructs every weight four ways, with no scikit-learn present: a plain sequential sum and a fused sum, each in first-seen order and in sorted order. Each fused multiply-add was emulated using exact rational arithmetic for the multiplication and the addition, followed by one rounding to binary64, and the rounded result was then used as the accumulator for the next operation; that is what the instruction does, one rounding per operation, not one rounding for the whole sum. Separately, scikit-learn 1.9.1 itself was run live on the same corpus and reproduced all 214 weights of the fixture.

The cosines carry the same signature through a different library. scikit-learn's cosine of a query and a document is scipy's sparse matrix product, and the double-precision product kernel in the scipy 1.18.1 wheel likewise holds one fused multiply-add in its inner loop. Of 192 query-document pairs in the fixture, scipy's live product reproduces every cosine; a plain dot in vocabulary order differs in 3 of the 17 pairs that share two or more terms; a fused dot differs in none.

Two consequences follow. The language gained an explicit fused multiply-add so that a held sum can say in its source what the reference's compiler decided silently. And the reference's own reproducibility is build-shaped: whether a kernel is contracted depends on the compiler, its floating-point options and the target the wheel was built for, so the same scikit-learn version may answer differently in the last digit from a wheel built for another architecture or with other flags. That second point is a hypothesis this work did not test; the verification here is for these wheels on one platform, and it is stated as such.

5. What else the last digit knew

• pandas's rolling kernels are compiled from Cython too, and on this platform their update steps contract to fused multiply-adds; a port that writes the same arithmetic without the fused operation drifts in the last digit across a window. The rolling variance also restarts a window whose running sum of squares has collapsed below a tolerance, a numerical guard that is semantics: a port without it is correct on paper and wrong on data.
• Sums are not one thing. numpy's plain sum is pairwise; pandas's group and window totals are compensated; a sequential sum matches neither. The frame module keeps each where pandas keeps it.
• Stopping rules are part of the model. Matching scipy's iteration and function-call counts, not only its points, provided evidence that the ports reproduced the reference algorithms' stopping behaviour rather than merely converging to similar answers; the same cases held the language's automatic differentiation to scipy's hand-written gradients.
• Tie-breaking decides real predictions. At tied splits scikit-learn's choice is a random feature order; a port that takes the first feature agrees on every leaf assignment and disagrees on importances, and the README says where.
• Determinism has an edge. PySR's deterministic parts, the evaluation and the scoring, match to 1e-12; its search is stochastic and is compared on what it recovers rather than on which equation it prints.
• The reference can be the artefact. In one golden-section case the two last points' values differ in the last bit of the exponential, and scipy picks the other; the port's answer is one step away and the README explains why it is left so.

6. Why it generalises

7. What it found underneath

A library held to a reference exercises its language harder than any unit test. The seven holdings found 38 unique defects beneath them; a defect found by several cases is counted once, and each has its own regression program in the parity suite that runs on both paths.

Three of them, each as the sequence a reviewer would want: the input that exposed it, the wrong behaviour, the root cause, the test that holds it.

• A function that was never emitted. A retrieval function was called only inside a collecting loop over a frame's column. The compiled program named a function it did not hold and did not build, while the interpreted path ran it. The lowering's scan had typed the function unknown because its only call site was in a loop body over a column, so the emitter skipped it. Fixed in the lowering; parity program 191 calls such a function on both paths.
• Sums a digit off. A group total and a column sum compared with pandas and numpy differed in the last digit. The runtime summed sequentially where pandas keeps compensated totals in groups and windows and numpy sums pairwise. Fixed by keeping each kind of sum where the reference keeps it; the 38 frame cases hold them.
• A cell read two ways. A timestamp column's cell, held in a variable and read back, came out as seconds on the compiled path and as text on the interpreted one. The lowering had forgotten the kind of a column's cells when the column was bound to a variable, at four separate sites. Fixed by keeping a column's cell kind through every binding; parity program 189 reads such a cell both ways and compares.

8. Cost and operation

The fixtures are small, and they are tracked: the retrieval reference writes 24 documents and 8 queries and the frame reference draws four tables from a seeded generator, so a re-run rewrites the same bytes. The holdings run inside the language's test gate with no Python present: 45 checks across the seven packages in 19.8 seconds and the 38 frame cases in 2.1 seconds on the installed binary of 2026-10-11. The one operational seam is version drift between the binary and the packages: a release binary three days old failed three of the checks today, two because the compiler had moved on and one because a test had been re-seeded to a new random stream, and the installed binary of the same morning passed all 45. A fixture pinned to a reference version is the point; a binary that is not pinned to its sources is the trap. The lesson generalises to any ML platform: reference version, fixtures, compiler revision and execution binary must form one reproducible unit, and a parity check that pins only the first two is pinned to nothing.

Toward a preprint

This is an engineering case study, not yet a research paper. The research question the holdings raise is a sharper one: how much of cross-language numerical disagreement is attributable to arithmetic contraction, to operation order, and to reference-specific stopping behaviour. The TF-IDF reconstruction is the candidate central experiment, since it isolates two of the three; the other packages are case studies for breadth. A preprint needs a reproducible artifact, the fixtures with their reference versions, the reconstruction code, minimal floating-point examples and the exact environment, which is the state this section's manifest moves toward.

Appendix: manifest

Scope

Seven packages, one language, one platform, one author, and the references' own source as the specification. What this documents is the arithmetic the references compute and do not state, the method that found it, the standard each package is held to, and the defects the method turned up beneath the packages; the central finding was reconstructed today from the fixture alone, confirmed in the reference wheel's object code, and the whole set re-run at one commit. What it does not claim is that every package is a bit-exact port, that every library should be held to the last digit, or that these findings hold for other builds. The ports were written by AI coding agents under the method above and the holdings were theirs to pass; the method, the refusal of tolerances and the decision to explain each remaining digit were mine.