diff --git a/qiskit/algorithms/gradients/__init__.py b/qiskit/algorithms/gradients/__init__.py
index f54f01390529..72aad2402c2d 100644
--- a/qiskit/algorithms/gradients/__init__.py
+++ b/qiskit/algorithms/gradients/__init__.py
@@ -30,6 +30,7 @@
    LinCombEstimatorGradient
    ParamShiftEstimatorGradient
    SPSAEstimatorGradient
+   ReverseEstimatorGradient
 
 Sampler Gradients
 =================
@@ -51,6 +52,7 @@
    BaseQGT
    LinCombQGT
    QFI
+   ReverseQGT
 
 Results
 =======
@@ -81,6 +83,8 @@
 from .sampler_gradient_result import SamplerGradientResult
 from .spsa_estimator_gradient import SPSAEstimatorGradient
 from .spsa_sampler_gradient import SPSASamplerGradient
+from .reverse_gradient.reverse_gradient import ReverseEstimatorGradient
+from .reverse_gradient.reverse_qgt import ReverseQGT
 
 __all__ = [
     "BaseEstimatorGradient",
@@ -101,4 +105,6 @@
     "SamplerGradientResult",
     "SPSAEstimatorGradient",
     "SPSASamplerGradient",
+    "ReverseEstimatorGradient",
+    "ReverseQGT",
 ]
diff --git a/qiskit/algorithms/gradients/base_estimator_gradient.py b/qiskit/algorithms/gradients/base_estimator_gradient.py
index 81a7bddc2041..6ceaa41a51aa 100644
--- a/qiskit/algorithms/gradients/base_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/base_estimator_gradient.py
@@ -48,21 +48,43 @@ def __init__(
         self,
         estimator: BaseEstimator,
         options: Options | None = None,
+        derivative_type: DerivativeType = DerivativeType.REAL,
     ):
-        """
+        r"""
         Args:
             estimator: The estimator 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
+            derivative_type: The type of derivative. Can be either ``DerivativeType.REAL``
+                ``DerivativeType.IMAG``, or ``DerivativeType.COMPLEX``.
+
+                    - ``DerivativeType.REAL`` computes :math:`2 \mathrm{Re}[⟨ψ(ω)|O(θ)|dω ψ(ω)〉]`.
+                    - ``DerivativeType.IMAG`` computes :math:`2 \mathrm{Im}[⟨ψ(ω)|O(θ)|dω ψ(ω)〉]`.
+                    - ``DerivativeType.COMPLEX`` computes :math:`2 ⟨ψ(ω)|O(θ)|dω ψ(ω)〉`.
+
+                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.
         """
         self._estimator: BaseEstimator = estimator
         self._default_options = Options()
         if options is not None:
             self._default_options.update_options(**options)
+        self._derivative_type = derivative_type
+
         self._gradient_circuit_cache: dict[QuantumCircuit, GradientCircuit] = {}
 
+    @property
+    def derivative_type(self) -> DerivativeType:
+        """Return the derivative type (real, imaginary or complex).
+
+        Returns:
+            The derivative type.
+        """
+        return self._derivative_type
+
     def run(
         self,
         circuits: Sequence[QuantumCircuit],
diff --git a/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py b/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
index 0df0e7b9d7b3..7c0e13ce19de 100644
--- a/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
@@ -87,8 +87,12 @@ def __init__(
                 Higher priority setting overrides lower priority setting.
         """
         self._lin_comb_cache = {}
+        super().__init__(estimator, options, derivative_type=derivative_type)
+
+    @BaseEstimatorGradient.derivative_type.setter
+    def derivative_type(self, derivative_type: DerivativeType) -> None:
+        """Set the derivative type."""
         self._derivative_type = derivative_type
