From fcc2c87f3de8db5a335d0d1218ae3845fc8f3a25 Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Fri, 6 Nov 2026 13:38:51 +0100 Subject: [PATCH 1/2] Marginalize two-qubit pauli error on lost partners --- crates/ppvm-pauli-sum/src/sum/noise.rs | 116 +++++++++--------- crates/ppvm-pauli-sum/tests/loss.rs | 42 +++++++ .../ppvm-tableau-sum/tests/sampler_vs_pure.rs | 10 +- crates/ppvm-tableau/src/noise.rs | 66 ++++++++-- crates/ppvm-tableau/src/tableau_like.rs | 43 +++++-- 5 files changed, 203 insertions(+), 74 deletions(-) diff --git a/crates/ppvm-pauli-sum/src/sum/noise.rs b/crates/ppvm-pauli-sum/src/sum/noise.rs index d6dc3aec1..713655b56 100644 --- a/crates/ppvm-pauli-sum/src/sum/noise.rs +++ b/crates/ppvm-pauli-sum/src/sum/noise.rs @@ -48,57 +48,63 @@ fn one_minus_two_sum(p: &[C; 15], indices: [usiz impl TwoQubitPauliError for PauliSum { fn two_qubit_pauli_error(&mut self, addr0: usize, addr1: usize, p: [::Coeff; 15]) { - self.scale(|k, v| match (k.get(addr0), k.get(addr1)) { - (Pauli::I, Pauli::I) => {} - (Pauli::I, Pauli::X) => { - *v *= one_minus_two_sum(&p, [1, 10, 13, 14, 2, 5, 6, 9]); - } - (Pauli::I, Pauli::Y) => { - *v *= one_minus_two_sum(&p, [0, 10, 12, 14, 2, 4, 6, 8]); - } - (Pauli::I, Pauli::Z) => { - *v *= one_minus_two_sum(&p, [0, 1, 12, 13, 4, 5, 8, 9]); - } - (Pauli::X, Pauli::I) => { - *v *= one_minus_two_sum(&p, [10, 11, 12, 13, 14, 7, 8, 9]); - } - (Pauli::X, Pauli::X) => { - *v *= one_minus_two_sum(&p, [1, 11, 12, 2, 5, 6, 7, 8]); - } - (Pauli::X, Pauli::Y) => { - *v *= one_minus_two_sum(&p, [0, 11, 13, 2, 4, 6, 7, 9]); - } - (Pauli::X, Pauli::Z) => { - *v *= one_minus_two_sum(&p, [0, 1, 10, 11, 14, 4, 5, 7]); - } - (Pauli::Y, Pauli::I) => { - *v *= one_minus_two_sum(&p, [11, 12, 13, 14, 3, 4, 5, 6]); - } - (Pauli::Y, Pauli::X) => { - *v *= one_minus_two_sum(&p, [1, 10, 11, 12, 2, 3, 4, 9]); - } - (Pauli::Y, Pauli::Y) => { - *v *= one_minus_two_sum(&p, [0, 10, 11, 13, 2, 3, 5, 8]); - } - (Pauli::Y, Pauli::Z) => { - *v *= one_minus_two_sum(&p, [0, 1, 11, 14, 3, 6, 8, 9]); - } - (Pauli::Z, Pauli::I) => { - *v *= one_minus_two_sum(&p, [10, 3, 4, 5, 6, 7, 8, 9]); - } - (Pauli::Z, Pauli::X) => { - *v *= one_minus_two_sum(&p, [1, 13, 14, 2, 3, 4, 7, 8]); - } - (Pauli::Z, Pauli::Y) => { - *v *= one_minus_two_sum(&p, [0, 12, 14, 2, 3, 5, 7, 9]); - } - (Pauli::Z, Pauli::Z) => { - *v *= one_minus_two_sum(&p, [0, 1, 10, 12, 13, 3, 6, 7]); - } - _ => { - // NOTE: if just one atom is lost, then there is no - // well-defined noise channel on the other atom - // so we don't apply any noise + self.scale(|k, v| { + // A lost qubit drops out: treating its slot as I applies the + // survivor's marginal channel. + let pair = match (k.get(addr0), k.get(addr1)) { + (Pauli::L, Pauli::L) => return, + (Pauli::L, p1) => (Pauli::I, p1), + (p0, Pauli::L) => (p0, Pauli::I), + pair => pair, + }; + match pair { + (Pauli::I, Pauli::I) => {} + (Pauli::I, Pauli::X) => { + *v *= one_minus_two_sum(&p, [1, 10, 13, 14, 2, 5, 6, 9]); + } + (Pauli::I, Pauli::Y) => { + *v *= one_minus_two_sum(&p, [0, 10, 12, 14, 2, 4, 6, 8]); + } + (Pauli::I, Pauli::Z) => { + *v *= one_minus_two_sum(&p, [0, 1, 12, 13, 4, 5, 8, 9]); + } + (Pauli::X, Pauli::I) => { + *v *= one_minus_two_sum(&p, [10, 11, 12, 13, 14, 7, 8, 9]); + } + (Pauli::X, Pauli::X) => { + *v *= one_minus_two_sum(&p, [1, 11, 12, 2, 5, 6, 7, 8]); + } + (Pauli::X, Pauli::Y) => { + *v *= one_minus_two_sum(&p, [0, 11, 13, 2, 4, 6, 7, 9]); + } + (Pauli::X, Pauli::Z) => { + *v *= one_minus_two_sum(&p, [0, 1, 10, 11, 14, 4, 5, 7]); + } + (Pauli::Y, Pauli::I) => { + *v *= one_minus_two_sum(&p, [11, 12, 13, 14, 3, 4, 5, 6]); + } + (Pauli::Y, Pauli::X) => { + *v *= one_minus_two_sum(&p, [1, 10, 11, 12, 2, 3, 4, 9]); + } + (Pauli::Y, Pauli::Y) => { + *v *= one_minus_two_sum(&p, [0, 10, 11, 13, 2, 3, 5, 8]); + } + (Pauli::Y, Pauli::Z) => { + *v *= one_minus_two_sum(&p, [0, 1, 11, 14, 3, 6, 8, 9]); + } + (Pauli::Z, Pauli::I) => { + *v *= one_minus_two_sum(&p, [10, 3, 4, 5, 6, 7, 8, 9]); + } + (Pauli::Z, Pauli::X) => { + *v *= one_minus_two_sum(&p, [1, 13, 14, 2, 3, 4, 7, 8]); + } + (Pauli::Z, Pauli::Y) => { + *v *= one_minus_two_sum(&p, [0, 12, 14, 2, 3, 5, 7, 9]); + } + (Pauli::Z, Pauli::Z) => { + *v *= one_minus_two_sum(&p, [0, 1, 10, 12, 13, 3, 6, 7]); + } + _ => unreachable!(), } }) } @@ -119,10 +125,10 @@ impl Depolarizing2 for PauliSum { fn depolarize2(&mut self, addr0: usize, addr1: usize, p: T::Coeff) { let factor = T::Coeff::from(1.0) - p * (16.0 / 15.0); self.scale(move |k, v| { - if k.get_lbit(addr0) || k.get_lbit(addr1) { - return; - } - if pauli_code(k, addr0) != 0 || pauli_code(k, addr1) != 0 { + // A lone survivor sees depolarize1(4p/5), whose factor is also 1 - 16p/15. + let active0 = !k.get_lbit(addr0) && pauli_code(k, addr0) != 0; + let active1 = !k.get_lbit(addr1) && pauli_code(k, addr1) != 0; + if active0 || active1 { *v *= factor.clone(); } }); diff --git a/crates/ppvm-pauli-sum/tests/loss.rs b/crates/ppvm-pauli-sum/tests/loss.rs index 924fe019e..35779c724 100644 --- a/crates/ppvm-pauli-sum/tests/loss.rs +++ b/crates/ppvm-pauli-sum/tests/loss.rs @@ -412,3 +412,45 @@ fn test_correlated_loss_ghz() { // XX should be affected by correlated loss assert!((x_exp - 0.9).abs() < 1e-8); } + +/// Coefficient of the single term `word` after applying `channel`. +fn scaled_coeff(word: &str, channel: impl FnOnce(&mut LossyPauliSum)) -> f64 { + let mut state = LossyPauliSum::builder().n_qubits(2).build(); + state += (word, 1.0); + channel(&mut state); + assert_eq!(state.data().iter().count(), 1, "{word} must not branch"); + state.data().iter().map(|(_, v)| *v).sum() +} + +#[test] +fn test_two_qubit_pauli_error_lost_slot_acts_as_identity() { + // Paulis act trivially on the lost level, so E†(|L⟩⟨L| ⊗ P) scales like + // E†(I ⊗ P): the survivor sees the marginal channel and nothing branches. + let p: [f64; 15] = std::array::from_fn(|k| 0.002 * (k + 1) as f64); + for pauli in ["X", "Y", "Z"] { + for (lost, ident) in [ + (format!("L{pauli}"), format!("I{pauli}")), + (format!("{pauli}L"), format!("{pauli}I")), + ] { + let c_lost = scaled_coeff(&lost, |s| s.two_qubit_pauli_error(0, 1, p)); + let c_ident = scaled_coeff(&ident, |s| s.two_qubit_pauli_error(0, 1, p)); + assert!(c_lost < 1.0, "{lost} must be damped"); + assert!( + (c_lost - c_ident).abs() < 1e-12, + "{lost}: {c_lost} vs {ident}: {c_ident}" + ); + + let d_lost = scaled_coeff(&lost, |s| s.depolarize2(0, 1, 0.1)); + let d_ident = scaled_coeff(&ident, |s| s.depolarize2(0, 1, 0.1)); + assert!( + (d_lost - d_ident).abs() < 1e-12, + "{lost}: {d_lost} vs {ident}: {d_ident}" + ); + } + } + assert_eq!( + scaled_coeff("LL", |s| s.two_qubit_pauli_error(0, 1, p)), + 1.0 + ); + assert_eq!(scaled_coeff("LL", |s| s.depolarize2(0, 1, 0.1)), 1.0); +} diff --git a/crates/ppvm-tableau-sum/tests/sampler_vs_pure.rs b/crates/ppvm-tableau-sum/tests/sampler_vs_pure.rs index 75ea3dd55..666ded5f8 100644 --- a/crates/ppvm-tableau-sum/tests/sampler_vs_pure.rs +++ b/crates/ppvm-tableau-sum/tests/sampler_vs_pure.rs @@ -1106,8 +1106,9 @@ fn two_qubit_pauli_error_nonuniform_on_ground_state() { } #[test] -fn two_qubit_pauli_error_on_lost_qubit_is_noop() { - // If either input qubit is lost the channel short-circuits to a no-op. +fn two_qubit_pauli_error_on_lost_qubit() { + // TableauSum still no-ops on a lost partner (marginalisation follows with + // its leakage support); the pure tableau applies the survivor's marginal. let shots = 1000; let p = [1.0 / 15.0_f64; 15]; let sum = run_sum(2, shots, 1e-12, |t| { @@ -1119,7 +1120,10 @@ fn two_qubit_pauli_error_on_lost_qubit_is_noop() { t.two_qubit_pauli_error(0, 1, p); }); assert!(sum.iter().all(|s| s == &vec![None, Some(false)])); - assert!(pure.iter().all(|s| s == &vec![None, Some(false)])); + // 8 of the 15 uniform-weight Paulis carry X or Y on the survivor. + assert!(pure.iter().all(|s| s[0].is_none())); + let flipped = pure.iter().filter(|s| s[1] == Some(true)).count() as f64 / shots as f64; + assert!((flipped - 8.0 / 15.0).abs() < 0.06, "flip rate {flipped}"); } #[test] diff --git a/crates/ppvm-tableau/src/noise.rs b/crates/ppvm-tableau/src/noise.rs index f09847222..ac6985889 100644 --- a/crates/ppvm-tableau/src/noise.rs +++ b/crates/ppvm-tableau/src/noise.rs @@ -608,15 +608,65 @@ mod tests { assert!(t.is_lost(1)); } - #[test] - fn two_qubit_pauli_error_first_lost_no_apply() { - // addr0 lost; p[0] = 1.0 (IX) → marginal p_x for addr1 = 1.0 - let mut t = tab(2); - t.qubit_status[0] = QubitStatus::Lost; - let mut p = [0.0f64; 15]; - p[0] = 1.0; // IX + /// Survivor's Z outcome for channel `p` with `dead` set to `status` + /// (or both live when `status` is `None`). + fn survivor_flipped(p: [f64; 15], dead: usize, status: Option, seed: u64) -> bool { + let mut t = TestTab::new_with_seed(2, 1e-12, seed); + if let Some(s) = status { + t.qubit_status[dead] = s; + } t.two_qubit_pauli_error(0, 1, p); - assert!(!t.measure(1).unwrap()); // nothing applied to addr1 + t.measure(1 - dead).unwrap() + } + + #[test] + fn two_qubit_pauli_error_survivor_independent_of_partner() { + // The survivor's reduced channel can't depend on whether its partner is + // present. Each `p` flips each qubit in Z deterministically, so one shot + // decides. Index of P_a ⊗ P_b is 4a + b - 1 (I=0, X=1, Y=2, Z=3). + let channels: [&[(usize, f64)]; 4] = [ + &[(0, 0.4), (12, 0.6)], // IX, ZX + &[(6, 0.5), (7, 0.5)], // XZ, YI + &[(3, 1.0)], // XI + &[(9, 0.2), (4, 0.3), (8, 0.5)], // YY, XX, YX + ]; + for entries in channels { + let mut p = [0.0; 15]; + for &(k, v) in entries { + p[k] = v; + } + for dead in [0, 1] { + let live = survivor_flipped(p, dead, None, 0); + for status in [QubitStatus::Lost, QubitStatus::Leaked] { + assert_eq!( + survivor_flipped(p, dead, Some(status), 0), + live, + "p = {entries:?}, dead = {dead}, {status:?}" + ); + } + } + } + } + + #[test] + fn depolarize2_lost_partner_survivor_flip_rate() { + // 8 of the 15 two-qubit Paulis carry X or Y on a given qubit, so its Z + // outcome flips with probability 8/15 whether or not the partner is there. + let shots = 4000; + for status in [None, Some(QubitStatus::Lost), Some(QubitStatus::Leaked)] { + let flips = (0..shots) + .filter(|&seed| { + let mut t = TestTab::new_with_seed(2, 1e-12, seed); + if let Some(s) = status { + t.qubit_status[0] = s; + } + t.depolarize2(0, 1, 1.0); + t.measure(1).unwrap() + }) + .count(); + let rate = flips as f64 / shots as f64; + assert!((rate - 8.0 / 15.0).abs() < 0.04, "{status:?}: rate {rate}"); + } } // === Depolarizing2 === diff --git a/crates/ppvm-tableau/src/tableau_like.rs b/crates/ppvm-tableau/src/tableau_like.rs index c72c0b540..9ac5c94fe 100644 --- a/crates/ppvm-tableau/src/tableau_like.rs +++ b/crates/ppvm-tableau/src/tableau_like.rs @@ -8,7 +8,7 @@ //! Depolarizing2, TwoQubitPauliError). New tableau backends only need to //! provide [`TableauLike::rng_mut`](crate::tableau_like::TableauLike::rng_mut) and (optionally) override //! [`TableauLike::is_qubit_lost`](crate::tableau_like::TableauLike::is_qubit_lost) -//! (which reports lost *or* leaked qubits, both skipped by the noise channels). +//! (which reports lost *or* leaked qubits; two-qubit channels marginalize them out). use num::Zero; use rand::{Rng, RngExt}; @@ -33,6 +33,21 @@ fn is_probability>(p: &C) -> bool { *p <= 1.0 } +/// Marginal X/Y/Z probabilities on one qubit of a two-qubit Pauli channel. +/// +/// `p` uses the [`TableauLike::two_qubit_pauli_error_impl`] layout, where +/// `p[4a + b - 1]` is the probability of `P_a ⊗ P_b` (0 = I, 1 = X, 2 = Y, +/// 3 = Z). Sums out the other qubit's Pauli: `q_s = Σ_o p(s, o)` when keeping +/// the first qubit, `q_s = Σ_o p(o, s)` when keeping the second. +fn marginal_pauli_probs(p: &[C; 15], keep_first: bool) -> [C; 3] { + core::array::from_fn(|i| { + let s = i + 1; + (0..4) + .map(|o| if keep_first { 4 * s + o } else { 4 * o + s }) + .fold(C::zero(), |acc, k| acc + p[k - 1].clone()) + }) +} + /// A stabilizer-tableau-like backend that supports Clifford gates and an RNG. /// /// Implementing this trait grants default implementations of the Pauli noise @@ -50,7 +65,7 @@ pub trait TableauLike: Clifford { fn rng_mut(&mut self) -> &mut Self::Rng; /// Whether the qubit at `addr` is inactive (lost or leaked), so the - /// Pauli noise channels skip it. Default: never inactive. + /// two-qubit Pauli channels act only on the surviving qubit. Default: never inactive. #[inline] fn is_qubit_lost(&self, _addr: usize) -> bool { false @@ -104,12 +119,27 @@ pub trait TableauLike: Clifford { } /// Two-qubit Pauli-error channel (15 non-identity Pauli combinations). + /// + /// If exactly one qubit is lost or leaked, the survivor gets the marginal + /// single-qubit channel (the dead qubit's Pauli is summed out). If both + /// are, nothing is applied. #[inline] fn two_qubit_pauli_error_impl(&mut self, addr0: usize, addr1: usize, p: [Self::Coeff; 15]) { - if self.is_qubit_lost(addr0) || self.is_qubit_lost(addr1) { + debug_assert!(p.iter().all(is_probability)); + + let is_lost0 = self.is_qubit_lost(addr0); + let is_lost1 = self.is_qubit_lost(addr1); + if is_lost0 && is_lost1 { return; } - debug_assert!(p.iter().all(is_probability)); + + if is_lost0 || is_lost1 { + let addr = if is_lost0 { addr1 } else { addr0 }; + let marginal_probs = marginal_pauli_probs(&p, !is_lost0); + self.pauli_error_impl(addr, marginal_probs); + return; + } + let r = self.rng_mut().random::(); let sum = Self::Coeff::zero(); let idx = p @@ -154,12 +184,9 @@ pub trait TableauLike: Clifford { } /// Two-qubit depolarizing channel: spreads `p` over the 15 non-identity - /// two-qubit Pauli errors. + /// two-qubit Pauli errors. A lone surviving qubit is depolarized with `4p/5`. #[inline] fn depolarize2_impl(&mut self, addr0: usize, addr1: usize, p: Self::Coeff) { - if self.is_qubit_lost(addr0) || self.is_qubit_lost(addr1) { - return; - } debug_assert!(is_probability(&p)); let p_arr: [Self::Coeff; 15] = core::array::from_fn(|_| p.clone() * (1.0 / 15.0)); self.two_qubit_pauli_error_impl(addr0, addr1, p_arr); From 3fd99341c99b22bb1338a54fa9b7b6ab72ab6a75 Mon Sep 17 00:00:00 2001 From: David Plankensteiner Date: Fri, 6 Nov 2026 13:52:04 +0100 Subject: [PATCH 2/2] fix(tableau): tolerate rounding in derived marginal probabilities Summed marginals of a normalized two-qubit channel can round to 1 + eps, tripping the debug probability check on a lost/leaked partner. Co-Authored-By: Claude Opus 5.5 (1M context) --- crates/ppvm-tableau/src/noise.rs | 18 ++++++++++++++++++ crates/ppvm-tableau/src/tableau_like.rs | 5 +++-- 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/crates/ppvm-tableau/src/noise.rs b/crates/ppvm-tableau/src/noise.rs index ac6985889..a26be0ac0 100644 --- a/crates/ppvm-tableau/src/noise.rs +++ b/crates/ppvm-tableau/src/noise.rs @@ -648,6 +648,24 @@ mod tests { } } + #[test] + fn two_qubit_pauli_error_saturated_marginal_on_lost_partner() { + // A normalized channel with X always on the survivor must flip it every + // shot, even though the summed marginal rounds to 1 + ε in f64. + let probs = [0.2, 0.4, 0.3, 0.1]; + for dead in [0, 1] { + let mut p = [0.0; 15]; + for (o, &v) in probs.iter().enumerate() { + // survivor = addr0: XI, XX, XY, XZ; survivor = addr1: IX, XX, YX, ZX + let k = if dead == 1 { 4 + o } else { 4 * o + 1 }; + p[k - 1] = v; + } + for seed in 0..20 { + assert!(survivor_flipped(p, dead, Some(QubitStatus::Lost), seed)); + } + } + } + #[test] fn depolarize2_lost_partner_survivor_flip_rate() { // 8 of the 15 two-qubit Paulis carry X or Y on a given qubit, so its Z diff --git a/crates/ppvm-tableau/src/tableau_like.rs b/crates/ppvm-tableau/src/tableau_like.rs index 9ac5c94fe..6005d97e6 100644 --- a/crates/ppvm-tableau/src/tableau_like.rs +++ b/crates/ppvm-tableau/src/tableau_like.rs @@ -18,7 +18,8 @@ use ppvm_traits::traits::{Clifford, Coefficient}; // `Zero` is used for `Self::Coeff::zero()` inside default method bodies; the // bound itself is redundant on `Self::Coeff` because `Coefficient: num::Zero`. -/// `true` iff `p` lies in the closed unit interval `[0, 1]`. +/// `true` iff `p` lies in the closed unit interval `[0, 1]`, up to a `1e-12` +/// rounding slack above 1 (e.g. summed marginals of a normalized channel). /// /// Spelled out with two explicit comparisons rather than /// `(0.0..=1.0).contains(p)` because the coefficient types here only implement @@ -30,7 +31,7 @@ fn is_probability>(p: &C) -> bool { if *p < 0.0 { return false; } - *p <= 1.0 + *p <= 1.0 + 1e-12 } /// Marginal X/Y/Z probabilities on one qubit of a two-qubit Pauli channel.