DEMs built with `qec.TwoQubitDepolarization` / `qec.TwoQubitBitFlip` silently omit all noise on the second qubit: the MSM disagrees with the sampler

[YaegakiY]

Summary

The QEC library's two-qubit noise channels (cudaq::qec::two_qubit_depolarization, cudaq::qec::two_qubit_bitflip) are documented as applying noise on both qubits. The C++ header promises "a bit flip on either qubit independently" / "a depolarization channel on either qubit independently" (libs/qec/include/cudaq/qec/noise_model.h:52-53, 76-77), the Python bindings say the same ("Models independent bit flip errors after a two-qubit operation.", libs/qec/python/bindings/py_code.cpp:258; analogously for TwoQubitDepolarization at :250), and the Kraus operators are explicitly constructed as the two-qubit product channel including the second-operand and both-operand terms (noise_model.h:58-67 for the bit flip: E1 = K0⊗K1, E2 = K1⊗K0, E3 = K1⊗K1; :92-97 for the depolarization: the full 16-term kron(k1, k2) double loop). The stim-target sampler applies exactly that. But the MSM / DEM machinery (compute_msm, dem_from_memory_circuit and friends) enumerates error mechanisms for the first operand only: no mechanism in the resulting detector error model ever flips anything caused by the second qubit's noise.

So every DEM built with these channels models roughly half of the noise that the same NoiseModel injects during sampling, with no error or warning. Decoders configured from such DEMs are calibrated against a substantially weaker error model than the circuit they decode.

The built-in cudaq.Depolarization2 is handled correctly (all 15 mechanisms), which isolates the problem to how these custom channels are interpreted.

Reproduction

Single CX with TwoQubitDepolarization(0.15) attached, both qubits measured. Compare the MSM mechanisms with the sampled flip statistics of the same noise model. Reproduced on cudaq-qec 0.6.0 / CUDA-Q 0.14.0, stim target (ghcr.io/nvidia/cudaqx image).

import cudaq
import cudaq_qec as qec

P = 0.15
SHOTS = 200_000
cudaq.set_target("stim")


@cudaq.kernel
def two_qubit_circuit():
    q = cudaq.qvector(2)
    x.ctrl(q[0], q[1])
    mz(q[0])
    mz(q[1])


def show(channel, label):
    noise = cudaq.NoiseModel()
    noise.add_all_qubit_channel("x", channel, num_controls=1)
    cudaq.set_noise(noise)
    msm, dims, probs, ids = qec.compute_msm(two_qubit_circuit, False)
    print(f"--- {label} ---")
    print(f"MSM error mechanisms ({dims[1]}):  (flip string = m0,m1)")
    for i, (p_, id_) in enumerate(zip(probs, ids)):
        print(f"  mechanism {i}: p={p_:.6f}  flips={msm[i]}  id={id_}")
    res = cudaq.sample(two_qubit_circuit, shots_count=SHOTS)
    counts = dict(res.items())
    m0 = sum(v for k, v in counts.items() if k[0] == "1") / SHOTS
    m1 = sum(v for k, v in counts.items() if k[1] == "1") / SHOTS
    p11 = counts.get("11", 0) / SHOTS
    print(f"sampled ({SHOTS} shots): P(m0 flips)={m0:.4f}  "
          f"P(m1 flips)={m1:.4f}  P(both)={p11:.4f}\n")
    cudaq.unset_noise()


print(f"expected for D1(p) (x) D1(p): per-qubit flip marginal 2p/3 = {2*P/3:.4f}, "
      f"P(both) = {(2*P/3)**2:.4f}\n")
show(qec.TwoQubitDepolarization(P), "qec.TwoQubitDepolarization(0.15)")
show(qec.TwoQubitBitFlip(P), "qec.TwoQubitBitFlip(0.15)")
show(cudaq.Depolarization2(P), "control: built-in cudaq.Depolarization2(0.15)")

Observed

expected for D1(p) (x) D1(p): per-qubit flip marginal 2p/3 = 0.1000, P(both) = 0.0100