-        super().__init__(estimator, options)
 
     def _run(
         self,
@@ -143,7 +147,7 @@ def _run_unique(
             )
             # If its derivative type is `DerivativeType.COMPLEX`, calculate the gradient
             # of the real and imaginary parts separately.
-            meta["derivative_type"] = self._derivative_type
+            meta["derivative_type"] = self.derivative_type
             metadata.append(meta)
             # Combine inputs into a single job to reduce overhead.
             if self._derivative_type == DerivativeType.COMPLEX:
@@ -173,10 +177,14 @@ def _run_unique(
         gradients = []
         partial_sum_n = 0
         for n in all_n:
-            if self._derivative_type == DerivativeType.COMPLEX:
+            # 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
diff --git a/qiskit/algorithms/gradients/lin_comb_qgt.py b/qiskit/algorithms/gradients/lin_comb_qgt.py
index 4d2374c6dd2b..3e3f3c721274 100644
--- a/qiskit/algorithms/gradients/lin_comb_qgt.py
+++ b/qiskit/algorithms/gradients/lin_comb_qgt.py
@@ -33,14 +33,14 @@
 
 
 class LinCombQGT(BaseQGT):
-    """Computes the Quantum Geometric Tensor (QGT) given a pure,
-    parameterized quantum state. This method employs a linear
-    combination of unitaries [1].
+    """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>`_
 
+        [1]: Schuld et al., "Evaluating analytic gradients on quantum hardware" (2018).
+             `arXiv:1811.11184 <https://arxiv.org/pdf/1811.11184.pdf>`_
     """
 
     SUPPORTED_GATES = [
diff --git a/qiskit/algorithms/gradients/reverse_gradient/__init__.py b/qiskit/algorithms/gradients/reverse_gradient/__init__.py
new file mode 100644
index 000000000000..fdb172d367f0
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/__init__.py
@@ -0,0 +1,11 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# 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.
diff --git a/qiskit/algorithms/gradients/reverse_gradient/bind.py b/qiskit/algorithms/gradients/reverse_gradient/bind.py
new file mode 100644
index 000000000000..5380090b7425
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/bind.py
@@ -0,0 +1,53 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022, 2023.
+#
+# 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.
+
+"""Bind values to a parametrized circuit, accepting binds for non-existing parameters in the circuit."""
+
+from __future__ import annotations
+from collections.abc import Iterable
+
+from qiskit.circuit import QuantumCircuit, Parameter
+
+# pylint: disable=inconsistent-return-statements
+def bind(
+    circuits: QuantumCircuit | Iterable[QuantumCircuit],
+    parameter_binds: dict[Parameter, float],
+    inplace: bool = False,
+) -> QuantumCircuit | Iterable[QuantumCircuit] | None:
+    """Bind parameters in a circuit (or list of circuits).
+
+    This method also allows passing parameter binds to parameters that are not in the circuit,
+    and thereby differs to :meth:`.QuantumCircuit.bind_parameters`.
+
+    Args:
+        circuits: Input circuit(s).
+        parameter_binds: A dictionary with ``{Parameter: float}`` pairs determining the values to
+            which the free parameters in the circuit(s) are bound.
+        inplace: If ``True``, bind the values in place, otherwise return circuit copies.
+
+    Returns:
+        The bound circuits, if ``inplace=False``, otherwise None.
+
+    """
+    if not isinstance(circuits, Iterable):
+        circuits = [circuits]
+        return_list = False
+    else:
+        return_list = True
+
+    bound = []
+    for circuit in circuits:
+        existing_parameter_binds = {p: parameter_binds[p] for p in circuit.parameters}
+        bound.append(circuit.assign_parameters(existing_parameter_binds, inplace=inplace))
+
+    if not inplace:
+        return bound if return_list else bound[0]
diff --git a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
new file mode 100644
index 000000000000..0932a5100b78
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
@@ -0,0 +1,156 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022, 2023.
+#
+# 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.
+
+"""Split a circuit into subcircuits, each containing a single parameterized gate."""
+
+from __future__ import annotations
+import itertools
+
+from qiskit.circuit import QuantumCircuit, Parameter, Gate
+from qiskit.circuit.library import RXGate, RYGate, RZGate, CRXGate, CRYGate, CRZGate
+
+
+def gradient_lookup(gate: Gate) -> list[tuple[complex, QuantumCircuit]]:
+    """Returns a circuit implementing the gradient of the input gate.
+
+    Args:
+        gate: The gate whose derivative is returned.
+
+    Returns:
+        The derivative of the input gate as list of ``(coeff, circuit)`` pairs,
+        where the sum of all ``coeff * circuit`` elements describes the full derivative.
+        The circuit is the unitary part of the derivative with a potential separate ``coeff``.
+        The output is a list as derivatives of e.g. controlled gates can only be described
+        as a sum of ``coeff * circuit`` pairs.
+
+    Raises:
+        NotImplementedError: If the derivative of ``gate`` is not implemented.
+    """
+
+    param = gate.params[0]
+    if isinstance(gate, RXGate):
+        derivative = QuantumCircuit(gate.num_qubits)
+        derivative.rx(param, 0)
+        derivative.x(0)
+        return [(-0.5j, derivative)]
+    if isinstance(gate, RYGate):
+        derivative = QuantumCircuit(gate.num_qubits)
+        derivative.ry(param, 0)
+        derivative.y(0)
+        return [(-0.5j, derivative)]
+    if isinstance(gate, RZGate):
+        derivative = QuantumCircuit(gate.num_qubits)
+        derivative.rz(param, 0)
+        derivative.z(0)
+        return [(-0.5j, derivative)]
+    if isinstance(gate, CRXGate):
+        proj1 = QuantumCircuit(gate.num_qubits)
+        proj1.rx(param, 1)
+        proj1.x(1)
+
+        proj2 = QuantumCircuit(gate.num_qubits)
+        proj2.z(0)
+        proj2.rx(param, 1)
+        proj2.x(1)
+
+        return [(-0.25j, proj1), (0.25j, proj2)]
+    if isinstance(gate, CRYGate):
+        proj1 = QuantumCircuit(gate.num_qubits)
+        proj1.ry(param, 1)
+        proj1.y(1)
+
+        proj2 = QuantumCircuit(gate.num_qubits)
+        proj2.z(0)
+        proj2.ry(param, 1)
+        proj2.y(1)
+
+        return [(-0.25j, proj1), (0.25j, proj2)]
+    if isinstance(gate, CRZGate):
+        proj1 = QuantumCircuit(gate.num_qubits)
+        proj1.rz(param, 1)
+        proj1.z(1)
+
+        proj2 = QuantumCircuit(gate.num_qubits)
+        proj2.z(0)
+        proj2.rz(param, 1)
+        proj2.z(1)
+
+        return [(-0.25j, proj1), (0.25j, proj2)]
+    raise NotImplementedError("Cannot implement gradient for", gate)
+
+
+def derive_circuit(
+    circuit: QuantumCircuit, parameter: Parameter
+) -> list[tuple[complex, QuantumCircuit]]:
+    """Return the analytic gradient expression of the input circuit wrt. a single parameter.
+
+    Returns a list of ``(coeff, gradient_circuit)`` tuples, where the derivative of the circuit is
+    given by the sum of the gradient circuits multiplied by their coefficient.
+
+    For example, the circuit::
+
+           ┌───┐┌───────┐┌─────┐
+        q: ┤ H ├┤ Rx(x) ├┤ Sdg ├
+           └───┘└───────┘└─────┘
+
+    returns the coefficient `-0.5j` and the circuit equivalent to::
+
+           ┌───┐┌───────┐┌───┐┌─────┐
+        q: ┤ H ├┤ Rx(x) ├┤ X ├┤ Sdg ├
+           └───┘└───────┘└───┘└─────┘
+
+    as the derivative of `Rx(x)` is `-0.5j Rx(x) X`.
+
+    Args:
+        circuit: The quantum circuit to derive.
+        parameter: The parameter with respect to which we derive.
+
+    Returns:
+        A list of ``(coeff, gradient_circuit)`` tuples.
+
+    Raises:
+        ValueError: If ``parameter`` is of the wrong type.
+        ValueError: If ``parameter`` is not in this circuit.
+        NotImplementedError: If a non-unique parameter is added, as the product rule is not yet
+            supported in this function.
+    """
+    # this is added as useful user-warning, since sometimes ``ParameterExpression``s are
+    # passed around instead of ``Parameter``s
+    if not isinstance(parameter, Parameter):
+        raise ValueError(f"parameter must be of type Parameter, not {type(parameter)}.")
+
+    if parameter not in circuit.parameters:
+        raise ValueError(f"The parameter {parameter} is not in this circuit.")
+
+    if len(circuit._parameter_table[parameter]) > 1:
+        raise NotImplementedError("No product rule support yet, circuit parameters must be unique.")
+
+    summands, op_context = [], []
+    for i, op in enumerate(circuit.data):
+        gate = op[0]
+        op_context += [op[1:]]
+        if parameter in gate.params:
+            coeffs_and_grads = gradient_lookup(gate)
+            summands += [coeffs_and_grads]
+        else:
+            summands += [[(1, gate)]]
+
+    gradient = []
+    for product_rule_term in itertools.product(*summands):
+        summand_circuit = QuantumCircuit(*circuit.qregs)
+        c = 1
+        for i, term in enumerate(product_rule_term):
+            c *= term[0]
+            summand_circuit.data.append([term[1], *op_context[i]])
+        gradient += [(c, summand_circuit.copy())]
+
+    return gradient
diff --git a/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py b/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
new file mode 100644
index 000000000000..004d09f1564c
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
@@ -0,0 +1,200 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# 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.
+
+"""Estimator gradients with the classically efficient reverse mode."""
+
+from __future__ import annotations
+from collections.abc import Sequence
+import logging
+
+import numpy as np
+
+from qiskit.circuit import QuantumCircuit, Parameter
+from qiskit.quantum_info.operators.base_operator import BaseOperator
+from qiskit.quantum_info import Statevector
+from qiskit.opflow import PauliSumOp
+from qiskit.primitives import Estimator
+
+from .bind import bind
+from .derive_circuit import derive_circuit
+from .split_circuits import split
+
+from ..base_estimator_gradient import BaseEstimatorGradient
+from ..estimator_gradient_result import EstimatorGradientResult
+from ..utils import DerivativeType
+
+logger = logging.getLogger(__name__)
+
+
+class ReverseEstimatorGradient(BaseEstimatorGradient):
+    """Estimator gradients with the classically efficient reverse mode.
+
+    .. note::
+
+        This gradient implementation is based on statevector manipulations and scales
+        exponentially with the number of qubits. However, for small system sizes it can be very fast
+        compared to circuit-based gradients.
+
+    This class implements the calculation of the expectation gradient as described in
+    [1]. By keeping track of two statevectors and iteratively sweeping through each parameterized
+    gate, this method scales only linearly with the number of parameters.
+
+    **References:**
+
+        [1]: Jones, T. and Gacon, J. "Efficient calculation of gradients in classical simulations
+             of variational quantum algorithms" (2020).
+             `arXiv:2009.02823 <https://arxiv.org/abs/2009.02823>`_.
+
+    """
+
+    SUPPORTED_GATES = ["rx", "ry", "rz", "cp", "crx", "cry", "crz"]
+
+    def __init__(self, derivative_type: DerivativeType = DerivativeType.REAL):
+        """
+        Args:
+            derivative_type: Defines whether the real, imaginary or real plus imaginary part
+                of the gradient is returned.
+        """
+        dummy_estimator = Estimator()  # this is required by the base class, but not used
+        super().__init__(dummy_estimator, derivative_type=derivative_type)
+
+    @BaseEstimatorGradient.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],
+        observables: Sequence[BaseOperator | PauliSumOp],
+        parameter_values: Sequence[Sequence[float]],
+        parameter_sets: Sequence[set[Parameter]],
+        **options,
+    ) -> EstimatorGradientResult:
+        """Compute the gradients of the expectation values by the parameter shift rule."""
+        g_circuits, g_parameter_values, g_parameter_sets = self._preprocess(
+            circuits, parameter_values, parameter_sets, self.SUPPORTED_GATES
+        )
+        results = self._run_unique(
+            g_circuits, observables, g_parameter_values, g_parameter_sets, **options
+        )
+        return self._postprocess(results, circuits, parameter_values, parameter_sets)
+
+    def _run_unique(
+        self,
+        circuits: Sequence[QuantumCircuit],
+        observables: Sequence[BaseOperator | PauliSumOp],
+        parameter_values: Sequence[Sequence[float]],
+        parameter_sets: Sequence[set[Parameter]],
+        **options,  # pylint: disable=unused-argument
+    ) -> EstimatorGradientResult:
+        num_gradients = len(circuits)
+        gradients = []
+        metadata = []
+
+        for i in range(num_gradients):
+            # temporary variables for easier access
+            circuit = circuits[i]
+            parameters = parameter_sets[i]
+            observable = observables[i]
+            values = parameter_values[i]
+
+            # the metadata only contains the parameters as there are no run configs here
+            metadata.append(
+                {
+                    "parameters": [p for p in circuits[i].parameters if p in parameter_sets[i]],
+                    "derivative_type": self.derivative_type,
+                }
+            )
+
+            # keep track of the parameter order of the circuit, as the circuit splitting might
+            # produce a list of unitaries in a different order
+            original_parameter_order = [p for p in circuit.parameters if p in parameters]
+
+            # split the circuit and generate lists of unitaries [U_1, U_2, ...] and
+            # parameters [p_1, p_2, ...] in these unitaries
+            unitaries, paramlist = split(circuit, parameters=parameters)
+
+            parameter_binds = dict(zip(circuit.parameters, values))
+            bound_circuit = bind(circuit, parameter_binds)
+
+            # initialize state variables -- we use the same naming as in the paper
+            phi = Statevector(bound_circuit)
+            lam = _evolve_by_operator(observable, phi)
+
+            # store gradients in a dictionary to return them in the correct order
+            grads = {param: 0j for param in original_parameter_order}
+
+            num_parameters = len(unitaries)
+            for j in reversed(range(num_parameters)):
+                unitary_j = unitaries[j]
+
+                # We currently only support gates with a single parameter -- which is reflected
+                # in self.SUPPORTED_GATES -- but generally we could also support gates with multiple
+                # parameters per gate
+                parameter_j = paramlist[j][0]
+
+                # get the analytic gradient d U_j / d p_j and bind the gate
+                deriv = derive_circuit(unitary_j, parameter_j)
+                for _, gate in deriv:
+                    bind(gate, parameter_binds, inplace=True)
+
+                # iterate the state variable
+                unitary_j_dagger = bind(unitary_j, parameter_binds).inverse()
+                phi = phi.evolve(unitary_j_dagger)
+
+                # compute current gradient
+                grad = sum(
+                    coeff * lam.conjugate().data.dot(phi.evolve(gate).data) for coeff, gate in deriv
+                )
+
+                # Compute the full gradient (real and complex parts) as all information is available.
+                # Later, based on the derivative type, cast to real/imag/complex.
+                grads[parameter_j] += grad
+
+                if j > 0:
+                    lam = lam.evolve(unitary_j_dagger)
+
+            gradient = np.array(list(grads.values()))
+            gradients.append(self._to_derivtype(gradient))
+
+        result = EstimatorGradientResult(gradients, metadata=metadata, options={})
+        return result
+
+    def _to_derivtype(self, gradient):
+        # 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.REAL:
+            return 2 * np.real(gradient)
+        if self.derivative_type == DerivativeType.IMAG:
+            return 2 * np.imag(gradient)
+
+        return 2 * gradient
+
+
+def _evolve_by_operator(operator, state):
+    """Evolve the Statevector state by operator."""
+
+    # try casting to sparse matrix and use sparse matrix-vector multiplication, which is
+    # a lot faster than using Statevector.evolve
+    if isinstance(operator, PauliSumOp):
+        operator = operator.primitive * operator.coeff
+
+    try:
+        spmatrix = operator.to_matrix(sparse=True)
+        evolved = spmatrix @ state.data
+        return Statevector(evolved)
+    except (TypeError, AttributeError):
+        logger.info("Operator is not castable to a sparse matrix, using Statevector.evolve.")
+
+    return state.evolve(operator)
diff --git a/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
new file mode 100644
index 000000000000..60e8c3cde2d6
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
@@ -0,0 +1,252 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2023.
+#
+# 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 with the classically efficient reverse mode."""
+
+from __future__ import annotations
+from collections.abc import Sequence
+import logging
+
+import numpy as np
+
+from qiskit.circuit import QuantumCircuit, Parameter
+from qiskit.quantum_info import Statevector
+from qiskit.providers import Options
+from qiskit.primitives import Estimator
+
+from ..base_qgt import BaseQGT
+from ..qgt_result import QGTResult
+from ..utils import DerivativeType
+
+from .split_circuits import split
+from .bind import bind
+from .derive_circuit import derive_circuit
+
+logger = logging.getLogger(__name__)
+
+
+class ReverseQGT(BaseQGT):
+    """QGT calculation with the classically efficient reverse mode.
+
+    .. note::
+
+        This QGT implementation is based on statevector manipulations and scales exponentially
+        with the number of qubits. However, for small system sizes it can be very fast
+        compared to circuit-based gradients.
+
+    This class implements the calculation of the QGT as described in [1].
+    By keeping track of three statevectors and iteratively sweeping through each parameterized
+    gate, this method scales only quadratically with the number of parameters.
+
+    **References:**
+
+        [1]: Jones, T. "Efficient classical calculation of the Quantum Natural Gradient" (2020).
+             `arXiv:2011.02991 <https://arxiv.org/abs/2011.02991>`_.
+
+    """
+
+    SUPPORTED_GATES = ["rx", "ry", "rz", "cp", "crx", "cry", "crz"]
+
+    def __init__(
+        self, phase_fix: bool = True, derivative_type: DerivativeType = DerivativeType.COMPLEX
+    ):
+        """
+        Args:
+            phase_fix: Whether or not to include the phase fix.
+            derivative_type: Determines whether the complex QGT or only the real or imaginary
+                parts are calculated.
+        """
+        dummy_estimator = Estimator()  # this method does not need an estimator
+        super().__init__(dummy_estimator, phase_fix, derivative_type)
+
+    @property
+    def options(self) -> Options:
+        """There are no options for the reverse QGT, returns an empty options dict.
+
+        Returns:
+            Empty options.
+        """
+        return Options()
+
+    def _run(
+        self,
+        circuits: Sequence[QuantumCircuit],
+        parameter_values: Sequence[Sequence[float]],
+        parameter_sets: Sequence[set[Parameter]],
+        **options,
+    ) -> QGTResult:
+        """Compute the QGT on the given circuits."""
+        g_circuits, g_parameter_values, g_parameter_sets = self._preprocess(
+            circuits, parameter_values, parameter_sets, self.SUPPORTED_GATES
+        )
+        results = self._run_unique(g_circuits, g_parameter_values, g_parameter_sets, **options)
+        return self._postprocess(results, circuits, parameter_values, parameter_sets)
+
+    def _run_unique(
+        self,
+        circuits: Sequence[QuantumCircuit],
+        parameter_values: Sequence[Sequence[float]],
+        parameter_sets: Sequence[set[Parameter]],
+        **options,  # pylint: disable=unused-argument
+    ) -> QGTResult:
+        num_qgts = len(circuits)
+        qgts = []
+        metadata = []
+
+        for k in range(num_qgts):
+            values = np.asarray(parameter_values[k])
+            circuit = circuits[k]
+            parameters = list(parameter_sets[k])
+
+            num_parameters = len(parameters)
+            original_parameter_order = [p for p in circuit.parameters if p in parameters]
+            metadata.append({"parameters": original_parameter_order})
+
+            unitaries, paramlist = split(circuit, parameters=parameters)
+
+            # initialize the phase fix vector and the hessian part ``metric``
+            num_parameters = len(unitaries)
+            phase_fixes = np.zeros(num_parameters, dtype=complex)
+            metric = np.zeros((num_parameters, num_parameters), dtype=complex)
+
+            # initialize the state variables -- naming convention is the same as the paper
+            parameter_binds = dict(zip(circuit.parameters, values))
+            bound_unitaries = bind(unitaries, parameter_binds)
+
+            chi = Statevector(bound_unitaries[0])
+            psi = chi.copy()
+            phi = Statevector.from_int(0, (2,) * circuit.num_qubits)
+
+            # Get the analytic gradient of the first unitary
+            # Note: We currently only support gates with a single parameter -- which is reflected
+            # in self.SUPPORTED_GATES -- but generally we could also support gates with multiple
+            # parameters per gate. This is the reason for the second 0-index.
+            deriv = derive_circuit(unitaries[0], paramlist[0][0])
+            for _, gate in deriv:
+                bind(gate, parameter_binds, inplace=True)
+
+            grad_coeffs = [coeff for coeff, _ in deriv]
+            grad_states = [phi.evolve(gate) for _, gate in deriv]
+
+            # compute phase fix (optional) and the hessian part
+            if self._phase_fix:
+                phase_fixes[0] = _phasefix_term(chi, grad_coeffs, grad_states)
+
+            metric[0, 0] = _l_term(grad_coeffs, grad_states, grad_coeffs, grad_states)
+
+            for j in range(1, num_parameters):
+                lam = psi.copy()
+                phi = psi.copy()
+
+                # get the analytic gradient d U_j / d p_j and apply it
+                deriv = derive_circuit(unitaries[j], paramlist[j][0])
+
+                for _, gate in deriv:
+                    bind(gate, parameter_binds, inplace=True)
+
+                # compute |phi> (in general it's a sum of states and coeffs)
+                grad_coeffs = [coeff for coeff, _ in deriv]
+                grad_states = [phi.evolve(gate) for _, gate in deriv]
+
+                # compute the digaonal element L_{j, j}
+                metric[j, j] += _l_term(grad_coeffs, grad_states, grad_coeffs, grad_states)
+
+                # compute the off diagonal elements L_{i, j}
+                for i in reversed(range(j)):
+                    # apply U_{i + 1}_dg
+                    unitary_ip_inv = bound_unitaries[i + 1].inverse()
+                    grad_states = [state.evolve(unitary_ip_inv) for state in grad_states]
+
+                    lam = lam.evolve(bound_unitaries[i].inverse())
+
+                    # get the gradient d U_i / d p_i and apply it
+                    deriv = derive_circuit(unitaries[i], paramlist[i][0])
+                    for _, gate in deriv:
+                        bind(gate, parameter_binds, inplace=True)
+
+                    grad_coeffs_mu = [coeff for coeff, _ in deriv]
+                    grad_states_mu = [lam.evolve(gate) for _, gate in deriv]
+
+                    metric[i, j] += _l_term(
+                        grad_coeffs_mu, grad_states_mu, grad_coeffs, grad_states
+                    )
+
+                if self._phase_fix:
+                    phase_fixes[j] += _phasefix_term(chi, grad_coeffs, grad_states)
+
+                psi = psi.evolve(bound_unitaries[j])
+
+            # The following code stacks the QGT together and maps the values into the
+            # correct original order of parameters
+
+            # map circuit parameter to global index in the circuit
+            param_to_circuit = {
+                param: index for index, param in enumerate(original_parameter_order)
+            }
+            # map global index to the local index used in the calculation, the new index can
+            # now be accessed by remap[index]
+            remap = {
+                index: param_to_circuit[_extract_parameter(plist[0])]
+                for index, plist in enumerate(paramlist)
+            }
+
+            qgt = np.zeros((num_parameters, num_parameters), dtype=complex)
+            for i in range(num_parameters):
+                iloc = remap[i]
+                for j in range(num_parameters):
+                    jloc = remap[j]
+                    if i <= j:
+                        qgt[iloc, jloc] += metric[i, j]
+                    else:
+                        qgt[iloc, jloc] += np.conj(metric[j, i])
+
+                    qgt[iloc, jloc] -= np.conj(phase_fixes[i]) * phase_fixes[j]
+
+            # append and cast to real/imag if required
+            qgts.append(self._to_derivtype(qgt))
+
+        result = QGTResult(qgts, self.derivative_type, metadata, options=None)
+        return result
+
+    def _to_derivtype(self, qgt):
+        if self.derivative_type == DerivativeType.REAL:
+            return np.real(qgt)
+        if self.derivative_type == DerivativeType.IMAG:
+            return np.imag(qgt)
+
+        return qgt
+
+
+def _l_term(coeffs_i, states_i, coeffs_j, states_j):
+    return sum(
+        sum(
+            np.conj(coeff_i) * coeff_j * np.conj(state_i.data).dot(state_j.data)
+            for coeff_i, state_i in zip(coeffs_i, states_i)
+        )
+        for coeff_j, state_j in zip(coeffs_j, states_j)
+    )
+
+
+def _phasefix_term(chi, coeffs, states):
+    return sum(
+        coeff_i * np.conj(chi.data).dot(state_i.data) for coeff_i, state_i in zip(coeffs, states)
+    )
+
+
+def _extract_parameter(expression):
+    if isinstance(expression, Parameter):
+        return expression
+
+    if len(expression.parameters) > 1:
+        raise ValueError("Expression has more than one parameter.")
+
+    return list(expression.parameters)[0]
diff --git a/qiskit/algorithms/gradients/reverse_gradient/split_circuits.py b/qiskit/algorithms/gradients/reverse_gradient/split_circuits.py
new file mode 100644
index 000000000000..822ee5853df1
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/split_circuits.py
@@ -0,0 +1,68 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# 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.
+
+"""Split a circuit into subcircuits, each containing a single parameterized gate."""
+
+from __future__ import annotations
+
+from collections.abc import Iterable
+from qiskit.circuit import QuantumCircuit, ParameterExpression, Parameter
+
+
+def split(
+    circuit: QuantumCircuit,
+    parameters: Iterable[Parameter] | None = None,
+) -> tuple[list[QuantumCircuit], list[list[Parameter]]]:
+    """Split the circuit at ParameterExpressions.
+
+    Args:
+        circuit: The circuit to split.
+        parameters: The parameters at which to split. If None, split at each parameter.
+
+    Returns:
+        A list of the split circuits along with a list of which parameters are in the subcircuits.
+    """
+    circuits = []
+    corresponding_parameters = []
+
+    sub = QuantumCircuit(*circuit.qregs, *circuit.cregs)
+    for op in circuit.data:
+        # check if new split must be created
+        if parameters is None:
+            params = [
+                param
+                for param in op[0].params
+                if isinstance(param, ParameterExpression) and len(param.parameters) > 0
+            ]
+        else:
+            if op[0].definition is not None:
+                free_op_params = op[0].definition.parameters
+            else:
+                free_op_params = {}
+
+            params = [p for p in parameters if p in free_op_params]
+
+        new_split = bool(len(params) > 0)
+
+        if new_split:
+            sub.data += [op]
+            circuits.append(sub)
+            corresponding_parameters.append(params)
+            sub = QuantumCircuit(*circuit.qregs, *circuit.cregs)
+        else:
+            sub.data += [op]
+
+    # handle leftover gates
+    if len(sub.data) > 0:
+        circuits[-1].compose(sub, inplace=True)
+
+    return circuits, corresponding_parameters
diff --git a/releasenotes/notes/turbo-gradients-5bebc6e665b900b2.yaml b/releasenotes/notes/turbo-gradients-5bebc6e665b900b2.yaml
new file mode 100644
index 000000000000..85609bd89cdc
--- /dev/null
+++ b/releasenotes/notes/turbo-gradients-5bebc6e665b900b2.yaml
@@ -0,0 +1,21 @@
+---
+features:
+  - |
+    Added the :class:`.ReverseEstimatorGradient` class for a classical, fast evaluation of
+    expectation value gradients based on backpropagation or reverse-mode gradients.
+    This class uses statevectors and thus provides exact gradients but scales 
+    exponentially in system size. It is designed for fast reference calculation of smaller system
+    sizes. It can for example be used as::
+
+      from qiskit.circuit.library import EfficientSU2
+      from qiskit.quantum_info import SparsePauliOp
+      from qiskit.algorithms.gradients import ReverseEstimatorGradient
+
+      observable = SparsePauliOp.from_sparse_list([("ZZ", [0, 1], 1)], num_qubits=10)
+      circuit = EfficientSU2(num_qubits=10)
+      values = [i / 100 for i in range(circuit.num_parameters)]
+      gradient = ReverseEstimatorGradient()
+
+      result = gradient.run([circuit], [observable], [values]).result()
+
+
diff --git a/test/python/algorithms/gradients/test_estimator_gradient.py b/test/python/algorithms/gradients/test_estimator_gradient.py
index 5ee28890f64e..00b9c8079fe3 100644
--- a/test/python/algorithms/gradients/test_estimator_gradient.py
+++ b/test/python/algorithms/gradients/test_estimator_gradient.py
@@ -24,6 +24,7 @@
     LinCombEstimatorGradient,
     ParamShiftEstimatorGradient,
     SPSAEstimatorGradient,
+    ReverseEstimatorGradient,
     DerivativeType,
 )
 from qiskit.circuit import Parameter
@@ -42,6 +43,7 @@
     lambda estimator: FiniteDiffEstimatorGradient(estimator, epsilon=1e-6, method="backward"),
     ParamShiftEstimatorGradient,
     LinCombEstimatorGradient,
+    lambda estimator: ReverseEstimatorGradient(),  # does not take an estimator!
 ]
 
 
