From 5b7c6b619910ca9bf35eb87c213bab7cb442a255 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 23 May 2024 10:00:56 -0500 Subject: [PATCH 01/43] gaussian elimination in rust --- crates/accelerate/src/lib.rs | 1 + crates/accelerate/src/linear_matrix.rs | 88 ++++++++++++++++++++++++++ 2 files changed, 89 insertions(+) create mode 100644 crates/accelerate/src/linear_matrix.rs diff --git a/crates/accelerate/src/lib.rs b/crates/accelerate/src/lib.rs index 0af8ea6a0fce..ede89b44f4bd 100644 --- a/crates/accelerate/src/lib.rs +++ b/crates/accelerate/src/lib.rs @@ -20,6 +20,7 @@ pub mod edge_collections; pub mod error_map; pub mod euler_one_qubit_decomposer; pub mod isometry; +pub mod linear_matrix; pub mod nlayout; pub mod optimize_1q_gates; pub mod pauli_exp_val; diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs new file mode 100644 index 000000000000..2737cc405a18 --- /dev/null +++ b/crates/accelerate/src/linear_matrix.rs @@ -0,0 +1,88 @@ +// This code is part of Qiskit. +// +// (C) Copyright IBM 2024 +// +// This code is licensed under the Apache License, Version 2.0. You may +// obtain a copy of this license in the LICENSE.txt file in the root directory +// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +// +// 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 ndarray::{s, Array2}; + +// Perform ROW operation on a matrix mat +fn _row_op(mat: &mut Array2, ctrl: usize, trgt: usize) { + let row0 = mat.row(ctrl).to_owned(); + let mut row1 = mat.row_mut(trgt); + row1.zip_mut_with(&row0, |x, &y| *x ^= y); +} + +// Perform COL operation on a matrix mat +fn _col_op(mat: &mut Array2, ctrl: usize, trgt: usize) { + let col0 = mat.column(ctrl).to_owned(); + let mut col1 = mat.column_mut(trgt); + col1.zip_mut_with(&col0, |x, &y| *x ^= y); +} + +// Gauss elimination of a matrix mat with m rows and n columns. +// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] +// Returns the matrix mat. + +fn _gauss_elimination( + mut mat: Array2, + ncols: Option, + full_elim: Option, +) -> Array2 { + let (m, mut n) = (mat.nrows(), mat.ncols()); // no. of rows and columns + if let Some(ncols_val) = ncols { + n = usize::min(n, ncols_val); // no. of active columns + } + + let mut r = 0; // current rank + let mut k = 0; // current pivot column + while (r < m) && (k < n) { + let mut is_non_zero = false; + let mut new_r = r; + for j in k..n { + for i in r..m { + if mat[(i, j)] { + is_non_zero = true; + k = j; + new_r = i; + break; + } + } + if is_non_zero { + break; + } + } + if !is_non_zero { + return mat; // A is in the canonical form + } + + if new_r != r { + let temp_r = mat.slice_mut(s![r, ..]).to_owned(); + let temp_new_r = mat.slice_mut(s![new_r, ..]).to_owned(); + mat.slice_mut(s![r, ..]).assign(&temp_new_r); + mat.slice_mut(s![new_r, ..]).assign(&temp_r); + } + + if full_elim.is_some() { + for i in 0..r { + if mat[(i, k)] { + _row_op(&mut mat, r, i); + } + } + } + + for i in r + 1..m { + if mat[(i, k)] { + _row_op(&mut mat, r, i); + } + } + r = r + 1; + } + mat +} From 10633a25efac8ce6505246428a0a4a0d17e3702c Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Sun, 26 May 2024 05:51:18 -0500 Subject: [PATCH 02/43] handle lint errors --- crates/accelerate/src/linear_matrix.rs | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 2737cc405a18..7aea5791e33e 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -41,15 +41,17 @@ fn _gauss_elimination( } let mut r = 0; // current rank - let mut k = 0; // current pivot column + let k = 0; // current pivot column + let mut new_k = 0; while (r < m) && (k < n) { let mut is_non_zero = false; let mut new_r = r; for j in k..n { + new_k = k; for i in r..m { if mat[(i, j)] { is_non_zero = true; - k = j; + new_k = j; new_r = i; break; } @@ -71,18 +73,18 @@ fn _gauss_elimination( if full_elim.is_some() { for i in 0..r { - if mat[(i, k)] { + if mat[(i, new_k)] { _row_op(&mut mat, r, i); } } } for i in r + 1..m { - if mat[(i, k)] { + if mat[(i, new_k)] { _row_op(&mut mat, r, i); } } - r = r + 1; + r += 1; } mat } From 352ec0b7fcfcb987847db83e904e58c3144841f7 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 30 May 2024 06:34:57 -0500 Subject: [PATCH 03/43] replace python function by rust function for gauss elimination --- crates/accelerate/src/linear_matrix.rs | 25 +++++++++++++++++-- crates/pyext/src/lib.rs | 11 ++++---- qiskit/__init__.py | 1 + .../synthesis/linear/linear_matrix_utils.py | 12 +++------ 4 files changed, 33 insertions(+), 16 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 7aea5791e33e..af9b7ef75da0 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -11,6 +11,8 @@ // that they have been altered from the originals. use ndarray::{s, Array2}; +use numpy::{AllowTypeChange, IntoPyArray, PyArray2, PyArrayLike2}; +use pyo3::prelude::*; // Perform ROW operation on a matrix mat fn _row_op(mat: &mut Array2, ctrl: usize, trgt: usize) { @@ -29,8 +31,7 @@ fn _col_op(mat: &mut Array2, ctrl: usize, trgt: usize) { // Gauss elimination of a matrix mat with m rows and n columns. // If full_elim = True, it allows full elimination of mat[:, 0 : ncols] // Returns the matrix mat. - -fn _gauss_elimination( +fn gauss_elimination( mut mat: Array2, ncols: Option, full_elim: Option, @@ -88,3 +89,23 @@ fn _gauss_elimination( } mat } + +#[pyfunction] +#[pyo3(signature = (mat, ncols=None, full_elim=false))] +fn _gauss_elimination( + py: Python, + mat: PyArrayLike2, + ncols: Option, + full_elim: Option, +) -> PyResult>> { + let view = mat.as_array().to_owned(); + Ok(gauss_elimination(view, ncols, full_elim) + .into_pyarray_bound(py) + .unbind()) +} + +#[pymodule] +pub fn linear_matrix(m: &Bound) -> PyResult<()> { + m.add_wrapped(wrap_pyfunction!(_gauss_elimination))?; + Ok(()) +} diff --git a/crates/pyext/src/lib.rs b/crates/pyext/src/lib.rs index a21b1307a88f..085b598421e3 100644 --- a/crates/pyext/src/lib.rs +++ b/crates/pyext/src/lib.rs @@ -16,11 +16,11 @@ use pyo3::wrap_pymodule; use qiskit_accelerate::{ convert_2q_block_matrix::convert_2q_block_matrix, dense_layout::dense_layout, error_map::error_map, euler_one_qubit_decomposer::euler_one_qubit_decomposer, - isometry::isometry, nlayout::nlayout, optimize_1q_gates::optimize_1q_gates, - pauli_exp_val::pauli_expval, results::results, sabre::sabre, sampled_exp_val::sampled_exp_val, - sparse_pauli_op::sparse_pauli_op, stochastic_swap::stochastic_swap, - two_qubit_decompose::two_qubit_decompose, uc_gate::uc_gate, utils::utils, - vf2_layout::vf2_layout, + isometry::isometry, linear_matrix::linear_matrix, nlayout::nlayout, + optimize_1q_gates::optimize_1q_gates, pauli_exp_val::pauli_expval, results::results, + sabre::sabre, sampled_exp_val::sampled_exp_val, sparse_pauli_op::sparse_pauli_op, + stochastic_swap::stochastic_swap, two_qubit_decompose::two_qubit_decompose, uc_gate::uc_gate, + utils::utils, vf2_layout::vf2_layout, }; #[pymodule] @@ -33,6 +33,7 @@ fn _accelerate(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pymodule!(error_map))?; m.add_wrapped(wrap_pymodule!(euler_one_qubit_decomposer))?; m.add_wrapped(wrap_pymodule!(isometry))?; + m.add_wrapped(wrap_pymodule!(linear_matrix))?; m.add_wrapped(wrap_pymodule!(nlayout))?; m.add_wrapped(wrap_pymodule!(optimize_1q_gates))?; m.add_wrapped(wrap_pymodule!(pauli_expval))?; diff --git a/qiskit/__init__.py b/qiskit/__init__.py index e4fbc1729e53..65ed57e0707e 100644 --- a/qiskit/__init__.py +++ b/qiskit/__init__.py @@ -70,6 +70,7 @@ sys.modules["qiskit._accelerate.euler_one_qubit_decomposer"] = ( qiskit._accelerate.euler_one_qubit_decomposer ) +sys.modules["qiskit._accelerate.linear_matrix"] = qiskit._accelerate.linear_matrix sys.modules["qiskit._accelerate.nlayout"] = qiskit._accelerate.nlayout sys.modules["qiskit._accelerate.optimize_1q_gates"] = qiskit._accelerate.optimize_1q_gates sys.modules["qiskit._accelerate.pauli_expval"] = qiskit._accelerate.pauli_expval diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index 7a5b60641475..afce1f3bf510 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -16,6 +16,9 @@ import numpy as np from qiskit.exceptions import QiskitError +# pylint: disable=unused-import +from qiskit._accelerate.linear_matrix import _gauss_elimination + def check_invertible_binary_matrix(mat: np.ndarray): """Check that a binary matrix is invertible. @@ -57,15 +60,6 @@ def random_invertible_binary_matrix( return mat -def _gauss_elimination(mat, ncols=None, full_elim=False): - """Gauss elimination of a matrix mat with m rows and n columns. - If full_elim = True, it allows full elimination of mat[:, 0 : ncols] - Returns the matrix mat.""" - - mat, _ = _gauss_elimination_with_perm(mat, ncols, full_elim) - return mat - - def _gauss_elimination_with_perm(mat, ncols=None, full_elim=False): """Gauss elimination of a matrix mat with m rows and n columns. If full_elim = True, it allows full elimination of mat[:, 0 : ncols] From 5b17b5870f868dd509e6d7f47af0901fed268861 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Sun, 2 Jun 2024 05:25:06 -0500 Subject: [PATCH 04/43] change matrix elements type from bool to i8 --- crates/accelerate/src/linear_matrix.rs | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index af9b7ef75da0..d708e71bd108 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -15,14 +15,14 @@ use numpy::{AllowTypeChange, IntoPyArray, PyArray2, PyArrayLike2}; use pyo3::prelude::*; // Perform ROW operation on a matrix mat -fn _row_op(mat: &mut Array2, ctrl: usize, trgt: usize) { +fn _row_op(mat: &mut Array2, ctrl: usize, trgt: usize) { let row0 = mat.row(ctrl).to_owned(); let mut row1 = mat.row_mut(trgt); row1.zip_mut_with(&row0, |x, &y| *x ^= y); } // Perform COL operation on a matrix mat -fn _col_op(mat: &mut Array2, ctrl: usize, trgt: usize) { +fn _col_op(mat: &mut Array2, ctrl: usize, trgt: usize) { let col0 = mat.column(ctrl).to_owned(); let mut col1 = mat.column_mut(trgt); col1.zip_mut_with(&col0, |x, &y| *x ^= y); @@ -32,10 +32,10 @@ fn _col_op(mat: &mut Array2, ctrl: usize, trgt: usize) { // If full_elim = True, it allows full elimination of mat[:, 0 : ncols] // Returns the matrix mat. fn gauss_elimination( - mut mat: Array2, + mut mat: Array2, ncols: Option, full_elim: Option, -) -> Array2 { +) -> Array2 { let (m, mut n) = (mat.nrows(), mat.ncols()); // no. of rows and columns if let Some(ncols_val) = ncols { n = usize::min(n, ncols_val); // no. of active columns @@ -50,7 +50,7 @@ fn gauss_elimination( for j in k..n { new_k = k; for i in r..m { - if mat[(i, j)] { + if mat[(i, j)] == 1 { is_non_zero = true; new_k = j; new_r = i; @@ -72,16 +72,16 @@ fn gauss_elimination( mat.slice_mut(s![new_r, ..]).assign(&temp_r); } - if full_elim.is_some() { + if full_elim.is_some() && full_elim == Some(true) { for i in 0..r { - if mat[(i, new_k)] { + if mat[(i, new_k)] == 1 { _row_op(&mut mat, r, i); } } } for i in r + 1..m { - if mat[(i, new_k)] { + if mat[(i, new_k)] == 1 { _row_op(&mut mat, r, i); } } @@ -94,10 +94,10 @@ fn gauss_elimination( #[pyo3(signature = (mat, ncols=None, full_elim=false))] fn _gauss_elimination( py: Python, - mat: PyArrayLike2, + mat: PyArrayLike2, ncols: Option, full_elim: Option, -) -> PyResult>> { +) -> PyResult>> { let view = mat.as_array().to_owned(); Ok(gauss_elimination(view, ncols, full_elim) .into_pyarray_bound(py) From 678ec0aed9b06c0d082875ee28451a8745dc2e0d Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Tue, 4 Jun 2024 08:51:10 -0500 Subject: [PATCH 05/43] add parallelization in row operations --- crates/accelerate/src/linear_matrix.rs | 32 ++++++++++++++++---------- 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index d708e71bd108..2305d49e56e8 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -10,9 +10,10 @@ // copyright notice, and modified files need to carry a notice indicating // that they have been altered from the originals. -use ndarray::{s, Array2}; +use ndarray::{s, Array2, Axis}; use numpy::{AllowTypeChange, IntoPyArray, PyArray2, PyArrayLike2}; use pyo3::prelude::*; +use rayon::prelude::*; // Perform ROW operation on a matrix mat fn _row_op(mat: &mut Array2, ctrl: usize, trgt: usize) { @@ -72,19 +73,26 @@ fn gauss_elimination( mat.slice_mut(s![new_r, ..]).assign(&temp_r); } - if full_elim.is_some() && full_elim == Some(true) { - for i in 0..r { - if mat[(i, new_k)] == 1 { - _row_op(&mut mat, r, i); - } - } + // Copy source row to avoid trying multiple borrows at once + let row0 = mat.row(r).to_owned(); + if full_elim == Some(true) { + mat.axis_iter_mut(Axis(0)) + .into_par_iter() + .enumerate() + .filter(|(i, row)| (*i < r) && row[new_k] == 1) + .for_each(|(_i, mut row)| { + row.zip_mut_with(&row0, |x, &y| *x ^= y); + }); } - for i in r + 1..m { - if mat[(i, new_k)] == 1 { - _row_op(&mut mat, r, i); - } - } + mat.axis_iter_mut(Axis(0)) + .into_par_iter() + .enumerate() + .filter(|(i, row)| (*i > r) && (*i < m) && row[new_k] == 1) + .for_each(|(_i, mut row)| { + row.zip_mut_with(&row0, |x, &y| *x ^= y); + }); + r += 1; } mat From 361b562934c82ea272297e9346d06c2f0d95110f Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Wed, 5 Jun 2024 03:12:12 -0500 Subject: [PATCH 06/43] update matrices in place --- crates/accelerate/src/linear_matrix.rs | 28 +++++++++----------------- 1 file changed, 10 insertions(+), 18 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 2305d49e56e8..6546196abc5f 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -10,20 +10,20 @@ // copyright notice, and modified files need to carry a notice indicating // that they have been altered from the originals. -use ndarray::{s, Array2, Axis}; -use numpy::{AllowTypeChange, IntoPyArray, PyArray2, PyArrayLike2}; +use ndarray::{s, ArrayViewMut2, Axis}; +use numpy::PyReadwriteArray2; use pyo3::prelude::*; use rayon::prelude::*; // Perform ROW operation on a matrix mat -fn _row_op(mat: &mut Array2, ctrl: usize, trgt: usize) { +fn _row_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { let row0 = mat.row(ctrl).to_owned(); let mut row1 = mat.row_mut(trgt); row1.zip_mut_with(&row0, |x, &y| *x ^= y); } // Perform COL operation on a matrix mat -fn _col_op(mat: &mut Array2, ctrl: usize, trgt: usize) { +fn _col_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { let col0 = mat.column(ctrl).to_owned(); let mut col1 = mat.column_mut(trgt); col1.zip_mut_with(&col0, |x, &y| *x ^= y); @@ -32,11 +32,7 @@ fn _col_op(mat: &mut Array2, ctrl: usize, trgt: usize) { // Gauss elimination of a matrix mat with m rows and n columns. // If full_elim = True, it allows full elimination of mat[:, 0 : ncols] // Returns the matrix mat. -fn gauss_elimination( - mut mat: Array2, - ncols: Option, - full_elim: Option, -) -> Array2 { +fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_elim: Option) { let (m, mut n) = (mat.nrows(), mat.ncols()); // no. of rows and columns if let Some(ncols_val) = ncols { n = usize::min(n, ncols_val); // no. of active columns @@ -63,7 +59,7 @@ fn gauss_elimination( } } if !is_non_zero { - return mat; // A is in the canonical form + return; // A is in the canonical form } if new_r != r { @@ -95,21 +91,17 @@ fn gauss_elimination( r += 1; } - mat } #[pyfunction] #[pyo3(signature = (mat, ncols=None, full_elim=false))] fn _gauss_elimination( - py: Python, - mat: PyArrayLike2, + mut mat: PyReadwriteArray2, ncols: Option, full_elim: Option, -) -> PyResult>> { - let view = mat.as_array().to_owned(); - Ok(gauss_elimination(view, ncols, full_elim) - .into_pyarray_bound(py) - .unbind()) +) { + let view = mat.as_array_mut(); + gauss_elimination(view, ncols, full_elim); } #[pymodule] From 19326bd27bedb9f3c435e033e2bc55d8ec1fcd9b Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Wed, 5 Jun 2024 04:15:53 -0500 Subject: [PATCH 07/43] change matrix type in rust code to bool --- crates/accelerate/src/linear_matrix.rs | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 6546196abc5f..d1e045024b64 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -16,14 +16,14 @@ use pyo3::prelude::*; use rayon::prelude::*; // Perform ROW operation on a matrix mat -fn _row_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { +fn _row_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { let row0 = mat.row(ctrl).to_owned(); let mut row1 = mat.row_mut(trgt); row1.zip_mut_with(&row0, |x, &y| *x ^= y); } // Perform COL operation on a matrix mat -fn _col_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { +fn _col_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { let col0 = mat.column(ctrl).to_owned(); let mut col1 = mat.column_mut(trgt); col1.zip_mut_with(&col0, |x, &y| *x ^= y); @@ -32,7 +32,7 @@ fn _col_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { // Gauss elimination of a matrix mat with m rows and n columns. // If full_elim = True, it allows full elimination of mat[:, 0 : ncols] // Returns the matrix mat. -fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_elim: Option) { +fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_elim: Option) { let (m, mut n) = (mat.nrows(), mat.ncols()); // no. of rows and columns if let Some(ncols_val) = ncols { n = usize::min(n, ncols_val); // no. of active columns @@ -47,7 +47,7 @@ fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_elim for j in k..n { new_k = k; for i in r..m { - if mat[(i, j)] == 1 { + if mat[(i, j)] { is_non_zero = true; new_k = j; new_r = i; @@ -75,7 +75,7 @@ fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_elim mat.axis_iter_mut(Axis(0)) .into_par_iter() .enumerate() - .filter(|(i, row)| (*i < r) && row[new_k] == 1) + .filter(|(i, row)| (*i < r) && row[new_k]) .for_each(|(_i, mut row)| { row.zip_mut_with(&row0, |x, &y| *x ^= y); }); @@ -84,7 +84,7 @@ fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_elim mat.axis_iter_mut(Axis(0)) .into_par_iter() .enumerate() - .filter(|(i, row)| (*i > r) && (*i < m) && row[new_k] == 1) + .filter(|(i, row)| (*i > r) && (*i < m) && row[new_k]) .for_each(|(_i, mut row)| { row.zip_mut_with(&row0, |x, &y| *x ^= y); }); @@ -96,7 +96,7 @@ fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_elim #[pyfunction] #[pyo3(signature = (mat, ncols=None, full_elim=false))] fn _gauss_elimination( - mut mat: PyReadwriteArray2, + mut mat: PyReadwriteArray2, ncols: Option, full_elim: Option, ) { From 0cf7502425b66426877fe845d6c6c6a90bacf531 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Wed, 5 Jun 2024 04:18:07 -0500 Subject: [PATCH 08/43] handle type in python code --- .../clifford/clifford_decompose_layers.py | 3 ++- qiskit/synthesis/linear/linear_matrix_utils.py | 14 ++++++++------ 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/qiskit/synthesis/clifford/clifford_decompose_layers.py b/qiskit/synthesis/clifford/clifford_decompose_layers.py index 2fc9ca5bdb29..9bb6bf8caa9d 100644 --- a/qiskit/synthesis/clifford/clifford_decompose_layers.py +++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py @@ -209,7 +209,8 @@ def _create_graph_state(cliff, validate=False): if rank < num_qubits: stab = cliff.stab[:, :-1] - stab = _gauss_elimination(stab, num_qubits) + stab = stab.astype(bool) + _gauss_elimination(stab, num_qubits) Cmat = stab[rank:num_qubits, num_qubits:] Cmat = np.transpose(Cmat) diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index afce1f3bf510..41d0da6a780d 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -118,7 +118,7 @@ def calc_inverse_matrix(mat: np.ndarray, verify: bool = False): verify: if True asserts that the multiplication of mat and its inverse is the identity matrix. Returns: - np.ndarray: the inverse matrix. + np.ndarray: the inverse boolean matrix. Raises: QiskitError: if the matrix is not square. @@ -130,8 +130,9 @@ def calc_inverse_matrix(mat: np.ndarray, verify: bool = False): n = mat.shape[0] # concatenate the matrix and identity - mat1 = np.concatenate((mat, np.eye(n, dtype=int)), axis=1) - mat1 = _gauss_elimination(mat1, None, full_elim=True) + mat1 = np.concatenate((mat.astype(bool), np.eye(n, dtype=bool)), axis=1) + _gauss_elimination(mat1, None, full_elim=True) + mat1 = mat1.astype(int) r = _compute_rank_after_gauss_elim(mat1[:, 0:n]) @@ -148,14 +149,15 @@ def calc_inverse_matrix(mat: np.ndarray, verify: bool = False): def _compute_rank_after_gauss_elim(mat): - """Given a matrix A after Gaussian elimination, computes its rank + """Given a boolean matrix A after Gaussian elimination, computes its rank (i.e. simply the number of nonzero rows)""" return np.sum(mat.any(axis=1)) def _compute_rank(mat): - """Given a matrix A computes its rank""" - mat = _gauss_elimination(mat) + """Given a boolean matrix A computes its rank""" + mat = mat.astype(bool) + _gauss_elimination(mat) return np.sum(mat.any(axis=1)) From c8aee3aa4945891c2c0b464838f469a8c6a2958f Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Wed, 5 Jun 2024 09:54:15 -0500 Subject: [PATCH 09/43] update filter following review --- crates/accelerate/src/linear_matrix.rs | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index d1e045024b64..0df19bbd598d 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -71,20 +71,13 @@ fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_el // Copy source row to avoid trying multiple borrows at once let row0 = mat.row(r).to_owned(); - if full_elim == Some(true) { - mat.axis_iter_mut(Axis(0)) - .into_par_iter() - .enumerate() - .filter(|(i, row)| (*i < r) && row[new_k]) - .for_each(|(_i, mut row)| { - row.zip_mut_with(&row0, |x, &y| *x ^= y); - }); - } - mat.axis_iter_mut(Axis(0)) .into_par_iter() .enumerate() - .filter(|(i, row)| (*i > r) && (*i < m) && row[new_k]) + .filter(|(i, row)| { + (full_elim == Some(true) && (*i < r) && row[new_k]) + || (*i > r && *i < m && row[new_k]) + }) .for_each(|(_i, mut row)| { row.zip_mut_with(&row0, |x, &y| *x ^= y); }); From 2f45606bb984e6ff3dde083457e5525a7c2e5a87 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 6 Jun 2024 07:02:09 -0500 Subject: [PATCH 10/43] remove parallelization using rayon --- crates/accelerate/src/linear_matrix.rs | 2 -- 1 file changed, 2 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 0df19bbd598d..6d526afc977c 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -13,7 +13,6 @@ use ndarray::{s, ArrayViewMut2, Axis}; use numpy::PyReadwriteArray2; use pyo3::prelude::*; -use rayon::prelude::*; // Perform ROW operation on a matrix mat fn _row_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { @@ -72,7 +71,6 @@ fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_el // Copy source row to avoid trying multiple borrows at once let row0 = mat.row(r).to_owned(); mat.axis_iter_mut(Axis(0)) - .into_par_iter() .enumerate() .filter(|(i, row)| { (full_elim == Some(true) && (*i < r) && row[new_k]) From 3bcda6d1a2cb7fd97bb3f7f3b7ab96535c85cbe7 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 13 Jun 2024 05:12:28 -0500 Subject: [PATCH 11/43] move _gauss_elimination_with_perm to rust --- crates/accelerate/src/linear_matrix.rs | 35 +++++++++++-- .../clifford/clifford_decompose_layers.py | 2 +- .../synthesis/linear/linear_matrix_utils.py | 52 +------------------ 3 files changed, 34 insertions(+), 55 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 6d526afc977c..f8d5314834f1 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -30,12 +30,18 @@ fn _col_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { // Gauss elimination of a matrix mat with m rows and n columns. // If full_elim = True, it allows full elimination of mat[:, 0 : ncols] -// Returns the matrix mat. -fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_elim: Option) { +// Returns the matrix mat, and the permutation perm that was done on the rows during the process. +// perm[0 : rank] represents the indices of linearly independent rows in the original matrix. +fn gauss_elimination_with_perm( + mut mat: ArrayViewMut2, + ncols: Option, + full_elim: Option, +) -> Vec { let (m, mut n) = (mat.nrows(), mat.ncols()); // no. of rows and columns if let Some(ncols_val) = ncols { n = usize::min(n, ncols_val); // no. of active columns } + let mut perm: Vec = Vec::from_iter(0..m); let mut r = 0; // current rank let k = 0; // current pivot column @@ -58,7 +64,7 @@ fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_el } } if !is_non_zero { - return; // A is in the canonical form + return perm; // A is in the canonical form } if new_r != r { @@ -66,6 +72,7 @@ fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_el let temp_new_r = mat.slice_mut(s![new_r, ..]).to_owned(); mat.slice_mut(s![r, ..]).assign(&temp_new_r); mat.slice_mut(s![new_r, ..]).assign(&temp_r); + perm.swap(r, new_r); } // Copy source row to avoid trying multiple borrows at once @@ -82,6 +89,27 @@ fn gauss_elimination(mut mat: ArrayViewMut2, ncols: Option, full_el r += 1; } + perm +} + +// Gauss elimination of a matrix mat with m rows and n columns. +// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] +// Returns the updated matrix mat. +fn gauss_elimination(mat: ArrayViewMut2, ncols: Option, full_elim: Option) { + let _perm = gauss_elimination_with_perm(mat, ncols, full_elim); +} + +#[pyfunction] +#[pyo3(signature = (mat, ncols=None, full_elim=false))] +fn _gauss_elimination_with_perm( + py: Python, + mut mat: PyReadwriteArray2, + ncols: Option, + full_elim: Option, +) -> PyResult { + let view = mat.as_array_mut(); + let perm = gauss_elimination_with_perm(view, ncols, full_elim); + Ok(perm.to_object(py)) } #[pyfunction] @@ -97,6 +125,7 @@ fn _gauss_elimination( #[pymodule] pub fn linear_matrix(m: &Bound) -> PyResult<()> { + m.add_wrapped(wrap_pyfunction!(_gauss_elimination_with_perm))?; m.add_wrapped(wrap_pyfunction!(_gauss_elimination))?; Ok(()) } diff --git a/qiskit/synthesis/clifford/clifford_decompose_layers.py b/qiskit/synthesis/clifford/clifford_decompose_layers.py index 9bb6bf8caa9d..4c66f8e6bb0c 100644 --- a/qiskit/synthesis/clifford/clifford_decompose_layers.py +++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py @@ -214,7 +214,7 @@ def _create_graph_state(cliff, validate=False): Cmat = stab[rank:num_qubits, num_qubits:] Cmat = np.transpose(Cmat) - Cmat, perm = _gauss_elimination_with_perm(Cmat) + perm = _gauss_elimination_with_perm(Cmat) perm = perm[0 : num_qubits - rank] # validate that the output matrix has the same rank diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index 41d0da6a780d..45c9ba591b91 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -17,7 +17,7 @@ from qiskit.exceptions import QiskitError # pylint: disable=unused-import -from qiskit._accelerate.linear_matrix import _gauss_elimination +from qiskit._accelerate.linear_matrix import _gauss_elimination, _gauss_elimination_with_perm def check_invertible_binary_matrix(mat: np.ndarray): @@ -60,56 +60,6 @@ def random_invertible_binary_matrix( return mat -def _gauss_elimination_with_perm(mat, ncols=None, full_elim=False): - """Gauss elimination of a matrix mat with m rows and n columns. - If full_elim = True, it allows full elimination of mat[:, 0 : ncols] - Returns the matrix mat, and the permutation perm that was done on the rows during the process. - perm[0 : rank] represents the indices of linearly independent rows in the original matrix.""" - - # Treat the matrix A as containing integer values - mat = np.array(mat, dtype=int, copy=True) - - m = mat.shape[0] # no. of rows - n = mat.shape[1] # no. of columns - if ncols is not None: - n = min(n, ncols) # no. of active columns - - perm = np.array(range(m)) # permutation on the rows - - r = 0 # current rank - k = 0 # current pivot column - while (r < m) and (k < n): - is_non_zero = False - new_r = r - for j in range(k, n): - for i in range(r, m): - if mat[i][j]: - is_non_zero = True - k = j - new_r = i - break - if is_non_zero: - break - if not is_non_zero: - return mat, perm # A is in the canonical form - - if new_r != r: - mat[[r, new_r]] = mat[[new_r, r]] - perm[r], perm[new_r] = perm[new_r], perm[r] - - if full_elim: - for i in range(0, r): - if mat[i][k]: - mat[i] = mat[i] ^ mat[r] - - for i in range(r + 1, m): - if mat[i][k]: - mat[i] = mat[i] ^ mat[r] - r += 1 - - return mat, perm - - def calc_inverse_matrix(mat: np.ndarray, verify: bool = False): """Given a square numpy(dtype=int) matrix mat, tries to compute its inverse. From b4f05bfdeb22f00fc99125d06944131b099e77bb Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 13 Jun 2024 07:11:01 -0500 Subject: [PATCH 12/43] fix fmt error --- crates/pyext/src/lib.rs | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/crates/pyext/src/lib.rs b/crates/pyext/src/lib.rs index 15f1f4a5f0a0..e8721b409ce2 100644 --- a/crates/pyext/src/lib.rs +++ b/crates/pyext/src/lib.rs @@ -16,11 +16,12 @@ use pyo3::wrap_pymodule; use qiskit_accelerate::{ convert_2q_block_matrix::convert_2q_block_matrix, dense_layout::dense_layout, error_map::error_map, euler_one_qubit_decomposer::euler_one_qubit_decomposer, - isometry::isometry, linear_matrix::linear_matrix, nlayout::nlayout, optimize_1q_gates::optimize_1q_gates, - pauli_exp_val::pauli_expval, permutation::permutation, results::results, sabre::sabre, - sampled_exp_val::sampled_exp_val, sparse_pauli_op::sparse_pauli_op, - stochastic_swap::stochastic_swap, two_qubit_decompose::two_qubit_decompose, uc_gate::uc_gate, - utils::utils, vf2_layout::vf2_layout, + isometry::isometry, linear_matrix::linear_matrix, nlayout::nlayout, + optimize_1q_gates::optimize_1q_gates, pauli_exp_val::pauli_expval, permutation::permutation, + results::results, sabre::sabre, sampled_exp_val::sampled_exp_val, + sparse_pauli_op::sparse_pauli_op, stochastic_swap::stochastic_swap, + two_qubit_decompose::two_qubit_decompose, uc_gate::uc_gate, utils::utils, + vf2_layout::vf2_layout, }; #[pymodule] From b79ea04e45006ba24233fd0433676da1b31f0148 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 13 Jun 2024 08:51:28 -0500 Subject: [PATCH 13/43] simplify _gauss_elimination function --- crates/accelerate/src/linear_matrix.rs | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index f8d5314834f1..77dd4db95aa8 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -92,13 +92,6 @@ fn gauss_elimination_with_perm( perm } -// Gauss elimination of a matrix mat with m rows and n columns. -// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] -// Returns the updated matrix mat. -fn gauss_elimination(mat: ArrayViewMut2, ncols: Option, full_elim: Option) { - let _perm = gauss_elimination_with_perm(mat, ncols, full_elim); -} - #[pyfunction] #[pyo3(signature = (mat, ncols=None, full_elim=false))] fn _gauss_elimination_with_perm( @@ -114,13 +107,16 @@ fn _gauss_elimination_with_perm( #[pyfunction] #[pyo3(signature = (mat, ncols=None, full_elim=false))] +// Gauss elimination of a matrix mat with m rows and n columns. +// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] +// Returns the updated matrix mat. fn _gauss_elimination( mut mat: PyReadwriteArray2, ncols: Option, full_elim: Option, ) { let view = mat.as_array_mut(); - gauss_elimination(view, ncols, full_elim); + let _perm = gauss_elimination_with_perm(view, ncols, full_elim); } #[pymodule] From 21f060e558c43aec6d9080c91edfdaf1b1493edd Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 13 Jun 2024 10:01:51 -0500 Subject: [PATCH 14/43] update _compute_rank_after_gauss_elim to rust --- crates/accelerate/src/linear_matrix.rs | 16 +++++++++++++++- qiskit/synthesis/linear/linear_matrix_utils.py | 14 ++++++-------- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 77dd4db95aa8..83f15aa2bd72 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -11,7 +11,7 @@ // that they have been altered from the originals. use ndarray::{s, ArrayViewMut2, Axis}; -use numpy::PyReadwriteArray2; +use numpy::{PyReadonlyArray2, PyReadwriteArray2}; use pyo3::prelude::*; // Perform ROW operation on a matrix mat @@ -119,9 +119,23 @@ fn _gauss_elimination( let _perm = gauss_elimination_with_perm(view, ncols, full_elim); } +#[pyfunction] +#[pyo3(signature = (mat))] +// Given a boolean matrix A after Gaussian elimination, computes its rank +// (i.e. simply the number of nonzero rows)""" +fn _compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyResult { + let view = mat.as_array(); + let rank: usize = view + .axis_iter(Axis(0)) + .map(|row| row.fold(false, |out, val| out | *val) as usize) + .sum(); + Ok(rank.to_object(py)) +} + #[pymodule] pub fn linear_matrix(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(_gauss_elimination_with_perm))?; m.add_wrapped(wrap_pyfunction!(_gauss_elimination))?; + m.add_wrapped(wrap_pyfunction!(_compute_rank_after_gauss_elim))?; Ok(()) } diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index 45c9ba591b91..8d17bda99421 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -17,7 +17,11 @@ from qiskit.exceptions import QiskitError # pylint: disable=unused-import -from qiskit._accelerate.linear_matrix import _gauss_elimination, _gauss_elimination_with_perm +from qiskit._accelerate.linear_matrix import ( + _gauss_elimination, + _gauss_elimination_with_perm, + _compute_rank_after_gauss_elim, +) def check_invertible_binary_matrix(mat: np.ndarray): @@ -82,13 +86,13 @@ def calc_inverse_matrix(mat: np.ndarray, verify: bool = False): # concatenate the matrix and identity mat1 = np.concatenate((mat.astype(bool), np.eye(n, dtype=bool)), axis=1) _gauss_elimination(mat1, None, full_elim=True) - mat1 = mat1.astype(int) r = _compute_rank_after_gauss_elim(mat1[:, 0:n]) if r < n: raise QiskitError("The matrix is not invertible.") + mat1 = mat1.astype(int) matinv = mat1[:, n : 2 * n] if verify: @@ -98,12 +102,6 @@ def calc_inverse_matrix(mat: np.ndarray, verify: bool = False): return matinv -def _compute_rank_after_gauss_elim(mat): - """Given a boolean matrix A after Gaussian elimination, computes its rank - (i.e. simply the number of nonzero rows)""" - return np.sum(mat.any(axis=1)) - - def _compute_rank(mat): """Given a boolean matrix A computes its rank""" mat = mat.astype(bool) From 598f91388d6164de4cd6a3472dfc9ed2c3861930 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Mon, 17 Jun 2024 02:45:30 -0500 Subject: [PATCH 15/43] update _row_op and _col_op --- crates/accelerate/src/linear_matrix.rs | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 83f15aa2bd72..c7c71f0791ca 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -16,15 +16,13 @@ use pyo3::prelude::*; // Perform ROW operation on a matrix mat fn _row_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { - let row0 = mat.row(ctrl).to_owned(); - let mut row1 = mat.row_mut(trgt); + let (row0, mut row1) = mat.multi_slice_mut((s![ctrl, ..], s![trgt, ..])); row1.zip_mut_with(&row0, |x, &y| *x ^= y); } // Perform COL operation on a matrix mat fn _col_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { - let col0 = mat.column(ctrl).to_owned(); - let mut col1 = mat.column_mut(trgt); + let (col0, mut col1) = mat.multi_slice_mut((s![.., ctrl], s![.., trgt])); col1.zip_mut_with(&col0, |x, &y| *x ^= y); } From 7510b2c93176c9ea0c6050cba2c906cf19b48714 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Tue, 18 Jun 2024 03:29:19 -0500 Subject: [PATCH 16/43] transfer _row_op and _col_op from python to rust --- crates/accelerate/src/linear_matrix.rs | 40 +++++++++++-------- qiskit/synthesis/linear/linear_depth_lnn.py | 2 +- .../synthesis/linear/linear_matrix_utils.py | 12 +----- 3 files changed, 27 insertions(+), 27 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index c7c71f0791ca..1451a8405be1 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -14,18 +14,6 @@ use ndarray::{s, ArrayViewMut2, Axis}; use numpy::{PyReadonlyArray2, PyReadwriteArray2}; use pyo3::prelude::*; -// Perform ROW operation on a matrix mat -fn _row_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { - let (row0, mut row1) = mat.multi_slice_mut((s![ctrl, ..], s![trgt, ..])); - row1.zip_mut_with(&row0, |x, &y| *x ^= y); -} - -// Perform COL operation on a matrix mat -fn _col_op(mat: &mut ArrayViewMut2, ctrl: usize, trgt: usize) { - let (col0, mut col1) = mat.multi_slice_mut((s![.., ctrl], s![.., trgt])); - col1.zip_mut_with(&col0, |x, &y| *x ^= y); -} - // Gauss elimination of a matrix mat with m rows and n columns. // If full_elim = True, it allows full elimination of mat[:, 0 : ncols] // Returns the matrix mat, and the permutation perm that was done on the rows during the process. @@ -98,8 +86,8 @@ fn _gauss_elimination_with_perm( ncols: Option, full_elim: Option, ) -> PyResult { - let view = mat.as_array_mut(); - let perm = gauss_elimination_with_perm(view, ncols, full_elim); + let matmut = mat.as_array_mut(); + let perm = gauss_elimination_with_perm(matmut, ncols, full_elim); Ok(perm.to_object(py)) } @@ -113,8 +101,8 @@ fn _gauss_elimination( ncols: Option, full_elim: Option, ) { - let view = mat.as_array_mut(); - let _perm = gauss_elimination_with_perm(view, ncols, full_elim); + let matmut = mat.as_array_mut(); + let _perm = gauss_elimination_with_perm(matmut, ncols, full_elim); } #[pyfunction] @@ -130,10 +118,30 @@ fn _compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> Py Ok(rank.to_object(py)) } +#[pyfunction] +#[pyo3(signature = (mat, ctrl, trgt))] +// Perform ROW operation on a matrix mat +fn _row_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { + let mut matmut = mat.as_array_mut(); + let (row0, mut row1) = matmut.multi_slice_mut((s![ctrl, ..], s![trgt, ..])); + row1.zip_mut_with(&row0, |x, &y| *x ^= y); +} + +#[pyfunction] +#[pyo3(signature = (mat, ctrl, trgt))] +// Perform COL operation on a matrix mat +fn _col_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { + let mut matmut = mat.as_array_mut(); + let (col0, mut col1) = matmut.multi_slice_mut((s![.., trgt], s![.., ctrl])); + col1.zip_mut_with(&col0, |x, &y| *x ^= y); +} + #[pymodule] pub fn linear_matrix(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(_gauss_elimination_with_perm))?; m.add_wrapped(wrap_pyfunction!(_gauss_elimination))?; m.add_wrapped(wrap_pyfunction!(_compute_rank_after_gauss_elim))?; + m.add_wrapped(wrap_pyfunction!(_row_op))?; + m.add_wrapped(wrap_pyfunction!(_col_op))?; Ok(()) } diff --git a/qiskit/synthesis/linear/linear_depth_lnn.py b/qiskit/synthesis/linear/linear_depth_lnn.py index 2d544f37ef94..9b1590291f7b 100644 --- a/qiskit/synthesis/linear/linear_depth_lnn.py +++ b/qiskit/synthesis/linear/linear_depth_lnn.py @@ -223,7 +223,7 @@ def _optimize_cx_circ_depth_5n_line(mat): # According to [1] the synthesis is done on the inverse matrix # so the matrix mat is inverted at this step mat_inv = mat.copy() - mat_cpy = calc_inverse_matrix(mat_inv) + mat_cpy = calc_inverse_matrix(mat_inv).astype(bool) n = len(mat_cpy) diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index 8d17bda99421..df4cb18a8a3e 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -21,6 +21,8 @@ _gauss_elimination, _gauss_elimination_with_perm, _compute_rank_after_gauss_elim, + _row_op, + _col_op, ) @@ -107,13 +109,3 @@ def _compute_rank(mat): mat = mat.astype(bool) _gauss_elimination(mat) return np.sum(mat.any(axis=1)) - - -def _row_op(mat, ctrl, trgt): - # Perform ROW operation on a matrix mat - mat[trgt] = mat[trgt] ^ mat[ctrl] - - -def _col_op(mat, ctrl, trgt): - # Perform COL operation on a matrix mat - mat[:, ctrl] = mat[:, trgt] ^ mat[:, ctrl] From 44296aefcbf809d990d5dc9c5959433200372aec Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Tue, 18 Jun 2024 05:01:36 -0500 Subject: [PATCH 17/43] fix code due to failing tests --- qiskit/synthesis/linear/linear_depth_lnn.py | 1 + qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py | 1 + 2 files changed, 2 insertions(+) diff --git a/qiskit/synthesis/linear/linear_depth_lnn.py b/qiskit/synthesis/linear/linear_depth_lnn.py index 9b1590291f7b..f12da92ad8b8 100644 --- a/qiskit/synthesis/linear/linear_depth_lnn.py +++ b/qiskit/synthesis/linear/linear_depth_lnn.py @@ -45,6 +45,7 @@ def _get_lower_triangular(n, mat, mat_inv): mat = mat.copy() mat_t = mat.copy() mat_inv_t = mat_inv.copy() + mat_inv_t = mat_inv_t.astype(bool) cx_instructions_rows = [] diff --git a/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py index 23f24e07eab3..119a65d6fe23 100644 --- a/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py +++ b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py @@ -245,6 +245,7 @@ def synth_cx_cz_depth_line_my(mat_x: np.ndarray, mat_z: np.ndarray) -> QuantumCi n = len(mat_x) mat_x = calc_inverse_matrix(mat_x) + mat_x = mat_x.astype(bool) cx_instructions_rows_m2nw, cx_instructions_rows_nw2id = _optimize_cx_circ_depth_5n_line(mat_x) # Meanwhile, also build the -CZ- circuit via Phase gate insertions as per Algorithm 2 [2] From 45fef0f5272dda48d9fc967f9d7013114cd7b315 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Tue, 18 Jun 2024 05:13:36 -0500 Subject: [PATCH 18/43] minor update of types --- qiskit/synthesis/linear/linear_depth_lnn.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/qiskit/synthesis/linear/linear_depth_lnn.py b/qiskit/synthesis/linear/linear_depth_lnn.py index f12da92ad8b8..21e8f3717b69 100644 --- a/qiskit/synthesis/linear/linear_depth_lnn.py +++ b/qiskit/synthesis/linear/linear_depth_lnn.py @@ -45,7 +45,6 @@ def _get_lower_triangular(n, mat, mat_inv): mat = mat.copy() mat_t = mat.copy() mat_inv_t = mat_inv.copy() - mat_inv_t = mat_inv_t.astype(bool) cx_instructions_rows = [] @@ -223,7 +222,7 @@ def _optimize_cx_circ_depth_5n_line(mat): # According to [1] the synthesis is done on the inverse matrix # so the matrix mat is inverted at this step - mat_inv = mat.copy() + mat_inv = mat.copy().astype(bool) mat_cpy = calc_inverse_matrix(mat_inv).astype(bool) n = len(mat_cpy) From 89479a5c4ba436dec971596e3ef445e2ee926c54 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Wed, 19 Jun 2024 06:59:13 -0500 Subject: [PATCH 19/43] move calc_inverse_matrix to rust, add _binary_matmul in rust --- crates/accelerate/src/linear_matrix.rs | 106 ++++++++++++++++-- .../clifford/clifford_decompose_layers.py | 5 +- .../synthesis/linear/linear_matrix_utils.py | 40 +------ .../python/synthesis/test_linear_synthesis.py | 4 +- 4 files changed, 107 insertions(+), 48 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 1451a8405be1..5fd5fe9cc4d2 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -10,10 +10,35 @@ // copyright notice, and modified files need to carry a notice indicating // that they have been altered from the originals. -use ndarray::{s, ArrayViewMut2, Axis}; -use numpy::{PyReadonlyArray2, PyReadwriteArray2}; +use ndarray::{concatenate, s, Array2, ArrayView2, ArrayViewMut2, Axis}; +use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2, PyReadwriteArray2}; use pyo3::prelude::*; +// Binary matrix multiplication +fn binary_matmul(mat1: ArrayView2, mat2: ArrayView2) -> Array2 { + let n1_rows = mat1.nrows(); + let n1_cols = mat1.ncols(); + let n2_rows = mat2.nrows(); + let n2_cols = mat2.ncols(); + if n1_cols != n2_rows { + panic!( + "Cannot multiply matrices with inappropriate dimensions {}, {}", + n1_cols, n2_rows + ); + } + + let mut result: Array2 = Array2::from_shape_fn((n1_rows, n2_cols), |(_i, _j)| false); + for i in 0..n1_rows { + for j in 0..n2_cols { + for k in 0..n2_rows { + result[(i, j)] ^= mat1[(i, k)] & mat2[(k, j)]; + } + } + } + + result +} + // Gauss elimination of a matrix mat with m rows and n columns. // If full_elim = True, it allows full elimination of mat[:, 0 : ncols] // Returns the matrix mat, and the permutation perm that was done on the rows during the process. @@ -78,6 +103,47 @@ fn gauss_elimination_with_perm( perm } +// Given a boolean matrix A after Gaussian elimination, computes its rank +// (i.e. simply the number of nonzero rows) +fn compute_rank_after_gauss_elim(mat: ArrayView2) -> usize { + let rank: usize = mat + .axis_iter(Axis(0)) + .map(|row| row.fold(false, |out, val| out | *val) as usize) + .sum(); + rank +} + +// Given a square boolean matrix mat, tries to compute its inverse. +fn calc_inverse_matrix_inner(mat: ArrayView2, verify: bool) -> Array2 { + if mat.shape()[0] != mat.shape()[1] { + panic!("Matrix to invert is a non-square matrix."); + } + let n = mat.shape()[0]; + + // concatenate the matrix and identity + let identity_matrix: Array2 = Array2::from_shape_fn((n, n), |(i, j)| i == j); + let mut mat1 = concatenate(Axis(1), &[mat.view(), identity_matrix.view()]).unwrap(); + + gauss_elimination_with_perm(mat1.view_mut(), None, Some(true)); + + let r = compute_rank_after_gauss_elim(mat1.slice(s![.., 0..n])); + if r < n { + panic!("The matrix is not invertible."); + } + + let invmat = mat1.slice(s![.., n..2 * n]).to_owned(); + + if verify { + let mat2 = binary_matmul(mat, (&invmat).into()); + let identity_matrix: Array2 = Array2::from_shape_fn((n, n), |(i, j)| i == j); + if mat2.ne(&identity_matrix) { + panic!("The inverse matrix is not correct."); + } + } + + invmat +} + #[pyfunction] #[pyo3(signature = (mat, ncols=None, full_elim=false))] fn _gauss_elimination_with_perm( @@ -108,16 +174,40 @@ fn _gauss_elimination( #[pyfunction] #[pyo3(signature = (mat))] // Given a boolean matrix A after Gaussian elimination, computes its rank -// (i.e. simply the number of nonzero rows)""" +// (i.e. simply the number of nonzero rows) fn _compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyResult { let view = mat.as_array(); - let rank: usize = view - .axis_iter(Axis(0)) - .map(|row| row.fold(false, |out, val| out | *val) as usize) - .sum(); + let rank = compute_rank_after_gauss_elim(view); Ok(rank.to_object(py)) } +#[pyfunction] +#[pyo3(signature = (mat, verify=false))] +// Calculates the inverse matrix +pub fn calc_inverse_matrix( + py: Python, + mat: PyReadonlyArray2, + verify: Option, +) -> PyResult>> { + let view = mat.as_array(); + let invmat = calc_inverse_matrix_inner(view, verify.is_some()); + Ok(invmat.into_pyarray_bound(py).unbind()) +} + +#[pyfunction] +#[pyo3(signature = (mat1, mat2))] +// Binary matrix multiplication +pub fn _binary_matmul( + py: Python, + mat1: PyReadonlyArray2, + mat2: PyReadonlyArray2, +) -> PyResult>> { + let view1 = mat1.as_array(); + let view2 = mat2.as_array(); + let result = binary_matmul(view1, view2); + Ok(result.into_pyarray_bound(py).unbind()) +} + #[pyfunction] #[pyo3(signature = (mat, ctrl, trgt))] // Perform ROW operation on a matrix mat @@ -141,7 +231,9 @@ pub fn linear_matrix(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(_gauss_elimination_with_perm))?; m.add_wrapped(wrap_pyfunction!(_gauss_elimination))?; m.add_wrapped(wrap_pyfunction!(_compute_rank_after_gauss_elim))?; + m.add_wrapped(wrap_pyfunction!(calc_inverse_matrix))?; m.add_wrapped(wrap_pyfunction!(_row_op))?; m.add_wrapped(wrap_pyfunction!(_col_op))?; + m.add_wrapped(wrap_pyfunction!(_binary_matmul))?; Ok(()) } diff --git a/qiskit/synthesis/clifford/clifford_decompose_layers.py b/qiskit/synthesis/clifford/clifford_decompose_layers.py index 4c66f8e6bb0c..4d4df02391b1 100644 --- a/qiskit/synthesis/clifford/clifford_decompose_layers.py +++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py @@ -36,6 +36,7 @@ _compute_rank, _gauss_elimination, _gauss_elimination_with_perm, + _binary_matmul, ) @@ -279,7 +280,7 @@ def _decompose_graph_state(cliff, validate, cz_synth_func): stabx = cliff.stab_x stabz = cliff.stab_z stabx_inv = calc_inverse_matrix(stabx, validate) - stabz_update = np.matmul(stabx_inv, stabz) % 2 + stabz_update = _binary_matmul(stabx_inv, stabz) # Assert that stabz_update is a symmetric matrix. if validate: @@ -341,7 +342,7 @@ def _decompose_hadamard_free( if not (stabx == np.zeros((num_qubits, num_qubits))).all(): raise QiskitError("The given Clifford is not Hadamard-free.") - destabz_update = np.matmul(calc_inverse_matrix(destabx), destabz) % 2 + destabz_update = _binary_matmul(calc_inverse_matrix(destabx), destabz) # Assert that destabz_update is a symmetric matrix. if validate: if (destabz_update != destabz_update.T).any(): diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index df4cb18a8a3e..6f3dfc2ec9d4 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -21,6 +21,8 @@ _gauss_elimination, _gauss_elimination_with_perm, _compute_rank_after_gauss_elim, + calc_inverse_matrix, + _binary_matmul, _row_op, _col_op, ) @@ -66,44 +68,6 @@ def random_invertible_binary_matrix( return mat -def calc_inverse_matrix(mat: np.ndarray, verify: bool = False): - """Given a square numpy(dtype=int) matrix mat, tries to compute its inverse. - - Args: - mat: a boolean square matrix. - verify: if True asserts that the multiplication of mat and its inverse is the identity matrix. - - Returns: - np.ndarray: the inverse boolean matrix. - - Raises: - QiskitError: if the matrix is not square. - QiskitError: if the matrix is not invertible. - """ - - if mat.shape[0] != mat.shape[1]: - raise QiskitError("Matrix to invert is a non-square matrix.") - - n = mat.shape[0] - # concatenate the matrix and identity - mat1 = np.concatenate((mat.astype(bool), np.eye(n, dtype=bool)), axis=1) - _gauss_elimination(mat1, None, full_elim=True) - - r = _compute_rank_after_gauss_elim(mat1[:, 0:n]) - - if r < n: - raise QiskitError("The matrix is not invertible.") - - mat1 = mat1.astype(int) - matinv = mat1[:, n : 2 * n] - - if verify: - mat2 = np.dot(mat, matinv) % 2 - assert np.array_equal(mat2, np.eye(n)) - - return matinv - - def _compute_rank(mat): """Given a boolean matrix A computes its rank""" mat = mat.astype(bool) diff --git a/test/python/synthesis/test_linear_synthesis.py b/test/python/synthesis/test_linear_synthesis.py index 98a49b6642f7..887ba94a451a 100644 --- a/test/python/synthesis/test_linear_synthesis.py +++ b/test/python/synthesis/test_linear_synthesis.py @@ -25,6 +25,7 @@ check_invertible_binary_matrix, calc_inverse_matrix, ) +from qiskit.synthesis.linear.linear_matrix_utils import _binary_matmul from qiskit.synthesis.linear.linear_circuits_utils import transpose_cx_circ, optimize_cx_4_options from test import QiskitTestCase # pylint: disable=wrong-import-order @@ -107,8 +108,9 @@ def test_invertible_matrix(self, n): """Test the functions for generating a random invertible matrix and inverting it.""" mat = random_invertible_binary_matrix(n, seed=1234) out = check_invertible_binary_matrix(mat) + mat = mat.astype(bool) mat_inv = calc_inverse_matrix(mat, verify=True) - mat_out = np.dot(mat, mat_inv) % 2 + mat_out = _binary_matmul(mat, mat_inv) self.assertTrue(np.array_equal(mat_out, np.eye(n))) self.assertTrue(out) From a333f3693ea4a7831b556c95d58ebf1abf55f45c Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 20 Jun 2024 02:33:13 -0500 Subject: [PATCH 20/43] fix failing tests, by changing mat type from int to bool --- qiskit/synthesis/stabilizer/stabilizer_decompose.py | 2 +- qiskit/transpiler/passes/synthesis/high_level_synthesis.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/qiskit/synthesis/stabilizer/stabilizer_decompose.py b/qiskit/synthesis/stabilizer/stabilizer_decompose.py index c43747105d04..0c667841130c 100644 --- a/qiskit/synthesis/stabilizer/stabilizer_decompose.py +++ b/qiskit/synthesis/stabilizer/stabilizer_decompose.py @@ -143,7 +143,7 @@ def _calc_pauli_diff_stabilizer(cliff, cliff_target): phase.extend(phase_stab) phase = np.array(phase, dtype=int) - A = cliff.symplectic_matrix.astype(int) + A = cliff.symplectic_matrix.astype(bool) Ainv = calc_inverse_matrix(A) # By carefully writing how X, Y, Z gates affect each qubit, all we need to compute diff --git a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py index 150874a84c7d..e52c5857dcbf 100644 --- a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py +++ b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py @@ -779,7 +779,7 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, ** use_inverted = options.get("use_inverted", False) use_transposed = options.get("use_transposed", False) - mat = high_level_object.linear.astype(int) + mat = high_level_object.linear.astype(bool) if use_transposed: mat = np.transpose(mat) @@ -831,7 +831,7 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, ** use_inverted = options.get("use_inverted", False) use_transposed = options.get("use_transposed", False) - mat = high_level_object.linear.astype(int) + mat = high_level_object.linear.astype(bool) if use_transposed: mat = np.transpose(mat) From 14d09abdb1835e4ad341338f245925baf8eec065 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 20 Jun 2024 07:18:05 -0500 Subject: [PATCH 21/43] update rust docstrings --- crates/accelerate/src/linear_matrix.rs | 40 +++++++++++++++++--------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 5fd5fe9cc4d2..81afaa126c74 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -14,7 +14,7 @@ use ndarray::{concatenate, s, Array2, ArrayView2, ArrayViewMut2, Axis}; use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2, PyReadwriteArray2}; use pyo3::prelude::*; -// Binary matrix multiplication +/// Binary matrix multiplication fn binary_matmul(mat1: ArrayView2, mat2: ArrayView2) -> Array2 { let n1_rows = mat1.nrows(); let n1_cols = mat1.ncols(); @@ -39,10 +39,10 @@ fn binary_matmul(mat1: ArrayView2, mat2: ArrayView2) -> Array2 result } -// Gauss elimination of a matrix mat with m rows and n columns. -// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] -// Returns the matrix mat, and the permutation perm that was done on the rows during the process. -// perm[0 : rank] represents the indices of linearly independent rows in the original matrix. +/// Gauss elimination of a matrix mat with m rows and n columns. +/// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] +/// Returns the matrix mat, and the permutation perm that was done on the rows during the process. +/// perm[0 : rank] represents the indices of linearly independent rows in the original matrix. fn gauss_elimination_with_perm( mut mat: ArrayViewMut2, ncols: Option, @@ -103,8 +103,8 @@ fn gauss_elimination_with_perm( perm } -// Given a boolean matrix A after Gaussian elimination, computes its rank -// (i.e. simply the number of nonzero rows) +/// Given a boolean matrix A after Gaussian elimination, computes its rank +/// (i.e. simply the number of nonzero rows) fn compute_rank_after_gauss_elim(mat: ArrayView2) -> usize { let rank: usize = mat .axis_iter(Axis(0)) @@ -113,7 +113,7 @@ fn compute_rank_after_gauss_elim(mat: ArrayView2) -> usize { rank } -// Given a square boolean matrix mat, tries to compute its inverse. +/// Given a square boolean matrix mat, tries to compute its inverse. fn calc_inverse_matrix_inner(mat: ArrayView2, verify: bool) -> Array2 { if mat.shape()[0] != mat.shape()[1] { panic!("Matrix to invert is a non-square matrix."); @@ -146,6 +146,10 @@ fn calc_inverse_matrix_inner(mat: ArrayView2, verify: bool) -> Array2, @@ -159,9 +163,9 @@ fn _gauss_elimination_with_perm( #[pyfunction] #[pyo3(signature = (mat, ncols=None, full_elim=false))] -// Gauss elimination of a matrix mat with m rows and n columns. -// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] -// Returns the updated matrix mat. +/// Gauss elimination of a matrix mat with m rows and n columns. +/// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] +/// Returns the updated matrix mat. fn _gauss_elimination( mut mat: PyReadwriteArray2, ncols: Option, @@ -173,8 +177,8 @@ fn _gauss_elimination( #[pyfunction] #[pyo3(signature = (mat))] -// Given a boolean matrix A after Gaussian elimination, computes its rank -// (i.e. simply the number of nonzero rows) +/// Given a boolean matrix A after Gaussian elimination, computes its rank +/// (i.e. simply the number of nonzero rows) fn _compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyResult { let view = mat.as_array(); let rank = compute_rank_after_gauss_elim(view); @@ -183,7 +187,15 @@ fn _compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> Py #[pyfunction] #[pyo3(signature = (mat, verify=false))] -// Calculates the inverse matrix +/// Given a boolean matrix mat, tries to calculate its inverse matrix +/// Args: +/// mat: a boolean square matrix. +/// verify: if True asserts that the multiplication of mat and its inverse is the identity matrix. +/// Returns: +/// the inverse matrix. +/// Raises: +/// QiskitError: if the matrix is not square. +/// QiskitError: if the matrix is not invertible. pub fn calc_inverse_matrix( py: Python, mat: PyReadonlyArray2, From 8c27cae5fb5f1136cc26c716c8c069d02d7128ef Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 20 Jun 2024 07:32:50 -0500 Subject: [PATCH 22/43] add function _add_row_or_col to rust code --- crates/accelerate/src/linear_matrix.rs | 27 +++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 81afaa126c74..142a253cecaa 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -206,6 +206,21 @@ pub fn calc_inverse_matrix( Ok(invmat.into_pyarray_bound(py).unbind()) } +/// Mutate a matrix inplace by adding the value of the ``ctrl`` row to the +/// ``target`` row. If ``add_cols`` is true, add columns instead of rows. +fn _add_row_or_col(mut mat: ArrayViewMut2, add_cols: &bool, ctrl: usize, trgt: usize) { + // get the two rows (or columns) + let info = if *add_cols { + (s![.., ctrl], s![.., trgt]) + } else { + (s![ctrl, ..], s![trgt, ..]) + }; + let (row0, mut row1) = mat.multi_slice_mut(info); + + // add them inplace + row1.zip_mut_with(&row0, |x, &y| *x ^= y); +} + #[pyfunction] #[pyo3(signature = (mat1, mat2))] // Binary matrix multiplication @@ -224,18 +239,16 @@ pub fn _binary_matmul( #[pyo3(signature = (mat, ctrl, trgt))] // Perform ROW operation on a matrix mat fn _row_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { - let mut matmut = mat.as_array_mut(); - let (row0, mut row1) = matmut.multi_slice_mut((s![ctrl, ..], s![trgt, ..])); - row1.zip_mut_with(&row0, |x, &y| *x ^= y); + let matmut = mat.as_array_mut(); + _add_row_or_col(matmut, &false, ctrl, trgt) } #[pyfunction] #[pyo3(signature = (mat, ctrl, trgt))] -// Perform COL operation on a matrix mat +// Perform COL operation on a matrix mat (in the inverse direction) fn _col_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { - let mut matmut = mat.as_array_mut(); - let (col0, mut col1) = matmut.multi_slice_mut((s![.., trgt], s![.., ctrl])); - col1.zip_mut_with(&col0, |x, &y| *x ^= y); + let matmut = mat.as_array_mut(); + _add_row_or_col(matmut, &true, trgt, ctrl) } #[pymodule] From d9b2b4b10805db2a978aa2a234d2f4c43767e745 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Mon, 24 Jun 2024 04:07:03 -0500 Subject: [PATCH 23/43] improve binary_matmul --- crates/accelerate/src/linear_matrix.rs | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/crates/accelerate/src/linear_matrix.rs b/crates/accelerate/src/linear_matrix.rs index 142a253cecaa..c9764d979051 100644 --- a/crates/accelerate/src/linear_matrix.rs +++ b/crates/accelerate/src/linear_matrix.rs @@ -27,16 +27,11 @@ fn binary_matmul(mat1: ArrayView2, mat2: ArrayView2) -> Array2 ); } - let mut result: Array2 = Array2::from_shape_fn((n1_rows, n2_cols), |(_i, _j)| false); - for i in 0..n1_rows { - for j in 0..n2_cols { - for k in 0..n2_rows { - result[(i, j)] ^= mat1[(i, k)] & mat2[(k, j)]; - } - } - } - - result + Array2::from_shape_fn((n1_rows, n2_cols), |(i, j)| { + (0..n2_rows) + .map(|k| mat1[[i, k]] & mat2[[k, j]]) + .fold(false, |acc, v| acc ^ v) + }) } /// Gauss elimination of a matrix mat with m rows and n columns. From e0d0e08d1c3483d607f54f197404bbe909f6d1aa Mon Sep 17 00:00:00 2001 From: AlexanderIvrii Date: Mon, 24 Jun 2024 14:02:29 +0300 Subject: [PATCH 24/43] proper error handling --- crates/accelerate/src/synthesis/linear/mod.rs | 13 +++++---- .../accelerate/src/synthesis/linear/utils.rs | 28 +++++++++++-------- 2 files changed, 25 insertions(+), 16 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index 9ca009ef468b..db640e7dd599 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -10,6 +10,7 @@ // copyright notice, and modified files need to carry a notice indicating // that they have been altered from the originals. +use crate::QiskitError; use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2, PyReadwriteArray2}; use pyo3::prelude::*; @@ -73,13 +74,14 @@ pub fn calc_inverse_matrix( verify: Option, ) -> PyResult>> { let view = mat.as_array(); - let invmat = utils::calc_inverse_matrix_inner(view, verify.is_some()); + let invmat = utils::calc_inverse_matrix_inner(view, verify.is_some()) + .map_err(|_| QiskitError::new_err("Inverse matrix computation failed."))?; Ok(invmat.into_pyarray_bound(py).unbind()) } #[pyfunction] #[pyo3(signature = (mat1, mat2))] -// Binary matrix multiplication +/// Binary matrix multiplication pub fn _binary_matmul( py: Python, mat1: PyReadonlyArray2, @@ -87,13 +89,14 @@ pub fn _binary_matmul( ) -> PyResult>> { let view1 = mat1.as_array(); let view2 = mat2.as_array(); - let result = utils::binary_matmul(view1, view2); + let result = utils::binary_matmul(view1, view2) + .map_err(|_| QiskitError::new_err("Binary matrix multiplication failed."))?; Ok(result.into_pyarray_bound(py).unbind()) } #[pyfunction] #[pyo3(signature = (mat, ctrl, trgt))] -// Perform ROW operation on a matrix mat +/// Perform ROW operation on a matrix mat fn _row_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { let matmut = mat.as_array_mut(); utils::_add_row_or_col(matmut, &false, ctrl, trgt) @@ -101,7 +104,7 @@ fn _row_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { #[pyfunction] #[pyo3(signature = (mat, ctrl, trgt))] -// Perform COL operation on a matrix mat (in the inverse direction) +/// Perform COL operation on a matrix mat (in the inverse direction) fn _col_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { let matmut = mat.as_array_mut(); utils::_add_row_or_col(matmut, &true, trgt, ctrl) diff --git a/crates/accelerate/src/synthesis/linear/utils.rs b/crates/accelerate/src/synthesis/linear/utils.rs index 93c24dbda6bd..c97966849c21 100644 --- a/crates/accelerate/src/synthesis/linear/utils.rs +++ b/crates/accelerate/src/synthesis/linear/utils.rs @@ -13,23 +13,26 @@ use ndarray::{concatenate, s, Array2, ArrayView2, ArrayViewMut2, Axis}; /// Binary matrix multiplication -pub fn binary_matmul(mat1: ArrayView2, mat2: ArrayView2) -> Array2 { +pub fn binary_matmul( + mat1: ArrayView2, + mat2: ArrayView2, +) -> Result, String> { let n1_rows = mat1.nrows(); let n1_cols = mat1.ncols(); let n2_rows = mat2.nrows(); let n2_cols = mat2.ncols(); if n1_cols != n2_rows { - panic!( + return Err(format!( "Cannot multiply matrices with inappropriate dimensions {}, {}", n1_cols, n2_rows - ); + )); } - Array2::from_shape_fn((n1_rows, n2_cols), |(i, j)| { + Ok(Array2::from_shape_fn((n1_rows, n2_cols), |(i, j)| { (0..n2_rows) .map(|k| mat1[[i, k]] & mat2[[k, j]]) .fold(false, |acc, v| acc ^ v) - }) + })) } /// Gauss elimination of a matrix mat with m rows and n columns. @@ -107,9 +110,12 @@ pub fn compute_rank_after_gauss_elim(mat: ArrayView2) -> usize { } /// Given a square boolean matrix mat, tries to compute its inverse. -pub fn calc_inverse_matrix_inner(mat: ArrayView2, verify: bool) -> Array2 { +pub fn calc_inverse_matrix_inner( + mat: ArrayView2, + verify: bool, +) -> Result, String> { if mat.shape()[0] != mat.shape()[1] { - panic!("Matrix to invert is a non-square matrix."); + return Err("Matrix to invert is a non-square matrix.".to_string()); } let n = mat.shape()[0]; @@ -121,20 +127,20 @@ pub fn calc_inverse_matrix_inner(mat: ArrayView2, verify: bool) -> Array2< let r = compute_rank_after_gauss_elim(mat1.slice(s![.., 0..n])); if r < n { - panic!("The matrix is not invertible."); + return Err("The matrix is not invertible.".to_string()); } let invmat = mat1.slice(s![.., n..2 * n]).to_owned(); if verify { - let mat2 = binary_matmul(mat, (&invmat).into()); + let mat2 = binary_matmul(mat, (&invmat).into())?; let identity_matrix: Array2 = Array2::from_shape_fn((n, n), |(i, j)| i == j); if mat2.ne(&identity_matrix) { - panic!("The inverse matrix is not correct."); + return Err("The inverse matrix is not correct.".to_string()); } } - invmat + Ok(invmat) } /// Mutate a matrix inplace by adding the value of the ``ctrl`` row to the From 5b1b3edadb0c3189911e43216b19fd27a45f9fe5 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Mon, 24 Jun 2024 06:08:01 -0500 Subject: [PATCH 25/43] unified format of function names --- crates/accelerate/src/synthesis/linear/mod.rs | 32 +++++++++---------- .../accelerate/src/synthesis/linear/utils.rs | 12 +++---- .../clifford/clifford_decompose_layers.py | 26 +++++++-------- qiskit/synthesis/linear/__init__.py | 1 + qiskit/synthesis/linear/linear_depth_lnn.py | 12 +++---- .../synthesis/linear/linear_matrix_utils.py | 21 ++++++------ .../python/synthesis/test_linear_synthesis.py | 4 +-- 7 files changed, 54 insertions(+), 54 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index 9ca009ef468b..7e538b5c9ff7 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -21,14 +21,14 @@ mod utils; /// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] /// Returns the matrix mat, and the permutation perm that was done on the rows during the process. /// perm[0 : rank] represents the indices of linearly independent rows in the original matrix. -fn _gauss_elimination_with_perm( +fn gauss_elimination_with_perm( py: Python, mut mat: PyReadwriteArray2, ncols: Option, full_elim: Option, ) -> PyResult { let matmut = mat.as_array_mut(); - let perm = utils::gauss_elimination_with_perm(matmut, ncols, full_elim); + let perm = utils::gauss_elimination_with_perm_inner(matmut, ncols, full_elim); Ok(perm.to_object(py)) } @@ -37,22 +37,22 @@ fn _gauss_elimination_with_perm( /// Gauss elimination of a matrix mat with m rows and n columns. /// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] /// Returns the updated matrix mat. -fn _gauss_elimination( +fn gauss_elimination( mut mat: PyReadwriteArray2, ncols: Option, full_elim: Option, ) { let matmut = mat.as_array_mut(); - let _perm = utils::gauss_elimination_with_perm(matmut, ncols, full_elim); + let _perm = utils::gauss_elimination_with_perm_inner(matmut, ncols, full_elim); } #[pyfunction] #[pyo3(signature = (mat))] /// Given a boolean matrix A after Gaussian elimination, computes its rank /// (i.e. simply the number of nonzero rows) -fn _compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyResult { +fn compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyResult { let view = mat.as_array(); - let rank = utils::compute_rank_after_gauss_elim(view); + let rank = utils::compute_rank_after_gauss_elim_inner(view); Ok(rank.to_object(py)) } @@ -80,21 +80,21 @@ pub fn calc_inverse_matrix( #[pyfunction] #[pyo3(signature = (mat1, mat2))] // Binary matrix multiplication -pub fn _binary_matmul( +pub fn binary_matmul( py: Python, mat1: PyReadonlyArray2, mat2: PyReadonlyArray2, ) -> PyResult>> { let view1 = mat1.as_array(); let view2 = mat2.as_array(); - let result = utils::binary_matmul(view1, view2); + let result = utils::binary_matmul_inner(view1, view2); Ok(result.into_pyarray_bound(py).unbind()) } #[pyfunction] #[pyo3(signature = (mat, ctrl, trgt))] // Perform ROW operation on a matrix mat -fn _row_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { +fn row_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { let matmut = mat.as_array_mut(); utils::_add_row_or_col(matmut, &false, ctrl, trgt) } @@ -102,19 +102,19 @@ fn _row_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { #[pyfunction] #[pyo3(signature = (mat, ctrl, trgt))] // Perform COL operation on a matrix mat (in the inverse direction) -fn _col_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { +fn col_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { let matmut = mat.as_array_mut(); utils::_add_row_or_col(matmut, &true, trgt, ctrl) } #[pymodule] pub fn linear(m: &Bound) -> PyResult<()> { - m.add_wrapped(wrap_pyfunction!(_gauss_elimination_with_perm))?; - m.add_wrapped(wrap_pyfunction!(_gauss_elimination))?; - m.add_wrapped(wrap_pyfunction!(_compute_rank_after_gauss_elim))?; + m.add_wrapped(wrap_pyfunction!(gauss_elimination_with_perm))?; + m.add_wrapped(wrap_pyfunction!(gauss_elimination))?; + m.add_wrapped(wrap_pyfunction!(compute_rank_after_gauss_elim))?; m.add_wrapped(wrap_pyfunction!(calc_inverse_matrix))?; - m.add_wrapped(wrap_pyfunction!(_row_op))?; - m.add_wrapped(wrap_pyfunction!(_col_op))?; - m.add_wrapped(wrap_pyfunction!(_binary_matmul))?; + m.add_wrapped(wrap_pyfunction!(row_op))?; + m.add_wrapped(wrap_pyfunction!(col_op))?; + m.add_wrapped(wrap_pyfunction!(binary_matmul))?; Ok(()) } diff --git a/crates/accelerate/src/synthesis/linear/utils.rs b/crates/accelerate/src/synthesis/linear/utils.rs index 93c24dbda6bd..9a16025e1dd8 100644 --- a/crates/accelerate/src/synthesis/linear/utils.rs +++ b/crates/accelerate/src/synthesis/linear/utils.rs @@ -13,7 +13,7 @@ use ndarray::{concatenate, s, Array2, ArrayView2, ArrayViewMut2, Axis}; /// Binary matrix multiplication -pub fn binary_matmul(mat1: ArrayView2, mat2: ArrayView2) -> Array2 { +pub fn binary_matmul_inner(mat1: ArrayView2, mat2: ArrayView2) -> Array2 { let n1_rows = mat1.nrows(); let n1_cols = mat1.ncols(); let n2_rows = mat2.nrows(); @@ -36,7 +36,7 @@ pub fn binary_matmul(mat1: ArrayView2, mat2: ArrayView2) -> Array2, ncols: Option, full_elim: Option, @@ -98,7 +98,7 @@ pub fn gauss_elimination_with_perm( /// Given a boolean matrix A after Gaussian elimination, computes its rank /// (i.e. simply the number of nonzero rows) -pub fn compute_rank_after_gauss_elim(mat: ArrayView2) -> usize { +pub fn compute_rank_after_gauss_elim_inner(mat: ArrayView2) -> usize { let rank: usize = mat .axis_iter(Axis(0)) .map(|row| row.fold(false, |out, val| out | *val) as usize) @@ -117,9 +117,9 @@ pub fn calc_inverse_matrix_inner(mat: ArrayView2, verify: bool) -> Array2< let identity_matrix: Array2 = Array2::from_shape_fn((n, n), |(i, j)| i == j); let mut mat1 = concatenate(Axis(1), &[mat.view(), identity_matrix.view()]).unwrap(); - gauss_elimination_with_perm(mat1.view_mut(), None, Some(true)); + gauss_elimination_with_perm_inner(mat1.view_mut(), None, Some(true)); - let r = compute_rank_after_gauss_elim(mat1.slice(s![.., 0..n])); + let r = compute_rank_after_gauss_elim_inner(mat1.slice(s![.., 0..n])); if r < n { panic!("The matrix is not invertible."); } @@ -127,7 +127,7 @@ pub fn calc_inverse_matrix_inner(mat: ArrayView2, verify: bool) -> Array2< let invmat = mat1.slice(s![.., n..2 * n]).to_owned(); if verify { - let mat2 = binary_matmul(mat, (&invmat).into()); + let mat2 = binary_matmul_inner(mat, (&invmat).into()); let identity_matrix: Array2 = Array2::from_shape_fn((n, n), |(i, j)| i == j); if mat2.ne(&identity_matrix) { panic!("The inverse matrix is not correct."); diff --git a/qiskit/synthesis/clifford/clifford_decompose_layers.py b/qiskit/synthesis/clifford/clifford_decompose_layers.py index 980a8c01b179..82ccadd4a1ce 100644 --- a/qiskit/synthesis/clifford/clifford_decompose_layers.py +++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py @@ -33,10 +33,10 @@ from qiskit.synthesis.linear_phase import synth_cz_depth_line_mr, synth_cx_cz_depth_line_my from qiskit.synthesis.linear.linear_matrix_utils import ( calc_inverse_matrix, - _compute_rank, - _gauss_elimination, - _gauss_elimination_with_perm, - _binary_matmul, + compute_rank, + gauss_elimination, + gauss_elimination_with_perm, + binary_matmul, ) @@ -204,25 +204,25 @@ def _create_graph_state(cliff, validate=False): """ num_qubits = cliff.num_qubits - rank = _compute_rank(cliff.stab_x) + rank = compute_rank(cliff.stab_x) H1_circ = QuantumCircuit(num_qubits, name="H1") cliffh = cliff.copy() if rank < num_qubits: stab = cliff.stab[:, :-1] stab = stab.astype(bool) - _gauss_elimination(stab, num_qubits) + gauss_elimination(stab, num_qubits) Cmat = stab[rank:num_qubits, num_qubits:] Cmat = np.transpose(Cmat) - perm = _gauss_elimination_with_perm(Cmat) + perm = gauss_elimination_with_perm(Cmat) perm = perm[0 : num_qubits - rank] # validate that the output matrix has the same rank if validate: - if _compute_rank(Cmat) != num_qubits - rank: + if compute_rank(Cmat) != num_qubits - rank: raise QiskitError("The matrix Cmat after Gauss elimination has wrong rank.") - if _compute_rank(stab[:, 0:num_qubits]) != rank: + if compute_rank(stab[:, 0:num_qubits]) != rank: raise QiskitError("The matrix after Gauss elimination has wrong rank.") # validate that we have a num_qubits - rank zero rows for i in range(rank, num_qubits): @@ -239,7 +239,7 @@ def _create_graph_state(cliff, validate=False): # validate that a layer of Hadamard gates and then appending cliff, provides a graph state. if validate: stabh = cliffh.stab_x - if _compute_rank(stabh) != num_qubits: + if compute_rank(stabh) != num_qubits: raise QiskitError("The state is not a graph state.") return H1_circ, cliffh @@ -269,7 +269,7 @@ def _decompose_graph_state(cliff, validate, cz_synth_func): """ num_qubits = cliff.num_qubits - rank = _compute_rank(cliff.stab_x) + rank = compute_rank(cliff.stab_x) cliff_cpy = cliff.copy() if rank < num_qubits: raise QiskitError("The stabilizer state is not a graph state.") @@ -280,7 +280,7 @@ def _decompose_graph_state(cliff, validate, cz_synth_func): stabx = cliff.stab_x stabz = cliff.stab_z stabx_inv = calc_inverse_matrix(stabx, validate) - stabz_update = _binary_matmul(stabx_inv, stabz) + stabz_update = binary_matmul(stabx_inv, stabz) # Assert that stabz_update is a symmetric matrix. if validate: @@ -342,7 +342,7 @@ def _decompose_hadamard_free( if not (stabx == np.zeros((num_qubits, num_qubits))).all(): raise QiskitError("The given Clifford is not Hadamard-free.") - destabz_update = _binary_matmul(calc_inverse_matrix(destabx), destabz) + destabz_update = binary_matmul(calc_inverse_matrix(destabx), destabz) # Assert that destabz_update is a symmetric matrix. if validate: if (destabz_update != destabz_update.T).any(): diff --git a/qiskit/synthesis/linear/__init__.py b/qiskit/synthesis/linear/__init__.py index 115fc557bfa4..f3537de9c3ff 100644 --- a/qiskit/synthesis/linear/__init__.py +++ b/qiskit/synthesis/linear/__init__.py @@ -18,6 +18,7 @@ random_invertible_binary_matrix, calc_inverse_matrix, check_invertible_binary_matrix, + binary_matmul, ) # This is re-import is kept for compatibility with Terra 0.23. Eligible for deprecation in 0.25+. diff --git a/qiskit/synthesis/linear/linear_depth_lnn.py b/qiskit/synthesis/linear/linear_depth_lnn.py index 33b1b0ae21e7..2eab7d3e8bff 100644 --- a/qiskit/synthesis/linear/linear_depth_lnn.py +++ b/qiskit/synthesis/linear/linear_depth_lnn.py @@ -28,15 +28,15 @@ from qiskit.synthesis.linear.linear_matrix_utils import ( calc_inverse_matrix, check_invertible_binary_matrix, - _col_op, - _row_op, + col_op, + row_op, ) def _row_op_update_instructions(cx_instructions, mat, a, b): # Add a cx gate to the instructions and update the matrix mat cx_instructions.append((a, b)) - _row_op(mat, a, b) + row_op(mat, a, b) def _get_lower_triangular(n, mat, mat_inv): @@ -62,7 +62,7 @@ def _get_lower_triangular(n, mat, mat_inv): first_j = j else: # cx_instructions_cols (L instructions) are not needed - _col_op(mat, j, first_j) + col_op(mat, j, first_j) # Use row operations directed upwards to zero out all "1"s above the remaining "1" in row i for k in reversed(range(0, i)): if mat[k, first_j]: @@ -70,8 +70,8 @@ def _get_lower_triangular(n, mat, mat_inv): # Apply only U instructions to get the permuted L for inst in cx_instructions_rows: - _row_op(mat_t, inst[0], inst[1]) - _col_op(mat_inv_t, inst[0], inst[1]) + row_op(mat_t, inst[0], inst[1]) + col_op(mat_inv_t, inst[0], inst[1]) return mat_t, mat_inv_t diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index dd1b11a28c6e..31cadee5db93 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -14,17 +14,16 @@ from typing import Optional, Union import numpy as np -from qiskit.exceptions import QiskitError # pylint: disable=unused-import from qiskit._accelerate.synthesis.linear import ( - _gauss_elimination, - _gauss_elimination_with_perm, - _compute_rank_after_gauss_elim, + gauss_elimination, + gauss_elimination_with_perm, + compute_rank_after_gauss_elim, calc_inverse_matrix, - _binary_matmul, - _row_op, - _col_op, + binary_matmul, + row_op, + col_op, ) @@ -40,7 +39,7 @@ def check_invertible_binary_matrix(mat: np.ndarray): if len(mat.shape) != 2 or mat.shape[0] != mat.shape[1]: return False - rank = _compute_rank(mat) + rank = compute_rank(mat) return rank == mat.shape[0] @@ -64,12 +63,12 @@ def random_invertible_binary_matrix( rank = 0 while rank != num_qubits: mat = rng.integers(2, size=(num_qubits, num_qubits)) - rank = _compute_rank(mat) + rank = compute_rank(mat) return mat -def _compute_rank(mat): +def compute_rank(mat): """Given a boolean matrix A computes its rank""" mat = mat.astype(bool) - _gauss_elimination(mat) + gauss_elimination(mat) return np.sum(mat.any(axis=1)) diff --git a/test/python/synthesis/test_linear_synthesis.py b/test/python/synthesis/test_linear_synthesis.py index 887ba94a451a..372fdab3523d 100644 --- a/test/python/synthesis/test_linear_synthesis.py +++ b/test/python/synthesis/test_linear_synthesis.py @@ -24,8 +24,8 @@ random_invertible_binary_matrix, check_invertible_binary_matrix, calc_inverse_matrix, + binary_matmul, ) -from qiskit.synthesis.linear.linear_matrix_utils import _binary_matmul from qiskit.synthesis.linear.linear_circuits_utils import transpose_cx_circ, optimize_cx_4_options from test import QiskitTestCase # pylint: disable=wrong-import-order @@ -110,7 +110,7 @@ def test_invertible_matrix(self, n): out = check_invertible_binary_matrix(mat) mat = mat.astype(bool) mat_inv = calc_inverse_matrix(mat, verify=True) - mat_out = _binary_matmul(mat, mat_inv) + mat_out = binary_matmul(mat, mat_inv) self.assertTrue(np.array_equal(mat_out, np.eye(n))) self.assertTrue(out) From a88ddba08ef7842102addb29061b2bde7622f091 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Mon, 24 Jun 2024 09:11:53 -0500 Subject: [PATCH 26/43] move compute_rank from python to rust, update errors --- crates/accelerate/src/synthesis/linear/mod.rs | 21 +++++++++++++------ .../accelerate/src/synthesis/linear/utils.rs | 7 +++++++ .../clifford/clifford_decompose_layers.py | 6 +++--- .../synthesis/linear/linear_matrix_utils.py | 15 +++++-------- 4 files changed, 30 insertions(+), 19 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index 137ed0a5d78c..d8a04174d5eb 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -49,7 +49,7 @@ fn gauss_elimination( #[pyfunction] #[pyo3(signature = (mat))] -/// Given a boolean matrix A after Gaussian elimination, computes its rank +/// Given a boolean matrix mat after Gaussian elimination, computes its rank /// (i.e. simply the number of nonzero rows) fn compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyResult { let view = mat.as_array(); @@ -57,6 +57,15 @@ fn compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyR Ok(rank.to_object(py)) } +#[pyfunction] +#[pyo3(signature = (mat))] +/// Given a boolean matrix mat computes its rank +fn compute_rank(py: Python, mut mat: PyReadwriteArray2) -> PyResult { + let matmut = mat.as_array_mut(); + let rank = utils::compute_rank_inner(matmut); + Ok(rank.to_object(py)) +} + #[pyfunction] #[pyo3(signature = (mat, verify=false))] /// Given a boolean matrix mat, tries to calculate its inverse matrix @@ -66,8 +75,7 @@ fn compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyR /// Returns: /// the inverse matrix. /// Raises: -/// QiskitError: if the matrix is not square. -/// QiskitError: if the matrix is not invertible. +/// QiskitError: if the matrix is not square or not invertible. pub fn calc_inverse_matrix( py: Python, mat: PyReadonlyArray2, @@ -75,7 +83,7 @@ pub fn calc_inverse_matrix( ) -> PyResult>> { let view = mat.as_array(); let invmat = utils::calc_inverse_matrix_inner(view, verify.is_some()) - .map_err(|_| QiskitError::new_err("Inverse matrix computation failed."))?; + .map_err(QiskitError::new_err)?; Ok(invmat.into_pyarray_bound(py).unbind()) } @@ -89,8 +97,8 @@ pub fn binary_matmul( ) -> PyResult>> { let view1 = mat1.as_array(); let view2 = mat2.as_array(); - let result = utils::binary_matmul_inner(view1, view2) - .map_err(|_| QiskitError::new_err("Binary matrix multiplication failed."))?; + let result = + utils::binary_matmul_inner(view1, view2).map_err(QiskitError::new_err)?; Ok(result.into_pyarray_bound(py).unbind()) } @@ -115,6 +123,7 @@ pub fn linear(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(gauss_elimination_with_perm))?; m.add_wrapped(wrap_pyfunction!(gauss_elimination))?; m.add_wrapped(wrap_pyfunction!(compute_rank_after_gauss_elim))?; + m.add_wrapped(wrap_pyfunction!(compute_rank))?; m.add_wrapped(wrap_pyfunction!(calc_inverse_matrix))?; m.add_wrapped(wrap_pyfunction!(row_op))?; m.add_wrapped(wrap_pyfunction!(col_op))?; diff --git a/crates/accelerate/src/synthesis/linear/utils.rs b/crates/accelerate/src/synthesis/linear/utils.rs index 405080333ae4..999aaa59212a 100644 --- a/crates/accelerate/src/synthesis/linear/utils.rs +++ b/crates/accelerate/src/synthesis/linear/utils.rs @@ -109,6 +109,13 @@ pub fn compute_rank_after_gauss_elim_inner(mat: ArrayView2) -> usize { rank } +/// Given a boolean matrix mat computes its rank +pub fn compute_rank_inner(mut mat: ArrayViewMut2) -> usize { + gauss_elimination_with_perm_inner(mat.view_mut(), None, Some(false)); + let rank = compute_rank_after_gauss_elim_inner(mat.view()); + rank +} + /// Given a square boolean matrix mat, tries to compute its inverse. pub fn calc_inverse_matrix_inner( mat: ArrayView2, diff --git a/qiskit/synthesis/clifford/clifford_decompose_layers.py b/qiskit/synthesis/clifford/clifford_decompose_layers.py index 82ccadd4a1ce..03b5ed216d9c 100644 --- a/qiskit/synthesis/clifford/clifford_decompose_layers.py +++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py @@ -204,7 +204,7 @@ def _create_graph_state(cliff, validate=False): """ num_qubits = cliff.num_qubits - rank = compute_rank(cliff.stab_x) + rank = compute_rank((cliff.stab_x).astype(bool)) H1_circ = QuantumCircuit(num_qubits, name="H1") cliffh = cliff.copy() @@ -238,7 +238,7 @@ def _create_graph_state(cliff, validate=False): # validate that a layer of Hadamard gates and then appending cliff, provides a graph state. if validate: - stabh = cliffh.stab_x + stabh = (cliffh.stab_x).astype(bool) if compute_rank(stabh) != num_qubits: raise QiskitError("The state is not a graph state.") @@ -269,7 +269,7 @@ def _decompose_graph_state(cliff, validate, cz_synth_func): """ num_qubits = cliff.num_qubits - rank = compute_rank(cliff.stab_x) + rank = compute_rank((cliff.stab_x).astype(bool)) cliff_cpy = cliff.copy() if rank < num_qubits: raise QiskitError("The stabilizer state is not a graph state.") diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index 31cadee5db93..4be1ef1d64b1 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -20,6 +20,7 @@ gauss_elimination, gauss_elimination_with_perm, compute_rank_after_gauss_elim, + compute_rank, calc_inverse_matrix, binary_matmul, row_op, @@ -39,8 +40,9 @@ def check_invertible_binary_matrix(mat: np.ndarray): if len(mat.shape) != 2 or mat.shape[0] != mat.shape[1]: return False - rank = compute_rank(mat) - return rank == mat.shape[0] + matcpy = mat.copy() + rank = compute_rank(matcpy) + return rank == matcpy.shape[0] def random_invertible_binary_matrix( @@ -62,13 +64,6 @@ def random_invertible_binary_matrix( rank = 0 while rank != num_qubits: - mat = rng.integers(2, size=(num_qubits, num_qubits)) + mat = rng.integers(2, size=(num_qubits, num_qubits), dtype=bool) rank = compute_rank(mat) return mat - - -def compute_rank(mat): - """Given a boolean matrix A computes its rank""" - mat = mat.astype(bool) - gauss_elimination(mat) - return np.sum(mat.any(axis=1)) From 2c0601f774b1bba78f931d734f0820e29de0ddb5 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Mon, 24 Jun 2024 10:32:19 -0500 Subject: [PATCH 27/43] update type of mat in compute_rank --- crates/accelerate/src/synthesis/linear/mod.rs | 12 +++++------- crates/accelerate/src/synthesis/linear/utils.rs | 7 ++++--- qiskit/synthesis/linear/linear_matrix_utils.py | 5 ++--- 3 files changed, 11 insertions(+), 13 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index d8a04174d5eb..695cd150ca36 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -60,9 +60,8 @@ fn compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyR #[pyfunction] #[pyo3(signature = (mat))] /// Given a boolean matrix mat computes its rank -fn compute_rank(py: Python, mut mat: PyReadwriteArray2) -> PyResult { - let matmut = mat.as_array_mut(); - let rank = utils::compute_rank_inner(matmut); +fn compute_rank(py: Python, mat: PyReadonlyArray2) -> PyResult { + let rank = utils::compute_rank_inner(mat.as_array()); Ok(rank.to_object(py)) } @@ -82,8 +81,8 @@ pub fn calc_inverse_matrix( verify: Option, ) -> PyResult>> { let view = mat.as_array(); - let invmat = utils::calc_inverse_matrix_inner(view, verify.is_some()) - .map_err(QiskitError::new_err)?; + let invmat = + utils::calc_inverse_matrix_inner(view, verify.is_some()).map_err(QiskitError::new_err)?; Ok(invmat.into_pyarray_bound(py).unbind()) } @@ -97,8 +96,7 @@ pub fn binary_matmul( ) -> PyResult>> { let view1 = mat1.as_array(); let view2 = mat2.as_array(); - let result = - utils::binary_matmul_inner(view1, view2).map_err(QiskitError::new_err)?; + let result = utils::binary_matmul_inner(view1, view2).map_err(QiskitError::new_err)?; Ok(result.into_pyarray_bound(py).unbind()) } diff --git a/crates/accelerate/src/synthesis/linear/utils.rs b/crates/accelerate/src/synthesis/linear/utils.rs index 999aaa59212a..151a6fbd0439 100644 --- a/crates/accelerate/src/synthesis/linear/utils.rs +++ b/crates/accelerate/src/synthesis/linear/utils.rs @@ -110,9 +110,10 @@ pub fn compute_rank_after_gauss_elim_inner(mat: ArrayView2) -> usize { } /// Given a boolean matrix mat computes its rank -pub fn compute_rank_inner(mut mat: ArrayViewMut2) -> usize { - gauss_elimination_with_perm_inner(mat.view_mut(), None, Some(false)); - let rank = compute_rank_after_gauss_elim_inner(mat.view()); +pub fn compute_rank_inner(mat: ArrayView2) -> usize { + let mut temp_mat = mat.to_owned(); + gauss_elimination_with_perm_inner(temp_mat.view_mut(), None, Some(false)); + let rank = compute_rank_after_gauss_elim_inner(temp_mat.view()); rank } diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index 4be1ef1d64b1..ee64a5b18617 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -40,9 +40,8 @@ def check_invertible_binary_matrix(mat: np.ndarray): if len(mat.shape) != 2 or mat.shape[0] != mat.shape[1]: return False - matcpy = mat.copy() - rank = compute_rank(matcpy) - return rank == matcpy.shape[0] + rank = compute_rank(mat) + return rank == mat.shape[0] def random_invertible_binary_matrix( From 26a95d210f8d17dea2fc1948a45d92985bcf50ee Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Tue, 25 Jun 2024 00:12:19 -0500 Subject: [PATCH 28/43] move random_invertible_binary_matrix and check_invertible_binary_matrix to rust --- crates/accelerate/src/synthesis/linear/mod.rs | 32 ++++++++++++++ .../accelerate/src/synthesis/linear/utils.rs | 33 +++++++++++++++ .../synthesis/linear/linear_matrix_utils.py | 42 +------------------ 3 files changed, 67 insertions(+), 40 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index 695cd150ca36..34c46f3a7090 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -116,6 +116,36 @@ fn col_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { utils::_add_row_or_col(matmut, &true, trgt, ctrl) } +#[pyfunction] +#[pyo3(signature = (num_qubits, seed))] +/// Generate a random invertible n x n binary matrix. +/// Args: +/// num_qubits: the matrix size. +/// seed: a random seed. +/// Returns: +/// np.ndarray: A random invertible binary matrix of size num_qubits. +fn random_invertible_binary_matrix( + py: Python, + num_qubits: usize, + seed: Option, +) -> PyResult>> { + let matrix = utils::random_invertible_binary_matrix_inner(num_qubits, seed); + Ok(matrix.into_pyarray_bound(py).unbind()) +} + +#[pyfunction] +#[pyo3(signature = (mat))] +/// Check that a binary matrix is invertible. +/// Args: +/// mat: a binary matrix. +/// Returns: +/// bool: True if mat in invertible and False otherwise. +fn check_invertible_binary_matrix(py: Python, mat: PyReadonlyArray2) -> PyResult { + let view = mat.as_array(); + let out = utils::check_invertible_binary_matrix_inner(view); + Ok(out.to_object(py)) +} + #[pymodule] pub fn linear(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(gauss_elimination_with_perm))?; @@ -126,5 +156,7 @@ pub fn linear(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(row_op))?; m.add_wrapped(wrap_pyfunction!(col_op))?; m.add_wrapped(wrap_pyfunction!(binary_matmul))?; + m.add_wrapped(wrap_pyfunction!(random_invertible_binary_matrix))?; + m.add_wrapped(wrap_pyfunction!(check_invertible_binary_matrix))?; Ok(()) } diff --git a/crates/accelerate/src/synthesis/linear/utils.rs b/crates/accelerate/src/synthesis/linear/utils.rs index 151a6fbd0439..b4dbf4993081 100644 --- a/crates/accelerate/src/synthesis/linear/utils.rs +++ b/crates/accelerate/src/synthesis/linear/utils.rs @@ -11,6 +11,8 @@ // that they have been altered from the originals. use ndarray::{concatenate, s, Array2, ArrayView2, ArrayViewMut2, Axis}; +use rand::{Rng, SeedableRng}; +use rand_pcg::Pcg64Mcg; /// Binary matrix multiplication pub fn binary_matmul_inner( @@ -165,3 +167,34 @@ pub fn _add_row_or_col(mut mat: ArrayViewMut2, add_cols: &bool, ctrl: usiz // add them inplace row1.zip_mut_with(&row0, |x, &y| *x ^= y); } + +/// Generate a random invertible n x n binary matrix. +pub fn random_invertible_binary_matrix_inner(num_qubits: usize, seed: Option) -> Array2 { + let mut rng = match seed { + Some(seed) => Pcg64Mcg::seed_from_u64(seed), + None => Pcg64Mcg::from_entropy(), + }; + + let mut matrix = Array2::from_elem((num_qubits, num_qubits), false); + + loop { + for value in matrix.iter_mut() { + *value = rng.gen_bool(0.5); + } + + let rank = compute_rank_inner(matrix.view()); + if rank == num_qubits { + break; + } + } + matrix +} + +/// Check that a binary matrix is invertible. +pub fn check_invertible_binary_matrix_inner(mat: ArrayView2) -> bool { + if mat.nrows() != mat.ncols() { + return false; + } + let rank = compute_rank_inner(mat); + rank == mat.nrows() +} diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index ee64a5b18617..eb47c0e20343 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -23,46 +23,8 @@ compute_rank, calc_inverse_matrix, binary_matmul, + random_invertible_binary_matrix, + check_invertible_binary_matrix, row_op, col_op, ) - - -def check_invertible_binary_matrix(mat: np.ndarray): - """Check that a binary matrix is invertible. - - Args: - mat: a binary matrix. - - Returns: - bool: True if mat in invertible and False otherwise. - """ - if len(mat.shape) != 2 or mat.shape[0] != mat.shape[1]: - return False - - rank = compute_rank(mat) - return rank == mat.shape[0] - - -def random_invertible_binary_matrix( - num_qubits: int, seed: Optional[Union[np.random.Generator, int]] = None -): - """Generates a random invertible n x n binary matrix. - - Args: - num_qubits: the matrix size. - seed: a random seed. - - Returns: - np.ndarray: A random invertible binary matrix of size num_qubits. - """ - if isinstance(seed, np.random.Generator): - rng = seed - else: - rng = np.random.default_rng(seed) - - rank = 0 - while rank != num_qubits: - mat = rng.integers(2, size=(num_qubits, num_qubits), dtype=bool) - rank = compute_rank(mat) - return mat From ec2ce02eceeef7788063ceeaf474d524375cdaf2 Mon Sep 17 00:00:00 2001 From: AlexanderIvrii Date: Tue, 25 Jun 2024 08:55:13 +0300 Subject: [PATCH 29/43] Updating HighLevelSynthesis tests that depend on the specific random number --- .../transpiler/test_high_level_synthesis.py | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/test/python/transpiler/test_high_level_synthesis.py b/test/python/transpiler/test_high_level_synthesis.py index ff54169374bc..a76ab08d90e1 100644 --- a/test/python/transpiler/test_high_level_synthesis.py +++ b/test/python/transpiler/test_high_level_synthesis.py @@ -534,22 +534,22 @@ def test_section_size(self): hls_config = HLSConfig(linear_function=[("pmh", {"section_size": 1})]) qct = HighLevelSynthesis(hls_config=hls_config)(qc) self.assertEqual(LinearFunction(qct), LinearFunction(qc)) - self.assertEqual(qct.size(), 22) - self.assertEqual(qct.depth(), 20) + self.assertEqual(qct.size(), 30) + self.assertEqual(qct.depth(), 27) with self.subTest("section_size_2"): hls_config = HLSConfig(linear_function=[("pmh", {"section_size": 2})]) qct = HighLevelSynthesis(hls_config=hls_config)(qc) self.assertEqual(LinearFunction(qct), LinearFunction(qc)) - self.assertEqual(qct.size(), 23) - self.assertEqual(qct.depth(), 19) + self.assertEqual(qct.size(), 27) + self.assertEqual(qct.depth(), 23) with self.subTest("section_size_3"): hls_config = HLSConfig(linear_function=[("pmh", {"section_size": 3})]) qct = HighLevelSynthesis(hls_config=hls_config)(qc) self.assertEqual(LinearFunction(qct), LinearFunction(qc)) - self.assertEqual(qct.size(), 23) - self.assertEqual(qct.depth(), 17) + self.assertEqual(qct.size(), 29) + self.assertEqual(qct.depth(), 23) def test_invert_and_transpose(self): """Test that the plugin takes the use_inverted and use_transposed arguments into account.""" @@ -623,7 +623,7 @@ def test_plugin_selection_all_with_metrix(self): # The seed is chosen so that we get different best circuits depending on whether we # want to minimize size or depth. - mat = random_invertible_binary_matrix(7, seed=37) + mat = random_invertible_binary_matrix(7, seed=38) qc = QuantumCircuit(7) qc.append(LinearFunction(mat), [0, 1, 2, 3, 4, 5, 6]) @@ -641,8 +641,8 @@ def test_plugin_selection_all_with_metrix(self): ) qct = HighLevelSynthesis(hls_config=hls_config)(qc) self.assertEqual(LinearFunction(qct), LinearFunction(qc)) - self.assertEqual(qct.size(), 20) - self.assertEqual(qct.depth(), 15) + self.assertEqual(qct.size(), 23) + self.assertEqual(qct.depth(), 19) with self.subTest("depth_fn"): # We want to minimize the "depth" (aka the number of layers) in the circuit @@ -658,8 +658,8 @@ def test_plugin_selection_all_with_metrix(self): ) qct = HighLevelSynthesis(hls_config=hls_config)(qc) self.assertEqual(LinearFunction(qct), LinearFunction(qc)) - self.assertEqual(qct.size(), 23) - self.assertEqual(qct.depth(), 12) + self.assertEqual(qct.size(), 24) + self.assertEqual(qct.depth(), 13) class TestKMSSynthesisLinearFunctionPlugin(QiskitTestCase): From 88e46dee423a8f361ecf78d1893e6b524eafcddc Mon Sep 17 00:00:00 2001 From: AlexanderIvrii Date: Tue, 25 Jun 2024 09:24:54 +0300 Subject: [PATCH 30/43] Updating LinearSynthesis tests to pass seeds --- test/python/synthesis/test_linear_synthesis.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/python/synthesis/test_linear_synthesis.py b/test/python/synthesis/test_linear_synthesis.py index 372fdab3523d..bbfab20a30fc 100644 --- a/test/python/synthesis/test_linear_synthesis.py +++ b/test/python/synthesis/test_linear_synthesis.py @@ -119,8 +119,9 @@ def test_synth_lnn_kms(self, num_qubits): """Test that synth_cnot_depth_line_kms produces the correct synthesis.""" rng = np.random.default_rng(1234) num_trials = 10 - for _ in range(num_trials): - mat = random_invertible_binary_matrix(num_qubits, seed=rng) + seeds = rng.integers(100000, size=num_trials, dtype=np.uint64) + for seed in seeds: + mat = random_invertible_binary_matrix(num_qubits, seed=seed) mat = np.array(mat, dtype=bool) qc = synth_cnot_depth_line_kms(mat) mat1 = LinearFunction(qc).linear From 634b92f943084970be3a15d7ae4446b85ad84f0f Mon Sep 17 00:00:00 2001 From: AlexanderIvrii Date: Tue, 25 Jun 2024 09:57:22 +0300 Subject: [PATCH 31/43] Updating tests in test_linear_function --- .../circuit/library/test_linear_function.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/test/python/circuit/library/test_linear_function.py b/test/python/circuit/library/test_linear_function.py index a2d868fbc1be..a3df1e9664a2 100644 --- a/test/python/circuit/library/test_linear_function.py +++ b/test/python/circuit/library/test_linear_function.py @@ -86,7 +86,8 @@ def random_linear_circuit( elif name == "linear": nqargs = rng.choice(range(1, num_qubits + 1)) qargs = rng.choice(range(num_qubits), nqargs, replace=False).tolist() - mat = random_invertible_binary_matrix(nqargs, seed=rng) + seed = rng.integers(100000, size=1, dtype=np.uint64)[0] + mat = random_invertible_binary_matrix(nqargs, seed=seed) circ.append(LinearFunction(mat), qargs) elif name == "permutation": nqargs = rng.choice(range(1, num_qubits + 1)) @@ -140,10 +141,11 @@ def test_conversion_to_matrix_and_back(self, num_qubits): and then synthesizing this linear function to a quantum circuit.""" rng = np.random.default_rng(1234) - for _ in range(10): - for num_gates in [0, 5, 5 * num_qubits]: + for num_gates in [0, 5, 5 * num_qubits]: + seeds = rng.integers(100000, size=10, dtype=np.uint64) + for seed in seeds: # create a random linear circuit - linear_circuit = random_linear_circuit(num_qubits, num_gates, seed=rng) + linear_circuit = random_linear_circuit(num_qubits, num_gates, seed=seed) self.assertIsInstance(linear_circuit, QuantumCircuit) # convert it to a linear function @@ -168,10 +170,11 @@ def test_conversion_to_linear_function_and_back(self, num_qubits): """Test correctness of first synthesizing a linear circuit from a linear function, and then converting this linear circuit to a linear function.""" rng = np.random.default_rng(5678) + seeds = rng.integers(100000, size=10, dtype=np.uint64) - for _ in range(10): + for seed in seeds: # create a random invertible binary matrix - binary_matrix = random_invertible_binary_matrix(num_qubits, seed=rng) + binary_matrix = random_invertible_binary_matrix(num_qubits, seed=seed) # create a linear function with this matrix linear_function = LinearFunction(binary_matrix, validate_input=True) From 7fdf87be764f054047b3863d9d9b271103ee3478 Mon Sep 17 00:00:00 2001 From: AlexanderIvrii Date: Tue, 25 Jun 2024 10:04:38 +0300 Subject: [PATCH 32/43] changing the matrix type in random dyhedral to be a matrix of ints rather than bools --- qiskit/quantum_info/operators/dihedral/random.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/qiskit/quantum_info/operators/dihedral/random.py b/qiskit/quantum_info/operators/dihedral/random.py index 4331d618d73d..79220c128673 100644 --- a/qiskit/quantum_info/operators/dihedral/random.py +++ b/qiskit/quantum_info/operators/dihedral/random.py @@ -53,7 +53,8 @@ def random_cnotdihedral(num_qubits, seed=None): random_invertible_binary_matrix, ) - linear = random_invertible_binary_matrix(num_qubits, seed=rng) + seed = rng.integers(100000, size=1, dtype=np.uint64)[0] + linear = random_invertible_binary_matrix(num_qubits, seed=seed).astype(int) elem.linear = linear # Random shift From 11ac8855321c694862532081d2b4c7d508c3087c Mon Sep 17 00:00:00 2001 From: AlexanderIvrii Date: Tue, 25 Jun 2024 10:06:54 +0300 Subject: [PATCH 33/43] updating cx_cz synthesis tests --- test/python/synthesis/test_cx_cz_synthesis.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/python/synthesis/test_cx_cz_synthesis.py b/test/python/synthesis/test_cx_cz_synthesis.py index 63353dab95df..28a26df181a2 100644 --- a/test/python/synthesis/test_cx_cz_synthesis.py +++ b/test/python/synthesis/test_cx_cz_synthesis.py @@ -39,8 +39,9 @@ def test_cx_cz_synth_lnn(self, num_qubits): rng = np.random.default_rng(seed) num_gates = 10 num_trials = 8 + seeds = rng.integers(100000, size=num_trials, dtype=np.uint64) - for _ in range(num_trials): + for seed in seeds: # Generate a random CZ circuit mat_z = np.zeros((num_qubits, num_qubits)) cir_z = QuantumCircuit(num_qubits) @@ -55,7 +56,7 @@ def test_cx_cz_synth_lnn(self, num_qubits): mat_z[j][i] = (mat_z[j][i] + 1) % 2 # Generate a random CX circuit - mat_x = random_invertible_binary_matrix(num_qubits, seed=rng) + mat_x = random_invertible_binary_matrix(num_qubits, seed=seed) mat_x = np.array(mat_x, dtype=bool) cir_x = synth_cnot_depth_line_kms(mat_x) From c2d63303b448d9d8afbf864e6a7daa6092be4aeb Mon Sep 17 00:00:00 2001 From: AlexanderIvrii Date: Tue, 25 Jun 2024 10:12:33 +0300 Subject: [PATCH 34/43] updating clifford tests --- .../quantum_info/operators/symplectic/test_clifford.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/python/quantum_info/operators/symplectic/test_clifford.py b/test/python/quantum_info/operators/symplectic/test_clifford.py index 043a9eca78b6..36d716d2d85a 100644 --- a/test/python/quantum_info/operators/symplectic/test_clifford.py +++ b/test/python/quantum_info/operators/symplectic/test_clifford.py @@ -424,9 +424,9 @@ def test_from_linear_function(self, num_qubits): """Test initialization from linear function.""" rng = np.random.default_rng(1234) samples = 50 - - for _ in range(samples): - mat = random_invertible_binary_matrix(num_qubits, seed=rng) + seeds = rng.integers(100000, size=samples, dtype=np.uint64) + for seed in seeds: + mat = random_invertible_binary_matrix(num_qubits, seed=seed) lin = LinearFunction(mat) cliff = Clifford(lin) self.assertTrue(Operator(cliff).equiv(Operator(lin))) From 4b8339547b58e9924e2ebfc5f8e1f8b52eb49003 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Tue, 25 Jun 2024 03:28:36 -0500 Subject: [PATCH 35/43] remove unused imports --- qiskit/synthesis/linear/linear_matrix_utils.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py index eb47c0e20343..a76efdbb8d72 100644 --- a/qiskit/synthesis/linear/linear_matrix_utils.py +++ b/qiskit/synthesis/linear/linear_matrix_utils.py @@ -12,9 +12,6 @@ """Utility functions for handling binary matrices.""" -from typing import Optional, Union -import numpy as np - # pylint: disable=unused-import from qiskit._accelerate.synthesis.linear import ( gauss_elimination, From 713d5935e79588f65d7d4eb1b57b9365392ddf65 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Tue, 25 Jun 2024 05:24:01 -0500 Subject: [PATCH 36/43] add option seed=None --- crates/accelerate/src/synthesis/linear/mod.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index 34c46f3a7090..3057c80328ae 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -117,7 +117,7 @@ fn col_op(mut mat: PyReadwriteArray2, ctrl: usize, trgt: usize) { } #[pyfunction] -#[pyo3(signature = (num_qubits, seed))] +#[pyo3(signature = (num_qubits, seed=None))] /// Generate a random invertible n x n binary matrix. /// Args: /// num_qubits: the matrix size. From 4406441886c2bc2c6c68fd9a7e050bfd1c97876a Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Tue, 25 Jun 2024 05:53:03 -0500 Subject: [PATCH 37/43] enhance rust docs --- crates/accelerate/src/synthesis/linear/mod.rs | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index 3057c80328ae..f0ede103bd75 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -22,6 +22,14 @@ mod utils; /// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] /// Returns the matrix mat, and the permutation perm that was done on the rows during the process. /// perm[0 : rank] represents the indices of linearly independent rows in the original matrix. +/// Args: +/// mat: a boolean matrix with n rows and m columns +/// ncols: the number of columns for the gaussian elimination, +/// if ncols=None, then the elimination is done over all the colums +/// full_elim: whether to do a full elimination, or partial (upper triangular form) +/// Returns: +/// mat: the matrix after the gaussian elimination +/// perm: the permutation perm that was done on the rows during the process fn gauss_elimination_with_perm( py: Python, mut mat: PyReadwriteArray2, @@ -38,6 +46,13 @@ fn gauss_elimination_with_perm( /// Gauss elimination of a matrix mat with m rows and n columns. /// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] /// Returns the updated matrix mat. +/// Args: +/// mat: a boolean matrix with n rows and m columns +/// ncols: the number of columns for the gaussian elimination, +/// if ncols=None, then the elimination is done over all the colums +/// full_elim: whether to do a full elimination, or partial (upper triangular form) +/// Returns: +/// mat: the matrix after the gaussian elimination fn gauss_elimination( mut mat: PyReadwriteArray2, ncols: Option, @@ -51,6 +66,10 @@ fn gauss_elimination( #[pyo3(signature = (mat))] /// Given a boolean matrix mat after Gaussian elimination, computes its rank /// (i.e. simply the number of nonzero rows) +/// Args: +/// mat: a boolean matrix after gaussian elimination +/// Returns: +/// rank: the rank of the matrix fn compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyResult { let view = mat.as_array(); let rank = utils::compute_rank_after_gauss_elim_inner(view); @@ -60,6 +79,10 @@ fn compute_rank_after_gauss_elim(py: Python, mat: PyReadonlyArray2) -> PyR #[pyfunction] #[pyo3(signature = (mat))] /// Given a boolean matrix mat computes its rank +/// Args: +/// mat: a boolean matrix +/// Returns: +/// rank: the rank of the matrix fn compute_rank(py: Python, mat: PyReadonlyArray2) -> PyResult { let rank = utils::compute_rank_inner(mat.as_array()); Ok(rank.to_object(py)) @@ -89,6 +112,13 @@ pub fn calc_inverse_matrix( #[pyfunction] #[pyo3(signature = (mat1, mat2))] /// Binary matrix multiplication +/// Args: +/// mat1: a boolean matrix +/// mat2: a boolean matrix +/// Returns: +/// a boolean matrix which is the multiplication of mat1 and mat2 +/// Raises: +/// QiskitError: if the dimensions of mat1 and mat2 do not match pub fn binary_matmul( py: Python, mat1: PyReadonlyArray2, From af7ac2a570fcad3717709b9247d60612754dc8e0 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Tue, 25 Jun 2024 05:53:18 -0500 Subject: [PATCH 38/43] add release notes --- .../linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 releasenotes/notes/linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml diff --git a/releasenotes/notes/linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml b/releasenotes/notes/linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml new file mode 100644 index 000000000000..86205834607c --- /dev/null +++ b/releasenotes/notes/linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml @@ -0,0 +1,8 @@ +--- +features_synthesis: + - | + Port internal binary matrix utils from Python to rust, including + binary matrix multiplication, gaussian elimination, rank calculation, + binary matrix inversion, and random invertible binary matrix generation. + These functions are not part of the Qiskit API, and porting them to rust + improves the performance of certain synthesis methods. From 3424807dbca68b0ba2307c5d3481476ec33897a4 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 27 Jun 2024 06:27:19 -0500 Subject: [PATCH 39/43] remove unnecessary copy in python --- qiskit/quantum_info/operators/dihedral/random.py | 2 +- qiskit/synthesis/clifford/clifford_decompose_layers.py | 6 +++--- qiskit/synthesis/linear/linear_depth_lnn.py | 4 ++-- qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py | 1 - qiskit/synthesis/stabilizer/stabilizer_decompose.py | 2 +- qiskit/transpiler/passes/synthesis/high_level_synthesis.py | 4 ++-- 6 files changed, 9 insertions(+), 10 deletions(-) diff --git a/qiskit/quantum_info/operators/dihedral/random.py b/qiskit/quantum_info/operators/dihedral/random.py index 79220c128673..8223f87e9a17 100644 --- a/qiskit/quantum_info/operators/dihedral/random.py +++ b/qiskit/quantum_info/operators/dihedral/random.py @@ -54,7 +54,7 @@ def random_cnotdihedral(num_qubits, seed=None): ) seed = rng.integers(100000, size=1, dtype=np.uint64)[0] - linear = random_invertible_binary_matrix(num_qubits, seed=seed).astype(int) + linear = random_invertible_binary_matrix(num_qubits, seed=seed).astype(int, copy=False) elem.linear = linear # Random shift diff --git a/qiskit/synthesis/clifford/clifford_decompose_layers.py b/qiskit/synthesis/clifford/clifford_decompose_layers.py index 03b5ed216d9c..d2ac1808129e 100644 --- a/qiskit/synthesis/clifford/clifford_decompose_layers.py +++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py @@ -204,7 +204,7 @@ def _create_graph_state(cliff, validate=False): """ num_qubits = cliff.num_qubits - rank = compute_rank((cliff.stab_x).astype(bool)) + rank = compute_rank((cliff.stab_x).astype(bool, copy=False)) H1_circ = QuantumCircuit(num_qubits, name="H1") cliffh = cliff.copy() @@ -238,7 +238,7 @@ def _create_graph_state(cliff, validate=False): # validate that a layer of Hadamard gates and then appending cliff, provides a graph state. if validate: - stabh = (cliffh.stab_x).astype(bool) + stabh = (cliffh.stab_x).astype(bool, copy=False) if compute_rank(stabh) != num_qubits: raise QiskitError("The state is not a graph state.") @@ -269,7 +269,7 @@ def _decompose_graph_state(cliff, validate, cz_synth_func): """ num_qubits = cliff.num_qubits - rank = compute_rank((cliff.stab_x).astype(bool)) + rank = compute_rank(np.asarray(cliff.stab_x, dtype=bool)) cliff_cpy = cliff.copy() if rank < num_qubits: raise QiskitError("The stabilizer state is not a graph state.") diff --git a/qiskit/synthesis/linear/linear_depth_lnn.py b/qiskit/synthesis/linear/linear_depth_lnn.py index 2eab7d3e8bff..7c7360915e0f 100644 --- a/qiskit/synthesis/linear/linear_depth_lnn.py +++ b/qiskit/synthesis/linear/linear_depth_lnn.py @@ -222,8 +222,8 @@ def _optimize_cx_circ_depth_5n_line(mat): # According to [1] the synthesis is done on the inverse matrix # so the matrix mat is inverted at this step - mat_inv = mat.copy().astype(bool) - mat_cpy = calc_inverse_matrix(mat_inv).astype(bool) + mat_inv = mat.astype(bool, copy=True) + mat_cpy = calc_inverse_matrix(mat_inv) n = len(mat_cpy) diff --git a/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py index 93e365c60249..c0956ea3bc7c 100644 --- a/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py +++ b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py @@ -245,7 +245,6 @@ def synth_cx_cz_depth_line_my(mat_x: np.ndarray, mat_z: np.ndarray) -> QuantumCi n = len(mat_x) mat_x = calc_inverse_matrix(mat_x) - mat_x = mat_x.astype(bool) cx_instructions_rows_m2nw, cx_instructions_rows_nw2id = _optimize_cx_circ_depth_5n_line(mat_x) # Meanwhile, also build the -CZ- circuit via Phase gate insertions as per Algorithm 2 [2] diff --git a/qiskit/synthesis/stabilizer/stabilizer_decompose.py b/qiskit/synthesis/stabilizer/stabilizer_decompose.py index a8ae402c3935..ef324bc3cad1 100644 --- a/qiskit/synthesis/stabilizer/stabilizer_decompose.py +++ b/qiskit/synthesis/stabilizer/stabilizer_decompose.py @@ -143,7 +143,7 @@ def _calc_pauli_diff_stabilizer(cliff, cliff_target): phase.extend(phase_stab) phase = np.array(phase, dtype=int) - A = cliff.symplectic_matrix.astype(bool) + A = cliff.symplectic_matrix.astype(bool, copy=False) Ainv = calc_inverse_matrix(A) # By carefully writing how X, Y, Z gates affect each qubit, all we need to compute diff --git a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py index 8157d2f13337..bbc986621050 100644 --- a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py +++ b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py @@ -779,7 +779,7 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, ** use_inverted = options.get("use_inverted", False) use_transposed = options.get("use_transposed", False) - mat = high_level_object.linear.astype(bool) + mat = high_level_object.linear.astype(bool, copy=False) if use_transposed: mat = np.transpose(mat) @@ -831,7 +831,7 @@ def run(self, high_level_object, coupling_map=None, target=None, qubits=None, ** use_inverted = options.get("use_inverted", False) use_transposed = options.get("use_transposed", False) - mat = high_level_object.linear.astype(bool) + mat = high_level_object.linear.astype(bool, copy=False) if use_transposed: mat = np.transpose(mat) From 3327ddc39a1bde4f00a0cffaa1d7a3747807b145 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 27 Jun 2024 06:56:18 -0500 Subject: [PATCH 40/43] another copy fix --- qiskit/synthesis/clifford/clifford_decompose_layers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qiskit/synthesis/clifford/clifford_decompose_layers.py b/qiskit/synthesis/clifford/clifford_decompose_layers.py index d2ac1808129e..b7e18f02e91b 100644 --- a/qiskit/synthesis/clifford/clifford_decompose_layers.py +++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py @@ -204,7 +204,7 @@ def _create_graph_state(cliff, validate=False): """ num_qubits = cliff.num_qubits - rank = compute_rank((cliff.stab_x).astype(bool, copy=False)) + rank = compute_rank(np.asarray(cliff.stab_x, dtype=bool)) H1_circ = QuantumCircuit(num_qubits, name="H1") cliffh = cliff.copy() From 63a114e573cbe3ea2d7194a481b45285c65d0137 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 27 Jun 2024 07:05:29 -0500 Subject: [PATCH 41/43] another copy fix --- qiskit/synthesis/clifford/clifford_decompose_layers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qiskit/synthesis/clifford/clifford_decompose_layers.py b/qiskit/synthesis/clifford/clifford_decompose_layers.py index b7e18f02e91b..8b745823dc10 100644 --- a/qiskit/synthesis/clifford/clifford_decompose_layers.py +++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py @@ -210,7 +210,7 @@ def _create_graph_state(cliff, validate=False): if rank < num_qubits: stab = cliff.stab[:, :-1] - stab = stab.astype(bool) + stab = stab.astype(bool, copy=True) gauss_elimination(stab, num_qubits) Cmat = stab[rank:num_qubits, num_qubits:] From 102cc906987f572674c045f8126ca00134f4c005 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 27 Jun 2024 08:16:01 -0500 Subject: [PATCH 42/43] update rust docstrings --- crates/accelerate/src/synthesis/linear/mod.rs | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index f0ede103bd75..2fa158ea761f 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -20,15 +20,15 @@ mod utils; #[pyo3(signature = (mat, ncols=None, full_elim=false))] /// Gauss elimination of a matrix mat with m rows and n columns. /// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] -/// Returns the matrix mat, and the permutation perm that was done on the rows during the process. -/// perm[0 : rank] represents the indices of linearly independent rows in the original matrix. +/// Modifies the matrix mat in-place, and returns the permutation perm that was done +/// on the rows during the process. perm[0 : rank] represents the indices of linearly +/// independent rows in the original matrix. /// Args: /// mat: a boolean matrix with n rows and m columns /// ncols: the number of columns for the gaussian elimination, -/// if ncols=None, then the elimination is done over all the colums +/// if ncols=None, then the elimination is done over all the columns /// full_elim: whether to do a full elimination, or partial (upper triangular form) /// Returns: -/// mat: the matrix after the gaussian elimination /// perm: the permutation perm that was done on the rows during the process fn gauss_elimination_with_perm( py: Python, @@ -45,14 +45,12 @@ fn gauss_elimination_with_perm( #[pyo3(signature = (mat, ncols=None, full_elim=false))] /// Gauss elimination of a matrix mat with m rows and n columns. /// If full_elim = True, it allows full elimination of mat[:, 0 : ncols] -/// Returns the updated matrix mat. +/// This function modifies the input matrix in-place. /// Args: /// mat: a boolean matrix with n rows and m columns /// ncols: the number of columns for the gaussian elimination, -/// if ncols=None, then the elimination is done over all the colums +/// if ncols=None, then the elimination is done over all the columns /// full_elim: whether to do a full elimination, or partial (upper triangular form) -/// Returns: -/// mat: the matrix after the gaussian elimination fn gauss_elimination( mut mat: PyReadwriteArray2, ncols: Option, From c1fc5599f84fee930fbeabb3f429743054e5ed73 Mon Sep 17 00:00:00 2001 From: ShellyGarion Date: Thu, 27 Jun 2024 08:16:19 -0500 Subject: [PATCH 43/43] update release notes --- .../notes/linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/releasenotes/notes/linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml b/releasenotes/notes/linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml index 86205834607c..a8e9ec743808 100644 --- a/releasenotes/notes/linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml +++ b/releasenotes/notes/linear-binary-matrix-utils-rust-c48b5577749c34ab.yaml @@ -1,7 +1,7 @@ --- features_synthesis: - | - Port internal binary matrix utils from Python to rust, including + Port internal binary matrix utils from Python to Rust, including binary matrix multiplication, gaussian elimination, rank calculation, binary matrix inversion, and random invertible binary matrix generation. These functions are not part of the Qiskit API, and porting them to rust