Repository: PennyLaneAI/pennylane
Pull request #1303 (2021-05-12 21:22:42)

[WIP] Adds decomposition to `DoubleExcitation` operator (and update gradient recipes)
**Context:** Currently `DoubleExcitation` is decomposed into `DoubleExcitationPlus` and `DoubleExcitationMinus`, neither of which are decomposed into elementary gates, such that `DoubleExcitation` cannot be implemented on all devices.

**Description of the Change:** Adds the decomposition into elementary gates based on [this paper](https://arxiv.org/abs/2104.05695) to `DoubleExcitation`. Furthermore, updates the gradient recipe of both `SingleExcitation` and `DoubleExcitation` to use the four-term parameter-shift rule of the aforementioned paper.

**Benefits:** `DoubleExcitation` can now run on any device; `Single` and `Double` excitation gradients are improved.

**Possible Drawbacks:** None

**Related GitHub Issues:** #1275 #1278 


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [.github/CHANGELOG.md](/PennyLaneAI/pennylane/blob/master/.github/CHANGELOG.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://codecov.io/gh/PennyLaneAI/pennylane/pull/1303?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1303](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1303?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (2d84e80) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/490f2d97a260834b6ffb4db6e99af39277fa08f4?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (490f2d9) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1303/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1303?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@           Coverage Diff           @@
##           master    #1303   +/-   ##
=======================================
  Coverage   98.15%   98.15%           
=======================================
  Files         147      147           
  Lines       11352    11354    +2     
=======================================
+ Hits        11142    11144    +2     
  Misses        210      210           
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1303?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/ops/qubit.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1303/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC5weQ==) | `98.36% <100.00%> (+<0.01%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1303?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1303?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [490f2d9...2d84e80](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1303?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index cc2f4ec7b..b4ad09c72 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -156,6 +156,11 @@ random_mat2 = rng.standard_normal(3, requires_grad=False)
 
 <h3>Improvements</h3>
 
+* The `qml.SingleExcitation` and `qml.DoubleExcitation` operations now
+  have decompositions over elementary gates, and their gradient recipes
+  have been updated to use the four-term parameter-shift rules.
+  [(#1303)](https://github.com/PennyLaneAI/pennylane/pull/1303)
+
 * The `qml.SingleExcitationPlus` and `qml.SingleExcitationMinus` operations now
   have decompositions over elementary gates.
   [(#1278)](https://github.com/PennyLaneAI/pennylane/pull/1278)
diff --git a/pennylane/ops/qubit.py b/pennylane/ops/qubit.py
index 9a5500918..9d981b604 100644
--- a/pennylane/ops/qubit.py
+++ b/pennylane/ops/qubit.py
@@ -1707,14 +1707,13 @@ class SingleExcitation(Operation):
 
     * Number of wires: 2
     * Number of parameters: 1
-    * Gradient recipe: Obtained from its decomposition in terms of the
-      :class:`~.SingleExcitationPlus` and :class:`~.SingleExcitationMinus` operations
+    * Gradient recipe: The ``SingleExcitation`` operator satisfies a four-term parameter-shift rule
+      (see Appendix F, https://arxiv.org/abs/2104.05695)
 
     Args:
         phi (float): rotation angle :math:`\phi`
         wires (Sequence[int]): the wires the operation acts on
 
-
     **Example**
 
     The following circuit performs the transformation :math:`|10\rangle\rightarrow \cos(
@@ -1737,6 +1736,7 @@ class SingleExcitation(Operation):
     num_wires = 2
     par_domain = "R"
     grad_method = "A"
+    grad_recipe = four_term_grad_recipe
     generator = [np.array([[0, 0, 0, 0], [0, 0, -1j, 0], [0, 1j, 0, 0], [0, 0, 0, 0]]), -1 / 2]
 
     @classmethod
@@ -1750,8 +1750,9 @@ class SingleExcitation(Operation):
     @staticmethod
     def decomposition(theta, wires):
         decomp_ops = [
-            SingleExcitationPlus(theta / 2, wires=wires),
-            SingleExcitationMinus(theta / 2, wires=wires),
+            qml.CNOT(wires=[wires[0], wires[1]]),
+            qml.CRY(theta, wires=[wires[1], wires[0]]),
+            qml.CNOT(wires=[wires[0], wires[1]]),
         ]
         return decomp_ops
 
@@ -2241,8 +2242,8 @@ class DoubleExcitation(Operation):
 
     * Number of wires: 4
     * Number of parameters: 1
-    * Gradient recipe: Obtained from its decomposition in terms of the
-      :class:`~.DoubleExcitationPlus` and :class:`~.DoubleExcitationMinus` operations
+    * Gradient recipe: The ``DoubleExcitation`` operator satisfies a four-term parameter-shift rule
+      (see Appendix F, https://arxiv.org/abs/2104.05695):
 
     Args:
         phi (float): rotation angle :math:`\phi`
@@ -2271,6 +2272,7 @@ class DoubleExcitation(Operation):
     num_wires = 4
     par_domain = "R"
     grad_method = "A"
+    grad_recipe = four_term_grad_recipe
 
     G = np.zeros((16, 16), dtype=np.complex64)
     G[3, 12] = -1j  # 3 (dec) = 0011 (bin)
@@ -2293,10 +2295,38 @@ class DoubleExcitation(Operation):
 
     @staticmethod
     def decomposition(theta, wires):
+        # This decomposition is the "upside down" version of that on p17 of https://arxiv.org/abs/2104.05695
         decomp_ops = [
-            DoubleExcitationPlus(theta / 2, wires=wires),
-            DoubleExcitationMinus(theta / 2, wires=wires),
+            qml.CNOT(wires=[wires[2], wires[3]]),
+            qml.CNOT(wires=[wires[0], wires[2]]),
+            qml.Hadamard(wires=wires[3]),
+            qml.Hadamard(wires=wires[0]),
+            qml.CNOT(wires=[wires[2], wires[3]]),
+            qml.CNOT(wires=[wires[0], wires[1]]),
+            qml.RY(theta / 8, wires=wires[1]),
+            qml.RY(-theta / 8, wires=wires[0]),
+            qml.CNOT(wires=[wires[0], wires[3]]),
+            qml.Hadamard(wires=wires[3]),
+            qml.CNOT(wires=[wires[3], wires[1]]),
+            qml.RY(theta / 8, wires=wires[1]),
+            qml.RY(-theta / 8, wires=wires[0]),
+            qml.CNOT(wires=[wires[2], wires[1]]),
+            qml.CNOT(wires=[wires[2], wires[0]]),
+            qml.RY(-theta / 8, wires=wires[1]),
+            qml.RY(theta / 8, wires=wires[0]),
+            qml.CNOT(wires=[wires[3], wires[1]]),
+            qml.Hadamard(wires=wires[3]),
+            qml.CNOT(wires=[wires[0], wires[3]]),
+            qml.RY(-theta / 8, wires=wires[1]),
+            qml.RY(theta / 8, wires=wires[0]),
+            qml.CNOT(wires=[wires[0], wires[1]]),
+            qml.CNOT(wires=[wires[2], wires[0]]),
+            qml.Hadamard(wires=wires[0]),
+            qml.Hadamard(wires=wires[3]),
+            qml.CNOT(wires=[wires[0], wires[2]]),
+            qml.CNOT(wires=[wires[2], wires[3]]),
         ]
+
         return decomp_ops
 
     def adjoint(self):


Test patch:
diff --git a/tests/gate_data.py b/tests/gate_data.py
index fe72a9e5b..d7119c7ee 100644
--- a/tests/gate_data.py
+++ b/tests/gate_data.py
@@ -14,6 +14,11 @@ Y = np.array([[0, -1j], [1j, 0]])  #: Pauli-Y matrix
 Z = np.array([[1, 0], [0, -1]])  #: Pauli-Z matrix
 
 H = np.array([[1, 1], [1, -1]]) / math.sqrt(2)  #: Hadamard gate
+
+# Single-qubit projectors
+StateZeroProjector = np.array([[1, 0], [0, 0]])
+StateOneProjector = np.array([[0, 0], [0, 1]])
+
 # Two qubit gates
 CNOT = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])  #: CNOT gate
 SWAP = np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]])  #: SWAP gate
