Repository: PennyLaneAI/pennylane
Pull request #1458 (2021-07-15 14:32:46)

Implement transform for full single-qubit gate fusion
**Context:** Addition of circuit optimization / compilation tools to PennyLane. 

**Description of the Change:** Adds a new transform that merges arbitrary-length sequences of adjacent single-qubit gates. Adds related tests, and an additional attribute to express single-qubit gates as three rotation angles.

**Benefits:** Can significantly reduce the length of circuits by squashing any number of single-qubit gates down to one (technically, three, but a single `Rot`).

**Possible Drawbacks:** None

**Related GitHub Issues:** None


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [.github/CHANGELOG.md](/PennyLaneAI/pennylane/blob/master/.github/CHANGELOG.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.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1458](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (20198e7) into [ch7104-implement-transforms-to-cancel-redundant](https://codecov.io/gh/PennyLaneAI/pennylane/commit/7ce969fcf8e0cf4dc2c00a80c702e99088856c4f?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (7ce969f) will **decrease** coverage by `0.02%`.
> The diff coverage is `87.50%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@                                 Coverage Diff                                 @@
##           ch7104-implement-transforms-to-cancel-redundant    #1458      +/-   ##
===================================================================================
- Coverage                                            98.24%   98.22%   -0.03%     
===================================================================================
  Files                                                  167      168       +1     
  Lines                                                12245    12308      +63     
===================================================================================
+ Hits                                                 12030    12089      +59     
- Misses                                                 215      219       +4     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/operation.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wZXJhdGlvbi5weQ==) | `95.77% <66.66%> (-0.19%)` | :arrow_down: |
| [pennylane/ops/qubit.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC5weQ==) | `98.37% <83.33%> (-0.32%)` | :arrow_down: |
| [...ane/transforms/optimization/single\_qubit\_fusion.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvb3B0aW1pemF0aW9uL3NpbmdsZV9xdWJpdF9mdXNpb24ucHk=) | `91.42% <91.42%> (ø)` | |
| [pennylane/transforms/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvX19pbml0X18ucHk=) | `100.00% <100.00%> (ø)` | |
| [pennylane/transforms/optimization/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvb3B0aW1pemF0aW9uL19faW5pdF9fLnB5) | `100.00% <100.00%> (ø)` | |
| [...lane/transforms/optimization/optimization\_utils.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvb3B0aW1pemF0aW9uL29wdGltaXphdGlvbl91dGlscy5weQ==) | `97.91% <0.00%> (+8.33%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [7ce969f...20198e7](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1458?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

glassnotes: > Looks great @glassnotes! I think the only thing still missing are tests for the single_qubit_rot_angles methods added to all the operations. Is there a clean way to easily test that they have the correct angles hardcoded?
>
> This is a small thing that could easily be hiding bugs

@josh146 thanks for pointing that out, have added a small test for this in the `test_qubit_ops` file! Think everything else should be covered now with the most recent updates to tests.

Patch:
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index c546152b6..5c885d81f 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -2,6 +2,37 @@
 
 <h3>New features since last release</h3>
 
+* A new quantum function transform has been added to perform full fusion of
+  arbitrary-length sequences of single-qubit gates.
+  [(#1458)](https://github.com/PennyLaneAI/pennylane/pull/1458)
+
+  The `single_qubit_fusion` transform acts on all sequences of
+  single-qubit operations in a quantum function, and converts each
+  sequence to a single `Rot` gate. For example given the circuit:
+
+  ```python
+  def circuit(x, y, z):
+      qml.Hadamard(wires=0)
+      qml.PauliZ(wires=1)
+      qml.RX(x, wires=0)
+      qml.RY(y, wires=1)
+      qml.CZ(wires=[1, 0])
+      qml.T(wires=0)
+      qml.SX(wires=0)
+      qml.Rot(x, y, z, wires=1)
+      qml.Rot(z, y, x, wires=1)
+      return qml.expval(qml.PauliX(wires=0))
+  ```
+
+  ```pycon
+  >>> optimized_circuit = qml.transforms.single_qubit_fusion()(circuit)
+  >>> dev = qml.device('default.qubit', wires=2)
+  >>> qnode = qml.QNode(optimized_circuit, dev)
+  >>> print(qml.draw(qnode)(0.1, 0.2, 0.3))
+   0: ──Rot(3.24, 1.57, 0)──╭Z──Rot(2.36, 1.57, -1.57)────┤ ⟨X⟩
+   1: ──Rot(3.14, 0.2, 0)───╰C──Rot(0.406, 0.382, 0.406)──┤
+  ```
+
 * Two new quantum function transforms have been added to enable the
   removal of redundant gates in quantum circuits.
   [(#1455)](https://github.com/PennyLaneAI/pennylane/pull/1455)
diff --git a/pennylane/operation.py b/pennylane/operation.py
index e246ca02c..1332de2c8 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -582,6 +582,18 @@ class Operation(Operator):
     transforms that merge adjacent rotations.
     """
 
+    @property
+    def single_qubit_rot_angles(self):
+        r"""The parameters required to implement a single-qubit gate as an
+        equivalent ``Rot`` gate, up to a global phase.
+
+        Returns:
+            tuple[float, float, float]: A list of values :math:`[\phi, \theta, \omega]`
+            such that :math:`RZ(\omega) RY(\theta) RZ(\phi)` is equivalent to the
+            original operation.
+        """
+        raise NotImplementedError
+
     def get_parameter_shift(self, idx, shift=np.pi / 2):
         """Multiplier and shift for the given parameter, based on its gradient recipe.
 
diff --git a/pennylane/ops/qubit.py b/pennylane/ops/qubit.py
index 22ae645cf..9434ddccc 100644
--- a/pennylane/ops/qubit.py
+++ b/pennylane/ops/qubit.py
@@ -98,6 +98,10 @@ class Hadamard(Observable, Operation):
     def adjoint(self):
         return Hadamard(wires=self.wires)
 
+    def single_qubit_rot_angles(self):
+        # H = RZ(\pi) RY(\pi/2) RZ(0)
+        return [np.pi, np.pi / 2, 0.0]
+
 
 class PauliX(Observable, Operation):
     r"""PauliX(wires)
@@ -157,6 +161,10 @@ class PauliX(Observable, Operation):
     def _controlled(self, wire):
         CNOT(wires=Wires(wire) + self.wires)
 
+    def single_qubit_rot_angles(self):
+        # X = RZ(-\pi/2) RY(\pi) RZ(\pi/2)
+        return [np.pi / 2, np.pi, -np.pi / 2]
+
 
 class PauliY(Observable, Operation):
     r"""PauliY(wires)
@@ -218,6 +226,10 @@ class PauliY(Observable, Operation):
     def _controlled(self, wire):
         CY(wires=Wires(wire) + self.wires)
 
+    def single_qubit_rot_angles(self):
+        # Y = RZ(0) RY(\pi) RZ(0)
+        return [0.0, np.pi, 0.0]
+
 
 class PauliZ(Observable, DiagonalOperation):
     r"""PauliZ(wires)
@@ -262,6 +274,10 @@ class PauliZ(Observable, DiagonalOperation):
     def _controlled(self, wire):
         CZ(wires=Wires(wire) + self.wires)
 
+    def single_qubit_rot_angles(self):
+        # Z = RZ(\pi) RY(0) RZ(0)
+        return [np.pi, 0.0, 0.0]
+
 
 class S(DiagonalOperation):
     r"""S(wires)
@@ -300,6 +316,10 @@ class S(DiagonalOperation):
     def adjoint(self):
         return S(wires=self.wires).inv()
 
+    def single_qubit_rot_angles(self):
+        # S = RZ(\pi/2) RY(0) RZ(0)
+        return [np.pi / 2, 0.0, 0.0]
+
 
 class T(DiagonalOperation):
     r"""T(wires)
@@ -338,6 +358,10 @@ class T(DiagonalOperation):
     def adjoint(self):
         return T(wires=self.wires).inv()
 
+    def single_qubit_rot_angles(self):
+        # T = RZ(\pi/4) RY(0) RZ(0)
+        return [np.pi / 4, 0.0, 0.0]
+
 
 class SX(Operation):
     r"""SX(wires)
@@ -381,6 +405,10 @@ class SX(Operation):
     def adjoint(self):
         return SX(wires=self.wires).inv()
 
+    def single_qubit_rot_angles(self):
+        # SX = RZ(-\pi/2) RY(\pi/2) RZ(\pi/2)
+        return [np.pi / 2, np.pi / 2, -np.pi / 2]
+
 
 class CNOT(Operation):
     r"""CNOT(wires)
@@ -773,6 +801,10 @@ class RX(Operation):
     def _controlled(self, wire):
         CRX(*self.parameters, wires=wire + self.wires)
 
+    def single_qubit_rot_angles(self):
+        # RX(\theta) = RZ(-\pi/2) RY(\theta) RZ(\pi/2)
+        return [np.pi / 2, self.data[0], -np.pi / 2]
+
 
 class RY(Operation):
     r"""RY(phi, wires)
@@ -815,6 +847,10 @@ class RY(Operation):
     def _controlled(self, wire):
         CRY(*self.parameters, wires=wire + self.wires)
 
+    def single_qubit_rot_angles(self):
+        # RY(\theta) = RZ(0) RY(\theta) RZ(0)
+        return [0.0, self.data[0], 0.0]
+
 
 class RZ(DiagonalOperation):
     r"""RZ(phi, wires)
@@ -863,6 +899,10 @@ class RZ(DiagonalOperation):
     def _controlled(self, wire):
         CRZ(*self.parameters, wires=wire + self.wires)
 
+    def single_qubit_rot_angles(self):
+        # RZ(\theta) = RZ(\theta) RY(0) RZ(0)
+        return [self.data[0], 0.0, 0.0]
+
 
 class PhaseShift(DiagonalOperation):
     r"""PhaseShift(phi, wires)
@@ -912,6 +952,10 @@ class PhaseShift(DiagonalOperation):
     def _controlled(self, wire):
         ControlledPhaseShift(*self.parameters, wires=wire + self.wires)
 
+    def single_qubit_rot_angles(self):
+        # PhaseShift(\theta) = RZ(\theta) RY(0) RZ(0)
+        return [self.data[0], 0.0, 0.0]
+
 
 class ControlledPhaseShift(DiagonalOperation):
     r"""ControlledPhaseShift(phi, wires)
@@ -1033,6 +1077,9 @@ class Rot(Operation):
     def _controlled(self, wire):
         CRot(*self.parameters, wires=wire + self.wires)
 
+    def single_qubit_rot_angles(self):
+        return self.data
+
 
 class MultiRZ(DiagonalOperation):
     r"""MultiRZ(theta, wires)
diff --git a/pennylane/transforms/__init__.py b/pennylane/transforms/__init__.py
index d643425d5..6817a1f48 100644
--- a/pennylane/transforms/__init__.py
+++ b/pennylane/transforms/__init__.py
@@ -47,6 +47,7 @@ containing quantum operations) that are used to construct QNodes.
     ~transforms.cancel_inverses
     ~ctrl
     ~transforms.merge_rotations
+    ~transforms.single_qubit_fusion
     ~transforms.unitary_to_rot
     ~transforms.invisible
     ~apply_controlled_Q
@@ -97,7 +98,7 @@ from .hamiltonian_expand import hamiltonian_expand
 from .invisible import invisible
 from .measurement_grouping import measurement_grouping
 from .metric_tensor import metric_tensor, metric_tensor_tape
-from .optimization import cancel_inverses, merge_rotations
+from .optimization import cancel_inverses, merge_rotations, single_qubit_fusion
 from .specs import specs
 from .qmc import apply_controlled_Q, quantum_monte_carlo
 from .unitary_to_rot import unitary_to_rot
diff --git a/pennylane/transforms/optimization/__init__.py b/pennylane/transforms/optimization/__init__.py
index adaff1307..67d4989ce 100644
--- a/pennylane/transforms/optimization/__init__.py
+++ b/pennylane/transforms/optimization/__init__.py
@@ -17,3 +17,4 @@ This subpackage contains quantum function transforms for optimizing quantum circ
 
 from .cancel_inverses import cancel_inverses
 from .merge_rotations import merge_rotations
+from .single_qubit_fusion import single_qubit_fusion
diff --git a/pennylane/transforms/optimization/optimization_utils.py b/pennylane/transforms/optimization/optimization_utils.py
index e62a133c6..870828cb0 100644
--- a/pennylane/transforms/optimization/optimization_utils.py
+++ b/pennylane/transforms/optimization/optimization_utils.py
@@ -41,7 +41,7 @@ def find_next_gate(wires, op_list):
 
 
 def _yzy_to_zyz(middle_yzy):
-    """Converts a set of angles representing a sequence of rotations RY, RZ, RY into
+    r"""Converts a set of angles representing a sequence of rotations RY, RZ, RY into
     an equivalent sequence of the form RZ, RY, RZ.
 
     Any rotation in 3-dimensional space (or, equivalently, any single-qubit unitary)
diff --git a/pennylane/transforms/optimization/single_qubit_fusion.py b/pennylane/transforms/optimization/single_qubit_fusion.py
new file mode 100644
index 000000000..8592fc725
--- /dev/null
+++ b/pennylane/transforms/optimization/single_qubit_fusion.py
@@ -0,0 +1,147 @@
+# Copyright 2018-2021 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 for fusing sequences of single-qubit gates."""
+# pylint: disable=too-many-branches
+from pennylane import apply
+from pennylane.transforms import qfunc_transform
+from pennylane.ops.qubit import Rot
+from pennylane.math import allclose, cast_like, stack, zeros
+
+from .optimization_utils import find_next_gate, fuse_rot_angles
+
+
+@qfunc_transform
+def single_qubit_fusion(tape, atol=1e-8, exclude_gates=None):
+    r"""Quantum function transform to fuse together groups of single-qubit
+    operations into a general single-qubit unitary operation (:class:`~.Rot`).
+
+    Fusion is performed only between gates that implement the property
+    ``single_qubit_rot_angles``. Any sequence of two or more single-qubit gates
+    (on the same qubit) with that property defined will be fused into one ``Rot``.
+
+    Args:
+        qfunc (function): A quantum function.
+        atol (float): An absolute tolerance for which to apply a rotation after
+            fusion. After fusion of gates, if the fused angles :math:`\theta` are such that
+            :math:`|\theta|\leq \text{atol}`, no rotation gate will be applied.
+        exclude_gates (None or list[str]): A list of gates that should be excluded
+            from full fusion. If set to ``None``, all single-qubit gates that can
+            be fused will be fused.
+
+    **Example**
+
+    Consider the following quantum function.
+
+    .. code-block:: python
+
+        def qfunc(r1, r2):
+            qml.Hadamard(wires=0)
+            qml.Rot(*r1, wires=0)
+            qml.Rot(*r2, wires=0)
+            qml.RZ(r1[0], wires=0)
+            qml.RZ(r2[0], wires=0)
+            return qml.expval(qml.PauliX(0))
+
+    The circuit before optimization:
+
+    >>> dev = qml.device('default.qubit', wires=1)
+    >>> qnode = qml.QNode(qfunc, dev)
+    >>> print(qml.draw(qnode)([0.1, 0.2, 0.3], [0.4, 0.5, 0.6]))
+    0: ──H──Rot(0.1, 0.2, 0.3)──Rot(0.4, 0.5, 0.6)──RZ(0.1)──RZ(0.4)──┤ ⟨X⟩
+
+    Full single-qubit gate fusion allows us to collapse this entire sequence into a
+    single ``qml.Rot`` rotation gate.
+
+    >>> optimized_qfunc = single_qubit_fusion()(qfunc)
+    >>> optimized_qnode = qml.QNode(optimized_qfunc, dev)
+    >>> print(qml.draw(optimized_qnode)([0.1, 0.2, 0.3], [0.4, 0.5, 0.6]))
+    0: ──Rot(3.57, 2.09, 2.05)──┤ ⟨X⟩
+
+    """
+    # Make a working copy of the list to traverse
+    list_copy = tape.operations.copy()
+
+    while len(list_copy) > 0:
+        current_gate = list_copy[0]
+
+        # If the gate should be excluded, queue it and move on regardless
+        # of fusion potential
+        if exclude_gates is not None:
+            if current_gate.name in exclude_gates:
+                apply(current_gate)
+                list_copy.pop(0)
+                continue
+
+        # Look for single_qubit_rot_angles; if not available, queue and move on.
+        # If available, grab the angles and try to fuse.
+        try:
+            cumulative_angles = stack(current_gate.single_qubit_rot_angles())
+        except (NotImplementedError, AttributeError):
+            apply(current_gate)
+            list_copy.pop(0)
+            continue
+
+        # Find the next gate that acts on the same wires
+        next_gate_idx = find_next_gate(current_gate.wires, list_copy[1:])
+
+        if next_gate_idx is None:
+            apply(current_gate)
+            list_copy.pop(0)
+            continue
+
+        # Before entering the loop, we check to make sure the next gate is not in the
+        # exclusion list. If it is, we should apply the original gate as-is, and not the
+        # Rot version (example in test test_single_qubit_fusion_exclude_gates).
+        if exclude_gates is not None:
+            next_gate = list_copy[next_gate_idx + 1]
+            if next_gate.name in exclude_gates:
+                apply(current_gate)
+                list_copy.pop(0)
+                continue
+
+        # Loop as long as a valid next gate exists
+        while next_gate_idx is not None:
+            next_gate = list_copy[next_gate_idx + 1]
+
+            # Check first if the next gate is in the exclusion list
+            if exclude_gates is not None:
+                if next_gate.name in exclude_gates:
+                    break
+
+            # Try to extract the angles; since the Rot angles are implemented
+            # solely for single-qubit gates, and we used find_next_gate to obtain
+            # the gate in question, only valid single-qubit gates on the same
+            # wire as the current gate will be fused.
+            try:
+                next_gate_angles = next_gate.single_qubit_rot_angles()
+            except (NotImplementedError, AttributeError):
+                break
+
+            cumulative_angles = fuse_rot_angles(
+                cumulative_angles, cast_like(stack(next_gate_angles), cumulative_angles)
+            )
+
+            list_copy.pop(next_gate_idx + 1)
+            next_gate_idx = find_next_gate(current_gate.wires, list_copy[1:])
+
+        # Only apply if the cumulative angle is not close to 0
+        if not allclose(cumulative_angles, zeros(3), atol=atol, rtol=0):
+            Rot(*cumulative_angles, wires=current_gate.wires)
+
+        # Remove the starting gate from the list
+        list_copy.pop(0)
+
+    # Queue the measurements normally
+    for m in tape.measurements:
+        apply(m)


Test patch:
diff --git a/tests/ops/test_qubit_ops.py b/tests/ops/test_qubit_ops.py
index cd2777f6d..4d50669eb 100644
--- a/tests/ops/test_qubit_ops.py
+++ b/tests/ops/test_qubit_ops.py
@@ -739,6 +739,35 @@ class TestOperations:
 
         assert np.allclose(decomposed_matrix, global_phase * op.matrix, atol=tol, rtol=0)
 
+    @pytest.mark.parametrize(
+        "op",
+        [
+            (qml.Hadamard(wires=0)),
+            (qml.PauliX(wires=0)),
+            (qml.PauliY(wires=0)),
+            (qml.PauliZ(wires=0)),
+            (qml.S(wires=0)),
+            (qml.T(wires=0)),
+            (qml.SX(wires=0)),
+            (qml.RX(0.3, wires=0)),
+            (qml.RY(0.3, wires=0)),
+            (qml.RZ(0.3, wires=0)),
+            (qml.PhaseShift(0.3, wires=0)),
+            (qml.Rot(0.3, 0.4, 0.5, wires=0)),
+        ],
+    )
+    def test_single_qubit_rot_angles(self, op):
+        """Tests that the Rot gates yielded by single_qubit_rot_angles
+        are equivalent to the true operations up to a global phase."""
+        angles = op.single_qubit_rot_angles()
+        obtained_mat = qml.Rot(*angles, wires=0).matrix
+
+        # Check whether the two matrices are each others conjugate transposes
+        mat_product = qml.math.dot(op.matrix, qml.math.conj(obtained_mat.T))
+        mat_product /= mat_product[0, 0]
+
+        assert qml.math.allclose(mat_product, I)
+
     def test_CY_decomposition(self, tol):
         """Tests that the decomposition of the CY gate is correct"""
         op = qml.CY(wires=[0, 1])
diff --git a/tests/transforms/test_optimization/test_single_qubit_fusion.py b/tests/transforms/test_optimization/test_single_qubit_fusion.py
new file mode 100644
index 000000000..5a718d31b
--- /dev/null
+++ b/tests/transforms/test_optimization/test_single_qubit_fusion.py
@@ -0,0 +1,292 @@
+# Copyright 2018-2021 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.
+
+import pytest
+from pennylane import numpy as np
+
+import pennylane as qml
+from pennylane.wires import Wires
+from pennylane.transforms.optimization import single_qubit_fusion
+
+from utils import *
+
+
+class TestSingleQubitFusion:
+    """Test that sequences of any single-qubit rotations are fully fused."""
+
+    def test_single_qubit_full_fusion(self):
+        """Test that a sequence of single-qubit gates all fuse."""
+
+        def qfunc():
+            qml.RZ(0.3, wires=0)
+            qml.Hadamard(wires=0)
+            qml.Rot(0.1, 0.2, 0.3, wires=0)
+            qml.RX(0.1, wires=0)
+            qml.SX(wires=0)
+            qml.T(wires=0)
+            qml.PauliX(wires=0)
+
+        transformed_qfunc = single_qubit_fusion()(qfunc)
+
+        original_ops = qml.transforms.make_tape(qfunc)().operations
+        transformed_ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        assert len(transformed_ops) == 1
+
+        # Compare matrices
+        matrix_expected = compute_matrix_from_ops_one_qubit(original_ops)
+        matrix_obtained = compute_matrix_from_ops_one_qubit(transformed_ops)
+        assert check_matrix_equivalence(matrix_expected, matrix_obtained)
+
+    def test_single_qubit_fusion_no_gates_after(self):
+        """Test that gates with nothing after are applied without modification."""
+
+        def qfunc():
+            qml.RZ(0.1, wires=0)
+            qml.Hadamard(wires=1)
+
+        transformed_qfunc = single_qubit_fusion()(qfunc)
+        transformed_ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["RZ", "Hadamard"]
+        wires_expected = [Wires(0), Wires(1)]
+        compare_operation_lists(transformed_ops, names_expected, wires_expected)
+
+    def test_single_qubit_cancelled_fusion(self):
+        """Test if a sequence of single-qubit gates that all cancel yields no operations."""
+
+        def qfunc():
+            qml.RZ(0.1, wires=0)
+            qml.RX(0.2, wires=0)
+            qml.RX(-0.2, wires=0)
+            qml.RZ(-0.1, wires=0)
+
+        transformed_qfunc = single_qubit_fusion()(qfunc)
+        transformed_ops = qml.transforms.make_tape(transformed_qfunc)().operations
+        assert len(transformed_ops) == 0
+
+    def test_single_qubit_fusion_not_implemented(self):
+        """Test that fusion is correctly skipped for single-qubit gates where
+        the rotation angles are not specified."""
+
+        def qfunc():
+            qml.RZ(0.1, wires=0)
+            qml.Hadamard(wires=0)
+            # No rotation angles specified for PauliRot since it is a gate that
+            # in principle acts on an arbitrary number of wires.
+            qml.PauliRot(0.2, "X", wires=0)
+            qml.RZ(0.1, wires=0)
+            qml.Hadamard(wires=0)
+
+        transformed_qfunc = single_qubit_fusion()(qfunc)
+        transformed_ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["Rot", "PauliRot", "Rot"]
+        wires_expected = [Wires(0)] * 3
+        compare_operation_lists(transformed_ops, names_expected, wires_expected)
+
+    def test_single_qubit_fusion_exclude_gates(self):
+        """Test that fusion is correctly skipped for gates explicitly on an
+        exclusion list."""
+
+        def qfunc():
+            # Excluded gate at the start
+            qml.RZ(0.1, wires=0)
+            qml.Hadamard(wires=0)
+            qml.PauliX(wires=0)
+            qml.RZ(0.1, wires=1)
+            qml.CNOT(wires=[0, 1])
+            qml.Hadamard(wires=0)
+            # Excluded gate after another gate
+            qml.RZ(0.1, wires=0)
+            qml.PauliX(wires=1)
+            qml.PauliZ(wires=1)
+            # Excluded gate after multiple others
+            qml.RZ(0.2, wires=1)
+
+        original_ops = qml.transforms.make_tape(qfunc)().operations
+
+        transformed_qfunc = single_qubit_fusion(exclude_gates=["RZ"])(qfunc)
+        transformed_ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["RZ", "Rot", "RZ", "CNOT", "Hadamard", "RZ", "Rot", "RZ"]
+        wires_expected = (
+            [Wires(0)] * 2 + [Wires(1)] + [Wires([0, 1])] + [Wires(0)] * 2 + [Wires(1)] * 2
+        )
+        compare_operation_lists(transformed_ops, names_expected, wires_expected)
+
+        # Compare matrices
+        matrix_expected = compute_matrix_from_ops_two_qubit(original_ops, wire_order=[0, 1])
+        matrix_obtained = compute_matrix_from_ops_two_qubit(transformed_ops, wire_order=[0, 1])
+        assert check_matrix_equivalence(matrix_expected, matrix_obtained)
+
+    def test_single_qubit_fusion_multiple_qubits(self):
+        """Test that all sequences of single-qubit gates across multiple qubits fuse properly."""
+
+        def qfunc():
+            qml.RZ(0.3, wires="a")
+            qml.RY(0.5, wires="a")
+            qml.Rot(0.1, 0.2, 0.3, wires="b")
+            qml.RX(0.1, wires="a")
+            qml.CNOT(wires=["b", "a"])
+            qml.SX(wires="b")
+            qml.S(wires="b")
+            qml.PhaseShift(0.3, wires="b")
+
+        transformed_qfunc = single_qubit_fusion()(qfunc)
+
+        original_ops = qml.transforms.make_tape(qfunc)().operations
+        transformed_ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["Rot", "Rot", "CNOT", "Rot"]
+        wires_expected = [Wires("a"), Wires("b"), Wires(["b", "a"]), Wires("b")]
+        compare_operation_lists(transformed_ops, names_expected, wires_expected)
+
+        # Check matrix representation
+        matrix_expected = compute_matrix_from_ops_two_qubit(original_ops, ["a", "b"])
+        matrix_obtained = compute_matrix_from_ops_two_qubit(transformed_ops, ["a", "b"])
+        assert check_matrix_equivalence(matrix_expected, matrix_obtained)
+
+
+# Example QNode and device for interface testing
+dev = qml.device("default.qubit", wires=3)
+
+# Test each of single-qubit, two-qubit, and Rot gates
+def qfunc(theta):
+    qml.Hadamard(wires=0)
+    qml.RZ(theta[0], wires=0)
+    qml.PauliY(wires=1)
+    qml.RZ(theta[1], wires=0)
+    qml.CNOT(wires=[1, 2])
+    qml.CRY(theta[2], wires=[1, 2])
+    qml.PauliZ(wires=0)
+    qml.CRY(theta[3], wires=[1, 2])
+    qml.Rot(theta[0], theta[1], theta[2], wires=1)
+    qml.Rot(theta[2], theta[3], theta[0], wires=1)
+    return qml.expval(qml.PauliX(0) @ qml.PauliX(2))
+
+
+transformed_qfunc = single_qubit_fusion()(qfunc)
+
+expected_op_list = ["Rot", "Rot", "CNOT", "CRY", "CRY", "Rot"]
+expected_wires_list = [Wires(0), Wires(1), Wires([1, 2]), Wires([1, 2]), Wires([1, 2]), Wires(1)]
+
+
+class TestSingleQubitFusionInterfaces:
+    """Test that rotation merging works in all interfaces."""
+
+    def test_single_qubit_fusion_autograd(self):
+        """Test QNode and gradient in autograd interface."""
+
+        original_qnode = qml.QNode(qfunc, dev)
+        transformed_qnode = qml.QNode(transformed_qfunc, dev)
+
+        input = np.array([0.1, 0.2, 0.3, 0.4], requires_grad=True)
+
+        # Check that the numerical output is the same
+        assert qml.math.allclose(original_qnode(input), transformed_qnode(input))
+
+        # Check that the gradient is the same
+        assert qml.math.allclose(
+            qml.grad(original_qnode)(input), qml.grad(transformed_qnode)(input)
+        )
+
+        # Check operation list
+        ops = transformed_qnode.qtape.operations
+        compare_operation_lists(ops, expected_op_list, expected_wires_list)
+
+    def test_single_qubit_fusion_torch(self):
+        """Test QNode and gradient in torch interface."""
+        torch = pytest.importorskip("torch", minversion="1.8")
+
+        original_qnode = qml.QNode(qfunc, dev, interface="torch")
+        transformed_qnode = qml.QNode(transformed_qfunc, dev, interface="torch")
+
+        original_input = torch.tensor([0.1, 0.2, 0.3, 0.4], requires_grad=True)
+        transformed_input = torch.tensor([0.1, 0.2, 0.3, 0.4], requires_grad=True)
+
+        original_result = original_qnode(original_input)
+        transformed_result = transformed_qnode(transformed_input)
+
+        # Check that the numerical output is the same
+        assert qml.math.allclose(original_result, transformed_result.detach().numpy())
+
+        # Check that the gradient is the same
+        original_result.backward()
+        transformed_result.backward()
+
+        assert qml.math.allclose(original_input.grad, transformed_input.grad)
+
+        # Check operation list
+        ops = transformed_qnode.qtape.operations
+        compare_operation_lists(ops, expected_op_list, expected_wires_list)
+
+    def test_single_qubit_fusion_tf(self):
+        """Test QNode and gradient in tensorflow interface."""
+        tf = pytest.importorskip("tensorflow")
+
+        original_qnode = qml.QNode(qfunc, dev, interface="tf")
+        transformed_qnode = qml.QNode(transformed_qfunc, dev, interface="tf")
+
+        original_input = tf.Variable([0.1, 0.2, 0.3, 0.4])
+        transformed_input = tf.Variable([0.1, 0.2, 0.3, 0.4])
+
+        original_result = original_qnode(original_input)
+        transformed_result = transformed_qnode(transformed_input)
+
+        # Check that the numerical output is the same
+        assert qml.math.allclose(original_result, transformed_result)
+
+        # Check that the gradient is the same
+        with tf.GradientTape() as tape:
+            loss = original_qnode(original_input)
+        original_grad = tape.gradient(loss, original_input)
+
+        with tf.GradientTape() as tape:
+            loss = transformed_qnode(transformed_input)
+        transformed_grad = tape.gradient(loss, transformed_input)
+
+        assert qml.math.allclose(original_grad, transformed_grad)
+
+        # Check operation list
+        ops = transformed_qnode.qtape.operations
+        compare_operation_lists(ops, expected_op_list, expected_wires_list)
+
+    def test_single_qubit_fusion_jax(self):
+        """Test QNode and gradient in JAX interface."""
+        jax = pytest.importorskip("jax")
+        from jax import numpy as jnp
+
+        # Enable float64 support
+        from jax.config import config
+
+        remember = config.read("jax_enable_x64")
+        config.update("jax_enable_x64", True)
+
+        original_qnode = qml.QNode(qfunc, dev, interface="jax")
+        transformed_qnode = qml.QNode(transformed_qfunc, dev, interface="jax")
+
+        input = jnp.array([0.1, 0.2, 0.3, 0.4], dtype=jnp.float64)
+
+        # Check that the numerical output is the same
+        assert qml.math.allclose(original_qnode(input), transformed_qnode(input))
+
+        # Check that the gradient is the same
+        assert qml.math.allclose(
+            jax.grad(original_qnode)(input), jax.grad(transformed_qnode)(input)
+        )
+
+        # Check operation list
+        ops = transformed_qnode.qtape.operations
+        compare_operation_lists(ops, expected_op_list, expected_wires_list)
diff --git a/tests/transforms/test_optimization/utils.py b/tests/transforms/test_optimization/utils.py
index f5995f427..2bf198a05 100644
--- a/tests/transforms/test_optimization/utils.py
+++ b/tests/transforms/test_optimization/utils.py
@@ -28,6 +28,23 @@ def compute_matrix_from_ops_one_qubit(ops):
     return mat
 
 
+def compute_matrix_from_ops_two_qubit(ops, wire_order):
+    """Given a list of two-qubit operations, construct its matrix representation."""
+
+    mat = qml.math.kron(I, I)
+
+    for op in ops:
+        if len(op.wires) == 1:
+            if op.wires == qml.wires.Wires(wire_order[0]):
+                mat = qml.math.dot(qml.math.kron(op.matrix, I), mat)
+            else:
+                mat = qml.math.dot(qml.math.kron(I, op.matrix), mat)
+        else:
+            mat = qml.math.dot(op.matrix, mat)
+
+    return mat
+
+
 def check_matrix_equivalence(matrix_expected, matrix_obtained):
     """Takes two matrices and checks if multiplying one by the conjugate
     transpose of the other gives the identity."""