Repository: PennyLaneAI/pennylane
Pull request #3450 (2022-11-30 21:40:54)

ControlledQubitUnitary inherits from Controlled
**Context:**
ControlledQubitUnitary was defined as a subclass of `QubitUnitary`, and implemented a variety of methods that are duplicated elsewhere.

**Description of the Change:**
- Reorganize so that `ControlledQubitUnitary` inherits from the Controlled class.
- Add broadcasting to the `Controlled` class. 
- Move `ControlledQubitUnitary` from `matrix_ops.py` to `controlled_ops.py` (avoids circular import). 

**Benefits:**
Less duplicate code. Inherits functionality it did not previously have from `Controlled` (methods `decomposition`, `expand` and `sparse_matrix` now return a result instead of raising an error). Future operators inheriting from `Controlled` will support parameter broadcasting as long as their base does.

**Possible Drawbacks:**
The hyperparameters dictionary now has the same format as other Controlled operators, so `op.hyperparameters['u_wires']` is no longer an option to access the wires the unitary matrix operates on. This information can be accessed via `op.base.wires` or `op.target_wires` instead. Otherwise operator behaviour should be the same as before reorganizing, so there shouldn't be any other issues.


codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#3450](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (e472548) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/b3b0bf8526d279bc88dc0a284f4dd1671eef236d?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (b3b0bf8) will **decrease** coverage by `0.01%`.
> The diff coverage is `95.16%`.

