Improved Clifford synthesis method (Qiskit#5779)

* add clifford greedy compiler * add documentation and handle parameter names * add spaces * remove unnecessary imports * updates following review * minor edits * add release notes * empty commit to rerun CI * empty commit to rerun CI * Update qiskit/quantum_info/synthesis/clifford_decompose.py Co-authored-by: Luciano Bello <bel@zurich.ibm.com> Co-authored-by: Christopher J. Wood <cjwood@us.ibm.com>
lcapelluto · Mar 5, 2021 · b549e01 · b549e01
1 parent 7c7a956
commit b549e01
Show file tree

Hide file tree

Showing 3 changed files with 335 additions and 6 deletions.
diff --git a/qiskit/quantum_info/synthesis/clifford_decompose.py b/qiskit/quantum_info/synthesis/clifford_decompose.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2017, 2020
+# (C) Copyright IBM 2017, 2021
 #
 # 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
@@ -19,18 +19,22 @@
 from qiskit.exceptions import QiskitError
 from qiskit.circuit import QuantumCircuit
 from qiskit.quantum_info.operators.symplectic.clifford_circuits import (
-    _append_z, _append_x, _append_h, _append_s, _append_v, _append_w, _append_cx, _append_swap)
+    _append_z, _append_x, _append_h, _append_s, _append_v, _append_w,
+    _append_cx, _append_swap)
 
 
-def decompose_clifford(clifford):
+def decompose_clifford(clifford, method=None):
     """Decompose a Clifford operator into a QuantumCircuit.
 
     For N <= 3 qubits this is based on optimal CX cost decomposition
     from reference [1]. For N > 3 qubits this is done using the general
-    non-optimal compilation routine from reference [2].
+    non-optimal greedy compilation routine from reference [3],
+    which typically yields better CX cost compared to the AG method in [2].
 
     Args:
         clifford (Clifford): a clifford operator.
+        method (str):  Optional, a synthesis method ('AG' or 'greedy').
+             If set this overrides optimal decomposition for N <=3 qubits.
 
     Return:
         QuantumCircuit: a circuit implementation of the Clifford.
@@ -43,13 +47,22 @@ def decompose_clifford(clifford):
         2. S. Aaronson, D. Gottesman, *Improved Simulation of Stabilizer Circuits*,
            Phys. Rev. A 70, 052328 (2004).
            `arXiv:quant-ph/0406196 <https://arxiv.org/abs/quant-ph/0406196>`_
+
+        3. Sergey Bravyi, Shaohan Hu, Dmitri Maslov, Ruslan Shaydulin,
+           *Clifford Circuit Optimization with Templates and Symbolic Pauli Gates*
     """
     num_qubits = clifford.num_qubits
 
+    if method == 'AG':
+        return decompose_clifford_ag(clifford)
+
+    if method == 'greedy':
+        return decompose_clifford_greedy(clifford)
+
     if num_qubits <= 3:
         return decompose_clifford_bm(clifford)
 
-    return decompose_clifford_ag(clifford)
+    return decompose_clifford_greedy(clifford)
 
 
 # ---------------------------------------------------------------------
@@ -413,3 +426,299 @@ def _set_row_z_zero(clifford, circuit, qubit):
             circuit.s(qubit)
         _append_h(clifford, qubit)
         circuit.h(qubit)
