Skip to content

Use nalgebra::Matrix4 as output for instructions_to_matrix #15871

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
alexanderivrii merged 2 commits into from
Mar 25, 2026
Merged

Use nalgebra::Matrix4 as output for instructions_to_matrix #15871

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
35f8d27
Use nalgebra::Matrix4 as output for instructions_to_matrix
mtreinish Mar 24, 2026
5ed8f82
Update safety comment
mtreinish Mar 25, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
107 changes: 82 additions & 25 deletions crates/quantum_info/src/convert_2q_block_matrix.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,12 +16,12 @@ use pyo3::prelude::*;

use num_complex::Complex64;
use numpy::PyReadonlyArray2;
use numpy::ndarray::{Array2, ArrayView2, ArrayViewMut2, arr2, aview2};
use numpy::nalgebra::{Matrix4, MatrixViewMut4};
use numpy::ndarray::{Array2, ArrayView2};
use rustworkx_core::petgraph::stable_graph::NodeIndex;

use qiskit_circuit::Qubit;
use qiskit_circuit::dag_circuit::DAGCircuit;
use qiskit_circuit::gate_matrix::TWO_QUBIT_IDENTITY;
use qiskit_circuit::imports::QI_OPERATOR;
use qiskit_circuit::interner::Interner;
use qiskit_circuit::operations::{ArrayType, OperationRef};
Expand All @@ -30,6 +30,28 @@ use qiskit_circuit::packed_instruction::PackedInstruction;
use crate::QiskitError;
use crate::versor_u2::{VersorSU2, VersorU2, VersorU2Error};

#[inline]
fn matrix4_from_pyreadonly(array: &PyReadonlyArray2<Complex64>) -> Matrix4<Complex64> {
Matrix4::new(
*array.get((0, 0)).unwrap(),
*array.get((0, 1)).unwrap(),
*array.get((0, 2)).unwrap(),
*array.get((0, 3)).unwrap(),
*array.get((1, 0)).unwrap(),
*array.get((1, 1)).unwrap(),
*array.get((1, 2)).unwrap(),
*array.get((1, 3)).unwrap(),
*array.get((2, 0)).unwrap(),
*array.get((2, 1)).unwrap(),
*array.get((2, 2)).unwrap(),
*array.get((2, 3)).unwrap(),
*array.get((3, 0)).unwrap(),
*array.get((3, 1)).unwrap(),
*array.get((3, 2)).unwrap(),
*array.get((3, 3)).unwrap(),
)
}

#[inline]
pub fn get_matrix_from_inst(inst: &PackedInstruction) -> PyResult<Array2<Complex64>> {
if let Some(mat) = inst.try_matrix() {
Expand Down Expand Up @@ -59,6 +81,34 @@ pub fn get_matrix_from_inst(inst: &PackedInstruction) -> PyResult<Array2<Complex
})
}

#[inline]
pub fn get_2q_matrix_from_inst(inst: &PackedInstruction) -> PyResult<Matrix4<Complex64>> {
if let Some(mat) = inst.try_matrix_as_nalgebra_2q() {
return Ok(mat);
}
if inst.op.try_standard_gate().is_some() {
return Err(QiskitError::new_err(
"Parameterized gates can't be consolidated",
));
}
let OperationRef::Gate(gate) = inst.op.view() else {
return Err(QiskitError::new_err(
"Can't compute matrix of non-unitary op",
));
};
Comment on lines +94 to +98
Copy link
Copy Markdown
Member

@alexanderivrii alexanderivrii Mar 25, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we also consider OperationRef::Operation variant here (and if so, then also in get_matrix_from_inst)?

We should also extend this to PauliProductRotation, but this can be done for all relevant functions in a separate PR.

Copy link
Copy Markdown
Member Author

@mtreinish mtreinish Mar 25, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just mirrored what was already there. I think we do want to expand this for any unitary operation. So this will be PPR and Operations that have a unitary matrix. Ideally we'd probably cover PPR and unitary OperationRef::Operations in PackedInstruction::try_matrix_as_nalgebra_2q this check is just the fallback for whether we call out to python and use quantum_info.Operator to compute the unitary from a gate's definition.

Copy link
Copy Markdown
Member Author

@mtreinish mtreinish Mar 25, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That being said I think we should fix this in a separate PR that fixes all the matrix methods

// If the operation is a custom python gate, we will acquire the gil and use an
// Operator. Otherwise, using op.matrix() should work. A user should not be
// able to reach this condition in Rust standalone mode.
Python::attach(|py| {
let res = QI_OPERATOR
.get_bound(py)
.call1((gate.instruction.clone_ref(py),))?
.getattr(intern!(py, "data"))?
.extract()?;
Ok(matrix4_from_pyreadonly(&res))
})
}

