Skip to content

Commit 3e80326

Browse files
mtreinishmtreinish
authored
Use nalgebra::Matrix4 as output for instructions_to_matrix (#15871)
* Use nalgebra::Matrix4 as output for instructions_to_matrix This commit switches the return type of convert_2q_block_matrix::instructions_to_matrix() to return a nalgebra Matrix4 rather than a dynamicly allocated ndarray Array2. The function explicitly returns a 4x4 Complex64 matrix. Using an nalgebra fixed size array is stack allocated and we avoid needing a dynamic allocation. This should also speedup matrix multiplication since nalgebra can leverage simd better (either directly or implicitly via the compiler) because it knows the fixed operations needed. We were already setup towards doing since #13649 which moved to using nalgebra internally for the 1q component, but it didn't update the whole path in that PR to use an nalgebra array for everything. Ideally we'd be using nalgebra data types throughout the two qubit decomposers too. We should still use faer for the more involved linear algebra operations in the module but we should be using Matrix4 and Matrix2 for fixed sized matrices where we know the size of the matrices in that module and generate faer MatRefs using qiskit_synthesis::linalg::nalgebra_to_faer() to do linear algebra with the matrix. This PR is an incremental step towards doing that. * Update safety comment
1 parent 03e8f5b commit 3e80326

4 files changed

Lines changed: 148 additions & 50 deletions

File tree

‎crates/quantum_info/src/convert_2q_block_matrix.rs‎

Lines changed: 82 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -16,12 +16,12 @@ use pyo3::prelude::*;
1616

1717
use num_complex::Complex64;
1818
use numpy::PyReadonlyArray2;
19-
use numpy::ndarray::{Array2, ArrayView2, ArrayViewMut2, arr2, aview2};
19+
use numpy::nalgebra::{Matrix4, MatrixViewMut4};
20+
use numpy::ndarray::{Array2, ArrayView2};
2021
use rustworkx_core::petgraph::stable_graph::NodeIndex;
2122

2223
use qiskit_circuit::Qubit;
2324
use qiskit_circuit::dag_circuit::DAGCircuit;
24-
use qiskit_circuit::gate_matrix::TWO_QUBIT_IDENTITY;
2525
use qiskit_circuit::imports::QI_OPERATOR;
2626
use qiskit_circuit::interner::Interner;
2727
use qiskit_circuit::operations::{ArrayType, OperationRef};
@@ -30,6 +30,28 @@ use qiskit_circuit::packed_instruction::PackedInstruction;
3030
use crate::QiskitError;
3131
use crate::versor_u2::{VersorSU2, VersorU2, VersorU2Error};
3232

33+
#[inline]
34+
fn matrix4_from_pyreadonly(array: &PyReadonlyArray2<Complex64>) -> Matrix4<Complex64> {
35+
Matrix4::new(
36+
*array.get((0, 0)).unwrap(),
37+
*array.get((0, 1)).unwrap(),
38+
*array.get((0, 2)).unwrap(),
39+
*array.get((0, 3)).unwrap(),
40+
*array.get((1, 0)).unwrap(),
41+
*array.get((1, 1)).unwrap(),
42+
*array.get((1, 2)).unwrap(),
43+
*array.get((1, 3)).unwrap(),
44+
*array.get((2, 0)).unwrap(),
45+
*array.get((2, 1)).unwrap(),
46+
*array.get((2, 2)).unwrap(),
47+
*array.get((2, 3)).unwrap(),
48+
*array.get((3, 0)).unwrap(),
49+
*array.get((3, 1)).unwrap(),
50+
*array.get((3, 2)).unwrap(),
51+
*array.get((3, 3)).unwrap(),
52+
)
53+
}
54+
3355
#[inline]
3456
pub fn get_matrix_from_inst(inst: &PackedInstruction) -> PyResult<Array2<Complex64>> {
3557
if let Some(mat) = inst.try_matrix() {
@@ -59,6 +81,34 @@ pub fn get_matrix_from_inst(inst: &PackedInstruction) -> PyResult<Array2<Complex
5981
})
6082
}
6183

