Where identical code gives different answers: a reproducibility budget for a protein-sequence pipeline

Abstract

We ran one small protein-sequence classification pipeline 192 times (153 scikit-learn runs, 39 PyTorch runs) on one Apple-silicon machine, changing one thing at a time: thread count, BLAS backend, CPU architecture (native arm64 against x86_64 under Rosetta), Python version, library version, input precision and random seed. We hashed the output of every stage and measured how far each difference travels, from bitwise identity to changed test predictions to changed conclusions. Floating-point noise near 1e-13 was invisible downstream for logistic regression and a neural network, but changed up to 4 of 1,998 random-forest predictions. One scikit-learn release changed 42 gradient-boosting predictions, about three quarters of what a different seed changes, and flipped a significance verdict. Seeds dominated everything else.

Motivation

Two of the nine software stacks we tested disagree about whether gradient boosting beats a small neural network on the same 1,998 held-out proteins. The code is identical, the seed is identical, the data are identical. In five stacks the paired bootstrap interval for the accuracy difference excludes zero. In the other four it does not. The only thing that separates the two groups is the version of scikit-learn.

Reproducibility advice in computational biology is mostly binary: pin the environment, record the seed, ship a container. That advice is right, but it gives a practitioner no way to decide which differences matter. A group that cannot pin everything (a new laptop, a cluster upgrade, a collaborator on a different operating system) needs a budget: for each stage of the pipeline and each thing that might change, how big is the difference, and does it reach the decision the pipeline exists to support?

We built a deliberately ordinary pipeline, because ordinary is what most people run. Proteins go in as sequences, become 3-mer count vectors, are weighted by TF-IDF, compressed by truncated SVD, clustered with k-means, searched for nearest neighbours, and used to train four classifiers that predict whether a protein has an annotated transmembrane segment. Every stage writes an array, and we hash every array.

Related work

Variance in deep-learning training has been measured carefully. Pham and colleagues attributed accuracy differences between identical training runs to algorithmic and implementation-level factors and found that most surveyed practitioners were unaware of the latter [1]. Zhuang and colleagues characterised how tooling choices contribute to non-determinism across hardware and showed that aggregate accuracy can hide large per-class sensitivity [2]. Summers and Dinneen traced run-to-run variability to training instability, so that even a one-bit change in initial weights produces models that diverge [3]. In neuroimaging, the operating system changed the output of commonly used analyses [4], and FreeSurfer version changes moved volume estimates by 8.8 percent on average between version 5.0.0 and earlier releases [5]. In single-cell analysis, Rich and colleagues quantified how package choice and version change results across a whole workflow [6].

What is new here is narrower than the topic. We study a classical protein-sequence pipeline rather than deep training, and we isolate factors in pairs of runs that differ in exactly one thing, including library versions at a fixed numpy and BLAS (scikit-learn 1.4.2, 1.5.2 and 1.9.1 on numpy 1.26.4) and BLAS backends at a fixed scikit-learn. We follow differences stage by stage and report them at three levels: bitwise identity, magnitude, and test-set decisions. We use the random seed as a yardstick, so that every effect is expressed in units of the noise the field already accepts. And we measure how each model class amplifies or damps the same small input perturbation.

Method

Data

We downloaded reviewed Swiss-Prot entries from the UniProt REST API (release 2026_03, accessed 2026-10-05) [7] for twelve organisms (human, mouse, budding yeast, fission yeast, Escherichia coli K-12, Bacillus subtilis, Arabidopsis, Drosophila, Caenorhabditis elegans, Mycobacterium tuberculosis, Dictyostelium and zebrafish), restricted to sequences of 60 to 700 residues made only of the twenty standard amino acids. We shuffled with a fixed seed, kept the first 6,000 proteins (mean length 337), and split them 4,002 for training and 1,998 for testing without any homology filtering. The label is whether UniProt records a transmembrane feature (23.5 percent of proteins; 468 of the 1,998 test proteins). The label is an annotation and not ground truth, which does not matter for this study because we compare runs with each other, not with biology.

Pipeline

The stages and their parameters were fixed in advance and never tuned.

