Skip to content

Marginalize two-qubit pauli error on lost partners #284

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
david-pl merged 2 commits into from
Nov 6, 2026
Merged

Marginalize two-qubit pauli error on lost partners #284

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
fcc2c87
Marginalize two-qubit pauli error on lost partners
david-pl Nov 6, 2026
3fd9934
fix(tableau): tolerate rounding in derived marginal probabilities
david-pl Nov 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
116 changes: 61 additions & 55 deletions crates/ppvm-pauli-sum/src/sum/noise.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,57 +48,63 @@ fn one_minus_two_sum<C: Coefficient, const N: usize>(p: &[C; 15], indices: [usiz

impl<T: Config> TwoQubitPauliError<T> for PauliSum<T> {
fn two_qubit_pauli_error(&mut self, addr0: usize, addr1: usize, p: [<T as Config>::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!(),
}
})
}
Expand All @@ -119,10 +125,10 @@ impl<T: Config> Depolarizing2<T> for PauliSum<T> {
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();
}
});
Expand Down
42 changes: 42 additions & 0 deletions crates/ppvm-pauli-sum/tests/loss.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
10 changes: 7 additions & 3 deletions crates/ppvm-tableau-sum/tests/sampler_vs_pure.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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| {
Expand All @@ -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]
Expand Down
84 changes: 76 additions & 8 deletions crates/ppvm-tableau/src/noise.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -608,15 +608,83 @@ 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<QubitStatus>, 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 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
// 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 ===
Expand Down
48 changes: 38 additions & 10 deletions crates/ppvm-tableau/src/tableau_like.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand All @@ -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
Expand All @@ -30,7 +31,22 @@ fn is_probability<C: PartialOrd<f64>>(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.
///
/// `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<C: Coefficient>(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.
Expand All @@ -50,7 +66,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
Expand Down Expand Up @@ -104,12 +120,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);
Comment thread
Copilot marked this conversation as resolved.
Resolved
return;
}

let r = self.rng_mut().random::<f64>();
let sum = Self::Coeff::zero();
let idx = p
Expand Down Expand Up @@ -154,12 +185,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);
Expand Down
Loading