From a47668b63fe93c7de723909b44a11cc68b95f954 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Tue, 13 Dec 2022 15:33:19 +0100
Subject: [PATCH 01/14] ckassically efficient gradients

---
 qiskit/algorithms/gradients/__init__.py       |   3 +
 .../gradients/reverse_gradient/__init__.py    |  11 ++
 .../gradients/reverse_gradient/bind.py        |  27 +++
 .../reverse_gradient/gradient_lookup.py       | 128 +++++++++++++
 .../reverse_gradient/reverse_gradient.py      | 179 ++++++++++++++++++
 .../reverse_gradient/split_circuits.py        |  73 +++++++
 .../gradients/test_estimator_gradient.py      |  19 ++
 7 files changed, 440 insertions(+)
 create mode 100644 qiskit/algorithms/gradients/reverse_gradient/__init__.py
 create mode 100644 qiskit/algorithms/gradients/reverse_gradient/bind.py
 create mode 100644 qiskit/algorithms/gradients/reverse_gradient/gradient_lookup.py
 create mode 100644 qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
 create mode 100644 qiskit/algorithms/gradients/reverse_gradient/split_circuits.py

diff --git a/qiskit/algorithms/gradients/__init__.py b/qiskit/algorithms/gradients/__init__.py
index 7c2e2af35efb..219db02b3527 100644
--- a/qiskit/algorithms/gradients/__init__.py
+++ b/qiskit/algorithms/gradients/__init__.py
@@ -30,6 +30,7 @@
    LinCombEstimatorGradient
    ParamShiftEstimatorGradient
    SPSAEstimatorGradient
+   ReverseEstimatorGradient
 
 Sampler Gradients
 =================
