binary matrices utils in rust (#12456)

* gaussian elimination in rust

* handle lint errors

* replace python function by rust function for gauss elimination

* change matrix elements type from bool to i8

* add parallelization in row operations

* update matrices in place

* change matrix type in rust code to bool

* handle type in python code

* update filter following review

* remove parallelization using rayon

* move _gauss_elimination_with_perm to rust

* fix fmt error

* simplify _gauss_elimination function

* update _compute_rank_after_gauss_elim to rust

* update _row_op and _col_op

* transfer _row_op and _col_op from python to rust

* fix code due to failing tests

* minor update of types

* move calc_inverse_matrix to rust, add _binary_matmul in rust

* fix failing tests, by changing mat type from int to bool

* update rust docstrings

* add function _add_row_or_col to rust code

* improve binary_matmul

* proper error handling

* unified format of function names

* move compute_rank from python to rust, update errors

* update type of mat in compute_rank

* move random_invertible_binary_matrix and check_invertible_binary_matrix to rust

* Updating HighLevelSynthesis tests that depend on the specific random number

* Updating LinearSynthesis tests to pass seeds

* Updating tests in test_linear_function

* changing the matrix type in random dyhedral to be a matrix of ints rather than bools

* updating cx_cz synthesis tests

* updating clifford tests

* remove unused imports

* add option seed=None

* enhance rust docs

* add release notes

* remove unnecessary copy in python

* another copy fix

* another copy fix

* update rust docstrings

* update release notes

---------

Co-authored-by: AlexanderIvrii <[email protected]>

crates/accelerate/src/synthesis/linear/mod.rs

@@ -0,0 +1,190 @@
+// 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 crate::QiskitError;
+use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2, PyReadwriteArray2};
+use pyo3::prelude::*;
+
+mod utils;
+
+#[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]
+/// 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 columns
+///     full_elim: whether to do a full elimination, or partial (upper triangular form)
+/// Returns:
+///     perm: the permutation perm that was done on the rows during the process
+fn gauss_elimination_with_perm(
+    py: Python,
+    mut mat: PyReadwriteArray2<bool>,
+    ncols: Option<usize>,
+    full_elim: Option<bool>,
+) -> PyResult<PyObject> {
+    let matmut = mat.as_array_mut();
+    let perm = utils::gauss_elimination_with_perm_inner(matmut, ncols, full_elim);
+    Ok(perm.to_object(py))
+}
+
+#[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]
+/// 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 columns
+///     full_elim: whether to do a full elimination, or partial (upper triangular form)
+fn gauss_elimination(
+    mut mat: PyReadwriteArray2<bool>,
+    ncols: Option<usize>,
+    full_elim: Option<bool>,
+) {
+    let matmut = mat.as_array_mut();
+    let _perm = utils::gauss_elimination_with_perm_inner(matmut, ncols, full_elim);
+}
+
+#[pyfunction]
+#[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<bool>) -> PyResult<PyObject> {
+    let view = mat.as_array();
+    let rank = utils::compute_rank_after_gauss_elim_inner(view);
+    Ok(rank.to_object(py))
+}
+
+#[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<bool>) -> PyResult<PyObject> {
+    let rank = utils::compute_rank_inner(mat.as_array());
+    Ok(rank.to_object(py))
+}
+
+#[pyfunction]
+#[pyo3(signature = (mat, verify=false))]
+/// 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 or not invertible.
+pub fn calc_inverse_matrix(
+    py: Python,
+    mat: PyReadonlyArray2<bool>,
+    verify: Option<bool>,
+) -> PyResult<Py<PyArray2<bool>>> {
+    let view = mat.as_array();
+    let invmat =
+        utils::calc_inverse_matrix_inner(view, verify.is_some()).map_err(QiskitError::new_err)?;
+    Ok(invmat.into_pyarray_bound(py).unbind())
+}
+
+#[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<bool>,
+    mat2: PyReadonlyArray2<bool>,
+) -> PyResult<Py<PyArray2<bool>>> {
+    let view1 = mat1.as_array();
+    let view2 = mat2.as_array();
+    let result = utils::binary_matmul_inner(view1, view2).map_err(QiskitError::new_err)?;
+    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<bool>, ctrl: usize, trgt: usize) {
+    let matmut = mat.as_array_mut();
+    utils::_add_row_or_col(matmut, &false, ctrl, trgt)
+}
+
+#[pyfunction]
+#[pyo3(signature = (mat, ctrl, trgt))]
+/// Perform COL operation on a matrix mat (in the inverse direction)
+fn col_op(mut mat: PyReadwriteArray2<bool>, ctrl: usize, trgt: usize) {
+    let matmut = mat.as_array_mut();
+    utils::_add_row_or_col(matmut, &true, trgt, ctrl)
+}
+
+#[pyfunction]
+#[pyo3(signature = (num_qubits, seed=None))]
+/// 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<u64>,
+) -> PyResult<Py<PyArray2<bool>>> {
+    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<bool>) -> PyResult<PyObject> {
+    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<PyModule>) -> 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))?;
+    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(())
+}

