Add pauli_twirl_2q_gates function

crates/accelerate/src/twirling.rs

@@ -12,7 +12,11 @@
 
 use std::f64::consts::PI;
 
-use hashbrown::HashSet;
+use hashbrown::{HashMap, HashSet};
+use ndarray::linalg::kron;
+use ndarray::prelude::*;
+use ndarray::ArrayView2;
+use num_complex::Complex64;
 use pyo3::intern;
 use pyo3::prelude::*;
 use pyo3::types::PyList;
@@ -23,9 +27,10 @@ use rand_pcg::Pcg64Mcg;
 use smallvec::SmallVec;
 
 use qiskit_circuit::circuit_data::CircuitData;
-use qiskit_circuit::circuit_instruction::ExtraInstructionAttributes;
+use qiskit_circuit::circuit_instruction::{ExtraInstructionAttributes, OperationFromPython};
 use qiskit_circuit::converters::dag_to_circuit;
 use qiskit_circuit::dag_circuit::DAGCircuit;
+use qiskit_circuit::gate_matrix::ONE_QUBIT_IDENTITY;
 use qiskit_circuit::imports::QUANTUM_CIRCUIT;
 use qiskit_circuit::operations::StandardGate::{IGate, XGate, YGate, ZGate};
 use qiskit_circuit::operations::{Operation, OperationRef, Param, PyInstruction, StandardGate};
@@ -123,12 +128,59 @@ const CZ_MASK: u8 = 4;
 const ECR_MASK: u8 = 2;
 const ISWAP_MASK: u8 = 1;
 
