Add Clifford-to-transvection decomposition from arXiv:2102.11380

paulimer/bindings/python/paulimer.pyi

@@ -582,6 +582,20 @@ class CliffordUnitary:
         """Get the symplectic matrix representation."""
         ...
 
+    def to_transvections(self) -> list[SparsePauli]:
+        """Decompose into an ordered product of Clifford transvections (pi/4 Pauli exponents).
+
+        Returns Pauli operators ``[P_1, ..., P_k]`` such that applying ``exp(i pi/4 P_1)``, then
+        ``exp(i pi/4 P_2)``, ..., then ``exp(i pi/4 P_k)`` reproduces this Clifford's symplectic
+        (conjugation) action, using a linear number of factors. Pauli-image signs and the global
+        phase are not reproduced.
+        """
+        ...
+
+    def centralizer(self) -> list[SparsePauli]:
+        """Generators of the centralizer: Paulis fixed up to sign under conjugation."""
+        ...
+
     def qubits(self) -> slice:
         """Return a slice representing the qubit indices."""
         ...

paulimer/bindings/python/src/py_clifford.rs

@@ -1,7 +1,7 @@
 use derive_more::{Deref, DerefMut, From, Into};
 use paulimer::clifford::{
-    group_encoding_clifford_of, split_phased_css, split_qubit_cliffords_and_css, Clifford, CliffordMutable,
-    CliffordUnitary, XOrZ,
+    clifford_centralizer, clifford_to_transvections, group_encoding_clifford_of, split_phased_css,
+    split_qubit_cliffords_and_css, Clifford, CliffordMutable, CliffordUnitary, XOrZ,
 };
 use paulimer::pauli::{as_sparse, DensePauli, SparsePauli};
 use pyo3::exceptions::PyValueError;
@@ -290,6 +290,29 @@ impl PyCliffordUnitary {
         self.inner.symplectic_matrix().into()
     }
 
