Allow multi controlled x decomposition without working wires (#5735)

**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]

---------

Co-authored-by: Mudit Pandey <mudit.pandey@xanadu.ai>
Co-authored-by: Christina Lee <christina@xanadu.ai>

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)

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,

pennylane/ops/op_math/controlled_decompositions.py

@@ -16,18 +16,26 @@
 """
 
 from copy import copy
-from typing import Tuple
 
 import numpy as np
 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,78 +158,79 @@ 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:
         raise ValueError(
             "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):

pennylane/ops/op_math/controlled_ops.py

@@ -1202,19 +1202,10 @@ def compute_decomposition(wires=None, work_wires=None, control_values=None, **kw
             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]