diff --git a/crates/accelerate/src/synthesis/linear/lnn.rs b/crates/accelerate/src/synthesis/linear/lnn.rs
new file mode 100644
index 000000000000..edb1ae0913a2
--- /dev/null
+++ b/crates/accelerate/src/synthesis/linear/lnn.rs
@@ -0,0 +1,325 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2025
+//
+// 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::synthesis::linear::utils::{_col_op, _row_op, _row_sum, calc_inverse_matrix_inner};
+use ndarray::{Array1, Array2, ArrayView1, ArrayView2, ArrayViewMut2};
+use numpy::PyReadonlyArray2;
+use smallvec::smallvec;
+
+use pyo3::prelude::*;
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::operations::{Param, StandardGate};
+use qiskit_circuit::Qubit;
+
+// Optimize the synthesis of an n-qubit circuit contains only CX gates for
+// linear nearest neighbor (LNN) connectivity.
+// The depth of the circuit is bounded by 5*n, while the gate count is approximately 2.5*n^2
+//
+// References:
+// [1]: Kutin, S., Moulton, D. P., Smithline, L. (2007).
+// Computation at a Distance.
+// `arXiv:quant-ph/0701194 <https://arxiv.org/abs/quant-ph/0701194>`_.
+
+type InstructionList = Vec<(usize, usize)>;
+
+/// Add a cx gate to the instructions and update the matrix mat
+fn _row_op_update_instructions(
+    cx_instructions: &mut InstructionList,
+    mat: ArrayViewMut2<bool>,
+    a: usize,
+    b: usize,
+) {
+    cx_instructions.push((a, b));
+    _row_op(mat, a, b);
+}
+
+/// Get the instructions for a lower triangular basis change of a matrix mat.
+/// See the proof of Proposition 7.3 in [1].
+/// mat_inv needs to be the inverted matrix of mat
+/// The outputs are the permuted versions of mat and mat_inv
+fn _get_lower_triangular<'a>(
+    n: usize,
+    mat: ArrayView2<bool>,
+    mut mat_inv: ArrayViewMut2<'a, bool>,
+) -> (Array2<bool>, ArrayViewMut2<'a, bool>) {
+    let mut mat = mat.to_owned();
+    let mut mat_t = mat.to_owned();
+
+    let mut cx_instructions_rows: InstructionList = Vec::new();
+    // Use the instructions in U, which contains only gates of the form cx(a,b) a>b
+    // to transform the matrix to a permuted lower-triangular matrix.
+    // The original Matrix mat is unchanged, but mat_inv is
+
+    for i in (0..n).rev() {
+        // Find the last "1" in row i, use COL operations to the left in order to
+        // zero out all other "1"s in that row.
+        let cols_to_update: Vec<usize> = (0..n).rev().filter(|&j| mat[[i, j]]).collect();
+        let (first_j, cols_to_update) = cols_to_update.split_first().unwrap();
+        cols_to_update.iter().for_each(|j| {
+            _col_op(mat.view_mut(), *first_j, *j);
+        });
+
+        // Use row operations directed upwards to zero out all "1"s above the remaining "1" in row i
+        let rows_to_update: Vec<usize> = (0..i).rev().filter(|k| mat[[*k, *first_j]]).collect();
+        rows_to_update.into_iter().for_each(|k| {
+            _row_op_update_instructions(&mut cx_instructions_rows, mat.view_mut(), i, k);
+        });
+    }
+    // Apply only U instructions to get the permuted L
+    for (ctrl, trgt) in cx_instructions_rows {
+        _row_op(mat_t.view_mut(), ctrl, trgt);
+        _col_op(mat_inv.view_mut(), trgt, ctrl); // performs an inverted col_op
+    }
+    (mat_t, mat_inv)
+}
+
+/// For each row in mat_t, save the column index of the last "1"
+fn _get_label_arr(n: usize, mat_t: ArrayView2<bool>) -> Vec<usize> {
+    (0..n)
+        .map(|i| (0..n).find(|&j| mat_t[[i, n - 1 - j]]).unwrap_or(n))
+        .collect()
+}
+
+/// Check if "row" is a linear combination of all rows in mat_inv_t not including the row labeled by k
+fn _in_linear_combination(
+    label_arr_t: &[usize],
+    mat_inv_t: ArrayView2<bool>,
+    row: ArrayView1<bool>,
+    k: usize,
+) -> bool {
+    // Find the linear combination of mat_t rows which produces "row"
+    !(0..row.len())
+        .filter(|&row_l| row[row_l])
+        .fold(Array1::from_elem(row.len(), false), |w_needed, row_l| {
+            _row_sum(w_needed.view(), mat_inv_t.row(row_l)).unwrap()
+        })[label_arr_t[k]]
+}
+
+/// Returns label_arr_t = label_arr^(-1)
+fn _get_label_arr_t(n: usize, label_arr: &[usize]) -> Vec<usize> {
+    let mut label_arr_t: Vec<usize> = vec![0; n];
+    (0..n).for_each(|i| label_arr_t[label_arr[i]] = i);
+    label_arr_t
+}
+
+/// Transform an arbitrary boolean invertible matrix to a north-west triangular matrix
+/// by Proposition 7.3 in [1]
+fn _matrix_to_north_west(
+    n: usize,
+    mut mat: ArrayViewMut2<bool>,
+    mut mat_inv: ArrayViewMut2<bool>,
+) -> InstructionList {
+    // The rows of mat_t hold all w_j vectors (see [1]). mat_inv_t is the inverted matrix of mat_t
+    // To save time on needless copying, we change mat_inv into mat_inv_t, since we won't need mat_inv anymore
+    let (mat_t, mat_inv_t) = _get_lower_triangular(n, mat.view(), mat_inv.view_mut());
+    // Get all pi(i) labels
+    let mut label_arr = _get_label_arr(n, mat_t.view());
+
+    // Save the original labels, exchange index <-> value
+    let label_arr_t = _get_label_arr_t(n, &label_arr);
+    let mut first_qubit = 0;
+    let mut empty_layers = 0;
+    let mut done = false;
+    let mut cx_instructions_rows: InstructionList = Vec::new();
+    while !done {
+        // At each iteration the values of i switch between even and odd
+        let mut at_least_one_needed = false;
+        for i in (first_qubit..n - 1).step_by(2) {
+            // "If j < k, we do nothing" (see [1])
+            // "If j > k, we swap the two labels, and we also perform a box" (see [1])
+            if label_arr[i] > label_arr[i + 1] {
+                at_least_one_needed = true;
+                // iterate on column indices, output rows as Vec<bool>
+                let row_sum = _row_sum(mat.row(i), mat.row(i + 1)).unwrap();
+                // "Let W be the span of all w_l for l!=k" (see [1])
+                // " We can perform a box on <i> and <i + 1> that writes a vector in W to wire <i + 1>."
+                // (see [1])
+                if _in_linear_combination(
+                    &label_arr_t,
+                    mat_inv_t.view(),
+                    mat.row(i + 1),
+                    label_arr[i + 1],
+                ) {
+                    // do nothing
+                } else if _in_linear_combination(
+                    &label_arr_t,
+                    mat_inv_t.view(),
+                    row_sum.view(),
+                    label_arr[i + 1],
+                ) {
+                    _row_op_update_instructions(
+                        &mut cx_instructions_rows,
+                        mat.view_mut(),
+                        i,
+                        i + 1,
+                    );
+                } else if _in_linear_combination(
+                    &label_arr_t,
+                    mat_inv_t.view(),
+                    mat.row(i),
+                    label_arr[i + 1],
+                ) {
+                    _row_op_update_instructions(
+                        &mut cx_instructions_rows,
+                        mat.view_mut(),
+                        i + 1,
+                        i,
+                    );
+                    _row_op_update_instructions(
+                        &mut cx_instructions_rows,
+                        mat.view_mut(),
+                        i,
+                        i + 1,
+                    );
+                }
+                (label_arr[i], label_arr[i + 1]) = (label_arr[i + 1], label_arr[i]);
+            }
+        }
+        if !at_least_one_needed {
+            empty_layers += 1;
+            if empty_layers > 1 {
+                // if nothing happened twice in a row, then finished.
+                done = true;
+            }
+        } else {
+            empty_layers = 0;
+        }
+        first_qubit = 1 - first_qubit;
+    }
+    cx_instructions_rows
+}
+
+/// Transform a north-west triangular matrix to identity in depth 3*n by Proposition 7.4 of [1]
+fn _north_west_to_identity(n: usize, mut mat: ArrayViewMut2<bool>) -> InstructionList {
+    // At start the labels are in reversed order
+    let mut label_arr: Vec<usize> = (0..n).rev().collect();
+    let mut first_qubit = 0;
+    let mut empty_layers = 0;
+    let mut done = false;
+    let mut cx_instructions_rows: InstructionList = Vec::new();
+    while !done {
+        let mut at_least_one_needed = false;
+        for i in (first_qubit..n - 1).step_by(2) {
+            // Exchange the labels if needed
+            if label_arr[i] > label_arr[i + 1] {
+                at_least_one_needed = true;
+                // If row i has "1" in column i+1, swap and remove the "1" (in depth 2)
+                // otherwise, only do a swap (in depth 3)
+                if !mat[[i, label_arr[i + 1]]] {
+                    // Adding this turns the operation to a SWAP
+                    _row_op_update_instructions(
+                        &mut cx_instructions_rows,
+                        mat.view_mut(),
+                        i + 1,
+                        i,
+                    );
+                }
+                _row_op_update_instructions(&mut cx_instructions_rows, mat.view_mut(), i, i + 1);
+                _row_op_update_instructions(&mut cx_instructions_rows, mat.view_mut(), i + 1, i);
+
+                (label_arr[i], label_arr[i + 1]) = (label_arr[i + 1], label_arr[i]);
+            }
+        }
+
+        if !at_least_one_needed {
+            empty_layers += 1;
+            if empty_layers > 1 {
+                // if nothing happened twice in a row, then finished.
+                done = true;
+            }
+        } else {
+            empty_layers = 0;
+        }
+        first_qubit = 1 - first_qubit;
+    }
+    cx_instructions_rows
+}
+
+/// Find instruction to synthesize CX circuit in depth bounded by 5n for LNN connectivity.
+/// The algorithm [1] has two steps:
+/// a) transform the original matrix to a north-west matrix (m2nw),
+/// b) transform the north-west matrix to identity (nw2id).
+///
+/// A square n-by-n matrix A is called north-west if A[i][j]=0 for all i+j>=n
+/// For example, the following matrix is north-west:
+/// [[0, 1, 0, 1]
+/// [1, 1, 1, 0]
+/// [0, 1, 0, 0]
+/// [1, 0, 0, 0]]
+
+/// According to [1] the synthesis is done on the inverse matrix
+/// so the matrix mat is inverted at this step
+
+/// References:
+/// [1]: Kutin, S., Moulton, D. P., Smithline, L. (2007).
+/// Computation at a Distance.
+/// `arXiv:quant-ph/0701194 <https://arxiv.org/abs/quant-ph/0701194>`_.
+fn _synth_cnot_lnn_instructions(arrayview: ArrayView2<bool>) -> (InstructionList, InstructionList) {
+    // According to [1] the synthesis is done on the inverse matrix
+    // so the matrix mat is inverted at this step
+    let mut mat_inv: Array2<bool> = arrayview.to_owned();
+    let mut mat_cpy = calc_inverse_matrix_inner(mat_inv.view(), false).unwrap();
+
+    let n = mat_cpy.nrows();
+
+    // Transform an arbitrary invertible matrix to a north-west triangular matrix
+    // by Proposition 7.3 of [1]
+
+    let cx_instructions_rows_m2nw =
+        _matrix_to_north_west(n, mat_cpy.view_mut(), mat_inv.view_mut());
+    // Transform a north-west triangular matrix to identity in depth 3*n
+    // by Proposition 7.4 of [1]
+
+    let cx_instructions_rows_nw2id = _north_west_to_identity(n, mat_cpy.view_mut());
+
+    (cx_instructions_rows_m2nw, cx_instructions_rows_nw2id)
+}
+
+/// Find instruction to synthesize CX circuit in depth bounded by 5n for LNN connectivity.
+/// Uses the algorithm by Kutin, Moulton, Smithline
+/// described in `arXiv:quant-ph/0701194 <https://arxiv.org/abs/quant-ph/0701194>`_.
+/// Returns: Tuple with two lists of instructions for CX gates
+/// Corresponding to the two parts of the algorithm
+#[pyfunction]
+#[pyo3(signature = (mat))]
+pub fn py_synth_cnot_lnn_instructions(
+    mat: PyReadonlyArray2<bool>,
+) -> PyResult<(InstructionList, InstructionList)> {
+    Ok(_synth_cnot_lnn_instructions(mat.as_array()))
+}
+
+/// Synthesize CX circuit in depth bounded by 5n for LNN connectivity.
+/// Uses the algorithm by Kutin, Moulton, Smithline
+/// described in `arXiv:quant-ph/0701194 <https://arxiv.org/abs/quant-ph/0701194>`_.
+/// Returns: The CircuitData of the synthesized circuit.
+#[pyfunction]
+#[pyo3(signature = (mat))]
+pub fn py_synth_cnot_depth_line_kms(
+    py: Python,
+    mat: PyReadonlyArray2<bool>,
+) -> PyResult<CircuitData> {
+    let num_qubits = mat.as_array().nrows(); // is a quadratic matrix
+    let (cx_instructions_rows_m2nw, cx_instructions_rows_nw2id) =
+        _synth_cnot_lnn_instructions(mat.as_array());
+
+    let instructions = cx_instructions_rows_m2nw
+        .into_iter()
+        .chain(cx_instructions_rows_nw2id)
+        .map(|(ctrl, target)| {
+            (
+                StandardGate::CXGate,
+                smallvec![],
+                smallvec![Qubit(ctrl as u32), Qubit(target as u32)],
+            )
+        });
+    CircuitData::from_standard_gates(py, num_qubits as u32, instructions, Param::Float(0.0))
+}
diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs
index a6406adabb48..688951864109 100644
--- a/crates/accelerate/src/synthesis/linear/mod.rs
+++ b/crates/accelerate/src/synthesis/linear/mod.rs
@@ -15,6 +15,7 @@ use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2, PyReadwriteArray2};
 use pyo3::prelude::*;
 use pyo3::IntoPyObjectExt;
 