1. 3-mer counts: an exact integer computation into a sparse 6,000 by 8,000 matrix.
2. TF-IDF weighting with L2 normalisation, in float64 or float32.
3. Truncated SVD to 64 components, once with the randomized algorithm (5 power iterations, the library default) [8] and once with ARPACK.
4. k-means with 30 clusters and 10 initialisations on the randomized embedding.
5. The 10 nearest neighbours of every protein by cosine similarity, with stable tie handling.
6. Four classifiers trained on the 64-dimensional embedding: logistic regression, a one-hidden-layer MLP (64 units, Adam, 60 epochs), histogram gradient boosting (100 iterations) and a random forest (200 trees).
7. A separate PyTorch MLP (64-256-256-1, 40 epochs, batch 64) trained on the reference embedding, so that PyTorch differences are not confounded with upstream differences.

Software stacks

Nine stacks, all on one Apple M4 machine with 16 GB of memory that was shared with other jobs throughout. Four use numpy 1.26.4 with its bundled OpenBLAS; five use numpy 2.x, which on this machine links Apple Accelerate. Two run as x86_64 binaries under Rosetta 2.

Each stack ran with 1, 2, 3, 4, 6, 8 and 10 threads (environment variables for OpenMP, OpenBLAS and Accelerate set together), each twice, plus float32 input at 1 and 4 threads. The reference stack also ran with seeds 1 to 9. All other runs use seed 0.

Factors and pairs

A factor effect is the difference between two runs that differ only in that factor: threads (same stack, 1 against n threads), repeat (same stack, same threads, run twice), BLAS backend (OpenBLAS against Accelerate at a fixed scikit-learn: 1.5.2 and 1.9.1, with numpy and scipy moving together because the numpy major version carries the BLAS), architecture (arm64 against Rosetta, at fixed versions: two pairs), Python (3.12.10 against 3.14.2), scikit-learn version at fixed numpy 1.26.4 (three pairs), precision (float64 against float32, nine pairs) and seed (reference seed 0 against seeds 1 to 9). In total 153 scikit-learn runs, 63 repeat pairs, 54 thread pairs, 2 BLAS pairs, 2 architecture pairs, 1 Python pair, 3 version pairs, 9 precision pairs and 9 seed pairs.

Measures

For every stage we store a SHA-256 prefix of the output bytes, so bitwise identity is a yes or no. For magnitude we use the relative Frobenius norm of the embedding difference after aligning column signs, the largest relative singular-value difference, the adjusted Rand index for k-means labels, the mean Jaccard index of nearest-neighbour sets, the largest absolute difference in predicted probability, and the number of test proteins whose predicted class changes at a threshold of 0.5 ("flips", out of 1,998). For decisions we ran a paired bootstrap on the test set (2,000 resamples with fixed indices, 95 percent percentile interval) for every pair of the four models, on accuracy and AUC, in every run.

Results

When nothing changes, nothing changes (almost)

Sixty-three pairs of runs repeated the same stack at the same thread count. Every array was bitwise identical, with one exception: the k-means inertia differed in 5 of 63 pairs, by at most 1.9e-16 relative. The cluster labels never differed. This is the null experiment for the whole audit: the hashing detects differences only where there are differences.

Threads, BLAS, architecture and Python: noise that does not travel

The table below summarises each factor. Bitwise identity is lost quickly and cheaply: any change to numpy's numerical backend changes the SVD embedding in its last digits. It does not change much else.

With OpenBLAS, every thread count other than one changed the embedding bits (24 of 24 pairs), by at most 1.6e-13 relative, and every downstream consumer except one was unmoved: k-means labels, neighbour lists, logistic-regression probabilities (to 1e-14), MLP probabilities (to 5.6e-14), boosting and forest predictions. With Accelerate, 29 of 30 thread pairs were bitwise identical. We do not read that as a finding about Accelerate: CPU time in a matrix-multiplication probe did not grow with the thread setting under Accelerate, as it did under OpenBLAS, so we cannot show that the setting changed anything.

The one downstream consumer that noticed was the random forest. Between native arm64 and x86_64 under Rosetta, with the same scikit-learn and numpy, the embedding differed by 8e-14 to 1.2e-13 and the forest changed 1 prediction in one pair and 4 in the other (largest probability difference 0.02, AUC change up to 8.4e-4). A forest chooses split thresholds from the data, and a 1e-13 shift can reorder near-ties. Gradient boosting, which bins features, did not notice noise of this size. The k-means labels and the neighbour lists were identical in all of these pairs.