diff --git a/tests/ops/test_qubit_ops.py b/tests/ops/test_qubit_ops.py
index 165a485ae..f7d017b97 100644
--- a/tests/ops/test_qubit_ops.py
+++ b/tests/ops/test_qubit_ops.py
@@ -30,6 +30,8 @@ from gate_data import (
     Y,
     Z,
     H,
+    StateZeroProjector,
+    StateOneProjector,
     CNOT,
     SWAP,
     CZ,
@@ -987,12 +989,30 @@ class TestSingleExcitation:
 
     @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4])
     def test_single_excitation_decomp(self, phi):
-        """Tests that the SingleExcitation operation calculates the correct decomposition"""
-        op = qml.SingleExcitation(phi, wires=[0, 1])
-        decomp = op.decomposition(phi, wires=[0, 1])
-        mats = [m.matrix for m in decomp]
-        decomposed_matrix = mats[0] @ mats[1]
+        """Tests that the SingleExcitation operation calculates the correct decomposition.
+
+        Need to consider the matrix of CRY separately, as the control is wire 1
+        and the target is wire 0 in the decomposition."""
+        decomp = qml.SingleExcitation.decomposition(phi, wires=[0, 1])
+
+        mats = []
+        for i in reversed(decomp):
+            if i.wires.tolist() == [1, 0] and isinstance(i, qml.CRY):
+                new_mat = np.array(
+                    [
+                        [1, 0, 0, 0],
+                        [0, np.cos(phi / 2), 0, -np.sin(phi / 2)],
+                        [0, 0, 1, 0],
+                        [0, np.sin(phi / 2), 0, np.cos(phi / 2)],
+                    ]
+                )
+                mats.append(new_mat)
+            else:
+                mats.append(i.matrix)
+
+        decomposed_matrix = np.linalg.multi_dot(mats)
         exp = SingleExcitation(phi)
