Efficient classical calculation of expectation gradients (Qiskit/qisk…

…it#9287)

* ckassically efficient gradients

* cleanup & reno

* Apply suggestions from code review

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>

* Remove support for Parameter = None

* Complete the docs

* QGT v0

* fix LCU tests

* final fixes

* Update after QGT merge

* print which parameter is not in the circuit

* Fix copyright

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>

* ``None`` is not actually supported

* only setter of derivative_type in LCU/Rev

* Update copyrights

Co-authored-by: Steve Wood <40241007+woodsp-ibm@users.noreply.github.com>

Co-authored-by: ElePT <57907331+ElePT@users.noreply.github.com>
Co-authored-by: Steve Wood <40241007+woodsp-ibm@users.noreply.github.com>
Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com>

qiskit_algorithms/gradients/__init__.py

@@ -30,6 +30,7 @@
    LinCombEstimatorGradient
    ParamShiftEstimatorGradient
    SPSAEstimatorGradient
+   ReverseEstimatorGradient
 
 Sampler Gradients
 =================
@@ -51,6 +52,7 @@
    BaseQGT
    LinCombQGT
    QFI
+   ReverseQGT
 
 Results
 =======
@@ -81,6 +83,8 @@
 from .sampler_gradient_result import SamplerGradientResult
 from .spsa_estimator_gradient import SPSAEstimatorGradient
 from .spsa_sampler_gradient import SPSASamplerGradient
+from .reverse_gradient.reverse_gradient import ReverseEstimatorGradient
+from .reverse_gradient.reverse_qgt import ReverseQGT
 
 __all__ = [
     "BaseEstimatorGradient",
@@ -101,4 +105,6 @@
     "SamplerGradientResult",
     "SPSAEstimatorGradient",
     "SPSASamplerGradient",
+    "ReverseEstimatorGradient",
+    "ReverseQGT",
 ]

qiskit_algorithms/gradients/base_estimator_gradient.py

@@ -48,21 +48,43 @@ def __init__(
         self,
         estimator: BaseEstimator,
         options: Options | None = None,
+        derivative_type: DerivativeType = DerivativeType.REAL,
     ):
-        """
+        r"""
         Args:
             estimator: The estimator used to compute the gradients.
             options: Primitive backend runtime options used for circuit execution.
                 The order of priority is: options in ``run`` method > gradient's
                 default options > primitive's default setting.
                 Higher priority setting overrides lower priority setting
+            derivative_type: The type of derivative. Can be either ``DerivativeType.REAL``
+                ``DerivativeType.IMAG``, or ``DerivativeType.COMPLEX``.
+
+                    - ``DerivativeType.REAL`` computes :math:`2 \mathrm{Re}[⟨ψ(ω)|O(θ)|dω ψ(ω)〉]`.
+                    - ``DerivativeType.IMAG`` computes :math:`2 \mathrm{Im}[⟨ψ(ω)|O(θ)|dω ψ(ω)〉]`.
+                    - ``DerivativeType.COMPLEX`` computes :math:`2 ⟨ψ(ω)|O(θ)|dω ψ(ω)〉`.
+
+                Defaults to ``DerivativeType.REAL``, as this yields e.g. the commonly-used energy
+                gradient and this type is the only supported type for function-level schemes like
+                finite difference.
         """
         self._estimator: BaseEstimator = estimator
         self._default_options = Options()
         if options is not None:
             self._default_options.update_options(**options)
+        self._derivative_type = derivative_type
+
         self._gradient_circuit_cache: dict[QuantumCircuit, GradientCircuit] = {}
 
+    @property
+    def derivative_type(self) -> DerivativeType:
+        """Return the derivative type (real, imaginary or complex).
+
+        Returns:
+            The derivative type.
+        """
+        return self._derivative_type
+
     def run(
         self,
         circuits: Sequence[QuantumCircuit],

qiskit_algorithms/gradients/lin_comb_estimator_gradient.py

@@ -87,8 +87,12 @@ def __init__(
                 Higher priority setting overrides lower priority setting.
         """
         self._lin_comb_cache = {}
+        super().__init__(estimator, options, derivative_type=derivative_type)
+
+    @BaseEstimatorGradient.derivative_type.setter
+    def derivative_type(self, derivative_type: DerivativeType) -> None:
+        """Set the derivative type."""
         self._derivative_type = derivative_type
-        super().__init__(estimator, options)
 
     def _run(
         self,
@@ -143,7 +147,7 @@ def _run_unique(
             )
             # If its derivative type is `DerivativeType.COMPLEX`, calculate the gradient
             # of the real and imaginary parts separately.
-            meta["derivative_type"] = self._derivative_type
+            meta["derivative_type"] = self.derivative_type
             metadata.append(meta)
             # Combine inputs into a single job to reduce overhead.
             if self._derivative_type == DerivativeType.COMPLEX:
@@ -173,10 +177,14 @@ def _run_unique(
         gradients = []
         partial_sum_n = 0
         for n in all_n:
-            if self._derivative_type == DerivativeType.COMPLEX:
+            # this disable is needed as Pylint does not understand derivative_type is a property if
+            # it is only defined in the base class and the getter is in the child
+            # pylint: disable=comparison-with-callable
+            if self.derivative_type == DerivativeType.COMPLEX:
                 gradient = np.zeros(n // 2, dtype="complex")
                 gradient.real = results.values[partial_sum_n : partial_sum_n + n // 2]
                 gradient.imag = results.values[partial_sum_n + n // 2 : partial_sum_n + n]
+
             else:
                 gradient = np.real(results.values[partial_sum_n : partial_sum_n + n])
             partial_sum_n += n

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 = [

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.

qiskit_algorithms/gradients/reverse_gradient/bind.py

@@ -0,0 +1,53 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022, 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Bind values to a parametrized circuit, accepting binds for non-existing parameters in the circuit."""
+
+from __future__ import annotations
+from collections.abc import Iterable
+
+from qiskit.circuit import QuantumCircuit, Parameter
+
+# pylint: disable=inconsistent-return-statements
+def bind(
+    circuits: QuantumCircuit | Iterable[QuantumCircuit],
+    parameter_binds: dict[Parameter, float],
+    inplace: bool = False,
+) -> QuantumCircuit | Iterable[QuantumCircuit] | None:
+    """Bind parameters in a circuit (or list of circuits).
+
+    This method also allows passing parameter binds to parameters that are not in the circuit,
+    and thereby differs to :meth:`.QuantumCircuit.bind_parameters`.
+
+    Args:
+        circuits: Input circuit(s).
+        parameter_binds: A dictionary with ``{Parameter: float}`` pairs determining the values to
+            which the free parameters in the circuit(s) are bound.
+        inplace: If ``True``, bind the values in place, otherwise return circuit copies.
+
+    Returns:
+        The bound circuits, if ``inplace=False``, otherwise None.
+
+    """
+    if not isinstance(circuits, Iterable):
+        circuits = [circuits]
+        return_list = False
+    else:
+        return_list = True
+
+    bound = []
+    for circuit in circuits:
+        existing_parameter_binds = {p: parameter_binds[p] for p in circuit.parameters}
+        bound.append(circuit.assign_parameters(existing_parameter_binds, inplace=inplace))
+
+    if not inplace:
+        return bound if return_list else bound[0]

qiskit_algorithms/gradients/reverse_gradient/derive_circuit.py

@@ -0,0 +1,156 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022, 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Split a circuit into subcircuits, each containing a single parameterized gate."""
+
+from __future__ import annotations
+import itertools
+
+from qiskit.circuit import QuantumCircuit, Parameter, Gate
+from qiskit.circuit.library import RXGate, RYGate, RZGate, CRXGate, CRYGate, CRZGate
+
+
+def gradient_lookup(gate: Gate) -> list[tuple[complex, QuantumCircuit]]:
+    """Returns a circuit implementing the gradient of the input gate.
+
+    Args:
+        gate: The gate whose derivative is returned.
+
+    Returns:
+        The derivative of the input gate as list of ``(coeff, circuit)`` pairs,
+        where the sum of all ``coeff * circuit`` elements describes the full derivative.
+        The circuit is the unitary part of the derivative with a potential separate ``coeff``.
+        The output is a list as derivatives of e.g. controlled gates can only be described
+        as a sum of ``coeff * circuit`` pairs.
+
+    Raises:
+        NotImplementedError: If the derivative of ``gate`` is not implemented.
+    """
+
+    param = gate.params[0]
+    if isinstance(gate, RXGate):
+        derivative = QuantumCircuit(gate.num_qubits)
+        derivative.rx(param, 0)
+        derivative.x(0)
+        return [(-0.5j, derivative)]
+    if isinstance(gate, RYGate):
+        derivative = QuantumCircuit(gate.num_qubits)
+        derivative.ry(param, 0)
+        derivative.y(0)
+        return [(-0.5j, derivative)]
+    if isinstance(gate, RZGate):
+        derivative = QuantumCircuit(gate.num_qubits)
+        derivative.rz(param, 0)
+        derivative.z(0)
+        return [(-0.5j, derivative)]
+    if isinstance(gate, CRXGate):
+        proj1 = QuantumCircuit(gate.num_qubits)
+        proj1.rx(param, 1)
+        proj1.x(1)
+
+        proj2 = QuantumCircuit(gate.num_qubits)
+        proj2.z(0)
+        proj2.rx(param, 1)
+        proj2.x(1)
+
+        return [(-0.25j, proj1), (0.25j, proj2)]
+    if isinstance(gate, CRYGate):
+        proj1 = QuantumCircuit(gate.num_qubits)
+        proj1.ry(param, 1)
+        proj1.y(1)
+
+        proj2 = QuantumCircuit(gate.num_qubits)
+        proj2.z(0)
+        proj2.ry(param, 1)
+        proj2.y(1)
+
+        return [(-0.25j, proj1), (0.25j, proj2)]
+    if isinstance(gate, CRZGate):
+        proj1 = QuantumCircuit(gate.num_qubits)
+        proj1.rz(param, 1)
+        proj1.z(1)
+
+        proj2 = QuantumCircuit(gate.num_qubits)
+        proj2.z(0)
+        proj2.rz(param, 1)
+        proj2.z(1)
+
+        return [(-0.25j, proj1), (0.25j, proj2)]
+    raise NotImplementedError("Cannot implement gradient for", gate)
+
+
+def derive_circuit(
+    circuit: QuantumCircuit, parameter: Parameter
+) -> list[tuple[complex, QuantumCircuit]]:
+    """Return the analytic gradient expression of the input circuit wrt. a single parameter.
+
+    Returns a list of ``(coeff, gradient_circuit)`` tuples, where the derivative of the circuit is
+    given by the sum of the gradient circuits multiplied by their coefficient.
+
+    For example, the circuit::
+
+           ┌───┐┌───────┐┌─────┐
+        q: ┤ H ├┤ Rx(x) ├┤ Sdg ├
+           └───┘└───────┘└─────┘
+
+    returns the coefficient `-0.5j` and the circuit equivalent to::
+
+           ┌───┐┌───────┐┌───┐┌─────┐
+        q: ┤ H ├┤ Rx(x) ├┤ X ├┤ Sdg ├
+           └───┘└───────┘└───┘└─────┘
+
+    as the derivative of `Rx(x)` is `-0.5j Rx(x) X`.
+
+    Args:
+        circuit: The quantum circuit to derive.
+        parameter: The parameter with respect to which we derive.
+
+    Returns:
+        A list of ``(coeff, gradient_circuit)`` tuples.
+
+    Raises:
+        ValueError: If ``parameter`` is of the wrong type.
+        ValueError: If ``parameter`` is not in this circuit.
+        NotImplementedError: If a non-unique parameter is added, as the product rule is not yet
+            supported in this function.
+    """
+    # this is added as useful user-warning, since sometimes ``ParameterExpression``s are
+    # passed around instead of ``Parameter``s
+    if not isinstance(parameter, Parameter):
+        raise ValueError(f"parameter must be of type Parameter, not {type(parameter)}.")
+
+    if parameter not in circuit.parameters:
+        raise ValueError(f"The parameter {parameter} is not in this circuit.")
+
+    if len(circuit._parameter_table[parameter]) > 1:
+        raise NotImplementedError("No product rule support yet, circuit parameters must be unique.")
+
+    summands, op_context = [], []
+    for i, op in enumerate(circuit.data):
+        gate = op[0]
+        op_context += [op[1:]]
+        if parameter in gate.params:
+            coeffs_and_grads = gradient_lookup(gate)
+            summands += [coeffs_and_grads]
+        else:
+            summands += [[(1, gate)]]
+
+    gradient = []
+    for product_rule_term in itertools.product(*summands):
+        summand_circuit = QuantumCircuit(*circuit.qregs)
+        c = 1
+        for i, term in enumerate(product_rule_term):
+            c *= term[0]
+            summand_circuit.data.append([term[1], *op_context[i]])
+        gradient += [(c, summand_circuit.copy())]
+
+    return gradient