Add support for arbitrary linear combination gradient recipes

.github/CHANGELOG.md

@@ -296,13 +296,22 @@
 
 <h3>Breaking changes</h3>
 
-* Updated how gradient recipes are stored for operations, allowing for
-  gradient recipe definitions involving custom multipliers and with more than
-  two terms.
+* Updated how parameter-shift gradient recipes are defined for operations, allowing for
+  gradient recipes that are specified as an arbitrary number of terms.
   [(#909)](https://github.com/PennyLaneAI/pennylane/pull/909)
 
-- The ``VQECost`` class has been renamed to ``ExpvalCost`` to reflect its general applicability
-  beyond VQE. Use of ``VQECost`` is still possible but will result in a deprecation warning.
+  Previously, `Operation.grad_recipe` was restricted to two-term parameter-shift formulas.
+  With this change, the gradient recipe now contains elements of the form
+  <img src="https://render.githubusercontent.com/render/math?math=[c_i, a_i, s_i]">, resulting in a gradient recipe of
+  <img src="https://render.githubusercontent.com/render/math?math=\frac{\partial}{\partial\phi_k}f(\phi_k) = \sum_{i} c_i * f(a_i * \phi_k+s_i)">.
+
+  As this is a breaking change, all custom operations with defined gradient recipes must be
+  updated to continue working with PennyLane 0.13. Note though that if `grad_recipe = None`, the
+  default gradient recipe remains unchanged, and corresponds to the two terms <img src="https://render.githubusercontent.com/render/math?math=[c_0, a_0, s_0]=[1/2, 1, \pi/2]">
+  and <img src="https://render.githubusercontent.com/render/math?math=[c_1, a_1, s_1]=[-1/2, 1, -\pi/2]"> for every parameter.
+
+- The `VQECost` class has been renamed to `ExpvalCost` to reflect its general applicability
+  beyond VQE. Use of `VQECost` is still possible but will result in a deprecation warning.
   [(#913)](https://github.com/PennyLaneAI/pennylane/pull/913)
 
 <h3>Documentation</h3>

doc/development/plugins.rst

@@ -469,17 +469,12 @@ where
 * :attr:`~.Operation.grad_recipe`: The gradient recipe for the analytic ``'A'``
   method. This is a tuple with one nested list per operation parameter. For
   parameter :math:`k`, the nested list contains elements of the form
-  :math:`[c_i, a_i, s_i]` where :math:`i \in I_{k}` is the index of the term,
-  resulting in a gradient recipe of
+  :math:`[c_i, a_i, s_i]`, resulting in a gradient recipe of
 
-  .. math:: \frac{\partial}{\partial\phi_k}O = \sum_{i \in I_{k}} c_i * O(a_i * \phi_k+s_i).
+  .. math:: \frac{\partial}{\partial\phi_k}f(\phi_k) = \sum_{i} c_i * f(a_i * \phi_k+s_i),
 
-  where :math:`f` is an expectation value that depends on :math:`O(\phi_k)`, an example being
-
-  .. math:: f(O(\phi_k)) = \braket{0 | O^{\dagger}(\phi_k) \hat{B} O(\phi_k) | 0}
-
-  which is the simple expectation value of the operator :math:`\hat{B}` evolved via the gate
-  :math:`O(\phi_k)`.
+  where :math:`f` is the expectation value of an observable on a circuit that has been evolved by
+  the operation being considered with parameter :math:`\phi_k`.
 
   Note that if ``grad_recipe = None``, the default gradient recipe containing
   the two terms :math:`[c_0, a_0, s_0]=[1/2, 1, \pi/2]` and :math:`[c_1, a_1,

pennylane/operation.py

@@ -65,17 +65,22 @@
     transformation on the quadrature operators.
 
 For gates that *are* supported via the analytic method, the gradient recipe
-(with multipliers :math:`c_i`, scaling factors :math:`a_i` and parameter shifts
-:math:`s_i` for parameter :math:`\phi_k` where :math:`i \in I_{k}`) works as
-follows:
+works as follows:
 
-.. math:: \frac{\partial}{\partial\phi_k}O = \sum_{i \in I_{k}} c_i * O(a_i * \phi_k+s_i).
+.. math:: \frac{\partial}{\partial\phi_k}f = \sum_{i \in I_{k}} c_i * f(a_i * \phi_k+s_i).
+
+where :math:`f` is the expectation value of an observable on a circuit that has
+been evolved by the operation being considered with parameter :math:`\phi_k`,
+there are multiple terms indexed with :math:`i` for each parameter :math:`\phi`
+and the :math:`[c_i, a_i, s_i]` are coefficients specific to the gate.
 
 The following specific case holds for example for qubit operations that are
 generated by one of the Pauli matrices and results in an overall positive and
 negative shift:
 
-.. math:: \frac{\partial}{\partial\phi_k}O = c_k\left[O(\phi_k+s_k)-O(\phi_k-s_k)\right].
+.. math:: \frac{\partial}{\partial\phi_k}f = \frac{1}{2}\left[f(\phi_k+\frac{\pi}{2})-f(\phi_k-\frac{\pi}{2})\right],
+
+i.e., so that :math:`[c_0, a_0, s_0]=[1/2, 1, \pi/2]` and :math:`[c_1, a_1, s_1]=[-1/2, 1, -\pi/2]`.
 
 CV Operation base classes
 ~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -652,7 +657,7 @@ def get_parameter_shift(self, idx, shift=np.pi / 2):
         a = 1
 
         # We set the default recipe following:
-        # ∂f(x) = c1*f(a1*x+s1) + c2*f(a2*x+s2)
+        # ∂f(x) = c*f(x+s) - c*f(x-s)
         # where we express a positive and a negative shift by default
         default_param_shift = [[multiplier, a, shift], [-multiplier, a, -shift]]
         param_shift = default_param_shift if recipe is None else recipe
@@ -1621,7 +1626,7 @@ def heisenberg_pd(self, idx):
         shift = np.pi / 2
 
         # We set the default recipe to as follows:
-        # ∂f(x) = c1*f(a1*x+s1) + c2*f(a2*x+s2)
+        # ∂f(x) = c*f(x+s) - c*f(x-s)
         default_param_shift = [[multiplier, a, shift], [-multiplier, a, -shift]]
         param_shift = default_param_shift if recipe is None else recipe