@@ -360,14 +362,18 @@ def test_gradient_random_parameters(self, grad):
 
     @data((DerivativeType.IMAG, -1.0), (DerivativeType.COMPLEX, -1.0j))
     @unpack
-    def test_lin_comb_imag_gradient(self, derivative_type, expected_gradient_value):
+    def test_complex_gradient(self, derivative_type, expected_gradient_value):
         """Tests if the ``LinCombEstimatorGradient`` has the correct value."""
         estimator = Estimator()
-        gradient = LinCombEstimatorGradient(estimator, derivative_type=derivative_type)
-        c = QuantumCircuit(1)
-        c.rz(Parameter("p"), 0)
-        result = gradient.run([c], [Pauli("I")], [[0.0]]).result()
-        self.assertAlmostEqual(result.gradients[0][0], expected_gradient_value)
+        lcu = LinCombEstimatorGradient(estimator, derivative_type=derivative_type)
+        reverse = ReverseEstimatorGradient(derivative_type=derivative_type)
+
+        for gradient in [lcu, reverse]:
+            with self.subTest(gradient=gradient):
+                c = QuantumCircuit(1)
+                c.rz(Parameter("p"), 0)
+                result = gradient.run([c], [Pauli("I")], [[0.0]]).result()
+                self.assertAlmostEqual(result.gradients[0][0], expected_gradient_value)
 
     @data(
         FiniteDiffEstimatorGradient,
diff --git a/test/python/algorithms/gradients/test_qfi.py b/test/python/algorithms/gradients/test_qfi.py
index ba112d85ecaf..b0d05ac7030f 100644
--- a/test/python/algorithms/gradients/test_qfi.py
+++ b/test/python/algorithms/gradients/test_qfi.py
@@ -11,27 +11,30 @@
 # that they have been altered from the originals.
 # =============================================================================
 
-""" Test QFI"""
+"""Test QFI."""
 
 import unittest
+from ddt import ddt, data
 
 import numpy as np
 
 from qiskit import QuantumCircuit
-from qiskit.algorithms.gradients import LinCombQGT, QFI
+from qiskit.algorithms.gradients import LinCombQGT, ReverseQGT, QFI, DerivativeType
 from qiskit.circuit import Parameter
 from qiskit.circuit.parametervector import ParameterVector
 from qiskit.primitives import Estimator
 from qiskit.test import QiskitTestCase
 
 
+@ddt
 class TestQFI(QiskitTestCase):
     """Test QFI"""
 
     def setUp(self):
         super().setUp()
         self.estimator = Estimator()
-        self.qgt = LinCombQGT(self.estimator)
+        self.lcu_qgt = LinCombQGT(self.estimator, derivative_type=DerivativeType.REAL)
+        self.reverse_qgt = ReverseQGT(derivative_type=DerivativeType.REAL)
 
     def test_qfi(self):
         """Test if the quantum fisher information calculation is correct for a simple test case.
@@ -47,7 +50,7 @@ def test_qfi(self):
         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.qgt)
+        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)
@@ -69,7 +72,8 @@ def test_qfi_phase_fix(self):
         qfis = qfi.run([qc], [param]).result().qfis
         np.testing.assert_allclose(qfis[0], correct_values, atol=1e-3)
 
-    def test_qfi_maxcut(self):
+    @data("lcu", "reverse")
+    def test_qfi_maxcut(self, qgt_kind):
         """Test the QFI for a simple MaxCut problem.
 
         This is interesting because it contains the same parameters in different gates.
@@ -101,7 +105,8 @@ def expiz(qubit0, qubit1):
         reference = np.array([[16.0, -5.551], [-5.551, 18.497]])
         param = [0.4, 0.69]
 
-        qfi = QFI(self.qgt)
+        qgt = self.lcu_qgt if qgt_kind == "lcu" else self.reverse_qgt
+        qfi = QFI(qgt)
         qfi_result = qfi.run([ansatz], [param]).result().qfis
         np.testing.assert_array_almost_equal(qfi_result[0], reference, decimal=3)
 
diff --git a/test/python/algorithms/gradients/test_qgt.py b/test/python/algorithms/gradients/test_qgt.py
index dd5a3d128aef..044fabdbedfa 100644
--- a/test/python/algorithms/gradients/test_qgt.py
+++ b/test/python/algorithms/gradients/test_qgt.py
@@ -11,14 +11,15 @@
 # that they have been altered from the originals.
 # =============================================================================
 
-""" Test QGT"""
+"""Test QGT."""
 
 import unittest
+from ddt import ddt, data
 
 import numpy as np
 
 from qiskit import QuantumCircuit
-from qiskit.algorithms.gradients import DerivativeType, LinCombQGT
+from qiskit.algorithms.gradients import DerivativeType, LinCombQGT, ReverseQGT
 from qiskit.circuit import Parameter
 from qiskit.circuit.library import RealAmplitudes
 from qiskit.primitives import Estimator
@@ -27,6 +28,7 @@
 from .logging_primitives import LoggingEstimator
 
 
+@ddt
 class TestQGT(QiskitTestCase):
     """Test QGT"""
 
@@ -34,8 +36,12 @@ def setUp(self):
         super().setUp()
         self.estimator = Estimator()
 
-    def test_qgt_derivative_type(self):
+    @data(LinCombQGT, ReverseQGT)
+    def test_qgt_derivative_type(self, qgt_type):
         """Test QGT derivative_type"""
+        args = () if qgt_type == ReverseQGT else (self.estimator,)
+        qgt = qgt_type(*args, derivative_type=DerivativeType.REAL)
+
         a, b = Parameter("a"), Parameter("b")
         qc = QuantumCircuit(1)
         qc.h(0)
@@ -43,37 +49,38 @@ def test_qgt_derivative_type(self):
         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 = LinCombQGT(self.estimator, derivative_type=DerivativeType.REAL)
-            correct_values = np.array([[[1, 0], [0, 0.5]], [[1, 0], [0, 1]]]) / 4
+            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], atol=1e-3)
+                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 = LinCombQGT(self.estimator, derivative_type=DerivativeType.IMAG)
-
-            correct_values = (
-                np.array([[[0, 0.707106781], [-0.707106781, 0]], [[0, 1], [-1, 0]]]) / 4
-            )
+            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], atol=1e-3)
+                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 = LinCombQGT(self.estimator, derivative_type=DerivativeType.COMPLEX)
-            correct_values = (
-                np.array([[[1, 0.707106781j], [-0.707106781j, 0.5]], [[1, 1j], [-1j, 1]]]) / 4
-            )
+            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)
 
