Fix extra CNOT inserted in Patel-Hayes-Markov algorithm (#5044)

* Fix for issue #2718

* Do inversion in graysynth rather than PMH

* Fixed tests

* Update qiskit/transpiler/synthesis/graysynth.py

* Changing some comments

* New test + out of bounds bug fix

* Typo

* linting

Co-authored-by: Luciano Bello <[email protected]>

qiskit/transpiler/synthesis/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)
+    qcir += cnot_synth(state, section_size).inverse()
     return qcir
 
 
@@ -241,6 +241,12 @@ def _lwr_cnot_synth(state, section_size):
     Patel, Ketan N., Igor L. Markov, and John P. Hayes.
     Quantum Information & Computation 8.3 (2008): 282-294.
 
+    Note:
+    This implementation tweaks the Patel, Markov, and Hayes algorithm by adding
+    a "back reduce" which adds rows below the pivot row with a high degree of
+    overlap back to it. The intuition is to avoid a high-weight pivot row
+    increasing the weight of lower rows.
+
     Args:
         state (ndarray): n x n matrix, describing a linear quantum circuit
         section_size (int): the section size the matrix columns are divided into
@@ -251,13 +257,10 @@ def _lwr_cnot_synth(state, section_size):
     """
     circuit = []
     num_qubits = state.shape[0]
+    cutoff = 1
 
-    # If the matrix is already an upper triangular one,
-    # there is no need for any transformations
-    if np.allclose(state, np.triu(state)):
-        return [state, circuit]
     # Iterate over column sections
-    for sec in range(1, int(np.ceil(num_qubits/section_size)+1)):
+    for sec in range(1, int(np.floor(num_qubits/section_size)+1)):
         # Remove duplicate sub-rows in section sec
         patt = {}
         for row in range((sec-1)*section_size, num_qubits):
@@ -284,6 +287,10 @@ def _lwr_cnot_synth(state, section_size):
                         diag_one = 1
                     state[row, :] ^= state[col, :]
                     circuit.append([col, row])
+                # Back reduce the pivot row using the current row
+                if sum(state[col, :] & state[row, :]) > cutoff:
+                    state[col, :] ^= state[row, :]
+                    circuit.append([row, col])
     return [state, circuit]
 
 

test/python/transpiler/test_synthesis.py

@@ -108,12 +108,12 @@ def test_paper_example(self):
         and only T gates as phase rotations,
 
         And should return the following circuit (or an equivalent one):
-                ┌───┐┌───┐     ┌───┐┌───┐┌───┐┌───┐┌───┐┌───┐┌───┐┌───┐┌───┐     ┌───┐
-        q_0: |0>┤ T ├┤ X ├─────┤ T ├┤ X ├┤ X ├┤ T ├┤ X ├┤ T ├┤ X ├┤ T ├┤ X ├─────┤ X ├
-                ├───┤└─┬─┘┌───┐└───┘└─┬─┘└─┬─┘└───┘└─┬─┘└───┘└─┬─┘└───┘└─┬─┘┌───┐└─┬─┘
-        q_1: |0>┤ X ├──┼──┤ T ├───────■────┼─────────┼─────────┼─────────■──┤ X ├──┼──
-                └─┬─┘  │  └───┘            │         │         │            └─┬─┘  │
-        q_2: |0>──■────┼───────────────────┼─────────■─────────┼──────────────■────┼──
+                ┌───┐┌───┐     ┌───┐┌───┐┌───┐┌───┐┌───┐┌───┐┌───┐┌───┐     ┌───┐┌───┐
+        q_0: |0>┤ T ├┤ X ├─────┤ T ├┤ X ├┤ X ├┤ T ├┤ X ├┤ T ├┤ X ├┤ T ├─────┤ X ├┤ X ├
+                ├───┤└─┬─┘┌───┐└───┘└─┬─┘└─┬─┘└───┘└─┬─┘└───┘└─┬─┘└───┘┌───┐└─┬─┘└─┬─┘
+        q_1: |0>┤ X ├──┼──┤ T ├───────■────┼─────────┼─────────┼───────┤ X ├──■────┼──
+                └─┬─┘  │  └───┘            │         │         │       └─┬─┘       │
+        q_2: |0>──■────┼───────────────────┼─────────■─────────┼─────────■─────────┼──
                        │                   │                   │                   │
         q_3: |0>───────■───────────────────■───────────────────■───────────────────■──
         """