--- qec.TwoQubitDepolarization(0.15) ---
MSM error mechanisms (3):  (flip string = m0,m1)
  mechanism 0: p=0.050000  flips=10  id=0
  mechanism 1: p=0.050000  flips=10  id=0
  mechanism 2: p=0.050000  flips=00  id=0
sampled (200000 shots): P(m0 flips)=0.0998  P(m1 flips)=0.0992  P(both)=0.0101

--- qec.TwoQubitBitFlip(0.15) ---
MSM error mechanisms (1):  (flip string = m0,m1)
  mechanism 0: p=0.150000  flips=10  id=0
sampled (200000 shots): P(m0 flips)=0.1508  P(m1 flips)=0.1498  P(both)=0.0225

--- control: built-in cudaq.Depolarization2(0.15) ---
MSM error mechanisms (15): all four flip patterns 01/10/11/00 present, p=0.01 each
sampled (200000 shots): P(m0 flips)=0.0794  P(m1 flips)=0.0803  P(both)=0.0398

The sampler flips m1 at the documented rate 2p/3 ≈ 0.10 (and P(both) = (2p/3)² confirms the independent product channel), but the MSM contains zero mechanisms that flip m1. For TwoQubitBitFlip the single retained mechanism also carries probability p instead of the product-channel value. The built-in Depolarization2 control shows the MSM machinery itself handles genuine two-qubit channels correctly.

(Note on probability semantics: MSM mechanism probabilities are per-case channel probabilities, i.e. mutually exclusive alternatives, as the control's 15 mechanisms at exactly p/15 show. They are not stim-DEM-style independent mechanism rates. Under that convention the expected TwoQubitBitFlip enumeration is three flip cases 10/01/11 at p(1−p), p(1−p), p², so the retained mechanism's expected value is p(1−p), and the expected TwoQubitDepolarization enumeration is the 15 non-identity products P1⊗P2 with probabilities p(P1)·p(P2), where p(I) = 1−p and p(X) = p(Y) = p(Z) = p/3.)

Expected

The mechanisms enumerated by compute_msm / dem_from_memory_circuit describe the same error process the sampler applies: for TwoQubitDepolarization(p) on a CX, single-qubit depolarizing mechanisms on both operands (and for TwoQubitBitFlip(p), X mechanisms on both operands).

Root cause (from reading the sources)

1. libs/qec/include/cudaq/qec/noise_model.h: two_qubit_depolarization builds the correct 16-operator product-channel Kraus set, but tags it noise_type = noise_model_type::depolarization_channel (the single-qubit depolarizing tag) with parameters = {p}. (two_qubit_bitflip analogously gets a single-qubit-style tag.)
2. runtime/nvqir/stim/StimCircuitSimulator.cpp, isValidStimNoiseChannel: dispatches on the tag only ("Check the old way first") and never inspects the Kraus operators. The single-qubit depolarizing tag maps to StimNoiseType{stim_name="DEPOLARIZE1", 3 mechanisms at p/3, num_targets=1}.
3. The two consumers of that StimNoiseType then diverge:
   • Sampling branch: noiseOps.safe_append_u(stim_name, qubits, parameters) appends DEPOLARIZE1(p) with the full two-qubit target list; under stim's multi-target semantics this is independent depolarizing on each qubit, which happens to match the documented channel.
   • MSM branch: the mechanism loop iterates t < res->num_targets (= 1), so each mechanism's flips are recorded on qubits[0] only. The second qubit's noise never enters the MSM, hence never enters any DEM.

So the same StimNoiseType is interpreted as "broadcast to all targets" by the sampler and "first target only" by the MSM. Any single-qubit-tagged channel attached to a multi-qubit gate hits this divergence; the QEC library's own two-qubit channels are the prominent case.

Impact