-    def test_qgt_phase_fix(self):
+    @data(LinCombQGT, ReverseQGT)
+    def test_qgt_phase_fix(self, qgt_type):
         """Test the phase-fix argument in a QGT calculation"""
+        args = () if qgt_type == ReverseQGT else (self.estimator,)
+        qgt = qgt_type(*args, phase_fix=False)
+
         # create the circuit
         a, b = Parameter("a"), Parameter("b")
         qc = QuantumCircuit(1)
@@ -82,43 +89,42 @@ def test_qgt_phase_fix(self):
         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 = LinCombQGT(self.estimator, phase_fix=False, derivative_type=DerivativeType.REAL)
-            correct_values = np.array([[[1, 0], [0, 1]], [[1, 0], [0, 1]]]) / 4
+            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], atol=1e-3)
+                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 = LinCombQGT(self.estimator, phase_fix=False, derivative_type=DerivativeType.IMAG)
-            correct_values = (
-                np.array([[[0, 0.707106781], [-0.707106781, 0]], [[0, 1], [-1, 0]]]) / 4
-            )
+            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], atol=1e-3)
+                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 = LinCombQGT(
-                self.estimator, phase_fix=False, derivative_type=DerivativeType.COMPLEX
-            )
-            correct_values = (
-                np.array([[[1, 0.707106781j], [-0.707106781j, 1]], [[1, 1j], [-1j, 1]]]) / 4
-            )
+            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)
 