Python 3.12.10 against 3.14.2 produced bitwise identical output at every stage.

Library versions: the largest effect that is not randomness

With numpy fixed at 1.26.4, going from scikit-learn 1.5.2 to 1.9.1 changed 42 of the 1,998 gradient-boosting predictions (largest probability difference 0.446), and nothing else. Embeddings, singular values, k-means labels, neighbour lists, logistic regression, MLP and forest outputs were bitwise identical. The scikit-learn 1.9 release notes list a fix to how histogram gradient boosting computes its bin edges [9], which is consistent with this pattern, although we did not trace the difference to that change. Boosting accuracy fell from 0.8984 to 0.8914, a net loss of 14 correct predictions.

Going from 1.4.2 to 1.5.2 changed 3 MLP predictions (largest probability difference 0.062), 2 boosting predictions and 12 forest predictions (largest probability difference 0.155), and left logistic regression untouched. We did not trace these. We also saw something that matters for anyone comparing raw embeddings: between 1.4.2 and 1.5.2, 24 of the 64 randomized-SVD columns and 26 of the 64 ARPACK columns changed sign. A downstream logistic regression does not care, but a script that diffs embeddings, or a plot that shows component 7, will.

Precision

Float32 input changed the embedding by up to 7.4e-5 relative, left logistic-regression and MLP predictions unchanged (MLP probabilities moved by at most 0.0022), and changed 43 boosting predictions in the median pair (maximum 47) and 7 forest predictions (maximum 8). Neighbour lists stayed identical for at least 99.57 percent of proteins, and k-means labels were identical in all nine pairs. Boosting is again the model that notices.

Seeds are the yardstick

A different seed changed the randomized-SVD embedding by a relative 0.63 to 0.73 after sign alignment, moved the singular values by up to 0.9 percent, and flipped the sign of 15 to 27 of the 64 columns. The ARPACK embedding, whose random starting vector does not survive convergence, moved by 2.4e-13. The seed effect on the embedding therefore comes from the randomized algorithm not having converged with 5 power iterations, and it propagates everywhere: k-means labels fell to an adjusted Rand index of 0.50 on average (range 0.40 to 0.54), the Jaccard index of 10-nearest-neighbour sets was 0.29, and test predictions flipped for a median of 13 proteins (logistic regression), 49 (MLP), 58 (boosting) and 21 (forest). Test accuracy across the ten seeds had a standard deviation of 0.0015 for logistic regression, 0.0027 for the MLP, 0.0028 for boosting and 0.0018 for the forest.

Against this yardstick, the budget reads as follows. Backend, threads and architecture move the embedding about twelve orders of magnitude less than a seed (1e-13 against 0.69) and change no prediction of logistic regression, the MLP or boosting; only the forest moved, by at most 4 predictions. Float32 input moves the embedding four orders of magnitude less than a seed, and the scikit-learn 1.9 release does not move it at all, yet each changes 42 to 43 boosting predictions, about three quarters of the seed effect (58).

Neural network on CPU and GPU

The PyTorch MLP gave bitwise identical probabilities for all 24 CPU runs (1, 2, 4 and 8 threads, deterministic algorithms on and off, three repeats each). The six runs on the Apple GPU (Metal Performance Shaders) were bitwise identical to each other, with and without the deterministic flag, and differed from the CPU: 5 of 1,998 predictions flipped, the largest probability difference was 0.135, and AUC changed by 3.0e-5. Nine other seeds on the CPU flipped between 85 and 107 predictions (largest probability difference above 0.99) and had an AUC standard deviation of 0.0026. The hardware path was therefore about twenty times smaller than the seed in flips. The network is small enough that its matrix products may not be split across threads at all, so identical CPU results do not generalise to larger models.

Do the conclusions change?

We compared the four models pairwise in every run, on accuracy and on AUC: twelve verdicts per run, each defined as whether the 95 percent paired bootstrap interval excludes zero. Eleven of the twelve were the same in all nine stacks. For example, the MLP is more accurate than logistic regression by about 1.4 percentage points in every stack, and the forest has lower AUC than each of the other three models in every stack. Across the ten seeds at the reference stack, nine of the twelve verdicts were the same in all ten runs; the other three (MLP against logistic regression on AUC, forest against MLP on accuracy, boosting against forest on accuracy) were significant in only 2, 8 and 9 of the ten seeds respectively.

