The Optimal Shot Allocation Rule Made Our Energies 50× Worse

The measurement problem

A VQE energy is not measured. It is assembled from pieces: E = Σᵢ cᵢ ⟨Pᵢ⟩, a weighted sum over the Pauli terms of the Hamiltonian. Benzene has 914 of them; H₆ has 919. Each ⟨Pᵢ⟩ has to be estimated by repeatedly preparing the state and measuring, and every one of those repetitions — every shot — costs real time on real hardware. The shot budget is the dominant cost of running VQE, and it is the reason VQE is hard to scale.

So: given a budget of N shots, how many should each term get? The default in most code is to divide evenly. That is obviously not optimal, because the terms are not equally uncertain. A Pauli string that happens to be close to an eigenoperator of the current state returns nearly the same answer every time — its σᵢ is nearly zero — and shots spent on it buy nothing at all.

Statistics has a name for the right answer here, and it is a century old. Neyman allocation gives each term a share proportional to |cᵢ| σᵢ. For a weighted sum of independent estimators that provably minimises the variance, and the improvement over uniform is exactly the spread in |cᵢ|σᵢ across terms — which in molecular Hamiltonians is large. PennyLane ships a related idea as Rosalin, which weights by |cᵢ| alone.

The distinction matters more than it looks. Weighting by |cᵢ| is blind to variance, so it happily pours shots into a large-coefficient term that is already known to four decimal places. That is the gap Neyman allocation closes — in principle.

What we measured

270 runs: 10 molecules from 20 to 919 Pauli terms, three parameter states each, three shot budgets, 200 repeats per scheme so the error distribution is actually resolved rather than guessed at.

The three states matter, and they are the part we got wrong the first time. Term variances depend on the quantum state, so an advantage measured at a random starting point tells you very little about the regime an optimiser actually spends its time in. We measured at a random start, 40 steps in, and fully converged — each found by exact statevector optimisation, so no sampling contaminates the choice of state.

Every scheme is charged the same total budget, and shots are counted per term so the total is exactly Σsᵢ. Neyman pays for its own pilot out of that budget. We also ran an oracle variant using the exact σᵢ from the statevector — not a usable method, but it marks the ceiling.

It looked like a triumph

Our first pass measured the standard deviation of each scheme, and Neyman allocation won everything. On LiH it cut the spread from 15.9 mHa to 1.8 — nearly nine times tighter than uniform, at identical cost. We had a headline.

It survived the obvious checks, too. The measured spreads matched the closed-form Neyman variance to within sampling noise, so the effect was not a fluke of one random draw; it was the theory working exactly as advertised.

The headline was false. Standard deviation measures how tightly the estimates cluster; it says nothing about where they cluster. Once we scored the same runs on RMSE, which counts bias as the error it is, Neyman allocation was not the best scheme in the study. It was by far the worst:

The tightest estimates in the study — a median spread half that of uniform — and centred six-tenths of a hartree away from the right answer. Chemical accuracy is 1.6 mHa.

The trap

A Pauli operator has eigenvalues ±1, so its variance is σ² = 1 − ⟨P⟩² and you estimate it by measuring ⟨P⟩ on a small pilot sample. Draw k heads out of n, set m = (2k − n)/n, take σ̂ = √(1 − m²).

When the pilot comes out all heads or all tails, m = ±1 and σ̂ is exactly zero. Not small — zero. That term is then allocated zero shots, is never measured, and contributes nothing to the energy. Its coefficient is simply deleted from the Hamiltonian.

And which terms come out all-heads? The near-deterministic ones — precisely the terms whose true σ is genuinely small. That much is self-consistent and sounds harmless. The problem is what those terms are: in a molecular Hamiltonian near a Hartree-Fock-like state, the near-deterministic Pauli strings are the Z-type diagonal terms, and those carry the largest coefficients in the whole operator.