+mod lnn;
 mod pmh;
 pub mod utils;
 
@@ -187,5 +188,7 @@ pub fn linear(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(random_invertible_binary_matrix))?;
     m.add_wrapped(wrap_pyfunction!(check_invertible_binary_matrix))?;
     m.add_wrapped(wrap_pyfunction!(pmh::synth_cnot_count_full_pmh))?;
+    m.add_wrapped(wrap_pyfunction!(lnn::py_synth_cnot_depth_line_kms))?;
+    m.add_wrapped(wrap_pyfunction!(lnn::py_synth_cnot_lnn_instructions))?;
     Ok(())
 }
diff --git a/crates/accelerate/src/synthesis/linear/utils.rs b/crates/accelerate/src/synthesis/linear/utils.rs
index 5f1fb99fa682..c620acb5eec0 100644
--- a/crates/accelerate/src/synthesis/linear/utils.rs
+++ b/crates/accelerate/src/synthesis/linear/utils.rs
@@ -10,7 +10,9 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-use ndarray::{azip, concatenate, s, Array2, ArrayView1, ArrayView2, ArrayViewMut2, Axis, Zip};
+use ndarray::{
+    azip, concatenate, s, Array1, Array2, ArrayView1, ArrayView2, ArrayViewMut2, Axis, Zip,
+};
 use rand::{Rng, SeedableRng};
 use rand_pcg::Pcg64Mcg;
 use rayon::iter::{IndexedParallelIterator, ParallelIterator};
