Rare-event estimation

The 10⁹-shot wall

To resolve a failure rate of p_L you need roughly 1/p_L shots just to see one failure, and many multiples of that for a tight estimate. At p_L ≈ 10⁻⁷ that is hundreds of millions of decodes — out of reach. Both estimators sidestep this by importance sampling: spend shots where the failures actually live, then correct the bias analytically.

The decode_fn contract

Both estimators are decoder-agnostic: you hand them a decode_fn(syndromes) that takes a list[list[int]] and returns an array-like of shape (len(syndromes), num_observables). Any of the three decoders works — wrap relay-BP to drop its iteration counts. That single contract is why one estimator harness serves BP+OSD, relay-BP, and BP+OTF alike.

glean.dss — Dynamical Subset Sampling

DSS (Heußen 2024) stratifies the independent faults by their weight w and writes p_L(p) = Σ_w A_w(p) · p_fail^(w). The weight probabilities A_w are exact and decoder-independent (a Poisson–binomial of the priors), so only the few contributing low-w classes (w ≳ d/2) need sampling — and the dominant w=0 class is never sampled at all.

import glean, numpy as np

g = glean.load_dem(dem_text, n_det, n_obs)
decode_fn = lambda s: np.asarray(g.decode_batch(s, osd_order=10, parallel=True))

# Stratify on fault weight; only the contributing low-w classes are sampled.
est = glean.dss.estimate(g, decode_fn, weights=range(1, 13), shots=50_000)
print(est.p_L, "+/-", est.err, "bracket", (est.p_L_lo, est.p_L_hi))

DssEstimate carries p_L, err, a rigorous truncation bracket (p_L_lo, p_L_hi), and the per-class breakdown. Even a class with zero observed failures gets a positive Wilson upper bound, which is what keeps the bracket honest.

✓ B8 — unbiased on Glean's pipeline, 0.0σ vs direct MC at p=3×10⁻³.

One campaign, a whole curve

Because p_fail^(w) does not depend on the physical error rate, you decode once and recombine with A_w(p) at every target rate to trace the whole p_L(p) curve — no re-decoding.

# Decode ONCE — p_fail^(w) is independent of the physical error rate p.
table = glean.dss.sample_pfail(g, decode_fn, weights=range(1, 13), shots=50_000)

# Recombine with the (decoder-independent, exact) weight probabilities A_w(p)
# at each target rate — no re-decoding.
for p in (1e-3, 2e-3, 3e-3):
    g_p = glean.load_dem(build_dem_at(p), n_det, n_obs)   # your DEM rebuilt at rate p
    A = glean.dss.weight_probs(g_p, w_max=12)
    print(p, glean.dss.combine(table, A).p_L)

glean.splitting — the deep-number pin

DSS hits a resolution wall when a single mode dominates the deepest classes. Splitting (Bravyi–Vargo 2013; Mayer 2025) has no such wall: it writes P_L(p_target) = P_L(p_1) · ∏_j R_j as a telescoping product of ratios, each estimated by Bennett's acceptance-ratio method over a Metropolis-in-F chain.

import glean, numpy as np

g = glean.load_dem(dem_text, n_det, n_obs)

# Telescope from a tractable rate down to the deep target. The decoder is fixed at
# the target rate, so the estimate is the standard matched logical error rate.
p_seq = [3e-3, 2.4e-3, 1.8e-3, 1.4e-3, 1e-3]
priors_per_rate = [build_priors(p) for p in p_seq]        # DEM priors rebuilt per rate

est = glean.splitting.estimate(
    g, p_seq, priors_per_rate,
    seed_log_pL=np.log(2.4e-3), seed_log_pL_err=0.1,
    n_chains=8, burn_in=50_000, n_samples=4000, thin=2000, seed=0)

print(est.pL, "band", (est.pL_lo, est.pL_hi), "rhat", est.rhat)

SplittingEstimate reports pL, a band (pL_lo, pL_hi), the per-rate ratios, and rhat — the Gelman–Rubin mixing diagnostic. Trust the estimate when rhat < 1.1. See benchmarks/splitting_gross.py for the full deep-point harness, including building the per-rate priors.

✓ B9/B10 — unbiased vs direct MC at p=3×10⁻³ (0.07σ) and p=2×10⁻³ (0.30σ); pins the Bravyi Nature deep point within convention+extrapolation uncertainty.

Which to use

• DSS — fast, gives the whole curve from one campaign, and is exact where a few low-weight classes dominate. Reach for it first.
• splitting — the deep-number pin. Use it when DSS's bracket is wide because one mode dominates, or to telescope two-plus decades into the rare regime.