```diff
@@            Coverage Diff             @@
##           master    #3450      +/-   ##
==========================================
- Coverage   99.87%   99.85%   -0.02%     
==========================================
  Files         314      315       +1     
  Lines       28019    28005      -14     
==========================================
- Hits        27983    27965      -18     
- Misses         36       40       +4     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/ops/qubit/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9fX2luaXRfXy5weQ==) | `100.00% <ø> (ø)` | |
| [pennylane/ops/op\_math/controlled\_class.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL2NvbnRyb2xsZWRfY2xhc3MucHk=) | `99.18% <93.33%> (-0.82%)` | :arrow_down: |
| [pennylane/ops/op\_math/controlled\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL2NvbnRyb2xsZWRfb3BzLnB5) | `96.29% <96.29%> (ø)` | |
| [pennylane/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL19faW5pdF9fLnB5) | `96.73% <100.00%> (ø)` | |
| [pennylane/ops/op\_math/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL19faW5pdF9fLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/qubit/matrix\_ops.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9tYXRyaXhfb3BzLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/drawer/tape\_mpl.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3450/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RyYXdlci90YXBlX21wbC5weQ==) | `98.83% <0.00%> (-1.17%)` | :arrow_down: |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

lillian542: Note: all the tests removed from `tests/ops/qubit/test_matrix_ops.py` were copy-pasted into `tests/ops/op_math/test_controlled_ops.py`. The tests were moved to reflect the new location of the operator they test (which had to be moved from the `qubit` folder into a file in the `math_op` folder because importing `Controlled` in the `qubit` folder created circular imports.
AlbertMitjans: [sc-31243]
lillian542: > It automatically imports everything in the `__all__` dunder list. I think it would be best to keep "ControlledQubitUnitary" in that list, and potentially even in `_qubit__ops__` as well.

Okay, thanks - I didn't know how that worked! 

I rearranged so there is a set in `ops.op_math.__init__` for this and any future `Controlled` ops, which will all then be included in `_qubit__ops__`,  ` _qubit__all__`, and the `__all__` dunder list for `ops`.


Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 6a11a13b1..4dcaa6a9a 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -132,6 +132,14 @@
 * Extended the `qml.equal` function to compare `Prod` and `Sum` operators.
   [(#3516)](https://github.com/PennyLaneAI/pennylane/pull/3516)
 
+* Reorganize `ControlledQubitUnitary` to inherit from `ControlledOp`. The class methods 
+  `decomposition`, `expand`, and `sparse_matrix` are now defined rather than raising an error.
+  [(#3450)](https://github.com/PennyLaneAI/pennylane/pull/3450)
+
+* Parameter broadcasting support is added for the `Controlled` class if the base operator supports 
+  broadcasting.
+  [(#3450)](https://github.com/PennyLaneAI/pennylane/pull/3450)
+
 * The `qml.generator` function now checks if the generator is hermitian, rather than whether it is a subclass of
   `Observable`, allowing it to return valid generators from `SymbolicOp` and `CompositeOp` classes.
  [(#3485)](https://github.com/PennyLaneAI/pennylane/pull/3485)
@@ -160,6 +168,10 @@
 
 <h3>Breaking changes</h3>
 
+* The target wires of the unitary for `ControlledQubitUnitary` are no longer available via `op.hyperparameters["u_wires"]`. 
+  Instead, they can be accesses via `op.base.wires` or `op.target_wires`.
+  [(#3450)](https://github.com/PennyLaneAI/pennylane/pull/3450)
+
 * The tape constructed by a QNode is no longer queued to surrounding contexts.
   [(#3509)](https://github.com/PennyLaneAI/pennylane/pull/3509)
 
diff --git a/pennylane/ops/__init__.py b/pennylane/ops/__init__.py
index a07e89109..8fbd75e26 100644
--- a/pennylane/ops/__init__.py
+++ b/pennylane/ops/__init__.py
@@ -29,16 +29,22 @@ from .cv import __all__ as _cv__all__
 from .cv import ops as _cv__ops__
 from .cv import obs as _cv__obs__
 
-from .qubit import __all__ as _qubit__all__
-from .qubit import ops as _qubit__ops__
-from .qubit import obs as _qubit__obs__
-
 from .qutrit import __all__ as _qutrit__all__
 from .qutrit import ops as _qutrit__ops__
 from .qutrit import obs as _qutrit__obs__
 
 from .channel import __all__ as _channel__ops__
 
+from .qubit import __all__ as _qubit__all__
+from .qubit import ops as _qubit__ops__
+from .qubit import obs as _qubit__obs__
+
+# we would like these to just live in .qubit, but can't because of circular imports
+from .op_math import controlled_qubit_ops as _controlled_qubit__ops__
+
+_qubit__ops__ = _qubit__ops__ | _controlled_qubit__ops__
+_qubit__all__ = _qubit__all__ + list(_controlled_qubit__ops__)
+
 
 __all__ = _cv__all__ + _qubit__all__ + _qutrit__all__ + _channel__ops__
 __all_ops__ = list(_cv__ops__ | _qubit__ops__ | _qutrit__ops__)
diff --git a/pennylane/ops/op_math/__init__.py b/pennylane/ops/op_math/__init__.py
index 0b912a715..7c988f339 100644
--- a/pennylane/ops/op_math/__init__.py
+++ b/pennylane/ops/op_math/__init__.py
@@ -51,11 +51,20 @@ Symbolic Classes
     ~SProd
     ~SymbolicOp
 
+Controlled Operator Classes
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. autosummary::
+    :toctree: api
+
+    ~ControlledQubitUnitary
+
 """
 
 from .adjoint import Adjoint, adjoint
 from .composite import CompositeOp
 from .controlled import Controlled, ControlledOp, ctrl
+from .controlled_ops import ControlledQubitUnitary
 from .dot import dot
 from .exp import Evolution, Exp, exp
 from .pow import Pow, pow
@@ -63,3 +72,8 @@ from .prod import Prod, prod
 from .sprod import SProd, s_prod
 from .sum import Sum, op_sum
 from .symbolicop import SymbolicOp
+
+
+controlled_qubit_ops = {
+    "ControlledQubitUnitary",
+}
diff --git a/pennylane/ops/op_math/controlled.py b/pennylane/ops/op_math/controlled.py
index ff9b02af6..d9741ee3f 100644
--- a/pennylane/ops/op_math/controlled.py
+++ b/pennylane/ops/op_math/controlled.py
@@ -140,6 +140,8 @@ class Controlled(SymbolicOp):
             length as ``control_wires``. Defaults to ``True`` for all control wires.
             Provided values are converted to `Bool` internally.
         work_wires (Any): Any auxiliary wires that can be used in the decomposition
+        do_queue(bool):  indicates whether the operator should be
+            recorded when created in a tape context
 
     .. note::
         This class, ``Controlled``, denotes a controlled version of any individual operation.
@@ -232,7 +234,7 @@ class Controlled(SymbolicOp):
 
     # pylint: disable=unused-argument
     def __new__(cls, base, *_, **__):
-        """If base is an ``Operation``, then the a ``ControlledOp`` should be used instead."""
+        """If base is an ``Operation``, then a ``ControlledOp`` should be used instead."""
         if isinstance(base, operation.Operation):
             return object.__new__(ControlledOp)
         return object.__new__(Controlled)
@@ -241,6 +243,7 @@ class Controlled(SymbolicOp):
     def __init__(
         self, base, control_wires, control_values=None, work_wires=None, do_queue=True, id=None
     ):
+
         control_wires = Wires(control_wires)
         work_wires = Wires([]) if work_wires is None else Wires(work_wires)
 
@@ -252,7 +255,8 @@ class Controlled(SymbolicOp):
                     "Specifying control values as a string is deprecated. Please use Sequence[Bool]",
                     UserWarning,
                 )
-                control_values = [(x == "1") for x in control_values]
+                # All values not 0 are cast as true. Assumes a string of 1s and 0s.
+                control_values = [(x != "0") for x in control_values]
 
             control_values = (
                 [bool(control_values)]
@@ -308,7 +312,19 @@ class Controlled(SymbolicOp):
     # pylint: disable=arguments-renamed, invalid-overridden-method
     @property
     def has_matrix(self):
-        return self.base.has_matrix if self.base.batch_size is None else False
+        return self.base.has_matrix
+
+    # pylint: disable=protected-access
+    def _check_batching(self, params):
+        self.base._check_batching(params)
+
+    @property
+    def batch_size(self):
+        return self.base.batch_size
+
+    @property
+    def ndim_params(self):
+        return self.base.ndim_params
 
     # Properties on the control values ######################
     @property
@@ -353,7 +369,7 @@ class Controlled(SymbolicOp):
         new_control_wires = Wires([wire_map.get(wire, wire) for wire in self.control_wires])
         new_work_wires = Wires([wire_map.get(wire, wire) for wire in self.work_wires])
 
-        return Controlled(
+        return self.__class__(
             base=new_base,
             control_wires=new_control_wires,
             control_values=self.control_values,
@@ -375,9 +391,6 @@ class Controlled(SymbolicOp):
 
     def matrix(self, wire_order=None):
 
-        if self.base.batch_size is not None:
-            raise qml.operation.MatrixUndefinedError
-
         base_matrix = self.base.matrix()
         interface = qmlmath.get_interface(base_matrix)
 
@@ -391,7 +404,13 @@ class Controlled(SymbolicOp):
         left_pad = qmlmath.cast_like(qmlmath.eye(padding_left, like=interface), 1j)
         right_pad = qmlmath.cast_like(qmlmath.eye(padding_right, like=interface), 1j)
 
-        canonical_matrix = qmlmath.block_diag([left_pad, base_matrix, right_pad])
+        shape = qml.math.shape(base_matrix)
+        if len(shape) == 3:  # stack if batching
+            canonical_matrix = qml.math.stack(
+                [qml.math.block_diag([left_pad, _U, right_pad]) for _U in base_matrix]
+            )
+        else:
+            canonical_matrix = qmlmath.block_diag([left_pad, base_matrix, right_pad])
 
         if wire_order is None or self.wires == Wires(wire_order):
             return qml.math.expand_matrix(
@@ -489,27 +508,35 @@ class Controlled(SymbolicOp):
         return self.base.has_adjoint
 
     def adjoint(self):
-        return Controlled(
-            self.base.adjoint(), self.control_wires, self.control_values, self.work_wires
+        return self.__class__(
+            self.base.adjoint(),
+            self.control_wires,
+            control_values=self.control_values,
+            work_wires=self.work_wires,
         )
 
     def pow(self, z):
         base_pow = self.base.pow(z)
         return [
-            Controlled(op, self.control_wires, self.control_values, self.work_wires)
+            self.__class__(
+                op,
+                self.control_wires,
+                control_values=self.control_values,
+                work_wires=self.work_wires,
+            )
             for op in base_pow
         ]
 
     def simplify(self) -> "Controlled":
         if isinstance(self.base, Controlled):
             base = self.base.base.simplify()
-            return Controlled(
+            return self.__class__(
                 base,
                 control_wires=self.control_wires + self.base.control_wires,
                 control_values=self.control_values + self.base.control_values,
                 work_wires=self.work_wires + self.base.work_wires,
             )
-        return Controlled(
+        return self.__class__(
             base=self.base.simplify(),
             control_wires=self.control_wires,
             control_values=self.control_values,
diff --git a/pennylane/ops/op_math/controlled_ops.py b/pennylane/ops/op_math/controlled_ops.py
new file mode 100644
index 000000000..e7f81899a
--- /dev/null
+++ b/pennylane/ops/op_math/controlled_ops.py
@@ -0,0 +1,138 @@
+# Copyright 2018-2022 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 submodule contains controlled operators based on the Controlled and ControlledOp class.
+"""
+
+import warnings
+from typing import Iterable
+
+from pennylane.operation import AnyWires
+from pennylane.ops.qubit.matrix_ops import QubitUnitary
+from .controlled import ControlledOp
+
+
+# pylint: disable=too-few-public-methods
+class ControlledQubitUnitary(ControlledOp):
+    r"""ControlledQubitUnitary(U, control_wires, wires, control_values)
+    Apply an arbitrary fixed unitary to ``wires`` with control from the ``control_wires``.
+
+    In addition to default ``Operation`` instance attributes, the following are
+    available for ``ControlledQubitUnitary``:
+
+    * ``control_wires``: wires that act as control for the operation
+    * ``control_values``: the state on which to apply the controlled operation (see below)
+    * ``target_wires``: the wires the unitary matrix will be applied to
+    * ``active_wires``: Wires modified by the operator. This is the control wires followed by the target wires.
+
+    **Details:**
+
+    * Number of wires: Any (the operation can act on any number of wires)
+    * Number of parameters: 1
+    * Number of dimensions per parameter: (2,)
+    * Gradient recipe: None
+
+    Args:
+        base (Union[array[complex], QubitUnitary]): square unitary matrix or a QubitUnitary operation. If passing a matrix,
+          this will be used to construct a QubitUnitary operator that will be used as the base operator. If providing
+          a ``qml.QubitUnitary``, this will be used as the base directly.
+        control_wires (Union[Wires, Sequence[int], or int]): the control wire(s)
+        wires (Union[Wires, Sequence[int], or int]): the wire(s) the unitary acts on (optional if U is provided as a QubitUnitary)
+        control_values (List[int, bool]): a list providing the state of the control qubits to control on (default is the all 1s state)
+        unitary_check (bool): whether to check whether an array U is unitary when creating the operator (default False)
+        do_queue (bool): indicates whether the operator should be recorded when created in a tape context
+
+    **Example**
+
+    The following shows how a single-qubit unitary can be applied to wire ``2`` with control on
+    both wires ``0`` and ``1``:
+
+    >>> U = np.array([[ 0.94877869,  0.31594146], [-0.31594146,  0.94877869]])
+    >>> qml.ControlledQubitUnitary(U, control_wires=[0, 1], wires=2)
+    Controlled(QubitUnitary(array([[ 0.94877869,  0.31594146],
+       [-0.31594146,  0.94877869]]), wires=[2]), control_wires=[0, 1])
+
+    Alternatively, the same operator can be constructed with a QubitUnitary:
+
+    >>> base = qml.QubitUnitary(U, wires=2)
+    >>> qml.ControlledQubitUnitary(base, control_wires=[0, 1])
+    Controlled(QubitUnitary(array([[ 0.94877869,  0.31594146],
+       [-0.31594146,  0.94877869]]), wires=[2]), control_wires=[0, 1])
+
+    Typically controlled operations apply a desired gate if the control qubits
+    are all in the state :math:`\vert 1\rangle`. However, there are some situations where
+    it is necessary to apply a gate conditioned on all qubits being in the
+    :math:`\vert 0\rangle` state, or a mix of the two.
+
+    The state on which to control can be changed by passing a string of bits to
+    `control_values`. For example, if we want to apply a single-qubit unitary to
+    wire ``3`` conditioned on three wires where the first is in state ``0``, the
+    second is in state ``1``, and the third in state ``1``, we can write:
+
+    >>> qml.ControlledQubitUnitary(U, control_wires=[0, 1, 2], wires=3, control_values=[0, 1, 1])
+
+    or
+
+    >>> qml.ControlledQubitUnitary(U, control_wires=[0, 1, 2], wires=3, control_values=[False, True, True])
+    """
+    num_wires = AnyWires
+    """int: Number of wires that the operator acts on."""
+
+    num_params = 1
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = (2,)
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    grad_method = None
+    """Gradient computation method."""
+
+    # pylint: disable= too-many-arguments
+    def __init__(
+        self,
+        base,
+        control_wires,
+        wires=None,
+        control_values=None,
+        do_queue=True,
+        unitary_check=False,
+        work_wires=None,
+    ):
+        if getattr(base, "wires", False) and wires is not None:
+            warnings.warn(
+                "base operator already has wires; values specified through wires kwarg will be ignored."
+            )
+
+        if isinstance(base, Iterable):
+            base = QubitUnitary(base, wires=wires, unitary_check=unitary_check)
+
+        super().__init__(
+            base,
+            control_wires,
+            control_values=control_values,
+            work_wires=work_wires,
+            do_queue=do_queue,
+        )
+        self._name = "ControlledQubitUnitary"
+
+    def _controlled(self, wire):
+        ctrl_wires = self.control_wires + wire
+        values = None if self.control_values is None else self.control_values + [True]
+        new_op = ControlledQubitUnitary(
+            self.base,
+            control_wires=ctrl_wires,
+            control_values=values,
+            work_wires=self.work_wires,
+        )
+        return new_op
diff --git a/pennylane/ops/qubit/__init__.py b/pennylane/ops/qubit/__init__.py
index 6e427471f..658fa8d4d 100644
--- a/pennylane/ops/qubit/__init__.py
+++ b/pennylane/ops/qubit/__init__.py
@@ -89,7 +89,6 @@ ops = {
     "QubitStateVector",
     "QubitDensityMatrix",
     "QubitUnitary",
-    "ControlledQubitUnitary",
     "MultiControlledX",
     "IntegerComparator",
     "DiagonalQubitUnitary",
diff --git a/pennylane/ops/qubit/matrix_ops.py b/pennylane/ops/qubit/matrix_ops.py
index ef54ee8a4..0383320c7 100644
--- a/pennylane/ops/qubit/matrix_ops.py
+++ b/pennylane/ops/qubit/matrix_ops.py
@@ -17,7 +17,6 @@ accept a hermitian or an unitary matrix as a parameter.
 """
 # pylint:disable=arguments-differ
 import warnings
-
 import numpy as np
 
 import pennylane as qml
@@ -74,36 +73,36 @@ class QubitUnitary(Operation):
     def __init__(
         self, U, wires, do_queue=True, id=None, unitary_check=False
     ):  # pylint: disable=too-many-arguments
-        wires = Wires(wires)
-
         # For pure QubitUnitary operations (not controlled), check that the number
         # of wires fits the dimensions of the matrix
-        if not isinstance(self, ControlledQubitUnitary):
-            U_shape = qml.math.shape(U)
 
-            dim = 2 ** len(wires)
+        wires = Wires(wires)
 
-            if len(U_shape) not in {2, 3} or U_shape[-2:] != (dim, dim):
-                raise ValueError(
-                    f"Input unitary must be of shape {(dim, dim)} or (batch_size, {dim}, {dim}) "
-                    f"to act on {len(wires)} wires."
-                )
+        U_shape = qml.math.shape(U)
 
-            # Check for unitarity; due to variable precision across the different ML frameworks,
-            # here we issue a warning to check the operation, instead of raising an error outright.
-            if unitary_check and not (
-                qml.math.is_abstract(U)
-                or qml.math.allclose(
-                    qml.math.einsum("...ij,...kj->...ik", U, qml.math.conj(U)),
-                    qml.math.eye(dim),
-                    atol=1e-6,
-                )
-            ):
-                warnings.warn(
-                    f"Operator {U}\n may not be unitary."
-                    "Verify unitarity of operation, or use a datatype with increased precision.",
-                    UserWarning,
-                )
+        dim = 2 ** len(wires)
+
+        if len(U_shape) not in {2, 3} or U_shape[-2:] != (dim, dim):
+            raise ValueError(
+                f"Input unitary must be of shape {(dim, dim)} or (batch_size, {dim}, {dim}) "
+                f"to act on {len(wires)} wires."
+            )
+
+        # Check for unitarity; due to variable precision across the different ML frameworks,
+        # here we issue a warning to check the operation, instead of raising an error outright.
+        if unitary_check and not (
+            qml.math.is_abstract(U)
+            or qml.math.allclose(
+                qml.math.einsum("...ij,...kj->...ik", U, qml.math.conj(U)),
+                qml.math.eye(dim),
+                atol=1e-6,
+            )
+        ):
+            warnings.warn(
+                f"Operator {U}\n may not be unitary."
+                "Verify unitarity of operation, or use a datatype with increased precision.",
+                UserWarning,
+            )
 
         super().__init__(U, wires=wires, do_queue=do_queue, id=id)
 
@@ -190,210 +189,13 @@ class QubitUnitary(Operation):
         return super().pow(z)
 
     def _controlled(self, wire):
-        new_op = ControlledQubitUnitary(*self.parameters, control_wires=wire, wires=self.wires)
+        new_op = qml.ControlledQubitUnitary(*self.parameters, control_wires=wire, wires=self.wires)
         return new_op.inv() if self.inverse else new_op
 
     def label(self, decimals=None, base_label=None, cache=None):
         return super().label(decimals=decimals, base_label=base_label or "U", cache=cache)
 
 
-class ControlledQubitUnitary(QubitUnitary):
-    r"""ControlledQubitUnitary(U, control_wires, wires, control_values)
-    Apply an arbitrary fixed unitary to ``wires`` with control from the ``control_wires``.
-
-    In addition to default ``Operation`` instance attributes, the following are
-    available for ``ControlledQubitUnitary``:
-
-    * ``control_wires``: wires that act as control for the operation
-    * ``U``: unitary applied to the target wires
-
-    **Details:**
-
-    * Number of wires: Any (the operation can act on any number of wires)
-    * Number of parameters: 1
-    * Number of dimensions per parameter: (2,)
-    * Gradient recipe: None
-
-    Args:
-        U (array[complex]): square unitary matrix
-        control_wires (Union[Wires, Sequence[int], or int]): the control wire(s)
-        wires (Union[Wires, Sequence[int], or int]): the wire(s) the unitary acts on
-        control_values (str): a string of bits representing the state of the control
-            qubits to control on (default is the all 1s state)
-
-    **Example**
-
-    The following shows how a single-qubit unitary can be applied to wire ``2`` with control on
-    both wires ``0`` and ``1``:
-
-    >>> U = np.array([[ 0.94877869,  0.31594146], [-0.31594146,  0.94877869]])
-    >>> qml.ControlledQubitUnitary(U, control_wires=[0, 1], wires=2)
-
-    Typically controlled operations apply a desired gate if the control qubits
-    are all in the state :math:`\vert 1\rangle`. However, there are some situations where
-    it is necessary to apply a gate conditioned on all qubits being in the
-    :math:`\vert 0\rangle` state, or a mix of the two.
-
-    The state on which to control can be changed by passing a string of bits to
-    `control_values`. For example, if we want to apply a single-qubit unitary to
-    wire ``3`` conditioned on three wires where the first is in state ``0``, the
-    second is in state ``1``, and the third in state ``1``, we can write:
-
-    >>> qml.ControlledQubitUnitary(U, control_wires=[0, 1, 2], wires=3, control_values='011')
-
-    """
-    num_wires = AnyWires
-    """int: Number of wires that the operator acts on."""
-
-    num_params = 1
-    """int: Number of trainable parameters that the operator depends on."""
-
-    ndim_params = (2,)
-    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
-
-    grad_method = None
-    """Gradient computation method."""
-
-    def __init__(
-        self,
-        *params,
-        control_wires=None,
-        wires=None,
-        control_values=None,
-        do_queue=True,
-    ):
-        if control_wires is None:
-            raise ValueError("Must specify control wires")
-
-        wires = Wires(wires)
-        control_wires = Wires(control_wires)
-
-        if Wires.shared_wires([wires, control_wires]):
-            raise ValueError(
-                "The control wires must be different from the wires specified to apply the unitary on."
-            )
-
-        self._hyperparameters = {
-            "u_wires": wires,
-            "control_wires": control_wires,
-            "control_values": control_values,
-        }
-
-        total_wires = control_wires + wires
-        super().__init__(*params, wires=total_wires, do_queue=do_queue)
-
-    @staticmethod
-    def compute_decomposition(*params, wires=None, **hyperparameters):
-        raise DecompositionUndefinedError
-
-    @staticmethod
-    def compute_matrix(
-        U, control_wires, u_wires, control_values=None
-    ):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.ControlledQubitUnitary.matrix`
-
-        Args:
-            U (tensor_like): unitary matrix
-            control_wires (Iterable): the control wire(s)
-            u_wires (Iterable): the wire(s) the unitary acts on
-            control_values (str or None): a string of bits representing the state of the control
-                qubits to control on (default is the all 1s state)
-
-        Returns:
-            tensor_like: canonical matrix
-
-        **Example**
-
-        >>> U = np.array([[ 0.94877869,  0.31594146], [-0.31594146,  0.94877869]])
-        >>> qml.ControlledQubitUnitary.compute_matrix(U, control_wires=[1], u_wires=[0], control_values="1")
-        [[ 1.        +0.j  0.        +0.j  0.        +0.j  0.        +0.j]
-         [ 0.        +0.j  1.        +0.j  0.        +0.j  0.        +0.j]
-         [ 0.        +0.j  0.        +0.j  0.94877869+0.j  0.31594146+0.j]
-         [ 0.        +0.j  0.        +0.j -0.31594146+0.j  0.94877869+0.j]]
-        """
-        target_dim = 2 ** len(u_wires)
-        shape = qml.math.shape(U)
-        if not (len(shape) in {2, 3} and shape[-2:] == (target_dim, target_dim)):
-            raise ValueError(
-                f"Input unitary must be of shape {(target_dim, target_dim)} or "
-                f"(batch_size, {target_dim}, {target_dim})."
-            )
-
-        # A multi-controlled operation is a block-diagonal matrix partitioned into
-        # blocks where the operation being applied sits in the block positioned at
-        # the integer value of the control string. For example, controlling a
-        # unitary U with 2 qubits will produce matrices with block structure
-        # (U, I, I, I) if the control is on bits '00', (I, U, I, I) if on bits '01',
-        # etc. The positioning of the block is controlled by padding the block diagonal
-        # to the left and right with the correct amount of identity blocks.
-
-        total_wires = qml.wires.Wires(control_wires) + qml.wires.Wires(u_wires)
-
-        # if control values unspecified, we control on the all-ones string
-        if not control_values:
-            control_values = "1" * len(control_wires)
-
-        if isinstance(control_values, str):
-            if len(control_values) != len(control_wires):
-                raise ValueError("Length of control bit string must equal number of control wires.")
-
-            # Make sure all values are either 0 or 1
-            if not set(control_values).issubset({"0", "1"}):
-                raise ValueError("String of control values can contain only '0' or '1'.")
-
-            control_int = int(control_values, 2)
-        else:
-            raise ValueError("Alternative control values must be passed as a binary string.")
-
-        padding_left = control_int * target_dim
-        padding_right = 2 ** len(total_wires) - target_dim - padding_left
-
-        interface = qml.math.get_interface(U)
-        left_pad = qml.math.cast_like(qml.math.eye(padding_left, like=interface), 1j)
-        right_pad = qml.math.cast_like(qml.math.eye(padding_right, like=interface), 1j)
-        if len(qml.math.shape(U)) == 3:
-            return qml.math.stack([qml.math.block_diag([left_pad, _U, right_pad]) for _U in U])
-        return qml.math.block_diag([left_pad, U, right_pad])
-
-    @property
-    def control_wires(self):
-        return self.hyperparameters["control_wires"]
-
-    @property
-    def control_values(self):
-        """str.  Specifies whether or not to control on zero "0" or one "1" for each
-        control wire."""
-        return self.hyperparameters["control_values"]
-
-    def pow(self, z):
-        if isinstance(z, int):
-            return [
-                ControlledQubitUnitary(
-                    qml.math.linalg.matrix_power(self.data[0], z),
-                    control_wires=self.control_wires,
-                    control_values=self.control_values,
-                    wires=self.hyperparameters["u_wires"],
-                )
-            ]
-        return super().pow(z)
-
-    def _controlled(self, wire):
-        ctrl_wires = self.control_wires + wire
-        values = None if self.control_values is None else f"{self.control_values}1"
-        new_op = ControlledQubitUnitary(
-            *self.parameters,
-            control_wires=ctrl_wires,
-            wires=self.hyperparameters["u_wires"],
-            control_values=values,
-        )
-        return new_op.inv() if self.inverse else new_op
-
-
 class DiagonalQubitUnitary(Operation):
     r"""DiagonalQubitUnitary(D, wires)
     Apply an arbitrary fixed diagonal unitary matrix.