/// Quaternion-based collect of two parallel runs of 1q gates.
#[derive(Clone, Debug)]
struct Separable1q {
Expand Down Expand Up @@ -86,7 +136,7 @@ impl Separable1q {

/// Construct the two-qubit gate matrix implied by these two runs.
#[inline]
fn matrix_into<'a>(&self, target: &'a mut [[Complex64; 4]; 4]) -> ArrayView2<'a, Complex64> {
fn matrix_into(&self, mut target: MatrixViewMut4<Complex64>) {
let q0 = VersorU2 {
phase: self.phase,
su2: self.qubits[0],
Expand All @@ -98,17 +148,16 @@ impl Separable1q {
let q1_row = [Complex64::new(scalar, iz), Complex64::new(iy, ix)];
for out_row in 0..2 {
for out_col in 0..4 {
target[out_row][out_col] = q1_row[(out_col & 2) >> 1] * q0[out_row][out_col & 1];
target[(out_row, out_col)] = q1_row[(out_col & 2) >> 1] * q0[out_row][out_col & 1];
}
}
let q1_row = [Complex64::new(-iy, ix), Complex64::new(scalar, -iz)];
for out_row in 0..2 {
for out_col in 0..4 {
target[out_row + 2][out_col] =
target[(out_row + 2, out_col)] =
q1_row[(out_col & 2) >> 1] * q0[out_row][out_col & 1];
}
}
aview2(target)
}
}

Expand Down Expand Up @@ -136,7 +185,7 @@ pub fn blocks_to_matrix(
dag: &DAGCircuit,
op_list: &[NodeIndex],
block_index_map: [Qubit; 2],
) -> PyResult<Array2<Complex64>> {
) -> PyResult<Matrix4<Complex64>> {
let inst_iter = op_list.iter().map(|node| dag[*node].unwrap_operation());
instructions_to_matrix(inst_iter, block_index_map, dag.qargs_interner())
}
Expand All @@ -145,7 +194,7 @@ pub fn instructions_to_matrix<'a>(
op_list: impl Iterator<Item = &'a PackedInstruction>,
block_index_map: [Qubit; 2],
qargs_interner: &'a Interner<[Qubit]>,
) -> PyResult<Array2<Complex64>> {
) -> PyResult<Matrix4<Complex64>> {
#[derive(Clone, Copy, PartialEq, Eq)]
enum Qarg {
Q0 = 0,
Expand All @@ -169,9 +218,9 @@ pub fn instructions_to_matrix<'a>(
.expect("not a 1q or 2q gate in the qargs block")]
};

let mut work: [[Complex64; 4]; 4] = Default::default();
let mut work: Matrix4<Complex64> = Matrix4::zeros();
let mut qubits_1q: Option<Separable1q> = None;
let mut output_matrix: Option<Array2<Complex64>> = None;
let mut output_matrix: Option<Matrix4<Complex64>> = None;
for inst in op_list {
let qarg = qarg_lookup(inst.qubits);
match qarg {
Expand All @@ -183,31 +232,36 @@ pub fn instructions_to_matrix<'a>(
};
}
Qarg::Q01 | Qarg::Q10 => {
let mut matrix = get_matrix_from_inst(inst)?;
let mut matrix = get_2q_matrix_from_inst(inst)?;
if qarg == Qarg::Q10 {
change_basis_inplace(matrix.view_mut());
change_basis_inplace(matrix.as_view_mut());
}
if let Some(sep) = qubits_1q.take() {
matrix = matrix.dot(&sep.matrix_into(&mut work));
sep.matrix_into(work.as_view_mut());
matrix *= work;
}
output_matrix = if let Some(state) = output_matrix {
Some(matrix.dot(&state))
Some(matrix * state)
} else {
Some(matrix)
};
}
}
}

