From 708c8fdedb78dddee0dafd20cf27b185f89df790 Mon Sep 17 00:00:00 2001 From: Gadi Aleksandrowicz Date: Mon, 13 Jan 2025 11:20:11 +0200 Subject: [PATCH 1/8] Port linear_depth_lnn to rust --- crates/accelerate/src/synthesis/linear/lnn.rs | 321 ++++++++++++++++++ crates/accelerate/src/synthesis/linear/mod.rs | 2 + qiskit/synthesis/linear/linear_depth_lnn.py | 215 +----------- .../synthesis/linear_phase/cx_cz_depth_lnn.py | 4 +- ...ust-linear-depth-lnn-360ae5fd4f284681.yaml | 5 + 5 files changed, 333 insertions(+), 214 deletions(-) create mode 100644 crates/accelerate/src/synthesis/linear/lnn.rs create mode 100644 releasenotes/notes/rust-linear-depth-lnn-360ae5fd4f284681.yaml diff --git a/crates/accelerate/src/synthesis/linear/lnn.rs b/crates/accelerate/src/synthesis/linear/lnn.rs new file mode 100644 index 000000000000..1bd23f92303b --- /dev/null +++ b/crates/accelerate/src/synthesis/linear/lnn.rs @@ -0,0 +1,321 @@ +// 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::{_add_row_or_col, calc_inverse_matrix_inner}; +use ndarray::{Array1, Array2, ArrayView1, ArrayView2, ArrayViewMut2}; +use numpy::PyReadonlyArray2; + +use pyo3::prelude::*; + +/// 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 `_. + +type InstructionList = Vec<(usize, usize)>; +fn _row_sum(row_1: ArrayView1, row_2: ArrayView1) -> Option> { + let n = row_1.len(); + if row_2.len() != n { + return None; + } + let result: Array1 = (0..n).map(|i| row_1[i] ^ row_2[i]).collect(); + Some(result) +} +/// Perform ROW operation on a matrix mat +fn _row_op(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { + _add_row_or_col(mat, &false, ctrl, trgt); +} + +/// Perform COL operation on a matrix mat (in the inverse direction) +fn _col_op(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { + _add_row_or_col(mat, &true, trgt, ctrl); +} + +/// Add a cx gate to the instructions and update the matrix mat +fn _row_op_update_instructions( + cx_instructions: &mut InstructionList, + mat: ArrayViewMut2, + 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]. +fn _get_lower_triangular( + n: usize, + mat: ArrayView2, + mat_inv: ArrayView2, +) -> (Array2, Array2) { + let mut mat = mat.to_owned(); + let mut mat_t = mat.to_owned(); + let mut mat_inv_t = mat_inv.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 is unchanged. + for i in (0..n).rev() { + let mut found_first = false; + let mut first_j = 0; + // 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 (0..n).rev() { + if mat[[i, j]] { + if !found_first { + found_first = true; + first_j = j; + } else { + // cx_instructions_cols (L instructions) are not needed + _col_op(mat.view_mut(), j, first_j); + } + } + } + // Use row operations directed upwards to zero out all "1"s above the remaining "1" in row i + for k in (0..i).rev() { + if mat[[k, first_j]] { + _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_t.view_mut(), ctrl, trgt); + } + (mat_t, mat_inv_t) +} + +/// For each row in mat_t, save the column index of the last "1" +fn _get_label_arr(n: usize, mat_t: ArrayView2) -> Vec { + let mut label_arr: Vec = Vec::new(); + for i in 0..n { + let mut j = 0; + while !mat_t[[i, n - 1 - j]] { + j += 1; + } + label_arr.push(j); + } + label_arr +} + +/// 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, + row: ArrayView1, + k: usize, +) -> bool { + let n = row.len(); + let indx_k = label_arr_t[k]; + let mut w_needed: Array1 = Array1::from_elem(n, false); + // Find the linear combination of mat_t rows which produces "row" + for row_l in 0..n { + 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 = _row_sum(w_needed.view(), mat_inv_t.row(row_l)).unwrap(); + } + } + // If the linear combination requires the row labeled by k + !w_needed[indx_k] +} + +/// Returns label_arr_t = label_arr^(-1) +fn _get_label_arr_t(n: usize, label_arr: &[usize]) -> Vec { + let mut label_err_t: Vec = vec![0; n]; + for i in 0..n { + label_err_t[label_arr[i]] = i + } + label_err_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, + mat_inv: ArrayView2, +) -> InstructionList { + // The rows of mat_t hold all w_j vectors (see [1]). mat_inv_t is the inverted matrix of mat_t + let (mat_t, mat_inv_t) = _get_lower_triangular(n, mat.view(), mat_inv.view()); + // 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 + 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 and that writes a vector in W to wire ." + // (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) -> InstructionList { + // At start the labels are in reversed order + let mut label_arr: Vec = (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 +} + +/// 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 +#[pyfunction] +#[pyo3(signature = (mat))] +pub fn optimize_cx_circ_depth_5n_line( + _py: Python, + mat: PyReadonlyArray2, +) -> PyResult<(InstructionList, InstructionList)> { + let arrayview = mat.as_array(); + + // According to [1] the synthesis is done on the inverse matrix + // so the matrix mat is inverted at this step + + let mat_inv: Array2 = 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()); + // 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()); + + let out = (cx_instructions_rows_m2nw, cx_instructions_rows_nw2id); + Ok(out) +} diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index 08a0b1e104b3..676ccb711a96 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -14,6 +14,7 @@ use crate::QiskitError; use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2, PyReadwriteArray2}; use pyo3::prelude::*; +mod lnn; mod pmh; pub mod utils; @@ -187,5 +188,6 @@ pub fn linear(m: &Bound) -> 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::optimize_cx_circ_depth_5n_line))?; Ok(()) } diff --git a/qiskit/synthesis/linear/linear_depth_lnn.py b/qiskit/synthesis/linear/linear_depth_lnn.py index 7c7360915e0f..516e79206da6 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 and that writes a vector in W to wire ." - # (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 optimize_cx_circ_depth_5n_line def synth_cnot_depth_line_kms(mat: np.ndarray[bool]) -> QuantumCircuit: @@ -267,7 +58,7 @@ def synth_cnot_depth_line_kms(mat: np.ndarray[bool]) -> QuantumCircuit: # 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) + 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]) diff --git a/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py index c0956ea3bc7c..629d1dd62f25 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 optimize_cx_circ_depth_5n_line 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 = optimize_cx_circ_depth_5n_line(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-360ae5fd4f284681.yaml b/releasenotes/notes/rust-linear-depth-lnn-360ae5fd4f284681.yaml new file mode 100644 index 000000000000..2cff39760d14 --- /dev/null +++ b/releasenotes/notes/rust-linear-depth-lnn-360ae5fd4f284681.yaml @@ -0,0 +1,5 @@ +--- +features_synthesis: + - | + The :method:`.optimize_cx_circ_depth_5n_line` synthesis method has been rewritten + in Rust for better performance. From 3f749111ae49b6b223a278e57032d1ea6b35e2ca Mon Sep 17 00:00:00 2001 From: Gadi Aleksandrowicz Date: Mon, 13 Jan 2025 14:54:18 +0200 Subject: [PATCH 2/8] Remove unnecessary release note --- .../notes/rust-linear-depth-lnn-360ae5fd4f284681.yaml | 5 ----- 1 file changed, 5 deletions(-) delete mode 100644 releasenotes/notes/rust-linear-depth-lnn-360ae5fd4f284681.yaml diff --git a/releasenotes/notes/rust-linear-depth-lnn-360ae5fd4f284681.yaml b/releasenotes/notes/rust-linear-depth-lnn-360ae5fd4f284681.yaml deleted file mode 100644 index 2cff39760d14..000000000000 --- a/releasenotes/notes/rust-linear-depth-lnn-360ae5fd4f284681.yaml +++ /dev/null @@ -1,5 +0,0 @@ ---- -features_synthesis: - - | - The :method:`.optimize_cx_circ_depth_5n_line` synthesis method has been rewritten - in Rust for better performance. From 4902815e1a194470e3897c74238b289ac40b23f2 Mon Sep 17 00:00:00 2001 From: Gadi Aleksandrowicz Date: Mon, 20 Jan 2025 12:43:46 +0200 Subject: [PATCH 3/8] Changes according to code review --- crates/accelerate/src/synthesis/linear/lnn.rs | 231 ++++++++++-------- crates/accelerate/src/synthesis/linear/mod.rs | 1 + .../accelerate/src/synthesis/linear/utils.rs | 10 + qiskit/synthesis/linear/linear_depth_lnn.py | 14 +- 4 files changed, 143 insertions(+), 113 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/lnn.rs b/crates/accelerate/src/synthesis/linear/lnn.rs index 1bd23f92303b..76239904a16d 100644 --- a/crates/accelerate/src/synthesis/linear/lnn.rs +++ b/crates/accelerate/src/synthesis/linear/lnn.rs @@ -10,40 +10,37 @@ // copyright notice, and modified files need to carry a notice indicating // that they have been altered from the originals. -use crate::synthesis::linear::utils::{_add_row_or_col, calc_inverse_matrix_inner}; +use crate::synthesis::linear::utils::{_col_op, _row_op, 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 `_. +// 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 `_. type InstructionList = Vec<(usize, usize)>; -fn _row_sum(row_1: ArrayView1, row_2: ArrayView1) -> Option> { - let n = row_1.len(); - if row_2.len() != n { - return None; +fn _row_sum(row_1: ArrayView1, row_2: ArrayView1) -> Result, 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()) } - let result: Array1 = (0..n).map(|i| row_1[i] ^ row_2[i]).collect(); - Some(result) -} -/// Perform ROW operation on a matrix mat -fn _row_op(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { - _add_row_or_col(mat, &false, ctrl, trgt); } - -/// Perform COL operation on a matrix mat (in the inverse direction) -fn _col_op(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { - _add_row_or_col(mat, &true, trgt, ctrl); -} - /// Add a cx gate to the instructions and update the matrix mat fn _row_op_update_instructions( cx_instructions: &mut InstructionList, @@ -55,43 +52,35 @@ fn _row_op_update_instructions( _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]. -fn _get_lower_triangular( +/// Get the instructions for a lower triangular basis change of a matrix mat. +/// See the proof of Proposition 7.3 in [1]. +fn _get_lower_triangular<'a>( n: usize, mat: ArrayView2, - mat_inv: ArrayView2, -) -> (Array2, Array2) { + mut mat_inv_t: ArrayViewMut2<'a, bool>, +) -> (Array2, ArrayViewMut2<'a, bool>) { let mut mat = mat.to_owned(); let mut mat_t = mat.to_owned(); - let mut mat_inv_t = mat_inv.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 is unchanged. + // The original Matrix mat is unchanged, but mat_inv_t is + for i in (0..n).rev() { - let mut found_first = false; - let mut first_j = 0; // 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 (0..n).rev() { - if mat[[i, j]] { - if !found_first { - found_first = true; - first_j = j; - } else { - // cx_instructions_cols (L instructions) are not needed - _col_op(mat.view_mut(), j, first_j); - } - } - } + let cols_to_update: Vec = (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(), *j, *first_j); + }); + // Use row operations directed upwards to zero out all "1"s above the remaining "1" in row i - for k in (0..i).rev() { - if mat[[k, first_j]] { - _row_op_update_instructions(&mut cx_instructions_rows, mat.view_mut(), i, k) - } - } + let rows_to_update: Vec = (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 { @@ -103,15 +92,9 @@ fn _get_lower_triangular( /// For each row in mat_t, save the column index of the last "1" fn _get_label_arr(n: usize, mat_t: ArrayView2) -> Vec { - let mut label_arr: Vec = Vec::new(); - for i in 0..n { - let mut j = 0; - while !mat_t[[i, n - 1 - j]] { - j += 1; - } - label_arr.push(j); - } - label_arr + (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 @@ -121,27 +104,18 @@ fn _in_linear_combination( row: ArrayView1, k: usize, ) -> bool { - let n = row.len(); - let indx_k = label_arr_t[k]; - let mut w_needed: Array1 = Array1::from_elem(n, false); // Find the linear combination of mat_t rows which produces "row" - for row_l in 0..n { - 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 = _row_sum(w_needed.view(), mat_inv_t.row(row_l)).unwrap(); - } - } - // If the linear combination requires the row labeled by k - !w_needed[indx_k] + !(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 { let mut label_err_t: Vec = vec![0; n]; - for i in 0..n { - label_err_t[label_arr[i]] = i - } + (0..n).for_each(|i| label_err_t[label_arr[i]] = i); label_err_t } @@ -150,10 +124,11 @@ fn _get_label_arr_t(n: usize, label_arr: &[usize]) -> Vec { fn _matrix_to_north_west( n: usize, mut mat: ArrayViewMut2, - mat_inv: ArrayView2, + mut mat_inv: ArrayViewMut2, ) -> InstructionList { // The rows of mat_t hold all w_j vectors (see [1]). mat_inv_t is the inverted matrix of mat_t - let (mat_t, mat_inv_t) = _get_lower_triangular(n, mat.view(), mat_inv.view()); + // 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()); @@ -231,7 +206,7 @@ fn _matrix_to_north_west( cx_instructions_rows } -/// Transform a north-west triangular matrix to identity in depth 3*n by Proposition 7.4 of [1] +/// 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) -> InstructionList { // At start the labels are in reversed order let mut label_arr: Vec = (0..n).rev().collect(); @@ -277,32 +252,13 @@ fn _north_west_to_identity(n: usize, mut mat: ArrayViewMut2) -> Instructio cx_instructions_rows } -/// 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 -#[pyfunction] -#[pyo3(signature = (mat))] -pub fn optimize_cx_circ_depth_5n_line( - _py: Python, - mat: PyReadonlyArray2, -) -> PyResult<(InstructionList, InstructionList)> { - let arrayview = mat.as_array(); - +/// Perform steps (a) and (b) of algorithm [1] +fn _optimize_cx_circ_depth_5n_line( + arrayview: ArrayView2, +) -> (InstructionList, InstructionList) { // According to [1] the synthesis is done on the inverse matrix // so the matrix mat is inverted at this step - - let mat_inv: Array2 = arrayview.to_owned(); + let mut mat_inv: Array2 = arrayview.to_owned(); let mut mat_cpy = calc_inverse_matrix_inner(mat_inv.view(), false).unwrap(); let n = mat_cpy.nrows(); @@ -310,12 +266,79 @@ pub fn optimize_cx_circ_depth_5n_line( // 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()); + 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()); - let out = (cx_instructions_rows_m2nw, cx_instructions_rows_nw2id); - Ok(out) + (cx_instructions_rows_m2nw, cx_instructions_rows_nw2id) +} + +/// 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 + +/// References: +/// [1]: Kutin, S., Moulton, D. P., Smithline, L. (2007). +/// Computation at a Distance. +/// `arXiv:quant-ph/0701194 `_. +#[pyfunction] +#[pyo3(signature = (mat))] +pub fn optimize_cx_circ_depth_5n_line( + _py: Python, + mat: PyReadonlyArray2, +) -> PyResult<(InstructionList, InstructionList)> { + Ok(_optimize_cx_circ_depth_5n_line(mat.as_array())) +} + +/// 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 + +/// References: +/// [1]: Kutin, S., Moulton, D. P., Smithline, L. (2007). +/// Computation at a Distance. +/// `arXiv:quant-ph/0701194 `_. +#[pyfunction] +#[pyo3(signature = (mat))] +pub fn synth_cnot_depth_line_kms(py: Python, mat: PyReadonlyArray2) -> PyResult { + let num_qubits = mat.as_array().nrows(); // is a quadratic matrix + let (cx_instructions_rows_m2nw, cx_instructions_rows_nw2id) = + _optimize_cx_circ_depth_5n_line(mat.as_array()); + + let instructions = cx_instructions_rows_m2nw + .into_iter() + .chain(cx_instructions_rows_nw2id.into_iter()) + .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 676ccb711a96..49e380a1543a 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -188,6 +188,7 @@ pub fn linear(m: &Bound) -> 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::synth_cnot_depth_line_kms))?; m.add_wrapped(wrap_pyfunction!(lnn::optimize_cx_circ_depth_5n_line))?; Ok(()) } diff --git a/crates/accelerate/src/synthesis/linear/utils.rs b/crates/accelerate/src/synthesis/linear/utils.rs index 5f1fb99fa682..e2eeb528003e 100644 --- a/crates/accelerate/src/synthesis/linear/utils.rs +++ b/crates/accelerate/src/synthesis/linear/utils.rs @@ -194,6 +194,16 @@ pub fn _add_row_or_col(mut mat: ArrayViewMut2, 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, ctrl: usize, trgt: usize) { + _add_row_or_col(mat, &false, ctrl, trgt); +} + +/// Perform COL operation on a matrix mat (in the inverse direction) +pub fn _col_op(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { + _add_row_or_col(mat, &true, trgt, ctrl); +} + /// 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 { diff --git a/qiskit/synthesis/linear/linear_depth_lnn.py b/qiskit/synthesis/linear/linear_depth_lnn.py index 516e79206da6..c7f4e5ec23b1 100644 --- a/qiskit/synthesis/linear/linear_depth_lnn.py +++ b/qiskit/synthesis/linear/linear_depth_lnn.py @@ -26,7 +26,7 @@ from qiskit.exceptions import QiskitError from qiskit.circuit import QuantumCircuit from qiskit.synthesis.linear.linear_matrix_utils import check_invertible_binary_matrix -from qiskit._accelerate.synthesis.linear import optimize_cx_circ_depth_5n_line +from qiskit._accelerate.synthesis.linear import synth_cnot_depth_line_kms as fast_kms def synth_cnot_depth_line_kms(mat: np.ndarray[bool]) -> QuantumCircuit: @@ -57,11 +57,7 @@ def synth_cnot_depth_line_kms(mat: np.ndarray[bool]) -> QuantumCircuit: # 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) From d29d7ee261a56d9b53d474be0fbd55e942d1db6d Mon Sep 17 00:00:00 2001 From: Gadi Aleksandrowicz Date: Mon, 20 Jan 2025 15:11:39 +0200 Subject: [PATCH 4/8] More corrections based on the review --- crates/accelerate/src/synthesis/linear/lnn.rs | 18 ++++-------------- .../accelerate/src/synthesis/linear/utils.rs | 19 +++++++++++++++++-- 2 files changed, 21 insertions(+), 16 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/lnn.rs b/crates/accelerate/src/synthesis/linear/lnn.rs index 76239904a16d..c065e93173ad 100644 --- a/crates/accelerate/src/synthesis/linear/lnn.rs +++ b/crates/accelerate/src/synthesis/linear/lnn.rs @@ -10,7 +10,7 @@ // 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, calc_inverse_matrix_inner}; +use crate::synthesis::linear::utils::{_col_op_inv, _row_op, _row_sum, calc_inverse_matrix_inner}; use ndarray::{Array1, Array2, ArrayView1, ArrayView2, ArrayViewMut2}; use numpy::PyReadonlyArray2; use smallvec::smallvec; @@ -30,17 +30,7 @@ use qiskit_circuit::Qubit; // `arXiv:quant-ph/0701194 `_. type InstructionList = Vec<(usize, usize)>; -fn _row_sum(row_1: ArrayView1, row_2: ArrayView1) -> Result, 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()) - } -} + /// Add a cx gate to the instructions and update the matrix mat fn _row_op_update_instructions( cx_instructions: &mut InstructionList, @@ -73,7 +63,7 @@ fn _get_lower_triangular<'a>( let cols_to_update: Vec = (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(), *j, *first_j); + _col_op_inv(mat.view_mut(), *j, *first_j); }); // Use row operations directed upwards to zero out all "1"s above the remaining "1" in row i @@ -85,7 +75,7 @@ fn _get_lower_triangular<'a>( // 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_t.view_mut(), ctrl, trgt); + _col_op_inv(mat_inv_t.view_mut(), ctrl, trgt); } (mat_t, mat_inv_t) } diff --git a/crates/accelerate/src/synthesis/linear/utils.rs b/crates/accelerate/src/synthesis/linear/utils.rs index e2eeb528003e..ba598ad44d4c 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}; @@ -200,10 +202,23 @@ pub fn _row_op(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { } /// Perform COL operation on a matrix mat (in the inverse direction) -pub fn _col_op(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { +pub fn _col_op_inv(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { _add_row_or_col(mat, &true, trgt, ctrl); } +/// Returns the element-wise sum of two boolean rows (i.e. addition modulo 2) +pub fn _row_sum(row_1: ArrayView1, row_2: ArrayView1) -> Result, 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) -> Array2 { let mut rng = match seed { From e3c69e26d92dd47c160424674ae427be93bcf0fd Mon Sep 17 00:00:00 2001 From: Gadi Aleksandrowicz Date: Thu, 23 Jan 2025 13:54:19 +0200 Subject: [PATCH 5/8] Changes according to code review --- crates/accelerate/src/synthesis/linear/lnn.rs | 40 ++++++++++--------- crates/accelerate/src/synthesis/linear/mod.rs | 4 +- .../accelerate/src/synthesis/linear/utils.rs | 6 +-- qiskit/synthesis/linear/linear_depth_lnn.py | 4 +- .../synthesis/linear_phase/cx_cz_depth_lnn.py | 4 +- ...ust-linear-depth-lnn-69532c9c6b86ffe8.yaml | 5 +++ 6 files changed, 34 insertions(+), 29 deletions(-) create mode 100644 releasenotes/notes/rust-linear-depth-lnn-69532c9c6b86ffe8.yaml diff --git a/crates/accelerate/src/synthesis/linear/lnn.rs b/crates/accelerate/src/synthesis/linear/lnn.rs index c065e93173ad..cfea34a7e4eb 100644 --- a/crates/accelerate/src/synthesis/linear/lnn.rs +++ b/crates/accelerate/src/synthesis/linear/lnn.rs @@ -10,7 +10,7 @@ // 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_inv, _row_op, _row_sum, calc_inverse_matrix_inner}; +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; @@ -44,10 +44,12 @@ fn _row_op_update_instructions( /// 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, - mut mat_inv_t: ArrayViewMut2<'a, bool>, + mut mat_inv: ArrayViewMut2<'a, bool>, ) -> (Array2, ArrayViewMut2<'a, bool>) { let mut mat = mat.to_owned(); let mut mat_t = mat.to_owned(); @@ -55,7 +57,7 @@ fn _get_lower_triangular<'a>( 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_t is + // 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 @@ -63,7 +65,7 @@ fn _get_lower_triangular<'a>( let cols_to_update: Vec = (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_inv(mat.view_mut(), *j, *first_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 @@ -75,9 +77,9 @@ fn _get_lower_triangular<'a>( // 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_inv(mat_inv_t.view_mut(), ctrl, trgt); + _col_op(mat_inv.view_mut(), trgt, ctrl); // performs an inverted col_op } - (mat_t, mat_inv_t) + (mat_t, mat_inv) } /// For each row in mat_t, save the column index of the last "1" @@ -104,9 +106,9 @@ fn _in_linear_combination( /// Returns label_arr_t = label_arr^(-1) fn _get_label_arr_t(n: usize, label_arr: &[usize]) -> Vec { - let mut label_err_t: Vec = vec![0; n]; - (0..n).for_each(|i| label_err_t[label_arr[i]] = i); - label_err_t + let mut label_arr_t: Vec = 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 @@ -243,9 +245,7 @@ fn _north_west_to_identity(n: usize, mut mat: ArrayViewMut2) -> Instructio } /// Perform steps (a) and (b) of algorithm [1] -fn _optimize_cx_circ_depth_5n_line( - arrayview: ArrayView2, -) -> (InstructionList, InstructionList) { +fn _synth_cnot_lnn_instructions(arrayview: ArrayView2) -> (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 = arrayview.to_owned(); @@ -266,7 +266,7 @@ fn _optimize_cx_circ_depth_5n_line( (cx_instructions_rows_m2nw, cx_instructions_rows_nw2id) } -/// Optimize CX circuit in depth bounded by 5n for LNN connectivity. +/// 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). @@ -287,14 +287,13 @@ fn _optimize_cx_circ_depth_5n_line( /// `arXiv:quant-ph/0701194 `_. #[pyfunction] #[pyo3(signature = (mat))] -pub fn optimize_cx_circ_depth_5n_line( - _py: Python, +pub fn py_synth_cnot_lnn_instructions( mat: PyReadonlyArray2, ) -> PyResult<(InstructionList, InstructionList)> { - Ok(_optimize_cx_circ_depth_5n_line(mat.as_array())) + Ok(_synth_cnot_lnn_instructions(mat.as_array())) } -/// Optimize CX circuit in depth bounded by 5n for LNN connectivity. +/// 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). @@ -315,10 +314,13 @@ pub fn optimize_cx_circ_depth_5n_line( /// `arXiv:quant-ph/0701194 `_. #[pyfunction] #[pyo3(signature = (mat))] -pub fn synth_cnot_depth_line_kms(py: Python, mat: PyReadonlyArray2) -> PyResult { +pub fn py_synth_cnot_depth_line_kms( + py: Python, + mat: PyReadonlyArray2, +) -> PyResult { let num_qubits = mat.as_array().nrows(); // is a quadratic matrix let (cx_instructions_rows_m2nw, cx_instructions_rows_nw2id) = - _optimize_cx_circ_depth_5n_line(mat.as_array()); + _synth_cnot_lnn_instructions(mat.as_array()); let instructions = cx_instructions_rows_m2nw .into_iter() diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index 49e380a1543a..a55a55769ae7 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -188,7 +188,7 @@ pub fn linear(m: &Bound) -> 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::synth_cnot_depth_line_kms))?; - m.add_wrapped(wrap_pyfunction!(lnn::optimize_cx_circ_depth_5n_line))?; + 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 ba598ad44d4c..c620acb5eec0 100644 --- a/crates/accelerate/src/synthesis/linear/utils.rs +++ b/crates/accelerate/src/synthesis/linear/utils.rs @@ -201,9 +201,9 @@ pub fn _row_op(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { _add_row_or_col(mat, &false, ctrl, trgt); } -/// Perform COL operation on a matrix mat (in the inverse direction) -pub fn _col_op_inv(mat: ArrayViewMut2, ctrl: usize, trgt: usize) { - _add_row_or_col(mat, &true, trgt, ctrl); +/// Perform COL operation on a matrix mat +pub fn _col_op(mat: ArrayViewMut2, 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) diff --git a/qiskit/synthesis/linear/linear_depth_lnn.py b/qiskit/synthesis/linear/linear_depth_lnn.py index c7f4e5ec23b1..fdbabd68f025 100644 --- a/qiskit/synthesis/linear/linear_depth_lnn.py +++ b/qiskit/synthesis/linear/linear_depth_lnn.py @@ -26,7 +26,7 @@ from qiskit.exceptions import QiskitError from qiskit.circuit import QuantumCircuit from qiskit.synthesis.linear.linear_matrix_utils import check_invertible_binary_matrix -from qiskit._accelerate.synthesis.linear import synth_cnot_depth_line_kms as fast_kms +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: @@ -55,8 +55,6 @@ 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 circuit_data = fast_kms(mat) # construct circuit from the data diff --git a/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py b/qiskit/synthesis/linear_phase/cx_cz_depth_lnn.py index 629d1dd62f25..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._accelerate.synthesis.linear 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 From c40c0e4183667996cab353924a9b205280484c0f Mon Sep 17 00:00:00 2001 From: Gadi Aleksandrowicz Date: Thu, 23 Jan 2025 15:25:01 +0200 Subject: [PATCH 6/8] Small fix --- crates/accelerate/src/synthesis/linear/lnn.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/accelerate/src/synthesis/linear/lnn.rs b/crates/accelerate/src/synthesis/linear/lnn.rs index cfea34a7e4eb..8049c3019844 100644 --- a/crates/accelerate/src/synthesis/linear/lnn.rs +++ b/crates/accelerate/src/synthesis/linear/lnn.rs @@ -324,7 +324,7 @@ pub fn py_synth_cnot_depth_line_kms( let instructions = cx_instructions_rows_m2nw .into_iter() - .chain(cx_instructions_rows_nw2id.into_iter()) + .chain(cx_instructions_rows_nw2id) .map(|(ctrl, target)| { ( StandardGate::CXGate, From 5d5ed8fd8f6e78437ae4a0234560f9b916b0be20 Mon Sep 17 00:00:00 2001 From: Gadi Aleksandrowicz Date: Mon, 27 Jan 2025 10:46:09 +0200 Subject: [PATCH 7/8] Docstring change according to PR review --- crates/accelerate/src/synthesis/linear/lnn.rs | 65 ++++++++----------- 1 file changed, 27 insertions(+), 38 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/lnn.rs b/crates/accelerate/src/synthesis/linear/lnn.rs index 8049c3019844..b14d975b4589 100644 --- a/crates/accelerate/src/synthesis/linear/lnn.rs +++ b/crates/accelerate/src/synthesis/linear/lnn.rs @@ -244,7 +244,25 @@ fn _north_west_to_identity(n: usize, mut mat: ArrayViewMut2) -> Instructio cx_instructions_rows } -/// Perform steps (a) and (b) of algorithm [1] +/// 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 `_. fn _synth_cnot_lnn_instructions(arrayview: ArrayView2) -> (InstructionList, InstructionList) { // According to [1] the synthesis is done on the inverse matrix // so the matrix mat is inverted at this step @@ -266,25 +284,11 @@ fn _synth_cnot_lnn_instructions(arrayview: ArrayView2) -> (InstructionList (cx_instructions_rows_m2nw, cx_instructions_rows_nw2id) } -/// 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 `_. +/// 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 `_. +/// 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( @@ -294,24 +298,9 @@ pub fn py_synth_cnot_lnn_instructions( } /// 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 `_. +/// Uses the algorithm by Kutin, Moulton, Smithline +/// described in `arXiv:quant-ph/0701194 `_. +/// Returns: The CircuitData of the synthesized circuit. #[pyfunction] #[pyo3(signature = (mat))] pub fn py_synth_cnot_depth_line_kms( From 3220dcd8b06ea04139a6597cfc21e54ec2866f97 Mon Sep 17 00:00:00 2001 From: Gadi Aleksandrowicz Date: Mon, 27 Jan 2025 12:12:25 +0200 Subject: [PATCH 8/8] Cargo fmt for version 1.79 --- crates/accelerate/src/synthesis/linear/lnn.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/crates/accelerate/src/synthesis/linear/lnn.rs b/crates/accelerate/src/synthesis/linear/lnn.rs index b14d975b4589..edb1ae0913a2 100644 --- a/crates/accelerate/src/synthesis/linear/lnn.rs +++ b/crates/accelerate/src/synthesis/linear/lnn.rs @@ -285,7 +285,7 @@ fn _synth_cnot_lnn_instructions(arrayview: ArrayView2) -> (InstructionList } /// Find instruction to synthesize CX circuit in depth bounded by 5n for LNN connectivity. -/// Uses the algorithm by Kutin, Moulton, Smithline +/// Uses the algorithm by Kutin, Moulton, Smithline /// described in `arXiv:quant-ph/0701194 `_. /// Returns: Tuple with two lists of instructions for CX gates /// Corresponding to the two parts of the algorithm @@ -298,7 +298,7 @@ pub fn py_synth_cnot_lnn_instructions( } /// Synthesize CX circuit in depth bounded by 5n for LNN connectivity. -/// Uses the algorithm by Kutin, Moulton, Smithline +/// Uses the algorithm by Kutin, Moulton, Smithline /// described in `arXiv:quant-ph/0701194 `_. /// Returns: The CircuitData of the synthesized circuit. #[pyfunction]