Efficient classical calculation of expectation gradients

qiskit/algorithms/gradients/__init__.py

@@ -52,6 +52,7 @@
    BaseQGT
    LinCombQGT
    QFI
+   ReverseQGT
 
 Results
 =======
@@ -83,6 +84,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",
@@ -104,4 +106,5 @@
     "SPSAEstimatorGradient",
     "SPSASamplerGradient",
     "ReverseEstimatorGradient",
+    "ReverseQGT",
 ]

qiskit/algorithms/gradients/base_estimator_gradient.py

@@ -212,9 +212,13 @@ def _postprocess(
         for idx, (circuit, parameter_values_, parameter_set) in enumerate(
             zip(circuits, parameter_values, parameter_sets)
         ):
-            dtype = complex if self.derivative_type == DerivativeType.COMPLEX else float
-            unique_gradient = np.zeros(len(parameter_set), dtype=dtype)
-
+            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.
+                gradient = gradient.astype("complex")
             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_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/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]

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)

qiskit/algorithms/gradients/reverse_gradient/reverse_qgt.py

@@ -0,0 +1,252 @@
+# 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 import Statevector
+from qiskit.providers import Options
+from qiskit.primitives import Estimator
+
+from ..base_qgt import BaseQGT
+from ..qgt_result import QGTResult
+from ..utils import DerivativeType
+
+from .split_circuits import split
+from .bind import bind
+from .derive_circuit import derive_circuit
+
+logger = logging.getLogger(__name__)
+
+
+class ReverseQGT(BaseQGT):
+    """QGT calculation with the classically efficient reverse mode.
+
+    .. note::
+
+        This QGT implementation is based on statevector manipulations and scales exponentially
+        with the number of qubits. However, for small system sizes it can be very fast
+        compared to circuit-based gradients.
+
+    This class implements the calculation of the QGT as described in [1].
+    By keeping track of three statevectors and iteratively sweeping through each parameterized
+    gate, this method scales only quadratically with the number of parameters.
+
+    **References:**
+
+        [1]: Jones, T. "Efficient classical calculation of the Quantum Natural Gradient" (2020).
+             `arXiv:2011.02991 <https://arxiv.org/abs/2011.02991>`_.
+
+    """
+
+    SUPPORTED_GATES = ["rx", "ry", "rz", "cp", "crx", "cry", "crz"]
+
+    def __init__(
+        self, phase_fix: bool = True, derivative_type: DerivativeType = DerivativeType.COMPLEX
+    ):
+        """
+        Args:
+            phase_fix: Whether or not to include the phase fix.
+            derivative_type: Determines whether the complex QGT or only the real or imaginary
+                parts are calculated.
+        """
+        dummy_estimator = Estimator()  # this method does not need an estimator
+        super().__init__(dummy_estimator, phase_fix, derivative_type)
+
+    @property
+    def options(self) -> Options:
+        """There are no options for the reverse QGT, returns an empty options dict.
+
+        Returns:
+            Empty options.
+        """
+        return Options()
+
+    def _run(
+        self,
+        circuits: Sequence[QuantumCircuit],
+        parameter_values: Sequence[Sequence[float]],
+        parameter_sets: Sequence[set[Parameter]],
+        **options,
+    ) -> QGTResult:
+        """Compute the QGT on the given circuits."""
+        g_circuits, g_parameter_values, g_parameter_sets = self._preprocess(
+            circuits, parameter_values, parameter_sets, self.SUPPORTED_GATES
+        )
+        results = self._run_unique(g_circuits, g_parameter_values, g_parameter_sets, **options)
+        return self._postprocess(results, circuits, parameter_values, parameter_sets)
+
+    def _run_unique(
+        self,
+        circuits: Sequence[QuantumCircuit],
+        parameter_values: Sequence[Sequence[float]],
+        parameter_sets: Sequence[set[Parameter]],
+        **options,  # pylint: disable=unused-argument
+    ) -> QGTResult:
+        num_qgts = len(circuits)
+        qgts = []
+        metadata = []
+
+        for k in range(num_qgts):
+            values = np.asarray(parameter_values[k])
+            circuit = circuits[k]
+            parameters = list(parameter_sets[k])
+
+            num_parameters = len(parameters)
+            original_parameter_order = [p for p in circuit.parameters if p in parameters]
+            metadata.append({"parameters": original_parameter_order})
+
+            unitaries, paramlist = split(circuit, parameters=parameters)
+
+            # initialize the phase fix vector and the hessian part ``metric``
+            num_parameters = len(unitaries)
+            phase_fixes = np.zeros(num_parameters, dtype=complex)
+            metric = np.zeros((num_parameters, num_parameters), dtype=complex)
+
+            # initialize the state variables -- naming convention is the same as the paper
+            parameter_binds = dict(zip(circuit.parameters, values))
+            bound_unitaries = bind(unitaries, parameter_binds)
+
+            chi = Statevector(bound_unitaries[0])
+            psi = chi.copy()
+            phi = Statevector.from_int(0, (2,) * circuit.num_qubits)
+
+            # Get the analytic gradient of the first unitary
+            # Note: We currently only support gates with a single parameter -- which is reflected
+            # in self.SUPPORTED_GATES -- but generally we could also support gates with multiple
+            # parameters per gate. This is the reason for the second 0-index.
+            deriv = derive_circuit(unitaries[0], paramlist[0][0])
+            for _, gate in deriv:
+                bind(gate, parameter_binds, inplace=True)
+
+            grad_coeffs = [coeff for coeff, _ in deriv]
+            grad_states = [phi.evolve(gate) for _, gate in deriv]
+
+            # compute phase fix (optional) and the hessian part
+            if self._phase_fix:
+                phase_fixes[0] = _phasefix_term(chi, grad_coeffs, grad_states)
+
+            metric[0, 0] = _l_term(grad_coeffs, grad_states, grad_coeffs, grad_states)
+
+            for j in range(1, num_parameters):
+                lam = psi.copy()
+                phi = psi.copy()
+
+                # get the analytic gradient d U_j / d p_j and apply it
+                deriv = derive_circuit(unitaries[j], paramlist[j][0])
+
+                for _, gate in deriv:
+                    bind(gate, parameter_binds, inplace=True)
+
+                # compute |phi> (in general it's a sum of states and coeffs)
+                grad_coeffs = [coeff for coeff, _ in deriv]
+                grad_states = [phi.evolve(gate) for _, gate in deriv]
+
+                # compute the digaonal element L_{j, j}
+                metric[j, j] += _l_term(grad_coeffs, grad_states, grad_coeffs, grad_states)
+
+                # compute the off diagonal elements L_{i, j}
+                for i in reversed(range(j)):
+                    # apply U_{i + 1}_dg
+                    unitary_ip_inv = bound_unitaries[i + 1].inverse()
+                    grad_states = [state.evolve(unitary_ip_inv) for state in grad_states]
+
+                    lam = lam.evolve(bound_unitaries[i].inverse())
+
+                    # get the gradient d U_i / d p_i and apply it
+                    deriv = derive_circuit(unitaries[i], paramlist[i][0])
+                    for _, gate in deriv:
+                        bind(gate, parameter_binds, inplace=True)
+
+                    grad_coeffs_mu = [coeff for coeff, _ in deriv]
+                    grad_states_mu = [lam.evolve(gate) for _, gate in deriv]
+
+                    metric[i, j] += _l_term(
+                        grad_coeffs_mu, grad_states_mu, grad_coeffs, grad_states
+                    )
+
+                if self._phase_fix:
+                    phase_fixes[j] += _phasefix_term(chi, grad_coeffs, grad_states)
+
+                psi = psi.evolve(bound_unitaries[j])
+
+            # The following code stacks the QGT together and maps the values into the
+            # correct original order of parameters
+
+            # map circuit parameter to global index in the circuit
+            param_to_circuit = {
+                param: index for index, param in enumerate(original_parameter_order)
+            }
+            # map global index to the local index used in the calculation, the new index can
+            # now be accessed by remap[index]
+            remap = {
+                index: param_to_circuit[_extract_parameter(plist[0])]
+                for index, plist in enumerate(paramlist)
+            }
+
+            qgt = np.zeros((num_parameters, num_parameters), dtype=complex)
+            for i in range(num_parameters):
+                iloc = remap[i]
+                for j in range(num_parameters):
+                    jloc = remap[j]
+                    if i <= j:
+                        qgt[iloc, jloc] += metric[i, j]
+                    else:
+                        qgt[iloc, jloc] += np.conj(metric[j, i])
+
+                    qgt[iloc, jloc] -= np.conj(phase_fixes[i]) * phase_fixes[j]
+
+            # append and cast to real/imag if required
+            qgts.append(self._to_derivtype(qgt))
+
+        result = QGTResult(qgts, self.derivative_type, metadata, options=None)
+        return result
+
+    def _to_derivtype(self, qgt):
+        if self.derivative_type == DerivativeType.REAL:
+            return np.real(qgt)
+        if self.derivative_type == DerivativeType.IMAG:
+            return np.imag(qgt)
+
+        return qgt
+
+
+def _l_term(coeffs_i, states_i, coeffs_j, states_j):
+    return sum(
+        sum(
+            np.conj(coeff_i) * coeff_j * np.conj(state_i.data).dot(state_j.data)
+            for coeff_i, state_i in zip(coeffs_i, states_i)
+        )
+        for coeff_j, state_j in zip(coeffs_j, states_j)
+    )
+
+
+def _phasefix_term(chi, coeffs, states):
+    return sum(
+        coeff_i * np.conj(chi.data).dot(state_i.data) for coeff_i, state_i in zip(coeffs, states)
+    )
+
+
+def _extract_parameter(expression):
+    if isinstance(expression, Parameter):
+        return expression
+
+    if len(expression.parameters) > 1:
+        raise ValueError("Expression has more than one parameter.")
+
+    return list(expression.parameters)[0]