Test patch:
diff --git a/pennylane/devices/tests/test_gates.py b/pennylane/devices/tests/test_gates.py
old mode 100755
new mode 100644
index c6606d57d..6fa3e9de5
--- a/pennylane/devices/tests/test_gates.py
+++ b/pennylane/devices/tests/test_gates.py
@@ -15,7 +15,10 @@
 Tests that application of gates and state preparations
 works correctly an a device.
 """
-# pylint:disable=no-self-use,too-many-arguments,pointless-statement,unnecessary-lambda-assignment
+# pylint: disable=no-self-use
+# pylint: disable=too-many-arguments
+# pylint: disable=pointless-statement
+# pylint: disable=unnecessary-lambda-assignment
 from cmath import exp
 from math import cos, sin, sqrt
 
diff --git a/tests/ops/op_math/test_controlled.py b/tests/ops/op_math/test_controlled.py
index 2ffece6f3..495e4f569 100644
--- a/tests/ops/op_math/test_controlled.py
+++ b/tests/ops/op_math/test_controlled.py
@@ -11,6 +11,7 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
+"""Unit tests for Controlled"""
 
 from copy import copy
 from functools import partial
@@ -33,6 +34,11 @@ from pennylane.tape import QuantumScript
 from pennylane.tape.tape import expand_tape
 from pennylane.wires import Wires
 
+# pylint: disable=too-few-public-methods
+# pylint: disable=protected-access
+# pylint: disable=pointless-statement
+# pylint: disable=expression-not-assigned
+
 base_num_control_mats = [
     (qml.PauliX("a"), 1, CNOT),
     (qml.PauliZ("a"), 1, CZ),
@@ -78,7 +84,7 @@ class TestControlledInheritance:
         base = CustomOp(1.234, wires=0)
         op = Controlled(base, 6.5)
 
-        assert type(op) is ControlledOp
+        assert type(op) is ControlledOp  # pylint: disable=unidiomatic-typecheck
 
         assert isinstance(op, Controlled)
         assert isinstance(op, Operator)
@@ -92,7 +98,7 @@ class TestControlledInheritance:
         base = TempOperator(1.234, wires="a")
         op = ControlledOp(base, "b")
 
-        assert type(op) is ControlledOp
+        assert type(op) is ControlledOp  # pylint: disable=unidiomatic-typecheck
 
 
 class TestInitialization:
@@ -217,13 +223,14 @@ class TestProperties:
         op = Controlled(DummyOp(1), 0)
         assert op.has_matrix is value
 
-    def test_no_matrix_if_batching(self):
-        """Test that has_matrix is false if the base operator has batching."""
+    @pytest.mark.parametrize(
+        "base", (qml.RX(1.23, 0), qml.Rot(1.2, 2.3, 3.4, 0), qml.QubitUnitary([[0, 1], [1, 0]], 0))
+    )
+    def test_ndim_params(self, base):
+        """Test that Controlled defers to base ndim_params"""
 
-        x = np.array([1.2, 3.4, 5.6])
-        base = qml.RX(x, wires=0)
         op = Controlled(base, 1)
-        assert op.has_matrix is False
+        assert op.ndim_params == base.ndim_params
 
     @pytest.mark.parametrize("cwires, cvalues", [(0, [0]), ([3, 0, 2], [1, 1, 0])])
     def test_has_decomposition_true_via_control_values(self, cwires, cvalues):
@@ -473,8 +480,8 @@ class TestOperationProperties:
 
         op.inv()
 
-        assert op.inverse == False
-        assert base.inverse == True
+        assert op.inverse is False
+        assert base.inverse is True
         assert op.name == "C(S.inv)"
         assert op.base_name == "C(S)"
 
@@ -488,8 +495,8 @@ class TestOperationProperties:
 
         op.inverse = True
 
-        assert op.inverse == False
-        assert base.inverse == True
+        assert op.inverse is False
+        assert base.inverse is True
         assert op.name == "C(T.inv)"
         assert op.base_name == "C(T)"
 
@@ -624,7 +631,7 @@ class TestQueuing:
 
         base = qml.PauliX(0)
         with qml.queuing.AnnotatedQueue() as q:
-            op = Controlled(base, 1, do_queue=False)
+            _ = Controlled(base, 1, do_queue=False)
 
         assert len(q) == 0
 
@@ -644,14 +651,15 @@ base_num_control_mats = [
 class TestMatrix:
     """Tests of Controlled.matrix and Controlled.sparse_matrix"""
 
-    def test_raises_error_if_batching(self):
-        """Test a MatrixUndefinedError is raised if the base is batched."""
+    def test_correct_matrix_dimenions_with_batching(self):
+        """Test batching returns a matrix of the correct dimensions"""
         x = np.array([1.0, 2.0, 3.0])
         base = qml.RX(x, 0)
         op = Controlled(base, 1)
 
-        with pytest.raises(qml.operation.MatrixUndefinedError):
-            op.matrix()
+        matrix = op.matrix()
+
+        assert matrix.shape == (3, 4, 4)
 
     @pytest.mark.parametrize("base, num_control, mat", base_num_control_mats)
     def test_matrix_compare_with_gate_data(self, base, num_control, mat):
@@ -819,7 +827,7 @@ class TestDecomposition:
         control_values = [True, False, False]
 
         base = TempOperator("a")
-        op = Controlled(base, (0, 1, 2), [True, False, False])
+        op = Controlled(base, control_wires, control_values)
 
         decomp = op.expand().circuit if test_expand else op.decomposition()
 
@@ -871,6 +879,7 @@ class TestDecomposition:
         op = Controlled(base, (0, 1, 2))
 
         with pytest.raises(DecompositionUndefinedError):
+            # pylint: disable=unused-variable
             decomp = op.expand().circuit if test_expand else op.decomposition()
 
 
@@ -930,10 +939,9 @@ class TestDifferentiation:
     @pytest.mark.autograd
     def test_autograd(self, diff_method):
         """Test differentiation using autograd"""
-        from pennylane import numpy as pnp
 
         dev = qml.device("default.qubit", wires=2)
-        init_state = pnp.array([1.0, -1.0], requires_grad=False) / np.sqrt(2)
+        init_state = np.array([1.0, -1.0], requires_grad=False) / np.sqrt(2)
 
         @qml.qnode(dev, diff_method=diff_method)
         def circuit(b):
@@ -941,7 +949,7 @@ class TestDifferentiation:
             Controlled(qml.RY(b, wires=1), control_wires=0)
             return qml.expval(qml.PauliX(0))
 
-        b = pnp.array(0.123, requires_grad=True)
+        b = np.array(0.123, requires_grad=True)
         res = qml.grad(circuit)(b)
         expected = np.sin(b / 2) / 2
 
@@ -1019,6 +1027,313 @@ class TestDifferentiation:
         assert np.allclose(res, expected)
 
 
+class TestControlledSupportsBroadcasting:
+    """Test that the Controlled version of qubit operations with the ``supports_broadcasting`` attribute
+    actually support broadcasting."""
+
+    single_scalar_single_wire_ops = [
+        "RX",
+        "RY",
+        "RZ",
+        "PhaseShift",
+        "U1",
+    ]
+
+    single_scalar_multi_wire_ops = [
+        "ControlledPhaseShift",
+        "CRX",
+        "CRY",
+        "CRZ",
+        "IsingXX",
+        "IsingYY",
+        "IsingZZ",
+        "IsingXY",
+        "SingleExcitation",
+        "SingleExcitationPlus",
+        "SingleExcitationMinus",
+        "DoubleExcitation",
+        "DoubleExcitationPlus",
+        "DoubleExcitationMinus",
+        "OrbitalRotation",
+        "FermionicSWAP",
+    ]
+
+    two_scalar_single_wire_ops = [
+        "U2",
+    ]
+
+    three_scalar_single_wire_ops = [
+        "Rot",
+        "U3",
+    ]
+
+    three_scalar_multi_wire_ops = [
+        "CRot",
+    ]
+
+    # When adding an operation to the following list, you
+    # actually need to write a new test!
+    separately_tested_ops = [
+        "QubitUnitary",
+        "ControlledQubitUnitary",
+        "DiagonalQubitUnitary",
+        "PauliRot",
+        "MultiRZ",
+        "QubitStateVector",
+        "AmplitudeEmbedding",
+        "AngleEmbedding",
+        "IQPEmbedding",
+        "QAOAEmbedding",
+    ]
+
+    @pytest.mark.parametrize("name", single_scalar_single_wire_ops)
+    def test_controlled_of_single_scalar_single_wire_ops(self, name):
+        """Test that a Controlled operation whose base is a single-scalar-parameter operations
+        on a single wire marked as supporting parameter broadcasting actually do support broadcasting."""
+        par = np.array([0.25, 2.1, -0.42])
+        wires = ["wire0"]
+
+        cls = getattr(qml, name)
+        base = cls(par, wires=wires)
+        op = Controlled(base, "wire1")
+
+        mat = op.matrix()
+        single_mats = [Controlled(cls(p, wires=wires), "wire1").matrix() for p in par]
+
+        assert qml.math.allclose(mat, single_mats)
+
+    @pytest.mark.parametrize("name", single_scalar_multi_wire_ops)
+    def test_controlled_single_scalar_multi_wire_ops(self, name):
+        """Test that a Controlled operation whose base is a single-scalar-parameter operations
+        on multiple wires marked as supporting parameter broadcasting actually do support broadcasting."""
+        par = np.array([0.25, 2.1, -0.42])
+        cls = getattr(qml, name)
+
+        # Provide up to 6 wires and take as many as the class requires
+        # This assumes that the class does *not* have `num_wires=qml.operation.AnyWires`
+        wires = ["wire0", 5, 41, "aux_wire", -1, 9][: cls.num_wires]
+        base = cls(par, wires=wires)
+        op = Controlled(base, "wire1")
+
+        mat = op.matrix()
+        single_mats = [Controlled(cls(p, wires=wires), "wire1").matrix() for p in par]
+
+        assert qml.math.allclose(mat, single_mats)
+
+    @pytest.mark.parametrize("name", two_scalar_single_wire_ops)
+    def test_controlled_two_scalar_single_wire_ops(self, name):
+        """Test that a Controlled operation whose base is a two-scalar-parameter operations
+        on a single wire marked as supporting parameter broadcasting actually do support broadcasting."""
+        par = (np.array([0.25, 2.1, -0.42]), np.array([-6.2, 0.12, 0.421]))
+        wires = ["wire0"]
+
+        cls = getattr(qml, name)
+        base = cls(*par, wires=wires)
+        op = Controlled(base, "wire1")
+
+        mat = op.matrix()
+        single_pars = [tuple(p[i] for p in par) for i in range(3)]
+        single_mats = [Controlled(cls(*p, wires=wires), "wire1").matrix() for p in single_pars]
+
+        assert qml.math.allclose(mat, single_mats)
+
+    @pytest.mark.parametrize("name", three_scalar_single_wire_ops)
+    def test_controlled_three_scalar_single_wire_ops(self, name):
+        """Test that a Controlled operation whose base is a three-scalar-parameter operations
+        on a single wire marked as supporting parameter broadcasting actually do support broadcasting."""
+        par = (
+            np.array([0.25, 2.1, -0.42]),
+            np.array([-6.2, 0.12, 0.421]),
+            np.array([0.2, 1.1, -5.2]),
+        )
+        wires = ["wire0"]
+
+        cls = getattr(qml, name)
+        base = cls(*par, wires=wires)
+        op = Controlled(base, "wire1")
+
+        mat = op.matrix()
+        single_pars = [tuple(p[i] for p in par) for i in range(3)]
+        single_mats = [Controlled(cls(*p, wires=wires), "wire1").matrix() for p in single_pars]
+
+        assert qml.math.allclose(mat, single_mats)
+
+    @pytest.mark.parametrize("name", three_scalar_multi_wire_ops)
+    def test_controlled_three_scalar_multi_wire_ops(self, name):
+        """Test that a Controlled operation whose base is a three-scalar-parameter operations
+        on multiple wires marked as supporting parameter broadcasting actually do support broadcasting."""
+        par = (
+            np.array([0.25, 2.1, -0.42]),
+            np.array([-6.2, 0.12, 0.421]),
+            np.array([0.2, 1.1, -5.2]),
+        )
+        wires = ["wire0", 214]
+
+        cls = getattr(qml, name)
+        base = cls(*par, wires=wires)
+        op = Controlled(base, "wire1")
+
+        mat = op.matrix()
+        single_pars = [tuple(p[i] for p in par) for i in range(3)]
+        single_mats = [Controlled(cls(*p, wires=wires), "wire1").matrix() for p in single_pars]
+
+        assert qml.math.allclose(mat, single_mats)
+
+    def test_controlled_diagonal_qubit_unitary(self):
+        """Test that a Controlled operation whose base is a DiagonalQubitUnitary, which is marked
+        as supporting parameter broadcasting, actually does support broadcasting."""
+        diag = np.array([[1j, 1, 1, -1j], [-1j, 1j, 1, -1], [1j, -1j, 1.0, -1]])
+        wires = ["a", 5]
+
+        base = qml.DiagonalQubitUnitary(diag, wires=wires)
+        op = Controlled(base, "wire1")
+
+        mat = op.matrix()
+        single_mats = [
+            Controlled(qml.DiagonalQubitUnitary(d, wires=wires), "wire1").matrix() for d in diag
+        ]
+
+        assert qml.math.allclose(mat, single_mats)
+
+    @pytest.mark.parametrize(
+        "pauli_word, wires", [("XYZ", [0, "4", 1]), ("II", [1, 5]), ("X", [7])]
+    )
+    def test_controlled_pauli_rot(self, pauli_word, wires):
+        """Test that a Controlled operation whose base is PauliRot, which is marked as supporting
+        parameter broadcasting, actually does support broadcasting."""
+        par = np.array([0.25, 2.1, -0.42])
+
+        base = qml.PauliRot(par, pauli_word, wires=wires)
+        op = Controlled(base, "wire1")
+
+        mat = op.matrix()
+        single_mats = [
+            Controlled(qml.PauliRot(p, pauli_word, wires=wires), "wire1").matrix() for p in par
+        ]
+
+        assert qml.math.allclose(mat, single_mats)
+
+    @pytest.mark.parametrize("wires", [[0, "4", 1], [1, 5], [7]])
+    def test_controlled_multi_rz(self, wires):
+        """Test that a Controlled operation whose base is MultiRZ, which is marked as supporting
+        parameter broadcasting, actually does support broadcasting."""
+        par = np.array([0.25, 2.1, -0.42])
+
+        base = qml.MultiRZ(par, wires=wires)
+        op = Controlled(base, "wire1")
+
+        mat = op.matrix()
+        single_mats = [Controlled(qml.MultiRZ(p, wires=wires), "wire1").matrix() for p in par]
+
+        assert qml.math.allclose(mat, single_mats)
+
+    @pytest.mark.parametrize(
+        "state_, num_wires",
+        [([1.0, 0.0], 1), ([0.5, -0.5j, 0.5, -0.5], 2), (np.ones(8) / np.sqrt(8), 3)],
+    )
+    def test_controlled_qubit_state_vector(self, state_, num_wires):
+        """Test that QubitStateVector, which is marked as supporting parameter broadcasting,
+        actually does support broadcasting."""
+
+        state = np.array([state_])
+        base = qml.QubitStateVector(state, wires=list(range(num_wires)))
+        op = Controlled(base, "wire1")
+
+        assert op.batch_size == 1
+        qml.QubitStateVector.compute_decomposition(state, list(range(num_wires)))
+        op.decomposition()
+
+        state = np.array([state_] * 3)
+        base = qml.QubitStateVector(state, wires=list(range(num_wires)))
+        op = Controlled(base, "wire1")
+        assert op.batch_size == 3
+        qml.QubitStateVector.compute_decomposition(state, list(range(num_wires)))
+        op.decomposition()
+
+    @pytest.mark.parametrize(
+        "state, num_wires",
+        [([1.0, 0.0], 1), ([0.5, -0.5j, 0.5, -0.5], 2), (np.ones(8) / np.sqrt(8), 3)],
+    )
+    def test_controlled_amplitude_embedding(self, state, num_wires):
+        """Test that AmplitudeEmbedding, which is marked as supporting parameter broadcasting,
+        actually does support broadcasting."""
+
+        features = np.array([state])
+        base = qml.AmplitudeEmbedding(features, wires=list(range(num_wires)))
+        op = Controlled(base, "wire1")
+        assert op.batch_size == 1
+        qml.AmplitudeEmbedding.compute_decomposition(features, list(range(num_wires)))
+        op.decomposition()
+
+        features = np.array([state] * 3)
+        base = qml.AmplitudeEmbedding(features, wires=list(range(num_wires)))
+        op = Controlled(base, "wire1")
+        assert op.batch_size == 3
+        qml.AmplitudeEmbedding.compute_decomposition(features, list(range(num_wires)))
+        op.decomposition()
+
+    @pytest.mark.parametrize(
+        "angles, num_wires",
+        [
+            (np.array([[0.5], [2.1]]), 1),
+            (np.array([[0.5, -0.5], [0.2, 1.5]]), 2),
+            (np.ones((2, 5)), 5),
+        ],
+    )
+    def test_controlled_angle_embedding(self, angles, num_wires):
+        """Test that AngleEmbedding, which is marked as supporting parameter broadcasting,
+        actually does support broadcasting."""
+
+        base = qml.AngleEmbedding(angles, wires=list(range(num_wires)))
+        op = Controlled(base, "wire1")
+        assert op.batch_size == 2
+        qml.AngleEmbedding.compute_decomposition(angles, list(range(num_wires)), rotation=qml.RX)
+        op.decomposition()
+
+    @pytest.mark.parametrize(
+        "features, num_wires",
+        [
+            (np.array([[0.5], [2.1]]), 1),
+            (np.array([[0.5, -0.5], [0.2, 1.5]]), 2),
+            (np.ones((2, 5)), 5),
+        ],
+    )
+    def test_controlled_iqp_embedding(self, features, num_wires):
+        """Test that IQPEmbedding, which is marked as supporting parameter broadcasting,
+        actually does support broadcasting."""
+
+        base = qml.IQPEmbedding(features, wires=list(range(num_wires)))
+        op = Controlled(base, "wire1")
+        assert op.batch_size == 2
+        qml.IQPEmbedding.compute_decomposition(
+            features,
+            list(range(num_wires)),
+            n_repeats=2,
+            pattern=op.base.hyperparameters["pattern"],
+        )
+        op.decomposition()
+
+    @pytest.mark.parametrize(
+        "features, weights, num_wires, batch_size",
+        [
+            (np.array([[0.5], [2.1]]), np.array([[0.61], [0.3]]), 1, 2),
+            (np.array([[0.5, -0.5], [0.2, 1.5]]), np.ones((2, 4, 3)), 2, 2),
+            (np.array([0.5, -0.5, 0.2]), np.ones((3, 2, 6)), 3, 3),
+        ],
+    )
+    def test_controlled_qaoa_embedding(self, features, weights, num_wires, batch_size):
+        """Test that QAOAEmbedding, which is marked as supporting parameter broadcasting,
+        actually does support broadcasting."""
+
+        base = qml.QAOAEmbedding(features, weights, wires=list(range(num_wires)))
+        op = Controlled(base, "wire1")
+        assert op.batch_size == batch_size
+        qml.QAOAEmbedding.compute_decomposition(
+            features, weights, wires=list(range(num_wires)), local_field=qml.RY
+        )
+        op.decomposition()
+
+
 ##### TESTS FOR THE ctrl TRANSFORM #####
 
 
@@ -1036,7 +1351,7 @@ def test_ctrl_sanity_check():
         qml.RX(0.123, wires=0)
         qml.RY(0.456, wires=2)
         qml.RX(0.789, wires=0)
-        qml.Rot(0.111, 0.222, 0.333, wires=2),
+        qml.Rot(0.111, 0.222, 0.333, wires=2)
         qml.PauliX(wires=2)
         qml.PauliY(wires=4)
         qml.PauliZ(wires=0)
@@ -1285,10 +1600,9 @@ class TestCtrlTransformDifferentiation:
     @pytest.mark.autograd
     def test_autograd(self, diff_method):
         """Test differentiation using autograd"""
-        from pennylane import numpy as pnp
 
         dev = qml.device("default.qubit", wires=2)
-        init_state = pnp.array([1.0, -1.0], requires_grad=False) / np.sqrt(2)
+        init_state = np.array([1.0, -1.0], requires_grad=False) / np.sqrt(2)
 
         @qml.qnode(dev, diff_method=diff_method)
         def circuit(b):
@@ -1296,7 +1610,7 @@ class TestCtrlTransformDifferentiation:
             qml.ctrl(qml.RY, control=0)(b, wires=[1])
             return qml.expval(qml.PauliX(0))
 
-        b = pnp.array(0.123, requires_grad=True)
+        b = np.array(0.123, requires_grad=True)
         res = qml.grad(circuit)(b)
         expected = np.sin(b / 2) / 2
 
@@ -1381,7 +1695,7 @@ def test_ctrl_values_sanity_check():
         qml.RX(0.123, wires=0)
         qml.RY(0.456, wires=2)
         qml.RX(0.789, wires=0)
-        qml.Rot(0.111, 0.222, 0.333, wires=2),
+        qml.Rot(0.111, 0.222, 0.333, wires=2)
         qml.PauliX(wires=2)
         qml.PauliY(wires=4)
         qml.PauliZ(wires=0)
diff --git a/tests/ops/op_math/test_controlled_ops.py b/tests/ops/op_math/test_controlled_ops.py
new file mode 100644
index 000000000..c4801b7cd
--- /dev/null
+++ b/tests/ops/op_math/test_controlled_ops.py
@@ -0,0 +1,433 @@
+# Copyright 2018-2020 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Unit tests for Operators inheriting from Controlled.
+"""
+
+import numpy as np
+import pytest
+from scipy.stats import unitary_group
+
+import pennylane as qml
+from pennylane.wires import Wires
+from pennylane.ops.qubit.matrix_ops import QubitUnitary
+
+
+X = np.array([[0, 1], [1, 0]])
+X_broadcasted = np.array([X] * 3)
+
+
+class TestControlledQubitUnitary:
+    """Tests specific to the ControlledQubitUnitary operation"""
+
+    def test_initialization_from_matrix_and_operator(self):
+        base_op = QubitUnitary(X, wires=1)
+
+        op1 = qml.ControlledQubitUnitary(X, control_wires=[0, 2], wires=1)
+        op2 = qml.ControlledQubitUnitary(base_op, control_wires=[0, 2])
+
+        assert qml.equal(op1, op2)
+
+    def test_no_control(self):
+        """Test if ControlledQubitUnitary raises an error if control wires are not specified"""
+        with pytest.raises(
+            TypeError, match="missing 1 required positional argument: 'control_wires'"
+        ):
+            qml.ControlledQubitUnitary(X, wires=2)
+
+    def test_shared_control(self):
+        """Test if ControlledQubitUnitary raises an error if control wires are shared with wires"""
+        with pytest.raises(
+            ValueError, match="The control wires must be different from the base operation wires"
+        ):
+            qml.ControlledQubitUnitary(X, control_wires=[0, 2], wires=2)
+
+    def test_wires_specified_twice_warning(self):
+        base = qml.QubitUnitary(X, 0)
+        with pytest.warns(
+            UserWarning,
+            match="base operator already has wires; values specified through wires kwarg will be ignored.",
+        ):
+            qml.ControlledQubitUnitary(base, control_wires=[1, 2], wires=3)
+
+    def test_wrong_shape(self):
+        """Test if ControlledQubitUnitary raises a ValueError if a unitary of shape inconsistent
+        with wires is provided"""
+        with pytest.raises(ValueError, match=r"Input unitary must be of shape \(2, 2\)"):
+            qml.ControlledQubitUnitary(np.eye(4), control_wires=[0, 1], wires=2).matrix()
+
+    @pytest.mark.parametrize("target_wire", range(3))
+    def test_toffoli(self, target_wire):
+        """Test if ControlledQubitUnitary acts like a Toffoli gate when the input unitary is a
+        single-qubit X. This test allows the target wire to be any of the three wires."""
+        control_wires = list(range(3))
+        del control_wires[target_wire]
+
+        # pick some random unitaries (with a fixed seed) to make the circuit less trivial
+        U1 = unitary_group.rvs(8, random_state=1)
+        U2 = unitary_group.rvs(8, random_state=2)
+
+        dev = qml.device("default.qubit", wires=3)
+
+        @qml.qnode(dev)
+        def f1():
+            qml.QubitUnitary(U1, wires=range(3))
+            qml.ControlledQubitUnitary(X, control_wires=control_wires, wires=target_wire)
+            qml.QubitUnitary(U2, wires=range(3))
+            return qml.state()
+
+        @qml.qnode(dev)
+        def f2():
+            qml.QubitUnitary(U1, wires=range(3))
+            qml.Toffoli(wires=control_wires + [target_wire])
+            qml.QubitUnitary(U2, wires=range(3))
+            return qml.state()
+
+        state_1 = f1()
+        state_2 = f2()
+
+        assert np.allclose(state_1, state_2)
+
+    @pytest.mark.parametrize("target_wire", range(3))
+    def test_toffoli_broadcasted(self, target_wire):
+        """Test if ControlledQubitUnitary acts like a Toffoli gate when the input unitary is a
+        broadcasted single-qubit X. Allows the target wire to be any of the three wires."""
+        control_wires = list(range(3))
+        del control_wires[target_wire]
+
+        # pick some random unitaries (with a fixed seed) to make the circuit less trivial
+        U1 = unitary_group.rvs(8, random_state=1)
+        U2 = unitary_group.rvs(8, random_state=2)
+
+        dev = qml.device("default.qubit", wires=3)
+
+        @qml.qnode(dev)
+        def f1():
+            qml.QubitUnitary(U1, wires=range(3))
+            qml.ControlledQubitUnitary(
+                X_broadcasted, control_wires=control_wires, wires=target_wire
+            )
+            qml.QubitUnitary(U2, wires=range(3))
+            return qml.state()
+
+        @qml.qnode(dev)
+        def f2():
+            qml.QubitUnitary(U1, wires=range(3))
+            qml.Toffoli(wires=control_wires + [target_wire])
+            qml.QubitUnitary(U2, wires=range(3))
+            return qml.state()
+
+        state_1 = f1()
+        state_2 = f2()
+
+        assert np.shape(state_1) == (3, 8)
+        assert np.allclose(state_1, state_1[0])  # Check that all broadcasted results are equal
+        assert np.allclose(state_1, state_2)
+
+    def test_arbitrary_multiqubit(self):
+        """Test if ControlledQubitUnitary applies correctly for a 2-qubit unitary with 2-qubit
+        control, where the control and target wires are not ordered."""
+        control_wires = [1, 3]
+        target_wires = [2, 0]
+
+        # pick some random unitaries (with a fixed seed) to make the circuit less trivial
+        U1 = unitary_group.rvs(16, random_state=1)
+        U2 = unitary_group.rvs(16, random_state=2)
+
+        # the two-qubit unitary
+        U = unitary_group.rvs(4, random_state=3)
+
+        # the 4-qubit representation of the unitary if the control wires were [0, 1] and the target
+        # wires were [2, 3]
+        U_matrix = np.eye(16, dtype=np.complex128)
+        U_matrix[12:16, 12:16] = U
+
+        # We now need to swap wires so that the control wires are [1, 3] and the target wires are
+        # [2, 0]
+        swap = qml.SWAP.compute_matrix()
+
+        # initial wire permutation: 0123
+        # target wire permutation: 1302
+        swap1 = np.kron(swap, np.eye(4))  # -> 1023
+        swap2 = np.kron(np.eye(4), swap)  # -> 1032
+        swap3 = np.kron(np.kron(np.eye(2), swap), np.eye(2))  # -> 1302
+        swap4 = np.kron(np.eye(4), swap)  # -> 1320
+
+        all_swap = swap4 @ swap3 @ swap2 @ swap1
+        U_matrix = all_swap.T @ U_matrix @ all_swap
+
+        dev = qml.device("default.qubit", wires=4)
+
+        @qml.qnode(dev)
+        def f1():
+            qml.QubitUnitary(U1, wires=range(4))
+            qml.ControlledQubitUnitary(U, control_wires=control_wires, wires=target_wires)
+            qml.QubitUnitary(U2, wires=range(4))
+            return qml.state()
+
+        @qml.qnode(dev)
+        def f2():
+            qml.QubitUnitary(U1, wires=range(4))
+            qml.QubitUnitary(U_matrix, wires=range(4))
+            qml.QubitUnitary(U2, wires=range(4))
+            return qml.state()
+
+        state_1 = f1()
+        state_2 = f2()
+
+        assert np.allclose(state_1, state_2)
+
+    def test_mismatched_control_value_length(self):
+        """Test if ControlledQubitUnitary properly handles invalid mixed-polarity
+        control values."""
+        control_wires = [0, 1]
+        wires = 2
+        control_values = "011"
+        target_wires = Wires(wires)
+
+        with pytest.raises(
+            ValueError, match="control_values should be the same length as control_wires"
+        ):
+            qml.ControlledQubitUnitary(
+                X, control_wires=control_wires, wires=target_wires, control_values=control_values
+            )
+
+    @pytest.mark.parametrize(
+        "control_wires,wires,control_values",
+        [
+            ([0], 1, "0"),
+            ([0, 1], 2, "00"),
+            ([0, 1], 2, "10"),
+            ([0, 1], 2, "11"),
+            ([1, 0], 2, "01"),
+            ([0, 1], [2, 3], "11"),
+            ([0, 2], [3, 1], "10"),
+            ([1, 2, 0], [3, 4], "100"),
+            ([1, 0, 2], [4, 3], "110"),
+        ],
+    )
+    def test_mixed_polarity_controls(self, control_wires, wires, control_values):
+        """Test if ControlledQubitUnitary properly applies mixed-polarity
+        control values."""
+        target_wires = Wires(wires)
+
+        dev = qml.device("default.qubit", wires=len(control_wires + target_wires))
+
+        # Pick a random unitary
+        U = unitary_group.rvs(2 ** len(target_wires), random_state=1967)
+
+        # Pick random starting state for the control and target qubits
+        control_state_weights = np.random.normal(size=(2 ** (len(control_wires) + 1) - 2))
+        target_state_weights = np.random.normal(size=(2 ** (len(target_wires) + 1) - 2))
+
+        @qml.qnode(dev)
+        def circuit_mixed_polarity():
+            qml.templates.ArbitraryStatePreparation(control_state_weights, wires=control_wires)
+            qml.templates.ArbitraryStatePreparation(target_state_weights, wires=target_wires)
+
+            qml.ControlledQubitUnitary(
+                U, control_wires=control_wires, wires=target_wires, control_values=control_values
+            )
+            return qml.state()
+
+        # The result of applying the mixed-polarity gate should be the same as
+        # if we conjugated the specified control wires with Pauli X and applied the
+        # "regular" ControlledQubitUnitary in between.
+
+        x_locations = [x for x in range(len(control_values)) if control_values[x] == "0"]
+
+        @qml.qnode(dev)
+        def circuit_pauli_x():
+            qml.templates.ArbitraryStatePreparation(control_state_weights, wires=control_wires)
+            qml.templates.ArbitraryStatePreparation(target_state_weights, wires=target_wires)
+
+            for wire in x_locations:
+                qml.PauliX(wires=control_wires[wire])
+
+            qml.ControlledQubitUnitary(U, control_wires=control_wires, wires=wires)
+
+            for wire in x_locations:
+                qml.PauliX(wires=control_wires[wire])
+
+            return qml.state()
+
+        mixed_polarity_state = circuit_mixed_polarity()
+        pauli_x_state = circuit_pauli_x()
+
+        assert np.allclose(mixed_polarity_state, pauli_x_state)
+
+    def test_same_as_Toffoli(self):
+        """Test if ControlledQubitUnitary returns the correct matrix for a control-control-X
+        (Toffoli) gate"""
+        mat = qml.ControlledQubitUnitary(X, control_wires=[0, 1], wires=2).matrix()
+        mat2 = qml.Toffoli(wires=[0, 1, 2]).matrix()
+        assert np.allclose(mat, mat2)
+
+    def test_matrix_representation(self, tol):
+        """Test that the matrix representation is defined correctly"""
+        U = np.array([[0.94877869, 0.31594146], [-0.31594146, 0.94877869]])
+        res_dynamic = qml.ControlledQubitUnitary(U, control_wires=[1], wires=0).matrix()
+        expected = np.array(
+            [
+                [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+                [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+                [0.0 + 0.0j, 0.0 + 0.0j, 0.94877869 + 0.0j, 0.31594146 + 0.0j],
+                [0.0 + 0.0j, 0.0 + 0.0j, -0.31594146 + 0.0j, 0.94877869 + 0.0j],
+            ]
+        )
+        assert np.allclose(res_dynamic, expected, atol=tol)
+
+    def test_matrix_representation_broadcasted(self, tol):
+        """Test that the matrix representation is defined correctly"""
+        U = np.array(
+            [
+                [[0.94877869, 0.31594146], [-0.31594146, 0.94877869]],
+                [[0.4125124, -0.91095199], [0.91095199, 0.4125124]],
+                [[0.31594146, 0.94877869j], [0.94877869j, 0.31594146]],
+            ]
+        )
+
+        res_dynamic = qml.ControlledQubitUnitary(U, control_wires=[1], wires=0).matrix()
+        expected = np.array(
+            [
+                [
+                    [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+                    [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+                    [0.0 + 0.0j, 0.0 + 0.0j, 0.94877869 + 0.0j, 0.31594146 + 0.0j],
+                    [0.0 + 0.0j, 0.0 + 0.0j, -0.31594146 + 0.0j, 0.94877869 + 0.0j],
+                ],
+                [
+                    [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+                    [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+                    [0.0 + 0.0j, 0.0 + 0.0j, 0.4125124 + 0.0j, -0.91095199 + 0.0j],
+                    [0.0 + 0.0j, 0.0 + 0.0j, 0.91095199 + 0.0j, 0.4125124 + 0.0j],
+                ],
+                [
+                    [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+                    [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+                    [0.0 + 0.0j, 0.0 + 0.0j, 0.31594146 + 0.0j, 0.0 + 0.94877869j],
+                    [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.94877869j, 0.31594146 + 0.0j],
+                ],
+            ]
+        )
+        assert np.allclose(res_dynamic, expected, atol=tol)
+
+    @pytest.mark.parametrize("n", (2, -1, -2))
+    def test_pow(self, n):
+        """Tests the metadata and unitary for a ControlledQubitUnitary raised to a power."""
+        U1 = np.array(
+            [
+                [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
+                [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
+            ]
+        )
+
+        op = qml.ControlledQubitUnitary(
+            U1, control_wires=("b", "c"), wires="a", control_values="01"
+        )
+
+        pow_ops = op.pow(n)
+        assert len(pow_ops) == 1
+
+        assert pow_ops[0].target_wires == op.target_wires
+        assert pow_ops[0].control_wires == op.control_wires
+        assert pow_ops[0].control_values == op.control_values
+
+        op_mat_to_pow = qml.math.linalg.matrix_power(op.data[0], n)
+        assert qml.math.allclose(pow_ops[0].data[0], op_mat_to_pow)
+
+    @pytest.mark.parametrize("n", (2, -1, -2))
+    def test_pow_broadcasted(self, n):
+        """Tests the metadata and unitary for a broadcasted
+        ControlledQubitUnitary raised to a power."""
+        U1 = np.tensordot(
+            np.array([1j, -1.0, 1j]),
+            np.array(
+                [
+                    [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
+                    [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
+                ]
+            ),
+            axes=0,
+        )
+
+        op = qml.ControlledQubitUnitary(U1, control_wires=("b", "c"), wires="a")
+
+        pow_ops = op.pow(n)
+        assert len(pow_ops) == 1
+
+        assert pow_ops[0].target_wires == op.target_wires
+        assert pow_ops[0].control_wires == op.control_wires
+
+        op_mat_to_pow = qml.math.linalg.matrix_power(op.data[0], n)
+        assert qml.math.allclose(pow_ops[0].data[0], op_mat_to_pow)
+
+    def test_noninteger_pow(self):
+        """Test that a ControlledQubitUnitary raised to a non-integer power raises an error."""
+        U1 = np.array(
+            [
+                [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
+                [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
+            ]
+        )
+
+        op = qml.ControlledQubitUnitary(U1, control_wires=("b", "c"), wires="a")
+
+        with pytest.raises(qml.operation.PowUndefinedError):
+            op.pow(0.12)
+
+    def test_noninteger_pow_broadcasted(self):
+        """Test that a ControlledQubitUnitary raised to a non-integer power raises an error."""
+        U1 = np.tensordot(
+            np.array([1j, -1.0, 1j]),
+            np.array(
+                [
+                    [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
+                    [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
+                ]
+            ),
+            axes=0,
+        )
+
+        op = qml.ControlledQubitUnitary(U1, control_wires=("b", "c"), wires="a")
+
+        with pytest.raises(qml.operation.PowUndefinedError):
+            op.pow(0.12)
+
+    @pytest.mark.parametrize("inverse", (True, False))
+    def test_controlled(self, inverse):
+        """Test the _controlled method for ControlledQubitUnitary."""
+
+        U = qml.PauliX(0).compute_matrix()
+
+        original = qml.ControlledQubitUnitary(U, control_wires=(0, 1), wires=4, control_values="01")
+        original.inverse = inverse
+        expected = qml.ControlledQubitUnitary(
+            U, control_wires=(0, 1, "a"), wires=4, control_values="011"
+        )
+        expected.inverse = inverse
+
+        out = original._controlled("a")  # pylint: disable=protected-access
+        assert qml.equal(out, expected)
+
+    def test_unitary_check(self):
+        unitary = np.array([[0.94877869j, 0.31594146], [-0.31594146, 0.94877869j]])
+        not_unitary = np.array([[0.94877869j, 0.31594146], [-5, 0.94877869j]])
+
+        qml.ControlledQubitUnitary(unitary, control_wires=[0, 2], wires=1, unitary_check=True)
+
+        with pytest.warns(UserWarning, match="may not be unitary"):
+            qml.ControlledQubitUnitary(
+                not_unitary, control_wires=[0, 2], wires=1, unitary_check=True
+            )
diff --git a/tests/ops/qubit/test_attributes.py b/tests/ops/qubit/test_attributes.py
index f99825ee0..2c3520239 100644
--- a/tests/ops/qubit/test_attributes.py
+++ b/tests/ops/qubit/test_attributes.py
@@ -21,6 +21,7 @@ from scipy.stats import unitary_group
 import pennylane as qml
 
 from pennylane.ops.qubit.attributes import Attribute
+from pennylane.ops import Controlled
 
 # Dummy attribute
 new_attribute = Attribute(["PauliX", "PauliY", "PauliZ", "Hadamard", "RZ"])
@@ -286,14 +287,12 @@ class TestSupportsBroadcasting:
         op = qml.ControlledQubitUnitary(U, wires=wires, control_wires=[1, "10"])
 
         mat1 = op.matrix()
-        mat2 = qml.ControlledQubitUnitary.compute_matrix(U, u_wires=wires, control_wires=[1, "10"])
         single_mats = [
             qml.ControlledQubitUnitary(_U, wires=wires, control_wires=[1, "10"]).matrix()
             for _U in U
         ]
 
         assert qml.math.allclose(mat1, single_mats)
-        assert qml.math.allclose(mat2, single_mats)
 
     def test_diagonal_qubit_unitary(self):
         """Test that DiagonalQubitUnitary, which is marked as supporting parameter broadcasting,
diff --git a/tests/ops/qubit/test_matrix_ops.py b/tests/ops/qubit/test_matrix_ops.py
index 5a0b3a512..4894fb6fe 100644
--- a/tests/ops/qubit/test_matrix_ops.py
+++ b/tests/ops/qubit/test_matrix_ops.py
@@ -519,399 +519,6 @@ class TestDiagonalQubitUnitary:
         assert np.allclose(grad, expected)
 
 
-X = np.array([[0, 1], [1, 0]])
-X_broadcasted = np.array([X] * 3)
-
-
-class TestControlledQubitUnitary:
-    """Tests for the ControlledQubitUnitary operation"""
-
-    def test_no_control(self):
-        """Test if ControlledQubitUnitary raises an error if control wires are not specified"""
-        with pytest.raises(ValueError, match="Must specify control wires"):
-            qml.ControlledQubitUnitary(X, wires=2)
-
-    def test_shared_control(self):
-        """Test if ControlledQubitUnitary raises an error if control wires are shared with wires"""
-        with pytest.raises(ValueError, match="The control wires must be different from the wires"):
-            qml.ControlledQubitUnitary(X, control_wires=[0, 2], wires=2)
-
-    def test_wrong_shape(self):
-        """Test if ControlledQubitUnitary raises a ValueError if a unitary of shape inconsistent
-        with wires is provided"""
-        with pytest.raises(ValueError, match=r"Input unitary must be of shape \(2, 2\)"):
-            qml.ControlledQubitUnitary(np.eye(4), control_wires=[0, 1], wires=2).matrix()
-
-    @pytest.mark.parametrize("target_wire", range(3))
-    def test_toffoli(self, target_wire):
-        """Test if ControlledQubitUnitary acts like a Toffoli gate when the input unitary is a
-        single-qubit X. This test allows the target wire to be any of the three wires."""
-        control_wires = list(range(3))
-        del control_wires[target_wire]
-
-        # pick some random unitaries (with a fixed seed) to make the circuit less trivial
-        U1 = unitary_group.rvs(8, random_state=1)
-        U2 = unitary_group.rvs(8, random_state=2)
-
-        dev = qml.device("default.qubit", wires=3)
-
-        @qml.qnode(dev)
-        def f1():
-            qml.QubitUnitary(U1, wires=range(3))
-            qml.ControlledQubitUnitary(X, control_wires=control_wires, wires=target_wire)
-            qml.QubitUnitary(U2, wires=range(3))
-            return qml.state()
-
-        @qml.qnode(dev)
-        def f2():
-            qml.QubitUnitary(U1, wires=range(3))
-            qml.Toffoli(wires=control_wires + [target_wire])
-            qml.QubitUnitary(U2, wires=range(3))
-            return qml.state()
-
-        state_1 = f1()
-        state_2 = f2()
-
-        assert np.allclose(state_1, state_2)
-
-    @pytest.mark.parametrize("target_wire", range(3))
-    def test_toffoli_broadcasted(self, target_wire):
-        """Test if ControlledQubitUnitary acts like a Toffoli gate when the input unitary is a
-        broadcasted single-qubit X. Allows the target wire to be any of the three wires."""
-        control_wires = list(range(3))
-        del control_wires[target_wire]
-
-        # pick some random unitaries (with a fixed seed) to make the circuit less trivial
-        U1 = unitary_group.rvs(8, random_state=1)
-        U2 = unitary_group.rvs(8, random_state=2)
-
-        dev = qml.device("default.qubit", wires=3)
-
-        @qml.qnode(dev)
-        def f1():
-            qml.QubitUnitary(U1, wires=range(3))
-            qml.ControlledQubitUnitary(
-                X_broadcasted, control_wires=control_wires, wires=target_wire
-            )
-            qml.QubitUnitary(U2, wires=range(3))
-            return qml.state()
-
-        @qml.qnode(dev)
-        def f2():
-            qml.QubitUnitary(U1, wires=range(3))
-            qml.Toffoli(wires=control_wires + [target_wire])
-            qml.QubitUnitary(U2, wires=range(3))
-            return qml.state()
-
-        state_1 = f1()
-        state_2 = f2()
-
-        assert np.shape(state_1) == (3, 8)
-        assert np.allclose(state_1, state_1[0])  # Check that all broadcasted results are equal
-        assert np.allclose(state_1, state_2)
-
-    def test_arbitrary_multiqubit(self):
-        """Test if ControlledQubitUnitary applies correctly for a 2-qubit unitary with 2-qubit
-        control, where the control and target wires are not ordered."""
-        control_wires = [1, 3]
-        target_wires = [2, 0]
-
-        # pick some random unitaries (with a fixed seed) to make the circuit less trivial
-        U1 = unitary_group.rvs(16, random_state=1)
-        U2 = unitary_group.rvs(16, random_state=2)
-
-        # the two-qubit unitary
-        U = unitary_group.rvs(4, random_state=3)
-
-        # the 4-qubit representation of the unitary if the control wires were [0, 1] and the target
-        # wires were [2, 3]
-        U_matrix = np.eye(16, dtype=np.complex128)
-        U_matrix[12:16, 12:16] = U
-
-        # We now need to swap wires so that the control wires are [1, 3] and the target wires are
-        # [2, 0]
-        swap = qml.SWAP.compute_matrix()
-
-        # initial wire permutation: 0123
-        # target wire permutation: 1302
-        swap1 = np.kron(swap, np.eye(4))  # -> 1023
-        swap2 = np.kron(np.eye(4), swap)  # -> 1032
-        swap3 = np.kron(np.kron(np.eye(2), swap), np.eye(2))  # -> 1302
-        swap4 = np.kron(np.eye(4), swap)  # -> 1320
-
-        all_swap = swap4 @ swap3 @ swap2 @ swap1
-        U_matrix = all_swap.T @ U_matrix @ all_swap
-
-        dev = qml.device("default.qubit", wires=4)
-
-        @qml.qnode(dev)
-        def f1():
-            qml.QubitUnitary(U1, wires=range(4))
-            qml.ControlledQubitUnitary(U, control_wires=control_wires, wires=target_wires)
-            qml.QubitUnitary(U2, wires=range(4))
-            return qml.state()
-
-        @qml.qnode(dev)
-        def f2():
-            qml.QubitUnitary(U1, wires=range(4))
-            qml.QubitUnitary(U_matrix, wires=range(4))
-            qml.QubitUnitary(U2, wires=range(4))
-            return qml.state()
-
-        state_1 = f1()
-        state_2 = f2()
-
-        assert np.allclose(state_1, state_2)
-
-    @pytest.mark.parametrize(
-        "control_wires,wires,control_values,expected_error_message",
-        [
-            ([0, 1], 2, "ab", "String of control values can contain only '0' or '1'."),
-            ([0, 1], 2, "011", "Length of control bit string must equal number of control wires."),
-            ([0, 1], 2, [0, 1], "Alternative control values must be passed as a binary string."),
-        ],
-    )
-    def test_invalid_mixed_polarity_controls(
-        self, control_wires, wires, control_values, expected_error_message
-    ):
-        """Test if ControlledQubitUnitary properly handles invalid mixed-polarity
-        control values."""
-        target_wires = Wires(wires)
-
-        with pytest.raises(ValueError, match=expected_error_message):
-            qml.ControlledQubitUnitary(
-                X, control_wires=control_wires, wires=target_wires, control_values=control_values
-            ).matrix()
-
-    @pytest.mark.parametrize(
-        "control_wires,wires,control_values",
-        [
-            ([0], 1, "0"),
-            ([0, 1], 2, "00"),
-            ([0, 1], 2, "10"),
-            ([0, 1], 2, "11"),
-            ([1, 0], 2, "01"),
-            ([0, 1], [2, 3], "11"),
-            ([0, 2], [3, 1], "10"),
-            ([1, 2, 0], [3, 4], "100"),
-            ([1, 0, 2], [4, 3], "110"),
-        ],
-    )
-    def test_mixed_polarity_controls(self, control_wires, wires, control_values):
-        """Test if ControlledQubitUnitary properly applies mixed-polarity
-        control values."""
-        target_wires = Wires(wires)
-
-        dev = qml.device("default.qubit", wires=len(control_wires + target_wires))
-
-        # Pick a random unitary
-        U = unitary_group.rvs(2 ** len(target_wires), random_state=1967)
-
-        # Pick random starting state for the control and target qubits
-        control_state_weights = np.random.normal(size=(2 ** (len(control_wires) + 1) - 2))
-        target_state_weights = np.random.normal(size=(2 ** (len(target_wires) + 1) - 2))
-
-        @qml.qnode(dev)
-        def circuit_mixed_polarity():
-            qml.templates.ArbitraryStatePreparation(control_state_weights, wires=control_wires)
-            qml.templates.ArbitraryStatePreparation(target_state_weights, wires=target_wires)
-
-            qml.ControlledQubitUnitary(
-                U, control_wires=control_wires, wires=target_wires, control_values=control_values
-            )
-            return qml.state()
-
-        # The result of applying the mixed-polarity gate should be the same as
-        # if we conjugated the specified control wires with Pauli X and applied the
-        # "regular" ControlledQubitUnitary in between.
-
-        x_locations = [x for x in range(len(control_values)) if control_values[x] == "0"]
-
-        @qml.qnode(dev)
-        def circuit_pauli_x():
-            qml.templates.ArbitraryStatePreparation(control_state_weights, wires=control_wires)
-            qml.templates.ArbitraryStatePreparation(target_state_weights, wires=target_wires)
-
-            for wire in x_locations:
-                qml.PauliX(wires=control_wires[wire])
-
-            qml.ControlledQubitUnitary(U, control_wires=control_wires, wires=wires)
-
-            for wire in x_locations:
-                qml.PauliX(wires=control_wires[wire])
-
-            return qml.state()
-
-        mixed_polarity_state = circuit_mixed_polarity()
-        pauli_x_state = circuit_pauli_x()
-
-        assert np.allclose(mixed_polarity_state, pauli_x_state)
-
-    def test_same_as_Toffoli(self):
-        """Test if ControlledQubitUnitary returns the correct matrix for a control-control-X
-        (Toffoli) gate"""
-        mat = qml.ControlledQubitUnitary(X, control_wires=[0, 1], wires=2).matrix()
-        mat2 = qml.Toffoli(wires=[0, 1, 2]).matrix()
-        assert np.allclose(mat, mat2)
-
-    def test_matrix_representation(self, tol):
-        """Test that the matrix representation is defined correctly"""
-        U = np.array([[0.94877869, 0.31594146], [-0.31594146, 0.94877869]])
-        res_static = qml.ControlledQubitUnitary.compute_matrix(U, control_wires=[1], u_wires=[0])
-        res_dynamic = qml.ControlledQubitUnitary(U, control_wires=[1], wires=0).matrix()
-        expected = np.array(
-            [
-                [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
-                [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
-                [0.0 + 0.0j, 0.0 + 0.0j, 0.94877869 + 0.0j, 0.31594146 + 0.0j],
-                [0.0 + 0.0j, 0.0 + 0.0j, -0.31594146 + 0.0j, 0.94877869 + 0.0j],
-            ]
-        )
-        assert np.allclose(res_static, expected, atol=tol)
-        assert np.allclose(res_dynamic, expected, atol=tol)
-
-    def test_matrix_representation_broadcasted(self, tol):
-        """Test that the matrix representation is defined correctly"""
-        U = np.array(
-            [
-                [[0.94877869, 0.31594146], [-0.31594146, 0.94877869]],
-                [[0.4125124, -0.91095199], [0.91095199, 0.4125124]],
-                [[0.31594146, 0.94877869j], [0.94877869j, 0.31594146]],
-            ]
-        )
-
-        res_static = qml.ControlledQubitUnitary.compute_matrix(U, control_wires=[1], u_wires=[0])
-        res_dynamic = qml.ControlledQubitUnitary(U, control_wires=[1], wires=0).matrix()
-        expected = np.array(
-            [
-                [
-                    [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
-                    [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
-                    [0.0 + 0.0j, 0.0 + 0.0j, 0.94877869 + 0.0j, 0.31594146 + 0.0j],
-                    [0.0 + 0.0j, 0.0 + 0.0j, -0.31594146 + 0.0j, 0.94877869 + 0.0j],
-                ],
-                [
-                    [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
-                    [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
-                    [0.0 + 0.0j, 0.0 + 0.0j, 0.4125124 + 0.0j, -0.91095199 + 0.0j],
-                    [0.0 + 0.0j, 0.0 + 0.0j, 0.91095199 + 0.0j, 0.4125124 + 0.0j],
-                ],
-                [
-                    [1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
-                    [0.0 + 0.0j, 1.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
-                    [0.0 + 0.0j, 0.0 + 0.0j, 0.31594146 + 0.0j, 0.0 + 0.94877869j],
-                    [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.94877869j, 0.31594146 + 0.0j],
-                ],
-            ]
-        )
-        assert np.allclose(res_static, expected, atol=tol)
-        assert np.allclose(res_dynamic, expected, atol=tol)
-
-    def test_no_decomp(self):
-        """Test that ControlledQubitUnitary raises a decomposition undefined
-        error."""
-        mat = qml.PauliX(0).matrix()
-        with pytest.raises(qml.operation.DecompositionUndefinedError):
-            qml.ControlledQubitUnitary(mat, wires=0, control_wires=1).decomposition()
-        with pytest.raises(qml.operation.DecompositionUndefinedError):
-            qml.ControlledQubitUnitary(X_broadcasted, wires=0, control_wires=1).decomposition()
-
-    @pytest.mark.parametrize("n", (2, -1, -2))
-    def test_pow(self, n):
-        """Tests the metadata and unitary for a ControlledQubitUnitary raised to a power."""
-        U1 = np.array(
-            [
-                [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
-                [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
-            ]
-        )
-
-        op = qml.ControlledQubitUnitary(
-            U1, control_wires=("b", "c"), wires="a", control_values="01"
-        )
-
-        pow_ops = op.pow(n)
-        assert len(pow_ops) == 1
-
-        assert pow_ops[0].hyperparameters["u_wires"] == op.hyperparameters["u_wires"]
-        assert pow_ops[0].control_wires == op.control_wires
-        assert pow_ops[0].control_values == op.control_values
-
-        op_mat_to_pow = qml.math.linalg.matrix_power(op.data[0], n)
-        assert qml.math.allclose(pow_ops[0].data[0], op_mat_to_pow)
-
-    @pytest.mark.parametrize("n", (2, -1, -2))
-    def test_pow_broadcasted(self, n):
-        """Tests the metadata and unitary for a broadcasted
-        ControlledQubitUnitary raised to a power."""
-        U1 = np.tensordot(
-            np.array([1j, -1.0, 1j]),
-            np.array(
-                [
-                    [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
-                    [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
-                ]
-            ),
-            axes=0,
-        )
-
-        op = qml.ControlledQubitUnitary(U1, control_wires=("b", "c"), wires="a")
-
-        pow_ops = op.pow(n)
-        assert len(pow_ops) == 1
-
-        assert pow_ops[0].hyperparameters["u_wires"] == op.hyperparameters["u_wires"]
-        assert pow_ops[0].control_wires == op.control_wires
-
-        op_mat_to_pow = qml.math.linalg.matrix_power(op.data[0], n)
-        assert qml.math.allclose(pow_ops[0].data[0], op_mat_to_pow)
-
-    def test_noninteger_pow(self):
-        """Test that a ControlledQubitUnitary raised to a non-integer power raises an error."""
-        U1 = np.array(
-            [
-                [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
-                [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
-            ]
-        )
-
-        op = qml.ControlledQubitUnitary(U1, control_wires=("b", "c"), wires="a")
-
-        with pytest.raises(qml.operation.PowUndefinedError):
-            op.pow(0.12)
-
-    def test_noninteger_pow_broadcasted(self):
-        """Test that a ControlledQubitUnitary raised to a non-integer power raises an error."""
-        U1 = np.array(
-            [
-                [0.73708696 + 0.61324932j, 0.27034258 + 0.08685028j],
-                [-0.24979544 - 0.1350197j, 0.95278437 + 0.1075819j],
-            ]
-            * 3
-        )
-
-        op = qml.ControlledQubitUnitary(U1, control_wires=("b", "c"), wires="a")
-
-        with pytest.raises(qml.operation.PowUndefinedError):
-            op.pow(0.12)
-
-    @pytest.mark.parametrize("inverse", (True, False))
-    def test_controlled(self, inverse):
-        """Test the _controlled method for ControlledQubitUnitary."""
-
-        U = qml.PauliX(0).compute_matrix()
-
-        original = qml.ControlledQubitUnitary(U, control_wires=(0, 1), wires=4, control_values="01")
-        original.inverse = inverse
-        expected = qml.ControlledQubitUnitary(
-            U, control_wires=(0, 1, "a"), wires=4, control_values="011"
-        )
-        expected.inverse = inverse
-
-        out = original._controlled("a")
-        assert qml.equal(out, expected)
-
-
 label_data = [
     (X, qml.QubitUnitary(X, wires=0)),
     (X, qml.ControlledQubitUnitary(X, control_wires=0, wires=1)),
diff --git a/tests/tests_passing_pylint.txt b/tests/tests_passing_pylint.txt
index aea07075b..1dc84ea3b 100644
--- a/tests/tests_passing_pylint.txt
+++ b/tests/tests_passing_pylint.txt
@@ -1,2 +1,4 @@
 tests/data/test_dataset.py
 tests/data/test_dataset_access.py
+tests/ops/op_math/test_controlled_ops.py
+tests/ops/op_math/test_controlled.py
\ No newline at end of file