Repository: PennyLaneAI/pennylane
Pull request #5735 (2024-05-23 21:56:25)

Allow multi controlled x decomposition without working wires
**Context:**
`MultiControlled.decomposition()` gives an error when no `work_wire` is provided.

**Description of the Change:**

Implement decomposition given in Lemma 7.5 of [1]  to decompose a general multi controlled unitary gate. `MultiControlledX`  being a particular case.

The function `ctrl_decomp_zyz` now returns the correct decomposition including a controlled global phase. It now checks that this decomposition only works for single control wire.

A decomposition for single controlled `GlobalPhase` is now available using Lemma 5.2 of [1], i.e. adds a `PhaseShift` to the control wire. Raise an exception when used with more than one control.

**Benefits:**

Decomposition of a MultiControlledX (and in general any multi controlled unitary) is now possible even if work wires are not provided. 

If controlled gate is SU(2), this resorts to the specialized algorithm from Lemma 7.9, yielding a decomposition with $\mathcal{O}(\text{len(control wires)})$ gates. For general unitaries, the decomposition yields $\mathcal{O}(\text{len(control wires)}^2)$ gates.


### Sources:

[1] Barenco et al. (1995). <https://arxiv.org/abs/quant-ph/9503016>

[SC-38954]


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5735?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 99.67%. Comparing base [(`b7b4b75`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/b7b4b7521dd76c7e9bcd7c323e303e6f1bb4a6e0?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`6713e63`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5735?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> Report is 4 commits behind head on master.


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #5735      +/-   ##
==========================================
- Coverage   99.68%   99.67%   -0.01%     
==========================================
  Files         415      416       +1     
  Lines       38898    38748     -150     
==========================================
- Hits        38774    38623     -151     
- Misses        124      125       +1     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5735?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

EmilianoG-byte: I think the code in line:
https://github.com/PennyLaneAI/pennylane/blob/7a216b65b559eaf87c5bb6377f9227c6be504844/pennylane/ops/op_math/controlled.py#L799
should be modified since `ctrl_decomp_zyz` can now be used for general unitaries. Is this necessary? Currently the code returns the controlled version of each of the operators in the `op.base` decomposition. In cases like:

```python
base = qml.SX(0)
op = qml.ops.Controlled(base, control_wires=1)
op.decomposition()
>>> [CRZ(1.5707963267948966, wires=<Wires = [1, 0]>),
 CRY(1.5707963267948966, wires=[1, 0])),
 CRZ(-3.141592653589793, wires=<Wires = [1, 0]>),
 ControlledPhaseShift(1.5707963267948966, wires=<Wires = [1, 0]>)]
```

shouldn't this be further decomposed since we know how to handle the decomposition of each of the controlled operators?

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 86f6889ea..255014a7c 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -174,6 +174,15 @@
 * Empty initialization of `PauliVSpace` is permitted.
   [(#5675)](https://github.com/PennyLaneAI/pennylane/pull/5675)
 
+* `MultiControlledX` can now be decomposed even when no `work_wires` are provided. The implementation returns $\mathcal{O}(\text{len(control\_wires)}^2)$ operations, and is applicable for any multi controlled unitary gate.
+  [(#5735)](https://github.com/PennyLaneAI/pennylane/pull/5735)
+
+* Single control unitary now includes the correct global phase.
+  [(#5735)](https://github.com/PennyLaneAI/pennylane/pull/5735)
+
+* Single control `GlobalPhase` has now a decomposition, i.e. relative phase on control wire.
+  [(#5735)](https://github.com/PennyLaneAI/pennylane/pull/5735)
+
 * `QuantumScript` properties are only calculated when needed, instead of on initialization. This decreases the classical overhead by >20%.
   `par_info`, `obs_sharing_wires`, and `obs_sharing_wires_id` are now public attributes.
   [(#5696)](https://github.com/PennyLaneAI/pennylane/pull/5696)
diff --git a/pennylane/ops/op_math/controlled.py b/pennylane/ops/op_math/controlled.py
index 7d901db0d..8fcc69e6a 100644
--- a/pennylane/ops/op_math/controlled.py
+++ b/pennylane/ops/op_math/controlled.py
@@ -764,6 +764,9 @@ def _decompose_custom_ops(op: Controlled) -> List["operation.Operator"]:
         # has some special case handling of its own for further decomposition
         return _decompose_pauli_x_based_no_control_values(op)
 
+    if isinstance(op.base, qml.GlobalPhase) and len(op.control_wires) == 1:
+        # use Lemma 5.2 from https://arxiv.org/pdf/quant-ph/9503016
+        return [qml.PhaseShift(phi=-op.data[0], wires=op.control_wires)]
     # A multi-wire controlled PhaseShift should be decomposed first using the decomposition
     # of ControlledPhaseShift. This is because the decomposition of PhaseShift contains a
     # GlobalPhase that we do not have a handling for.
@@ -802,9 +805,9 @@ def _decompose_no_control_values(op: Controlled) -> List["operation.Operator"]:
         return None
 
     base_decomp = op.base.decomposition()
-    if len(base_decomp) == 0 and isinstance(op.base, qml.GlobalPhase):
+    if len(base_decomp) == 0 and isinstance(op.base, qml.GlobalPhase) and len(op.control_wires) > 1:
         warnings.warn(
-            "Controlled-GlobalPhase currently decomposes to nothing, and this will likely "
+            "Multi-Controlled-GlobalPhase currently decomposes to nothing, and this will likely "
             "produce incorrect results. Consider implementing your circuit with a different set "
             "of operations, or use a device that natively supports GlobalPhase.",
             UserWarning,
diff --git a/pennylane/ops/op_math/controlled_decompositions.py b/pennylane/ops/op_math/controlled_decompositions.py
index 307bbe71b..a72815a4a 100644
--- a/pennylane/ops/op_math/controlled_decompositions.py
+++ b/pennylane/ops/op_math/controlled_decompositions.py
@@ -16,7 +16,6 @@ This submodule defines functions to decompose controlled operations
 """
 
 from copy import copy
-from typing import Tuple
 
 import numpy as np
 import numpy.linalg as npl
@@ -24,10 +23,19 @@ import numpy.linalg as npl
 import pennylane as qml
 from pennylane import math
 from pennylane.operation import Operation, Operator
+from pennylane.ops.op_math.decompositions.single_qubit_unitary import (
+    _get_single_qubit_rot_angles_via_matrix,
+)
 from pennylane.wires import Wires
 
 
-def _convert_to_su2(U):
+def _is_single_qubit_special_unitary(op):
+    mat = op.matrix()
+    det = mat[0, 0] * mat[1, 1] - mat[0, 1] * mat[1, 0]
+    return qml.math.allclose(det, 1)
+
+
+def _convert_to_su2(U, return_global_phase=False):
     r"""Convert a 2x2 unitary matrix to :math:`SU(2)`.
 
     Args:
@@ -43,11 +51,12 @@ def _convert_to_su2(U):
         :math:`SU(2)` equivalent and the second, the global phase.
     """
     # Compute the determinants
-    dets = math.linalg.det(U)
+    with np.errstate(divide="ignore", invalid="ignore"):
+        dets = math.linalg.det(U)
 
-    exp_angles = math.cast_like(math.angle(dets), 1j) / 2
-    U_SU2 = math.cast_like(U, dets) * math.exp(-1j * exp_angles)
-    return U_SU2
+    global_phase = math.cast_like(math.angle(dets), 1j) / 2
+    U_SU2 = math.cast_like(U, dets) * math.exp(-1j * global_phase)
+    return (U_SU2, global_phase) if return_global_phase else U_SU2
 
 
 def _convert_to_real_diagonal(q: np.ndarray) -> np.ndarray:
@@ -128,16 +137,13 @@ def _bisect_compute_b(u: np.ndarray):
 def ctrl_decomp_zyz(target_operation: Operator, control_wires: Wires):
     """Decompose the controlled version of a target single-qubit operation
 
-    This function decomposes a controlled single-qubit target operation using the
-    decomposition defined in section 5 of
+    This function decomposes a controlled single-qubit target operation with one
+    single control using the decomposition defined in Lemma 4.3 and Lemma 5.1 of
     `Barenco et al. (1995) <https://arxiv.org/abs/quant-ph/9503016>`_.
 
-    .. warning:: This method will add a global phase for target operations that do not
-        belong to the SU(2) group.
-
     Args:
         target_operation (~.operation.Operator): the target operation to decompose
-        control_wires (~.wires.Wires): the control wires of the operation
+        control_wires (~.wires.Wires): the control wires of the operation.
 
     Returns:
         list[Operation]: the decomposed operations
@@ -152,31 +158,28 @@ def ctrl_decomp_zyz(target_operation: Operator, control_wires: Wires):
 
     .. code-block:: python
 
-        dev = qml.device("default.qubit", wires=3)
+        dev = qml.device("default.qubit", wires=2)
 
         @qml.qnode(dev)
         def expected_circuit(op):
             qml.Hadamard(wires=0)
-            qml.Hadamard(wires=1)
-            qml.ctrl(op, [0,1])
+            qml.ctrl(op, [0])
             return qml.probs()
 
         @qml.qnode(dev)
         def decomp_circuit(op):
             qml.Hadamard(wires=0)
-            qml.Hadamard(wires=1)
-            qml.ops.ctrl_decomp_zyz(op, [0,1])
+            qml.ops.ctrl_decomp_zyz(op, [0])
             return qml.probs()
 
     Measurements on both circuits will give us the same results:
 
-    >>> op = qml.RX(0.123, wires=2)
+    >>> op = qml.RX(0.123, wires=1)
     >>> expected_circuit(op)
-    tensor([0.25      , 0.        , 0.25      , 0.        , 0.25      ,
-        0.        , 0.24905563, 0.00094437], requires_grad=True)
+    tensor([0.5       , 0.        , 0.49811126, 0.00188874], requires_grad=True)
+
     >>> decomp_circuit(op)
-    tensor([0.25      , 0.        , 0.25      , 0.        , 0.25      ,
-        0.        , 0.24905563, 0.00094437], requires_grad=True)
+    tensor([0.5       , 0.        , 0.49811126, 0.00188874], requires_grad=True)
 
     """
     if len(target_operation.wires) != 1:
@@ -184,46 +187,50 @@ def ctrl_decomp_zyz(target_operation: Operator, control_wires: Wires):
             "The target operation must be a single-qubit operation, instead "
             f"got {target_operation.__class__.__name__}."
         )
+    control_wires = Wires(control_wires)
+    if len(control_wires) > 1:
+        raise ValueError(
+            f"The control_wires should be a single wire, instead got: {len(control_wires)}-wires"
+        )
 
     target_wire = target_operation.wires
 
-    def get_single_qubit_rot_angles_via_matrix() -> Tuple[float, float, float]:
-        """Returns a triplet of angles representing the single-qubit decomposition
-        of the matrix of the target operation using ZYZ rotations.
-        """
-        with qml.QueuingManager.stop_recording():
-            zyz_decomp = qml.ops.one_qubit_decomposition(
-                qml.matrix(target_operation),
-                wire=target_wire,
-                rotations="ZYZ",
-            )
-        return tuple(gate.parameters[0] for gate in zyz_decomp)  # type: ignore
-
     if isinstance(target_operation, Operation):
         try:
             phi, theta, omega = target_operation.single_qubit_rot_angles()
         except NotImplementedError:
-            phi, theta, omega = get_single_qubit_rot_angles_via_matrix()
+            phi, theta, omega = _get_single_qubit_rot_angles_via_matrix(
+                qml.matrix(target_operation)
+            )
     else:
-        phi, theta, omega = get_single_qubit_rot_angles_via_matrix()
+        phi, theta, omega = _get_single_qubit_rot_angles_via_matrix(qml.matrix(target_operation))
 
-    decomp = []
+    _, global_phase = _convert_to_su2(qml.matrix(target_operation), return_global_phase=True)
 
+    # We use the conditional statements to account when decomposition is ran within a queue
+    decomp = []
+    # Add negative of global phase. Compare definition of qml.GlobalPhase and Ph(delta) from section 4.1 of Barenco et al.
+    if not qml.math.allclose(0.0, global_phase, atol=1e-8, rtol=0):
+        decomp.append(
+            qml.ctrl(qml.GlobalPhase(phi=-global_phase, wires=target_wire), control=control_wires)
+        )
+    # Add A operator
     if not qml.math.allclose(0.0, phi, atol=1e-8, rtol=0):
         decomp.append(qml.RZ(phi, wires=target_wire))
     if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
-        decomp.extend(
-            [
-                qml.RY(theta / 2, wires=target_wire),
-                qml.ctrl(qml.X(target_wire), control=control_wires),
-                qml.RY(-theta / 2, wires=target_wire),
-            ]
-        )
-    else:
-        decomp.append(qml.ctrl(qml.X(target_wire), control=control_wires))
+        decomp.append(qml.RY(theta / 2, wires=target_wire))
+
+    decomp.append(qml.ctrl(qml.X(target_wire), control=control_wires))
+
+    # Add B operator
+    if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.RY(-theta / 2, wires=target_wire))
     if not qml.math.allclose(0.0, -(phi + omega) / 2, atol=1e-6, rtol=0):
         decomp.append(qml.RZ(-(phi + omega) / 2, wires=target_wire))
+
     decomp.append(qml.ctrl(qml.PauliX(wires=target_wire), control=control_wires))
+
+    # Add C operator
     if not qml.math.allclose(0.0, (omega - phi) / 2, atol=1e-8, rtol=0):
         decomp.append(qml.RZ((omega - phi) / 2, wires=target_wire))
 
@@ -458,14 +465,91 @@ def ctrl_decomp_bisect(
 
 
 def decompose_mcx(control_wires, target_wire, work_wires):
-    """Decomposes the multi-controlled PauliX gate"""
-
+    """Decomposes the multi-controlled PauliX gate using decompositions from
+    `Barenco et al. (1995) <https://arxiv.org/abs/quant-ph/9503016>`_"""
     num_work_wires_needed = len(control_wires) - 2
 
+    if len(control_wires) == 1:
+        return [qml.CNOT(wires=control_wires + Wires(target_wire))]
+    if len(control_wires) == 2:
+        return qml.Toffoli.compute_decomposition(wires=control_wires + Wires(target_wire))
+
     if len(work_wires) >= num_work_wires_needed:
+        # Lemma 7.2
         return _decompose_mcx_with_many_workers(control_wires, target_wire, work_wires)
+    if len(work_wires) >= 1:
+        # Lemma 7.3
+        return _decompose_mcx_with_one_worker(control_wires, target_wire, work_wires[0])
+
+    # Lemma 7.5
+    with qml.QueuingManager.stop_recording():
+        op = qml.X(target_wire)
+    return _decompose_multicontrolled_unitary(op, control_wires)
 
-    return _decompose_mcx_with_one_worker(control_wires, target_wire, work_wires[0])
+
+def _decompose_multicontrolled_unitary(op, control_wires):
+    """Decomposes general multi controlled unitary with no work wires
+    Follows approach from Lemma 7.5 combined with 7.3 and 7.2 of
+    https://arxiv.org/abs/quant-ph/9503016.
+
+    We are assuming this decomposition is used only in the general cases
+    """
+    if not op.has_matrix or len(op.wires) != 1:
+        raise ValueError(
+            "The target operation must be a single-qubit operation with a matrix representation"
+        )
+
+    target_wire = op.wires
+    if len(control_wires) == 0:
+        return [op]
+    if len(control_wires) == 1:
+        return ctrl_decomp_zyz(op, control_wires)
+    if _is_single_qubit_special_unitary(op):
+        return ctrl_decomp_bisect(op, control_wires)
+    # use recursive decomposition of general gate
+    return _decompose_recursive(op, 1.0, control_wires, target_wire, Wires([]))
+
+
+def _decompose_recursive(op, power, control_wires, target_wire, work_wires):
+    """Decompose multicontrolled operator recursively using lemma 7.5
+    Number of gates in decomposition are: O(len(control_wires)^2)
+    """
+    if len(control_wires) == 1:
+        with qml.QueuingManager.stop_recording():
+            powered_op = qml.pow(op, power, lazy=True)
+        return ctrl_decomp_zyz(powered_op, control_wires)
+
+    with qml.QueuingManager.stop_recording():
+        cnots = decompose_mcx(
+            control_wires=control_wires[:-1],
+            target_wire=control_wires[-1],
+            work_wires=work_wires + target_wire,
+        )
+    with qml.QueuingManager.stop_recording():
+        powered_op = qml.pow(op, 0.5 * power, lazy=True)
+        powered_op_adj = qml.adjoint(powered_op, lazy=True)
+
+    if qml.QueuingManager.recording():
+        decomposition = [
+            *ctrl_decomp_zyz(powered_op, control_wires[-1]),
+            *(qml.apply(o) for o in cnots),
+            *ctrl_decomp_zyz(powered_op_adj, control_wires[-1]),
+            *(qml.apply(o) for o in cnots),
+            *_decompose_recursive(
+                op, 0.5 * power, control_wires[:-1], target_wire, control_wires[-1] + work_wires
+            ),
+        ]
+    else:
+        decomposition = [
+            *ctrl_decomp_zyz(powered_op, control_wires[-1]),
+            *cnots,
+            *ctrl_decomp_zyz(powered_op_adj, control_wires[-1]),
+            *cnots,
+            *_decompose_recursive(
+                op, 0.5 * power, control_wires[:-1], target_wire, control_wires[-1] + work_wires
+            ),
+        ]
+    return decomposition
 
 
 def _decompose_mcx_with_many_workers(control_wires, target_wire, work_wires):
diff --git a/pennylane/ops/op_math/controlled_ops.py b/pennylane/ops/op_math/controlled_ops.py
index 34b13535d..c1953ba6d 100644
--- a/pennylane/ops/op_math/controlled_ops.py
+++ b/pennylane/ops/op_math/controlled_ops.py
@@ -1202,19 +1202,10 @@ class MultiControlledX(ControlledOp):
             control_values = [True] * len(control_wires)
 
         work_wires = work_wires or []
-        if len(control_wires) > 2 and len(work_wires) == 0:
-            raise ValueError(
-                "At least one work wire is required to decompose operation: MultiControlledX"
-            )
 
         flips1 = [qml.X(w) for w, val in zip(control_wires, control_values) if not val]
 
-        if len(control_wires) == 1:
-            decomp = [qml.CNOT(wires=wires)]
-        elif len(control_wires) == 2:
-            decomp = qml.Toffoli.compute_decomposition(wires=wires)
-        else:
-            decomp = decompose_mcx(control_wires, target_wire, work_wires)
+        decomp = decompose_mcx(control_wires, target_wire, work_wires)
 
         flips2 = [qml.X(w) for w, val in zip(control_wires, control_values) if not val]
 
diff --git a/pennylane/ops/op_math/decompositions/single_qubit_unitary.py b/pennylane/ops/op_math/decompositions/single_qubit_unitary.py
index 7497e1570..956082b7c 100644
--- a/pennylane/ops/op_math/decompositions/single_qubit_unitary.py
+++ b/pennylane/ops/op_math/decompositions/single_qubit_unitary.py
@@ -15,6 +15,8 @@
 operations into elementary gates.
 """
 
+from typing import Tuple
+
 import numpy as np
 
 import pennylane as qml
@@ -39,7 +41,8 @@ def _convert_to_su2(U, return_global_phase=False):
 
     # Compute the determinants
     U = qml.math.cast(U, "complex128")
-    determinants = math.linalg.det(U)
+    with np.errstate(divide="ignore", invalid="ignore"):
+        determinants = math.linalg.det(U)
     phase = math.angle(determinants) / 2
     U = math.cast_like(U, determinants) * math.exp(-1j * math.cast_like(phase, 1j))[:, None, None]
 
@@ -158,6 +161,26 @@ def _rot_decomposition(U, wire, return_global_phase=False):
     return operations
 
 
+def _get_single_qubit_rot_angles_via_matrix(
+    U, return_global_phase=False
+) -> Tuple[float, float, float]:
+    """Returns a triplet of angles representing the single-qubit decomposition
+    of the matrix of the target operation using ZYZ rotations.
+    """
+    # Cast to batched format for more consistent code
+    U = math.expand_dims(U, axis=0) if len(U.shape) == 2 else U
+
+    # Convert to SU(2) format and extract global phase
+    U_su2, global_phase = _convert_to_su2(U, return_global_phase=True)
+
+    # Compute the zyz rotation angles
+    phis, thetas, omegas = _zyz_get_rotation_angles(U_su2)
+    angles = (phis, thetas, omegas)
+    if return_global_phase:
+        angles += (global_phase,)
+    return angles
+
+
 def _zyz_decomposition(U, wire, return_global_phase=False):
     r"""Compute the decomposition of a single-qubit matrix :math:`U` in terms
     of elementary operations, as a product of Z and Y rotations in the form
@@ -194,20 +217,14 @@ def _zyz_decomposition(U, wire, return_global_phase=False):
      GlobalPhase(1.1759220332464762, wires=[])]
 
     """
-
-    # Cast to batched format for more consistent code
-    U = math.expand_dims(U, axis=0) if len(U.shape) == 2 else U
-
-    # Convert to SU(2) format and extract global phase
-    U_det1, alphas = _convert_to_su2(U, return_global_phase=True)
-
-    # Compute the zyz rotation angles
-    phis, thetas, omegas = _zyz_get_rotation_angles(U_det1)
+    phis, thetas, omegas, *global_phase = _get_single_qubit_rot_angles_via_matrix(
+        U, return_global_phase=True
+    )
 
     operations = [qml.RZ(phis, wire), qml.RY(thetas, wire), qml.RZ(omegas, wire)]
     if return_global_phase:
-        alphas = math.squeeze(alphas)
-        operations.append(qml.GlobalPhase(-alphas))
+        global_phase = math.squeeze(global_phase)
+        operations.append(qml.GlobalPhase(-global_phase))
 
     return operations
 


Test patch:
diff --git a/tests/ops/op_math/test_controlled.py b/tests/ops/op_math/test_controlled.py
index 225bfb76b..831d23074 100644
--- a/tests/ops/op_math/test_controlled.py
+++ b/tests/ops/op_math/test_controlled.py
@@ -1051,16 +1051,15 @@ class TestDecomposition:
         assert qml.math.allclose(op.matrix(), decomp_mat)
 
     def test_differentiable_one_qubit_special_unitary(self):
-        """Assert that a differentiable qubit speical unitary uses the zyz decomposition."""
+        """Assert that a differentiable qubit special unitary uses the zyz decomposition."""
 
-        op = qml.ctrl(qml.RZ(qml.numpy.array(1.2), 0), (1, 2, 3, 4))
+        op = qml.ctrl(qml.RZ(qml.numpy.array(1.2), 0), (1))
         decomp = op.decomposition()
 
-        assert qml.equal(decomp[0], qml.RZ(qml.numpy.array(1.2), 0))
-        assert qml.equal(decomp[1], qml.MultiControlledX(wires=(1, 2, 3, 4, 0)))
-        assert qml.equal(decomp[2], qml.RZ(qml.numpy.array(-0.6), wires=0))
-        assert qml.equal(decomp[3], qml.MultiControlledX(wires=(1, 2, 3, 4, 0)))
-        assert qml.equal(decomp[4], qml.RZ(qml.numpy.array(-0.6), wires=0))
+        assert qml.equal(decomp[0], qml.PhaseShift(qml.numpy.array(1.2 / 2), 0))
+        assert qml.equal(decomp[1], qml.CNOT(wires=(1, 0)))
+        assert qml.equal(decomp[2], qml.PhaseShift(qml.numpy.array(-1.2 / 2), 0))
+        assert qml.equal(decomp[3], qml.CNOT(wires=(1, 0)))
 
         decomp_mat = qml.matrix(op.decomposition, wire_order=op.wires)()
         assert qml.math.allclose(op.matrix(), decomp_mat)
@@ -1153,8 +1152,8 @@ class TestDecomposition:
             op.decomposition()
 
     def test_global_phase_decomp_raises_warning(self):
-        """Test that ctrl(GlobalPhase).decomposition() raises a warning."""
-        op = qml.ctrl(qml.GlobalPhase(1.23), control=[0])
+        """Test that ctrl(GlobalPhase).decomposition() raises a warning with more than one control."""
+        op = qml.ctrl(qml.GlobalPhase(1.23), control=[0, 1])
         with pytest.warns(
             UserWarning, match="Controlled-GlobalPhase currently decomposes to nothing"
         ):
@@ -1731,6 +1730,7 @@ class TestCtrl:
 
         if isinstance(op, qml.QubitUnitary):
             pytest.skip("ControlledQubitUnitary can accept any number of control wires.")
+            expected = None  # to pass pylint(possibly-used-before-assignment) error
         elif isinstance(op, Controlled):
             expected = Controlled(
                 op.base,
diff --git a/tests/ops/op_math/test_controlled_decompositions.py b/tests/ops/op_math/test_controlled_decompositions.py
index b13017bea..f80c45082 100644
--- a/tests/ops/op_math/test_controlled_decompositions.py
+++ b/tests/ops/op_math/test_controlled_decompositions.py
@@ -23,7 +23,7 @@ import pytest
 import pennylane as qml
 from pennylane import math
 from pennylane.ops import ctrl_decomp_zyz
-from pennylane.ops.op_math.controlled import Controlled
+from pennylane.ops.op_math.controlled import Controlled, _is_single_qubit_special_unitary
 from pennylane.ops.op_math.controlled_decompositions import (
     _bisect_compute_a,
     _bisect_compute_b,
@@ -33,7 +33,10 @@ from pennylane.ops.op_math.controlled_decompositions import (
     _ctrl_decomp_bisect_od,
     _decompose_mcx_with_many_workers,
     _decompose_mcx_with_one_worker,
+    _decompose_multicontrolled_unitary,
+    _decompose_recursive,
     ctrl_decomp_bisect,
+    decompose_mcx,
 )
 from pennylane.wires import Wires
 
@@ -84,6 +87,11 @@ class TestControlledDecompositionZYZ:
         ):
             _ = ctrl_decomp_zyz(qml.CNOT([0, 1]), [2])
 
+    def test_invalid_num_controls(self):
+        """Tests that an error is raised when an invalid number of control wires is passed"""
+        with pytest.raises(ValueError, match="The control_wires should be a single wire"):
+            _ = ctrl_decomp_zyz(qml.X([1]), [0, 1])
+
     su2_ops = [
         qml.RX(0.123, wires=0),
         qml.RY(0.123, wires=0),
@@ -109,7 +117,7 @@ class TestControlledDecompositionZYZ:
     ]
 
     @pytest.mark.parametrize("op", su2_ops + unitary_ops)
-    @pytest.mark.parametrize("control_wires", ([1], [1, 2], [1, 2, 3]))
+    @pytest.mark.parametrize("control_wires", ([1], [2], [3]))
     def test_decomposition_circuit(self, op, control_wires, tol):
         """Tests that the controlled decomposition of a single-qubit operation
         behaves as expected in a quantum circuit"""
@@ -132,7 +140,7 @@ class TestControlledDecompositionZYZ:
 
         assert np.allclose(res, expected, atol=tol, rtol=0)
 
-    @pytest.mark.parametrize("control_wires", ([1], [1, 2], [1, 2, 3]))
+    @pytest.mark.parametrize("control_wires", ([1], [2], [3]))
     def test_decomposition_circuit_gradient(self, control_wires, tol):
         """Tests that the controlled decomposition of a single-qubit operation
         behaves as expected in a quantum circuit"""
@@ -168,7 +176,7 @@ class TestControlledDecompositionZYZ:
         assert np.allclose(jac_ad, jac_bp, atol=tol, rtol=0)
 
     @pytest.mark.parametrize("op", su2_ops)
-    @pytest.mark.parametrize("control_wires", ([1], [1, 2], [1, 2, 3]))
+    @pytest.mark.parametrize("control_wires", ([1], [2], [3]))
     def test_decomposition_matrix(self, op, control_wires, tol):
         """Tests that the matrix representation of the controlled ZYZ decomposition
         of a single-qubit operation is correct"""
@@ -182,16 +190,16 @@ class TestControlledDecompositionZYZ:
         """Test that the operations in the decomposition are correct."""
         phi, theta, omega = 0.123, 0.456, 0.789
         op = qml.Rot(phi, theta, omega, wires=0)
-        control_wires = [1, 2, 3]
+        control_wires = [1]
         decomps = ctrl_decomp_zyz(op, Wires(control_wires))
 
         expected_ops = [
             qml.RZ(0.123, wires=0),
             qml.RY(0.456 / 2, wires=0),
-            qml.MultiControlledX(wires=control_wires + [0]),
+            qml.CNOT(wires=control_wires + [0]),
             qml.RY(-0.456 / 2, wires=0),
             qml.RZ(-(0.123 + 0.789) / 2, wires=0),
-            qml.MultiControlledX(wires=control_wires + [0]),
+            qml.CNOT(wires=control_wires + [0]),
             qml.RZ((0.789 - 0.123) / 2, wires=0),
         ]
         assert all(
@@ -201,7 +209,7 @@ class TestControlledDecompositionZYZ:
         assert len(decomps) == 7
 
     @pytest.mark.parametrize("op", su2_ops + unitary_ops)
-    @pytest.mark.parametrize("control_wires", ([1], [1, 2], [1, 2, 3]))
+    @pytest.mark.parametrize("control_wires", ([1], [2], [3]))
     def test_decomp_queues_correctly(self, op, control_wires, tol):
         """Test that any incorrect operations aren't queued when using
         ``ctrl_decomp_zyz``."""
@@ -268,6 +276,7 @@ class TestControlledDecompositionZYZ:
                     + qml.s_prod(1 / np.sqrt(2), qml.PauliX(0))
                 ),
                 [
+                    qml.ctrl(qml.GlobalPhase(phi=-np.pi / 2, wires=0), control=1),
                     qml.RZ(np.pi / 2, wires=0),
                     qml.RY(np.pi / 2, wires=0),
                     qml.CNOT(wires=[1, 0]),
@@ -758,6 +767,207 @@ class TestControlledBisectGeneral:
         assert np.allclose(res, expected, atol=tol, rtol=tol)
 
 
+class TestMultiControlledUnitary:
+    """tests for qml.ops._ops_math.controlled_decompositions._decompose_multicontrolled_unitary"""
+
+    def test_invalid_op_size_error(self):
+        """Tests that an error is raised when op acts on more than one wire"""
+        with pytest.raises(
+            ValueError, match="The target operation must be a single-qubit operation"
+        ):
+            _ = _decompose_multicontrolled_unitary(qml.CNOT([0, 1]), [2])
+
+    def test_invalid_op_matrix(self):
+        """Tests that an error is raised when op does not define a matrix"""
+
+        # pylint: disable=too-few-public-methods
+        class MyOp(qml.operation.Operator):
+            num_wires = 1
+
+        with pytest.raises(
+            ValueError, match="The target operation must be a single-qubit operation"
+        ):
+            _ = _decompose_multicontrolled_unitary(MyOp, [1])
+
+    su2_gen_ops = [
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [0, 1],
+                    [-1, 0],
+                ]
+            ),
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [0, 1j],
+                    [1j, 0],
+                ]
+            ),
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [1j, 1j],
+                    [1j, -1j],
+                ]
+            )
+            * 2**-0.5,
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [1, 1],
+                    [-1, 1],
+                ]
+            )
+            * 2**-0.5,
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [1 + 2j, -3 + 4j],
+                    [3 + 4j, 1 - 2j],
+                ]
+            )
+            * 30**-0.5,
+            wires=0,
+        ),
+    ]
+
+    gen_ops = [
+        qml.PauliX(0),
+        qml.PauliZ(0),
+        qml.Hadamard(0),
+        qml.Rot(0.123, 0.456, 0.789, wires=0),
+    ]
+
+    @pytest.mark.parametrize("op", gen_ops + su2_gen_ops)
+    @pytest.mark.parametrize("control_wires", cw5)
+    def test_decomposition_circuit(self, op, control_wires, tol):
+        """Tests that the controlled decomposition of a single-qubit operation
+        behaves as expected in a quantum circuit"""
+        dev = qml.device("default.qubit", wires=max(control_wires) + 1)
+
+        @qml.qnode(dev)
+        def decomp_circuit():
+            qml.broadcast(unitary=qml.Hadamard, pattern="single", wires=control_wires)
+            _decompose_multicontrolled_unitary(op, Wires(control_wires))
+            return qml.probs()
+
+        @qml.qnode(dev)
+        def expected_circuit():
+            qml.broadcast(unitary=qml.Hadamard, pattern="single", wires=control_wires)
+            qml.ctrl(op, control_wires)
+            return qml.probs()
+
+        res = decomp_circuit()
+        expected = expected_circuit()
+        assert np.allclose(res, expected, atol=tol, rtol=tol)
+
+    controlled_wires = tuple(list(range(1, 1 + n)) for n in range(0, 2))
+
+    @pytest.mark.parametrize("op", gen_ops + su2_gen_ops)
+    @pytest.mark.parametrize("control_wires", controlled_wires)
+    def test_auto_select_wires(self, op, control_wires):
+        """
+        Test that the auto selection is correct and optimal.
+        """
+        if len(control_wires) == 0:
+            expected = [op]
+        else:  # we only have zero or one control wires
+            expected = ctrl_decomp_zyz(op, Wires(control_wires))
+
+        res = _decompose_multicontrolled_unitary(op, Wires(control_wires))
+        assert equal_list(res, expected)
+
+    controlled_wires = tuple(list(range(2, 1 + n)) for n in range(3, 7))
+
+    @pytest.mark.parametrize("op", gen_ops + su2_gen_ops)
+    @pytest.mark.parametrize("control_wires", controlled_wires)
+    def test_auto_select_su2(self, op, control_wires):
+        """
+        Test that the auto selection is correct and optimal.
+        """
+        if _is_single_qubit_special_unitary(op):
+            expected = ctrl_decomp_bisect(op, Wires(control_wires))
+        else:
+            target_wire = op.wires
+            expected = _decompose_recursive(op, 1.0, Wires(control_wires), target_wire, Wires([]))
+
+        res = _decompose_multicontrolled_unitary(op, Wires(control_wires))
+        assert equal_list(res, expected)
+
+    @pytest.mark.parametrize("op", gen_ops + su2_gen_ops)
+    @pytest.mark.parametrize("control_wires", cw5)
+    def test_decomposition_matrix(self, op, control_wires, tol):
+        """Tests that the matrix representation of the controlled decomposition
+        of a single-qubit operation is correct"""
+
+        expected_op = qml.ctrl(op, control_wires)
+        res = qml.matrix(
+            record_from_list(_decompose_multicontrolled_unitary), wire_order=control_wires + [0]
+        )(op, Wires(control_wires))
+        expected = expected_op.matrix()
+
+        assert np.allclose(res, expected, atol=tol, rtol=tol)
+
+
+class TestControlledUnitaryRecursive:
+    """tests for qml.ops._decompose_recursive"""
+
+    gen_ops = [
+        qml.PauliX(0),
+        qml.PauliZ(0),
+        qml.Hadamard(0),
+    ]
+    controlled_wires = tuple(list(range(1, 1 + n)) for n in range(1, 6))
+
+    @pytest.mark.parametrize("op", gen_ops)
+    @pytest.mark.parametrize("control_wires", controlled_wires)
+    def test_decomposition_circuit(self, op, control_wires, tol):
+        """Tests that the controlled decomposition of a single-qubit operation
+        behaves as expected in a quantum circuit"""
+        dev = qml.device("default.qubit", wires=max(control_wires) + 1)
+
+        @qml.qnode(dev)
+        def decomp_circuit():
+            qml.broadcast(unitary=qml.Hadamard, pattern="single", wires=control_wires)
+            record_from_list(_decompose_recursive)(
+                op, 1.0, Wires(control_wires), op.wires, Wires([])
+            )
+            return qml.probs()
+
+        @qml.qnode(dev)
+        def expected_circuit():
+            qml.broadcast(unitary=qml.Hadamard, pattern="single", wires=control_wires)
+            qml.ctrl(op, control_wires)
+            return qml.probs()
+
+        res = decomp_circuit()
+        expected = expected_circuit()
+        assert np.allclose(res, expected, atol=tol, rtol=tol)
+
+    @pytest.mark.parametrize("op", gen_ops)
+    @pytest.mark.parametrize("control_wires", controlled_wires)
+    def test_decomposition_matrix(self, op, control_wires, tol):
+        """Tests that the matrix representation of the controlled decomposition
+        of a single-qubit operation is correct"""
+
+        expected_op = qml.ctrl(op, control_wires)
+        res = qml.matrix(record_from_list(_decompose_recursive), wire_order=control_wires + [0])(
+            op, 1.0, Wires(control_wires), op.wires, Wires([])
+        )
+        expected = expected_op.matrix()
+
+        assert np.allclose(res, expected, atol=tol, rtol=tol)
+
+
 class TestMCXDecomposition:
     @pytest.mark.parametrize("n_ctrl_wires", range(3, 6))
     def test_decomposition_with_many_workers(self, n_ctrl_wires):
@@ -820,6 +1030,35 @@ class TestMCXDecomposition:
         ).T
         assert np.allclose(u, np.eye(2 ** (n_ctrl_wires + 1)))
 
+    @pytest.mark.parametrize("n_ctrl_wires", range(3, 6))
+    def test_decomposition_with_no_workers(self, n_ctrl_wires):
+        """Test that the decomposed MultiControlledX gate performs the same unitary as the
+        matrix-based version by checking if U^dagger U applies the identity to each basis
+        state. This test focuses on the case where there is no work wires."""
+
+        # pylint: disable=protected-access
+        control_wires = Wires(range(n_ctrl_wires))
+        target_wire = n_ctrl_wires
+
+        dev = qml.device("default.qubit", wires=n_ctrl_wires + 1)
+        with qml.queuing.AnnotatedQueue() as q:
+            decompose_mcx(control_wires, Wires(target_wire), work_wires=Wires([]))
+        tape = qml.tape.QuantumScript.from_queue(q)
+        tape = tape.expand(depth=1)
+
+        @qml.qnode(dev)
+        def f(bitstring):
+            qml.BasisState(bitstring, wires=range(n_ctrl_wires + 1))
+            qml.MultiControlledX(wires=list(control_wires) + [target_wire])
+            for op in tape.operations:
+                op.queue()
+            return qml.probs(wires=range(n_ctrl_wires + 1))
+
+        u = np.array(
+            [f(np.array(b)) for b in itertools.product(range(2), repeat=n_ctrl_wires + 1)]
+        ).T
+        assert np.allclose(u, np.eye(2 ** (n_ctrl_wires + 1)))
+
 
 def test_ControlledQubitUnitary_has_decomposition_correct():
     """Test that ControlledQubitUnitary reports has_decomposition=False if it is False"""
diff --git a/tests/ops/qubit/test_non_parametric_ops.py b/tests/ops/qubit/test_non_parametric_ops.py
index 3ca810fc5..39936c56e 100644
--- a/tests/ops/qubit/test_non_parametric_ops.py
+++ b/tests/ops/qubit/test_non_parametric_ops.py
@@ -717,7 +717,7 @@ class TestMultiControlledX:
         assert qml.math.allclose(mpmct_state, pauli_x_state)
 
     def test_decomposition_not_enough_wires(self):
-        """Test that the decomposition raises an error if the number of wires"""
+        """Test that the decomposition raises an error if the number of wires is lower than two"""
         with pytest.raises(ValueError, match="Wrong number of wires"):
             qml.MultiControlledX.compute_decomposition((0,), control_values=[1])
 
@@ -792,16 +792,6 @@ class TestMultiControlledX:
 
         assert qml.math.allclose(mpmct_state, pauli_x_state)
 
-    def test_not_enough_workers(self):
-        """Test that a ValueError is raised when more than 2 control wires are to be decomposed with
-        no work wires supplied"""
-        control_target_wires = range(4)
-        op = qml.MultiControlledX(wires=control_target_wires)
-
-        match = "At least one work wire is required to decompose operation: MultiControlledX"
-        with pytest.raises(ValueError, match=match):
-            op.decomposition()
-
     def test_not_unique_wires(self):
         """Test that a ValueError is raised when work_wires is not complementary to control_wires"""
         control_target_wires = range(4)