Repository: PennyLaneAI/pennylane
Pull request #1932 (2021-11-24 09:23:43)

Adds utility functions for the general parameter-shift rule
**Context:** In #1788, a utility function to generate parameter-shift rules, given an operators spectrum, was added. As per the recent gradient ADR, the decision was made to use this pipeline as the default method for the parameter-shift rule in PennyLane, and deprecate the `Operator.grad_recipe` attribute.

This PR introduces several changes, and new utility function, in preparation for this move.

**Description of the Change:**

- The function `get_shift_rule` was renamed to `_get_shift_rule`, and wrapped by a new user-facing wrapper function `generate_shift_rule`, to emphasize that the shift rule is computed, rather than retrieved.

- As the `grad_recipe` is due to be removed from the parameter-shift logic (kept only as a fallback option), `generate_shift_rule` has been modified to simply return an array of shifts and coefficients.

- `generate_shift_rule` automatically removes negative frequencies, so that it can be used with `qml.fourier.spectrum` (which returns both positive and negative frequencies).

- `_process_gradient_recipe` has been added, which processes the shifts and coefficients to remove redundant shifts. In particular:

  * All small coefficients and shifts (up to a certain tolerance) are set to 0

  * Terms where the coefficients are 0 are removed

  * The terms are sorted as per `abs(shift)`, to ensure that any unshifted terms are first.

  * Terms with the same shift value are combined, with corresponding coefficients summed.

- A utility function, `eigvals_to_frequencies` is provided, allowing users to easily convert generator eigenvalues to the corresponding frequency spectrum.

- A utility function, `frequencies_to_period`, is provided, allowing users to determine the period  of `f(t) = <0|U(t)† O U(t)|0>` given the frequencies of `G` (where `U(t) = exp(iGt)`).