test/python/algorithms/gradients/test_qfi.py

@@ -11,27 +11,30 @@
 # that they have been altered from the originals.
 # =============================================================================
 
-""" Test QFI"""
+"""Test QFI."""
 
 import unittest
+from ddt import ddt, data
 
 import numpy as np
 
 from qiskit import QuantumCircuit
-from qiskit.algorithms.gradients import LinCombQGT, QFI
+from qiskit.algorithms.gradients import LinCombQGT, ReverseQGT, QFI, DerivativeType
 from qiskit.circuit import Parameter
 from qiskit.circuit.parametervector import ParameterVector
 from qiskit.primitives import Estimator
 from qiskit.test import QiskitTestCase
 
 
+@ddt
 class TestQFI(QiskitTestCase):
     """Test QFI"""
 
     def setUp(self):
         super().setUp()
         self.estimator = Estimator()
-        self.qgt = LinCombQGT(self.estimator)
+        self.lcu_qgt = LinCombQGT(self.estimator, derivative_type=DerivativeType.REAL)
+        self.reverse_qgt = ReverseQGT(derivative_type=DerivativeType.REAL)
 
     def test_qfi(self):
         """Test if the quantum fisher information calculation is correct for a simple test case.
@@ -47,7 +50,7 @@ def test_qfi(self):
         param_list = [[np.pi / 4, 0.1], [np.pi, 0.1], [np.pi / 2, 0.1]]
         correct_values = [[[1, 0], [0, 0.5]], [[1, 0], [0, 0]], [[1, 0], [0, 1]]]
 
-        qfi = QFI(self.qgt)
+        qfi = QFI(self.lcu_qgt)
         for i, param in enumerate(param_list):
             qfis = qfi.run([qc], [param]).result().qfis
             np.testing.assert_allclose(qfis[0], correct_values[i], atol=1e-3)
@@ -69,7 +72,8 @@ def test_qfi_phase_fix(self):
         qfis = qfi.run([qc], [param]).result().qfis
         np.testing.assert_allclose(qfis[0], correct_values, atol=1e-3)
 
-    def test_qfi_maxcut(self):
+    @data("lcu", "reverse")
+    def test_qfi_maxcut(self, qgt_kind):
         """Test the QFI for a simple MaxCut problem.
 
         This is interesting because it contains the same parameters in different gates.
@@ -101,7 +105,8 @@ def expiz(qubit0, qubit1):
         reference = np.array([[16.0, -5.551], [-5.551, 18.497]])
         param = [0.4, 0.69]
 
-        qfi = QFI(self.qgt)
+        qgt = self.lcu_qgt if qgt_kind == "lcu" else self.reverse_qgt
+        qfi = QFI(qgt)
         qfi_result = qfi.run([ansatz], [param]).result().qfis
         np.testing.assert_array_almost_equal(qfi_result[0], reference, decimal=3)