Adding linear synthesis algorithm for LNN (#9098)

* update init file

* add a test for kms synthesis for lnn

* add a synthesis algorithm based on kms for lnn in depth 5n.

Co-authored-by: Ben Zindorf <benzindorf@gmail.com>

* minor fixes in linear_depth_lnn

* fix import

* add release notes

* updates following review

* more updates following review

* add synth_cnot_depth_line_kms to docs

* updates following review comments

* fix init

* add Patel-Markov-Hayes to docs

* update docstring

* add Patel-Markov-Hayes to release notes

Co-authored-by: Ben Zindorf <benzindorf@gmail.com>
Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com>

qiskit/synthesis/__init__.py

@@ -29,6 +29,14 @@ Evolution Synthesis
    SuzukiTrotter
    MatrixExponential
 
+Linear Function Synthesis
+=========================
+.. autosummary::
+   :toctree: ../stubs/
+
+    synth_cnot_count_full_pmh
+    synth_cnot_depth_line_kms
+
 Permutation Synthesis
 =====================
 
@@ -48,4 +56,5 @@ from .evolution import (
     QDrift,
 )
 
+from .linear import synth_cnot_count_full_pmh, synth_cnot_depth_line_kms
 from .permutation import synth_permutation_depth_lnn_kms

qiskit/synthesis/linear/__init__.py

@@ -14,6 +14,7 @@
 
 
 from .graysynth import graysynth, synth_cnot_count_full_pmh
+from .linear_depth_lnn import synth_cnot_depth_line_kms
 from .linear_matrix_utils import (
     random_invertible_binary_matrix,
     calc_inverse_matrix,

qiskit/synthesis/linear/graysynth.py

@@ -182,27 +182,30 @@ def graysynth(cnots, angles, section_size=2):
 
 def synth_cnot_count_full_pmh(state, section_size=2):
     """
-    This function is an implementation of the Patel–Markov–Hayes algorithm
-    for optimal synthesis of linear reversible circuits, as specified by an
-    n x n matrix.
+    Synthesize linear reversible circuits for all-to-all architecture
+    using Patel, Markov and Hayes method.
 
-    The algorithm is described in detail in the following paper:
-    "Optimal synthesis of linear reversible circuits."
-    Patel, Ketan N., Igor L. Markov, and John P. Hayes.
-    Quantum Information & Computation 8.3 (2008): 282-294.
+    This function is an implementation of the Patel, Markov and Hayes algorithm from [1]
+    for optimal synthesis of linear reversible circuits for all-to-all architecture,
+    as specified by an n x n matrix.
 
     Args:
-        state (list[list] or ndarray): n x n matrix, describing the state
+        state (list[list] or ndarray): n x n boolean invertible matrix, describing the state
             of the input circuit
-        section_size (int): the size of each section, used in _lwr_cnot_synth(), in the
-            Patel–Markov–Hayes algorithm. section_size must be a factor of num_qubits.
+        section_size (int): the size of each section, used in the
+            Patel–Markov–Hayes algorithm [1]. section_size must be a factor of num_qubits.
 
     Returns:
-        QuantumCircuit: a CNOT-only circuit implementing the
-            desired linear transformation
+        QuantumCircuit: a CX-only circuit implementing the linear transformation.
 
     Raises:
-        QiskitError: when variable "state" isn't of type numpy.matrix
+        QiskitError: when variable "state" isn't of type numpy.ndarray
+
+    References:
+        1. Patel, Ketan N., Igor L. Markov, and John P. Hayes,
+           *Optimal synthesis of linear reversible circuits*,
+           Quantum Information & Computation 8.3 (2008): 282-294.
+           `arXiv:quant-ph/0302002 [quant-ph] <https://arxiv.org/abs/quant-ph/0302002>`_
     """
     if not isinstance(state, (list, np.ndarray)):
         raise QiskitError(

qiskit/synthesis/linear/linear_depth_lnn.py

@@ -0,0 +1,275 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022
+#
+# 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.
+
+"""
+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>`_.
+"""
+
+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 originl 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.copy()
+    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
+
+
+def synth_cnot_depth_line_kms(mat):
+    """
+    Synthesize linear reversible circuit for linear nearest-neighbor architectures using
+    Kutin, Moulton, Smithline method.
+
+    Synthesis algorithm for linear reversible circuits from [1], Chapter 7.
+    Synthesizes any linear reversible circuit of n qubits over linear nearest-neighbor
+    architecture using CX gates with depth at most 5*n.
+
+    Args:
+        mat(np.ndarray]): A boolean invertible matrix.
+
+    Returns:
+        QuantumCircuit: the synthesized quantum circuit.
+
+    Raises:
+        QiskitError: if mat is not invertible.
+
+    References:
+        1. Kutin, S., Moulton, D. P., Smithline, L.,
+           *Computation at a distance*, Chicago J. Theor. Comput. Sci., vol. 2007, (2007),
+           `arXiv:quant-ph/0701194 <https://arxiv.org/abs/quant-ph/0701194>`_
+    """
+    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

qiskit/synthesis/linear/linear_matrix_utils.py

@@ -147,3 +147,13 @@ 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 _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]

releasenotes/notes/add-linear-synthesis-lnn-depth-5n-36c1aeda02b8bc6f.yaml

@@ -0,0 +1,15 @@
+---
+features:
+  - |
+    Added a depth-efficient synthesis algorithm
+    :func:`qiskit.synthesis.linear.linear_depth_lnn.synth_cnot_depth_line_kms`
+    for linear reversible circuits :class:`~qiskit.circuit.library.LinearFunction`
+    over the linear nearest-neighbor architecture,
+    following the paper <https://arxiv.org/abs/quant-ph/0701194>`__.
+upgrade:
+  - |
+    Added to the Qiskit documentation the synthesis algorithm
+    :func:`qiskit.synthesis.linear.linear_depth_lnn.synth_cnot_count_full_pmh`
+    for linear reversible circuits :class:`~qiskit.circuit.library.LinearFunction`
+    for all-to-all architecture, following the paper
+    <https://arxiv.org/abs/quant-ph/0302002>`__.

test/python/synthesis/test_linear_synthesis.py

@@ -15,10 +15,12 @@
 import unittest
 
 import numpy as np
+from ddt import ddt, data
 from qiskit import QuantumCircuit
 from qiskit.circuit.library import LinearFunction
 from qiskit.synthesis.linear import (
     synth_cnot_count_full_pmh,
+    synth_cnot_depth_line_kms,
     random_invertible_binary_matrix,
     check_invertible_binary_matrix,
     calc_inverse_matrix,
@@ -27,6 +29,7 @@ from qiskit.synthesis.linear.linear_circuits_utils import transpose_cx_circ, opt
 from qiskit.test import QiskitTestCase
 
 
+@ddt
 class TestLinearSynth(QiskitTestCase):
     """Test the linear reversible circuit synthesis functions."""
 
@@ -99,9 +102,9 @@ class TestLinearSynth(QiskitTestCase):
         self.assertEqual(optimized_qc.depth(), 15)
         self.assertEqual(optimized_qc.count_ops()["cx"], 23)
 
-    def test_invertible_matrix(self):
+    @data(5, 6)
+    def test_invertible_matrix(self, n):
         """Test the functions for generating a random invertible matrix and inverting it."""
-        n = 5
         mat = random_invertible_binary_matrix(n, seed=1234)
         out = check_invertible_binary_matrix(mat)
         mat_inv = calc_inverse_matrix(mat, verify=True)
@@ -109,6 +112,30 @@ class TestLinearSynth(QiskitTestCase):
         self.assertTrue(np.array_equal(mat_out, np.eye(n)))
         self.assertTrue(out)
 
+    @data(5, 6)
+    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)
+            mat = np.array(mat, dtype=bool)
+            qc = synth_cnot_depth_line_kms(mat)
+            mat1 = LinearFunction(qc).linear
+            self.assertTrue((mat == mat1).all())
+
+            # Check that the circuit depth is bounded by 5*num_qubits
+            depth = qc.depth()
+            self.assertTrue(depth <= 5 * num_qubits)
+
+            # Check that the synthesized circuit qc fits LNN connectivity
+            for inst in qc.data:
+                self.assertEqual(inst.operation.name, "cx")
+                q0 = qc.find_bit(inst.qubits[0]).index
+                q1 = qc.find_bit(inst.qubits[1]).index
+                dist = abs(q0 - q1)
+                self.assertEqual(dist, 1)
+
 
 if __name__ == "__main__":
     unittest.main()