@@ -77,6 +78,7 @@
 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
 
 __all__ = [
     "BaseEstimatorGradient",
@@ -95,4 +97,5 @@
     "SamplerGradientResult",
     "SPSAEstimatorGradient",
     "SPSASamplerGradient",
+    "ReverseEstimatorGradient",
 ]
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..1b35b133e51b
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/bind.py
@@ -0,0 +1,27 @@
+# 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.
+
+"""Bind a parameters to a circuit, accepting parameters not existing in the circuit."""
+
+# pylint: disable=inconsistent-return-statements
+def bind(circuits, parameter_binds, inplace=False):
+    if not isinstance(circuits, list):
+        existing_parameter_binds = {p: parameter_binds[p] for p in circuits.parameters}
+        return circuits.assign_parameters(existing_parameter_binds, inplace=inplace)
+
+    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
diff --git a/qiskit/algorithms/gradients/reverse_gradient/gradient_lookup.py b/qiskit/algorithms/gradients/reverse_gradient/gradient_lookup.py
new file mode 100644
index 000000000000..f4a08cf25544
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/gradient_lookup.py
@@ -0,0 +1,128 @@
+# 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."""
+
+import itertools
+from qiskit.circuit import QuantumCircuit, Parameter
+from qiskit.circuit.library import RXGate, RYGate, RZGate, CRXGate, CRYGate, CRZGate
+
+
+def extract_single_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]
+
+
+def gradient_lookup(gate):
+    """Returns a circuit implementing the gradient of the input gate."""
+
+    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 for", gate)
+
+
+def analytic_gradient(circuit, parameter=None, return_parameter=False):
+    """Return the analytic gradient of the input circuit."""
+
+    if parameter is not None:
+        single = extract_single_parameter(parameter)
+
+        if single not in circuit.parameters:
+            raise ValueError("Parameter not in this circuit.")
+
+        if len(circuit._parameter_table[single]) > 1:
+            raise NotImplementedError("No product rule support yet, params must be unique.")
+
+    summands, op_context = [], []
+    for i, op in enumerate(circuit.data):
+        gate = op[0]
+        op_context += [op[1:]]
+        if (parameter is None and len(gate.params) > 0) or parameter in gate.params:
+            coeffs_and_grads = gradient_lookup(gate)
+            if not isinstance(parameter, Parameter):
+                # is not a fully decomposed parameter
+                if len(parameter.parameters) > 1:
+                    raise NotImplementedError("Cannot support multiple parameters in one gate yet.")
+                single_parameter = list(parameter.parameters)[0]
+
+                # multiply coefficient with parameter derivative
+                for k, coeff_and_grad in enumerate(coeffs_and_grads):
+                    coeffs_and_grads[k][0] = (
+                        parameter.gradient(single_parameter) * coeff_and_grad[0]
+                    )
+
+            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, a in enumerate(product_rule_term):
+            c *= a[0]
+            summand_circuit.data.append([a[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..5ac20c8c4ebd
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
@@ -0,0 +1,179 @@
+# 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 .gradient_lookup import analytic_gradient
+from .split_circuits import split
+
+from ..base_estimator_gradient import BaseEstimatorGradient
+from ..estimator_gradient_result import EstimatorGradientResult
+from ..lin_comb_estimator_gradient import DerivativeType
+
+logger = logging.getLogger(__name__)
+
+
+class ReverseEstimatorGradient(BaseEstimatorGradient):
+    """Estimator gradients with the classically efficient reverse mode."""
+
+    SUPPORTED_GATES = ["rx", "ry", "rz", "cp", "crx", "cry", "crz"]
+
+    def __init__(self, derivative_type: DerivativeType = DerivativeType.REAL):
+        """
+
+        .. note::
+
+            These gradients are calculated directly on the statevectors. This is
+            inefficient in the number of qubits, but very fast in the number of gates.
+
+        Args:
+            derivative_type: Selects whether the real, imaginary or real + imaginary part
+                of the gradient is returned.
+
+        """
+        alibi_estimator = Estimator()  # this is never used
+        super().__init__(alibi_estimator)
+        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]]}
+            )
+
+            # 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 = analytic_gradient(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):
+        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 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/split_circuits.py b/qiskit/algorithms/gradients/reverse_gradient/split_circuits.py
new file mode 100644
index 000000000000..9f35430bf2d4
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/split_circuits.py
@@ -0,0 +1,73 @@
+# 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 collections.abc import Iterable
+from qiskit.circuit import QuantumCircuit, ParameterExpression, Parameter
+
+
+def split(
+    circuit: QuantumCircuit,
+    parameters: Iterable[Parameter] | None = None,
+):
+    """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
+            ]
+        elif isinstance(parameters, Iterable):
+            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]
+        else:
+            raise NotImplementedError("Unsupported type of parameters:", parameters)
+
+        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]
+
+    if len(sub.data) > 0:  # handle leftover gates
+        # if separate_parameterized_gates or len(circuits) == 0:
+        #     corresponding_parameters.append(params)
+        #     circuits.append(
+        #         sub
+        #     )  # create new if parameterized gates should be separated
+        # else:
+        circuits[-1].compose(sub, inplace=True)  # or attach to last
+
+    return circuits, corresponding_parameters
diff --git a/test/python/algorithms/gradients/test_estimator_gradient.py b/test/python/algorithms/gradients/test_estimator_gradient.py
index feea14423ea0..aac8cd1065b4 100644
--- a/test/python/algorithms/gradients/test_estimator_gradient.py
+++ b/test/python/algorithms/gradients/test_estimator_gradient.py
@@ -25,6 +25,7 @@
     LinCombEstimatorGradient,
     ParamShiftEstimatorGradient,
     SPSAEstimatorGradient,
+    ReverseEstimatorGradient,
     DerivativeType,
 )
 from qiskit.circuit import Parameter
@@ -41,6 +42,7 @@
     lambda estimator: FiniteDiffEstimatorGradient(estimator, epsilon=1e-6, method="backward"),
     ParamShiftEstimatorGradient,
     LinCombEstimatorGradient,
+    ReverseEstimatorGradient,
 ]
 
 
@@ -125,6 +127,23 @@ def test_gradient_u(self, grad):
             for j, value in enumerate(gradients):
                 self.assertAlmostEqual(value, correct_results[i][j], 3)
 
+    @combine(grad=gradient_factories)
+    def test_gradient_shuffled_order(self, grad):
+        """Test computing the gradients out-of-order."""
+        estimator = Estimator()
+        a = Parameter("a")
+        b = Parameter("b")
+        qc = QuantumCircuit(1)
+        qc.rx(0.12, 0)
+        qc.ry(a, 0)
+        qc.rz(b, 0)
+        gradient = grad(estimator)
+        op = SparsePauliOp.from_list([("Z", 1)])
+        values = [1, 2]
+
+        print(gradient.run([qc], [op], [values]).result())
+        print(gradient.run([qc], [op], [values], parameters=[[b, a]]).result())
+
     @combine(grad=gradient_factories)
     def test_gradient_efficient_su2(self, grad):
         """Test the estimator gradient for EfficientSU2"""

From fa85220b10c0b5dc752250c4bde64bd9c21523d5 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Wed, 14 Dec 2022 16:59:04 +0100
Subject: [PATCH 02/14] cleanup & reno

---
 .../gradients/base_estimator_gradient.py      | 21 +++--
 .../gradients/lin_comb_estimator_gradient.py  | 14 ++-
 .../gradients/reverse_gradient/bind.py        | 24 ++++-
 .../{gradient_lookup.py => derive_circuit.py} | 90 +++++++++++--------
 .../reverse_gradient/reverse_gradient.py      | 56 ++++++++----
 .../reverse_gradient/split_circuits.py        | 19 ++--
 .../turbo-gradients-5bebc6e665b900b2.yaml     | 21 +++++
 .../gradients/test_estimator_gradient.py      | 35 +++-----
 8 files changed, 178 insertions(+), 102 deletions(-)
 rename qiskit/algorithms/gradients/reverse_gradient/{gradient_lookup.py => derive_circuit.py} (54%)
 create mode 100644 releasenotes/notes/turbo-gradients-5bebc6e665b900b2.yaml

diff --git a/qiskit/algorithms/gradients/base_estimator_gradient.py b/qiskit/algorithms/gradients/base_estimator_gradient.py
index 2b1d169ec1a0..82abe32eb953 100644
--- a/qiskit/algorithms/gradients/base_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/base_estimator_gradient.py
@@ -63,6 +63,17 @@ def __init__(
             self._default_options.update_options(**options)
         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.
+        """
+        # the default case is real, as this yields e.g. the energy gradient and this type
+        # is also supported by function-level schemes like finite difference or SPSA
+        return DerivativeType.REAL
+
     def run(
         self,
         circuits: Sequence[QuantumCircuit],
@@ -204,13 +215,9 @@ def _postprocess(
         for idx, (circuit, parameter_values_, parameter_set) in enumerate(
             zip(circuits, parameter_values, parameter_sets)
         ):
-            unique_gradient = np.zeros(len(parameter_set))
-            if (
-                "derivative_type" in results.metadata[idx]
-                and results.metadata[idx]["derivative_type"] == DerivativeType.COMPLEX
-            ):
-                # If the derivative type is complex, cast the gradient to complex.
-                unique_gradient = unique_gradient.astype("complex")
+            dtype = complex if self.derivative_type == DerivativeType.COMPLEX else float
+            unique_gradient = np.zeros(len(parameter_set), dtype=dtype)
+
             gradient_circuit = self._gradient_circuit_cache[_circuit_key(circuit)]
             g_parameter_set = _make_gradient_parameter_set(gradient_circuit, parameter_set)
             # Make a map from the gradient parameter to the respective index in the gradient.
diff --git a/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py b/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
index 2797ce48fe54..fb9f47927858 100644
--- a/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
@@ -91,6 +91,16 @@ def __init__(
         self._derivative_type = derivative_type
         super().__init__(estimator, options)
 
+    @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],
@@ -144,7 +154,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:
@@ -174,7 +184,7 @@ def _run_unique(
         gradients = []
         partial_sum_n = 0
         for n in all_n:
-            if self._derivative_type == DerivativeType.COMPLEX:
+            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]
diff --git a/qiskit/algorithms/gradients/reverse_gradient/bind.py b/qiskit/algorithms/gradients/reverse_gradient/bind.py
index 1b35b133e51b..24202f3baf7b 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/bind.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/bind.py
@@ -12,11 +12,27 @@
 
 """Bind a parameters to a circuit, accepting parameters not existing 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, parameter_binds, inplace=False):
+def bind(
+    circuits: QuantumCircuit | Iterable[QuantumCircuit],
+    parameter_binds: dict[Parameter, float],
+    inplace: bool = False,
+) -> QuantumCircuit | Iterable[QuantumCircuit] | None:
+    """Bind parameters to 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`.
+    """
     if not isinstance(circuits, list):
-        existing_parameter_binds = {p: parameter_binds[p] for p in circuits.parameters}
-        return circuits.assign_parameters(existing_parameter_binds, inplace=inplace)
+        circuits = [circuits]
+        return_list = False
+    else:
+        return_list = True
 
     bound = []
     for circuit in circuits:
@@ -24,4 +40,4 @@ def bind(circuits, parameter_binds, inplace=False):
         bound.append(circuit.assign_parameters(existing_parameter_binds, inplace=inplace))
 
     if not inplace:
-        return bound
+        return bound if return_list else bound[0]
diff --git a/qiskit/algorithms/gradients/reverse_gradient/gradient_lookup.py b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
similarity index 54%
rename from qiskit/algorithms/gradients/reverse_gradient/gradient_lookup.py
rename to qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
index f4a08cf25544..ff027067fa57 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/gradient_lookup.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
@@ -12,22 +12,14 @@
 
 """Split a circuit into subcircuits, each containing a single parameterized gate."""
 
+from __future__ import annotations
 import itertools
-from qiskit.circuit import QuantumCircuit, Parameter
-from qiskit.circuit.library import RXGate, RYGate, RZGate, CRXGate, CRYGate, CRZGate
-
-
-def extract_single_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]
+from qiskit.circuit import QuantumCircuit, Parameter, Gate
+from qiskit.circuit.library import RXGate, RYGate, RZGate, CRXGate, CRYGate, CRZGate
 
 
-def gradient_lookup(gate):
+def gradient_lookup(gate: Gate) -> list[tuple[complex, QuantumCircuit]]:
     """Returns a circuit implementing the gradient of the input gate."""
 
     param = gate.params[0]
@@ -35,17 +27,17 @@ def gradient_lookup(gate):
         derivative = QuantumCircuit(gate.num_qubits)
         derivative.rx(param, 0)
         derivative.x(0)
-        return [[-0.5j, derivative]]
+        return [(-0.5j, derivative)]
     if isinstance(gate, RYGate):
         derivative = QuantumCircuit(gate.num_qubits)
         derivative.ry(param, 0)
         derivative.y(0)
-        return [[-0.5j, derivative]]
+        return [(-0.5j, derivative)]
     if isinstance(gate, RZGate):
         derivative = QuantumCircuit(gate.num_qubits)
         derivative.rz(param, 0)
         derivative.z(0)
-        return [[-0.5j, derivative]]
+        return [(-0.5j, derivative)]
     if isinstance(gate, CRXGate):
         proj1 = QuantumCircuit(gate.num_qubits)
         proj1.rx(param, 1)
@@ -56,7 +48,7 @@ def gradient_lookup(gate):
         proj2.rx(param, 1)
         proj2.x(1)
 
-        return [[-0.25j, proj1], [0.25j, proj2]]
+        return [(-0.25j, proj1), (0.25j, proj2)]
     if isinstance(gate, CRYGate):
         proj1 = QuantumCircuit(gate.num_qubits)
         proj1.ry(param, 1)
@@ -67,7 +59,7 @@ def gradient_lookup(gate):
         proj2.ry(param, 1)
         proj2.y(1)
 
-        return [[-0.25j, proj1], [0.25j, proj2]]
+        return [(-0.25j, proj1), (0.25j, proj2)]
     if isinstance(gate, CRZGate):
         proj1 = QuantumCircuit(gate.num_qubits)
         proj1.rz(param, 1)
@@ -78,20 +70,56 @@ def gradient_lookup(gate):
         proj2.rz(param, 1)
         proj2.z(1)
 
-        return [[-0.25j, proj1], [0.25j, proj2]]
+        return [(-0.25j, proj1), (0.25j, proj2)]
     raise NotImplementedError("Cannot implement for", gate)
 
 
-def analytic_gradient(circuit, parameter=None, return_parameter=False):
-    """Return the analytic gradient of the input circuit."""
+def derive_circuit(
+    circuit: QuantumCircuit, parameter: Parameter | None = None
+) -> 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.
+    """
 
     if parameter is not None:
-        single = extract_single_parameter(parameter)
+        # 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("``parameter`` must be  None or of type ``Parameter``.")
 
-        if single not in circuit.parameters:
+        if parameter not in circuit.parameters:
             raise ValueError("Parameter not in this circuit.")
 
-        if len(circuit._parameter_table[single]) > 1:
+        if len(circuit._parameter_table[parameter]) > 1:
             raise NotImplementedError("No product rule support yet, params must be unique.")
 
     summands, op_context = [], []
@@ -100,21 +128,9 @@ def analytic_gradient(circuit, parameter=None, return_parameter=False):
         op_context += [op[1:]]
         if (parameter is None and len(gate.params) > 0) or parameter in gate.params:
             coeffs_and_grads = gradient_lookup(gate)
-            if not isinstance(parameter, Parameter):
-                # is not a fully decomposed parameter
-                if len(parameter.parameters) > 1:
-                    raise NotImplementedError("Cannot support multiple parameters in one gate yet.")
-                single_parameter = list(parameter.parameters)[0]
-
-                # multiply coefficient with parameter derivative
-                for k, coeff_and_grad in enumerate(coeffs_and_grads):
-                    coeffs_and_grads[k][0] = (
-                        parameter.gradient(single_parameter) * coeff_and_grad[0]
-                    )
-
             summands += [coeffs_and_grads]
         else:
-            summands += [[[1, gate]]]
+            summands += [[(1, gate)]]
 
     gradient = []
     for product_rule_term in itertools.product(*summands):
@@ -123,6 +139,6 @@ def analytic_gradient(circuit, parameter=None, return_parameter=False):
         for i, a in enumerate(product_rule_term):
             c *= a[0]
             summand_circuit.data.append([a[1], *op_context[i]])
-        gradient += [[c, summand_circuit.copy()]]
+        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
index 5ac20c8c4ebd..f254f1b93ed9 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
@@ -18,7 +18,6 @@
 
 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
@@ -26,7 +25,7 @@
 from qiskit.primitives import Estimator
 
 from .bind import bind
-from .gradient_lookup import analytic_gradient
+from .derive_circuit import derive_circuit
 from .split_circuits import split
 
 from ..base_estimator_gradient import BaseEstimatorGradient
@@ -37,26 +36,47 @@
 
 
 class ReverseEstimatorGradient(BaseEstimatorGradient):
-    """Estimator gradients with the classically efficient reverse mode."""
+    """Estimator gradients with the classically efficient reverse mode.
 
-    SUPPORTED_GATES = ["rx", "ry", "rz", "cp", "crx", "cry", "crz"]
+    .. note::
 
-    def __init__(self, derivative_type: DerivativeType = DerivativeType.REAL):
-        """
+        This gradient implementation is based on statevector manipulations and scales
+        exponentially in numbers 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 in 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>`_.
 
-        .. note::
+    """
 
-            These gradients are calculated directly on the statevectors. This is
-            inefficient in the number of qubits, but very fast in the number of gates.
+    SUPPORTED_GATES = ["rx", "ry", "rz", "cp", "crx", "cry", "crz"]
 
+    def __init__(self, derivative_type: DerivativeType = DerivativeType.REAL):
+        """
         Args:
-            derivative_type: Selects whether the real, imaginary or real + imaginary part
+            derivative_type: Selects whether the real, imaginary or real plus imaginary part
                 of the gradient is returned.
-
         """
         alibi_estimator = Estimator()  # this is never used
         super().__init__(alibi_estimator)
-        self.derivative_type = derivative_type
+        self._derivative_type = derivative_type
+
+    @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,
@@ -96,7 +116,10 @@ def _run_unique(
 
             # 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]]}
+                {
+                    "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
@@ -112,7 +135,7 @@ def _run_unique(
 
             # initialize state variables -- we use the same naming as in the paper
             phi = Statevector(bound_circuit)
-            lam = evolve_by_operator(observable, phi)
+            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}
@@ -127,7 +150,7 @@ def _run_unique(
                 parameter_j = paramlist[j][0]
 
                 # get the analytic gradient d U_j / d p_j and bind the gate
-                deriv = analytic_gradient(unitary_j, parameter_j)
+                deriv = derive_circuit(unitary_j, parameter_j)
                 for _, gate in deriv:
                     bind(gate, parameter_binds, inplace=True)
 
@@ -158,10 +181,11 @@ def _to_derivtype(self, gradient):
             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):
+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
diff --git a/qiskit/algorithms/gradients/reverse_gradient/split_circuits.py b/qiskit/algorithms/gradients/reverse_gradient/split_circuits.py
index 9f35430bf2d4..822ee5853df1 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/split_circuits.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/split_circuits.py
@@ -12,6 +12,8 @@
 
 """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
 
@@ -19,7 +21,7 @@
 def split(
     circuit: QuantumCircuit,
     parameters: Iterable[Parameter] | None = None,
-):
+) -> tuple[list[QuantumCircuit], list[list[Parameter]]]:
     """Split the circuit at ParameterExpressions.
 
     Args:
@@ -41,15 +43,13 @@ def split(
                 for param in op[0].params
                 if isinstance(param, ParameterExpression) and len(param.parameters) > 0
             ]
-        elif isinstance(parameters, Iterable):
+        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]
-        else:
-            raise NotImplementedError("Unsupported type of parameters:", parameters)
 
         new_split = bool(len(params) > 0)
 
@@ -61,13 +61,8 @@ def split(
         else:
             sub.data += [op]
 
-    if len(sub.data) > 0:  # handle leftover gates
-        # if separate_parameterized_gates or len(circuits) == 0:
-        #     corresponding_parameters.append(params)
-        #     circuits.append(
-        #         sub
-        #     )  # create new if parameterized gates should be separated
-        # else:
-        circuits[-1].compose(sub, inplace=True)  # or attach to last
+    # 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..f93c29f52812
--- /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 is based on 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 aac8cd1065b4..8ae7d2695b06 100644
--- a/test/python/algorithms/gradients/test_estimator_gradient.py
+++ b/test/python/algorithms/gradients/test_estimator_gradient.py
@@ -42,7 +42,7 @@
     lambda estimator: FiniteDiffEstimatorGradient(estimator, epsilon=1e-6, method="backward"),
     ParamShiftEstimatorGradient,
     LinCombEstimatorGradient,
-    ReverseEstimatorGradient,
+    lambda estimator: ReverseEstimatorGradient(),  # does not take an estimator!
 ]
 
 
@@ -127,23 +127,6 @@ def test_gradient_u(self, grad):
             for j, value in enumerate(gradients):
                 self.assertAlmostEqual(value, correct_results[i][j], 3)
 
-    @combine(grad=gradient_factories)
-    def test_gradient_shuffled_order(self, grad):
-        """Test computing the gradients out-of-order."""
-        estimator = Estimator()
-        a = Parameter("a")
-        b = Parameter("b")
-        qc = QuantumCircuit(1)
-        qc.rx(0.12, 0)
-        qc.ry(a, 0)
-        qc.rz(b, 0)
-        gradient = grad(estimator)
-        op = SparsePauliOp.from_list([("Z", 1)])
-        values = [1, 2]
-
-        print(gradient.run([qc], [op], [values]).result())
-        print(gradient.run([qc], [op], [values], parameters=[[b, a]]).result())
-
     @combine(grad=gradient_factories)
     def test_gradient_efficient_su2(self, grad):
         """Test the estimator gradient for EfficientSU2"""
@@ -380,14 +363,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)
 
     @combine(
         grad=[

From fc375abc1faed989c470ee105e9f5f1080e42aed Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Tue, 10 Jan 2023 16:41:54 +0100
Subject: [PATCH 03/14] Apply suggestions from code review

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>
---
 qiskit/algorithms/gradients/reverse_gradient/bind.py |  4 ++--
 .../gradients/reverse_gradient/derive_circuit.py     | 10 +++++-----
 .../gradients/reverse_gradient/reverse_gradient.py   | 12 ++++++------
 .../notes/turbo-gradients-5bebc6e665b900b2.yaml      |  2 +-
 4 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/qiskit/algorithms/gradients/reverse_gradient/bind.py b/qiskit/algorithms/gradients/reverse_gradient/bind.py
index 24202f3baf7b..06dcf56befee 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/bind.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/bind.py
@@ -10,7 +10,7 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-"""Bind a parameters to a circuit, accepting parameters not existing in the circuit."""
+"""Bind parameter values to a parametrized circuit, accepting binds for parameters not existing in the circuit."""
 
 from __future__ import annotations
 from collections.abc import Iterable
@@ -23,7 +23,7 @@ def bind(
     parameter_binds: dict[Parameter, float],
     inplace: bool = False,
 ) -> QuantumCircuit | Iterable[QuantumCircuit] | None:
-    """Bind parameters to a circuit (or list of circuits).
+    """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`.
diff --git a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
index ff027067fa57..a8959d0a48eb 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
@@ -71,7 +71,7 @@ def gradient_lookup(gate: Gate) -> list[tuple[complex, QuantumCircuit]]:
         proj2.z(1)
 
         return [(-0.25j, proj1), (0.25j, proj2)]
-    raise NotImplementedError("Cannot implement for", gate)
+    raise NotImplementedError("Cannot implement gradient for", gate)
 
 
 def derive_circuit(
@@ -114,7 +114,7 @@ def derive_circuit(
         # 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("``parameter`` must be  None or of type ``Parameter``.")
+            raise ValueError("parameter  must be  None or of type Parameter.")
 
         if parameter not in circuit.parameters:
             raise ValueError("Parameter not in this circuit.")
@@ -136,9 +136,9 @@ def derive_circuit(
     for product_rule_term in itertools.product(*summands):
         summand_circuit = QuantumCircuit(*circuit.qregs)
         c = 1
-        for i, a in enumerate(product_rule_term):
-            c *= a[0]
-            summand_circuit.data.append([a[1], *op_context[i]])
+        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
index f254f1b93ed9..29af8964d489 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
@@ -30,7 +30,7 @@
 
 from ..base_estimator_gradient import BaseEstimatorGradient
 from ..estimator_gradient_result import EstimatorGradientResult
-from ..lin_comb_estimator_gradient import DerivativeType
+from ..utils import DerivativeType
 
 logger = logging.getLogger(__name__)
 
@@ -41,12 +41,12 @@ class ReverseEstimatorGradient(BaseEstimatorGradient):
     .. note::
 
         This gradient implementation is based on statevector manipulations and scales
-        exponentially in numbers of qubits. However, for small system sizes it can be very fast
+        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 in the number of parameters.
+    gate, this method scales only linearly with the number of parameters.
 
     **References:**
 
@@ -61,11 +61,11 @@ class ReverseEstimatorGradient(BaseEstimatorGradient):
     def __init__(self, derivative_type: DerivativeType = DerivativeType.REAL):
         """
         Args:
-            derivative_type: Selects whether the real, imaginary or real plus imaginary part
+            derivative_type: Defines whether the real, imaginary or real plus imaginary part
                 of the gradient is returned.
         """
-        alibi_estimator = Estimator()  # this is never used
-        super().__init__(alibi_estimator)
+        dummy_estimator = Estimator()  # this is required by the base class, but not used
+        super().__init__(dummy_estimator)
         self._derivative_type = derivative_type
 
     @property
diff --git a/releasenotes/notes/turbo-gradients-5bebc6e665b900b2.yaml b/releasenotes/notes/turbo-gradients-5bebc6e665b900b2.yaml
index f93c29f52812..85609bd89cdc 100644
--- a/releasenotes/notes/turbo-gradients-5bebc6e665b900b2.yaml
+++ b/releasenotes/notes/turbo-gradients-5bebc6e665b900b2.yaml
@@ -3,7 +3,7 @@ features:
   - |
     Added the :class:`.ReverseEstimatorGradient` class for a classical, fast evaluation of
     expectation value gradients based on backpropagation or reverse-mode gradients.
-    This class is based on statevectors and thus provides exact gradients but scales 
+    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 bcc3d8d2452b0718fe7a292efbb7ec5075660d51 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Tue, 10 Jan 2023 16:48:51 +0100
Subject: [PATCH 04/14] Remove support for Parameter = None

---
 .../gradients/reverse_gradient/bind.py        |  4 ++--
 .../reverse_gradient/derive_circuit.py        | 22 +++++++++----------
 2 files changed, 12 insertions(+), 14 deletions(-)

diff --git a/qiskit/algorithms/gradients/reverse_gradient/bind.py b/qiskit/algorithms/gradients/reverse_gradient/bind.py
index 06dcf56befee..b65442cc334d 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/bind.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/bind.py
@@ -10,7 +10,7 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-"""Bind parameter values to a parametrized circuit, accepting binds for parameters not existing in the circuit."""
+"""Bind values to a parametrized circuit, accepting binds for non-existing parameters in the circuit."""
 
 from __future__ import annotations
 from collections.abc import Iterable
@@ -28,7 +28,7 @@ def bind(
     This method also allows passing parameter binds to parameters that are not in the circuit,
     and thereby differs to :meth:`.QuantumCircuit.bind_parameters`.
     """
-    if not isinstance(circuits, list):
+    if not isinstance(circuits, Iterable):
         circuits = [circuits]
         return_list = False
     else:
diff --git a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
index a8959d0a48eb..4a4e5b9f4619 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
@@ -75,7 +75,7 @@ def gradient_lookup(gate: Gate) -> list[tuple[complex, QuantumCircuit]]:
 
 
 def derive_circuit(
-    circuit: QuantumCircuit, parameter: Parameter | None = None
+    circuit: QuantumCircuit, parameter: Parameter
 ) -> list[tuple[complex, QuantumCircuit]]:
     """Return the analytic gradient expression of the input circuit wrt. a single parameter.
 
@@ -109,24 +109,22 @@ def derive_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("parameter  must be  None or of type Parameter.")
 
-    if parameter is not None:
-        # 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("parameter  must be  None or of type Parameter.")
+    if parameter not in circuit.parameters:
+        raise ValueError("Parameter not in this circuit.")
 
-        if parameter not in circuit.parameters:
-            raise ValueError("Parameter not in this circuit.")
-
-        if len(circuit._parameter_table[parameter]) > 1:
-            raise NotImplementedError("No product rule support yet, params must be unique.")
+    if len(circuit._parameter_table[parameter]) > 1:
+        raise NotImplementedError("No product rule support yet, params must be unique.")
 
     summands, op_context = [], []
     for i, op in enumerate(circuit.data):
         gate = op[0]
         op_context += [op[1:]]
-        if (parameter is None and len(gate.params) > 0) or parameter in gate.params:
+        if parameter in gate.params:
             coeffs_and_grads = gradient_lookup(gate)
             summands += [coeffs_and_grads]
         else:

From d112def56cbd9a36ebe743d7e37f5d75888061ab Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Wed, 11 Jan 2023 09:29:33 +0100
Subject: [PATCH 05/14] Complete the docs

---
 .../gradients/reverse_gradient/bind.py        | 10 ++++++++++
 .../reverse_gradient/derive_circuit.py        | 19 +++++++++++++++----
 2 files changed, 25 insertions(+), 4 deletions(-)

diff --git a/qiskit/algorithms/gradients/reverse_gradient/bind.py b/qiskit/algorithms/gradients/reverse_gradient/bind.py
index b65442cc334d..0236a30ff7e0 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/bind.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/bind.py
@@ -27,6 +27,16 @@ def bind(
 
     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]
diff --git a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
index 4a4e5b9f4619..266b02abaa0d 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
@@ -20,7 +20,18 @@
 
 
 def gradient_lookup(gate: Gate) -> list[tuple[complex, QuantumCircuit]]:
-    """Returns a circuit implementing the gradient of the input gate."""
+    """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.
+    """
 
     param = gate.params[0]
     if isinstance(gate, RXGate):
@@ -112,13 +123,13 @@ def derive_circuit(
     # 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("parameter  must be  None or of type Parameter.")
+        raise ValueError(f"parameter must be None or of type Parameter, not {type(parameter)}.")
 
     if parameter not in circuit.parameters:
-        raise ValueError("Parameter not in this circuit.")
+        raise ValueError("The parameter is not in this circuit.")
 
     if len(circuit._parameter_table[parameter]) > 1:
-        raise NotImplementedError("No product rule support yet, params must be unique.")
+        raise NotImplementedError("No product rule support yet, circuit parameters must be unique.")
 
     summands, op_context = [], []
     for i, op in enumerate(circuit.data):

From 1a4f4c44956426840b8a47efcfab5f3b414e5bdc Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Fri, 13 Jan 2023 15:58:15 +0100
Subject: [PATCH 06/14] QGT v0

---
 qiskit/algorithms/gradients/__init__.py       |   3 +
 .../gradients/reverse_gradient/reverse_qgt.py | 261 ++++++++++++++++++
 test/python/algorithms/gradients/test_qfi.py  |  11 +-
 3 files changed, 272 insertions(+), 3 deletions(-)
 create mode 100644 qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py

diff --git a/qiskit/algorithms/gradients/__init__.py b/qiskit/algorithms/gradients/__init__.py
index 219db02b3527..15b307c4c3a9 100644
--- a/qiskit/algorithms/gradients/__init__.py
+++ b/qiskit/algorithms/gradients/__init__.py
@@ -51,6 +51,7 @@
 
    BaseQFI
    LinCombQFI
+   ReverseQGT
 
 Results
 =======
@@ -79,6 +80,7 @@
 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",
@@ -98,4 +100,5 @@
     "SPSAEstimatorGradient",
     "SPSASamplerGradient",
     "ReverseEstimatorGradient",
+    "ReverseQGT",
 ]
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..636e7dfadd44
--- /dev/null
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
@@ -0,0 +1,261 @@
+# 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.
+
+"""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.operators.base_operator import BaseOperator
+from qiskit.quantum_info import Statevector
+from qiskit.opflow import PauliSumOp
+from qiskit.quantum_info import Statevector
+
+from ..base_qfi import BaseQFI
+from ..qfi_result import QFIResult
+from ..utils import DerivativeType
+
+logger = logging.getLogger(__name__)
+
+
+from surfer.tools.gradient_lookup import extract_single_parameter
+
+from .split_circuits import split
+from .bind import bind
+from .derive_circuit import derive_circuit
+
+
+class ReverseQGT(BaseQFI):
+    """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.
+        """
+        super().__init__()
+        self.phase_fix = phase_fix
+        self._derivative_type = derivative_type
+
+    # TODO this should be in the base class of QGT
+    @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],
+        **options,
+    ) -> QFIResult:
+        # cast to array, if values are a list
+        num_qgts = len(circuits)
+        qgts = []
+        metadata = []
+
+        for k in range(num_qgts):
+
+            # TODO unrolling done by QGT
+            # circuit = self.unroller(circuit)
+
+            values = np.asarray(parameter_values[k])
+            circuit = circuits[k]
+
+            # TODO parameters is None captured by QGT
+            if parameters[k] is None:
+                parameters[k] = circuit.parameters
+
+            num_parameters = len(parameters[k])
+            original_parameter_order = [p for p in circuit.parameters if p in parameters]
+
+            metadata.append(
+                {
+                    "parameters": original_parameter_order,
+                    "derivative_type": self.derivative_type,
+                }
+            )
+
+            unitaries, paramlist = split(circuit, parameters=parameters)
+
+            # initialize the phase fix vector and the hessian part ``lis``
+            num_parameters = len(unitaries)
+            phase_fixes = np.zeros(num_parameters, dtype=complex)
+            lis = 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
+            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)
+
+            lis[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.decompose()) for _, gate in deriv]
+
+                # compute the digaonal element L_{j, j}
+                lis[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]
+
+                    lis[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])
+
+            # stack quantum geometric tensor together and take into account the original
+            # order of parameters
+            param_to_circuit = {
+                param: index for index, param in enumerate(original_parameter_order)
+            }
+            remap = {
+                index: param_to_circuit[extract_single_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] += lis[i, j]
+                    else:
+                        qgt[iloc, jloc] += np.conj(lis[j, i])
+
+                    qgt[iloc, jloc] -= np.conj(phase_fixes[i]) * phase_fixes[j]
+
+            if self.derivative_type == DerivativeType.REAL:
+                qgt = np.real(qgt)
+            elif self.derivative_type == DerivativeType.IMAG:
+                qgt = np.imag(qgt)
+
+            qgts.append(self._to_derivtype(qgt))
+
+        result = QFIResult(qgts, metadata, options=None)
+        return result
+
+    def _to_derivtype(self, qgt):
+        # TODO remove factor 4 once the QGT interface is there
+        if self.derivative_type == DerivativeType.REAL:
+            return 4 * np.real(qgt)
+        if self.derivative_type == DerivativeType.IMAG:
+            return 4 * np.imag(qgt)
+
+        return 4 * gradient
+
+
+def _l_term(coeffs_i, states_i, coeffs_j, states_j):
+    return sum(
+        sum(
+            np.conj(c_i) * c_j * np.conj(state_i.data).dot(state_j.data)
+            for c_i, state_i in zip(coeffs_i, states_i)
+        )
+        for c_j, state_j in zip(coeffs_j, states_j)
+    )
+
+
+def _phasefix_term(chi, coeffs, states):
+    return sum(c_i * np.conj(chi.data).dot(state_i.data) for c_i, state_i in zip(coeffs, states))
+
+
+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 AttributeError:
+        logger.info("Operator is not castable to a sparse matrix, using Statevector.evolve.")
+
+    return state.evolve(operator)
diff --git a/test/python/algorithms/gradients/test_qfi.py b/test/python/algorithms/gradients/test_qfi.py
index 276b1c7c370b..19cc379f9d1e 100644
--- a/test/python/algorithms/gradients/test_qfi.py
+++ b/test/python/algorithms/gradients/test_qfi.py
@@ -14,12 +14,13 @@
 """ Test QFI"""
 
 import unittest
+from ddt import ddt, data
 
 import numpy as np
 
 from qiskit import QuantumCircuit
 from qiskit.algorithms.gradients.lin_comb_estimator_gradient import DerivativeType
-from qiskit.algorithms.gradients.lin_comb_qfi import LinCombQFI
+from qiskit.algorithms.gradients import LinCombQFI, ReverseQGT
 from qiskit.circuit import Parameter
 from qiskit.circuit.library import RealAmplitudes
 from qiskit.circuit.parametervector import ParameterVector
@@ -27,6 +28,7 @@
 from qiskit.test import QiskitTestCase
 
 
+@ddt
 class TestQFI(QiskitTestCase):
     """Test QFI"""
 
@@ -34,11 +36,15 @@ def setUp(self):
         super().setUp()
         self.estimator = Estimator()
 
-    def test_qfi_simple(self):
+    @data(LinCombQFI, ReverseQGT)
+    def test_qfi_simple(self, qgt_type):
         """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)]]
         """
+        args = (self.estimator,) if qgt_type != ReverseQGT else (,)
+        qfi = qgt_type(*args)
+
         # create the circuit
         a, b = Parameter("a"), Parameter("b")
         qc = QuantumCircuit(1)
@@ -49,7 +55,6 @@ def test_qfi_simple(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 = LinCombQFI(self.estimator)
         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)

From 62e10dfd2fefcde32ea33fc79ce03f4f5e2ba23f Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Fri, 13 Jan 2023 16:49:50 +0100
Subject: [PATCH 07/14] fix LCU tests

---
 qiskit/algorithms/gradients/lin_comb_qfi.py   | 35 +++++++++++++------
 .../gradients/reverse_gradient/reverse_qgt.py | 30 ++++++++++------
 test/python/algorithms/gradients/test_qfi.py  | 34 +++++++++++-------
 3 files changed, 66 insertions(+), 33 deletions(-)

diff --git a/qiskit/algorithms/gradients/lin_comb_qfi.py b/qiskit/algorithms/gradients/lin_comb_qfi.py
index 60d18da562bf..e85b511ce6d2 100644
--- a/qiskit/algorithms/gradients/lin_comb_qfi.py
+++ b/qiskit/algorithms/gradients/lin_comb_qfi.py
@@ -95,6 +95,16 @@ def __init__(
         )
         self._qfi_circuit_cache = {}
 
+    @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],
@@ -123,7 +133,7 @@ def _run(
                         result_map[i] = -1
 
             meta = {"parameters": [p for p in circuit.parameters if p in param_set]}
-            meta["derivative_type"] = self._derivative_type
+            meta["derivative_type"] = self.derivative_type
             metadata_.append(meta)
 
             observable = SparsePauliOp.from_list([("I" * circuit.num_qubits, 1)])
@@ -241,20 +251,25 @@ def _run(
                     qfi_[idx] += complex(0, coeff * grad_)
             else:
                 for grad_, idx, coeff in zip(result.values, result_indices_all[i], coeffs_all[i]):
-                    qfi_[idx] += coeff * grad_
-                qfi_ = qfi_.real
+                    if metadata_[i]["derivative_type"] == DerivativeType.REAL:
+                        qfi_[idx] += coeff * grad_
+                    else:
+                        qfi_[idx] += complex(0, coeff * grad_)
 
-            if metadata_[i]["derivative_type"] == DerivativeType.REAL:
-                phase_fixes[i] = phase_fixes[i].real
-            elif metadata_[i]["derivative_type"] == DerivativeType.IMAG:
-                phase_fixes[i] = phase_fixes[i].imag
-            qfi = qfi_ - phase_fixes[i]
-            qfi += np.triu(qfi_, k=1).T
-            qfis.append(qfi)
+            qfi = qfi_ + np.conj(np.triu(qfi_, k=1).T) - phase_fixes[i]
+            qfis.append(self._to_derivtype(qfi))
 
         run_opt = self._get_local_options(options)
         return QFIResult(qfis=qfis, metadata=metadata_, options=run_opt)
 
+    def _to_derivtype(self, qfi):
+        if self.derivative_type == DerivativeType.REAL:
+            return np.real(qfi)
+        if self.derivative_type == DerivativeType.IMAG:
+            return np.imag(qfi)
+
+        return qfi
+
     @property
     def options(self) -> Options:
         """Return the union of estimator options setting and QFI default options,
diff --git a/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
index 636e7dfadd44..434fa3547251 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
@@ -23,6 +23,7 @@
 from qiskit.quantum_info import Statevector
 from qiskit.opflow import PauliSumOp
 from qiskit.quantum_info import Statevector
+from qiskit.providers import Options
 
 from ..base_qfi import BaseQFI
 from ..qfi_result import QFIResult
@@ -60,8 +61,9 @@ class ReverseQGT(BaseQFI):
 
     SUPPORTED_GATES = ["rx", "ry", "rz", "cp", "crx", "cry", "crz"]
 
+    # TODO: should the default for QGT be complex?
     def __init__(
-        self, phase_fix: bool = True, derivative_type: DerivativeType = DerivativeType.COMPLEX
+        self, phase_fix: bool = True, derivative_type: DerivativeType = DerivativeType.REAL
     ):
         """
         Args:
@@ -73,6 +75,15 @@ def __init__(
         self.phase_fix = phase_fix
         self._derivative_type = 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()
+
     # TODO this should be in the base class of QGT
     @property
     def derivative_type(self) -> DerivativeType:
@@ -106,10 +117,12 @@ def _run(
 
             # TODO parameters is None captured by QGT
             if parameters[k] is None:
-                parameters[k] = circuit.parameters
+                parameters_ = list(circuit.parameters)
+            else:
+                parameters_ = list(parameters[k])
 
-            num_parameters = len(parameters[k])
-            original_parameter_order = [p for p in circuit.parameters if p in parameters]
+            num_parameters = len(parameters_)
+            original_parameter_order = [p for p in circuit.parameters if p in parameters_]
 
             metadata.append(
                 {
@@ -118,7 +131,7 @@ def _run(
                 }
             )
 
-            unitaries, paramlist = split(circuit, parameters=parameters)
+            unitaries, paramlist = split(circuit, parameters=parameters_)
 
             # initialize the phase fix vector and the hessian part ``lis``
             num_parameters = len(unitaries)
@@ -209,11 +222,6 @@ def _run(
 
                     qgt[iloc, jloc] -= np.conj(phase_fixes[i]) * phase_fixes[j]
 
-            if self.derivative_type == DerivativeType.REAL:
-                qgt = np.real(qgt)
-            elif self.derivative_type == DerivativeType.IMAG:
-                qgt = np.imag(qgt)
-
             qgts.append(self._to_derivtype(qgt))
 
         result = QFIResult(qgts, metadata, options=None)
@@ -226,7 +234,7 @@ def _to_derivtype(self, qgt):
         if self.derivative_type == DerivativeType.IMAG:
             return 4 * np.imag(qgt)
 
-        return 4 * gradient
+        return 4 * qgt
 
 
 def _l_term(coeffs_i, states_i, coeffs_j, states_j):
diff --git a/test/python/algorithms/gradients/test_qfi.py b/test/python/algorithms/gradients/test_qfi.py
index 19cc379f9d1e..0cf34f35a73e 100644
--- a/test/python/algorithms/gradients/test_qfi.py
+++ b/test/python/algorithms/gradients/test_qfi.py
@@ -42,7 +42,7 @@ def test_qfi_simple(self, qgt_type):
 
         QFI = [[1, 0], [0, 1]] - [[0, 0], [0, cos^2(a)]]
         """
-        args = (self.estimator,) if qgt_type != ReverseQGT else (,)
+        args = (self.estimator,) if qgt_type != ReverseQGT else ()
         qfi = qgt_type(*args)
 
         # create the circuit
@@ -59,11 +59,15 @@ def test_qfi_simple(self, qgt_type):
             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):
+    @data(LinCombQFI, ReverseQGT)
+    def test_qfi_phase_fix(self, qgt_type):
         """Test the phase-fix argument in a QFI calculation
 
         QFI = [[1, 0], [0, 1]].
         """
+        args = (self.estimator,) if qgt_type != ReverseQGT else ()
+        qfi = qgt_type(*args, phase_fix=False)
+
         # create the circuit
         a, b = Parameter("a"), Parameter("b")
         qc = QuantumCircuit(1)
@@ -74,11 +78,11 @@ def test_qfi_phase_fix(self):
         param = [np.pi / 4, 0.1]
         # test for different values
         correct_values = [[1, 0], [0, 1]]
-        qfi = LinCombQFI(self.estimator, phase_fix=False)
         qfi_result = qfi.run([qc], [param]).result().qfis
         np.testing.assert_allclose(qfi_result[0], correct_values, atol=1e-3)
 
-    def test_qfi_maxcut(self):
+    @data(LinCombQFI, ReverseQGT)
+    def test_qfi_maxcut(self, qgt_type):
         """Test the QFI for a simple MaxCut problem.
 
         This is interesting because it contains the same parameters in different gates.
@@ -110,11 +114,14 @@ def expiz(qubit0, qubit1):
         reference = np.array([[16.0, -5.551], [-5.551, 18.497]])
         param = [0.4, 0.69]
 
-        qfi = LinCombQFI(self.estimator)
+        args = (self.estimator,) if qgt_type != ReverseQGT else ()
+        qfi = qgt_type(*args)
+
         qfi_result = qfi.run([ansatz], [param]).result().qfis
         np.testing.assert_array_almost_equal(qfi_result[0], reference, decimal=3)
 
-    def test_qfi_derivative_type(self):
+    @data(LinCombQFI, ReverseQGT)
+    def test_qfi_derivative_type(self, qgt_type):
         """Test QFI derivative_type"""
         a, b = Parameter("a"), Parameter("b")
         qc = QuantumCircuit(1)
@@ -122,19 +129,22 @@ def test_qfi_derivative_type(self):
         qc.rz(a, 0)
         qc.rx(b, 0)
 
+        args = (self.estimator,) if qgt_type != ReverseQGT else ()
+        qfi = qgt_type(*args)
+
+        param_list = [[np.pi / 4, 0], [np.pi / 2, np.pi / 4]]
         # test imaginary derivative
         with self.subTest("Test with DerivativeType.IMAG"):
-            qfi = LinCombQFI(self.estimator, derivative_type=DerivativeType.IMAG)
-            param_list = [[np.pi / 4, 0], [np.pi / 2, np.pi / 4]]
-            correct_values = [[[0, 0.707106781], [0.707106781, 0]], [[0, 1], [1, 0]]]
+            qfi.derivative_type = DerivativeType.IMAG
+            correct_values = [[[0, 0.707106781], [-0.707106781, 0]], [[0, 1], [-1, 0]]]
             for i, param in enumerate(param_list):
                 qfi_result = qfi.run([qc], [param]).result().qfis
                 np.testing.assert_allclose(qfi_result[0], correct_values[i], atol=1e-3)
 
         # test real + imaginary derivative
-        with self.subTest("Test with DerivativeType.IMAG"):
-            qfi = LinCombQFI(self.estimator, derivative_type=DerivativeType.COMPLEX)
-            correct_values = [[[1, 0.707106781j], [0.707106781j, 0.5]], [[1, 1j], [1j, 1]]]
+        with self.subTest("Test with DerivativeType.COMPLEX"):
+            qfi.derivative_type = DerivativeType.COMPLEX
+            correct_values = [[[1, 0.707106781j], [-0.707106781j, 0.5]], [[1, 1j], [-1j, 1]]]
             for i, param in enumerate(param_list):
                 qfi_result = qfi.run([qc], [param]).result().qfis
                 np.testing.assert_allclose(qfi_result[0], correct_values[i], atol=1e-3)

From e47ca94e0995cea5c1279616fa0e34f24038f45f Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Fri, 13 Jan 2023 19:53:48 +0100
Subject: [PATCH 08/14] final fixes

---
 .../gradients/reverse_gradient/reverse_qgt.py | 81 ++++++++++---------
 test/python/algorithms/gradients/test_qfi.py  | 22 +++--
 2 files changed, 57 insertions(+), 46 deletions(-)

diff --git a/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
index 434fa3547251..4a2ca10dac17 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
@@ -19,24 +19,21 @@
 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.quantum_info import Statevector
 from qiskit.providers import Options
+from qiskit.transpiler.passes import TranslateParameterizedGates
 
 from ..base_qfi import BaseQFI
 from ..qfi_result import QFIResult
 from ..utils import DerivativeType
 
-logger = logging.getLogger(__name__)
-
-
-from surfer.tools.gradient_lookup import extract_single_parameter
-
 from .split_circuits import split
 from .bind import bind
 from .derive_circuit import derive_circuit
+from .reverse_gradient import _evolve_by_operator
+
+logger = logging.getLogger(__name__)
 
 
 class ReverseQGT(BaseQFI):
@@ -108,13 +105,13 @@ def _run(
         metadata = []
 
         for k in range(num_qgts):
-
-            # TODO unrolling done by QGT
-            # circuit = self.unroller(circuit)
-
             values = np.asarray(parameter_values[k])
             circuit = circuits[k]
 
+            # TODO unrolling done by QGT
+            translator = TranslateParameterizedGates(self.SUPPORTED_GATES)
+            circuit = translator(circuit)
+
             # TODO parameters is None captured by QGT
             if parameters[k] is None:
                 parameters_ = list(circuit.parameters)
@@ -133,10 +130,10 @@ def _run(
 
             unitaries, paramlist = split(circuit, parameters=parameters_)
 
-            # initialize the phase fix vector and the hessian part ``lis``
+            # initialize the phase fix vector and the hessian part ``metric``
             num_parameters = len(unitaries)
             phase_fixes = np.zeros(num_parameters, dtype=complex)
-            lis = np.zeros((num_parameters, 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))
@@ -146,7 +143,10 @@ def _run(
             psi = chi.copy()
             phi = Statevector.from_int(0, (2,) * circuit.num_qubits)
 
-            # get the analytic gradient of the first unitary
+            # 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)
@@ -158,7 +158,7 @@ def _run(
             if self.phase_fix:
                 phase_fixes[0] = _phasefix_term(chi, grad_coeffs, grad_states)
 
-            lis[0, 0] = _l_term(grad_coeffs, grad_states, 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()
@@ -175,7 +175,7 @@ def _run(
                 grad_states = [phi.evolve(gate.decompose()) for _, gate in deriv]
 
                 # compute the digaonal element L_{j, j}
-                lis[j, j] += _l_term(grad_coeffs, grad_states, grad_coeffs, grad_states)
+                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)):
@@ -193,20 +193,26 @@ def _run(
                     grad_coeffs_mu = [coeff for coeff, _ in deriv]
                     grad_states_mu = [lam.evolve(gate) for _, gate in deriv]
 
-                    lis[i, j] += _l_term(grad_coeffs_mu, grad_states_mu, grad_coeffs, grad_states)
+                    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])
 
-            # stack quantum geometric tensor together and take into account the original
-            # order of parameters
+            # 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_single_parameter(plist[0])]
+                index: param_to_circuit[_extract_parameter(plist[0])]
                 for index, plist in enumerate(paramlist)
             }
 
@@ -216,12 +222,13 @@ def _run(
                 for j in range(num_parameters):
                     jloc = remap[j]
                     if i <= j:
-                        qgt[iloc, jloc] += lis[i, j]
+                        qgt[iloc, jloc] += metric[i, j]
                     else:
-                        qgt[iloc, jloc] += np.conj(lis[j, i])
+                        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 = QFIResult(qgts, metadata, options=None)
@@ -240,30 +247,24 @@ def _to_derivtype(self, qgt):
 def _l_term(coeffs_i, states_i, coeffs_j, states_j):
     return sum(
         sum(
-            np.conj(c_i) * c_j * np.conj(state_i.data).dot(state_j.data)
-            for c_i, state_i in zip(coeffs_i, states_i)
+            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 c_j, state_j in zip(coeffs_j, states_j)
+        for coeff_j, state_j in zip(coeffs_j, states_j)
     )
 
 
 def _phasefix_term(chi, coeffs, states):
-    return sum(c_i * np.conj(chi.data).dot(state_i.data) for c_i, state_i in zip(coeffs, states))
-
+    return sum(
+        coeff_i * np.conj(chi.data).dot(state_i.data) for coeff_i, state_i in zip(coeffs, states)
+    )
 
-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
+def _extract_parameter(expression):
+    if isinstance(expression, Parameter):
+        return expression
 
-    try:
-        spmatrix = operator.to_matrix(sparse=True)
-        evolved = spmatrix @ state.data
-        return Statevector(evolved)
-    except AttributeError:
-        logger.info("Operator is not castable to a sparse matrix, using Statevector.evolve.")
+    if len(expression.parameters) > 1:
+        raise ValueError("Expression has more than one parameter.")
 
-    return state.evolve(operator)
+    return list(expression.parameters)[0]
diff --git a/test/python/algorithms/gradients/test_qfi.py b/test/python/algorithms/gradients/test_qfi.py
index 34b95ac00ceb..7790ec89f46d 100644
--- a/test/python/algorithms/gradients/test_qfi.py
+++ b/test/python/algorithms/gradients/test_qfi.py
@@ -151,13 +151,16 @@ def test_qfi_derivative_type(self, qgt_type):
                 qfi_result = qfi.run([qc], [param]).result().qfis
                 np.testing.assert_allclose(qfi_result[0], correct_values[i], atol=1e-3)
 
-    def test_qfi_coefficients(self):
+    @data(LinCombQFI, ReverseQGT)
+    def test_qfi_coefficients(self, qgt_type):
         """Test the derivative option of QFI"""
         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)
 
-        qfi = LinCombQFI(self.estimator)
+        args = (self.estimator,) if qgt_type != ReverseQGT else ()
+        qfi = qgt_type(*args)
+
         # test imaginary derivative
         param_list = [
             [np.pi / 4 for param in qc.parameters],
@@ -181,19 +184,24 @@ def test_qfi_coefficients(self):
             qfi_result = qfi.run([qc], [param]).result().qfis
             np.testing.assert_allclose(qfi_result[0], correct_values[i], atol=1e-3)
 
-    def test_qfi_specify_parameters(self):
+    @data(LinCombQFI, ReverseQGT)
+    def test_qfi_specify_parameters(self, qgt_type):
         """Test the QFI with specified parameters"""
         a = Parameter("a")
         b = Parameter("b")
         qc = QuantumCircuit(1)
         qc.rx(a, 0)
         qc.ry(b, 0)
-        qfi = LinCombQFI(self.estimator)
+
+        args = (self.estimator,) if qgt_type != ReverseQGT else ()
+        qfi = qgt_type(*args)
+
         param_list = [np.pi / 4, np.pi / 4]
         qfi_result = qfi.run([qc], [param_list], [[a]]).result().qfis
         np.testing.assert_allclose(qfi_result[0], [[1]], atol=1e-3)
 
-    def test_qfi_multi_arguments(self):
+    @data(LinCombQFI, ReverseQGT)
+    def test_qfi_multi_arguments(self, qgt_type):
         """Test the QFI for multiple arguments"""
         a = Parameter("a")
         b = Parameter("b")
@@ -203,7 +211,9 @@ def test_qfi_multi_arguments(self):
         qc2 = QuantumCircuit(1)
         qc2.rx(a, 0)
         qc2.ry(b, 0)
-        qfi = LinCombQFI(self.estimator)
+
+        args = (self.estimator,) if qgt_type != ReverseQGT else ()
+        qfi = qgt_type(*args)
 
         param_list = [[np.pi / 4], [np.pi / 2]]
         correct_values = [

From 88c5f2587447fca3581d4645f222bda2a310b04c Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Tue, 17 Jan 2023 17:13:11 +0100
Subject: [PATCH 09/14] Update after QGT merge

---
 qiskit/algorithms/gradients/lin_comb_qgt.py   | 10 +-
 .../reverse_gradient/derive_circuit.py        |  3 +
 .../reverse_gradient/reverse_gradient.py      |  2 +-
 .../gradients/reverse_gradient/reverse_qgt.py | 88 +++++++-----------
 test/python/algorithms/gradients/test_qfi.py  | 15 ++-
 test/python/algorithms/gradients/test_qgt.py  | 93 +++++++++++--------
 6 files changed, 108 insertions(+), 103 deletions(-)

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/derive_circuit.py b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
index 266b02abaa0d..f97335c7e56b 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
@@ -31,6 +31,9 @@ def gradient_lookup(gate: Gate) -> list[tuple[complex, QuantumCircuit]]:
         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]
diff --git a/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py b/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
index 29af8964d489..9abe2792e178 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
@@ -197,7 +197,7 @@ def _evolve_by_operator(operator, state):
         spmatrix = operator.to_matrix(sparse=True)
         evolved = spmatrix @ state.data
         return Statevector(evolved)
-    except AttributeError:
+    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
index 4a2ca10dac17..e423e0bcbe95 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
@@ -19,24 +19,22 @@
 import numpy as np
 
 from qiskit.circuit import QuantumCircuit, Parameter
-from qiskit.opflow import PauliSumOp
 from qiskit.quantum_info import Statevector
 from qiskit.providers import Options
-from qiskit.transpiler.passes import TranslateParameterizedGates
+from qiskit.primitives import Estimator
 
-from ..base_qfi import BaseQFI
-from ..qfi_result import QFIResult
+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
-from .reverse_gradient import _evolve_by_operator
 
 logger = logging.getLogger(__name__)
 
 
-class ReverseQGT(BaseQFI):
+class ReverseQGT(BaseQGT):
     """QGT calculation with the classically efficient reverse mode.
 
     .. note::
@@ -58,9 +56,8 @@ class ReverseQGT(BaseQFI):
 
     SUPPORTED_GATES = ["rx", "ry", "rz", "cp", "crx", "cry", "crz"]
 
-    # TODO: should the default for QGT be complex?
     def __init__(
-        self, phase_fix: bool = True, derivative_type: DerivativeType = DerivativeType.REAL
+        self, phase_fix: bool = True, derivative_type: DerivativeType = DerivativeType.COMPLEX
     ):
         """
         Args:
@@ -68,9 +65,8 @@ def __init__(
             derivative_type: Determines whether the complex QGT or only the real or imaginary
                 parts are calculated.
         """
-        super().__init__()
-        self.phase_fix = phase_fix
-        self._derivative_type = derivative_type
+        dummy_estimator = Estimator()  # this method does not need an estimator
+        super().__init__(dummy_estimator, phase_fix, derivative_type)
 
     @property
     def options(self) -> Options:
@@ -81,25 +77,27 @@ def options(self) -> Options:
         """
         return Options()
 
-    # TODO this should be in the base class of QGT
-    @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],
+        parameter_sets: Sequence[set[Parameter]],
         **options,
-    ) -> QFIResult:
-        # cast to array, if values are a list
+    ) -> 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 = []
@@ -107,28 +105,13 @@ def _run(
         for k in range(num_qgts):
             values = np.asarray(parameter_values[k])
             circuit = circuits[k]
+            parameters = list(parameter_sets[k])
 
-            # TODO unrolling done by QGT
-            translator = TranslateParameterizedGates(self.SUPPORTED_GATES)
-            circuit = translator(circuit)
-
-            # TODO parameters is None captured by QGT
-            if parameters[k] is None:
-                parameters_ = list(circuit.parameters)
-            else:
-                parameters_ = list(parameters[k])
-
-            num_parameters = len(parameters_)
-            original_parameter_order = [p for p in circuit.parameters if p in parameters_]
-
-            metadata.append(
-                {
-                    "parameters": original_parameter_order,
-                    "derivative_type": self.derivative_type,
-                }
-            )
+            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_)
+            unitaries, paramlist = split(circuit, parameters=parameters)
 
             # initialize the phase fix vector and the hessian part ``metric``
             num_parameters = len(unitaries)
@@ -155,7 +138,7 @@ def _run(
             grad_states = [phi.evolve(gate) for _, gate in deriv]
 
             # compute phase fix (optional) and the hessian part
-            if self.phase_fix:
+            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)
@@ -172,7 +155,7 @@ def _run(
 
                 # compute |phi> (in general it's a sum of states and coeffs)
                 grad_coeffs = [coeff for coeff, _ in deriv]
-                grad_states = [phi.evolve(gate.decompose()) for _, gate 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)
@@ -197,7 +180,7 @@ def _run(
                         grad_coeffs_mu, grad_states_mu, grad_coeffs, grad_states
                     )
 
-                if self.phase_fix:
+                if self._phase_fix:
                     phase_fixes[j] += _phasefix_term(chi, grad_coeffs, grad_states)
 
                 psi = psi.evolve(bound_unitaries[j])
@@ -231,17 +214,16 @@ def _run(
             # append and cast to real/imag if required
             qgts.append(self._to_derivtype(qgt))
 
-        result = QFIResult(qgts, metadata, options=None)
+        result = QGTResult(qgts, self.derivative_type, metadata, options=None)
         return result
 
     def _to_derivtype(self, qgt):
-        # TODO remove factor 4 once the QGT interface is there
         if self.derivative_type == DerivativeType.REAL:
-            return 4 * np.real(qgt)
+            return np.real(qgt)
         if self.derivative_type == DerivativeType.IMAG:
-            return 4 * np.imag(qgt)
+            return np.imag(qgt)
 
-        return 4 * qgt
+        return qgt
 
 
 def _l_term(coeffs_i, states_i, coeffs_j, states_j):
diff --git a/test/python/algorithms/gradients/test_qfi.py b/test/python/algorithms/gradients/test_qfi.py
index 2b0c83e32a08..b0d05ac7030f 100644
--- a/test/python/algorithms/gradients/test_qfi.py
+++ b/test/python/algorithms/gradients/test_qfi.py
@@ -14,24 +14,27 @@
 """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):

From 53bffc22300e3722f7986120f42d9105621cc15c Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Tue, 17 Jan 2023 17:19:19 +0100
Subject: [PATCH 10/14] print which parameter is not in the circuit

---
 qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
index f97335c7e56b..ca9550f1ceef 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
@@ -129,7 +129,7 @@ def derive_circuit(
         raise ValueError(f"parameter must be None or of type Parameter, not {type(parameter)}.")
 
     if parameter not in circuit.parameters:
-        raise ValueError("The parameter is not in this circuit.")
+        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.")

From 33587723c487f63677a51ce568d54adac98ffdef Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Tue, 17 Jan 2023 17:33:16 +0100
Subject: [PATCH 11/14] Fix copyright

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>
---
 qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
index e423e0bcbe95..60e8c3cde2d6 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2022.
+# (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

From 0d4f7ffcf2be6c4783a272d11393c2fab16e9a2f Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Tue, 17 Jan 2023 17:34:09 +0100
Subject: [PATCH 12/14] ``None`` is not actually supported

---
 qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
index ca9550f1ceef..e31db92265d4 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
@@ -126,7 +126,7 @@ def derive_circuit(
     # 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 None or of type Parameter, not {type(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.")

From 6fac2ec2ad64ace48863db82f6f18caf8297e1b3 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Tue, 17 Jan 2023 18:56:45 +0100
Subject: [PATCH 13/14] only setter of derivative_type in LCU/Rev

---
 .../gradients/base_estimator_gradient.py      | 19 +++++++++++++++----
 .../gradients/lin_comb_estimator_gradient.py  | 14 ++++++--------
 .../reverse_gradient/reverse_gradient.py      | 13 +++++--------
 3 files changed, 26 insertions(+), 20 deletions(-)

diff --git a/qiskit/algorithms/gradients/base_estimator_gradient.py b/qiskit/algorithms/gradients/base_estimator_gradient.py
index 67b5d5328693..6ceaa41a51aa 100644
--- a/qiskit/algorithms/gradients/base_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/base_estimator_gradient.py
@@ -48,19 +48,32 @@ 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
@@ -70,9 +83,7 @@ def derivative_type(self) -> DerivativeType:
         Returns:
             The derivative type.
         """
-        # the default case is real, as this yields e.g. the energy gradient and this type
-        # is also supported by function-level schemes like finite difference or SPSA
-        return DerivativeType.REAL
+        return self._derivative_type
 
     def run(
         self,
diff --git a/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py b/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
index e2047e9cf3ec..7c0e13ce19de 100644
--- a/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
+++ b/qiskit/algorithms/gradients/lin_comb_estimator_gradient.py
@@ -87,15 +87,9 @@ def __init__(
                 Higher priority setting overrides lower priority setting.
         """
         self._lin_comb_cache = {}
-        self._derivative_type = derivative_type
-        super().__init__(estimator, options)
-
-    @property
-    def derivative_type(self) -> DerivativeType:
-        """The derivative type."""
-        return self._derivative_type
+        super().__init__(estimator, options, derivative_type=derivative_type)
 
-    @derivative_type.setter
+    @BaseEstimatorGradient.derivative_type.setter
     def derivative_type(self, derivative_type: DerivativeType) -> None:
         """Set the derivative type."""
         self._derivative_type = derivative_type
@@ -183,10 +177,14 @@ def _run_unique(
         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
diff --git a/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py b/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
index 9abe2792e178..004d09f1564c 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/reverse_gradient.py
@@ -65,15 +65,9 @@ def __init__(self, derivative_type: DerivativeType = DerivativeType.REAL):
                 of the gradient is returned.
         """
         dummy_estimator = Estimator()  # this is required by the base class, but not used
-        super().__init__(dummy_estimator)
-        self._derivative_type = derivative_type
-
-    @property
-    def derivative_type(self) -> DerivativeType:
-        """The derivative type."""
-        return self._derivative_type
+        super().__init__(dummy_estimator, derivative_type=derivative_type)
 
-    @derivative_type.setter
+    @BaseEstimatorGradient.derivative_type.setter
     def derivative_type(self, derivative_type: DerivativeType) -> None:
         """Set the derivative type."""
         self._derivative_type = derivative_type
@@ -177,6 +171,9 @@ def _run_unique(
         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:

From c7954bfbdad7d9b48a7e05e9685cf82ab77ae359 Mon Sep 17 00:00:00 2001
From: Julien Gacon <gaconju@gmail.com>
Date: Tue, 17 Jan 2023 21:42:10 +0100
Subject: [PATCH 14/14] Update copyrights

Co-authored-by: Steve Wood <40241007+woodsp-ibm@users.noreply.github.com>
---
 qiskit/algorithms/gradients/reverse_gradient/bind.py           | 2 +-
 qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/qiskit/algorithms/gradients/reverse_gradient/bind.py b/qiskit/algorithms/gradients/reverse_gradient/bind.py
index 0236a30ff7e0..5380090b7425 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/bind.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/bind.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2022.
+# (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
diff --git a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
index e31db92265d4..0932a5100b78 100644
--- a/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
+++ b/qiskit/algorithms/gradients/reverse_gradient/derive_circuit.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2022.
+# (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