Optimisations: Physics Tutorial 09 of 10

The earlier pages built estimators that spend shots wisely: the noise sampler weights trajectories, the stratified sampler cancels the weight magnitude before it forms a number. This page works one level below, on the cost of a single shot. Each of the four moves below is the simulator refusing to compute something it already knows. They all live in the Rust core (src/tableau.rs), reached through qliff.noise and qliff.qec.

The Pauli frame sampler

Sampling a noisy Clifford circuit the obvious way builds a fresh tableau for every shot, replays every gate, and rolls the coins. But the gates never change between shots, and in many circuits the measurements never roll a coin. Almost all of that work is redundant.

Split the run in two. First a single noiseless reference pass, replaying the gates once and recording each measurement bit. Then, over the shots, propagate only the difference the noise makes: a Pauli frame, a phase-free two-bit-per-qubit operator (an X-frame and a Z-frame) that says which Pauli each shot carries relative to the reference. The frame drops the Pauli's phase, the overall $\pm 1$ or $\pm i$ in front of the operator, so propagating it becomes pure XOR with no i-phases and no row sums. And because 64 shots pack into one u64, a whole cohort moves through a gate in a handful of word operations.

Dropping the phase is safe because a measurement outcome depends only on whether the error commutes or anticommutes with the measured operator, which is set by its X and Z bits, not by its phase. This is the operator's phase, a different thing from the signed quasiprobability weight of the noise decomposition. That weight's sign is kept, not dropped, and the signed-weights section below is where it earns its keep.

The frame rules are the gate bit-algebra with the phase term deleted:

Every gate acting on a phase-free Pauli frame

• fx: the X-frame (which shots carry an X-part)
• fz: the Z-frame (which shots carry a Z-part)

gate	action on the frame	why
H(a)	swap fx[a], fz[a]	X-frame and Z-frame trade places
S / S†(a)	fz[a] ^= fx[a]	one map for both: the phase that split them is gone
X / Y / Z(a)	(nothing)	a Pauli only changes the phase, which the frame dropped
CX(a,b)	fx[b] ^= fx[a] fz[a] ^= fz[b]	X copies forward, Z copies back
CZ(a,b)	fz[a] ^= fx[b] fz[b] ^= fx[a]	symmetric, one pass
SWAP(a,b)	swap frames a, b	exchange the two columns

Copied from frame_planes in src/tableau.rs. Compare the CX and CZ rows against the signed versions on the gates page: identical bit moves, no phase term. S and S-dagger collapse to the same rule once the phase is dropped, so the frame cannot tell them apart and does not need to.

A deterministic measurement is then free. Its noiseless value is a known reference bit $b_{\text{ref}}$, and the noise only flips it when the frame carries an X-part on that qubit. So the recorded outcome for all 64 shots in the word is one XOR:

$$\text{record}=b_{\text{ref}}\oplus \text{fx}[q],$$

a single word operation for the whole cohort, no tableau consulted. A random measurement stays on the frame path too. The reference pass forces its outcome to 0 and records a kick: the stabilizer the measured Pauli anticommuted with, keyed by the instruction's position. Before that readout every shot word XORs the kick in times a fresh per-shot coin, which reproduces the coin flip and its correlations with later measurements (Gidney 2021). So frame_reference always returns (ref_bits, kicks), and every Pauli-noise circuit samples on frames.

Bit-packed frames through one syndrome check

Legend:

• coloured cell = shot carries an X-part on that qubit
• fx q0, fx q1 = X faults on the two data qubits
• record = ancilla word, reference XOR fx[a]
• toggle off = random final measurement (the reference records a kick, and the frame path still covers every shot)

Take a memory experiment. There the speedup is a large constant factor over building a tableau per shot, because the per-shot cost falls from "replay every gate on a tableau" to "XOR a few words":

Frame sampler vs per-shot tableau, same repetition-code memory

• sample_batch_signed: one full tableau per shot
• frame_run_folded: bit-packed frames, 64 shots per word

distance	rounds	shots	per-shot tableau	frame	ratio
5	5	100k	0.18 s	0.001 s	~160x
11	11	200k	1.11 s	0.010 s	~110x

Measured on this machine, 2026-07-23, X_ERROR at p = 0.02, best of several runs. Both paths return statistically identical detector records, and the ratio is the stable quantity while the absolute wall-clock moves with machine load and core count. The reference pass is the frame method's one serial cost, and a sampler computes it once and reuses it for every sample() call.

