Repository: PennyLaneAI/pennylane
Pull request #1455 (2021-07-08 17:09:40)

Implement transforms to cancel redundant single-qubit operations
**Context:** Addition of circuit optimization / compilation tools to PennyLane.

**Description of the Change:** Creates a subdirectory `transforms/optimization` with new transforms, `cancel_inverses`, and `merge_rotations`. The `Operator` class, and relevant single-qubit operations, are augmented with additional attributes to enable this.

**Benefits:** Enables reduction of single-qubit gates in a circuit.

**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/1455?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1455](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (76947be) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/aa75f41636b60627397d846ef1c9e72cd780c0e8?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (aa75f41) will **decrease** coverage by `1.12%`.
> The diff coverage is `90.60%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/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/1455?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #1455      +/-   ##
==========================================
- Coverage   99.26%   98.14%   -1.13%     
==========================================
  Files         158      167       +9     
  Lines       11554    12231     +677     
==========================================
+ Hits        11469    12004     +535     
- Misses         85      227     +142     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [...lane/transforms/optimization/optimization\_utils.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvb3B0aW1pemF0aW9uL29wdGltaXphdGlvbl91dGlscy5weQ==) | `75.00% <75.00%> (ø)` | |
| [...nnylane/transforms/optimization/cancel\_inverses.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvb3B0aW1pemF0aW9uL2NhbmNlbF9pbnZlcnNlcy5weQ==) | `96.87% <96.87%> (ø)` | |
| [pennylane/math/single\_dispatch.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL21hdGgvc2luZ2xlX2Rpc3BhdGNoLnB5) | `99.42% <100.00%> (+0.01%)` | :arrow_up: |
| [pennylane/operation.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wZXJhdGlvbi5weQ==) | `95.94% <100.00%> (-3.83%)` | :arrow_down: |
| [pennylane/ops/qubit.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC5weQ==) | `98.68% <100.00%> (-0.16%)` | :arrow_down: |
| [pennylane/transforms/optimization/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/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%> (ø)` | |
| [...nnylane/transforms/optimization/merge\_rotations.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvb3B0aW1pemF0aW9uL21lcmdlX3JvdGF0aW9ucy5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/fourier/coefficients.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ZvdXJpZXIvY29lZmZpY2llbnRzLnB5) | `54.83% <0.00%> (-45.17%)` | :arrow_down: |
| [pennylane/tape/operation\_recorder.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RhcGUvb3BlcmF0aW9uX3JlY29yZGVyLnB5) | `64.51% <0.00%> (-35.49%)` | :arrow_down: |
| [pennylane/measure.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL21lYXN1cmUucHk=) | `91.86% <0.00%> (-8.14%)` | :arrow_down: |
| ... and [36 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455?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/1455?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [aa75f41...76947be](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1455?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).

anthayes92: This looks fantastic @glassnotes ! Some really useful utils functions for compilation too :grin: 

Patch:
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 5407dae5b..c3b3a0696 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -2,6 +2,63 @@
 
 <h3>New features since last release</h3>
 
+* 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)
+
+  The `cancel_inverses` transform loops through a list of operations,
+  and removes adjacent pairs of operations that cancel out. For example,
+
+  ```python
+  def circuit():
+      qml.Hadamard(wires=0)
+      qml.PauliZ(wires=1)
+      qml.Hadamard(wires=0)
+      qml.T(wires=0)
+      qml.CZ(wires=[0, 1])
+      qml.CZ(wires=[1, 0])
+      return qml.expval(qml.PauliX(wires=0))
+  ```
+
+  ```pycon
+  >>> dev = qml.device('default.qubit', wires=2)
+  >>> qnode = qml.QNode(circuit, dev)
+  >>> print(qml.draw(qnode)())
+   0: ──H──H──T──╭C──╭Z──┤ ⟨X⟩
+   1: ──Z────────╰Z──╰C──┤
+  >>> optimized_circuit = qml.transforms.cancel_inverses(circuit)
+  >>> optimized_qnode = qml.QNode(optimized_circuit, dev)
+  >>> print(qml.draw(optimized_qnode)())
+   0: ──T──┤ ⟨X⟩
+   1: ──Z──┤
+  ```
+
+  The `merge_rotations` transform combines adjacent rotation gates of
+  the same type into a single gate, including controlled rotations.
+
+  ```python
+  def circuit(x, y, z):
+      qml.RX(x, wires=0)
+      qml.RX(x, wires=0)
+      qml.Rot(x, y, z, wires=1)
+      qml.Rot(y, z, x, wires=1)
+      qml.CRY(y, wires=[0, 1])
+      qml.CRY(y + z, wires=[0, 1])
+      return qml.expval(qml.PauliX(wires=0))
+  ```
+
+  ```pycon
+  >>> qnode = qml.QNode(circuit, dev)
+  >>> print(qml.draw(qnode)(0.1, 0.2, 0.3))
+   0: ──RX(0.1)─────────────RX(0.1)─────────────╭C────────╭C────────┤ ⟨X⟩
+   1: ──Rot(0.1, 0.2, 0.3)──Rot(0.2, 0.3, 0.1)──╰RY(0.2)──╰RY(0.5)──┤
+  >>> optimized_circuit = qml.transforms.merge_rotations()(circuit)
+  >>> optimized_qnode = qml.QNode(optimized_circuit, dev)
+  >>> print(qml.draw(optimized_qnode)(0.1, 0.2, 0.3))
+   0: ──RX(0.2)───────────────────╭C────────┤ ⟨X⟩
+   1: ──Rot(0.409, 0.485, 0.306)──╰RY(0.7)──┤
+  ```
+
 * A decomposition has been added to ``QubitUnitary`` that makes the
   single-qubit case fully differentiable in all interfaces. Furthermore,
   a quantum function transform, ``unitary_to_rot()``, has been added to decompose all
diff --git a/pennylane/math/single_dispatch.py b/pennylane/math/single_dispatch.py
index cbf56110d..312daa9c2 100644
--- a/pennylane/math/single_dispatch.py
+++ b/pennylane/math/single_dispatch.py
@@ -117,12 +117,14 @@ ar.autoray._SUBMODULE_ALIASES["tensorflow", "angle"] = "tensorflow.math"
 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", "diag"] = "tensorflow.linalg"
 
 
 ar.autoray._FUNC_ALIASES["tensorflow", "arcsin"] = "asin"
 ar.autoray._FUNC_ALIASES["tensorflow", "arccos"] = "acos"
 ar.autoray._FUNC_ALIASES["tensorflow", "arctan"] = "atan"
+ar.autoray._FUNC_ALIASES["tensorflow", "arctan2"] = "atan2"
 ar.autoray._FUNC_ALIASES["tensorflow", "diag"] = "diag"
 
 
@@ -243,6 +245,9 @@ ar.register_function(
     ),
 )
 