The oracle version confirms the diagnosis. Given the exact σᵢ it funds every term, has essentially zero bias, and does beat uniform. So the allocation rule is fine. What is broken is the estimator of σ — and nothing about that is specific to quantum computing. It is a small-sample estimation failure that happens to be catastrophic here because the consequence of σ̂ = 0 is not a bad estimate but no estimate at all.

The fix

Two repairs, both cheap, and neither of them clever:

1. Shrink the variance estimate. Use the Agresti-Coull proportion p = (k + 1)/(n + 2) instead of k/n. It is never exactly 0 or 1, so σ̂ is never exactly zero, so no term is ever starved. Add a floor of one shot per term for good measure.

2. Pool the pilot with the main pass. The pilot already measured every term and you already paid for it. Discarding it is pure waste. Combine the two samples per term, weighted by shot count.

Scored on RMSE, over the same 270 runs:

Because sampling variance falls as 1/N, a 1.53× reduction in RMSE is 2.33× fewer shots for the same accuracy. The pilot costs only 6% above the oracle ceiling — knowing the variances exactly would buy you almost nothing more.

And it is not carried by one lucky molecule. It holds in every one of the ten, from 20 terms to 919:

It also holds across the optimisation trajectory — 1.46× at a random start, 2.07× mid-optimisation, 1.32× converged — and across budgets from 10⁴ to 10⁶ shots.

A correction

A post we published earlier this month, and have since withdrawn, reported 8.9× for Neyman allocation on LiH. That number was wrong in three ways, and this study was built to find out whether it was:

• It was measured at a random starting point, not where an optimiser actually works. The same molecule gives 2.0× once converged.
• It compared standard deviation only, so the bias described above was invisible. Scored on RMSE the scheme was fifty times worse than uniform, not nine times better.
• It rested on two molecules and one parameter point.

The validated figure is ~1.5× in RMSE, ~2.3× in shots. That is a smaller and much duller number than 8.9×, and it is the one we can defend. We would rather publish the dull number and the failure that produced it than keep the exciting one.

What this does not show

Stated plainly, because a result like this invites overreading:

• Shots are counted per term. Grouping commuting terms into joint measurements is a separate and larger saving; it composes with any of these schemes and is not modelled here. No claim above depends on it.
• This measures the energy estimator, not a full optimisation run. We have since run that experiment separately (1,200 optimisations, in experiments/shot_allocation_opt/) and the answer is sometimes, and not for the reason you would guess: at the budget where our archived LiH failure happened, better allocation moves the result by 0.7 mHa out of 707 — with perfect arithmetic the same run still fails at 652 mHa, so that collapse was evaluation starvation, not measurement. Given 10× the evaluations it matters a great deal, up to 27×, but by stopping optimisers from mistaking sampling noise for convergence and quitting early, not by taking better steps.
• Statevector simulation throughout — no gate noise, no readout error, no device topology. Sampling noise is zero-mean; gate noise is not.

One methodological note. Terms are sampled from Binomial(s, (1+⟨P⟩)/2) rather than by executing the circuit s times. For a ±1 observable that is the exact sampling distribution, not an approximation, and it is what makes 200 repeats across 919 terms affordable. We tested that claim rather than asserting it: both paths agree within sampling noise. The check is in the repository.

The general lesson

The bug here was not in the allocation rule, which is correct, nor really in the variance estimator, which is standard. It was in the metric. We measured spread, spread is what an allocation scheme is designed to reduce, and a scheme that reduces spread by quietly deleting half the Hamiltonian scored beautifully.

Any optimality argument that assumes a quantity is known — and every optimal-allocation argument assumes σ is known — becomes a claim about your estimator the moment you implement it. That is where it can fail, and it can fail in a direction the optimality proof never mentions.

All 270 runs, the raw cluster logs, the tooling and the sampling-model check are in the repository under experiments/shot_allocation/. Failed and superseded runs are included, as always.