Repository: PennyLaneAI/pennylane
Pull request #1486 (2021-07-30 13:42:47)

 (non)Differentiable CV parameter-shift gradient transform
**Context:** As part of the roadmap for supporting differentiable batch execution, gradient logic will be moved out of the subclasses and into a module of pure functions. This is the third PR after #1476 and #1479; here, we move the CV parameter-shift logic out of `CVParamShiftTape` and into a new gradients package.

**Description of the Change:**

- Adds a function `second_order_param_shift` for computing the second-order CV parameter-shift gradient of a tape terminating in expectation values. Port of `CVQubitParamShift.parameter_shift_second_order`.

- Adds a function `var_param_shift` for computing the CV parameter-shift gradient of a tape terminating in one or more variances. Directly equivalent to `CVQubitParamShift.parameter_shift_var`.

- Adds a wrapper function `param_shift_cv`, directly analagous to `CVQubitParamShift.parameter_shift`, which does the following:

  - Basic input validation
  - Performs static analysis of the tape parameters, to determine which support parameter-shift rules.
  - A fallback mode; `finite_diff` is called for any unsupported parameters
  - Dispatches to one of `expval_param_shift`, `var_param_shift`, or `second_order_param_shift` depending on the structure of the tape
  - Provides a processing function to combine the gradients computed via fallback and the parameter-shift methods.

**Benefits:**

- The parameter-shift logic is now much more user and dev accessible

**Possible Drawbacks:**

- See drawbacks in #1476 
- This transforms is _not_ differentiable, as it depends on a lot of legacy CV code (such as Heisenberg transform functions) which are written in vanilla NumPy and not differentiable.
- The size of the PR - there is very little new code here, it is all a port of `CVParamShiftTape`, including the tests!

