diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs
new file mode 100644
index 000000000000..2fa158ea761f
--- /dev/null
+++ b/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(())
+}
diff --git a/crates/accelerate/src/synthesis/linear/utils.rs b/crates/accelerate/src/synthesis/linear/utils.rs
new file mode 100644
index 000000000000..b4dbf4993081
--- /dev/null
+++ b/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()
+}
diff --git a/crates/accelerate/src/synthesis/mod.rs b/crates/accelerate/src/synthesis/mod.rs
index f1a720459211..db28751437f6 100644
--- a/crates/accelerate/src/synthesis/mod.rs
+++ b/crates/accelerate/src/synthesis/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.
 
+mod linear;
 mod permutation;
 
 use pyo3::prelude::*;
@@ -18,5 +19,6 @@ use pyo3::wrap_pymodule;
 #[pymodule]
 pub fn synthesis(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pymodule!(permutation::permutation))?;
+    m.add_wrapped(wrap_pymodule!(linear::linear))?;
     Ok(())
 }
diff --git a/qiskit/__init__.py b/qiskit/__init__.py
index 5b8505654428..aca555da8cb8 100644
--- a/qiskit/__init__.py
+++ b/qiskit/__init__.py
@@ -81,6 +81,7 @@
 sys.modules["qiskit._accelerate.two_qubit_decompose"] = _accelerate.two_qubit_decompose
 sys.modules["qiskit._accelerate.vf2_layout"] = _accelerate.vf2_layout
 sys.modules["qiskit._accelerate.synthesis.permutation"] = _accelerate.synthesis.permutation
+sys.modules["qiskit._accelerate.synthesis.linear"] = _accelerate.synthesis.linear
 
 from qiskit.exceptions import QiskitError, MissingOptionalLibraryError
 
diff --git a/qiskit/quantum_info/operators/dihedral/random.py b/qiskit/quantum_info/operators/dihedral/random.py
index 4331d618d73d..8223f87e9a17 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, 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 21bea89b6575..8b745823dc10 100644
--- a/qiskit/synthesis/clifford/clifford_decompose_layers.py
+++ b/qiskit/synthesis/clifford/clifford_decompose_layers.py
@@ -33,9 +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,
+    compute_rank,
+    gauss_elimination,
+    gauss_elimination_with_perm,
+    binary_matmul,
 )
 
 
@@ -203,24 +204,25 @@ def _create_graph_state(cliff, validate=False):
     """
 
     num_qubits = cliff.num_qubits
-    rank = _compute_rank(cliff.stab_x)
+    rank = compute_rank(np.asarray(cliff.stab_x, dtype=bool))
     H1_circ = QuantumCircuit(num_qubits, name="H1")
     cliffh = cliff.copy()
 
     if rank < num_qubits:
         stab = cliff.stab[:, :-1]
-        stab = _gauss_elimination(stab, num_qubits)
+        stab = stab.astype(bool, copy=True)
+        gauss_elimination(stab, num_qubits)
 
         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
         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):
@@ -236,8 +238,8 @@ 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:
+            stabh = (cliffh.stab_x).astype(bool, copy=False)
+            if compute_rank(stabh) != num_qubits:
                 raise QiskitError("The state is not a graph state.")
 
     return H1_circ, cliffh
@@ -267,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(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.")
@@ -278,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:
@@ -340,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/__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 2811b755fa42..7c7360915e0f 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
 
 
@@ -222,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.astype(bool, copy=True)
     mat_cpy = calc_inverse_matrix(mat_inv)
 
     n = len(mat_cpy)
diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py
index 7a5b60641475..a76efdbb8d72 100644
--- a/qiskit/synthesis/linear/linear_matrix_utils.py
+++ b/qiskit/synthesis/linear/linear_matrix_utils.py
@@ -12,164 +12,16 @@
 
 """Utility functions for handling binary matrices."""
 
-from typing import Optional, Union
-import numpy as np
-from qiskit.exceptions import QiskitError
-
-
-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))
-        rank = _compute_rank(mat)
-    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]
-    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.
-
-    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 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, np.eye(n, dtype=int)), axis=1)
-    mat1 = _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.")
-
-    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_after_gauss_elim(mat):
-    """Given a 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)
-    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]
+# pylint: disable=unused-import
+from qiskit._accelerate.synthesis.linear import (
+    gauss_elimination,
+    gauss_elimination_with_perm,
+    compute_rank_after_gauss_elim,
+    compute_rank,
+    calc_inverse_matrix,
+    binary_matmul,
+    random_invertible_binary_matrix,
+    check_invertible_binary_matrix,
+    row_op,
+    col_op,
+)
diff --git a/qiskit/synthesis/stabilizer/stabilizer_decompose.py b/qiskit/synthesis/stabilizer/stabilizer_decompose.py
index ecdc1b3257ef..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(int)
+    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 668d10f32bcc..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(int)
+        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(int)
+        mat = high_level_object.linear.astype(bool, copy=False)
 
         if use_transposed:
             mat = np.transpose(mat)
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..a8e9ec743808
--- /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.
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)
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)))
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)
 
diff --git a/test/python/synthesis/test_linear_synthesis.py b/test/python/synthesis/test_linear_synthesis.py
index 98a49b6642f7..bbfab20a30fc 100644
--- a/test/python/synthesis/test_linear_synthesis.py
+++ b/test/python/synthesis/test_linear_synthesis.py
@@ -24,6 +24,7 @@
     random_invertible_binary_matrix,
     check_invertible_binary_matrix,
     calc_inverse_matrix,
+    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)
 
@@ -117,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
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):