Repository: Qiskit/qiskit-terra (python)

Issue:
Solovay-Kitaev, ValueError: Input matrix is not orthogonal.
### Environment

- **Commit hash**:  dc8e232a
- **Python version**: 3.8.7
- **Operating system**: Ubuntu 20.04.5 LTS
- **Qiskit Terra version**: 0.24.0

### What is happening?

I'm trying to run the newly added SolovayKitaev (SK) algorithm. For some simple qiskit circuits the SK works, but most of the time it produces the following error:

```
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[13], line 15
     12 # Time computation
     13 start_time = time.time()
---> 15 discretized = skd(transpiled_circ_qasm)
     17 print("\t Time of conversion:", (time.time() - start_time) / 60, " min")

File ~/.pyenv/versions/circuit_comp_qiskit_dev/lib/python3.8/site-packages/qiskit/transpiler/basepasses.py:121, in BasePass.__call__(self, circuit, property_set)
    118 if isinstance(property_set_, PropertySet):
    119     self.property_set = property_set_
--> 121 result = self.run(circuit_to_dag(circuit))
    123 result_circuit = circuit
    125 if isinstance(property_set, dict):  # this includes (dict, PropertySet)

File ~/.pyenv/versions/circuit_comp_qiskit_dev/lib/python3.8/site-packages/qiskit/transpiler/passes/synthesis/solovay_kitaev_synthesis.py:158, in SolovayKitaev.run(self, dag)
    155 matrix = node.op.to_matrix()
    157 # call solovay kitaev
--> 158 approximation = self._sk.run(matrix, self.recursion_degree, return_dag=True)
    160 # convert to a dag and replace the gate by the approximation
    161 dag.substitute_node_with_dag(node, approximation)

File ~/.pyenv/versions/circuit_comp_qiskit_dev/lib/python3.8/site-packages/qiskit/synthesis/discrete_basis/solovay_kitaev.py:107, in SolovayKitaevDecomposition.run(self, gate_matrix, recursion_degree, return_dag)
    104 global_phase = np.arctan2(np.imag(z), np.real(z))
    106 # get the decompositon as GateSequence type
--> 107 decomposition = self._recurse(gate_matrix_su2, recursion_degree)
    109 # simplify
    110 _remove_identities(decomposition)

File ~/.pyenv/versions/circuit_comp_qiskit_dev/lib/python3.8/site-packages/qiskit/synthesis/discrete_basis/solovay_kitaev.py:146, in SolovayKitaevDecomposition._recurse(self, sequence, n)
    142 u_n1 = self._recurse(sequence, n - 1)
    144 v_n, w_n = commutator_decompose(sequence.dot(u_n1.adjoint()).product)
--> 146 v_n1 = self._recurse(v_n, n - 1)
    147 w_n1 = self._recurse(w_n, n - 1)
    148 return v_n1.dot(w_n1).dot(v_n1.adjoint()).dot(w_n1.adjoint()).dot(u_n1)

File ~/.pyenv/versions/circuit_comp_qiskit_dev/lib/python3.8/site-packages/qiskit/synthesis/discrete_basis/solovay_kitaev.py:144, in SolovayKitaevDecomposition._recurse(self, sequence, n)
    140     return self.find_basic_approximation(sequence)
    142 u_n1 = self._recurse(sequence, n - 1)
--> 144 v_n, w_n = commutator_decompose(sequence.dot(u_n1.adjoint()).product)
    146 v_n1 = self._recurse(v_n, n - 1)
    147 w_n1 = self._recurse(w_n, n - 1)

File ~/.pyenv/versions/circuit_comp_qiskit_dev/lib/python3.8/site-packages/qiskit/synthesis/discrete_basis/commutator_decompose.py:224, in commutator_decompose(u_so3, check_input)
    220     identity = np.identity(3)
    221     if not (
    222         np.allclose(u_so3.dot(u_so3.T), identity) and np.allclose(u_so3.T.dot(u_so3), identity)
    223     ):
--> 224         raise ValueError("Input matrix is not orthogonal.")
    226 angle = _solve_decomposition_angle(u_so3)
    228 # Compute rotation about x-axis with angle 'angle'

ValueError: Input matrix is not orthogonal.
```

### How can we reproduce the issue?

Since I can't share the qiskit circuit, I will upload its QASM file which contains about 100 (out of few millions) gates  of the original circuit. This QASM file is enough to replicate the error.

