Improve LinearFunction synthesis

qiskit/circuit/library/generalized_gates/linear_function.py

@@ -16,6 +16,7 @@
 import numpy as np
 from qiskit.circuit import QuantumCircuit, Gate
 from qiskit.circuit.exceptions import CircuitError
+from qiskit.synthesis.linear import check_invertible_binary_matrix
 
 
 class LinearFunction(Gate):
@@ -110,8 +111,7 @@ def __init__(
 
             # Optionally, check that the matrix is invertible
             if validate_input:
-                det = np.linalg.det(linear) % 2
-                if not np.allclose(det, 1):
+                if not check_invertible_binary_matrix(linear):
                     raise CircuitError(
                         "A linear function must be represented by an invertible matrix."
                     )
@@ -134,7 +134,7 @@ def synthesize(self):
         Returns:
             QuantumCircuit: A circuit implementing the evolution.
         """
-        from qiskit.transpiler.synthesis import cnot_synth
+        from qiskit.synthesis.linear import cnot_synth
 
         return cnot_synth(self.linear)
 

qiskit/quantum_info/operators/dihedral/random.py

@@ -48,10 +48,9 @@ def random_cnotdihedral(num_qubits, seed=None):
 
     # Random affine function
     # Random invertible binary matrix
-    det = 0
-    while np.allclose(det, 0) or np.allclose(det, 2):
-        linear = rng.integers(2, size=(num_qubits, num_qubits))
-        det = np.linalg.det(linear) % 2
+    from qiskit.synthesis.linear import random_invertible_binary_matrix
+
+    linear = random_invertible_binary_matrix(num_qubits, seed=rng)
     elem.linear = linear
 
     # Random shift

qiskit/synthesis/linear/__init__.py

@@ -0,0 +1,21 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2018.
+#
+# 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.
+
+"""Module containing cnot circuits and cnot-phase circuit synthesize."""
+
+
+from .graysynth import graysynth, cnot_synth
+from .linear_matrix_utils import (
+    random_invertible_binary_matrix,
+    calc_inverse_matrix,
+    check_invertible_binary_matrix,
+)

qiskit/transpiler/synthesis/graysynth.py → qiskit/synthesis/linear/graysynth.py

@@ -176,7 +176,7 @@ def graysynth(cnots, angles, section_size=2):
         else:
             sta.append([cnots1, list(set(ilist).difference([j])), qubit])
         sta.append([cnots0, list(set(ilist).difference([j])), qubit])
-    qcir += cnot_synth(state, section_size).inverse()
+    qcir = qcir.compose(cnot_synth(state, section_size).inverse())
     return qcir
 
 

qiskit/synthesis/linear/linear_circuits_utils.py

@@ -0,0 +1,66 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 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.
+
+"""Utility functions for handling linear reversible circuits."""
+
+import copy
+import numpy as np
+from . import calc_inverse_matrix
+
+
+def transpose_cx_circ(qc):
+    """Transpose all cx gates in a circuit."""
+    data = qc.data
+    for i, _ in enumerate(data):
+        if data[i][0].name == "cx" and data[i].operation.num_qubits == 2:
+            data[i][1][0], data[i][1][1] = data[i][1][1], data[i][1][0]
+
+
+def optimize_cx_4_options(function, mat, optimize_count=True):
+    """Get best implementation of CX, implementing M,M^(-1),M^T,M^(-1)^T"""
+    qc = function(mat)
+    best_qc = qc
+    best_depth = qc.depth()
+    best_count = qc.count_ops()["cx"]
+
+    for i in range(1, 4):
+        mat_cpy = copy.deepcopy(mat)
+        # i=1 inverse, i=2 transpose, i=3 transpose and inverse
+        if i == 1:
+            mat_cpy = calc_inverse_matrix(mat_cpy)
+            qc = function(mat_cpy)
+            qc = qc.inverse()
+        elif i == 2:
+            mat_cpy = np.transpose(mat_cpy)
+            qc = function(mat_cpy)
+            transpose_cx_circ(qc)
+            qc = qc.inverse()
+        elif i == 3:
+            mat_cpy = calc_inverse_matrix(np.transpose(mat_cpy))
+            qc = function(mat_cpy)
+            transpose_cx_circ(qc)
+
+        new_depth = qc.depth()
+        new_count = qc.count_ops()["cx"]
+        better_count = (optimize_count and best_count > new_count) or (
+            not optimize_count and best_depth == new_depth and best_count > new_count
+        )
+        better_depth = (not optimize_count and best_depth > new_depth) or (
+            optimize_count and best_count == new_count and best_depth > new_depth
+        )
+
+        if better_count or better_depth:
+            best_count = new_count
+            best_depth = new_depth
+            best_qc = qc
+
+    return best_qc

qiskit/synthesis/linear/linear_matrix_utils.py

@@ -0,0 +1,119 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 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.
+
+"""Utility functions for handling binary matrices."""
+
+import numpy as np
+from qiskit.exceptions import QiskitError
+
+
+def check_invertible_binary_matrix(mat):
+    """Check that a binary matrix is invertible."""
+    if len(mat.shape) != 2 or mat.shape[0] != mat.shape[1]:
+        return False
+
+    mat = _gauss_elimination(mat)
+    rank = _compute_rank_after_gauss_elim(mat)
+    return rank == mat.shape[0]
+
+
+def random_invertible_binary_matrix(num_qubits, seed=None):
+    """Generates a random invertible n x n binary matrix."""
+    if isinstance(seed, np.random.Generator):
+        rng = seed
+    else:
+        rng = np.random.default_rng(seed)
+
+    rank = 0
+    while rank != num_qubits:
+        mat = rng.integers(2, size=(num_qubits, num_qubits))
+        mat = _gauss_elimination(mat)
+        rank = _compute_rank_after_gauss_elim(mat)
+    return mat
+
+
+def _gauss_elimination(mat, ncols=None, full_elim=False):
+    """Gauss elimination of a matrix mat with m rows and n columns.
+    If full_elim = True, it allows full elimination of mat[:, 0 : ncols]
+    Returns the matrix mat."""
+
+    # Treat the matrix A as containing integer values
+    mat = np.array(mat, dtype=int, copy=False)
+
+    m = mat.shape[0]  # no. of rows
+    n = mat.shape[1]  # no. of columns
+    if ncols is not None:
+        n = min(n, ncols)  # no. of active columns
+
+    r = 0  # current rank
+    k = 0  # current pivot column
+    while (r < m) and (k < n):
+        is_non_zero = False
+        new_r = r
+        for j in range(k, n):
+            for i in range(r, m):
+                if mat[i][j]:
+                    is_non_zero = True
+                    k = j
+                    new_r = i
+                    break
+            if is_non_zero:
+                break
+        if not is_non_zero:
+            return mat  # A is in the canonical form
+
+        if new_r != r:
+            mat[[r, new_r]] = mat[[new_r, r]]
+
+        if full_elim:
+            for i in range(0, r):
+                if mat[i][k]:
+                    mat[i] = mat[i] ^ mat[r]
+
+        for i in range(r + 1, m):
+            if mat[i][k]:
+                mat[i] = mat[i] ^ mat[r]
+        r += 1
+
+    return mat
+
+
+def calc_inverse_matrix(mat, verify=False):
+    """Given a square numpy(dtype=int) matrix mat, tries to compute its inverse.
+    Returns None if the matrix is not invertible, and the inverse otherwise."""
+
+    if mat.shape[0] != mat.shape[1]:
+        raise QiskitError("Matrix to invert is a non-square matrix.")
+
+    n = mat.shape[0]
+    # concatenate the matrix and identity
+    mat1 = np.concatenate((mat, np.eye(n, dtype=int)), axis=1)
+    mat1, _ = _gauss_elimination(mat1, None, full_elim=True)
+
+    r = _compute_rank_after_gauss_elim(mat1[:, 0:n])
+
+    if r < n:
+        raise QiskitError("The matrix is not invertible")
+
+    matinv = mat1[:, n : 2 * n]
+
+    if verify:
+        mat2 = np.dot(mat, matinv) % 2
+        assert np.array_equal(mat2, np.eye(n))
+
+    return matinv
+
+
+def _compute_rank_after_gauss_elim(mat):
+    """Given a matrix A after Gaussian elimination, computes its rank
+    (i.e. simply the number of nonzero rows)"""
+    return np.sum(mat.any(axis=1))

qiskit/transpiler/synthesis/__init__.py

@@ -11,6 +11,3 @@
 # that they have been altered from the originals.
 
 """Module containing transpiler synthesize."""
-
-
-from .graysynth import graysynth, cnot_synth

test/python/circuit/library/test_linear_function.py

@@ -22,6 +22,7 @@
 from qiskit.circuit.library.standard_gates import CXGate, SwapGate
 from qiskit.circuit.library.generalized_gates import LinearFunction
 from qiskit.circuit.exceptions import CircuitError
+from qiskit.synthesis.linear import random_invertible_binary_matrix
 
 from qiskit.quantum_info.operators import Operator
 
@@ -51,23 +52,6 @@ def random_linear_circuit(num_qubits, num_gates, seed=None):
     return circ
 
 
-def random_invertible_binary_matrix(num_qubits, seed=None):
-    """Generates a random invertible n x n binary matrix."""
-
-    # This code is adapted from random_cnotdihedral
-    if isinstance(seed, np.random.Generator):
-        rng = seed
-    else:
-        rng = np.random.default_rng(seed)
-
-    det = 0
-    while np.allclose(det, 0) or np.allclose(det, 2):
-        binary_matrix = rng.integers(2, size=(num_qubits, num_qubits))
-        det = np.linalg.det(binary_matrix) % 2
-
-    return binary_matrix
-
-
 @ddt
 class TestLinearFunctions(QiskitTestCase):
     """Tests for clifford append gate functions."""
@@ -132,7 +116,7 @@ def test_conversion_to_linear_function_and_back(self, num_qubits):
     def test_patel_markov_hayes(self):
         """Checks the explicit example from Patel-Markov-Hayes's paper."""
 
-        # This code is adapted from test_synthesis.py
+        # This code is adapted from test_gray_synthesis.py
         binary_matrix = [
             [1, 1, 0, 0, 0, 0],
             [1, 0, 0, 1, 1, 0],

test/python/synthesis/__init__.py

@@ -0,0 +1,13 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 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.
+
+"""Qiskit synthesis unit tests."""

test/python/transpiler/test_synthesis.py → test/python/synthesis/test_gray_synthesis.py

@@ -10,12 +10,14 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-"""Test synthesis algorithms"""
+"""Test cnot circuit and cnot-phase circuit synthesis algorithms"""
+
+import unittest
 
 from qiskit.circuit import QuantumCircuit, QuantumRegister
 from qiskit.quantum_info.operators import Operator
 from qiskit.extensions.unitary import UnitaryGate
-from qiskit.transpiler.synthesis import graysynth, cnot_synth
+from qiskit.synthesis.linear import graysynth, cnot_synth
 from qiskit.test import QiskitTestCase
 
 
@@ -262,3 +264,7 @@ def test_patel_markov_hayes(self):
 
         # Check if the two circuits are equivalent
         self.assertEqual(unitary_patel, unitary_compare)
+
+
+if __name__ == "__main__":
+    unittest.main()