@@ -140,14 +140,66 @@ def test_paper_example(self):
         c_compare.t(q[0])
         c_compare.cx(q[3], q[0])
         c_compare.t(q[0])
-        c_compare.cx(q[1], q[0])
         c_compare.cx(q[2], q[1])
+        c_compare.cx(q[1], q[0])
         c_compare.cx(q[3], q[0])
         unitary_compare = UnitaryGate(Operator(c_compare))
 
         # Check if the two circuits are equivalent
         self.assertEqual(unitary_gray, unitary_compare)
 
+    def test_ccz(self):
+        """Test synthesis of the doubly-controlled Z gate.
+
+        The diagonal operator in Example 4.3
+            U|x> = e^(2.pi.i.f(x))|x>,
+        where
+            f(x) = 1/8*(x0 + x1 + x2 - x0^x1 - x0^x2 - x1^x2 + x0^x1^x2)
+
+        The algorithm should take the following matrix as an input:
+        S = [[1, 0, 0, 1, 1, 0, 1],
+             [0, 1, 0, 1, 0, 1, 1],
+             [0, 0, 1, 0, 1, 1, 1]]
+
+        and only T and T* gates as phase rotations,
+
+        And should return the following circuit (or an equivalent one):
+                ┌───┐
+        q_0: |0>┤ T ├───────■──────────────■───────────────────■──────────────■──
+                └───┘┌───┐┌─┴─┐┌───┐       │                   │            ┌─┴─┐
+        q_1: |0>─────┤ T ├┤ X ├┤ T*├───────┼─────────■─────────┼─────────■──┤ X ├
+                     └───┘└───┘└───┘┌───┐┌─┴─┐┌───┐┌─┴─┐┌───┐┌─┴─┐┌───┐┌─┴─┐└───┘
+        q_2: |0>────────────────────┤ T ├┤ X ├┤ T*├┤ X ├┤ T*├┤ X ├┤ T ├┤ X ├─────
+                                    └───┘└───┘└───┘└───┘└───┘└───┘└───┘└───┘
+        """
+        cnots = [[1, 0, 0, 1, 1, 0, 1],
+                 [0, 1, 0, 1, 0, 1, 1],
+                 [0, 0, 1, 0, 1, 1, 1]]
+        angles = ['t', 't', 't', 'tdg', 'tdg', 'tdg', 't']
+        c_gray = graysynth(cnots, angles)
+        unitary_gray = UnitaryGate(Operator(c_gray))
+
+        # Create the circuit displayed above:
+        q = QuantumRegister(3, 'q')
+        c_compare = QuantumCircuit(q)
+        c_compare.t(q[0])
+        c_compare.t(q[1])
+        c_compare.cx(q[0], q[1])
+        c_compare.tdg(q[1])
+        c_compare.t(q[2])
+        c_compare.cx(q[0], q[2])
+        c_compare.tdg(q[2])
+        c_compare.cx(q[1], q[2])
+        c_compare.tdg(q[2])
+        c_compare.cx(q[0], q[2])
+        c_compare.t(q[2])
+        c_compare.cx(q[1], q[2])
+        c_compare.cx(q[0], q[1])
+        unitary_compare = UnitaryGate(Operator(c_compare))
+
+        # Check if the two circuits are equivalent
+        self.assertEqual(unitary_gray, unitary_compare)
+
 
 class TestPatelMarkovHayes(QiskitTestCase):
     """Test the Patel-Markov-Hayes algorithm for synthesizing linear
@@ -187,6 +239,7 @@ def test_patel_markov_hayes(self):
                  [1, 1, 0, 1, 1, 1],
                  [0, 0, 1, 1, 1, 0]]
         c_patel = cnot_synth(state)
+        unitary_patel = UnitaryGate(Operator(c_patel))
 
         # Create the circuit displayed above:
         q = QuantumRegister(6, 'q')
@@ -206,6 +259,7 @@ def test_patel_markov_hayes(self):
         c_compare.cx(q[0], q[1])
         c_compare.cx(q[0], q[4])
         c_compare.cx(q[0], q[3])
+        unitary_compare = UnitaryGate(Operator(c_compare))
 
         # Check if the two circuits are equivalent
-        self.assertEqual(c_patel, c_compare)
+        self.assertEqual(unitary_patel, unitary_compare)