Refactor PauliProductRotation matrix support to reuse SparsePauliOp logic

Author: Adithyaphani

crates/circuit/src/operations.rs

@@ -3122,23 +3122,6 @@ impl PyInstruction {
         })
     }
 
-    pub fn matrix(&self) -> Option<Array2<Complex64>> {
-        Python::attach(|py| -> Option<Array2<Complex64>> {
-            match self.instruction.getattr(py, intern!(py, "to_matrix")) {
-                Ok(to_matrix) => {
-                    let res: Option<Py<PyAny>> = to_matrix.call0(py).ok()?.extract(py).ok();
-                    match res {
-                        Some(x) => {
-                            let array: PyReadonlyArray2<Complex64> = x.extract(py).ok()?;
-                            Some(array.as_array().to_owned())
-                        }
-                        None => None,
-                    }
-                }
-                Err(_) => None,
-            }
-        })
-    }
 
     pub fn definition(&self) -> Option<CircuitData> {
         Python::attach(|py| -> Option<CircuitData> {
@@ -3170,6 +3153,25 @@ impl PyInstruction {
             Some([[arr[[0, 0]], arr[[0, 1]]], [arr[[1, 0]], arr[[1, 1]]]])
         })
     }
+
+    pub fn matrix(&self) -> Option<Array2<Complex64>> {
+        Python::attach(|py| -> Option<Array2<Complex64>> {
+            match self.instruction.getattr(py, intern!(py, "to_matrix")) {
+                Ok(to_matrix) => {
+                    let res: Option<Py<PyAny>> = to_matrix.call0(py).ok()?.extract(py).ok();
+                    match res {
+                        Some(x) => {
+                            let array: PyReadonlyArray2<Complex64> = x.extract(py).ok()?;
+                            Some(array.as_array().to_owned())
+                        }
+                        None => None,
+                    }
+                }
+                Err(_) => None,
+            }
+        })
+    }
+
 }
 
 #[derive(Clone, Debug)]
@@ -3336,19 +3338,6 @@ pub struct PauliProductRotation {
 }
 
 
-// Helper function to convert x/z bool vectors into a Pauli string in "IXYZ" notation.
-fn pauli_product_to_string(zs: &[bool], xs: &[bool]) -> String {
-    zs.iter()
-        .zip(xs.iter())
-        .map(|(&z, &x)| match (z, x) {
-            (false, false) => 'I',
-            (false, true) => 'X',
-            (true, false) => 'Z',
-            (true, true) => 'Y',
-        })
-        .collect()
-}
-
 impl PauliProductRotation {
     pub fn create_py_op(&self, py: Python, label: Option<&str>) -> PyResult<Py<PyAny>> {
         let z = self.z.to_pyarray(py);
@@ -3368,18 +3357,25 @@ impl PauliProductRotation {
             )?;
         Ok(gate.unbind())
     }
-
     pub fn matrix(&self) -> Option<Array2<Complex64>> {
         let angle = match self.angle {
             Param::Float(f) => f,
             _ => return None,
         };
-        Some(gate_matrix::pauli_product_rotation_matrix(
-            angle,
+        let pauli_mat = gate_matrix::pauli_zx_to_dense_matrix(
             &self.z,
             &self.x,
-        ))
+        );
+        let cos_val = (angle / 2.0).cos();
+        let sin_val = (angle / 2.0).sin();
+        let dim = pauli_mat.shape()[0];
+        let identity = Array2::<Complex64>::eye(dim);
+        Some(
+            identity.mapv(|v| Complex64::new(v.re * cos_val, 0.))
+                + pauli_mat.mapv(|v| Complex64::new(0., -sin_val) * v)
+        )
     }
+
 }
 
 impl Operation for PauliProductRotation {

crates/quantum_info/src/sparse_pauli_op.rs

@@ -992,6 +992,62 @@ fn to_matrix_dense_inner(paulis: &MatrixCompressedPaulis, parallel: bool) -> Vec
     out
 }
 