84+
#[inline]
85+
pub fn get_2q_matrix_from_inst(inst: &PackedInstruction) -> PyResult<Matrix4<Complex64>> {
86+
if let Some(mat) = inst.try_matrix_as_nalgebra_2q() {
87+
return Ok(mat);
88+
}
89+
if inst.op.try_standard_gate().is_some() {
90+
return Err(QiskitError::new_err(
91+
"Parameterized gates can't be consolidated",
92+
));
93+
}
94+
let OperationRef::Gate(gate) = inst.op.view() else {
95+
return Err(QiskitError::new_err(
96+
"Can't compute matrix of non-unitary op",
97+
));
98+
};
99+
// If the operation is a custom python gate, we will acquire the gil and use an
100+
// Operator. Otherwise, using op.matrix() should work. A user should not be
101+
// able to reach this condition in Rust standalone mode.
102+
Python::attach(|py| {
103+
let res = QI_OPERATOR
104+
.get_bound(py)
105+
.call1((gate.instruction.clone_ref(py),))?
106+
.getattr(intern!(py, "data"))?
107+
.extract()?;
108+
Ok(matrix4_from_pyreadonly(&res))
109+
})
110+
}
111+
62112
/// Quaternion-based collect of two parallel runs of 1q gates.
63113
#[derive(Clone, Debug)]
64114
struct Separable1q {
@@ -86,7 +136,7 @@ impl Separable1q {
86136

87137
/// Construct the two-qubit gate matrix implied by these two runs.
88138
#[inline]
89-
fn matrix_into<'a>(&self, target: &'a mut [[Complex64; 4]; 4]) -> ArrayView2<'a, Complex64> {
139+
fn matrix_into(&self, mut target: MatrixViewMut4<Complex64>) {
90140
let q0 = VersorU2 {
91141
phase: self.phase,
92142
su2: self.qubits[0],
@@ -98,17 +148,16 @@ impl Separable1q {
98148
let q1_row = [Complex64::new(scalar, iz), Complex64::new(iy, ix)];
99149
for out_row in 0..2 {
100150
for out_col in 0..4 {
101-
target[out_row][out_col] = q1_row[(out_col & 2) >> 1] * q0[out_row][out_col & 1];
151+
target[(out_row, out_col)] = q1_row[(out_col & 2) >> 1] * q0[out_row][out_col & 1];
102152
}
103153
}
104154
let q1_row = [Complex64::new(-iy, ix), Complex64::new(scalar, -iz)];
105155
for out_row in 0..2 {
106156
for out_col in 0..4 {
107-
target[out_row + 2][out_col] =
157+
target[(out_row + 2, out_col)] =
108158
q1_row[(out_col & 2) >> 1] * q0[out_row][out_col & 1];
109159
}
110160
}
111-
aview2(target)
112161
}
113162
}
114163

@@ -136,7 +185,7 @@ pub fn blocks_to_matrix(
136185
dag: &DAGCircuit,
137186
op_list: &[NodeIndex],
138187
block_index_map: [Qubit; 2],
139-
) -> PyResult<Array2<Complex64>> {
188+
) -> PyResult<Matrix4<Complex64>> {
140189
let inst_iter = op_list.iter().map(|node| dag[*node].unwrap_operation());
141190
instructions_to_matrix(inst_iter, block_index_map, dag.qargs_interner())
142191
}
@@ -145,7 +194,7 @@ pub fn instructions_to_matrix<'a>(
145194
op_list: impl Iterator<Item = &'a PackedInstruction>,
146195
block_index_map: [Qubit; 2],
147196
qargs_interner: &'a Interner<[Qubit]>,
148-
) -> PyResult<Array2<Complex64>> {
197+
) -> PyResult<Matrix4<Complex64>> {
149198
#[derive(Clone, Copy, PartialEq, Eq)]
150199
enum Qarg {
151200
Q0 = 0,
@@ -169,9 +218,9 @@ pub fn instructions_to_matrix<'a>(
169218
.expect("not a 1q or 2q gate in the qargs block")]
170219
};
171220

172-
let mut work: [[Complex64; 4]; 4] = Default::default();
221+
let mut work: Matrix4<Complex64> = Matrix4::zeros();
173222
let mut qubits_1q: Option<Separable1q> = None;
174-
let mut output_matrix: Option<Array2<Complex64>> = None;
223+
let mut output_matrix: Option<Matrix4<Complex64>> = None;
175224
for inst in op_list {
176225
let qarg = qarg_lookup(inst.qubits);
177226
match qarg {
@@ -183,31 +232,36 @@ pub fn instructions_to_matrix<'a>(
183232
};
184233
}
185234
Qarg::Q01 | Qarg::Q10 => {
186-
let mut matrix = get_matrix_from_inst(inst)?;
235+
let mut matrix = get_2q_matrix_from_inst(inst)?;
187236
if qarg == Qarg::Q10 {
188-
change_basis_inplace(matrix.view_mut());
237+
change_basis_inplace(matrix.as_view_mut());
189238
}
190239
if let Some(sep) = qubits_1q.take() {
191-
matrix = matrix.dot(&sep.matrix_into(&mut work));
240+
sep.matrix_into(work.as_view_mut());
241+
matrix *= work;
192242
}
193243
output_matrix = if let Some(state) = output_matrix {
194-
Some(matrix.dot(&state))
244+
Some(matrix * state)
195245
} else {
196246
Some(matrix)
197247
};
198248
}
199249
}
200250
}
251+
201252
match (
202-
qubits_1q.map(|sep| sep.matrix_into(&mut work)),
253+
qubits_1q.map(|sep| {
254+
sep.matrix_into(work.as_view_mut());
255+
work
256+
}),
203257
output_matrix,
204258
) {
205-
(Some(sep), Some(state)) => Ok(sep.dot(&state)),
259+
(Some(sep), Some(state)) => Ok(sep * state),
206260
(None, Some(state)) => Ok(state),
207-
(Some(sep), None) => Ok(sep.to_owned()),
261+
(Some(sep), None) => Ok(sep),
208262
// This shouldn't actually ever trigger, because we expect blocks to be non-empty, but it's
209263
// trivial to handle anyway.
210-
(None, None) => Ok(arr2(&TWO_QUBIT_IDENTITY)),
264+
(None, None) => Ok(Matrix4::identity()),
211265
}
212266
}
213267

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