The twelfth verdict is the one from the opening paragraph. Boosting minus MLP accuracy was 0.93 to 0.97 percentage points, with a bootstrap interval excluding zero (lower end 0.05 to 0.10 points), in the five stacks with scikit-learn 1.4.2, 1.5.2 or 1.6.1. In the four stacks with 1.9.1 it was 0.27 percentage points and the interval ran from minus 0.70 to plus 1.25. Across ten seeds at the reference stack the same difference ranged from minus 0.29 to plus 0.65 percentage points, changed sign, and was significant in none. The bootstrap ignores sequence homology among test proteins and is therefore too narrow; the companion paper in this series measures by how much. That does not change the point, which is that this verdict depends on a version number and on a seed.

The cheapest reduction in the budget

Logistic regression is deterministic given its inputs, so the 13 flips it shows across seeds come entirely from the embedding. We asked how fast they shrink with more power iterations in the randomized SVD (six seeds per setting, reference stack, flips counted against seed 0):

Fifty iterations cost about six and a half times the default and cut the logistic-regression flips from 13.6 to 1.0, while raising the neighbour-set agreement from 0.29 to 0.81. The ARPACK solver, which is stable to 2.4e-13 across seeds, would remove the seed from the embedding entirely; we did not run the downstream models on it.

Falsification and limits

We tried to break the audit in four ways and report what survived.

The null held: repeating identical runs gave bitwise identical outputs apart from a 1e-16 effect in the k-means inertia. The hash-based measure is therefore not producing false differences.

The attribution of version effects is partial. We isolated scikit-learn at fixed numpy, but the old and new stacks also differ in scipy (1.11.4, 1.13.1), and the backend pairs move numpy and scipy together. We name the library whose version string changed and we do not claim to have identified the line of code responsible. The bin-edge change in the 1.9 notes is a candidate, not a finding.

The hardware is one machine. "Architecture" here is x86_64 code executed by Rosetta translation, not Intel or AMD silicon, "GPU" is one Apple GPU, and we had no Linux, no MKL and no CUDA. The Accelerate thread setting may not have been honoured, as noted above. Operating-system and compiler effects, which the neuroimaging literature found large [4, 5], are untested here because the operating system was constant.

The pipeline is small. Matrices of 6,000 by 64 may fall below the size where BLAS libraries split work across threads, so the thread-count results are probably an underestimate of what a larger pipeline shows. We would not extrapolate the 1e-13 scale to a model with millions of parameters.

There are only ten seeds and nine stacks. The decision analysis rests on a single test set of 1,998 proteins, flips are counted at one threshold, and the bootstrap ignores homology. The split of five stacks against four in the one verdict that flipped is a count of stacks, and the stacks are not independent draws from anything. The percentile intervals for that verdict barely exclude zero in the old stacks, so it is a borderline verdict that a library release pushed across the line, not a large effect reversed.

What this means in practice

If you cannot pin everything, pin in this order of importance: the random seed and, within it, the convergence of any randomized solver; the versions of tree-based learners, which can change under a minor release; the numeric precision of the input. Thread count, BLAS backend, Python version and CPU architecture came in far lower. The one visible effect was that a random forest changed up to about one prediction in 500 when its inputs moved by 1e-13.

If you compare two models, log the verdict under a second seed and a second library version before you report it. Eleven of our twelve verdicts were stable across stacks; the twelfth, a 0.3 to 1 percentage-point accuracy difference on 1,998 proteins, was not. A margin of that size on a test set of that size is inside the budget.

If you publish an embedding, publish it with a sign convention. Component signs flipped for 24 of 64 columns between two library releases and for 15 to 27 columns between seeds.

If you need an inexpensive win, replace a randomized SVD with a converged solver, or raise its iteration count, and rerun the ten-seed comparison. In our pipeline that was the largest reduction in seed-driven change that we found.

Reproducibility appendix

