Repository: PennyLaneAI/pennylane
Pull request #4801 (2023-11-07 22:26:05)

Solovay-Kitaev Decomposition Algorithm
**Context:** Implements the `Solovay-Kitaev` decomposition algorithm using the description mentioned in [arXiv:0505030](https://arxiv.org/abs/quant-ph/0505030)

**Description of the Change:** This adds a Solovay-Kitaev decomposition for the Clifford+T workflow in the `solovay_kitaev.py` file. 

1. `sk_decomposition` technique is the parent method, which takes in any single-qubit operator and uses its `SU(2)` representation for finding its decomposition in the specified basis up to a global phase via recursive implementation contained in `_solovay_kitaev` function.
2. `_approximate_set` method allows us to build an approximate set of Clifford sequences that are used in the base case return for the above `_solovay_kitaev` recursive call.  
3. `_group_commutator_decompose` allows us to perform the group commutation decomposition and uses the recipe mentioned in Sec 4.1 [arXiv:0505030](https://arxiv.org/abs/quant-ph/0505030).

**Benefits:** Should allow one to do the Clifford+T decomposition.

**Possible Drawbacks:** Better algorithms exist in the literature, specifically for decomposing rotation gates. 

**Related GitHub Issues:**



Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 0706193a1..e444885f8 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -6,6 +6,10 @@
 
 <h4>Decompose circuits into the Clifford+T gateset 🧩</h4>
 
+* `qml.transforms.decompositions.sk_decomposition` method implements the Solovay-Kitaev algorithm for
+  approximately decomposing any single-qubit operation to Clifford+T basis.
+  [(#4801)](https://github.com/PennyLaneAI/pennylane/pull/4801)
+
 <h4>Transforms (TODO: better title) 🤖</h4>
 
 <h4>Use an iterative approach for quantum phase estimation 🔄</h4>
diff --git a/pennylane/ops/functions/map_wires.py b/pennylane/ops/functions/map_wires.py
index f13eb40a1..36133c505 100644
--- a/pennylane/ops/functions/map_wires.py
+++ b/pennylane/ops/functions/map_wires.py
@@ -100,17 +100,17 @@ def map_wires(
             return new_op
         return input.map_wires(wire_map=wire_map)
     if isinstance(input, (QuantumScript, QNode)) or callable(input):
-        return _map_wires_transform(input, wire_map=wire_map)
+        return _map_wires_transform(input, wire_map=wire_map, queue=queue)
 
     raise ValueError(f"Cannot map wires of object {input} of type {type(input)}.")
 
 
 @partial(transform)
 def _map_wires_transform(
-    tape: qml.tape.QuantumTape, wire_map=None
+    tape: qml.tape.QuantumTape, wire_map=None, queue=False
 ) -> (Sequence[qml.tape.QuantumTape], Callable):
-    ops = [map_wires(op, wire_map) for op in tape.operations]
-    measurements = [map_wires(m, wire_map) for m in tape.measurements]
+    ops = [map_wires(op, wire_map, queue=queue) for op in tape.operations]
+    measurements = [map_wires(m, wire_map, queue=queue) for m in tape.measurements]
 
     out = tape.__class__(ops=ops, measurements=measurements, shots=tape.shots)
     out.trainable_params = tape.trainable_params
diff --git a/pennylane/transforms/__init__.py b/pennylane/transforms/__init__.py
index 3f2c4626a..d6512ed4e 100644
--- a/pennylane/transforms/__init__.py
+++ b/pennylane/transforms/__init__.py
@@ -74,6 +74,16 @@ There are also utility functions that take a circuit and return a DAG.
     ~transforms.CommutationDAG
     ~transforms.CommutationDAGNode
 
+Transform for Clifford+T compilation
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The following functions assist in decomposing operations to the Clifford+T basis.
+
+.. autosummary::
+    :toctree: api
+
+    ~transforms.sk_decomposition
+
 Transform for circuit cutting
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -230,6 +240,7 @@ from .compile import compile
 from .decompositions import (
     one_qubit_decomposition,
     two_qubit_decomposition,
+    sk_decomposition,
 )
 from .defer_measurements import defer_measurements
 from .sign_expand import sign_expand
diff --git a/pennylane/transforms/decompositions/__init__.py b/pennylane/transforms/decompositions/__init__.py
index 1eef12881..ab18a931d 100644
--- a/pennylane/transforms/decompositions/__init__.py
+++ b/pennylane/transforms/decompositions/__init__.py
@@ -17,3 +17,5 @@ unitary operations into sequences of elementary operations.
 
 from .single_qubit_unitary import one_qubit_decomposition
 from .two_qubit_unitary import two_qubit_decomposition
+
+from .clifford_t import sk_decomposition
diff --git a/pennylane/transforms/decompositions/clifford_t/__init__.py b/pennylane/transforms/decompositions/clifford_t/__init__.py
new file mode 100644
index 000000000..bc9bb0e07
--- /dev/null
+++ b/pennylane/transforms/decompositions/clifford_t/__init__.py
@@ -0,0 +1,16 @@
+# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""This module contains methods related to Clifford+T decomposition"""
+
+from .solovay_kitaev import sk_decomposition
diff --git a/pennylane/transforms/decompositions/clifford_t/solovay_kitaev.py b/pennylane/transforms/decompositions/clifford_t/solovay_kitaev.py
new file mode 100644
index 000000000..87d3fa655
--- /dev/null
+++ b/pennylane/transforms/decompositions/clifford_t/solovay_kitaev.py
@@ -0,0 +1,347 @@
+# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Solovay-Kitaev implementation for approximate single-qubit unitary decomposition."""
+import math
+import warnings
+from functools import lru_cache
+
+from scipy.spatial import KDTree
+
+import pennylane as qml
+from pennylane.queuing import QueuingManager
+from pennylane.tape import QuantumScript
+from pennylane.transforms.optimization import cancel_inverses
+
+# Defining Clifford+T basis
+_CLIFFORD_T_BASIS = {
+    "I": qml.Identity(0),
+    "X": qml.PauliX(0),
+    "Y": qml.PauliY(0),
+    "Z": qml.PauliZ(0),
+    "H": qml.Hadamard(0),
+    "T": qml.T(0),
+    "T*": qml.adjoint(qml.T(0)),
+    "S": qml.S(0),
+    "S*": qml.adjoint(qml.S(0)),
+}
+
+
+def _SU2_transform(matrix):
+    r"""Perform a U(2) to SU(2) transformation via a global phase addition.
+
+    A general element of :math:`\text{SU}_2(\mathbb{C})` has the following form:
+
+    .. math::
+        \text{SU}_{2} = \begin{bmatrix} \alpha & \beta \\ -\beta^{*} & \alpha^{*} \end{bmatrix},
+
+    where :math:`\alpha, \beta \in \mathbb{C}` and :math:`|\alpha|^2 + |\beta|^2 = 1`.
+    """
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", RuntimeWarning)
+        factor = qml.math.linalg.det(matrix)
+
+    # We put the phase first in [0, \pi] and then convert it to [0, \pi)
+    gphase = qml.math.mod(qml.math.angle(factor), 2 * math.pi) / 2
+    rphase = (-1) ** qml.math.isclose(gphase, math.pi)
+
+    # Get the final matrix form using the phase information
+    s2_mat = rphase * matrix * qml.math.exp(-1j * qml.math.cast_like(gphase, 1j))
+    return s2_mat, gphase if rphase == 1 else 0.0
+
+
+def _quaternion_transform(matrix):
+    r"""Perform a SU(2) to quaternion transformation.
+
+    Any element :math:`\text{SU}_2(\mathbb{C})` can be be written as a unique quaternion
+    :math:`\alpha_0\mathbb{I} + \alpha_1\mathbf{i} + \alpha_2\mathbf{j} + \alpha_3\mathbf{k}`,
+    where :math:`\mathbf{i}=-i\mathbf{X},\ \mathbf{j}=-i\mathbf{Y},\ \text{and}\ \mathbf{k}=-i\mathbf{Z}`.
+    """
+    return qml.math.array(
+        [
+            qml.math.real(matrix[0, 0]),
+            -qml.math.imag(matrix[0, 1]),
+            -qml.math.real(matrix[0, 1]),
+            -qml.math.imag(matrix[0, 0]),
+        ],
+        dtype=float,
+    )
+
+
+def _contains_SU2(op_mat, ops_vecs, tol=1e-8):
+    r"""Checks if a given SU(2) matrix is contained in a list of quaternions for a given tolerance.
+
+    Args:
+        op_mat (TensorLike): SU(2) matrix for the operation to be searched
+        op_vecs (list(TensorLike)): List of quaternion for the operations that makes the search space.
+        tol (float): Tolerance for the match to be considered ``True``.
+
+    Returns:
+        Tuple(bool, TensorLike): A bool that shows whether an operation similar to the given operations
+        was found, and the quaternion representation of the searched operation.
+    """
+    node_points = qml.math.array(ops_vecs)
+    gate_points = qml.math.array([_quaternion_transform(op_mat)])
+
+    tree = KDTree(node_points)
+    dist = tree.query(gate_points, workers=-1)[0][0]
+
+    return (dist < tol, gate_points[0])
+
+
+@lru_cache()
+def _approximate_set(basis_gates, max_length=10):
+    r"""Builds an approximate unitary set required for the `Solovay-Kitaev algorithm <https://arxiv.org/abs/quant-ph/0505030>`_.
+
+    Args:
+        basis_gates (list(str)): Basis set to be used for Solovay-Kitaev decomposition build using
+            following terms, ``['X', 'Y', 'Z', 'H', 'T', 'T*', 'S', 'S*']``, where `*` refers
+            to the gate adjoint.
+        max_length (int): Maximum expansion length of Clifford+T sequences in the approximation set. Default is `10`
+
+    Returns:
+        Tuple(list[list[~pennylane.operation.Operation]], list[TensorLike], list[float], list[TensorLike]): A tuple containing the list of
+        Clifford+T sequences that will be used for approximating a matrix in the base case of recursive implementation of
+        Solovay-Kitaev algorithm, with their corresponding SU(2) representations, global phases, and quaternion representations.
+    """
+    # Maintains the basis gates
+    basis = [_CLIFFORD_T_BASIS[gate.upper()] for gate in basis_gates]
+
+    # Get the SU(2) data for the gates in basis set
+    basis_mat, basis_gph = {}, {}
+    for gate in basis:
+        su2_mat, su2_gph = _SU2_transform(gate.matrix())
+        basis_mat.update({gate: su2_mat})
+        basis_gph.update({gate: su2_gph})
+
+    # Maintains a trie-like structure for each depth
+    gtrie_ids = [[[gate] for gate in basis]]
+    gtrie_mat = [list(basis_mat.values())]
+    gtrie_gph = [list(basis_gph.values())]
+
+    # Maintains the approximate set for gates' names, SU(2)s, global phases and quaternions
+    approx_set_ids = list(gtrie_ids[0])
+    approx_set_mat = list(gtrie_mat[0])
+    approx_set_gph = list(gtrie_gph[0])
+    approx_set_qat = [_quaternion_transform(mat) for mat in approx_set_mat]
+
+    # We will perform a breadth-first search (BFS) style set building for the set
+    for depth in range(max_length - 1):
+        # Add the containers for next depth while we explore the current
+        gtrie_id, gtrie_mt, gtrie_gp = [], [], []
+        for node, su2m, gphase in zip(gtrie_ids[depth], gtrie_mat[depth], gtrie_gph[depth]):
+            # Get the last operation in the current node
+            last_op = qml.adjoint(node[-1], lazy=False) if node else None
+
+            # Now attempt extending the current node for each basis gate
+            for op in basis:
+                # If basis gate is adjoint of last op in the node, skip.
+                if qml.equal(op, last_op):
+                    continue
+
+                # Extend and check if the node already exists in the approximate set.
+                su2_gp = basis_gph[op] + gphase
+                su2_op = (-1.0) ** bool(su2_gp >= math.pi) * (basis_mat[op] @ su2m)
+                exists, quaternion = _contains_SU2(su2_op, approx_set_qat)
+                if not exists:
+                    approx_set_ids.append(node + [op])
+                    approx_set_mat.append(su2_op)
+                    approx_set_qat.append(quaternion)
+
+                    # Add to the containers for next depth
+                    gtrie_id.append(node + [op])
+                    gtrie_mt.append(su2_op)
+
+                    # Add the global phase data
+                    global_phase = qml.math.mod(su2_gp, math.pi)
+                    approx_set_gph.append(global_phase)
+                    gtrie_gp.append(global_phase)
+
+        # Add to the next depth for next iteration
+        gtrie_ids.append(gtrie_id)
+        gtrie_mat.append(gtrie_mt)
+        gtrie_gph.append(gtrie_gp)
+
+    return approx_set_ids, approx_set_mat, approx_set_gph, approx_set_qat
+
+
+def _group_commutator_decompose(matrix, tol=1e-5):
+    r"""Performs a group commutator decomposition :math:`U = V' \times W' \times V'^{\dagger} \times W'^{\dagger}`
+    as given in the Section 4.1 of `arXiv:0505030 <https://arxiv.org/abs/quant-ph/0505030>`_."""
+    # Use the quaternion form to get the rotation axis and angle on the Bloch sphere,
+    # while using clipping for dealing with floating point precision errors.
+    quaternion = _quaternion_transform(matrix)
+    theta, axis = 2 * qml.math.arccos(qml.math.clip(quaternion[0], -1.0, 1.0)), quaternion[1:]
+
+    # Early return for the case where matrix is I or -I, where I is Identity.
+    if qml.math.allclose(axis, 0.0, atol=tol) and qml.math.isclose(theta % math.pi, 0.0, atol=tol):
+        return qml.math.eye(2, dtype=complex), qml.math.eye(2, dtype=complex)
+
+    # The angle phi comes from the Eq. 10 in the Solovay-Kitaev algorithm paper (arXiv:0505030).
+    phi = 2.0 * qml.math.arcsin(qml.math.sqrt(qml.math.sqrt((0.5 - 0.5 * qml.math.cos(theta / 2)))))
+
+    # Begin decomposition by computing the rotation operations V and W.
+    v = qml.RX(phi, [0])
+    w = qml.RY(2 * math.pi - phi, [0]) if axis[2] > 0 else qml.RY(phi, [0])
+
+    # Get the similarity transormation matrices S and S.adjoint().
+    ud = qml.math.linalg.eig(matrix)[1]
+    vwd = qml.math.linalg.eig(qml.matrix(v @ w @ v.adjoint() @ w.adjoint()))[1]
+    s = ud @ qml.math.conj(qml.math.transpose(vwd))
+    sdg = vwd @ qml.math.conj(qml.math.transpose(ud))
+
+    # Get the required matrices V' and W'.
+    v_hat = s @ v.matrix() @ sdg
+    w_hat = s @ w.matrix() @ sdg
+
+    return w_hat, v_hat
+
+
+def sk_decomposition(op, epsilon, *, max_depth=5, basis_set=("T", "T*", "H"), basis_length=10):
+    r"""Approximate an arbitrary single-qubit gate in the Clifford+T basis using the `Solovay-Kitaev algorithm <https://arxiv.org/abs/quant-ph/0505030>`_.
+
+    This method implements the Solovay-Kitaev decomposition algorithm that approximates any single-qubit
+    operation with :math:`\epsilon > 0` error. The procedure exits when the approximation error
+    becomes less than :math:`\epsilon`, or when ``max_depth`` approximation passes have been made. In the
+    latter case, the approximation error could be :math:`\geq \epsilon`.
+
+    This algorithm produces a decomposition with :math:`O(\text{log}^{3.97}(1/\epsilon))` operations.
+
+    Args:
+        op (~pennylane.operation.Operation): A single-qubit gate operation.
+        epsilon (float): The maximum permissible error.
+
+    Keyword Args:
+        max_depth (int): The maximum number of approximation passes. A smaller :math:`\epsilon` would generally require
+            a greater number of passes. Default is ``5``.
+        basis_set (list[str]): Basis set to be used for the decomposition and building an approximate set internally.
+            It accepts the following gate terms: ``['X', 'Y', 'Z', 'H', 'T', 'T*', 'S', 'S*']``, where ``*`` refers
+            to the gate adjoint. Default value is ``['T', 'T*', 'H']``.
+        basis_length (int): Maximum expansion length of Clifford+T sequences in the internally-built approximate set.
+            Default is ``10``.
+
+    Returns:
+        list[~pennylane.operation.Operation]: A list of gates in the Clifford+T basis set that approximates the given
+        operation along with a final global phase operation. The operations are in the circuit-order.
+
+    Raises:
+        ValueError: If the given operator acts on more than one wires.
+
+    **Example**
+
+    Suppose one would like to decompose :class:`~.RZ` with rotation angle :math:`\phi = \pi/3`:
+
+    .. code-block:: python3
+
+        import numpy as np
+        import pennylane as qml
+
+        op  = qml.RZ(np.pi/3, wires=0)
+
+        # Get the gate decomposition in ['T', 'T*', 'H']
+        ops = qml.transforms.decompositions.sk_decomposition(op, epsilon=1e-3)
+
+        # Get the approximate matrix from the ops
+        matrix_sk = qml.prod(*reversed(ops)).matrix()
+
+    When the function is run for a sufficient ``depth`` with a good enough approximate set,
+    the output gate sequence should implement the same operation approximately.
+
+    >>> qml.math.allclose(op.matrix(), matrix_sk, atol=1e-3)
+    True
+
+    """
+    # Check for length of wires in the operation
+    if len(op.wires) != 1:
+        raise ValueError(
+            f"Operator must be a single qubit operation, got {op} acting on {op.wires} wires."
+        )
+
+    with QueuingManager.stop_recording():
+        # Build the approximate set with caching
+        approx_set_ids, approx_set_mat, approx_set_gph, approx_set_qat = _approximate_set(
+            basis_set, max_length=basis_length
+        )
+
+        # Build the k-d tree with the current approximation set for querying in the base case
+        kd_tree = KDTree(qml.math.array(approx_set_qat))
+
+        # Obtain the SU(2) and quaternion for the operation
+        gate_mat, gate_gph = _SU2_transform(op.matrix())
+        gate_qat = _quaternion_transform(gate_mat)
+
+        def _solovay_kitaev(umat, n, u_n1_ids, u_n1_mat):
+            """Recursive method as given in the Section 3 of arXiv:0505030"""
+
+            if not n:
+                # Check the approximate gate in our approximate set
+                seq_node = qml.math.array([_quaternion_transform(umat)])
+                _, [index] = kd_tree.query(seq_node, workers=-1)
+                return approx_set_ids[index], approx_set_mat[index]
+
+            # Get the decomposition for the remaining unitary: U @ U'.adjoint()
+            v_n, w_n = _group_commutator_decompose(
+                umat @ qml.math.conj(qml.math.transpose(u_n1_mat))
+            )
+
+            # Get the approximation for the residual commutator unitaries: V and W
+            c_ids_mats = []
+            for c_n in [v_n, w_n]:
+                # Get the approximation for each commutator iteratively: C --> C'
+                c_n1_ids, c_n1_mat = None, None
+                for i in range(n):
+                    c_n1_ids, c_n1_mat = _solovay_kitaev(c_n, i, c_n1_ids, c_n1_mat)
+
+                # Get the adjoints C' --> C'.adjoint()
+                c_n1_ids_adj = [qml.adjoint(gate, lazy=False) for gate in reversed(c_n1_ids)]
+                c_n1_mat_adj = qml.math.conj(qml.math.transpose(c_n1_mat))
+
+                # Store the decompositions and matrices for C'
+                c_ids_mats.append([c_n1_ids, c_n1_mat, c_n1_ids_adj, c_n1_mat_adj])
+
+            # Get the V' and W'
+            v_n1_ids, v_n1_mat, v_n1_ids_adj, v_n1_mat_adj = c_ids_mats[0]
+            w_n1_ids, w_n1_mat, w_n1_ids_adj, w_n1_mat_adj = c_ids_mats[1]
+
+            # Build the operations and their SU(2) for return
+            approx_ids = u_n1_ids + w_n1_ids_adj + v_n1_ids_adj + w_n1_ids + v_n1_ids
+            approx_mat = v_n1_mat @ w_n1_mat @ v_n1_mat_adj @ w_n1_mat_adj @ u_n1_mat
+
+            return approx_ids, approx_mat
+
+        # If we have it already, use that, otherwise proceed for decomposition
+        _, [index] = kd_tree.query(qml.math.array([gate_qat]), workers=-1)
+        decomposition, u_prime = approx_set_ids[index], approx_set_mat[index]
+
+        # Iterate until max_depth while doing an epsilon-error comparision
+        for depth in range(max_depth):
+            # For a SU(2) matrix Hilbert-Schmidt norm is √(|α|^2 + |β|^2),
+            # which is simply the L2-norm for the first row of the matrix.
+            if qml.math.norm(gate_mat[0] - u_prime[0]) <= epsilon:
+                break
+            # Approximate the residual with the approximation from previous iteration
+            decomposition, u_prime = _solovay_kitaev(gate_mat, depth + 1, decomposition, u_prime)
+
+        # Remove inverses if any in the decomposition and handle trivial case
+        [new_tape], _ = cancel_inverses(QuantumScript(decomposition or [qml.Identity(0)]))
+
+    # Map the wires to that of the operation and queue
+    [map_tape], _ = qml.map_wires(new_tape, wire_map={0: op.wires[0]}, queue=True)
+
+    # Get phase information based on the decomposition effort
+    global_phase = qml.GlobalPhase(
+        approx_set_gph[index] - gate_gph if depth or qml.math.isclose(gate_gph, 0.0) else 0.0
+    )
+
+    # Return the gates from the mapped tape and global phase
+    return map_tape.operations + [global_phase]


Test patch:
diff --git a/tests/transforms/test_clifford_t/test_solovay_kitaev.py b/tests/transforms/test_clifford_t/test_solovay_kitaev.py
new file mode 100644
index 000000000..1a5c98572
--- /dev/null
+++ b/tests/transforms/test_clifford_t/test_solovay_kitaev.py
@@ -0,0 +1,222 @@
+# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for Solovay-Kitaev decomposition's implementation."""
+
+import math
+
+import pytest
+import scipy as sp
+import pennylane as qml
+
+from pennylane.transforms.decompositions.clifford_t.solovay_kitaev import (
+    _SU2_transform,
+    _quaternion_transform,
+    _contains_SU2,
+    _approximate_set,
+    _group_commutator_decompose,
+    sk_decomposition,
+)
+
+
+@pytest.mark.parametrize(
+    ("op", "matrix", "phase"),
+    [
+        (
+            qml.Identity(0),  # computed manually
+            qml.math.array([[1.0 + 0.0j, 0.0 + 0.0j], [0.0 + 0.0j, 1.0 + 0.0j]]),
+            0.0,
+        ),
+        (
+            qml.T(0),  # computed manually
+            qml.math.array(
+                [[0.92387953 - 0.38268343j, 0.0 + 0.0j], [0.0 + 0.0j, 0.92387953 + 0.38268343j]]
+            ),
+            0.3926990817,
+        ),
+        (
+            qml.Hadamard(0),  # computed manually
+            qml.math.array(
+                [
+                    [-0.70710678j, -0.70710678j],
+                    [-0.70710678j, 0.70710678j],
+                ]
+            ),
+            1.5707963268,
+        ),
+    ],
+)
+def test_su2_transform(op, matrix, phase):
+    """Test the functionality to create SU(2) representation"""
+    su2_matrix, su2_phase = _SU2_transform(op.matrix())
+    assert qml.math.allclose(su2_matrix, matrix)
+    assert qml.math.isclose(su2_phase, phase)
+
+
+@pytest.mark.parametrize(
+    "op,quaternion",
+    [
+        (qml.Identity(0), qml.math.array([1.0, 0.0, 0.0, 0.0])),
+        (qml.T(0), qml.math.array([0.92387953, 0.0, 0.0, 0.38268343])),
+        (qml.Hadamard(0), qml.math.array([0, 0.707106781, 0, 0.707106781])),
+    ],
+)
+def test_quaternion_transform(op, quaternion):
+    """Test the functionality to create quaternion representation"""
+    su2_matrix, _ = _SU2_transform(op.matrix())
+    op_quaternion = _quaternion_transform(su2_matrix)
+    assert qml.math.allclose(op_quaternion, quaternion)
+
+
+def test_contains_SU2():
+    """Test the functionality to check an operation is in an existing approximate set"""
+    op = qml.prod(qml.T(0), qml.T(0), qml.T(0), qml.T(0), qml.T(0))
+    target = _SU2_transform(op.matrix())[0]
+
+    approx_ids, _, _, approx_vec = _approximate_set(("T", "T*", "H"), max_length=3)
+
+    exists, quaternion = _contains_SU2(target, approx_vec)
+    assert exists
+
+    result = [qml.adjoint(qml.T(0)), qml.adjoint(qml.T(0)), qml.adjoint(qml.T(0))]
+
+    node_points = qml.math.array(approx_vec)
+    seq_node = qml.math.array([quaternion])
+    _, [index] = sp.spatial.KDTree(node_points).query(seq_node, workers=-1)
+
+    assert approx_ids[index] == result
+
+
+def test_approximate_sets():
+    """Test the functionality to create approximate set"""
+    # Check length of approximated sets with gate depths
+    approx_set1, _, _, _ = _approximate_set((), max_length=0)
+    approx_set2, _, _, _ = _approximate_set(("t", "h"), max_length=1)
+    assert len(approx_set1) == 0
+    assert len(approx_set2) == 2
+
+    # Verify exist-check reduce redundancy via GP-like summation
+    approx_set3, _, _, _ = _approximate_set(("t", "t*", "h"), max_length=3)
+    approx_set4, _, _, _ = _approximate_set(("t", "t*", "h"), max_length=5)
+    assert len(approx_set3) == 21  # should be <  39 = 3 + 3x3 + 3x3x3
+    assert len(approx_set4) == 84  # should be < 263 = 3 + ... + 3x3x3x3x3
+
+
+@pytest.mark.parametrize(
+    ("op"),
+    [
+        qml.U3(0.0, 0.0, 0.0, wires=[0]),
+        qml.U3(1.0, 2.0, 3.0, wires=[1]),
+        qml.U3(-1.0, 2.0, 3.0, wires=["b"]),
+        qml.U3(1.0, 2.0, -3.0, wires=[0]),
+    ],
+)
+def test_group_commutator_decompose(op):
+    """Test group commutator deocompose method for SU(2)"""
+
+    su2_matrix, _ = _SU2_transform(qml.matrix(op))
+
+    v_hat, w_hat = _group_commutator_decompose(su2_matrix)
+    decomposed_matrix = (
+        v_hat
+        @ w_hat
+        @ qml.math.conj(qml.math.transpose(v_hat))
+        @ qml.math.conj(qml.math.transpose(w_hat))
+    )
+
+    assert qml.math.allclose(decomposed_matrix, su2_matrix)
+
+
+@pytest.mark.parametrize(
+    ("op"),
+    [
+        qml.RX(math.pi / 42, wires=[1]),
+        qml.RY(math.pi / 7, wires=["a"]),
+        qml.prod(*[qml.RX(1.0, "a"), qml.T("a")]),
+        qml.prod(*[qml.T(0), qml.Hadamard(0)] * 5),
+    ],
+)
+def test_solovay_kitaev(op):
+    """Test Solovay-Kitaev decomposition method"""
+
+    with qml.queuing.AnnotatedQueue() as q:
+        gates = sk_decomposition(op, epsilon=1e-4, max_depth=5, basis_set=("T", "T*", "H"))
+    assert q.queue == gates
+
+    matrix_sk = qml.prod(*reversed(gates)).matrix()
+
+    assert qml.math.allclose(qml.matrix(op), matrix_sk, atol=1e-2)
+    assert qml.prod(*gates, lazy=False).wires == op.wires
+
+
+@pytest.mark.parametrize(
+    ("basis_length", "basis_set"),
+    [(10, ("T*", "T", "H")), (8, ("H", "T")), (10, ("H", "S", "T"))],
+)
+def test_solovay_kitaev_with_basis_gates(basis_length, basis_set):
+    """Test Solovay-Kitaev decomposition method with additional basis information provided."""
+    op = qml.PhaseShift(math.pi / 4, 0)
+    gates = sk_decomposition(op, 1e-4, max_depth=3, basis_set=basis_set, basis_length=basis_length)
+
+    matrix_sk = qml.prod(*reversed(gates)).matrix()
+
+    assert qml.math.allclose(qml.matrix(op), matrix_sk, atol=1e-5)
+
+
+@pytest.mark.parametrize(
+    "op,query_count",
+    [
+        (qml.Hadamard(0), 1),
+        (qml.prod(qml.Hadamard(0), qml.T(0)), 1),
+        (qml.prod(qml.Hadamard(0), qml.T(0), qml.RX(1e-9, 0)), 1),
+        (qml.RZ(1.23, 0), 27),
+        (qml.prod(*[qml.T(0), qml.Hadamard(0)] * 5), 1),
+    ],
+)
+def test_close_approximations_do_not_go_deep(op, query_count, mocker):
+    """Test that the recursive solver is only used when necessary."""
+    basis_set = ("H", "T", "T*")
+    basis_length = 10
+    _ = _approximate_set(basis_set, max_length=basis_length)  # pre-compute so spy is accurate
+
+    spy = mocker.spy(sp.spatial.KDTree, "query")
+    _ = sk_decomposition(op, 1e-4, max_depth=3, basis_set=basis_set, basis_length=basis_length)
+    assert spy.call_count == query_count
+
+
+@pytest.mark.parametrize("epsilon", [2e-2, 3e-2, 7e-2])
+def test_epsilon_value_respected(epsilon):
+    """Test that resulting decompositions are within an epsilon value."""
+    op = qml.RX(1.234, 0)
+    gates = sk_decomposition(op, epsilon, max_depth=5)
+    matrix_sk = qml.prod(*reversed(gates)).matrix()
+    assert qml.math.norm(qml.matrix(op)[0] - matrix_sk[0]) < epsilon
+
+
+def test_epsilon_value_effect():
+    """Test that different epsilon values create different decompositions."""
+    op = qml.RZ(math.pi / 5, 0)
+    decomp_with_error = sk_decomposition(op, 9e-2, max_depth=5)
+    decomp_less_error = sk_decomposition(op, 1e-2, max_depth=5)
+    assert len(decomp_with_error) < len(decomp_less_error)
+
+
+def test_exception():
+    """Test operation wire exception in Solovay-Kitaev"""
+    op = qml.SingleExcitation(1.0, wires=[1, 2])
+
+    with pytest.raises(
+        ValueError,
+        match=r"Operator must be a single qubit operation",
+    ):
+        sk_decomposition(op, epsilon=1e-4, max_depth=1)