Efficient classical calculation of expectation gradients

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",
 ]

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.

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]

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,43 @@
+# 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 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`.
+    """
+    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,142 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Split a circuit into subcircuits, each containing a single parameterized gate."""
+
+from __future__ import annotations
+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."""
+
+    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("parameter  must be  None or of type Parameter.")
+
+    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.")
+
+    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