Repository: PennyLaneAI/pennylane
Pull request #3851 (2023-03-01 19:30:33)

Implement multi-controlled SU(2) decompositions based on Vale et al. (2023)
**Context:**

Currently only the ZYZ multi-controlled decomposition from Barenco et al. (1995) is implemented. More decomposition algorithms would provide more options which can potentially run on more devices or be more efficient.

**Description of the Change:**

Adds the specialized real off-diagonal and real main-diagonal decompositions which use 16n+O(1) CNOTs, and the generalized decomposition which uses 20n+O(1) CNOTs. n here is the number of control wires. All of them share the following properties:

* Works for any number of control wires n >= 2 (n=0 and n=1 are trivial cases, please do not use these functions)
* Target operation must be on 1 wire and have a matrix which is SU(2)
* Does not need or use auxiliary/work wires
* Outputs `MultiControlledX`, `QubitUnitary`, and `Hadamard`

These are described in https://arxiv.org/abs/2302.06377

**Benefits:**

It uses less CNOTs than the method from Barenco et al. (1995). As a candidate algorithm for auto-selecting decompositions, it can reduce the cost of some multi-controlled operations.

Since this algorithm is applicable to Pauli X as well and requires no auxiliary wires, it can be used as a potential decomposition for `MultiControlledX` if there are no auxiliary wires, however, this is significantly more expensive than if auxiliary wires are available.

**Possible Drawbacks:**

Relies on the ability to directly run or further decompose arbitrary single qubit unitary operations. If the single qubit unitary operations are horrible for a device, this method can actually be worse than Barenco et al. (1995), which requires only some constant number of arbitrary Y and Z rotations.

Only works for SU(2), not U(2). Operations not in SU(2) have their global phase removed to coerce them to SU(2).

Not differentiable due to the method of calculating the single qubit unitary matrices. As implemented, will only select these algorithms when using numpy.

**Related GitHub Issues:**

Fixes bug https://github.com/PennyLaneAI/pennylane/issues/3878

**Notes:**

Decomposition hints will be added in a future PR. For now the general algorithm outputs 6 plain `MultiControlledX` which are naively decomposed to 24n+O(1) CNOTs. With decomposition hints, we could tell the `MultiControlledX` to decompose in a certain way that makes the cancellation of some CNOTs more obvious. After fully decomposing and applying the optimizer, there will be only 20n+O(1) CNOTs, matching the complexity claim from the paper.

mudit2812: > If the single qubit unitary operations are horrible for a device, this method can actually be worse than Barenco et al. (1995)

Just wanted to add that realistically, any decomposition for multi-controlled operations is likely going to be most relevant for:
* hardware devices that have a limited gate set to work with, or
* simulators that have been given compilation constraints to limit the allowed gates.

In either case, I wouldn't worry too much about using `QubitUnitary` for the operations applied to the target wire as it can further be decomposed using the single-qubit ZYZ decomposition during compilation, and for a fully decomposed multicontrolled operation, the number of CNOT operations would be far more relevant due to two-qubit operations contributing to decoherence much more significantly.

Obviously, optimizations should be added where feasible, but worst case there will be 10 single-qubit operations to apply to the target wire if the qubit unitaries are decomposed and your implementation uses the $\tilde{H}$ gate from Vale et al. (assuming the gates supported by the device include `Hadamard`, `RX`, `RY`, `RZ`). That would be somewhat insignificant compared to potentially 60+ CNOTs for 3 or more control wires.

Great work so far!
komiamiko: Tests are passing. All 3 algorithms should be implemented correctly now.
komiamiko: Actually, there is one more thing. As the paper noted, the algorithm should achieve 20n+O(1) in the general case, but as implemented it is 24n+O(1). This is because there are some more cancellations - part of the first MCX2 cancels with part of the second MCX2. I don't expect the optimizer to be smart enough to figure this out, so either we should decompose those explicitly with the cancelling applied, or add some kind of hint to `MultiControlledX` so it decomposes in a way that the optimizer can more easily see what will cancel out.
komiamiko: Decomposition hints will be left to future work, and is out of scope for this specific PR.
Work I do need to do though is not use magic to compute the B matrix.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3851?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#3851](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3851?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (0a1a548) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/f5e2ee90f06b2cf007549e43b6e0310902d4504e?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (f5e2ee9) will **decrease** coverage by `0.02%`.
> The diff coverage is `97.48%`.