+ar.autoray._SUBMODULE_ALIASES["torch", "arctan2"] = "torch"
+ar.autoray._FUNC_ALIASES["torch", "arctan2"] = "atan2"
+
 
 def _take_torch(tensor, indices, axis=None):
     """Torch implementation of np.take"""
diff --git a/pennylane/operation.py b/pennylane/operation.py
index 45d947249..e246ca02c 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -539,6 +539,49 @@ class Operation(Operator):
         s_1]=[-1/2, 1, -\pi/2]` is assumed for every parameter.
     """
 
+    # Attributes for compilation transforms
+    is_self_inverse = None
+    """bool or None: ``True`` if the operation is its own inverse.
+
+    If ``None``, all instances of the given operation will be ignored during
+    compilation transforms involving inverse cancellation.
+    """
+
+    is_symmetric_over_all_wires = None
+    """bool or None: ``True`` if the operation is the same if you exchange the order
+    of wires.
+
+    For example, ``qml.CZ(wires=[0, 1])`` has the same effect as ``qml.CZ(wires=[1,
+    0])`` due to symmetry of the operation.
+
+    If ``None``, all instances of the operation will be ignored during
+    compilation transforms that check for wire symmetry.
+    """
+
+    is_symmetric_over_control_wires = None
+    """bool or None: ``True`` if the operation is the same if you exchange the order
+    of all but the last wire.
+
+    For example, ``qml.Toffoli(wires=[0, 1, 2])`` has the same effect as
+    ``qml.Toffoli(wires=[1, 0, 2])``, but neither are the same as
+    ``qml.Toffoli(wires=[0, 2, 1])``.
+
+    If ``None``, all instances of the operation will be ignored during
+    compilation transforms that check for control-wire symmetry.
+    """
+
+    is_composable_rotation = None
+    """bool or None: ``True`` if composing multiple copies of the operation
+    results in an addition (or alternative accumulation) of parameters.
+
+    For example, ``qml.RZ`` is a composable rotation. Applying ``qml.RZ(0.1,
+    wires=0)`` followed by ``qml.RZ(0.2, wires=0)`` is equivalent to performing
+    a single rotation ``qml.RZ(0.3, wires=0)``.
+
+    If set to ``None``, the operation will be ignored during compilation
+    transforms that merge adjacent rotations.
+    """
+
     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 b2442c12d..22ae645cf 100644
--- a/pennylane/ops/qubit.py
+++ b/pennylane/ops/qubit.py
@@ -58,6 +58,7 @@ class Hadamard(Observable, Operation):
     num_params = 0
     num_wires = 1
     par_domain = None
+    is_self_inverse = True
     eigvals = pauli_eigs(1)
     matrix = np.array([[INV_SQRT2, INV_SQRT2], [INV_SQRT2, -INV_SQRT2]])
 
@@ -115,6 +116,7 @@ class PauliX(Observable, Operation):
     num_params = 0
     num_wires = 1
     par_domain = None
+    is_self_inverse = True
     eigvals = pauli_eigs(1)
     matrix = np.array([[0, 1], [1, 0]])
 
@@ -173,6 +175,7 @@ class PauliY(Observable, Operation):
     num_params = 0
     num_wires = 1
     par_domain = None
+    is_self_inverse = True
     eigvals = pauli_eigs(1)
     matrix = np.array([[0, -1j], [1j, 0]])
 
@@ -233,6 +236,7 @@ class PauliZ(Observable, DiagonalOperation):
     num_params = 0
     num_wires = 1
     par_domain = None
+    is_self_inverse = True
     eigvals = pauli_eigs(1)
     matrix = np.array([[1, 0], [0, -1]])
 
@@ -402,6 +406,7 @@ class CNOT(Operation):
     num_params = 0
     num_wires = 2
     par_domain = None
+    is_self_inverse = True
     matrix = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
 
     @classmethod
@@ -439,6 +444,8 @@ class CZ(DiagonalOperation):
     num_params = 0
     num_wires = 2
     par_domain = None
+    is_self_inverse = True
+    is_symmetric_over_all_wires = True
     eigvals = np.array([1, 1, 1, -1])
     matrix = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]])
 
@@ -478,6 +485,7 @@ class CY(Operation):
     num_params = 0
     num_wires = 2
     par_domain = None
+    is_self_inverse = True
     matrix = np.array(
         [
             [1, 0, 0, 0],
@@ -522,6 +530,8 @@ class SWAP(Operation):
     num_params = 0
     num_wires = 2
     par_domain = None
+    is_self_inverse = True
+    is_symmetric_over_all_wires = True
     matrix = np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]])
 
     @classmethod
@@ -678,6 +688,8 @@ class Toffoli(Operation):
     num_params = 0
     num_wires = 3
     par_domain = None
+    is_self_inverse = True
+    is_symmetric_over_control_wires = True
     matrix = np.array(
         [
             [1, 0, 0, 0, 0, 0, 0, 0],
@@ -743,6 +755,7 @@ class RX(Operation):
     num_params = 1
     num_wires = 1
     par_domain = "R"
+    is_composable_rotation = True
     grad_method = "A"
     generator = [PauliX, -1 / 2]
 
@@ -784,6 +797,7 @@ class RY(Operation):
     num_params = 1
     num_wires = 1
     par_domain = "R"
+    is_composable_rotation = True
     grad_method = "A"
     generator = [PauliY, -1 / 2]
 
@@ -825,6 +839,7 @@ class RZ(DiagonalOperation):
     num_params = 1
     num_wires = 1
     par_domain = "R"
+    is_composable_rotation = True
     grad_method = "A"
     generator = [PauliZ, -1 / 2]
 
@@ -872,6 +887,7 @@ class PhaseShift(DiagonalOperation):
     num_params = 1
     num_wires = 1
     par_domain = "R"
+    is_composable_rotation = True
     grad_method = "A"
     generator = [np.array([[0, 0], [0, 1]]), 1]
 
@@ -924,6 +940,7 @@ class ControlledPhaseShift(DiagonalOperation):
     num_params = 1
     num_wires = 2
     par_domain = "R"
+    is_composable_rotation = True
     grad_method = "A"
     generator = [np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1]]), 1]
 
@@ -988,6 +1005,7 @@ class Rot(Operation):
     num_params = 3
     num_wires = 1
     par_domain = "R"
+    is_composable_rotation = True
     grad_method = "A"
 
     @classmethod
@@ -1341,6 +1359,7 @@ class CRX(Operation):
     num_params = 1
     num_wires = 2
     par_domain = "R"
+    is_composable_rotation = True
     grad_method = "A"
     grad_recipe = four_term_grad_recipe
 
@@ -1410,6 +1429,7 @@ class CRY(Operation):
     num_params = 1
     num_wires = 2
     par_domain = "R"
+    is_composable_rotation = True
     grad_method = "A"
     grad_recipe = four_term_grad_recipe
 
@@ -1480,6 +1500,7 @@ class CRZ(DiagonalOperation):
     num_params = 1
     num_wires = 2
     par_domain = "R"
+    is_composable_rotation = True
     grad_method = "A"
     grad_recipe = four_term_grad_recipe
 
diff --git a/pennylane/transforms/__init__.py b/pennylane/transforms/__init__.py
index 9a18760fb..d643425d5 100644
--- a/pennylane/transforms/__init__.py
+++ b/pennylane/transforms/__init__.py
@@ -44,7 +44,9 @@ containing quantum operations) that are used to construct QNodes.
     :toctree: api
 
     ~adjoint
+    ~transforms.cancel_inverses
     ~ctrl
+    ~transforms.merge_rotations
     ~transforms.unitary_to_rot
     ~transforms.invisible
     ~apply_controlled_Q
@@ -95,6 +97,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 .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
new file mode 100644
index 000000000..adaff1307
--- /dev/null
+++ b/pennylane/transforms/optimization/__init__.py
@@ -0,0 +1,19 @@
+# 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.
+"""
+This subpackage contains quantum function transforms for optimizing quantum circuits.
+"""
+
+from .cancel_inverses import cancel_inverses
+from .merge_rotations import merge_rotations
diff --git a/pennylane/transforms/optimization/cancel_inverses.py b/pennylane/transforms/optimization/cancel_inverses.py
new file mode 100644
index 000000000..8e1536dd4
--- /dev/null
+++ b/pennylane/transforms/optimization/cancel_inverses.py
@@ -0,0 +1,144 @@
+# 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 cancelling adjacent inverse gates in quantum circuits."""
+# pylint: disable=too-many-branches
+from pennylane import apply
+from pennylane.wires import Wires
+from pennylane.transforms import qfunc_transform
+
+from .optimization_utils import find_next_gate
+
+
+@qfunc_transform
+def cancel_inverses(tape):
+    """Quantum function transform to remove any operations that are applied next to their
+    (self-)inverse.
+
+    Args:
+        qfunc (function): A quantum function.
+
+    **Example**
+
+    Consider the following quantum function:
+
+    .. code-block:: python
+
+        def qfunc(x, y, z):
+            qml.Hadamard(wires=0)
+            qml.Hadamard(wires=1)
+            qml.Hadamard(wires=0)
+            qml.RX(x, wires=2)
+            qml.RY(y, wires=1)
+            qml.PauliX(wires=1)
+            qml.RZ(z, wires=0)
+            qml.RX(y, wires=2)
+            qml.CNOT(wires=[0, 2])
+            qml.PauliX(wires=1)
+            return qml.expval(qml.PauliZ(0))
+
+    The circuit before optimization:
+
+    >>> dev = qml.device('default.qubit', wires=3)
+    >>> qnode = qml.QNode(qfunc, dev)
+    >>> print(qml.draw(qnode)(1, 2, 3))
+    0: ──H──────H──────RZ(3)─────╭C──┤ ⟨Z⟩
+    1: ──H──────RY(2)──X──────X──│───┤
+    2: ──RX(1)──RX(2)────────────╰X──┤
+
+    We can see that there are two adjacent Hadamards on the first qubit that
+    should cancel each other out. Similarly, there are two Pauli-X gates on the
+    second qubit that should cancel. We can obtain a simplified circuit by running
+    the ``cancel_inverses`` transform:
+
+    >>> optimized_qfunc = cancel_inverses(qfunc)
+    >>> optimized_qnode = qml.QNode(optimized_qfunc, dev)
+    >>> print(qml.draw(optimized_qnode)(1, 2, 3))
+    0: ──RZ(3)─────────╭C──┤ ⟨Z⟩
+    1: ──H──────RY(2)──│───┤
+    2: ──RX(1)──RX(2)──╰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]
+
+        # Find the next gate that acts on at least one of the same wires
+        next_gate_idx = find_next_gate(current_gate.wires, list_copy[1:])
+
+        # If no such gate is found queue the operation and move on
+        if next_gate_idx is None:
+            apply(current_gate)
+            list_copy.pop(0)
+            continue
+
+        # Otherwise, get the next gate
+        next_gate = list_copy[next_gate_idx + 1]
+
+        # There are then three possibilities that may lead to inverse cancellation. For a gate U,
+        # 1. U is self-inverse, and the next gate is also U
+        # 2. The current gate is U.inv and the next gate is U
+        # 3. The current gate is U and the next gate is U.inv
+
+        # Case 1
+        are_self_inverses = current_gate.is_self_inverse and current_gate.name == next_gate.name
+
+        # Cases 2 and 3
+        are_inverses = False
+        name_set = set([current_gate.name, next_gate.name])
+        shortest_name = min(name_set)
+        if set([shortest_name, shortest_name + ".inv"]) == name_set:
+            are_inverses = True
+
+        # If either of the two flags is true, we can potentially cancel the gates
+        if are_self_inverses or are_inverses:
+            # If the wires are the same, then we can safely remove both
+            if current_gate.wires == next_gate.wires:
+                list_copy.pop(next_gate_idx + 1)
+            # If wires are not equal, there are two things that can happen
+            else:
+                # There is not full overlap in the wires; we cannot cancel
+                if len(Wires.shared_wires([current_gate.wires, next_gate.wires])) != len(
+                    current_gate.wires
+                ):
+                    apply(current_gate)
+                # There is full overlap, but the wires are in a different order
+                else:
+                    # If the wires are in a different order, gates that are "symmetric"
+                    # over all wires (e.g., CZ), can be cancelled.
+                    if current_gate.is_symmetric_over_all_wires:
+                        list_copy.pop(next_gate_idx + 1)
+                    # For other gates, as long as the control wires are the same, we can still
+                    # cancel (e.g., the Toffoli gate).
+                    elif current_gate.is_symmetric_over_control_wires:
+                        if (
+                            len(Wires.shared_wires([current_gate.wires[:-1], next_gate.wires[:-1]]))
+                            == len(current_gate.wires) - 1
+                        ):
+                            list_copy.pop(next_gate_idx + 1)
+                        else:
+                            apply(current_gate)
+                    # Apply gate any cases where there is no wire symmetry
+                    else:
+                        apply(current_gate)
+        # If neither of the flags are true, queue and move on to the next item
+        else:
+            apply(current_gate)
+
+        # Remove this gate from the working list
+        list_copy.pop(0)
+
+    # Queue the measurements normally
+    for m in tape.measurements:
+        apply(m)
diff --git a/pennylane/transforms/optimization/merge_rotations.py b/pennylane/transforms/optimization/merge_rotations.py
new file mode 100644
index 000000000..f78a9d5ef
--- /dev/null
+++ b/pennylane/transforms/optimization/merge_rotations.py
@@ -0,0 +1,153 @@
+# 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 merging adjacent rotations of the same type in a quantum circuit."""
+
+from pennylane import apply
+from pennylane.transforms import qfunc_transform
+from pennylane.math import allclose, stack, cast_like, zeros
+
+from .optimization_utils import find_next_gate, fuse_rot_angles
+
+
+@qfunc_transform
+def merge_rotations(tape, atol=1e-8, include_gates=None):
+    r"""Quantum function transform to combine rotation gates of the same type
+    that act sequentially.
+
+    If the combination of two rotation produces an angle that is close to 0,
+    neither gate will be applied.
+
+    Args:
+        qfunc (function): A quantum function.
+        atol (float): After fusion of gates, if the fused angle :math:`\theta` is such that
+            :math:`|\theta|\leq \text{atol}`, no rotation gate will be applied.
+        include_gates (None or list[str]): A list of specific operations to merge. If
+            set to ``None`` (default), all operations with the ``is_composable_rotation``
+            attribute set to ``True`` will be merged. Otherwise, only the operations whose
+            names match those in the list will undergo merging.
+
+    **Example**
+
+    Consider the following quantum function.
+
+    .. code-block:: python
+
+        def qfunc(x, y, z):
+            qml.RX(x, wires=0)
+            qml.RX(y, wires=0)
+            qml.CNOT(wires=[1, 2])
+            qml.RY(y, wires=1)
+            qml.Hadamard(wires=2)
+            qml.CRZ(z, wires=[2, 0])
+            qml.RY(-y, wires=1)
+            return qml.expval(qml.PauliZ(0))
+
+    The circuit before optimization:
+
+    >>> dev = qml.device('default.qubit', wires=3)
+    >>> qnode = qml.QNode(qfunc, dev)
+    >>> print(qml.draw(qnode)(1, 2, 3))
+     0: ───RX(1)──RX(2)──────────╭RZ(3)──┤ ⟨Z⟩
+     1: ──╭C──────RY(2)──RY(-2)──│───────┤
+     2: ──╰X──────H──────────────╰C──────┤
+
+    By inspection, we can combine the two ``RX`` rotations on the first qubit.
+    On the second qubit, we have a cumulative angle of 0, and the gates will cancel.
+
+    >>> optimized_qfunc = merge_rotations()(qfunc)
+    >>> optimized_qnode = qml.QNode(optimized_qfunc, dev)
+    >>> print(qml.draw(optimized_qnode)(1, 2, 3))
+     0: ───RX(3)─────╭RZ(3)──┤ ⟨Z⟩
+     1: ──╭C─────────│───────┤
+     2: ──╰X──────H──╰C──────┤
+
+    It is also possible to explicitly specify which rotations ``merge_rotations`` should
+    be merged using the ``include_gates`` argument. For example, if in the above
+    circuit we wanted only to merge the "RX" gates, we could do so as follows:
+
+    >>> optimized_qfunc = merge_rotations(include_gates=["RX"])(qfunc)
+    >>> optimized_qnode = qml.QNode(optimized_qfunc, dev)
+    >>> print(qml.draw(optimized_qnode)(1, 2, 3))
+     0: ───RX(3)─────────────────╭RZ(3)──┤ ⟨Z⟩
+     1: ──╭C──────RY(2)──RY(-2)──│───────┤
+     2: ──╰X──────H──────────────╰C──────┤
+    """
+    # 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 a specific list of operations is specified, check and see if our
+        # op is in it, then try to merge. If not, queue and move on.
+        if include_gates is not None:
+            if current_gate.name not in include_gates:
+                apply(current_gate)
+                list_copy.pop(0)
+                continue
+
+        # Check if the rotation is composable; if it is not, move on.
+        if not current_gate.is_composable_rotation:
+            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 no such gate is found (either there simply is none, or there are other gates
+        # "in the way", queue the operation and move on
+        if next_gate_idx is None:
+            apply(current_gate)
+            list_copy.pop(0)
+            continue
+
+        # We need to use stack to get this to work and be differentiable in all interfaces
+        cumulative_angles = stack(current_gate.parameters)
+
+        # 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
+            next_gate = list_copy[next_gate_idx + 1]
+
+            # If next gate is of the same type, we can merge the angles
+            if current_gate.name == next_gate.name and current_gate.wires == next_gate.wires:
+                list_copy.pop(next_gate_idx + 1)
+
+                # The Rot gate must be treated separately
+                if current_gate.name == "Rot":
+                    cumulative_angles = fuse_rot_angles(
+                        cumulative_angles, cast_like(stack(next_gate.parameters), cumulative_angles)
+                    )
+                # Other, single-parameter rotation gates just have the angle summed
+                else:
+                    cumulative_angles = cumulative_angles + cast_like(
+                        stack(next_gate.parameters), cumulative_angles
+                    )
+            # If it is not, we need to stop
+            else:
+                break
+
+            # If we did merge, look now at the next gate
+            next_gate_idx = find_next_gate(current_gate.wires, list_copy[1:])
+
+        if not allclose(cumulative_angles, zeros(len(cumulative_angles)), atol=atol, rtol=0):
+            current_gate.__class__(*cumulative_angles, wires=current_gate.wires)
+
+        # Remove the first gate gate from the working list
+        list_copy.pop(0)
+
+    # Queue the measurements normally
+    for m in tape.measurements:
+        apply(m)
diff --git a/pennylane/transforms/optimization/optimization_utils.py b/pennylane/transforms/optimization/optimization_utils.py
new file mode 100644
index 000000000..e62a133c6
--- /dev/null
+++ b/pennylane/transforms/optimization/optimization_utils.py
@@ -0,0 +1,158 @@
+# 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.
+"""Utility functions for circuit optimization."""
+# pylint: disable=too-many-return-statements
+from pennylane.math import allclose, sin, cos, arccos, arctan2, zeros, cast_like, stack
+from pennylane.wires import Wires
+
+
+def find_next_gate(wires, op_list):
+    """Given a list of operations, finds the next operation that acts on at least one of
+    the same set of wires, if present.
+
+    Args:
+        wires (Wires): A set of wires acted on by a quantum operation.
+        op_list (list[Operation]): A list of operations that are implemented after the
+            operation that acts on ``wires``.
+
+    Returns:
+        int or None: The index, in ``op_list``, of the earliest gate that uses one or more
+        of the same wires, or ``None`` if no such gate is present.
+    """
+    next_gate_idx = None
+
+    for op_idx, op in enumerate(op_list):
+        if len(Wires.shared_wires([wires, op.wires])) > 0:
+            next_gate_idx = op_idx
+            break
+
+    return next_gate_idx
+
+
+def _yzy_to_zyz(middle_yzy):
+    """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)
+    can be expressed as a sequence of rotations about 3 axes in 12 different ways.
+    Typically, the arbitrary single-qubit rotation is expressed as RZ(a) RY(b) RZ(c),
+    but there are some situations, e.g., composing two such rotations, where we need
+    to convert between representations. This function converts the angles of a sequence
+
+    .. math::
+
+       RY(y_1) RZ(z) RY(y_2)
+
+    into the form
+
+    .. math::
+
+       RZ(z_1) RY(y) RZ(z_2)
+
+    This is accomplished by first converting the rotation to quaternion form, and then
+    extracting the desired set of angles.
+
+    Args:
+        y1 (float): The angle of the first ``RY`` rotation.
+        z (float): The angle of the inner ``RZ`` rotation.
+        y2 (float): The angle of the second ``RY`` rotation.
+
+    Returns:
+        tuple[float, float, float]: A list of rotation angles in the ZYZ representation.
+    """
+    if allclose(stack(middle_yzy), cast_like(zeros(3), stack(middle_yzy))):
+        return stack([0.0, 0.0, 0.0])
+
+    y1, z, y2 = middle_yzy[0], middle_yzy[1], middle_yzy[2]
+
+    # First, compute the quaternion representation
+    # https://ntrs.nasa.gov/api/citations/19770024290/downloads/19770024290.pdf
+    qw = cos(z / 2) * cos(0.5 * (y1 + y2))
+    qx = sin(z / 2) * sin(0.5 * (y1 - y2))
+    qy = cos(z / 2) * sin(0.5 * (y1 + y2))
+    qz = sin(z / 2) * cos(0.5 * (y1 - y2))
+
+    # Now convert from YZY Euler angles to ZYZ angles
+    # Source: http://bediyap.com/programming/convert-quaternion-to-euler-rotations/
+    z1_arg1 = 2 * (qy * qz - qw * qx)
+    z1_arg2 = 2 * (qx * qz + qw * qy)
+    z1 = arctan2(z1_arg1, z1_arg2)
+
+    y = arccos(qw ** 2 - qx ** 2 - qy ** 2 + qz ** 2)
+
+    z2_arg1 = 2 * (qy * qz + qw * qx)
+    z2_arg2 = -2 * (qx * qz - qw * qy)
+    z2 = arctan2(z2_arg1, z2_arg2)
+
+    return stack([z1, y, z2])
+
+
+def fuse_rot_angles(angles_1, angles_2):
+    """Computed the set of rotation angles that is obtained when composing
+    two ``qml.Rot`` operations.
+
+    The ``qml.Rot`` operation represents the most general single-qubit operation.
+    Two such operations can be fused into a new operation, however the angular dependence
+    is non-trivial.
+
+    Args:
+        angles_1 (float): A set of three angles for the first ``qml.Rot`` operation.
+        angles_2 (float): A set of three angles for the second ``qml.Rot`` operation.
+
+    Returns:
+        array[float]: A tuple of rotation angles for a single ``qml.Rot`` operation
+        that implements the same operation as the two sets of input angles.
+    """
+    # Check for all-zero instances; if there are some, we can just return the sum.
+    are_angles_1_zero = allclose(angles_1, zeros(3))
+    are_angles_2_zero = allclose(angles_2, zeros(3))
+
+    if are_angles_1_zero or are_angles_2_zero:
+        return stack([angles_1[i] + angles_2[i] for i in range(3)])
+
+    # RZ(a) RY(b) RZ(c) fused with RZ(d) RY(e) RZ(f)
+    # first produces RZ(a) RY(b) RZ(c+d) RY(e) RZ(f)
+    left_z = angles_1[0]
+    middle_yzy = angles_1[1], angles_1[2] + angles_2[0], angles_2[1]
+    right_z = angles_2[2]
+
+    # There are a few other cases to consider where things can be 0 and
+    # avoid having to use the quaternion conversion routine
+    # If b = 0, then we have RZ(a + c + d) RY(e) RZ(f)
+    if allclose(middle_yzy[0], 0.0):
+        # Then if e is close to zero, return a single rotation RZ(a + c + d + f)
+        if allclose(middle_yzy[2], 0.0):
+            return [left_z + middle_yzy[1] + right_z, 0.0, 0.0]
+        return stack([left_z + middle_yzy[1], middle_yzy[2], right_z])
+
+    # If c + d is close to 0, then we have the case RZ(a) RY(b + e) RZ(f)
+    if allclose(middle_yzy[1], 0.0):
+        # If b + e is 0, we have RZ(a + f)
+        if allclose(middle_yzy[0] + middle_yzy[2], 0.0):
+            return [left_z + right_z, 0.0, 0.0]
+        return stack([left_z, middle_yzy[0] + middle_yzy[2], right_z])
+
+    # If e is close to 0, then we have the case RZ(a) RY(b) RZ(c + d + f)
+    # The case where b is 0 actually already covered in the first loop,
+    # so only one case here.
+    if allclose(middle_yzy[2], 0.0):
+        return stack([left_z, middle_yzy[0], middle_yzy[1] + right_z])
+
+    # Otherwise, we need to turn the RY(b) RZ(c+d) RY(e) into something
+    # of the form RZ(u) RY(v) RZ(w)
+    u, v, w = _yzy_to_zyz(middle_yzy)
+
+    # Then we can combine to create
+    # RZ(a + u) RY(v) RZ(w + f)
+    return stack([left_z + u, v, w + right_z])


Test patch:
diff --git a/tests/transforms/test_optimization/test_cancel_inverses.py b/tests/transforms/test_optimization/test_cancel_inverses.py
new file mode 100644
index 000000000..16e597cc4
--- /dev/null
+++ b/tests/transforms/test_optimization/test_cancel_inverses.py
@@ -0,0 +1,312 @@
+# 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 cancel_inverses
+from utils import compare_operation_lists
+
+
+class TestCancelInverses:
+    """Test that adjacent inverse gates are cancelled."""
+
+    def test_one_qubit_cancel_adjacent_self_inverse(self):
+        """Test that a single-qubit circuit with adjacent self-inverse gate cancels."""
+
+        def qfunc():
+            qml.Hadamard(wires=0)
+            qml.Hadamard(wires=0)
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        new_tape = qml.transforms.make_tape(transformed_qfunc)()
+
+        assert len(new_tape.operations) == 0
+
+    def test_one_qubit_cancel_adjacent_inverse(self):
+        """Test that a single-qubit circuit with adjacent inverse gate cancels."""
+
+        def qfunc():
+            qml.S(wires=0)
+            qml.adjoint(qml.S)(wires=0)
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        new_tape = qml.transforms.make_tape(transformed_qfunc)()
+
+        assert len(new_tape.operations) == 0
+
+    def test_one_qubit_no_inverse(self):
+        """Test that a one-qubit circuit with a gate in the way does not cancel the inverses."""
+
+        def qfunc():
+            qml.Hadamard(wires=0)
+            qml.RZ(0.3, wires=0)
+            qml.Hadamard(wires=0)
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["Hadamard", "RZ", "Hadamard"]
+        wires_expected = [Wires(0)] * 3
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+    def test_two_qubits_no_inverse(self):
+        """Test that a two-qubit circuit self-inverse on each qubit does not cancel."""
+
+        def qfunc():
+            qml.Hadamard(wires=0)
+            qml.Hadamard(wires=1)
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["Hadamard"] * 2
+        wires_expected = [Wires(0), Wires(1)]
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+    def test_three_qubits_inverse_after_cnot(self):
+        """Test that a three-qubit circuit with a CNOT still allows cancellation."""
+
+        def qfunc():
+            qml.Hadamard(wires=0)
+            qml.PauliX(wires=1)
+            qml.CNOT(wires=[0, 2])
+            qml.RZ(0.5, wires=2)
+            qml.PauliX(wires=1)
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["Hadamard", "CNOT", "RZ"]
+        wires_expected = [Wires(0), Wires([0, 2]), Wires(2)]
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+    def test_three_qubits_blocking_cnot(self):
+        """Test that a three-qubit circuit with a blocking CNOT causes no cancellation."""
+
+        def qfunc():
+            qml.Hadamard(wires=0)
+            qml.PauliX(wires=1)
+            qml.CNOT(wires=[0, 1])
+            qml.RZ(0.5, wires=2)
+            qml.PauliX(wires=1)
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["Hadamard", "PauliX", "CNOT", "RZ", "PauliX"]
+        wires_expected = [Wires(0), Wires(1), Wires([0, 1]), Wires(2), Wires(1)]
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+    def test_three_qubits_toffolis(self):
+        """Test that Toffolis on different permutations of wires cancel correctly."""
+
+        def qfunc():
+            # These two will cancel
+            qml.Toffoli(wires=["a", "b", "c"])
+            qml.Toffoli(wires=["b", "a", "c"])
+            # These three will not cancel
+            qml.Toffoli(wires=["a", "b", "c"])
+            qml.Toffoli(wires=["a", "c", "b"])
+            qml.Toffoli(wires=["a", "c", "d"])
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["Toffoli"] * 3
+        wires_expected = [Wires(["a", "b", "c"]), Wires(["a", "c", "b"]), Wires(["a", "c", "d"])]
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+    def test_two_qubits_cnot_same_direction(self):
+        """Test that two adjacent CNOTs cancel."""
+
+        def qfunc():
+            qml.CNOT(wires=[0, 1])
+            qml.CNOT(wires=[0, 1])
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        assert len(ops) == 0
+
+    def test_two_qubits_cnot_opposite_direction(self):
+        """Test that two adjacent CNOTs with the control/target flipped do NOT cancel."""
+
+        def qfunc():
+            qml.CNOT(wires=[0, 1])
+            qml.CNOT(wires=[1, 0])
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["CNOT"] * 2
+        wires_expected = [Wires([0, 1]), Wires([1, 0])]
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+    def test_two_qubits_cz_opposite_direction(self):
+        """Test that two adjacent CZ with the control/target flipped do cancel due to symmetry."""
+
+        def qfunc():
+            qml.CZ(wires=[0, 1])
+            qml.CZ(wires=[1, 0])
+
+        transformed_qfunc = cancel_inverses(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        assert len(ops) == 0
+
+
+# Example QNode and device for interface testing
+dev = qml.device("default.qubit", wires=3)
+
+
+def qfunc(theta):
+    qml.Hadamard(wires=0)
+    qml.PauliX(wires=1)
+    qml.S(wires=1)
+    qml.adjoint(qml.S)(wires=1)
+    qml.Hadamard(wires=0)
+    qml.CNOT(wires=[0, 1])
+    qml.RZ(theta[0], wires=2)
+    qml.PauliX(wires=1)
+    qml.CZ(wires=[1, 0])
+    qml.RY(theta[1], wires=2)
+    qml.CZ(wires=[0, 1])
+    return qml.expval(qml.PauliX(0) @ qml.PauliX(2))
+
+
+transformed_qfunc = cancel_inverses(qfunc)
+
+expected_op_list = ["PauliX", "CNOT", "RZ", "PauliX", "RY"]
+expected_wires_list = [Wires(1), Wires([0, 1]), Wires(2), Wires(1), Wires(2)]
+
+
+class TestCancelInversesInterfaces:
+    """Test that adjacent inverse gates are cancelled in all interfaces."""
+
+    def test_cancel_inverses_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], 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_cancel_inverses_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], requires_grad=True)
+        transformed_input = torch.tensor([0.1, 0.2], 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)
+
+        # 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_cancel_inverses_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])
+        transformed_input = tf.Variable([0.1, 0.2])
+
+        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_cancel_inverses_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], 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/test_merge_rotations.py b/tests/transforms/test_optimization/test_merge_rotations.py
new file mode 100644
index 000000000..8020bb643
--- /dev/null
+++ b/tests/transforms/test_optimization/test_merge_rotations.py
@@ -0,0 +1,400 @@
+# 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 merge_rotations
+from utils import compare_operation_lists
+
+
+class TestMergeRotations:
+    """Test that adjacent rotation gates of the same type will add the angles."""
+
+    @pytest.mark.parametrize(
+        ("theta_1", "theta_2", "expected_ops"),
+        [
+            (0.3, -0.2, [qml.RZ(0.1, wires=0)]),
+            (0.15, -0.15, []),
+        ],
+    )
+    def test_one_qubit_rotation_merge(self, theta_1, theta_2, expected_ops):
+        """Test that a single-qubit circuit with adjacent rotation along the same
+        axis either merge, or cancel if the angles sum to 0."""
+
+        def qfunc():
+            qml.RZ(theta_1, wires=0)
+            qml.RZ(theta_2, wires=0)
+
+        transformed_qfunc = merge_rotations()(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        assert len(ops) == len(expected_ops)
+
+        # Check that all operations and parameter values are as expected
+        for op_obtained, op_expected in zip(ops, expected_ops):
+            assert op_obtained.name == op_expected.name
+            assert np.allclose(op_obtained.parameters, op_expected.parameters)
+
+    @pytest.mark.parametrize(
+        ("theta_1", "theta_2", "expected_ops"),
+        [
+            (0.3, -0.2, [qml.RZ(0.3, wires=0), qml.RZ(-0.2, wires=1)]),
+            (0.15, -0.15, [qml.RZ(0.15, wires=0), qml.RZ(-0.15, wires=1)]),
+        ],
+    )
+    def test_two_qubits_rotation_no_merge(self, theta_1, theta_2, expected_ops):
+        """Test that a two-qubit circuit with rotations on different qubits
+        do not get merged."""
+
+        def qfunc():
+            qml.RZ(theta_1, wires=0)
+            qml.RZ(theta_2, wires=1)
+
+        transformed_qfunc = merge_rotations()(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        assert len(ops) == len(expected_ops)
+
+        for op_obtained, op_expected in zip(ops, expected_ops):
+            assert op_obtained.name == op_expected.name
+            assert np.allclose(op_obtained.parameters, op_expected.parameters)
+
+    def test_two_qubits_rotation_merge_tolerance(self):
+        """Test whether tolerance argument is respected for merging."""
+
+        def qfunc():
+            qml.RZ(1e-7, wires=0)
+            qml.RZ(-2e-7, wires=0)
+
+        # Try with default tolerance; these ops should still be applied
+        transformed_qfunc = merge_rotations()(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        assert len(ops) == 1
+        assert ops[0].name == "RZ"
+        assert ops[0].parameters[0] == -1e-7
+
+        # Now try with higher tolerance threshold; the ops should cancel
+        transformed_qfunc = merge_rotations(atol=1e-5)(qfunc)
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+        assert len(ops) == 0
+
+    @pytest.mark.parametrize(
+        ("theta_11", "theta_12", "theta_21", "theta_22", "expected_ops"),
+        [
+            (0.3, -0.2, 0.5, -0.8, [qml.RX(0.1, wires=0), qml.RY(-0.3, wires=1)]),
+            (0.3, -0.3, 0.7, -0.1, [qml.RY(0.6, wires=1)]),
+        ],
+    )
+    def test_two_qubits_rotation_no_merge(
+        self, theta_11, theta_12, theta_21, theta_22, expected_ops
+    ):
+        """Test that a two-qubit circuit with rotations on different qubits get merged."""
+
+        def qfunc():
+            qml.RX(theta_11, wires=0)
+            qml.RY(theta_21, wires=1)
+            qml.RX(theta_12, wires=0)
+            qml.RY(theta_22, wires=1)
+
+        transformed_qfunc = merge_rotations()(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        assert len(ops) == len(expected_ops)
+
+        for op_obtained, op_expected in zip(ops, expected_ops):
+            assert op_obtained.name == op_expected.name
+            assert np.allclose(op_obtained.parameters, op_expected.parameters)
+
+    @pytest.mark.parametrize(
+        ("theta_11", "theta_12", "theta_21", "theta_22", "expected_ops"),
+        [
+            (0.3, -0.2, 0.5, -0.8, [qml.CRX(0.5, wires=[0, 1]), qml.RY(-1.3, wires=1)]),
+            (0.3, 0.3, 0.7, -0.1, [qml.RY(-0.8, wires=1)]),
+        ],
+    )
+    def test_two_qubits_merge_with_adjoint(
+        self, theta_11, theta_12, theta_21, theta_22, expected_ops
+    ):
+        """Test that adjoint rotations on different qubits get merged."""
+
+        def qfunc():
+            qml.CRX(theta_11, wires=[0, 1])
+            qml.adjoint(qml.RY)(theta_21, wires=2)
+            qml.adjoint(qml.CRX)(theta_12, wires=[0, 1])
+            qml.RY(theta_22, wires=2)
+
+        transformed_qfunc = merge_rotations()(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        assert len(ops) == len(expected_ops)
+
+        for op_obtained, op_expected in zip(ops, expected_ops):
+            assert op_obtained.name == op_expected.name
+            assert np.allclose(op_obtained.parameters, op_expected.parameters)
+
+    def test_two_qubits_merge_gate_subset(self):
+        """Test that specifying a subset of operations to include merges correctly."""
+
+        def qfunc():
+            qml.CRX(0.1, wires=[0, 1])
+            qml.CRX(0.2, wires=[0, 1])
+            qml.RY(0.3, wires=["a"])
+            qml.RY(0.5, wires=["a"])
+            qml.RX(-0.5, wires=[2])
+            qml.RX(0.2, wires=[2])
+
+        transformed_qfunc = merge_rotations(include_gates=["RX", "CRX"])(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["CRX", "RY", "RY", "RX"]
+        wires_expected = [Wires([0, 1]), Wires("a"), Wires("a"), Wires(2)]
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+        assert qml.math.isclose(ops[0].parameters[0], 0.3)
+        assert qml.math.isclose(ops[1].parameters[0], 0.3)
+        assert qml.math.isclose(ops[2].parameters[0], 0.5)
+        assert qml.math.isclose(ops[3].parameters[0], -0.3)
+
+    def test_one_qubit_rotation_blocked(self):
+        """Test that rotations on one-qubit separated by a "blocking" operation don't merge."""
+
+        def qfunc():
+            qml.RX(0.5, wires=0)
+            qml.Hadamard(wires=0)
+            qml.RX(0.4, wires=0)
+
+        transformed_qfunc = merge_rotations()(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["RX", "Hadamard", "RX"]
+        wires_expected = [Wires(0)] * 3
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+        assert ops[0].parameters[0] == 0.5
+        assert ops[2].parameters[0] == 0.4
+
+    def test_two_qubits_rotation_blocked(self):
+        """Test that rotations on a two-qubit system separated by a "blocking" operation
+        don't merge."""
+
+        def qfunc():
+            qml.RX(-0.42, wires=0)
+            qml.CNOT(wires=[0, 1])
+            qml.RX(0.8, wires=0)
+
+        transformed_qfunc = merge_rotations()(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["RX", "CNOT", "RX"]
+        wires_expected = [Wires(0), Wires([0, 1]), Wires(0)]
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+        assert ops[0].parameters[0] == -0.42
+        assert ops[2].parameters[0] == 0.8
+
+    @pytest.mark.parametrize(
+        ("theta_1", "theta_2", "expected_ops"),
+        [
+            (0.3, -0.2, [qml.CRY(0.1, wires=["w1", "w2"])]),
+            (0.15, -0.15, []),
+        ],
+    )
+    def test_controlled_rotation_merge(self, theta_1, theta_2, expected_ops):
+        """Test that adjacent controlled rotations on the same wires in same order get merged."""
+
+        def qfunc():
+            qml.CRY(theta_1, wires=["w1", "w2"])
+            qml.CRY(theta_2, wires=["w1", "w2"])
+
+        transformed_qfunc = merge_rotations()(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        assert len(ops) == len(expected_ops)
+
+        # Check that all operations and parameter values are as expected
+        for op_obtained, op_expected in zip(ops, expected_ops):
+            assert op_obtained.name == op_expected.name
+            assert np.allclose(op_obtained.parameters, op_expected.parameters)
+
+    def test_controlled_rotation_no_merge(self):
+        """Test that adjacent controlled rotations on the same wires in different order don't merge."""
+
+        def qfunc():
+            qml.CRX(0.2, wires=["w1", "w2"])
+            qml.CRX(0.3, wires=["w2", "w1"])
+
+        transformed_qfunc = merge_rotations()(qfunc)
+
+        ops = qml.transforms.make_tape(transformed_qfunc)().operations
+
+        names_expected = ["CRX", "CRX"]
+        wires_expected = [Wires(["w1", "w2"]), Wires(["w2", "w1"])]
+        compare_operation_lists(ops, names_expected, wires_expected)
+
+        assert ops[0].parameters[0] == 0.2
+        assert ops[1].parameters[0] == 0.3
+
+
+# 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)
+    qml.Rot(0.0, 0.0, 0.0, wires=1)
+    return qml.expval(qml.PauliX(0) @ qml.PauliX(2))
+
+
+transformed_qfunc = merge_rotations()(qfunc)
+
+expected_op_list = ["Hadamard", "RZ", "PauliY", "CNOT", "CRY", "PauliZ", "Rot"]
+expected_wires_list = [
+    Wires(0),
+    Wires(0),
+    Wires(1),
+    Wires([1, 2]),
+    Wires([1, 2]),
+    Wires(0),
+    Wires(1),
+]
+
+
+class TestMergeRotationsInterfaces:
+    """Test that rotation merging works in all interfaces."""
+
+    def test_merge_rotations_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_merge_rotations_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_merge_rotations_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_merge_rotations_jax(self):
+        """Test QNode and gradient in JAX interface."""
+        jax = pytest.importorskip("jax")
+        from jax import numpy as jnp
+
+        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/test_optimization_utils.py b/tests/transforms/test_optimization/test_optimization_utils.py
new file mode 100644
index 000000000..49d04be68
--- /dev/null
+++ b/tests/transforms/test_optimization/test_optimization_utils.py
@@ -0,0 +1,104 @@
+# 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
+
+import pennylane as qml
+from pennylane.transforms.optimization.optimization_utils import (
+    find_next_gate,
+    _yzy_to_zyz,
+    fuse_rot_angles,
+)
+
+from utils import check_matrix_equivalence, compute_matrix_from_ops_one_qubit
+
+
+sample_op_list = [
+    qml.Hadamard(wires="a"),
+    qml.CNOT(wires=["a", "b"]),
+    qml.Toffoli(wires=[0, 1, "b"]),
+    qml.Hadamard(wires="c"),
+]
+
+
+class TestFindNextGate:
+    @pytest.mark.parametrize(
+        ("wires,op_list,next_gate_idx"),
+        [
+            ("a", sample_op_list, 0),
+            ("b", sample_op_list, 1),
+            ("e", sample_op_list, None),
+            ([0, 2], sample_op_list, 2),
+            ([0, 1], sample_op_list, 2),
+            ("c", sample_op_list, 3),
+        ],
+    )
+    def test_find_next_gate(self, wires, op_list, next_gate_idx):
+        """Test find_next_gate correctly identifies the next gate in a list of operations that share any number of wires."""
+        assert find_next_gate(qml.wires.Wires(wires), op_list) == next_gate_idx
+
+
+class TestRotGateFusion:
+    """Test that utility functions for fusing two qml.Rot gates function as expected."""
+
+    @pytest.mark.parametrize(
+        ("angles"),
+        [([0.15, 0.25, -0.90]), ([0.0, 0.0, 0.0]), ([0.15, 0.25, -0.90]), ([0.05, -1.34, 4.12])],
+    )
+    def test_yzy_to_zyz(self, angles):
+        """Test that a set of rotations of the form YZY is correctly converted
+        to a sequence of the form ZYZ."""
+
+        original_ops = [
+            qml.RY(angles[0], wires=0),
+            qml.RZ(angles[1], wires=0),
+            qml.RY(angles[2], wires=0),
+        ]
+        product_yzy = compute_matrix_from_ops_one_qubit(original_ops)
+
+        z1, y, z2 = _yzy_to_zyz(angles)
+        transformed_ops = [qml.RZ(z1, wires=0), qml.RY(y, wires=0), qml.RZ(z2, wires=0)]
+        product_zyz = compute_matrix_from_ops_one_qubit(transformed_ops)
+
+        assert check_matrix_equivalence(product_yzy, product_zyz)
+
+    @pytest.mark.parametrize(
+        ("angles_1", "angles_2"),
+        [
+            ([0.15, 0.25, -0.90], [-0.5, 0.25, 1.3]),
+            ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0]),
+            ([0.1, 0.2, 0.3], [0.0, 0.0, 0.0]),
+            ([0.0, 0.0, 0.0], [0.2, 0.4, -0.5]),
+            ([0.15, 0.25, -0.90], [-0.15, -0.25, 0.9]),
+            ([0.05, -1.34, 4.12], [-0.8, 0.2, 0.12]),
+            ([0.05, 0.0, 4.12], [-0.8, 0.2, 0.12]),
+            ([0.05, -1.34, 0.0], [-0.8, 0.2, 0.12]),
+            ([0.05, 0.0, 0.1], [-0.2, 0.0, 0.12]),
+            ([0.05, 0.0, 0.0], [0.0, 0.0, 0.12]),
+            ([0.05, 0.2, 0.0], [0.0, -0.6, 0.12]),
+            ([0.05, -1.34, 4.12], [0.0, 0.2, 0.12]),
+            ([0.05, -1.34, 4.12], [0.3, 0.0, 0.12]),
+        ],
+    )
+    def test_full_rot_fusion_autograd(self, angles_1, angles_2):
+        """Test that the fusion of two Rot gates has the same effect as
+        applying the Rots sequentially."""
+
+        original_ops = [qml.Rot(*angles_1, wires=0), qml.Rot(*angles_2, wires=0)]
+        matrix_expected = compute_matrix_from_ops_one_qubit(original_ops)
+
+        fused_angles = fuse_rot_angles(angles_1, angles_2)
+        matrix_obtained = qml.Rot(*fused_angles, wires=0).matrix
+
+        assert check_matrix_equivalence(matrix_expected, matrix_obtained)
diff --git a/tests/transforms/test_optimization/utils.py b/tests/transforms/test_optimization/utils.py
new file mode 100644
index 000000000..f5995f427
--- /dev/null
+++ b/tests/transforms/test_optimization/utils.py
@@ -0,0 +1,44 @@
+# 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.
+"""Convenient utility functions for testing optimization transforms."""
+
+import pennylane as qml
+
+from gate_data import I
+
+
+def compute_matrix_from_ops_one_qubit(ops):
+    """Given a list of single-qubit operations, construct its matrix representation."""
+
+    mat = I
+
+    for op in ops:
+        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."""
+
+    mat_product = qml.math.dot(qml.math.conj(matrix_obtained.T), matrix_expected)
+    mat_product /= mat_product[0, 0]
+    return qml.math.allclose(mat_product, qml.math.eye(matrix_expected.shape[0]))
+
+
+def compare_operation_lists(ops_obtained, names_expected, wires_expected):
+    """Compare two lists of operations."""
+    assert len(ops_obtained) == len(names_expected)
+    assert all([op.name == exp_name for (op, exp_name) in zip(ops_obtained, names_expected)])
+    assert all([op.wires == exp_wires for (op, exp_wires) in zip(ops_obtained, wires_expected)])