First, I will give the code:
```
# Load the QASM file
transpiled_circ_qasm = QuantumCircuit.from_qasm_file(filename)

# Do SK pass
depth = 4
gate_approx_library = generate_basic_approximations(
    basis_gates=["h", "t", "tdg", "s", "z"], depth=depth
)

skd = SolovayKitaev(recursion_degree=2, basic_approximations=gate_approx_library)
discretized = skd(transpiled_circ_qasm)

``` 

Here is the partial QASM file generated from  a qiskit circuit:
```
OPENQASM 2.0;
include "qelib1.inc";
qreg eigenvalue_ancilla_register[1];
qreg eigenvalue_register[4];
qreg coin_register[3];
qreg coin_ancilla_register[1];
qreg power_register[8];
qreg state_register[8];
creg classical_register[8];
u3(pi/2,-pi,-pi) eigenvalue_ancilla_register[0];
u3(pi/2,pi/2,-pi) eigenvalue_register[0];
u3(pi/2,pi/2,-pi) eigenvalue_register[1];
u3(pi/2,pi/2,-pi) eigenvalue_register[2];
u3(pi/2,pi/2,-pi) eigenvalue_register[3];
u3(pi/2,pi/2,-pi) coin_register[0];
cx eigenvalue_register[0],coin_register[0];
u3(-pi/4,0,-pi/2) coin_register[0];
cx eigenvalue_register[0],coin_register[0];
u3(pi/4,-pi/16,0) coin_register[0];
u3(0,0,pi/16) eigenvalue_register[0];
cx eigenvalue_register[0],coin_register[0];
u3(0,0,pi/16) coin_register[0];
u3(0,pi/4,pi/4) eigenvalue_register[0];
u3(pi/2,pi/2,-pi) coin_register[1];
cx eigenvalue_register[0],coin_register[1];
u3(-pi/4,0,-pi/2) coin_register[1];
cx eigenvalue_register[0],coin_register[1];
u3(pi/4,-pi/16,0) coin_register[1];
u3(0,0,pi/16) eigenvalue_register[0];
cx eigenvalue_register[0],coin_register[1];
u3(pi/2,0,-15*pi/16) coin_register[1];
cx coin_register[0],coin_register[1];
u1(-pi/4) coin_register[1];
cx eigenvalue_register[0],coin_register[1];
u1(pi/4) coin_register[1];
cx coin_register[0],coin_register[1];
u1(pi/4) coin_register[0];
u1(-pi/4) coin_register[1];
cx eigenvalue_register[0],coin_register[1];
u3(pi/2,0,-3*pi/4) coin_register[1];
cx eigenvalue_register[0],coin_register[0];
u1(-pi/4) coin_register[0];
cx eigenvalue_register[0],coin_register[0];
u3(0,0,-pi/16) eigenvalue_register[0];
cx eigenvalue_register[0],coin_register[1];
u3(0,0,pi/16) coin_register[1];
cx eigenvalue_register[0],coin_register[1];
u3(pi/2,0,15*pi/16) coin_register[1];
cx coin_register[0],coin_register[1];
u1(-pi/4) coin_register[1];
cx eigenvalue_register[0],coin_register[1];
u1(pi/4) coin_register[1];
cx coin_register[0],coin_register[1];
u1(pi/4) coin_register[0];
u1(-pi/4) coin_register[1];
cx eigenvalue_register[0],coin_register[1];
u3(pi/2,0,-3*pi/4) coin_register[1];
cx eigenvalue_register[0],coin_register[0];
u1(-pi/4) coin_register[0];
cx eigenvalue_register[0],coin_register[0];
u3(0,pi/4,pi/4) eigenvalue_register[0];
u3(pi/2,pi/2,-pi) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u3(-pi/4,0,-pi/2) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u3(pi/4,-pi/16,0) coin_register[2];
u3(0,0,pi/16) eigenvalue_register[0];
cx eigenvalue_register[0],coin_register[2];
u3(pi/2,0,-15*pi/16) coin_register[2];
cx coin_register[1],coin_register[2];
u1(-pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u1(pi/4) coin_register[2];
cx coin_register[1],coin_register[2];
u1(pi/4) coin_register[1];
u1(-pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u3(pi/2,0,-3*pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[1];
u1(-pi/4) coin_register[1];
cx eigenvalue_register[0],coin_register[1];
u3(0,0,-pi/16) eigenvalue_register[0];
cx eigenvalue_register[0],coin_register[2];
u3(0,0,pi/16) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u3(pi/2,0,15*pi/16) coin_register[2];
cx coin_register[0],coin_register[2];
u1(-pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u1(pi/4) coin_register[2];
cx coin_register[0],coin_register[2];
u1(pi/4) coin_register[0];
u1(-pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u3(pi/2,0,-3*pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[0];
u1(-pi/4) coin_register[0];
cx eigenvalue_register[0],coin_register[0];
u3(0,0,pi/16) eigenvalue_register[0];
cx eigenvalue_register[0],coin_register[2];
u3(0,0,-pi/16) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u3(pi/2,0,-15*pi/16) coin_register[2];
cx coin_register[1],coin_register[2];
u1(-pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u1(pi/4) coin_register[2];
cx coin_register[1],coin_register[2];
u1(pi/4) coin_register[1];
u1(-pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u3(pi/2,0,-3*pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[1];
u1(-pi/4) coin_register[1];
cx eigenvalue_register[0],coin_register[1];
u3(0,0,-pi/16) eigenvalue_register[0];
cx eigenvalue_register[0],coin_register[2];
u3(0,0,pi/16) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u3(pi/2,0,15*pi/16) coin_register[2];
cx coin_register[0],coin_register[2];
u1(-pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u1(pi/4) coin_register[2];
cx coin_register[0],coin_register[2];
u1(pi/4) coin_register[0];
u1(-pi/4) coin_register[2];
cx eigenvalue_register[0],coin_register[2];
u3(pi/2,0,-3*pi/4) coin_register[2];
```

