Repository: PennyLaneAI/pennylane
Pull request #4802 (2023-11-07 22:47:13)

Clifford+T Decomposition Transform
Context: Transform to perform Clifford+T decomposition

Description of the Change: Unrolls the circuit into Clifford+Rotation basis and then uses the Solovay-Kitaev algorithm for approximating arbitrary z-rotations.

Benefits: Allows one to transform any random circuit to a Clifford+Phase basis.

Possible Drawbacks: Might be sub-optimal

Related GitHub Issues:

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 268e85a86..0f8392706 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.clifford_t_decomposition` method implements a transform for
+  approximately decomposing any circuit to the Clifford+T basis.
+  [(#4802)](https://github.com/PennyLaneAI/pennylane/pull/4802)
+
 * `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)
diff --git a/pennylane/__init__.py b/pennylane/__init__.py
index f9a650475..8a62587a7 100644
--- a/pennylane/__init__.py
+++ b/pennylane/__init__.py
@@ -98,6 +98,7 @@ from pennylane.transforms import (
     commutation_dag,
     pattern_matching,
     pattern_matching_optimization,
+    clifford_t_decomposition,
 )
 from pennylane.ops.functions import (
     dot,
diff --git a/pennylane/math/single_dispatch.py b/pennylane/math/single_dispatch.py
index 536581d74..6e1a24c4c 100644
--- a/pennylane/math/single_dispatch.py
+++ b/pennylane/math/single_dispatch.py
@@ -235,6 +235,7 @@ ar.autoray._SUBMODULE_ALIASES["tensorflow", "arcsin"] = "tensorflow.math"
 ar.autoray._SUBMODULE_ALIASES["tensorflow", "arccos"] = "tensorflow.math"
 ar.autoray._SUBMODULE_ALIASES["tensorflow", "arctan"] = "tensorflow.math"
 ar.autoray._SUBMODULE_ALIASES["tensorflow", "arctan2"] = "tensorflow.math"
+ar.autoray._SUBMODULE_ALIASES["tensorflow", "mod"] = "tensorflow.math"
 ar.autoray._SUBMODULE_ALIASES["tensorflow", "diag"] = "tensorflow.linalg"
 ar.autoray._SUBMODULE_ALIASES["tensorflow", "kron"] = "tensorflow.experimental.numpy"
 ar.autoray._SUBMODULE_ALIASES["tensorflow", "moveaxis"] = "tensorflow.experimental.numpy"
@@ -527,6 +528,7 @@ ar.register_function("torch", "expand_dims", lambda x, axis: _i("torch").unsquee
 ar.register_function("torch", "shape", lambda x: tuple(x.shape))
 ar.register_function("torch", "gather", lambda x, indices: x[indices])
 ar.register_function("torch", "equal", lambda x, y: _i("torch").eq(x, y))
+ar.register_function("torch", "mod", lambda x, y: x % y)
 
 ar.register_function(
     "torch",
diff --git a/pennylane/transforms/__init__.py b/pennylane/transforms/__init__.py
index b3dce3869..993b6f007 100644
--- a/pennylane/transforms/__init__.py
+++ b/pennylane/transforms/__init__.py
@@ -74,10 +74,17 @@ There are also utility functions that take a circuit and return a DAG.
     ~transforms.CommutationDAG
     ~transforms.CommutationDAGNode
 
-Transform for Clifford+T compilation
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Transform for Clifford+T decomposition
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-The following functions assist in decomposing operations to the Clifford+T basis.
+This transform accepts quantum circuits and decomposes them to the Clifford+T basis.
+
+.. autosummary::
+    :toctree: api
+
+    ~clifford_t_decomposition
+
+The following function assists in decomposing operations to the Clifford+T basis.
 
 .. autosummary::
     :toctree: api
@@ -199,6 +206,7 @@ from .decompositions import (
     one_qubit_decomposition,
     two_qubit_decomposition,
     sk_decomposition,
+    clifford_t_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 ab18a931d..706898f5c 100644
--- a/pennylane/transforms/decompositions/__init__.py
+++ b/pennylane/transforms/decompositions/__init__.py
@@ -18,4 +18,4 @@ 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
+from .clifford_t import clifford_t_decomposition, sk_decomposition
diff --git a/pennylane/transforms/decompositions/clifford_t/__init__.py b/pennylane/transforms/decompositions/clifford_t/__init__.py
index bc9bb0e07..fffd399b2 100644
--- a/pennylane/transforms/decompositions/clifford_t/__init__.py
+++ b/pennylane/transforms/decompositions/clifford_t/__init__.py
@@ -13,4 +13,6 @@
 # limitations under the License.
 """This module contains methods related to Clifford+T decomposition"""
 
+from .clifford_t_transform import clifford_t_decomposition
+
 from .solovay_kitaev import sk_decomposition
diff --git a/pennylane/transforms/decompositions/clifford_t/clifford_t_transform.py b/pennylane/transforms/decompositions/clifford_t/clifford_t_transform.py
new file mode 100644
index 000000000..7989353e9
--- /dev/null
+++ b/pennylane/transforms/decompositions/clifford_t/clifford_t_transform.py
@@ -0,0 +1,500 @@
+# 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.
+"""Transform function for the Clifford+T decomposition."""
+
+import math
+import warnings
+from itertools import product
+from typing import Sequence, Callable
+
+import pennylane as qml
+from pennylane.queuing import QueuingManager
+from pennylane.transforms.core import transform
+from pennylane.tape import QuantumTape
+from pennylane.transforms.optimization import (
+    cancel_inverses,
+    commute_controlled,
+    merge_rotations,
+    remove_barrier,
+)
+from pennylane.transforms.optimization.optimization_utils import find_next_gate, _fuse_global_phases
+from .solovay_kitaev import sk_decomposition
+
+# Single qubits Clifford+T gates in PL
+_CLIFFORD_T_ONE_GATES = [
+    qml.Identity,
+    qml.PauliX,
+    qml.PauliY,
+    qml.PauliZ,
+    qml.Hadamard,
+    qml.S,
+    qml.SX,
+    qml.T,
+]
+
+# Two qubits Clifford+T gates in PL
+_CLIFFORD_T_TWO_GATES = [
+    qml.CNOT,
+    qml.CY,
+    qml.CZ,
+    qml.SWAP,
+    qml.ISWAP,
+]
+
+# Single-parameter gates in PL
+# note: _simplify_param makes use of their periodic nature,
+# any additions to this, should be reflected there as well.
+_PARAMETER_GATES = (qml.RX, qml.RY, qml.RZ, qml.Rot, qml.PhaseShift)
+
+# Clifford+T gate set
+_CLIFFORD_T_GATES = tuple(_CLIFFORD_T_ONE_GATES + _CLIFFORD_T_TWO_GATES) + (qml.GlobalPhase,)
+
+# Gates to be skipped during decomposition
+_SKIP_OP_TYPES = (qml.Barrier, qml.Snapshot, qml.WireCut)
+
+
+def _check_clifford_op(op, use_decomposition=False):
+    r"""Checks if an operator is Clifford or not.
+
+    For a given unitary operator :math:`U` acting on :math:`N` qubits, this method checks that the
+    transformation :math:`UPU^{\dagger}` maps the Pauli tensor products :math:`P = {I, X, Y, Z}^{\otimes N}`
+    to Pauli tensor products using the decomposition of the matrix for :math:`U` in the Pauli basis.
+
+    Args:
+        op (~pennylane.operation.Operation): the operator that needs to be tested
+        use_decomposition (bool): if ``True``, use operator's decomposition to compute the matrix, in case
+            it doesn't define a ``compute_matrix`` method. Default is ``False``.
+
+    Returns:
+        Bool that represents whether the provided operator is Clifford or not.
+    """
+
+    # Check if matrix can be calculated for the operator
+    if (not op.has_matrix and not use_decomposition) or (
+        use_decomposition and not op.expand().wires
+    ):
+        return False
+
+    # Compute the LCUs for the operator in Pauli basis
+    pauli_terms = qml.pauli_decompose(qml.matrix(op), wire_order=op.wires, check_hermitian=False)
+    pauli_terms_adj = qml.Hamiltonian(qml.math.conj(pauli_terms.coeffs), pauli_terms.ops)
+
+    # Build Pauli tensor products P in the Hamiltonian form
+    def pauli_group(x):
+        return [qml.Identity(x), qml.PauliX(x), qml.PauliY(x), qml.PauliZ(x)]
+
+    # Build PauliSentence and Hamiltonians representations for set P
+    pauli_sens = [
+        qml.pauli.pauli_sentence(qml.prod(*pauli))
+        for pauli in product(*(pauli_group(idx) for idx in op.wires))
+    ]
+    pauli_hams = (pauli_sen.hamiltonian(wire_order=op.wires) for pauli_sen in pauli_sens)
+
+    # Perform U@P@U^\dagger and check if the result exists in set P
+    for pauli_prod in product([pauli_terms], pauli_hams, [pauli_terms_adj]):
+        # hopefully op_math.prod scales better than matrix multiplication, i.e., O((2^N)^3)
+        upu = qml.pauli.pauli_sentence(qml.prod(*pauli_prod))
+        upu.simplify()
+        # Pauli sum always lie outside set P
+        if len(upu) != 1:
+            return False
+
+    return True
+
+
+def check_clifford_t(op, use_decomposition=False):
+    r"""Checks whether the gate is in the standard Clifford+T basis.
+
+    For a given unitary operator :math:`U` acting on :math:`N` qubits, which is not a T-gate,
+    this method checks that the transformation :math:`UPU^{\dagger}` maps the Pauli tensor products
+    :math:`P = {I, X, Y, Z}^{\otimes N}` to Pauli tensor products using the decomposition of the
+    matrix for :math:`U` in the Pauli basis.
+
+    Args:
+        op (~pennylane.operation.Operation): the operator that needs to be checked
+        use_decomposition (bool): if ``True``, use operator's decomposition to compute the matrix, in case
+            it doesn't define a ``compute_matrix`` method. Default is ``False``.
+
+    Returns:
+        Bool that represents whether the provided operator is Clifford+T or not.
+    """
+    # Get the base operation for an adjointed operation
+    base = getattr(op, "base", None)
+
+    # Save time and check from the pre-computed list
+    if isinstance(op, _CLIFFORD_T_GATES) or isinstance(base, _CLIFFORD_T_GATES):
+        return True
+
+    # Save time and check from the parameter of rotation gates
+    if isinstance(op, _PARAMETER_GATES) or isinstance(base, _PARAMETER_GATES):
+        theta = op.data[0]
+        return (
+            False
+            if qml.math.is_abstract(theta)
+            else qml.math.allclose(qml.math.mod(theta, math.pi), 0.0)
+        )
+
+    return _check_clifford_op(op, use_decomposition=use_decomposition)
+
+
+def _simplify_param(theta, gate):
+    r"""Check if the parameter allows simplification for the rotation gate.
+
+    For the cases where theta is an integer multiple of π: (a) returns a global phase
+    when even, and (b) returns combination of provided gate with global phase when odd.
+    In rest of the other cases it would return None.
+    """
+    if qml.math.is_abstract(theta):  # pragma: no cover
+        return None
+
+    if qml.math.allclose(theta, 0.0, atol=1e-6):
+        return [qml.GlobalPhase(0.0)]
+
+    rem_, mod_ = qml.math.divide(theta, math.pi), qml.math.mod(theta, math.pi)
+    if qml.math.allclose(mod_, 0.0, atol=1e-6):
+        ops = [qml.GlobalPhase(theta / 2)]
+        if qml.math.allclose(qml.math.mod(rem_, 2), 1.0, atol=1e-6):
+            ops.append(gate)
+        return ops
+
+    return None
+
+
+# pylint: disable= too-many-branches,
+def _rot_decompose(op):
+    r"""Decomposes a rotation operation: :class:`~.Rot`, :class:`~.RX`, :class:`~.RY`, :class:`~.RZ`,
+    :class:`~.PhaseShift` into a basis composed of :class:`~.RZ`, :class:`~.S`, and :class:`~.Hadamard`.
+    """
+    d_ops = []
+    # Extend for Rot operation with Rz.Ry.Rz decompositions
+    if isinstance(op, qml.Rot):
+        (phi, theta, omega), wires = op.parameters, op.wires
+        for dec in [qml.RZ(phi, wires), qml.RY(theta, wires), qml.RZ(omega, wires)]:
+            d_ops.extend(_rot_decompose(dec))
+        return d_ops
+
+    (theta,), wires = op.data, op.wires
+    if isinstance(op, qml.ops.Adjoint):  # pylint: disable=no-member
+        ops_ = _rot_decompose(op.base.adjoint())
+    elif isinstance(op, qml.RX):
+        ops_ = _simplify_param(theta, qml.PauliX(wires=wires))
+        if ops_ is None:  # Use Rx = H @ Rz @ H
+            ops_ = [qml.Hadamard(wires), qml.RZ(theta, wires), qml.Hadamard(wires)]
+    elif isinstance(op, qml.RY):
+        ops_ = _simplify_param(theta, qml.PauliY(wires=wires))
+        if ops_ is None:  # Use Ry = S @ H @ Rz @ H @ S.adjoint()
+            ops_ = [
+                qml.S(wires),
+                qml.Hadamard(wires),
+                qml.RZ(theta, wires),
+                qml.Hadamard(wires),
+                qml.adjoint(qml.S(wires)),
+            ][::-1]
+    elif isinstance(op, qml.RZ):
+        ops_ = _simplify_param(theta, qml.PauliZ(wires=wires))
+        if ops_ is None:
+            ops_ = [qml.RZ(theta, wires)]
+    elif isinstance(op, qml.PhaseShift):
+        ops_ = _simplify_param(theta, qml.PauliZ(wires=wires))
+        if ops_ is None:
+            ops_ = [qml.RZ(theta, wires=wires), qml.GlobalPhase(theta / 2)]
+        else:
+            ops_.append(qml.GlobalPhase(-theta / 2))
+
+    else:
+        raise ValueError(
+            f"Operation {op} is not a supported Pauli rotation: qml.RX, qml.RY, qml.RZ, qml.Rot and qml.PhaseShift."
+        )
+
+    d_ops.extend(ops_)
+    return d_ops
+
+
+def _one_qubit_decompose(op):
+    r"""Decomposition for single qubit operations using combination of :class:`~.RZ`, :class:`~.S`, and
+    :class:`~.Hadamard`."""
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", RuntimeWarning)
+        sd_ops = qml.transforms.one_qubit_decomposition(
+            qml.matrix(op), op.wires, "ZXZ", return_global_phase=True
+        )
+    # Get the global phase
+    gphase_op = sd_ops.pop()
+
+    # Decompose the rotation gates
+    d_ops = []
+    for sd_op in sd_ops:
+        d_ops.extend(_rot_decompose(sd_op) if sd_op.num_params else [sd_op])
+
+    return d_ops, gphase_op
+
+
+def _two_qubit_decompose(op):
+    r"""Decomposition for two qubit operations using combination of :class:`~.RZ`, :class:`~.S`,
+    :class:`~.Hadamard`, and :class:`~.CNOT`."""
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", RuntimeWarning)
+        td_ops = qml.transforms.two_qubit_decomposition(qml.matrix(op), op.wires)
+
+    d_ops = []
+    for td_op in td_ops:
+        d_ops.extend(
+            _rot_decompose(td_op) if td_op.num_params and td_op.num_wires == 1 else [td_op]
+        )
+
+    return d_ops
+
+
+def _merge_param_gates(operations, merge_ops=None):
+    """Merge the provided parameterized gates on the same wires that are adjacent to each other"""
+
+    copied_ops = operations.copy()
+    merged_ops, number_ops = [], 0
+
+    while len(copied_ops) > 0:
+        curr_gate = copied_ops.pop(0)
+
+        # If gate is not to be merged, let it go
+        if merge_ops is not None and curr_gate.name not in merge_ops:
+            merged_ops.append(curr_gate)
+            continue
+
+        # If gate is in the merge_ops, update counter
+        if curr_gate.name in merge_ops:
+            number_ops += 1
+
+        # Find the next gate that acts on the same wires
+        next_gate_idx = find_next_gate(curr_gate.wires, copied_ops)
+        if next_gate_idx is None:
+            merged_ops.append(curr_gate)
+            continue
+
+        # Initialize the current angle and iterate until end of merge
+        curr_params = curr_gate.parameters
+        curr_intrfc = qml.math.get_deep_interface(curr_gate.parameters)
+        cumulative_angles = qml.math.array(curr_params, dtype=float, like=curr_intrfc)
+        next_gate = copied_ops[next_gate_idx]
+        while curr_gate.name == next_gate.name and curr_gate.wires == next_gate.wires:
+            cumulative_angles += qml.math.array(next_gate.parameters, like=curr_intrfc)
+            # Check if the subsequent gate exists in the vicinity
+            copied_ops.pop(next_gate_idx)
+            next_gate_idx = find_next_gate(curr_gate.wires, copied_ops)
+            if next_gate_idx is None:
+                break
+            next_gate = copied_ops[next_gate_idx]
+
+        # Replace the current gate, add it to merged list and pop it from orginal one
+        merged_ops.append(curr_gate.__class__(*cumulative_angles, wires=curr_gate.wires))
+
+    return merged_ops, number_ops
+
+
+# pylint: disable= too-many-nested-blocks, too-many-branches, too-many-statements, unnecessary-lambda-assignment
+@transform
+def clifford_t_decomposition(
+    tape: QuantumTape,
+    epsilon=1e-4,
+    max_expansion=6,
+    method="sk",
+    **method_kwargs,
+) -> (Sequence[QuantumTape], Callable):
+    r"""Decomposes a circuit into the Clifford+T basis.
+
+    This method first decomposes the gate operations to a basis comprised of Clifford, :class:`~.T`, :class:`~.RZ` and
+    :class:`~.GlobalPhase` operations (and their adjoints). The Clifford gates include the following PennyLane operations:
+
+    - Single qubit gates - :class:`~.Identity`, :class:`~.PauliX`, :class:`~.PauliY`, :class:`~.PauliZ`,
+      :class:`~.SX`, :class:`~.S`, and :class:`~.Hadamard`.
+    - Two qubit gates - :class:`~.CNOT`, :class:`~.CY`, :class:`~.CZ`, :class:`~.SWAP`, and :class:`~.ISWAP`.
+
+    Then, the leftover single qubit :class:`~.RZ` operations are approximated in the Clifford+T basis with
+    :math:`\epsilon > 0` error. By default, we use the Solovay-Kitaev algorithm described in
+    `Dawson and Nielsen (2005) <https://arxiv.org/abs/quant-ph/0505030>`_ for this.
+
+    Args:
+        tape (QNode or QuantumTape or Callable): The quantum circuit to be decomposed.
+        epsilon (float): The maximum permissible operator norm error of the complete circuit decomposition. Defaults to ``0.0001``.
+        max_expansion (int): The depth to be used for tape expansion before manual decomposition to Clifford+T basis is applied.
+        method (str): Method to be used for Clifford+T decomposition. Default value is ``"sk"`` for Solovay-Kitaev.
+        **method_kwargs: Keyword argument to pass options for the ``method`` used for decompositions.
+
+    Returns:
+        qnode (QNode) or quantum function (Callable) or tuple[List[QuantumTape], function]: The transformed circuit as described
+        in the :func:`qml.transform <pennylane.transform>`.
+
+    **Keyword Arguments**
+
+    - Solovay-Kitaev decomposition --
+        **max_depth** (int), **basis_set** (list[str]), **basis_length** (int) -- arguments for the ``"sk"`` method,
+        where the decomposition is performed using the :func:`~.sk_decomposition` method.
+
+    Raises:
+        ValueError: If a gate operation does not have a decomposition when required.
+        NotImplementedError: If chosen decomposition ``method`` is not supported.
+
+    .. seealso:: :func:`~.sk_decomposition` for Solovay-Kitaev decomposition.
+
+    **Example**
+
+    .. code-block:: python3
+
+        @qml.qnode(qml.device("default.qubit"))
+        def circuit(x, y):
+            qml.RX(x, 0)
+            qml.CNOT([0, 1])
+            qml.RY(y, 0)
+            return qml.expval(qml.PauliZ(0))
+
+        x, y = 1.1, 2.2
+        decomposed_circuit = qml.transforms.clifford_t_decomposition(circuit)
+        result = circuit(x, y)
+        approx = decomposed_circuit(x, y)
+
+    >>> qml.math.allclose(result, approx, atol=1e-4)
+    True
+    """
+    with QueuingManager.stop_recording():
+        # Build the basis set and the pipeline for intial compilation pass
+        basis_set = [op.__name__ for op in _PARAMETER_GATES + _CLIFFORD_T_GATES]
+        pipelines = [remove_barrier, commute_controlled, cancel_inverses, merge_rotations]
+
+        # Compile the tape according to depth provided by the user and expand it
+        [compiled_tape], _ = qml.compile(
+            tape, pipelines, basis_set=basis_set, expand_depth=max_expansion
+        )
+
+        # Now iterate over the expanded tape operations
+        decomp_ops, gphase_ops = [], []
+        for op in compiled_tape.operations:
+            # Check whether operation is to be skipped
+            if isinstance(op, _SKIP_OP_TYPES):
+                decomp_ops.append(op)
+
+            # Check whether the operation is a global phase
+            elif isinstance(op, qml.GlobalPhase):
+                gphase_ops.append(op)
+
+            # Check whether the operation is a Clifford or a T-gate
+            elif op.name in basis_set and check_clifford_t(op):
+                if op.num_params:
+                    decomp_ops.extend(_rot_decompose(op))
+                else:
+                    decomp_ops.append(op)
+
+            # Decompose and then iteratively go deeper via DFS
+            else:
+                # Single qubit unitary decomposition with ZXZ rotations
+                if op.num_wires == 1:
+                    if op.name in basis_set:
+                        d_ops = _rot_decompose(op)
+                    else:
+                        d_ops, g_op = _one_qubit_decompose(op)
+                        gphase_ops.append(g_op)
+                    decomp_ops.extend(d_ops)
+
+                # Two qubit unitary decomposition with SU(4) rotations
+                elif op.num_wires == 2:
+                    d_ops = _two_qubit_decompose(op)
+                    decomp_ops.extend(d_ops)
+
+                # For special multi-qubit gates and ones constructed from matrix
+                else:
+                    try:
+                        # Attempt decomposing the operation
+                        md_ops = op.decomposition()
+                        idx = 0  # might not be fast but at least is not recursive
+                        while idx < len(md_ops):
+                            md_op = md_ops[idx]
+                            if md_op.name not in basis_set or not check_clifford_t(md_op):
+                                # For the gates acting on one qubit
+                                if len(md_op.wires) == 1:
+                                    if md_op.name in basis_set:  # For known recipe
+                                        d_ops = _rot_decompose(md_op)
+                                    else:  # Resort to decomposing manually
+                                        d_ops, g_op = _one_qubit_decompose(md_op)
+                                        gphase_ops.append(g_op)
+
+                                # For the gates acting on two qubits
+                                elif len(md_op.wires) == 2:
+                                    # Resort to decomposing manually
+                                    d_ops = _two_qubit_decompose(md_op)
+
+                                # Final resort (should not enter in an ideal situtation)
+                                else:
+                                    d_ops = md_op.decomposition()
+
+                                # Expand the list and iterate over
+                                del md_ops[idx]
+                                md_ops[idx:idx] = d_ops
+                            idx += 1
+
+                        decomp_ops.extend(md_ops)
+
+                    # If we don't know how to decompose the operation
+                    except Exception as exc:
+                        raise ValueError(
+                            f"Cannot unroll {op} into the Clifford+T basis as no rule exists for its decomposition"
+                        ) from exc
+
+        # Merge RZ rotations together
+        merged_ops, number_ops = _merge_param_gates(decomp_ops, merge_ops=["RZ"])
+
+        # Squeeze global phases into a single global phase
+        new_operations = _fuse_global_phases(merged_ops + gphase_ops)
+
+        # Compute the per-gate epsilon value
+        epsilon /= number_ops or 1
+
+        # Every decomposition implementation should have the following shape:
+        # def decompose_fn(op: Operator, epsilon: float, **method_kwargs) -> List[Operator]
+        # note: the last operator in the decomposition must be a GlobalPhase
+
+        # Build the approximation set for Solovay-Kitaev decomposition
+        if method == "sk":
+            decompose_fn = sk_decomposition
+
+        else:
+            raise NotImplementedError(
+                f"Currently we only support Solovay-Kitaev ('sk') decompostion, got {method}"
+            )
+
+        decomp_ops = []
+        phase = new_operations.pop().data[0]
+        for op in new_operations:
+            if isinstance(op, qml.RZ):
+                clifford_ops = decompose_fn(op, epsilon, **method_kwargs)
+                phase += qml.math.convert_like(clifford_ops.pop().data[0], phase)
+                decomp_ops.extend(clifford_ops)
+            else:
+                decomp_ops.append(op)
+
+        # check if phase is non-zero for non jax-jit cases
+        if qml.math.is_abstract(phase) or not qml.math.allclose(phase, 0.0):
+            decomp_ops.append(qml.GlobalPhase(phase))
+
+    # Construct a new tape with the expanded set of operations
+    new_tape = type(tape)(decomp_ops, compiled_tape.measurements, shots=tape.shots)
+
+    # Perform a final attempt of simplification before return
+    [new_tape], _ = cancel_inverses(new_tape)
+
+    def null_postprocessing(results):
+        """A postprocesing function returned by a transform that only converts the batch of results
+        into a result for a single ``QuantumTape``.
+        """
+        return results[0]
+
+    return [new_tape], null_postprocessing
diff --git a/pennylane/transforms/decompositions/clifford_t/solovay_kitaev.py b/pennylane/transforms/decompositions/clifford_t/solovay_kitaev.py
index 87d3fa655..4e00731a0 100644
--- a/pennylane/transforms/decompositions/clifford_t/solovay_kitaev.py
+++ b/pennylane/transforms/decompositions/clifford_t/solovay_kitaev.py
@@ -277,7 +277,9 @@ def sk_decomposition(op, epsilon, *, max_depth=5, basis_set=("T", "T*", "H"), ba
         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())
+        op_matrix = op.matrix()
+        interface = qml.math.get_deep_interface(op_matrix)
+        gate_mat, gate_gph = _SU2_transform(qml.math.unwrap(op_matrix))
         gate_qat = _quaternion_transform(gate_mat)
 
         def _solovay_kitaev(umat, n, u_n1_ids, u_n1_mat):
@@ -339,9 +341,8 @@ def sk_decomposition(op, epsilon, *, max_depth=5, basis_set=("T", "T*", "H"), ba
     [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
-    )
+    phase = approx_set_gph[index] - gate_gph if depth or qml.math.allclose(gate_gph, 0.0) else 0.0
+    global_phase = qml.GlobalPhase(qml.math.array(phase, like=interface))
 
     # Return the gates from the mapped tape and global phase
     return map_tape.operations + [global_phase]
diff --git a/pennylane/transforms/optimization/merge_rotations.py b/pennylane/transforms/optimization/merge_rotations.py
index 91cb67efa..883c47be0 100644
--- a/pennylane/transforms/optimization/merge_rotations.py
+++ b/pennylane/transforms/optimization/merge_rotations.py
@@ -17,7 +17,7 @@ from typing import Sequence, Callable
 
 from pennylane.tape import QuantumTape
 from pennylane.transforms import transform
-from pennylane.math import allclose, stack, cast_like, zeros, is_abstract
+from pennylane.math import allclose, stack, cast_like, zeros, is_abstract, get_interface
 from pennylane.queuing import QueuingManager
 
 from pennylane.ops.qubit.attributes import composable_rotations
@@ -149,6 +149,7 @@ def merge_rotations(
 
         # We need to use stack to get this to work and be differentiable in all interfaces
         cumulative_angles = stack(current_gate.parameters)
+        interface = get_interface(cumulative_angles)
         # As long as there is a valid next gate, check if we can merge the angles
         while next_gate_idx is not None:
             # Get the next gate
@@ -165,7 +166,9 @@ def merge_rotations(
                     else:
                         cumulative_angles = fuse_rot_angles(
                             cumulative_angles,
-                            cast_like(stack(next_gate.parameters), cumulative_angles),
+                            cast_like(
+                                stack(next_gate.parameters, like=interface), cumulative_angles
+                            ),
                         )
                 # Other, single-parameter rotation gates just have the angle summed
                 else:
@@ -174,7 +177,7 @@ def merge_rotations(
                         cumulative_angles = cumulative_angles + stack(next_gate.parameters)
                     else:
                         cumulative_angles = cumulative_angles + cast_like(
-                            stack(next_gate.parameters), cumulative_angles
+                            stack(next_gate.parameters, like=interface), cumulative_angles
                         )
             # If it is not, we need to stop
             else:
diff --git a/pennylane/transforms/optimization/optimization_utils.py b/pennylane/transforms/optimization/optimization_utils.py
index 11e5af4de..b6727b507 100644
--- a/pennylane/transforms/optimization/optimization_utils.py
+++ b/pennylane/transforms/optimization/optimization_utils.py
@@ -18,8 +18,10 @@ from functools import partial
 import numpy as np
 
 from pennylane.math import abs as math_abs
+from pennylane.math import sum as math_sum
 from pennylane.math import allclose, arccos, arctan2, cos, get_interface, is_abstract, sin, stack
 from pennylane.wires import Wires
+from pennylane.ops.identity import GlobalPhase
 
 
 def find_next_gate(wires, op_list):
@@ -185,3 +187,24 @@ def fuse_rot_angles(angles_1, angles_2):
         return _no_fuse(angles_1, angles_2)
 
     return _fuse(angles_1, angles_2)
+
+
+def _fuse_global_phases(operations):
+    """Fuse all the global phase operations into single one
+
+    Args:
+        operations (list[Operation]): list of operations to be iterated over
+
+    Returns:
+        transformed list with a single :func:`~.pennylane.GlobalPhase` operation.
+    """
+
+    fused_ops, global_ops = [], []
+    for op in operations:
+        if isinstance(op, GlobalPhase):
+            global_ops.append(op)
+        else:
+            fused_ops.append(op)
+
+    fused_ops.append(GlobalPhase(math_sum(op.data[0] for op in global_ops)))
+    return fused_ops


Test patch:
diff --git a/tests/transforms/test_clifford_t/test_cliffordt_transform.py b/tests/transforms/test_clifford_t/test_cliffordt_transform.py
new file mode 100644
index 000000000..399c985f1
--- /dev/null
+++ b/tests/transforms/test_clifford_t/test_cliffordt_transform.py
@@ -0,0 +1,377 @@
+# 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.
+"""Unit tests for the Clifford+T transform."""
+
+import math
+from functools import reduce
+
+import pytest
+import pennylane as qml
+
+from pennylane.transforms.decompositions.clifford_t.clifford_t_transform import (
+    check_clifford_t,
+    clifford_t_decomposition,
+    _rot_decompose,
+    _one_qubit_decompose,
+    _two_qubit_decompose,
+    _CLIFFORD_T_GATES,
+)
+
+from pennylane.transforms.optimization.optimization_utils import _fuse_global_phases
+
+_SKIP_GATES = (qml.Barrier, qml.Snapshot, qml.WireCut)
+_CLIFFORD_PHASE_GATES = _CLIFFORD_T_GATES + _SKIP_GATES
+
+INVSQ2 = 1 / math.sqrt(2)
+PI = math.pi
+
+
+def circuit_1():
+    """Circuit 1 with quantum chemistry gates"""
+    qml.RZ(1.0, wires=[0])
+    qml.PhaseShift(1.0, wires=[1])
+    qml.SingleExcitation(2.0, wires=[1, 2])
+    qml.PauliX(0)
+    return qml.expval(qml.PauliZ(1))
+
+
+def circuit_2():
+    """Circuit 2 without chemistry gates"""
+    qml.CRX(1, wires=[0, 1])
+    qml.ISWAP(wires=[0, 1])
+    qml.CSWAP(wires=[0, 1, 2])
+    return qml.expval(qml.PauliZ(0))
+
+
+def circuit_3():
+    """Circuit 3 with Clifford gates"""
+    qml.GlobalPhase(PI)
+    qml.CNOT(wires=[0, 1])
+    qml.PauliX(wires=[1])
+    qml.ISWAP(wires=[0, 1])
+    qml.Hadamard(wires=[0])
+    qml.WireCut(wires=[1])
+    qml.RZ(PI, wires=[0])
+    return qml.expval(qml.PauliZ(0))
+
+
+def circuit_4():
+    """Circuit 4 with a Template"""
+    qml.RandomLayers(weights=qml.math.array([[0.1, -2.1, 1.4]]), wires=range(2))
+    return qml.expval(qml.PauliZ(0))
+
+
+def circuit_5():
+    """Circuit 5 with Qubit Unitaries"""
+    qml.DiagonalQubitUnitary(
+        [qml.math.exp(1j * 0.1), qml.math.exp(1j * PI), INVSQ2 * (1 + 1j), INVSQ2 * (1 - 1j)],
+        wires=[0, 1],
+    )
+    return qml.expval(qml.PauliZ(0))
+
+
+class TestCliffordCompile:
+    """Unit tests for clifford compilation function."""
+
+    @pytest.mark.parametrize(
+        "op, res",
+        [
+            (qml.DoubleExcitation(2.0, wires=[0, 1, 2, 3]), False),
+            (qml.PauliX(wires=[1]), True),
+            (qml.RX(3 * PI, wires=[1]), True),
+            (qml.PhaseShift(2 * PI, wires=["a"]), True),
+            (qml.ECR(wires=["e", "f"]), True),
+            (qml.CH(wires=["a", "b"]), False),
+            (qml.WireCut(0), False),
+        ],
+    )
+    def test_clifford_checker(self, op, res):
+        """Test Clifford checker operation for gate"""
+        assert check_clifford_t(op) == res
+        assert check_clifford_t(op, use_decomposition=True) == res
+
+    @pytest.mark.parametrize(
+        ("circuit, max_expansion"),
+        [(circuit_1, 1), (circuit_2, 0), (circuit_3, 0), (circuit_4, 1), (circuit_5, 0)],
+    )
+    def test_decomposition(self, circuit, max_expansion):
+        """Test decomposition for the Clifford transform."""
+
+        old_tape = qml.tape.make_qscript(circuit)()
+
+        [new_tape], tape_fn = clifford_t_decomposition(
+            old_tape, max_expansion=max_expansion, max_depth=3
+        )
+
+        assert all(
+            isinstance(op, _CLIFFORD_PHASE_GATES)
+            or isinstance(getattr(op, "base", None), _CLIFFORD_PHASE_GATES)
+            for op in new_tape.operations
+        )
+
+        dev = qml.device("default.qubit")
+        transform_program, _ = dev.preprocess()
+        res1, res2 = qml.execute(
+            [old_tape, new_tape], device=dev, transform_program=transform_program
+        )
+        qml.math.isclose(res1, tape_fn([res2]), atol=1e-2)
+
+    def test_qnode_decomposition(self):
+        """Test decomposition for the Clifford transform."""
+
+        dev = qml.device("default.qubit")
+
+        def qfunc():
+            qml.PhaseShift(1.0, wires=[0])
+            qml.PhaseShift(2.0, wires=[1])
+            qml.ISWAP(wires=[0, 1])
+            return qml.expval(qml.PauliZ(0))
+
+        original_qnode = qml.QNode(qfunc, dev)
+        transfmd_qnode = qml.QNode(
+            clifford_t_decomposition(qfunc, max_depth=3, basis_length=10), dev
+        )
+
+        res1, res2 = original_qnode(), transfmd_qnode()
+        assert qml.math.isclose(res1, res2, atol=1e-2)
+
+        assert all(
+            isinstance(op, _CLIFFORD_PHASE_GATES)
+            or isinstance(getattr(op, "base", None), _CLIFFORD_PHASE_GATES)
+            for op in transfmd_qnode.tape.operations
+        )
+
+    @pytest.mark.parametrize("epsilon", [2e-2, 5e-2, 9e-2])
+    @pytest.mark.parametrize("circuit", [circuit_3, circuit_4, circuit_5])
+    def test_total_error(self, epsilon, circuit):
+        """Ensure that given a certain epsilon, the total operator error is below the threshold."""
+        dev = qml.device("default.qubit")
+
+        qnode_basic = qml.QNode(circuit, dev)
+        qnode_transformed = clifford_t_decomposition(qnode_basic, epsilon=epsilon, max_depth=10)
+        mat_exact = qml.matrix(qnode_basic, wire_order=[0, 1])()
+        mat_approx = qml.matrix(qnode_transformed, wire_order=[0, 1])()
+        phase = qml.math.divide(
+            mat_exact,
+            mat_approx,
+            out=qml.math.zeros_like(mat_exact, dtype=complex),
+            where=mat_exact != 0,
+        )[qml.math.nonzero(qml.math.round(mat_exact, 10))][0]
+        mat_exact /= phase
+        diff = mat_exact - mat_approx
+        error = qml.math.sqrt(qml.math.real(qml.math.trace(qml.math.conj(diff).T @ diff)) / 2)
+        assert error < epsilon
+
+    @pytest.mark.parametrize(
+        "op", [qml.RX(1.0, wires="a"), qml.U3(1, 2, 3, wires=[1]), qml.PhaseShift(1.0, wires=[2])]
+    )
+    def test_one_qubit_decomposition(self, op):
+        """Test decomposition for the Clifford transform."""
+
+        decomp_ops, global_ops = _one_qubit_decompose(op)
+        decomp_ops = _fuse_global_phases(decomp_ops + [global_ops])
+
+        valid_gates = _CLIFFORD_PHASE_GATES + (qml.RZ,)
+        assert all(
+            isinstance(op, valid_gates) or isinstance(getattr(op, "base", None), valid_gates)
+            for op in decomp_ops
+        )
+
+        global_op = decomp_ops.pop()
+        assert isinstance(global_op, qml.GlobalPhase)
+
+        matrix_op = reduce(
+            lambda x, y: x @ y, [qml.matrix(op) for op in decomp_ops][::-1]
+        ) * qml.matrix(global_op)
+
+        # check for matrice equivalence up to global phase
+        phase = qml.math.divide(
+            matrix_op,
+            qml.matrix(op),
+            out=qml.math.zeros_like(matrix_op, dtype=complex),
+            where=matrix_op != 0,
+        )[qml.math.nonzero(qml.math.round(matrix_op, 10))]
+        assert qml.math.allclose(
+            phase / phase[0], qml.math.ones(qml.math.shape(phase)[0]), atol=1e-5
+        )
+
+    @pytest.mark.parametrize(
+        "op",
+        [
+            qml.PSWAP(1.0, wires=["a", "b"]),
+            qml.SingleExcitation(1, wires=[1, 2]),
+            qml.IsingXX(1.0, wires=[2, 3]),
+        ],
+    )
+    def test_two_qubit_decomposition(self, op):
+        """Test decomposition for the Clifford transform."""
+
+        decomp_ops = _fuse_global_phases(_two_qubit_decompose(op))
+
+        valid_gates = _CLIFFORD_PHASE_GATES + (qml.RZ,)
+        assert all(
+            isinstance(op, valid_gates) or isinstance(getattr(op, "base", None), valid_gates)
+            for op in decomp_ops
+        )
+
+        global_op = decomp_ops.pop()
+        wire_map = {wire: idx for idx, wire in enumerate(op.wires)}
+        mapped_op = [qml.map_wires(op, wire_map=wire_map) for op in decomp_ops][::-1]
+        matrix_op = reduce(
+            lambda x, y: x @ y, [qml.matrix(op, wire_order=[0, 1]) for op in mapped_op]
+        ) * qml.matrix(global_op, wire_order=[0, 1])
+
+        # check for matrice equivalence up to global phase
+        phase = qml.math.divide(
+            matrix_op,
+            qml.matrix(op),
+            out=qml.math.zeros_like(matrix_op, dtype=complex),
+            where=matrix_op != 0,
+        )[qml.math.nonzero(qml.math.round(matrix_op, 10))]
+        assert qml.math.allclose(phase / phase[0], qml.math.ones(qml.math.shape(phase)[0]))
+
+    @pytest.mark.parametrize(
+        "op",
+        [
+            qml.adjoint(qml.RX(1.0, wires=["b"])),
+            qml.Rot(1, 2, 3, wires=[2]),
+            qml.PhaseShift(1.0, wires=[0]),
+            qml.PhaseShift(3 * PI, wires=[0]),
+        ],
+    )
+    def test_rot_decomposition(self, op):
+        """Test decomposition for the Clifford transform."""
+
+        decomp_ops = _fuse_global_phases(_rot_decompose(op))
+
+        valid_gates = _CLIFFORD_PHASE_GATES + (qml.RZ,)
+        assert all(
+            isinstance(op, valid_gates) or isinstance(getattr(op, "base", None), valid_gates)
+            for op in decomp_ops
+        )
+
+        decomp_ops, global_ops = decomp_ops[:-1], decomp_ops[-1]
+        matrix_op = reduce(
+            lambda x, y: x @ y, [qml.matrix(op) for op in decomp_ops][::-1]
+        ) * qml.matrix(global_ops)
+
+        # check for matrice equivalence up to global phase
+        phase = qml.math.divide(
+            matrix_op,
+            qml.matrix(op),
+            out=qml.math.zeros_like(matrix_op, dtype=complex),
+            where=matrix_op != 0,
+        )[qml.math.nonzero(qml.math.round(matrix_op, 10))]
+        assert qml.math.allclose(phase / phase[0], qml.math.ones(qml.math.shape(phase)[0]))
+
+    def test_raise_with_cliffordt_decomposition(self):
+        """Test that exception is correctly raise when decomposing gates without any decomposition"""
+
+        tape = qml.tape.QuantumScript([qml.QubitUnitary(qml.math.eye(8), wires=[0, 1, 2])])
+
+        with pytest.raises(ValueError, match="Cannot unroll"):
+            clifford_t_decomposition(tape, max_expansion=0)
+
+    @pytest.mark.parametrize("op", [qml.U1(1.0, wires=["b"])])
+    def test_raise_with_rot_decomposition(self, op):
+        """Test that exception is correctly raise when decomposing parameterized gates for which we already don't have a recipe"""
+
+        with pytest.raises(
+            ValueError,
+            match="qml.RX, qml.RY, qml.RZ, qml.Rot and qml.PhaseShift",
+        ):
+            _rot_decompose(op)
+
+    def test_zero_global_phase(self):
+        """Test that global phase operation is added only when it is non-zero"""
+
+        with qml.tape.QuantumTape() as tape:
+            qml.CNOT([0, 1])
+
+        [tape], _ = qml.clifford_t_decomposition(tape)
+
+        assert not sum([isinstance(op, qml.GlobalPhase) for op in tape.operations])
+
+    def test_raise_with_decomposition_method(self):
+        """Test that exception is correctly raise when using incorrect decomposing method"""
+
+        def qfunc():
+            qml.RX(1.0, wires=[0])
+            return qml.expval(qml.PauliZ(0))
+
+        decomposed_qfunc = clifford_t_decomposition(qfunc, method="fast")
+
+        with pytest.raises(
+            NotImplementedError,
+            match=r"Currently we only support Solovay-Kitaev \('sk'\) decompostion",
+        ):
+            decomposed_qfunc()
+
+    # pylint: disable= import-outside-toplevel
+    @pytest.mark.all_interfaces
+    def test_clifford_decompose_interfaces(self):
+        """Test that unwrap converts lists to lists and interface variables to numpy."""
+
+        dev = qml.device("default.qubit", wires=3)
+
+        def circuit(x):
+            qml.RZ(x[0], wires=[0])
+            qml.PhaseShift(x[1], wires=[1])
+            qml.SingleExcitation(x[2], wires=[1, 2])
+            qml.PauliX(0)
+            return qml.expval(qml.PauliZ(1))
+
+        original_qnode = qml.QNode(circuit, dev)
+        transfmd_qnode = qml.QNode(
+            clifford_t_decomposition(circuit, max_depth=3, basis_length=10), dev
+        )
+
+        import jax
+        import torch
+        import tensorflow as tf
+
+        funres = []
+        igrads = []
+        coeffs = [1.0, 2.0, 3.0]
+        for qcirc in [original_qnode, transfmd_qnode]:
+            # Autograd Interface
+            A = qml.numpy.array(coeffs)
+            fres_numpy = qcirc(A)
+            grad_numpy = qml.grad(qcirc, argnum=0)(A)
+
+            # Jax Interface
+            A = jax.numpy.array(coeffs)
+            fres_jax = qcirc(A)
+            grad_jax = jax.grad(qcirc, argnums=0)(A)
+
+            # Tensorflow Interface
+            A = tf.Variable(qml.numpy.array(coeffs))
+            with tf.GradientTape() as tape:
+                loss = qcirc(A)
+            grad_tflow = tape.gradient(loss, A)
+            fres_tflow = loss
+
+            # PyTorch Interface
+            A = torch.tensor(coeffs, requires_grad=True)
+            result = qcirc(A)
+            result.backward()
+            grad_torch = A.grad
+            fres_torch = result
+
+            funres.append([fres_numpy, fres_jax, fres_torch, fres_tflow])
+            igrads.append([grad_numpy, grad_jax, grad_torch, grad_tflow])
+
+        # Compare results
+        assert all(qml.math.allclose(res1, res2, atol=1e-2) for res1, res2 in zip(*funres))
+        assert all(qml.math.allclose(res1, res2, atol=1e-2) for res1, res2 in zip(*igrads))