Sampler-based VQE for diagonal operators, plus QAOA

qiskit/algorithms/minimum_eigensolvers/__init__.py

@@ -26,6 +26,9 @@
    NumPyMinimumEigensolver
    VQE
    AdaptVQE
+   SamplingMinimumEigensolver
+   SamplingVQE
+   QAOA
 
 .. autosummary::
    :toctree: ../stubs/
@@ -34,12 +37,17 @@
    NumPyMinimumEigensolverResult
    VQEResult
    AdaptVQEResult
+   SamplingMinimumEigensolverResult
+   SamplingVQEResult
 """
 
 from .adapt_vqe import AdaptVQE, AdaptVQEResult
 from .minimum_eigensolver import MinimumEigensolver, MinimumEigensolverResult
 from .numpy_minimum_eigensolver import NumPyMinimumEigensolver, NumPyMinimumEigensolverResult
 from .vqe import VQE, VQEResult
+from .sampling_mes import SamplingMinimumEigensolver, SamplingMinimumEigensolverResult
+from .sampling_vqe import SamplingVQE, SamplingVQEResult
+from .qaoa import QAOA
 
 __all__ = [
     "AdaptVQE",
@@ -50,4 +58,9 @@
     "NumPyMinimumEigensolverResult",
     "VQE",
     "VQEResult",
+    "SamplingMinimumEigensolver",
+    "SamplingMinimumEigensolverResult",
+    "SamplingVQE",
+    "SamplingVQEResult",
+    "QAOA",
 ]

qiskit/algorithms/minimum_eigensolvers/diagonal_estimator.py

@@ -0,0 +1,202 @@
+# 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.
+
+"""Expectation value for a diagonal observable using a sampler primitive."""
+
+from __future__ import annotations
+
+from collections.abc import Callable, Sequence, Mapping
+from typing import Any
+
+from dataclasses import dataclass
+
+import numpy as np
+from qiskit.algorithms.algorithm_job import AlgorithmJob
+from qiskit.circuit import QuantumCircuit
+from qiskit.primitives import BaseSampler, BaseEstimator, EstimatorResult
+from qiskit.primitives.utils import init_observable, _circuit_key
+from qiskit.opflow import PauliSumOp
+from qiskit.quantum_info import SparsePauliOp
+from qiskit.quantum_info.operators.base_operator import BaseOperator
+
+
+@dataclass(frozen=True)
+class _DiagonalEstimatorResult(EstimatorResult):
+    """A result from an expectation of a diagonal observable."""
+
+    # TODO make each measurement a dataclass rather than a dict
+    best_measurements: Sequence[Mapping[str, Any]] | None = None
+
+
+class _DiagonalEstimator(BaseEstimator):
+    """An estimator for diagonal observables."""
+
+    def __init__(
+        self,
+        sampler: BaseSampler,
+        aggregation: float | Callable[[Sequence[tuple[float, float]]], float] | None = None,
+        callback: Callable[[Sequence[Mapping[str, Any]]], None] | None = None,
+        **options,
+    ) -> None:
+        r"""Evaluate the expectation of quantum state with respect to a diagonal operator.
+
+        Args:
+            sampler: The sampler used to evaluate the circuits.
+            aggregation: The aggregation function to aggregate the measurement outcomes. If a float
+                this specified the CVaR :math:`\alpha` parameter.
+            callback: A callback which is given the best measurements of all circuits in each
+                evaluation.
+            run_options: Options for the sampler.
+
+        """
+        super().__init__(options=options)
+        self.sampler = sampler
+        if not callable(aggregation):
+            aggregation = _get_cvar_aggregation(aggregation)
+
+        self.aggregation = aggregation
+        self.callback = callback
+
+    def _run(
+        self,
+        circuits: Sequence[QuantumCircuit],
+        observables: Sequence[BaseOperator | PauliSumOp],
+        parameter_values: Sequence[Sequence[float]],
+        **run_options,
+    ) -> AlgorithmJob:
+        circuit_indices = []
+        for circuit in circuits:
+            key = _circuit_key(circuit)
+            index = self._circuit_ids.get(key)
+            if index is not None:
+                circuit_indices.append(index)
+            else:
+                circuit_indices.append(len(self._circuits))
+                self._circuit_ids[key] = len(self._circuits)
+                self._circuits.append(circuit)
+                self._parameters.append(circuit.parameters)
+        observable_indices = []
+        for observable in observables:
+            index = self._observable_ids.get(id(observable))
+            if index is not None:
+                observable_indices.append(index)
+            else:
+                observable_indices.append(len(self._observables))
+                self._observable_ids[id(observable)] = len(self._observables)
+                converted_observable = init_observable(observable)
+                _check_observable_is_diagonal(converted_observable)  # check it's diagonal
+                self._observables.append(converted_observable)
+        job = AlgorithmJob(
+            self._call, circuit_indices, observable_indices, parameter_values, **run_options
+        )
+        job.submit()
+        return job
+
+    def _call(
+        self,
+        circuits: Sequence[int],
+        observables: Sequence[int],
+        parameter_values: Sequence[Sequence[float]],
+        **run_options,
+    ) -> _DiagonalEstimatorResult:
+        job = self.sampler.run(
+            [self._circuits[i] for i in circuits],
+            parameter_values,
+            **run_options,
+        )
+        sampler_result = job.result()
+        samples = sampler_result.quasi_dists
+
+        # a list of dictionaries containing: {state: (measurement probability, value)}
+        evaluations = [
+            {
+                state: (probability, _evaluate_sparsepauli(state, self._observables[i]))
+                for state, probability in sampled.items()
+            }
+            for i, sampled in zip(observables, samples)
+        ]
+
+        results = np.array([self.aggregation(evaluated.values()) for evaluated in evaluations])
+
+        # get the best measurements
+        best_measurements = []
+        num_qubits = self._circuits[0].num_qubits
+        for evaluated in evaluations:
+            best_result = min(evaluated.items(), key=lambda x: x[1][1])
+            best_measurements.append(
+                {
+                    "state": best_result[0],
+                    "bitstring": bin(best_result[0])[2:].zfill(num_qubits),
+                    "value": best_result[1][1],
+                    "probability": best_result[1][0],
+                }
+            )
+
+        if self.callback is not None:
+            self.callback(best_measurements)
+
+        return _DiagonalEstimatorResult(
+            values=results, metadata=sampler_result.metadata, best_measurements=best_measurements
+        )
+
+
+def _get_cvar_aggregation(alpha):
+    """Get the aggregation function for CVaR with confidence level ``alpha``."""
+    if alpha is None:
+        alpha = 1
+    elif not 0 <= alpha <= 1:
+        raise ValueError(f"alpha must be in [0, 1] but was {alpha}")
+
+    # if alpha is close to 1 we can avoid the sorting
+    if np.isclose(alpha, 1):
+
+        def aggregate(measurements):
+            return sum(probability * value for probability, value in measurements)
+
+    else:
+
+        def aggregate(measurements):
+            # sort by values
+            sorted_measurements = sorted(measurements, key=lambda x: x[1])
+
+            accumulated_percent = 0  # once alpha is reached, stop
+            cvar = 0
+            for probability, value in sorted_measurements:
+                cvar += value * max(probability, alpha - accumulated_percent)
+                accumulated_percent += probability
+                if accumulated_percent >= alpha:
+                    break
+
+            return cvar / alpha
+
+    return aggregate
+
+
+def _evaluate_sparsepauli(state: int, observable: SparsePauliOp) -> complex:
+    return sum(
+        coeff * _evaluate_bitstring(state, paulistring)
+        for paulistring, coeff in observable.label_iter()
+    )
+
+
+def _evaluate_bitstring(state: int, paulistring: str) -> float:
+    """Evaluate a bitstring on a Pauli label."""
+    n = len(paulistring) - 1
+    return np.prod(
+        [-1 if state & (1 << (n - i)) else 1 for i, pauli in enumerate(paulistring) if pauli == "Z"]
+    )
+
+
+def _check_observable_is_diagonal(observable: SparsePauliOp) -> bool:
+    is_diagonal = not np.any(observable.paulis.x)
+    if not is_diagonal:
+        raise ValueError("The observable must be diagonal.")

qiskit/algorithms/minimum_eigensolvers/qaoa.py

@@ -0,0 +1,149 @@
+# 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.
+
+"""The quantum approximate optimization algorithm. """
+
+from __future__ import annotations
+
+from typing import Callable, Any
+import numpy as np
+
+from qiskit.algorithms.optimizers import Minimizer, Optimizer
+from qiskit.circuit import QuantumCircuit
+from qiskit.circuit.library.n_local.qaoa_ansatz import QAOAAnsatz
+from qiskit.quantum_info.operators.base_operator import BaseOperator
+from qiskit.opflow import PauliSumOp, PrimitiveOp
+from qiskit.primitives import BaseSampler
+from qiskit.utils.validation import validate_min
+
+from ..exceptions import AlgorithmError
+from .sampling_vqe import SamplingVQE
+
+
+class QAOA(SamplingVQE):
+    r"""
+    The Quantum Approximate Optimization Algorithm (QAOA).
+
+    QAOA is a well-known algorithm for finding approximate solutions to combinatorial-optimization
+    problems [1].
+
+    The QAOA implementation directly extends :class:`.SamplingVQE` and inherits its optimization
+    structure. However, unlike VQE, which can be configured with arbitrary ansatzes, QAOA uses its
+    own fine-tuned ansatz, which comprises :math:`p` parameterized global :math:`x` rotations and
+    :math:`p` different parameterizations of the problem hamiltonian. QAOA is thus principally
+    configured  by the single integer parameter, ``reps``, which dictates the depth of the ansatz,
+    and thus affects the approximation quality.
+
+    An optional array of :math:`2p` parameter values, as the :attr:`initial_point`, may be provided
+    as the starting :math:`\beta` and :math:`\gamma` parameters for the QAOA ansatz [1].
+
+    An operator or a parameterized quantum circuit may optionally also be provided as a custom
+    :attr:`mixer` Hamiltonian. This allows in the case of quantum annealing [2] and QAOA [3], to run
+    constrained optimization problems where the mixer constrains the evolution to a feasible
+    subspace of the full Hilbert space.
+
+    The following attributes can be set via the initializer but can also be read and updated once
+    the QAOA object has been constructed.
+
+    Attributes:
+        sampler (BaseSampler): The sampler primitive to sample the circuits.
+        optimizer (Optimizer | Minimizer): A classical optimizer to find the minimum energy. This
+            can either be a Qiskit :class:`.Optimizer` or a callable implementing the
+            :class:`.Minimizer` protocol.
+        reps (int): The integer parameter :math:`p`. Has a minimum valid value of 1.
+        initial_state: An optional initial state to prepend the QAOA circuit with.
+        mixer (QuantumCircuit | BaseOperator | PauliSumOp): The mixer Hamiltonian to evolve with or
+            a custom quantum circuit. Allows support of optimizations in constrained subspaces [2,
+            3] as well as warm-starting the optimization [4].
+        aggregation (float | Callable[[list[float]], float] | None): A float or callable to specify
+            how the objective function evaluated on the basis states should be aggregated. If a
+            float, this specifies the :math:`\alpha \in [0,1]` parameter for a CVaR expectation
+            value.
+        callback (Callable[[int, np.ndarray, float, dict[str, Any]], None] | None): A callback
+            that can access the intermediate data at each optimization step. These data are: the
+            evaluation count, the optimizer parameters for the ansatz, the evaluated value, the
+            the metadata dictionary, and the best measurement.
+
+    References:
+        [1]: Farhi, E., Goldstone, J., Gutmann, S., "A Quantum Approximate Optimization Algorithm"
+            `arXiv:1411.4028 <https://arxiv.org/abs/1411.4028>`__
+        [2]: Hen, I., Spedalieri, F. M.,  "Quantum Annealing for Constrained Optimization"
+            `PhysRevApplied.5.034007 <https://doi.org/10.1103/PhysRevApplied.5.034007>`__
+        [3]: Hadfield, S. et al, "From the Quantum Approximate Optimization Algorithm to a Quantum
+            Alternating Operator Ansatz" `arXiv:1709.03489 <https://arxiv.org/abs/1709.03489>`__
+        [4]: Egger, D. J., Marecek, J., Woerner, S., "Warm-starting quantum optimization"
+            `arXiv: 2009.10095 <https://arxiv.org/abs/2009.10095>`__
+    """
+
+    def __init__(
+        self,
+        sampler: BaseSampler,
+        optimizer: Optimizer | Minimizer,
+        *,
+        reps: int = 1,
+        initial_state: QuantumCircuit | None = None,
+        mixer: QuantumCircuit | BaseOperator | PauliSumOp = None,
+        initial_point: np.ndarray | None = None,
+        aggregation: float | Callable[[list[float]], float] | None = None,
+        callback: Callable[[int, np.ndarray, float, dict[str, Any]], None] | None = None,
+    ) -> None:
+        r"""
+        Args:
+            sampler: The sampler primitive to sample the circuits.
+            optimizer: A classical optimizer to find the minimum energy. This can either be a
+                Qiskit :class:`.Optimizer` or a callable implementing the :class:`.Minimizer`
+                protocol.
+            reps: The integer parameter :math:`p`. Has a minimum valid value of 1.
+            initial_state: An optional initial state to prepend the QAOA circuit with.
+            mixer: The mixer Hamiltonian to evolve with or a custom quantum circuit. Allows support
+                of optimizations in constrained subspaces [2, 3] as well as warm-starting the
+                optimization [4].
+            initial_point: An optional initial point (i.e. initial parameter values) for the
+                optimizer. The length of the initial point must match the number of :attr:`ansatz`
+                parameters. If ``None``, a random point will be generated within certain parameter
+                bounds. ``QAOA`` will look to the ansatz for these bounds. If the ansatz does not
+                specify bounds, bounds of :math:`-2\pi`, :math:`2\pi` will be used.
+            aggregation: A float or callable to specify how the objective function evaluated on the
+                basis states should be aggregated. If a float, this specifies the :math:`\alpha \in
+                [0,1]` parameter for a CVaR expectation value.
+            callback: A callback that can access the intermediate data at each optimization step.
+                These data are: the evaluation count, the optimizer parameters for the ansatz, the
+                evaluated value, the the metadata dictionary.
+        """
+        validate_min("reps", reps, 1)
+
+        self.reps = reps
+        self.mixer = mixer
+        self.initial_state = initial_state
+        self._cost_operator = None
+
+        super().__init__(
+            sampler=sampler,
+            ansatz=None,
+            optimizer=optimizer,
+            initial_point=initial_point,
+            aggregation=aggregation,
+            callback=callback,
+        )
+
+    def _check_operator_ansatz(self, operator: BaseOperator | PauliSumOp):
+        if isinstance(operator, BaseOperator):
+            try:
+                operator = PrimitiveOp(operator)
+            except TypeError as error:
+                raise AlgorithmError(
+                    f"Unsupported operator type {type(operator)} passed to QAOA."
+                ) from error
+        # Recreates a circuit based on operator parameter.
+        self.ansatz = QAOAAnsatz(
+            operator, self.reps, initial_state=self.initial_state, mixer_operator=self.mixer
+        ).decompose()  # TODO remove decompose once #6674 is fixed