+#[inline(always)]
+fn diff_frob_norm(array: ArrayView2<Complex64>, gate_matrix: ArrayView2<Complex64>) -> f64 {
+    let mut res: f64 = 0.;
+    for i in 0..4 {
+        for j in 0..4 {
+            let gate = gate_matrix[[i, j]];
+            let twirled = array[[i, j]];
+            let diff = twirled - gate;
+            res += (diff.conj() * diff).re;
+        }
+    }
+    res.sqrt()
+}
+
+fn generate_twirling_set(gate_matrix: ArrayView2<Complex64>) -> Vec<([StandardGate; 4], f64)> {
+    let mut out_vec = Vec::with_capacity(16);
+    let i_matrix = aview2(&ONE_QUBIT_IDENTITY);
+    let x_matrix = aview2(&qiskit_circuit::gate_matrix::X_GATE);
+    let y_matrix = aview2(&qiskit_circuit::gate_matrix::Y_GATE);
+    let z_matrix = aview2(&qiskit_circuit::gate_matrix::Z_GATE);
+    let iter_set = [
+        (IGate, i_matrix),
+        (XGate, x_matrix),
+        (YGate, y_matrix),
+        (ZGate, z_matrix),
+    ];
+    for (i, i_mat) in &iter_set {
+        for (j, j_mat) in &iter_set {
+            let before_matrix = kron(j_mat, i_mat);
+            let half_twirled_matrix = gate_matrix.dot(&before_matrix);
+            for (k, k_mat) in &iter_set {
+                for (l, l_mat) in &iter_set {
+                    let after_matrix = kron(l_mat, k_mat);
+                    let twirled_matrix = after_matrix.dot(&half_twirled_matrix);
+                    let norm: f64 = diff_frob_norm(twirled_matrix.view(), gate_matrix);
+                    if norm.abs() < 1e-15 {
+                        out_vec.push(([*i, *j, *k, *l], 0.));
+                    } else if (norm - 4.).abs() < 1e-15 {
+                        out_vec.push(([*i, *j, *k, *l], PI));
+                    }
+                }
+            }
+        }
+    }
+    out_vec
+}
+
 fn twirl_gate(
     py: Python,
     circ: &CircuitData,
     rng: &mut Pcg64Mcg,
     out_circ: &mut CircuitData,
-    twirl_set: &[([StandardGate; 4], f64); 16],
+    twirl_set: &[([StandardGate; 4], f64)],
     inst: &PackedInstruction,
 ) -> PyResult<()> {
     let qubits = circ.get_qargs(inst.qubits);
@@ -198,6 +250,7 @@ fn generate_twirled_circuit(
     circ: &CircuitData,
     rng: &mut Pcg64Mcg,
     twirling_mask: u8,
+    custom_gate_map: Option<&HashMap<String, Vec<([StandardGate; 4], f64)>>>,
     run_pass: bool,
     optimizer_target: Option<&Target>,
     optimizer_global_decomposer: Option<&HashSet<String>>,
@@ -206,6 +259,12 @@ fn generate_twirled_circuit(
     let mut out_circ = CircuitData::clone_empty_like(circ, None);
 
     for inst in circ.data() {
+        if let Some(custom_gate_map) = custom_gate_map {
+            if let Some(twirling_set) = custom_gate_map.get(inst.op.name()) {
+                twirl_gate(py, circ, rng, &mut out_circ, twirling_set.as_slice(), inst)?;
+                continue;
+            }
+        }
         match inst.op.view() {
             OperationRef::Standard(gate) => match gate {
                 StandardGate::CXGate => {
@@ -249,6 +308,7 @@ fn generate_twirled_circuit(
                                 block,
                                 rng,
                                 twirling_mask,
+                                custom_gate_map,
                                 run_pass,
                                 optimizer_target,
                                 optimizer_global_decomposer,
@@ -324,11 +384,12 @@ fn generate_twirled_circuit(
 
 #[allow(clippy::too_many_arguments)]
 #[pyfunction]
-#[pyo3(signature=(circ, twirled_gate, seed=None, num_twirls=1, run_pass=false, optimizer_target=None, optimizer_global_decomposer=None, optimizer_basis_gates=None))]
+#[pyo3(signature=(circ, twirled_gate=None, custom_twirled_gates=None, seed=None, num_twirls=1, run_pass=false, optimizer_target=None, optimizer_global_decomposer=None, optimizer_basis_gates=None))]
 pub(crate) fn twirl_circuit(
     py: Python,
     circ: &CircuitData,
     twirled_gate: Option<Vec<StandardGate>>,
+    custom_twirled_gates: Option<Vec<OperationFromPython>>,
     seed: Option<u64>,
     num_twirls: usize,
     run_pass: bool,
@@ -359,16 +420,43 @@ pub(crate) fn twirl_circuit(
             }
             out_mask
         }
-        None => 15,
+        None => {
+            if custom_twirled_gates.is_none() {
+                15
+            } else {
+                0
+            }
+        }
     };
-
+    let custom_gate_twirling_sets: Option<HashMap<String, Vec<([StandardGate; 4], f64)>>> =
+        custom_twirled_gates.map(|gates| {
+            gates
+                .into_iter()
+                .filter_map(|gate| {
+                    let matrix = gate.operation.matrix(&gate.params);
+                    if let Some(matrix) = matrix {
+                        let twirl_set = generate_twirling_set(matrix.view());
+                        if twirl_set.is_empty() {
+                            None
+                        } else {
+                            Some(Ok((gate.operation.name().to_string(), twirl_set)))
+                        }
+                    } else {
+                        return Some(Err(QiskitError::new_err(
+                            format!("Provided gate to twirl, {}, does not have a matrix defined and can't be twirled", gate.operation.name())
+                        )));
+                    }
+                })
+                .collect()
+        }).transpose()?;
     (0..num_twirls)
         .map(|_| {
             generate_twirled_circuit(
                 py,
                 circ,
                 &mut rng,
                 twirling_mask,
+                custom_gate_twirling_sets.as_ref(),
                 run_pass,
                 optimizer_target,
                 optimizer_global_decomposer.as_ref(),

qiskit/circuit/twirling.py

@@ -35,7 +35,7 @@
 
 def twirl_circuit(
     circuit: QuantumCircuit,
-    twirling_gate: None | str | type[Gate] | list[str] | list[type[Gate]] = None,
+    twirling_gate: None | str | Gate | list[str] | list[Gate] = None,
     seed: int | None = None,
     num_twirls: int | None = None,
     optimize_pass: Optimize1qGatesDecomposition | None = None,
@@ -45,10 +45,14 @@ def twirl_circuit(
 
     Args:
         circuit: The circuit to twirl
-        twirling_gate: The gate to twirl, defaults to `None` which means twirl all supported gates.
-            If supplied it can either be a single gate or a list of gates either as a gate class
-            or it's string name. Currently only :class:`.CXGate` (`"cx"`), :class:`.CZGate` (`"cz"`),
-            :class:`.ECRGate` (`"ecr"`), and :class:`.iSwapGate` (`"iswap"`) are supported.
+        twirling_gate: The gate to twirl, defaults to `None` which means twirl all default gates:
+            :class:`.CXGate`, :class:`.CZGate`, :class:`.ECRGate`, and :class:`.iSwapGate`.
+            If supplied it can either be a single gate or a list of gates either as either a gate
+            object or it's string name. Currently only the names `"cx"`, `"cz"`, `"ecr"`,  and
+            `"iswap"` are supported. If a gate object is provided outside the default gates it must
+            have a matrix defined from it's :class:`~.Gate.to_matrix` method for the gate to potentially
+            be twirled. If a valid twirling configuration can't be computed that particular gate will
+            be silently ignored and not twirled.
         seed: An integer seed for the random number generator used internally.
             If specified this must be between 0 and 18,446,744,073,709,551,615.
         num_twirls: The number of twirling circuits to build. This defaults to None and will return
@@ -62,12 +66,14 @@ def twirl_circuit(
         A copy of the given circuit with Pauli twirling applied to each
         instance of the specified twirling gate.
     """
+    custom_gates = None
     if isinstance(twirling_gate, str):
         gate = NAME_TO_CLASS.get(twirling_gate, None)
         if gate is None:
             raise QiskitError(f"The specified gate name {twirling_gate} is not supported")
         twirling_std_gate = [gate]
     elif isinstance(twirling_gate, list):
+        custom_gates = []
         twirling_std_gate = []
         for gate in twirling_gate:
             if isinstance(gate, str):
@@ -77,14 +83,29 @@ def twirl_circuit(
                 twirling_std_gate.append(gate)
             else:
                 twirling_gate = getattr(gate, "_standard_gate", None)
+
                 if twirling_gate is None:
-                    raise QiskitError("This function can only be used with standard gates")
-                twirling_std_gate.append(twirling_gate)
+                    custom_gates.append(gate)
+                else:
+                    if twirling_gate in NAME_TO_CLASS.values():
+                        twirling_std_gate.append(twirling_gate)
+                    else:
+                        custom_gates.append(gate)
+        if not custom_gates:
+            custom_gates = None
+        if not twirling_std_gate:
+            twirling_std_gate = None
     elif twirling_gate is not None:
-        twirling_std_gate = getattr(twirling_gate, "_standard_gate", None)
-        if twirling_std_gate is None:
-            raise QiskitError("This function can only be used with standard gates")
-        twirling_std_gate = [twirling_std_gate]
+        std_gate = getattr(twirling_gate, "_standard_gate", None)
+        if std_gate is None:
+            twirling_std_gate = None
+            custom_gates = [twirling_gate]
+        else:
+            if std_gate in NAME_TO_CLASS.values():
+                twirling_std_gate = [std_gate]
+            else:
+                twirling_std_gate = None
+                custom_gates = [twirling_gate]
     else:
         twirling_std_gate = twirling_gate
     out_twirls = num_twirls
@@ -101,6 +122,7 @@ def twirl_circuit(
     new_data = twirl_rs(
         circuit._data,
         twirling_std_gate,
+        custom_gates,
         seed,
         out_twirls,
         run_pass,

test/python/circuit/test_twirling.py

@@ -67,7 +67,7 @@ def test_twirl_circuit_list(self):
         qc.cz(0, 1)
         qc.ecr(0, 1)
         qc.iswap(0, 1)
-        res = twirl_circuit(qc, twirling_gate=["cx", iSwapGate], seed=12345)
+        res = twirl_circuit(qc, twirling_gate=["cx", iSwapGate()], seed=12345)
         np.testing.assert_allclose(
             Operator(qc), Operator(res), err_msg=f"{qc}\nnot equiv to\n{res}"
         )
@@ -99,14 +99,17 @@ def __init__(self):
         qc.append(MyGate(), (0, 1))
 
         with self.assertRaises(QiskitError):
-            twirl_circuit(qc, twirling_gate=MyGate)
+            twirl_circuit(qc, twirling_gate=MyGate())
 
-    def test_invalid_standard_gate(self):
+    def test_custom_standard_gate(self):
         """Test an error is raised with an unsupported standard gate."""
         qc = QuantumCircuit(2)
         qc.swap(0, 1)
-        with self.assertRaises(QiskitError):
-            twirl_circuit(qc, twirling_gate=SwapGate)
+        res = twirl_circuit(qc, twirling_gate=SwapGate())
+        np.testing.assert_allclose(
+            Operator(qc), Operator(res), err_msg=f"gate: {qc} not equiv to\n{res}"
+        )
+        self.assertNotEqual(qc, res)
 
     def test_invalid_string(self):
         """Test an error is raised with an unsupported standard gate."""
@@ -126,8 +129,11 @@ def test_invalid_class_entry_in_list(self):
         """Test an error is raised with an unsupported string gate in list."""
         qc = QuantumCircuit(2)
         qc.swap(0, 1)
-        with self.assertRaises(QiskitError):
-            twirl_circuit(qc, twirling_gate=[SwapGate, "cx"])
+        res = twirl_circuit(qc, twirling_gate=[SwapGate(), "cx"])
+        np.testing.assert_allclose(
+            Operator(qc), Operator(res), err_msg=f"gate: {qc} not equiv to\n{res}"
+        )
+        self.assertNotEqual(qc, res)
 
     @ddt.data(CXGate, ECRGate, CZGate, iSwapGate)
     def test_full_circuit(self, gate):