+    /// Decomposes this Clifford into an ordered product of Clifford transvections (pi/4 Pauli
+    /// exponents), reproducing its symplectic action with a linear number of factors.
+    ///
+    /// Returns Pauli operators ``[P_1, ..., P_k]`` such that applying ``exp(i pi/4 P_1)``, then
+    /// ``exp(i pi/4 P_2)``, ..., then ``exp(i pi/4 P_k)`` reproduces the conjugation action of this
+    /// Clifford. Pauli-image signs and the global phase are not reproduced; see
+    /// :meth:`to_pauli_exponents` for the sign-exact (but ``O(n^2)``) decomposition.
+    fn to_transvections(&self) -> Vec<PySparsePauli> {
+        clifford_to_transvections(&self.inner)
+            .into_iter()
+            .map(PySparsePauli::from)
+            .collect()
+    }
+
+    /// Returns generators of this Clifford's centralizer: the Pauli operators fixed up to sign under
+    /// conjugation (``clifford * P * clifford_dagger == +/- P``).
+    fn centralizer(&self) -> Vec<PySparsePauli> {
+        clifford_centralizer(&self.inner)
+            .into_iter()
+            .map(PySparsePauli::from)
+            .collect()
+    }
+
     #[allow(clippy::needless_pass_by_value)]
     fn left_mul(&mut self, unitary_op: PyUnitaryOp, support: Vec<usize>) {
         self.inner.left_mul(unitary_op.into(), &support);

paulimer/bindings/python/tests/transvection_test.py

@@ -0,0 +1,120 @@
+"""Tests for the Clifford -> transvection decomposition bindings (arXiv:2102.11380).
+
+The decomposition reproduces a Clifford's symplectic (conjugation) action with a linear number of
+pi/4 Pauli exponents, ignoring Pauli-image signs and the global phase.
+"""
+
+from hypothesis import given, settings
+from hypothesis import strategies as st
+
+from paulimer import CliffordUnitary, SparsePauli, UnitaryOpcode
+
+
+def _rebuild_from_transvections(transvections, qubit_count):
+    rebuilt = CliffordUnitary.identity(qubit_count)
+    for pauli in transvections:
+        rebuilt.left_mul_pauli_exp(pauli)
+    return rebuilt
+
+
+def _residue_rank(clifford):
+    return 2 * clifford.qubit_count - len(clifford.centralizer())
+
+
+def _is_conjugation_fixed(clifford, pauli):
+    image = SparsePauli.from_dense(clifford.image_of(pauli))
+    return (image * pauli).weight == 0
+
+
+def _assert_valid_decomposition(clifford):
+    qubit_count = clifford.qubit_count
+    transvections = clifford.to_transvections()
+
+    rebuilt = _rebuild_from_transvections(transvections, qubit_count)
+    assert rebuilt.symplectic_matrix == clifford.symplectic_matrix
+
+    for pauli in transvections:
+        assert pauli.weight > 0
+
+    minimum = _residue_rank(clifford)
+    assert len(transvections) >= minimum
+    assert len(transvections) <= 4 * qubit_count + 2
+
+
+def test_identity_has_no_transvections():
+    for qubit_count in range(5):
+        identity = CliffordUnitary.identity(qubit_count)
+        assert identity.to_transvections() == []
+        assert len(identity.centralizer()) == 2 * qubit_count
+
+
+def test_single_qubit_gate_lengths():
+    s_gate = CliffordUnitary.from_name("SqrtZ", [0], 1)
+    _assert_valid_decomposition(s_gate)
+    assert len(s_gate.to_transvections()) == 1
+
+    hadamard = CliffordUnitary.from_name("Hadamard", [0], 1)
+    _assert_valid_decomposition(hadamard)
+    assert len(hadamard.to_transvections()) == 1
+
+
+def test_swap_hyperbolic_branch():
+    swap = CliffordUnitary.from_name("Swap", [0, 1], 2)
+    _assert_valid_decomposition(swap)
+    assert _residue_rank(swap) == 2
+    assert len(swap.to_transvections()) == 3
+    assert len(swap.centralizer()) == 2
+
+
+def test_two_qubit_gates():
+    for name in ("ControlledX", "ControlledZ"):
+        clifford = CliffordUnitary.from_name(name, [0, 1], 2)
+        _assert_valid_decomposition(clifford)
+
+
+def test_centralizer_generators_are_conjugation_fixed():
+    clifford = CliffordUnitary.identity(3)
+    clifford.left_mul(UnitaryOpcode.Hadamard, [0])
+    clifford.left_mul(UnitaryOpcode.ControlledX, [0, 1])
+    clifford.left_mul(UnitaryOpcode.SqrtZ, [2])
+    centralizer = clifford.centralizer()
+    assert all(_is_conjugation_fixed(clifford, pauli) for pauli in centralizer)
+    assert all(pauli.weight > 0 for pauli in centralizer)
+
+
+_SINGLE_QUBIT_GATES = ["Hadamard", "SqrtZ", "SqrtX", "X", "Y", "Z"]
+_TWO_QUBIT_GATES = ["ControlledX", "ControlledZ", "Swap"]
+
+
+@st.composite
+def _random_clifford(draw, max_qubits=5):
+    qubit_count = draw(st.integers(min_value=1, max_value=max_qubits))
+    gate_count = draw(st.integers(min_value=0, max_value=3 * qubit_count))
+    clifford = CliffordUnitary.identity(qubit_count)
+    for _ in range(gate_count):
+        if qubit_count >= 2 and draw(st.booleans()):
+            name = draw(st.sampled_from(_TWO_QUBIT_GATES))
+            first = draw(st.integers(min_value=0, max_value=qubit_count - 1))
+            second = draw(
+                st.integers(min_value=0, max_value=qubit_count - 1).filter(lambda q: q != first)
+            )
+            clifford.left_mul(getattr(UnitaryOpcode, name), [first, second])
+        else:
+            name = draw(st.sampled_from(_SINGLE_QUBIT_GATES))
+            qubit = draw(st.integers(min_value=0, max_value=qubit_count - 1))
+            clifford.left_mul(getattr(UnitaryOpcode, name), [qubit])
+    return clifford
+
+
+@settings(max_examples=200)
+@given(_random_clifford())
+def test_random_cliffords_reproduce_symplectic_action(clifford):
+    _assert_valid_decomposition(clifford)
+
+
+@settings(max_examples=200)
+@given(_random_clifford())
+def test_random_centralizers_are_conjugation_fixed(clifford):
+    for generator in clifford.centralizer():
+        assert _is_conjugation_fixed(clifford, generator)
+        assert generator.weight > 0

paulimer/src/clifford.rs

@@ -308,5 +308,8 @@ pub use clifford_impl::{
     z_images_partition_transform,
 };
 
+mod transvection;
+pub use transvection::{clifford_centralizer, clifford_to_transvections};
+
 #[derive(Debug, PartialEq, Eq, Default)]
 pub struct CliffordStringParsingError;

paulimer/src/clifford/transvection.rs

@@ -0,0 +1,206 @@
+//! Minimal decomposition of Clifford unitaries into Clifford transvections (`π/4` Pauli exponents).
+//!
+//! A *Clifford transvection* is the `π/4` Pauli exponent `exp(iπ/4·P_v)`, whose conjugation action
+//! on Pauli operators is the *symplectic transvection*
+//!
+//! ```text
+//! x ↦ x + ⟨x, v⟩ v,
+//! ```
+//!
+//! where `⟨·,·⟩` is the symplectic (commutation) form. This module follows the transvection
+//! framework of [arXiv:2102.11380](https://arxiv.org/abs/2102.11380) (Pllaha, Volanto & Tirkkonen,
+//! *Decomposition of Clifford Gates*): every Clifford is a product of transvections, and the
+//! *minimal* number of factors is `r = 2n − dim Fix(F)` (or `r + 1` when the symplectic action `F`
+//! is hyperbolic), where `Fix(F)` is the space of Pauli operators fixed by conjugation.
+//!
+//! The decomposition here uses a greedy O'Meara-style reduction: it always produces a **linear
+//! number of factors** (`O(n)`), reproducing the symplectic action exactly, but it is **not
+//! guaranteed to hit the strict `r`/`r + 1` minimum** — intermediate maps can become hyperbolic,
+//! adding an occasional extra factor. In practice it stays within a small additive constant of the
+//! minimum. The strict-minimum variant (via the paper's congruence-triangulation machinery) is
+//! tracked as a follow-up.
+//!
+//! Unlike [`clifford_to_pauli_exponents`](super::clifford_to_pauli_exponents), which reproduces the
+//! full signed tableau (and hence an exact global phase when replayed on a phased operator), this
+//! decomposition reproduces only the **symplectic action** — it ignores Pauli-image signs and the
+//! global phase. Its advantage is the linear factor count `O(n)`, versus `O(n²)` for the
+//! Gaussian-elimination decomposition.
+
+use binar::matrix::{AlignedBitMatrix, kernel_basis_matrix};
+use binar::{Bitwise, IndexSet};
+
+use crate::clifford::{Clifford, CliffordMutable, CliffordUnitary};
+use crate::pauli::DensePauli;
+use crate::{Pauli, PauliBinaryOps, PauliMutable, SparsePauli, anti_commutes_with};
+
+/// Decomposes `clifford` into an ordered product of Clifford transvections.
+///
+/// Returns a list of Pauli operators `[P₁, …, P_k]` such that left-multiplying the identity by the
+/// transvections `exp(iπ/4·P₁)`, then `exp(iπ/4·P₂)`, …, then `exp(iπ/4·P_k)` reproduces the
+/// **symplectic action** of `clifford` (its conjugation map on Pauli operators). The Pauli-image
+/// signs and the global phase are *not* reproduced; see the module docs for the contrast with
+/// [`clifford_to_pauli_exponents`](super::clifford_to_pauli_exponents).
+///
+/// The number of factors is **linear** in the qubit count (`O(n)`). It is close to, but not
+/// guaranteed to equal, the strict minimum `r = 2n − dim Fix(clifford)` (`r + 1` when the
+/// symplectic action is hyperbolic) of [arXiv:2102.11380](https://arxiv.org/abs/2102.11380); the
+/// greedy reduction here can add an occasional extra factor when an intermediate map becomes
+/// hyperbolic. The count is always at least `r`.
+///
+/// Every factor is returned with phase exponent `0`; the sign of a transvection does not affect its
+/// symplectic action, so `exp(iπ/4·P)` and `exp(−iπ/4·P)` are interchangeable here.
+///
+/// # Examples
+///
+/// ```
+/// use paulimer::CliffordUnitary;
+/// use paulimer::clifford::{clifford_to_transvections, Clifford, CliffordMutable};
+///
+/// let mut clifford = CliffordUnitary::identity(2);
+/// clifford.left_mul_hadamard(0);
+/// clifford.left_mul_cx(0, 1);
+///
+/// let transvections = clifford_to_transvections(&clifford);
+///
+/// let mut rebuilt = CliffordUnitary::identity(2);
+/// for pauli in &transvections {
+///     rebuilt.left_mul_pauli_exp(pauli);
+/// }
+/// // The symplectic actions agree (signs and global phase may differ).
+/// assert_eq!(rebuilt.symplectic_matrix(), clifford.symplectic_matrix());
+/// ```
+#[must_use]
+pub fn clifford_to_transvections(clifford: &CliffordUnitary) -> Vec<SparsePauli> {
+    let qubit_count = clifford.num_qubits();
+    let mut working = clifford.clone();
+    let mut recorded = Vec::new();
+    // Reduce the symplectic action to the identity by left-multiplying transvections `T_{v₁}, …,
+    // T_{v_k}`, so that `T_{v_k} ⋯ T_{v₁} · F = I` and hence `F = T_{v₁} ⋯ T_{v_k}`. Replaying the
+    // factors in reverse order rebuilds `F` from the identity.
+    while let Some(transvection) = next_transvection(&working) {
+        working.left_mul_pauli_exp(&transvection);
+        recorded.push(transvection);
+        debug_assert!(
+            recorded.len() <= 4 * qubit_count + 2,
+            "transvection reduction exceeded its linear termination bound"
+        );
+    }
+    recorded.reverse();
+    recorded
+}
+
+/// Returns generators of the Clifford's centralizer: the Pauli operators fixed (up to sign) by
+/// conjugation, i.e. the `P` with `clifford · P · clifford† = ±P`.
+///
+/// This is `Fix(F)`, the kernel of the residue map `P ↦ conj(P) · P`, computed as the left null
+/// space of the residue matrix over GF(2). The returned Paulis are independent generators (with
+/// phase exponent `0`); the centralizer they span has dimension `dim Fix(F) = 2n − r`, where `r` is
+/// the number of factors returned by [`clifford_to_transvections`] for a non-hyperbolic action.
+///
+/// # Examples
+///
+/// ```
+/// use paulimer::{CliffordUnitary, Pauli};
+/// use paulimer::clifford::{clifford_centralizer, Clifford, CliffordMutable};
+///
+/// let mut clifford = CliffordUnitary::identity(1);
+/// clifford.left_mul_root_z(0); // S fixes Z, sends X -> Y
+///
+/// let generators = clifford_centralizer(&clifford);
+/// // Every generator is fixed (up to sign) under conjugation.
+/// assert!(generators.iter().all(|pauli| {
+///     let image = clifford.image(pauli);
+///     image.x_bits() == pauli.x_bits() && image.z_bits() == pauli.z_bits()
+/// }));
+/// ```
+#[must_use]
+pub fn clifford_centralizer(clifford: &CliffordUnitary) -> Vec<SparsePauli> {
+    let qubit_count = clifford.num_qubits();
+    let dimension = 2 * qubit_count;
+    let mut residue = AlignedBitMatrix::zeros(dimension, dimension);
+    for (row, basis) in symplectic_basis(qubit_count).enumerate() {
+        let vector = residue_vector(&basis, &clifford.image(&basis));
+        for qubit in 0..qubit_count {
+            if vector.x_bits().index(qubit) {
+                residue.set((row, qubit), true);
+            }
+            if vector.z_bits().index(qubit) {
+                residue.set((row, qubit_count + qubit), true);
+            }
+        }
+    }
+    let kernel = kernel_basis_matrix(&residue.transposed());
+    (0..kernel.row_count())
+        .map(|row| {
+            let x_bits: IndexSet = (0..qubit_count).filter(|&qubit| kernel[(row, qubit)]).collect();
+            let z_bits: IndexSet = (0..qubit_count)
+                .filter(|&qubit| kernel[(row, qubit_count + qubit)])
+                .collect();
+            SparsePauli::from_bits(x_bits, z_bits, 0)
+        })
+        .collect()
+}
+
+/// The `2n` standard basis Pauli operators `X₀, …, X_{n−1}, Z₀, …, Z_{n−1}`.
+fn symplectic_basis(qubit_count: usize) -> impl Iterator<Item = SparsePauli> {
+    (0..qubit_count)
+        .map(move |qubit| SparsePauli::x(qubit, qubit_count))
+        .chain((0..qubit_count).map(move |qubit| SparsePauli::z(qubit, qubit_count)))
+}
+
+/// The next transvection `T_v` reducing the residue of `working`, or `None` if `working` already
+/// acts as the identity on Pauli operators (up to sign).
+///
+/// Following the O'Meara strategy of [arXiv:2102.11380](https://arxiv.org/abs/2102.11380): find a
+/// vector `x` with `⟨x, conj(x)⟩ = 1` (`x` anticommutes with its own image) and set `v = x + conj(x)`
+/// — a residue vector — which lowers the residue rank by one. If no such `x` exists but `working`
+/// is non-trivial (the hyperbolic case), any nonzero residue vector `v` makes the action
+/// non-hyperbolic while preserving the residue space, costing one extra transvection.
+fn next_transvection(working: &CliffordUnitary) -> Option<SparsePauli> {
+    let qubit_count = working.num_qubits();
+    let basis: Vec<SparsePauli> = symplectic_basis(qubit_count).collect();
+    let images: Vec<DensePauli> = basis.iter().map(|pauli| working.image(pauli)).collect();
+
+    for (pauli, image) in basis.iter().zip(&images) {
+        if anti_commutes_with(pauli, image) {
+            return Some(residue_vector(pauli, image));
+        }
+    }
+
+    let dimension = basis.len();
+    for first in 0..dimension {
+        for second in (first + 1)..dimension {
+            let anticommuting =
+                anti_commutes_with(&basis[first], &images[second]) ^ anti_commutes_with(&basis[second], &images[first]);
+            if anticommuting {
+                let mut sum = basis[first].clone();
+                sum.mul_assign_left(&basis[second]);
+                let mut image = images[first].clone();
+                image.mul_assign_left(&images[second]);
+                return Some(residue_vector(&sum, &image));
+            }
+        }
+    }
+
+    basis
+        .iter()
+        .zip(&images)
+        .find(|(pauli, image)| !acts_trivially_on(pauli, image))
+        .map(|(pauli, image)| residue_vector(pauli, image))
+}
+
+/// The residue vector `v = x + conj(x)` as a phaseless Pauli (its symplectic vector is the product
+/// `x · conj(x)`).
+fn residue_vector(pauli: &SparsePauli, image: &DensePauli) -> SparsePauli {
+    let mut vector: SparsePauli = image.clone().into();
+    vector.mul_assign_left(pauli);
+    vector.assign_phase_exp(0);
+    vector
+}
+
+/// Whether `image` equals `pauli` as a symplectic vector (i.e. conjugation fixes `pauli` up to sign).
+fn acts_trivially_on(pauli: &SparsePauli, image: &DensePauli) -> bool {
+    let mut difference: SparsePauli = image.clone().into();
+    difference.mul_assign_left(pauli);
+    difference.x_bits().is_zero() && difference.z_bits().is_zero()
+}