Worked example: a deterministic reference, then a kicked one

A distance-5 bit-flip repetition memory over 5 rounds emits 20 ancilla parities plus a final 5-qubit data readout, 25 measurements. In the noiseless run every data qubit stays in $|0⟩$ and every Z parity is even, so all 25 read 0 and no kick is needed.

from qliff.noise import Sampler
from qliff.qec.codes import repetition_code

code = repetition_code(distance=5, rounds=5, p=0.02, channel="X_ERROR")
bits = Sampler(code).reference_bits()
print("measurements:", len(bits))
print("all read 0 in the reference:", not any(bits))
# -> measurements: 25
# -> all read 0 in the reference: True

Now a coin. A GHZ state measured in the Z basis is a fair coin on the first qubit. The reference forces it to 0 and records one kick at the first M (instruction position 3) spanning $XXX$, so every shot flips all three bits together or none:

from qliff import Circuit
from qliff.noise import Sampler

ghz = Circuit(3)
ghz.append("H", 0)
ghz.append("CX", [0, 1])
ghz.append("CX", [1, 2])
ghz.append("M", [0, 1, 2])
sampler = Sampler(ghz)
print(sampler._reference)
# -> ([False, False, False], [(3, [0, 1, 2], [])])

recs = sampler.sample(10_000, seed=1)
print(sorted({tuple(r) for r in recs.tolist()}), recs[:, 0].mean())
# -> [(0, 0, 0), (1, 1, 1)] 0.509

Result: the frame path is not a special case of a special code. A deterministic measurement costs one XOR against its reference bit; a random one adds one kick, a coin word XORed into the frames it spans.

Rare-error noise: skip to the faulty shots

Inside the frame engine there is a second refusal to do redundant work. A single-qubit noise location fires with some per-shot probability $\phi$, and at QEC error rates $\phi \sim {10}^{-3}$. Visiting all 64 shots in a word to draw a fault that almost never happens is wasteful: 999 of every 1000 draws come back "no fault".

The fix uses the fact that the faulty shots are Bernoulli($\phi$), so the gaps between them are geometric. Instead of drawing a coin per shot, draw the gap directly and jump:

$$\text{skip}=\lfloor \frac{\ln (1-U)}{\ln (1-\phi )}\rfloor ,U\sim \text{Uniform}[0,1),$$

which lands on the next faulty shot without touching the quiet ones in between. The loop does one skip and one branch draw per faulty shot, so it touches $\sim \phi \cdot \text{shots}$ shots and the work drops by $\sim 1/\phi$. The two edge cases fall out cleanly: a location with no fault mass ($\phi =0$) is skipped entirely, and $\phi =1$ makes $\ln (1-\phi )$ non-finite, so every skip is zero and every shot faults.

How few shots a rare location touches

Legend:

• lit cell = a shot the location faulted on (the only shots the loop visits)
• uniform draws spent = one geometric skip per faulty shot
• work ratio = shots / faults, tracking 1 / phi

Counting the faults confirms the touch rate: a lone X_ERROR location, measured, records a 1 on the shots it flipped.

from qliff import Circuit

for phi in (0.001, 0.01, 0.1):
    c = Circuit(1)
    c.append("X_ERROR", [0], phi)
    c.append("M", [0])
    faulted = int(c.sample(100_000, seed=1).sum())
    print(f"phi={phi:<6} faulted {faulted:5d} of 100000")
# -> phi=0.001  faulted    97 of 100000
# -> phi=0.01   faulted   989 of 100000
# -> phi=0.1    faulted  9916 of 100000

At $\phi ={10}^{-3}$ the sampler wrote 97 faults into 100000 shots and did skip-arithmetic 97 times. The other 99903 shots cost nothing but the XOR that carries their (empty) frame forward.

Signed weights for non-Pauli noise

Pauli noise is a coin. Amplitude damping and coherent rotation are not. Both are still handled by one tableau per trajectory rather than a state vector, because qliff writes a non-Pauli channel as a signed quasiprobability mix of Clifford branches. The noise page builds that decomposition and the stratified page drives its variance down. This section's optimisation is only the sampling loop, which lives in Rust (ColTableau::estimate, one trajectory per shot), off the Python GIL.