**Related GitHub Issues:**


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [.github/CHANGELOG.md](/PennyLaneAI/pennylane/blob/master/.github/CHANGELOG.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1486](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (9378d59) into [parameter-shift](https://codecov.io/gh/PennyLaneAI/pennylane/commit/0c5666c7bdb3bb390ca967660d272953973eed1b?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (0c5666c) will **decrease** coverage by `0.00%`.
> The diff coverage is `98.18%`.

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

```diff
@@                 Coverage Diff                 @@
##           parameter-shift    #1486      +/-   ##
===================================================
- Coverage            98.32%   98.32%   -0.01%     
===================================================
  Files                  173      181       +8     
  Lines                12737    13027     +290     
===================================================
+ Hits                 12524    12809     +285     
- Misses                 213      218       +5     
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/tape/jacobian\_tape.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RhcGUvamFjb2JpYW5fdGFwZS5weQ==) | `98.14% <ø> (-0.02%)` | :arrow_down: |
| [pennylane/gradients/parameter\_shift\_cv.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyYWRpZW50cy9wYXJhbWV0ZXJfc2hpZnRfY3YucHk=) | `98.15% <98.15%> (ø)` | |
| [pennylane/gradients/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/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%> (ø)` | |
| [pennylane/gradients/finite\_difference.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyYWRpZW50cy9maW5pdGVfZGlmZmVyZW5jZS5weQ==) | `100.00% <100.00%> (ø)` | |
| [pennylane/ops/\_\_init\_\_.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9fX2luaXRfXy5weQ==) | `100.00% <0.00%> (ø)` | |
| [pennylane/gradients/parameter\_shift.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL2dyYWRpZW50cy9wYXJhbWV0ZXJfc2hpZnQucHk=) | `99.35% <0.00%> (ø)` | |
| [pennylane/transforms/hamiltonian\_expand.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RyYW5zZm9ybXMvaGFtaWx0b25pYW5fZXhwYW5kLnB5) | `100.00% <0.00%> (ø)` | |
| [pennylane/ops/qubit.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC5weQ==) | | |
| [pennylane/ops/qubit/observables.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9xdWJpdC9vYnNlcnZhYmxlcy5weQ==) | `100.00% <0.00%> (ø)` | |
| ... and [9 more](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486/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/1486?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/1486?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [0c5666c...9378d59](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1486?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: [ch7888]

Patch:
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 121168f70..66ceb552f 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -75,10 +75,14 @@
 * A new gradients module `qml.gradients` has been added, which provides
   differentiable quantum gradient transforms.
   [(#1476)](https://github.com/PennyLaneAI/pennylane/pull/1476)
+  [(#1479)](https://github.com/PennyLaneAI/pennylane/pull/1479)
+  [(#1486)](https://github.com/PennyLaneAI/pennylane/pull/1486)
 
   Available quantum gradient transforms include:
 
   - `qml.gradients.finite_diff`
+  - `qml.gradients.param_shift`
+  - `qml.gradients.param_shift_cv`
 
   For example,
 
diff --git a/pennylane/gradients/__init__.py b/pennylane/gradients/__init__.py
index 105201a17..2db236e9d 100644
--- a/pennylane/gradients/__init__.py
+++ b/pennylane/gradients/__init__.py
@@ -16,6 +16,8 @@ import pennylane as qml
 
 from . import finite_difference
 from . import parameter_shift
+from . import parameter_shift_cv
 
 from .finite_difference import finite_diff, finite_diff_coeffs, generate_shifted_tapes
 from .parameter_shift import param_shift
+from .parameter_shift_cv import param_shift_cv
diff --git a/pennylane/gradients/finite_difference.py b/pennylane/gradients/finite_difference.py
index cc91ba7fa..f823d761f 100644
--- a/pennylane/gradients/finite_difference.py
+++ b/pennylane/gradients/finite_difference.py
@@ -161,7 +161,7 @@ def generate_shifted_tapes(tape, idx, shifts, multipliers=None):
         ``idx`` shifted by consecutive values of ``shift``. The length
         of the returned list of tapes will match the length of ``shifts``.
     """
-    params = tape.get_parameters()
+    params = list(tape.get_parameters())
     tapes = []
 
     for i, s in enumerate(shifts):
@@ -300,7 +300,7 @@ def finite_diff(tape, argnum=None, h=1e-7, approx_order=1, n=1, strategy="forwar
         # In the future, we might want to change this so that only tuples
         # of arrays are returned.
         for i, g in enumerate(grads):
-            g = qml.math.convert_like(g, res[0])
+            g = qml.math.convert_like(g, results[0])
             if hasattr(g, "dtype") and g.dtype is np.dtype("object"):
                 grads[i] = qml.math.hstack(g)
 
diff --git a/pennylane/gradients/parameter_shift_cv.py b/pennylane/gradients/parameter_shift_cv.py
new file mode 100644
index 000000000..5db6b6b45
--- /dev/null
+++ b/pennylane/gradients/parameter_shift_cv.py
@@ -0,0 +1,693 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+This module contains functions for computing the parameter-shift gradient
+of a CV-based quantum tape.
+"""
+# pylint: disable=protected-access,too-many-arguments,too-many-statements,too-many-branches
+import itertools
+import warnings
+
+import numpy as np
+
+import pennylane as qml
+
+from .finite_difference import finite_diff, generate_shifted_tapes
+from .parameter_shift import expval_param_shift, _get_operation_recipe, _process_gradient_recipe
+
+
+def _grad_method(tape, idx):
+    """Determine the best CV parameter-shift gradient recipe for a given
+    parameter index of a tape.
+
+    Args:
+        tape (.QuantumTape): input tape
+        idx (int): positive integer corresponding to the parameter location
+            on the tape to inspect
+
+    Returns:
+        str: a string containing either ``"A"`` (for first-order analytic method),
+            ``"A2"`` (second-order analytic method), ``"F"`` (finite differences),
+            or ``"0"`` (constant parameter).
+    """
+
+    op = tape._par_info[idx]["op"]
+
+    if op.grad_method in (None, "F"):
+        return op.grad_method
+
+    if op.grad_method != "A":
+        raise ValueError(f"Operation {op} has unknown gradient method {op.grad_method}")
+
+    # Operation supports the CV parameter-shift rule.
+    # Create an empty list to store the 'best' partial derivative method
+    # for each observable
+    best = []
+
+    for m in tape.measurements:
+
+        if (m.return_type is qml.operation.Probability) or (m.obs.ev_order not in (1, 2)):
+            # Higher-order observables (including probability) only support finite differences.
+            best.append("F")
+            continue
+
+        # get the set of operations betweens the operation and the observable
+        ops_between = tape.graph.nodes_between(op, m.obs)
+
+        if not ops_between:
+            # if there is no path between the operation and the observable,
+            # the operator has a zero gradient.
+            best.append("0")
+            continue
+
+        # For parameter-shift compatible CV gates, we need to check both the
+        # intervening gates, and the type of the observable.
+        best_method = "A"
+
+        if any(not k.supports_heisenberg for k in ops_between):
+            # non-Gaussian operators present in-between the operation
+            # and the observable. Must fallback to numeric differentiation.
+            best_method = "F"
+
+        elif m.obs.ev_order == 2:
+
+            if m.return_type is qml.operation.Expectation:
+                # If the observable is second-order, we must use the second-order
+                # CV parameter shift rule
+                best_method = "A2"
+
+            elif m.return_type is qml.operation.Variance:
+                # we only support analytic variance gradients for
+                # first-order observables
+                best_method = "F"
+
+        best.append(best_method)
+
+    if all(k == "0" for k in best):
+        # if the operation is independent of *all* observables
+        # in the circuit, the gradient will be 0
+        return "0"
+
+    if "F" in best:
+        # one non-analytic observable path makes the whole operation
+        # gradient method fallback to finite-difference
+        return "F"
+
+    if "A2" in best:
+        # one second-order observable makes the whole operation gradient
+        # require the second-order parameter-shift rule
+        return "A2"
+
+    return "A"
+
+
+def _gradient_analysis(tape):
+    """Update the parameter information dictionary of the tape with
+    gradient information of each parameter."""
+
+    if getattr(tape, "_gradient_fn", None) is param_shift_cv:
+        # gradient analysis has already been performed on this tape
+        return
+
+    tape._gradient_fn = param_shift_cv
+
+    for idx, info in tape._par_info.items():
+        info["grad_method"] = _grad_method(tape, idx)
+
+
+def _transform_observable(obs, Z, device_wires):
+    """Apply a Gaussian linear transformation to an observable.
+
+    Args:
+        obs (.Observable): observable to transform
+        Z (array[float]): Heisenberg picture representation of the linear transformation
+        device_wires (.Wires): wires on the device the transformed observable is to be
+            measured on
+
+    Returns:
+        .Observable: the transformed observable
+    """
+    # Get the Heisenberg representation of the observable
+    # in the position/momentum basis. The returned matrix/vector
+    # will have been expanded to act on the entire device.
+    if obs.ev_order > 2:
+        raise NotImplementedError("Transforming observables of order > 2 not implemented.")
+
+    A = obs.heisenberg_obs(device_wires)
+
+    if A.ndim != obs.ev_order:
+        raise ValueError(
+            "Mismatch between the polynomial order of observable and its Heisenberg representation"
+        )
+
+    # transform the observable by the linear transformation Z
+    A = A @ Z
+
+    if A.ndim == 2:
+        A = A + A.T
+
+    # TODO: if the A matrix corresponds to a known observable in PennyLane,
+    # for example qml.X, qml.P, qml.NumberOperator, we should return that
+    # instead. This will allow for greater device compatibility.
+    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):
+    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.
+
+    Expectation values may be of first- or second-order observables,
+    but variances can only be taken of first-order variables.
+
+    .. warning::
+
+        This method can only be executed on devices that support the
+        :class:`~.PolyXP` observable.
+
+    Args:
+        tape (.QuantumTape): quantum tape to differentiate
+        dev_wires (.Wires): wires on the device the parameter-shift method is computed on
+        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``.
+        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.
+        f0 (tensor_like[float] or None): Output of the evaluated input tape. If provided,
+            and the gradient recipe contains an unshifted term, this value is used,
+            saving a quantum evaluation.
+
+    Returns:
+        tuple[list[QuantumTape], function]: A tuple containing a
+        list of generated tapes, in addition to a post-processing
+        function to be applied to the evaluated tapes.
+    """
+    argnum = argnum or tape.trainable_params
+
+    # Determine the locations of any variance measurements in the measurement queue.
+    var_mask = [m.return_type is qml.operation.Variance for m in tape.measurements]
+    var_idx = np.where(var_mask)[0]
+
+    # Get <A>, the expectation value of the tape with unshifted parameters.
+    expval_tape = tape.copy(copy_operations=True)
+
+    # Convert all variance measurements on the tape into expectation values
+    for i in var_idx:
+        obs = expval_tape._measurements[i].obs
+        expval_tape._measurements[i] = qml.measure.MeasurementProcess(
+            qml.operation.Expectation, obs=obs
+        )
+
+    gradient_tapes = [expval_tape]
+
+    # evaluate the analytic derivative of <A>
+    pdA_tapes, pdA_fn = expval_param_shift(expval_tape, argnum, shift, gradient_recipes, f0)
+    gradient_tapes.extend(pdA_tapes)
+
+    # Store the number of first derivative tapes, so that we know
+    # the number of results to post-process later.
+    tape_boundary = len(pdA_tapes) + 1
+    expval_sq_tape = tape.copy(copy_operations=True)
+
+    for i in var_idx:
+        # We need to calculate d<A^2>/dp; to do so, we replace the
+        # observables A in the queue with A^2.
+        obs = expval_sq_tape._measurements[i].obs
+
+        # CV first-order observable
+        # get the heisenberg representation
+        # This will be a real 1D vector representing the
+        # first-order observable in the basis [I, x, p]
+        A = obs._heisenberg_rep(obs.parameters)
+
+        # take the outer product of the heisenberg representation
+        # with itself, to get a square symmetric matrix representing
+        # the square of the observable
+        obs = qml.PolyXP(np.outer(A, A), wires=obs.wires)
+        expval_sq_tape._measurements[i] = qml.measure.MeasurementProcess(
+            qml.operation.Expectation, obs=obs
+        )
+
+    # Non-involutory observables are present; the partial derivative of <A^2>
+    # 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
+    )
+    gradient_tapes.extend(pdA2_tapes)
+
+    def processing_fn(results):
+        mask = qml.math.convert_like(qml.math.reshape(var_mask, [-1, 1]), results[0])
+        f0 = qml.math.expand_dims(results[0], -1)
+
+        pdA = pdA_fn(results[1:tape_boundary])
+        pdA2 = pdA2_fn(results[tape_boundary:])
+
+        # return d(var(A))/dp = d<A^2>/dp -2 * <A> * d<A>/dp for the variances (mask==True)
+        # d<A>/dp for plain expectations (mask==False)
+        return qml.math.where(mask, pdA2 - 2 * f0 * pdA, pdA)
+
+    return gradient_tapes, processing_fn
+
+
+def second_order_param_shift(tape, dev_wires, argnum=None, shift=np.pi / 2, 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.
+
+    .. note::
+
+        The 2nd order method can handle also first-order observables, but
+        1st order method may be more efficient unless it's really easy to
+        experimentally measure arbitrary 2nd order observables.
+
+    .. warning::
+
+        The 2nd order method can only be executed on devices that support the
+        :class:`~.PolyXP` observable.
+
+    Args:
+        tape (.QuantumTape): quantum tape to differentiate
+        dev_wires (.Wires): wires on the device the parameter-shift method is computed on
+        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``.
+        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.
+
+    Returns:
+        tuple[list[QuantumTape], function]: A tuple containing a
+        list of generated tapes, in addition to a post-processing
+        function to be applied to the evaluated tapes.
+    """
+    argnum = argnum or list(tape.trainable_params)
+    gradient_recipes = gradient_recipes or [None] * len(argnum)
+
+    gradient_tapes = []
+    shapes = []
+    obs_indices = []
+    gradient_values = []
+
+    for idx, _ in enumerate(tape.trainable_params):
+        t_idx = list(tape.trainable_params)[idx]
+        op = tape._par_info[t_idx]["op"]
+
+        if idx not in argnum:
+            # parameter has zero gradient
+            shapes.append(0)
+            obs_indices.append([])
+            gradient_values.append([])
+            continue
+
+        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:
+            # The 2nd order CV parameter-shift rule only accepts two-term shifts
+            raise NotImplementedError(
+                "Taking the analytic gradient for order-2 operators is "
+                f"unsupported for operation {op} which has a "
+                "gradient recipe of more than two terms."
+            )
+
+        shifted_tapes = generate_shifted_tapes(tape, idx, shifts, multipliers)
+
+        # evaluate transformed observables at the original parameter point
+        # first build the Heisenberg picture transformation matrix Z
+        Z0 = op.heisenberg_tr(dev_wires, inverse=True)
+        Z2 = shifted_tapes[0]._par_info[t_idx]["op"].heisenberg_tr(dev_wires)
+        Z1 = shifted_tapes[1]._par_info[t_idx]["op"].heisenberg_tr(dev_wires)
+
+        # derivative of the operation
+        Z = Z2 * coeffs[0] + Z1 * coeffs[1]
+        Z = Z @ Z0
+
+        # conjugate Z with all the descendant operations
+        B = np.eye(1 + 2 * len(dev_wires))
+        B_inv = B.copy()
+
+        succ = tape.graph.descendants_in_order((op,))
+        operation_descendents = itertools.filterfalse(qml.circuit_graph._is_observable, succ)
+        observable_descendents = filter(qml.circuit_graph._is_observable, succ)
+
+        for BB in operation_descendents:
+            if not BB.supports_heisenberg:
+                # if the descendant gate is non-Gaussian in parameter-shift differentiation
+                # mode, then there must be no observable following it.
+                continue
+
+            B = BB.heisenberg_tr(dev_wires) @ B
+            B_inv = B_inv @ BB.heisenberg_tr(dev_wires, inverse=True)
+
+        Z = B @ Z @ B_inv  # conjugation
+
+        g_tape = tape.copy(copy_operations=True)
+        constants = []
+
+        # transform the descendant observables into their derivatives using Z
+        transformed_obs_idx = []
+
+        for obs in observable_descendents:
+            # get the index of the descendent observable
+            idx = tape.observables.index(obs)
+            transformed_obs_idx.append(idx)
+
+            transformed_obs = _transform_observable(obs, Z, dev_wires)
+
+            A = transformed_obs.parameters[0]
+            constant = None
+
+            # Check if the transformed observable corresponds to a constant term.
+            if len(A.nonzero()[0]) == 1:
+                if A.ndim == 2 and A[0, 0] != 0:
+                    constant = A[0, 0]
+
+                elif A.ndim == 1 and A[0] != 0:
+                    constant = A[0]
+
+            constants.append(constant)
+
+            g_tape._measurements[idx] = qml.measure.MeasurementProcess(
+                qml.operation.Expectation, _transform_observable(obs, Z, dev_wires)
+            )
+
+        if not any(i is None for i in constants):
+            # Check if *all* transformed observables corresponds to a constant term.
+            # term. If this is the case for all transformed observables on the tape,
+            # then <psi|A|psi> = A<psi|psi> = A,
+            # and we can avoid the device execution.
+            shapes[-1] = 0
+            obs_indices.append(transformed_obs_idx)
+            gradient_values.append(constants)
+            continue
+
+        gradient_tapes.append(g_tape)
+        obs_indices.append(transformed_obs_idx)
+        gradient_values.append(None)
+
+    def processing_fn(results):
+        grads = []
+        start = 0
+
+        if not results:
+            results = [np.zeros([tape.output_dim])]
+
+        interface = qml.math.get_interface(results[0])
+        iterator = enumerate(zip(shapes, gradient_values, obs_indices))
+
+        for i, (shape, grad_value, obs_ind) in iterator:
+
+            if shape == 0:
+                # parameter has zero gradient
+                g = qml.math.zeros_like(results[0], like=interface)
+
+                if grad_value:
+                    g = qml.math.scatter_element_add(g, obs_ind, grad_value, like=interface)
+
+                grads.append(g)
+                continue
+
+            obs_result = results[start : start + shape]
+            start = start + shape
+
+            # compute the linear combination of results and coefficients
+            obs_result = qml.math.stack(obs_result[0])
+            g = qml.math.zeros_like(obs_result, like=interface)
+
+            if qml.math.get_interface(g) not in ("tensorflow", "autograd"):
+                obs_ind = (obs_ind,)
+
+            g = qml.math.scatter_element_add(g, obs_ind, obs_result[obs_ind], like=interface)
+            grads.append(g)
+
+        # The following is for backwards compatibility; currently,
+        # the device stacks multiple measurement arrays, even if not the same
+        # size, resulting in a ragged array.
+        # In the future, we might want to change this so that only tuples
+        # of arrays are returned.
+        for i, g in enumerate(grads):
+            g = qml.math.convert_like(g, results[0])
+            if hasattr(g, "dtype") and g.dtype is np.dtype("object"):
+                grads[i] = qml.math.hstack(g)
+
+        return qml.math.T(qml.math.stack(grads))
+
+    return gradient_tapes, processing_fn
+
+
+def param_shift_cv(
+    tape,
+    dev,
+    argnum=None,
+    shift=np.pi / 2,
+    gradient_recipes=None,
+    fallback_fn=finite_diff,
+    f0=None,
+    force_order2=False,
+):
+    r"""Generate the CV parameter-shift tapes and postprocessing methods required
+    to compute the gradient of a gate parameter with respect to the CV output.
+
+    Args:
+        tape (.QuantumTape): quantum tape to differentiate
+        dev (.Device): device the parameter-shift method is to be computed on
+        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``.
+        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.
+
+            This is a tuple with one nested list per 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
+
+            .. math:: \frac{\partial}{\partial\phi_k}f = \sum_{i} c_i f(a_i \phi_k + s_i).
+
+            If ``None``, the default gradient recipe containing the two terms
+            :math:`[c_0, a_0, s_0]=[1/2, 1, \pi/2]` and :math:`[c_1, a_1,
+            s_1]=[-1/2, 1, -\pi/2]` is assumed for every parameter.
+        fallback_fn (None or Callable): a fallback grdient function to use for
+            any parameters that do not support the parameter-shift rule.
+        f0 (tensor_like[float] or None): Output of the evaluated input tape. If provided,
+            and the gradient recipe contains an unshifted term, this value is used,
+            saving a quantum evaluation.
+        force_order2 (bool): if True, use the order-2 method even if not necessary
+
+    Returns:
+        tuple[list[QuantumTape], function]: A tuple containing a
+        list of generated tapes, in addition to a post-processing
+        function to be applied to the evaluated tapes.
+
+    This transform supports analytic gradients of Gaussian CV operations using
+    the parameter-shift rule. This gradient method returns *exact* gradients,
+    and can be computed directly on quantum hardware.
+
+    Analytic gradients of photonic circuits that satisfy
+    the following constraints with regards to measurements are supported:
+
+    * Expectation values are restricted to observables that are first- and
+      second-order in :math:`\hat{x}` and :math:`\hat{p}` only.
+      This includes :class:`~.X`, :class:`~.P`, :class:`~.QuadOperator`,
+      :class:`~.PolyXP`, and :class:`~.NumberOperator`.
+
+      For second-order observables, the device **must support** :class:`~.PolyXP`.
+
+    * Variances are restricted to observables that are first-order
+      in :math:`\hat{x}` and :math:`\hat{p}` only. This includes :class:`~.X`, :class:`~.P`,
+      :class:`~.QuadOperator`, and *some* parameter values of :class:`~.PolyXP`.
+
+      The device **must support** :class:`~.PolyXP`.
+
+    .. warning::
+
+        Fock state probabilities (tapes that return :func:`~pennylane.probs` or
+        expectation values of :class:`~.FockStateProjector`) are not supported.
+
+    In addition, the tape operations must fulfill the following requirements:
+
+    * Only Gaussian operations are differentiable.
+
+    * Non-differentiable Fock states and Fock operations may *precede* all differentiable Gaussian,
+      operations. For example, the following is permissible:
+
+      .. code-block:: python
+
+          with qml.tape.JacobianTape() as tape:
+              # Non-differentiable Fock operations
+              qml.FockState(2, wires=0)
+              qml.Kerr(0.654, wires=1)
+
+              # differentiable Gaussian operations
+              qml.Displacement(0.6, 0.5, wires=0)
+              qml.Beamsplitter(0.5, 0.1, wires=[0, 1])
+              qml.expval(qml.NumberOperator(0))
+
+          tape.trainable_params = {2, 3, 4}
+
+    * If a Fock operation succeeds a Gaussian operation, the Fock operation must
+      not contribute to any measurements. For example, the following is allowed:
+
+      .. code-block:: python
+
+          with qml.tape.JacobianTape() as tape:
+              qml.Displacement(0.6, 0.5, wires=0)
+              qml.Beamsplitter(0.5, 0.1, wires=[0, 1])
+              qml.Kerr(0.654, wires=1)  # there is no measurement on wire 1
+              qml.expval(qml.NumberOperator(0))
+
+          tape.trainable_params = {0, 1, 2}
+
+    If any of the above constraints are not followed, the tape cannot be differentiated
+    via the CV parameter-shift rule. Please use numerical differentiation instead.
+
+    **Example**
+
+    >>> r0, phi0, r1, phi1 = [0.4, -0.3, -0.7, 0.2]
+    >>> dev = qml.device("default.gaussian", wires=1)
+    >>> with qml.tape.JacobianTape() as tape:
+    ...     qml.Squeezing(r0, phi0, wires=[0])
+    ...     qml.Squeezing(r1, phi1, wires=[0])
+    ...     qml.expval(qml.NumberOperator(0))  # second-order
+    >>> tape.trainable_params = {0, 2}
+    >>> gradient_tapes, fn = qml.gradients.param_shift_cv(tape, dev)
+    >>> res = dev.batch_execute(gradient_tapes)
+    >>> fn(res)
+    array([[-0.32487113, -0.87049853]])
+    """
+
+    # 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 = []
+    shapes = []
+    fns = []
+
+    def _update(data):
+        """Utility function to update the list of gradient tapes,
+        the corresponding number of gradient tapes, and the processing functions"""
+        gradient_tapes.extend(data[0])
+        shapes.append(len(data[0]))
+        fns.append(data[1])
+
+    # TODO: replace the JacobianTape._grad_method_validation
+    # functionality before deprecation.
+    diff_methods = tape._grad_method_validation("analytic" if fallback_fn is None else "best")
+    all_params_grad_method_zero = all(g == "0" for g in diff_methods)
+
+    if not tape.trainable_params or all_params_grad_method_zero:
+        return gradient_tapes, 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))
+    var_present = any(m.return_type is qml.operation.Variance for m in tape.measurements)
+
+    unsupported_params = []
+    first_order_params = []
+    second_order_params = []
+
+    for idx, g in method_map.items():
+        if g == "F":
+            unsupported_params.append(idx)
+
+        elif g == "A":
+            first_order_params.append(idx)
+
+        elif g == "A2":
+            second_order_params.append(idx)
+
+    if force_order2:
+        # all analytic parameters should be computed using the second-order method
+        second_order_params += first_order_params
+        first_order_params = []
+
+    if "PolyXP" not in dev.observables and (second_order_params or var_present):
+        warnings.warn(
+            f"The device {dev.short_name} does not support "
+            "the PolyXP observable. The analytic parameter-shift cannot be used for "
+            "second-order observables; falling back to finite-differences.",
+            UserWarning,
+        )
+
+        if var_present:
+            unsupported_params += first_order_params
+            first_order_params = []
+
+        unsupported_params += second_order_params
+        second_order_params = []
+
+    # If there are unsupported operations, call the fallback gradient function
+    if unsupported_params:
+        _update(fallback_fn(tape, argnum=unsupported_params))
+
+    # collect all the analytic parameters
+    argnum = first_order_params + second_order_params
+
+    if not argnum:
+        # No analytic parameters. Return the existing fallback tapes/fn
+        return gradient_tapes, fns[-1]
+
+    gradient_recipes = gradient_recipes or [None] * len(argnum)
+
+    if var_present:
+        _update(var_param_shift(tape, dev.wires, argnum, shift, 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))
+
+        if second_order_params:
+            _update(
+                second_order_param_shift(
+                    tape, dev.wires, second_order_params, shift, gradient_recipes
+                )
+            )
+
+    def processing_fn(results):
+        start = 0
+        grads = []
+
+        for s, f in zip(shapes, fns):
+            grads.append(f(results[start : start + s]))
+            start += s
+
+        return sum(grads)
+
+    return gradient_tapes, processing_fn
diff --git a/pennylane/operation.py b/pennylane/operation.py
index 99168eac5..5fdc762f3 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -1692,7 +1692,7 @@ class CVOperation(CV, Operation):
         Returns:
             array[float]: :math:`\tilde{U}`, the Heisenberg picture representation of the linear transformation
         """
-        p = self.parameters
+        p = [qml.math.toarray(a) for a in self.parameters]
         if inverse:
             if self.par_domain == "A":
                 # TODO: expand this for the new par domain class, for non-unitary matrices.
diff --git a/pennylane/tape/jacobian_tape.py b/pennylane/tape/jacobian_tape.py
index 4054f9ba3..77ce708b2 100644
--- a/pennylane/tape/jacobian_tape.py
+++ b/pennylane/tape/jacobian_tape.py
@@ -399,7 +399,8 @@ class JacobianTape(QuantumTape):
         """
         raise NotImplementedError
 
-    def _choose_params_with_methods(self, diff_methods, argnum):
+    @staticmethod
+    def _choose_params_with_methods(diff_methods, argnum):
         """Chooses the trainable parameters to use for computing the Jacobian
         by returning a map of their indices and differentiation methods.
 
@@ -423,11 +424,6 @@ class JacobianTape(QuantumTape):
         if isinstance(argnum, int):
             argnum = [argnum]
 
-        if not all(ind in self.trainable_params for ind in argnum):
-            raise ValueError(
-                "Incorrect trainable parameters were specified for the argnum argument."
-            )
-
         num_params = len(argnum)
 
         if num_params == 0:


Test patch:
diff --git a/tests/gradients/test_parameter_shift_cv.py b/tests/gradients/test_parameter_shift_cv.py
new file mode 100644
index 000000000..c43548a6b
--- /dev/null
+++ b/tests/gradients/test_parameter_shift_cv.py
@@ -0,0 +1,1105 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for the gradients.parameter_shift_cv module."""
+import pytest
+
+import pennylane as qml
+from pennylane import numpy as np
+from pennylane.gradients import param_shift_cv
+from pennylane.gradients.parameter_shift_cv import (
+    _grad_method,
+    _gradient_analysis,
+    _transform_observable,
+)
+
+
+hbar = 2
+
+
+class TestGradAnalysis:
+    """Tests for parameter gradient methods"""
+
+    def test_non_differentiable(self):
+        """Test that a non-differentiable parameter is
+        correctly marked"""
+
+        with qml.tape.JacobianTape() as tape:
+            qml.FockState(1, wires=0)
+            qml.Displacement(0.543, 0, wires=[1])
+            qml.Beamsplitter(0, 0, wires=[0, 1])
+            qml.expval(qml.X(wires=[0]))
+
+        assert _grad_method(tape, 0) is None
+        assert _grad_method(tape, 1) == "A"
+        assert _grad_method(tape, 2) == "A"
+        assert _grad_method(tape, 3) == "A"
+        assert _grad_method(tape, 4) == "A"
+
+        _gradient_analysis(tape)
+
+        assert tape._par_info[0]["grad_method"] is None
+        assert tape._par_info[1]["grad_method"] == "A"
+        assert tape._par_info[2]["grad_method"] == "A"
+        assert tape._par_info[3]["grad_method"] == "A"
+        assert tape._par_info[4]["grad_method"] == "A"
+
+        _gradient_analysis(tape)
+
+    def test_independent(self):
+        """Test that an independent variable is properly marked
+        as having a zero gradient"""
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(0.543, wires=[0])
+            qml.Rotation(-0.654, wires=[1])
+            qml.expval(qml.P(0))
+
+        assert _grad_method(tape, 0) == "A"
+        assert _grad_method(tape, 1) == "0"
+
+        _gradient_analysis(tape)
+
+        assert tape._par_info[0]["grad_method"] == "A"
+        assert tape._par_info[1]["grad_method"] == "0"
+
+    def test_finite_diff(self, monkeypatch):
+        """If an op has grad_method=F, this should be respected
+        by the qml.tape.JacobianTape"""
+        monkeypatch.setattr(qml.Rotation, "grad_method", "F")
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(0.543, wires=[0])
+            qml.Squeezing(0.543, 0, wires=[0])
+            qml.expval(qml.P(0))
+
+        assert _grad_method(tape, 0) == "F"
+        assert _grad_method(tape, 1) == "A"
+        assert _grad_method(tape, 2) == "A"
+
+    def test_non_gaussian_operation(self):
+        """Test that a non-Gaussian operation succeeding
+        a differentiable Gaussian operation results in
+        numeric differentiation."""
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(1.0, wires=[0])
+            qml.Rotation(1.0, wires=[1])
+            # Non-Gaussian
+            qml.Kerr(1.0, wires=[1])
+            qml.expval(qml.P(0))
+            qml.expval(qml.X(1))
+
+        # First rotation gate has no succeeding non-Gaussian operation
+        assert _grad_method(tape, 0) == "A"
+        # Second rotation gate does no succeeding non-Gaussian operation
+        assert _grad_method(tape, 1) == "F"
+        # Kerr gate does not support the parameter-shift rule
+        assert _grad_method(tape, 2) == "F"
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(1.0, wires=[0])
+            qml.Rotation(1.0, wires=[1])
+            # entangle the modes
+            qml.Beamsplitter(1.0, 0.0, wires=[0, 1])
+            # Non-Gaussian
+            qml.Kerr(1.0, wires=[1])
+            qml.expval(qml.P(0))
+            qml.expval(qml.X(1))
+
+        # After entangling the modes, the Kerr gate now succeeds
+        # both initial rotations
+        assert _grad_method(tape, 0) == "F"
+        assert _grad_method(tape, 1) == "F"
+        assert _grad_method(tape, 2) == "F"
+
+    def test_probability(self):
+        """Probability is the expectation value of a
+        higher order observable, and thus only supports numerical
+        differentiation"""
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(0.543, wires=[0])
+            qml.Squeezing(0.543, 0, wires=[0])
+            qml.probs(wires=0)
+
+        assert _grad_method(tape, 0) == "F"
+        assert _grad_method(tape, 1) == "F"
+        assert _grad_method(tape, 2) == "F"
+
+    def test_variance(self):
+        """If the variance of the observable is first order, then
+        parameter-shift is supported. If the observable is second order,
+        however, only finite-differences is supported."""
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(1.0, wires=[0])
+            qml.var(qml.P(0))  # first order
+
+        assert _grad_method(tape, 0) == "A"
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(1.0, wires=[0])
+            qml.var(qml.NumberOperator(0))  # second order
+
+        assert _grad_method(tape, 0) == "F"
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(1.0, wires=[0])
+            qml.Rotation(1.0, wires=[1])
+            qml.Beamsplitter(0.5, 0.0, wires=[0, 1])
+            qml.var(qml.NumberOperator(0))  # fourth order
+            qml.expval(qml.NumberOperator(1))
+
+        assert _grad_method(tape, 0) == "F"
+        assert _grad_method(tape, 1) == "F"
+        assert _grad_method(tape, 2) == "F"
+        assert _grad_method(tape, 3) == "F"
+
+    def test_second_order_expectation(self):
+        """Test that the expectation of a second-order observable forces
+        the gradient method to use the second-order parameter-shift rule"""
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(1.0, wires=[0])
+            qml.expval(qml.NumberOperator(0))  # second order
+
+        assert _grad_method(tape, 0) == "A2"
+
+    def test_unknown_op_grad_method(self, monkeypatch):
+        """Test that an exception is raised if an operator has a
+        grad method defined that the CV parameter-shift tape
+        doesn't recognize"""
+        monkeypatch.setattr(qml.Rotation, "grad_method", "B")
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(1.0, wires=0)
+            qml.expval(qml.X(0))
+
+        with pytest.raises(ValueError, match="unknown gradient method"):
+            _grad_method(tape, 0)
+
+
+class TestTransformObservable:
+    """Tests for the _transform_observable method"""
+
+    def test_incorrect_heisenberg_size(self, monkeypatch):
+        """The number of dimensions of a CV observable Heisenberg representation does
+        not match the ev_order attribute."""
+        monkeypatch.setattr(qml.P, "ev_order", 2)
+
+        with pytest.raises(ValueError, match="Mismatch between the polynomial order"):
+            _transform_observable(qml.P(0), np.identity(3), device_wires=[0])
+
+    def test_higher_order_observable(self, monkeypatch):
+        """An exception should be raised if the observable is higher than 2nd order."""
+        monkeypatch.setattr(qml.P, "ev_order", 3)
+
+        with pytest.raises(NotImplementedError, match="order > 2 not implemented"):
+            _transform_observable(qml.P(0), np.identity(3), device_wires=[0])
+
+    def test_first_order_transform(self, tol):
+        """Test that a first order observable is transformed correctly"""
+        # create a symmetric transformation
+        Z = np.arange(3 ** 2).reshape(3, 3)
+        Z = Z.T + Z
+
+        obs = qml.X(0)
+        res = _transform_observable(obs, Z, device_wires=[0])
+
+        # The Heisenberg representation of the X
+        # operator is simply... X
+        expected = np.array([0, 1, 0]) @ Z
+
+        assert isinstance(res, qml.PolyXP)
+        assert res.wires.labels == (0,)
+        assert np.allclose(res.data[0], expected, atol=tol, rtol=0)
+
+    def test_second_order_transform(self, tol):
+        """Test that a second order observable is transformed correctly"""
+        # create a symmetric transformation
+        Z = np.arange(3 ** 2).reshape(3, 3)
+        Z = Z.T + Z
+
+        obs = qml.NumberOperator(0)
+        res = _transform_observable(obs, Z, device_wires=[0])
+
+        # The Heisenberg representation of the number operator
+        # is (X^2 + P^2) / (2*hbar) - 1/2
+        A = np.array([[-0.5, 0, 0], [0, 0.25, 0], [0, 0, 0.25]])
+        expected = A @ Z + Z @ A
+
+        assert isinstance(res, qml.PolyXP)
+        assert res.wires.labels == (0,)
+        assert np.allclose(res.data[0], expected, atol=tol, rtol=0)
+
+    def test_device_wire_expansion(self, tol):
+        """Test that the transformation works correctly
+        for the case where the transformation applies to more wires
+        than the observable."""
+
+        # create a 3-mode symmetric transformation
+        wires = qml.wires.Wires([0, "a", 2])
+        ndim = 1 + 2 * len(wires)
+
+        Z = np.arange(ndim ** 2).reshape(ndim, ndim)
+        Z = Z.T + Z
+
+        obs = qml.NumberOperator(0)
+        res = _transform_observable(obs, Z, device_wires=wires)
+
+        # The Heisenberg representation of the number operator
+        # is (X^2 + P^2) / (2*hbar) - 1/2. We use the ordering
+        # I, X0, Xa, X2, P0, Pa, P2.
+        A = np.diag([-0.5, 0.25, 0.25, 0, 0, 0, 0])
+        expected = A @ Z + Z @ A
+
+        assert isinstance(res, qml.PolyXP)
+        assert res.wires == wires
+        assert np.allclose(res.data[0], expected, atol=tol, rtol=0)
+
+
+class TestParameterShiftLogic:
+    """Test for the dispatching logic of the parameter shift method"""
+
+    def test_state_non_differentiable_error(self):
+        """Test error raised if attempting to differentiate with
+        respect to a state"""
+        with qml.tape.JacobianTape() as tape:
+            qml.state()
+
+        with pytest.raises(ValueError, match=r"return the state is not supported"):
+            qml.gradients.param_shift_cv(tape, None)
+
+    def test_force_order2(self, mocker):
+        """Test that if the force_order2 keyword argument is provided,
+        the second order parameter shift rule is forced"""
+        spy = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
+
+        dev = qml.device("default.gaussian", wires=1)
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(1.0, 0.0, wires=[0])
+            qml.Rotation(2.0, wires=[0])
+            qml.expval(qml.X(0))
+
+        tape.trainable_params = {0, 1, 2}
+
+        qml.gradients.param_shift_cv(tape, dev, force_order2=False)
+        spy.assert_not_called()
+
+        qml.gradients.param_shift_cv(tape, dev, force_order2=True)
+        spy.assert_called()
+
+    def test_no_poly_xp_support(self, mocker, monkeypatch, caplog):
+        """Test that if a device does not support PolyXP
+        and the second-order parameter-shift rule is required,
+        we fallback to finite differences."""
+        spy_second_order = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
+
+        dev = qml.device("default.gaussian", wires=1)
+
+        monkeypatch.delitem(dev._observable_map, "PolyXP")
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(1.0, wires=[0])
+            qml.expval(qml.NumberOperator(0))
+
+        tape.trainable_params = {0}
+
+        with pytest.warns(UserWarning, match="does not support the PolyXP observable"):
+            qml.gradients.param_shift_cv(tape, dev)
+
+        spy_second_order.assert_not_called()
+
+    def test_no_poly_xp_support_variance(self, mocker, monkeypatch, caplog):
+        """Test that if a device does not support PolyXP
+        and the variance parameter-shift rule is required,
+        we fallback to finite differences."""
+        spy = mocker.spy(qml.gradients.parameter_shift_cv, "var_param_shift")
+        dev = qml.device("default.gaussian", wires=1)
+
+        monkeypatch.delitem(dev._observable_map, "PolyXP")
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Rotation(1.0, wires=[0])
+            qml.var(qml.X(0))
+
+        tape.trainable_params = {0}
+
+        with pytest.warns(UserWarning, match="does not support the PolyXP observable"):
+            qml.gradients.param_shift_cv(tape, dev)
+
+        spy.assert_not_called()
+
+    def test_independent_parameters_analytic(self):
+        """Test the case where expectation values are independent of some parameters. For those
+        parameters, the gradient should be evaluated to zero without executing the device."""
+        dev = qml.device("default.gaussian", wires=2)
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(1, 0, wires=[1])
+            qml.Displacement(1, 0, wires=[0])
+            qml.expval(qml.X(0))
+
+        tape.trainable_params = {0, 2}
+        tapes, fn = qml.gradients.param_shift_cv(tape, dev)
+
+        # We should only be executing the device to differentiate 1 parameter
+        # (first order, so 2 executions)
+        assert len(tapes) == 2
+
+        res = fn(dev.batch_execute(tapes))
+        assert np.allclose(res, [0, 2])
+
+        tape.trainable_params = {0, 2}
+        tapes, fn = qml.gradients.param_shift_cv(tape, dev, force_order2=True)
+
+        # We should only be executing the device to differentiate 1 parameter
+        # (second order, so 0 executions)
+        assert len(tapes) == 0
+
+        res = fn(dev.batch_execute(tapes))
+        assert np.allclose(res, [0, 2])
+
+    def test_all_independent(self):
+        """Test the case where expectation values are independent of all parameters."""
+        dev = qml.device("default.gaussian", wires=2)
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(1, 0, wires=[1])
+            qml.Displacement(1, 0, wires=[1])
+            qml.expval(qml.X(0))
+
+        tape.trainable_params = {0, 2}
+        tapes, fn = qml.gradients.param_shift_cv(tape, dev)
+        assert len(tapes) == 0
+
+        grad = fn(dev.batch_execute(tapes))
+        assert np.allclose(grad, [0, 0])
+
+
+class TestExpectationQuantumGradients:
+    """Tests for the quantum gradients of various gates
+    with expectation value output"""
+
+    def test_rotation_gradient(self, mocker, tol):
+        """Test the gradient of the rotation gate"""
+        dev = qml.device("default.gaussian", wires=2, hbar=hbar)
+
+        alpha = 0.5643
+        theta = 0.23354
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(alpha, 0.0, wires=[0])
+            qml.Rotation(theta, wires=[0])
+            qml.expval(qml.X(0))
+
+        tape.trainable_params = {2}
+
+        spy2 = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
+
+        tapes, fn = param_shift_cv(tape, dev)
+        grad_A = fn(dev.batch_execute(tapes))
+        spy2.assert_not_called()
+
+        tapes, fn = param_shift_cv(tape, dev, force_order2=True)
+        grad_A2 = fn(dev.batch_execute(tapes))
+        spy2.assert_called()
+
+        expected = -hbar * alpha * np.sin(theta)
+        assert np.allclose(grad_A, expected, atol=tol, rtol=0)
+        assert np.allclose(grad_A2, expected, atol=tol, rtol=0)
+
+    def test_beamsplitter_gradient(self, mocker, tol):
+        """Test the gradient of the beamsplitter gate"""
+        dev = qml.device("default.gaussian", wires=2, hbar=hbar)
+
+        alpha = 0.5643
+        theta = 0.23354
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(alpha, 0.0, wires=[0])
+            qml.Beamsplitter(theta, 0.0, wires=[0, 1])
+            qml.expval(qml.X(0))
+
+        tape.trainable_params = {2}
+
+        spy2 = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
+
+        tapes, fn = param_shift_cv(tape, dev)
+        grad_A = fn(dev.batch_execute(tapes))
+        spy2.assert_not_called()
+
+        tapes, fn = param_shift_cv(tape, dev, force_order2=True)
+        grad_A2 = fn(dev.batch_execute(tapes))
+        spy2.assert_called()
+
+        expected = -hbar * alpha * np.sin(theta)
+        assert np.allclose(grad_A, expected, atol=tol, rtol=0)
+        assert np.allclose(grad_A2, expected, atol=tol, rtol=0)
+
+    def test_displacement_gradient(self, mocker, tol):
+        """Test the gradient of the displacement gate"""
+        dev = qml.device("default.gaussian", wires=2, hbar=hbar)
+
+        r = 0.5643
+        phi = 0.23354
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(r, phi, wires=[0])
+            qml.expval(qml.X(0))
+
+        tape.trainable_params = {0, 1}
+
+        spy2 = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
+
+        tapes, fn = param_shift_cv(tape, dev)
+        grad_A = fn(dev.batch_execute(tapes))
+        spy2.assert_not_called()
+
+        tapes, fn = param_shift_cv(tape, dev, force_order2=True)
+        grad_A2 = fn(dev.batch_execute(tapes))
+        spy2.assert_called()
+
+        expected = [hbar * np.cos(phi), -hbar * r * np.sin(phi)]
+        assert np.allclose(grad_A, expected, atol=tol, rtol=0)
+        assert np.allclose(grad_A2, expected, atol=tol, rtol=0)
+
+    def test_squeezed_gradient(self, mocker, tol):
+        """Test the gradient of the squeezed gate. We also
+        ensure that the gradient is correct even when an operation
+        with no Heisenberg representation is a descendent."""
+        dev = qml.device("default.gaussian", wires=2, hbar=hbar)
+
+        class Rotation(qml.operation.CVOperation):
+            """Dummy operation that does not support
+            heisenberg representation"""
+
+            num_wires = 1
+            num_params = 1
+            par_domain = "R"
+            grad_method = "A"
+
+        alpha = 0.5643
+        r = 0.23354
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(alpha, 0.0, wires=[0])
+            qml.Squeezing(r, 0.0, wires=[0])
+
+            # The following two gates have no effect
+            # on the circuit gradient and expectation value
+            qml.Beamsplitter(0.0, 0.0, wires=[0, 1])
+            Rotation(0.543, wires=[1])
+
+            qml.expval(qml.X(0))
+
+        tape.trainable_params = {2}
+
+        spy2 = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
+
+        tapes, fn = param_shift_cv(tape, dev)
+        grad_A = fn(dev.batch_execute(tapes))
+        spy2.assert_not_called()
+
+        tapes, fn = param_shift_cv(tape, dev, force_order2=True)
+        grad_A2 = fn(dev.batch_execute(tapes))
+        spy2.assert_called()
+
+        expected = -np.exp(-r) * hbar * alpha
+        assert np.allclose(grad_A, expected, atol=tol, rtol=0)
+        assert np.allclose(grad_A2, expected, atol=tol, rtol=0)
+
+    def test_squeezed_number_state_gradient(self, mocker, tol):
+        """Test the numerical gradient of the squeeze gate with
+        with number state expectation is correct"""
+        dev = qml.device("default.gaussian", wires=2, hbar=hbar)
+
+        r = 0.23354
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Squeezing(r, 0.0, wires=[0])
+            # the fock state projector is a 'non-Gaussian' observable
+            qml.expval(qml.FockStateProjector(np.array([2, 0]), wires=[0, 1]))
+
+        tape.trainable_params = {0}
+
+        spy = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
+
+        tapes, fn = param_shift_cv(tape, dev)
+        grad = fn(dev.batch_execute(tapes))
+        assert tape._par_info[0]["grad_method"] == "F"
+
+        spy.assert_not_called()
+
+        # (d/dr) |<2|S(r)>|^2 = 0.5 tanh(r)^3 (2 csch(r)^2 - 1) sech(r)
+        expected = 0.5 * np.tanh(r) ** 3 * (2 / (np.sinh(r) ** 2) - 1) / np.cosh(r)
+        assert np.allclose(grad, expected, atol=tol, rtol=0)
+
+    def test_multiple_squeezing_gradient(self, mocker, tol):
+        """Test that the gradient of a circuit with two squeeze
+        gates is correct."""
+        dev = qml.device("default.gaussian", wires=2, hbar=hbar)
+
+        r0, phi0, r1, phi1 = [0.4, -0.3, -0.7, 0.2]
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Squeezing(r0, phi0, wires=[0])
+            qml.Squeezing(r1, phi1, wires=[0])
+            qml.expval(qml.NumberOperator(0))  # second order
+
+        spy2 = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
+        tapes, fn = param_shift_cv(tape, dev, force_order2=True)
+        grad_A2 = fn(dev.batch_execute(tapes))
+        spy2.assert_called()
+
+        # check against the known analytic formula
+        expected = np.zeros([4])
+        expected[0] = np.cosh(2 * r1) * np.sinh(2 * r0) + np.cos(phi0 - phi1) * np.cosh(
+            2 * r0
+        ) * np.sinh(2 * r1)
+        expected[1] = -0.5 * np.sin(phi0 - phi1) * np.sinh(2 * r0) * np.sinh(2 * r1)
+        expected[2] = np.cos(phi0 - phi1) * np.cosh(2 * r1) * np.sinh(2 * r0) + np.cosh(
+            2 * r0
+        ) * np.sinh(2 * r1)
+        expected[3] = 0.5 * np.sin(phi0 - phi1) * np.sinh(2 * r0) * np.sinh(2 * r1)
+
+        assert np.allclose(grad_A2, expected, atol=tol, rtol=0)
+
+    def test_multiple_second_order_observables(self, mocker, tol):
+        """Test that the gradient of a circuit with multiple
+        second order observables is correct"""
+
+        dev = qml.device("default.gaussian", wires=2, hbar=hbar)
+        r = [0.4, -0.7, 0.1, 0.2]
+        p = [0.1, 0.2, 0.3, 0.4]
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Squeezing(r[0], p[0], wires=[0])
+            qml.Squeezing(r[1], p[1], wires=[0])
+            qml.Squeezing(r[2], p[2], wires=[1])
+            qml.Squeezing(r[3], p[3], wires=[1])
+            qml.expval(qml.NumberOperator(0))  # second order
+            qml.expval(qml.NumberOperator(1))  # second order
+
+        spy2 = mocker.spy(qml.gradients.parameter_shift_cv, "second_order_param_shift")
+        tapes, fn = param_shift_cv(tape, dev)
+        grad_A2 = fn(dev.batch_execute(tapes))
+        spy2.assert_called()
+
+        # check against the known analytic formula
+
+        def expected_grad(r, p):
+            return np.array(
+                [
+                    np.cosh(2 * r[1]) * np.sinh(2 * r[0])
+                    + np.cos(p[0] - p[1]) * np.cosh(2 * r[0]) * np.sinh(2 * r[1]),
+                    -0.5 * np.sin(p[0] - p[1]) * np.sinh(2 * r[0]) * np.sinh(2 * r[1]),
+                    np.cos(p[0] - p[1]) * np.cosh(2 * r[1]) * np.sinh(2 * r[0])
+                    + np.cosh(2 * r[0]) * np.sinh(2 * r[1]),
+                    0.5 * np.sin(p[0] - p[1]) * np.sinh(2 * r[0]) * np.sinh(2 * r[1]),
+                ]
+            )
+
+        expected = np.zeros([2, 8])
+        expected[0, :4] = expected_grad(r[:2], p[:2])
+        expected[1, 4:] = expected_grad(r[2:], p[2:])
+
+        assert np.allclose(grad_A2, expected, atol=tol, rtol=0)
+
+    cv_ops = [getattr(qml, name) for name in qml.ops._cv__ops__]
+    analytic_cv_ops = [cls for cls in cv_ops if cls.supports_parameter_shift]
+
+    @pytest.mark.parametrize("obs", [qml.X, qml.P, qml.NumberOperator, qml.Identity])
+    @pytest.mark.parametrize("op", analytic_cv_ops)
+    def test_gradients_gaussian_circuit(self, op, obs, mocker, tol):
+        """Tests that the gradients of circuits of gaussian gates match between the
+        finite difference and analytic methods."""
+        tol = 1e-2
+
+        args = np.linspace(0.2, 0.5, op.num_params)
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(0.5, 0, wires=0)
+            op(*args, wires=range(op.num_wires))
+            qml.Beamsplitter(1.3, -2.3, wires=[0, 1])
+            qml.Displacement(-0.5, 0.1, wires=0)
+            qml.Squeezing(0.5, -1.5, wires=0)
+            qml.Rotation(-1.1, wires=0)
+            qml.expval(obs(wires=0))
+
+        dev = qml.device("default.gaussian", wires=2)
+        res = tape.execute(dev)
+
+        tape.trainable_params = set(range(2, 2 + op.num_params))
+
+        tapes, fn = qml.gradients.finite_diff(tape)
+        grad_F = fn(dev.batch_execute(tapes))
+
+        tapes, fn = param_shift_cv(tape, dev, force_order2=True)
+        grad_A2 = fn(dev.batch_execute(tapes))
+
+        # check that every parameter is analytic
+        for i in range(op.num_params):
+            assert tape._par_info[2 + i]["grad_method"][0] == "A"
+
+        assert np.allclose(grad_A2, grad_F, atol=tol, rtol=0)
+
+        if obs.ev_order == 1:
+            tapes, fn = param_shift_cv(tape, dev)
+            grad_A = fn(dev.batch_execute(tapes))
+            assert np.allclose(grad_A, grad_F, atol=tol, rtol=0)
+
+    @pytest.mark.parametrize("t", [0, 1])
+    def test_interferometer(self, t, tol):
+        """An integration test for CV gates that support analytic differentiation
+        if succeeding the gate to be differentiated, but cannot be differentiated
+        themselves (for example, they may be Gaussian but accept no parameters,
+        or may accept a numerical array parameter.).
+
+        This ensures that, assuming their _heisenberg_rep is defined, the quantum
+        gradient analytic method can still be used, and returns the correct result.
+
+        Currently, the only such operation is qml.Interferometer. In the future,
+        we may consider adding a qml.GaussianTransfom operator.
+        """
+
+        if t == 1:
+            pytest.xfail(
+                "There is a bug in the second order CV parameter-shift rule; "
+                "phase arguments return the incorrect derivative."
+            )
+
+            # Note: this bug currently affects PL core as well:
+            #
+            # dev = qml.device("default.gaussian", wires=2)
+            #
+            # U = np.array([[ 0.51310276+0.81702166j,  0.13649626+0.22487759j],
+            #         [ 0.26300233+0.00556194j, -0.96414101-0.03508489j]])
+            #
+            # @qml.qnode(dev)
+            # def circuit(r, phi):
+            #     qml.Displacement(r, phi, wires=0)
+            #     qml.Interferometer(U, wires=[0, 1])
+            #     return qml.expval(qml.X(0))
+            #
+            # r = 0.543
+            # phi = 0.
+            #
+            # >>> print(circuit.jacobian([r, phi], options={"force_order2":False}))
+            # [[ 1.02620552 0.14823494]]
+            # >>> print(circuit.jacobian([r, phi], options={"force_order2":True}))
+            # [[ 1.02620552 -0.88728552]]
+
+        U = np.array(
+            [
+                [0.51310276 + 0.81702166j, 0.13649626 + 0.22487759j],
+                [0.26300233 + 0.00556194j, -0.96414101 - 0.03508489j],
+            ]
+        )
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(0.543, 0, wires=0)
+            qml.Interferometer(U, wires=[0, 1])
+            qml.expval(qml.X(0))
+
+        tape.trainable_params = {t}
+
+        dev = qml.device("default.gaussian", wires=2)
+
+        tapes, fn = qml.gradients.finite_diff(tape)
+        grad_F = fn(dev.batch_execute(tapes))
+
+        tapes, fn = param_shift_cv(tape, dev)
+        grad_A = fn(dev.batch_execute(tapes))
+
+        tapes, fn = param_shift_cv(tape, dev, force_order2=True)
+        grad_A2 = fn(dev.batch_execute(tapes))
+
+        assert tape._par_info[0]["grad_method"] == "A"
+        assert tape._par_info[1]["grad_method"] == "A"
+
+        # the different methods agree
+        assert np.allclose(grad_A, grad_F, atol=tol, rtol=0)
+        assert np.allclose(grad_A2, grad_F, atol=tol, rtol=0)
+
+
+class TestVarianceQuantumGradients:
+    """Tests for the quantum gradients of various gates
+    with variance measurements"""
+
+    def test_first_order_observable(self, tol):
+        """Test variance of a first order CV observable"""
+        dev = qml.device("default.gaussian", wires=1)
+
+        r = 0.543
+        phi = -0.654
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Squeezing(r, 0, wires=0)
+            qml.Rotation(phi, wires=0)
+            qml.var(qml.X(0))
+
+        tape.trainable_params = {0, 2}
+
+        res = tape.execute(dev)
+        expected = np.exp(2 * r) * np.sin(phi) ** 2 + np.exp(-2 * r) * np.cos(phi) ** 2
+        assert np.allclose(res, expected, atol=tol, rtol=0)
+
+        # circuit jacobians
+        tapes, fn = qml.gradients.finite_diff(tape)
+        grad_F = fn(dev.batch_execute(tapes))
+
+        tapes, fn = param_shift_cv(tape, dev)
+        grad_A = fn(dev.batch_execute(tapes))
+
+        expected = np.array(
+            [
+                [
+                    2 * np.exp(2 * r) * np.sin(phi) ** 2 - 2 * np.exp(-2 * r) * np.cos(phi) ** 2,
+                    2 * np.sinh(2 * r) * np.sin(2 * phi),
+                ]
+            ]
+        )
+        assert np.allclose(grad_A, expected, atol=tol, rtol=0)
+        assert np.allclose(grad_F, expected, atol=tol, rtol=0)
+
+    def test_second_order_cv(self, tol):
+        """Test variance of a second order CV expectation value"""
+        dev = qml.device("default.gaussian", wires=1)
+
+        n = 0.12
+        a = 0.765
+
+        with qml.tape.JacobianTape() as tape:
+            qml.ThermalState(n, wires=0)
+            qml.Displacement(a, 0, wires=0)
+            qml.var(qml.NumberOperator(0))
+
+        tape.trainable_params = {0, 1}
+
+        res = tape.execute(dev)
+        expected = n ** 2 + n + np.abs(a) ** 2 * (1 + 2 * n)
+        assert np.allclose(res, expected, atol=tol, rtol=0)
+
+        # circuit jacobians
+        tapes, fn = qml.gradients.finite_diff(tape)
+        grad_F = fn(dev.batch_execute(tapes))
+
+        expected = np.array([[2 * a ** 2 + 2 * n + 1, 2 * a * (2 * n + 1)]])
+        assert np.allclose(grad_F, expected, atol=tol, rtol=0)
+
+    def test_expval_and_variance(self, tol):
+        """Test that the gradient works for a combination of CV expectation
+        values and variances"""
+        dev = qml.device("default.gaussian", wires=3)
+
+        a, b = [0.54, -0.423]
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(0.5, 0, wires=0)
+            qml.Squeezing(a, 0, wires=0)
+            qml.Squeezing(b, 0, wires=1)
+            qml.Beamsplitter(0.6, -0.3, wires=[0, 1])
+            qml.Squeezing(-0.3, 0, wires=2)
+            qml.Beamsplitter(1.4, 0.5, wires=[1, 2])
+            qml.var(qml.X(0))
+            qml.expval(qml.X(1))
+            qml.var(qml.X(2))
+
+        tape.trainable_params = {2, 4}
+
+        # jacobians must match
+        tapes, fn = qml.gradients.finite_diff(tape)
+        grad_F = fn(dev.batch_execute(tapes))
+
+        tapes, fn = param_shift_cv(tape, dev)
+        grad_A = fn(dev.batch_execute(tapes))
+
+        assert np.allclose(grad_A, grad_F, atol=tol, rtol=0)
+
+    def test_error_analytic_second_order(self):
+        """Test exception raised if attempting to use a second
+        order observable to compute the variance derivative analytically"""
+        dev = qml.device("default.gaussian", wires=1)
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(1.0, 0, wires=0)
+            qml.var(qml.NumberOperator(0))
+
+        tape.trainable_params = {0}
+
+        with pytest.raises(ValueError, match=r"cannot be used with the argument\(s\) \{0\}"):
+            param_shift_cv(tape, dev, fallback_fn=None)
+
+    def test_error_unsupported_grad_recipe(self, monkeypatch):
+        """Test exception raised if attempting to use the second order rule for
+        computing the gradient analytically of an expectation value that
+        contains an operation with more than two terms in the gradient recipe"""
+
+        class DummyOp(qml.operation.CVOperation):
+            num_wires = 1
+            num_params = 1
+            par_domain = "R"
+            grad_method = "A"
+            grad_recipe = ([[1, 1, 1], [1, 1, 1], [1, 1, 1]],)
+
+        dev = qml.device("default.gaussian", wires=1)
+
+        dev.operations.add(DummyOp)
+
+        with qml.tape.JacobianTape() as tape:
+            DummyOp(1, wires=[0])
+            qml.expval(qml.X(0))
+
+        with monkeypatch.context() as m:
+            m.setattr(tape, "_grad_method_validation", lambda *args: ("A",))
+            tape._par_info[0]["grad_method"] = "A"
+            tape.trainable_params = {0}
+
+            with pytest.raises(
+                NotImplementedError, match=r"analytic gradient for order-2 operators is unsupported"
+            ):
+                param_shift_cv(tape, dev, force_order2=True)
+
+    cv_ops = [getattr(qml, name) for name in qml.ops._cv__ops__]
+    analytic_cv_ops = [cls for cls in cv_ops if cls.supports_parameter_shift]
+
+    @pytest.mark.parametrize("obs", [qml.X, qml.P, qml.NumberOperator, qml.Identity])
+    @pytest.mark.parametrize("op", analytic_cv_ops)
+    def test_gradients_gaussian_circuit(self, op, obs, tol):
+        """Tests that the gradients of circuits of gaussian gates match between the
+        finite difference and analytic methods."""
+        tol = 1e-2
+
+        args = np.linspace(0.2, 0.5, op.num_params)
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Displacement(0.5, 0, wires=0)
+            op(*args, wires=range(op.num_wires))
+            qml.Beamsplitter(1.3, -2.3, wires=[0, 1])
+            qml.Displacement(-0.5, 0.1, wires=0)
+            qml.Squeezing(0.5, -1.5, wires=0)
+            qml.Rotation(-1.1, wires=0)
+            qml.var(obs(wires=0))
+
+        dev = qml.device("default.gaussian", wires=2)
+        res = tape.execute(dev)
+
+        tape.trainable_params = set(range(2, 2 + op.num_params))
+
+        # jacobians must match
+        tapes, fn = qml.gradients.finite_diff(tape)
+        grad_F = fn(dev.batch_execute(tapes))
+
+        tapes, fn = param_shift_cv(tape, dev)
+        grad_A = fn(dev.batch_execute(tapes))
+
+        tapes, fn = param_shift_cv(tape, dev, force_order2=True)
+        grad_A2 = fn(dev.batch_execute(tapes))
+
+        assert np.allclose(grad_A2, grad_F, atol=tol, rtol=0)
+        assert np.allclose(grad_A, grad_F, atol=tol, rtol=0)
+
+        # check that every parameter is analytic
+        if obs != qml.NumberOperator:
+            for i in range(op.num_params):
+                assert tape._par_info[2 + i]["grad_method"][0] == "A"
+
+    def test_squeezed_mean_photon_variance(self, tol):
+        """Test gradient of the photon variance of a displaced thermal state"""
+        dev = qml.device("default.gaussian", wires=1)
+
+        r = 0.12
+        phi = 0.105
+
+        with qml.tape.JacobianTape() as tape:
+            qml.Squeezing(r, 0, wires=0)
+            qml.Rotation(phi, wires=0)
+            qml.var(qml.X(wires=[0]))
+
+        tape.trainable_params = {0, 2}
+        tapes, fn = param_shift_cv(tape, dev)
+        grad = fn(dev.batch_execute(tapes))
+        expected = np.array(
+            [
+                2 * np.exp(2 * r) * np.sin(phi) ** 2 - 2 * np.exp(-2 * r) * np.cos(phi) ** 2,
+                2 * np.sinh(2 * r) * np.sin(2 * phi),
+            ]
+        )
+        assert np.allclose(grad, expected, atol=tol, rtol=0)
+
+    def test_displaced_thermal_mean_photon_variance(self, tol):
+        """Test gradient of the photon variance of a displaced thermal state"""
+        dev = qml.device("default.gaussian", wires=1)
+
+        n = 0.12
+        a = 0.105
+
+        with qml.tape.JacobianTape() as tape:
+            qml.ThermalState(n, wires=0)
+            qml.Displacement(a, 0, wires=0)
+            qml.var(qml.TensorN(wires=[0]))
+
+        tape.trainable_params = {0, 1}
+        tapes, fn = param_shift_cv(tape, dev)
+        grad = fn(dev.batch_execute(tapes))
+        expected = np.array([2 * a ** 2 + 2 * n + 1, 2 * a * (2 * n + 1)])
+        assert np.allclose(grad, expected, atol=tol, rtol=0)
+
+
+class TestParamShiftInterfaces:
+    """Test that the transform is differentiable"""
+
+    def test_autograd_gradient(self, tol):
+        """Tests that the output of the parameter-shift CV transform
+        can be differentiated using autograd, yielding second derivatives."""
+        dev = qml.device("default.gaussian", wires=1)
+        from pennylane.interfaces.autograd import AutogradInterface
+
+        r = 0.12
+        phi = 0.105
+
+        def cost_fn(x):
+            with AutogradInterface.apply(qml.tape.CVParamShiftTape()) as tape:
+                qml.Squeezing(x[0], 0, wires=0)
+                qml.Rotation(x[1], wires=0)
+                qml.var(qml.X(wires=[0]))
+
+            tapes, fn = param_shift_cv(tape, dev)
+            return fn([t.execute(dev) for t in tapes])[0, 1]
+
+        params = np.array([r, phi], requires_grad=True)
+        grad = qml.jacobian(cost_fn)(params)
+        expected = np.array(
+            [4 * np.cosh(2 * r) * np.sin(2 * phi), 4 * np.cos(2 * phi) * np.sinh(2 * r)]
+        )
+        assert np.allclose(grad, expected, atol=tol, rtol=0)
+
+    def test_tf(self, tol):
+        """Tests that the output of the parameter-shift CV transform
+        can be executed using TF"""
+        tf = pytest.importorskip("tensorflow")
+        from pennylane.interfaces.tf import TFInterface
+
+        dev = qml.device("default.gaussian", wires=1)
+        params = tf.Variable([0.543, -0.654], dtype=tf.float64)
+
+        with tf.GradientTape() as t:
+            with TFInterface.apply(qml.tape.CVParamShiftTape()) as tape:
+                qml.Squeezing(params[0], 0, wires=0)
+                qml.Rotation(params[1], wires=0)
+                qml.var(qml.X(wires=[0]))
+
+            tapes, fn = qml.gradients.param_shift_cv(tape, dev)
+            jac = fn([tp.execute(dev) for tp in tapes])
+            res = jac[0, 1]
+
+        r, phi = 1.0 * params
+
+        expected = np.array(
+            [
+                2 * np.exp(2 * r) * np.sin(phi) ** 2 - 2 * np.exp(-2 * r) * np.cos(phi) ** 2,
+                2 * np.sinh(2 * r) * np.sin(2 * phi),
+            ]
+        )
+        assert np.allclose(jac, expected, atol=tol, rtol=0)
+
+        grad = t.jacobian(res, params)
+        expected = np.array(
+            [4 * np.cosh(2 * r) * np.sin(2 * phi), 4 * np.cos(2 * phi) * np.sinh(2 * r)]
+        )
+        assert np.allclose(grad, expected, atol=tol, rtol=0)
+
+    def test_torch(self, tol):
+        """Tests that the output of the parameter-shift CV transform
+        can be executed using Torch."""
+        torch = pytest.importorskip("torch")
+        from pennylane.interfaces.torch import TorchInterface
+
+        dev = qml.device("default.gaussian", wires=1)
+        params = torch.tensor([0.543, -0.654], dtype=torch.float64, requires_grad=True)
+
+        with TorchInterface.apply(qml.tape.CVParamShiftTape()) as tape:
+            qml.Squeezing(params[0], 0, wires=0)
+            qml.Rotation(params[1], wires=0)
+            qml.var(qml.X(wires=[0]))
+
+        tapes, fn = qml.gradients.param_shift_cv(tape, dev)
+        jac = fn([t.execute(dev) for t in tapes])
+
+        r, phi = params.detach().numpy()
+
+        expected = np.array(
+            [
+                2 * np.exp(2 * r) * np.sin(phi) ** 2 - 2 * np.exp(-2 * r) * np.cos(phi) ** 2,
+                2 * np.sinh(2 * r) * np.sin(2 * phi),
+            ]
+        )
+        assert np.allclose(jac.detach().numpy(), expected, atol=tol, rtol=0)
+
+        cost = jac[0, 1]
+        cost.backward()
+        hess = params.grad
+        expected = np.array(
+            [4 * np.cosh(2 * r) * np.sin(2 * phi), 4 * np.cos(2 * phi) * np.sinh(2 * r)]
+        )
+        assert np.allclose(hess.detach().numpy(), expected, atol=0.1, rtol=0)
+
+    def test_jax(self, tol):
+        """Tests that the output of the parameter-shift CV transform
+        can be differentiated using JAX, yielding second derivatives."""
+        jax = pytest.importorskip("jax")
+        from jax import numpy as jnp
+        from pennylane.interfaces.jax import JAXInterface
+        from jax.config import config
+
+        config.update("jax_enable_x64", True)
+
+        dev = qml.device("default.gaussian", wires=2)
+        params = jnp.array([0.543, -0.654])
+
+        def cost_fn(x):
+            with JAXInterface.apply(qml.tape.CVParamShiftTape()) as tape:
+                qml.Squeezing(params[0], 0, wires=0)
+                qml.Rotation(params[1], wires=0)
+                qml.var(qml.X(wires=[0]))
+
+            tape.trainable_params = {0, 2}
+            tapes, fn = qml.gradients.param_shift_cv(tape, dev)
+            jac = fn([t.execute(dev) for t in tapes])
+            return jac
+
+        r, phi = params
+        res = cost_fn(params)
+        expected = np.array(
+            [
+                2 * np.exp(2 * r) * np.sin(phi) ** 2 - 2 * np.exp(-2 * r) * np.cos(phi) ** 2,
+                2 * np.sinh(2 * r) * np.sin(2 * phi),
+            ]
+        )
+        assert np.allclose(res, expected, atol=tol, rtol=0)
+
+        pytest.xfail("The CV Operation methods have not been updated to support autodiff")
+
+        res = jax.jacobian(cost_fn)(params)
+        expected = np.array(
+            [
+                [
+                    4 * np.exp(-2 * r) * (np.cos(phi) ** 2 + np.exp(4 * r) * np.sin(phi) ** 2),
+                    4 * np.cosh(2 * r) * np.sin(2 * phi),
+                ],
+                [4 * np.cosh(2 * r) * np.sin(2 * phi), 4 * np.cos(2 * phi) * np.sinh(2 * r)],
+            ]
+        )
+        assert np.allclose(res, expected, atol=tol, rtol=0)
diff --git a/tests/tape/test_jacobian_tape.py b/tests/tape/test_jacobian_tape.py
index 528d25039..983b96de3 100644
--- a/tests/tape/test_jacobian_tape.py
+++ b/tests/tape/test_jacobian_tape.py
@@ -456,19 +456,6 @@ class TestJacobian:
 
         assert len(res) == num_params
 
-    @pytest.mark.parametrize("argnum", [1, 2, 3, -1])
-    def test_choose_params_and_methods_raises(self, argnum):
-        """Test that the _choose_params_and_methods helper method raises an
-        error if incorrect trainable parameters are specified."""
-        tape = JacobianTape()
-        tape.trainable_params = [0]
-        diff_methods = ["F"]
-        with pytest.raises(
-            ValueError,
-            match="Incorrect trainable parameters",
-        ):
-            res = tape._choose_params_with_methods(diff_methods, argnum)
-
     def test_choose_params_and_methods_warns_no_params(self):
         """Test that the _choose_params_and_methods helper method warns if an
         empty list was passed as argnum."""