match (
qubits_1q.map(|sep| sep.matrix_into(&mut work)),
qubits_1q.map(|sep| {
sep.matrix_into(work.as_view_mut());
work
}),
output_matrix,
) {
(Some(sep), Some(state)) => Ok(sep.dot(&state)),
(Some(sep), Some(state)) => Ok(sep * state),
(None, Some(state)) => Ok(state),
(Some(sep), None) => Ok(sep.to_owned()),
(Some(sep), None) => Ok(sep),
// This shouldn't actually ever trigger, because we expect blocks to be non-empty, but it's
// trivial to handle anyway.
(None, None) => Ok(arr2(&TWO_QUBIT_IDENTITY)),
(None, None) => Ok(Matrix4::identity()),
}
}

Expand All @@ -227,11 +281,14 @@ pub fn change_basis(matrix: ArrayView2<Complex64>) -> Array2<Complex64> {

/// Change the qubit order of a 2q matrix in place.
#[inline]
pub fn change_basis_inplace(mut matrix: ArrayViewMut2<Complex64>) {
matrix.swap((0, 1), (0, 2));
matrix.swap((3, 1), (3, 2));
matrix.swap((1, 0), (2, 0));
matrix.swap((1, 3), (2, 3));
matrix.swap((1, 1), (2, 2));
matrix.swap((1, 2), (2, 1));
pub fn change_basis_inplace(mut matrix: MatrixViewMut4<Complex64>) {
// SAFETY: The bounds are all valid for a 4x4 matrix
unsafe {
matrix.swap_unchecked((0, 1), (0, 2));
matrix.swap_unchecked((3, 1), (3, 2));
matrix.swap_unchecked((1, 0), (2, 0));
matrix.swap_unchecked((1, 3), (2, 3));
matrix.swap_unchecked((1, 1), (2, 2));
matrix.swap_unchecked((1, 2), (2, 1));
}
}
12 changes: 11 additions & 1 deletion crates/synthesis/src/linalg/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@

use approx::{abs_diff_eq, relative_ne};
use faer::MatRef;
use nalgebra::{DMatrix, DMatrixView, Dim, Dyn, MatrixView, ViewStorage};
use nalgebra::{DMatrix, DMatrixView, Dim, Dyn, MatrixView, Scalar, ViewStorage};
use ndarray::ArrayView2;
use ndarray::ShapeBuilder;
use num_complex::Complex64;
Expand All @@ -22,6 +22,16 @@ pub mod cos_sin_decomp;
const ATOL_DEFAULT: f64 = 1e-8;
const RTOL_DEFAULT: f64 = 1e-5;

#[inline]
pub fn nalgebra_array_view<T: Scalar, R: Dim, C: Dim>(mat: MatrixView<T, R, C>) -> ArrayView2<T> {
let dim = ndarray::Dim(mat.shape());
let strides = ndarray::Dim(mat.strides());
// SAFETY: We know the array is a 2d array from nalgebra and we get the pointer and memory
// layout description from nalgebra so we don't need to check for invalid format as nalgebra
// has already validated this
unsafe { ArrayView2::from_shape_ptr(dim.strides(strides), mat.get_unchecked(0)) }
}

#[inline]
pub fn ndarray_to_faer<T>(array: ArrayView2<'_, T>) -> MatRef<'_, T> {
// This function's code is based on faer-ext's IntoFaer::into_faer implementation for ArrayView2:
Expand Down
33 changes: 21 additions & 12 deletions crates/synthesis/src/qsd.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ use std::sync::OnceLock;

use approx::abs_diff_eq;
use hashbrown::HashMap;
use nalgebra::{DMatrix, DMatrixView, DVector, Matrix4, QR, stack};
use nalgebra::{DMatrix, DMatrixView, DVector, Matrix4, QR, U4, stack};
use ndarray::prelude::*;
use num_complex::Complex64;
use numpy::PyReadonlyArray2;
Expand All @@ -28,7 +28,9 @@ use thiserror::Error;
use crate::euler_one_qubit_decomposer::{
EulerBasis, EulerBasisSet, unitary_to_gate_sequence_inner,
};
use crate::linalg::{closest_unitary, is_hermitian_matrix, svd_decomposition, verify_unitary};
use crate::linalg::{
closest_unitary, is_hermitian_matrix, nalgebra_array_view, svd_decomposition, verify_unitary,
};
use crate::two_qubit_decompose::{TwoQubitBasisDecomposer, two_qubit_decompose_up_to_diagonal};
use qiskit_circuit::bit::ShareableQubit;
use qiskit_circuit::circuit_data::{CircuitData, CircuitDataError, PyCircuitData};
Expand Down Expand Up @@ -716,19 +718,26 @@ fn apply_a2(
.collect();
for ind in ind2q.windows(2) {
let mat1 = match diagonal_rollover.get(&ind[0]) {
Some(circ) => instructions_to_matrix(
circ.data().iter(),
[Qubit(0), Qubit(1)],
circ.qargs_interner(),
)?,
Some(circ) => {
let mat: Matrix4<Complex64> = instructions_to_matrix(
circ.data().iter(),
[Qubit(0), Qubit(1)],
circ.qargs_interner(),
)?;
nalgebra_array_view::<Complex64, U4, U4>(mat.as_view()).to_owned()
}
None => new_matrices[&ind[0]].to_owned(),
};
let mat2 = match diagonal_rollover.get(&ind[1]) {
Some(circ) => instructions_to_matrix(
circ.data().iter(),
[Qubit(0), Qubit(1)],
circ.qargs_interner(),
)?,
Some(circ) => nalgebra_array_view::<Complex64, U4, U4>(
instructions_to_matrix(
circ.data().iter(),
[Qubit(0), Qubit(1)],
circ.qargs_interner(),
)?
.as_view(),
)
.to_owned(),
None => new_matrices[&ind[1]].to_owned(),
};
let (diagonal_mat, qc2cx) = two_qubit_decompose_up_to_diagonal(mat1.view()).unwrap();
Expand Down
46 changes: 34 additions & 12 deletions crates/transpiler/src/passes/consolidate_blocks.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@

use super::optimize_1q_gates_decomposition::matmul_1q;
use hashbrown::{HashMap, HashSet};
use nalgebra::Matrix2;
use nalgebra::{Matrix2, Matrix4, U4};
use ndarray::ArrayView2;
use ndarray::{Array2, aview2};
use num_complex::Complex64;
Expand All @@ -25,14 +25,14 @@ use qiskit_circuit::circuit_data::CircuitData;
use qiskit_circuit::dag_circuit::{DAGCircuit, NodeType};
use qiskit_circuit::gate_matrix::{
CH_GATE, CX_GATE, CY_GATE, CZ_GATE, DCX_GATE, ECR_GATE, ISWAP_GATE, ONE_QUBIT_IDENTITY,
TWO_QUBIT_IDENTITY,
};
use qiskit_circuit::imports::QI_OPERATOR;
use qiskit_circuit::interner::Interned;
use qiskit_circuit::operations::StandardGate;
use qiskit_circuit::operations::{ArrayType, Operation, Param, UnitaryGate};
use qiskit_circuit::packed_instruction::PackedOperation;
use qiskit_quantum_info::convert_2q_block_matrix::{blocks_to_matrix, get_matrix_from_inst};
use qiskit_synthesis::linalg::nalgebra_array_view;
use qiskit_synthesis::two_qubit_decompose::RXXEquivalent;
use qiskit_synthesis::two_qubit_decompose::{
TwoQubitBasisDecomposer, TwoQubitControlledUDecomposer,
Expand All @@ -44,6 +44,29 @@ use crate::passes::unitary_synthesis::{PARAM_SET, TWO_QUBIT_BASIS_SET};
use crate::target::{Qargs, Target};
use qiskit_circuit::PhysicalQubit;

static IDENTITY_2Q: Matrix4<Complex64> = Matrix4::new(
// Row 1
Complex64::ONE,
Complex64::ZERO,
Complex64::ZERO,
Complex64::ZERO,
// Row 2
Complex64::ZERO,
Complex64::ONE,
Complex64::ZERO,
Complex64::ZERO,
// Row 3
Complex64::ZERO,
Complex64::ZERO,
Complex64::ONE,
Complex64::ZERO,
// Row 4
Complex64::ZERO,
Complex64::ZERO,
Complex64::ZERO,
Complex64::ONE,
);
Comment on lines +47 to +68
Copy link
Copy Markdown
Member

@alexanderivrii alexanderivrii Mar 25, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like this matrix: no need to think whether the elements are passed in row-major or column-major order 😄

Copy link
Copy Markdown
Member Author

@mtreinish mtreinish Mar 25, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ideally I would have been able to use Matrix4::identity() here, but it's not a const function so I couldn't use it in a static context. This was the best way I could come up with, but I also do a double take every time I work with nalgebra about row major vs col major because I've gotten it backwards too many times.


#[allow(clippy::large_enum_variant)]
#[derive(Clone, Debug)]
pub enum DecomposerType {
Expand Down Expand Up @@ -367,12 +390,13 @@ fn py_run_consolidate_blocks(
let matrix = blocks_to_matrix(dag, &block, block_index_map).ok();
if let Some(matrix) = matrix {
let num_basis_gates = match decomposer {
DecomposerType::TwoQubitBasis(ref decomp) => {
decomp.num_basis_gates_inner(matrix.view())?
}
DecomposerType::TwoQubitControlledU(ref decomp) => {
decomp.num_basis_gates_inner(matrix.view())?
}
DecomposerType::TwoQubitBasis(ref decomp) => decomp.num_basis_gates_inner(
nalgebra_array_view::<Complex64, U4, U4>(matrix.as_view()),
)?,
DecomposerType::TwoQubitControlledU(ref decomp) => decomp
.num_basis_gates_inner(nalgebra_array_view::<Complex64, U4, U4>(
matrix.as_view(),
))?,
};

if force_consolidate
Expand All @@ -381,15 +405,13 @@ fn py_run_consolidate_blocks(
|| (basis_gates.is_some() && outside_basis)
|| (target.is_some() && outside_basis)
{
if approx::abs_diff_eq!(aview2(&TWO_QUBIT_IDENTITY), matrix) {
if approx::abs_diff_eq!(IDENTITY_2Q, matrix) {
for node in block {
dag.remove_op_node(node);
}
} else {
// TODO: Use Matrix4/ArrayType::TwoQ when we're using nalgebra
// for consolidation
let unitary_gate = UnitaryGate {
array: ArrayType::NDArray(matrix),
array: ArrayType::TwoQ(matrix),
};
let qubit_pos_map = block_index_map
.into_iter()
Expand Down
⚡