Fix SolovayKitaev for roundoff errors in Qiskit gates (#9441) (#9444)

* fix tol and skip test for Qiskit gates

* Add regression test

Co-authored-by: Shelly Garion <46566946+ShellyGarion@users.noreply.github.com>

* add reno

* flip the check_input!

* properly include ``check_input`` in all ``_recurse``

* rm reno -- SK not released yet

Co-authored-by: Shelly Garion <46566946+ShellyGarion@users.noreply.github.com>
Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com>
(cherry picked from commit 088a84f)

Co-authored-by: Julien Gacon <gaconju@gmail.com>

qiskit/synthesis/discrete_basis/commutator_decompose.py

@@ -16,6 +16,7 @@
 
 import math
 import numpy as np
+from qiskit.quantum_info.operators.predicates import is_identity_matrix
 from .gate_sequence import _check_is_so3, GateSequence
 
 
@@ -217,10 +218,7 @@ def commutator_decompose(
         # assert that the input matrix is really SO(3)
         _check_is_so3(u_so3)
 
-        identity = np.identity(3)
-        if not (
-            np.allclose(u_so3.dot(u_so3.T), identity) and np.allclose(u_so3.T.dot(u_so3), identity)
-        ):
+        if not is_identity_matrix(u_so3.dot(u_so3.T)):
             raise ValueError("Input matrix is not orthogonal.")
 
     angle = _solve_decomposition_angle(u_so3)

qiskit/synthesis/discrete_basis/solovay_kitaev.py

@@ -86,14 +86,20 @@ def load_basic_approximations(self, data: list | str | dict) -> list[GateSequenc
         return sequences
 
     def run(
-        self, gate_matrix: np.ndarray, recursion_degree: int, return_dag: bool = False
+        self,
+        gate_matrix: np.ndarray,
+        recursion_degree: int,
+        return_dag: bool = False,
+        check_input: bool = True,
     ) -> "QuantumCircuit" | "DAGCircuit":
         r"""Run the algorithm.
 
         Args:
-            gate_matrix: The 2x2 matrix representing the gate. Does not need to be SU(2).
+            gate_matrix: The 2x2 matrix representing the gate. This matrix has to be SU(2)
+                up to global phase.
             recursion_degree: The recursion degree, called :math:`n` in the paper.
             return_dag: If ``True`` return a :class:`.DAGCircuit`, else a :class:`.QuantumCircuit`.
+            check_input: If ``True`` check that the input matrix is valid for the decomposition.
 
         Returns:
             A one-qubit circuit approximating the ``gate_matrix`` in the specified discrete basis.
@@ -104,7 +110,7 @@ def run(
         global_phase = np.arctan2(np.imag(z), np.real(z))
 
         # get the decompositon as GateSequence type
-        decomposition = self._recurse(gate_matrix_su2, recursion_degree)
+        decomposition = self._recurse(gate_matrix_su2, recursion_degree, check_input=check_input)
 
         # simplify
         _remove_identities(decomposition)
@@ -120,31 +126,35 @@ def run(
 
         return out
 
-    def _recurse(self, sequence: GateSequence, n: int) -> GateSequence:
+    def _recurse(self, sequence: GateSequence, n: int, check_input: bool = True) -> GateSequence:
         """Performs ``n`` iterations of the Solovay-Kitaev algorithm on ``sequence``.
 
         Args:
-            sequence: GateSequence to which the Solovay-Kitaev algorithm is applied.
-            n: number of iterations that the algorithm needs to run.
+            sequence: ``GateSequence`` to which the Solovay-Kitaev algorithm is applied.
+            n: The number of iterations that the algorithm needs to run.
+            check_input: If ``True`` check that the input matrix represented by ``GateSequence``
+                is valid for the decomposition.
 
         Returns:
             GateSequence that approximates ``sequence``.
 
         Raises:
-            ValueError: if ``u`` does not represent an SO(3)-matrix.
+            ValueError: If the matrix in ``GateSequence`` does not represent an SO(3)-matrix.
         """
         if sequence.product.shape != (3, 3):
             raise ValueError("Shape of U must be (3, 3) but is", sequence.shape)
 
         if n == 0:
             return self.find_basic_approximation(sequence)
 
-        u_n1 = self._recurse(sequence, n - 1)
+        u_n1 = self._recurse(sequence, n - 1, check_input=check_input)
 
-        v_n, w_n = commutator_decompose(sequence.dot(u_n1.adjoint()).product)
+        v_n, w_n = commutator_decompose(
+            sequence.dot(u_n1.adjoint()).product, check_input=check_input
+        )
 
-        v_n1 = self._recurse(v_n, n - 1)
-        w_n1 = self._recurse(w_n, n - 1)
+        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)
 
     def find_basic_approximation(self, sequence: GateSequence) -> Gate:

qiskit/transpiler/passes/synthesis/solovay_kitaev_synthesis.py

@@ -18,6 +18,7 @@
 import numpy as np
 
 from qiskit.converters import circuit_to_dag
+from qiskit.circuit.gate import Gate
 from qiskit.dagcircuit import DAGCircuit
 from qiskit.synthesis.discrete_basis.solovay_kitaev import SolovayKitaevDecomposition
 from qiskit.synthesis.discrete_basis.generate_basis_approximations import (
@@ -146,16 +147,21 @@ def run(self, dag: DAGCircuit) -> DAGCircuit:
             if not node.op.num_qubits == 1:
                 continue  # ignore all non-single qubit gates
 
+            # we do not check the input matrix as we know it comes from a Qiskit gate, as this
+            # we know it will generate a valid SU(2) matrix
+            check_input = not isinstance(node.op, Gate)
+
             if not hasattr(node.op, "to_matrix"):
                 raise TranspilerError(
-                    "SolovayKitaev does not support gate without "
-                    f"to_matrix method: {node.op.name}"
+                    f"SolovayKitaev does not support gate without to_matrix method: {node.op.name}"
                 )
 
             matrix = node.op.to_matrix()
 
             # call solovay kitaev
-            approximation = self._sk.run(matrix, self.recursion_degree, return_dag=True)
+            approximation = self._sk.run(
+                matrix, self.recursion_degree, return_dag=True, check_input=check_input
+            )
 
             # convert to a dag and replace the gate by the approximation
             dag.substitute_node_with_dag(node, approximation)

test/python/transpiler/test_solovay_kitaev.py

@@ -202,6 +202,35 @@ def test_approximation_on_qft(self):
             discretized = skd(transpiled)
             self.assertLess(_trace_distance(transpiled, discretized), 7)
 
+    def test_u_gates_work(self):
+        """Test SK works on Qiskit's UGate.
+
+        Regression test of Qiskit/qiskit-terra#9437.
+        """
+        circuit = QuantumCircuit(1)
+        circuit.u(np.pi / 2, -np.pi, -np.pi, 0)
+        circuit.u(np.pi / 2, np.pi / 2, -np.pi, 0)
+        circuit.u(-np.pi / 4, 0, -np.pi / 2, 0)
+        circuit.u(np.pi / 4, -np.pi / 16, 0, 0)
+        circuit.u(0, 0, np.pi / 16, 0)
+        circuit.u(0, np.pi / 4, np.pi / 4, 0)
+        circuit.u(np.pi / 2, 0, -15 * np.pi / 16, 0)
+        circuit.p(-np.pi / 4, 0)
+        circuit.p(np.pi / 4, 0)
+        circuit.u(np.pi / 2, 0, -3 * np.pi / 4, 0)
+        circuit.u(0, 0, -np.pi / 16, 0)
+        circuit.u(np.pi / 2, 0, 15 * np.pi / 16, 0)
+
+        depth = 4
+        basis_gates = ["h", "t", "tdg", "s", "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)
+
 
 @ddt
 class TestGateSequence(QiskitTestCase):