+
+pub fn pauli_zx_to_matrix(
+    z: &[bool],
+    x: &[bool],
+) -> ndarray::Array2<num_complex::Complex64> {
+    use ndarray::Array2;
+    use num_complex::Complex64;
+
+    assert_eq!(z.len(), x.len());
+
+    let n = z.len();
+    let dim = 1usize << n;
+    let mut out = Array2::<Complex64>::zeros((dim, dim));
+
+    for col in 0..dim {
+        let mut row = col;
+        let mut phase = Complex64::new(1.0, 0.0);
+
+        for i in 0..n {
+            let bit_index = n - 1 - i;
+            let bit = (col >> bit_index) & 1;
+
+            match (z[i], x[i]) {
+                (false, false) => {
+                    // I
+                }
+                (false, true) => {
+                    // X
+                    row ^= 1 << bit_index;
+                }
+                (true, false) => {
+                    // Z
+                    if bit == 1 {
+                        phase = -phase;
+                    }
+                }
+                (true, true) => {
+                    // Y = iXZ
+                    row ^= 1 << bit_index;
+                    phase *= if bit == 0 {
+                        Complex64::new(0.0, 1.0)
+                    } else {
+                        Complex64::new(0.0, -1.0)
+                    };
+                }
+            }
+        }
+
+        out[(row, col)] = phase;
+    }
+
+    out
+}
+
+
+
 type CSRData<T> = (Vec<Complex64>, Vec<T>, Vec<T>);
 type ToCSRData<T> = fn(&MatrixCompressedPaulis) -> CSRData<T>;
 

test/python/circuit/test_pauli_product_rotation_matrix.py

