Repository: PennyLaneAI/pennylane
Pull request #5097 (2024-01-23 23:49:48)

Remove `qml.ExpValCost`
Remove` qml.ExpValCost` after completing deprecation cycle started here: https://github.com/PennyLaneAI/pennylane/pull/2571



codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5097?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Comparison is base [(`4e2e344`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/4e2e344f33e9d289f0f1ec6be8c32627dce515cb?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.68% compared to head [(`ecce86b`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5097?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.63%.
> Report is 1 commits behind head on master.


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #5097      +/-   ##
==========================================
- Coverage   99.68%   99.63%   -0.06%     
==========================================
  Files         392      390       -2     
  Lines       35844    35519     -325     
==========================================
- Hits        35732    35389     -343     
- Misses        112      130      +18     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5097?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).


Patch:
diff --git a/doc/code/qml_qaoa.rst b/doc/code/qml_qaoa.rst
index 20daf7b35..01f6e01a9 100644
--- a/doc/code/qml_qaoa.rst
+++ b/doc/code/qml_qaoa.rst
@@ -87,22 +87,25 @@ computational basis states, and then repeatedly apply QAOA layers with the
 .. code-block:: python3
 
     # Repeatedly applies layers of the QAOA ansatz
-    def circuit(params, **kwargs):
+    def circuit(params):
 
         for w in wires:
             qml.Hadamard(wires=w)
 
         qml.layer(qaoa_layer, 2, params[0], params[1])
 
-With the circuit defined, we call the device on which QAOA will be executed, as well as the ``qml.ExpvalCost``, which
-creates the QAOA cost function: the expected value of the cost Hamiltonian with respect to the parametrized output
+With the circuit defined, we call the device on which QAOA will be executed and use ``qml.expval()`` to
+create the QAOA cost function: the expected value of the cost Hamiltonian with respect to the parametrized output
 of the QAOA circuit.
 
 .. code-block:: python3
 
     # Defines the device and the QAOA cost function
     dev = qml.device('default.qubit', wires=len(wires))
-    cost_function = qml.ExpvalCost(circuit, cost_h, dev)
+    @qml.qnode(dev)
+    def cost_function(params):
+        circuit(params)
+        return qml.expval(cost_h)
 
 >>> print(cost_function([[1, 1], [1, 1]]))
 -1.8260274380964299
diff --git a/doc/development/deprecations.rst b/doc/development/deprecations.rst
index 26f4fb514..0f7b15ff1 100644
--- a/doc/development/deprecations.rst
+++ b/doc/development/deprecations.rst
@@ -62,21 +62,6 @@ Pending deprecations
   - Deprecated in v0.34
   - Will be removed in v0.36
 
-* ``qml.ExpvalCost`` has been deprecated, and usage will now raise a warning.
-
-  - Deprecated in v0.24
-  - Will be removed in v0.35
-
-  Instead, it is recommended to simply
-  pass Hamiltonians to the ``qml.expval`` function inside QNodes:
-
-  .. code-block:: python
-
-    @qml.qnode(dev)
-    def ansatz(params):
-        some_qfunc(params)
-        return qml.expval(Hamiltonian)
-
 * ``PauliWord`` and ``PauliSentence`` no longer use ``*`` for matrix and tensor products,
   but instead use ``@`` to conform with the PennyLane convention.
 
@@ -140,8 +125,8 @@ Completed deprecation cycles
   - Deprecated in v0.34
   - Removed in v0.35
 
-* ``QuantumScript.is_sampled`` and ``QuantumScript.all_sampled`` have been removed. Users should now validate
-  these properties manually.
+* ``QuantumScript.is_sampled`` and ``QuantumScript.all_sampled`` have been removed.
+  Users should now validate these properties manually.
 
   .. code-block:: python
 
@@ -154,6 +139,18 @@ Completed deprecation cycles
   - Deprecated in v0.34
   - Removed in v0.35
 
+* ``qml.ExpvalCost`` has been removed. Users should use ``qml.expval()`` instead.
+
+  .. code-block:: python
+
+    @qml.qnode(dev)
+    def cost_function(params):
+        some_qfunc(params)
+        return qml.expval(Hamiltonian)
+
+  - Deprecated in v0.24
+  - Removed in v0.35
+
 * Specifying ``control_values`` passed to ``qml.ctrl`` as a string is no longer supported.
 
   - Deprecated in v0.25
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index d37156fde..1ace79487 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -173,6 +173,9 @@
   you should use `qml.ops.one_qubit_decomposition` and `qml.ops.two_qubit_decomposition`.
   [(#5091)](https://github.com/PennyLaneAI/pennylane/pull/5091)
 
+* `qml.ExpvalCost` has been removed. Users should use `qml.expval()` moving forward.
+  [(#5097)](https://github.com/PennyLaneAI/pennylane/pull/5097)
+
 <h3>Deprecations 👋</h3>
 
 * `Operator.validate_subspace(subspace)` has been relocated to the `qml.ops.qutrit.parametric_ops`
diff --git a/pennylane/__init__.py b/pennylane/__init__.py
index a000a50e4..05e64279b 100644
--- a/pennylane/__init__.py
+++ b/pennylane/__init__.py
@@ -116,7 +116,6 @@ from pennylane.ops.functions import (
     comm,
 )
 from pennylane.optimize import *
-from pennylane.vqe import ExpvalCost
 from pennylane.debugging import snapshots
 from pennylane.shadows import ClassicalShadow
 from pennylane.qcut import cut_circuit, cut_circuit_mc
diff --git a/pennylane/optimize/qng.py b/pennylane/optimize/qng.py
index 8ba0f02e4..98c75f8d0 100644
--- a/pennylane/optimize/qng.py
+++ b/pennylane/optimize/qng.py
@@ -178,11 +178,11 @@ class QNGOptimizer(GradientDescentOptimizer):
             prior to the step
         """
         # pylint: disable=arguments-differ
-        if not isinstance(qnode, (qml.QNode, qml.ExpvalCost)) and metric_tensor_fn is None:
+        if not isinstance(qnode, qml.QNode) and metric_tensor_fn is None:
             raise ValueError(
-                "The objective function must either be encoded as a single QNode or "
-                "an ExpvalCost object for the natural gradient to be automatically computed. "
-                "Otherwise, metric_tensor_fn must be explicitly provided to the optimizer."
+                "The objective function must be encoded as a single QNode for the natural gradient "
+                "to be automatically computed. Otherwise, metric_tensor_fn must be explicitly "
+                "provided to the optimizer."
             )
 
         if recompute_tensor or self.metric_tensor is None:
diff --git a/pennylane/optimize/shot_adaptive.py b/pennylane/optimize/shot_adaptive.py
index 78f649f52..948f5ce59 100644
--- a/pennylane/optimize/shot_adaptive.py
+++ b/pennylane/optimize/shot_adaptive.py
@@ -23,8 +23,8 @@ from .gradient_descent import GradientDescentOptimizer
 
 
 class ShotAdaptiveOptimizer(GradientDescentOptimizer):
-    r"""Optimizer where the shot rate is adaptively calculated using the variances of the parameter-shift
-    gradient.
+    r"""Optimizer where the shot rate is adaptively calculated using the variances of the
+    parameter-shift gradient.
 
     By keeping a running average of the parameter-shift gradient and the *variance*
     of the parameter-shift gradient, this optimizer frugally distributes a shot
@@ -33,33 +33,17 @@ class ShotAdaptiveOptimizer(GradientDescentOptimizer):
     In addition, weighted random sampling can be used to further distribute the
     shot budget across the local terms from which the Hamiltonian is constructed.
 
-    .. warning::
-
-        ``ExpvalCost`` is deprecated. Instead, it is recommended to simply
-        pass Hamiltonians to the :func:`~pennylane.expval` function inside QNodes.
-
-        .. code-block:: python
-
-            @qml.qnode(dev)
-            def ansatz(params):
-                some_qfunc(params)
-                return qml.expval(Hamiltonian)
-
-        You should also set ``term_sampling='weighted_random_sampling'`` to continue using weighted
-        random sampling in ``ShotAdaptiveOptimizer``.
-
     .. note::
 
-        The shot adaptive optimizer only supports single QNodes or :class:`~.ExpvalCost` objects as
-        objective functions. The bound device must also be instantiated with a finite number
-        of shots.
+        The shot adaptive optimizer only supports single QNode objects as objective functions.
+        The bound device must also be instantiated with a finite number of shots.
 
     Args:
         min_shots (int): The minimum number of shots used to estimate the expectations
             of each term in the Hamiltonian. Note that this must be larger than 2 for the variance
             of the gradients to be computed.
-        mu (float): The running average constant :math:`\mu \in [0, 1]`. Used to control how quickly the
-            number of shots recommended for each gradient component changes.
+        mu (float): The running average constant :math:`\mu \in [0, 1]`. Used to control how
+            quickly the number of shots recommended for each gradient component changes.
         b (float): Regularization bias. The bias should be kept small, but non-zero.
         term_sampling (str): The random sampling algorithm to multinomially distribute the shot
             budget across terms in the Hamiltonian expectation value. Currently, only
@@ -219,66 +203,6 @@ class ShotAdaptiveOptimizer(GradientDescentOptimizer):
 
         super().__init__(stepsize=stepsize)
 
-    @staticmethod
-    def weighted_random_sampling(qnodes, coeffs, shots, argnums, *args, **kwargs):
-        """Returns an array of length ``shots`` containing single-shot estimates
-        of the Hamiltonian gradient. The shots are distributed randomly over
-        the terms in the Hamiltonian, as per a multinomial distribution.
-
-        Args:
-            qnodes (Sequence[.QNode]): Sequence of QNodes, each one when evaluated
-                returning the corresponding expectation value of a term in the Hamiltonian.
-            coeffs (Sequence[float]): Sequences of coefficients corresponding to
-                each term in the Hamiltonian. Must be the same length as ``qnodes``.
-            shots (int): The number of shots used to estimate the Hamiltonian expectation
-                value. These shots are distributed over the terms in the Hamiltonian,
-                as per a Multinomial distribution.
-            argnums (Sequence[int]): the QNode argument indices which are trainable
-            *args: Arguments to the QNodes
-            **kwargs: Keyword arguments to the QNodes
-
-        Returns:
-            array[float]: the single-shot gradients of the Hamiltonian expectation value
-        """
-
-        # determine the shot probability per term
-        prob_shots = np.abs(coeffs) / np.sum(np.abs(coeffs))
-
-        # construct the multinomial distribution, and sample
-        # from it to determine how many shots to apply per term
-        si = multinomial(n=shots, p=prob_shots)
-        shots_per_term = si.rvs()[0]
-
-        grads = []
-
-        for h, c, p, s in zip(qnodes, coeffs, prob_shots, shots_per_term):
-            # if the number of shots is 0, do nothing
-            if s == 0:
-                continue
-
-            # set the QNode device shots
-            new_shots = 1 if s == 1 else [(1, int(s))]
-
-            if s > 1:
-
-                def cost(*args, **kwargs):
-                    # pylint: disable=cell-var-from-loop
-                    return qml.math.stack(h(*args, **kwargs))
-
-            else:
-                cost = h
-
-            jacs = qml.jacobian(cost, argnum=argnums)(*args, **kwargs, shots=new_shots)
-
-            if s == 1:
-                jacs = [np.expand_dims(j, 0) for j in jacs]
-
-            # Divide each term by the probability per shot. This is
-            # because we are sampling one at a time.
-            grads.append([(c * j / p) for j in jacs])
-
-        return [np.concatenate(i) for i in zip(*grads)]
-
     @staticmethod
     def qnode_weighted_random_sampling(qnode, coeffs, observables, shots, argnums, *args, **kwargs):
         """Returns an array of length ``shots`` containing single-shot estimates
@@ -379,38 +303,6 @@ class ShotAdaptiveOptimizer(GradientDescentOptimizer):
         if self.stepsize > 2 / self.lipschitz:
             raise ValueError(f"The learning rate must be less than {2 / self.lipschitz}")
 
-    def _single_shot_expval_gradients(self, expval_cost, args, kwargs):
-        """Compute the single shot gradients of an ExpvalCost object"""
-        qnodes = expval_cost.qnodes
-        coeffs = expval_cost.hamiltonian.coeffs
-
-        self.check_device(qnodes[0].device)
-
-        if self.lipschitz is None:
-            self.check_learning_rate(coeffs)
-
-        if self.term_sampling == "weighted_random_sampling":
-            grads = self.weighted_random_sampling(
-                qnodes, coeffs, self.max_shots, self.trainable_args, *args, **kwargs
-            )
-        elif self.term_sampling is None:
-            new_shots = [(1, int(self.max_shots))]
-            # We iterate over each trainable argument, rather than using
-            # qml.jacobian(expval_cost), to take into account the edge case where
-            # different arguments have different shapes and cannot be stacked.
-            grads = [
-                qml.jacobian(expval_cost, argnum=i)(*args, **kwargs, shots=new_shots)
-                for i in self.trainable_args
-            ]
-        else:
-            raise ValueError(
-                f"Unknown Hamiltonian term sampling method {self.term_sampling}. "
-                "Only term_sampling='weighted_random_sampling' and "
-                "term_sampling=None currently supported."
-            )
-
-        return grads
-
     def _single_shot_qnode_gradients(self, qnode, args, kwargs):
         """Compute the single shot gradients of a QNode."""
         self.check_device(qnode.device)
@@ -460,14 +352,12 @@ class ShotAdaptiveOptimizer(GradientDescentOptimizer):
             tuple[array[float], array[float]]: a tuple of NumPy arrays containing the gradient
             :math:`\nabla f(x^{(t)})` and the variance of the gradient
         """
-        if isinstance(objective_fn, qml.ExpvalCost):
-            grads = self._single_shot_expval_gradients(objective_fn, args, kwargs)
-        elif isinstance(objective_fn, qml.QNode) or hasattr(objective_fn, "device"):
+        if isinstance(objective_fn, qml.QNode) or hasattr(objective_fn, "device"):
             grads = self._single_shot_qnode_gradients(objective_fn, args, kwargs)
         else:
             raise ValueError(
-                "The objective function must either be encoded as a single QNode or "
-                "an ExpvalCost object for the Shot adaptive optimizer. "
+                "The objective function must be encoded as a single QNode object for the shot "
+                "adaptive optimizer. "
             )
 
         # grads will have dimension [max(self.s), *params.shape]
diff --git a/pennylane/templates/layers/particle_conserving_u1.py b/pennylane/templates/layers/particle_conserving_u1.py
index fc0787ec4..3b6430e3d 100644
--- a/pennylane/templates/layers/particle_conserving_u1.py
+++ b/pennylane/templates/layers/particle_conserving_u1.py
@@ -216,10 +216,13 @@ class ParticleConservingU1(Operation):
             dev = qml.device('default.qubit', wires=qubits)
 
             # Define the ansatz
-            ansatz = partial(qml.ParticleConservingU1, init_state=ref_state)
+            ansatz = partial(qml.ParticleConservingU1, init_state=ref_state, wires=dev.wires)
 
             # Define the cost function
-            cost_fn = qml.ExpvalCost(ansatz, h, dev)
+            @qml.qnode(dev)
+            def cost_fn(params):
+                ansatz(params)
+                return qml.expval(h)
 
             # Compute the expectation value of 'h'
             layers = 2
diff --git a/pennylane/templates/layers/particle_conserving_u2.py b/pennylane/templates/layers/particle_conserving_u2.py
index 8518b6dcd..70fa3f12c 100644
--- a/pennylane/templates/layers/particle_conserving_u2.py
+++ b/pennylane/templates/layers/particle_conserving_u2.py
@@ -100,7 +100,6 @@ class ParticleConservingU2(Operation):
     .. details::
         :title: Usage Details
 
-
         #. The number of wires has to be equal to the number of spin orbitals included in
            the active space.
 
@@ -126,10 +125,13 @@ class ParticleConservingU2(Operation):
             dev = qml.device('default.qubit', wires=qubits)
 
             # Define the ansatz
-            ansatz = partial(qml.ParticleConservingU2, init_state=ref_state)
+            ansatz = partial(qml.ParticleConservingU2, init_state=ref_state, wires=dev.wires)
 
             # Define the cost function
-            cost_fn = qml.ExpvalCost(ansatz, h, dev)
+            @qml.qnode(dev)
+            def cost_fn(params):
+                ansatz(params)
+                return qml.expval(h)
 
             # Compute the expectation value of 'h' for a given set of parameters
             layers = 1
diff --git a/pennylane/transforms/batch_transform.py b/pennylane/transforms/batch_transform.py
index 60bdc55fb..97a1d9f79 100644
--- a/pennylane/transforms/batch_transform.py
+++ b/pennylane/transforms/batch_transform.py
@@ -360,7 +360,7 @@ class batch_transform:
             # tapes, fn = some_transform(tape, *transform_args)
             return self._tape_wrapper(*targs, **tkwargs)(qnode)
 
-        if isinstance(qnode, (qml.QNode, qml.ExpvalCost)):
+        if isinstance(qnode, qml.QNode):
             # Input is a QNode:
             # result = some_transform(qnode, *transform_args)(*qnode_args)
             wrapper = self.qnode_wrapper(qnode, targs, tkwargs)
diff --git a/pennylane/vqe/__init__.py b/pennylane/vqe/__init__.py
deleted file mode 100644
index b0119b021..000000000
--- a/pennylane/vqe/__init__.py
+++ /dev/null
@@ -1,18 +0,0 @@
-# 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 package contains functionality for running Variational Quantum Eigensolver (VQE)
-computations using PennyLane.
-"""
-from .vqe import ExpvalCost
diff --git a/pennylane/vqe/vqe.py b/pennylane/vqe/vqe.py
deleted file mode 100644
index 328618086..000000000
--- a/pennylane/vqe/vqe.py
+++ /dev/null
@@ -1,225 +0,0 @@
-# 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 submodule contains functionality for running Variational Quantum Eigensolver (VQE)
-computations using PennyLane.
-"""
-import warnings
-
-# pylint: disable=too-many-arguments, too-few-public-methods
-from collections.abc import Sequence
-
-import pennylane as qml
-
-
-class ExpvalCost:
-    """
-    Create a cost function that gives the expectation value of an input Hamiltonian.
-
-    This cost function is useful for a range of problems including VQE and QAOA.
-
-    .. warning::
-        ``ExpvalCost`` is deprecated. Instead, it is recommended to simply
-        pass Hamiltonians to the :func:`~pennylane.expval` function inside QNodes.
-
-        .. code-block:: python
-
-            @qml.qnode(dev)
-            def ansatz(params):
-                some_qfunc(params)
-                return qml.expval(Hamiltonian)
-
-        In order to optimize the Hamiltonian evaluation taking into account commuting terms, use the ``grouping_type`` keyword in :class:`~.Hamiltonian`.
-
-    Args:
-        ansatz (callable): The ansatz for the circuit before the final measurement step.
-            Note that the ansatz **must** have the following signature:
-
-            .. code-block:: python
-
-                ansatz(params, wires, **kwargs)
-
-            where ``params`` are the trainable weights of the variational circuit,
-            ``wires`` is the wires the circuit acts on, and ``kwargs`` are any additional
-            keyword arguments that need to be passed to the template.
-        hamiltonian (~.Hamiltonian): Hamiltonian operator whose expectation value should be measured
-        device (Device, Sequence[Device]): Corresponding device(s) where the resulting
-            cost function should be executed. This can either be a single device, or a list
-            of devices of length matching the number of terms in the Hamiltonian.
-        interface (str, None): Which interface to use.
-            This affects the types of objects that can be passed to/returned to the cost function.
-            Supports all interfaces supported by the :func:`~.qnode` decorator.
-        diff_method (str, None): The method of differentiation to use with the created cost function.
-            Supports all differentiation methods supported by the :func:`~.qnode` decorator.
-        optimize (bool): Whether to optimize the observables composing the Hamiltonian by
-            separating them into qubit-wise commuting groups. Each group can then be executed
-            within a single QNode, resulting in fewer QNodes to evaluate.
-
-    Returns:
-        callable: a cost function with signature ``cost_fn(params, **kwargs)`` that evaluates
-        the expectation of the Hamiltonian on the provided device(s)
-
-    .. seealso:: :class:`~.Hamiltonian`, :func:`~.molecular_hamiltonian`, :func:`~.map`, :func:`~.dot`
-
-    **Example:**
-
-    To construct an ``ExpvalCost`` cost function, we require a Hamiltonian to measure, and an ansatz
-    for our variational circuit.
-
-    We can construct a Hamiltonian manually,
-
-    .. code-block:: python
-
-        coeffs = [0.2, -0.543]
-        obs = [
-            qml.PauliX(0) @ qml.PauliZ(1) @ qml.PauliY(3),
-            qml.PauliZ(0) @ qml.Hadamard(2)
-        ]
-        H = qml.Hamiltonian(coeffs, obs)
-
-    Alternatively, the :func:`~.molecular_hamiltonian` function from the
-    :doc:`/introduction/chemistry` module can be used to generate a molecular Hamiltonian.
-
-    Once we have our Hamiltonian, we can select an ansatz and construct
-    the cost function.
-
-    >>> ansatz = qml.templates.StronglyEntanglingLayers
-    >>> dev = qml.device("default.qubit", wires=4)
-    >>> cost = qml.ExpvalCost(ansatz, H, dev, interface="torch")
-    >>> params = torch.rand([2, 4, 3])
-    >>> cost(params)
-    tensor(-0.2316, dtype=torch.float64)
-
-    The cost function can then be minimized using any gradient descent-based
-    :doc:`optimizer </introduction/interfaces>`.
-
-    .. details::
-        :title: Usage Details
-
-        **Optimizing observables:**
-
-        Setting ``optimize=True`` can be used to decrease the number of device executions. The
-        observables composing the Hamiltonian can be separated into groups that are qubit-wise
-        commuting using the :mod:`~.grouping` module. These groups can be executed together on a
-        *single* qnode, resulting in a lower device overhead:
-
-        .. code-block:: python
-
-            commuting_obs = [qml.PauliX(0), qml.PauliX(0) @ qml.PauliZ(1)]
-            H = qml.Hamiltonian([1, 1], commuting_obs)
-
-            dev = qml.device("default.qubit", wires=2)
-            ansatz = qml.templates.StronglyEntanglingLayers
-
-            cost_opt = qml.ExpvalCost(ansatz, H, dev, optimize=True)
-            cost_no_opt = qml.ExpvalCost(ansatz, H, dev, optimize=False)
-
-            shape = qml.templates.StronglyEntanglingLayers.shape(n_layers=3, n_wires=2)
-            params = np.random.random(shape)
-
-        Grouping these commuting observables leads to fewer device executions:
-
-        >>> with qml.Tracker(dev) as tracker:
-        ...     cost_opt(params)
-        >>> ex_opt = tracker.totals["executions"]
-        >>> with tracker:
-        ...     cost_no_opt(params)
-        >>> ex_no_opt = tracker.totals["executions"]
-        >>> print("Number of executions:", ex_no_opt)
-        Number of executions: 2
-        >>> print("Number of executions (optimized):", ex_opt)
-        Number of executions (optimized): 1
-    """
-
-    def __init__(
-        self,
-        ansatz,
-        hamiltonian,
-        device,
-        interface="autograd",
-        diff_method="best",
-        optimize=False,
-        **kwargs,
-    ):
-        warnings.warn(
-            "ExpvalCost is deprecated, use qml.expval() instead. "
-            "For optimizing Hamiltonian measurements with measuring commuting "
-            "terms in parallel, use the grouping_type keyword in qml.Hamiltonian.",
-            qml.PennyLaneDeprecationWarning,
-        )
-
-        if kwargs.get("measure", "expval") != "expval":
-            raise ValueError("ExpvalCost can only be used to construct sums of expectation values.")
-
-        if not callable(ansatz):
-            raise ValueError("Could not create QNodes. The template is not a callable function.")
-
-        self.hamiltonian = hamiltonian
-        """Hamiltonian: the input Hamiltonian."""
-
-        self._multiple_devices = isinstance(device, Sequence)
-        """Bool: Records if multiple devices are input"""
-
-        self.qnodes = []
-        """The QNodes to be evaluated."""
-
-        if optimize:
-            if self._multiple_devices:
-                raise ValueError("Using multiple devices is not supported when optimize=True")
-
-            hamiltonian.compute_grouping()
-
-            @qml.qnode(device, interface=interface, diff_method=diff_method, **kwargs)
-            def circuit(params, **circuit_kwargs):
-                ansatz(params, wires=device.wires, **circuit_kwargs)
-                return qml.expval(hamiltonian)
-
-            self.qnodes.append(circuit)
-            self.cost_fn = self.qnodes[0]
-
-        else:
-            coeffs, observables = hamiltonian.terms()
-            if not self._multiple_devices:
-                device = [device] * len(coeffs)
-
-            for obs, dev in zip(observables, device):
-                wires = dev.wires
-
-                @qml.qnode(
-                    dev,  # pylint: disable=cell-var-from-loop
-                    interface=interface,
-                    diff_method=diff_method,
-                    **kwargs,
-                )
-                def circuit(params, _wires=wires, _obs=obs, **circuit_kwargs):
-                    ansatz(params, wires=_wires, **circuit_kwargs)
-                    return qml.expval(_obs)
-
-                self.qnodes.append(circuit)
-
-            def cost_fn(*args, **kwargs):
-                res = [q(*args, **kwargs) for q in self.qnodes]
-                # pylint: disable=no-member
-                res = [
-                    qml.math.stack(r)
-                    if isinstance(r, (tuple, qml.numpy.builtins.SequenceBox))
-                    else r
-                    for r in res
-                ]
-                return sum(c * q for c, q in zip(coeffs, res))
-
-            self.cost_fn = cost_fn
-
-    def __call__(self, *args, **kwargs):
-        return self.cost_fn(*args, **kwargs)


Test patch:
diff --git a/tests/optimize/test_optimize_shot_adaptive.py b/tests/optimize/test_optimize_shot_adaptive.py
index 06d327459..10dbb7e9c 100644
--- a/tests/optimize/test_optimize_shot_adaptive.py
+++ b/tests/optimize/test_optimize_shot_adaptive.py
@@ -20,14 +20,6 @@ import pennylane as qml
 from pennylane import numpy as np
 
 
-def catch_warn_ExpvalCost(ansatz, hamiltonian, device, **kwargs):
-    """Computes the ExpvalCost and catches the initial deprecation warning."""
-
-    with pytest.warns(UserWarning, match="is deprecated,"):
-        res = qml.ExpvalCost(ansatz, hamiltonian, device, **kwargs)
-    return res
-
-
 class TestExceptions:
     """Test exceptions are raised for incorrect usage"""
 
@@ -36,23 +28,18 @@ class TestExceptions:
         H = qml.Hamiltonian([0.3, 0.1], [qml.PauliX(0), qml.PauliZ(0)])
         dev = qml.device("default.qubit", wires=1, shots=None)
 
-        def ansatz(x):
+        @qml.qnode(dev)
+        def expval_cost(x):
             qml.RX(x, wires=0)
-
-        expval_cost = catch_warn_ExpvalCost(ansatz, H, dev)
+            return qml.expval(H)
 
         opt = qml.ShotAdaptiveOptimizer(min_shots=10)
 
         x = np.array(0.5, requires_grad=True)
 
-        # test expval cost
         with pytest.raises(ValueError, match="can only be used with devices that"):
             opt.step(expval_cost, x)
 
-        # test qnode cost
-        with pytest.raises(ValueError, match="can only be used with devices that"):
-            opt.step(expval_cost.qnodes[0], x)
-
     def test_learning_error(self):
         """Test that an exception is raised if the learning rate is beyond the
         lipschitz bound"""
@@ -60,10 +47,10 @@ class TestExceptions:
         H = qml.Hamiltonian(coeffs, [qml.PauliX(0), qml.PauliZ(0)])
         dev = qml.device("default.qubit", wires=1, shots=100)
 
-        def ansatz(x, **kwargs):
+        @qml.qnode(dev)
+        def expval_cost(x):
             qml.RX(x, wires=0)
-
-        expval_cost = catch_warn_ExpvalCost(ansatz, H, dev)
+            return qml.expval(H)
 
         opt = qml.ShotAdaptiveOptimizer(min_shots=10, stepsize=100.0)
 
@@ -77,10 +64,20 @@ class TestExceptions:
         ):
             opt.step(expval_cost, np.array(0.5, requires_grad=True))
 
-        # for a single QNode, the lipschitz constant is simply 1
-        opt = qml.ShotAdaptiveOptimizer(min_shots=10, stepsize=100.0)
-        with pytest.raises(ValueError, match=f"The learning rate must be less than {2 / 1}"):
-            opt.step(expval_cost.qnodes[0], np.array(0.5, requires_grad=True))
+    def test_compute_grad_no_qnode_error(self):
+        """Test that an exception is raised if a cost_function is not encoded
+        as a QNode Object for compute_grad()"""
+
+        def cost_fn():
+            return None
+
+        opt = qml.ShotAdaptiveOptimizer(min_shots=10)
+        x_init = np.array(0.5, requires_grad=True)
+
+        with pytest.raises(
+            ValueError, match="The objective function must be encoded as a single QNode object"
+        ):
+            opt.compute_grad(cost_fn, [x_init], {})
 
 
 def ansatz0(x, **kwargs):
@@ -105,42 +102,32 @@ class TestSingleShotGradientIntegration:
     for a variety of argument types."""
 
     dev = qml.device("default.qubit", wires=1, shots=100)
-    H = qml.Hamiltonian([1.0], [qml.PauliZ(0)])
-
-    expval_cost = catch_warn_ExpvalCost(ansatz0, H, dev)
 
+    @staticmethod
     @qml.qnode(dev)
-    def qnode(x):
-        # pylint: disable=no-self-argument
+    def cost_fn0(x):
         qml.RX(x, wires=0)
         return qml.expval(qml.PauliZ(0))
 
     @flaky(max_runs=3)
-    @pytest.mark.parametrize("cost_fn", [qnode, expval_cost])
-    def test_single_argument_step(self, cost_fn, mocker, monkeypatch):
+    def test_single_argument_step(self, mocker, monkeypatch):
         """Test that a simple QNode with a single argument correctly performs an optimization step,
         and that the single-shot gradients generated have the correct shape"""
         # pylint: disable=protected-access
 
         opt = qml.ShotAdaptiveOptimizer(min_shots=10)
-        spy_single_shot_expval = mocker.spy(opt, "_single_shot_expval_gradients")
         spy_single_shot_qnodes = mocker.spy(opt, "_single_shot_qnode_gradients")
         spy_grad = mocker.spy(opt, "compute_grad")
 
         x_init = np.array(0.5, requires_grad=True)
-        new_x = opt.step(cost_fn, x_init)
+        new_x = opt.step(self.cost_fn0, x_init)
 
         assert isinstance(new_x, np.tensor)
         assert new_x != x_init
 
         spy_grad.assert_called_once()
-
-        if isinstance(cost_fn, qml.ExpvalCost):
-            spy_single_shot_expval.assert_called_once()
-            single_shot_grads = opt._single_shot_expval_gradients(cost_fn, [x_init], {})
-        else:
-            spy_single_shot_qnodes.assert_called_once()
-            single_shot_grads = opt._single_shot_qnode_gradients(cost_fn, [x_init], {})
+        spy_single_shot_qnodes.assert_called_once()
+        single_shot_grads = opt._single_shot_qnode_gradients(self.cost_fn0, [x_init], {})
 
         # assert single shot gradients are computed correctly
         assert len(single_shot_grads) == 1
@@ -151,15 +138,12 @@ class TestSingleShotGradientIntegration:
         monkeypatch.setattr(
             opt, "_single_shot_qnode_gradients", lambda *args, **kwargs: single_shot_grads
         )
-        monkeypatch.setattr(
-            opt, "_single_shot_expval_gradients", lambda *args, **kwargs: single_shot_grads
-        )
 
         # reset the shot budget
         opt.s = [np.array(10)]
 
         # check that the gradient and variance are computed correctly
-        grad, grad_variance = opt.compute_grad(cost_fn, [x_init], {})
+        grad, grad_variance = opt.compute_grad(self.cost_fn0, [x_init], {})
         assert len(grad) == 1
         assert len(grad_variance) == 1
         assert np.allclose(grad, np.mean(single_shot_grads))
@@ -168,38 +152,34 @@ class TestSingleShotGradientIntegration:
         # check that the gradient and variance are computed correctly
         # with a different shot budget
         opt.s = [np.array(5)]
-        grad, grad_variance = opt.compute_grad(cost_fn, [x_init], {})
+        grad, grad_variance = opt.compute_grad(self.cost_fn0, [x_init], {})
         assert len(grad) == 1
         assert len(grad_variance) == 1
         assert np.allclose(grad, np.mean(single_shot_grads[0][:5]))
         assert np.allclose(grad_variance, np.var(single_shot_grads[0][:5], ddof=1))
 
-    expval_cost = catch_warn_ExpvalCost(ansatz1, H, dev)
-    qnode = expval_cost.qnodes[0]
+    @staticmethod
+    @qml.qnode(dev)
+    def cost_fn1(params):
+        ansatz1(params)
+        return qml.expval(qml.PauliZ(0))
 
-    @pytest.mark.parametrize("cost_fn", [qnode, expval_cost])
-    def test_single_array_argument_step(self, cost_fn, mocker, monkeypatch):
+    def test_single_array_argument_step(self, mocker, monkeypatch):
         """Test that a simple QNode with a single array argument correctly performs an optimization step,
         and that the single-shot gradients generated have the correct shape"""
         # pylint: disable=protected-access
         opt = qml.ShotAdaptiveOptimizer(min_shots=10)
-        spy_single_shot_expval = mocker.spy(opt, "_single_shot_expval_gradients")
         spy_single_shot_qnodes = mocker.spy(opt, "_single_shot_qnode_gradients")
         spy_grad = mocker.spy(opt, "compute_grad")
 
         x_init = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
-        new_x = opt.step(cost_fn, x_init)
+        new_x = opt.step(self.cost_fn1, x_init)
 
         assert isinstance(new_x, np.ndarray)
         assert not np.allclose(new_x, x_init)
 
-        if isinstance(cost_fn, qml.ExpvalCost):
-            spy_single_shot_expval.assert_called_once()
-            single_shot_grads = opt._single_shot_expval_gradients(cost_fn, [x_init], {})
-        else:
-            spy_single_shot_qnodes.assert_called_once()
-            single_shot_grads = opt._single_shot_qnode_gradients(cost_fn, [x_init], {})
-
+        spy_single_shot_qnodes.assert_called_once()
+        single_shot_grads = opt._single_shot_qnode_gradients(self.cost_fn1, [x_init], {})
         spy_grad.assert_called_once()
 
         # assert single shot gradients are computed correctly
@@ -211,15 +191,12 @@ class TestSingleShotGradientIntegration:
         monkeypatch.setattr(
             opt, "_single_shot_qnode_gradients", lambda *args, **kwargs: single_shot_grads
         )
-        monkeypatch.setattr(
-            opt, "_single_shot_expval_gradients", lambda *args, **kwargs: single_shot_grads
-        )
 
         # reset the shot budget
         opt.s = [10 * np.ones([2, 3], dtype=np.int64)]
 
         # check that the gradient and variance are computed correctly
-        grad, grad_variance = opt.compute_grad(cost_fn, [x_init], {})
+        grad, grad_variance = opt.compute_grad(self.cost_fn1, [x_init], {})
         assert len(grad) == 1
         assert len(grad_variance) == 1
         assert grad[0].shape == x_init.shape
@@ -233,7 +210,7 @@ class TestSingleShotGradientIntegration:
         opt.s[0] = opt.s[0] // 2  # all array elements have a shot budget of 5
         opt.s[0][0, 0] = 8  # set the shot budget of the zeroth element to 8
 
-        grad, grad_variance = opt.compute_grad(cost_fn, [x_init], {})
+        grad, grad_variance = opt.compute_grad(self.cost_fn1, [x_init], {})
         assert len(grad) == 1
         assert len(grad_variance) == 1
 
@@ -247,34 +224,30 @@ class TestSingleShotGradientIntegration:
 
     dev = qml.device("default.qubit", wires=2, shots=100)
 
-    expval_cost = catch_warn_ExpvalCost(ansatz2, H, dev)
-    qnode = expval_cost.qnodes[0]
+    @staticmethod
+    @qml.qnode(dev)
+    def cost_fn2(params):
+        ansatz2(params)
+        return qml.expval(qml.PauliZ(0))
 
-    @pytest.mark.parametrize("cost_fn", [qnode, expval_cost])
-    def test_padded_single_array_argument_step(self, cost_fn, mocker, monkeypatch):
+    def test_padded_single_array_argument_step(self, mocker, monkeypatch):
         """Test that a simple QNode with a single array argument with extra dimensions correctly
         performs an optimization step, and that the single-shot gradients generated have the
         correct shape"""
         # pylint: disable=protected-access
         opt = qml.ShotAdaptiveOptimizer(min_shots=10)
-        spy_single_shot_expval = mocker.spy(opt, "_single_shot_expval_gradients")
         spy_single_shot_qnodes = mocker.spy(opt, "_single_shot_qnode_gradients")
         spy_grad = mocker.spy(opt, "compute_grad")
 
         shape = qml.StronglyEntanglingLayers.shape(n_layers=1, n_wires=2)
         x_init = np.ones(shape) * 0.5
-        new_x = opt.step(cost_fn, x_init)
+        new_x = opt.step(self.cost_fn2, x_init)
 
         assert isinstance(new_x, np.ndarray)
         assert not np.allclose(new_x, x_init)
 
-        if isinstance(cost_fn, qml.ExpvalCost):
-            spy_single_shot_expval.assert_called_once()
-            single_shot_grads = opt._single_shot_expval_gradients(cost_fn, [x_init], {})
-        else:
-            spy_single_shot_qnodes.assert_called_once()
-            single_shot_grads = opt._single_shot_qnode_gradients(cost_fn, [x_init], {})
-
+        spy_single_shot_qnodes.assert_called_once()
+        single_shot_grads = opt._single_shot_qnode_gradients(self.cost_fn2, [x_init], {})
         spy_grad.assert_called_once()
 
         # assert single shot gradients are computed correctly
@@ -286,15 +259,12 @@ class TestSingleShotGradientIntegration:
         monkeypatch.setattr(
             opt, "_single_shot_qnode_gradients", lambda *args, **kwargs: single_shot_grads
         )
-        monkeypatch.setattr(
-            opt, "_single_shot_expval_gradients", lambda *args, **kwargs: single_shot_grads
-        )
 
         # reset the shot budget
         opt.s = [10 * np.ones(shape, dtype=np.int64)]
 
         # check that the gradient and variance are computed correctly
-        grad, grad_variance = opt.compute_grad(cost_fn, [x_init], {})
+        grad, grad_variance = opt.compute_grad(self.cost_fn2, [x_init], {})
         assert len(grad) == 1
         assert len(grad_variance) == 1
         assert grad[0].shape == x_init.shape
@@ -308,7 +278,7 @@ class TestSingleShotGradientIntegration:
         opt.s[0] = opt.s[0] // 2  # all array elements have a shot budget of 5
         opt.s[0][0, 0, 0] = 8  # set the shot budget of the zeroth element to 8
 
-        grad, grad_variance = opt.compute_grad(cost_fn, [x_init], {})
+        grad, grad_variance = opt.compute_grad(self.cost_fn2, [x_init], {})
         assert len(grad) == 1
         assert len(grad_variance) == 1
 
@@ -327,7 +297,7 @@ class TestSingleShotGradientIntegration:
         # Step twice to ensure that `opt.s` does not get reshaped.
         # If it was reshaped, its shape would not match `new_x`
         # and an error would get raised.
-        _ = opt.step(cost_fn, new_x)
+        _ = opt.step(self.cost_fn2, new_x)
 
     def test_multiple_argument_step(self, mocker, monkeypatch):
         """Test that a simple QNode with multiple scalar arguments correctly performs an optimization step,
@@ -477,103 +447,6 @@ class TestSingleShotGradientIntegration:
         assert np.allclose(grad_variance[1][1], np.var(single_shot_grads[1][:2, 1], ddof=1))
 
 
-class TestWeightedRandomSampling:
-    """Tests for weighted random Hamiltonian term sampling"""
-
-    def test_wrs_expval_cost(self, mocker):
-        """Checks that cost functions that are expval costs can
-        make use of weighted random sampling"""
-        coeffs = [0.2, 0.1]
-        dev = qml.device("default.qubit", wires=2, shots=100)
-        H = qml.Hamiltonian(coeffs, [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliZ(1)])
-
-        expval_cost = catch_warn_ExpvalCost(qml.templates.StronglyEntanglingLayers, H, dev)
-        weights = np.random.random(qml.templates.StronglyEntanglingLayers.shape(3, 2))
-
-        opt = qml.ShotAdaptiveOptimizer(min_shots=10, term_sampling="weighted_random_sampling")
-        spy = mocker.spy(opt, "weighted_random_sampling")
-
-        _ = opt.step(expval_cost, weights)
-        spy.assert_called_once()
-
-        grads = opt.weighted_random_sampling(expval_cost.qnodes, coeffs, 10, [0], weights)
-        assert len(grads) == 1
-        assert grads[0].shape == (10, *weights.shape)
-
-    def test_wrs_disabled(self, mocker):
-        """Checks that cost functions that are expval costs can
-        disable use of weighted random sampling"""
-        coeffs = [0.2, 0.1]
-        dev = qml.device("default.qubit", wires=2, shots=100)
-        H = qml.Hamiltonian(coeffs, [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliZ(1)])
-
-        expval_cost = catch_warn_ExpvalCost(qml.templates.StronglyEntanglingLayers, H, dev)
-        weights = np.random.random(qml.templates.StronglyEntanglingLayers.shape(3, 2))
-
-        opt = qml.ShotAdaptiveOptimizer(min_shots=10, term_sampling=None)
-        spy = mocker.spy(opt, "weighted_random_sampling")
-
-        opt.step(expval_cost, weights)
-        spy.assert_not_called()
-
-    def test_unknown_term_sampling_method(self):
-        """Checks that an exception is raised if the term sampling method is unknown"""
-        coeffs = [0.2, 0.1]
-        dev = qml.device("default.qubit", wires=2, shots=100)
-        H = qml.Hamiltonian(coeffs, [qml.PauliZ(0), qml.PauliZ(0) @ qml.PauliZ(1)])
-
-        expval_cost = catch_warn_ExpvalCost(qml.templates.StronglyEntanglingLayers, H, dev)
-        weights = np.random.random(qml.templates.StronglyEntanglingLayers.shape(3, 2))
-
-        opt = qml.ShotAdaptiveOptimizer(min_shots=10, term_sampling="uniform_random_sampling")
-
-        with pytest.raises(ValueError, match="Unknown Hamiltonian term sampling method"):
-            opt.step(expval_cost, weights)
-
-    def test_zero_shots(self, mocker):
-        """Test that, if the shot budget for a single term is 0,
-        that the jacobian computation is skipped"""
-        coeffs = [0.2, 0.1, 0.1]
-        dev = qml.device("default.qubit", wires=2, shots=100)
-        H = qml.Hamiltonian(coeffs, [qml.PauliZ(0), qml.PauliX(1), qml.PauliZ(0) @ qml.PauliZ(1)])
-
-        expval_cost = catch_warn_ExpvalCost(qml.templates.StronglyEntanglingLayers, H, dev)
-        weights = np.random.random(qml.templates.StronglyEntanglingLayers.shape(3, 2))
-
-        opt = qml.ShotAdaptiveOptimizer(min_shots=10)
-        spy = mocker.spy(qml, "jacobian")
-        mocker.patch(
-            "scipy.stats._multivariate.multinomial_gen.rvs", return_value=np.array([[4, 0, 6]])
-        )
-        grads = opt.weighted_random_sampling(expval_cost.qnodes, coeffs, 10, [0], weights)
-
-        assert len(spy.call_args_list) == 2
-        assert len(grads) == 1
-        assert grads[0].shape == (10, *weights.shape)
-
-    def test_single_shots(self, mocker):
-        """Test that, if the shot budget for a single term is 1,
-        that the number of dimensions for the returned Jacobian is expanded"""
-        coeffs = [0.2, 0.1, 0.1]
-        dev = qml.device("default.qubit", wires=2, shots=100)
-        H = qml.Hamiltonian(coeffs, [qml.PauliZ(0), qml.PauliX(1), qml.PauliZ(0) @ qml.PauliZ(1)])
-
-        expval_cost = catch_warn_ExpvalCost(qml.templates.StronglyEntanglingLayers, H, dev)
-        weights = np.random.random(qml.templates.StronglyEntanglingLayers.shape(3, 2))
-
-        opt = qml.ShotAdaptiveOptimizer(min_shots=10)
-
-        spy = mocker.spy(qml, "jacobian")
-        mocker.patch(
-            "scipy.stats._multivariate.multinomial_gen.rvs", return_value=np.array([[4, 1, 5]])
-        )
-        grads = opt.weighted_random_sampling(expval_cost.qnodes, coeffs, 10, [0], weights)
-
-        assert len(spy.call_args_list) == 3
-        assert len(grads) == 1
-        assert grads[0].shape == (10, *weights.shape)
-
-
 class TestQNodeWeightedRandomSampling:
     """Tests for weighted random Hamiltonian term sampling"""
 
@@ -770,7 +643,11 @@ class TestOptimization:
             qml.Rot(*x[3], wires=1)
             qml.CNOT(wires=[0, 1])
 
-        cost = catch_warn_ExpvalCost(ansatz, H, dev)
+        @qml.qnode(dev)
+        def cost(params):
+            ansatz(params)
+            return qml.expval(H)
+
         params = np.random.random((4, 3), requires_grad=True)
         initial_loss = cost(params)
 
@@ -788,6 +665,7 @@ class TestOptimization:
 
 
 class TestStepAndCost:
+    # pylint: disable=too-few-public-methods
     """Tests for the step_and_cost method"""
 
     @flaky(max_runs=3)
@@ -812,24 +690,3 @@ class TestStepAndCost:
             params, res = opt.step_and_cost(circuit, params)
 
         assert np.allclose(res, -1, atol=tol, rtol=0)
-
-    @flaky(max_runs=3)
-    def test_expval_cost(self, tol):
-        """Test that the cost is correctly returned
-        when using a QNode as the cost function"""
-        np.random.seed(0)
-        dev = qml.device("default.qubit", wires=1, shots=10)
-
-        def ansatz(x, **kwargs):
-            qml.RX(x[0], wires=0)
-            qml.RY(x[1], wires=0)
-
-        H = qml.Hamiltonian([1.0], [qml.PauliZ(0)])
-        circuit = catch_warn_ExpvalCost(ansatz, H, dev)
-        params = np.array([0.1, 0.3], requires_grad=True)
-        opt = qml.ShotAdaptiveOptimizer(min_shots=10)
-
-        for _ in range(100):
-            params, res = opt.step_and_cost(circuit, params)
-
-        assert np.allclose(res, -1, atol=tol, rtol=0)
diff --git a/tests/optimize/test_qng.py b/tests/optimize/test_qng.py
index b541b6609..4afc0e543 100644
--- a/tests/optimize/test_qng.py
+++ b/tests/optimize/test_qng.py
@@ -42,7 +42,7 @@ class TestExceptions:
 
         with pytest.raises(
             ValueError,
-            match="The objective function must either be encoded as a single QNode or an ExpvalCost object",
+            match="The objective function must be encoded as a single QNode.",
         ):
             opt.step(cost, params)
 
@@ -114,15 +114,18 @@ class TestOptimize:
         opt = qml.QNGOptimizer(stepsize=0.01)
 
         # With autograd gradient function
-        grad_fn = qml.grad(circuit)
-        step1, cost1 = opt.step_and_cost(circuit, var, grad_fn=grad_fn)
-        step2 = opt.step(circuit, var, grad_fn=grad_fn)
+        grad_fn1 = qml.grad(circuit)
+        step1, cost1 = opt.step_and_cost(circuit, var, grad_fn=grad_fn1)
+        step2 = opt.step(circuit, var, grad_fn=grad_fn1)
 
         # With more custom gradient function, forward has to be computed explicitly.
-        grad_fn = lambda param: np.array(qml.grad(circuit)(param))
-        step3, cost2 = opt.step_and_cost(circuit, var, grad_fn=grad_fn)
-        opt.step(circuit, var, grad_fn=grad_fn)
-        expected_step = var - opt.stepsize * 4 * grad_fn(var)
+        def grad_fn2(param):
+            return np.array(qml.grad(circuit)(param))
+
+        # grad_fn = lambda param: np.array(qml.grad(circuit)(param))
+        step3, cost2 = opt.step_and_cost(circuit, var, grad_fn=grad_fn2)
+        opt.step(circuit, var, grad_fn=grad_fn2)
+        expected_step = var - opt.stepsize * 4 * grad_fn2(var)
         expected_cost = circuit(var)
 
         for step in [step1, step2, step3, step3]:
@@ -147,15 +150,17 @@ class TestOptimize:
         opt = qml.QNGOptimizer(stepsize=0.01)
 
         # With autograd gradient function
-        grad_fn = qml.grad(circuit)
-        step1, cost1 = opt.step_and_cost(circuit, *var, grad_fn=grad_fn)
-        step2 = opt.step(circuit, *var, grad_fn=grad_fn)
+        grad_fn1 = qml.grad(circuit)
+        step1, cost1 = opt.step_and_cost(circuit, *var, grad_fn=grad_fn1)
+        step2 = opt.step(circuit, *var, grad_fn=grad_fn1)
 
         # With more custom gradient function, forward has to be computed explicitly.
-        grad_fn = lambda params_0, params_1: np.array(qml.grad(circuit)(params_0, params_1))
-        step3, cost2 = opt.step_and_cost(circuit, *var, grad_fn=grad_fn)
-        opt.step(circuit, *var, grad_fn=grad_fn)
-        expected_step = var - opt.stepsize * 4 * grad_fn(*var)
+        def grad_fn2(params_0, params_1):
+            return np.array(qml.grad(circuit)(params_0, params_1))
+
+        step3, cost2 = opt.step_and_cost(circuit, *var, grad_fn=grad_fn2)
+        opt.step(circuit, *var, grad_fn=grad_fn2)
+        expected_step = var - opt.stepsize * 4 * grad_fn2(*var)
         expected_cost = circuit(*var)
 
         for step in [step1, step2, step3, step3]:
diff --git a/tests/qchem/of_tests/test_convert.py b/tests/qchem/of_tests/test_convert.py
index 95b942fb9..1b9786d50 100644
--- a/tests/qchem/of_tests/test_convert.py
+++ b/tests/qchem/of_tests/test_convert.py
@@ -34,14 +34,6 @@ pauli_ops_and_prod = (qml.PauliX, qml.PauliY, qml.PauliZ, qml.Identity, qml.ops.
 pauli_ops_and_tensor = (qml.PauliX, qml.PauliY, qml.PauliZ, qml.Identity, qml.operation.Tensor)
 
 
-def catch_warn_ExpvalCost(ansatz, hamiltonian, device, **kwargs):
-    """Computes the ExpvalCost and catches the initial deprecation warning."""
-
-    with pytest.warns(UserWarning, match="is deprecated,"):
-        res = qml.ExpvalCost(ansatz, hamiltonian, device, **kwargs)
-    return res
-
-
 @pytest.fixture(
     scope="module",
     params=[
@@ -678,7 +670,7 @@ def test_pennylane_to_openfermion_no_decomp():
 def test_integration_observable_to_vqe_cost(
     monkeypatch, _, terms_ref, expected_cost, custom_wires, tol
 ):
-    r"""Test if `import_operator()` integrates with `ExpvalCost()` in pennylane"""
+    r"""Test if `import_operator()` integrates with `QNode` in pennylane"""
 
     qOp = openfermion.QubitOperator()
     if terms_ref is not None:
@@ -686,7 +678,6 @@ def test_integration_observable_to_vqe_cost(
     vqe_observable = qml.qchem.convert.import_operator(qOp, "openfermion", custom_wires)
 
     num_qubits = len(vqe_observable.wires)
-    assert vqe_observable.terms.__repr__()  # just to satisfy codecov
 
     if custom_wires is None:
         wires = num_qubits
@@ -694,14 +685,15 @@ def test_integration_observable_to_vqe_cost(
         wires = qml.qchem.convert._process_wires(custom_wires)
     else:
         wires = custom_wires[:num_qubits]
+
     dev = qml.device("default.qubit", wires=wires)
 
-    # can replace the ansatz with more suitable ones later.
-    def dummy_ansatz(phis, wires):
-        for phi, w in zip(phis, wires):
+    @qml.qnode(dev)
+    def dummy_cost(params):
+        for phi, w in zip(params, dev.wires):
             qml.RX(phi, wires=w)
+        return qml.expval(vqe_observable)
 
-    dummy_cost = catch_warn_ExpvalCost(dummy_ansatz, vqe_observable, dev)
     params = [0.1 * i for i in range(num_qubits)]
     res = dummy_cost(params)
 
@@ -776,7 +768,7 @@ def test_integration_mol_file_to_vqe_cost(
     name, core, active, mapping, expected_cost, custom_wires, tol
 ):
     r"""Test if the output of `decompose()` works with `import_operator()`
-    to generate `ExpvalCost()`"""
+    to generate `QNode`"""
     ref_dir = os.path.join(os.path.dirname(os.path.realpath(__file__)), "test_ref_files")
     hf_file = os.path.join(ref_dir, name)
     qubit_hamiltonian = qchem.decompose(
@@ -800,16 +792,16 @@ def test_integration_mol_file_to_vqe_cost(
         wires = qml.qchem.convert._process_wires(custom_wires)
     else:
         wires = custom_wires[:num_qubits]
+
     dev = qml.device("default.qubit", wires=wires)
+    phis = np.load(os.path.join(ref_dir, "dummy_ansatz_parameters.npy"))
 
-    # can replace the ansatz with more suitable ones later.
-    def dummy_ansatz(phis, wires):
-        for phi, w in zip(phis, wires):
+    @qml.qnode(dev)
+    def dummy_cost(params):
+        for phi, w in zip(params, dev.wires):
             qml.RX(phi, wires=w)
+        return qml.expval(vqe_hamiltonian)
 
-    phis = np.load(os.path.join(ref_dir, "dummy_ansatz_parameters.npy"))
-
-    dummy_cost = catch_warn_ExpvalCost(dummy_ansatz, vqe_hamiltonian, dev)
     res = dummy_cost(phis)
 
     assert np.abs(res - expected_cost) < tol["atol"]
diff --git a/tests/qchem/of_tests/test_dipole_of.py b/tests/qchem/of_tests/test_dipole_of.py
index 22b5a201f..aee3b527d 100644
--- a/tests/qchem/of_tests/test_dipole_of.py
+++ b/tests/qchem/of_tests/test_dipole_of.py
@@ -242,18 +242,16 @@ def test_dipole_obs(symbols, coords, charge, core, active, mapping, coeffs, ops,
 def test_dipole(symbols, coords, charge, hf_state, exp_dipole, tol, tmpdir):
     r"""Tests the correctness of the computed dipole moment."""
 
-    n_qubits = len(hf_state)
-
-    dev = qml.device("default.qubit", wires=n_qubits)
+    dev = qml.device("default.qubit")
 
     dip_obs = qml.qchem.dipole_of(symbols, coords, charge=charge, outpath=tmpdir.strpath)
 
-    def circuit(params, wires):
-        # pylint: disable=unused-argument
-        qml.BasisState(hf_state, wires=wires)
+    @qml.qnode(dev)
+    def circuit(hf_state, obs):
+        qml.BasisState(hf_state, wires=range(len(hf_state)))
+        return qml.expval(obs)
 
-    with pytest.warns(UserWarning, match="is deprecated,"):
-        dipole = np.array([qml.ExpvalCost(circuit, obs, dev)(None) for obs in dip_obs])
+    dipole = np.array([circuit(hf_state, obs) for obs in dip_obs])
 
     assert np.allclose(dipole, exp_dipole, **tol)
 
diff --git a/tests/test_qaoa.py b/tests/test_qaoa.py
index 2913bcd0d..2854d79ed 100644
--- a/tests/test_qaoa.py
+++ b/tests/test_qaoa.py
@@ -74,14 +74,6 @@ b_rx.add_nodes_from(["b", 1, 0.3])
 b_rx.add_edges_from([(0, 1, ""), (1, 2, ""), (0, 2, "")])
 
 
-def catch_warn_ExpvalCost(ansatz, hamiltonian, device, **kwargs):
-    """Computes the ExpvalCost and catches the initial deprecation warning."""
-
-    with pytest.warns(UserWarning, match="is deprecated,"):
-        res = qml.ExpvalCost(ansatz, hamiltonian, device, **kwargs)
-    return res
-
-
 def decompose_hamiltonian(hamiltonian):
     coeffs = list(qml.math.toarray(hamiltonian.coeffs))
     ops = [i.name for i in hamiltonian.ops]
@@ -1206,7 +1198,11 @@ class TestIntegration:
 
         # Defines the device and the QAOA cost function
         dev = qml.device("default.qubit", wires=len(wires))
-        cost_function = catch_warn_ExpvalCost(circuit, cost_h, dev)
+
+        @qml.qnode(dev)
+        def cost_function(params):
+            circuit(params)
+            return qml.expval(cost_h)
 
         res = cost_function([[1, 1], [1, 1]])
         expected = -1.8260274380964299
@@ -1241,7 +1237,11 @@ class TestIntegration:
 
         # Defines the device and the QAOA cost function
         dev = qml.device("default.qubit", wires=len(wires))
-        cost_function = catch_warn_ExpvalCost(circuit, cost_h, dev)
+
+        @qml.qnode(dev)
+        def cost_function(params):
+            circuit(params)
+            return qml.expval(cost_h)
 
         res = cost_function([[1, 1], [1, 1]])
         expected = -1.8260274380964299
@@ -1873,7 +1873,10 @@ class TestCycles:
         def states(basis_state, **kwargs):
             qml.BasisState(basis_state, wires=range(wires))
 
-        cost = catch_warn_ExpvalCost(states, h, dev, optimize=True)
+        @qml.qnode(dev)
+        def cost(params):
+            states(params)
+            return qml.expval(h)
 
         # Calculate the set of all bitstrings
         bitstrings = itertools.product([0, 1], repeat=wires)
@@ -1920,11 +1923,10 @@ class TestCycles:
         # We use PL to find the energies corresponding to each possible bitstring
         dev = qml.device("default.qubit", wires=wires)
 
-        # pylint: disable=unused-argument
-        def energy(basis_state, **kwargs):
+        @qml.qnode(dev)
+        def cost(basis_state):
             qml.BasisState(basis_state, wires=range(wires))
-
-        cost = catch_warn_ExpvalCost(energy, h, dev, optimize=True)
+            return qml.expval(h)
 
         # Calculate the set of all bitstrings
         states = itertools.product([0, 1], repeat=wires)
@@ -1999,7 +2001,10 @@ class TestCycles:
         def states(basis_state, **kwargs):
             qml.BasisState(basis_state, wires=range(wires))
 
-        cost = catch_warn_ExpvalCost(states, h, dev, optimize=True)
+        @qml.qnode(dev)
+        def cost(params):
+            states(params)
+            return qml.expval(h)
 
         # Calculate the set of all bitstrings
         bitstrings = itertools.product([0, 1], repeat=wires)
diff --git a/tests/test_vqe.py b/tests/test_vqe.py
index 0b1834fb1..1f274c298 100644
--- a/tests/test_vqe.py
+++ b/tests/test_vqe.py
@@ -29,11 +29,15 @@ def seed():
     np.random.seed(0)
 
 
-def catch_warn_ExpvalCost(ansatz, hamiltonian, device, **kwargs):
-    """Computes the ExpvalCost and catches the initial deprecation warning."""
+def generate_cost_fn(ansatz, hamiltonian, device, **kwargs):
+    """Generates a QNode and computes the expectation value of a cost Hamiltonian with respect
+    to the parameters provided to an ansatz"""
+
+    @qml.qnode(device, **kwargs)
+    def res(params):
+        ansatz(params, wires=device.wires)
+        return qml.expval(hamiltonian)
 
-    with pytest.warns(UserWarning, match="is deprecated,"):
-        res = qml.ExpvalCost(ansatz, hamiltonian, device, **kwargs)
     return res
 
 
@@ -60,8 +64,6 @@ OBSERVABLES_NO_HERMITIAN = [
     (qml.PauliX(0) @ qml.PauliZ(1), qml.PauliY(0) @ qml.PauliZ(1), qml.PauliZ(1)),
 ]
 
-JUNK_INPUTS = [None, [], tuple(), 5.0, {"junk": -1}]
-
 hamiltonians_with_expvals = [
     ((-0.6,), (qml.PauliZ(0),), [-0.6 * 1.0]),
     ((1.0,), (qml.PauliX(0),), [0.0]),
@@ -273,7 +275,7 @@ class TestVQE:
         """Tests that the cost function evaluates properly"""
         hamiltonian = qml.Hamiltonian(coeffs, observables)
         dev = qml.device("default.qubit", wires=3)
-        expval = catch_warn_ExpvalCost(ansatz, hamiltonian, dev)
+        expval = generate_cost_fn(ansatz, hamiltonian, dev)
         assert expval(params).dtype == np.float64
         assert np.shape(expval(params)) == ()  # expval should be scalar
 
@@ -284,60 +286,38 @@ class TestVQE:
         """Tests that the cost function returns correct expectation values"""
         dev = qml.device("default.qubit", wires=2)
         hamiltonian = qml.Hamiltonian(coeffs, observables)
-        cost = catch_warn_ExpvalCost(lambda params, **kwargs: None, hamiltonian, dev)
+        cost = generate_cost_fn(lambda params, **kwargs: None, hamiltonian, dev)
         assert cost([]) == sum(expected)
 
-    @pytest.mark.parametrize("ansatz", JUNK_INPUTS)
-    def test_cost_invalid_ansatz(self, ansatz, mock_device):
-        """Tests that the cost function raises an exception if the ansatz is not valid"""
-        hamiltonian = qml.Hamiltonian((1.0,), [qml.PauliZ(0)])
-        with pytest.raises(ValueError, match="not a callable function."):
-            catch_warn_ExpvalCost(4, hamiltonian, mock_device())
-
-    @pytest.mark.autograd
-    @pytest.mark.parametrize("coeffs, observables, expected", hamiltonians_with_expvals)
-    def test_passing_kwargs(self, coeffs, observables, expected):
-        """Test that the step size and order used for the finite differences
-        differentiation method were passed to the QNode instances using the
-        keyword arguments."""
-        dev = qml.device("default.qubit", wires=2)
-        hamiltonian = qml.Hamiltonian(coeffs, observables)
-        cost = catch_warn_ExpvalCost(
-            lambda params, **kwargs: None, hamiltonian, dev, h=123, order=2
-        )
-
-        # Checking that the qnodes contain the step size and order
-        for qnode in cost.qnodes:
-            assert qnode.gradient_kwargs["h"] == 123
-            assert qnode.gradient_kwargs["order"] == 2
-
     # pylint: disable=protected-access
     @pytest.mark.torch
     @pytest.mark.slow
     @pytest.mark.parametrize("dev_name", ["default.qubit", "default.qubit.legacy"])
     @pytest.mark.parametrize("shots", [None, [(8000, 5)], [(8000, 5), (9000, 4)]])
     def test_optimize_torch(self, dev_name, shots):
-        """Test that an ExpvalCost with observable optimization gives the same result as another
-        ExpvalCost without observable optimization."""
+        """Test that a Hamiltonian cost function is the same with and without
+        grouping optimization when using the Torch interface."""
+
+        if dev_name == "default.qubit.legacy" and shots is None:
+            pytest.xfail(reason="DQ legacy does not count hardware executions in analytic mode")
 
         dev = qml.device(dev_name, wires=4, shots=shots)
 
         hamiltonian1 = copy.copy(big_hamiltonian)
         hamiltonian2 = copy.copy(big_hamiltonian)
+        hamiltonian1.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
+        cost = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian1,
             dev,
-            optimize=True,
             interface="torch",
             diff_method="parameter-shift",
         )
-        cost2 = catch_warn_ExpvalCost(
+        cost2 = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian2,
             dev,
-            optimize=False,
             interface="torch",
             diff_method="parameter-shift",
         )
@@ -366,27 +346,26 @@ class TestVQE:
     @pytest.mark.slow
     @pytest.mark.parametrize("shots", [None, [(8000, 5)], [(8000, 5), (9000, 4)]])
     def test_optimize_tf(self, shots):
-        """Test that an ExpvalCost with observable optimization gives the same result as another
-        ExpvalCost without observable optimization."""
+        """Test that a Hamiltonian cost function is the same with and without
+        grouping optimization when using the TensorFlow interface."""
 
         dev = qml.device("default.qubit", wires=4, shots=shots)
 
         hamiltonian1 = copy.copy(big_hamiltonian)
         hamiltonian2 = copy.copy(big_hamiltonian)
+        hamiltonian1.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
+        cost = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian1,
             dev,
-            optimize=True,
             interface="tf",
             diff_method="parameter-shift",
         )
-        cost2 = catch_warn_ExpvalCost(
+        cost2 = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian2,
             dev,
-            optimize=False,
             interface="tf",
             diff_method="parameter-shift",
         )
@@ -413,27 +392,26 @@ class TestVQE:
     @pytest.mark.slow
     @pytest.mark.parametrize("shots", [None, [(8000, 5)], [(8000, 5), (9000, 4)]])
     def test_optimize_autograd(self, shots):
-        """Test that an ExpvalCost with observable optimization gives the same result as another
-        ExpvalCost without observable optimization."""
+        """Test that a Hamiltonian cost function is the same with and without
+        grouping optimization when using the autograd interface."""
 
         dev = qml.device("default.qubit", wires=4, shots=shots)
 
         hamiltonian1 = copy.copy(big_hamiltonian)
         hamiltonian2 = copy.copy(big_hamiltonian)
+        hamiltonian1.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
+        cost = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian1,
             dev,
-            optimize=True,
             interface="autograd",
             diff_method="parameter-shift",
         )
-        cost2 = catch_warn_ExpvalCost(
+        cost2 = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian2,
             dev,
-            optimize=False,
             interface="autograd",
             diff_method="parameter-shift",
         )
@@ -458,9 +436,9 @@ class TestVQE:
     # pylint: disable=protected-access
     @pytest.mark.autograd
     def test_optimize_multiple_terms_autograd(self):
-        """Test that an ExpvalCost with observable optimization gives the same
-        result as another ExpvalCost without observable optimization even when there
-        are non-unique Hamiltonian terms."""
+        """Test that a Hamiltonian cost function is the same with and without
+        grouping optimization when using the autograd interface, even when
+        there are non-unique Hamiltonian terms."""
 
         dev = qml.device("default.qubit", wires=5)
         obs = [
@@ -477,20 +455,19 @@ class TestVQE:
         coefs = (np.random.rand(len(obs)) - 0.5) * 2
         hamiltonian1 = qml.Hamiltonian(coefs, obs)
         hamiltonian2 = qml.Hamiltonian(coefs, obs)
+        hamiltonian1.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
+        cost = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian1,
             dev,
-            optimize=True,
             interface="autograd",
             diff_method="parameter-shift",
         )
-        cost2 = catch_warn_ExpvalCost(
+        cost2 = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian2,
             dev,
-            optimize=False,
             interface="autograd",
             diff_method="parameter-shift",
         )
@@ -515,8 +492,8 @@ class TestVQE:
     # pylint: disable=protected-access
     @pytest.mark.torch
     def test_optimize_multiple_terms_torch(self):
-        """Test that an ExpvalCost with observable optimization gives the same
-        result as another ExpvalCost without observable optimization even when there
+        """Test that a Hamiltonian cost function is the same with and without
+        grouping optimization when using the Torch interface, even when there
         are non-unique Hamiltonian terms."""
 
         dev = qml.device("default.qubit", wires=5)
@@ -534,20 +511,19 @@ class TestVQE:
         coefs = (np.random.rand(len(obs)) - 0.5) * 2
         hamiltonian1 = qml.Hamiltonian(coefs, obs)
         hamiltonian2 = qml.Hamiltonian(coefs, obs)
+        hamiltonian1.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
+        cost = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian1,
             dev,
-            optimize=True,
             interface="torch",
             diff_method="parameter-shift",
         )
-        cost2 = catch_warn_ExpvalCost(
+        cost2 = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian2,
             dev,
-            optimize=False,
             interface="torch",
             diff_method="parameter-shift",
         )
@@ -573,9 +549,9 @@ class TestVQE:
     # pylint: disable=protected-access
     @pytest.mark.tf
     def test_optimize_multiple_terms_tf(self):
-        """Test that an ExpvalCost with observable optimization gives the same
-        result as another ExpvalCost without observable optimization even when there
-        are non-unique Hamiltonian terms."""
+        """Test that a Hamiltonian cost function is the same with and without
+        grouping optimization when using the TensorFlow interface, even when
+        there are non-unique Hamiltonian terms."""
 
         dev = qml.device("default.qubit", wires=5)
         obs = [
@@ -592,20 +568,19 @@ class TestVQE:
         coefs = (np.random.rand(len(obs)) - 0.5) * 2
         hamiltonian1 = qml.Hamiltonian(coefs, obs)
         hamiltonian2 = qml.Hamiltonian(coefs, obs)
+        hamiltonian1.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
+        cost = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian1,
             dev,
-            optimize=True,
             interface="tf",
             diff_method="parameter-shift",
         )
-        cost2 = catch_warn_ExpvalCost(
+        cost2 = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian2,
             dev,
-            optimize=False,
             interface="tf",
             diff_method="parameter-shift",
         )
@@ -631,25 +606,25 @@ class TestVQE:
     # pylint: disable=protected-access
     @pytest.mark.autograd
     def test_optimize_grad(self):
-        """Test that the gradient of ExpvalCost is accessible and correct when using observable
-        optimization and the autograd interface."""
+        """Test that the gradient of a Hamiltonian cost function is accessible
+        and correct when using observable grouping optimization and the
+        autograd interface."""
         dev = qml.device("default.qubit", wires=4)
 
         hamiltonian1 = copy.copy(big_hamiltonian)
         hamiltonian2 = copy.copy(big_hamiltonian)
+        hamiltonian1.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
+        cost = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian1,
             dev,
-            optimize=True,
             diff_method="parameter-shift",
         )
-        cost2 = catch_warn_ExpvalCost(
+        cost2 = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian2,
             dev,
-            optimize=False,
             diff_method="parameter-shift",
         )
 
@@ -672,16 +647,18 @@ class TestVQE:
     @pytest.mark.autograd
     @pytest.mark.parametrize("opt", [True, False])
     def test_grad_zero_hamiltonian(self, opt):
-        """Test that the gradient of ExpvalCost is accessible and correct when using observable
-        optimization and the autograd interface with a zero Hamiltonian."""
+        """Test that the gradient of a Hamiltonian cost function is accessible
+        and correct when using observable grouping optimization and the
+        autograd interface, with a zero Hamiltonian."""
         dev = qml.device("default.qubit", wires=4)
         hamiltonian = qml.Hamiltonian([0], [qml.PauliX(0)])
+        if opt:
+            hamiltonian.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
+        cost = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian,
             dev,
-            optimize=opt,
             diff_method="parameter-shift",
         )
 
@@ -695,18 +672,19 @@ class TestVQE:
     @pytest.mark.torch
     @pytest.mark.slow
     def test_optimize_grad_torch(self):
-        """Test that the gradient of ExpvalCost is accessible and correct when using observable
-        optimization and the Torch interface."""
+        """Test that the gradient of a Hamiltonian cost function is accessible
+        and correct when using observable grouping optimization and the Torch
+        interface."""
         import torch
 
         dev = qml.device("default.qubit", wires=4)
         hamiltonian = big_hamiltonian
+        hamiltonian.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
+        cost = generate_cost_fn(
             qml.templates.StronglyEntanglingLayers,
             hamiltonian,
             dev,
-            optimize=True,
             interface="torch",
         )
 
@@ -724,15 +702,17 @@ class TestVQE:
     @pytest.mark.tf
     @pytest.mark.slow
     def test_optimize_grad_tf(self):
-        """Test that the gradient of ExpvalCost is accessible and correct when using observable
-        optimization and the TensorFlow interface."""
+        """Test that the gradient of a Hamiltonian cost function is accessible
+        and correct when using observable grouping optimization and the
+        TensorFlow interface."""
         import tensorflow as tf
 
         dev = qml.device("default.qubit", wires=4)
         hamiltonian = big_hamiltonian
+        hamiltonian.compute_grouping()
 
-        cost = catch_warn_ExpvalCost(
-            qml.templates.StronglyEntanglingLayers, hamiltonian, dev, optimize=True, interface="tf"
+        cost = generate_cost_fn(
+            qml.templates.StronglyEntanglingLayers, hamiltonian, dev, interface="tf"
         )
 
         np.random.seed(1967)
@@ -747,26 +727,6 @@ class TestVQE:
 
         assert np.allclose(dc, big_hamiltonian_grad)
 
-    def test_multiple_devices_opt_true(self):
-        """Test if a ValueError is raised when multiple devices are passed when optimize=True."""
-        dev = [qml.device("default.qubit", wires=2), qml.device("default.qubit", wires=2)]
-
-        h = qml.Hamiltonian([1, 1], [qml.PauliZ(0), qml.PauliZ(1)])
-
-        with pytest.raises(ValueError, match="Using multiple devices is not supported when"):
-            catch_warn_ExpvalCost(qml.templates.StronglyEntanglingLayers, h, dev, optimize=True)
-
-    def test_variance_error(self):
-        """Test that an error is raised if attempting to use ExpvalCost to measure
-        variances"""
-        dev = qml.device("default.qubit", wires=4)
-        hamiltonian = big_hamiltonian
-
-        with pytest.raises(ValueError, match="sums of expectation values"):
-            catch_warn_ExpvalCost(
-                qml.templates.StronglyEntanglingLayers, hamiltonian, dev, measure="var"
-            )
-
 
 # Test data
 np.random.seed(1967)
@@ -1091,7 +1051,7 @@ class TestInterfaces:
         a, b = 0.54, 0.123
         params = np.array([a, b])
 
-        cost = catch_warn_ExpvalCost(ansatz, H, dev, interface=interface)
+        cost = generate_cost_fn(ansatz, H, dev, interface=interface)
         dcost = qml.grad(cost, argnum=[0])
         res = dcost(params)
 
@@ -1120,7 +1080,7 @@ class TestInterfaces:
         a, b = 0.54, 0.123
         params = torch.autograd.Variable(torch.tensor([a, b]), requires_grad=True)
 
-        cost = catch_warn_ExpvalCost(ansatz, H, dev, interface="torch")
+        cost = generate_cost_fn(ansatz, H, dev, interface="torch")
         loss = cost(params)
         loss.backward()
 
@@ -1150,7 +1110,7 @@ class TestInterfaces:
         H = qml.Hamiltonian(coeffs, observables)
         a, b = 0.54, 0.123
         params = tf.Variable([a, b], dtype=tf.float64)
-        cost = catch_warn_ExpvalCost(ansatz, H, dev, interface="tf")
+        cost = generate_cost_fn(ansatz, H, dev, interface="tf")
 
         with tf.GradientTape() as tape:
             loss = cost(params)
@@ -1185,7 +1145,7 @@ class TestInterfaces:
         w = tf.Variable(params)
         ansatz = qml.templates.layers.StronglyEntanglingLayers
 
-        cost = catch_warn_ExpvalCost(ansatz, H, dev, interface="tf")
+        cost = generate_cost_fn(ansatz, H, dev, interface="tf")
 
         with tf.GradientTape() as tape:
             loss = cost(w)
@@ -1196,7 +1156,7 @@ class TestInterfaces:
         w = torch.autograd.Variable(w, requires_grad=True)
         ansatz = qml.templates.layers.StronglyEntanglingLayers
 
-        cost = catch_warn_ExpvalCost(ansatz, H, dev, interface="torch")
+        cost = generate_cost_fn(ansatz, H, dev, interface="torch")
         loss = cost(w)
         loss.backward()
         res_torch = w.grad.numpy()
@@ -1204,7 +1164,7 @@ class TestInterfaces:
         # NumPy interface
         w = params
         ansatz = qml.templates.layers.StronglyEntanglingLayers
-        cost = catch_warn_ExpvalCost(ansatz, H, dev, interface="autograd")
+        cost = generate_cost_fn(ansatz, H, dev, interface="autograd")
         dcost = qml.grad(cost, argnum=[0])
         res = dcost(w)