```diff
@@            Coverage Diff             @@
##           master    #3851      +/-   ##
==========================================
- Coverage   99.71%   99.70%   -0.02%     
==========================================
  Files         338      338              
  Lines       29537    29693     +156     
==========================================
+ Hits        29452    29604     +152     
- Misses         85       89       +4     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3851?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [.../transforms/decompositions/single\_qubit\_unitary.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3851?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvZGVjb21wb3NpdGlvbnMvc2luZ2xlX3F1Yml0X3VuaXRhcnkucHk=) | `100.00% <ø> (ø)` | |
| [pennylane/ops/op\_math/controlled\_decompositions.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3851?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL2NvbnRyb2xsZWRfZGVjb21wb3NpdGlvbnMucHk=) | `97.80% <97.43%> (-2.20%)` | :arrow_down: |
| [pennylane/ops/op\_math/controlled.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/3851?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL2NvbnRyb2xsZWQucHk=) | `100.00% <100.00%> (ø)` | |

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)

komiamiko: For future PR: make MultiControlledX use bisect when there are 0 work wires
mudit2812: Now that the bisect decomposition is the go-to for 2 or more control wires, the only case where the ZYZ decomposition would be used is for a single control wire. This is actually a good thing, because for a single control wire, the ZYZ decomposition can be updated as follows:
* Replace `MultiControlledX` with `CNOT`.
* Add a `PhaseShift` to the control wire with the global phase of U(2) target operations, which would eliminate the phase that the ZYZ decomposition was adding. This will make the decomposition exact rather than having a global phase.

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 06ebf21a1..b657c7b17 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -57,6 +57,10 @@
   to be set as trainable), and therefore the derivatives are computed more efficiently.
   [(#3697)](https://github.com/PennyLaneAI/pennylane/pull/3697)
 
+* 3 new decomposition algorithms are added for n-controlled operations with single-qubit target,
+  and are selected automatically when they produce a better result. They can be accessed via
+  `ops.op_math.ctrl_decomp_bisect`.
+  [(#3851)](https://github.com/PennyLaneAI/pennylane/pull/3851)
 
 <h3>Breaking changes</h3>
 
@@ -110,6 +114,7 @@
 
 This release contains contributions from (in alphabetical order):
 
+Komi Amiko
 Utkarsh Azad
 Lillian M. A. Frederiksen
 Soran Jahangiri
diff --git a/pennylane/ops/op_math/__init__.py b/pennylane/ops/op_math/__init__.py
index 511f59841..ec314dc90 100644
--- a/pennylane/ops/op_math/__init__.py
+++ b/pennylane/ops/op_math/__init__.py
@@ -68,6 +68,7 @@ Decompositions
     :toctree: api
 
     ~ctrl_decomp_zyz
+    ~ctrl_decomp_bisect
 
 """
 
@@ -82,7 +83,7 @@ from .prod import Prod, prod
 from .sprod import SProd, s_prod
 from .sum import Sum, op_sum, sum
 from .symbolicop import ScalarSymbolicOp, SymbolicOp
-from .controlled_decompositions import ctrl_decomp_zyz
+from .controlled_decompositions import ctrl_decomp_zyz, ctrl_decomp_bisect
 
 controlled_qubit_ops = {
     "ControlledQubitUnitary",
diff --git a/pennylane/ops/op_math/controlled.py b/pennylane/ops/op_math/controlled.py
index 24d5fac02..851622bb9 100644
--- a/pennylane/ops/op_math/controlled.py
+++ b/pennylane/ops/op_math/controlled.py
@@ -30,7 +30,10 @@ from pennylane.operation import Operator
 from pennylane.wires import Wires
 
 from .symbolicop import SymbolicOp
-from .controlled_decompositions import ctrl_decomp_zyz
+from .controlled_decompositions import (
+    ctrl_decomp_zyz,
+    ctrl_decomp_bisect,
+)
 
 
 def ctrl(op, control, control_values=None, work_wires=None):
@@ -556,11 +559,23 @@ def _decompose_no_control_values(op: "operation.Operator") -> List["operation.Op
     no decomposition.
     """
     if len(op.control_wires) == 1 and hasattr(op.base, "_controlled"):
-        return [op.base._controlled(op.control_wires[0])]
+        result = op.base._controlled(op.control_wires[0])
+        # disallow decomposing to itself
+        # pylint: disable=unidiomatic-typecheck
+        if type(result) != type(op):
+            return [result]
+        qml.QueuingManager.remove(result)
     if isinstance(op.base, qml.PauliX):
-        if len(op.control_wires) == 2:
-            return [qml.Toffoli(op.active_wires)]
+        # has some special case handling of its own for further decomposition
         return [qml.MultiControlledX(wires=op.active_wires, work_wires=op.work_wires)]
+    if (
+        len(op.base.wires) == 1
+        and len(op.control_wires) >= 2
+        and getattr(op.base, "has_matrix", False)
+        and qmlmath.get_interface(*op.data) == "numpy"  # as implemented, not differentiable
+    ):
+        # Bisect algorithms use CNOTs and single qubit unitary
+        return ctrl_decomp_bisect(op.base, op.control_wires)
     if len(op.base.wires) == 1 and getattr(op.base, "has_matrix", False):
         return ctrl_decomp_zyz(op.base, op.control_wires)
 
diff --git a/pennylane/ops/op_math/controlled_decompositions.py b/pennylane/ops/op_math/controlled_decompositions.py
index 78232e07d..b414c70d6 100644
--- a/pennylane/ops/op_math/controlled_decompositions.py
+++ b/pennylane/ops/op_math/controlled_decompositions.py
@@ -15,9 +15,111 @@
 This submodule defines functions to decompose controlled operations
 """
 
+from copy import copy
+import numpy as np
+import numpy.linalg as npl
 import pennylane as qml
 from pennylane.operation import Operator
 from pennylane.wires import Wires
+from pennylane import math
+
+
+def _convert_to_su2(U):
+    r"""Convert a 2x2 unitary matrix to :math:`SU(2)`.
+
+    Args:
+        U (array[complex]): A matrix, presumed to be :math:`2 \times 2` and unitary.
+        return_global_phase (bool): If `True`, the return will include
+        the global phase. If `False`, only the :math:`SU(2)` representative
+        is returned.
+
+    Returns:
+        array[complex]: A :math:`2 \times 2` matrix in :math:`SU(2)` that is
+        equivalent to U up to a global phase. If ``return_global_phase=True``,
+        a 2-element tuple is returned, with the first element being the
+        :math:`SU(2)` equivalent and the second, the global phase.
+    """
+    # Compute the determinants
+    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
+
+
+def _convert_to_real_diagonal(q: np.ndarray) -> np.ndarray:
+    """
+    Change the phases of Q so the main diagonal is real, and return the modified Q.
+    """
+    exp_angles = np.angle(np.diag(q))
+    return q * np.exp(-1j * exp_angles).reshape((1, 2))
+
+
+def _param_su2(ar: float, ai: float, br: float, bi: float):
+    """
+    Create a matrix in the SU(2) form from complex parameters a, b.
+    The resulting matrix is not guaranteed to be in SU(2), unless |a|^2 + |b|^2 = 1.
+    """
+    return np.array([[ar + 1j * ai, -br + 1j * bi], [br + 1j * bi, ar + 1j * -ai]])
+
+
+def _bisect_compute_a(u: np.ndarray):
+    """
+    Given the U matrix, compute the A matrix such that
+    At x A x At x A x = U
+    where At is the adjoint of A
+    and x is the Pauli X matrix.
+    """
+    x = np.real(u[0, 1])
+    z = u[1, 1]
+    zr = np.real(z)
+    zi = np.imag(z)
+    if np.isclose(zr, -1):
+        # special case [[-1, 0], [0, -1]]
+        # would cause divide by 0 with the other formula, so we use hardcoded solution
+        return np.array([[1, -1], [1, 1]]) * 2**-0.5
+    ar = np.sqrt((np.sqrt((zr + 1) / 2) + 1) / 2)
+    mul = 1 / (2 * np.sqrt((zr + 1) * (np.sqrt((zr + 1) / 2) + 1)))
+    ai = zi * mul
+    br = x * mul
+    bi = 0
+    return _param_su2(ar, ai, br, bi)
+
+
+def _bisect_compute_b(u: np.ndarray):
+    """
+    Given the U matrix, compute the B matrix such that
+    H Bt x B x H = U
+    where Bt is the adjoint of B,
+    H is the Hadamard matrix,
+    and x is the Pauli X matrix.
+    """
+    sqrt = np.sqrt
+    Abs = np.abs
+    w = np.real(u[0, 0])
+    s = np.real(u[1, 0])
+    t = np.imag(u[1, 0])
+    if np.isclose(s, 0):
+        b = 0
+        if np.isclose(t, 0):
+            if w < 0:
+                c = 0
+                d = sqrt(-w)
+            else:
+                c = sqrt(w)
+                d = 0
+        else:
+            c = sqrt(2 - 2 * w) * (-w / 2 - 1 / 2) / t
+            d = sqrt(2 - 2 * w) / 2
+    elif np.isclose(t, 0):
+        b = (1 / 2 - w / 2) * sqrt(2 * w + 2) / s
+        c = sqrt(2 * w + 2) / 2
+        d = 0
+    else:
+        b = sqrt(2) * s * sqrt((1 - w) / (s**2 + t**2)) * Abs(t) / (2 * t)
+        c = sqrt(2) * sqrt((1 - w) / (s**2 + t**2)) * (w + 1) * Abs(t) / (2 * t)
+        d = -sqrt(2) * sqrt((1 - w) / (s**2 + t**2)) * Abs(t) / 2
+    return _param_su2(c, d, b, 0)
 
 
 def ctrl_decomp_zyz(target_operation: Operator, control_wires: Wires):
@@ -112,3 +214,230 @@ def ctrl_decomp_zyz(target_operation: Operator, control_wires: Wires):
         decomp.append(qml.RZ((omega - phi) / 2, wires=target_wire))
 
     return decomp
+
+
+def _ctrl_decomp_bisect_od(
+    u: np.ndarray,
+    target_wire: Wires,
+    control_wires: Wires,
+):
+    """Decompose the controlled version of a target single-qubit operation
+
+    Not backpropagation compatible (as currently implemented). Use only with numpy.
+
+    This function decomposes a controlled single-qubit target operation using the
+    decomposition defined in section 3.1, Theorem 1 of
+    `Vale et al. (2023) <https://arxiv.org/abs/2302.06377>`_.
+
+    The target operation's matrix must have a real off-diagonal for this specialized method to work.
+
+    .. warning:: This method will add a global phase for target operations that do not
+        belong to the SU(2) group.
+
+    Args:
+        u (np.ndarray): the target operation's matrix
+        target_wire (~.wires.Wires): the target wire of the controlled operation
+        control_wires (~.wires.Wires): the control wires of the operation
+
+    Returns:
+        list[Operation]: the decomposed operations
+
+    Raises:
+        ValueError: if ``u`` does not have a real off-diagonal
+
+    """
+    ui = np.imag(u)
+    if not np.isclose(ui[1, 0], 0) or not np.isclose(ui[0, 1], 0):
+        raise ValueError(f"Target operation's matrix must have real off-diagonal, but it is {u}")
+
+    a = _bisect_compute_a(u)
+
+    mid = (len(control_wires) + 1) // 2  # for odd n, make control_k1 bigger
+    control_k1 = control_wires[:mid]
+    control_k2 = control_wires[mid:]
+
+    def component():
+        return [
+            qml.MultiControlledX(wires=control_k1 + target_wire, work_wires=control_k2),
+            qml.QubitUnitary(a, target_wire),
+            qml.MultiControlledX(wires=control_k2 + target_wire, work_wires=control_k1),
+            qml.adjoint(qml.QubitUnitary(a, target_wire)),
+        ]
+
+    return component() + component()
+
+
+def _ctrl_decomp_bisect_md(
+    u: np.ndarray,
+    target_wire: Wires,
+    control_wires: Wires,
+):
+    """Decompose the controlled version of a target single-qubit operation
+
+    Not backpropagation compatible (as currently implemented). Use only with numpy.
+
+    This function decomposes a controlled single-qubit target operation using the
+    decomposition defined in section 3.1, Theorem 2 of
+    `Vale et al. (2023) <https://arxiv.org/abs/2302.06377>`_.
+
+    The target operation's matrix must have a real main-diagonal for this specialized method to work.
+
+    .. warning:: This method will add a global phase for target operations that do not
+        belong to the SU(2) group.
+
+    Args:
+        u (np.ndarray): the target operation's matrix
+        target_wire (~.wires.Wires): the target wire of the controlled operation
+        control_wires (~.wires.Wires): the control wires of the operation
+
+    Returns:
+        list[Operation]: the decomposed operations
+
+    Raises:
+        ValueError: if ``u`` does not have a real main-diagonal
+
+    """
+    ui = np.imag(u)
+    if not np.isclose(ui[0, 0], 0) or not np.isclose(ui[1, 1], 0):
+        raise ValueError(f"Target operation's matrix must have real main-diagonal, but it is {u}")
+
+    h_matrix = qml.Hadamard.compute_matrix()
+    mod_u = h_matrix @ u @ h_matrix
+
+    decomposition = [qml.Hadamard(target_wire)]
+    decomposition += _ctrl_decomp_bisect_od(mod_u, target_wire, control_wires)
+    decomposition.append(qml.Hadamard(target_wire))
+
+    return decomposition
+
+
+def _ctrl_decomp_bisect_general(
+    u: np.ndarray,
+    target_wire: Wires,
+    control_wires: Wires,
+):
+    """Decompose the controlled version of a target single-qubit operation
+
+    Not backpropagation compatible (as currently implemented). Use only with numpy.
+
+    This function decomposes a controlled single-qubit target operation using the
+    decomposition defined in section 3.2 of
+    `Vale et al. (2023) <https://arxiv.org/abs/2302.06377>`_.
+
+    .. warning:: This method will add a global phase for target operations that do not
+        belong to the SU(2) group.
+
+    Args:
+        u (np.ndarray): the target operation's matrix
+        target_wire (~.wires.Wires): the target wire of the controlled operation
+        control_wires (~.wires.Wires): the control wires of the operation
+
+    Returns:
+        list[Operation]: the decomposed operations
+    """
+    x_matrix = qml.PauliX.compute_matrix()
+    h_matrix = qml.Hadamard.compute_matrix()
+    alternate_h_matrix = x_matrix @ h_matrix @ x_matrix
+
+    d, q = npl.eig(u)
+    d = np.diag(d)
+    q = _convert_to_real_diagonal(q)
+    b = _bisect_compute_b(q)
+    c1 = b @ alternate_h_matrix
+    c2t = b @ h_matrix
+
+    mid = (len(control_wires) + 1) // 2  # for odd n, make control_k1 bigger
+    control_k1 = control_wires[:mid]
+    control_k2 = control_wires[mid:]
+
+    component = [
+        qml.QubitUnitary(c2t, target_wire),
+        qml.MultiControlledX(wires=control_k2 + target_wire, work_wires=control_k1),
+        qml.adjoint(qml.QubitUnitary(c1, target_wire)),
+        qml.MultiControlledX(wires=control_k1 + target_wire, work_wires=control_k2),
+    ]
+
+    od_decomp = _ctrl_decomp_bisect_od(d, target_wire, control_wires)
+
+    # cancel two identical multicontrolled x gates
+    qml.QueuingManager.remove(component[3])
+    qml.QueuingManager.remove(od_decomp[0])
+
+    adjoint_component = [qml.adjoint(copy(op), lazy=False) for op in reversed(component)]
+
+    return component[0:3] + od_decomp[1:] + adjoint_component
+
+
+def ctrl_decomp_bisect(
+    target_operation: Operator,
+    control_wires: Wires,
+):
+    """Decompose the controlled version of a target single-qubit operation
+
+    Not backpropagation compatible (as currently implemented). Use only with numpy.
+
+    Automatically selects the best algorithm based on the matrix (uses specialized more efficient
+    algorithms if the matrix has a certain form, otherwise falls back to the general algorithm).
+    These algorithms are defined in section 3.1 and 3.2 of
+    `Vale et al. (2023) <https://arxiv.org/abs/2302.06377>`_.
+
+    .. 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
+
+    Returns:
+        list[Operation]: the decomposed operations
+
+    Raises:
+        ValueError: if ``target_operation`` is not a single-qubit operation
+
+    **Example:**
+
+    >>> op = qml.T(0) # uses OD algorithm
+    >>> print(qml.draw(ctrl_decomp_bisect, wire_order=(0,1,2,3,4,5))(op, (1,2,3,4,5)))
+    0: ─╭X──U(M0)─╭X──U(M0)†─╭X──U(M0)─╭X──U(M0)†─┤
+    1: ─├●────────│──────────├●────────│──────────┤
+    2: ─├●────────│──────────├●────────│──────────┤
+    3: ─╰●────────│──────────╰●────────│──────────┤
+    4: ───────────├●───────────────────├●─────────┤
+    5: ───────────╰●───────────────────╰●─────────┤
+    >>> op = qml.QubitUnitary([[0,1j],[1j,0]], 0) # uses MD algorithm
+    >>> print(qml.draw(ctrl_decomp_bisect, wire_order=(0,1,2,3,4,5))(op, (1,2,3,4,5)))
+    0: ──H─╭X──U(M0)─╭X──U(M0)†─╭X──U(M0)─╭X──U(M0)†──H─┤
+    1: ────├●────────│──────────├●────────│─────────────┤
+    2: ────├●────────│──────────├●────────│─────────────┤
+    3: ────╰●────────│──────────╰●────────│─────────────┤
+    4: ──────────────├●───────────────────├●────────────┤
+    5: ──────────────╰●───────────────────╰●────────────┤
+    >>> op = qml.Hadamard(0) # uses general algorithm
+    >>> print(qml.draw(ctrl_decomp_bisect, wire_order=(0,1,2,3,4,5))(op, (1,2,3,4,5)))
+    0: ──U(M0)─╭X──U(M1)†──U(M2)─╭X──U(M2)†─╭X──U(M2)─╭X──U(M2)†─╭X──U(M1)─╭X──U(M0)─┤
+    1: ────────│─────────────────│──────────├●────────│──────────├●────────│─────────┤
+    2: ────────│─────────────────│──────────├●────────│──────────├●────────│─────────┤
+    3: ────────│─────────────────│──────────╰●────────│──────────╰●────────│─────────┤
+    4: ────────├●────────────────├●───────────────────├●───────────────────├●────────┤
+    5: ────────╰●────────────────╰●───────────────────╰●───────────────────╰●────────┤
+
+    """
+    if len(target_operation.wires) > 1:
+        raise ValueError(
+            "The target operation must be a single-qubit operation, instead "
+            f"got {target_operation}."
+        )
+    target_matrix = target_operation.matrix()
+    target_wire = target_operation.wires
+
+    target_matrix = _convert_to_su2(target_matrix)
+    target_matrix_imag = np.imag(target_matrix)
+
+    if np.isclose(target_matrix_imag[1, 0], 0) and np.isclose(target_matrix_imag[0, 1], 0):
+        # Real off-diagonal specialized algorithm - 16n+O(1) CNOTs
+        return _ctrl_decomp_bisect_od(target_matrix, target_wire, control_wires)
+    if np.isclose(target_matrix_imag[0, 0], 0) and np.isclose(target_matrix_imag[1, 1], 0):
+        # Real main-diagonal specialized algorithm - 16n+O(1) CNOTs
+        return _ctrl_decomp_bisect_md(target_matrix, target_wire, control_wires)
+    # General algorithm - 20n+O(1) CNOTs
+    return _ctrl_decomp_bisect_general(target_matrix, target_wire, control_wires)
diff --git a/pennylane/transforms/decompositions/single_qubit_unitary.py b/pennylane/transforms/decompositions/single_qubit_unitary.py
index 91f7b5bc9..df8e1a432 100644
--- a/pennylane/transforms/decompositions/single_qubit_unitary.py
+++ b/pennylane/transforms/decompositions/single_qubit_unitary.py
@@ -20,16 +20,17 @@ from pennylane import math
 
 
 def _convert_to_su2(U, return_global_phase=False):
-    r"""Convert a 2x2 unitary matrix to :math:`SU(2)`.
+    r"""Convert a 2x2 unitary matrix to :math:`SU(2)`. (batched operation)
 
     Args:
-        U (array[complex]): A matrix, presumed to be :math:`2 \times 2` and unitary.
+        U (array[complex]): A matrix with a batch dimension, presumed to be
+        of shape :math:`n \times 2 \times 2` and unitary for any positive integer n.
         return_global_phase (bool): If `True`, the return will include
         the global phase. If `False`, only the :math:`SU(2)` representative
         is returned.
 
     Returns:
-        array[complex]: A :math:`2 \times 2` matrix in :math:`SU(2)` that is
+        array[complex]: A :math:`n \times 2 \times 2` matrix in :math:`SU(2)` that is
         equivalent to U up to a global phase. If ``return_global_phase=True``,
         a 2-element tuple is returned, with the first element being the
         :math:`SU(2)` equivalent and the second, the global phase.


Test patch:
diff --git a/tests/ops/op_math/test_controlled.py b/tests/ops/op_math/test_controlled.py
index a0e6a32a7..4a1c500c0 100644
--- a/tests/ops/op_math/test_controlled.py
+++ b/tests/ops/op_math/test_controlled.py
@@ -39,6 +39,15 @@ from pennylane.wires import Wires
 # pylint: disable=pointless-statement
 # pylint: disable=expression-not-assigned
 
+
+def equal_list(lhs, rhs):
+    if not isinstance(lhs, list):
+        lhs = [lhs]
+    if not isinstance(rhs, list):
+        rhs = [rhs]
+    return len(lhs) == len(rhs) and all(qml.equal(l, r) for l, r in zip(lhs, rhs))
+
+
 base_num_control_mats = [
     (qml.PauliX("a"), 1, CNOT),
     (qml.PauliZ("a"), 1, CZ),
@@ -760,7 +769,7 @@ class TestHelperMethod:
         op = Controlled(qml.PauliX("c"), ("a", 2))
         decomp = _decompose_no_control_values(op)
         assert len(decomp) == 1
-        assert qml.equal(decomp[0], qml.Toffoli(("a", 2, "c")))
+        assert equal_list(decomp, qml.MultiControlledX(wires=("a", 2, "c")))
 
     def test_multicontrolledx(self):
         """Test case when PauliX has many controls."""
@@ -821,9 +830,8 @@ class TestDecomposition:
         op = Controlled(base, (0, 1), (True, False))
 
         decomp = op.expand().circuit if test_expand else op.decomposition()
-        expected = [qml.PauliX(1), qml.Toffoli((0, 1, 2)), qml.PauliX(1)]
-        for op1, op2 in zip(decomp, expected):
-            assert qml.equal(op1, op2)
+        expected = [qml.PauliX(1), qml.MultiControlledX(wires=(0, 1, 2)), qml.PauliX(1)]
+        assert equal_list(decomp, expected)
 
     def test_no_control_values_special_decomp(self, test_expand):
         """Test a case with no control values but a special decomposition."""
@@ -856,36 +864,6 @@ class TestDecomposition:
             # pylint: disable=unused-variable
             decomp = op.expand().circuit if test_expand else op.decomposition()
 
-    def test_zyz_decomp_no_control_values(self, test_expand):
-        """Test that the ZYZ decomposition is used for single qubit target operations
-        when other decompositions aren't available."""
-
-        base = qml.RX(0.123, wires="a")
-        op = Controlled(base, (0, 1, 2))
-
-        assert op.has_decomposition
-        decomp = op.expand().circuit if test_expand else op.decomposition()
-        expected = qml.ops.ctrl_decomp_zyz(base, (0, 1, 2))
-        for op, exp_op in zip(decomp, expected):
-            assert qml.equal(op, exp_op)
-
-    def test_zyz_decomp_control_values(self, test_expand):
-        """Test that the ZYZ decomposition is used for single qubit target operations
-        when other decompositions aren't available and control values are present."""
-
-        base = qml.RX(0.123, wires="a")
-        op = Controlled(base, (0, 1, 2), control_values=[True, False, False])
-
-        assert op.has_decomposition
-        decomp = op.expand().circuit if test_expand else op.decomposition()
-        assert qml.equal(qml.PauliX(1), decomp[0])
-        assert qml.equal(qml.PauliX(2), decomp[1])
-        assert qml.equal(qml.PauliX(1), decomp[-2])
-        assert qml.equal(qml.PauliX(2), decomp[-1])
-        expected = qml.ops.ctrl_decomp_zyz(base, (0, 1, 2))
-        for op, exp_op in zip(decomp[2:-2], expected):
-            assert qml.equal(op, exp_op)
-
 
 class TestArithmetic:
     """Test arithmetic decomposition methods."""