@@ -194,6 +196,29 @@ pub fn _add_row_or_col(mut mat: ArrayViewMut2<bool>, add_cols: &bool, ctrl: usiz
     row1.zip_mut_with(&row0, |x, &y| *x ^= y);
 }
 
+/// Perform ROW operation on a matrix mat
+pub fn _row_op(mat: ArrayViewMut2<bool>, ctrl: usize, trgt: usize) {
+    _add_row_or_col(mat, &false, ctrl, trgt);
+}
+
+/// Perform COL operation on a matrix mat
+pub fn _col_op(mat: ArrayViewMut2<bool>, ctrl: usize, trgt: usize) {
+    _add_row_or_col(mat, &true, ctrl, trgt);
+}
+
+/// Returns the element-wise sum of two boolean rows (i.e. addition modulo 2)
+pub fn _row_sum(row_1: ArrayView1<bool>, row_2: ArrayView1<bool>) -> Result<Array1<bool>, String> {
+    if row_1.len() != row_2.len() {
+        Err(format!(
+            "row sum performed on rows with different lengths ({} and {})",
+            row_1.len(),
+            row_2.len()
+        ))
+    } else {
+        Ok((0..row_1.len()).map(|i| row_1[i] ^ row_2[i]).collect())
+    }
+}
+
 /// 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 {