Seeds and settings. Data: Swiss-Prot 2026_03, 12 organisms, 1,000 proteins sampled per organism with seed 20261005, first 6,000 of a shuffled pool used. Pipeline seed 0 for all stack comparisons, seeds 0 to 9 on the reference stack. Bootstrap: 2,000 resamples, rng seed 12345. Run times: median 11.9 s per scikit-learn run (range 4.5 to 342 s on a loaded machine), 55 minutes of summed wall time for 153 runs; PyTorch runs 2 to 36 s. Each factor-pair statistic in this paper is computed from the stored hashes and arrays of the runs listed above.

The core of the pipeline:

import numpy as np, scipy.sparse as sp, hashlib
from sklearn.feature_extraction.text import TfidfTransformer
from sklearn.decomposition import TruncatedSVD
from sklearn.cluster import KMeans
from sklearn.linear_model import LogisticRegression
from sklearn.neural_network import MLPClassifier
from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier

AA = 'ACDEFGHIKLMNPQRSTVWY'
LUT = np.full(256, -1, dtype=np.int64)
for i, c in enumerate(AA):
    LUT[ord(c)] = i

def kmer_counts(seqs, k=3):
    rows, cols, vals = [], [], []
    for r, s in enumerate(seqs):
        x = LUT[np.frombuffer(s.encode(), dtype=np.uint8)]
        code = np.zeros(len(x) - k + 1, dtype=np.int64)
        for j in range(k):
            code = code * 20 + x[j:len(x) - k + 1 + j]
        u, c = np.unique(code, return_counts=True)
        rows.append(np.full(len(u), r)); cols.append(u); vals.append(c)
    return sp.csr_matrix((np.concatenate(vals).astype(np.float64),
                          (np.concatenate(rows), np.concatenate(cols))),
                         shape=(len(seqs), 20 ** k))

def digest(a):
    return hashlib.sha256(np.ascontiguousarray(a).tobytes()).hexdigest()[:16]

def run(seqs, y, ntr, seed, dtype):
    X = TfidfTransformer().fit_transform(kmer_counts(seqs)).astype(dtype)
    E = TruncatedSVD(64, algorithm='randomized', n_iter=5, random_state=seed).fit_transform(X)
    km = KMeans(30, n_init=10, random_state=seed).fit(E)
    En = E / np.linalg.norm(E, axis=1, keepdims=True)
    S = En @ En.T
    np.fill_diagonal(S, -2)
    nn = np.argsort(-S, axis=1, kind='stable')[:, :10]
    tr, te = slice(0, ntr), slice(ntr, None)
    models = {
        'lr': LogisticRegression(max_iter=1000),
        'mlp': MLPClassifier(hidden_layer_sizes=(64,), max_iter=60, random_state=seed),
        'hgb': HistGradientBoostingClassifier(random_state=seed, max_iter=100),
        'rf': RandomForestClassifier(n_estimators=200, random_state=seed),
    }
    probs = {k: m.fit(E[tr], y[:ntr]).predict_proba(E[te])[:, 1] for k, m in models.items()}
    return {'emb': E, 'km': km.labels_, 'nn': nn, **{'p_' + k: v for k, v in probs.items()}}

The comparison of two runs:

from sklearn.metrics import adjusted_rand_score

def compare(a, b):
    out = {k: digest(a[k]) == digest(b[k]) for k in a}
    sgn = lambda E: np.sign(E[np.argmax(np.abs(E), axis=0), np.arange(E.shape[1])])
    Eb = b['emb'] * sgn(b['emb']) * sgn(a['emb'])
    out['emb_rel'] = np.linalg.norm(a['emb'] - Eb) / np.linalg.norm(a['emb'])
    out['km_ari'] = adjusted_rand_score(a['km'], b['km'])
    for k in ('lr', 'mlp', 'hgb', 'rf'):
        pa, pb = a['p_' + k], b['p_' + k]
        out['flips_' + k] = int(((pa > 0.5) != (pb > 0.5)).sum())
        out['maxdp_' + k] = float(np.max(np.abs(pa - pb)))
    return out

Each stack ran under OMP_NUM_THREADS, OPENBLAS_NUM_THREADS and VECLIB_MAXIMUM_THREADS set to the same value; x86_64 stacks ran as arch -x86_64 python.