Add Rust matrix support for PauliProductRotation

crates/circuit/src/operations.rs

@@ -9,7 +9,6 @@
 // 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.
-
 use approx::relative_eq;
 use std::f64::consts::PI;
 use std::fmt::Debug;
@@ -3262,24 +3261,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> {
             match self.instruction.getattr(py, intern!(py, "definition")) {
@@ -3311,6 +3292,24 @@ 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 matrix_as_static_2q(&self) -> Option<[[Complex64; 4]; 4]> {
         if self.num_qubits() != 2 {
             return None;
@@ -3550,24 +3549,6 @@ pub struct PauliProductRotation {
     pub angle: Param,
 }
 
-impl Operation for PauliProductRotation {
-    fn name(&self) -> &str {
-        "pauli_product_rotation"
-    }
-    fn num_qubits(&self) -> u32 {
-        self.z.len() as u32
-    }
-    fn num_clbits(&self) -> u32 {
-        0
-    }
-    fn num_params(&self) -> u32 {
-        1
-    }
-    fn directive(&self) -> bool {
-        false
-    }
-}
-
 impl PauliProductRotation {
     pub fn create_py_op(&self, py: Python, label: Option<&str>) -> PyResult<Py<PyAny>> {
         let z = self.z.to_pyarray(py);
@@ -3587,6 +3568,43 @@ impl PauliProductRotation {
             )?;
         Ok(gate.unbind())
     }
+    pub fn matrix(&self) -> Option<Array2<Complex64>> {
+        let angle = match self.angle {
+            Param::Float(f) => f,
+            _ => return None,
+        };
+        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 {
+    fn name(&self) -> &str {
+        "pauli_product_rotation"
+    }
+
+    fn num_qubits(&self) -> u32 {
+        self.z.len() as u32
+    }
+
+    fn num_clbits(&self) -> u32 {
+        0
+    }
+
+    fn num_params(&self) -> u32 {
+        1
+    }
+
+    fn directive(&self) -> bool {
+        false
+    }
 }
 
 impl PartialEq for PauliProductRotation {

crates/circuit/src/packed_instruction.rs

@@ -23,6 +23,7 @@ use crate::operations::{
     PyOperationTypes, PythonOperation, StandardGate, StandardInstruction, UnitaryGate,
 };
 use crate::{Block, Clbit, Qubit};
+
 use hashbrown::HashMap;
 use nalgebra::{Matrix2, Matrix4};
 use ndarray::{Array2, CowArray, Ix2};

crates/quantum_info/src/sparse_pauli_op.rs

@@ -992,6 +992,56 @@ 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>;
 

releasenotes/notes/ppr-rust-matrix-support-15869.rst

@@ -0,0 +1,6 @@
+---
+features:
+  - |
+    Added Rust-side matrix support for ``PauliProductRotationGate`` by reusing
+    the existing Pauli matrix construction logic from the Rust
+    ``SparsePauliOp`` implementation.

test/python/circuit/library/test_pauli_product_rotation.py

@@ -24,7 +24,7 @@
     PauliEvolutionGate,
     PauliProductMeasurement,
 )
-from qiskit.quantum_info import Pauli, Operator, SparsePauliOp
+from qiskit.quantum_info import Pauli, Operator, SparsePauliOp, Statevector
 from qiskit.qpy import dump, load
 from qiskit.compiler import transpile
 from test import QiskitTestCase  # pylint: disable=wrong-import-order
@@ -223,6 +223,119 @@ def test_matrix_conventions(self):
             evo = PauliEvolutionGate(pauli, time=angle / 2)
             self.assertTrue(np.allclose(evo.to_matrix(), ppr.to_matrix()))
 
+    def test_single_qubit_x_rotation_matrix(self):
+        """Test matrix support for a single-qubit X rotation."""
+        theta = np.pi / 3
+        gate = PauliProductRotationGate(Pauli("X"), theta)
+        mat = np.asarray(Operator(gate).data)
+
+        expected = np.array(
+            [
+                [np.cos(theta / 2), -1j * np.sin(theta / 2)],
+                [-1j * np.sin(theta / 2), np.cos(theta / 2)],
+            ],
+            dtype=complex,
+        )
+
+        self.assertEqual(mat.shape, (2, 2))
+        self.assertTrue(np.allclose(mat, expected))
+
+    def test_single_qubit_z_rotation_matrix(self):
+        """Test matrix support for a single-qubit Z rotation."""
+        theta = np.pi / 4
+        gate = PauliProductRotationGate(Pauli("Z"), theta)
+        mat = np.asarray(Operator(gate).data)
+
+        expected = np.array(
+            [
+                [np.exp(-1j * theta / 2), 0],
+                [0, np.exp(1j * theta / 2)],
+            ],
+            dtype=complex,
+        )
+
+        self.assertEqual(mat.shape, (2, 2))
+        self.assertTrue(np.allclose(mat, expected))
+
+    def test_two_qubit_zz_matrix(self):
+        """Test matrix support for a two-qubit ZZ rotation."""
+        theta = np.pi / 5
+        gate = PauliProductRotationGate(Pauli("ZZ"), theta)
+        mat = np.asarray(Operator(gate).data)
+
+        self.assertEqual(mat.shape, (4, 4))
+        self.assertTrue(np.allclose(mat.conj().T @ mat, np.eye(4)))
+
+    def test_three_qubit_xyz_matrix(self):
+        """Test matrix support for a three-qubit XYZ rotation."""
+        theta = 0.7
+        gate = PauliProductRotationGate(Pauli("XYZ"), theta)
+        mat = np.asarray(Operator(gate).data)
+
+        self.assertEqual(mat.shape, (8, 8))
+        self.assertTrue(np.allclose(mat.conj().T @ mat, np.eye(8)))
+
+    def test_identity_pauli_is_global_phase(self):
+        """Test that an all-identity Pauli only contributes a global phase."""
+        theta = np.pi / 6
+        gate = PauliProductRotationGate(Pauli("II"), theta)
+        mat = np.asarray(Operator(gate).data)
+        expected = np.exp(-1j * theta / 2) * np.eye(4)
+
+        self.assertTrue(np.allclose(mat, expected))
+
+    def test_theta_zero_is_identity(self):
+        """Test that zero angle gives the identity."""
+        gate = PauliProductRotationGate(Pauli("XZ"), 0.0)
+        mat = np.asarray(Operator(gate).data)
+
+        self.assertTrue(np.allclose(mat, np.eye(4)))
+
+    def test_theta_2pi_is_negative_identity(self):
+        """Test that a 2π rotation gives -I up to global phase."""
+        gate = PauliProductRotationGate(Pauli("Z"), 2 * np.pi)
+        mat = np.asarray(Operator(gate).data)
+
+        self.assertTrue(np.allclose(mat, -np.eye(2)))
+
+    def test_matrix_consistent_with_simulation(self):
+        """Test matrix action is consistent with statevector simulation."""
+        theta = np.pi / 3
+        gate = PauliProductRotationGate(Pauli("ZZ"), theta)
+
+        qc = QuantumCircuit(2)
+        qc.append(gate, [0, 1])
+
+        psi0 = Statevector.from_label("00")
+        evolved = psi0.evolve(qc)
+        mat = np.asarray(Operator(gate).data)
+        expected = Statevector(mat @ psi0.data)
+
+        self.assertTrue(np.allclose(evolved.data, expected.data))
+
+    @data(
+        ("X", 1),
+        ("Y", 1),
+        ("Z", 1),
+        ("XX", 2),
+        ("YY", 2),
+        ("ZZ", 2),
+        ("XY", 2),
+        ("YZ", 2),
+        ("XXX", 3),
+        ("ZZZ", 3),
+    )
+    def test_matrix_unitary_for_common_paulis(self, case):
+        """Test matrix support for common Pauli strings."""
+        pauli, n_qubits = case
+        theta = 1.23
+
+        gate = PauliProductRotationGate(Pauli(pauli), theta)
+        mat = np.asarray(Operator(gate).data)
+
+        self.assertEqual(mat.shape, (2**n_qubits, 2**n_qubits))
+        self.assertTrue(np.allclose(mat.conj().T @ mat, np.eye(2**n_qubits)))
+
     @data(0, 1, 2)
     def test_control(self, num_ctrl_qubits):
         """Check calling the control method."""