diff --git a/qiskit/synthesis/linear/linear_depth_lnn.py b/qiskit/synthesis/linear/linear_depth_lnn.py
index 7c7360915e0f..fdbabd68f025 100644
--- a/qiskit/synthesis/linear/linear_depth_lnn.py
+++ b/qiskit/synthesis/linear/linear_depth_lnn.py
@@ -25,217 +25,8 @@
 import numpy as np
 from qiskit.exceptions import QiskitError
 from qiskit.circuit import QuantumCircuit
-from qiskit.synthesis.linear.linear_matrix_utils import (
-    calc_inverse_matrix,
-    check_invertible_binary_matrix,
-    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)
-
-
-def _get_lower_triangular(n, mat, mat_inv):
-    # Get the instructions for a lower triangular basis change of a matrix mat.
-    # See the proof of Proposition 7.3 in [1].
-    mat = mat.copy()
-    mat_t = mat.copy()
-    mat_inv_t = mat_inv.copy()
-
-    cx_instructions_rows = []
-
-    # Use the instructions in U, which contains only gates of the form cx(a,b) a>b
-    # to transform the matrix to a permuted lower-triangular matrix.
-    # The original Matrix is unchanged.
-    for i in reversed(range(0, n)):
-        found_first = False
-        # Find the last "1" in row i, use COL operations to the left in order to
-        # zero out all other "1"s in that row.
-        for j in reversed(range(0, n)):
-            if mat[i, j]:
-                if not found_first:
-                    found_first = True
-                    first_j = j
-                else:
-                    # cx_instructions_cols (L instructions) are not needed
-                    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]:
-                _row_op_update_instructions(cx_instructions_rows, mat, i, k)
-
-    # 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])
-    return mat_t, mat_inv_t
-
-
-def _get_label_arr(n, mat_t):
-    # For each row in mat_t, save the column index of the last "1"
-    label_arr = []
-    for i in range(n):
-        j = 0
-        while not mat_t[i, n - 1 - j]:
-            j += 1
-        label_arr.append(j)
-    return label_arr
-
-
-def _in_linear_combination(label_arr_t, mat_inv_t, row, k):
-    # Check if "row" is a linear combination of all rows in mat_inv_t not including the row labeled by k
-    indx_k = label_arr_t[k]
-    w_needed = np.zeros(len(row), dtype=bool)
-    # Find the linear combination of mat_t rows which produces "row"
-    for row_l, _ in enumerate(row):
-        if row[row_l]:
-            # mat_inv_t can be thought of as a set of instructions. Row l in mat_inv_t
-            # indicates which rows from mat_t are necessary to produce the elementary vector e_l
-            w_needed = w_needed ^ mat_inv_t[row_l]
-    # If the linear combination requires the row labeled by k
-    if w_needed[indx_k]:
-        return False
-    return True
-
-
-def _get_label_arr_t(n, label_arr):
-    # Returns label_arr_t = label_arr^(-1)
-    label_arr_t = [None] * n
-    for i in range(n):
-        label_arr_t[label_arr[i]] = i
-    return label_arr_t
-
-
-def _matrix_to_north_west(n, mat, mat_inv):
-    # Transform an arbitrary boolean invertible matrix to a north-west triangular matrix
-    # by Proposition 7.3 in [1]
-
-    # The rows of mat_t hold all w_j vectors (see [1]). mat_inv_t is the inverted matrix of mat_t
-    mat_t, mat_inv_t = _get_lower_triangular(n, mat, mat_inv)
-
-    # Get all pi(i) labels
-    label_arr = _get_label_arr(n, mat_t)
-
-    # Save the original labels, exchange index <-> value
-    label_arr_t = _get_label_arr_t(n, label_arr)
-
-    first_qubit = 0
-    empty_layers = 0
-    done = False
-    cx_instructions_rows = []
-
-    while not done:
-        # At each iteration the values of i switch between even and odd
-        at_least_one_needed = False
-
-        for i in range(first_qubit, n - 1, 2):
-            # "If j < k, we do nothing" (see [1])
-            # "If j > k, we swap the two labels, and we also perform a box" (see [1])
-            if label_arr[i] > label_arr[i + 1]:
-                at_least_one_needed = True
-                # "Let W be the span of all w_l for l!=k" (see [1])
-                # " We can perform a box on <i> and <i + 1> that writes a vector in W to wire <i + 1>."
-                # (see [1])
-                if _in_linear_combination(label_arr_t, mat_inv_t, mat[i + 1], label_arr[i + 1]):
-                    pass
-
-                elif _in_linear_combination(
-                    label_arr_t, mat_inv_t, mat[i + 1] ^ mat[i], label_arr[i + 1]
-                ):
-                    _row_op_update_instructions(cx_instructions_rows, mat, i, i + 1)
-
-                elif _in_linear_combination(label_arr_t, mat_inv_t, mat[i], label_arr[i + 1]):
-                    _row_op_update_instructions(cx_instructions_rows, mat, i + 1, i)
-                    _row_op_update_instructions(cx_instructions_rows, mat, i, i + 1)
-
-                label_arr[i], label_arr[i + 1] = label_arr[i + 1], label_arr[i]
-
-        if not at_least_one_needed:
-            empty_layers += 1
-            if empty_layers > 1:  # if nothing happened twice in a row, then finished.
-                done = True
-        else:
-            empty_layers = 0
-
-        first_qubit = int(not first_qubit)
-
-    return cx_instructions_rows
-
-
-def _north_west_to_identity(n, mat):
-    # Transform a north-west triangular matrix to identity in depth 3*n by Proposition 7.4 of [1]
-
-    # At start the labels are in reversed order
-    label_arr = list(reversed(range(n)))
-    first_qubit = 0
-    empty_layers = 0
-    done = False
-    cx_instructions_rows = []
-
-    while not done:
-        at_least_one_needed = False
-
-        for i in range(first_qubit, n - 1, 2):
-            # Exchange the labels if needed
-            if label_arr[i] > label_arr[i + 1]:
-                at_least_one_needed = True
-
-                # If row i has "1" in column i+1, swap and remove the "1" (in depth 2)
-                # otherwise, only do a swap (in depth 3)
-                if not mat[i, label_arr[i + 1]]:
-                    # Adding this turns the operation to a SWAP
-                    _row_op_update_instructions(cx_instructions_rows, mat, i + 1, i)
-
-                _row_op_update_instructions(cx_instructions_rows, mat, i, i + 1)
-                _row_op_update_instructions(cx_instructions_rows, mat, i + 1, i)
-
-                label_arr[i], label_arr[i + 1] = label_arr[i + 1], label_arr[i]
-
-        if not at_least_one_needed:
-            empty_layers += 1
-            if empty_layers > 1:  # if nothing happened twice in a row, then finished.
-                done = True
-        else:
-            empty_layers = 0
-
-        first_qubit = int(not first_qubit)
-
-    return cx_instructions_rows
-
-
-def _optimize_cx_circ_depth_5n_line(mat):
-    # Optimize CX circuit in depth bounded by 5n for LNN connectivity.
-    # The algorithm [1] has two steps:
-    # a) transform the original matrix to a north-west matrix (m2nw),
-    # b) transform the north-west matrix to identity (nw2id).
-    #
-    # A square n-by-n matrix A is called north-west if A[i][j]=0 for all i+j>=n
-    # For example, the following matrix is north-west:
-    # [[0, 1, 0, 1]
-    #  [1, 1, 1, 0]
-    #  [0, 1, 0, 0]
-    #  [1, 0, 0, 0]]
-
-    # According to [1] the synthesis is done on the inverse matrix
-    # so the matrix mat is inverted at this step
-    mat_inv = mat.astype(bool, copy=True)
-    mat_cpy = calc_inverse_matrix(mat_inv)
-
-    n = len(mat_cpy)
-
-    # Transform an arbitrary invertible matrix to a north-west triangular matrix
-    # by Proposition 7.3 of [1]
-    cx_instructions_rows_m2nw = _matrix_to_north_west(n, mat_cpy, mat_inv)
-
-    # Transform a north-west triangular matrix to identity in depth 3*n
-    # by Proposition 7.4 of [1]
-    cx_instructions_rows_nw2id = _north_west_to_identity(n, mat_cpy)
-
-    return cx_instructions_rows_m2nw, cx_instructions_rows_nw2id
+from qiskit.synthesis.linear.linear_matrix_utils import check_invertible_binary_matrix
+from qiskit._accelerate.synthesis.linear import py_synth_cnot_depth_line_kms as fast_kms
 
 
 def synth_cnot_depth_line_kms(mat: np.ndarray[bool]) -> QuantumCircuit:
@@ -264,13 +55,7 @@ def synth_cnot_depth_line_kms(mat: np.ndarray[bool]) -> QuantumCircuit:
     if not check_invertible_binary_matrix(mat):
         raise QiskitError("The input matrix is not invertible.")
 
-    # Returns the quantum circuit constructed from the instructions
-    # that we got in _optimize_cx_circ_depth_5n_line
-    num_qubits = len(mat)
-    cx_inst = _optimize_cx_circ_depth_5n_line(mat)
-    qc = QuantumCircuit(num_qubits)
-    for pair in cx_inst[0]:
-        qc.cx(pair[0], pair[1])
-    for pair in cx_inst[1]:
-        qc.cx(pair[0], pair[1])
-    return qc
+    circuit_data = fast_kms(mat)
+
+    # construct circuit from the data
+    return QuantumCircuit._from_circuit_data(circuit_data, add_regs=True)
diff --git a/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py
index c0956ea3bc7c..9773d1d63f84 100644
--- a/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py
+++ b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py
@@ -33,7 +33,7 @@
 
 from qiskit.circuit import QuantumCircuit
 from qiskit.synthesis.linear.linear_matrix_utils import calc_inverse_matrix
-from qiskit.synthesis.linear.linear_depth_lnn import _optimize_cx_circ_depth_5n_line
+from qiskit._accelerate.synthesis.linear import py_synth_cnot_lnn_instructions
 
 
 def _initialize_phase_schedule(mat_z):
@@ -245,7 +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)
 
-    cx_instructions_rows_m2nw, cx_instructions_rows_nw2id = _optimize_cx_circ_depth_5n_line(mat_x)
+    cx_instructions_rows_m2nw, cx_instructions_rows_nw2id = py_synth_cnot_lnn_instructions(mat_x)
 
     # Meanwhile, also build the -CZ- circuit via Phase gate insertions as per Algorithm 2 [2]
     phase_schedule = _initialize_phase_schedule(mat_z)
diff --git a/releasenotes/notes/rust-linear-depth-lnn-69532c9c6b86ffe8.yaml b/releasenotes/notes/rust-linear-depth-lnn-69532c9c6b86ffe8.yaml
new file mode 100644
index 000000000000..24f7f0b710e2
--- /dev/null
+++ b/releasenotes/notes/rust-linear-depth-lnn-69532c9c6b86ffe8.yaml
@@ -0,0 +1,5 @@
+---
+features_transpiler:
+  - |
+    The :meth:`.synth_cnot_depth_line_kms` pass was ported into rust, with preliminary
+    benchmarks pointing at a factor of 20x speedup.
\ No newline at end of file