Add support for multi-controlled zyz

Makefile

@@ -70,17 +70,17 @@ coverage:
 .PHONY:format
 format:
 ifdef check
-	isort --py 311 --profile black -l 100 -o autoray -p ./pennylane --skip __init__.py --filter-files ./pennylane ./tests --check
-	black -t py39 -t py310 -t py311 -l 100 ./pennylane ./tests --check
+	$(PYTHON) -m isort --py 311 --profile black -l 100 -o autoray -p ./pennylane --skip __init__.py --filter-files ./pennylane ./tests --check
+	$(PYTHON) -m black -t py39 -t py310 -t py311 -l 100 ./pennylane ./tests --check
 else
-	isort --py 311 --profile black -l 100 -o autoray -p ./pennylane --skip __init__.py --filter-files ./pennylane ./tests
-	black -t py39 -t py310 -t py311 -l 100 ./pennylane ./tests
+	$(PYTHON) -m isort --py 311 --profile black -l 100 -o autoray -p ./pennylane --skip __init__.py --filter-files ./pennylane ./tests
+	$(PYTHON) -m black -t py39 -t py310 -t py311 -l 100 ./pennylane ./tests
 endif
 
 .PHONY: lint
 lint:
-	pylint pennylane --rcfile .pylintrc
+	$(PYTHON) -m pylint pennylane --rcfile .pylintrc
 
 .PHONY: lint-test
 lint-test:
-	pylint tests pennylane/devices/tests --rcfile tests/.pylintrc
+	$(PYTHON) -m pylint tests pennylane/devices/tests --rcfile tests/.pylintrc

doc/releases/changelog-dev.md

@@ -216,6 +216,9 @@
 
 <h4>Other improvements</h4>
 
+* Added the decomposition of zyz for special unitaries with multiple control wires.
+  [(#6042)](https://github.com/PennyLaneAI/pennylane/pull/6042)
+
 * A new method `process_density_matrix` has been added to the `ProbabilityMP` and `DensityMatrixMP`
   classes, allowing for more efficient handling of quantum density matrices, particularly with batch
   processing support. This method simplifies the calculation of probabilities from quantum states
@@ -395,6 +398,7 @@ This release contains contributions from (in alphabetical order):
 
 Tarun Kumar Allamsetty,
 Guillermo Alonso,
+Ali Asadi,
 Utkarsh Azad,
 Tonmoy T. Bhattacharya,
 Gabriel Bottrill,

pennylane/ops/op_math/controlled.py

@@ -864,7 +864,7 @@ def _decompose_custom_ops(op: Controlled) -> list["operation.Operator"]:
     return None
 
 
-def _decompose_no_control_values(op: Controlled) -> list["operation.Operator"]:
+def _decompose_no_control_values(op: Controlled) -> Optional[list["operation.Operator"]]:
     """Decompose without considering control values. Returns None if no decomposition."""
 
     decomp = _decompose_custom_ops(op)
@@ -874,7 +874,7 @@ def _decompose_no_control_values(op: Controlled) -> list["operation.Operator"]:
     if _is_single_qubit_special_unitary(op.base):
         if len(op.control_wires) >= 2 and qmlmath.get_interface(*op.data) == "numpy":
             return ctrl_decomp_bisect(op.base, op.control_wires)
-        return ctrl_decomp_zyz(op.base, op.control_wires)
+        return ctrl_decomp_zyz(op.base, op.control_wires, work_wires=op.work_wires)
 
     if not op.base.has_decomposition:
         return None

pennylane/ops/op_math/controlled_decompositions.py

@@ -16,6 +16,7 @@
 """
 
 from copy import copy
+from typing import Optional
 
 import numpy as np
 import numpy.linalg as npl
@@ -134,12 +135,95 @@ def _bisect_compute_b(u: np.ndarray):
     return _param_su2(c, d, b, 0)
 
 
-def ctrl_decomp_zyz(target_operation: Operator, control_wires: Wires):
+def _multi_controlled_zyz(
+    rot_angles,
+    global_phase,
+    target_wire: Wires,
+    control_wires: Wires,
+    work_wires: Optional[Wires] = None,
+) -> list[Operator]:
+    # The decomposition of zyz for special unitaries with multiple control wires
+    # defined in Lemma 7.9 of https://arxiv.org/pdf/quant-ph/9503016
+
+    if not qml.math.allclose(0.0, global_phase, atol=1e-6, rtol=0):
+        raise ValueError(f"The global_phase should be zero, instead got: {global_phase}.")
+
+    # Unpack the rotation angles
+    phi, theta, omega = rot_angles
+
+    # We use the conditional statements to account when decomposition is ran within a queue
+    decomp = []
+
+    cop_wires = (control_wires[-1], target_wire[0])
+
+    # Add A operator
+    if not qml.math.allclose(0.0, phi, atol=1e-8, rtol=0):
+        decomp.append(qml.CRZ(phi, wires=cop_wires))
+    if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.CRY(theta / 2, wires=cop_wires))
+
+    decomp.append(qml.ctrl(qml.X(target_wire), control=control_wires[:-1], work_wires=work_wires))
+
+    # Add B operator
+    if not qml.math.allclose(0.0, theta / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.CRY(-theta / 2, wires=cop_wires))
+    if not qml.math.allclose(0.0, -(phi + omega) / 2, atol=1e-6, rtol=0):
+        decomp.append(qml.CRZ(-(phi + omega) / 2, wires=cop_wires))
+
+    decomp.append(qml.ctrl(qml.X(target_wire), control=control_wires[:-1], work_wires=work_wires))
+
+    # Add C operator
+    if not qml.math.allclose(0.0, (omega - phi) / 2, atol=1e-8, rtol=0):
+        decomp.append(qml.CRZ((omega - phi) / 2, wires=cop_wires))
+
+    return decomp
+
+
+def _single_control_zyz(rot_angles, global_phase, target_wire, control_wires: Wires):
+    # The zyz decomposition of a general unitary with single control wire
+    # defined in Lemma 7.9 of https://arxiv.org/pdf/quant-ph/9503016
+
+    # Unpack the rotation angles
+    phi, theta, omega = rot_angles
+    # 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.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.X(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))
+
+    return decomp
+
+
+def ctrl_decomp_zyz(
+    target_operation: Operator, control_wires: Wires, work_wires: Optional[Wires] = None
+) -> list[Operator]:
     """Decompose the controlled version of a target single-qubit operation
 
-    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>`_.
+    This function decomposes both single and multiple controlled single-qubit
+    target operations using the decomposition defined in Lemma 4.3 and Lemma 5.1
+    for single `controlled_wires`, and Lemma 7.9 for multiple `controlled_wires`
+    from `Barenco et al. (1995) <https://arxiv.org/abs/quant-ph/9503016>`_.
 
     Args:
         target_operation (~.operation.Operator): the target operation to decompose
@@ -190,53 +274,24 @@ def decomp_circuit(op):
             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
 
     if isinstance(target_operation, Operation):
         try:
-            phi, theta, omega = target_operation.single_qubit_rot_angles()
+            rot_angles = target_operation.single_qubit_rot_angles()
         except NotImplementedError:
-            phi, theta, omega = _get_single_qubit_rot_angles_via_matrix(
-                qml.matrix(target_operation)
-            )
+            rot_angles = _get_single_qubit_rot_angles_via_matrix(qml.matrix(target_operation))
     else:
-        phi, theta, omega = _get_single_qubit_rot_angles_via_matrix(qml.matrix(target_operation))
+        rot_angles = _get_single_qubit_rot_angles_via_matrix(qml.matrix(target_operation))
 
     _, 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.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))
