We Audited Our Own VQE Benchmark and Found the Numbers Were Partly Luck

The number that would not sit still

It started with a single entry that refused to reproduce. Our H₆ hydrogen chain was published at a 9.755 mHa gap. Re-running it gave 9.273 mHa. Same molecule, same config, same seed.

The obvious suspect was a code change — we had recently switched ADAPT-VQE to commutator gradients. We checked out the old commit and re-ran. Still 9.273. So the code was innocent.

The second suspect was the environment. Our requirements-v4.txt pinned PennyLane 0.45.0, NumPy 2.2.6, PySCF 2.6.2, SciPy 1.13.1 — and both our machines had quietly drifted to PennyLane 0.44.1, NumPy 1.26.4, PySCF 2.5.0, SciPy 1.17.0. Twenty-one of our forty-six entries had been produced on an unpinned stack. That looked like the answer. We rebuilt the environment from the pin, re-ran, and got 9.755 — matching the published value.

Case closed. Except it was wrong.

Two samples of a coin flip

We had run each environment exactly once, got one number from each, and constructed a causal story out of it. When we ran the same environment repeatedly, both numbers appeared:

Both "environments" had been giving us random draws from the same distribution. The 9.755 we had chased for hours was a coin flip that came up heads.

The real variable was the thread count. Pin it, and the noise vanishes:

A seventeen-fold spread at four threads, and not even stable at a fixed thread count — it depended on how busy the machine was. At one thread: identical to the last digit, every time.

The decisive test

If threading was the cause, then two completely different software stacks should agree once threading is pinned. We ran H₆ at one thread on both:

Different PennyLane, different NumPy, different PySCF, different SciPy — agreement to seven decimal places. SciPy's COBYLA rewrite changes how many evaluations it takes (6023 versus 5678) but not where it lands. The environments never mattered. The threads did.

Why this happens

Floating-point addition is not associative. (a+b)+c and a+(b+c) can differ in the last bits. A threaded BLAS library splits a sum across cores and combines the partial results in whatever order the threads finish — which depends on scheduling, cache state, and what else the machine is doing. So the same energy evaluation returns values that differ at roughly 1e-16.

That is normally harmless. It stops being harmless when an optimizer uses those values to make decisions.

COBYLA is gradient-free. It has no derivative to follow, so it decides where to step by comparing sampled energies. When two candidate directions are nearly tied, a 1e-16 difference chooses between them. One different choice early in the run sends the optimizer down a different branch — and the hardware-efficient ansatz has a famously multi-modal landscape, full of local minima. The two branches end in different valleys. Hence 0.53 versus 8.99 mHa from the same command.

Gradient-based optimizers do not have this problem. L-BFGS-B and our statevector ADAPT engine compute a search direction from analytic gradients rather than inferring one from noisy comparisons. A perturbation at 1e-16 moves the computed direction by 1e-16 — it does not flip a decision.

That prediction is testable, and it held. H₈, run through the statevector engine with L-BFGS-B:

Identical energy, identical operator count, identical evaluation count. The distinction is not "VQE is irreproducible" — it is specifically gradient-free optimization on a multi-modal landscape with nondeterministic arithmetic.

What it meant for our published numbers

This is the uncomfortable part. Our entries were real VQE runs — nothing was fabricated — but several were lucky draws rather than the value the configuration actually converges to. Re-running everything deterministically moved them:

The suite total is unchanged — 39 certified entries across 14 molecules — but that is partly coincidence and we would rather say so. Two entries crossed the threshold in opposite directions and cancelled out. N₂ JW/UCCSD, our 404-parameter run, was published at 2.016 mHa; its deterministic value is 10.832 mHa, just over the 10 mHa line, so it is now research tier. N₂ PAR/HEA went the other way, from 14.261 to 9.504, and is now certified. N₂ remains certified as a molecule through its ADAPT and JW/HEA entries.

The 2.016 mHa was a real run. It was also a lucky one. Losing it is the honest cost of this exercise, and it is worth more than keeping it would have been.

