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 deterministic 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_chunk 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. The frame cannot handle a random measurement, because a coin flip is not reference-XOR-anything. So the reference pass returns None the moment it meets one, and the sampler falls back to the per-shot tableau path (sample_batch).

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 (knocks out the frame path)

Take a circuit the frame path can handle, a deterministic 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: one full tableau per shot
• frame_run: bit-packed frames, 64 shots per word

distance	rounds	shots	sample_batch	frame_run	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 CompiledSampler computes it once and reuses it across a whole logical-error-rate sweep.

Worked example: why the repetition memory is frame-eligible

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: a fully deterministic reference.

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

code = repetition_code(distance=5, rounds=5, p=0.02, channel="X_ERROR")
sampler = CompiledSampler(code)
instrs, _tables = sampler._compiled
ref = sampler._core.frame_reference(instrs)
print("deterministic measurements:", len(ref))
print("all read 0 in the reference:", not any(ref))
# -> deterministic measurements: 25
# -> all read 0 in the reference: True

Now break it. A GHZ state measured in the Z basis is a fair coin on the first qubit, so the reference pass gives up and the sampler uses the per-shot tableau:

from qliff import Circuit
from qliff.noise import CompiledSampler

ghz = Circuit(3)
ghz.append("H", 0)
ghz.append("CX", [0, 1])
ghz.append("CX", [1, 2])
ghz.append("M", [0, 1, 2])
gs = CompiledSampler(ghz)
print(gs._core.frame_reference(gs._compiled[0]))
# -> None

Result: the frame path is not a special case of a special code; it is precisely the circuits whose measurements are all deterministic in the noiseless run, which is most memory experiments and no logical-coin experiment.

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
from qliff.noise import CompiledSampler