• Any DEM produced by dem_from_memory_circuit / x_dem... / z_dem... / compute_msm under these channels omits all second-operand error mechanisms. On a surface_code d=3, 3-round memory DEM with TwoQubitDepolarization(0.01), 70 of the 160 distinct (detector set, observable flip) symptom classes of the exact channel (44%) are entirely absent; 62 of those 70 are unreachable by first-operand-only enumeration regardless of any later DEM processing, the remaining 8 drop out at the canonicalization stage. (Method: independent frame-propagation enumeration of all 144 single-qubit depolarizing instances using the DEM's own row convention. Method validation: emulating the first-operand-only enumeration reproduces the DEM's symptom-class structure exactly, all 90 of its classes with no class in the DEM that the emulation cannot produce, and bit-matches the probabilities on 70 of the 90 classes (median deviation 0); the residual deviations trace to canonicalization-stage merge behavior downstream of mechanism enumeration, not to the enumeration comparison itself.)
• Decoders configured from these DEMs are calibrated against roughly half the physical noise; logical error rates and thresholds derived from them are silently optimistic.
• Sampling results are unaffected (the sampler matches the channel's documentation), which makes the mismatch hard to notice: the tests that exercise these channels (libs/qec/python/tests/test_code.py, libs/qec/unittests/backend-specific/stim/test_qec_stim.cpp) all go through the sampling path and assert only syndrome shapes and the presence of nonzero entries; nothing compares the MSM/DEM against the sampled statistics.
• The shipped examples use the built-in cudaq.Depolarization2 and are unaffected, but the two channels are public, documented API of the QEC library.

Suggested fix

Either layer works; both would be more robust:

1. In noise_model.h, construct these channels as cudaq::pauli2 with the exact 15 product probabilities (for the depolarizing product: P(P1⊗P2) = p(P1)·p(P2) with p(I) = 1−p, p(X)=p(Y)=p(Z) = p/3; for the bit-flip product: P(X⊗I) = p(1−p), P(I⊗X) = (1−p)p, P(X⊗X) = p²). pauli2 is already handled correctly by the stim target, so no backend change is needed. Verified on 0.6.0: attaching cudaq.Pauli2 with these 15 product probabilities to the same single-CX circuit yields the correct 15-mechanism MSM covering both operands (per-mechanism probabilities exact to 1e-12, flip signatures correct for every Pauli pair) with sampler statistics unchanged (all within 1.2σ of 2p/3 and (2p/3)² at 200k shots). Until the library is fixed, this construction also works as a user-side workaround.
2. In StimCircuitSimulator.cpp, make the MSM branch mirror the sampling branch's broadcast semantics (enumerate the 1-target mechanisms once per listed target), or reject loudly when a 1-target StimNoiseType is applied to more than one qubit.

Environment

cudaq-qec 0.6.0 (cudaqx 84d18ca), CUDA-Q 0.14.0 (cuda-quantum d845683), stim target, ghcr.io/nvidia/cudaqx:latest Docker image (arm64-cu13-0.14.0), CPU only. File and line references throughout are to these two commits.

[bmhowe23]

@YaegakiY - thank you for the very detailed bug report. I posted a draft PR with a fix for CUDA-Q based on your analysis above, but as mentioned in that PR, this area of the code is undergoing some changes/improvements right now (switching to a new "dem" mode vice "msm" mode, and this bug only appears in the "msm" mode). Those ongoing changes might eliminate the need for that CUDA-Q PR, but we will keep you posted in the coming days.

For now, I am hesitant to make the proposed change to noise_model.h because the current description of those noise classes are describing independent noise mechanisms on each qubit (shown below), suggesting only a total of 6 errors should be added (X/Y/Z for each of 2 qubits) instead of the 15 that would come from pauli2 or depolarization2. E.g., from two_qubit_depolarization:

  /// @brief Construct a two qubit kraus channel that applies a depolarization
  /// channel on either qubit independently.

One could, of course, argue that the depolarization2 is the more preferred noise class for most cases, and it is certainly confusing that depolarization2 differs in behavior from two_qubit_depolarization. Do you have any opinions here?