Repository: PennyLaneAI/pennylane
Pull request #2208 (2022-02-17 09:19:08)

Use par frequencies to master
This PR pulls the changes from #2182 that were accidentally merged into `op-refactor` through to `master`.

codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#2208](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (6044a75) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/88a7b36eba72b6314fc045fc288e75f0c1019b40?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (88a7b36) will **decrease** coverage by `2.75%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/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/2208?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@            Coverage Diff             @@
##           master    #2208      +/-   ##
==========================================
- Coverage   99.25%   96.50%   -2.76%     
==========================================
  Files         231      231              
  Lines       18093    18116      +23     
==========================================
- Hits        17959    17483     -476     
- Misses        134      633     +499     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/gradients/general\_shift\_rules.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyYWRpZW50cy9nZW5lcmFsX3NoaWZ0X3J1bGVzLnB5) | `100.00% <100.00%> (ø)` | |
| [pennylane/gradients/parameter\_shift.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyYWRpZW50cy9wYXJhbWV0ZXJfc2hpZnQucHk=) | `100.00% <100.00%> (ø)` | |
| [pennylane/gradients/parameter\_shift\_cv.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyYWRpZW50cy9wYXJhbWV0ZXJfc2hpZnRfY3YucHk=) | `99.15% <100.00%> (+0.01%)` | :arrow_up: |
| [pennylane/operation.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wZXJhdGlvbi5weQ==) | `95.91% <100.00%> (+0.05%)` | :arrow_up: |
| [pennylane/interfaces/batch/jax.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvYmF0Y2gvamF4LnB5) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/interfaces/batch/jax\_jit.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvYmF0Y2gvamF4X2ppdC5weQ==) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [pennylane/interfaces/jax.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2ludGVyZmFjZXMvamF4LnB5) | `0.00% <0.00%> (-97.44%)` | :arrow_down: |
| [pennylane/devices/default\_qubit\_jax.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2RldmljZXMvZGVmYXVsdF9xdWJpdF9qYXgucHk=) | `0.00% <0.00%> (-96.63%)` | :arrow_down: |
| [pennylane/qnn/cost.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3Fubi9jb3N0LnB5) | `42.85% <0.00%> (-57.15%)` | :arrow_down: |
| [pennylane/transforms/classical\_jacobian.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvY2xhc3NpY2FsX2phY29iaWFuLnB5) | `77.08% <0.00%> (-22.92%)` | :arrow_down: |
| ... and [17 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208?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/2208?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [88a7b36...6044a75](https://codecov.io/gh/PennyLaneAI/pennylane/pull/2208?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/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 3c3926b94..61d42cdc8 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -24,8 +24,10 @@
 
 * Parametric operations now have the `parameter_frequencies`
   method that returns the frequencies with which a parameter
-  enters a circuit when using the operation.
+  enters a circuit. In addition, the general parameter-shift
+  rule is now automatically used by `qml.gradients.param_shift`.
   [(#2180)](https://github.com/PennyLaneAI/pennylane/pull/2180)
+  [(#2182)](https://github.com/PennyLaneAI/pennylane/pull/2182)
 
   The frequencies can be used for circuit analysis, optimization
   via the `RotosolveOptimizer` and differentiation with the
@@ -33,6 +35,23 @@
   expectation values or probabilities, for a variance
   measurement the frequencies will differ.
 
+  By default, the frequencies will be obtained from the
+  `generator` property (if it is defined).
+
+  When using `qml.gradients.param_shift`, the parameter frequencies
+  are used to obtain the shift rule for the operation.
+
+  For operations that are registered to have an analytic gradient
+  method but that do not provide parameter frequencies, the
+  `grad_recipe` of the operation will be used for differentiation
+  instead. If there is no `grad_recipe`, the standard two-term shift
+  rule will be used.
+
+  See [Vidal and Theis (2018)](https://arxiv.org/abs/1812.06323)
+  and [Wierichs et al. (2021)](https://arxiv.org/abs/2107.12390)
+  for theoretical background information on the general
+  parameter-shift rule.
+
 * Continued development of the circuit-cutting compiler:
 
   A method for converting a quantum tape to a directed multigraph that is amenable
diff --git a/pennylane/gradients/general_shift_rules.py b/pennylane/gradients/general_shift_rules.py
index efcfdca54..fc06ae21f 100644
--- a/pennylane/gradients/general_shift_rules.py
+++ b/pennylane/gradients/general_shift_rules.py
@@ -20,22 +20,25 @@ import numpy as np
 import pennylane as qml
 
 
-def _process_shifts(rule, tol=1e-10):
+def process_shifts(rule, tol=1e-10, check_duplicates=True):
     """Utility function to process gradient rules.
 
     Args:
-        rule (array): a ``(N, M)`` array corresponding to ``M`` terms
-            with ``N-1`` simultaneous function shifts. The first row of
+        rule (array): a ``(2, M)`` array corresponding to ``M`` terms
+            with function shifts. The first row of
             the array corresponds to the linear combination coefficients;
-            subsequent rows correspond to parameter shifts for these coefficients.
+            the second row contains the shift values.
         tol (float): floating point tolerance used when comparing shifts/coefficients
+            Terms with coefficients below ``tol`` will be removed.
+        check_duplicates (bool): whether to check the input ``rule`` for duplicate
+            shift values in its second row.
 
     This utility function accepts coefficients and shift values, and performs the following
     processing:
 
-    - Removes all small (within absolute tolerance ``tol``) coefficients and shifts
+    - Sets all small (within absolute tolerance ``tol``) coefficients and shifts to 0
 
-    - Removes terms with the coefficients are 0
+    - Removes terms with the coefficients are 0 (including the ones set to 0 in the previous step)
 
     - Terms with the same shift value are combined into a single term.
 
@@ -48,25 +51,21 @@ def _process_shifts(rule, tol=1e-10):
     # remove columns where the coefficients are 0
     rule = rule[:, ~(rule[0] == 0)]
 
-    # sort columns according to abs(shift)
-    rule = rule[:, np.argsort(np.abs(rule)[-1])]
-
-    # determine unique shifts
-    round_decimals = int(-np.log10(tol))
-    rounded_rule = np.round(rule[-1], round_decimals)
-    unique_shifts = np.unique(rounded_rule)
+    if check_duplicates:
+        # determine unique shifts
+        round_decimals = int(-np.log10(tol))
+        rounded_rule = np.round(rule[-1], round_decimals)
+        unique_shifts = np.unique(rounded_rule)
 
-    if rule.shape[-1] != len(unique_shifts):
-        # sum columns that have the same shift value
-        coeffs = [
-            np.sum(rule[:, np.nonzero(rounded_rule == s)[0]], axis=1)[0] for s in unique_shifts
-        ]
-        rule = np.stack([np.stack(coeffs), unique_shifts])
+        if rule.shape[-1] != len(unique_shifts):
+            # sum columns that have the same shift value
+            coeffs = [
+                np.sum(rule[:, np.nonzero(rounded_rule == s)[0]], axis=1)[0] for s in unique_shifts
+            ]
+            rule = np.stack([np.stack(coeffs), unique_shifts])
 
-        # sort columns according to abs(shift)
-        rule = rule[:, np.argsort(np.abs(rule)[-1])]
-
-    return rule
+    # sort columns according to abs(shift)
+    return rule[:, np.argsort(np.abs(rule)[-1])]
 
 
 @functools.lru_cache(maxsize=None)
@@ -255,7 +254,7 @@ def generate_shift_rule(frequencies, shifts=None, order=1):
         # array, with coefficients on the first row and shifts on the second row.
         rule = qml.math.stack(combined_rules).T
 
-    return _process_shifts(rule, tol=1e-10)
+    return process_shifts(rule, tol=1e-10)
 
 
 def generate_multi_shift_rule(frequencies, shifts=None, orders=None):
@@ -313,7 +312,7 @@ def generate_multi_shift_rule(frequencies, shifts=None, orders=None):
 
     for f, s, o in zip(frequencies, shifts, orders):
         rule = generate_shift_rule(f, shifts=s, order=o)
-        rules.append(_process_shifts(rule).T)
+        rules.append(process_shifts(rule).T)
 
     combined_rules = []
 
diff --git a/pennylane/gradients/parameter_shift.py b/pennylane/gradients/parameter_shift.py
index 4db333ea2..419b03547 100644
--- a/pennylane/gradients/parameter_shift.py
+++ b/pennylane/gradients/parameter_shift.py
@@ -23,6 +23,7 @@ import pennylane as qml
 
 from .gradient_transform import gradient_transform
 from .finite_difference import finite_diff, generate_shifted_tapes
+from .general_shift_rules import process_shifts
 
 
 NONINVOLUTORY_OBS = {
@@ -61,27 +62,49 @@ def _square_observable(obs):
     return NONINVOLUTORY_OBS[obs.name](obs)
 
 
-def _get_operation_recipe(tape, t_idx, shift=np.pi / 2):
+def _get_operation_recipe(tape, t_idx, shifts):
     """Utility function to return the parameter-shift rule
     of the operation corresponding to trainable parameter
     t_idx on tape.
 
-    If the corresponding operation has grad_recipe=None, then
-    the default two-term parameter-shift rule is assumed.
+    This function performs multiple attempts to obtain the recipe:
+
+    - If ``parameter_frequencies`` yield a result, the frequencies are
+      used to construct the general parameter-shift rule via
+      ``qml.gradients.generate_shift_rule``
+
+    - If ``parameter_frequencies`` raises an error (because the operation
+      has no custom ``parameter_frequencies`` and no ``generator`` is defined),
+      the :meth:`.Operator.get_parameter_shift` method is used.
+      If in turn no custom version of ``get_parameter_shift`` is defined
+      and the operation does not have a ``grad_recipe``, the two-term
+      parameter-shift rule is assumed.
+
+    That is, a default to the two-term rule only is returned if no custom
+    ``parameter_frequencies``, ``generator``, ``get_parameter_shift``,
+    and ``grad_recipe`` are defined.
     """
     op, p_idx = tape.get_operation(t_idx)
-    return op.get_parameter_shift(p_idx, shift=shift)
+    try:
+        # Obtain frequencies, either via custom implementation or from generator eigvals
+        frequencies = op.parameter_frequencies[p_idx]
+        # Create shift rule from frequencies with given shifts
+        coeffs, shifts = qml.gradients.generate_shift_rule(frequencies, shifts=shifts, order=1)
+        # The shift rules do not include a rescaling of the parameter, only shifts.
+        mults = np.ones_like(coeffs)
+        return coeffs, mults, shifts
+
+    except qml.operation.OperatorPropertyUndefined:
+        # if no frequencies could be obtained, use the operation's shift rule if present.
+        # Note that the keyword argument ``shift`` only is used for the default two-term rule
+        if shifts is None:
+            use_shift = None
+        else:
+            use_shift = shifts[0]
 
+        recipe = np.array(op.get_parameter_shift(p_idx, shift=use_shift)).T
 
-def _process_gradient_recipe(gradient_recipe, tol=1e-10):
-    """Utility function to process gradient recipes."""
-    gradient_recipe = np.array(gradient_recipe).T
-    # remove all small coefficients, shifts, and multipliers
-    gradient_recipe[np.abs(gradient_recipe) < tol] = 0
-    # remove columns where the coefficients are 0
-    gradient_recipe = gradient_recipe[:, ~(gradient_recipe[0] == 0)]
-    # sort columns according to abs(shift)
-    return gradient_recipe[:, np.argsort(np.abs(gradient_recipe)[-1])]
+        return process_shifts(recipe, check_duplicates=False)
 
 
 def _gradient_analysis(tape, use_graph=True):
@@ -97,19 +120,25 @@ def _gradient_analysis(tape, use_graph=True):
     for idx, info in tape._par_info.items():
 
         if idx not in tape.trainable_params:
+            # non-trainable parameters do not require a grad_method
             info["grad_method"] = None
         else:
             op = tape._par_info[idx]["op"]
 
-            if op.grad_method == "F":
-                info["grad_method"] = "F"
-            else:
-                info["grad_method"] = tape._grad_method(
-                    idx, use_graph=use_graph, default_method="A"
-                )
+            if not qml.operation.has_grad_method(op):
+                # no differentiation method is registered for this operation
+                info["grad_method"] = None
 
+            elif (tape._graph is not None) or use_graph:
+                if not any(tape.graph.has_path(op, ob) for ob in tape.observables):
+                    # there is no influence of this operation on any of the observables
+                    info["grad_method"] = "0"
+                    continue
 
-def expval_param_shift(tape, argnum=None, shift=np.pi / 2, gradient_recipes=None, f0=None):
+            info["grad_method"] = op.grad_method
+
+
+def expval_param_shift(tape, argnum=None, shifts=None, gradient_recipes=None, f0=None):
     r"""Generate the parameter-shift tapes and postprocessing methods required
     to compute the gradient of a gate parameter with respect to an
     expectation value.
@@ -119,10 +148,10 @@ def expval_param_shift(tape, argnum=None, shift=np.pi / 2, gradient_recipes=None
         argnum (int or list[int] or None): Trainable parameter indices to differentiate
             with respect to. If not provided, the derivatives with respect to all
             trainable indices are returned.
-        shift (float): The shift value to use for the two-term parameter-shift formula.
-            Only valid if the operation in question supports the two-term parameter-shift
-            rule (that is, it has two distinct eigenvalues) and ``gradient_recipes``
-            is ``None``.
+        shifts (list[tuple[int or float]]): List containing tuples of shift values.
+            If provided, one tuple of shifts should be given per trainable parameter
+            and the tuple should match the number of frequencies for that parameter.
+            If unspecified, equidistant shifts are assumed.
         gradient_recipes (tuple(list[list[float]] or None)): List of gradient recipes
             for the parameter-shift method. One gradient recipe must be provided
             per trainable parameter.
@@ -171,13 +200,18 @@ def expval_param_shift(tape, argnum=None, shift=np.pi / 2, gradient_recipes=None
             continue
 
         # get the gradient recipe for the trainable parameter
-        recipe = gradient_recipes[argnum.index(idx)]
-        recipe = recipe or _get_operation_recipe(tape, idx, shift=shift)
-        recipe = _process_gradient_recipe(recipe)
-        coeffs, multipliers, shifts = recipe
+        arg_idx = argnum.index(idx)
+        recipe = gradient_recipes[arg_idx]
+        if recipe is not None:
+            recipe = process_shifts(np.array(recipe).T)
+        else:
+            op_shifts = None if shifts is None else shifts[arg_idx]
+            recipe = _get_operation_recipe(tape, idx, shifts=op_shifts)
+        coeffs, multipliers, op_shifts = recipe
         fns.append(None)
 
-        if shifts[0] == 0 and multipliers[0] == 1:
+        # Extract zero-shift term if present (if so, it will always be the first)
+        if op_shifts[0] == 0 and multipliers[0] == 1:
             # Gradient recipe includes a term with zero shift.
 
             if not unshifted_coeffs and f0 is None:
@@ -188,11 +222,11 @@ def expval_param_shift(tape, argnum=None, shift=np.pi / 2, gradient_recipes=None
             # Store the unshifted coefficient. We know that
             # it will always be the first coefficient due to processing.
             unshifted_coeffs.append(coeffs[0])
-            coeffs, multipliers, shifts = recipe[:, 1:]
+            coeffs, multipliers, op_shifts = recipe[:, 1:]
 
         # generate the gradient tapes
         gradient_coeffs.append(coeffs)
-        g_tapes = generate_shifted_tapes(tape, idx, shifts, multipliers)
+        g_tapes = generate_shifted_tapes(tape, idx, op_shifts, multipliers)
 
         gradient_tapes.extend(g_tapes)
         shapes.append(len(g_tapes))
@@ -241,7 +275,7 @@ def expval_param_shift(tape, argnum=None, shift=np.pi / 2, gradient_recipes=None
     return gradient_tapes, processing_fn
 
 
-def var_param_shift(tape, argnum, shift=np.pi / 2, gradient_recipes=None, f0=None):
+def var_param_shift(tape, argnum, shifts=None, gradient_recipes=None, f0=None):
     r"""Generate the parameter-shift tapes and postprocessing methods required
     to compute the gradient of a gate parameter with respect to a
     variance value.
@@ -251,10 +285,10 @@ def var_param_shift(tape, argnum, shift=np.pi / 2, gradient_recipes=None, f0=Non
         argnum (int or list[int] or None): Trainable parameter indices to differentiate
             with respect to. If not provided, the derivative with respect to all
             trainable indices are returned.
-        shift (float): The shift value to use for the two-term parameter-shift formula.
-            Only valid if the operation in question supports the two-term parameter-shift
-            rule (that is, it has two distinct eigenvalues) and ``gradient_recipes``
-            is ``None``.
+        shifts (list[tuple[int or float]]): List containing tuples of shift values.
+            If provided, one tuple of shifts should be given per trainable parameter
+            and the tuple should match the number of frequencies for that parameter.
+            If unspecified, equidistant shifts are assumed.
         gradient_recipes (tuple(list[list[float]] or None)): List of gradient recipes
             for the parameter-shift method. One gradient recipe must be provided
             per trainable parameter.
@@ -286,7 +320,7 @@ def var_param_shift(tape, argnum, shift=np.pi / 2, gradient_recipes=None, f0=Non
     gradient_tapes = [expval_tape]
 
     # evaluate the analytic derivative of <A>
-    pdA_tapes, pdA_fn = expval_param_shift(expval_tape, argnum, shift, gradient_recipes, f0)
+    pdA_tapes, pdA_fn = expval_param_shift(expval_tape, argnum, shifts, gradient_recipes, f0)
     gradient_tapes.extend(pdA_tapes)
 
     # Store the number of first derivative tapes, so that we know
@@ -327,7 +361,7 @@ def var_param_shift(tape, argnum, shift=np.pi / 2, gradient_recipes=None, f0=Non
         # may be non-zero. Here, we calculate the analytic derivatives of the <A^2>
         # observables.
         pdA2_tapes, pdA2_fn = expval_param_shift(
-            expval_sq_tape, argnum, shift, gradient_recipes, f0
+            expval_sq_tape, argnum, shifts, gradient_recipes, f0
         )
         gradient_tapes.extend(pdA2_tapes)
 
@@ -382,7 +416,7 @@ def var_param_shift(tape, argnum, shift=np.pi / 2, gradient_recipes=None, f0=Non
 
 @gradient_transform
 def param_shift(
-    tape, argnum=None, shift=np.pi / 2, gradient_recipes=None, fallback_fn=finite_diff, f0=None
+    tape, argnum=None, shifts=None, gradient_recipes=None, fallback_fn=finite_diff, f0=None
 ):
     r"""Transform a QNode to compute the parameter-shift gradient of all gate
     parameters with respect to its inputs.
@@ -392,10 +426,10 @@ def param_shift(
         argnum (int or list[int] or None): Trainable parameter indices to differentiate
             with respect to. If not provided, the derivative with respect to all
             trainable indices are returned.
-        shift (float): The shift value to use for the two-term parameter-shift formula.
-            Only valid if the operation in question supports the two-term parameter-shift
-            rule (that is, it has two distinct eigenvalues) and ``gradient_recipes``
-            is ``None``.
+        shifts (list[tuple[int or float]]): List containing tuples of shift values.
+            If provided, one tuple of shifts should be given per trainable parameter
+            and the tuple should match the number of frequencies for that parameter.
+            If unspecified, equidistant shifts are assumed.
         gradient_recipes (tuple(list[list[float]] or None)): List of gradient recipes
             for the parameter-shift method. One gradient recipe must be provided
             per trainable parameter.
@@ -539,36 +573,35 @@ def param_shift(
          [ 0.69916862  0.34072424  0.69202359]]
     """
 
-    # perform gradient method validation
     if any(m.return_type is qml.operation.State for m in tape.measurements):
         raise ValueError(
             "Computing the gradient of circuits that return the state is not supported."
         )
 
-    _gradient_analysis(tape)
-    gradient_tapes = []
-
     if argnum is None and not tape.trainable_params:
         warnings.warn(
             "Attempted to compute the gradient of a tape with no trainable parameters. "
             "If this is unintended, please mark trainable parameters in accordance with the "
             "chosen auto differentiation framework, or via the 'tape.trainable_params' property."
         )
-        return gradient_tapes, lambda _: np.zeros([tape.output_dim, len(tape.trainable_params)])
+        return [], lambda _: np.zeros([tape.output_dim, len(tape.trainable_params)])
+
+    _gradient_analysis(tape)
+    method = "analytic" if fallback_fn is None else "best"
 
     # TODO: replace the JacobianTape._grad_method_validation
     # functionality before deprecation.
-    method = "analytic" if fallback_fn is None else "best"
     diff_methods = tape._grad_method_validation(method)
-    all_params_grad_method_zero = all(g == "0" for g in diff_methods)
-    if all_params_grad_method_zero:
-        return gradient_tapes, lambda _: np.zeros([tape.output_dim, len(tape.trainable_params)])
+
+    if all(g == "0" for g in diff_methods):
+        return [], lambda _: np.zeros([tape.output_dim, len(tape.trainable_params)])
 
     # TODO: replace the JacobianTape._choose_params_with_methods
     # functionality before deprecation.
     method_map = dict(tape._choose_params_with_methods(diff_methods, argnum))
 
     # If there are unsupported operations, call the fallback gradient function
+    gradient_tapes = []
     unsupported_params = {idx for idx, g in method_map.items() if g == "F"}
     argnum = [i for i, dm in method_map.items() if dm == "A"]
 
@@ -589,9 +622,9 @@ def param_shift(
         gradient_recipes = [None] * len(argnum)
 
     if any(m.return_type is qml.operation.Variance for m in tape.measurements):
-        g_tapes, fn = var_param_shift(tape, argnum, shift, gradient_recipes, f0)
+        g_tapes, fn = var_param_shift(tape, argnum, shifts, gradient_recipes, f0)
     else:
-        g_tapes, fn = expval_param_shift(tape, argnum, shift, gradient_recipes, f0)
+        g_tapes, fn = expval_param_shift(tape, argnum, shifts, gradient_recipes, f0)
 
     gradient_tapes.extend(g_tapes)
 
diff --git a/pennylane/gradients/parameter_shift_cv.py b/pennylane/gradients/parameter_shift_cv.py
index f06fecf9f..dde10dfd2 100644
--- a/pennylane/gradients/parameter_shift_cv.py
+++ b/pennylane/gradients/parameter_shift_cv.py
@@ -25,7 +25,8 @@ import pennylane as qml
 
 from .gradient_transform import gradient_transform
 from .finite_difference import finite_diff, generate_shifted_tapes
-from .parameter_shift import expval_param_shift, _get_operation_recipe, _process_gradient_recipe
+from .parameter_shift import expval_param_shift, _get_operation_recipe
+from .general_shift_rules import process_shifts
 
 
 def _grad_method(tape, idx):
@@ -164,7 +165,7 @@ def _transform_observable(obs, Z, device_wires):
     return qml.PolyXP(A, wires=device_wires)
 
 
-def var_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, gradient_recipes=None, f0=None):
+def var_param_shift(tape, dev_wires, argnum=None, shifts=None, gradient_recipes=None, f0=None):
     r"""Partial derivative using the first-order or second-order parameter-shift rule of a tape
     consisting of a mixture of expectation values and variances of observables.
 
@@ -182,10 +183,10 @@ def var_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, gradient_reci
         argnum (int or list[int] or None): Trainable parameter indices to differentiate
             with respect to. If not provided, the derivative with respect to all
             trainable indices are returned.
-        shift (float): The shift value to use for the two-term parameter-shift formula.
-            Only valid if the operation in question supports the two-term parameter-shift
-            rule (that is, it has two distinct eigenvalues) and ``gradient_recipes``
-            is ``None``.
+        shifts (list[tuple[int or float]]): List containing tuples of shift values.
+            If provided, one tuple of shifts should be given per trainable parameter
+            and the tuple should match the number of frequencies for that parameter.
+            If unspecified, equidistant shifts are assumed.
         gradient_recipes (tuple(list[list[float]] or None)): List of gradient recipes
             for the parameter-shift method. One gradient recipe must be provided
             per trainable parameter.
@@ -217,7 +218,7 @@ def var_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, gradient_reci
     gradient_tapes = [expval_tape]
 
     # evaluate the analytic derivative of <A>
-    pdA_tapes, pdA_fn = expval_param_shift(expval_tape, argnum, shift, gradient_recipes, f0)
+    pdA_tapes, pdA_fn = expval_param_shift(expval_tape, argnum, shifts, gradient_recipes, f0)
     gradient_tapes.extend(pdA_tapes)
 
     # Store the number of first derivative tapes, so that we know
@@ -248,7 +249,7 @@ def var_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, gradient_reci
     # may be non-zero. Here, we calculate the analytic derivatives of the <A^2>
     # observables.
     pdA2_tapes, pdA2_fn = second_order_param_shift(
-        expval_sq_tape, dev_wires, argnum, shift, gradient_recipes
+        expval_sq_tape, dev_wires, argnum, shifts, gradient_recipes
     )
     gradient_tapes.extend(pdA2_tapes)
 
@@ -266,7 +267,7 @@ def var_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, gradient_reci
     return gradient_tapes, processing_fn
 
 
-def second_order_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, gradient_recipes=None):
+def second_order_param_shift(tape, dev_wires, argnum=None, shifts=None, gradient_recipes=None):
     r"""Generate the second-order CV parameter-shift tapes and postprocessing methods required
     to compute the gradient of a gate parameter with respect to an
     expectation value.
@@ -288,10 +289,10 @@ def second_order_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, grad
         argnum (int or list[int] or None): Trainable parameter indices to differentiate
             with respect to. If not provided, the derivative with respect to all
             trainable indices are returned.
-        shift (float): The shift value to use for the two-term parameter-shift formula.
-            Only valid if the operation in question supports the two-term parameter-shift
-            rule (that is, it has two distinct eigenvalues) and ``gradient_recipes``
-            is ``None``.
+        shifts (list[tuple[int or float]]): List containing tuples of shift values.
+            If provided, one tuple of shifts should be given per trainable parameter
+            and the tuple should match the number of frequencies for that parameter.
+            If unspecified, equidistant shifts are assumed.
         gradient_recipes (tuple(list[list[float]] or None)): List of gradient recipes
             for the parameter-shift method. One gradient recipe must be provided
             per trainable parameter.
@@ -323,12 +324,16 @@ def second_order_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, grad
         shapes.append(1)
 
         # get the gradient recipe for the trainable parameter
-        recipe = gradient_recipes[argnum.index(idx)]
-        recipe = recipe or _get_operation_recipe(tape, idx, shift=shift)
-        recipe = _process_gradient_recipe(recipe)
-        coeffs, multipliers, shifts = recipe
-
-        if len(shifts) != 2:
+        arg_idx = argnum.index(idx)
+        recipe = gradient_recipes[arg_idx]
+        if recipe is not None:
+            recipe = process_shifts(np.array(recipe).T)
+        else:
+            op_shifts = None if shifts is None else shifts[arg_idx]
+            recipe = _get_operation_recipe(tape, idx, shifts=op_shifts)
+        coeffs, multipliers, op_shifts = recipe
+
+        if len(op_shifts) != 2:
             # The 2nd order CV parameter-shift rule only accepts two-term shifts
             raise NotImplementedError(
                 "Taking the analytic gradient for order-2 operators is "
@@ -336,7 +341,7 @@ def second_order_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, grad
                 "gradient recipe of more than two terms."
             )
 
-        shifted_tapes = generate_shifted_tapes(tape, idx, shifts, multipliers)
+        shifted_tapes = generate_shifted_tapes(tape, idx, op_shifts, multipliers)
 
         # evaluate transformed observables at the original parameter point
         # first build the Heisenberg picture transformation matrix Z
@@ -466,7 +471,7 @@ def param_shift_cv(
     tape,
     dev,
     argnum=None,
-    shift=np.pi / 2,
+    shifts=None,
     gradient_recipes=None,
     fallback_fn=finite_diff,
     f0=None,
@@ -481,10 +486,10 @@ def param_shift_cv(
         argnum (int or list[int] or None): Trainable parameter indices to differentiate
             with respect to. If not provided, the derivative with respect to all
             trainable indices are returned.
-        shift (float): The shift value to use for the two-term parameter-shift formula.
-            Only valid if the operation in question supports the two-term parameter-shift
-            rule (that is, it has two distinct eigenvalues) and ``gradient_recipes``
-            is ``None``.
+        shifts (list[tuple[int or float]]): List containing tuples of shift values.
+            If provided, one tuple of shifts should be given per trainable parameter
+            and the tuple should match the number of frequencies for that parameter.
+            If unspecified, equidistant shifts are assumed.
         gradient_recipes (tuple(list[list[float]] or None)): List of gradient recipes
             for the parameter-shift method. One gradient recipe must be provided
             per trainable parameter.
@@ -721,17 +726,17 @@ def param_shift_cv(
     gradient_recipes = gradient_recipes or [None] * len(argnum)
 
     if var_present:
-        _update(var_param_shift(tape, dev.wires, argnum, shift, gradient_recipes, f0))
+        _update(var_param_shift(tape, dev.wires, argnum, shifts, gradient_recipes, f0))
 
     else:
         # Only expectation values were specified
         if first_order_params:
-            _update(expval_param_shift(tape, first_order_params, shift, gradient_recipes, f0))
+            _update(expval_param_shift(tape, first_order_params, shifts, gradient_recipes, f0))
 
         if second_order_params:
             _update(
                 second_order_param_shift(
-                    tape, dev.wires, second_order_params, shift, gradient_recipes
+                    tape, dev.wires, second_order_params, shifts, gradient_recipes
                 )
             )
 
diff --git a/pennylane/operation.py b/pennylane/operation.py
index 29d113937..ccbd66c9b 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -1083,7 +1083,13 @@ class Operation(Operator):
 
         Default is ``'F'``, or ``None`` if the Operation has zero parameters.
         """
-        return None if self.num_params == 0 else "F"
+        if self.num_params == 0:
+            return None
+        try:
+            self.parameter_frequencies  # pylint:disable=pointless-statement
+            return "A"
+        except OperatorPropertyUndefined:
+            return "F"
 
     grad_recipe = None
     r"""tuple(Union(list[list[float]], None)) or None: Gradient recipe for the
@@ -1134,28 +1140,38 @@ class Operation(Operator):
         """
         raise NotImplementedError
 
-    def get_parameter_shift(self, idx, shift=np.pi / 2):
-        """Multiplier and shift for the parameter at position idx, based on its gradient recipe.
+    def get_parameter_shift(self, idx, shift=None):
+        r"""Multiplier and shift for the given parameter, based on its gradient recipe.
 
         Args:
             idx (int): parameter index
+            shift (float or None): The shift value to use for the two-term
+                parameter-shift rule. This is only used if the operation
+                does not have :attr:`Operator.grad_recipe` defined. If ``None``,
+                a shift value of :math:`\pi/2` is used.
 
         Returns:
             list[[float, float, float]]: list of multiplier, coefficient, shift for each term in the gradient recipe
         """
         # get the gradient recipe for this parameter
         recipe = self.grad_recipe[idx]
+        if recipe is not None:
+            return recipe
 
         # Default values
+        if shift is None:
+            shift = np.pi / 2
         multiplier = 0.5 / np.sin(shift)
         a = 1
 
         # We set the following default recipe:
         # ∂f(x) = c*f(a*x+s) - c*f(a*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
-        return param_shift
+        default_param_shift = [
+            [multiplier, a, shift],
+            [-multiplier, a, -shift],  # pylint: disable=invalid-unary-operand-type
+        ]
+        return default_param_shift
 
     @property
     def parameter_frequencies(self):
@@ -1327,8 +1343,6 @@ class Operation(Operator):
                 assert (
                     len(self.grad_recipe) == self.num_params
                 ), "Gradient recipe must have one entry for each parameter!"
-        else:
-            assert self.grad_recipe is None, "Gradient recipe is only used by the A method!"
 
 
 class Channel(Operation, abc.ABC):
@@ -1632,7 +1646,7 @@ class Tensor(Observable):
 
     def __copy__(self):
         cls = self.__class__
-        copied_op = cls.__new__(cls)
+        copied_op = cls.__new__(cls)  # pylint: disable=no-value-for-parameter
         copied_op.obs = self.obs.copy()
         copied_op._eigvals_cache = self._eigvals_cache
         return copied_op


Test patch:
diff --git a/tests/gradients/test_gradient_transform.py b/tests/gradients/test_gradient_transform.py
index 2171ce37e..d86610988 100644
--- a/tests/gradients/test_gradient_transform.py
+++ b/tests/gradients/test_gradient_transform.py
@@ -76,7 +76,8 @@ class TestGradientTransformIntegration:
             qml.CNOT(wires=[0, 1])
             return qml.expval(qml.PauliZ(0)), qml.var(qml.PauliX(1))
 
-        grad_fn = qml.gradients.param_shift(circuit, shift=np.pi / 4)
+        shifts = [(np.pi / 4,), (np.pi / 3,)]
+        grad_fn = qml.gradients.param_shift(circuit, shifts=shifts)
 
         w = np.array([0.543, -0.654], requires_grad=True)
         res = grad_fn(w)
@@ -85,7 +86,7 @@ class TestGradientTransformIntegration:
         expected = np.array([[-np.sin(x), 0], [0, -2 * np.cos(y) * np.sin(y)]])
         assert np.allclose(res, expected, atol=tol, rtol=0)
 
-        assert spy.call_args[0][2] == np.pi / 4
+        assert spy.call_args[0][2] == shifts
 
     def test_expansion(self, mocker, tol):
         """Test that a gradient transform correctly
diff --git a/tests/gradients/test_parameter_shift.py b/tests/gradients/test_parameter_shift.py
index 3633037f7..ef5d14187 100644
--- a/tests/gradients/test_parameter_shift.py
+++ b/tests/gradients/test_parameter_shift.py
@@ -52,7 +52,7 @@ class TestGradAnalysis:
             qml.CNOT(wires=[0, 1])
             qml.probs(wires=[0, 1])
 
-        spy = mocker.spy(tape, "_grad_method")
+        spy = mocker.spy(qml.operation, "has_grad_method")
         _gradient_analysis(tape)
         spy.assert_called()
 
@@ -60,7 +60,7 @@ class TestGradAnalysis:
         assert tape._par_info[1]["grad_method"] == "A"
         assert tape._par_info[2]["grad_method"] == "A"
 
-        spy = mocker.spy(tape, "_grad_method")
+        spy = mocker.spy(qml.operation, "has_grad_method")
         _gradient_analysis(tape)
         spy.assert_not_called()
 
@@ -239,7 +239,7 @@ class TestParamShift:
         assert res.size == 0
         assert np.all(res == np.array([[]]))
 
-    def test_y0(self, mocker):
+    def test_y0(self):
         """Test that if the gradient recipe has a zero-shift component, then
         the tape is executed only once using the current parameter
         values."""
@@ -256,7 +256,7 @@ class TestParamShift:
         # one tape per parameter, plus one global call
         assert len(tapes) == tape.num_params + 1
 
-    def test_y0_provided(self, mocker):
+    def test_y0_provided(self):
         """Test that if the original tape output is provided, then
         the tape is executed only once using the current parameter
         values."""
@@ -311,6 +311,10 @@ class TestParamShift:
         shift rule"""
 
         class RX(qml.RX):
+            @property
+            def parameter_frequencies(self):
+                raise qml.operation.OperatorPropertyUndefined
+
             @property
             def grad_recipe(self):
                 # The gradient is given by [f(2x) - f(0)] / (2 sin(x)), by subsituting
@@ -335,6 +339,39 @@ class TestParamShift:
         grad = fn(dev.batch_execute(tapes))
         assert np.allclose(grad, -np.sin(x))
 
+    @pytest.mark.parametrize("shifts", [None, [(0.1, 0.3, 0.5)]])
+    def test_default_two_term_rule(self, shifts):
+        """Test that the default two-term shift rule is used if no
+        other information is available to obtain the shift rule."""
+
+        class RX(qml.RX):
+            @property
+            def parameter_frequencies(self):
+                raise qml.operation.OperatorPropertyUndefined
+
+            def generator(self):
+                raise GeneratorUndefinedError(f"Operation {self.name} does not have a generator")
+
+        x = np.array(0.654, requires_grad=True)
+        dev = qml.device("default.qubit", wires=2)
+
+        with qml.tape.JacobianTape() as tape:
+            RX(x, wires=0)
+            qml.expval(qml.PauliZ(0))
+
+        tapes, fn = qml.gradients.param_shift(tape, shifts=shifts)
+
+        assert len(tapes) == 2
+        if shifts is None:
+            assert tapes[0].operations[0].data[0] == x + np.pi / 2
+            assert tapes[1].operations[0].data[0] == x - np.pi / 2
+        else:
+            assert tapes[0].operations[0].data[0] == x + shifts[0][0]
+            assert tapes[1].operations[0].data[0] == x - shifts[0][0]
+
+        grad = fn(dev.batch_execute(tapes))
+        assert np.allclose(grad, -np.sin(x))
+
 
 class TestParameterShiftRule:
     """Tests for the parameter shift implementation"""
@@ -354,7 +391,7 @@ class TestParameterShiftRule:
 
         tape.trainable_params = {1}
 
-        tapes, fn = qml.gradients.param_shift(tape, shift=shift)
+        tapes, fn = qml.gradients.param_shift(tape, shifts=[(shift,)])
         assert len(tapes) == 2
 
         autograd_val = fn(dev.batch_execute(tapes))
@@ -364,7 +401,7 @@ class TestParameterShiftRule:
         ) / 2
         assert np.allclose(autograd_val, manualgrad_val, atol=tol, rtol=0)
 
-        assert spy.call_args[1]["shift"] == shift
+        assert spy.call_args[1]["shifts"] == (shift,)
 
         # compare to finite differences
         tapes, fn = qml.gradients.finite_diff(tape)
@@ -386,7 +423,7 @@ class TestParameterShiftRule:
 
         tape.trainable_params = {1, 2, 3}
 
-        tapes, fn = qml.gradients.param_shift(tape, shift=shift)
+        tapes, fn = qml.gradients.param_shift(tape, shifts=[(shift,)] * 3)
         assert len(tapes) == 2 * len(tape.trainable_params)
 
         autograd_val = fn(dev.batch_execute(tapes))
@@ -402,7 +439,7 @@ class TestParameterShiftRule:
             manualgrad_val[0, idx] = (forward - backward) / 2
 
         assert np.allclose(autograd_val, manualgrad_val, atol=tol, rtol=0)
-        assert spy.call_args[1]["shift"] == shift
+        assert spy.call_args[1]["shifts"] == (shift,)
 
         # compare to finite differences
         tapes, fn = qml.gradients.finite_diff(tape)
diff --git a/tests/gradients/test_parameter_shift_cv.py b/tests/gradients/test_parameter_shift_cv.py
index d58054223..6a27b415f 100644
--- a/tests/gradients/test_parameter_shift_cv.py
+++ b/tests/gradients/test_parameter_shift_cv.py
@@ -434,7 +434,11 @@ class TestExpectationQuantumGradients:
     """Tests for the quantum gradients of various gates
     with expectation value output"""
 
-    def test_rotation_gradient(self, mocker, tol):
+    @pytest.mark.parametrize(
+        "gradient_recipes",
+        [None, ([[1 / np.sqrt(2), 1, np.pi / 4], [-1 / np.sqrt(2), 1, -np.pi / 4]],)],
+    )
+    def test_rotation_gradient(self, gradient_recipes, mocker, tol):
         """Test the gradient of the rotation gate"""
         dev = qml.device("default.gaussian", wires=2, hbar=hbar)
 
@@ -450,11 +454,11 @@ class TestExpectationQuantumGradients:
 
         spy2 = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
 
-        tapes, fn = param_shift_cv(tape, dev)
+        tapes, fn = param_shift_cv(tape, dev, gradient_recipes=gradient_recipes)
         grad_A = fn(dev.batch_execute(tapes))
         spy2.assert_not_called()
 
-        tapes, fn = param_shift_cv(tape, dev, force_order2=True)
+        tapes, fn = param_shift_cv(tape, dev, gradient_recipes=gradient_recipes, force_order2=True)
         grad_A2 = fn(dev.batch_execute(tapes))
         spy2.assert_called()
 
diff --git a/tests/interfaces/test_batch_autograd.py b/tests/interfaces/test_batch_autograd.py
index d40625b36..f98eba3ab 100644
--- a/tests/interfaces/test_batch_autograd.py
+++ b/tests/interfaces/test_batch_autograd.py
@@ -66,13 +66,16 @@ class TestAutogradExecuteUnitTests:
                 qml.expval(qml.PauliZ(0))
 
             return execute(
-                [tape], device, gradient_fn=param_shift, gradient_kwargs={"shift": np.pi / 4}
+                [tape],
+                device,
+                gradient_fn=param_shift,
+                gradient_kwargs={"shifts": [(np.pi / 4,)] * 2},
             )[0]
 
         res = qml.jacobian(cost)(a, device=dev)
 
         for args in spy.call_args_list:
-            assert args[1]["shift"] == np.pi / 4
+            assert args[1]["shifts"] == [(np.pi / 4,)] * 2
 
     def test_incorrect_mode(self):
         """Test that an error is raised if a gradient transform
diff --git a/tests/interfaces/test_batch_jax.py b/tests/interfaces/test_batch_jax.py
index c38026817..aef2bbed3 100644
--- a/tests/interfaces/test_batch_jax.py
+++ b/tests/interfaces/test_batch_jax.py
@@ -48,14 +48,14 @@ class TestJaxExecuteUnitTests:
                 [tape],
                 device,
                 gradient_fn=param_shift,
-                gradient_kwargs={"shift": np.pi / 4},
+                gradient_kwargs={"shifts": [(np.pi / 4,)] * 2},
                 interface=interface,
             )[0][0]
 
         res = jax.grad(cost)(a, device=dev)
 
         for args in spy.call_args_list:
-            assert args[1]["shift"] == np.pi / 4
+            assert args[1]["shifts"] == [(np.pi / 4,)] * 2
 
     def test_incorrect_mode(self, interface):
         """Test that an error is raised if an gradient transform
diff --git a/tests/interfaces/test_batch_tensorflow.py b/tests/interfaces/test_batch_tensorflow.py
index 1a2a56f11..837e2063a 100644
--- a/tests/interfaces/test_batch_tensorflow.py
+++ b/tests/interfaces/test_batch_tensorflow.py
@@ -46,14 +46,14 @@ class TestTensorFlowExecuteUnitTests:
                 [tape],
                 dev,
                 gradient_fn=param_shift,
-                gradient_kwargs={"shift": np.pi / 4},
+                gradient_kwargs={"shifts": [(np.pi / 4,)] * 2},
                 interface="tf",
             )[0]
 
         res = t.jacobian(res, a)
 
         for args in spy.call_args_list:
-            assert args[1]["shift"] == np.pi / 4
+            assert args[1]["shifts"] == [(np.pi / 4,)] * 2
 
     def test_incorrect_mode(self):
         """Test that an error is raised if a gradient transform
diff --git a/tests/interfaces/test_batch_torch.py b/tests/interfaces/test_batch_torch.py
index 5a0db39a8..b926d6083 100644
--- a/tests/interfaces/test_batch_torch.py
+++ b/tests/interfaces/test_batch_torch.py
@@ -44,14 +44,14 @@ class TestTorchExecuteUnitTests:
             [tape],
             dev,
             gradient_fn=param_shift,
-            gradient_kwargs={"shift": np.pi / 4},
+            gradient_kwargs={"shifts": [(np.pi / 4,)] * 2},
             interface="torch",
         )[0]
 
         res.backward()
 
         for args in spy.call_args_list:
-            assert args[1]["shift"] == np.pi / 4
+            assert args[1]["shift"] == [(np.pi / 4,)] * 2
 
     def test_incorrect_mode(self):
         """Test that an error is raised if a gradient transform
diff --git a/tests/optimize/test_lie_algebra.py b/tests/optimize/test_lie_algebra.py
index 80e127eed..99105c494 100644
--- a/tests/optimize/test_lie_algebra.py
+++ b/tests/optimize/test_lie_algebra.py
@@ -323,8 +323,7 @@ def test_docstring_example():
 
     for step in range(6):
         circuit, cost = opt.step_and_cost()
-        print(f"Step {step} - cost {cost}")
-    print(circuit())
+    circuit()
 
     assert np.isclose(cost, -2.23, atol=1e-2)
 
@@ -349,7 +348,6 @@ def test_docstring_example_exact():
 
     for step in range(6):
         _, cost = opt.step_and_cost()
-        print(f"Step {step} - cost {cost}")
     assert np.isclose(cost, -2.23, atol=1e-2)
 
 
@@ -372,4 +370,3 @@ def test_example_shots():
 
     for step in range(3):
         _, cost = opt.step_and_cost()
-        print(f"Step {step} - cost {cost}")
diff --git a/tests/test_operation.py b/tests/test_operation.py
index b63313d5c..90d4b4c7a 100644
--- a/tests/test_operation.py
+++ b/tests/test_operation.py
@@ -183,6 +183,64 @@ class TestOperationConstruction:
         op = DummyOp(x, wires=0)
         assert op.grad_recipe == ([[1.0, 1.0, x], [1.0, 0.0, -x]],)
 
+    def test_default_grad_method_with_frequencies(self):
+        """Test that the correct ``grad_method`` is returned by default
+        if ``parameter_frequencies`` are present.
+        """
+
+        class DummyOp(qml.operation.Operation):
+            r"""Dummy custom operation"""
+            num_wires = 1
+
+            @property
+            def parameter_frequencies(self):
+                return [(0.4, 1.2)]
+
+        x = 0.654
+        op = DummyOp(x, wires=0)
+        assert op.grad_method == "A"
+
+    def test_default_grad_method_with_generator(self):
+        """Test that the correct ``grad_method`` is returned by default
+        if a generator is present to determine parameter_frequencies from.
+        """
+
+        class DummyOp(qml.operation.Operation):
+            r"""Dummy custom operation"""
+            num_wires = 1
+
+            def generator(self):
+                return -0.2 * qml.PauliX(wires=self.wires)
+
+        x = 0.654
+        op = DummyOp(x, wires=0)
+        assert op.grad_method == "A"
+
+    def test_default_grad_method_numeric(self):
+        """Test that the correct ``grad_method`` is returned by default
+        if no information is present to deduce an analytic gradient method.
+        """
+
+        class DummyOp(qml.operation.Operation):
+            r"""Dummy custom operation"""
+            num_wires = 1
+
+        x = 0.654
+        op = DummyOp(x, wires=0)
+        assert op.grad_method == "F"
+
+    def test_default_grad_no_param(self):
+        """Test that the correct ``grad_method`` is returned by default
+        if an operation does not have a parameter.
+        """
+
+        class DummyOp(qml.operation.Operation):
+            r"""Dummy custom operation"""
+            num_wires = 1
+
+        op = DummyOp(wires=0)
+        assert op.grad_method is None
+
     def test_frequencies_default_single_param(self):
         """Test that an operation with default parameter frequencies
         and a single parameter works correctly."""
diff --git a/tests/transforms/test_draw.py b/tests/transforms/test_draw.py
index 3dbfd9765..c378bd2bb 100644
--- a/tests/transforms/test_draw.py
+++ b/tests/transforms/test_draw.py
@@ -219,7 +219,10 @@ class TestLayering:
 
 @pytest.mark.parametrize(
     "transform",
-    [qml.gradients.param_shift(shift=0.2), partial(qml.gradients.param_shift, shift=0.2)],
+    [
+        qml.gradients.param_shift(shifts=[(0.2,)]),
+        partial(qml.gradients.param_shift, shifts=[(0.2,)]),
+    ],
 )
 def test_draw_batch_transform(transform):
     """Test that drawing a batch transform works correctly."""
diff --git a/tests/transforms/test_draw_old.py b/tests/transforms/test_draw_old.py
index 3acba1d45..c197081a1 100644
--- a/tests/transforms/test_draw_old.py
+++ b/tests/transforms/test_draw_old.py
@@ -15,7 +15,7 @@
 """
 Unit tests for the old draw transform.
 """
-import functools
+from functools import partial
 import pytest
 
 import pennylane as qml
@@ -237,7 +237,10 @@ def test_invalid_wires():
 
 @pytest.mark.parametrize(
     "transform",
-    [qml.gradients.param_shift(shift=0.2), functools.partial(qml.gradients.param_shift, shift=0.2)],
+    [
+        qml.gradients.param_shift(shifts=[(0.2,)]),
+        partial(qml.gradients.param_shift, shifts=[(0.2,)]),
+    ],
 )
 def test_draw_batch_transform(transform):
     """Test that drawing a batch transform works correctly"""