crates/accelerate/src/synthesis/linear/utils.rs

@@ -0,0 +1,200 @@
+// 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::{concatenate, s, Array2, ArrayView2, ArrayViewMut2, Axis};
+use rand::{Rng, SeedableRng};
+use rand_pcg::Pcg64Mcg;
+
+/// Binary matrix multiplication
+pub fn binary_matmul_inner(
+    mat1: ArrayView2<bool>,
+    mat2: ArrayView2<bool>,
+) -> Result<Array2<bool>, 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 {
+        return Err(format!(
+            "Cannot multiply matrices with inappropriate dimensions {}, {}",
+            n1_cols, n2_rows
+        ));
+    }
+
+    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.
+/// 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.
+pub fn gauss_elimination_with_perm_inner(
+    mut mat: ArrayViewMut2<bool>,
+    ncols: Option<usize>,
+    full_elim: Option<bool>,
+) -> Vec<usize> {
+    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<usize> = Vec::from_iter(0..m);
+
+    let mut r = 0; // current rank
+    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;
+                    new_k = j;
+                    new_r = i;
+                    break;
+                }
+            }
+            if is_non_zero {
+                break;
+            }
+        }
+        if !is_non_zero {
+            return perm; // 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);
+            perm.swap(r, new_r);
+        }
+
+        // Copy source row to avoid trying multiple borrows at once
+        let row0 = mat.row(r).to_owned();
+        mat.axis_iter_mut(Axis(0))
+            .enumerate()
+            .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);
+            });
+
+        r += 1;
+    }
+    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_inner(mat: ArrayView2<bool>) -> usize {
+    let rank: usize = mat
+        .axis_iter(Axis(0))
+        .map(|row| row.fold(false, |out, val| out | *val) as usize)
+        .sum();
+    rank
+}
+
+/// Given a boolean matrix mat computes its rank
+pub fn compute_rank_inner(mat: ArrayView2<bool>) -> 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
+}
+
+/// Given a square boolean matrix mat, tries to compute its inverse.
+pub fn calc_inverse_matrix_inner(
+    mat: ArrayView2<bool>,
+    verify: bool,
+) -> Result<Array2<bool>, String> {
+    if mat.shape()[0] != mat.shape()[1] {
+        return Err("Matrix to invert is a non-square matrix.".to_string());
+    }
+    let n = mat.shape()[0];
+
+    // concatenate the matrix and identity
+    let identity_matrix: Array2<bool> = 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_inner(mat1.view_mut(), None, Some(true));
+
+    let r = compute_rank_after_gauss_elim_inner(mat1.slice(s![.., 0..n]));
+    if r < n {
+        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_inner(mat, (&invmat).into())?;
+        let identity_matrix: Array2<bool> = Array2::from_shape_fn((n, n), |(i, j)| i == j);
+        if mat2.ne(&identity_matrix) {
+            return Err("The inverse matrix is not correct.".to_string());
+        }
+    }
+
+    Ok(invmat)
+}
+
+/// 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.
+pub fn _add_row_or_col(mut mat: ArrayViewMut2<bool>, 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);
+}
+
+/// Generate a random invertible n x n binary matrix.
+pub fn random_invertible_binary_matrix_inner(num_qubits: usize, seed: Option<u64>) -> Array2<bool> {
+    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>) -> bool {
+    if mat.nrows() != mat.ncols() {
+        return false;
+    }
+    let rank = compute_rank_inner(mat);
+    rank == mat.nrows()
+}