@@ -0,0 +1,141 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2026.
+#
+# This code is licensed under the Apache License, Version 2.0.
+"""Tests for PauliProductRotation matrix methods in Rust (Issue #15869)."""
+
+import numpy as np
+import pytest
+from qiskit.circuit import QuantumCircuit
+from qiskit.circuit.library.generalized_gates.pauli_product_rotation import (
+    PauliProductRotationGate,
+)
+from qiskit.quantum_info import Pauli
+
+class TestPauliProductRotationMatrix:
+    """Test that PauliProductRotation returns correct unitary matrices."""
+
+    def test_single_qubit_x_rotation(self):
+        """PPR with 'X' should match RX gate matrix."""
+        theta = np.pi / 3
+        gate = PauliProductRotationGate(Pauli("X"), theta)
+        mat = gate.to_matrix()
+
+        # RX(theta) = cos(θ/2)*I - i*sin(θ/2)*X
+        expected = np.array([
+            [np.cos(theta / 2), -1j * np.sin(theta / 2)],
+            [-1j * np.sin(theta / 2), np.cos(theta / 2)],
+        ])
+        np.testing.assert_allclose(mat, expected, atol=1e-10,
+            err_msg="PPR('X', θ) must match RX(θ)")
+
+    def test_single_qubit_z_rotation(self):
+        """PPR with 'Z' should match RZ gate (up to global phase)."""
+        theta = np.pi / 4
+        gate = PauliProductRotationGate(Pauli("Z"),theta)
+        mat = gate.to_matrix()
+
+        # exp(-i*θ/2*Z) = diag(exp(-iθ/2), exp(+iθ/2))
+        expected = np.diag([
+            np.exp(-1j * theta / 2),
+            np.exp(+1j * theta / 2),
+        ])
+        np.testing.assert_allclose(mat, expected, atol=1e-10,
+            err_msg="PPR('Z', θ) must match RZ(θ)")
+
+    def test_two_qubit_zz(self):
+        """PPR 'ZZ' must be a 4×4 unitary."""
+        theta = np.pi / 5
+        gate = PauliProductRotationGate(Pauli("ZZ"),theta)
+        mat = gate.to_matrix()
+        assert mat.shape == (4, 4), "ZZ gate must be 4×4"
+
+        # Must be unitary: U† U = I
+        ident = mat.conj().T @ mat
+        np.testing.assert_allclose(ident, np.eye(4), atol=1e-10,
+            err_msg="PPR('ZZ', θ) matrix must be unitary")
+
+    def test_three_qubit_xyz(self):
+        """PPR 'XYZ' must be an 8×8 unitary."""
+        theta = 0.7
+        gate = PauliProductRotationGate(Pauli("XYZ"),theta)
+        mat = gate.to_matrix()
+        assert mat.shape == (8, 8), "XYZ gate must be 8×8"
+
+        ident = mat.conj().T @ mat
+        np.testing.assert_allclose(ident, np.eye(8), atol=1e-10,
+            err_msg="PPR('XYZ', θ) matrix must be unitary")
+
+    def test_identity_pauli_is_global_phase(self):
+        """PPR with all-identity Pauli 'II' reduces to a global phase gate."""
+        theta = np.pi / 6
+        gate = PauliProductRotationGate(Pauli("II"),theta)
+        mat = gate.to_matrix()
+
+        # exp(-i*θ/2*I⊗I) = exp(-i*θ/2) * I4
+        expected = np.exp(-1j * theta / 2) * np.eye(4)
+        np.testing.assert_allclose(mat, expected, atol=1e-10)
+
+    def test_theta_zero_is_identity(self):
+        """At θ=0, PPR must be the identity (cos(0)=1, sin(0)=0)."""
+        gate = PauliProductRotationGate(Pauli("XZ"),0.0)
+        mat = gate.to_matrix()
+        np.testing.assert_allclose(mat, np.eye(4), atol=1e-10,
+            err_msg="PPR at θ=0 must be identity")
+
+    def test_theta_2pi_is_negative_identity(self):
+        """At θ=2π, PPR = -I (a global phase of -1)."""
+        gate = PauliProductRotationGate(Pauli("Z"),2 * np.pi)
+        mat = gate.to_matrix()
+        np.testing.assert_allclose(mat, -np.eye(2), atol=1e-10)
+
+    def test_packed_instruction_try_matrix(self):
+        """PackedInstruction.try_matrix must work for PauliProductRotation."""
+        from qiskit._accelerate.circuit import CircuitData
+
+        theta = np.pi / 7
+        gate = PauliProductRotationGate(Pauli("XX"),theta)
+        qc = QuantumCircuit(2)
+        qc.append(gate, [0, 1])
+
+        # Access inner packed instruction
+        packed = qc._data[0]
+        mat = packed.operation.to_matrix()  # goes through Rust try_matrix
+        assert mat is not None, "try_matrix must not return None for PPR"
+        assert mat.shape == (4, 4)
+
+        # Must match direct gate matrix
+        np.testing.assert_allclose(mat, gate.to_matrix(), atol=1e-10)
+
+    def test_matrix_consistent_with_simulation(self):
+        """Statevector evolved by PPR must match matrix multiplication."""
+        from qiskit.quantum_info import Statevector
+
+        theta = np.pi / 3
+        gate = PauliProductRotationGate(Pauli("ZZ"),theta)
+        qc = QuantumCircuit(2)
+        qc.h(0)
+        qc.cx(0, 1)
+        qc.append(gate, [0, 1])
+
+        sv = Statevector(qc)
+        mat = gate.to_matrix()
+        assert sv.is_valid(), "Statevector must remain normalized"
+
+    @pytest.mark.parametrize("pauli,n_qubits", [
+        ("X", 1), ("Y", 1), ("Z", 1),
+        ("XX", 2), ("YY", 2), ("ZZ", 2), ("XY", 2), ("YZ", 2),
+        ("XXX", 3), ("ZZZ", 3),
+    ])
+    def test_unitary_parametrized(self, pauli, n_qubits):
+        """All common Pauli strings must produce valid unitaries."""
+        theta = 1.23
+        gate = PauliProductRotationGate(Pauli(pauli), theta)
+        mat = gate.to_matrix()
+        dim = 2**n_qubits
+        assert mat.shape == (dim, dim)
+        np.testing.assert_allclose(
+            mat.conj().T @ mat, np.eye(dim), atol=1e-10,
+            err_msg=f"PPR('{pauli}') matrix not unitary"
+        )