Repository: qiskit-community/qiskit-machine-learning (language: python)

Issue:
Cleanup and bugfix to support different primitives. (#55)
### Summary
V2 primitives from different dependencies such as qiskit_aer etc. do not require isa_circuits, hence, removed and cleaned up the gradients. 
- PassManager object is not a strict requirement for V2 primitives anymore when using gradients and state_fidelities.
- Updated tests for gradients.
- Removed unsupported functionalities from gradients.
Cleanup and bugfix to support different primitives. (#55)
### Summary
V2 primitives from different dependencies such as qiskit_aer etc. do not require isa_circuits, hence, removed and cleaned up the gradients. 
- PassManager object is not a strict requirement for V2 primitives anymore when using gradients and state_fidelities.
- Updated tests for gradients.
- Removed unsupported functionalities from gradients.


Discussion / hints:
My concern with this is that for a stablebackport it looks like a lot of breaking changes if you were using this code, i.e. for any existing code you may have E.g. some classes are removed, some classes the parameter order is changed e.g. because passmanager was added or changed in order. One parameter has a default added which is different but ok as it being mandatory before any code would have had to state it so that code will still work unlike the other changes I see here.
## Pull Request Test Coverage Report for [Build 11891912016](https://coveralls.io/builds/70917641)


### Details

* **49** of **65**   **(75.38%)**  changed or added relevant lines in **9** files are covered.
* **10** unchanged lines in **4** files lost coverage.
* Overall coverage decreased (**-0.7%**) to **90.496%**

---

|  Changes Missing Coverage | Covered Lines | Changed/Added Lines | % |
| :-----|--------------|--------|---: |
| [qiskit_machine_learning/gradients/param_shift/param_shift_estimator_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fparam_shift%2Fparam_shift_estimator_gradient.py#L128) | 5 | 7 | 71.43%
| [qiskit_machine_learning/neural_networks/estimator_qnn.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fneural_networks%2Festimator_qnn.py#L175) | 8 | 12 | 66.67%
| [qiskit_machine_learning/neural_networks/sampler_qnn.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fneural_networks%2Fsampler_qnn.py#L198) | 7 | 11 | 63.64%
| [qiskit_machine_learning/gradients/param_shift/param_shift_sampler_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fparam_shift%2Fparam_shift_sampler_gradient.py#L104) | 1 | 7 | 14.29%
<!-- | **Total:** | **49** | **65** | **75.38%** | -->

|  Files with Coverage Reduction | New Missed Lines | % |
| :-----|--------------|--: |
| [qiskit_machine_learning/gradients/base/base_estimator_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fbase%2Fbase_estimator_gradient.py#L251) | 1 | 95.45% |
| [qiskit_machine_learning/gradients/param_shift/param_shift_sampler_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fparam_shift%2Fparam_shift_sampler_gradient.py#L110) | 2 | 71.01% |
| [qiskit_machine_learning/gradients/base/base_sampler_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fbase%2Fbase_sampler_gradient.py#L72) | 3 | 92.93% |
| [qiskit_machine_learning/gradients/utils.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Futils.py#L199) | 4 | 92.48% |
<!-- | **Total:** | **10** |  | -->

|  Totals | [![Coverage Status](https://coveralls.io/builds/70917641/badge)](https://coveralls.io/builds/70917641) |
| :-- | --: |
| Change from base [Build 11829406769](https://coveralls.io/builds/70860703): |  -0.7% |
| Covered Lines: | 4304 |
| Relevant Lines: | 4756 |

---
##### 💛  - [Coveralls](https://coveralls.io)

Hi @woodsp-ibm, I preserved the old formatting as you suggested. The changes in this PR might seem major, but, there are minor bug fixes to allow other types of V2 primitives for different gradients, and, a few fixes for docs. I think we should consider this for the stable branch if people need these changes.
My concern with this is that for a stablebackport it looks like a lot of breaking changes if you were using this code, i.e. for any existing code you may have E.g. some classes are removed, some classes the parameter order is changed e.g. because passmanager was added or changed in order. One parameter has a default added which is different but ok as it being mandatory before any code would have had to state it so that code will still work unlike the other changes I see here.
## Pull Request Test Coverage Report for [Build 11891912016](https://coveralls.io/builds/70917641)


### Details

* **49** of **65**   **(75.38%)**  changed or added relevant lines in **9** files are covered.
* **10** unchanged lines in **4** files lost coverage.
* Overall coverage decreased (**-0.7%**) to **90.496%**

---

|  Changes Missing Coverage | Covered Lines | Changed/Added Lines | % |
| :-----|--------------|--------|---: |
| [qiskit_machine_learning/gradients/param_shift/param_shift_estimator_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fparam_shift%2Fparam_shift_estimator_gradient.py#L128) | 5 | 7 | 71.43%
| [qiskit_machine_learning/neural_networks/estimator_qnn.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fneural_networks%2Festimator_qnn.py#L175) | 8 | 12 | 66.67%
| [qiskit_machine_learning/neural_networks/sampler_qnn.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fneural_networks%2Fsampler_qnn.py#L198) | 7 | 11 | 63.64%
| [qiskit_machine_learning/gradients/param_shift/param_shift_sampler_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fparam_shift%2Fparam_shift_sampler_gradient.py#L104) | 1 | 7 | 14.29%
<!-- | **Total:** | **49** | **65** | **75.38%** | -->

|  Files with Coverage Reduction | New Missed Lines | % |
| :-----|--------------|--: |
| [qiskit_machine_learning/gradients/base/base_estimator_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fbase%2Fbase_estimator_gradient.py#L251) | 1 | 95.45% |
| [qiskit_machine_learning/gradients/param_shift/param_shift_sampler_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fparam_shift%2Fparam_shift_sampler_gradient.py#L110) | 2 | 71.01% |
| [qiskit_machine_learning/gradients/base/base_sampler_gradient.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Fbase%2Fbase_sampler_gradient.py#L72) | 3 | 92.93% |
| [qiskit_machine_learning/gradients/utils.py](https://coveralls.io/builds/70917641/source?filename=qiskit_machine_learning%2Fgradients%2Futils.py#L199) | 4 | 92.48% |
<!-- | **Total:** | **10** |  | -->

|  Totals | [![Coverage Status](https://coveralls.io/builds/70917641/badge)](https://coveralls.io/builds/70917641) |
| :-- | --: |
| Change from base [Build 11829406769](https://coveralls.io/builds/70860703): |  -0.7% |
| Covered Lines: | 4304 |
| Relevant Lines: | 4756 |

---
##### 💛  - [Coveralls](https://coveralls.io)

Hi @woodsp-ibm, I preserved the old formatting as you suggested. The changes in this PR might seem major, but, there are minor bug fixes to allow other types of V2 primitives for different gradients, and, a few fixes for docs. I think we should consider this for the stable branch if people need these changes.

Fix patch:
diff --git a/.pylintdict b/.pylintdict
index 70e2bb17f..c059dfefa 100644
--- a/.pylintdict
+++ b/.pylintdict
@@ -390,6 +390,7 @@ platt
 polyfit
 postprocess
 powell
+pragma
 pre
 precompute
 precomputed
diff --git a/qiskit_machine_learning/gradients/__init__.py b/qiskit_machine_learning/gradients/__init__.py
index 9f24daa53..dcc426ae3 100644
--- a/qiskit_machine_learning/gradients/__init__.py
+++ b/qiskit_machine_learning/gradients/__init__.py
@@ -25,11 +25,9 @@
    :nosignatures:
 
    BaseEstimatorGradient
-   BaseQGT
    BaseSamplerGradient
    EstimatorGradientResult
    SamplerGradientResult
-   QGTResult
 
 Linear Combination of Unitaries
 -------------------------------
@@ -40,7 +38,6 @@
 
    LinCombEstimatorGradient
    LinCombSamplerGradient
-   LinCombQGT
 
 Parameter Shift Rules
 ---------------------
@@ -52,16 +49,6 @@
    ParamShiftEstimatorGradient
    ParamShiftSamplerGradient
 
-Quantum Fisher Information
---------------------------
-
-.. autosummary::
-   :toctree: ../stubs/
-   :nosignatures:
-
-   QFIResult
-   QFI
-
 Simultaneous Perturbation Stochastic Approximation
 --------------------------------------------------
 
@@ -74,35 +61,25 @@
 """
 
 from .base.base_estimator_gradient import BaseEstimatorGradient
-from .base.base_qgt import BaseQGT
 from .base.base_sampler_gradient import BaseSamplerGradient
 from .base.estimator_gradient_result import EstimatorGradientResult
 from .lin_comb.lin_comb_estimator_gradient import DerivativeType, LinCombEstimatorGradient
-from .lin_comb.lin_comb_qgt import LinCombQGT
 from .lin_comb.lin_comb_sampler_gradient import LinCombSamplerGradient
 from .param_shift.param_shift_estimator_gradient import ParamShiftEstimatorGradient
 from .param_shift.param_shift_sampler_gradient import ParamShiftSamplerGradient
-from .qfi import QFI
-from .qfi_result import QFIResult
-from .base.qgt_result import QGTResult
 from .base.sampler_gradient_result import SamplerGradientResult
 from .spsa.spsa_estimator_gradient import SPSAEstimatorGradient
 from .spsa.spsa_sampler_gradient import SPSASamplerGradient
 
 __all__ = [
     "BaseEstimatorGradient",
-    "BaseQGT",
     "BaseSamplerGradient",
     "DerivativeType",
     "EstimatorGradientResult",
     "LinCombEstimatorGradient",
-    "LinCombQGT",
     "LinCombSamplerGradient",
     "ParamShiftEstimatorGradient",
     "ParamShiftSamplerGradient",
-    "QFI",
-    "QFIResult",
-    "QGTResult",
     "SamplerGradientResult",
     "SPSAEstimatorGradient",
     "SPSASamplerGradient",
diff --git a/qiskit_machine_learning/gradients/base/base_estimator_gradient.py b/qiskit_machine_learning/gradients/base/base_estimator_gradient.py
index bb85cd179..2bb0c6735 100644
--- a/qiskit_machine_learning/gradients/base/base_estimator_gradient.py
+++ b/qiskit_machine_learning/gradients/base/base_estimator_gradient.py
@@ -56,7 +56,6 @@ def __init__(
         r"""
         Args:
             estimator: The estimator used to compute the gradients.
-            pass_manager: pass manager for isa_circuit transpilation.
             options: Primitive backend runtime options used for circuit execution.
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
@@ -71,6 +70,8 @@ def __init__(
                 Defaults to ``DerivativeType.REAL``, as this yields e.g. the commonly-used energy
                 gradient and this type is the only supported type for function-level schemes like
                 finite difference.
+            pass_manager: The pass manager to transpile the circuits if necessary.
+            Defaults to ``None``, as some primitives do not need transpiled circuits.
         """
         if isinstance(estimator, BaseEstimatorV1):
             issue_deprecation_msg(
diff --git a/qiskit_machine_learning/gradients/base/base_qgt.py b/qiskit_machine_learning/gradients/base/base_qgt.py
deleted file mode 100644
index 9094a26a5..000000000
--- a/qiskit_machine_learning/gradients/base/base_qgt.py
+++ /dev/null
@@ -1,388 +0,0 @@
-# This code is part of a Qiskit project.
-#
-# (C) Copyright IBM 2022, 2024.
-#
-# This code is licensed under the Apache License, Version 2.0. You may
-# obtain a copy of this license in the LICENSE.txt file in the root directory
-# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
-#
-# Any modifications or derivative works of this code must retain this
-# copyright notice, and modified files need to carry a notice indicating
-# that they have been altered from the originals.
-
-"""
-Abstract base class of the Quantum Geometric Tensor (QGT).
-"""
-
-from __future__ import annotations
-
-from abc import ABC, abstractmethod
-from collections.abc import Sequence
-from copy import copy
-
-import numpy as np
-
-from qiskit.circuit import Parameter, ParameterExpression, QuantumCircuit
-from qiskit.primitives import BaseEstimator
-from qiskit.primitives.utils import _circuit_key
-from qiskit.providers import Options
-from qiskit.transpiler.passes import TranslateParameterizedGates
-
-from .qgt_result import QGTResult
-from ..utils import (
-    DerivativeType,
-    GradientCircuit,
-    _assign_unique_parameters,
-    _make_gradient_parameters,
-    _make_gradient_parameter_values,
-)
-
-from ...algorithm_job import AlgorithmJob
-
-
-class BaseQGT(ABC):
-    r"""Base class to computes the Quantum Geometric Tensor (QGT) given a pure,
-    parameterized quantum state. QGT is defined as:
-
-    .. math::
-
-        \mathrm{QGT}_{ij}= \langle \partial_i \psi | \partial_j \psi \rangle
-            - \langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle.
-    """
-
-    def __init__(
-        self,
-        estimator: BaseEstimator,
-        phase_fix: bool = True,
-        derivative_type: DerivativeType = DerivativeType.COMPLEX,
-        options: Options | None = None,
-    ):
-        r"""
-        Args:
-            estimator: The estimator used to compute the QGT.
-            phase_fix: Whether to calculate the second term (phase fix) of the QGT, which is
-                :math:`\langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle`.
-                Defaults to ``True``.
-            derivative_type: The type of derivative. Can be either ``DerivativeType.REAL``
-                ``DerivativeType.IMAG``, or ``DerivativeType.COMPLEX``. Defaults to
-                ``DerivativeType.REAL``.
-
-                - ``DerivativeType.REAL`` computes
-
-                .. math::
-
-                    \mathrm{Re(QGT)}_{ij}= \mathrm{Re}[\langle \partial_i \psi | \partial_j \psi \rangle
-                        - \langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle].
-
-                - ``DerivativeType.IMAG`` computes
-
-                .. math::
-
-                    \mathrm{Im(QGT)}_{ij}= \mathrm{Im}[\langle \partial_i \psi | \partial_j \psi \rangle
-                        - \langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle].
-
-                - ``DerivativeType.COMPLEX`` computes
-
-                .. math::
-
-                    \mathrm{QGT}_{ij}= [\langle \partial_i \psi | \partial_j \psi \rangle
-                        - \langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle].
-
-            options: Backend runtime options used for circuit execution. The order of priority is:
-                options in ``run`` method > QGT's default options > primitive's default
-                setting. Higher priority setting overrides lower priority setting.
-        """
-        self._estimator: BaseEstimator = estimator
-        self._phase_fix: bool = phase_fix
-        self._derivative_type: DerivativeType = derivative_type
-        self._default_options = Options()
-        if options is not None:
-            self._default_options.update_options(**options)
-        self._qgt_circuit_cache: dict[tuple, GradientCircuit] = {}
-        self._gradient_circuit_cache: dict[tuple, GradientCircuit] = {}
-
-    @property
-    def derivative_type(self) -> DerivativeType:
-        """The derivative type."""
-        return self._derivative_type
-
-    @derivative_type.setter
-    def derivative_type(self, derivative_type: DerivativeType) -> None:
-        """Set the derivative type."""
-        self._derivative_type = derivative_type
-
-    def run(
-        self,
-        circuits: Sequence[QuantumCircuit],
-        parameter_values: Sequence[Sequence[float]],
-        parameters: Sequence[Sequence[Parameter] | None] | None = None,
-        **options,
-    ) -> AlgorithmJob:
-        """Run the job of the QGTs on the given circuits.
-
-        Args:
-            circuits: The list of quantum circuits to compute the QGTs.
-            parameter_values: The list of parameter values to be bound to the circuit.
-            parameters: The sequence of parameters to calculate only the QGTs of
-                the specified parameters. Each sequence of parameters corresponds to a circuit in
-                ``circuits``. Defaults to None, which means that the QGTs of all parameters in
-                each circuit are calculated.
-            options: Primitive backend runtime options used for circuit execution.
-                The order of priority is: options in ``run`` method > QGT's
-                default options > primitive's default setting.
-                Higher priority setting overrides lower priority setting.
-
-        Returns:
-            The job object of the QGTs of the expectation values. The i-th result corresponds to
-            ``circuits[i]`` evaluated with parameters bound as ``parameter_values[i]``.
-
-        Raises:
-            ValueError: Invalid arguments are given.
-        """
-        if isinstance(circuits, QuantumCircuit):
-            # Allow a single circuit to be passed in.
-            circuits = (circuits,)
-
-        if parameters is None:
-            # If parameters is None, we calculate the gradients of all parameters in each circuit.
-            parameters = [circuit.parameters for circuit in circuits]
-        else:
-            # If parameters is not None, we calculate the gradients of the specified parameters.
-            # None in parameters means that the gradients of all parameters in the corresponding
-            # circuit are calculated.
-            parameters = [
-                params if params is not None else circuits[i].parameters
-                for i, params in enumerate(parameters)
-            ]
-        # Validate the arguments.
-        self._validate_arguments(circuits, parameter_values, parameters)
-        # The priority of run option is as follows:
-        # options in ``run`` method > QGT's default options > primitive's default setting.
-        opts = copy(self._default_options)
-        opts.update_options(**options)
-        job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **opts.__dict__)
-        job.submit()
-        return job
-
-    @abstractmethod
-    def _run(
-        self,
-        circuits: Sequence[QuantumCircuit],
-        parameter_values: Sequence[Sequence[float]],
-        parameters: Sequence[Sequence[Parameter]],
-        **options,
-    ) -> QGTResult:
-        """Compute the QGTs on the given circuits."""
-        raise NotImplementedError()
-
-    def _preprocess(
-        self,
-        circuits: Sequence[QuantumCircuit],
-        parameter_values: Sequence[Sequence[float]],
-        parameters: Sequence[Sequence[Parameter]],
-        supported_gates: Sequence[str],
-    ) -> tuple[Sequence[QuantumCircuit], Sequence[Sequence[float]], Sequence[Sequence[Parameter]]]:
-        """Preprocess the gradient. This makes a gradient circuit for each circuit. The gradient
-        circuit is a transpiled circuit by using the supported gates, and has unique parameters.
-        ``parameter_values`` and ``parameters`` are also updated to match the gradient circuit.
-
-        Args:
-            circuits: The list of quantum circuits to compute the gradients.
-            parameter_values: The list of parameter values to be bound to the circuit.
-            parameters: The sequence of parameters to calculate only the gradients of the specified
-                parameters.
-            supported_gates: The supported gates used to transpile the circuit.
-
-        Returns:
-            The list of gradient circuits, the list of parameter values, and the list of parameters.
-            parameter_values and parameters are updated to match the gradient circuit.
-        """
-        translator = TranslateParameterizedGates(supported_gates)
-        g_circuits: list[QuantumCircuit] = []
-        g_parameter_values: list[Sequence[float]] = []
-        g_parameters: list[Sequence[Parameter]] = []
-        for circuit, parameter_value_, parameters_ in zip(circuits, parameter_values, parameters):
-            circuit_key = _circuit_key(circuit)
-            if circuit_key not in self._gradient_circuit_cache:
-                unrolled = translator(circuit)
-                self._gradient_circuit_cache[circuit_key] = _assign_unique_parameters(unrolled)
-            gradient_circuit = self._gradient_circuit_cache[circuit_key]
-            g_circuits.append(gradient_circuit.gradient_circuit)
-            g_parameter_values.append(
-                _make_gradient_parameter_values(  # type: ignore[arg-type]
-                    circuit, gradient_circuit, parameter_value_
-                )
-            )
-            g_parameters_ = [
-                g_param
-                for g_param in gradient_circuit.gradient_circuit.parameters
-                if g_param in _make_gradient_parameters(gradient_circuit, parameters_)
-            ]
-            g_parameters.append(g_parameters_)
-        return g_circuits, g_parameter_values, g_parameters
-
-    def _postprocess(
-        self,
-        results: QGTResult,
-        circuits: Sequence[QuantumCircuit],
-        parameter_values: Sequence[Sequence[float]],
-        parameters: Sequence[Sequence[Parameter]],
-    ) -> QGTResult:
-        """Postprocess the QGTs. This method computes the QGTs of the original circuits
-        by applying the chain rule to the QGTs of the circuits with unique parameters.
-
-        Args:
-            results: The computed QGT for the circuits with unique parameters.
-            circuits: The list of original circuits submitted for gradient computation.
-            parameter_values: The list of parameter values to be bound to the circuits.
-            parameters: The sequence of parameters to calculate only the gradients of the specified
-                parameters.
-
-        Returns:
-            The QGTs of the original circuits.
-        """
-        qgts, metadata = [], []
-        for idx, (circuit, parameter_values_, parameters_) in enumerate(
-            zip(circuits, parameter_values, parameters)
-        ):
-            dtype = complex if self.derivative_type == DerivativeType.COMPLEX else float
-            qgt: np.ndarray = np.zeros((len(parameters_), len(parameters_)), dtype=dtype)
-
-            gradient_circuit = self._gradient_circuit_cache[_circuit_key(circuit)]
-            g_parameters = _make_gradient_parameters(gradient_circuit, parameters_)
-            # Make a map from the gradient parameter to the respective index in the gradient.
-            # parameters_ = [param for param in circuit.parameters if param in parameters_]
-            g_parameter_indices = [
-                param
-                for param in gradient_circuit.gradient_circuit.parameters
-                if param in g_parameters
-            ]
-            g_parameter_indices_d = {param: i for i, param in enumerate(g_parameter_indices)}
-            rows, cols = np.triu_indices(len(parameters_))
-            for row, col in zip(rows, cols):
-                for g_parameter1, coeff1 in gradient_circuit.parameter_map[parameters_[row]]:
-                    for g_parameter2, coeff2 in gradient_circuit.parameter_map[parameters_[col]]:
-                        if isinstance(coeff1, ParameterExpression):
-                            local_map = {
-                                p: parameter_values_[circuit.parameters.data.index(p)]
-                                for p in coeff1.parameters
-                            }
-                            bound_coeff1 = coeff1.bind(local_map)
-                        else:
-                            bound_coeff1 = coeff1
-                        if isinstance(coeff2, ParameterExpression):
-                            local_map = {
-                                p: parameter_values_[circuit.parameters.data.index(p)]
-                                for p in coeff2.parameters
-                            }
-                            bound_coeff2 = coeff2.bind(local_map)
-                        else:
-                            bound_coeff2 = coeff2
-                        qgt[row, col] += (
-                            float(bound_coeff1)
-                            * float(bound_coeff2)
-                            * results.qgts[idx][
-                                g_parameter_indices_d[g_parameter1],
-                                g_parameter_indices_d[g_parameter2],
-                            ]
-                        )
-
-            if self.derivative_type == DerivativeType.IMAG:
-                qgt += -1 * np.triu(qgt, k=1).T
-            else:
-                qgt += np.triu(qgt, k=1).conjugate().T
-            qgts.append(qgt)
-            metadata.append([{"parameters": parameters_}])
-        return QGTResult(
-            qgts=qgts,
-            derivative_type=self.derivative_type,
-            metadata=metadata,
-            options=results.options,
-        )
-
-    @staticmethod
-    def _validate_arguments(
-        circuits: Sequence[QuantumCircuit],
-        parameter_values: Sequence[Sequence[float]],
-        parameters: Sequence[Sequence[Parameter]],
-    ) -> None:
-        """Validate the arguments of the ``run`` method.
-
-        Args:
-            circuits: The list of quantum circuits to compute the QGTs.
-            parameter_values: The list of parameter values to be bound to the circuits.
-            parameters: The sequence of parameters with respect to which the QGTs should be
-                computed.
-
-        Raises:
-            ValueError: Invalid arguments are given.
-        """
-        if len(circuits) != len(parameter_values):
-            raise ValueError(
-                f"The number of circuits ({len(circuits)}) does not match "
-                f"the number of parameter values ({len(parameter_values)})."
-            )
-
-        if len(circuits) != len(parameters):
-            raise ValueError(
-                f"The number of circuits ({len(circuits)}) does not match "
-                f"the number of the specified parameter sets ({len(parameters)})."
-            )
-
-        for i, (circuit, parameter_value) in enumerate(zip(circuits, parameter_values)):
-            if not circuit.num_parameters:
-                raise ValueError(f"The {i}-th circuit is not parameterised.")
-            if len(parameter_value) != circuit.num_parameters:
-                raise ValueError(
-                    f"The number of values ({len(parameter_value)}) does not match "
-                    f"the number of parameters ({circuit.num_parameters}) for the {i}-th circuit."
-                )
-
-        if len(circuits) != len(parameters):
-            raise ValueError(
-                f"The number of circuits ({len(circuits)}) does not match "
-                f"the number of the list of specified parameters ({len(parameters)})."
-            )
-
-        for i, (circuit, parameters_) in enumerate(zip(circuits, parameters)):
-            if not set(parameters_).issubset(circuit.parameters):
-                raise ValueError(
-                    f"The {i}-th parameters contains parameters not present in the "
-                    f"{i}-th circuit."
-                )
-
-    @property
-    def options(self) -> Options:
-        """Return the union of estimator options setting and QGT default options,
-        where, if the same field is set in both, the QGT's default options override
-        the primitive's default setting.
-
-        Returns:
-            The QGT default + estimator options.
-        """
-        return self._get_local_options(self._default_options.__dict__)
-
-    def update_default_options(self, **options):
-        """Update the gradient's default options setting.
-
-        Args:
-            **options: The fields to update the default options.
-        """
-
-        self._default_options.update_options(**options)
-
-    def _get_local_options(self, options: Options) -> Options:
-        """Return the union of the primitive's default setting,
-        the QGT default options, and the options in the ``run`` method.
-        The order of priority is: options in ``run`` method > QGT's default options > primitive's
-        default setting.
-
-        Args:
-            options: The fields to update the options
-
-        Returns:
-            The QGT default + estimator + run options.
-        """
-        opts = copy(self._estimator.options)
-        opts.update_options(**options)
-        return opts
diff --git a/qiskit_machine_learning/gradients/base/base_sampler_gradient.py b/qiskit_machine_learning/gradients/base/base_sampler_gradient.py
index 3db0c3e31..ea8ad98e4 100644
--- a/qiskit_machine_learning/gradients/base/base_sampler_gradient.py
+++ b/qiskit_machine_learning/gradients/base/base_sampler_gradient.py
@@ -46,7 +46,6 @@ def __init__(
         self,
         sampler: BaseSampler,
         options: Options | None = None,
-        len_quasi_dist: int | None = None,
         pass_manager: BasePassManager | None = None,
     ):
         """
@@ -56,6 +55,8 @@ def __init__(
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
                 Higher priority setting overrides lower priority setting
+            pass_manager: The pass manager to transpile the circuits if necessary.
+            Defaults to ``None``, as some primitives do not need transpiled circuits.
         """
         if isinstance(sampler, BaseSamplerV1):
             issue_deprecation_msg(
@@ -66,7 +67,6 @@ def __init__(
             )
         self._sampler: BaseSampler = sampler
         self._pass_manager = pass_manager
-        self._len_quasi_dist = len_quasi_dist
         self._default_options = Options()
         if options is not None:
             self._default_options.update_options(**options)
diff --git a/qiskit_machine_learning/gradients/base/qgt_result.py b/qiskit_machine_learning/gradients/base/qgt_result.py
deleted file mode 100644
index acdb6710e..000000000
--- a/qiskit_machine_learning/gradients/base/qgt_result.py
+++ /dev/null
@@ -1,39 +0,0 @@
-# This code is part of a Qiskit project.
-#
-# (C) Copyright IBM 2022, 2024.
-#
-# This code is licensed under the Apache License, Version 2.0. You may
-# obtain a copy of this license in the LICENSE.txt file in the root directory
-# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
-#
-# Any modifications or derivative works of this code must retain this
-# copyright notice, and modified files need to carry a notice indicating
-# that they have been altered from the originals.
-"""
-QGT result class
-"""
-
-from __future__ import annotations
-
-from dataclasses import dataclass
-from typing import Any
-
-import numpy as np
-
-from qiskit.providers import Options
-
-from ..utils import DerivativeType
-
-
-@dataclass(frozen=True)
-class QGTResult:
-    """Result of QGT."""
-
-    qgts: list[np.ndarray]
-    """The QGT."""
-    derivative_type: DerivativeType
-    """The type of derivative."""
-    metadata: list[dict[str, Any]] | list[list[dict[str, Any]]]
-    """Additional information about the job."""
-    options: Options
-    """Primitive runtime options for the execution of the job."""
diff --git a/qiskit_machine_learning/gradients/lin_comb/lin_comb_estimator_gradient.py b/qiskit_machine_learning/gradients/lin_comb/lin_comb_estimator_gradient.py
index e70876a26..9c0a0e336 100644
--- a/qiskit_machine_learning/gradients/lin_comb/lin_comb_estimator_gradient.py
+++ b/qiskit_machine_learning/gradients/lin_comb/lin_comb_estimator_gradient.py
@@ -12,7 +12,6 @@
 """
 Gradient of probabilities with linear combination of unitaries (LCU)
 """
-
 from __future__ import annotations
 
 from collections.abc import Sequence
@@ -20,7 +19,10 @@
 import numpy as np
 
 from qiskit.circuit import Parameter, QuantumCircuit
-from qiskit.primitives import BaseEstimator
+from qiskit.primitives.base import BaseEstimatorV2
+from qiskit.primitives import BaseEstimator, BaseEstimatorV1
+from qiskit.transpiler.passmanager import BasePassManager
+
 from qiskit.primitives.utils import init_observable, _circuit_key
 from qiskit.providers import Options
 from qiskit.quantum_info.operators.base_operator import BaseOperator
@@ -69,6 +71,7 @@ def __init__(
         estimator: BaseEstimator,
         derivative_type: DerivativeType = DerivativeType.REAL,
         options: Options | None = None,
+        pass_manager: BasePassManager | None = None,
     ):
         r"""
         Args:
@@ -85,9 +88,13 @@ def __init__(
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
                 Higher priority setting overrides lower priority setting.
+            pass_manager: The pass manager to transpile the circuits if necessary.
+                Defaults to ``None``, as some primitives do not need transpiled circuits.
         """
         self._lin_comb_cache: dict[tuple, dict[Parameter, QuantumCircuit]] = {}
-        super().__init__(estimator, options, derivative_type=derivative_type)
+        super().__init__(
+            estimator, options=options, derivative_type=derivative_type, pass_manager=pass_manager
+        )
 
     @BaseEstimatorGradient.derivative_type.setter  # type: ignore[attr-defined]
     def derivative_type(self, derivative_type: DerivativeType) -> None:
@@ -118,7 +125,7 @@ def _run_unique(
         parameter_values: Sequence[Sequence[float]],
         parameters: Sequence[Sequence[Parameter]],
         **options,
-    ) -> EstimatorGradientResult:
+    ) -> EstimatorGradientResult:  # pragma: no cover
         """Compute the estimator gradients on the given circuits."""
         job_circuits, job_observables, job_param_values, metadata = [], [], [], []
         all_n = []
@@ -161,34 +168,79 @@ def _run_unique(
                 job_param_values.extend([parameter_values_] * n)
                 all_n.append(n)
 
-        # Run the single job with all circuits.
-        job = self._estimator.run(
-            job_circuits,
-            job_observables,
-            job_param_values,
-            **options,
-        )
-        try:
-            results = job.result()
-        except AlgorithmError as exc:
-            raise AlgorithmError("Estimator job failed.") from exc
-
-        # Compute the gradients.
-        gradients = []
-        partial_sum_n = 0
-        for n in all_n:
-            # this disable is needed as Pylint does not understand derivative_type is a property if
-            # it is only defined in the base class and the getter is in the child
-            # pylint: disable=comparison-with-callable
-            if self.derivative_type == DerivativeType.COMPLEX:
-                gradient = np.zeros(n // 2, dtype="complex")
-                gradient.real = results.values[partial_sum_n : partial_sum_n + n // 2]
-                gradient.imag = results.values[partial_sum_n + n // 2 : partial_sum_n + n]
-
+        if isinstance(self._estimator, BaseEstimatorV1):
+            # Run the single job with all circuits.
+            job = self._estimator.run(
+                job_circuits,
+                job_observables,
+                job_param_values,
+                **options,
+            )
+            try:
+                results = job.result()
+            except Exception as exc:
+                raise AlgorithmError("Estimator job failed.") from exc
+            # Compute the gradients.
+            gradients = []
+            partial_sum_n = 0
+            for n in all_n:
+                # this disable is needed as Pylint does not understand derivative_type is a property if
+                # it is only defined in the base class and the getter is in the child
+                # pylint: disable=comparison-with-callable
+                if self.derivative_type == DerivativeType.COMPLEX:
+                    gradient = np.zeros(n // 2, dtype="complex")
+                    gradient.real = results.values[partial_sum_n : partial_sum_n + n // 2]
+                    gradient.imag = results.values[partial_sum_n + n // 2 : partial_sum_n + n]
+
+                else:
+                    gradient = np.real(results.values[partial_sum_n : partial_sum_n + n])
+                partial_sum_n += n
+                gradients.append(gradient)
+
+            opt = self._get_local_options(options)
+        elif isinstance(self._estimator, BaseEstimatorV2):
+            if self._pass_manager is None:
+                circs = job_circuits
+                observables = job_observables
             else:
-                gradient = np.real(results.values[partial_sum_n : partial_sum_n + n])
-            partial_sum_n += n
-            gradients.append(gradient)
-
-        opt = self._get_local_options(options)
+                circs = self._pass_manager.run(job_circuits)
+                observables = [
+                    op.apply_layout(circs[i].layout) for i, op in enumerate(job_observables)
+                ]
+            # Prepare circuit-observable-parameter tuples (PUBs)
+            circuit_observable_params = []
+            for pub in zip(circs, observables, job_param_values):
+                circuit_observable_params.append(pub)
+
+            # For BaseEstimatorV2, run the estimator using PUBs and specified precision
+            job = self._estimator.run(circuit_observable_params)
+            try:
+                results = job.result()
+            except Exception as exc:
+                raise AlgorithmError("Estimator job failed.") from exc
+            results = np.array([float(r.data.evs) for r in results])
+            opt = Options(**options)
+            # Compute the gradients.
+            gradients = []
+            partial_sum_n = 0
+            for n in all_n:
+                # this disable is needed as Pylint does not understand derivative_type is a property if
+                # it is only defined in the base class and the getter is in the child
+                # pylint: disable=comparison-with-callable
+                if self.derivative_type == DerivativeType.COMPLEX:
+                    gradient = np.zeros(n // 2, dtype="complex")
+                    gradient.real = results[partial_sum_n : partial_sum_n + n // 2]
+                    gradient.imag = results[partial_sum_n + n // 2 : partial_sum_n + n]
+
+                else:
+                    gradient = np.real(results[partial_sum_n : partial_sum_n + n])
+                partial_sum_n += n
+                gradients.append(gradient)
+
+        else:
+            raise AlgorithmError(
+                "The accepted estimators are BaseEstimatorV1 and BaseEstimatorV2; got "
+                + f"{type(self._estimator)} instead. Note that BaseEstimatorV1 is deprecated in"
+                + "Qiskit and removed in Qiskit IBM Runtime."
+            )
         return EstimatorGradientResult(gradients=gradients, metadata=metadata, options=opt)
diff --git a/qiskit_machine_learning/gradients/lin_comb/lin_comb_qgt.py b/qiskit_machine_learning/gradients/lin_comb/lin_comb_qgt.py
deleted file mode 100644
index afa452ae5..000000000
--- a/qiskit_machine_learning/gradients/lin_comb/lin_comb_qgt.py
+++ /dev/null
@@ -1,258 +0,0 @@
-# This code is part of a Qiskit project.
-#
-# (C) Copyright IBM 2022, 2024.
-#
-# This code is licensed under the Apache License, Version 2.0. You may
-# obtain a copy of this license in the LICENSE.txt file in the root directory
-# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
-#
-# Any modifications or derivative works of this code must retain this
-# copyright notice, and modified files need to carry a notice indicating
-# that they have been altered from the originals.
-"""
-A class for the Linear Combination Quantum Gradient Tensor.
-"""
-
-from __future__ import annotations
-
-from collections.abc import Sequence
-
-import numpy as np
-
-from qiskit.circuit import Parameter, QuantumCircuit
-from qiskit.primitives import BaseEstimator
-from qiskit.primitives.utils import _circuit_key
-from qiskit.providers import Options
-from qiskit.quantum_info import SparsePauliOp
-
-from ..base.base_qgt import BaseQGT
-from .lin_comb_estimator_gradient import LinCombEstimatorGradient
-from ..base.qgt_result import QGTResult
-from ..utils import DerivativeType, _make_lin_comb_qgt_circuit, _make_lin_comb_observables
-
-from ...exceptions import AlgorithmError
-
-
-class LinCombQGT(BaseQGT):
-    """Computes the Quantum Geometric Tensor (QGT) given a pure, parameterized quantum state.
-
-    This method employs a linear combination of unitaries [1].
-
-    **Reference:**
-
-        [1]: Schuld et al., "Evaluating analytic gradients on quantum hardware" (2018).
-             `arXiv:1811.11184 <https://arxiv.org/pdf/1811.11184.pdf>`_
-    """
-
-    SUPPORTED_GATES = [
-        "rx",
-        "ry",
-        "rz",
-        "rzx",
-        "rzz",
-        "ryy",
-        "rxx",
-        "cx",
-        "cy",
-        "cz",
-        "ccx",
-        "swap",
-        "iswap",
-        "h",
-        "t",
-        "s",
-        "sdg",
-        "x",
-        "y",
-        "z",
-    ]
-
-    def __init__(
-        self,
-        estimator: BaseEstimator,
-        phase_fix: bool = True,
-        derivative_type: DerivativeType = DerivativeType.COMPLEX,
-        options: Options | None = None,
-    ):
-        r"""
-        Args:
-            estimator: The estimator used to compute the QGT.
-            phase_fix: Whether to calculate the second term (phase fix) of the QGT, which is
-                :math:`\langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle`.
-                Default to ``True``.
-            derivative_type: The type of derivative. Can be either ``DerivativeType.REAL``
-                ``DerivativeType.IMAG``, or ``DerivativeType.COMPLEX``. Defaults to
-                ``DerivativeType.REAL``.
-
-                - ``DerivativeType.REAL`` computes
-
-                .. math::
-
-                    \mathrm{Re(QGT)}_{ij}= \mathrm{Re}[\langle \partial_i \psi | \partial_j \psi \rangle
-                        - \langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle].
-
-                - ``DerivativeType.IMAG`` computes
-
-                .. math::
-
-                    \mathrm{Re(QGT)}_{ij}= \mathrm{Im}[\langle \partial_i \psi | \partial_j \psi \rangle
-                        - \langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle].
-
-                - ``DerivativeType.COMPLEX`` computes
-
-                .. math::
-
-                    \mathrm{QGT}_{ij}= [\langle \partial_i \psi | \partial_j \psi \rangle
-                        - \langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle].
-
-            options: Backend runtime options used for circuit execution. The order of priority is:
-                options in ``run`` method > QGT's default options > primitive's default
-                setting. Higher priority setting overrides lower priority setting.
-        """
-        super().__init__(estimator, phase_fix, derivative_type, options=options)
-        self._gradient = LinCombEstimatorGradient(
-            estimator, derivative_type=DerivativeType.COMPLEX, options=options
-        )
-        self._lin_comb_qgt_circuit_cache: dict[
-            tuple, dict[tuple[Parameter, Parameter], QuantumCircuit]
-        ] = {}
-
-    def _run(
-        self,
-        circuits: Sequence[QuantumCircuit],
-        parameter_values: Sequence[Sequence[float]],
-        parameters: Sequence[Sequence[Parameter]],
-        **options,
-    ) -> QGTResult:
-        """Compute the QGT on the given circuits."""
-        g_circuits, g_parameter_values, g_parameters = self._preprocess(
-            circuits, parameter_values, parameters, self.SUPPORTED_GATES
-        )
-        results = self._run_unique(g_circuits, g_parameter_values, g_parameters, **options)
-        return self._postprocess(results, circuits, parameter_values, parameters)
-
-    def _run_unique(
-        self,
-        circuits: Sequence[QuantumCircuit],
-        parameter_values: Sequence[Sequence[float]],
-        parameters: Sequence[Sequence[Parameter]],
-        **options,
-    ) -> QGTResult:
-        """Compute the QGTs on the given circuits."""
-        job_circuits, job_observables, job_param_values, metadata = [], [], [], []
-        all_n, all_m = [], []
-        phase_fixes: list[int | np.ndarray] = []
-
-        for circuit, parameter_values_, parameters_ in zip(circuits, parameter_values, parameters):
-            # Prepare circuits for the gradient of the specified parameters.
-            parameters_ = [p for p in circuit.parameters if p in parameters_]
-            meta = {"parameters": parameters_}
-            metadata.append(meta)
-
-            # Compute the first term in the QGT
-            circuit_key = _circuit_key(circuit)
-            if circuit_key not in self._lin_comb_qgt_circuit_cache:
-                # generate the all of the circuits for the first term in the QGT and cache them.
-                # Only the circuit related to specified parameters will be executed.
-                # In the future, we can generate the specified circuits on demand.
-                self._lin_comb_qgt_circuit_cache[circuit_key] = _make_lin_comb_qgt_circuit(circuit)
-            lin_comb_qgt_circuits = self._lin_comb_qgt_circuit_cache[circuit_key]
-
-            qgt_circuits = []
-            rows, cols = np.triu_indices(len(parameters_))
-            for row, col in zip(rows, cols):
-                param_i = parameters_[row]
-                param_j = parameters_[col]
-                qgt_circuits.append(lin_comb_qgt_circuits[(param_i, param_j)])
-
-            observable = SparsePauliOp.from_list([("I" * circuit.num_qubits, 1)])
-            observable_1, observable_2 = _make_lin_comb_observables(
-                observable, self._derivative_type
-            )
-
-            n = len(qgt_circuits)
-            if self._derivative_type == DerivativeType.COMPLEX:
-                job_circuits.extend(qgt_circuits * 2)
-                job_observables.extend([observable_1] * n + [observable_2] * n)
-                job_param_values.extend([parameter_values_] * 2 * n)
-                all_m.append(len(parameters_))
-                all_n.append(2 * n)
-            else:
-                job_circuits.extend(qgt_circuits)
-                job_observables.extend([observable_1] * n)
-                job_param_values.extend([parameter_values_] * n)
-                all_m.append(len(parameters_))
-                all_n.append(n)
-
-        # Run the single job with all circuits.
-        job = self._estimator.run(
-            job_circuits,
-            job_observables,
-            job_param_values,
-            **options,
-        )
-
-        if self._phase_fix:
-            # Compute the second term in the QGT if phase fix is enabled.
-            phase_fix_obs = [
-                SparsePauliOp.from_list([("I" * circuit.num_qubits, 1)]) for circuit in circuits
-            ]
-            phase_fix_job = self._gradient.run(
-                circuits=circuits,
-                observables=phase_fix_obs,
-                parameter_values=parameter_values,
-                parameters=parameters,
-                **options,
-            )
-
-        try:
-            results = job.result()
-            if self._phase_fix:
-                gradient_results = phase_fix_job.result()
-        except AlgorithmError as exc:
-            raise AlgorithmError("Estimator job or gradient job failed.") from exc
-
-        # Compute the phase fix
-        if self._phase_fix:
-            for gradient in gradient_results.gradients:
-                phase_fix = np.outer(np.conjugate(gradient), gradient)
-                # Select the real or imaginary part of the phase fix if needed
-                if self.derivative_type == DerivativeType.REAL:
-                    phase_fix = np.real(phase_fix)
-                elif self.derivative_type == DerivativeType.IMAG:
-                    phase_fix = np.imag(phase_fix)
-                phase_fixes.append(phase_fix)
-        else:
-            phase_fixes = [0 for i in range(len(circuits))]
-        # Compute the QGT
-        qgts = []
-        partial_sum_n = 0
-        for i, (n, m) in enumerate(zip(all_n, all_m)):
-            qgt = np.zeros((m, m), dtype="complex")
-            # Compute the first term in the QGT
-            if self.derivative_type == DerivativeType.COMPLEX:
-                qgt[np.triu_indices(m)] = results.values[partial_sum_n : partial_sum_n + n // 2]
-                qgt[np.triu_indices(m)] += (
-                    1j * results.values[partial_sum_n + n // 2 : partial_sum_n + n]
-                )
-            elif self.derivative_type == DerivativeType.REAL:
-                qgt[np.triu_indices(m)] = results.values[partial_sum_n : partial_sum_n + n]
-            elif self.derivative_type == DerivativeType.IMAG:
-                qgt[np.triu_indices(m)] = 1j * results.values[partial_sum_n : partial_sum_n + n]
-
-            # Add the conjugate of the upper triangle to the lower triangle
-            qgt += np.triu(qgt, k=1).conjugate().T
-            if self.derivative_type == DerivativeType.REAL:
-                qgt = np.real(qgt)
-            elif self.derivative_type == DerivativeType.IMAG:
-                qgt = np.imag(qgt)
-
-            # Subtract the phase fix from the QGT
-            qgt = qgt - phase_fixes[i]
-            partial_sum_n += n
-            qgts.append(qgt / 4)
-
-        opt = self._get_local_options(options)
-        return QGTResult(
-            qgts=qgts, derivative_type=self.derivative_type, metadata=metadata, options=opt
-        )
diff --git a/qiskit_machine_learning/gradients/lin_comb/lin_comb_sampler_gradient.py b/qiskit_machine_learning/gradients/lin_comb/lin_comb_sampler_gradient.py
index cbbe8bf45..96e4a65d5 100644
--- a/qiskit_machine_learning/gradients/lin_comb/lin_comb_sampler_gradient.py
+++ b/qiskit_machine_learning/gradients/lin_comb/lin_comb_sampler_gradient.py
@@ -19,9 +19,13 @@
 from collections.abc import Sequence
 
 from qiskit.circuit import Parameter, QuantumCircuit
-from qiskit.primitives import BaseSampler
 from qiskit.primitives.utils import _circuit_key
+
+from qiskit.primitives import BaseSampler, BaseSamplerV1
+from qiskit.primitives.base import BaseSamplerV2
+from qiskit.result import QuasiDistribution
 from qiskit.providers import Options
+from qiskit.transpiler.passmanager import BasePassManager
 
 from ..base.base_sampler_gradient import BaseSamplerGradient
 from ..base.sampler_gradient_result import SamplerGradientResult
@@ -62,17 +66,24 @@ class LinCombSamplerGradient(BaseSamplerGradient):
         "z",
     ]
 
-    def __init__(self, sampler: BaseSampler, options: Options | None = None):
+    def __init__(
+        self,
+        sampler: BaseSampler,
+        options: Options | None = None,
+        pass_manager: BasePassManager | None = None,
+    ):
         """
         Args:
             sampler: The sampler used to compute the gradients.
             options: Primitive backend runtime options used for circuit execution.
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
-                Higher priority setting overrides lower priority setting
+                Higher priority setting overrides lower priority setting.
+            pass_manager: The pass manager to transpile the circuits if necessary.
+                Defaults to ``None``, as some primitives do not need transpiled circuits.
         """
         self._lin_comb_cache: dict[tuple, dict[Parameter, QuantumCircuit]] = {}
-        super().__init__(sampler, options)
+        super().__init__(sampler, options, pass_manager=pass_manager)
 
     def _run(
         self,
@@ -94,7 +105,7 @@ def _run_unique(
         parameter_values: Sequence[Sequence[float]],
         parameters: Sequence[Sequence[Parameter]],
         **options,
-    ) -> SamplerGradientResult:
+    ) -> SamplerGradientResult:  # pragma: no cover
         """Compute the sampler gradients on the given circuits."""
         job_circuits, job_param_values, metadata = [], [], []
         all_n = []
@@ -119,8 +130,25 @@ def _run_unique(
             job_param_values.extend([parameter_values_] * n)
             all_n.append(n)
 
+        opt = options
         # Run the single job with all circuits.
-        job = self._sampler.run(job_circuits, job_param_values, **options)
+        if isinstance(self._sampler, BaseSamplerV1):
+            job = self._sampler.run(job_circuits, job_param_values, **options)
+            opt = self._get_local_options(options)
+        elif isinstance(self._sampler, BaseSamplerV2):
+            if self._pass_manager is None:
+                circs = job_circuits
+                _len_quasi_dist = 2 ** job_circuits[0].num_qubits
+            else:
+                circs = self._pass_manager.run(job_circuits)
+                _len_quasi_dist = 2 ** circs[0].layout._input_qubit_count
+            circ_params = [(circs[i], job_param_values[i]) for i in range(len(job_param_values))]
+            job = self._sampler.run(circ_params)
+        else:
+            raise AlgorithmError(
+                "The accepted estimators are BaseSamplerV1 (deprecated) and BaseSamplerV2; got "
+                + f"{type(self._sampler)} instead."
+            )
         try:
             results = job.result()
         except Exception as exc:
@@ -131,7 +159,21 @@ def _run_unique(
         partial_sum_n = 0
         for i, n in enumerate(all_n):
             gradient = []
-            result = results.quasi_dists[partial_sum_n : partial_sum_n + n]
+            if isinstance(self._sampler, BaseSamplerV1):
+                result = results.quasi_dists[partial_sum_n : partial_sum_n + n]
+
+            elif isinstance(self._sampler, BaseSamplerV2):
+                result = []
+                for x in range(partial_sum_n, partial_sum_n + n):
+                    bitstring_counts = results[x].data.meas.get_counts()
+
+                    # Normalize the counts to probabilities
+                    total_shots = sum(bitstring_counts.values())
+                    probabilities = {k: v / total_shots for k, v in bitstring_counts.items()}
+
+                    # Convert to quasi-probabilities
+                    counts = QuasiDistribution(probabilities)
+                    result.append({k: v for k, v in counts.items() if int(k) < _len_quasi_dist})
             m = 2 ** circuits[i].num_qubits
             for dist in result:
                 grad_dist: dict[int, float] = defaultdict(float)
@@ -144,5 +186,4 @@ def _run_unique(
             gradients.append(gradient)
             partial_sum_n += n
 
-        opt = self._get_local_options(options)
         return SamplerGradientResult(gradients=gradients, metadata=metadata, options=opt)
diff --git a/qiskit_machine_learning/gradients/param_shift/param_shift_estimator_gradient.py b/qiskit_machine_learning/gradients/param_shift/param_shift_estimator_gradient.py
index 8bbe5f051..fe65c8a70 100644
--- a/qiskit_machine_learning/gradients/param_shift/param_shift_estimator_gradient.py
+++ b/qiskit_machine_learning/gradients/param_shift/param_shift_estimator_gradient.py
@@ -28,7 +28,6 @@
 from ..base.base_estimator_gradient import BaseEstimatorGradient
 from ..base.estimator_gradient_result import EstimatorGradientResult
 from ..utils import _make_param_shift_parameter_values
-from ...exceptions import QiskitMachineLearningError
 
 
 class ParamShiftEstimatorGradient(BaseEstimatorGradient):
@@ -101,6 +100,7 @@ def _run_unique(
             job_param_values.extend(param_shift_parameter_values)
             all_n.append(n)
 
+        opt = Options(**options)
         # Determine how to run the estimator based on its version
         if isinstance(self._estimator, BaseEstimatorV1):
             # Run the single job with all circuits.
@@ -124,13 +124,17 @@ def _run_unique(
             opt = self._get_local_options(options)
 
         elif isinstance(self._estimator, BaseEstimatorV2):
-            isa_g_circs = self._pass_manager.run(job_circuits)
-            isa_g_observables = [
-                op.apply_layout(isa_g_circs[i].layout) for i, op in enumerate(job_observables)
-            ]
+            if self._pass_manager is None:
+                circs_ = job_circuits
+                observables_ = job_observables
+            else:
+                circs_ = self._pass_manager.run(job_circuits)
+                observables_ = [
+                    op.apply_layout(circs_[i].layout) for i, op in enumerate(job_observables)
+                ]
             # Prepare circuit-observable-parameter tuples (PUBs)
             circuit_observable_params = []
-            for pub in zip(isa_g_circs, isa_g_observables, job_param_values):
+            for pub in zip(circs_, observables_, job_param_values):
                 circuit_observable_params.append(pub)
 
             # For BaseEstimatorV2, run the estimator using PUBs and specified precision
@@ -147,13 +151,4 @@ def _run_unique(
                 gradients.append(gradient_)
                 partial_sum_n += n
 
-            opt = Options(**options)
-
-        else:
-            raise QiskitMachineLearningError(
-                "The accepted estimators are BaseEstimatorV1 and BaseEstimatorV2; got "
-                + f"{type(self._estimator)} instead. Note that BaseEstimatorV1 is deprecated in"
-                + "Qiskit and removed in Qiskit IBM Runtime."
-            )
-
         return EstimatorGradientResult(gradients=gradients, metadata=metadata, options=opt)
diff --git a/qiskit_machine_learning/gradients/param_shift/param_shift_sampler_gradient.py b/qiskit_machine_learning/gradients/param_shift/param_shift_sampler_gradient.py
index f327b6453..89efe6ec8 100644
--- a/qiskit_machine_learning/gradients/param_shift/param_shift_sampler_gradient.py
+++ b/qiskit_machine_learning/gradients/param_shift/param_shift_sampler_gradient.py
@@ -27,7 +27,7 @@
 from ..base.base_sampler_gradient import BaseSamplerGradient
 from ..base.sampler_gradient_result import SamplerGradientResult
 from ..utils import _make_param_shift_parameter_values
-from ...exceptions import AlgorithmError, QiskitMachineLearningError
+from ...exceptions import AlgorithmError
 
 
 class ParamShiftSamplerGradient(BaseSamplerGradient):
@@ -78,7 +78,10 @@ def _run_unique(
         parameters: Sequence[Sequence[Parameter]],
         **options,
     ) -> SamplerGradientResult:
-        """Compute the sampler gradients on the given circuits."""
+        """Compute the sampler gradients on the given circuits.
+        Raises:
+            AlgorithmError: If an invalid ``sampler``provided or if sampler job failed.
+        """
         job_circuits, job_param_values, metadata = [], [], []
         all_n = []
         for circuit, parameter_values_, parameters_ in zip(circuits, parameter_values, parameters):
@@ -98,21 +101,18 @@ def _run_unique(
             job = self._sampler.run(job_circuits, job_param_values, **options)
         elif isinstance(self._sampler, BaseSamplerV2):
             if self._pass_manager is None:
-                raise QiskitMachineLearningError(
-                    "To use ParameterShifSamplerGradient with SamplerV2 you "
-                    + "must pass a gradient with a pass manager"
-                )
-            isa_g_circs = self._pass_manager.run(job_circuits)
-            circ_params = [
-                (isa_g_circs[i], job_param_values[i]) for i in range(len(job_param_values))
-            ]
+                _circs = job_circuits
+                _len_quasi_dist = 2 ** job_circuits[0].num_qubits
+            else:
+                _circs = self._pass_manager.run(job_circuits)
+                _len_quasi_dist = 2 ** _circs[0].layout._input_qubit_count
+            circ_params = [(_circs[i], job_param_values[i]) for i in range(len(job_param_values))]
             job = self._sampler.run(circ_params)
         else:
             raise AlgorithmError(
                 "The accepted estimators are BaseSamplerV1 (deprecated) and BaseSamplerV2; got "
                 + f"{type(self._sampler)} instead."
             )
-
         try:
             results = job.result()
         except Exception as exc:
@@ -140,9 +140,7 @@ def _run_unique(
 
                     # Convert to quasi-probabilities
                     counts = QuasiDistribution(probabilities)
-                    result.append(
-                        {k: v for k, v in counts.items() if int(k) < self._len_quasi_dist}
-                    )
+                    result.append({k: v for k, v in counts.items() if int(k) < _len_quasi_dist})
                     opt = options
 
             for dist_plus, dist_minus in zip(result[: n // 2], result[n // 2 :]):
diff --git a/qiskit_machine_learning/gradients/qfi.py b/qiskit_machine_learning/gradients/qfi.py
deleted file mode 100644
index ad0c83e85..000000000
--- a/qiskit_machine_learning/gradients/qfi.py
+++ /dev/null
@@ -1,171 +0,0 @@
-# This code is part of a Qiskit project.
-#
-# (C) Copyright IBM 2022, 2024.
-#
-# This code is licensed under the Apache License, Version 2.0. You may
-# obtain a copy of this license in the LICENSE.txt file in the root directory
-# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
-#
-# Any modifications or derivative works of this code must retain this
-# copyright notice, and modified files need to carry a notice indicating
-# that they have been altered from the originals.
-"""
-A class for the Quantum Fisher Information.
-"""
-
-from __future__ import annotations
-
-from abc import ABC
-from collections.abc import Sequence
-from copy import copy
-
-from qiskit.circuit import Parameter, QuantumCircuit
-from qiskit.providers import Options
-
-from .base.base_qgt import BaseQGT
-from .lin_comb.lin_comb_estimator_gradient import DerivativeType
-from .qfi_result import QFIResult
-
-from ..algorithm_job import AlgorithmJob
-from ..exceptions import AlgorithmError
-
-
-class QFI(ABC):
-    r"""Computes the Quantum Fisher Information (QFI) given a pure,
-    parameterized quantum state. QFI is defined as:
-
-    .. math::
-
-        \mathrm{QFI}_{ij}= 4 \mathrm{Re}[\langle \partial_i \psi | \partial_j \psi \rangle
-        - \langle\partial_i \psi | \psi \rangle \langle\psi | \partial_j \psi \rangle].
-    """
-
-    def __init__(
-        self,
-        qgt: BaseQGT,
-        options: Options | None = None,
-    ):
-        r"""
-        Args:
-            qgt: The quantum geometric tensor used to compute the QFI.
-            options: Backend runtime options used for circuit execution. The order of priority is:
-                options in ``run`` method > QFI's default options > primitive's default
-                setting. Higher priority setting overrides lower priority setting.
-        """
-        self._qgt: BaseQGT = qgt
-        self._default_options = Options()
-        if options is not None:
-            self._default_options.update_options(**options)
-
-    def run(
-        self,
-        circuits: Sequence[QuantumCircuit],
-        parameter_values: Sequence[Sequence[float]],
-        parameters: Sequence[Sequence[Parameter] | None] | None = None,
-        **options,
-    ) -> AlgorithmJob:
-        """Run the job of the QFIs on the given circuits.
-
-        Args:
-            circuits: The list of quantum circuits to compute the QFIs.
-            parameter_values: The list of parameter values to be bound to the circuit.
-            parameters: The sequence of parameters to calculate only the QFIs of
-                the specified parameters. Each sequence of parameters corresponds to a circuit in
-                ``circuits``. Defaults to None, which means that the QFIs of all parameters in
-                each circuit are calculated.
-            options: Primitive backend runtime options used for circuit execution.
-                The order of priority is: options in ``run`` method > QFI's
-                default options > QGT's default setting.
-                Higher priority setting overrides lower priority setting.
-
-        Returns:
-            The job object of the QFIs of the expectation values. The i-th result corresponds to
-            ``circuits[i]`` evaluated with parameters bound as ``parameter_values[i]``.
-        """
-
-        if isinstance(circuits, QuantumCircuit):
-            # Allow a single circuit to be passed in.
-            circuits = (circuits,)
-
-        if parameters is None:
-            # If parameters is None, we calculate the gradients of all parameters in each circuit.
-            parameters = [circuit.parameters for circuit in circuits]
-        else:
-            # If parameters is not None, we calculate the gradients of the specified parameters.
-            # None in parameters means that the gradients of all parameters in the corresponding
-            # circuit are calculated.
-            parameters = [
-                params if params is not None else circuits[i].parameters
-                for i, params in enumerate(parameters)
-            ]
-        # The priority of run option is as follows:
-        # options in ``run`` method > QFI's default options > QGT's default setting.
-        opts = copy(self._default_options)
-        opts.update_options(**options)
-        job = AlgorithmJob(self._run, circuits, parameter_values, parameters, **opts.__dict__)
-        job.submit()
-        return job
-
-    def _run(
-        self,
-        circuits: Sequence[QuantumCircuit],
-        parameter_values: Sequence[Sequence[float]],
-        parameters: Sequence[Sequence[Parameter]],
-        **options,
-    ) -> QFIResult:
-        """Compute the QFI on the given circuits."""
-        # Set the derivative type to real
-        temp_derivative_type, self._qgt.derivative_type = (
-            self._qgt.derivative_type,
-            DerivativeType.REAL,
-        )
-        job = self._qgt.run(circuits, parameter_values, parameters, **options)
-
-        try:
-            result = job.result()
-        except AlgorithmError as exc:
-            raise AlgorithmError("Estimator job or gradient job failed.") from exc
-
-        self._qgt.derivative_type = temp_derivative_type
-
-        return QFIResult(
-            qfis=[4 * qgt.real for qgt in result.qgts],
-            metadata=result.metadata,
-            options=result.options,
-        )
-
-    @property
-    def options(self) -> Options:
-        """Return the union of QGT's options setting and QFI's default options,
-        where, if the same field is set in both, the QFI's default options override
-        the QGT's default setting.
-
-        Returns:
-            The QFI default + QGT options.
-        """
-        return self._get_local_options(self._default_options.__dict__)
-
-    def update_default_options(self, **options):
-        """Update the gradient's default options setting.
-
-        Args:
-            **options: The fields to update the default options.
-        """
-
-        self._default_options.update_options(**options)
-
-    def _get_local_options(self, options: Options) -> Options:
-        """Return the union of the QFI default setting,
-        the QGT default options, and the options in the ``run`` method.
-        The order of priority is: options in ``run`` method > QFI's default options > QGT's
-        default setting.
-
-        Args:
-            options: The fields to update the options
-
-        Returns:
-            The QFI default + QGT default + run options.
-        """
-        opts = copy(self._qgt.options)
-        opts.update_options(**options)
-        return opts
diff --git a/qiskit_machine_learning/gradients/qfi_result.py b/qiskit_machine_learning/gradients/qfi_result.py
deleted file mode 100644
index 57aeeb932..000000000
--- a/qiskit_machine_learning/gradients/qfi_result.py
+++ /dev/null
@@ -1,35 +0,0 @@
-# This code is part of a Qiskit project.
-#
-# (C) Copyright IBM 2022, 2024.
-#
-# This code is licensed under the Apache License, Version 2.0. You may
-# obtain a copy of this license in the LICENSE.txt file in the root directory
-# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
-#
-# Any modifications or derivative works of this code must retain this
-# copyright notice, and modified files need to carry a notice indicating
-# that they have been altered from the originals.
-"""
-QFI result class
-"""
-
-from __future__ import annotations
-
-from dataclasses import dataclass
-from typing import Any
-
-import numpy as np
-
-from qiskit.providers import Options
-
-
-@dataclass(frozen=True)
-class QFIResult:
-    """Result of QFI."""
-
-    qfis: list[np.ndarray]
-    """The QFI."""
-    metadata: list[dict[str, Any]]
-    """Additional information about the job."""
-    options: Options
-    """Primitive runtime options for the execution of the job."""
diff --git a/qiskit_machine_learning/gradients/spsa/spsa_estimator_gradient.py b/qiskit_machine_learning/gradients/spsa/spsa_estimator_gradient.py
index 8f524a0bf..801e48182 100644
--- a/qiskit_machine_learning/gradients/spsa/spsa_estimator_gradient.py
+++ b/qiskit_machine_learning/gradients/spsa/spsa_estimator_gradient.py
@@ -19,9 +19,11 @@
 import numpy as np
 
 from qiskit.circuit import Parameter, QuantumCircuit
-from qiskit.primitives import BaseEstimator
 from qiskit.providers import Options
 from qiskit.quantum_info.operators.base_operator import BaseOperator
+from qiskit.primitives.base import BaseEstimatorV2
+from qiskit.primitives import BaseEstimator, BaseEstimatorV1
+from qiskit.transpiler.passmanager import BasePassManager
 
 from ..base.base_estimator_gradient import BaseEstimatorGradient
 from ..base.estimator_gradient_result import EstimatorGradientResult
@@ -44,10 +46,11 @@ class SPSAEstimatorGradient(BaseEstimatorGradient):
     def __init__(
         self,
         estimator: BaseEstimator,
-        epsilon: float,
+        epsilon: float = 1e-6,
         batch_size: int = 1,
         seed: int | None = None,
         options: Options | None = None,
+        pass_manager: BasePassManager | None = None,
     ):
         """
         Args:
@@ -58,7 +61,9 @@ def __init__(
             options: Primitive backend runtime options used for circuit execution.
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
-                Higher priority setting overrides lower priority setting
+                Higher priority setting overrides lower priority setting.
+            pass_manager: The pass manager to transpile the circuits if necessary.
+                Defaults to ``None``, as some primitives do not need transpiled circuits.
 
         Raises:
             ValueError: If ``epsilon`` is not positive.
@@ -69,7 +74,7 @@ def __init__(
         self._batch_size = batch_size
         self._seed = np.random.default_rng(seed)
 
-        super().__init__(estimator, options)
+        super().__init__(estimator, options=options, pass_manager=pass_manager)
 
     def _run(
         self,
@@ -78,7 +83,7 @@ def _run(
         parameter_values: Sequence[Sequence[float]] | np.ndarray,
         parameters: Sequence[Sequence[Parameter]],
         **options,
-    ) -> EstimatorGradientResult:
+    ) -> EstimatorGradientResult:  # pragma: no cover
         """Compute the estimator gradients on the given circuits."""
         job_circuits, job_observables, job_param_values, metadata, offsets = [], [], [], [], []
         all_n = []
@@ -102,34 +107,84 @@ def _run(
             job_observables.extend([observable] * 2 * self._batch_size)
             job_param_values.extend(plus + minus)
             all_n.append(2 * self._batch_size)
+        if isinstance(self._estimator, BaseEstimatorV1):
+            # Run the single job with all circuits.
+            job = self._estimator.run(
+                job_circuits,
+                job_observables,
+                job_param_values,
+                **options,
+            )
+            try:
+                results = job.result()
+            except Exception as exc:
+                raise AlgorithmError("Estimator job failed.") from exc
+            # Compute the gradients.
+            gradients = []
+            partial_sum_n = 0
+            for i, n in enumerate(all_n):
+                result = results.values[partial_sum_n : partial_sum_n + n]
+                partial_sum_n += n
+                n = len(result) // 2
+                diffs = (result[:n] - result[n:]) / (2 * self._epsilon)
+                # Calculate the gradient for each batch.
+                # Note that (``diff`` / ``offset``) is the gradient
+                # since ``offset`` is a perturbation vector of 1s and -1s.
+                batch_gradients = np.array(
+                    [diff / offset for diff, offset in zip(diffs, offsets[i])]
+                )
+                # Take the average of the batch gradients.
+                gradient = np.mean(batch_gradients, axis=0)
+                indices = [circuits[i].parameters.data.index(p) for p in metadata[i]["parameters"]]
+                gradients.append(gradient[indices])
+            opt = self._get_local_options(options)
+        elif isinstance(self._estimator, BaseEstimatorV2):
+            if self._pass_manager is None:
+                circs = job_circuits
+                observables = job_observables
+            else:
+                circs = self._pass_manager.run(job_circuits)
+                observables = [
+                    op.apply_layout(circs[x].layout) for x, op in enumerate(job_observables)
+                ]
+            # Prepare circuit-observable-parameter tuples (PUBs)
+            circuit_observable_params = []
+            for pub in zip(circs, observables, job_param_values):
+                circuit_observable_params.append(pub)
+
+            # For BaseEstimatorV2, run the estimator using PUBs and specified precision
+            job = self._estimator.run(circuit_observable_params)
+            try:
+                results = job.result()
+            except Exception as exc:
+                raise AlgorithmError("Estimator job failed.") from exc
+            results = np.array([float(r.data.evs) for r in results])
+            opt = Options(**options)
+
+            # Compute the gradients.
+            gradients = []
+            partial_sum_n = 0
+            for i, n in enumerate(all_n):
+                result = results[partial_sum_n : partial_sum_n + n]
+                partial_sum_n += n
+                n = len(result) // 2
+                diffs = (result[:n] - result[n:]) / (2 * self._epsilon)
+                # Calculate the gradient for each batch.
+                # Note that (``diff`` / ``offset``) is the gradient
+                # since ``offset`` is a perturbation vector of 1s and -1s.
+                batch_gradients = np.array(
+                    [diff / offset for diff, offset in zip(diffs, offsets[i])]
+                )
+                # Take the average of the batch gradients.
+                gradient = np.mean(batch_gradients, axis=0)
+                indices = [circuits[i].parameters.data.index(p) for p in metadata[i]["parameters"]]
+                gradients.append(gradient[indices])
+
+        else:
+            raise AlgorithmError(
+                "The accepted estimators are BaseEstimatorV1 and BaseEstimatorV2; got "
+                + f"{type(self._estimator)} instead. Note that BaseEstimatorV1 is deprecated in"
+                + "Qiskit and removed in Qiskit IBM Runtime."
+            )
 
-        # Run the single job with all circuits.
-        job = self._estimator.run(
-            job_circuits,
-            job_observables,
-            job_param_values,
-            **options,
-        )
-        try:
-            results = job.result()
-        except Exception as exc:
-            raise AlgorithmError("Estimator job failed.") from exc
-
-        # Compute the gradients.
-        gradients = []
-        partial_sum_n = 0
-        for i, n in enumerate(all_n):
-            result = results.values[partial_sum_n : partial_sum_n + n]
-            partial_sum_n += n
-            n = len(result) // 2
-            diffs = (result[:n] - result[n:]) / (2 * self._epsilon)
-            # Calculate the gradient for each batch. Note that (``diff`` / ``offset``) is the gradient
-            # since ``offset`` is a perturbation vector of 1s and -1s.
-            batch_gradients = np.array([diff / offset for diff, offset in zip(diffs, offsets[i])])
-            # Take the average of the batch gradients.
-            gradient = np.mean(batch_gradients, axis=0)
-            indices = [circuits[i].parameters.data.index(p) for p in metadata[i]["parameters"]]
-            gradients.append(gradient[indices])
-
-        opt = self._get_local_options(options)
         return EstimatorGradientResult(gradients=gradients, metadata=metadata, options=opt)
diff --git a/qiskit_machine_learning/gradients/spsa/spsa_sampler_gradient.py b/qiskit_machine_learning/gradients/spsa/spsa_sampler_gradient.py
index 1c25b8aaa..922e3d68c 100644
--- a/qiskit_machine_learning/gradients/spsa/spsa_sampler_gradient.py
+++ b/qiskit_machine_learning/gradients/spsa/spsa_sampler_gradient.py
@@ -20,8 +20,12 @@
 import numpy as np
 
 from qiskit.circuit import Parameter, QuantumCircuit
-from qiskit.primitives import BaseSampler
+
+from qiskit.primitives import BaseSampler, BaseSamplerV1
+from qiskit.primitives.base import BaseSamplerV2
+from qiskit.result import QuasiDistribution
 from qiskit.providers import Options
+from qiskit.transpiler.passmanager import BasePassManager
 
 from ..base.base_sampler_gradient import BaseSamplerGradient
 from ..base.sampler_gradient_result import SamplerGradientResult
@@ -44,10 +48,11 @@ class SPSASamplerGradient(BaseSamplerGradient):
     def __init__(
         self,
         sampler: BaseSampler,
-        epsilon: float,
+        epsilon: float = 1e-6,
         batch_size: int = 1,
         seed: int | None = None,
         options: Options | None = None,
+        pass_manager: BasePassManager | None = None,
     ):
         """
         Args:
@@ -59,6 +64,8 @@ def __init__(
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
                 Higher priority setting overrides lower priority setting
+            pass_manager: The pass manager to transpile the circuits if necessary.
+                Defaults to ``None``, as some primitives do not need transpiled circuits.
 
         Raises:
             ValueError: If ``epsilon`` is not positive.
@@ -69,7 +76,7 @@ def __init__(
         self._epsilon = epsilon
         self._seed = np.random.default_rng(seed)
 
-        super().__init__(sampler, options)
+        super().__init__(sampler, options, pass_manager=pass_manager)
 
     def _run(
         self,
@@ -77,7 +84,7 @@ def _run(
         parameter_values: Sequence[Sequence[float]],
         parameters: Sequence[Sequence[Parameter]],
         **options,
-    ) -> SamplerGradientResult:
+    ) -> SamplerGradientResult:  # pragma: no cover
         """Compute the sampler gradients on the given circuits."""
         job_circuits, job_param_values, metadata, offsets = [], [], [], []
         all_n = []
@@ -101,8 +108,25 @@ def _run(
             job_param_values.extend(plus + minus)
             all_n.append(n)
 
+        opt = options
         # Run the single job with all circuits.
-        job = self._sampler.run(job_circuits, job_param_values, **options)
+        if isinstance(self._sampler, BaseSamplerV1):
+            job = self._sampler.run(job_circuits, job_param_values, **options)
+            opt = self._get_local_options(options)
+        elif isinstance(self._sampler, BaseSamplerV2):
+            if self._pass_manager is None:
+                _circs = job_circuits
+                _len_quasi_dist = 2 ** job_circuits[0].num_qubits
+            else:
+                _circs = self._pass_manager.run(job_circuits)
+                _len_quasi_dist = 2 ** _circs[0].layout._input_qubit_count
+            _circ_params = [(_circs[i], job_param_values[i]) for i in range(len(job_param_values))]
+            job = self._sampler.run(_circ_params)
+        else:
+            raise AlgorithmError(
+                "The accepted estimators are BaseSamplerV1 (deprecated) and BaseSamplerV2; got "
+                + f"{type(self._sampler)} instead."
+            )
         try:
             results = job.result()
         except Exception as exc:
@@ -110,10 +134,24 @@ def _run(
 
         # Compute the gradients.
         gradients = []
+        result = []
         partial_sum_n = 0
         for i, n in enumerate(all_n):
             dist_diffs = {}
-            result = results.quasi_dists[partial_sum_n : partial_sum_n + n]
+            if isinstance(self._sampler, BaseSamplerV1):
+                result = results.quasi_dists[partial_sum_n : partial_sum_n + n]
+            elif isinstance(self._sampler, BaseSamplerV2):
+                _result = []
+                for m in range(partial_sum_n, partial_sum_n + n):
+                    _bitstring_counts = results[m].data.meas.get_counts()
+                    # Normalize the counts to probabilities
+                    _total_shots = sum(_bitstring_counts.values())
+                    _probabilities = {k: v / _total_shots for k, v in _bitstring_counts.items()}
+                    # Convert to quasi-probabilities
+                    _counts = QuasiDistribution(_probabilities)
+                    _result.append({k: v for k, v in _counts.items() if int(k) < _len_quasi_dist})
+                    result = [{key: d[key] for key in sorted(d)} for d in _result]
+
             for j, (dist_plus, dist_minus) in enumerate(zip(result[: n // 2], result[n // 2 :])):
                 dist_diff: dict[int, float] = defaultdict(float)
                 for key, value in dist_plus.items():
@@ -121,6 +159,7 @@ def _run(
                 for key, value in dist_minus.items():
                     dist_diff[key] -= value / (2 * self._epsilon)
                 dist_diffs[j] = dist_diff
+
             gradient = []
             indices = [circuits[i].parameters.data.index(p) for p in metadata[i]["parameters"]]
             for j in indices:
@@ -133,5 +172,4 @@ def _run(
             gradients.append(gradient)
             partial_sum_n += n
 
-        opt = self._get_local_options(options)
         return SamplerGradientResult(gradients=gradients, metadata=metadata, options=opt)
diff --git a/qiskit_machine_learning/gradients/utils.py b/qiskit_machine_learning/gradients/utils.py
index 53ef7fcc2..572815a54 100644
--- a/qiskit_machine_learning/gradients/utils.py
+++ b/qiskit_machine_learning/gradients/utils.py
@@ -43,7 +43,6 @@
     RZGate,
     RZXGate,
     RZZGate,
-    XGate,
 )
 from qiskit.quantum_info import SparsePauliOp
 
@@ -109,7 +108,7 @@ def _make_param_shift_parameter_values(  # pylint: disable=invalid-name
 
 
 ################################################################################
-## Linear combination gradient and Linear combination QGT
+## Linear combination gradient
 ################################################################################
 def _make_lin_comb_gradient_circuit(
     circuit: QuantumCircuit, add_measurement: bool = False
@@ -178,78 +177,6 @@ def _gate_gradient(gate: Gate) -> Instruction:
     raise TypeError(f"Unrecognized parameterized gate, {gate}")
 
 
-def _make_lin_comb_qgt_circuit(
-    circuit: QuantumCircuit, add_measurement: bool = False
-) -> dict[tuple[Parameter, Parameter], QuantumCircuit]:
-    """Makes a circuit that computes the linear combination of the QGT circuits."""
-    circuit_temp = circuit.copy()
-    qr_aux = QuantumRegister(1, "aux")
-    circuit_temp.add_register(qr_aux)
-    if add_measurement:
-        cr_aux = ClassicalRegister(1, "aux")
-        circuit_temp.add_bits(cr_aux)
-    circuit_temp.h(qr_aux)
-    circuit_temp.data.insert(0, circuit_temp.data.pop())
-
-    lin_comb_qgt_circuits = {}
-    for i, instruction_i in enumerate(circuit_temp.data):
-        if not instruction_i.operation.is_parameterized():
-            continue
-        for j, instruction_j in enumerate(circuit_temp.data):
-            if not instruction_j.operation.is_parameterized():
-                continue
-            # Calculate the QGT of the i-th gate with respect to the j-th gate.
-            param_i = instruction_i.operation.params[0]
-            param_j = instruction_j.operation.params[0]
-
-            for p_i in param_i.parameters:
-                for p_j in param_j.parameters:
-                    if circuit_temp.parameters.data.index(p_i) > circuit_temp.parameters.data.index(
-                        p_j
-                    ):
-                        continue
-                    gate_i = _gate_gradient(instruction_i.operation)
-                    gate_j = _gate_gradient(instruction_j.operation)
-                    lin_comb_qgt_circuit = circuit_temp.copy()
-                    if i < j:
-                        # insert gate_j to j-th position
-                        lin_comb_qgt_circuit.append(
-                            gate_j, [qr_aux[0]] + list(instruction_j.qubits), []
-                        )
-                        lin_comb_qgt_circuit.data.insert(j, lin_comb_qgt_circuit.data.pop())
-                        # insert gate_i to i-th position with two X gates at its sides
-                        lin_comb_qgt_circuit.append(XGate(), [qr_aux[0]], [])
-                        lin_comb_qgt_circuit.data.insert(i, lin_comb_qgt_circuit.data.pop())
-                        lin_comb_qgt_circuit.append(
-                            gate_i, [qr_aux[0]] + list(instruction_i.qubits), []
-                        )
-                        lin_comb_qgt_circuit.data.insert(i, lin_comb_qgt_circuit.data.pop())
-                        lin_comb_qgt_circuit.append(XGate(), [qr_aux[0]], [])
-                        lin_comb_qgt_circuit.data.insert(i, lin_comb_qgt_circuit.data.pop())
-                    else:
-                        # insert gate_i to i-th position
-                        lin_comb_qgt_circuit.append(
-                            gate_i, [qr_aux[0]] + list(instruction_i.qubits), []
-                        )
-                        lin_comb_qgt_circuit.data.insert(i, lin_comb_qgt_circuit.data.pop())
-                        # insert gate_j to j-th position with two X gates at its sides
-                        lin_comb_qgt_circuit.append(XGate(), [qr_aux[0]], [])
-                        lin_comb_qgt_circuit.data.insert(j, lin_comb_qgt_circuit.data.pop())
-                        lin_comb_qgt_circuit.append(
-                            gate_j, [qr_aux[0]] + list(instruction_j.qubits), []
-                        )
-                        lin_comb_qgt_circuit.data.insert(j, lin_comb_qgt_circuit.data.pop())
-                        lin_comb_qgt_circuit.append(XGate(), [qr_aux[0]], [])
-                        lin_comb_qgt_circuit.data.insert(j, lin_comb_qgt_circuit.data.pop())
-
-                    lin_comb_qgt_circuit.h(qr_aux)
-                    if add_measurement:
-                        lin_comb_qgt_circuit.measure(qr_aux, cr_aux)
-                    lin_comb_qgt_circuits[(p_i, p_j)] = lin_comb_qgt_circuit
-
-    return lin_comb_qgt_circuits
-
-
 def _make_lin_comb_observables(
     observable: SparsePauliOp,
     derivative_type: DerivativeType,
diff --git a/qiskit_machine_learning/neural_networks/estimator_qnn.py b/qiskit_machine_learning/neural_networks/estimator_qnn.py
index 36b3d92ce..ca0da35cc 100644
--- a/qiskit_machine_learning/neural_networks/estimator_qnn.py
+++ b/qiskit_machine_learning/neural_networks/estimator_qnn.py
@@ -15,7 +15,6 @@
 from __future__ import annotations
 
 import logging
-import warnings
 from copy import copy
 from typing import Sequence
 import numpy as np
@@ -25,6 +24,7 @@
 from qiskit.primitives import BaseEstimator, BaseEstimatorV1, Estimator, EstimatorResult
 from qiskit.quantum_info import SparsePauliOp
 from qiskit.quantum_info.operators.base_operator import BaseOperator
+from qiskit.transpiler.passmanager import BasePassManager
 
 
 from ..gradients import (
@@ -115,8 +115,8 @@ def __init__(
         weight_params: Sequence[Parameter] | None = None,
         gradient: BaseEstimatorGradient | None = None,
         input_gradients: bool = False,
-        num_virtual_qubits: int | None = None,
         default_precision: float = 0.015625,
+        pass_manager: BasePassManager | None = None,
     ):
         r"""
         Args:
@@ -147,11 +147,12 @@ def __init__(
                 Note that this parameter is ``False`` by default, and must be explicitly set to
                 ``True`` for a proper gradient computation when using
                 :class:`~qiskit_machine_learning.connectors.TorchConnector`.
-            num_virtual_qubits: Number of virtual qubits.
             default_precision: The default precision for the estimator if not specified during run.
-
+            pass_manager: The pass manager to transpile the circuits, if necessary.
+                Defaults to ``None``, as some primitives do not need transpiled circuits.
         Raises:
             QiskitMachineLearningError: Invalid parameter values.
+            QiskitMachineLearningError: Gradient is required if
         """
         if estimator is None:
             estimator = Estimator()
@@ -164,19 +165,17 @@ def __init__(
                 period="4 months",
             )
         self.estimator = estimator
-        self._org_circuit = circuit
 
-        if num_virtual_qubits is None:
-            self.num_virtual_qubits = circuit.num_qubits
-            warnings.warn(
-                f"No number of qubits was not specified ({num_virtual_qubits}) and was retrieved from "
-                + f"`circuit` ({self.num_virtual_qubits:d}). If `circuit` is transpiled, this may cause "
-                + "unstable behaviour.",
-                UserWarning,
-                stacklevel=2,
-            )
+        if hasattr(circuit.layout, "_input_qubit_count"):
+            self.num_virtual_qubits = circuit.layout._input_qubit_count
         else:
-            self.num_virtual_qubits = num_virtual_qubits
+            if pass_manager is None:
+                self.num_virtual_qubits = circuit.num_qubits
+            else:
+                circuit = pass_manager.run(circuit)
+                self.num_virtual_qubits = circuit.layout._input_qubit_count
+
+        self._org_circuit = circuit
 
         if observables is None:
             observables = SparsePauliOp.from_sparse_list(
@@ -196,14 +195,18 @@ def __init__(
             self._input_params = list(input_params) if input_params is not None else []
             self._weight_params = list(weight_params) if weight_params is not None else []
 
+        # set gradient
         if gradient is None:
-            if isinstance(self.estimator, BaseEstimatorV2):
-                raise QiskitMachineLearningError(
-                    "Please provide a gradient with pass manager initialised."
+            if isinstance(estimator, BaseEstimatorV1):
+                gradient = ParamShiftEstimatorGradient(estimator=self.estimator)
+            else:
+                logger.warning(
+                    "No gradient function provided, creating a gradient function."
+                    " If your Estimator requires transpilation, please provide a pass manager."
+                )
+                gradient = ParamShiftEstimatorGradient(
+                    estimator=self.estimator, pass_manager=pass_manager
                 )
-
-            gradient = ParamShiftEstimatorGradient(self.estimator)
-
         self._default_precision = default_precision
         self.gradient = gradient
         self._input_gradients = input_gradients
diff --git a/qiskit_machine_learning/neural_networks/sampler_qnn.py b/qiskit_machine_learning/neural_networks/sampler_qnn.py
index bb5ca4023..0596d6030 100644
--- a/qiskit_machine_learning/neural_networks/sampler_qnn.py
+++ b/qiskit_machine_learning/neural_networks/sampler_qnn.py
@@ -24,6 +24,7 @@
 from qiskit.circuit import Parameter, QuantumCircuit
 from qiskit.primitives import BaseSampler, SamplerResult, Sampler
 from qiskit.result import QuasiDistribution
+from qiskit.transpiler.passmanager import BasePassManager
 
 import qiskit_machine_learning.optionals as _optionals
 
@@ -132,7 +133,6 @@ def __init__(
         self,
         *,
         circuit: QuantumCircuit,
-        num_virtual_qubits: int | None = None,
         sampler: BaseSampler | None = None,
         input_params: Sequence[Parameter] | None = None,
         weight_params: Sequence[Parameter] | None = None,
@@ -141,6 +141,7 @@ def __init__(
         output_shape: int | tuple[int, ...] | None = None,
         gradient: BaseSamplerGradient | None = None,
         input_gradients: bool = False,
+        pass_manager: BasePassManager | None = None,
     ):
         """
         Args: sampler: The sampler primitive used to compute the neural network's results. If
@@ -170,7 +171,10 @@ def __init__(
         input_gradients: Determines whether to compute gradients with respect to input data. Note
         that this parameter is ``False`` by default, and must be explicitly set to ``True`` for a
         proper gradient computation when using
-        :class:`~qiskit_machine_learning.connectors.TorchConnector`. Raises:
+        :class:`~qiskit_machine_learning.connectors.TorchConnector`.
+        pass_manager: The pass manager to transpile the circuits, if necessary.
+        Defaults to ``None``, as some primitives do not need transpiled circuits.
+        Raises:
         QiskitMachineLearningError: Invalid parameter values.
         """
         # set primitive, provide default
@@ -185,11 +189,14 @@ def __init__(
                 period="4 months",
             )
         self.sampler = sampler
-
-        if num_virtual_qubits is None:
-            # print statement
-            num_virtual_qubits = circuit.num_qubits
-        self.num_virtual_qubits = num_virtual_qubits
+        if hasattr(circuit.layout, "_input_qubit_count"):
+            self.num_virtual_qubits = circuit.layout._input_qubit_count
+        else:
+            if pass_manager is None:
+                self.num_virtual_qubits = circuit.num_qubits
+            else:
+                circuit = pass_manager.run(circuit)
+                self.num_virtual_qubits = circuit.layout._input_qubit_count
 
         self._org_circuit = circuit
 
@@ -204,10 +211,18 @@ def __init__(
             _optionals.HAS_SPARSE.require_now("DOK")
 
         self.set_interpret(interpret, output_shape)
-
         # set gradient
         if gradient is None:
-            gradient = ParamShiftSamplerGradient(sampler=self.sampler)
+            if isinstance(sampler, BaseSamplerV1):
+                gradient = ParamShiftSamplerGradient(sampler=self.sampler)
+            else:
+                logger.warning(
+                    "No gradient function provided, creating a gradient function."
+                    " If your Sampler requires transpilation, please provide a pass manager."
+                )
+                gradient = ParamShiftSamplerGradient(
+                    sampler=self.sampler, pass_manager=pass_manager
+                )
         self.gradient = gradient
 
         self._input_gradients = input_gradients
@@ -270,7 +285,11 @@ def _compute_output_shape(
         interpret: Callable[[int], int | tuple[int, ...]] | None = None,
         output_shape: int | tuple[int, ...] | None = None,
     ) -> tuple[int, ...]:
-        """Validate and compute the output shape."""
+        """Validate and compute the output shape.
+        Raises:
+            QiskitMachineLearningError: If no output shape is given.
+            QiskitMachineLearningError: If an invalid ``sampler``provided.
+        """
 
         # this definition is required by mypy
         output_shape_: tuple[int, ...] = (-1,)
diff --git a/qiskit_machine_learning/state_fidelities/compute_uncompute.py b/qiskit_machine_learning/state_fidelities/compute_uncompute.py
index 03a9d7354..dd96f9731 100644
--- a/qiskit_machine_learning/state_fidelities/compute_uncompute.py
+++ b/qiskit_machine_learning/state_fidelities/compute_uncompute.py
@@ -18,7 +18,7 @@
 from copy import copy
 
 from qiskit import QuantumCircuit
-from qiskit.primitives import BaseSampler, BaseSamplerV1, SamplerResult, StatevectorSampler
+from qiskit.primitives import BaseSampler, BaseSamplerV1, SamplerResult
 from qiskit.primitives.base import BaseSamplerV2
 from qiskit.transpiler.passmanager import PassManager
 from qiskit.result import QuasiDistribution
@@ -59,10 +59,9 @@ def __init__(
         self,
         sampler: BaseSampler | BaseSamplerV2,
         *,
-        num_virtual_qubits: int | None = None,
-        pass_manager: PassManager | None = None,
         options: Options | None = None,
         local: bool = False,
+        pass_manager: PassManager | None = None,
     ) -> None:
         r"""
         Args:
@@ -82,7 +81,8 @@ def __init__(
                 This coincides with the standard (global) fidelity in the limit of
                 the fidelity approaching 1. Might be used to increase the variance
                 to improve trainability in algorithms such as :class:`~.time_evolvers.PVQD`.
-
+            pass_manager: The pass manager to transpile the circuits, if necessary.
+                Defaults to ``None``, as some primitives do not need transpiled circuits.
         Raises:
             ValueError: If the sampler is not an instance of ``BaseSampler``.
         """
@@ -91,16 +91,7 @@ def __init__(
                 f"The sampler should be an instance of BaseSampler or BaseSamplerV2, "
                 f"but got {type(sampler)}"
             )
-        if (
-            isinstance(sampler, BaseSamplerV2)
-            and (pass_manager is None)
-            and not isinstance(sampler, StatevectorSampler)
-        ):
-            raise ValueError(f"A pass_manager should be provided for {type(sampler)}.")
-        if (pass_manager is not None) and (num_virtual_qubits is None):
-            raise ValueError(
-                f"Number of virtual qubits should be provided for {type(pass_manager)}."
-            )
+
         if isinstance(sampler, BaseSamplerV1):
             issue_deprecation_msg(
                 msg="V1 Primitives are deprecated",
@@ -109,8 +100,7 @@ def __init__(
                 period="4 months",
             )
         self._sampler: BaseSampler = sampler
-        self.num_virtual_qubits = num_virtual_qubits
-        self.pass_manager = pass_manager
+        self._pass_manager = pass_manager
         self._local = local
         self._default_options = Options()
         if options is not None:
@@ -138,8 +128,8 @@ def create_fidelity_circuit(
 
         circuit = circuit_1.compose(circuit_2.inverse())
         circuit.measure_all()
-        if self.pass_manager is not None:
-            circuit = self.pass_manager.run(circuit)
+        if self._pass_manager is not None:
+            circuit = self._pass_manager.run(circuit)
         return circuit
 
     def _run(
@@ -171,8 +161,9 @@ def _run(
         Raises:
             ValueError: At least one pair of circuits must be defined.
             AlgorithmError: If the sampler job is not completed successfully.
+            QiskitMachineLearningError: If the sampler is not an instance
+                of ``BaseSamplerV1`` or ``BaseSamplerV2``.
         """
-
         circuits = self._construct_circuits(circuits_1, circuits_2)
         if len(circuits) == 0:
             raise ValueError(
@@ -190,11 +181,13 @@ def _run(
             sampler_job = self._sampler.run(
                 circuits=circuits, parameter_values=values, **opts.__dict__
             )
+            _len_quasi_dist = circuits[0].num_qubits
             local_opts = self._get_local_options(opts.__dict__)
         elif isinstance(self._sampler, BaseSamplerV2):
             sampler_job = self._sampler.run(
                 [(circuits[i], values[i]) for i in range(len(circuits))], **opts.__dict__
             )
+            _len_quasi_dist = circuits[0].layout._input_qubit_count
             local_opts = opts.__dict__
         else:
             raise QiskitMachineLearningError(
@@ -209,7 +202,7 @@ def _run(
             local_opts,
             self._sampler,
             self._post_process_v2,
-            self.num_virtual_qubits,
+            _len_quasi_dist,
         )
 
     @staticmethod
@@ -230,7 +223,7 @@ def _call(
         if isinstance(_sampler, BaseSamplerV1):
             quasi_dists = result.quasi_dists
         elif isinstance(_sampler, BaseSamplerV2):
-            quasi_dists = _post_process_v2(result)
+            quasi_dists = _post_process_v2(result, num_virtual_qubits)
 
         if local:
             raw_fidelities = [
@@ -293,7 +286,7 @@ def _get_local_options(self, options: Options) -> Options:
         opts.update_options(**options)
         return opts
 
-    def _post_process_v2(self, result: SamplerResult):
+    def _post_process_v2(self, result: SamplerResult, num_virtual_qubits: int):
         quasis = []
         for i in range(len(result)):
             bitstring_counts = result[i].data.meas.get_counts()
@@ -304,7 +297,7 @@ def _post_process_v2(self, result: SamplerResult):
 
             # Convert to quasi-probabilities
             counts = QuasiDistribution(probabilities)
-            quasi_probs = {k: v for k, v in counts.items() if int(k) < 2**self.num_virtual_qubits}
+            quasi_probs = {k: v for k, v in counts.items() if int(k) < 2**num_virtual_qubits}
             quasis.append(quasi_probs)
         return quasis
 


Test patch:
diff --git a/test/gradients/test_estimator_gradient.py b/test/gradients/test_estimator_gradient.py
index ab2a97fac..f075d2900 100644
--- a/test/gradients/test_estimator_gradient.py
+++ b/test/gradients/test_estimator_gradient.py
@@ -25,6 +25,10 @@
 from qiskit.circuit.library.standard_gates import RXXGate, RYYGate, RZXGate, RZZGate
 from qiskit.primitives import Estimator
 from qiskit.quantum_info import SparsePauliOp
+from qiskit.providers.fake_provider import GenericBackendV2
+from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
+
+from qiskit_ibm_runtime import Session, EstimatorV2
 
 from qiskit_machine_learning.gradients import (
     LinCombEstimatorGradient,
@@ -44,16 +48,19 @@
 class TestEstimatorGradient(QiskitAlgorithmsTestCase):
     """Test Estimator Gradient"""
 
+    def __init__(self, TestCase):
+        self.estimator = Estimator()
+        super().__init__(TestCase)
+
     @data(*gradient_factories)
     def test_gradient_operators(self, grad):
         """Test the estimator gradient for different operators"""
-        estimator = Estimator()
         a = Parameter("a")
         qc = QuantumCircuit(1)
         qc.h(0)
         qc.p(a, 0)
         qc.h(0)
-        gradient = grad(estimator)
+        gradient = grad(self.estimator)
         op = SparsePauliOp.from_list([("Z", 1)])
         correct_result = -1 / np.sqrt(2)
         param = [np.pi / 4]
@@ -66,13 +73,13 @@ def test_gradient_operators(self, grad):
     @data(*gradient_factories)
     def test_single_circuit_observable(self, grad):
         """Test the estimator gradient for a single circuit and observable"""
-        estimator = Estimator()
+
         a = Parameter("a")
         qc = QuantumCircuit(1)
         qc.h(0)
         qc.p(a, 0)
         qc.h(0)
-        gradient = grad(estimator)
+        gradient = grad(self.estimator)
         op = SparsePauliOp.from_list([("Z", 1)])
         correct_result = -1 / np.sqrt(2)
         param = [np.pi / 4]
@@ -82,13 +89,13 @@ def test_single_circuit_observable(self, grad):
     @data(*gradient_factories)
     def test_gradient_p(self, grad):
         """Test the estimator gradient for p"""
-        estimator = Estimator()
+
         a = Parameter("a")
         qc = QuantumCircuit(1)
         qc.h(0)
         qc.p(a, 0)
         qc.h(0)
-        gradient = grad(estimator)
+        gradient = grad(self.estimator)
         op = SparsePauliOp.from_list([("Z", 1)])
         param_list = [[np.pi / 4], [0], [np.pi / 2]]
         correct_results = [[-1 / np.sqrt(2)], [0], [-1]]
@@ -100,7 +107,7 @@ def test_gradient_p(self, grad):
     @data(*gradient_factories)
     def test_gradient_u(self, grad):
         """Test the estimator gradient for u"""
-        estimator = Estimator()
+
         a = Parameter("a")
         b = Parameter("b")
         c = Parameter("c")
@@ -108,7 +115,7 @@ def test_gradient_u(self, grad):
         qc.h(0)
         qc.u(a, b, c, 0)
         qc.h(0)
-        gradient = grad(estimator)
+        gradient = grad(self.estimator)
         op = SparsePauliOp.from_list([("Z", 1)])
 
         param_list = [[np.pi / 4, 0, 0], [np.pi / 4, np.pi / 4, np.pi / 4]]
@@ -121,10 +128,10 @@ def test_gradient_u(self, grad):
     @data(*gradient_factories)
     def test_gradient_efficient_su2(self, grad):
         """Test the estimator gradient for EfficientSU2"""
-        estimator = Estimator()
+
         qc = EfficientSU2(2, reps=1)
         op = SparsePauliOp.from_list([("ZI", 1)])
-        gradient = grad(estimator)
+        gradient = grad(self.estimator)
         param_list = [
             [np.pi / 4 for param in qc.parameters],
             [np.pi / 2 for param in qc.parameters],
@@ -149,7 +156,7 @@ def test_gradient_efficient_su2(self, grad):
     @data(*gradient_factories)
     def test_gradient_2qubit_gate(self, grad):
         """Test the estimator gradient for 2 qubit gates"""
-        estimator = Estimator()
+
         for gate in [RXXGate, RYYGate, RZZGate, RZXGate]:
             param_list = [[np.pi / 4], [np.pi / 2]]
             correct_results = [
@@ -160,7 +167,7 @@ def test_gradient_2qubit_gate(self, grad):
             for i, param in enumerate(param_list):
                 a = Parameter("a")
                 qc = QuantumCircuit(2)
-                gradient = grad(estimator)
+                gradient = grad(self.estimator)
 
                 if gate is RZZGate:
                     qc.h([0, 1])
@@ -174,14 +181,14 @@ def test_gradient_2qubit_gate(self, grad):
     @data(*gradient_factories)
     def test_gradient_parameter_coefficient(self, grad):
         """Test the estimator gradient for parameter variables with coefficients"""
-        estimator = Estimator()
+
         qc = RealAmplitudes(num_qubits=2, reps=1)
         qc.rz(qc.parameters[0].exp() + 2 * qc.parameters[1], 0)
         qc.rx(3.0 * qc.parameters[0] + qc.parameters[1].sin(), 1)
         qc.u(qc.parameters[0], qc.parameters[1], qc.parameters[3], 1)
         qc.p(2 * qc.parameters[0] + 1, 0)
         qc.rxx(qc.parameters[0] + 2, 0, 1)
-        gradient = grad(estimator)
+        gradient = grad(self.estimator)
         param_list = [[np.pi / 4 for _ in qc.parameters], [np.pi / 2 for _ in qc.parameters]]
         correct_results = [
             [-0.7266653, -0.4905135, -0.0068606, -0.9228880],
@@ -195,13 +202,13 @@ def test_gradient_parameter_coefficient(self, grad):
     @data(*gradient_factories)
     def test_gradient_parameters(self, grad):
         """Test the estimator gradient for parameters"""
-        estimator = Estimator()
+
         a = Parameter("a")
         b = Parameter("b")
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
         qc.rx(b, 0)
-        gradient = grad(estimator)
+        gradient = grad(self.estimator)
         param_list = [[np.pi / 4, np.pi / 2]]
         correct_results = [
             [-0.70710678],
@@ -229,7 +236,7 @@ def test_gradient_parameters(self, grad):
             param = [[a, b, c], [c, b, a], [a, c], [c, a]]
             op = SparsePauliOp.from_list([("Z", 1)])
             for i, p in enumerate(param):  # pylint: disable=invalid-name
-                gradient = grad(estimator)
+                gradient = grad(self.estimator)
                 gradients = (
                     gradient.run([qc], [op], param_list, parameters=[p]).result().gradients[0]
                 )
@@ -238,14 +245,14 @@ def test_gradient_parameters(self, grad):
     @data(*gradient_factories)
     def test_gradient_multi_arguments(self, grad):
         """Test the estimator gradient for multiple arguments"""
-        estimator = Estimator()
+
         a = Parameter("a")
         b = Parameter("b")
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
         qc2 = QuantumCircuit(1)
         qc2.rx(b, 0)
-        gradient = grad(estimator)
+        gradient = grad(self.estimator)
         param_list = [[np.pi / 4], [np.pi / 2]]
         correct_results = [
             [-0.70710678],
@@ -277,11 +284,11 @@ def test_gradient_multi_arguments(self, grad):
     @data(*gradient_factories)
     def test_gradient_validation(self, grad):
         """Test estimator gradient's validation"""
-        estimator = Estimator()
+
         a = Parameter("a")
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
-        gradient = grad(estimator)
+        gradient = grad(self.estimator)
         param_list = [[np.pi / 4], [np.pi / 2]]
         op = SparsePauliOp.from_list([("Z", 1)])
         with self.assertRaises(ValueError):
@@ -295,9 +302,9 @@ def test_gradient_validation(self, grad):
 
     def test_spsa_gradient(self):
         """Test the SPSA estimator gradient"""
-        estimator = Estimator()
+
         with self.assertRaises(ValueError):
-            _ = SPSAEstimatorGradient(estimator, epsilon=-0.1)
+            _ = SPSAEstimatorGradient(self.estimator, epsilon=-0.1)
         a = Parameter("a")
         b = Parameter("b")
         qc = QuantumCircuit(2)
@@ -306,13 +313,13 @@ def test_spsa_gradient(self):
         param_list = [[1, 1]]
         correct_results = [[-0.84147098, 0.84147098]]
         op = SparsePauliOp.from_list([("ZI", 1)])
-        gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, seed=123)
+        gradient = SPSAEstimatorGradient(self.estimator, epsilon=1e-6, seed=123)
         gradients = gradient.run([qc], [op], param_list).result().gradients
         np.testing.assert_allclose(gradients, correct_results, atol=1e-3)
 
         # multi parameters
         with self.subTest(msg="Multiple parameters"):
-            gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, seed=123)
+            gradient = SPSAEstimatorGradient(self.estimator, epsilon=1e-6, seed=123)
             param_list2 = [[1, 1], [1, 1], [3, 3]]
             gradients2 = (
                 gradient.run([qc] * 3, [op] * 3, param_list2, parameters=[None, [b], None])
@@ -326,13 +333,13 @@ def test_spsa_gradient(self):
         # batch size
         with self.subTest(msg="Batch size"):
             correct_results = [[-0.84147098, 0.1682942]]
-            gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, batch_size=5, seed=123)
+            gradient = SPSAEstimatorGradient(self.estimator, epsilon=1e-6, batch_size=5, seed=123)
             gradients = gradient.run([qc], [op], param_list).result().gradients
             np.testing.assert_allclose(gradients, correct_results, atol=1e-3)
 
         # parameter order
         with self.subTest(msg="The order of gradients"):
-            gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, seed=123)
+            gradient = SPSAEstimatorGradient(self.estimator, epsilon=1e-6, seed=123)
             c = Parameter("c")
             qc = QuantumCircuit(1)
             qc.rx(a, 0)
@@ -348,7 +355,7 @@ def test_spsa_gradient(self):
                 [0.3535525, -0.3535525],
             ]
             for i, p in enumerate(param):  # pylint: disable=invalid-name
-                gradient = SPSAEstimatorGradient(estimator, epsilon=1e-6, seed=123)
+                gradient = SPSAEstimatorGradient(self.estimator, epsilon=1e-6, seed=123)
                 gradients = (
                     gradient.run([qc], [op], param_list3, parameters=[p]).result().gradients[0]
                 )
@@ -447,5 +454,458 @@ def operations_callback(op):
             self.assertAlmostEqual(result.gradients[0].item(), expect, places=5)
 
 
+@ddt
+class TestEstimatorGradientV2(QiskitAlgorithmsTestCase):
+    """Test Estimator Gradient"""
+
+    def __init__(self, TestCase):
+        backend = GenericBackendV2(num_qubits=3, seed=123)
+        session = Session(backend=backend)
+        self.estimator = EstimatorV2(mode=session)
+        self.pass_manager = generate_preset_pass_manager(optimization_level=1, backend=backend)
+        super().__init__(TestCase)
+
+    @data(*gradient_factories)
+    def test_gradient_operators(self, grad):
+        """Test the estimator gradient for different operators"""
+
+        a = Parameter("a")
+        qc = QuantumCircuit(1)
+        qc.h(0)
+        qc.p(a, 0)
+        qc.h(0)
+        gradient = grad(self.estimator, pass_manager=self.pass_manager)
+        op = SparsePauliOp.from_list([("Z", 1)])
+        correct_result = -1 / np.sqrt(2)
+        param = [np.pi / 4]
+        value = gradient.run([qc], [op], [param]).result().gradients[0]
+        self.assertAlmostEqual(value[0], correct_result, 1)
+        op = SparsePauliOp.from_list([("Z", 1)])
+        value = gradient.run([qc], [op], [param]).result().gradients[0]
+        self.assertAlmostEqual(value[0], correct_result, 1)
+
+    @data(*gradient_factories)
+    def test_single_circuit_observable(self, grad):
+        """Test the estimator gradient for a single circuit and observable"""
+
+        a = Parameter("a")
+        qc = QuantumCircuit(1)
+        qc.h(0)
+        qc.p(a, 0)
+        qc.h(0)
+        gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+        op = SparsePauliOp.from_list([("Z", 1)])
+        correct_result = -1 / np.sqrt(2)
+        param = [np.pi / 4]
+        value = gradient.run(qc, op, [param]).result().gradients[0]
+        self.assertAlmostEqual(value[0], correct_result, 1)
+
+    @data(*gradient_factories)
+    def test_gradient_p(self, grad):
+        """Test the estimator gradient for p"""
+
+        a = Parameter("a")
+        qc = QuantumCircuit(1)
+        qc.h(0)
+        qc.p(a, 0)
+        qc.h(0)
+        gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+        op = SparsePauliOp.from_list([("Z", 1)])
+        param_list = [[np.pi / 4], [0], [np.pi / 2]]
+        correct_results = [[-1 / np.sqrt(2)], [0], [-1]]
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [op], [param]).result().gradients[0]
+            for j, value in enumerate(gradients):
+                self.assertAlmostEqual(value, correct_results[i][j], 1)
+
+    @data(*gradient_factories)
+    def test_gradient_u(self, grad):
+        """Test the estimator gradient for u"""
+
+        a = Parameter("a")
+        b = Parameter("b")
+        c = Parameter("c")
+        qc = QuantumCircuit(1)
+        qc.h(0)
+        qc.u(a, b, c, 0)
+        qc.h(0)
+        gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+        op = SparsePauliOp.from_list([("Z", 1)])
+
+        param_list = [[np.pi / 4, 0, 0], [np.pi / 4, np.pi / 4, np.pi / 4]]
+        correct_results = [[-0.70710678, 0.0, 0.0], [-0.35355339, -0.85355339, -0.85355339]]
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [op], [param]).result().gradients[0]
+            for j, value in enumerate(gradients):
+                self.assertAlmostEqual(value, correct_results[i][j], 1)
+
+    @data(*gradient_factories)
+    def test_gradient_efficient_su2(self, grad):
+        """Test the estimator gradient for EfficientSU2"""
+
+        qc = EfficientSU2(2, reps=1)
+        op = SparsePauliOp.from_list([("ZI", 1)])
+        gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+        param_list = [
+            [np.pi / 4 for param in qc.parameters],
+            [np.pi / 2 for param in qc.parameters],
+        ]
+        correct_results = [
+            [
+                -0.35355339,
+                -0.70710678,
+                0,
+                0.35355339,
+                0,
+                -0.70710678,
+                0,
+                0,
+            ],
+            [0, 0, 0, 1, 0, 0, 0, 0],
+        ]
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [op], [param]).result().gradients[0]
+            np.testing.assert_allclose(gradients, correct_results[i], atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    def test_gradient_2qubit_gate(self, grad):
+        """Test the estimator gradient for 2 qubit gates"""
+
+        for gate in [RXXGate, RYYGate, RZZGate, RZXGate]:
+            param_list = [[np.pi / 4], [np.pi / 2]]
+            correct_results = [
+                [-0.70710678],
+                [-1],
+            ]
+            op = SparsePauliOp.from_list([("ZI", 1)])
+            for i, param in enumerate(param_list):
+                a = Parameter("a")
+                qc = QuantumCircuit(2)
+                gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+
+                if gate is RZZGate:
+                    qc.h([0, 1])
+                    qc.append(gate(a), [qc.qubits[0], qc.qubits[1]], [])
+                    qc.h([0, 1])
+                else:
+                    qc.append(gate(a), [qc.qubits[0], qc.qubits[1]], [])
+                gradients = gradient.run([qc], [op], [param]).result().gradients[0]
+                np.testing.assert_allclose(gradients, correct_results[i], atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    def test_gradient_parameter_coefficient(self, grad):
+        """Test the estimator gradient for parameter variables with coefficients"""
+
+        qc = RealAmplitudes(num_qubits=2, reps=1)
+        qc.rz(qc.parameters[0].exp() + 2 * qc.parameters[1], 0)
+        qc.rx(3.0 * qc.parameters[0] + qc.parameters[1].sin(), 1)
+        qc.u(qc.parameters[0], qc.parameters[1], qc.parameters[3], 1)
+        qc.p(2 * qc.parameters[0] + 1, 0)
+        qc.rxx(qc.parameters[0] + 2, 0, 1)
+        gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+        param_list = [[np.pi / 4 for _ in qc.parameters], [np.pi / 2 for _ in qc.parameters]]
+        correct_results = [
+            [-0.7266653, -0.4905135, -0.0068606, -0.9228880],
+            [-3.5972095, 0.10237173, -0.3117748, 0],
+        ]
+        op = SparsePauliOp.from_list([("ZI", 1)])
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [op], [param]).result().gradients[0]
+            np.testing.assert_allclose(gradients, correct_results[i], atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    def test_gradient_parameters(self, grad):
+        """Test the estimator gradient for parameters"""
+
+        a = Parameter("a")
+        b = Parameter("b")
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        qc.rx(b, 0)
+        gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+        param_list = [[np.pi / 4, np.pi / 2]]
+        correct_results = [
+            [-0.70710678],
+        ]
+        op = SparsePauliOp.from_list([("Z", 1)])
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [op], [param], parameters=[[a]]).result().gradients[0]
+            np.testing.assert_allclose(gradients, correct_results[i], atol=1e-1, rtol=1e-1)
+
+        # parameter order
+        with self.subTest(msg="The order of gradients"):
+            c = Parameter("c")
+            qc = QuantumCircuit(1)
+            qc.rx(a, 0)
+            qc.rz(b, 0)
+            qc.rx(c, 0)
+
+            param_list = [[np.pi / 4, np.pi / 2, np.pi / 3]]
+            correct_results = [
+                [-0.35355339, 0.61237244, -0.61237244],
+                [-0.61237244, 0.61237244, -0.35355339],
+                [-0.35355339, -0.61237244],
+                [-0.61237244, -0.35355339],
+            ]
+            param = [[a, b, c], [c, b, a], [a, c], [c, a]]
+            op = SparsePauliOp.from_list([("Z", 1)])
+            for i, p in enumerate(param):  # pylint: disable=invalid-name
+                gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+                gradients = (
+                    gradient.run([qc], [op], param_list, parameters=[p]).result().gradients[0]
+                )
+                np.testing.assert_allclose(gradients, correct_results[i], atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    def test_gradient_multi_arguments(self, grad):
+        """Test the estimator gradient for multiple arguments"""
+
+        a = Parameter("a")
+        b = Parameter("b")
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        qc2 = QuantumCircuit(1)
+        qc2.rx(b, 0)
+        gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+        param_list = [[np.pi / 4], [np.pi / 2]]
+        correct_results = [
+            [-0.70710678],
+            [-1],
+        ]
+        op = SparsePauliOp.from_list([("Z", 1)])
+        gradients = gradient.run([qc, qc2], [op] * 2, param_list).result().gradients
+        np.testing.assert_allclose(gradients, correct_results, atol=1e-1, rtol=1e-1)
+
+        c = Parameter("c")
+        qc3 = QuantumCircuit(1)
+        qc3.rx(c, 0)
+        qc3.ry(a, 0)
+        param_list2 = [[np.pi / 4], [np.pi / 4, np.pi / 4], [np.pi / 4, np.pi / 4]]
+        correct_results2 = [
+            [-0.70710678],
+            [-0.5],
+            [-0.5, -0.5],
+        ]
+        gradients2 = (
+            gradient.run([qc, qc3, qc3], [op] * 3, param_list2, parameters=[[a], [c], None])
+            .result()
+            .gradients
+        )
+        np.testing.assert_allclose(gradients2[0], correct_results2[0], atol=1e-1, rtol=1e-1)
+        np.testing.assert_allclose(gradients2[1], correct_results2[1], atol=1e-1, rtol=1e-1)
+        np.testing.assert_allclose(gradients2[2], correct_results2[2], atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    def test_gradient_validation(self, grad):
+        """Test estimator gradient's validation"""
+
+        a = Parameter("a")
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        gradient = grad(estimator=self.estimator, pass_manager=self.pass_manager)
+        param_list = [[np.pi / 4], [np.pi / 2]]
+        op = SparsePauliOp.from_list([("Z", 1)])
+        with self.assertRaises(ValueError):
+            gradient.run([qc], [op], param_list)
+        with self.assertRaises(ValueError):
+            gradient.run([qc, qc], [op, op], param_list, parameters=[[a]])
+        with self.assertRaises(ValueError):
+            gradient.run([qc, qc], [op], param_list, parameters=[[a]])
+        with self.assertRaises(ValueError):
+            gradient.run([qc], [op], [[np.pi / 4, np.pi / 4]])
+
+    @unittest.skip("Skipping due to noise.")
+    def test_spsa_gradient(self):
+        """Test the SPSA estimator gradient"""
+
+        with self.assertRaises(ValueError):
+            _ = SPSAEstimatorGradient(
+                estimator=self.estimator, pass_manager=self.pass_manager, epsilon=-0.1
+            )
+        a = Parameter("a")
+        b = Parameter("b")
+        qc = QuantumCircuit(2)
+        qc.rx(b, 0)
+        qc.rx(a, 1)
+        param_list = [[1, 1]]
+        correct_results = [[-0.84147098, 0.84147098]]
+        op = SparsePauliOp.from_list([("ZI", 1)])
+        gradient = SPSAEstimatorGradient(
+            estimator=self.estimator, pass_manager=self.pass_manager, epsilon=1e-6, seed=123
+        )
+        gradients = gradient.run([qc], [op], param_list).result().gradients
+        np.testing.assert_allclose(gradients, correct_results, atol=1e-1, rtol=1e-1)
+
+        # multi parameters
+        with self.subTest(msg="Multiple parameters"):
+            gradient = SPSAEstimatorGradient(
+                estimator=self.estimator, pass_manager=self.pass_manager, epsilon=1e-6, seed=123
+            )
+            param_list2 = [[1, 1], [1, 1], [3, 3]]
+            gradients2 = (
+                gradient.run([qc] * 3, [op] * 3, param_list2, parameters=[None, [b], None])
+                .result()
+                .gradients
+            )
+            correct_results2 = [[-0.84147098, 0.84147098], [0.84147098], [-0.14112001, 0.14112001]]
+            for grad, correct in zip(gradients2, correct_results2):
+                np.testing.assert_allclose(grad, correct, atol=1e-1, rtol=1e-1)
+
+        # batch size
+        with self.subTest(msg="Batch size"):
+            correct_results = [[-0.84147098, 0.1682942]]
+            gradient = SPSAEstimatorGradient(
+                estimator=self.estimator,
+                pass_manager=self.pass_manager,
+                epsilon=1e-6,
+                batch_size=5,
+                seed=123,
+            )
+            gradients = gradient.run([qc], [op], param_list).result().gradients
+            np.testing.assert_allclose(gradients, correct_results, atol=1e-1, rtol=1e-1)
+
+        # parameter order
+        with self.subTest(msg="The order of gradients"):
+            gradient = SPSAEstimatorGradient(
+                estimator=self.estimator, pass_manager=self.pass_manager, epsilon=1e-6, seed=123
+            )
+            c = Parameter("c")
+            qc = QuantumCircuit(1)
+            qc.rx(a, 0)
+            qc.rz(b, 0)
+            qc.rx(c, 0)
+            op = SparsePauliOp.from_list([("Z", 1)])
+            param_list3 = [[np.pi / 4, np.pi / 2, np.pi / 3]]
+            expected = [
+                [-0.3535525, 0.3535525, 0.3535525],
+                [0.3535525, 0.3535525, -0.3535525],
+                [-0.3535525, 0.3535525],
+                [0.3535525, -0.3535525],
+            ]
+            param = [[a, b, c], [c, b, a], [a, c], [c, a]]
+            for i, p in enumerate(param):  # pylint: disable=invalid-name
+                gradient = SPSAEstimatorGradient(
+                    estimator=self.estimator, pass_manager=self.pass_manager, epsilon=1e-6, seed=123
+                )
+                gradients = (
+                    gradient.run([qc], [op], param_list3, parameters=[p]).result().gradients[0]
+                )
+                np.testing.assert_allclose(gradients, expected[i], atol=1e-1, rtol=1e-1)
+
+    @data(
+        ParamShiftEstimatorGradient,
+        LinCombEstimatorGradient,
+        SPSAEstimatorGradient,
+    )
+    @unittest.skip("Options needs to be added for V2.")
+    def test_options(self, grad):
+        """Test estimator gradient's run options"""
+        a = Parameter("a")
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        op = SparsePauliOp.from_list([("Z", 1)])
+        estimator = EstimatorV2(options={"shots": 100})
+        with self.subTest("estimator"):
+            if grad is SPSAEstimatorGradient:
+                gradient = grad(estimator=estimator, pass_manager=self.pass_manager, epsilon=1e-6)
+            else:
+                gradient = grad(estimator=estimator, pass_manager=self.pass_manager)
+            options = gradient.options
+            result = gradient.run([qc], [op], [[1]]).result()
+            self.assertEqual(result.options.get("shots"), 100)
+            self.assertEqual(options.get("shots"), 100)
+
+        with self.subTest("gradient init"):
+            if grad is SPSAEstimatorGradient:
+                gradient = grad(
+                    estimator=self.estimator,
+                    pass_manager=self.pass_manager,
+                    epsilon=1e-6,
+                    options={"shots": 200},
+                )
+            else:
+                gradient = grad(
+                    estimator=self.estimator, pass_manager=self.pass_manager, options={"shots": 200}
+                )
+            options = gradient.options
+            result = gradient.run([qc], [op], [[1]]).result()
+            self.assertEqual(result.options.get("shots"), 200)
+            self.assertEqual(options.get("shots"), 200)
+
+        with self.subTest("gradient update"):
+            if grad is SPSAEstimatorGradient:
+                gradient = grad(
+                    estimator=self.estimator,
+                    pass_manager=self.pass_manager,
+                    epsilon=1e-6,
+                    options={"shots": 200},
+                )
+            else:
+                gradient = grad(
+                    estimator=self.estimator, pass_manager=self.pass_manager, options={"shots": 200}
+                )
+            gradient.update_default_options(shots=100)
+            options = gradient.options
+            result = gradient.run([qc], [op], [[1]]).result()
+            self.assertEqual(result.options.get("shots"), 100)
+            self.assertEqual(options.get("shots"), 100)
+
+        with self.subTest("gradient run"):
+            if grad is SPSAEstimatorGradient:
+                gradient = grad(
+                    estimator=self.estimator,
+                    pass_manager=self.pass_manager,
+                    epsilon=1e-6,
+                    options={"shots": 200},
+                )
+            else:
+                gradient = grad(
+                    estimator=self.estimator, pass_manager=self.pass_manager, options={"shots": 200}
+                )
+            options = gradient.options
+            result = gradient.run([qc], [op], [[1]], shots=300).result()
+            self.assertEqual(result.options.get("shots"), 300)
+            # Only default + estimator options. Not run.
+            self.assertEqual(options.get("shots"), 200)
+
+    @data(
+        ParamShiftEstimatorGradient,
+        LinCombEstimatorGradient,
+        SPSAEstimatorGradient,
+    )
+    def test_operations_preserved(self, gradient_cls):
+        """Test non-parameterized instructions are preserved and not unrolled."""
+        x = Parameter("x")
+        circuit = QuantumCircuit(2)
+        circuit.initialize([0.5, 0.5, 0.5, 0.5])  # this should remain as initialize
+        circuit.crx(x, 0, 1)  # this should get unrolled
+
+        values = [np.pi / 2]
+        expect = -1 / (2 * np.sqrt(2))
+
+        observable = SparsePauliOp(["XX"])
+
+        ops = []
+
+        def operations_callback(op):
+            ops.append(op)
+
+        estimator = LoggingEstimator(operations_callback=operations_callback)
+
+        if gradient_cls in [SPSAEstimatorGradient]:
+            gradient = gradient_cls(estimator, epsilon=0.01)
+        else:
+            gradient = gradient_cls(estimator)
+
+        job = gradient.run([circuit], [observable], [values])
+        result = job.result()
+
+        with self.subTest(msg="assert initialize is preserved"):
+            self.assertTrue(all("initialize" in ops_i[0].keys() for ops_i in ops))
+
+        with self.subTest(msg="assert result is correct"):
+            self.assertAlmostEqual(result.gradients[0].item(), expect, places=5)
+
+
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/gradients/test_qfi.py b/test/gradients/test_qfi.py
deleted file mode 100644
index 8acbcaf09..000000000
--- a/test/gradients/test_qfi.py
+++ /dev/null
@@ -1,151 +0,0 @@
-# This code is part of a Qiskit project.
-#
-# (C) Copyright IBM 2022, 2024.
-#
-# This code is licensed under the Apache License, Version 2.0. You may
-# obtain a copy of this license in the LICENSE.txt file in the root directory
-# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
-#
-# Any modifications or derivative works of this code must retain this
-# copyright notice, and modified files need to carry a notice indicating
-# that they have been altered from the originals.
-# =============================================================================
-
-"""Test QFI."""
-
-import unittest
-from test import QiskitAlgorithmsTestCase
-
-from ddt import ddt, data
-import numpy as np
-
-from qiskit import QuantumCircuit
-from qiskit.circuit import Parameter
-from qiskit.circuit.parametervector import ParameterVector
-from qiskit.primitives import Estimator
-
-from qiskit_machine_learning.gradients import LinCombQGT, QFI, DerivativeType
-
-
-@ddt
-class TestQFI(QiskitAlgorithmsTestCase):
-    """Test QFI"""
-
-    def setUp(self):
-        super().setUp()
-        self.estimator = Estimator()
-        self.lcu_qgt = LinCombQGT(self.estimator, derivative_type=DerivativeType.REAL)
-
-    def test_qfi(self):
-        """Test if the quantum fisher information calculation is correct for a simple test case.
-        QFI = [[1, 0], [0, 1]] - [[0, 0], [0, cos^2(a)]]
-        """
-        # create the circuit
-        a, b = Parameter("a"), Parameter("b")
-        qc = QuantumCircuit(1)
-        qc.h(0)
-        qc.rz(a, 0)
-        qc.rx(b, 0)
-
-        param_list = [[np.pi / 4, 0.1], [np.pi, 0.1], [np.pi / 2, 0.1]]
-        correct_values = [[[1, 0], [0, 0.5]], [[1, 0], [0, 0]], [[1, 0], [0, 1]]]
-
-        qfi = QFI(self.lcu_qgt)
-        for i, param in enumerate(param_list):
-            qfis = qfi.run([qc], [param]).result().qfis
-            np.testing.assert_allclose(qfis[0], correct_values[i], atol=1e-3)
-
-    def test_qfi_phase_fix(self):
-        """Test the phase-fix argument in the QFI calculation"""
-        # create the circuit
-        a, b = Parameter("a"), Parameter("b")
-        qc = QuantumCircuit(1)
-        qc.h(0)
-        qc.rz(a, 0)
-        qc.rx(b, 0)
-
-        param = [np.pi / 4, 0.1]
-        # test for different values
-        correct_values = [[1, 0], [0, 1]]
-        qgt = LinCombQGT(self.estimator, phase_fix=False)
-        qfi = QFI(qgt)
-        qfis = qfi.run([qc], [param]).result().qfis
-        np.testing.assert_allclose(qfis[0], correct_values, atol=1e-3)
-
-    @data("lcu")
-    def test_qfi_maxcut(self, qgt_kind):  # pylint: disable=unused-argument
-        """Test the QFI for a simple MaxCut problem.
-
-        This is interesting because it contains the same parameters in different gates.
-        """
-        # create maxcut circuit for the hamiltonian
-        # H = (I ^ I ^ Z ^ Z) + (I ^ Z ^ I ^ Z) + (Z ^ I ^ I ^ Z) + (I ^ Z ^ Z ^ I)
-
-        x = ParameterVector("x", 2)
-        ansatz = QuantumCircuit(4)
-
-        # initial hadamard layer
-        ansatz.h(ansatz.qubits)
-
-        # e^{iZZ} layers
-        def expiz(qubit0, qubit1):
-            ansatz.cx(qubit0, qubit1)
-            ansatz.rz(2 * x[0], qubit1)
-            ansatz.cx(qubit0, qubit1)
-
-        expiz(2, 1)
-        expiz(3, 0)
-        expiz(2, 0)
-        expiz(1, 0)
-
-        # mixer layer with RX gates
-        for i in range(ansatz.num_qubits):
-            ansatz.rx(2 * x[1], i)
-
-        reference = np.array([[16.0, -5.551], [-5.551, 18.497]])
-        param = [0.4, 0.69]
-
-        qgt = self.lcu_qgt
-        qfi = QFI(qgt)
-        qfi_result = qfi.run([ansatz], [param]).result().qfis
-        np.testing.assert_array_almost_equal(qfi_result[0], reference, decimal=3)
-
-    def test_options(self):
-        """Test QFI's options"""
-        a = Parameter("a")
-        qc = QuantumCircuit(1)
-        qc.rx(a, 0)
-        qgt = LinCombQGT(estimator=self.estimator, options={"shots": 100})
-
-        with self.subTest("QGT"):
-            qfi = QFI(qgt=qgt)
-            options = qfi.options
-            result = qfi.run([qc], [[1]]).result()
-            self.assertEqual(result.options.get("shots"), 100)
-            self.assertEqual(options.get("shots"), 100)
-
-        with self.subTest("QFI init"):
-            qfi = QFI(qgt=qgt, options={"shots": 200})
-            result = qfi.run([qc], [[1]]).result()
-            options = qfi.options
-            self.assertEqual(result.options.get("shots"), 200)
-            self.assertEqual(options.get("shots"), 200)
-
-        with self.subTest("QFI update"):
-            qfi = QFI(qgt, options={"shots": 200})
-            qfi.update_default_options(shots=100)
-            options = qfi.options
-            result = qfi.run([qc], [[1]]).result()
-            self.assertEqual(result.options.get("shots"), 100)
-            self.assertEqual(options.get("shots"), 100)
-
-        with self.subTest("QFI run"):
-            qfi = QFI(qgt=qgt, options={"shots": 200})
-            result = qfi.run([qc], [[0]], shots=300).result()
-            options = qfi.options
-            self.assertEqual(result.options.get("shots"), 300)
-            self.assertEqual(options.get("shots"), 200)
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/test/gradients/test_qgt.py b/test/gradients/test_qgt.py
deleted file mode 100644
index cdab18353..000000000
--- a/test/gradients/test_qgt.py
+++ /dev/null
@@ -1,310 +0,0 @@
-# This code is part of a Qiskit project.
-#
-# (C) Copyright IBM 2022, 2024.
-#
-# This code is licensed under the Apache License, Version 2.0. You may
-# obtain a copy of this license in the LICENSE.txt file in the root directory
-# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
-#
-# Any modifications or derivative works of this code must retain this
-# copyright notice, and modified files need to carry a notice indicating
-# that they have been altered from the originals.
-# =============================================================================
-
-"""Test QGT."""
-
-import unittest
-from test import QiskitAlgorithmsTestCase
-
-from ddt import ddt, data
-import numpy as np
-
-from qiskit import QuantumCircuit
-from qiskit.circuit import Parameter
-from qiskit.circuit.library import RealAmplitudes
-from qiskit.primitives import Estimator
-
-from qiskit_machine_learning.gradients import DerivativeType, LinCombQGT
-
-from .logging_primitives import LoggingEstimator
-
-
-@ddt
-class TestQGT(QiskitAlgorithmsTestCase):
-    """Test QGT"""
-
-    def setUp(self):
-        super().setUp()
-        self.estimator = Estimator()
-
-    @data(LinCombQGT)
-    def test_qgt_derivative_type(self, qgt_type):
-        """Test QGT derivative_type"""
-        args = (self.estimator,)
-        qgt = qgt_type(*args, derivative_type=DerivativeType.REAL)
-
-        a, b = Parameter("a"), Parameter("b")
-        qc = QuantumCircuit(1)
-        qc.h(0)
-        qc.rz(a, 0)
-        qc.rx(b, 0)
-
-        param_list = [[np.pi / 4, 0], [np.pi / 2, np.pi / 4]]
-        correct_values = [
-            np.array([[1, 0.707106781j], [-0.707106781j, 0.5]]) / 4,
-            np.array([[1, 1j], [-1j, 1]]) / 4,
-        ]
-
-        # test real derivative
-        with self.subTest("Test with DerivativeType.REAL"):
-            qgt.derivative_type = DerivativeType.REAL
-            for i, param in enumerate(param_list):
-                qgt_result = qgt.run([qc], [param]).result().qgts
-                np.testing.assert_allclose(qgt_result[0], correct_values[i].real, atol=1e-3)
-
-        # test imaginary derivative
-        with self.subTest("Test with DerivativeType.IMAG"):
-            qgt.derivative_type = DerivativeType.IMAG
-            for i, param in enumerate(param_list):
-                qgt_result = qgt.run([qc], [param]).result().qgts
-                np.testing.assert_allclose(qgt_result[0], correct_values[i].imag, atol=1e-3)
-
-        # test real + imaginary derivative
-        with self.subTest("Test with DerivativeType.COMPLEX"):
-            qgt.derivative_type = DerivativeType.COMPLEX
-            for i, param in enumerate(param_list):
-                qgt_result = qgt.run([qc], [param]).result().qgts
-                np.testing.assert_allclose(qgt_result[0], correct_values[i], atol=1e-3)
-
-    @data(LinCombQGT)
-    def test_qgt_phase_fix(self, qgt_type):
-        """Test the phase-fix argument in a QGT calculation"""
-        args = (self.estimator,)
-        qgt = qgt_type(*args, phase_fix=False)
-
-        # create the circuit
-        a, b = Parameter("a"), Parameter("b")
-        qc = QuantumCircuit(1)
-        qc.h(0)
-        qc.rz(a, 0)
-        qc.rx(b, 0)
-
-        param_list = [[np.pi / 4, 0], [np.pi / 2, np.pi / 4]]
-        correct_values = [
-            np.array([[1, 0.707106781j], [-0.707106781j, 1]]) / 4,
-            np.array([[1, 1j], [-1j, 1]]) / 4,
-        ]
-
-        # test real derivative
-        with self.subTest("Test phase fix with DerivativeType.REAL"):
-            qgt.derivative_type = DerivativeType.REAL
-            for i, param in enumerate(param_list):
-                qgt_result = qgt.run([qc], [param]).result().qgts
-                np.testing.assert_allclose(qgt_result[0], correct_values[i].real, atol=1e-3)
-
-        # test imaginary derivative
-        with self.subTest("Test phase fix with DerivativeType.IMAG"):
-            qgt.derivative_type = DerivativeType.IMAG
-            for i, param in enumerate(param_list):
-                qgt_result = qgt.run([qc], [param]).result().qgts
-                np.testing.assert_allclose(qgt_result[0], correct_values[i].imag, atol=1e-3)
-
-        # test real + imaginary derivative
-        with self.subTest("Test phase fix with DerivativeType.COMPLEX"):
-            qgt.derivative_type = DerivativeType.COMPLEX
-            for i, param in enumerate(param_list):
-                qgt_result = qgt.run([qc], [param]).result().qgts
-                np.testing.assert_allclose(qgt_result[0], correct_values[i], atol=1e-3)
-
-    @data(LinCombQGT)
-    def test_qgt_coefficients(self, qgt_type):
-        """Test the derivative option of QGT"""
-        args = (self.estimator,)
-        qgt = qgt_type(*args, derivative_type=DerivativeType.REAL)
-
-        qc = RealAmplitudes(num_qubits=2, reps=1)
-        qc.rz(qc.parameters[0].exp() + 2 * qc.parameters[1], 0)
-        qc.rx(3.0 * qc.parameters[2] + qc.parameters[3].sin(), 1)
-
-        # test imaginary derivative
-        param_list = [
-            [np.pi / 4 for param in qc.parameters],
-            [np.pi / 2 for param in qc.parameters],
-        ]
-        correct_values = (
-            np.array(
-                [
-                    [
-                        [5.707309, 4.2924833, 1.5295868, 0.1938604],
-                        [4.2924833, 4.9142136, 0.75, 0.8838835],
-                        [1.5295868, 0.75, 3.4430195, 0.0758252],
-                        [0.1938604, 0.8838835, 0.0758252, 1.1357233],
-                    ],
-                    [
-                        [1.0, 0.0, 1.0, 0.0],
-                        [0.0, 1.0, 0.0, 0.0],
-                        [1.0, 0.0, 10.0, -0.0],
-                        [0.0, 0.0, -0.0, 1.0],
-                    ],
-                ]
-            )
-            / 4
-        )
-        for i, param in enumerate(param_list):
-            qgt_result = qgt.run([qc], [param]).result().qgts
-            np.testing.assert_allclose(qgt_result[0], correct_values[i], atol=1e-3)
-
-    @data(LinCombQGT)
-    def test_qgt_parameters(self, qgt_type):
-        """Test the QGT with specified parameters"""
-        args = (self.estimator,)
-        qgt = qgt_type(*args, derivative_type=DerivativeType.REAL)
-
-        a = Parameter("a")
-        b = Parameter("b")
-        qc = QuantumCircuit(1)
-        qc.rx(a, 0)
-        qc.ry(b, 0)
-        param_values = [np.pi / 4, np.pi / 4]
-        qgt_result = qgt.run([qc], [param_values], [[a]]).result().qgts
-        np.testing.assert_allclose(qgt_result[0], [[1 / 4]], atol=1e-3)
-
-        with self.subTest("Test with different parameter orders"):
-            c = Parameter("c")
-            qc2 = QuantumCircuit(1)
-            qc2.rx(a, 0)
-            qc2.rz(b, 0)
-            qc2.rx(c, 0)
-            param_values = [np.pi / 4, np.pi / 4, np.pi / 4]
-            params = [[a, b, c], [c, b, a], [a, c], [b, a]]
-            expected = [
-                np.array(
-                    [
-                        [0.25, 0.0, 0.1767767],
-                        [0.0, 0.125, -0.08838835],
-                        [0.1767767, -0.08838835, 0.1875],
-                    ]
-                ),
-                np.array(
-                    [
-                        [0.1875, -0.08838835, 0.1767767],
-                        [-0.08838835, 0.125, 0.0],
-                        [0.1767767, 0.0, 0.25],
-                    ]
-                ),
-                np.array([[0.25, 0.1767767], [0.1767767, 0.1875]]),
-                np.array([[0.125, 0.0], [0.0, 0.25]]),
-            ]
-            for i, param in enumerate(params):
-                qgt_result = qgt.run([qc2], [param_values], [param]).result().qgts
-                np.testing.assert_allclose(qgt_result[0], expected[i], atol=1e-3)
-
-    @data(LinCombQGT)
-    def test_qgt_multi_arguments(self, qgt_type):
-        """Test the QGT for multiple arguments"""
-        args = (self.estimator,)
-        qgt = qgt_type(*args, derivative_type=DerivativeType.REAL)
-
-        a = Parameter("a")
-        b = Parameter("b")
-        qc = QuantumCircuit(1)
-        qc.rx(a, 0)
-        qc.ry(b, 0)
-        qc2 = QuantumCircuit(1)
-        qc2.rx(a, 0)
-        qc2.ry(b, 0)
-
-        param_list = [[np.pi / 4], [np.pi / 2]]
-        correct_values = [[[1 / 4]], [[1 / 4, 0], [0, 0]]]
-        param_list = [[np.pi / 4, np.pi / 4], [np.pi / 2, np.pi / 2]]
-        qgt_results = qgt.run([qc, qc2], param_list, [[a], None]).result().qgts
-        for i, _ in enumerate(param_list):
-            np.testing.assert_allclose(qgt_results[i], correct_values[i], atol=1e-3)
-
-    @data(LinCombQGT)
-    def test_qgt_validation(self, qgt_type):
-        """Test estimator QGT's validation"""
-        args = (self.estimator,)
-        qgt = qgt_type(*args)
-
-        a = Parameter("a")
-        qc = QuantumCircuit(1)
-        qc.rx(a, 0)
-        parameter_values = [[np.pi / 4]]
-        with self.subTest("assert number of circuits does not match"):
-            with self.assertRaises(ValueError):
-                qgt.run([qc, qc], parameter_values)
-        with self.subTest("assert number of parameter values does not match"):
-            with self.assertRaises(ValueError):
-                qgt.run([qc], [[np.pi / 4], [np.pi / 2]])
-        with self.subTest("assert number of parameters does not match"):
-            with self.assertRaises(ValueError):
-                qgt.run([qc], parameter_values, parameters=[[a], [a]])
-
-    def test_options(self):
-        """Test QGT's options"""
-        a = Parameter("a")
-        qc = QuantumCircuit(1)
-        qc.rx(a, 0)
-        estimator = Estimator(options={"shots": 100})
-
-        with self.subTest("estimator"):
-            qgt = LinCombQGT(estimator)
-            options = qgt.options
-            result = qgt.run([qc], [[1]]).result()
-            self.assertEqual(result.options.get("shots"), 100)
-            self.assertEqual(options.get("shots"), 100)
-
-        with self.subTest("QGT init"):
-            qgt = LinCombQGT(estimator, options={"shots": 200})
-            result = qgt.run([qc], [[1]]).result()
-            options = qgt.options
-            self.assertEqual(result.options.get("shots"), 200)
-            self.assertEqual(options.get("shots"), 200)
-
-        with self.subTest("QGT update"):
-            qgt = LinCombQGT(estimator, options={"shots": 200})
-            qgt.update_default_options(shots=100)
-            options = qgt.options
-            result = qgt.run([qc], [[1]]).result()
-            self.assertEqual(result.options.get("shots"), 100)
-            self.assertEqual(options.get("shots"), 100)
-
-        with self.subTest("QGT run"):
-            qgt = LinCombQGT(estimator, options={"shots": 200})
-            result = qgt.run([qc], [[0]], shots=300).result()
-            options = qgt.options
-            self.assertEqual(result.options.get("shots"), 300)
-            self.assertEqual(options.get("shots"), 200)
-
-    def test_operations_preserved(self):
-        """Test non-parameterized instructions are preserved and not unrolled."""
-        x, y = Parameter("x"), Parameter("y")
-        circuit = QuantumCircuit(2)
-        circuit.initialize([0.5, 0.5, 0.5, 0.5])  # this should remain as initialize
-        circuit.crx(x, 0, 1)  # this should get unrolled
-        circuit.ry(y, 0)
-
-        values = [np.pi / 2, np.pi]
-        expect = np.diag([0.25, 0.5]) / 4
-
-        ops = []
-
-        def operations_callback(op):
-            ops.append(op)
-
-        estimator = LoggingEstimator(operations_callback=operations_callback)
-        qgt = LinCombQGT(estimator, derivative_type=DerivativeType.REAL)
-
-        job = qgt.run([circuit], [values])
-        result = job.result()
-
-        with self.subTest(msg="assert initialize is preserved"):
-            self.assertTrue(all("initialize" in ops_i[0].keys() for ops_i in ops))
-
-        with self.subTest(msg="assert result is correct"):
-            np.testing.assert_allclose(result.qgts[0], expect, atol=1e-5)
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/test/gradients/test_sampler_gradient.py b/test/gradients/test_sampler_gradient.py
index d586dbe4b..da2408a94 100644
--- a/test/gradients/test_sampler_gradient.py
+++ b/test/gradients/test_sampler_gradient.py
@@ -16,7 +16,6 @@
 import unittest
 from test import QiskitAlgorithmsTestCase
 from typing import List
-
 import numpy as np
 from ddt import ddt, data
 
@@ -26,6 +25,10 @@
 from qiskit.circuit.library.standard_gates import RXXGate
 from qiskit.primitives import Sampler
 from qiskit.result import QuasiDistribution
+from qiskit.providers.fake_provider import GenericBackendV2
+from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
+
+from qiskit_ibm_runtime import Session, SamplerV2
 
 from qiskit_machine_learning.gradients import (
     LinCombSamplerGradient,
@@ -45,17 +48,21 @@
 class TestSamplerGradient(QiskitAlgorithmsTestCase):
     """Test Sampler Gradient"""
 
+    def __init__(self, TestCase):
+        self.sampler = Sampler()
+        super().__init__(TestCase)
+
     @data(*gradient_factories)
     def test_single_circuit(self, grad):
         """Test the sampler gradient for a single circuit"""
-        sampler = Sampler()
+
         a = Parameter("a")
         qc = QuantumCircuit(1)
         qc.h(0)
         qc.p(a, 0)
         qc.h(0)
         qc.measure_all()
-        gradient = grad(sampler)
+        gradient = grad(self.sampler)
         param_list = [[np.pi / 4], [0], [np.pi / 2]]
         expected = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
@@ -71,14 +78,14 @@ def test_single_circuit(self, grad):
     @data(*gradient_factories)
     def test_gradient_p(self, grad):
         """Test the sampler gradient for p"""
-        sampler = Sampler()
+
         a = Parameter("a")
         qc = QuantumCircuit(1)
         qc.h(0)
         qc.p(a, 0)
         qc.h(0)
         qc.measure_all()
-        gradient = grad(sampler)
+        gradient = grad(self.sampler)
         param_list = [[np.pi / 4], [0], [np.pi / 2]]
         expected = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
@@ -94,7 +101,7 @@ def test_gradient_p(self, grad):
     @data(*gradient_factories)
     def test_gradient_u(self, grad):
         """Test the sampler gradient for u"""
-        sampler = Sampler()
+
         a = Parameter("a")
         b = Parameter("b")
         c = Parameter("c")
@@ -103,7 +110,7 @@ def test_gradient_u(self, grad):
         qc.u(a, b, c, 0)
         qc.h(0)
         qc.measure_all()
-        gradient = grad(sampler)
+        gradient = grad(self.sampler)
         param_list = [[np.pi / 4, 0, 0], [np.pi / 4, np.pi / 4, np.pi / 4]]
         expected = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}, {0: 0, 1: 0}, {0: 0, 1: 0}],
@@ -118,10 +125,10 @@ def test_gradient_u(self, grad):
     @data(*gradient_factories)
     def test_gradient_efficient_su2(self, grad):
         """Test the sampler gradient for EfficientSU2"""
-        sampler = Sampler()
+
         qc = EfficientSU2(2, reps=1)
         qc.measure_all()
-        gradient = grad(sampler)
+        gradient = grad(self.sampler)
         param_list = [
             [np.pi / 4 for param in qc.parameters],
             [np.pi / 2 for param in qc.parameters],
@@ -212,7 +219,7 @@ def test_gradient_efficient_su2(self, grad):
     @data(*gradient_factories)
     def test_gradient_2qubit_gate(self, grad):
         """Test the sampler gradient for 2 qubit gates"""
-        sampler = Sampler()
+
         for gate in [RXXGate]:
             param_list = [[np.pi / 4], [np.pi / 2]]
             correct_results = [
@@ -224,7 +231,7 @@ def test_gradient_2qubit_gate(self, grad):
                 qc = QuantumCircuit(2)
                 qc.append(gate(a), [qc.qubits[0], qc.qubits[1]], [])
                 qc.measure_all()
-                gradient = grad(sampler)
+                gradient = grad(self.sampler)
                 gradients = gradient.run([qc], [param]).result().gradients[0]
                 array1 = _quasi2array(gradients, num_qubits=2)
                 array2 = _quasi2array(correct_results[i], num_qubits=2)
@@ -233,7 +240,7 @@ def test_gradient_2qubit_gate(self, grad):
     @data(*gradient_factories)
     def test_gradient_parameter_coefficient(self, grad):
         """Test the sampler gradient for parameter variables with coefficients"""
-        sampler = Sampler()
+
         qc = RealAmplitudes(num_qubits=2, reps=1)
         qc.rz(qc.parameters[0].exp() + 2 * qc.parameters[1], 0)
         qc.rx(3.0 * qc.parameters[0] + qc.parameters[1].sin(), 1)
@@ -241,7 +248,7 @@ def test_gradient_parameter_coefficient(self, grad):
         qc.p(2 * qc.parameters[0] + 1, 0)
         qc.rxx(qc.parameters[0] + 2, 0, 1)
         qc.measure_all()
-        gradient = grad(sampler)
+        gradient = grad(self.sampler)
         param_list = [[np.pi / 4 for _ in qc.parameters], [np.pi / 2 for _ in qc.parameters]]
         correct_results = [
             [
@@ -307,14 +314,14 @@ def test_gradient_parameter_coefficient(self, grad):
     @data(*gradient_factories)
     def test_gradient_parameters(self, grad):
         """Test the sampler gradient for parameters"""
-        sampler = Sampler()
+
         a = Parameter("a")
         b = Parameter("b")
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
         qc.rz(b, 0)
         qc.measure_all()
-        gradient = grad(sampler)
+        gradient = grad(self.sampler)
         param_list = [[np.pi / 4, np.pi / 2]]
         expected = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
@@ -364,7 +371,7 @@ def test_gradient_parameters(self, grad):
     @data(*gradient_factories)
     def test_gradient_multi_arguments(self, grad):
         """Test the sampler gradient for multiple arguments"""
-        sampler = Sampler()
+
         a = Parameter("a")
         b = Parameter("b")
         qc = QuantumCircuit(1)
@@ -373,7 +380,7 @@ def test_gradient_multi_arguments(self, grad):
         qc2 = QuantumCircuit(1)
         qc2.rx(b, 0)
         qc2.measure_all()
-        gradient = grad(sampler)
+        gradient = grad(self.sampler)
         param_list = [[np.pi / 4], [np.pi / 2]]
         correct_results = [
             [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
@@ -411,12 +418,12 @@ def test_gradient_multi_arguments(self, grad):
     @data(*gradient_factories)
     def test_gradient_validation(self, grad):
         """Test sampler gradient's validation"""
-        sampler = Sampler()
+
         a = Parameter("a")
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
         qc.measure_all()
-        gradient = grad(sampler)
+        gradient = grad(self.sampler)
         param_list = [[np.pi / 4], [np.pi / 2]]
         with self.assertRaises(ValueError):
             gradient.run([qc], param_list)
@@ -427,9 +434,9 @@ def test_gradient_validation(self, grad):
 
     def test_spsa_gradient(self):
         """Test the SPSA sampler gradient"""
-        sampler = Sampler()
+
         with self.assertRaises(ValueError):
-            _ = SPSASamplerGradient(sampler, epsilon=-0.1)
+            _ = SPSASamplerGradient(self.sampler, epsilon=-0.1)
 
         a = Parameter("a")
         b = Parameter("b")
@@ -445,7 +452,7 @@ def test_spsa_gradient(self):
                 {0: -0.2273244, 1: 0.6480598, 2: -0.2273244, 3: -0.1934111},
             ],
         ]
-        gradient = SPSASamplerGradient(sampler, epsilon=1e-6, seed=123)
+        gradient = SPSASamplerGradient(self.sampler, epsilon=1e-6, seed=123)
         for i, param in enumerate(param_list):
             gradients = gradient.run([qc], [param]).result().gradients[0]
             array1 = _quasi2array(gradients, num_qubits=2)
@@ -468,7 +475,7 @@ def test_spsa_gradient(self):
                     {0: 0.0141129, 1: 0.0564471, 2: 0.3642884, 3: -0.4348484},
                 ],
             ]
-            gradient = SPSASamplerGradient(sampler, epsilon=1e-6, seed=123)
+            gradient = SPSASamplerGradient(self.sampler, epsilon=1e-6, seed=123)
             gradients = (
                 gradient.run([qc] * 3, param_list2, parameters=[None, [b], None]).result().gradients
             )
@@ -480,7 +487,7 @@ def test_spsa_gradient(self):
         # batch size
         with self.subTest(msg="Batch size"):
             param_list = [[1, 1]]
-            gradient = SPSASamplerGradient(sampler, epsilon=1e-6, batch_size=4, seed=123)
+            gradient = SPSASamplerGradient(self.sampler, epsilon=1e-6, batch_size=4, seed=123)
             gradients = gradient.run([qc], param_list).result().gradients
             correct_results3 = [
                 [
@@ -533,7 +540,7 @@ def test_spsa_gradient(self):
                 ],
             ]
             for i, p in enumerate(param):  # pylint: disable=invalid-name
-                gradient = SPSASamplerGradient(sampler, epsilon=1e-6, seed=123)
+                gradient = SPSASamplerGradient(self.sampler, epsilon=1e-6, seed=123)
                 gradients = gradient.run([qc], param_list, parameters=[p]).result().gradients[0]
                 array1 = _quasi2array(gradients, num_qubits=1)
                 array2 = _quasi2array(correct_results[i], num_qubits=1)
@@ -544,54 +551,583 @@ def test_spsa_gradient(self):
         LinCombSamplerGradient,
         SPSASamplerGradient,
     )
-    def test_options(self, grad):
-        """Test sampler gradient's run options"""
+    def test_operations_preserved(self, gradient_cls):
+        """Test non-parameterized instructions are preserved and not unrolled."""
+        x = Parameter("x")
+        circuit = QuantumCircuit(2)
+        circuit.initialize(np.array([1, 1, 0, 0]) / np.sqrt(2))  # this should remain as initialize
+        circuit.crx(x, 0, 1)  # this should get unrolled
+        circuit.measure_all()
+
+        values = [np.pi / 2]
+        expect = [{0: 0, 1: -0.25, 2: 0, 3: 0.25}]
+
+        ops = []
+
+        def operations_callback(op):
+            ops.append(op)
+
+        sampler = LoggingSampler(operations_callback=operations_callback)
+
+        if gradient_cls in [SPSASamplerGradient]:
+            gradient = gradient_cls(sampler, epsilon=0.01)
+        else:
+            gradient = gradient_cls(sampler)
+
+        job = gradient.run([circuit], [values])
+        result = job.result()
+
+        with self.subTest(msg="assert initialize is preserved"):
+            self.assertTrue(all("initialize" in ops_i[0].keys() for ops_i in ops))
+
+        with self.subTest(msg="assert result is correct"):
+            array1 = _quasi2array(result.gradients[0], num_qubits=2)
+            array2 = _quasi2array(expect, num_qubits=2)
+            np.testing.assert_allclose(array1, array2, atol=1e-5)
+
+
+@ddt
+class TestSamplerGradientV2(QiskitAlgorithmsTestCase):
+    """Test Sampler Gradient"""
+
+    def __init__(self, TestCase):
+        backend = GenericBackendV2(num_qubits=3, seed=123)
+        session = Session(backend=backend)
+        self.sampler = SamplerV2(mode=session)
+        self.pass_manager = generate_preset_pass_manager(optimization_level=1, backend=backend)
+        super().__init__(TestCase)
+
+    @data(*gradient_factories)
+    @unittest.skip("Skipping due to noise sensitivity.")
+    def test_single_circuit(self, grad):
+        """Test the sampler gradient for a single circuit"""
         a = Parameter("a")
         qc = QuantumCircuit(1)
-        qc.rx(a, 0)
+        qc.h(0)
+        qc.p(a, 0)
+        qc.h(0)
+        qc.measure_all()
+
+        gradient = grad(sampler=self.sampler, pass_manager=self.pass_manager)
+        param_list = [[np.pi / 4], [0], [np.pi / 2]]
+        expected = [
+            [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
+            [{0: 0, 1: 0}],
+            [{0: -0.499999, 1: 0.499999}],
+        ]
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [param]).result().gradients[0]
+            array1 = _quasi2array(gradients, num_qubits=1)
+            array2 = _quasi2array(expected[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=0.5, rtol=0.5)
+
+    @data(*gradient_factories)
+    @unittest.skip("Skipping due to noise sensitivity.")
+    def test_gradient_p(self, grad):
+        """Test the sampler gradient for p"""
+
+        a = Parameter("a")
+        qc = QuantumCircuit(1)
+        qc.h(0)
+        qc.p(a, 0)
+        qc.h(0)
+        qc.measure_all()
+
+        gradient = grad(
+            sampler=self.sampler,
+            pass_manager=self.pass_manager,
+        )
+        param_list = [[np.pi / 4], [0], [np.pi / 2]]
+        expected = [
+            [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
+            [{0: 0, 1: 0}],
+            [{0: -0.499999, 1: 0.499999}],
+        ]
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [param]).result().gradients[0]
+            array1 = _quasi2array(gradients, num_qubits=1)
+            array2 = _quasi2array(expected[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    @unittest.skip("Skipping due to noise sensitivity.")
+    def test_gradient_u(self, grad):
+        """Test the sampler gradient for u"""
+
+        a = Parameter("a")
+        b = Parameter("b")
+        c = Parameter("c")
+        qc = QuantumCircuit(1)
+        qc.h(0)
+        qc.u(a, b, c, 0)
+        qc.h(0)
+        qc.measure_all()
+        gradient = grad(
+            sampler=self.sampler,
+            pass_manager=self.pass_manager,
+        )
+        param_list = [[np.pi / 4, 0, 0], [np.pi / 4, np.pi / 4, np.pi / 4]]
+        expected = [
+            [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}, {0: 0, 1: 0}, {0: 0, 1: 0}],
+            [{0: -0.176777, 1: 0.176777}, {0: -0.426777, 1: 0.426777}, {0: -0.426777, 1: 0.426777}],
+        ]
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [param]).result().gradients[0]
+            array1 = _quasi2array(gradients, num_qubits=1)
+            array2 = _quasi2array(expected[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    @unittest.skip("Skipping due to noise sensitivity.")
+    def test_gradient_efficient_su2(self, grad):
+        """Test the sampler gradient for EfficientSU2"""
+
+        qc = EfficientSU2(2, reps=1)
+        qc.measure_all()
+        gradient = grad(
+            sampler=self.sampler,
+            pass_manager=self.pass_manager,
+        )
+        param_list = [
+            [np.pi / 4 for param in qc.parameters],
+            [np.pi / 2 for param in qc.parameters],
+        ]
+        expected = [
+            [
+                {
+                    0: -0.11963834764831836,
+                    1: -0.05713834764831845,
+                    2: -0.21875000000000003,
+                    3: 0.39552669529663675,
+                },
+                {
+                    0: -0.32230339059327373,
+                    1: -0.031250000000000014,
+                    2: 0.2339150429449554,
+                    3: 0.11963834764831843,
+                },
+                {
+                    0: 0.012944173824159189,
+                    1: -0.01294417382415923,
+                    2: 0.07544417382415919,
+                    3: -0.07544417382415919,
+                },
+                {
+                    0: 0.2080266952966367,
+                    1: -0.03125000000000002,
+                    2: -0.11963834764831842,
+                    3: -0.057138347648318405,
+                },
+                {
+                    0: -0.11963834764831838,
+                    1: 0.11963834764831838,
+                    2: -0.21875000000000003,
+                    3: 0.21875,
+                },
+                {
+                    0: -0.2781092167691146,
+                    1: -0.0754441738241592,
+                    2: 0.27810921676911443,
+                    3: 0.07544417382415924,
+                },
+                {0: 0.0, 1: 0.0, 2: 0.0, 3: 0.0},
+                {0: 0.0, 1: 0.0, 2: 0.0, 3: 0.0},
+            ],
+            [
+                {
+                    0: -4.163336342344337e-17,
+                    1: 2.7755575615628914e-17,
+                    2: -4.163336342344337e-17,
+                    3: 0.0,
+                },
+                {0: 0.0, 1: -1.3877787807814457e-17, 2: 4.163336342344337e-17, 3: 0.0},
+                {
+                    0: -0.24999999999999994,
+                    1: 0.24999999999999994,
+                    2: 0.24999999999999994,
+                    3: -0.24999999999999994,
+                },
+                {
+                    0: 0.24999999999999994,
+                    1: 0.24999999999999994,
+                    2: -0.24999999999999994,
+                    3: -0.24999999999999994,
+                },
+                {
+                    0: -4.163336342344337e-17,
+                    1: 4.163336342344337e-17,
+                    2: -4.163336342344337e-17,
+                    3: 5.551115123125783e-17,
+                },
+                {
+                    0: -0.24999999999999994,
+                    1: 0.24999999999999994,
+                    2: 0.24999999999999994,
+                    3: -0.24999999999999994,
+                },
+                {0: 0.0, 1: 2.7755575615628914e-17, 2: 0.0, 3: 2.7755575615628914e-17},
+                {0: 0.0, 1: 0.0, 2: 0.0, 3: 0.0},
+            ],
+        ]
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [param]).result().gradients[0]
+            array1 = _quasi2array(gradients, num_qubits=2)
+            array2 = _quasi2array(expected[i], num_qubits=2)
+            np.testing.assert_allclose(array1, array2, atol=0.5, rtol=0.5)
+
+    @data(*gradient_factories)
+    @unittest.skip("Skipping due to noise sensitivity.")
+    def test_gradient_2qubit_gate(self, grad):
+        """Test the sampler gradient for 2 qubit gates"""
+
+        for gate in [RXXGate]:
+            param_list = [[np.pi / 4], [np.pi / 2]]
+            correct_results = [
+                [{0: -0.5 / np.sqrt(2), 1: 0, 2: 0, 3: 0.5 / np.sqrt(2)}],
+                [{0: -0.5, 1: 0, 2: 0, 3: 0.5}],
+            ]
+            for i, param in enumerate(param_list):
+                a = Parameter("a")
+                qc = QuantumCircuit(2)
+                qc.append(gate(a), [qc.qubits[0], qc.qubits[1]], [])
+                qc.measure_all()
+                gradient = grad(sampler=self.sampler, pass_manager=self.pass_manager)
+                gradients = gradient.run([qc], [param]).result().gradients[0]
+                array1 = _quasi2array(gradients, num_qubits=2)
+                array2 = _quasi2array(correct_results[i], num_qubits=2)
+                np.testing.assert_allclose(array1, array2, atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    @unittest.skip("Skipping due to noise sensitivity.")
+    def test_gradient_parameter_coefficient(self, grad):
+        """Test the sampler gradient for parameter variables with coefficients"""
+
+        qc = RealAmplitudes(num_qubits=2, reps=1)
+        qc.rz(qc.parameters[0].exp() + 2 * qc.parameters[1], 0)
+        qc.rx(3.0 * qc.parameters[0] + qc.parameters[1].sin(), 1)
+        qc.u(qc.parameters[0], qc.parameters[1], qc.parameters[3], 1)
+        qc.p(2 * qc.parameters[0] + 1, 0)
+        qc.rxx(qc.parameters[0] + 2, 0, 1)
         qc.measure_all()
-        sampler = Sampler(options={"shots": 100})
-        with self.subTest("sampler"):
-            if grad is SPSASamplerGradient:
-                gradient = grad(sampler, epsilon=1e-6)
-            else:
-                gradient = grad(sampler)
-            options = gradient.options
-            result = gradient.run([qc], [[1]]).result()
-            self.assertEqual(result.options.get("shots"), 100)
-            self.assertEqual(options.get("shots"), 100)
-
-        with self.subTest("gradient init"):
-            if grad is SPSASamplerGradient:
-                gradient = grad(sampler, epsilon=1e-6, options={"shots": 200})
-            else:
-                gradient = grad(sampler, options={"shots": 200})
-            options = gradient.options
-            result = gradient.run([qc], [[1]]).result()
-            self.assertEqual(result.options.get("shots"), 200)
-            self.assertEqual(options.get("shots"), 200)
-
-        with self.subTest("gradient update"):
-            if grad is SPSASamplerGradient:
-                gradient = grad(sampler, epsilon=1e-6, options={"shots": 200})
-            else:
-                gradient = grad(sampler, options={"shots": 200})
-            gradient.update_default_options(shots=100)
-            options = gradient.options
-            result = gradient.run([qc], [[1]]).result()
-            self.assertEqual(result.options.get("shots"), 100)
-            self.assertEqual(options.get("shots"), 100)
-
-        with self.subTest("gradient run"):
-            if grad is SPSASamplerGradient:
-                gradient = grad(sampler, epsilon=1e-6, options={"shots": 200})
-            else:
-                gradient = grad(sampler, options={"shots": 200})
-            options = gradient.options
-            result = gradient.run([qc], [[1]], shots=300).result()
-            self.assertEqual(result.options.get("shots"), 300)
-            # Only default + sampler options. Not run.
-            self.assertEqual(options.get("shots"), 200)
+
+        gradient = grad(
+            sampler=self.sampler,
+            pass_manager=self.pass_manager,
+        )
+        param_list = [[np.pi / 4 for _ in qc.parameters], [np.pi / 2 for _ in qc.parameters]]
+        correct_results = [
+            [
+                {
+                    0: 0.30014831912265927,
+                    1: -0.6634809704357856,
+                    2: 0.343589357193753,
+                    3: 0.019743294119373426,
+                },
+                {
+                    0: 0.16470607453981906,
+                    1: -0.40996282450610577,
+                    2: 0.08791803062881773,
+                    3: 0.15733871933746948,
+                },
+                {
+                    0: 0.27036068339663866,
+                    1: -0.273790986018701,
+                    2: 0.12752010079553433,
+                    3: -0.12408979817347202,
+                },
+                {
+                    0: -0.2098616294167757,
+                    1: -0.2515823946449894,
+                    2: 0.21929102305386305,
+                    3: 0.24215300100790207,
+                },
+            ],
+            [
+                {
+                    0: -1.844810060881004,
+                    1: 0.04620532700836027,
+                    2: 1.6367366426074323,
+                    3: 0.16186809126521057,
+                },
+                {
+                    0: 0.07296073407769421,
+                    1: -0.021774869186331716,
+                    2: 0.02177486918633173,
+                    3: -0.07296073407769456,
+                },
+                {
+                    0: -0.07794369186049102,
+                    1: -0.07794369186049122,
+                    2: 0.07794369186049117,
+                    3: 0.07794369186049112,
+                },
+                {
+                    0: 0.0,
+                    1: 0.0,
+                    2: 0.0,
+                    3: 0.0,
+                },
+            ],
+        ]
+
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [param]).result().gradients[0]
+            array1 = _quasi2array(gradients, num_qubits=2)
+            array2 = _quasi2array(correct_results[i], num_qubits=2)
+            np.testing.assert_allclose(array1, array2, atol=0.5, rtol=0.5)
+
+    @data(*gradient_factories)
+    @unittest.skip("Skipping due to noise sensitivity.")
+    def test_gradient_parameters(self, grad):
+        """Test the sampler gradient for parameters"""
+
+        a = Parameter("a")
+        b = Parameter("b")
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        qc.rz(b, 0)
+        qc.measure_all()
+
+        gradient = grad(
+            sampler=self.sampler,
+            pass_manager=self.pass_manager,
+        )
+        param_list = [[np.pi / 4, np.pi / 2]]
+        expected = [
+            [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
+        ]
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [param], parameters=[[a]]).result().gradients[0]
+            array1 = _quasi2array(gradients, num_qubits=1)
+            array2 = _quasi2array(expected[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=0.5, rtol=0.5)
+
+        # parameter order
+        with self.subTest(msg="The order of gradients"):
+            c = Parameter("c")
+            qc = QuantumCircuit(1)
+            qc.rx(a, 0)
+            qc.rz(b, 0)
+            qc.rx(c, 0)
+            qc.measure_all()
+
+            param_values = [[np.pi / 4, np.pi / 2, np.pi / 3]]
+            params = [[a, b, c], [c, b, a], [a, c], [c, a]]
+            expected = [
+                [
+                    {0: -0.17677666583387008, 1: 0.17677666583378482},
+                    {0: 0.3061861668168149, 1: -0.3061861668167012},
+                    {0: -0.3061861668168149, 1: 0.30618616681678645},
+                ],
+                [
+                    {0: -0.3061861668168149, 1: 0.30618616681678645},
+                    {0: 0.3061861668168149, 1: -0.3061861668167012},
+                    {0: -0.17677666583387008, 1: 0.17677666583378482},
+                ],
+                [
+                    {0: -0.17677666583387008, 1: 0.17677666583378482},
+                    {0: -0.3061861668168149, 1: 0.30618616681678645},
+                ],
+                [
+                    {0: -0.3061861668168149, 1: 0.30618616681678645},
+                    {0: -0.17677666583387008, 1: 0.17677666583378482},
+                ],
+            ]
+            for i, p in enumerate(params):  # pylint: disable=invalid-name
+                gradients = gradient.run([qc], param_values, parameters=[p]).result().gradients[0]
+                array1 = _quasi2array(gradients, num_qubits=1)
+                array2 = _quasi2array(expected[i], num_qubits=1)
+                np.testing.assert_allclose(array1, array2, atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    @unittest.skip("Skipping due to noise sensitivity.")
+    def test_gradient_multi_arguments(self, grad):
+        """Test the sampler gradient for multiple arguments"""
+
+        a = Parameter("a")
+        b = Parameter("b")
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        qc.measure_all()
+
+        qc2 = QuantumCircuit(1)
+        qc2.rx(b, 0)
+        qc2.measure_all()
+        gradient = grad(
+            sampler=self.sampler,
+            pass_manager=self.pass_manager,
+        )
+        param_list = [[np.pi / 4], [np.pi / 2]]
+        correct_results = [
+            [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
+            [{0: -0.499999, 1: 0.499999}],
+        ]
+        gradients = gradient.run([qc, qc2], param_list).result().gradients
+        for i, q_dists in enumerate(gradients):
+            array1 = _quasi2array(q_dists, num_qubits=1)
+            array2 = _quasi2array(correct_results[i], num_qubits=1)
+            np.testing.assert_allclose(array1, array2, atol=1e-1, rtol=1e-1)
+
+        # parameters
+        with self.subTest(msg="Different parameters"):
+            c = Parameter("c")
+            qc3 = QuantumCircuit(1)
+            qc3.rx(c, 0)
+            qc3.ry(a, 0)
+            qc3.measure_all()
+            param_list2 = [[np.pi / 4], [np.pi / 4, np.pi / 4], [np.pi / 4, np.pi / 4]]
+            gradients = (
+                gradient.run([qc, qc3, qc3], param_list2, parameters=[[a], [c], None])
+                .result()
+                .gradients
+            )
+            correct_results = [
+                [{0: -0.5 / np.sqrt(2), 1: 0.5 / np.sqrt(2)}],
+                [{0: -0.25, 1: 0.25}],
+                [{0: -0.25, 1: 0.25}, {0: -0.25, 1: 0.25}],
+            ]
+            for i, q_dists in enumerate(gradients):
+                array1 = _quasi2array(q_dists, num_qubits=1)
+                array2 = _quasi2array(correct_results[i], num_qubits=1)
+                np.testing.assert_allclose(array1, array2, atol=1e-1, rtol=1e-1)
+
+    @data(*gradient_factories)
+    def test_gradient_validation(self, grad):
+        """Test sampler gradient's validation"""
+
+        a = Parameter("a")
+        qc = QuantumCircuit(1)
+        qc.rx(a, 0)
+        qc.measure_all()
+
+        gradient = grad(
+            sampler=self.sampler,
+            pass_manager=self.pass_manager,
+        )
+        param_list = [[np.pi / 4], [np.pi / 2]]
+        with self.assertRaises(ValueError):
+            gradient.run([qc], param_list)
+        with self.assertRaises(ValueError):
+            gradient.run([qc, qc], param_list, parameters=[[a]])
+        with self.assertRaises(ValueError):
+            gradient.run([qc], [[np.pi / 4, np.pi / 4]])
+
+    @unittest.skip("Skipping due to noise sensitivity.")
+    def test_spsa_gradient(self):
+        """Test the SPSA sampler gradient"""
+
+        with self.assertRaises(ValueError):
+            _ = SPSASamplerGradient(self.sampler, epsilon=-0.01)
+
+        a = Parameter("a")
+        b = Parameter("b")
+        c = Parameter("c")
+        qc = QuantumCircuit(2)
+        qc.rx(b, 0)
+        qc.rx(a, 1)
+        qc.measure_all()
+        param_list = [[1, 2]]
+        correct_results = [
+            [
+                {0: 0.2273244, 1: -0.6480598, 2: 0.2273244, 3: 0.1934111},
+                {0: -0.2273244, 1: 0.6480598, 2: -0.2273244, 3: -0.1934111},
+            ],
+        ]
+        gradient = SPSASamplerGradient(
+            sampler=self.sampler, pass_manager=self.pass_manager, epsilon=1e-6, seed=123
+        )
+        for i, param in enumerate(param_list):
+            gradients = gradient.run([qc], [param]).result().gradients[0]
+            array1 = _quasi2array(gradients, num_qubits=2)
+            array2 = _quasi2array(correct_results[i], num_qubits=2)
+            np.testing.assert_allclose(array1, array2, atol=1e-3)
+
+        # multi parameters
+        with self.subTest(msg="Multiple parameters"):
+            param_list2 = [[1, 2], [1, 2], [3, 4]]
+            correct_results2 = [
+                [
+                    {0: 0.2273244, 1: -0.6480598, 2: 0.2273244, 3: 0.1934111},
+                    {0: -0.2273244, 1: 0.6480598, 2: -0.2273244, 3: -0.1934111},
+                ],
+                [
+                    {0: -0.2273244, 1: 0.6480598, 2: -0.2273244, 3: -0.1934111},
+                ],
+                [
+                    {0: -0.0141129, 1: -0.0564471, 2: -0.3642884, 3: 0.4348484},
+                    {0: 0.0141129, 1: 0.0564471, 2: 0.3642884, 3: -0.4348484},
+                ],
+            ]
+            gradient = SPSASamplerGradient(self.sampler, epsilon=1e-6, seed=123)
+            gradients = (
+                gradient.run([qc] * 3, param_list2, parameters=[None, [b], None]).result().gradients
+            )
+            for i, result in enumerate(gradients):
+                array1 = _quasi2array(result, num_qubits=2)
+                array2 = _quasi2array(correct_results2[i], num_qubits=2)
+                np.testing.assert_allclose(array1, array2, atol=1e-3)
+
+        # batch size
+        with self.subTest(msg="Batch size"):
+            param_list = [[1, 1]]
+            gradient = SPSASamplerGradient(self.sampler, epsilon=1e-6, batch_size=4, seed=123)
+            gradients = gradient.run([qc], param_list).result().gradients
+            correct_results3 = [
+                [
+                    {
+                        0: -0.1620149622932887,
+                        1: -0.25872053011771756,
+                        2: 0.3723827084675668,
+                        3: 0.04835278392088804,
+                    },
+                    {
+                        0: -0.1620149622932887,
+                        1: 0.3723827084675668,
+                        2: -0.25872053011771756,
+                        3: 0.04835278392088804,
+                    },
+                ]
+            ]
+            for i, q_dists in enumerate(gradients):
+                array1 = _quasi2array(q_dists, num_qubits=2)
+                array2 = _quasi2array(correct_results3[i], num_qubits=2)
+                np.testing.assert_allclose(array1, array2, atol=1e-3)
+
+        # parameter order
+        with self.subTest(msg="The order of gradients"):
+            qc = QuantumCircuit(1)
+            qc.rx(a, 0)
+            qc.rz(b, 0)
+            qc.rx(c, 0)
+            qc.measure_all()
+            param_list = [[np.pi / 4, np.pi / 2, np.pi / 3]]
+            param = [[a, b, c], [c, b, a], [a, c], [c, a]]
+            correct_results = [
+                [
+                    {0: -0.17677624757590138, 1: 0.17677624757590138},
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                ],
+                [
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                    {0: -0.17677624757590138, 1: 0.17677624757590138},
+                ],
+                [
+                    {0: -0.17677624757590138, 1: 0.17677624757590138},
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                ],
+                [
+                    {0: 0.17677624757590138, 1: -0.17677624757590138},
+                    {0: -0.17677624757590138, 1: 0.17677624757590138},
+                ],
+            ]
+            for i, p in enumerate(param):  # pylint: disable=invalid-name
+                gradient = SPSASamplerGradient(self.sampler, epsilon=1e-6, seed=123)
+                gradients = gradient.run([qc], param_list, parameters=[p]).result().gradients[0]
+                array1 = _quasi2array(gradients, num_qubits=1)
+                array2 = _quasi2array(correct_results[i], num_qubits=1)
+                np.testing.assert_allclose(array1, array2, atol=1e-3)
 
     @data(
         ParamShiftSamplerGradient,
diff --git a/test/neural_networks/test_estimator_qnn_v2.py b/test/neural_networks/test_estimator_qnn_v2.py
index b8fad6557..f5c3539aa 100644
--- a/test/neural_networks/test_estimator_qnn_v2.py
+++ b/test/neural_networks/test_estimator_qnn_v2.py
@@ -191,8 +191,10 @@ def __init__(
         TestCase,
     ):
         self.estimator = EstimatorV2(mode=self.session, options={"default_shots": 1e3})
-        self.pm = generate_preset_pass_manager(backend=self.backend, optimization_level=0)
-        self.gradient = ParamShiftEstimatorGradient(estimator=self.estimator, pass_manager=self.pm)
+        self.pass_manager = generate_preset_pass_manager(backend=self.backend, optimization_level=0)
+        self.gradient = ParamShiftEstimatorGradient(
+            estimator=self.estimator, pass_manager=self.pass_manager
+        )
         super().__init__(TestCase)
 
     def _test_network_passes(
@@ -247,7 +249,7 @@ def test_estimator_qnn_1_1(self):
         qc.h(0)
         qc.ry(params[0], 0)
         qc.rx(params[1], 0)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op = SparsePauliOp.from_list([("Z", 1), ("X", 1)])
         isa_ob = op.apply_layout(isa_qc.layout)
         estimator_qnn = EstimatorQNN(
@@ -257,7 +259,6 @@ def test_estimator_qnn_1_1(self):
             weight_params=[params[1]],
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
 
         self._test_network_passes(estimator_qnn, CASE_DATA["shape_1_1"])
@@ -276,7 +277,7 @@ def test_estimator_qnn_2_1(self):
         qc.ry(params[1], 1)
         qc.rx(params[2], 0)
         qc.rx(params[3], 1)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op = SparsePauliOp.from_list([("ZZ", 1), ("XX", 1)])
         op = op.apply_layout(isa_qc.layout)
         estimator_qnn = EstimatorQNN(
@@ -286,7 +287,6 @@ def test_estimator_qnn_2_1(self):
             weight_params=params[2:],
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
 
         self._test_network_passes(estimator_qnn, CASE_DATA["shape_2_1"])
@@ -299,7 +299,7 @@ def test_estimator_qnn_1_2(self):
         qc.ry(params[0], 0)
         qc.rx(params[1], 0)
 
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op1 = SparsePauliOp.from_list([("Z", 1), ("X", 1)])
         op1 = op1.apply_layout(isa_qc.layout)
         op2 = SparsePauliOp.from_list([("Z", 2), ("X", 2)])
@@ -313,7 +313,6 @@ def test_estimator_qnn_1_2(self):
             weight_params=[params[1]],
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
 
         self._test_network_passes(estimator_qnn, CASE_DATA["shape_1_2"])
@@ -332,7 +331,7 @@ def test_estimator_qnn_2_2(self):
         qc.ry(params[1], 1)
         qc.rx(params[2], 0)
         qc.rx(params[3], 1)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op1 = SparsePauliOp.from_list([("ZZ", 1)])
         op1 = op1.apply_layout(isa_qc.layout)
         op2 = SparsePauliOp.from_list([("XX", 1)])
@@ -345,7 +344,6 @@ def test_estimator_qnn_2_2(self):
             weight_params=params[2:],
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
 
         self._test_network_passes(estimator_qnn, CASE_DATA["shape_2_2"])
@@ -357,7 +355,7 @@ def test_no_input_parameters(self):
         qc.h(0)
         qc.ry(params[0], 0)
         qc.rx(params[1], 0)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op = SparsePauliOp.from_list([("Z", 1), ("X", 1)])
         op = op.apply_layout(isa_qc.layout)
         estimator_qnn = EstimatorQNN(
@@ -367,7 +365,6 @@ def test_no_input_parameters(self):
             weight_params=params,
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
         self._test_network_passes(estimator_qnn, CASE_DATA["no_input_parameters"])
 
@@ -378,7 +375,7 @@ def test_no_weight_parameters(self):
         qc.h(0)
         qc.ry(params[0], 0)
         qc.rx(params[1], 0)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op = SparsePauliOp.from_list([("Z", 1), ("X", 1)])
         op = op.apply_layout(isa_qc.layout)
         estimator_qnn = EstimatorQNN(
@@ -388,7 +385,6 @@ def test_no_weight_parameters(self):
             weight_params=None,
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
         self._test_network_passes(estimator_qnn, CASE_DATA["no_weight_parameters"])
 
@@ -396,7 +392,7 @@ def test_no_parameters(self):
         """Test Estimator QNN with no parameters."""
         qc = QuantumCircuit(1)
         qc.h(0)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op = SparsePauliOp.from_list([("Z", 1), ("X", 1)])
         op = op.apply_layout(isa_qc.layout)
         estimator_qnn = EstimatorQNN(
@@ -406,7 +402,6 @@ def test_no_parameters(self):
             weight_params=None,
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
         self._test_network_passes(estimator_qnn, CASE_DATA["no_parameters"])
 
@@ -417,14 +412,13 @@ def test_default_observables(self):
         qc.h(0)
         qc.ry(params[0], 0)
         qc.rx(params[1], 0)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         estimator_qnn = EstimatorQNN(
             circuit=isa_qc,
             input_params=[params[0]],
             weight_params=[params[1]],
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
         self._test_network_passes(estimator_qnn, CASE_DATA["default_observables"])
 
@@ -435,7 +429,7 @@ def test_single_observable(self):
         qc.h(0)
         qc.ry(params[0], 0)
         qc.rx(params[1], 0)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op = SparsePauliOp.from_list([("Z", 1), ("X", 1)])
         op = op.apply_layout(isa_qc.layout)
         estimator_qnn = EstimatorQNN(
@@ -445,7 +439,6 @@ def test_single_observable(self):
             weight_params=[params[1]],
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=isa_qc.num_qubits,
         )
         self._test_network_passes(estimator_qnn, CASE_DATA["single_observable"])
 
@@ -456,7 +449,7 @@ def test_setters_getters(self):
         qc.h(0)
         qc.ry(params[0], 0)
         qc.rx(params[1], 0)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op = SparsePauliOp.from_list([("Z", 1), ("X", 1)])
         op = op.apply_layout(isa_qc.layout)
         estimator_qnn = EstimatorQNN(
@@ -466,7 +459,6 @@ def test_setters_getters(self):
             weight_params=[params[1]],
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
         with self.subTest("Test circuit getter."):
             self.assertEqual(estimator_qnn.circuit, isa_qc)
@@ -491,7 +483,7 @@ def test_qnn_qc_circuit_construction(self):
         qc = QuantumCircuit(num_qubits)
         qc.compose(feature_map, inplace=True)
         qc.compose(ansatz, inplace=True)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
 
         estimator_qc = EstimatorQNN(
             circuit=isa_qc,
@@ -500,11 +492,10 @@ def test_qnn_qc_circuit_construction(self):
             input_gradients=True,
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
 
         qnn_qc = QNNCircuit(num_qubits=num_qubits, feature_map=feature_map, ansatz=ansatz)
-        isa_qnn_qc = self.pm.run(qnn_qc)
+        isa_qnn_qc = self.pass_manager.run(qnn_qc)
         estimator_qnn_qc = EstimatorQNN(
             circuit=isa_qnn_qc,
             input_params=qnn_qc.feature_map.parameters,
@@ -512,7 +503,6 @@ def test_qnn_qc_circuit_construction(self):
             input_gradients=True,
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
 
         input_data = [1, 2]
@@ -544,7 +534,7 @@ def test_binding_order(self):
         weight = Parameter("weight")
         for i in range(qc.num_qubits):
             qc.rx(weight, i)
-        isa_qc = self.pm.run(qc)
+        isa_qc = self.pass_manager.run(qc)
         op = SparsePauliOp.from_list([("Z" * isa_qc.num_qubits, 1)])
         op = op.apply_layout(isa_qc.layout)
         estimator_qnn = EstimatorQNN(
@@ -554,7 +544,6 @@ def test_binding_order(self):
             weight_params=[weight],
             estimator=self.estimator,
             gradient=self.gradient,
-            num_virtual_qubits=qc.num_qubits,
         )
 
         estimator_qnn_weights = [3]
@@ -562,7 +551,7 @@ def test_binding_order(self):
         res = estimator_qnn.forward(estimator_qnn_input, estimator_qnn_weights)
         # When parameters were used in circuit order, before being assigned correctly, so inputs
         # went to input params, weights to weight params, this gave 0.00613403
-        self.assertAlmostEqual(res[0][0], 0.00040017, delta=0.05)
+        self.assertAlmostEqual(res[0][0], 0.00040017, delta=0.1)
 
 
 if __name__ == "__main__":
diff --git a/test/neural_networks/test_sampler_qnn.py b/test/neural_networks/test_sampler_qnn.py
index 9651a93d4..e2b821cc8 100644
--- a/test/neural_networks/test_sampler_qnn.py
+++ b/test/neural_networks/test_sampler_qnn.py
@@ -78,7 +78,6 @@ def setUp(self):
         self.qc.append(feature_map, range(2))
         self.qc.append(var_form, range(2))
         self.qc.measure_all()
-        self.num_virtual_qubits = num_qubits
 
         # store params
         self.input_params = list(feature_map.parameters)
@@ -106,7 +105,7 @@ def interpret_2d(x):
         self.backend = GenericBackendV2(num_qubits=8)
         self.session = Session(backend=self.backend)
         self.sampler_v2 = SamplerV2(mode=self.session)
-        self.pm = None
+        self.pass_manager = None
         self.array_type = {True: SparseArray, False: np.ndarray}
 
     # pylint: disable=too-many-positional-arguments
@@ -134,12 +133,13 @@ def _get_qnn(
             sampler = self.sampler_v2
 
             if self.qc.layout is None:
-                self.pm = generate_preset_pass_manager(optimization_level=1, backend=self.backend)
-                self.qc = self.pm.run(self.qc)
+                self.pass_manager = generate_preset_pass_manager(
+                    optimization_level=1, backend=self.backend
+                )
+                self.qc = self.pass_manager.run(self.qc)
             gradient = ParamShiftSamplerGradient(
                 sampler=self.sampler,
-                len_quasi_dist=2**self.num_virtual_qubits,
-                pass_manager=self.pm,
+                pass_manager=self.pass_manager,
             )
         else:
             sampler = None
@@ -158,7 +158,6 @@ def _get_qnn(
         qnn = SamplerQNN(
             sampler=sampler,
             circuit=self.qc,
-            num_virtual_qubits=self.num_virtual_qubits,
             input_params=input_params,
             weight_params=weight_params,
             sparse=sparse,
diff --git a/test/state_fidelities/test_compute_uncompute_v2.py b/test/state_fidelities/test_compute_uncompute_v2.py
index 819b206fc..30d276333 100644
--- a/test/state_fidelities/test_compute_uncompute_v2.py
+++ b/test/state_fidelities/test_compute_uncompute_v2.py
@@ -63,7 +63,7 @@ def setUp(self):
         )
         self.session = Session(backend=self.backend)
         self._sampler = SamplerV2(mode=self.session)
-        self.pm = generate_preset_pass_manager(optimization_level=0, backend=self.backend)
+        self.pass_manager = generate_preset_pass_manager(optimization_level=0, backend=self.backend)
 
         self._left_params = np.array([[0, 0], [np.pi / 2, 0], [0, np.pi / 2], [np.pi, np.pi]])
         self._right_params = np.array([[0, 0], [0, 0], [np.pi / 2, 0], [0, 0]])
@@ -71,7 +71,8 @@ def setUp(self):
     def test_1param_pair(self):
         """test for fidelity with one pair of parameters"""
         fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[0].num_qubits
+            self._sampler,
+            pass_manager=self.pass_manager,
         )
         job = fidelity.run(
             self._circuit[0], self._circuit[1], self._left_params[0], self._right_params[0]
@@ -84,8 +85,7 @@ def test_1param_pair_local(self):
         fidelity = ComputeUncompute(
             self._sampler,
             local=True,
-            pass_manager=self.pm,
-            num_virtual_qubits=self._circuit[0].num_qubits,
+            pass_manager=self.pass_manager,
         )
         job = fidelity.run(
             self._circuit[0], self._circuit[1], self._left_params[0], self._right_params[0]
@@ -98,14 +98,12 @@ def test_local(self):
         fidelity_global = ComputeUncompute(
             self._sampler,
             local=False,
-            pass_manager=self.pm,
-            num_virtual_qubits=self._circuit[2].num_qubits,
+            pass_manager=self.pass_manager,
         )
         fidelity_local = ComputeUncompute(
             self._sampler,
             local=True,
-            pass_manager=self.pm,
-            num_virtual_qubits=self._circuit[2].num_qubits,
+            pass_manager=self.pass_manager,
         )
         fidelities = []
         for fidelity in [fidelity_global, fidelity_local]:
@@ -116,9 +114,7 @@ def test_local(self):
 
     def test_4param_pairs(self):
         """test for fidelity with four pairs of parameters"""
-        fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[0].num_qubits
-        )
+        fidelity = ComputeUncompute(self._sampler, pass_manager=self.pass_manager)
         n = len(self._left_params)
         job = fidelity.run(
             [self._circuit[0]] * n, [self._circuit[1]] * n, self._left_params, self._right_params
@@ -130,9 +126,7 @@ def test_4param_pairs(self):
 
     def test_symmetry(self):
         """test for fidelity with the same circuit"""
-        fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[0].num_qubits
-        )
+        fidelity = ComputeUncompute(self._sampler, pass_manager=self.pass_manager)
         n = len(self._left_params)
         job_1 = fidelity.run(
             [self._circuit[0]] * n, [self._circuit[0]] * n, self._left_params, self._right_params
@@ -148,7 +142,8 @@ def test_symmetry(self):
     def test_no_params(self):
         """test for fidelity without parameters"""
         fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[2].num_qubits
+            self._sampler,
+            pass_manager=self.pass_manager,
         )
         job = fidelity.run([self._circuit[2]], [self._circuit[3]])
         results = job.result()
@@ -161,7 +156,8 @@ def test_no_params(self):
     def test_left_param(self):
         """test for fidelity with only left parameters"""
         fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[1].num_qubits
+            self._sampler,
+            pass_manager=self.pass_manager,
         )
         n = len(self._left_params)
         job = fidelity.run(
@@ -175,7 +171,8 @@ def test_left_param(self):
     def test_right_param(self):
         """test for fidelity with only right parameters"""
         fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[1].num_qubits
+            self._sampler,
+            pass_manager=self.pass_manager,
         )
         n = len(self._left_params)
         job = fidelity.run(
@@ -188,9 +185,7 @@ def test_right_param(self):
 
     def test_not_set_circuits(self):
         """test for fidelity with no circuits."""
-        fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[0].num_qubits
-        )
+        fidelity = ComputeUncompute(self._sampler, pass_manager=self.pass_manager)
         with self.assertRaises(TypeError):
             job = fidelity.run(
                 circuits_1=None,
@@ -202,9 +197,7 @@ def test_not_set_circuits(self):
 
     def test_circuit_mismatch(self):
         """test for fidelity with different number of left/right circuits."""
-        fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[0].num_qubits
-        )
+        fidelity = ComputeUncompute(self._sampler, pass_manager=self.pass_manager)
         n = len(self._left_params)
         with self.assertRaises(ValueError):
             job = fidelity.run(
@@ -219,9 +212,7 @@ def test_asymmetric_params(self):
         """test for fidelity when the 2 circuits have different number of
         left/right parameters."""
 
-        fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[0].num_qubits
-        )
+        fidelity = ComputeUncompute(self._sampler, pass_manager=self.pass_manager)
         n = len(self._left_params)
         right_params = [[p] for p in self._right_params[:, 0]]
         job = fidelity.run(
@@ -236,9 +227,7 @@ def test_input_format(self):
         """test for different input format variations"""
 
         circuit = RealAmplitudes(2)
-        fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=circuit.num_qubits
-        )
+        fidelity = ComputeUncompute(self._sampler, pass_manager=self.pass_manager)
         values = np.random.random(circuit.num_parameters)
         shift = np.ones_like(values) * 0.01
 
@@ -266,9 +255,7 @@ def test_input_format(self):
 
     def test_input_measurements(self):
         """test for fidelity with measurements on input circuits"""
-        fidelity = ComputeUncompute(
-            self._sampler, pass_manager=self.pm, num_virtual_qubits=self._circuit[0].num_qubits
-        )
+        fidelity = ComputeUncompute(self._sampler, pass_manager=self.pass_manager)
         circuit_1 = self._circuit[0]
         circuit_1.measure_all()
         circuit_2 = self._circuit[1]
@@ -284,9 +271,7 @@ def test_options(self):
 
         with self.subTest("sampler"):
             # Only options in sampler
-            fidelity = ComputeUncompute(
-                sampler_shots, pass_manager=self.pm, num_virtual_qubits=self._circuit[2].num_qubits
-            )
+            fidelity = ComputeUncompute(sampler_shots, pass_manager=self.pass_manager)
             options = fidelity.options
             job = fidelity.run(self._circuit[2], self._circuit[3])
             result = job.result()
@@ -299,8 +284,7 @@ def test_options(self):
             fidelity = ComputeUncompute(
                 sampler_shots,
                 options={"shots": 2048, "dummy": 100},
-                pass_manager=self.pm,
-                num_virtual_qubits=self._circuit[2].num_qubits,
+                pass_manager=self.pass_manager,
             )
             options = fidelity.options
             job = fidelity.run(self._circuit[2], self._circuit[3])
@@ -313,8 +297,7 @@ def test_options(self):
             fidelity = ComputeUncompute(
                 sampler_shots,
                 options={"shots": 2048, "dummy": 100},
-                pass_manager=self.pm,
-                num_virtual_qubits=self._circuit[2].num_qubits,
+                pass_manager=self.pass_manager,
             )
             fidelity.update_default_options(shots=100)
             options = fidelity.options
@@ -328,8 +311,7 @@ def test_options(self):
             fidelity = ComputeUncompute(
                 sampler_shots,
                 options={"shots": 2048, "dummy": 100},
-                pass_manager=self.pm,
-                num_virtual_qubits=self._circuit[2].num_qubits,
+                pass_manager=self.pass_manager,
             )
             job = fidelity.run(self._circuit[2], self._circuit[3], shots=50, dummy=None)
             options = fidelity.options