Per trajectory: at each noise location, draw one branch with probability $|w_k|/\gamma$ (where $\gamma =\sum _k|w_k|$ is the location's negativity), apply that branch's Clifford ops, and multiply the running weight by $\mathrm{sign}(w_k)\gamma$. After the last gate, read the observable off the final tableau in one word-parallel pass. The estimate is the mean of $w\cdot ⟨O⟩$ over shots, unbiased for any channel:

$$⟨O⟩\approx \frac{1}{N}\sum _{\text{shots}}w⟨O{⟩}_{\text{shot}},w=\prod _i\mathrm{sign}(w_{k_i}){\gamma }_i.$$

Two details make it cheap. The branch ops are plain Clifford opcodes plus one extra, opcode 9 = reset, which is how amplitude damping's $R$ branch collapses a qubit to $|0⟩$ without leaving the stabilizer formalism. And every weight has the same magnitude $\mathrm{\Gamma }=\prod _i{\gamma }_i$, so only its sign varies from shot to shot. The whole trajectory payload is that one bit: the quasiprobability sign, the one the frame sampler could discard but this loop must keep.

One branch per location, weight = sign x Gamma

Legend:

• grey chip I = identity (no-fault) branch
• green chip + = positive-weight fault branch
• red chip - = negative-weight fault branch
• +Gamma / -Gamma = trajectory weight, sign set by the red-chip parity

Run it against a known state and the estimate lands where the density matrix says it should. Exciting a qubit and damping it at rate $p$ leaves $\rho =(1-p)|1⟩⟨1|+p|0⟩⟨0|$, so $⟨Z⟩=2p-1$:

from qliff import Circuit
from qliff.noise import Sampler

c = Circuit(1)
c.append("X", [0])                      # |1>, so damping has population to act on
c.append("AMPLITUDE_DAMP", [0], 0.3)
print(round(Sampler(c).expect("Z", shots=200_000, seed=3), 4))
# -> -0.3975       (true rho gives <Z> = 2p - 1 = -0.4)

The whole trajectory loop of draw, apply, weight, and evaluate ran in the Rust core. Python only handed it the compiled branch tables and read back one float.

Word-parallel measurement

The last optimisation is under the gate algebra, in the one place a stabilizer simulator cannot avoid a global operation: measuring a Pauli that anticommutes with the state needs a rowsum, combining two length-$2n$ Pauli rows and tracking the power of $i$ that accumulates. The gates page covers the per-gate sign bookkeeping and the one-pass CZ. What it does not cover is how that rowsum is done 64 lanes at a time.

The phase contribution of multiplying two Paulis, lane by lane, is a function $g\in \{-1,0,+1\}$ of the four bits $(x_i,z_i,x_h,z_h)$. Written per qubit it is a branch; written per word it is two masks and two popcounts. The six Pauli-pair cases that give $+1$ live on disjoint lanes, as do the six that give $-1$, so OR them together and count:

rowsum_words: the rowsum phase, 64 lanes at once

• +1 lanes: (X,Y) (Z,X) (Y,Z)
• -1 lanes: (X,Z) (Z,Y) (Y,X)

// inside rowsum_words: row h *= row i, one word k of the row at a time
let (xi, zi, xh, zh) = (xs[ib + k], zs[ib + k], xs[hb + k], zs[hb + k]);
let plus = (xi & !zi & xh & zh) | (!xi & zi & xh & !zh) | (xi & zi & !xh & zh);
let minus = (xi & !zi & !xh & zh) | (!xi & zi & xh & zh) | (xi & zi & xh & !zh);
plus_pc += plus.count_ones() as i64;
minus_pc += minus.count_ones() as i64;

One machine popcount per mask replaces 64 per-qubit branches. rowsum_words accumulates the phase 2*r_h + 2*r_i + plus_pc - minus_pc over the row's words, and the same trick drives the expectation used by the signed estimator above. Tail lanes past n are all-zero (identity), which contributes 0, so no masking is needed.

On a full word, plus - minus equals the qubit-by-qubit sum of the scalar phase $g$. It holds because the plus and minus Pauli-pair cases never share a lane, and the core's expectation_matches_scalar test checks the word-parallel sign against a qubit-by-qubit rebuild with the scalar $g$.

Four refusals to do redundant work, and they compound. The frame sampler pays for the gates once and carries 64 shots per word. The rare-error skip stops it visiting shots that do nothing. The signed estimator keeps non-Pauli noise inside one tableau. And rowsum_words measures a whole word in two popcounts. None of them change a single number the sampler returns. They change how little it has to compute to return it, which is the entire budget for a scalable noise simulation. The logical error rate is what all of this is spent on. The next page, Optimisations: Memory, is where the tableau itself gets smaller.