@@ -1440,7 +1418,7 @@ def test_multi_ctrl():
     op = tape.operations[0]
     assert isinstance(op, Controlled)
     new_tape = tape.expand(depth=1)
-    assert qml.equal(new_tape[0], qml.Toffoli(wires=[3, 7, 0]))
+    assert qml.equal(new_tape[0], qml.MultiControlledX(wires=[3, 7, 0]))
 
 
 def test_ctrl_with_qnode():
@@ -1512,31 +1490,63 @@ def test_diagonal_ctrl():
     )
 
 
-def test_qubit_unitary():
+@pytest.mark.parametrize(
+    "M",
+    [
+        qml.PauliX.compute_matrix(),
+        qml.PauliY.compute_matrix(),
+        qml.PauliZ.compute_matrix(),
+        qml.Hadamard.compute_matrix(),
+        np.array(
+            [
+                [1 + 2j, -3 + 4j],
+                [3 + 4j, 1 - 2j],
+            ]
+        )
+        * 30**-0.5,
+    ],
+)
+def test_qubit_unitary(M):
     """Test ctrl on QubitUnitary and ControlledQubitUnitary"""
-    M = np.array([[1.0, 1.0], [1.0, -1.0]]) / np.sqrt(2.0)
     with qml.queuing.AnnotatedQueue() as q_tape:
         ctrl(qml.QubitUnitary, 1)(M, wires=0)
 
     tape = QuantumScript.from_queue(q_tape)
