Add unitary gate representation to rust

This commit expands the rust circuit data model to include a definition
of UnitaryGate. This is a more abstract operation than what we've
defined so far in Rust, a gate represented solely by it's unitary matrix
but during normal transpilation we create and interact with these
operations for optimization purposes so having a native representation
in rust is necessary to avoid calling to Python when working with them.

This introduces a new UnitaryGate struct which represents the unpacked
operation. It has 3 internal storage variants based on either an ndarray
arbitrary sized array, a 2x2 nalgebra fixed sized array, or a 4x4
nalgebra fixed sized array. From the python perspective these all look
the same, but being able to internally work with all 3 variants lets us
optimize the code paths for 1q and 2q unitary gates (which are by far
more common) and avoid one layer of pointer indirection.

When stored in a circuit the packed representation is just a pointer to
the actual UnitaryGate which we put in a Box. This is necessary because
the struct size is too large to fit in our compact representation of
operations. Even the ndarray which is a heap allocated type requires more
than our allotted space in a PackedOperation so we need to reduce it's
size by putting it in a Box.

The one major difference here from the previous python based
representation the unitary matrix was stored as an object type parameter
in the PackedInstruction.params field, which now it is stored internally
in the operation itself. There is arguably a behavior change around
this because it's no longer possible to mutate the array of a
UnitaryGate in place once it's inserted into the circuit. While doing
this was horribly unsound, because there was no guardrails for doing it
a release note is added because there is a small use case where it would
have worked and it wasn't explicitly documented.

Closes Qiskit#13272

Cargo.toml

@@ -20,6 +20,7 @@ indexmap.version = "2.7.1"
 hashbrown.version = "0.14.5"
 num-bigint = "0.4"
 num-complex = "0.4"
+nalgebra = "0.33"
 ndarray = "0.15"
 numpy = "0.23"
 smallvec = "1.13"

crates/accelerate/src/target_transpiler/mod.rs

@@ -783,6 +783,9 @@ impl Target {
                     OperationRef::Gate(gate) => gate.gate.clone_ref(py),
                     OperationRef::Instruction(instruction) => instruction.instruction.clone_ref(py),
                     OperationRef::Operation(operation) => operation.operation.clone_ref(py),
+                    OperationRef::Unitary(unitary) => unitary
+                        .create_py_op(py, &ExtraInstructionAttributes::default())?
+                        .into_any(),
                 },
                 TargetOperation::Variadic(op_cls) => op_cls.clone_ref(py),
             };

crates/circuit/Cargo.toml

@@ -20,10 +20,10 @@ bytemuck.workspace = true
 bitfield-struct.workspace = true
 num-complex.workspace = true
 ndarray.workspace = true
-numpy.workspace = true
 thiserror.workspace = true
 approx.workspace = true
 itertools.workspace = true
+nalgebra.workspace = true
 
 [dependencies.pyo3]
 workspace = true
@@ -41,6 +41,10 @@ features = ["rayon"]
 workspace = true
 features = ["union"]
 
+[dependencies.numpy]
+workspace = true
+features = ["nalgebra"]
+
 [features]
 cache_pygates = []
 

crates/circuit/src/circuit_instruction.rs

@@ -13,23 +13,24 @@
 #[cfg(feature = "cache_pygates")]
 use std::sync::OnceLock;
 
-use numpy::{IntoPyArray, PyArray2};
+use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2};
 use pyo3::basic::CompareOp;
 use pyo3::exceptions::{PyDeprecationWarning, PyTypeError};
 use pyo3::prelude::*;
 use pyo3::types::{PyBool, PyList, PyString, PyTuple, PyType};
 use pyo3::IntoPyObjectExt;
 use pyo3::{intern, PyObject, PyResult};
 
+use nalgebra::{MatrixView2, MatrixView4};
 use num_complex::Complex64;
 use smallvec::SmallVec;
 
 use crate::imports::{
     CONTROLLED_GATE, CONTROL_FLOW_OP, GATE, INSTRUCTION, OPERATION, WARNINGS_WARN,
 };
 use crate::operations::{
-    Operation, OperationRef, Param, PyGate, PyInstruction, PyOperation, StandardGate,
-    StandardInstruction, StandardInstructionType,
+    ArrayType, Operation, OperationRef, Param, PyGate, PyInstruction, PyOperation, StandardGate,
+    StandardInstruction, StandardInstructionType, UnitaryGate,
 };
 use crate::packed_instruction::PackedOperation;
 
