Bug fixes in Solovay Kitaev Decomposition (backport #14217) (#14294)

* Bug fixes in Solovay Kitaev Decomposition (#14217)

* Fix _compute_rotation_axis method

The method used to return [0, 0, 0] for rotation axis when theta is very close to 0 (and yet not within 1e010).
This commit fixes this by noting that the rotation axis of an SO(3) matrix is simply the eigenvector of this
matrix corresponding to eigenvalue 1.

* Fixes related to SU(2) vs. SO(3)

The implemented algorithm has a global phase uncertainty of +-1, due to approximating not the
original SU(2) matrix but its projection onto SO(3). This fixes the global phase of the
computed approximation.

In addition, the product_su2 matrix in GateSequence was not correctly computed in many cases,
leading to incorrect values when using this attribute. Since this is not needed for the actual
algorithm (and spends unnecessary effort for computing it), this also completely removes this
attribute.

Co-authored-by: Almudena Carrera Vazquez <[email protected]>

* Fix the tests.

One of the tests was simply wrong because the reference circuit has the wrong global phase.
The tests that check the gates in the final decomposition should not assume that SGate is included and SdgGate is not.

* reno

* fix to the global phase

* bug fix (forgetting that gate_matrix_su2 is not SU(2) was a sequence)

* adding test

* Restoring the previous fast code for _compute_rotation_axis based on Rodrigues formula, while additionally handling the case of 180-degree rotation

* removing old code

* review comments

* switching back to using _to_dag and _to_circuit

I actually no longer think that these methods would ever return a circuit with unitary gates, so it does not really matter

* Addressing review comments

* pass over release notes combining the original and the suggested wordings

* Update releasenotes/notes/fix-sk-decomposition-23da3ee4b6a10d62.yaml

Co-authored-by: Julien Gacon <[email protected]>

---------

Co-authored-by: Almudena Carrera Vazquez <[email protected]>
Co-authored-by: Julien Gacon <[email protected]>
(cherry picked from commit 6892608)

# Conflicts:
#	test/python/transpiler/test_solovay_kitaev.py

* Removing tests added in PR #13690 that was not backported to 1.4

---------

Co-authored-by: Alexander Ivrii <[email protected]>
Co-authored-by: Julien Gacon <[email protected]>
Co-authored-by: Elena Peña Tapia <[email protected]>

Author: 4 people

qiskit/synthesis/discrete_basis/commutator_decompose.py

@@ -53,16 +53,40 @@ def _compute_rotation_axis(matrix: np.ndarray) -> np.ndarray:
     """
     _check_is_so3(matrix)
 
+    # If theta represents the rotation angle, then trace = 1 + 2cos(theta).
     trace = _compute_trace_so3(matrix)
-    theta = math.acos(0.5 * (trace - 1))
-    if math.sin(theta) > 1e-10:
-        x = 1 / (2 * math.sin(theta)) * (matrix[2][1] - matrix[1][2])
-        y = 1 / (2 * math.sin(theta)) * (matrix[0][2] - matrix[2][0])
-        z = 1 / (2 * math.sin(theta)) * (matrix[1][0] - matrix[0][1])
-    else:
+
+    if trace >= 3 - 1e-10:
+        # The matrix is the identity (rotation by 0)
         x = 1.0
         y = 0.0
         z = 0.0
+
+    elif trace <= -1 + 1e-10:
+        # The matrix is the 180-degree rotation
+        squares = (1 + np.diagonal(matrix)) / 2
+        index_of_max = np.argmax(squares)
+
+        if index_of_max == 0:
+            x = math.sqrt(squares[0])
+            y = matrix[0][1] / (2 * x)
+            z = matrix[0][2] / (2 * x)
+        elif index_of_max == 1:
+            y = math.sqrt(squares[1])
+            x = matrix[0][1] / (2 * y)
+            z = matrix[1][2] / (2 * y)
+        else:
+            z = math.sqrt(squares[2])
+            x = matrix[0][2] / (2 * z)
+            y = matrix[1][2] / (2 * z)
+
+    else:
+        # The matrix is the rotation by theta with sin(theta)!=0
+        theta = math.acos(0.5 * (trace - 1))
+        x = 1 / (2 * math.sin(theta)) * (matrix[2][1] - matrix[1][2])
+        y = 1 / (2 * math.sin(theta)) * (matrix[0][2] - matrix[2][0])
+        z = 1 / (2 * math.sin(theta)) * (matrix[1][0] - matrix[0][1])
+
     return np.array([x, y, z])
 
 

qiskit/synthesis/discrete_basis/gate_sequence.py

@@ -55,7 +55,6 @@ def __init__(self, gates: Sequence[Gate] = ()) -> None:
         self.name = " ".join(self.labels)
         self.global_phase = global_phase
         self.product = so3_matrix
-        self.product_su2 = su2_matrix
 
     def remove_cancelling_pair(self, indices: Sequence[int]) -> None:
         """Remove a pair of indices that cancel each other and *do not* change the matrices."""
@@ -106,6 +105,16 @@ def to_circuit(self):
 
         return circuit
 
+    def _to_u2(self):
+        """Creates the U2 matrix corresponding to the stored sequence of gates
+        and the global phase.
+        """
+        u2 = np.eye(2, dtype=complex)
+        for mat in self.gates:
+            u2 = mat.to_matrix().dot(u2)
+        u2 = np.exp(1j * self.global_phase) * u2
+        return u2
+
     def to_dag(self):
         """Convert to a :class:`.DAGCircuit`.
 
@@ -149,7 +158,6 @@ def append(self, gate: Gate) -> "GateSequence":
         so3 = _convert_su2_to_so3(su2)
 
         self.product = so3.dot(self.product)
-        self.product_su2 = su2.dot(self.product_su2)
         self.global_phase = self.global_phase + phase
 
         self.gates.append(gate)
@@ -172,7 +180,6 @@ def adjoint(self) -> "GateSequence":
         adjoint.labels = [inv.name for inv in adjoint.gates]
         adjoint.name = " ".join(adjoint.labels)
         adjoint.product = np.conj(self.product).T
-        adjoint.product_su2 = np.conj(self.product_su2).T
         adjoint.global_phase = -self.global_phase
 
         return adjoint
@@ -190,7 +197,6 @@ def copy(self) -> "GateSequence":
         out.matrices = self.matrices.copy()
         out.global_phase = self.global_phase
         out.product = self.product.copy()
-        out.product_su2 = self.product_su2.copy()
         out.name = self.name
         out._eulers = self._eulers
         return out

qiskit/synthesis/discrete_basis/generate_basis_approximations.py

@@ -76,7 +76,7 @@ def _check_candidate_greedy(candidate, existing_sequences, tol=1e-10):
         return False
 
     for existing in existing_sequences:
-        if matrix_equal(existing.product_su2, candidate.product_su2, ignore_phase=True, atol=tol):
+        if matrix_equal(existing.product, candidate.product, ignore_phase=True, atol=tol):
             # is the new sequence less or more efficient?
             return len(candidate.gates) < len(existing.gates)
     return True

qiskit/synthesis/discrete_basis/solovay_kitaev.py

@@ -24,7 +24,7 @@
 class SolovayKitaevDecomposition:
     """The Solovay Kitaev discrete decomposition algorithm.
 
-    This class is called recursively by the transpiler pass, which is why it is separeted.
+    This class is called recursively by the transpiler pass, which is why it is separated.
     See :class:`qiskit.transpiler.passes.SolovayKitaev` for more information.
     """
 
@@ -33,7 +33,7 @@ def __init__(
     ) -> None:
         """
         Args:
-            basic_approximations: A specification of the basic SU(2) approximations in terms
+            basic_approximations: A specification of the basic SO(3) approximations in terms
                 of discrete gates. At each iteration this algorithm, the remaining error is
                 approximated with the closest sequence of gates in this set.
                 If a ``str``, this specifies a ``.npy`` filename from which to load the
@@ -116,23 +116,33 @@ def run(
         """
         # make input matrix SU(2) and get the according global phase
         z = 1 / np.sqrt(np.linalg.det(gate_matrix))
-        gate_matrix_su2 = GateSequence.from_matrix(z * gate_matrix)
+
+        gate_matrix_su2 = z * gate_matrix
+        gate_matrix_as_sequence = GateSequence.from_matrix(gate_matrix_su2)
         global_phase = np.arctan2(np.imag(z), np.real(z))
 
         # get the decomposition as GateSequence type
-        decomposition = self._recurse(gate_matrix_su2, recursion_degree, check_input=check_input)
+        decomposition = self._recurse(
+            gate_matrix_as_sequence, recursion_degree, check_input=check_input
+        )
 
         # simplify
         _remove_identities(decomposition)
         _remove_inverse_follows_gate(decomposition)
 
+        # adjust to the correct SU(2) phase
+        adjust_phase = (
+            np.pi if _should_adjust_phase(decomposition._to_u2(), gate_matrix_su2) else 0.0
+        )
+
         # convert to a circuit and attach the right phases
         if return_dag:
             out = decomposition.to_dag()
         else:
             out = decomposition.to_circuit()
 
-        out.global_phase = decomposition.global_phase - global_phase
+        out.global_phase += adjust_phase
+        out.global_phase -= global_phase
 
         return out
 
@@ -155,17 +165,20 @@ def _recurse(self, sequence: GateSequence, n: int, check_input: bool = True) ->
             raise ValueError("Shape of U must be (3, 3) but is", sequence.shape)
 
         if n == 0:
-            return self.find_basic_approximation(sequence)
+            res = self.find_basic_approximation(sequence)
 
-        u_n1 = self._recurse(sequence, n - 1, check_input=check_input)
+        else:
+            u_n1 = self._recurse(sequence, n - 1, check_input=check_input)
 
-        v_n, w_n = commutator_decompose(
-            sequence.dot(u_n1.adjoint()).product, check_input=check_input
-        )
+            v_n, w_n = commutator_decompose(
+                sequence.dot(u_n1.adjoint()).product, check_input=check_input
+            )
+
+            v_n1 = self._recurse(v_n, n - 1, check_input=check_input)
+            w_n1 = self._recurse(w_n, n - 1, check_input=check_input)
+            res = v_n1.dot(w_n1).dot(v_n1.adjoint()).dot(w_n1.adjoint()).dot(u_n1)
 
-        v_n1 = self._recurse(v_n, n - 1, check_input=check_input)
-        w_n1 = self._recurse(w_n, n - 1, check_input=check_input)
-        return v_n1.dot(w_n1).dot(v_n1.adjoint()).dot(w_n1.adjoint()).dot(u_n1)
+        return res
 
     def find_basic_approximation(self, sequence: GateSequence) -> GateSequence:
         """Find ``GateSequence`` in ``self._basic_approximations`` that approximates ``sequence``.
@@ -215,3 +228,13 @@ def _remove_identities(sequence):
             sequence.gates.pop(index)
         else:
             index += 1
+
+
+def _should_adjust_phase(computed: np.ndarray, target: np.ndarray) -> bool:
+    """
+    The implemented SolovayKitaevDecomposition has a global phase uncertainty of +-1,
+    due to approximating not the original SU(2) matrix but its projection onto SO(3).
+    This function returns ``True`` if the global phase of the computed approximation
+    should be adjusted (by adding pi) to better much the target.
+    """
+    return np.linalg.norm(-computed - target) < np.linalg.norm(computed - target)

releasenotes/notes/fix-sk-decomposition-23da3ee4b6a10d62.yaml

@@ -0,0 +1,12 @@
+---
+fixes:
+  - |
+    Fixed a problem in the :class:`.SolovayKitaev` transpiler pass where the pass could
+    crash due to encountering a 180 degree rotation in the internal recursion,
+    which was not handled correctly.
+  - |
+    Fixed a problem in the :class:`.SolovayKitaev` transpiler pass where the generated
+    approximation could have a phase that differs by :math:`\pi` from the correct value.
+    This resulted due to the internal :math:`SO(3)` representation, which requires additional
+    handling to obtain the correct sign of the qubit gate matrix.
+    Fixed `#9552 <https://github.com/Qiskit/qiskit-terra/issues/9552>`__

test/python/transpiler/test_solovay_kitaev.py

@@ -98,9 +98,7 @@ def test_unitary_synthesis(self):
         passes = PassManager([_1q, _cons, _synth])
         compiled = passes.run(circuit)
 
-        diff = np.linalg.norm(Operator(compiled) - Operator(circuit))
-        self.assertLess(diff, 1)
-        self.assertEqual(set(compiled.count_ops().keys()), {"h", "s", "cx"})
+        self.assertLessEqual(set(compiled.count_ops().keys()), {"h", "s", "sdg", "cx"})
 
     def test_plugin(self):
         """Test calling the plugin directly."""
@@ -112,12 +110,7 @@ def test_plugin(self):
         plugin = SolovayKitaevSynthesis()
         out = plugin.run(unitary, basis_gates=["h", "s"])
 
-        reference = QuantumCircuit(1, global_phase=3 * np.pi / 4)
-        reference.h(0)
-        reference.s(0)
-        reference.h(0)
-
-        self.assertEqual(dag_to_circuit(out), reference)
+        self.assertLessEqual(set(out.count_ops().keys()), {"h", "s", "sdg", "cx"})
 
     def test_multiple_plugins(self):
         """Test calling the plugins directly but with different instances of basis set."""
@@ -199,11 +192,14 @@ def test_str_basis_gates(self):
         dag = circuit_to_dag(circuit)
         discretized = dag_to_circuit(synth.run(dag))
 
-        reference = QuantumCircuit(1, global_phase=7 * np.pi / 8)
+        reference = QuantumCircuit(1, global_phase=15 * np.pi / 8)
         reference.h(0)
         reference.t(0)
         reference.h(0)
 
+        # Make sure that the discretized circuit gives a valid approximation
+        diff = _trace_distance(circuit, discretized)
+        self.assertLess(diff, 0.01)
         self.assertEqual(discretized, reference)
 
     def test_approximation_on_qft(self):
@@ -242,14 +238,14 @@ def test_u_gates_work(self):
         circuit.u(np.pi / 2, 0, 15 * np.pi / 16, 0)
 
         depth = 4
-        basis_gates = ["h", "t", "tdg", "s", "z"]
+        basis_gates = ["h", "t", "tdg", "s", "sdg", "z"]
         gate_approx_library = generate_basic_approximations(basis_gates=basis_gates, depth=depth)
 
         skd = SolovayKitaev(recursion_degree=2, basic_approximations=gate_approx_library)
         discretized = skd(circuit)
 
         included_gates = set(discretized.count_ops().keys())
-        self.assertEqual(set(basis_gates), included_gates)
+        self.assertLessEqual(included_gates, set(basis_gates))
 
     def test_load_from_file(self):
         """Test loading basic approximations from a file works.
@@ -262,24 +258,32 @@ def test_load_from_file(self):
             fullpath = os.path.join(tmp_dir, filename)
 
             # dump approximations to file
-            generate_basic_approximations(basis_gates=["h", "s", "sdg"], depth=3, filename=fullpath)
+            gate_approx_library = generate_basic_approximations(
+                basis_gates=["h", "s", "sdg"], depth=3, filename=fullpath
+            )
 
             # circuit to decompose and reference decomp
             circuit = QuantumCircuit(1)
             circuit.rx(0.8, 0)
 
-            reference = QuantumCircuit(1, global_phase=3 * np.pi / 4)
-            reference.h(0)
-            reference.s(0)
-            reference.h(0)
+            # Run SK pass using gate_approx_library
+            reference = SolovayKitaev(basic_approximations=gate_approx_library)(circuit)
 
-            # load the decomp and compare to reference
-            skd = SolovayKitaev(basic_approximations=fullpath)
-            # skd = SolovayKitaev(basic_approximations=filename)
-            discretized = skd(circuit)
+            # Run SK pass using stored basis_approximations
+            discretized = SolovayKitaev(basic_approximations=fullpath)(circuit)
 
+        # Check that both flows produce the same result
         self.assertEqual(discretized, reference)
 
+    def test_y_gate(self):
+        """Test the Solovay-Kitaev decomposition on the circuit with a Y-gate (see issue #9552)."""
+        circuit = QuantumCircuit(1)
+        circuit.y(0)
+
+        transpiled = SolovayKitaev()(circuit)
+        diff = _trace_distance(circuit, transpiled)
+        self.assertLess(diff, 1e-6)
+
 
 @ddt
 class TestGateSequence(QiskitTestCase):