-    tape = tape.expand(3, stop_at=lambda op: not isinstance(op, Controlled))
-
     expected = qml.ControlledQubitUnitary(M, control_wires=1, wires=0)
+    assert equal_list(list(tape), expected)
 
-    assert qml.equal(tape[0], expected)
-
-
-def test_controlledqubitunitary():
+    # causes decomposition into more basic operators
+    tape = tape.expand(3, stop_at=lambda op: not isinstance(op, Controlled))
+    assert not equal_list(list(tape), expected)
+
+
+@pytest.mark.parametrize(
+    "M",
+    [
+        qml.PauliX.compute_matrix(),
+        qml.PauliY.compute_matrix(),
+        qml.PauliZ.compute_matrix(),
+        qml.Hadamard.compute_matrix(),
+        np.array(
+            [
+                [1 + 2j, -3 + 4j],
+                [3 + 4j, 1 - 2j],
+            ]
+        )
+        * 30**-0.5,
+    ],
+)
+def test_controlledqubitunitary(M):
     """Test ctrl on ControlledQubitUnitary."""
-    M = np.array([[1.0, 1.0], [1.0, -1.0]]) / np.sqrt(2.0)
     with qml.queuing.AnnotatedQueue() as q_tape:
         ctrl(qml.ControlledQubitUnitary, 1)(M, control_wires=2, wires=0)
 
     tape = QuantumScript.from_queue(q_tape)