228282
/// Change the qubit order of a 2q matrix in place.
229283
#[inline]
230-
pub fn change_basis_inplace(mut matrix: ArrayViewMut2<Complex64>) {
231-
matrix.swap((0, 1), (0, 2));
232-
matrix.swap((3, 1), (3, 2));
233-
matrix.swap((1, 0), (2, 0));
234-
matrix.swap((1, 3), (2, 3));
235-
matrix.swap((1, 1), (2, 2));
236-
matrix.swap((1, 2), (2, 1));
284+
pub fn change_basis_inplace(mut matrix: MatrixViewMut4<Complex64>) {
285+
// SAFETY: The bounds are all valid for a 4x4 matrix
286+
unsafe {
287+
matrix.swap_unchecked((0, 1), (0, 2));
288+
matrix.swap_unchecked((3, 1), (3, 2));
289+
matrix.swap_unchecked((1, 0), (2, 0));
290+
matrix.swap_unchecked((1, 3), (2, 3));
291+
matrix.swap_unchecked((1, 1), (2, 2));
292+
matrix.swap_unchecked((1, 2), (2, 1));
293+
}
237294
}

‎crates/synthesis/src/linalg/mod.rs‎

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@
1212

1313
use approx::{abs_diff_eq, relative_ne};
1414
use faer::MatRef;
15-
use nalgebra::{DMatrix, DMatrixView, Dim, Dyn, MatrixView, ViewStorage};
15+
use nalgebra::{DMatrix, DMatrixView, Dim, Dyn, MatrixView, Scalar, ViewStorage};
1616
use ndarray::ArrayView2;
1717
use ndarray::ShapeBuilder;
1818
use num_complex::Complex64;
@@ -22,6 +22,16 @@ pub mod cos_sin_decomp;
2222
const ATOL_DEFAULT: f64 = 1e-8;
2323
const RTOL_DEFAULT: f64 = 1e-5;
2424

25+
#[inline]
26+
pub fn nalgebra_array_view<T: Scalar, R: Dim, C: Dim>(mat: MatrixView<T, R, C>) -> ArrayView2<T> {
27+
let dim = ndarray::Dim(mat.shape());
28+
let strides = ndarray::Dim(mat.strides());
29+
// SAFETY: We know the array is a 2d array from nalgebra and we get the pointer and memory
30+
// layout description from nalgebra so we don't need to check for invalid format as nalgebra
31+
// has already validated this
32+
unsafe { ArrayView2::from_shape_ptr(dim.strides(strides), mat.get_unchecked(0)) }
33+
}
34+
2535
#[inline]
2636
pub fn ndarray_to_faer<T>(array: ArrayView2<'_, T>) -> MatRef<'_, T> {
2737
// This function's code is based on faer-ext's IntoFaer::into_faer implementation for ArrayView2:

‎crates/synthesis/src/qsd.rs‎

Lines changed: 21 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ use std::sync::OnceLock;
1515

1616
use approx::abs_diff_eq;
1717
use hashbrown::HashMap;
18-
use nalgebra::{DMatrix, DMatrixView, DVector, Matrix4, QR, stack};
18+
use nalgebra::{DMatrix, DMatrixView, DVector, Matrix4, QR, U4, stack};
1919
use ndarray::prelude::*;
2020
use num_complex::Complex64;
2121
use numpy::PyReadonlyArray2;
@@ -28,7 +28,9 @@ use thiserror::Error;
2828
use crate::euler_one_qubit_decomposer::{
2929
EulerBasis, EulerBasisSet, unitary_to_gate_sequence_inner,
3030
};
31-
use crate::linalg::{closest_unitary, is_hermitian_matrix, svd_decomposition, verify_unitary};
31+
use crate::linalg::{
32+
closest_unitary, is_hermitian_matrix, nalgebra_array_view, svd_decomposition, verify_unitary,
33+
};
3234
use crate::two_qubit_decompose::{TwoQubitBasisDecomposer, two_qubit_decompose_up_to_diagonal};
3335
use qiskit_circuit::bit::ShareableQubit;
3436
use qiskit_circuit::circuit_data::{CircuitData, CircuitDataError, PyCircuitData};
@@ -716,19 +718,26 @@ fn apply_a2(
716718
.collect();
717719
for ind in ind2q.windows(2) {
718720
let mat1 = match diagonal_rollover.get(&ind[0]) {
719-
Some(circ) => instructions_to_matrix(
720-
circ.data().iter(),
721-
[Qubit(0), Qubit(1)],
722-
circ.qargs_interner(),
723-
)?,
721+
Some(circ) => {
722+
let mat: Matrix4<Complex64> = instructions_to_matrix(
723+
circ.data().iter(),
724+
[Qubit(0), Qubit(1)],
725+
circ.qargs_interner(),
726+
)?;
727+
nalgebra_array_view::<Complex64, U4, U4>(mat.as_view()).to_owned()
728+
}
724729
None => new_matrices[&ind[0]].to_owned(),
725730
};
726731
let mat2 = match diagonal_rollover.get(&ind[1]) {
727-
Some(circ) => instructions_to_matrix(
728-
circ.data().iter(),
729-
[Qubit(0), Qubit(1)],
730-
circ.qargs_interner(),
731-
)?,
732+
Some(circ) => nalgebra_array_view::<Complex64, U4, U4>(
733+
instructions_to_matrix(
734+
circ.data().iter(),
735+
[Qubit(0), Qubit(1)],
736+
circ.qargs_interner(),
737+
)?
738+
.as_view(),
739+
)
740+
.to_owned(),
732741
None => new_matrices[&ind[1]].to_owned(),
733742
};
734743
let (diagonal_mat, qc2cx) = two_qubit_decompose_up_to_diagonal(mat1.view()).unwrap();

‎crates/transpiler/src/passes/consolidate_blocks.rs‎

Lines changed: 34 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@
1212

1313
use super::optimize_1q_gates_decomposition::matmul_1q;
1414
use hashbrown::{HashMap, HashSet};
15-
use nalgebra::Matrix2;
15+
use nalgebra::{Matrix2, Matrix4, U4};
1616
use ndarray::ArrayView2;
1717
use ndarray::{Array2, aview2};
1818
use num_complex::Complex64;
@@ -25,14 +25,14 @@ use qiskit_circuit::circuit_data::CircuitData;
2525
use qiskit_circuit::dag_circuit::{DAGCircuit, NodeType};
2626
use qiskit_circuit::gate_matrix::{
2727
CH_GATE, CX_GATE, CY_GATE, CZ_GATE, DCX_GATE, ECR_GATE, ISWAP_GATE, ONE_QUBIT_IDENTITY,
28-
TWO_QUBIT_IDENTITY,
2928
};
3029
use qiskit_circuit::imports::QI_OPERATOR;
3130
use qiskit_circuit::interner::Interned;
3231
use qiskit_circuit::operations::StandardGate;
3332
use qiskit_circuit::operations::{ArrayType, Operation, Param, UnitaryGate};
3433
use qiskit_circuit::packed_instruction::PackedOperation;
3534
use qiskit_quantum_info::convert_2q_block_matrix::{blocks_to_matrix, get_matrix_from_inst};
35+
use qiskit_synthesis::linalg::nalgebra_array_view;
3636
use qiskit_synthesis::two_qubit_decompose::RXXEquivalent;
3737
use qiskit_synthesis::two_qubit_decompose::{
3838
TwoQubitBasisDecomposer, TwoQubitControlledUDecomposer,
@@ -44,6 +44,29 @@ use crate::passes::unitary_synthesis::{PARAM_SET, TWO_QUBIT_BASIS_SET};
4444
use crate::target::{Qargs, Target};
4545
use qiskit_circuit::PhysicalQubit;
4646

47+
static IDENTITY_2Q: Matrix4<Complex64> = Matrix4::new(
48+
// Row 1
49+
Complex64::ONE,
50+
Complex64::ZERO,
51+
Complex64::ZERO,
52+
Complex64::ZERO,
53+
// Row 2
54+
Complex64::ZERO,
55+
Complex64::ONE,
56+
Complex64::ZERO,
57+
Complex64::ZERO,
58+
// Row 3
59+
Complex64::ZERO,
60+
Complex64::ZERO,
61+
Complex64::ONE,
62+
Complex64::ZERO,
63+
// Row 4
64+
Complex64::ZERO,
65+
Complex64::ZERO,
66+
Complex64::ZERO,
67+
Complex64::ONE,
68+
);
69+
4770
#[allow(clippy::large_enum_variant)]
4871
#[derive(Clone, Debug)]
4972
pub enum DecomposerType {
@@ -367,12 +390,13 @@ fn py_run_consolidate_blocks(
367390
let matrix = blocks_to_matrix(dag, &block, block_index_map).ok();
368391
if let Some(matrix) = matrix {
369392
let num_basis_gates = match decomposer {
370-
DecomposerType::TwoQubitBasis(ref decomp) => {
371-
decomp.num_basis_gates_inner(matrix.view())?
372-
}
373-
DecomposerType::TwoQubitControlledU(ref decomp) => {
374-
decomp.num_basis_gates_inner(matrix.view())?
375-
}
393+
DecomposerType::TwoQubitBasis(ref decomp) => decomp.num_basis_gates_inner(
394+
nalgebra_array_view::<Complex64, U4, U4>(matrix.as_view()),
395+
)?,
396+
DecomposerType::TwoQubitControlledU(ref decomp) => decomp
397+
.num_basis_gates_inner(nalgebra_array_view::<Complex64, U4, U4>(
398+
matrix.as_view(),
399+
))?,
376400
};
377401

378402
if force_consolidate
@@ -381,15 +405,13 @@ fn py_run_consolidate_blocks(
381405
|| (basis_gates.is_some() && outside_basis)
382406
|| (target.is_some() && outside_basis)
383407
{
384-
if approx::abs_diff_eq!(aview2(&TWO_QUBIT_IDENTITY), matrix) {
408+
if approx::abs_diff_eq!(IDENTITY_2Q, matrix) {
385409
for node in block {
386410
dag.remove_op_node(node);
387411
}
388412
} else {
389-
// TODO: Use Matrix4/ArrayType::TwoQ when we're using nalgebra
390-
// for consolidation
391413
let unitary_gate = UnitaryGate {
392-
array: ArrayType::NDArray(matrix),
414+
array: ArrayType::TwoQ(matrix),
393415
};
394416
let qubit_pos_map = block_index_map
395417
.into_iter()

0 commit comments

Comments
 (0)
⚡