diff --git a/crates/circuit/src/util.rs b/crates/circuit/src/util.rs index 76840b551a7f..112044bf8cfd 100644 --- a/crates/circuit/src/util.rs +++ b/crates/circuit/src/util.rs @@ -47,3 +47,12 @@ pub const C_ONE: Complex64 = c64(1., 0.); pub const C_M_ONE: Complex64 = c64(-1., 0.); pub const IM: Complex64 = c64(0., 1.); pub const M_IM: Complex64 = c64(0., -1.); + +use std::f64::consts::{FRAC_PI_2, FRAC_PI_4, FRAC_PI_8, SQRT_2}; +pub const C_FRAC_PI_2: Complex64 = c64(FRAC_PI_2, 0.0); +pub const C_FRAC_PI_4: Complex64 = c64(FRAC_PI_4, 0.0); +pub const C_FRAC_PI_8: Complex64 = c64(FRAC_PI_8, 0.0); +pub const C_FRAC_PI_2_SQRT_2: Complex64 = c64(FRAC_PI_2 / SQRT_2, 0.0); +pub const C_M_FRAC_PI_4: Complex64 = c64(-FRAC_PI_4, 0.0); +pub const C_M_FRAC_PI_8: Complex64 = c64(-FRAC_PI_8, 0.0); +pub const C_M_FRAC_PI_2_SQRT_2: Complex64 = c64(-FRAC_PI_2 / SQRT_2, 0.0); diff --git a/crates/pyext/src/lib.rs b/crates/pyext/src/lib.rs index 27b8c64c2124..7cb7fcca01ee 100644 --- a/crates/pyext/src/lib.rs +++ b/crates/pyext/src/lib.rs @@ -68,6 +68,7 @@ fn _accelerate(m: &Bound) -> PyResult<()> { add_submodule(m, ::qiskit_transpiler::passes::sabre::sabre, "sabre")?; add_submodule(m, ::qiskit_accelerate::sampled_exp_val::sampled_exp_val, "sampled_exp_val")?; add_submodule(m, ::qiskit_quantum_info::sparse_observable::sparse_observable, "sparse_observable")?; + add_submodule(m, ::qiskit_quantum_info::sparse_observable::standard_generators::standard_generators, "standard_generators")?; add_submodule(m, ::qiskit_quantum_info::sparse_pauli_op::sparse_pauli_op, "sparse_pauli_op")?; add_submodule(m, ::qiskit_quantum_info::unitary_sim::unitary_sim, "unitary_sim")?; add_submodule(m, ::qiskit_transpiler::passes::split_2q_unitaries_mod, "split_2q_unitaries")?; diff --git a/crates/quantum_info/src/sparse_observable/mod.rs b/crates/quantum_info/src/sparse_observable/mod.rs index fcd4f113c7f4..77f94842037d 100644 --- a/crates/quantum_info/src/sparse_observable/mod.rs +++ b/crates/quantum_info/src/sparse_observable/mod.rs @@ -10,6 +10,8 @@ // copyright notice, and modified files need to carry a notice indicating // that they have been altered from the originals. +pub mod standard_generators; + mod lookup; use hashbrown::HashSet; @@ -4261,8 +4263,19 @@ fn coerce_to_observable<'py>( } } } +#[pyfunction(name = "_generator_observable")] +#[pyo3(signature = (gate, params = None))] +pub fn generator_observable_py( + gate: qiskit_circuit::operations::StandardGate, + params: Option>, +) -> Option { + let params = params.unwrap_or_default(); + standard_generators::generator_observable(gate, ¶ms) +} + pub fn sparse_observable(m: &Bound) -> PyResult<()> { m.add_class::()?; + m.add_function(wrap_pyfunction!(generator_observable_py, m)?)?; Ok(()) } diff --git a/crates/quantum_info/src/sparse_observable/standard_generators.rs b/crates/quantum_info/src/sparse_observable/standard_generators.rs new file mode 100644 index 000000000000..a8c639bc36fe --- /dev/null +++ b/crates/quantum_info/src/sparse_observable/standard_generators.rs @@ -0,0 +1,459 @@ +// This code is part of Qiskit. +// +// (C) Copyright IBM 2026 +// +// This code is licensed under the Apache License, Version 2.0. You may +// obtain a copy of this license in the LICENSE.txt file in the root directory +// of this source tree or at https://www.apache.org/licenses/LICENSE-2.0. +// +// Any modifications or derivative works of this code must retain this +// copyright notice, and modified files need to carry a notice indicating +// that they have been altered from the originals. +// src/sparse_observable/standard_generators.rs +// +// This module maps standard quantum gates to their Hamiltonian generators H_gate such that: +// +// gate = exp(-i * H_gate) (up to global phase) +// +// The generator is returned as a SparseObservable (a sum of Pauli tensor products). +// +// For single-qubit gates: X = exp(-i*(pi/2)*X), H = exp(-i*(pi/2)*(X+Z)/sqrt(2)), etc. +// For multi-qubit controlled gates: CX = exp(-i*(pi/4)*(ZX - ZI - IX)), etc. +// For rotation gates: RX(t) = exp(-i*(t/2)*X), RY(t) = exp(-i*(t/2)*Y), etc. +// +// The SoA (Struct-of-Arrays) layout is used: `bit_terms`, `indices`, and `boundaries` are +// flattened arrays. Term `i` uses bit_terms[boundaries[i]..boundaries[i+1]] and +// indices[boundaries[i]..boundaries[i+1]] for its Pauli operators and qubit targets. + +use super::BitTerm; +use super::SparseObservable; + +use pyo3::prelude::*; +use pyo3::wrap_pyfunction; +use qiskit_circuit::operations::{Operation, Param, StandardGate}; +use qiskit_circuit::util::{ + C_FRAC_PI_2, C_FRAC_PI_2_SQRT_2, C_FRAC_PI_4, C_FRAC_PI_8, C_M_FRAC_PI_2_SQRT_2, C_M_FRAC_PI_4, + C_M_FRAC_PI_8, C_ZERO, c64, +}; + +const BETA_TOLERANCE: f64 = 1e-10; + +/// For parametric gates (e.g., `RX(theta)`), the generator $H$ depends on the +/// gate parameters (e.g., $H = (\theta/2)X$). This function extracts parameter values +/// from the `params` slice to compute the concrete coefficients for the returned +/// `SparseObservable`. +/// +/// If parameters are missing or symbolic (non-`Float`), it defaults to a coefficient +/// of 1.0. This allows the commutation checker to still prove commutation in cases +/// where the generator's Pauli structure alone is sufficient (e.g., `[theta*X, X] = 0` +/// for any `theta`), but it avoids attempting to store parametric expressions, which +/// `SparseObservable` does not currently support. +pub fn generator_observable(gate: StandardGate, params: &[Param]) -> Option { + let _params = params; + let num_qubits = gate.num_qubits(); + + use BitTerm::*; + + // Global phase gate + if num_qubits == 0 { + if let StandardGate::GlobalPhase = gate { + // Global Phase is exp(-i * theta * I) wait... + // Qiskit GlobalPhaseGate(theta) matrix is exp(i * theta). + // Generator H such that exp(-i * H) = exp(i * theta) -> H = -theta. + // A 0-qubit Identity operator has 1 term (the empty string). + let mut theta = 1.0; + if let [Param::Float(t)] = _params { + theta = *t; + } + return Some( + SparseObservable::new( + 0, + vec![c64(-theta, 0.0)], + vec![], // no paulis -> Identity + vec![], // no target qubits + vec![0, 0], // 1 term of length 0 + ) + .expect("invalid 0-qubit generator layout"), + ); + } + return None; + } + + // Single-qubit gates + if num_qubits == 1 { + let (coeffs, terms, indices) = match gate { + // H = exp(-i*(pi/2)*(X + Z)/sqrt(2)) + // => H_gen = (pi/2)*(X + Z)/sqrt(2) = (pi / (2*sqrt(2))) * X + (pi / (2*sqrt(2))) * Z + // Numerically: pi / (2*sqrt(2)) ≈ 1.1107... + // Note: the sign must be negative so H_gate uses coefficients +1/sqrt(2) each. + StandardGate::H => { + let c = C_FRAC_PI_2_SQRT_2.re; + (vec![c64(c, 0.0), c64(c, 0.0)], vec![X, Z], vec![0, 0]) + } + // X = exp(-i*(pi/2)*X), Y = exp(-i*(pi/2)*Y), Z = exp(-i*(pi/2)*Z) + StandardGate::X => (vec![C_FRAC_PI_2], vec![X], vec![0]), + StandardGate::Y => (vec![C_FRAC_PI_2], vec![Y], vec![0]), + StandardGate::Z => (vec![C_FRAC_PI_2], vec![Z], vec![0]), + // Identity: exp(-i * 0) = I. + StandardGate::I => (vec![C_ZERO], vec![], vec![]), + // S = exp(-i*(pi/4)*Z), Sdg = exp(-i*(-pi/4)*Z) + StandardGate::S => (vec![C_FRAC_PI_4], vec![Z], vec![0]), + StandardGate::Sdg => (vec![C_M_FRAC_PI_4], vec![Z], vec![0]), + // T = exp(-i*(pi/8)*Z), Tdg = exp(-i*(-pi/8)*Z) + StandardGate::T => (vec![C_FRAC_PI_8], vec![Z], vec![0]), + StandardGate::Tdg => (vec![C_M_FRAC_PI_8], vec![Z], vec![0]), + // SX = exp(-i*(pi/4)*X), SXdg = exp(-i*(-pi/4)*X) + StandardGate::SX => (vec![C_FRAC_PI_4], vec![X], vec![0]), + StandardGate::SXdg => (vec![C_M_FRAC_PI_4], vec![X], vec![0]), + // RX(t) = exp(-i*(t/2)*X), RY(t) = exp(-i*(t/2)*Y) + // RZ(t) = exp(-i*(t/2)*Z), Phase(t) = exp(-i*(t/2)*Z) (same generator) + StandardGate::RX | StandardGate::RY | StandardGate::RZ | StandardGate::Phase => { + // Qiskit's `_generator_observable` falls back to `1.0` if no parameters are available or the parameter is an unbound expression. + // This corresponds effectively to returning the base operator for the Pauli (e.g. `X.generator() == X`). + // In normal workflows the `params` tuple is fully concrete during commutation logic (i.e., `Float`). + let theta = if let [Param::Float(t)] = _params { + *t + } else { + 1.0 + }; + let term = match gate { + StandardGate::RX => X, + StandardGate::RY => Y, + _ => Z, + }; + (vec![c64(theta / 2.0, 0.0)], vec![term], vec![0]) + } + _ => return None, + }; + + // For 1-qubit gates each term has exactly 1 Pauli operator, so + // boundaries = [0, 1, 2, ..., N]. + // Exception: Identity gate has 1 term of length 0. + let boundaries: Vec = if gate == StandardGate::I { + vec![0, 0] + } else { + (0..=terms.len()).collect() + }; + + return Some( + SparseObservable::new(num_qubits, coeffs, terms, indices, boundaries) + .expect("invalid 1-qubit generator layout"), + ); + } + + // Multi-qubit gates + // Returns (coeffs, bit_terms, indices, boundaries) in SoA layout. + // Qubit ordering convention: index 0 = qubit 0 = LEAST significant (rightmost in Qiskit strings). + // For a 2q gate: control = qubit 0, target = qubit 1. + let (coeffs, bit_terms, indices, boundaries) = match gate { + // CX (CNOT): CX = exp(-i*(pi/4)*(Z0*X1 - Z0 - X1)) + // Generator H = (pi/4)*(Z0*X1 - Z0 - X1) + // Terms: [Z0X1 coeff=+pi/4], [Z0 coeff=-pi/4], [X1 coeff=-pi/4] + StandardGate::CX => ( + vec![C_M_FRAC_PI_4, C_FRAC_PI_4, C_FRAC_PI_4], + // Term 0: Z(q0) X(q1); Term 1: Z(q0); Term 2: X(q1) + vec![Z, X, Z, X], + vec![0, 1, 0, 1], + vec![0, 2, 3, 4], + ), + // CY: CY = exp(-i*(pi/4)*(Z0*Y1 - Z0 - Y1)) + StandardGate::CY => ( + vec![C_M_FRAC_PI_4, C_FRAC_PI_4, C_FRAC_PI_4], + vec![Z, Y, Z, Y], + vec![0, 1, 0, 1], + vec![0, 2, 3, 4], + ), + // CZ: CZ = exp(-i*(pi/4)*(Z0*Z1 - Z0 - Z1)) + StandardGate::CZ => ( + vec![C_M_FRAC_PI_4, C_FRAC_PI_4, C_FRAC_PI_4], + vec![Z, Z, Z, Z], + vec![0, 1, 0, 1], + vec![0, 2, 3, 4], + ), + // CS (sqrt(CZ)): CS = exp(-i*(pi/8)*(Z0 + Z1 - Z0*Z1)) + // Generator H = (pi/8)*(Z0 + Z1 - Z0*Z1) + // Wait, Shelly prefers factor -pi/8, so H = -pi/8 * (ZZ - ZI - IZ) = -pi/8*ZZ + pi/8*ZI + pi/8*IZ + // That matches exp(-i * (-pi/8) * (Z0 + Z1 - Z0*Z1)). + StandardGate::CS => ( + vec![C_FRAC_PI_8, C_FRAC_PI_8, C_M_FRAC_PI_8], + vec![Z, Z, Z, Z], + vec![0, 1, 0, 1], + vec![0, 1, 2, 4], + ), + // CSdg (inv sqrt(CZ)): factor +pi/8 + StandardGate::CSdg => ( + vec![C_M_FRAC_PI_8, C_M_FRAC_PI_8, C_FRAC_PI_8], + vec![Z, Z, Z, Z], + vec![0, 1, 0, 1], + vec![0, 1, 2, 4], + ), + // CSX (sqrt(CX)/controlled-SX): factor -pi/8 + StandardGate::CSX => ( + vec![C_FRAC_PI_8, C_FRAC_PI_8, C_M_FRAC_PI_8], + vec![Z, X, Z, X], + vec![0, 1, 0, 1], + vec![0, 1, 2, 4], + ), + // CRX(t): CRX(t) = exp(-i*(t/4)*(-Z0*X1 + X1)) + // Generator = -t/4 * Z0*X1 + t/4 * X1 + StandardGate::CRX => { + let t = if let [Param::Float(theta)] = _params { + *theta + } else { + 1.0 + }; + ( + vec![c64(-t / 4.0, 0.0), c64(t / 4.0, 0.0)], + vec![Z, X, X], + vec![0, 1, 1], + vec![0, 2, 3], + ) + } + // CRY(t): Generator = -t/4 * Z0*Y1 + t/4 * Y1 + StandardGate::CRY => { + let t = if let [Param::Float(theta)] = _params { + *theta + } else { + 1.0 + }; + ( + vec![c64(-t / 4.0, 0.0), c64(t / 4.0, 0.0)], + vec![Z, Y, Y], + vec![0, 1, 1], + vec![0, 2, 3], + ) + } + // CRZ(t): Generator = -t/4 * Z0*Z1 + t/4 * Z1 + StandardGate::CRZ => { + let t = if let [Param::Float(theta)] = _params { + *theta + } else { + 1.0 + }; + ( + vec![c64(-t / 4.0, 0.0), c64(t / 4.0, 0.0)], + vec![Z, Z, Z], + vec![0, 1, 1], + vec![0, 2, 3], + ) + } + // CPhase(t): Generator = -t/4 * Z0*Z1 + t/4 * Z0 + t/4 * Z1 + // (same as CRZ but shifts both Z0 and Z1, not just Z1) + StandardGate::CPhase => { + let t = if let [Param::Float(theta)] = _params { + *theta + } else { + 1.0 + }; + ( + vec![c64(-t / 4.0, 0.0), c64(t / 4.0, 0.0), c64(t / 4.0, 0.0)], + vec![Z, Z, Z, Z], + vec![0, 1, 0, 1], + vec![0, 2, 3, 4], + ) + } + // Swap: Swap = exp(-i*(pi/4)*(X0*X1 + Y0*Y1 + Z0*Z1)) + // (note: Swap = exp(-i*pi/4*(XX+YY+ZZ)) treats Swap as "swap up to phase for each sector") + StandardGate::Swap => ( + vec![C_FRAC_PI_4, C_FRAC_PI_4, C_FRAC_PI_4], + vec![X, X, Y, Y, Z, Z], + vec![0, 1, 0, 1, 0, 1], + vec![0, 2, 4, 6], + ), + // ISwap: ISwap = exp(-i*(-pi/4)*(X0*X1 + Y0*Y1)) + StandardGate::ISwap => ( + vec![C_M_FRAC_PI_4, C_M_FRAC_PI_4], + vec![X, X, Y, Y], + vec![0, 1, 0, 1], + vec![0, 2, 4], + ), + // RXX(t) = exp(-i*(t/2)*X0*X1) + StandardGate::RXX => { + let t = if let [Param::Float(theta)] = _params { + *theta + } else { + 1.0 + }; + (vec![c64(t / 2.0, 0.0)], vec![X, X], vec![0, 1], vec![0, 2]) + } + // RYY(t) = exp(-i*(t/2)*Y0*Y1) + StandardGate::RYY => { + let t = if let [Param::Float(theta)] = _params { + *theta + } else { + 1.0 + }; + (vec![c64(t / 2.0, 0.0)], vec![Y, Y], vec![0, 1], vec![0, 2]) + } + // RZZ(t) = exp(-i*(t/2)*Z0*Z1) + StandardGate::RZZ => { + let t = if let [Param::Float(theta)] = _params { + *theta + } else { + 1.0 + }; + (vec![c64(t / 2.0, 0.0)], vec![Z, Z], vec![0, 1], vec![0, 2]) + } + // RZX(t) = exp(-i*(t/2)*Z0*X1) + StandardGate::RZX => { + let t = if let [Param::Float(theta)] = _params { + *theta + } else { + 1.0 + }; + (vec![c64(t / 2.0, 0.0)], vec![Z, X], vec![0, 1], vec![0, 2]) + } + // XXPlusYY(theta, beta): Generator = (theta/4)*(X0*X1 + Y0*Y1) + // (the beta angle just rotates the YY axis; for the commutation check only XX+YY matters) + StandardGate::XXPlusYY | StandardGate::XXMinusYY => { + let t = if let [Param::Float(theta), ..] = _params { + *theta + } else { + 1.0 + }; + + // The beta angle rotates the YY axis. Ensure beta=0 (or assert it) so commutation + // is strictly XX +/- YY. Otherwise, fallback to matrix checking. + let beta = if let [_, Param::Float(b)] = _params { + *b + } else { + return None; + }; + + if beta.abs() > BETA_TOLERANCE { + return None; + } + + match gate { + StandardGate::XXPlusYY => ( + vec![c64(t / 4.0, 0.0), c64(t / 4.0, 0.0)], + vec![X, X, Y, Y], + vec![0, 1, 0, 1], + vec![0, 2, 4], + ), + StandardGate::XXMinusYY => ( + vec![c64(t / 4.0, 0.0), c64(-t / 4.0, 0.0)], + vec![X, X, Y, Y], + vec![0, 1, 0, 1], + vec![0, 2, 4], + ), + _ => unreachable!(), + } + } + // CCX (Toffoli): CCX = exp(-i*(pi/8)*(Z0*Z1*X2 - Z0*X2 - Z1*X2 - Z0*Z1 + Z0 + Z1 + X2)) + // 7 terms in total. + // qubit ordering: q0=ctrl0, q1=ctrl1, q2=target + StandardGate::CCX => ( + vec![ + C_FRAC_PI_8, + C_M_FRAC_PI_8, + C_M_FRAC_PI_8, + C_M_FRAC_PI_8, + C_FRAC_PI_8, + C_FRAC_PI_8, + C_FRAC_PI_8, + ], + vec![Z, Z, X, Z, X, Z, X, Z, Z, Z, Z, X], + vec![0, 1, 2, 0, 2, 1, 2, 0, 1, 0, 1, 2], + vec![0, 3, 5, 7, 9, 10, 11, 12], + ), + // CCZ: CCZ = exp(-i*(pi/8)*(Z0*Z1*Z2 - Z0*Z2 - Z1*Z2 - Z0*Z1 + Z0 + Z1 + Z2)) + StandardGate::CCZ => ( + vec![ + C_M_FRAC_PI_8, + C_FRAC_PI_8, + C_FRAC_PI_8, + C_FRAC_PI_8, + C_M_FRAC_PI_8, + C_M_FRAC_PI_8, + C_M_FRAC_PI_8, + ], + vec![Z, Z, Z, Z, Z, Z, Z, Z, Z, Z, Z, Z], + vec![0, 1, 2, 0, 2, 1, 2, 0, 1, 0, 1, 2], + vec![0, 3, 5, 7, 9, 10, 11, 12], + ), + // CSwap (Fredkin): CSwap = exp(-i*(pi/8)*(Z0 - Z0*X1*X2 - Z0*Y1*Y2 - Z0*Z1*Z2 + X1*X2 + Y1*Y2 + Z1*Z2)) + // 7 terms: Z0(-ZXX=ZYY=ZZZ), +XX, +YY, +ZZ + StandardGate::CSwap => ( + vec![ + C_FRAC_PI_8, + C_M_FRAC_PI_8, + C_M_FRAC_PI_8, + C_M_FRAC_PI_8, + C_FRAC_PI_8, + C_FRAC_PI_8, + C_FRAC_PI_8, + ], + vec![Z, Z, X, X, Z, Y, Y, Z, Z, Z, X, X, Y, Y, Z, Z], + vec![0, 0, 1, 2, 0, 1, 2, 0, 1, 2, 1, 2, 1, 2, 1, 2], + vec![0, 1, 4, 7, 10, 12, 14, 16], + ), + // ECR: ECR = exp(-i * (pi/2/sqrt(2)) * (IX - XY)) + // Terms: X1 with coeff +pi/(2*sqrt(2)), X1Y0 with coeff -pi/(2*sqrt(2)) + StandardGate::ECR => ( + vec![C_FRAC_PI_2_SQRT_2, C_M_FRAC_PI_2_SQRT_2], + vec![X, Y, X], + vec![0, 0, 1], + vec![0, 1, 3], + ), + _ => return None, + }; + + Some( + SparseObservable::new(num_qubits, coeffs, bit_terms, indices, boundaries) + .expect("invalid multi-qubit generator layout"), + ) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn rx_has_some_generator() { + let obs = generator_observable(StandardGate::RX, &[]).expect("Rx should have a generator"); + assert!(!obs.bit_terms().is_empty()); + } + + #[test] + fn cx_has_some_generator() { + let obs = generator_observable(StandardGate::CX, &[]).expect("CX should have a generator"); + assert_eq!(obs.num_terms(), 3); + } + + #[test] + fn ccz_has_seven_terms() { + let obs = + generator_observable(StandardGate::CCZ, &[]).expect("CCZ should have a generator"); + assert_eq!(obs.num_terms(), 7); + } + + #[test] + fn cswap_has_seven_terms() { + let obs = + generator_observable(StandardGate::CSwap, &[]).expect("CSwap should have a generator"); + assert_eq!(obs.num_terms(), 7); + } + + #[test] + fn ecr_has_two_terms() { + let obs = + generator_observable(StandardGate::ECR, &[]).expect("ECR should have a generator"); + assert_eq!(obs.num_terms(), 2); + } +} + +#[pyfunction(name = "generator_observable")] +#[pyo3(signature = (gate, params=None))] +pub fn generator_observable_py( + gate: StandardGate, + params: Option>, +) -> Option { + let params = params.unwrap_or_default(); + generator_observable(gate, ¶ms) +} + +pub fn standard_generators(m: &Bound) -> PyResult<()> { + m.add_function(wrap_pyfunction!(generator_observable_py, m)?)?; + Ok(()) +} diff --git a/crates/quantum_info/src/sparse_pauli_op.rs b/crates/quantum_info/src/sparse_pauli_op.rs index 1709cde418dd..f26140961356 100644 --- a/crates/quantum_info/src/sparse_pauli_op.rs +++ b/crates/quantum_info/src/sparse_pauli_op.rs @@ -1035,7 +1035,8 @@ pub fn to_matrix_sparse( // Pessimistic estimation of whether we can fit in `i32`. If there's any risk of overflowing // `i32`, we use `i64`, but Scipy will always try to downcast to `i32`, so we try to match it. - let max_entries_per_row = (paulis.num_ops() as u64).min(1u64 << (paulis.num_qubits() - 1)); + let max_entries_per_row = + (paulis.num_ops() as u64).min(1u64 << paulis.num_qubits().saturating_sub(1)); let use_32_bit = max_entries_per_row.saturating_mul(1u64 << paulis.num_qubits()) <= (i32::MAX as u64); if use_32_bit { diff --git a/crates/transpiler/src/commutation_checker.rs b/crates/transpiler/src/commutation_checker.rs index ea8685255c71..4bad2f784fc6 100644 --- a/crates/transpiler/src/commutation_checker.rs +++ b/crates/transpiler/src/commutation_checker.rs @@ -16,8 +16,8 @@ use ndarray::linalg::kron; use num_complex::Complex64; use num_complex::ComplexFloat; use qiskit_circuit::object_registry::PyObjectAsKey; -use qiskit_quantum_info::sparse_observable::PySparseObservable; -use qiskit_quantum_info::sparse_observable::SparseObservable; +use qiskit_quantum_info::sparse_observable::standard_generators::generator_observable; +use qiskit_quantum_info::sparse_observable::{PySparseObservable, SparseObservable}; use smallvec::SmallVec; use std::fmt::Debug; @@ -91,6 +91,14 @@ const fn build_supported_ops() -> [bool; STANDARD_GATE_SIZE] { lut[StandardGate::ISwap as usize] = true; lut[StandardGate::ECR as usize] = true; lut[StandardGate::CCX as usize] = true; + lut[StandardGate::CCZ as usize] = true; + lut[StandardGate::CS as usize] = true; + lut[StandardGate::CSdg as usize] = true; + lut[StandardGate::CSX as usize] = true; + lut[StandardGate::I as usize] = true; + lut[StandardGate::GlobalPhase as usize] = true; + lut[StandardGate::XXPlusYY as usize] = true; + lut[StandardGate::XXMinusYY as usize] = true; lut[StandardGate::CSwap as usize] = true; lut } @@ -236,9 +244,17 @@ fn try_extract_op_from_ppr( fn try_pauli_generator( operation: &OperationRef, + params: &[Param], qubits: &[Qubit], num_qubits: u32, ) -> Option { + if let OperationRef::StandardGate(gate) = operation { + if let Some(local) = generator_observable(*gate, params) { + let out = SparseObservable::identity(num_qubits); + return Some(out.compose_map(&local, |i| qubits[i as usize].0)); + } + } + match operation.name() { "pauli" => try_extract_op_from_pauli_gate(operation, qubits, num_qubits), "PauliEvolution" => try_extract_op_from_pauli_evo(operation, qubits, num_qubits), @@ -256,8 +272,7 @@ where T: From + Copy + Debug, u32: From, { - // Using `PyObjectAsKey` here is a total hack, but this is a short-term workaround before a - // larger refactor of the commutation checker. + // Use `PyObjectAsKey` as a workaround before a larger refactor of the commutation checker. let mut registry: ObjectRegistry = ObjectRegistry::new(); for bit in py_bits1.iter().chain(py_bits2.iter()) { @@ -278,6 +293,11 @@ where /// It handles the actual commutation checking, cache management, and library /// lookups. It's not meant to be a public facing Python object though and only used /// internally by the Python class. +/// +/// This implementation now supports efficient commutation checks for complex gates +/// (including multi-qubit gates like CCX, CSwap) by analyzing their generators using +/// the `sparse_observable::standard_generators` module. For efficiency, generator data +/// is stored and processed using a flattened Struct-of-Arrays (SoA) layout. #[pyclass(module = "qiskit._accelerate.commutation_checker")] pub struct CommutationChecker { library: CommutationLibrary, @@ -523,9 +543,9 @@ impl CommutationChecker { // Handle commutations in between Pauli-based gates, like PauliGate or PauliEvolutionGate let size = qargs1.iter().chain(qargs2.iter()).max().unwrap().0 + 1; - if let Some(obs1) = try_pauli_generator(op1, qargs1, size) { - if let Some(obs2) = try_pauli_generator(op2, qargs2, size) { - return Ok(obs1.commutes(&obs2, tol)); + if let Some(pauli1) = try_pauli_generator(op1, params1, qargs1, size) { + if let Some(pauli2) = try_pauli_generator(op2, params2, qargs2, size) { + return Ok(pauli1.commutes(&pauli2, tol)); } } diff --git a/releasenotes/notes/pauli-evolution-commutation-8f7caab8.yaml b/releasenotes/notes/pauli-evolution-commutation-8f7caab8.yaml new file mode 100644 index 000000000000..228316baa200 --- /dev/null +++ b/releasenotes/notes/pauli-evolution-commutation-8f7caab8.yaml @@ -0,0 +1,15 @@ +--- +features_transpiler: + - | + The :class:`.CommutationChecker` can now efficiently evaluate commutativity + between explicitly defined :class:`~qiskit.quantum_info.Pauli`, :class:`~qiskit.quantum_info.SparsePauliOp`, and + :class:`~qiskit.circuit.library.PauliEvolutionGate` operations with many other standard gates + (including multi-qubit gates like :class:`~qiskit.circuit.library.CCXGate`). This avoids converting + operators to their full matrix representation when finding commutation + relations, making commutation checking significantly faster and avoiding + out-of-memory errors for operations spanning large numbers of qubits. +fixes: + - | + Fixes an issue where the :class:`.LightCone` pass would fail with index errors or + out-of-memory panics when encountering large or complex operations. + See `#15021 `__ for context. diff --git a/test/python/circuit/test_commutation_checker.py b/test/python/circuit/test_commutation_checker.py index 65438b8bc5ec..95b85c8759a5 100644 --- a/test/python/circuit/test_commutation_checker.py +++ b/test/python/circuit/test_commutation_checker.py @@ -13,10 +13,12 @@ """Test commutation checker class .""" import unittest -from test import QiskitTestCase +import pickle +from test import QiskitTestCase # pylint: disable=wrong-import-order import numpy as np from ddt import idata, ddt, data, unpack +from qiskit.quantum_info import Operator from qiskit.circuit import ( AnnotatedOperation, @@ -32,11 +34,13 @@ from qiskit.circuit.library import ( Barrier, CCXGate, + CHGate, CPhaseGate, CRXGate, CRYGate, CRZGate, CXGate, + DCXGate, CUGate, LinearFunction, MCXGate, @@ -53,17 +57,39 @@ RZXGate, RZZGate, SGate, + SdgGate, XGate, YGate, ZGate, HGate, + IGate, UnitaryGate, + CSGate, + CSdgGate, + CZGate, + CYGate, + CSXGate, + SwapGate, + CCZGate, + CSwapGate, + SXGate, + SXdgGate, + TGate, + TdgGate, + iSwapGate, + GlobalPhaseGate, + ECRGate, UGate, PauliEvolutionGate, + XXPlusYYGate, + XXMinusYYGate, PauliProductMeasurement, ) from qiskit.dagcircuit import DAGOpNode from qiskit.quantum_info import SparseObservable, SparsePauliOp, Pauli +from qiskit._accelerate.circuit import StandardGate +from qiskit._accelerate.sparse_observable import _generator_observable +from qiskit._accelerate import standard_generators ROTATION_GATES = [ RXGate, @@ -302,6 +328,16 @@ def test_barrier(self): # Though in this case, it probably makes sense to say that barrier and gate can be swapped. self.assertTrue(scc.commute(Barrier(4), [0, 1, 2, 3], [], CXGate(), [5, 6], [])) + @data( + (CHGate(), [0, 1], HGate(), [1], True), + (UGate(np.pi, 0, np.pi), [0], XGate(), [0], True), + (RGate(np.pi, 0), [0], XGate(), [0], True), + ) + @unpack + def test_unsupported_gates_commute_fallback(self, gate1, q1, gate2, q2, expected): + """Verify that gates without Rust generators still commute correctly via matrix fallback.""" + self.assertEqual(expected, scc.commute(gate1, q1, [], gate2, q2, [])) + def test_reset(self): """Check commutativity involving resets.""" # A gate should not commute with reset when the qubits intersect. @@ -377,7 +413,6 @@ def test_wide_gates_over_disjoint_qubits(self): def test_serialization(self): """Test that the commutation checker is correctly serialized""" - import pickle cx_like = CUGate(np.pi, 0, np.pi, 0) @@ -397,8 +432,6 @@ def test_serialization(self): def test_cutoff_angles(self, gate_cls): """Check rotations with a small enough angle are cut off.""" max_power = 30 - from qiskit.circuit.library import DCXGate - generic_gate = DCXGate() # gate that does not commute with any rotation gate # the cutoff angle depends on the average gate fidelity; i.e. it is the angle @@ -598,5 +631,305 @@ def build_pauli_gate(pauli_string: str, gate_type: str) -> Gate: raise ValueError(f"Invalid gate type: {gate_type}") +@ddt +class TestGeneratorObservableCommutation(QiskitTestCase): + """Test that `_generator_observable` produces correct Pauli generators + for use in commutation checking. + + The function returns H such that gate ≈ exp(-i H) (up to global phase). + Two gates A and B commute iff [H_A, H_B] = 0. We verify here that the + generators correctly predict commutation vs. non-commutation for a + representative set of gate pairs. + """ + + def _commute(self, gate_a, gate_b, qargs_a, qargs_b): + """Compute commutation via the full CommutationChecker.""" + return scc.commute(gate_a, qargs_a, [], gate_b, qargs_b, []) + + def test_cs_commutation(self): + """CS and CZ share the same Z0 Z1 generator subspace, so they commute.""" + self.assertTrue( + self._commute(CSGate(), CZGate(), [0, 1], [0, 1]), + "CS and CZ should commute (both are ZZ-type generators)", + ) + + def test_cx_cy_noncommute(self): + """CX and CY do not commute on overlapping qubits (different generators).""" + self.assertFalse( + self._commute(CXGate(), CYGate(), [0, 1], [0, 1]), + "CX and CY should not commute on the same qubits", + ) + + def test_csx_commutes_with_cx(self): + """CSX and CX share the ZX generator subspace (same up to scale), so they commute.""" + self.assertTrue( + self._commute(CSXGate(), CXGate(), [0, 1], [0, 1]), + "CSX and CX should commute (both ZX-type generators)", + ) + + @data( + StandardGate.CX, + StandardGate.CY, + StandardGate.CZ, + StandardGate.CS, + StandardGate.CSdg, + StandardGate.CSX, + StandardGate.Swap, + StandardGate.CCX, + StandardGate.CCZ, + StandardGate.CSwap, + StandardGate.X, + StandardGate.Y, + StandardGate.Z, + StandardGate.S, + StandardGate.Sdg, + StandardGate.T, + StandardGate.Tdg, + StandardGate.H, + StandardGate.SX, + StandardGate.SXdg, + StandardGate.ISwap, + StandardGate.ECR, + StandardGate.CCZ, + StandardGate.I, + StandardGate.CPhase, + StandardGate.CRX, + StandardGate.CRY, + StandardGate.CRZ, + StandardGate.GlobalPhase, + StandardGate.RXX, + StandardGate.RYY, + StandardGate.RZZ, + StandardGate.RZX, + StandardGate.XXPlusYY, + StandardGate.XXMinusYY, + ) + def test_clifford_gates_have_generators(self, std_gate): + """ + `_generator_observable` should return a SparseObservable for all + the standard gates that PR #15488 adds commutation support for. + """ + # If a gate requires parameters (like XXPlusYY), provide defaults. + # Otherwise, fall back to empty params. + params = [0.1, 0.0] if std_gate.name in ("xx_plus_yy", "xx_minus_yy") else [] + obs = _generator_observable(std_gate, params) + self.assertIsNotNone( + obs, + f"{std_gate.name} should have a generator", + ) + + @data( + (StandardGate.RX, 0.7), + (StandardGate.RY, 0.7), + (StandardGate.RZ, 0.7), + (StandardGate.Phase, 0.7), + (StandardGate.RXX, 0.7), + (StandardGate.RYY, 0.7), + (StandardGate.RZZ, 0.7), + (StandardGate.RZX, 0.7), + (StandardGate.CPhase, 0.7), + (StandardGate.CRX, 0.7), + (StandardGate.CRY, 0.7), + (StandardGate.CRZ, 0.7), + ) + @unpack + def test_rotation_gates_have_generators(self, std_gate, angle): + """Parametric rotation gates should return generators proportional to the angle.""" + obs = _generator_observable(std_gate, [angle]) + self.assertIsNotNone(obs) + + # Controlled rotations typically use angle/4 in their generator expansion + # (e.g. CRX(t) ~ exp(-i * t/4 * (X - ZX))) + expected_scale = 4.0 if std_gate.name in ("cp", "crx", "cry", "crz") else 2.0 + + self.assertAlmostEqual( + abs(obs.coeffs[0]), + abs(angle / expected_scale), + places=10, + msg=f"{std_gate.name} coefficient should be angle/{expected_scale}", + ) + + @data( + (StandardGate.X, XGate()), + (StandardGate.Y, YGate()), + (StandardGate.Z, ZGate()), + (StandardGate.I, IGate()), + (StandardGate.H, HGate()), + (StandardGate.S, SGate()), + (StandardGate.Sdg, SdgGate()), + (StandardGate.T, TGate()), + (StandardGate.Tdg, TdgGate()), + (StandardGate.SX, SXGate()), + (StandardGate.SXdg, SXdgGate()), + (StandardGate.CX, CXGate()), + (StandardGate.CY, CYGate()), + (StandardGate.CZ, CZGate()), + (StandardGate.CS, CSGate()), + (StandardGate.CSdg, CSdgGate()), + (StandardGate.CSX, CSXGate()), + (StandardGate.Swap, SwapGate()), + (StandardGate.ISwap, iSwapGate()), + (StandardGate.ECR, ECRGate()), + (StandardGate.CCX, CCXGate()), + (StandardGate.CCZ, CCZGate()), + (StandardGate.CSwap, CSwapGate()), + (StandardGate.RX, RXGate(0.5)), + (StandardGate.RY, RYGate(0.5)), + (StandardGate.RZ, RZGate(0.5)), + (StandardGate.Phase, PhaseGate(0.5)), + (StandardGate.RXX, RXXGate(0.5)), + (StandardGate.RYY, RYYGate(0.5)), + (StandardGate.RZZ, RZZGate(0.5)), + (StandardGate.RZX, RZXGate(0.5)), + (StandardGate.CPhase, CPhaseGate(0.5)), + (StandardGate.CRX, CRXGate(0.5)), + (StandardGate.CRY, CRYGate(0.5)), + (StandardGate.CRZ, CRZGate(0.5)), + (StandardGate.XXPlusYY, XXPlusYYGate(0.5, 0.0)), + (StandardGate.XXMinusYY, XXMinusYYGate(0.5, 0.0)), + (StandardGate.GlobalPhase, GlobalPhaseGate(0.5)), + ) + @unpack + def test_all_gates_operator_equivalence_ddt(self, std_gate, gate_obj): + """Verify that gate ≈ exp(-i * H) for all supported standard gates.""" + obs = _generator_observable(std_gate, gate_obj.params) + self.assertIsNotNone(obs, f"{std_gate.name} should have a generator") + + # Convert generator SparseObservable to Operator via PauliEvolutionGate + # Convention: gate = exp(-i * H). So time=1.0. + evo = PauliEvolutionGate(obs, time=1.0) + gen_op = Operator(evo) + target_op = Operator(gate_obj) + + # Operator.equiv() checks for equivalence up to global phase + self.assertTrue( + gen_op.equiv(target_op), + f"Generator for {std_gate.name} is not equivalent to the gate unitary. " + f"Dist: {np.linalg.norm(gen_op.data - target_op.data)}", + ) + + def test_large_pauli_evolution_commutation(self): + """Test that a large PauliEvolutionGate correctly commutes with standard gates. + This verifies that matrix-based fallback is avoided for large commutative Pauli strings.""" + pauli_op = SparsePauliOp.from_list([("X" * 50, 1.0)]) + evo_gate = PauliEvolutionGate(pauli_op) + + # Create a circuit with the evolution gate and a standard gate + qr = QuantumRegister(51) + qc = QuantumCircuit(qr) + qc.append(evo_gate, range(50)) + qc.x(50) + + evo_op = qc.data[0] + x_op = qc.data[1] + + cc = scc + + # They act on disjoint qubits, so they should commute + commutes = cc.commute( + evo_op.operation, evo_op.qubits, evo_op.clbits, x_op.operation, x_op.qubits, x_op.clbits + ) + self.assertTrue(commutes) + + # Test commuting with an X gate that overlaps (on qubit 0). + # "X...X" and X commute! + qc.x(0) + x_op2 = qc.data[2] + commutes_overlap = cc.commute( + evo_op.operation, + evo_op.qubits, + evo_op.clbits, + x_op2.operation, + x_op2.qubits, + x_op2.clbits, + ) + self.assertTrue(commutes_overlap) + + # Test non-commuting with Y gate on qubit 0 + qc.y(0) + y_op = qc.data[3] + commutes_y = cc.commute( + evo_op.operation, evo_op.qubits, evo_op.clbits, y_op.operation, y_op.qubits, y_op.clbits + ) + self.assertFalse(commutes_y) + + @data( + StandardGate.CH, + StandardGate.R, + StandardGate.U, + StandardGate.U1, + StandardGate.U2, + StandardGate.U3, + StandardGate.CU, + StandardGate.CU3, + ) + def test_unsupported_gates_return_none(self, std_gate): + """Standard gates not yet implemented in Rust should return None.""" + obs = _generator_observable(std_gate, [0.1] * std_gate.num_params) + self.assertIsNone(obs, f"{std_gate.name} should not have a generator implementation yet") + + @data( + (StandardGate.XXPlusYY, [0.5, 0.1]), + (StandardGate.XXMinusYY, [0.5, 0.1]), + ) + @unpack + def test_unsupported_params_return_none(self, std_gate, params): + """Standard gates with unsupported parameter values should return None.""" + obs = _generator_observable(std_gate, params) + self.assertIsNone( + obs, f"{std_gate.name} with beta={params[1]} should not have a generator implementation" + ) + + def test_pauli_product_rotation_commutation(self): + """Test commutation between standard gates and PauliProductRotationGate.""" + from qiskit.circuit.library import PauliProductRotationGate + + # CX(0,1) commutes with exp(-i*theta*XI) and exp(-i*theta*IZ) + # but not exp(-i*theta*IX) + ppr_xi = PauliProductRotationGate(Pauli("XI"), 0.5) + ppr_iz = PauliProductRotationGate(Pauli("IZ"), 0.5) + ppr_ix = PauliProductRotationGate(Pauli("IX"), 0.5) + + self.assertTrue(scc.commute(CXGate(), [0, 1], [], ppr_xi, [0, 1], [])) + self.assertTrue(scc.commute(CXGate(), [0, 1], [], ppr_iz, [0, 1], [])) + self.assertFalse(scc.commute(CXGate(), [0, 1], [], ppr_ix, [0, 1], [])) + + def test_large_qubit_scaling(self): + """Test commutation with a very large Number of qubits.""" + num_qubits = 100 + # PauliEvolutionGate on qubits 0 and 99 + op = SparsePauliOp.from_list([("Z" + "I" * (num_qubits - 2) + "Z", 1.0)]) + evo = PauliEvolutionGate(op, time=0.5) + + # X gate on qubit 50 should commute + self.assertTrue(scc.commute(XGate(), [50], [], evo, list(range(num_qubits)), [])) + + # X gate on qubit 0 should not commute + self.assertFalse(scc.commute(XGate(), [0], [], evo, list(range(num_qubits)), [])) + + def test_exposed_generator_observable(self): + """Test the exposed Python binding for generator_observable.""" + # CCXGate generator should be exposed + obs = standard_generators.generator_observable(StandardGate.CCX) + self.assertIsInstance(obs, SparseObservable) + self.assertEqual(obs.num_terms, 7) + + @data( + (RXXGate(0.5), [0, 1], XGate(), [0], True), + (RXXGate(0.5), [0, 1], ZGate(), [0], False), + (RYYGate(0.5), [0, 1], YGate(), [1], True), + (RZZGate(0.5), [0, 1], ZGate(), [1], True), + (RZXGate(0.5), [0, 1], ZGate(), [0], True), + (RZXGate(0.5), [0, 1], XGate(), [1], True), + (XXPlusYYGate(0.5), [0, 1], ZGate(), [0], False), + (GlobalPhaseGate(0.5), [], XGate(), [0], True), + ) + @unpack + def test_more_standard_gates_commutation(self, gate1, q1, gate2, q2, expected): + """Test more standard gates requested by reviewers.""" + res = scc.commute(gate1, q1, [], gate2, q2, []) + self.assertEqual(res, expected, f"Commutation of {gate1.name} and {gate2.name} failed") + + if __name__ == "__main__": unittest.main() diff --git a/test/python/transpiler/test_light_cone.py b/test/python/transpiler/test_light_cone.py index fb5ccaab0b0e..cecdc34da6f2 100644 --- a/test/python/transpiler/test_light_cone.py +++ b/test/python/transpiler/test_light_cone.py @@ -13,6 +13,7 @@ """Test the LightCone pass""" import unittest +import numpy as np from test import QiskitTestCase import ddt @@ -21,7 +22,7 @@ Parameter, QuantumCircuit, ) -from qiskit.circuit.library import real_amplitudes +from qiskit.circuit.library import RXXGate, real_amplitudes from qiskit.circuit.library.n_local.efficient_su2 import efficient_su2 from qiskit.converters import circuit_to_dag from qiskit.quantum_info import SparsePauliOp, SparseObservable @@ -359,6 +360,41 @@ def test_raise_error_when_circuit_measurements_and_observable_present(self): ): light_cone.run(dag) + def test_rxx_commuting(self): + """Test for a commuting RXX gate in the LightCone pass.""" + # Simple test: X commutes with RXX(0, 1), so it should be dropped + qc_simple = QuantumCircuit(2) + qc_simple.append(RXXGate(0.5), [0, 1]) + + light_cone = LightCone(bit_terms="X", indices=[0]) + pm = PassManager([light_cone]) + new_circuit = pm.run(qc_simple) + self.assertEqual(len(new_circuit), 0) + + # Non-commuting test: Z does not commute with RXX(0, 1) + lc_z = LightCone(bit_terms="Z", indices=[0]) + pm_z = PassManager([lc_z]) + new_z = pm_z.run(qc_simple) + self.assertEqual(len(new_z), 1) + + def test_original_bug_15021(self): + """Test for the scattered-qubit bug reported in issue #15021.""" + qc = QuantumCircuit(18) + qc.rx(np.pi / 3, range(qc.num_qubits)) + bit_terms = "ZZZZZZZZZ" + indices = [0, 1, 2, 3, 4, 9, 10, 12, 13] + pm = PassManager( + [ + LightCone( + bit_terms=bit_terms, + indices=indices, + ) + ] + ) + # This used to fail with Index results / UnsortedIndices + # Verify the reduced circuit contains exactly 9 gates on the expected light cone + self.assertEqual(len(pm.run(qc)), 9) + if __name__ == "__main__": unittest.main()