-
-    return decomp
+    return (
+        _multi_controlled_zyz(rot_angles, global_phase, target_wire, control_wires, work_wires)
+        if len(control_wires) > 1
+        else _single_control_zyz(rot_angles, global_phase, target_wire, control_wires)
+    )
 
 
 def _ctrl_decomp_bisect_od(

tests/ops/op_math/test_controlled.py

@@ -1050,20 +1050,39 @@ def test_non_differentiable_one_qubit_special_unitary(self):
         decomp_mat = qml.matrix(op.decomposition, wire_order=op.wires)()
         assert qml.math.allclose(op.matrix(), decomp_mat)
 
-    def test_differentiable_one_qubit_special_unitary(self):
-        """Assert that a differentiable qubit special unitary uses the zyz decomposition."""
+    def test_differentiable_one_qubit_special_unitary_single_ctrl(self):
+        """
+        Assert that a differentiable qubit special unitary uses the zyz decomposition with a single controlled wire.
+        """
 
-        op = qml.ctrl(qml.RZ(qml.numpy.array(1.2), 0), (1))
+        theta = 1.2
+        op = qml.ctrl(qml.RZ(qml.numpy.array(theta), 0), (1))
         decomp = op.decomposition()
 
-        qml.assert_equal(decomp[0], qml.PhaseShift(qml.numpy.array(1.2 / 2), 0))
+        qml.assert_equal(decomp[0], qml.PhaseShift(qml.numpy.array(theta / 2), 0))
         qml.assert_equal(decomp[1], qml.CNOT(wires=(1, 0)))
-        qml.assert_equal(decomp[2], qml.PhaseShift(qml.numpy.array(-1.2 / 2), 0))
+        qml.assert_equal(decomp[2], qml.PhaseShift(qml.numpy.array(-theta / 2), 0))
         qml.assert_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)
 
+    def test_differentiable_one_qubit_special_unitary_multiple_ctrl(self):
+        """Assert that a differentiable qubit special unitary uses the zyz decomposition with multiple controlled wires."""
+
+        theta = 1.2
+        op = qml.ctrl(qml.RZ(qml.numpy.array(theta), 0), (1, 2, 3, 4))
+        decomp = op.decomposition()
+
+        assert qml.equal(decomp[0], qml.CRZ(qml.numpy.array(theta), [4, 0]))
+        assert qml.equal(decomp[1], qml.MultiControlledX(wires=[1, 2, 3, 0]))
+        assert qml.equal(decomp[2], qml.CRZ(qml.numpy.array(-theta / 2), wires=[4, 0]))
+        assert qml.equal(decomp[3], qml.MultiControlledX(wires=[1, 2, 3, 0]))
+        assert qml.equal(decomp[4], qml.CRZ(qml.numpy.array(-theta / 2), wires=[4, 0]))
+
+        decomp_mat = qml.matrix(op.decomposition, wire_order=op.wires)()
+        assert qml.math.allclose(op.matrix(), decomp_mat)
+
     @pytest.mark.parametrize(
         "base_cls, params, base_wires, ctrl_wires, custom_ctrl_cls, expected",
         custom_ctrl_op_decomps,
@@ -1730,7 +1749,6 @@ def test_custom_controlled_ops_wrong_wires(self, op, ctrl_wires, _):
 
         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,