- `generate_shift_rule` can now be used to return iterated shift rules that return higher order derivatives of an operator (e.g., `∂^2 f/∂t^2`

- `generate_multi_shift_rule` can now be used to return iterated shift rules that return higher order derivatives of _multiple_ operators (e.g., `∂^(n+m) f / ∂x^n ∂y^m`). This can be used to generate shift rules for the off-diagonal elements of Hessians.

- `CommutingEvolution` has been updated to take into account these changes.

**Benefits:**

- It is now possible to query for shift rules of multiple operations, to arbitrary order.

- It is easier to understand the output of `generate_shift_rules`.

**Possible Drawbacks:** n/a

**Related GitHub Issues:** n/a


codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1932?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1932](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1932?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (5a094a7) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/d705f0e3cc51bd876c8b910cd25b85bfc08694b0?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (d705f0e) will **decrease** coverage by `0.14%`.
> The diff coverage is `51.92%`.

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

```diff
@@            Coverage Diff             @@
##           master    #1932      +/-   ##
==========================================
- Coverage   98.84%   98.69%   -0.15%     
==========================================
  Files         224      224              
  Lines       17101    17147      +46     
==========================================
+ Hits        16903    16924      +21     
- Misses        198      223      +25     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1932?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/1932/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyYWRpZW50cy9nZW5lcmFsX3NoaWZ0X3J1bGVzLnB5) | `67.10% <46.80%> (-32.90%)` | :arrow_down: |
| [pennylane/gradients/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1932/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyYWRpZW50cy9fX2luaXRfXy5weQ==) | `100.00% <100.00%> (ø)` | |
| [...ylane/templates/subroutines/commuting\_evolution.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1932/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RlbXBsYXRlcy9zdWJyb3V0aW5lcy9jb21tdXRpbmdfZXZvbHV0aW9uLnB5) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1932?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/1932?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [d705f0e...5a094a7](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1932?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).

josh146: [sc-11839]
josh146: Thanks @mariaschuld @dwierichs! Your comments were super helpful. I've implemented every one. Two of the biggest changes are:

- I've removed the term 'gradient recipe' everywhere, as I realize this is technical jargon that makes it harder for external developers. Instead, I've replaced it with 'gradient rule', which should be more familiar due to 'parameter-shift rule'.

- I've updated `frequencies_to_period` to use `np.gcd`, and added tests for non-equidistant frequencies 🙂 

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index d00d0e867..91c2bd766 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -265,38 +265,55 @@
   operation on our quantum circuits for both qubit and CV devices.
   [(#1829)](https://github.com/PennyLaneAI/pennylane/pull/1829)
 
-* For Hamiltonians whose eigenvalue frequency spectrum is known, `qml.gradients.get_shift_rule` is
-  a function that computes the generalized parameter shift rules for the time evolution.
+* Given an operator of the form :math:`U=e^{iHt}`, where :math:`H` has
+  commuting terms and known eigenvalues,
+  `qml.gradients.generate_shift_rule` computes the generalized parameter shift rules for determining
+  the gradient of the expectation value :math:`f(t) = \langle 0|U(t)^\dagger \hat{O} U(t)|0\rangle` on
+  hardware.
   [(#1788)](https://github.com/PennyLaneAI/pennylane/pull/1788)
+  [(#1932)](https://github.com/PennyLaneAI/pennylane/pull/1932)
 
-  Given a Hamiltonian's frequency spectrum of `R` unique frequencies, `qml.gradients.get_shift_rule`
-  returns the parameter shift rules to compute expectation value gradients of the Hamiltonian's
-  time parameter using `2R` shifted cost function evaluations. This becomes cheaper than
-  the standard application of the chain rule and two-term shift rule when `R` is less than the
-  number of Pauli words in the Hamiltonian generator.
+  Given
 
-  For example, a four-term shift rule is generated for the frequency spectrum `[1, 2]`, which
-  corresponds to a generator eigenspectrum of e.g., `[-1, 0, 1]`:
+  .. math:: H = \sum_i a_i h_i,
+
+  where the eigenvalues of :math:`H` are known and all :math:`h_i` commute, we can compute
+  the *frequencies* (the unique positive differences of any two eigenvalues) using
+  `qml.gradients.eigvals_to_frequencies`.
+
+  `qml.gradients.generate_shift_rule` can then be used to compute the parameter
+  shift rules to compute :math:`f'(t)` using `2R` shifted cost function evaluations.
+  This becomes cheaper than the standard application of the chain rule and
+  two-term shift rule when `R` is less than the
+  number of Pauli words in the generator.
+
+  For example, consider the case where :math:`H` has eigenspectrum ``(-1, 0, 1)``:
 
   ```pycon
-  >>> frequencies = (1,2)
-  >>> grad_recipe = qml.gradients.get_shift_rule(frequencies)
-  >>> grad_recipe
-  ([[0.8535533905932737, 1, 0.7853981633974483], [-0.14644660940672624, 1, 2.356194490192345],
-    [-0.8535533905932737, 1, -0.7853981633974483], [0.14644660940672624, 1, -2.356194490192345]],)
+  >>> frequencies = qml.gradients.eigvals_to_frequencies((-1, 0, 1))
+  >>> frequencies
+  (1, 2)
+  >>> coeffs, shifts = qml.gradients.generate_shift_rule(frequencies)
+  >>> coeffs
+  array([ 0.85355339, -0.85355339, -0.14644661,  0.14644661])
+  >>> shifts
+  array([ 0.78539816, -0.78539816,  2.35619449, -2.35619449])
   ```
 
-  As we can see, `get_shift_rule` returns a tuple containing a list of four nested lists for the
-  four term parameter shift rule. Each term :math:`[c_i, a_i, s_i]` specifies a term in the
-  gradient reconstructed via parameter shifts as
+  As we can see, `generate_shift_rule` returns four coefficients :math:`c_i` and shifts
+  :math:`s_i` corresponding to a four term parameter shift rule. The gradient can then
+  be reconstructed via:
+
+  .. math:: \frac{\partial}{\partial\phi}f = \sum_{i} c_i f(\phi + s_i),
 
-  .. math:: \frac{\partial}{\partial\phi_k}f = \sum_{i} c_i f(a_i \phi_k + s_i).
+  where :math:`f(\phi) = \langle 0|U(\phi)^\dagger \hat{O} U(\phi)|0\rangle`
+  for some observable :math:`\hat{O}` and the unitary :math:`U(\phi)=e^{iH\phi}`.
 
 * A circuit template for time evolution under a commuting Hamiltonian utilizing generalized
   parameter shift rules for cost function gradients is available as `qml.CommutingEvolution`.
   [(#1788)](https://github.com/PennyLaneAI/pennylane/pull/1788)
 
-  If the template is handed a frequency spectrum during its instantiation, then `get_shift_rule`
+  If the template is handed a frequency spectrum during its instantiation, then `generate_shift_rule`
   is internally called to obtain the general parameter shift rules with respect to
   `CommutingEvolution`'s :math:`t` parameter, otherwise the shift rule for a decomposition of
   `CommutingEvolution` will be used.
diff --git a/pennylane/gradients/__init__.py b/pennylane/gradients/__init__.py
index 6dc141209..79142392c 100644
--- a/pennylane/gradients/__init__.py
+++ b/pennylane/gradients/__init__.py
@@ -62,7 +62,9 @@ Utility functions
 
     finite_diff_coeffs
     generate_shifted_tapes
-    get_shift_rule
+    generate_shift_rule
+    generate_multi_shift_rule
+    eigvals_to_frequencies
     compute_vjp
     batch_vjp
     vjp
@@ -288,4 +290,8 @@ from .parameter_shift import param_shift
 from .parameter_shift_cv import param_shift_cv
 from .vjp import compute_vjp, batch_vjp, vjp
 from .hamiltonian_grad import hamiltonian_grad
-from .general_shift_rules import get_shift_rule
+from .general_shift_rules import (
+    eigvals_to_frequencies,
+    generate_shift_rule,
+    generate_multi_shift_rule,
+)
diff --git a/pennylane/gradients/general_shift_rules.py b/pennylane/gradients/general_shift_rules.py
index 1531e9ff2..ab554e7c9 100644
--- a/pennylane/gradients/general_shift_rules.py
+++ b/pennylane/gradients/general_shift_rules.py
@@ -13,78 +13,126 @@
 # limitations under the License.
 """Contains a function for generating generalized parameter shift rules."""
 import functools
+import itertools
 import warnings
+
 import numpy as np
 import pennylane as qml
 
 
-@functools.lru_cache(maxsize=None)
-def get_shift_rule(frequencies, shifts=None):
-    r"""Computes the parameter shift rule for a unitary based on its generator's eigenvalue frequency
-    spectrum.
+def _process_shifts(rule, tol=1e-10):
+    """Utility function to process gradient rules.
 
-    To compute gradients of circuit parameters in variational quantum algorithms, expressions for
-    cost function first derivatives with respect to the variational parameters can be cast into
-    linear combinations of expectation values at shifted parameter values. These "gradient recipes"
-    can be obtained from the unitary generator's eigenvalue frequency spectrum. Details can be
-    found in `arXiv:2107.12390 <https://arxiv.org/abs/2107.12390>`__.
+    Args:
+        rule (array): a ``(N, M)`` array corresponding to ``M`` terms
+            with ``N-1`` simultaneous function shifts. The first row of
+            the array corresponds to the linear combination coefficients;
+            subsequent rows correspond to parameter shifts for these coefficients.
+        tol (float): floating point tolerance used when comparing shifts/coefficients
+
+    This utility function accepts coefficients and shift values, and performs the following
+    processing:
+
+    - Removes all small (within absolute tolerance ``tol``) coefficients and shifts
+
+    - Removes terms with the coefficients are 0
+
+    - Terms with the same shift value are combined into a single term.
+
+    - Finally, the terms are sorted according to the absolute value of ``shift``,
+      ensuring that, if there is a zero-shift term, this is returned first.
+    """
+    # remove all small coefficients and shifts
+    rule[np.abs(rule) < tol] = 0
+
+    # 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 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
+
+
+@functools.lru_cache(maxsize=None)
+def eigvals_to_frequencies(eigvals):
+    r"""Convert an eigenvalue spectrum to frequency values, defined
+    as the the set of positive, unique differences of the eigenvalues in the spectrum.
 
     Args:
-        frequencies (tuple[int or float]): the tuple of eigenvalue frequencies. Eigenvalue
-            frequencies are defined as the unique positive differences obtained from a set of
-            eigenvalues.
-        shifts (tuple[int or float]): the tuple of shift values. If unspecified, the equidistant
-            shifts `(2 * mu - 1) * pi / (2 * n_freqs * freq_min)` are assumed, where `n_freqs` is
-            the number of frequencies, `freq_min` is the smallest frequency, and
-            `mu = 1, ..., n_freqs`. If supplied, the length of this tuple should match the number
-            of given frequencies.
+        eigvals (tuple[int, float]): eigenvalue spectra
 
     Returns:
-         tuple: a tuple of one nested list describing the gradient recipe
-         for the parameter-shift method.
-         This is a tuple with one nested list per operation parameter. For
-         parameter :math:`\phi_k`, the nested list contains elements of the form
-         :math:`[c_i, a_i, s_i]` where :math:`i` is the index of the
-         term, resulting in a gradient recipe of
+        tuple[int, float]: frequencies
 
-            .. math:: \frac{\partial}{\partial\phi_k}f = \sum_{i} c_i f(a_i \phi_k + s_i).
+    **Example**
 
+    >>> eigvals = (-0.5, 0, 0, 0.5)
+    >>> eigvals_to_frequencies(eigvals)
+    (0.5, 1.0)
+    """
+    unique_eigvals = sorted(set(eigvals))
+    return tuple({j - i for i, j in itertools.combinations(unique_eigvals, 2)})
 
-    Raises:
-        ValueError: if ``frequencies`` is not a list of unique positive values, or if ``shifts``
-            (if specified) is not a list of unique values the same length as ``frequencies``.
 
-    **Examples**
+@functools.lru_cache(maxsize=None)
+def frequencies_to_period(frequencies, decimals=5):
+    r"""Returns the period of a Fourier series as defined
+    by a set of frequencies.
 
-    An example of obtaining the frequencies from a set of unique eigenvalues and obtaining the
-    parameter shift rule:
+    The period is simply :math:`2\pi/gcd(frequencies)`,
+    where :math:`\text{gcd}` is the greatest common divisor.
 
-    >>> unique_eigenvals = sorted([1, -1, 0])
-    >>> from itertools import combinations
-    >>> frequencies = {abs(i - j) for i, j in combinations(unique_eigenvals, 2)}
-    >>> get_shift_rule(tuple(frequencies))
-    ([[0.8535533905932737, 1, 0.7853981633974483], [-0.14644660940672624, 1, 2.356194490192345],
-    [-0.8535533905932737, 1, -0.7853981633974483], [0.14644660940672624, 1, -2.356194490192345]],)
+    Args:
+        spectra (tuple[int, float]): frequency spectra
+        decimals (int): Number of decimal places to round to
+            if there are non-integral frequencies.
 
-    An example with explicitly specified shift values:
+    Returns:
+        tuple[int, float]: frequencies
 
-    >>> frequencies = (1, 2, 4)
-    >>> shifts = (np.pi / 3, 2 * np.pi / 3, np.pi / 4)
-    >>> get_shift_rule(frequencies, shifts)
-    ([[-2.0907702751760278, 1, 1.0471975511965976], [0.2186308015824754, 1, 2.0943951023931953],
-    [3.0000000000000004, 1, 0.7853981633974483], [2.0907702751760278, 1, -1.0471975511965976],
-    [-0.2186308015824754, 1, -2.0943951023931953], [-3.0000000000000004, 1, -0.7853981633974483]],)
+    **Example**
+
+    >>> frequencies = (0.5, 1.0)
+    >>> frequencies_to_period(frequencies)
+    12.566370614359172
     """
+    try:
+        gcd = np.gcd.reduce(frequencies)
+
+    except TypeError:
+        # np.gcd only support integer frequencies
+        exponent = 10 ** decimals
+        frequencies = np.round(frequencies, decimals) * exponent
+        gcd = np.gcd.reduce(np.int64(frequencies)) / exponent
+
+    return 2 * np.pi / gcd
 
+
+@functools.lru_cache(maxsize=None)
+def _get_shift_rule(frequencies, shifts=None):
     n_freqs = len(frequencies)
     frequencies = qml.math.sort(qml.math.stack(frequencies))
     freq_min = frequencies[0]
 
     if len(set(frequencies)) != n_freqs or freq_min <= 0:
         raise ValueError(
-            "Expected frequencies to be a list of unique positive values, instead got {}.".format(
-                frequencies
-            )
+            f"Expected frequencies to be a list of unique positive values, instead got {frequencies}."
         )
 
     mu = np.arange(1, n_freqs + 1)
@@ -96,12 +144,10 @@ def get_shift_rule(frequencies, shifts=None):
         shifts = qml.math.sort(qml.math.stack(shifts))
         if len(shifts) != n_freqs:
             raise ValueError(
-                "Expected number of shifts to equal the number of frequencies ({}), instead got {}.".format(
-                    n_freqs, shifts
-                )
+                f"Expected number of shifts to equal the number of frequencies ({n_freqs}), instead got {shifts}."
             )
         if len(set(shifts)) != n_freqs:
-            raise ValueError("Shift values must be unique, instead got {}".format(shifts))
+            raise ValueError(f"Shift values must be unique, instead got {shifts}")
 
         equ_shifts = all(np.isclose(shifts, (2 * mu - 1) * np.pi / (2 * n_freqs * freq_min)))
 
@@ -117,13 +163,163 @@ def get_shift_rule(frequencies, shifts=None):
         det_sin_matr = np.linalg.det(sin_matr)
         if abs(det_sin_matr) < 1e-6:
             warnings.warn(
-                "Solving linear problem with near zero determinant ({}) may give unstable results for the parameter shift rules.".format(
-                    det_sin_matr
-                )
+                f"Solving linear problem with near zero determinant ({det_sin_matr}) "
+                "may give unstable results for the parameter shift rules."
             )
 
         coeffs = -2 * np.linalg.solve(sin_matr.T, frequencies)
 
     coeffs = np.concatenate((coeffs, -coeffs))
     shifts = np.concatenate((shifts, -shifts))
-    return ([[coeffs[mu], 1, shifts[mu]] for mu in range(0, 2 * n_freqs)],)
+    return np.stack([coeffs, shifts])
+
+
+@functools.lru_cache()
+def generate_shift_rule(frequencies, shifts=None, order=1):
+    r"""Computes the parameter shift rule for a unitary based on its generator's eigenvalue
+    frequency spectrum.
+
+    To compute gradients of circuit parameters in variational quantum algorithms, expressions for
+    cost function first derivatives with respect to the variational parameters can be cast into
+    linear combinations of expectation values at shifted parameter values. The coefficients and
+    shifts defining the linear combination can be obtained from the unitary generator's eigenvalue
+    frequency spectrum. Details can be found in https://arxiv.org/abs/2107.12390.
+
+    Args:
+        frequencies (tuple[int or float]): The tuple of eigenvalue frequencies. Eigenvalue
+            frequencies are defined as the unique positive differences obtained from a set of
+            eigenvalues.
+        shifts (tuple[int or float]): the tuple of shift values. If unspecified,
+            equidistant shifts are assumed. If supplied, the length of this tuple should match the
+            number of given frequencies.
+        order (int): the order of differentiation to compute the shift rule for
+
+    Returns:
+        tuple: a tuple of coefficients and shifts describing the gradient rule for the
+        parameter-shift method. For parameter :math:`\phi`, the coefficients :math:`c_i` and the
+        shifts :math:`s_i` combine to give a gradient rule of the following form:
+
+        .. math:: \frac{\partial}{\partial\phi}f = \sum_{i} c_i f(\phi + s_i).
+
+        where :math:`f(\phi) = \langle 0|U(\phi)^\dagger \hat{O} U(\phi)|0\rangle`
+        for some observable :math:`\hat{O}` and the unitary :math:`U(\phi)=e^{iH\phi}`.
+
+    Raises:
+        ValueError: if ``frequencies`` is not a list of unique positive values, or if ``shifts``
+            (if specified) is not a list of unique values the same length as ``frequencies``.
+
+    **Examples**
+
+    An example of obtaining the frequencies from generator eigenvalues, and obtaining the parameter
+    shift rule:
+
+    >>> eigvals = (-0.5, 0, 0, 0.5)
+    >>> frequencies = eigvals_to_frequencies(eigvals)
+    >>> generate_shift_rule(frequencies)
+    [[ 0.85355339, -0.85355339, -0.14644661,  0.14644661],
+     [0.78539816,  -0.78539816,  2.35619449, -2.35619449]]
+
+    An example with explicitly specified shift values:
+
+    >>> frequencies = (1, 2, 4)
+    >>> shifts = (np.pi / 3, 2 * np.pi / 3, np.pi / 4)
+    >>> generate_shift_rule(frequencies, shifts)
+    [[ 3.        , -3.        , -2.09077028,  2.09077028, 0.2186308, -0.2186308 ],
+     [ 0.78539816, -0.78539816,  1.04719755, -1.04719755, 2.0943951, -2.0943951 ]]
+
+    Higher order shift rules (corresponding to the :math:`n`-th derivative of the parameter) can be
+    requested via the ``order`` argument. For example, to extract the second order shift rule for a
+    gate with generator :math:`X/2`:
+
+    >>> eigvals = (0.5, -0.5) frequencies = eigvals_to_frequencies(eigvals)
+    >>> generate_shift_rule(frequencies, order=2)
+    [[-0.5       ,   0.5       ],
+     [ 0.        ,  -3.14159265]]
+
+    This corresponds to :math:`\frac{\partial^2 f}{\partial phi^2} = \frac{1}{2} \left[f(\phi) -
+    f(\phi-\pi)\right]`.
+    """
+    frequencies = tuple(f for f in frequencies if f > 0)
+    rule = _get_shift_rule(frequencies, shifts=shifts)
+
+    if order > 1:
+        combined_rules = []
+        T = frequencies_to_period(frequencies)
+
+        for partial_rule in itertools.product(rule.T, repeat=order):
+            c, s = np.stack(partial_rule).T
+            new_shift = np.mod(sum(s) + 0.5 * T, T) - 0.5 * T
+            combined_rules.append(np.stack([np.prod(c), new_shift]))
+
+        # combine all terms in the linear combination into a single
+        # 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)
+
+
+def generate_multi_shift_rule(frequencies, shifts=None, orders=None):
+    r"""Computes the parameter shift rule with respect to two parametrized unitaries,
+    given their generator's eigenvalue frequency spectrum. This corresponds to a
+    shift rule that computes off-diagonal elements of higher order derivative tensors.
+    For the second order, this corresponds to the Hessian.
+
+    Args:
+        frequencies (list[tuple[int or float]]): List of eigenvalue frequencies corresponding
+            to the each parametrized unitary.
+        shifts (list[tuple[int or float]]): List of shift values corresponding to each parametrized
+            unitary. If unspecified, equidistant shifts are assumed. If supplied, the length
+            of each tuple in the list must be the same as the length of the corresponding tuple in
+            ``frequencies``.
+        orders (list[int]): the order of differentiation for each parametrized unitary.
+            If unspecified, the first order derivative shift rule is computed for each parametrized
+            unitary.
+
+    Returns:
+        tuple: a tuple of coefficients, shifts for the first parameter, and shifts for the
+        second parameter, describing the gradient rule
+        for the parameter-shift method.
+
+        For parameters :math:`\phi_a` and :math:`\phi_b`, the
+        coefficients :math:`c_i` and the shifts :math:`s^{(a)}_i`, :math:`s^{(b)}_i`,
+        combine to give a gradient rule of the following form:
+
+        .. math::
+
+            \frac{\partial^2}{\partial\phi_a \partial\phi_b}f
+            = \sum_{i} c_i f(\phi_a + s^{(a)}_i, \phi_b + s^{(b)}_i).
+
+        where :math:`f(\phi_a, \phi_b) = \langle 0|U(\phi_a)^\dagger V(\phi_b)^\dagger \hat{O} V(\phi_b) U(\phi_a)|0\rangle`
+        for some observable :math:`\hat{O}` and unitaries :math:`U(\phi_a)=e^{iH_a\phi_a}` and :math:`V(\phi_b)=e^{iH_b\phi_b}`.
+
+    **Example**
+
+    >>> generate_multi_shift_rule([(1,), (1,)])
+    [[ 0.25      , -0.25      , -0.25      ,  0.25      ],
+     [ 1.57079633,  1.57079633, -1.57079633, -1.57079633],
+     [ 1.57079633, -1.57079633,  1.57079633, -1.57079633]])
+
+    This corresponds to the gradient rule
+
+    .. math::
+
+        \frac{\partial^2 f}{\partial x\partial y} &= \frac{1}{4}
+        \left[f(x+\pi/2, y+\pi/2) - f(x+\pi/2, y-\pi/2)\\
+        &~~~- f(x-\pi/2, y+\pi/2) + f(x-\pi/2, y-\pi/2) \right].
+    """
+    rules = []
+    shifts = shifts or [None] * len(frequencies)
+    orders = orders or [1] * len(frequencies)
+
+    for f, s, o in zip(frequencies, shifts, orders):
+        rule = generate_shift_rule(f, shifts=s, order=o)
+        rules.append(_process_shifts(rule).T)
+
+    combined_rules = []
+
+    for partial_rules in itertools.product(*rules):
+        c, s = np.stack(partial_rules).T
+        combined = np.concatenate([[np.prod(c)], s])
+        combined_rules.append(np.stack(combined))
+
+    return np.stack(combined_rules).T
diff --git a/pennylane/math/single_dispatch.py b/pennylane/math/single_dispatch.py
index 1a55c9648..f2b86e49e 100644
--- a/pennylane/math/single_dispatch.py
+++ b/pennylane/math/single_dispatch.py
@@ -187,6 +187,16 @@ ar.register_function(
 )
 
 
+def _round_tf(tensor, decimals=0):
+    """Implement a TensorFlow version of np.round"""
+    tf = _i("tf")
+    tol = 10 ** decimals
+    return tf.round(tensor * tol) / tol
+
+
+ar.register_function("tensorflow", "round", _round_tf)
+
+
 def _take_tf(tensor, indices, axis=None):
     """Implement a TensorFlow version of np.take"""
     tf = _i("tf")
@@ -323,6 +333,16 @@ ar.autoray._SUBMODULE_ALIASES["torch", "arctan2"] = "torch"
 ar.autoray._FUNC_ALIASES["torch", "arctan2"] = "atan2"
 
 
+def _round_torch(tensor, decimals=0):
+    """Implement a Torch version of np.round"""
+    torch = _i("torch")
+    tol = 10 ** decimals
+    return torch.round(tensor * tol) / tol
+
+
+ar.register_function("torch", "round", _round_torch)
+
+
 def _take_torch(tensor, indices, axis=None):
     """Torch implementation of np.take"""
     torch = _i("torch")
diff --git a/pennylane/operation.py b/pennylane/operation.py
index 0dc624e8f..483091137 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -236,10 +236,10 @@ def classproperty(func):
 
 def _process_data(op):
     if op.name in ("RX", "RY", "RZ", "PhaseShift", "Rot"):
-        return str([d % (2 * np.pi) for d in op.data])
+        return str([qml.math.round(d % (2 * np.pi), 10) for d in op.data])
 
     if op.name in ("CRX", "CRY", "CRZ", "CRot"):
-        return str([d % (4 * np.pi) for d in op.data])
+        return str([qml.math.round(d % (4 * np.pi), 10) for d in op.data])
 
     return str(op.data)
 
diff --git a/pennylane/templates/subroutines/commuting_evolution.py b/pennylane/templates/subroutines/commuting_evolution.py
index f0d46dc89..094af7c9b 100644
--- a/pennylane/templates/subroutines/commuting_evolution.py
+++ b/pennylane/templates/subroutines/commuting_evolution.py
@@ -108,11 +108,9 @@ class CommutingEvolution(Operation):
     grad_method = None
 
     def __init__(self, hamiltonian, time, frequencies=None, shifts=None, do_queue=True, id=None):
-        from pennylane.gradients.general_shift_rules import get_shift_rule
-
         # pylint: disable=import-outside-toplevel
         from pennylane.gradients.general_shift_rules import (
-            get_shift_rule,
+            generate_shift_rule,
         )
 
         if not isinstance(hamiltonian, qml.Hamiltonian):
@@ -121,9 +119,9 @@ class CommutingEvolution(Operation):
 
         trainable_hamiltonian = qml.math.requires_grad(hamiltonian.coeffs)
         if frequencies is not None and not trainable_hamiltonian:
-            self.grad_recipe = (get_shift_rule(frequencies, shifts)[0],) + (None,) * len(
-                hamiltonian.data
-            )
+            c, s = generate_shift_rule(frequencies, shifts)
+            recipe = qml.math.stack([c, qml.math.ones_like(c), s]).T
+            self.grad_recipe = (recipe,) + (None,) * len(hamiltonian.data)
             self.grad_method = "A"
 
         self.hamiltonian = hamiltonian


Test patch:
diff --git a/tests/gradients/test_general_shift_rules.py b/tests/gradients/test_general_shift_rules.py
index d042c169b..06eeed002 100644
--- a/tests/gradients/test_general_shift_rules.py
+++ b/tests/gradients/test_general_shift_rules.py
@@ -17,11 +17,15 @@ import pytest
 import numpy as np
 import pennylane as qml
 from pennylane import numpy as np
-from pennylane.gradients.general_shift_rules import get_shift_rule
+from pennylane.gradients.general_shift_rules import (
+    generate_shift_rule,
+    _get_shift_rule,
+    generate_multi_shift_rule,
+)
 
 
-class TestGetShiftRule:
-    """Tests of input validation and output correctness of function `get_shift_rule`."""
+class TestGenerateShiftRule:
+    """Tests of input validation and output correctness of function `generate_shift_rule`."""
 
     def test_invalid_frequency_spectrum(self):
         """Tests ValueError is raised if input frequency spectrum is non positive or non unique."""
@@ -29,8 +33,8 @@ class TestGetShiftRule:
         non_positive_frequency_spectrum = (-1, 1)
         non_unique_frequency_spectrum = (1, 2, 2, 3)
 
-        assert pytest.raises(ValueError, get_shift_rule, non_positive_frequency_spectrum)
-        assert pytest.raises(ValueError, get_shift_rule, non_unique_frequency_spectrum)
+        assert pytest.raises(ValueError, _get_shift_rule, non_positive_frequency_spectrum)
+        assert pytest.raises(ValueError, _get_shift_rule, non_unique_frequency_spectrum)
 
     def test_invalid_shifts(self):
         """Tests ValueError is raised if specified shifts is not of the same length as
@@ -40,8 +44,8 @@ class TestGetShiftRule:
         invalid_shifts_num = (np.pi / 8, 3 * np.pi / 8, 5 * np.pi / 8)
         non_unique_shifts = (np.pi / 8, 3 * np.pi / 8, 5 * np.pi / 8, np.pi / 8)
 
-        assert pytest.raises(ValueError, get_shift_rule, frequencies, invalid_shifts_num)
-        assert pytest.raises(ValueError, get_shift_rule, frequencies, non_unique_shifts)
+        assert pytest.raises(ValueError, generate_shift_rule, frequencies, invalid_shifts_num)
+        assert pytest.raises(ValueError, generate_shift_rule, frequencies, non_unique_shifts)
 
     def test_two_term_rule_default_shifts(self):
         """Tests the correct two term equidistant rule is generated using default shift pi/2.
@@ -50,11 +54,10 @@ class TestGetShiftRule:
         frequencies = (1,)
 
         n_terms = 2
-        correct_terms = [[0.5, 1.0, np.pi / 2], [-0.5, 1.0, -np.pi / 2]]
+        correct_terms = [[0.5, -0.5], [np.pi / 2, -np.pi / 2]]
+        generated_terms = generate_shift_rule(frequencies)
 
-        generated_terms = get_shift_rule(frequencies)[0]
-
-        assert all([all(np.isclose(generated_terms[i], correct_terms[i])) for i in range(n_terms)])
+        assert np.allclose(generated_terms, correct_terms)
 
     def test_four_term_rule_default_shifts(self):
         """Tests the correct two term equidistant rule is generated using the default shifts [pi/4, 3*pi/4].
@@ -64,15 +67,12 @@ class TestGetShiftRule:
 
         n_terms = 4
         correct_terms = [
-            [0.8535533905932737, 1.0, np.pi / 4],
-            [-0.14644660940672624, 1.0, 3 * np.pi / 4],
-            [-0.8535533905932737, 1.0, -np.pi / 4],
-            [0.14644660940672624, 1.0, -3 * np.pi / 4],
+            [0.8535533905932737, -0.8535533905932737, -0.14644660940672624, 0.14644660940672624],
+            [np.pi / 4, -np.pi / 4, 3 * np.pi / 4, -3 * np.pi / 4],
         ]
+        generated_terms = generate_shift_rule(frequencies)
 
-        generated_terms = get_shift_rule(frequencies)[0]
-
-        assert all([all(np.isclose(generated_terms[i], correct_terms[i])) for i in range(n_terms)])
+        assert np.allclose(generated_terms, correct_terms)
 
     def test_eight_term_rule_non_equidistant_default_shifts(self):
         """Tests the correct non-equidistant eight term shift rule is generated given the
@@ -82,20 +82,21 @@ class TestGetShiftRule:
         frequencies = (1, 4, 5, 6)
 
         n_terms = 8
+
         correct_terms = [
-            [2.8111804455102014, 1.0, np.pi / 8],
-            [0.31327576445128014, 1.0, 3 * np.pi / 8],
-            [-0.8080445791083615, 1.0, 5 * np.pi / 8],
-            [-0.3101398980494395, 1.0, 7 * np.pi / 8],
-            [-2.8111804455102014, 1.0, -np.pi / 8],
-            [-0.31327576445128014, 1.0, -3 * np.pi / 8],
-            [0.8080445791083615, 1.0, -5 * np.pi / 8],
-            [0.3101398980494395, 1.0, -7 * np.pi / 8],
+            [2.8111804455102014, np.pi / 8],
+            [-2.8111804455102014, -np.pi / 8],
+            [0.31327576445128014, 3 * np.pi / 8],
+            [-0.31327576445128014, -3 * np.pi / 8],
+            [-0.8080445791083615, 5 * np.pi / 8],
+            [0.8080445791083615, -5 * np.pi / 8],
+            [-0.3101398980494395, 7 * np.pi / 8],
+            [0.3101398980494395, -7 * np.pi / 8],
         ]
 
-        generated_terms = get_shift_rule(frequencies)[0]
+        generated_terms = generate_shift_rule(frequencies)
 
-        assert all([all(np.isclose(generated_terms[i], correct_terms[i])) for i in range(n_terms)])
+        assert np.allclose(generated_terms.T, correct_terms)
 
     def test_eight_term_rule_non_equidistant_custom_shifts(self):
         """Tests the correct non-equidistant eight term shift rule is generated given the
@@ -107,19 +108,19 @@ class TestGetShiftRule:
 
         n_terms = 8
         correct_terms = [
-            [2.709571194594805, 1.0, np.pi / 13],
-            [-0.12914139932030527, 1.0, 3 * np.pi / 7],
-            [-0.3820906256032637, 1.0, 2 * np.pi / 3],
-            [0.436088184940856, 1.0, 3 * np.pi / 4],
-            [-2.709571194594805, 1.0, -np.pi / 13],
-            [0.12914139932030527, 1.0, -3 * np.pi / 7],
-            [0.3820906256032637, 1.0, -2 * np.pi / 3],
-            [-0.436088184940856, 1.0, -3 * np.pi / 4],
+            [2.709571194594805, np.pi / 13],
+            [-2.709571194594805, -np.pi / 13],
+            [-0.12914139932030527, 3 * np.pi / 7],
+            [0.12914139932030527, -3 * np.pi / 7],
+            [-0.3820906256032637, 2 * np.pi / 3],
+            [0.3820906256032637, -2 * np.pi / 3],
+            [0.436088184940856, 3 * np.pi / 4],
+            [-0.436088184940856, -3 * np.pi / 4],
         ]
 
-        generated_terms = get_shift_rule(frequencies, custom_shifts)[0]
+        generated_terms = generate_shift_rule(frequencies, custom_shifts)
 
-        assert all([all(np.isclose(generated_terms[i], correct_terms[i])) for i in range(n_terms)])
+        assert np.allclose(generated_terms.T, correct_terms)
 
     def test_non_integer_frequency_default_shifts(self):
         """Tests the correct four term shift rule is generated given non-integer frequencies."""
@@ -128,15 +129,15 @@ class TestGetShiftRule:
         n_terms = 4
 
         correct_terms = [
-            [0.2845177968644246, 1, 3 * np.pi / 4],
-            [-0.048815536468908745, 1, 9 * np.pi / 4],
-            [-0.2845177968644246, 1, -3 * np.pi / 4],
-            [0.048815536468908745, 1, -9 * np.pi / 4],
+            [0.2845177968644246, 3 * np.pi / 4],
+            [-0.2845177968644246, -3 * np.pi / 4],
+            [-0.048815536468908745, 9 * np.pi / 4],
+            [0.048815536468908745, -9 * np.pi / 4],
         ]
 
-        generated_terms = get_shift_rule(frequencies)[0]
+        generated_terms = generate_shift_rule(frequencies)
 
-        assert all([all(np.isclose(generated_terms[i], correct_terms[i])) for i in range(n_terms)])
+        assert np.allclose(generated_terms.T, correct_terms)
 
     def test_non_integer_frequency_custom_shifts(self):
         """Tests the correct four term shift rule is generated given non-integer frequencies using
@@ -151,17 +152,17 @@ class TestGetShiftRule:
         n_terms = 6
 
         correct_terms = [
-            [1.7548361197453346, 1, 0.7853981633974483],
-            [-0.8720240894718643, 1, 1.0471975511965976],
-            [0.016695190986336428, 1, 2.0943951023931953],
-            [-1.7548361197453346, 1, -0.7853981633974483],
-            [0.8720240894718643, 1, -1.0471975511965976],
-            [-0.016695190986336428, 1, -2.0943951023931953],
+            [1.7548361197453346, 0.7853981633974483],
+            [-1.7548361197453346, -0.7853981633974483],
+            [-0.8720240894718643, 1.0471975511965976],
+            [0.8720240894718643, -1.0471975511965976],
+            [0.016695190986336428, 2.0943951023931953],
+            [-0.016695190986336428, -2.0943951023931953],
         ]
 
-        generated_terms = get_shift_rule(frequencies, custom_shifts)[0]
+        generated_terms = generate_shift_rule(frequencies, custom_shifts)
 
-        assert all([all(np.isclose(generated_terms[i], correct_terms[i])) for i in range(n_terms)])
+        assert np.allclose(generated_terms.T, correct_terms)
 
     def test_near_singular_warning(self):
         """Tests a warning is raised if the determinant of the matrix to be inverted is near zero
@@ -170,7 +171,7 @@ class TestGetShiftRule:
         frequencies = (1, 2, 3, 4, 5, 67)
 
         with pytest.warns(None) as warnings:
-            get_shift_rule(frequencies)
+            generate_shift_rule(frequencies)
 
         raised_warning = False
         for warning in warnings:
@@ -178,3 +179,137 @@ class TestGetShiftRule:
                 raised_warning = True
 
         assert raised_warning
+
+    def test_second_order_two_term_shift_rule(self):
+        """Test that the second order shift rule is correct and
+        properly simplified"""
+        frequencies = (1,)
+        generated_terms = generate_shift_rule(frequencies, order=2)
+        correct_terms = [[-0.5, 0], [0.5, -np.pi]]
+        assert np.allclose(generated_terms.T, correct_terms)
+
+    def test_second_order_two_term_shift_rule_custom_shifts(self):
+        """Test that the second order shift rule is correct and
+        properly simplified when custom shift values are provided"""
+        frequencies = (1,)
+        generated_terms = generate_shift_rule(frequencies, shifts=(np.pi / 4,), order=2)
+        correct_terms = [[-1, 0], [0.5, -np.pi / 2], [0.5, np.pi / 2]]
+        assert np.allclose(generated_terms.T, correct_terms)
+
+    def test_second_order_four_term_shift_rule(self):
+        """Test that the second order shift rule is correct and
+        properly simplified for generators with 4-term rules"""
+        frequencies = (0.5, 1)
+        generated_terms = generate_shift_rule(frequencies, order=2)
+        correct_terms = [
+            [-0.375, 0],
+            [0.25, -np.pi],
+            [0.25, np.pi],
+            [-0.125, -2 * np.pi],
+        ]
+        assert np.allclose(generated_terms.T, correct_terms)
+
+    def test_second_order_non_equidistant_shift_rule(self):
+        """Test that the second order shift rule is correct and
+        properly simplified for generators with non-equidistant frequencies"""
+        frequencies = (2, 3)
+        generated_terms = generate_shift_rule(frequencies, order=2)
+        correct_terms = [
+            [-6, 0],
+            [3.91421356, -np.pi / 4],
+            [3.91421356, np.pi / 4],
+            [-1, -np.pi / 2],
+            [-1, np.pi / 2],
+            [0.08578644, -3 * np.pi / 4],
+            [0.08578644, 3 * np.pi / 4],
+        ]
+        assert np.allclose(generated_terms.T, correct_terms)
+
+
+class TestMultiShiftRule:
+    """Tests for the generate_multi_shift_rule function"""
+
+    def test_single_parameter(self):
+        """Test that the generate_multi_shift_rule function
+        correctly returns a single-parameter shift rule"""
+        res = generate_multi_shift_rule([(1,)]).T
+        expected = [[0.5, np.pi / 2], [-0.5, -np.pi / 2]]
+        assert np.allclose(res, expected)
+
+        res = generate_multi_shift_rule([(1,)], orders=[2]).T
+        expected = [[-0.5, 0], [0.5, -np.pi]]
+        assert np.allclose(res, expected)
+
+        res = generate_multi_shift_rule([(1,)], orders=[2], shifts=[(np.pi / 4,)]).T
+        expected = [[-1, 0], [0.5, -np.pi / 2], [0.5, np.pi / 2]]
+        assert np.allclose(res, expected)
+
+    def test_two_single_frequency(self):
+        """Test that two independent single-frequency parameters
+        are correctly combined."""
+        res = generate_multi_shift_rule([(1,), (1,)]).T
+        expected = [
+            [0.25, np.pi / 2, np.pi / 2],
+            [-0.25, np.pi / 2, -np.pi / 2],
+            [-0.25, -np.pi / 2, np.pi / 2],
+            [0.25, -np.pi / 2, -np.pi / 2],
+        ]
+        assert np.allclose(res, expected)
+
+    def test_three_single_frequency(self):
+        """Test that three independent single-frequency parameters
+        are correctly combined."""
+        res = generate_multi_shift_rule([(1,), (1,), (1,)]).T
+        expected = [
+            [0.125, np.pi / 2, np.pi / 2, np.pi / 2],
+            [-0.125, np.pi / 2, np.pi / 2, -np.pi / 2],
+            [-0.125, np.pi / 2, -np.pi / 2, np.pi / 2],
+            [0.125, np.pi / 2, -np.pi / 2, -np.pi / 2],
+            [-0.125, -np.pi / 2, np.pi / 2, np.pi / 2],
+            [0.125, -np.pi / 2, np.pi / 2, -np.pi / 2],
+            [0.125, -np.pi / 2, -np.pi / 2, np.pi / 2],
+            [-0.125, -np.pi / 2, -np.pi / 2, -np.pi / 2],
+        ]
+        assert np.allclose(res, expected)
+
+    def test_two_frequency(self):
+        """Test that two independent 2-frequency parameters
+        are correctly combined."""
+        c1 = (np.sqrt(2) + 1) / (4 * np.sqrt(2))
+        c2 = (np.sqrt(2) - 1) / (4 * np.sqrt(2))
+        f = [(1,), (0.5, 1)]
+
+        res = generate_multi_shift_rule(f).T
+        expected = [
+            [c1 * 0.5, np.pi / 2, np.pi / 2],
+            [-c1 * 0.5, np.pi / 2, -np.pi / 2],
+            [-c2 * 0.5, np.pi / 2, 3 * np.pi / 2],
+            [c2 * 0.5, np.pi / 2, -3 * np.pi / 2],
+            [-c1 * 0.5, -np.pi / 2, np.pi / 2],
+            [c1 * 0.5, -np.pi / 2, -np.pi / 2],
+            [c2 * 0.5, -np.pi / 2, 3 * np.pi / 2],
+            [-c2 * 0.5, -np.pi / 2, -3 * np.pi / 2],
+        ]
+        assert np.allclose(res, expected)
+
+
+class TestEigvalsToFrequency:
+    """Tests for the eigvals_to_frequencies function"""
+
+    def test_two_eigvals(self):
+        """Test the case of two eigenvalues"""
+        res = qml.gradients.eigvals_to_frequencies((-0.5, 0.5))
+        expected = (1,)
+        assert res == expected
+
+    def test_four_eigvals(self):
+        """Test the case of four eigenvalues"""
+        res = qml.gradients.eigvals_to_frequencies((0.5, -0.5, 0, 0))
+        expected = (0.5, 1)
+        assert res == expected
+
+    def test_nonequidistant_eigvals(self):
+        """Test the case of non-equidistant eigenvalues"""
+        res = qml.gradients.eigvals_to_frequencies((0.453, 0.65, -1.2, 0))
+        expected = (0.453, 1.2, 1.85, 1.653, 0.65, 0.197)
+        assert res == expected
diff --git a/tests/interfaces/test_batch_autograd.py b/tests/interfaces/test_batch_autograd.py
index 63e3818c3..0b0f72ad1 100644
--- a/tests/interfaces/test_batch_autograd.py
+++ b/tests/interfaces/test_batch_autograd.py
@@ -300,7 +300,7 @@ class TestCaching:
 
         expected_runs_ideal = 1  # forward pass
         expected_runs_ideal += 2 * N  # Jacobian
-        expected_runs_ideal += 2 * N + 1  # Hessian diagonal
+        expected_runs_ideal += N + 1  # Hessian diagonal
         expected_runs_ideal += 4 * N * (N - 1) // 2  # Hessian off-diagonal
         assert dev.num_executions == expected_runs_ideal
         assert expected_runs_ideal < expected_runs
diff --git a/tests/interfaces/test_batch_tensorflow.py b/tests/interfaces/test_batch_tensorflow.py
index ea1b09586..ba77353ae 100644
--- a/tests/interfaces/test_batch_tensorflow.py
+++ b/tests/interfaces/test_batch_tensorflow.py
@@ -267,7 +267,7 @@ class TestCaching:
 
         expected_runs_ideal = 1  # forward pass
         expected_runs_ideal += 2 * N  # Jacobian
-        expected_runs_ideal += 2 * N + 1  # Hessian diagonal
+        expected_runs_ideal += N + 1  # Hessian diagonal
         expected_runs_ideal += 4 * N * (N - 1) // 2  # Hessian off-diagonal
         assert dev.num_executions == expected_runs_ideal
         assert expected_runs_ideal < nonideal_runs
diff --git a/tests/interfaces/test_batch_torch.py b/tests/interfaces/test_batch_torch.py
index 72d8240e1..d27522049 100644
--- a/tests/interfaces/test_batch_torch.py
+++ b/tests/interfaces/test_batch_torch.py
@@ -278,7 +278,7 @@ class TestCaching:
 
         expected_runs_ideal = 1  # forward pass
         expected_runs_ideal += 2 * N  # Jacobian
-        expected_runs_ideal += 2 * N + 1  # Hessian diagonal
+        expected_runs_ideal += N + 1  # Hessian diagonal
         expected_runs_ideal += 4 * N * (N - 1) // 2  # Hessian off-diagonal
         assert dev.num_executions == expected_runs_ideal
         assert expected_runs_ideal < expected_runs