+
+
+# ---------------------------------------------------------------------
+# Synthesis based on Bravyi et. al. greedy clifford compiler
+# ---------------------------------------------------------------------
+
+def decompose_clifford_greedy(clifford):
+    """Decompose a Clifford operator into a QuantumCircuit.
+
+    Args:
+        clifford (Clifford): a clifford operator.
+
+    Return:
+        QuantumCircuit: a circuit implementation of the Clifford.
+
+    Raises:
+        QiskitError: if symplectic Gaussian elimination fails.
+    """
+
+    num_qubits = clifford.num_qubits
+    circ = QuantumCircuit(num_qubits, name=str(clifford))
+    qubit_list = list(range(num_qubits))
+    clifford_cpy = clifford.copy()
+
+    # Reducing the original Clifford to identity
+    # via symplectic Gaussian elimination
+    while len(qubit_list) > 0:
+        clifford_cpy_inv = clifford_cpy.adjoint()
+
+        list_greedy_cost = []
+        for qubit in qubit_list:
+            cliff_ox = clifford_cpy.copy()
+            _append_x(cliff_ox, qubit)
+            cliff_ox = cliff_ox.compose(clifford_cpy_inv)
+
+            cliff_oz = clifford_cpy.copy()
+            _append_z(cliff_oz, qubit)
+            cliff_oz = cliff_oz.compose(clifford_cpy_inv)
+
+            list_pairs = []
+            pauli_count = 0
+
+            # Compute the CNOT cost in order to find the qubit with the minimal cost
+            for i in qubit_list:
+                typeq = _from_pair_cliffs_to_type(cliff_ox, cliff_oz, i)
+                list_pairs.append(typeq)
+                pauli_count += 1
+            cost = _compute_greedy_cost(list_pairs)
+            list_greedy_cost.append([cost, qubit])
+
+        _, min_qubit = (sorted(list_greedy_cost))[0]
+
+        # Gaussian elimination step for the qubit with minimal CNOT cost
+        cliff_ox = clifford_cpy.copy()
+        _append_x(cliff_ox, min_qubit)
+        cliff_ox = cliff_ox.compose(clifford_cpy_inv)
+
+        cliff_oz = clifford_cpy.copy()
+        _append_z(cliff_oz, min_qubit)
+        cliff_oz = cliff_oz.compose(clifford_cpy_inv)
+
+        # Compute the decoupling operator of cliff_ox and cliff_oz
+        decouple_circ, decouple_cliff = _calc_decoupling(cliff_ox, cliff_oz, qubit_list,
+                                                         min_qubit, num_qubits)
+        circ = circ.compose(decouple_circ)
+
+        # Now the clifford acts trivially on min_qubit
+        clifford_cpy = decouple_cliff.adjoint().compose(clifford_cpy)
+        qubit_list.remove(min_qubit)
+
+    # Add the phases (Pauli gates) to the Clifford circuit
+    for qubit in range(num_qubits):
+        stab = clifford_cpy.stabilizer.phase[qubit]
+        destab = clifford_cpy.destabilizer.phase[qubit]
+        if destab and stab:
+            circ.y(qubit)
+        elif not destab and stab:
+            circ.x(qubit)
+        elif destab and not stab:
+            circ.z(qubit)
+
+    return circ
+
+
+# ---------------------------------------------------------------------
+# Helper functions for Bravyi et. al. greedy clifford compiler
+# ---------------------------------------------------------------------
+
+# Global arrays of the 16 pairs of Pauli operators
+# divided into 5 equivalence classes under the action of single-qubit Cliffords
+
+# Class A - canonical representative is 'XZ'
+A_class = [[[False, True], [True, True]],  # 'XY'
+           [[False, True], [True, False]],  # 'XZ'
+           [[True, True], [False, True]],  # 'YX'
+           [[True, True], [True, False]],  # 'YZ'
+           [[True, False], [False, True]],  # 'ZX'
+           [[True, False], [True, True]]]  # 'ZY'
+
+# Class B - canonical representative is 'XX'
+B_class = [[[True, False], [True, False]],  # 'ZZ'
+           [[False, True], [False, True]],  # 'XX'
+           [[True, True], [True, True]]]  # 'YY'
+
+# Class C - canonical representative is 'XI'
+C_class = [[[True, False], [False, False]],  # 'ZI'
+           [[False, True], [False, False]],  # 'XI'
+           [[True, True], [False, False]]]  # 'YI'
+
+# Class D - canonical representative is 'IZ'
+D_class = [[[False, False], [False, True]],  # 'IX'
+           [[False, False], [True, False]],  # 'IZ'
+           [[False, False], [True, True]]]  # 'IY'
+
+# Class E - only 'II'
+E_class = [[[False, False], [False, False]]]  # 'II'
+
+
+def _from_pair_cliffs_to_type(cliff_ox, cliff_oz, qubit):
+    """Converts a pair of Paulis Ox and Oz into a type"""
+
+    type_ox = [cliff_ox.destabilizer.phase[qubit], cliff_ox.stabilizer.phase[qubit]]
+    type_oz = [cliff_oz.destabilizer.phase[qubit], cliff_oz.stabilizer.phase[qubit]]
+    return [type_ox, type_oz]
+
+
+def _compute_greedy_cost(list_pairs):
+    """Compute the CNOT cost of one step of the algorithm"""
+
+    A_num = 0
+    B_num = 0
+    C_num = 0
+    D_num = 0
+
+    for pair in list_pairs:
+        if pair in A_class:
+            A_num += 1
+        elif pair in B_class:
+            B_num += 1
+        elif pair in C_class:
+            C_num += 1
+        elif pair in D_class:
+            D_num += 1
+
+    if (A_num % 2) == 0:
+        raise QiskitError("Symplectic Gaussian elimination fails.")
+
+    # Calculate the CNOT cost
+    cost = 3 * (A_num - 1) / 2 + (B_num + 1) * (B_num > 0) + C_num + D_num
+    if list_pairs[0] not in A_class:  # additional SWAP
+        cost += 3
+
+    return cost
+
+
+def _calc_decoupling(cliff_ox, cliff_oz, qubit_list, min_qubit, num_qubits):
+    """Calculate a decoupling operator D:
+    D^{-1} * Ox * D = x1
+    D^{-1} * Oz * D = z1
+    and reduce the clifford such that it will act trivially on min_qubit
+    """
+
+    circ = QuantumCircuit(num_qubits)
+
+    # decouple_cliff is initialized to an identity clifford
+    decouple_cliff = cliff_ox.copy()
+    num_qubits = decouple_cliff.num_qubits
+    decouple_cliff.table.phase = np.zeros(2 * num_qubits)
+    if (decouple_cliff.table.array != np.eye(2 * num_qubits)).any():
+        raise QiskitError("Symplectic Gaussian elimination fails.")
+
+    qubit0 = min_qubit  # The qubit for the symplectic Gaussian elimination
+
+    # Reduce the pair of Paulis to a representative in the equivalence class
+    # ['XZ', 'XX', 'XI', 'IZ', 'II'] by adding single-qubit gates
+    for qubit in qubit_list:
+
+        typeq = _from_pair_cliffs_to_type(cliff_ox, cliff_oz, qubit)
+
+        if typeq in [[[True, True], [False, False]],  # 'YI'
+                     [[True, True], [True, True]],  # 'YY'
+                     [[True, True], [True, False]]]:  # 'YZ':
+            circ.s(qubit)
+            _append_s(decouple_cliff, qubit)
+
+        elif typeq in [[[True, False], [False, False]],  # 'ZI'
+                       [[True, False], [True, False]],   # 'ZZ'
+                       [[True, False], [False, True]],  # 'ZX'
+                       [[False, False], [False, True]]]:  # 'IX'
+            circ.h(qubit)
+            _append_h(decouple_cliff, qubit)
+
+        elif typeq in [[[False, False], [True, True]],  # 'IY'
+                       [[True, False], [True, True]]]:  # 'ZY'
+            circ.s(qubit)
+            circ.h(qubit)
+            _append_s(decouple_cliff, qubit)
+            _append_h(decouple_cliff, qubit)
+
+        elif typeq == [[True, True], [False, True]]:  # 'YX'
+            circ.h(qubit)
+            circ.s(qubit)
+            _append_h(decouple_cliff, qubit)
+            _append_s(decouple_cliff, qubit)
+
+        elif typeq == [[False, True], [True, True]]:  # 'XY'
+            circ.s(qubit)
+            circ.h(qubit)
+            circ.s(qubit)
+            _append_s(decouple_cliff, qubit)
+            _append_h(decouple_cliff, qubit)
+            _append_s(decouple_cliff, qubit)
+
+    # Reducing each pair of Paulis (except of qubit0) to 'II'
+    # by adding two-qubit gates and single-qubit gates
+    A_qubits = []
+    B_qubits = []
+    C_qubits = []
+    D_qubits = []
+
+    for qubit in qubit_list:
+        typeq = _from_pair_cliffs_to_type(cliff_ox, cliff_oz, qubit)
+        if typeq in A_class:
+            A_qubits.append(qubit)
+        elif typeq in B_class:
+            B_qubits.append(qubit)
+        elif typeq in C_class:
+            C_qubits.append(qubit)
+        elif typeq in D_class:
+            D_qubits.append(qubit)
+
+    if len(A_qubits) % 2 != 1:
+        raise QiskitError("Symplectic Gaussian elimination fails.")
+
+    if qubit0 not in A_qubits:  # SWAP qubit0 and qubitA
+        qubitA = A_qubits[0]
+        circ.swap(qubit0, qubitA)
+        _append_swap(decouple_cliff, qubit0, qubitA)
+        if qubit0 in B_qubits:
+            B_qubits.remove(qubit0)
+            B_qubits.append(qubitA)
+            A_qubits.remove(qubitA)
+            A_qubits.append(qubit0)
+        elif qubit0 in C_qubits:
+            C_qubits.remove(qubit0)
+            C_qubits.append(qubitA)
+            A_qubits.remove(qubitA)
+            A_qubits.append(qubit0)
+        elif qubit0 in D_qubits:
+            D_qubits.remove(qubit0)
+            D_qubits.append(qubitA)
+            A_qubits.remove(qubitA)
+            A_qubits.append(qubit0)
+        else:
+            A_qubits.remove(qubitA)
+            A_qubits.append(qubit0)
+
+    # Reduce pairs in Class C to 'II'
+    for qubit in C_qubits:
+        circ.cx(qubit0, qubit)
+        _append_cx(decouple_cliff, qubit0, qubit)
+
+    # Reduce pairs in Class D to 'II'
+    for qubit in D_qubits:
+        circ.cx(qubit, qubit0)
+        _append_cx(decouple_cliff, qubit, qubit0)
+
+    # Reduce pairs in Class B to 'II'
+    if len(B_qubits) > 1:
+        for qubit in B_qubits[1:]:
+            qubitB = B_qubits[0]
+            circ.cx(qubitB, qubit)
+            _append_cx(decouple_cliff, qubitB, qubit)
+
+    if len(B_qubits) > 0:
+        qubitB = B_qubits[0]
+        circ.cx(qubit0, qubitB)
+        circ.h(qubitB)
+        circ.cx(qubitB, qubit0)
+        _append_cx(decouple_cliff, qubit0, qubitB)
+        _append_h(decouple_cliff, qubitB)
+        _append_cx(decouple_cliff, qubitB, qubit0)
+
+    # Reduce pairs in Class A (except of qubit0) to 'II'
+    Alen = int((len(A_qubits) - 1) / 2)
+    if Alen > 0:
+        A_qubits.remove(qubit0)
+    for qubit in range(Alen):
+        circ.cx(A_qubits[2 * qubit + 1], A_qubits[2 * qubit])
+        circ.cx(A_qubits[2 * qubit], qubit0)
+        circ.cx(qubit0, A_qubits[2 * qubit + 1])
+        _append_cx(decouple_cliff, A_qubits[2 * qubit + 1], A_qubits[2 * qubit])
+        _append_cx(decouple_cliff, A_qubits[2 * qubit], qubit0)
+        _append_cx(decouple_cliff, qubit0, A_qubits[2 * qubit + 1])
+
+    return circ, decouple_cliff
diff --git a/releasenotes/notes/improve-clifford-synthesis-6bddeefa78a2cac1.yaml b/releasenotes/notes/improve-clifford-synthesis-6bddeefa78a2cac1.yaml
@@ -0,0 +1,9 @@
+---
+features:
+  - |
+    The :meth:`~Clifford.to_circuit` method now uses a non-optimal
+    greedy compilation routine for Clifford elements synthesis, by Bravyi et. al.,
+    which typically yields better CX cost compared to the previousely
+    used Aaronson-Gottesman method (for more than two qubits).
+    This is now the default method for Clifford circuit synthesis of more
+    than three qubits, unless another method is specified.
diff --git a/test/python/quantum_info/operators/symplectic/test_clifford.py b/test/python/quantum_info/operators/symplectic/test_clifford.py
@@ -28,7 +28,7 @@
 from qiskit.quantum_info.operators import Clifford, Operator
 from qiskit.quantum_info.operators.symplectic.clifford_circuits import _append_circuit
 from qiskit.quantum_info.synthesis.clifford_decompose import (
-    decompose_clifford_ag, decompose_clifford_bm)
+    decompose_clifford_ag, decompose_clifford_bm, decompose_clifford_greedy)
 
 
 class VGate(Gate):
@@ -413,6 +413,17 @@ def test_decompose_2q_ag(self, num_qubits):
             value = Clifford(decompose_clifford_ag(target))
             self.assertEqual(value, target)
 
+    @combine(num_qubits=[1, 2, 3, 4, 5])
+    def test_decompose_2q_greedy(self, num_qubits):
+        """Test greedy synthesis for set of {num_qubits}-qubit Cliffords"""
+        rng = np.random.default_rng(1234)
+        samples = 50
+        for _ in range(samples):
+            circ = random_clifford_circuit(num_qubits, 5 * num_qubits, seed=rng)
+            target = Clifford(circ)
+            value = Clifford(decompose_clifford_greedy(target))
+            self.assertEqual(value, target)
+
 
 @ddt
 class TestCliffordDecomposition(QiskitTestCase):