+    # will immediately decompose according to selected decomposition algorithm
     tape = tape.expand(3, stop_at=lambda op: not isinstance(op, Controlled))
 
-    expected = qml.ControlledQubitUnitary(M, control_wires=[2, 1], wires=0)
-    qml.equal(tape[0], expected)
+    expected = qml.ControlledQubitUnitary(M, control_wires=[2, 1], wires=0).decomposition()
+    assert equal_list(list(tape), expected)
 
 
 def test_no_control_defined():
@@ -1742,7 +1752,7 @@ def test_multi_ctrl_values(ctrl_values):
         for i, j in enumerate(ctrl_val):
             if not bool(j):
                 exp_op.append(qml.PauliX(ctrl_wires[i]))
-        exp_op.append(qml.Toffoli(wires=[3, 7, 0]))
+        exp_op.append(qml.MultiControlledX(wires=[3, 7, 0]))
         for i, j in enumerate(ctrl_val):
             if not bool(j):
                 exp_op.append(qml.PauliX(ctrl_wires[i]))
@@ -1757,5 +1767,4 @@ def test_multi_ctrl_values(ctrl_values):
     op = tape.operations[0]
     assert isinstance(op, Controlled)
     new_tape = expand_tape(tape, 1)
-    for op1, op2 in zip(new_tape, expected_ops(ctrl_values)):
-        assert qml.equal(op1, op2)
+    assert equal_list(list(new_tape), expected_ops(ctrl_values))
diff --git a/tests/ops/op_math/test_controlled_decompositions.py b/tests/ops/op_math/test_controlled_decompositions.py
index 7700a5afa..4875ca1a4 100644
--- a/tests/ops/op_math/test_controlled_decompositions.py
+++ b/tests/ops/op_math/test_controlled_decompositions.py
@@ -20,6 +20,55 @@ import numpy as np
 import pennylane as qml
 from pennylane.ops import ctrl_decomp_zyz
 from pennylane.wires import Wires