-    def test_qgt_coefficients(self):
+    @data(LinCombQGT, ReverseQGT)
+    def test_qgt_coefficients(self, qgt_type):
         """Test the derivative option of QGT"""
+        args = () if qgt_type == ReverseQGT else (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)
 
-        qgt = LinCombQGT(self.estimator, derivative_type=DerivativeType.REAL)
         # test imaginary derivative
         param_list = [
             [np.pi / 4 for param in qc.parameters],
@@ -147,20 +153,27 @@ def test_qgt_coefficients(self):
             qgt_result = qgt.run([qc], [param]).result().qgts
             np.testing.assert_allclose(qgt_result[0], correct_values[i], atol=1e-3)
 
-    def test_qgt_specify_parameters(self):
+    @data(LinCombQGT, ReverseQGT)
+    def test_qgt_specify_parameters(self, qgt_type):
         """Test the QGT with specified parameters"""
+        args = () if qgt_type == ReverseQGT else (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)
-        qgt = LinCombQGT(self.estimator, derivative_type=DerivativeType.REAL)
         param_list = [np.pi / 4, np.pi / 4]
         qgt_result = qgt.run([qc], [param_list], [[a]]).result().qgts
         np.testing.assert_allclose(qgt_result[0], [[1 / 4]], atol=1e-3)
 
-    def test_qgt_multi_arguments(self):
+    @data(LinCombQGT, ReverseQGT)
+    def test_qgt_multi_arguments(self, qgt_type):
         """Test the QGT for multiple arguments"""
+        args = () if qgt_type == ReverseQGT else (self.estimator,)
+        qgt = qgt_type(*args, derivative_type=DerivativeType.REAL)
+
         a = Parameter("a")
         b = Parameter("b")
         qc = QuantumCircuit(1)
@@ -169,7 +182,6 @@ def test_qgt_multi_arguments(self):
         qc2 = QuantumCircuit(1)
         qc2.rx(a, 0)
         qc2.ry(b, 0)
-        qgt = LinCombQGT(self.estimator, derivative_type=DerivativeType.REAL)
 
         param_list = [[np.pi / 4], [np.pi / 2]]
         correct_values = [[[1 / 4]], [[1 / 4, 0], [0, 0]]]
@@ -178,12 +190,15 @@ def test_qgt_multi_arguments(self):
         for i, _ in enumerate(param_list):
             np.testing.assert_allclose(qgt_results[i], correct_values[i], atol=1e-3)
 
-    def test_qgt_validation(self):
+    @data(LinCombQGT, ReverseQGT)
+    def test_qgt_validation(self, qgt_type):
         """Test estimator QGT's validation"""
+        args = () if qgt_type == ReverseQGT else (self.estimator,)
+        qgt = qgt_type(*args)
+
         a = Parameter("a")
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
-        qgt = LinCombQGT(self.estimator)
         parameter_values = [[np.pi / 4]]
         with self.subTest("assert number of circuits does not match"):
             with self.assertRaises(ValueError):