### What should happen?

I would like to emphasize that even if I feed the circuit directly into SK (without loading it as QASM), the SK algorithm still produces the error. The original circuit was created exclusively with qiskit's default toolset and it runs on qiskit simulators without problems.

The expectation is that the following lines of code should produce a new quantum circuit which is comprised of $S, T, T^*, Z, CNOT$ and $H$ gates. Instead there is an error `ValueError: Input matrix is not orthogonal.`

```
skd = SolovayKitaev(recursion_degree=2, basic_approximations=gate_approx_library)
discretized = skd(transpiled_circ_qasm)
```

### Any suggestions?

I can only kindly ask to fix this issue ASAP.

Pull request description:
Fix `SolovayKitaev` for roundoff errors in Qiskit gates
<!--
⚠️ If you do not respect this template, your pull request will be closed.
⚠️ Your pull request title should be short detailed and understandable for all.
⚠️ Also, please add a release note file using reno if the change needs to be
  documented in the release notes.
⚠️ If your pull request fixes an open issue, please link to the issue.

- [ ] I have added the tests to cover my changes.
- [ ] I have updated the documentation accordingly.
- [ ] I have read the CONTRIBUTING document.
-->

### Summary

Fixes #9437.




Patch:
diff --git a/qiskit/synthesis/discrete_basis/commutator_decompose.py b/qiskit/synthesis/discrete_basis/commutator_decompose.py
index c50a19dfa..10622e76e 100644
--- a/qiskit/synthesis/discrete_basis/commutator_decompose.py
+++ b/qiskit/synthesis/discrete_basis/commutator_decompose.py
@@ -16,6 +16,7 @@ from __future__ import annotations
 
 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)
diff --git a/qiskit/synthesis/discrete_basis/solovay_kitaev.py b/qiskit/synthesis/discrete_basis/solovay_kitaev.py
index cd470265f..62ad50582 100644
--- a/qiskit/synthesis/discrete_basis/solovay_kitaev.py
+++ b/qiskit/synthesis/discrete_basis/solovay_kitaev.py
@@ -86,14 +86,20 @@ class SolovayKitaevDecomposition:
         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 @@ class SolovayKitaevDecomposition:
         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,18 +126,20 @@ class SolovayKitaevDecomposition:
 
         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)
@@ -139,12 +147,14 @@ class SolovayKitaevDecomposition:
         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:
diff --git a/qiskit/transpiler/passes/synthesis/solovay_kitaev_synthesis.py b/qiskit/transpiler/passes/synthesis/solovay_kitaev_synthesis.py
index feb40ad53..43ae079b9 100644
--- a/qiskit/transpiler/passes/synthesis/solovay_kitaev_synthesis.py
+++ b/qiskit/transpiler/passes/synthesis/solovay_kitaev_synthesis.py
@@ -18,6 +18,7 @@ from __future__ import annotations
 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 @@ class SolovayKitaev(TransformationPass):
             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 patch:
diff --git a/test/python/transpiler/test_solovay_kitaev.py b/test/python/transpiler/test_solovay_kitaev.py
index 95dfd0676..d3d85b8cf 100644
--- a/test/python/transpiler/test_solovay_kitaev.py
+++ b/test/python/transpiler/test_solovay_kitaev.py
@@ -202,6 +202,35 @@ class TestSolovayKitaev(QiskitTestCase):
             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):