@@ -341,6 +342,9 @@ impl CircuitInstruction {
             OperationRef::Gate(gate) => gate.gate.clone_ref(py),
             OperationRef::Instruction(instruction) => instruction.instruction.clone_ref(py),
             OperationRef::Operation(operation) => operation.operation.clone_ref(py),
+            OperationRef::Unitary(unitary) => {
+                unitary.create_py_op(py, &self.extra_attrs)?.into_any()
+            }
         };
 
         #[cfg(feature = "cache_pygates")]
@@ -762,6 +766,45 @@ impl<'py> FromPyObject<'py> for OperationFromPython {
             });
         }
 
+        // We need to check by name here to avoid a circular import during initial loading
+        if ob.getattr(intern!(py, "name"))?.extract::<String>()? == "unitary" {
+            let params = extract_params()?;
+            if let Param::Obj(data) = &params[0] {
+                let py_matrix: PyReadonlyArray2<Complex64> = data.extract(py)?;
+                let matrix: Option<MatrixView2<Complex64>> = py_matrix.try_as_matrix();
+                if let Some(x) = matrix {
+                    let unitary_gate = UnitaryGate {
+                        array: ArrayType::OneQ(x.into_owned()),
+                    };
+                    return Ok(OperationFromPython {
+                        operation: PackedOperation::from_unitary(Box::new(unitary_gate)),
+                        params: SmallVec::new(),
+                        extra_attrs: extract_extra()?,
+                    });
+                }
+                let matrix: Option<MatrixView4<Complex64>> = py_matrix.try_as_matrix();
+                if let Some(x) = matrix {
+                    let unitary_gate = UnitaryGate {
+                        array: ArrayType::TwoQ(x.into_owned()),
+                    };
+                    return Ok(OperationFromPython {
+                        operation: PackedOperation::from_unitary(Box::new(unitary_gate)),
+                        params: SmallVec::new(),
+                        extra_attrs: extract_extra()?,
+                    });
+                } else {
+                    let unitary_gate = UnitaryGate {
+                        array: ArrayType::NDArray(py_matrix.as_array().to_owned()),
+                    };
+                    return Ok(OperationFromPython {
+                        operation: PackedOperation::from_unitary(Box::new(unitary_gate)),
+                        params: SmallVec::new(),
+                        extra_attrs: extract_extra()?,
+                    });
+                };
+            }
+        }
+
         if ob_type.is_subclass(GATE.get_bound(py))? {
             let params = extract_params()?;
             let gate = Box::new(PyGate {

crates/circuit/src/dag_circuit.rs

@@ -13,6 +13,7 @@
 use std::hash::Hash;
 
 use ahash::RandomState;
+use approx::relative_eq;
 use smallvec::SmallVec;
 
 use crate::bit_data::BitData;
@@ -26,7 +27,7 @@ use crate::dot_utils::build_dot;
 use crate::error::DAGCircuitError;
 use crate::imports;
 use crate::interner::{Interned, Interner};
-use crate::operations::{Operation, OperationRef, Param, PyInstruction, StandardGate};
+use crate::operations::{ArrayType, Operation, OperationRef, Param, PyInstruction, StandardGate};
 use crate::packed_instruction::{PackedInstruction, PackedOperation};
 use crate::rustworkx_core_vnext::isomorphism;
 use crate::{BitType, Clbit, Qubit, TupleLikeArg};
@@ -2631,6 +2632,96 @@ def _format(operand):
                         | [OperationRef::Instruction(_), OperationRef::StandardInstruction(_)] => {
                             Ok(inst1.py_op_eq(py, inst2)? && check_args() && check_conditions()?)
                         }
+                        [OperationRef::Unitary(op_a), OperationRef::Unitary(op_b)] => {
+                            match [&op_a.array, &op_b.array] {
+                                [ArrayType::NDArray(a), ArrayType::NDArray(b)] => {
+                                    Ok(relative_eq!(a, b, max_relative = 1e-5, epsilon = 1e-8))
+                                }
+                                [ArrayType::NDArray(a), ArrayType::OneQ(b)] => {
+                                    if a.shape()[0] == 2 {
+                                        for i in 0..2 {
+                                            for j in 0..2 {
+                                                if !relative_eq!(
+                                                    a[[i, j]],
+                                                    b[(i, j)],
+                                                    max_relative = 1e-5,
+                                                    epsilon = 1e-8
+                                                ) {
+                                                    return Ok(false);
+                                                }
+                                            }
+                                        }
+                                        Ok(true)
+                                    } else {
+                                        Ok(false)
+                                    }
+                                }
+                                [ArrayType::NDArray(a), ArrayType::TwoQ(b)] => {
+                                    if a.shape()[0] == 4 {
+                                        for i in 0..4 {
+                                            for j in 0..4 {
+                                                if !relative_eq!(
+                                                    a[[i, j]],
+                                                    b[(i, j)],
+                                                    max_relative = 1e-5,
+                                                    epsilon = 1e-8
+                                                ) {
+                                                    return Ok(false);
+                                                }
+                                            }
+                                        }
+                                        Ok(true)
+                                    } else {
+                                        Ok(false)
+                                    }
+                                }
+                                [ArrayType::OneQ(a), ArrayType::NDArray(b)] => {
+                                    if b.shape()[0] == 2 {
+                                        for i in 0..2 {
+                                            for j in 0..2 {
+                                                if !relative_eq!(
+                                                    b[[i, j]],
+                                                    a[(i, j)],
+                                                    max_relative = 1e-5,
+                                                    epsilon = 1e-8
+                                                ) {
+                                                    return Ok(false);
+                                                }
+                                            }
+                                        }
+                                        Ok(true)
+                                    } else {
+                                        Ok(false)
+                                    }
+                                }
+                                [ArrayType::TwoQ(a), ArrayType::NDArray(b)] => {
+                                    if b.shape()[0] == 4 {
+                                        for i in 0..4 {
+                                            for j in 0..4 {
+                                                if !relative_eq!(
+                                                    b[[i, j]],
+                                                    a[(i, j)],
+                                                    max_relative = 1e-5,
+                                                    epsilon = 1e-8
+                                                ) {
+                                                    return Ok(false);
+                                                }
+                                            }
+                                        }
+                                        Ok(true)
+                                    } else {
+                                        Ok(false)
+                                    }
+                                }
+                                [ArrayType::OneQ(a), ArrayType::OneQ(b)] => {
+                                    Ok(relative_eq!(a, b, max_relative = 1e-5, epsilon = 1e-8))
+                                }
+                                [ArrayType::TwoQ(a), ArrayType::TwoQ(b)] => {
+                                    Ok(relative_eq!(a, b, max_relative = 1e-5, epsilon = 1e-8))
+                                }
+                                _ => Ok(false),
+                            }
+                        }
                         _ => Ok(false),
                     }
                 }
@@ -3295,7 +3386,8 @@ def _format(operand):
                                     py_op.operation.setattr(py, "condition", new_condition)?;
                                 }
                                 OperationRef::StandardGate(_)
-                                | OperationRef::StandardInstruction(_) => {}
+                                | OperationRef::StandardInstruction(_)
+                                | OperationRef::Unitary(_) => {}
                             }
                         }
                     }
@@ -6245,9 +6337,9 @@ impl DAGCircuit {
             };
             #[cfg(feature = "cache_pygates")]
             let py_op = match new_op.operation.view() {
-                OperationRef::StandardGate(_) | OperationRef::StandardInstruction(_) => {
-                    OnceLock::new()
-                }
+                OperationRef::StandardGate(_)
+                | OperationRef::StandardInstruction(_)
+                | OperationRef::Unitary(_) => OnceLock::new(),
                 OperationRef::Gate(gate) => OnceLock::from(gate.gate.clone_ref(py)),
                 OperationRef::Instruction(instruction) => {
                     OnceLock::from(instruction.instruction.clone_ref(py))