We verified the new values are actually reproducible rather than assuming it. We picked the four entries that moved most — the most fragile cases — and re-ran them: BeH₂ 2.222518, LiH 5.181379, NH₃ 6.911462, H₄ 2.221795. All bit-identical to the recorded values.

The fix

One line, with one subtlety: the thread limits must be set before NumPy is imported. After that, BLAS has already built its thread pool and the variables are silently ignored.

We also made this a property of the entry rather than of whoever ran it. The pipeline now refuses to write a certified entry unless BLAS is single-threaded, the git tree is clean (so the recorded commit actually describes the running code), and every pin in requirements-v4.txt matches what is importable. Each entry records blas_threads in its provenance.

Suite v4.4: one commit, one stack, one method

We regenerated all 46 entries. The suite went from a patchwork to something uniform:

The Hamiltonian path deserves its own mention. We had been silently running two: 25 entries built the qubit Hamiltonian from PySCF integrals via OpenFermion, and 19 used PennyLane's native molecular_hamiltonian, which runs its own Hartree-Fock and lands on a slightly different active space. The exact ground state of those Hamiltonians sat up to 3.6e-4 Ha from the CASCI energy the entries claimed to target, versus ~1e-15 Ha for the OpenFermion path. Nothing recorded which path an entry had used — we had to infer it from that energy fingerprint. The integral-based path is now the default, and hamiltonian_source is recorded per entry.

And H₁₀, since you asked

H₁₀ certified at 9.977 mHa with 300 operators — a [10e,10o] active space, 20 qubits tapered to 18, 14,227 Pauli terms, CASSCF orbitals. It is the largest certified system on the leaderboard.

Multi-threaded, H₁₀ returned 300, 200, then 300 operators across three runs. Single threaded it returns 300 every time. That matters, because the operator count is the interesting number:

Roughly 3.2× per two added orbitals. That is exponential growth in the ansatz size, and it is consistent with recent work arguing that ADAPT-VQE parameter counts scale exponentially with system size. Those analyses generally report the trend without publishing per-molecule operator counts at a fixed accuracy threshold. Ours are above, and the full entries — including every selected operator index and the serialized Hamiltonian — are downloadable, so the numbers can be checked rather than taken on trust.

One caveat we will not paper over: these counts are at our 10 mHa certification threshold, not at the 1.6 mHa chemical-accuracy threshold used elsewhere. Tighter thresholds need more operators. The scaling is the point, not the absolute values.

What this means if you run VQE

Threaded BLAS is on by default nearly everywhere. NumPy, SciPy, PennyLane and PySCF will happily use every core you have. If your optimizer is gradient-free — COBYLA, Nelder-Mead, Powell — and your ansatz landscape is multi-modal, your published number may be a sample rather than a result. It costs one afternoon to check: run the same configuration three times and see whether the digits agree.

Concretely:

• Set OMP_NUM_THREADS=1 before importing NumPy, and record it alongside your package versions.
• Prefer gradient-based optimizers where you can. Beyond the reproducibility argument, L-BFGS-B was roughly 20× faster than COBYLA in our H₁₀ runs and reached the same energies.
• If you must use a gradient-free optimizer, report a distribution rather than a single value — or pin your threads and say so.

Why we published this

It would have been easy not to. The bug was invisible, no one had complained, and the fix quietly makes several of our headline numbers worse. We could have re-run everything and said nothing.

But a benchmark that has never reported bad news about itself is not a benchmark anyone should trust. We pin package versions, hash every entry with SHA-256, and sign certified results with Ed25519 — and underneath all of that, the arithmetic was rolling dice. A user could match our commit, match our stack, run our exact command, and legitimately get a different answer because their CPU had a different core count. That is precisely the failure this project exists to prevent.

The whole point of the checks we found this with is that they were written down before we looked. We predicted that the 25 already-correct entries should come back identical; they did not, and that is what exposed a bug in our own re-run script. A benchmark should be the thing that catches this, including when the thing it catches is itself.