+
         assert np.allclose(decomposed_matrix, exp)
 
     @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4])
@@ -1247,6 +1267,53 @@ class TestDoubleExcitation:
 
         assert np.allclose(res, exp)
 
+    @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
+    def test_double_excitation_decomp(self, phi):
+        """Tests that the DoubleExcitation operation calculates the correct decomposition.
+
+        The decomposition has already been expressed in terms of single-qubit rotations
+        and CNOTs. For each term in the decomposition we need to construct the appropriate
+        four-qubit tensor product matrix and then multiply them together.
+        """
+        decomp = qml.DoubleExcitation.decomposition(phi, wires=[0, 1, 2, 3])
+
+        from functools import reduce
+
+        # To compute the matrix for CX on an arbitrary number of qubits, use the fact that
+        # CU  = |0><0| \otimes I + |1><1| \otimes U
+        def cnot_four_qubits(wires):
+            proj_0_term = [StateZeroProjector if idx == wires[0] else np.eye(2) for idx in range(4)]
+
+            proj_1_term = [np.eye(2) for idx in range(4)]
+            proj_1_term[wires[0]] = StateOneProjector
+            proj_1_term[wires[1]] = X
+
+            proj_0_kron = reduce(np.kron, proj_0_term)
+            proj_1_kron = reduce(np.kron, proj_1_term)
+
+            return proj_0_kron + proj_1_kron
+
+        # Inserts a single-qubit matrix into a four-qubit matrix at the right place
+        def single_mat_four_qubits(mat, wire):
+            individual_mats = [mat if idx == wire else np.eye(2) for idx in range(4)]
+            return reduce(np.kron, individual_mats)
+
+        mats = []
+        for i in reversed(decomp):
+            # Single-qubit gate
+            if len(i.wires.tolist()) == 1:
+                mat = single_mat_four_qubits(i.matrix, i.wires.tolist()[0])
+                mats.append(mat)
+            # Two-qubit gate
+            else:
+                mat = cnot_four_qubits(i.wires.tolist())
+                mats.append(mat)
+
+        decomposed_matrix = np.linalg.multi_dot(mats)
+        exp = DoubleExcitation(phi)
+
+        assert np.allclose(decomposed_matrix, exp)
+
     @pytest.mark.parametrize("phi", [-0.1, 0.2, np.pi / 4])
     def test_double_excitation_plus_matrix(self, phi):
         """Tests that the DoubleExcitationPlus operation calculates the correct matrix"""