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,100 @@
+# 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
+from typing import Callable
+import numpy as np
+from qiskit import QuantumCircuit
+from qiskit.exceptions import QiskitError
+from qiskit.circuit.exceptions import CircuitError
+from . import calc_inverse_matrix, check_invertible_binary_matrix
+
+
+def transpose_cx_circ(qc: QuantumCircuit):
+    """Takes a circuit having only CX gates, and calculates its transpose.
+    This is done by recursively replacing CX(i, j) with CX(j, i) in all instructions.
+
+    Args:
+        qc: a QuantumCircuit containing only CX gates.
+
+    Returns:
+        QuantumCircuit: the transposed circuit.
+
+    Raises:
+        CircuitError: if qc has a non-CX gate.
+    """
+    transposed_circ = QuantumCircuit(qc.qubits, qc.clbits, name=qc.name + "_transpose")
+    for instruction in reversed(qc.data):
+        if instruction.operation.name != "cx":
+            raise CircuitError("The circuit contains non-CX gates.")
+        transposed_circ._append(instruction.replace(qubits=reversed(instruction.qubits)))
+    return transposed_circ
+
+
+def optimize_cx_4_options(function: Callable, mat: np.ndarray, optimize_count: bool = True):
+    """Get the best implementation of a circuit implementing a binary invertible matrix M,
+    by considering all four options: M,M^(-1),M^T,M^(-1)^T.
+    Optimizing either the CX count or the depth.
+
+    Args:
+        function: the synthesis function.
+        mat: a binary invertible matrix.
+        optimize_count: True if the number of CX gates in optimize, False if the depth is optimized.
+
+    Returns:
+        QuantumCircuit: an optimized QuantumCircuit, has the best depth or CX count of the four options.
+
+    Raises:
+        QiskitError: if mat is not an invertible matrix.
+    """
+    if not check_invertible_binary_matrix(mat):
+        raise QiskitError("The matrix is not invertible.")
+
+    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)
+            qc = transpose_cx_circ(qc)
+        elif i == 3:
+            mat_cpy = calc_inverse_matrix(np.transpose(mat_cpy))
+            qc = function(mat_cpy)
+            qc = transpose_cx_circ(qc)
+            qc = qc.inverse()
+
+        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,149 @@
+# 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."""
+
+from typing import Optional, Union
+import numpy as np
+from qiskit.exceptions import QiskitError
+
+
+def check_invertible_binary_matrix(mat: np.ndarray):
+    """Check that a binary matrix is invertible.
+
+    Args:
+        mat: a binary matrix.
+
+    Returns:
+        bool: True if mat in invertible and False otherwise.
+    """
+    if len(mat.shape) != 2 or mat.shape[0] != mat.shape[1]:
+        return False
+
+    mat = _gauss_elimination(mat)
+    rank = _compute_rank_after_gauss_elim(mat)
+    return rank == mat.shape[0]
+
+
+def random_invertible_binary_matrix(
+    num_qubits: int, seed: Optional[Union[np.random.Generator, int]] = None
+):
+    """Generates a random invertible n x n binary matrix.
+
+    Args:
+        num_qubits: the matrix size.
+        seed: a random seed.
+
+    Returns:
+        np.ndarray: A random invertible binary matrix of size num_qubits.
+    """
+    if isinstance(seed, np.random.Generator):
+        rng = seed
+    else:
+        rng = np.random.default_rng(seed)
+
+    rank = 0
+    while rank != num_qubits:
+        mat = rng.integers(2, size=(num_qubits, num_qubits))
+        mat_gauss = mat.copy()
+        mat_gauss = _gauss_elimination(mat_gauss)
+        rank = _compute_rank_after_gauss_elim(mat_gauss)
+    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: np.ndarray, verify: bool = False):
+    """Given a square numpy(dtype=int) matrix mat, tries to compute its inverse.
+
+    Args:
+        mat: a boolean square matrix.
+        verify: if True asserts that the multiplication of mat and its inverse is the identity matrix.
+
+    Returns:
+        np.ndarray: the inverse matrix.
+
+    Raises:
+         QiskitError: if the matrix is not square.
+         QiskitError: if the matrix is not invertible.
+    """
+
+    if mat.shape[0] != mat.shape[1]:
+        raise QiskitError("Matrix to invert is a non-square matrix.")
+
+    n = mat.shape[0]
+    # concatenate the matrix and identity
+    mat1 = np.concatenate((mat, np.eye(n, dtype=int)), axis=1)
+    mat1 = _gauss_elimination(mat1, None, full_elim=True)
+
+    r = _compute_rank_after_gauss_elim(mat1[:, 0:n])
+
+    if r < n:
+        raise QiskitError("The matrix is not invertible.")
+
+    matinv = mat1[:, n : 2 * n]
+
+    if verify:
+        mat2 = np.dot(mat, matinv) % 2
+        assert np.array_equal(mat2, np.eye(n))
+
+    return matinv
+
+
+def _compute_rank_after_gauss_elim(mat):
+    """Given a matrix A after Gaussian elimination, computes its rank
+    (i.e. simply the number of nonzero rows)"""
+    return np.sum(mat.any(axis=1))

qiskit/transpiler/passes/synthesis/high_level_synthesis.py

@@ -19,7 +19,7 @@
 from qiskit.dagcircuit.dagcircuit import DAGCircuit
 from qiskit.transpiler.exceptions import TranspilerError
 from qiskit.quantum_info import decompose_clifford
-from qiskit.transpiler.synthesis import cnot_synth
+from qiskit.synthesis.linear import cnot_synth
 from .plugin import HighLevelSynthesisPluginManager, HighLevelSynthesisPlugin
 
 

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

releasenotes/notes/linear_function_synthesis_utils-f2f96924ca45e1fb.yaml

@@ -0,0 +1,23 @@
+---
+features:
+  - |
+    Add new internal utility functions to handle linear functions and circuits:
+    A utility function that checks whether a NxN binary matrix is invertible and a utility function
+    that computes an inverse of a NxN invertible binary matrix.
+    A utility function that calculates the transpose of a linear reversible circuit.
+    Implement the following "meta" idea for synthesizing linear functions:
+    given an invertible binary matrix A and a synthesis algorithm f: A -> QuantumCircuit,
+    in addition to applying f to A, we can also apply f to A^{-1} and invert the computed circuit,
+    apply f to A^{t} and transpose the computed circuit, and to apply f to A^t^{-1} and to invert
+    and transpose the computed circuit, thus leading to 4 generally different synthesized circuits for A,
+    allowing to pick the best one in terms of depth or the number of gates.
+upgrade:
+  - |
+    Improve code structure by moving the internal code graysynth.py from qiskit/transpiler/synthesis
+    to qiskit/synthesis/linear.
+fixes:
+  - |
+    Provide an internal utility function to generate a random NxN invertible binary matrix.
+    This function is completely robust, as opposed to the previously implemented method based on
+    computing determinant (via np.linalg.det) and checking that it's close to 1,
+    where we already saw cases of floating point errors leading to an incorrect result (for N>100).