for phi in (0.001, 0.01, 0.1):
    c = Circuit(1)
    c.append("X_ERROR", [0], phi)
    c.append("M", [0])
    faulted = int(CompiledSampler(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 / estimate_chunk), 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:

word_g: the rowsum phase, 64 lanes at once

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

fn word_g(xi: u64, zi: u64, xh: u64, zh: u64) -> i32 {
    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.count_ones() as i32 - minus.count_ones() as i32
}

One machine popcount per mask replaces 64 per-qubit branches. rowsum accumulates 2*r_h + 2*r_i + sum_k word_g(...) a word at a time, 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.

word_g on a full word equals the qubit-by-qubit sum of the scalar phase $g$. That is the property the core's test asserts, and it holds because the plus and minus Pauli-pair cases never share a lane.

Implementation

The tabs write the four moves down three ways. Pseudocode is the frame engine's inner loop, the one that carries a whole cohort. Std Python reproduces each optimisation in NumPy, small enough to read: the phase-free XOR propagation of one check, the geometric-skip fault list, and the bit-parallel word_g checked against a lane-by-lane sum. Qliff runs the real thing: the frame path, its fallback, and the signed estimator.

Pseudocode

# --- frame sampler: one reference pass, then bit-packed frames --------------
ref_bits = run the circuit ONCE, noiseless, recording each measurement
if any measurement was random:
    fall back to a per-shot tableau (sample_batch)      # frame can't roll a coin

for each word of 64 shots:                              # 64 shots share one u64
    fx[*] = fz[*] = 0                                    # frame = identity
    for instruction in circuit:
        gate:   XOR the frame per the phase-free table   # H swap, S xor, CX/CZ propagate
        noise:  s = 0                                     # rare-error geometric skip
                loop:
                    s += floor(ln(1-U) / ln(1-phi))      # jump to the next faulty shot
                    if s >= 64: break
                    pick a fault branch ~ |weight|
                    flip bit s of fx / fz per the Pauli
                    s += 1
        measure(q): record_word = ref_bit XOR fx[q]      # whole word, one XOR

Std Python

import numpy as np

# 1. A Pauli FRAME is phase-free, so propagation is pure XOR. One Z-parity check,
#    64 shots packed as a boolean per shot; X faults on the two data qubits copy
#    forward through CX into the ancilla frame (fx[a] ^= fx[q0] ^ fx[q1]).
rng = np.random.default_rng(0)
phi = 0.08
fx_q0 = rng.random(64) < phi
fx_q1 = rng.random(64) < phi
fx_anc = fx_q0 ^ fx_q1
record = False ^ fx_anc                 # reference bit 0, so record = fx[a]
print("shots faulted:", int((fx_q0 | fx_q1).sum()))
print("syndrome ones:", int(record.sum()))
# -> shots faulted: 12
# -> syndrome ones: 11

# 2. Rare-error skips: draw the GAP to the next faulty shot, not a coin per shot.
def faulty_shots(phi, shots, seed):
    r = np.random.default_rng(seed)
    ln1m = np.log(1 - phi)
    out, s = [], 0
    while True:
        s += int(np.log(1 - r.random()) / ln1m)   # skip = floor(ln(1-U)/ln(1-phi))
        if s >= shots:
            break
        out.append(s)
        s += 1
    return out

hits = faulty_shots(0.01, 100_000, 1)
print("faulted shots:", len(hits), " uniform draws:", len(hits) + 1)
# -> faulted shots: 982  uniform draws: 983

# 3. word_g: the measurement phase, bit-parallel. The +1 and -1 Pauli-pair cases
#    are disjoint masks, so one popcount each replaces 64 per-lane branches.
MASK = (1 << 64) - 1

def word_g(xi, zi, xh, zh):
    plus = (xi & ~zi & xh & zh) | (~xi & zi & xh & ~zh) | (xi & zi & ~xh & zh)
    minus = (xi & ~zi & ~xh & zh) | (~xi & zi & xh & zh) | (xi & zi & xh & ~zh)
    return bin(plus & MASK).count("1") - bin(minus & MASK).count("1")

def g_scalar(x1, z1, x2, z2):
    if not x1 and not z1:
        return 0
    if x1 and not z1:
        return z2 * (2 * x2 - 1)
    if not x1 and z1:
        return x2 * (1 - 2 * z2)
    return z2 - x2

r = np.random.default_rng(7)
xi, zi, xh, zh = (int(r.integers(0, 1 << 63)) for _ in range(4))
lanes = sum(
    g_scalar((xi >> k) & 1, (zi >> k) & 1, (xh >> k) & 1, (zh >> k) & 1)
    for k in range(64)
)
print("word_g:", word_g(xi, zi, xh, zh), " lane-by-lane:", lanes)   # they agree
# -> word_g: 3  lane-by-lane: 3

Qliff

from qliff import Circuit
from qliff.noise import CompiledSampler, Sampler
from qliff.qec.codes import repetition_code

# Frame path: a deterministic reference -> bit-packed frames over the shots.
mem = repetition_code(distance=5, rounds=5, p=0.02, channel="X_ERROR")
s = CompiledSampler(mem)
ref = s._core.frame_reference(s._compiled[0])
print("frame path:", ref is not None, " deterministic measurements:", len(ref))
print("detector rate:", round(float(s.sample(50_000, seed=1).mean()), 4))
# -> frame path: True  deterministic measurements: 25
# -> detector rate: 0.103

# Fallback: a random measurement is a coin the frame can't model.
gh = Circuit(2)
gh.append("H", 0)
gh.append("CX", [0, 1])
gh.append("M", [0, 1])
print("random-measurement reference:", CompiledSampler(gh)._core.frame_reference(
    CompiledSampler(gh)._compiled[0]
))
# -> random-measurement reference: None

# Signed estimator: unbiased for non-Pauli noise, whole loop in the Rust core.
q = Circuit(1)
q.append("X", [0])
q.append("AMPLITUDE_DAMP", [0], 0.3)
print("estimate <Z>:", round(Sampler(q).expect("Z", shots=200_000, seed=3), 4))
# -> estimate <Z>: -0.3975

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 word_g 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.