+from pennylane.ops.op_math.controlled_decompositions import (
+    _ctrl_decomp_bisect_od,
+    _ctrl_decomp_bisect_md,
+    _ctrl_decomp_bisect_general,
+    ctrl_decomp_bisect,
+    _convert_to_su2,
+    _bisect_compute_a,
+    _bisect_compute_b,
+)
+from pennylane.ops.op_math.controlled import (
+    Controlled,
+)
+from pennylane import math
+
+cw5 = tuple(list(range(1, 1 + n)) for n in range(2, 6))
+
+
+def _matrix_adjoint(matrix: np.ndarray):
+    return math.transpose(math.conj(matrix))
+
+
+def record_from_list(func):
+    """
+    Decorates a function to
+    - not record any operators instantiated during the function
+    - record the operators returned by the function, in order
+    - not return the return value
+
+    Not returning the return value is intentional, as the contexts where this is used
+    (qnode, matrix) only care about recorded operators and not the return value.
+    """
+
+    def irecord_from_list(*args, **kwargs):
+        with qml.QueuingManager.stop_recording():
+            decomposition = func(*args, **kwargs)
+
+        if qml.QueuingManager.recording():
+            for iop in decomposition:
+                qml.apply(iop)
+
+    return irecord_from_list
+
+
+def equal_list(lhs, rhs):
+    if not isinstance(lhs, list):
+        lhs = [lhs]
+    if not isinstance(rhs, list):
+        rhs = [rhs]
+    return len(lhs) == len(rhs) and all(qml.equal(l, r) for l, r in zip(lhs, rhs))
 
 
 class TestControlledDecompositionZYZ:
@@ -151,3 +200,462 @@ class TestControlledDecompositionZYZ:
 
         assert len(decomp) == 5
         assert all(qml.equal(o, e) for o, e in zip(decomp, expected))
+
+    @pytest.mark.parametrize("test_expand", [False, True])
+    def test_zyz_decomp_no_control_values(self, test_expand):
+        """Test that the ZYZ decomposition is used for single qubit target operations
+        when other decompositions aren't available."""
+
+        base = qml.QubitUnitary(
+            np.array(
+                [
+                    [1, 1],
+                    [-1, 1],
+                ]
+            )
+            * 2**-0.5,
+            wires="a",
+        )
+        op = Controlled(base, (0,))
+
+        assert op.has_decomposition
+        decomp = (
+            op.expand().expand().circuit if test_expand else op.decomposition()[0].decomposition()
+        )
+        expected = qml.ops.ctrl_decomp_zyz(base, (0,))
+        assert equal_list(decomp, expected)
+
+    @pytest.mark.parametrize("test_expand", [False, True])
+    def test_zyz_decomp_control_values(self, test_expand):
+        """Test that the ZYZ decomposition is used for single qubit target operations
+        when other decompositions aren't available and control values are present."""
+
+        base = qml.QubitUnitary(
+            np.array(
+                [
+                    [1, 1],
+                    [-1, 1],
+                ]
+            )
+            * 2**-0.5,
+            wires="a",
+        )
+        op = Controlled(base, (0,), control_values=[False])
+
+        assert op.has_decomposition
+        decomp = op.expand().circuit if test_expand else op.decomposition()
+        assert len(decomp) == 3
+        assert qml.equal(qml.PauliX(0), decomp[0])
+        assert qml.equal(qml.PauliX(0), decomp[-1])
+        decomp = decomp[1]
+        decomp = decomp.expand().circuit if test_expand else decomp.decomposition()
+        expected = qml.ops.ctrl_decomp_zyz(base, (0,))
+        assert equal_list(decomp, expected)
+
+
+class TestControlledBisectOD:
+    """tests for qml.ops._ctrl_decomp_bisect_od"""
+
+    def test_invalid_op_error(self):
+        """Tests that an error is raised when an invalid operation is passed"""
+        with pytest.raises(ValueError, match="Target operation's matrix must have real"):
+            _ = _ctrl_decomp_bisect_od(_convert_to_su2(qml.Hadamard.compute_matrix()), 0, [1, 2])
+
+    su2_od_ops = [
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [0, 1],
+                    [-1, 0],
+                ]
+            ),
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [1, 1],
+                    [-1, 1],
+                ]
+            )
+            * 2**-0.5,
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [1j, 0],
+                    [0, -1j],
+                ]
+            ),
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [1, 0],
+                    [0, 1],
+                ]
+            ),
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [-1, 0],
+                    [0, -1],
+                ]
+            ),
+            wires=0,
+        ),
+    ]
+
+    od_ops = [
+        qml.PauliZ(0),
+    ]
+
+    @pytest.mark.parametrize("op", su2_od_ops + od_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)
+            record_from_list(_ctrl_decomp_bisect_od)(
+                _convert_to_su2(op.matrix()), op.wires, 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)
+
+    @pytest.mark.parametrize("op", su2_od_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"""
+        assert np.allclose(op.matrix(), _convert_to_su2(op.matrix()), atol=tol, rtol=tol)
+
+        expected_op = qml.ctrl(op, control_wires)
+        res = qml.matrix(record_from_list(_ctrl_decomp_bisect_od), wire_order=control_wires + [0])(
+            op.matrix(), op.wires, Wires(control_wires)
+        )
+        expected = expected_op.matrix()
+
+        assert np.allclose(res, expected, atol=tol, rtol=tol)
+
+    @pytest.mark.parametrize("op", su2_od_ops)
+    def test_decomposed_operators(self, op, tol):
+        """Tests that the operators in the decomposition match expectations."""
+        control_wires = [1, 2, 3, 4, 5]
+
+        su = op.matrix()
+        sx = qml.PauliX.compute_matrix()
+        op_seq = _ctrl_decomp_bisect_od(op.matrix(), op.wires, Wires(control_wires))
+
+        assert len(op_seq) == 8
+
+        mcx1 = qml.MultiControlledX(
+            control_wires=Wires([1, 2, 3]), wires=Wires(0), work_wires=Wires([4, 5])
+        )
+        assert qml.equal(mcx1, op_seq[0])
+        assert qml.equal(mcx1, op_seq[4])
+
+        mcx2 = qml.MultiControlledX(
+            control_wires=Wires([4, 5]), wires=Wires(0), work_wires=Wires([1, 2, 3])
+        )
+        assert qml.equal(mcx2, op_seq[2])
+        assert qml.equal(mcx2, op_seq[6])
+
+        a = op_seq[1].matrix()
+        at = op_seq[3].matrix()
+        a2 = op_seq[5].matrix()
+        at2 = op_seq[7].matrix()
+        assert np.array_equal(a, a2)
+        assert np.array_equal(at, at2)
+
+        i2 = np.identity(2)
+        assert np.allclose(a @ at, i2, atol=tol, rtol=tol)
+        assert np.allclose(at @ a, i2, atol=tol, rtol=tol)
+
+        assert np.allclose(at @ sx @ a @ sx @ at @ sx @ a @ sx, su, atol=tol, rtol=tol)
+
+    @pytest.mark.parametrize("op", su2_od_ops)
+    def test_a_matrix(self, op, tol):
+        """Tests that the A matrix subroutine returns a correct A matrix."""
+        su = op.matrix()
+        sx = qml.PauliX.compute_matrix()
+        a = _bisect_compute_a(su)
+        at = _matrix_adjoint(a)
+        assert np.allclose(at @ sx @ a @ sx @ at @ sx @ a @ sx, su, atol=tol, rtol=tol)
+
+
+class TestControlledBisectMD:
+    """tests for qml.ops._ctrl_decomp_bisect_md"""
+
+    def test_invalid_op_error(self):
+        """Tests that an error is raised when an invalid operation is passed"""
+        with pytest.raises(ValueError, match="Target operation's matrix must have real"):
+            _ = _ctrl_decomp_bisect_md(_convert_to_su2(qml.Hadamard.compute_matrix()), 0, [1, 2])
+
+    su2_md_ops = [
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [0, 1j],
+                    [1j, 0],
+                ]
+            ),
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [0, 1],
+                    [-1, 0],
+                ]
+            ),
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [1, 1],
+                    [-1, 1],
+                ]
+            )
+            * 2**-0.5,
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [-1, 0],
+                    [0, -1],
+                ]
+            ),
+            wires=0,
+        ),
+    ]
+
+    md_ops = [
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [0, 1],
+                    [1, 0],
+                ]
+            ),
+            wires=0,
+        ),
+        qml.QubitUnitary(
+            np.array(
+                [
+                    [1j, 1j],
+                    [-1j, 1j],
+                ]
+            )
+            * 2**-0.5,
+            wires=0,
+        ),
+    ]
+
+    @pytest.mark.parametrize("op", su2_md_ops + md_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)
+            record_from_list(_ctrl_decomp_bisect_md)(
+                _convert_to_su2(op.matrix()), op.wires, 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)
+
+    @pytest.mark.parametrize("op", su2_md_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"""
+        assert np.allclose(op.matrix(), _convert_to_su2(op.matrix()), atol=tol, rtol=tol)
+
+        expected_op = qml.ctrl(op, control_wires)
+        res = qml.matrix(record_from_list(_ctrl_decomp_bisect_md), wire_order=control_wires + [0])(
+            op.matrix(), op.wires, control_wires
+        )
+        expected = expected_op.matrix()
+
+        assert np.allclose(res, expected, atol=tol, rtol=tol)
+
+    @pytest.mark.parametrize("op", su2_md_ops)
+    def test_b_matrix(self, op, tol):
+        """Tests that the B matrix subroutine returns a correct A matrix."""
+        su = op.matrix()
+        sx = qml.PauliX.compute_matrix()
+        sh = qml.Hadamard.compute_matrix()
+        b = _bisect_compute_b(su)
+        bt = _matrix_adjoint(b)
+        assert np.allclose(sh @ bt @ sx @ b @ sx @ sh, su, atol=tol, rtol=tol)
+
+
+class TestControlledBisectGeneral:
+    """tests for qml.ops._ctrl_decomp_bisect_general"""
+
+    def test_invalid_op_error(self):
+        """Tests that an error is raised when an invalid operation is passed"""
+        with pytest.raises(
+            ValueError, match="The target operation must be a single-qubit operation"
+        ):
+            _ = ctrl_decomp_bisect(qml.CNOT([0, 1]), [2])
+
+    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.PauliY(0),
+        qml.PauliZ(0),
+        qml.Hadamard(0),
+        qml.Rot(0.123, 0.456, 0.789, wires=0),
+    ]
+
+    gen_ops_best = [
+        _ctrl_decomp_bisect_md,
+        _ctrl_decomp_bisect_od,
+        _ctrl_decomp_bisect_od,
+        _ctrl_decomp_bisect_general,
+        _ctrl_decomp_bisect_general,
+    ]
+
+    @pytest.mark.parametrize("op", su2_gen_ops + gen_ops)
+    @pytest.mark.parametrize("control_wires", cw5)
+    @pytest.mark.parametrize("auto", [False, True])
+    def test_decomposition_circuit(self, op, control_wires, auto, 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)
+            if auto:
+                ctrl_decomp_bisect(op, Wires(control_wires))
+            else:
+                record_from_list(_ctrl_decomp_bisect_general)(
+                    _convert_to_su2(op.matrix()), op.wires, 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)
+
+    @pytest.mark.parametrize("op", zip(gen_ops, gen_ops_best))
+    @pytest.mark.parametrize("control_wires", cw5)
+    @pytest.mark.parametrize("all_the_way_from_ctrl", [False, True])
+    def test_auto_select(self, op, control_wires, all_the_way_from_ctrl):
+        """
+        Test that the auto selection is correct and optimal.
+        """
+        op, best = op
+        if all_the_way_from_ctrl:
+            if isinstance(op, qml.PauliX):
+                # X has its own special case
+                pytest.skip()
+            res = qml.ctrl(op, control_wires).decomposition()
+        else:
+            res = ctrl_decomp_bisect(op, Wires(control_wires))
+        expected = best(_convert_to_su2(op.matrix()), op.wires, Wires(control_wires))
+        assert equal_list(res, expected)
+
+    @pytest.mark.parametrize("op", 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"""
+        assert np.allclose(op.matrix(), _convert_to_su2(op.matrix()), atol=tol, rtol=tol)
+
+        expected_op = qml.ctrl(op, control_wires)
+        res = qml.matrix(
+            record_from_list(_ctrl_decomp_bisect_general), wire_order=control_wires + [0]
+        )(op.matrix(), op.wires, control_wires)
+        expected = expected_op.matrix()
+
+        assert np.allclose(res, expected, atol=tol, rtol=tol)