Update optimizers to allow new primitives-based algorithms (#436)

* let MinimumEigenOptimizer allow new primitive-based algorithms

* update other optimizers

* update other tests

* add sampler

* added unittests

* finalized the codes and unittests

* update docstring

* add reno

* update reno, requirements.txt, and tests

- Catch deprecation messages in tests

* Updates tests to suppress warnings

- catch PendingDeprecationWarnings of legacy tests
- replace a deprecated networkx function in test_sk_model
- replace a deprecated matplotlib function in bin_packing

* revert update of matplotlib and networkx

* add prelude

* Update releasenotes/notes/add-primitives-support-31af39549b5e66e3.yaml

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

* update reno

* revert grover reno file to avoid sphinx warning

* Update releasenotes/notes/add-primitives-support-31af39549b5e66e3.yaml

* update sample_most_likely

* update reno

* Apply Steve's suggestions from code review

Co-authored-by: a-matsuo <MATSUOA@jp.ibm.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>

Author: 4 people

.pylintdict

@@ -126,6 +126,7 @@ nosignatures
 np
 num
 numpy
+numpyminimumeigensolver
 october
 optimality
 optimizationresult

README.md

@@ -68,9 +68,9 @@ from docplex.mp.model import Model
 from qiskit_optimization.algorithms import MinimumEigenOptimizer
 from qiskit_optimization.translators import from_docplex_mp
 
-from qiskit.utils import algorithm_globals, QuantumInstance
-from qiskit import BasicAer
-from qiskit.algorithms import QAOA
+from qiskit.utils import algorithm_globals
+from qiskit.primitives import Sampler
+from qiskit.algorithms.minimum_eigensolvers import QAOA
 from qiskit.algorithms.optimizers import SPSA
 
 # Generate a graph of 4 nodes
@@ -97,9 +97,8 @@ seed = 1234
 algorithm_globals.random_seed = seed
 
 spsa = SPSA(maxiter=250)
-backend = BasicAer.get_backend('qasm_simulator')
-q_i = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)
-qaoa = QAOA(optimizer=spsa, reps=5, quantum_instance=q_i)
+sampler = Sampler()
+qaoa = QAOA(sampler=sampler, optimizer=spsa, reps=5)
 algorithm = MinimumEigenOptimizer(qaoa)
 result = algorithm.solve(problem)
 print(result.prettyprint())  # prints solution, x=[1, 0, 1, 0], the cost, fval=4

qiskit_optimization/__init__.py

@@ -26,10 +26,11 @@
 A uniform interface as well as automatic conversion between different problem representations
 allows users to solve problems using a large set of algorithms, from variational quantum algorithms,
 such as the Quantum Approximate Optimization Algorithm
-(:class:`~qiskit.algorithms.QAOA`), to
+(:class:`~qiskit.algorithms.minimum_eigensolver.QAOA`), to
 `Grover Adaptive Search <https://arxiv.org/abs/quant-ph/9607014>`_
 (:class:`~algorithms.GroverOptimizer`), leveraging
-fundamental :mod:`~qiskit.algorithms` provided by Qiskit Terra. Furthermore, the modular design
+fundamental :mod:`~qiskit.algorithms.minimum_eigensolver` provided by Qiskit Terra.
+Furthermore, the modular design
 of the optimization module allows it to be easily extended and facilitates rapid development and
 testing of new algorithms. Compatible classical optimizers are also provided for testing,
 validation, and benchmarking.

qiskit_optimization/algorithms/admm_optimizer.py

@@ -17,21 +17,21 @@
 from typing import List, Optional, Tuple, cast
 
 import numpy as np
-from qiskit.algorithms import NumPyMinimumEigensolver
+from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver
 
+from ..converters import MaximizeToMinimize
+from ..problems.constraint import Constraint
+from ..problems.linear_constraint import LinearConstraint
+from ..problems.linear_expression import LinearExpression
+from ..problems.quadratic_program import QuadraticProgram
+from ..problems.variable import Variable, VarType
 from .minimum_eigen_optimizer import MinimumEigenOptimizer
 from .optimization_algorithm import (
-    OptimizationResultStatus,
     OptimizationAlgorithm,
     OptimizationResult,
+    OptimizationResultStatus,
 )
 from .slsqp_optimizer import SlsqpOptimizer
-from ..problems.constraint import Constraint
-from ..problems.linear_constraint import LinearConstraint
-from ..problems.linear_expression import LinearExpression
-from ..problems.quadratic_program import QuadraticProgram
-from ..problems.variable import VarType, Variable
-from ..converters import MaximizeToMinimize
 
 UPDATE_RHO_BY_TEN_PERCENT = 0
 UPDATE_RHO_BY_RESIDUALS = 1

qiskit_optimization/algorithms/grover_optimizer.py

@@ -16,6 +16,7 @@
 import math
 from copy import deepcopy
 from typing import Optional, Dict, Union, List, cast
+import warnings
 
 import numpy as np
 
@@ -24,6 +25,7 @@
 from qiskit.utils import QuantumInstance, algorithm_globals
 from qiskit.algorithms.amplitude_amplifiers.grover import Grover
 from qiskit.circuit.library import QuadraticForm
+from qiskit.primitives import BaseSampler
 from qiskit.providers import Backend
 from qiskit.quantum_info import partial_trace
 from .optimization_algorithm import (
@@ -36,6 +38,7 @@
     QuadraticProgramToQubo,
     QuadraticProgramConverter,
 )
+from ..exceptions import QiskitOptimizationError
 from ..problems import Variable
 from ..problems.quadratic_program import QuadraticProgram
 
@@ -54,6 +57,7 @@ def __init__(
             Union[QuadraticProgramConverter, List[QuadraticProgramConverter]]
         ] = None,
         penalty: Optional[float] = None,
+        sampler: Optional[BaseSampler] = None,
     ) -> None:
         """
         Args:
@@ -66,20 +70,35 @@ def __init__(
                 :class:`~qiskit_optimization.converters.QuadraticProgramToQubo` will be used.
             penalty: The penalty factor used in the default
                 :class:`~qiskit_optimization.converters.QuadraticProgramToQubo` converter
+            sampler: A Sampler to use for sampling the results of the circuits.
 
         Raises:
+            ValueError: If both a quantum instance and sampler are set.
             TypeError: When there one of converters is an invalid type.
         """
         self._num_value_qubits = num_value_qubits
         self._num_key_qubits = 0
         self._n_iterations = num_iterations
-        self._quantum_instance = None  # type: Optional[QuantumInstance]
         self._circuit_results = {}  # type: dict
+        self._converters = self._prepare_converters(converters, penalty)
 
-        if quantum_instance is not None:
-            self.quantum_instance = quantum_instance
+        if quantum_instance is not None and sampler is not None:
+            raise ValueError("Only one of quantum_instance or sampler can be passed, not both!")
 
-        self._converters = self._prepare_converters(converters, penalty)
+        self._quantum_instance = None  # type: Optional[QuantumInstance]
+        if quantum_instance is not None:
+            warnings.warn(
+                "The quantum_instance argument has been superseded by the sampler argument. "
+                "This argument will be deprecated in a future release and subsequently "
+                "removed after that.",
+                category=PendingDeprecationWarning,
+                stacklevel=2,
+            )
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore", category=PendingDeprecationWarning)
+                self.quantum_instance = quantum_instance
+
+        self._sampler = sampler
 
     @property
     def quantum_instance(self) -> QuantumInstance:
@@ -88,6 +107,13 @@ def quantum_instance(self) -> QuantumInstance:
         Returns:
             The quantum instance used in the algorithm.
         """
+        warnings.warn(
+            "The quantum_instance argument has been superseded by the sampler argument. "
+            "This argument will be deprecated in a future release and subsequently "
+            "removed after that.",
+            category=PendingDeprecationWarning,
+            stacklevel=2,
+        )
         return self._quantum_instance
 
     @quantum_instance.setter
@@ -97,6 +123,13 @@ def quantum_instance(self, quantum_instance: Union[Backend, QuantumInstance]) ->
         Args:
             quantum_instance: The quantum instance to be used in the algorithm.
         """
+        warnings.warn(
+            "The GroverOptimizer.quantum_instance setter is pending deprecation. "
+            "This property will be deprecated in a future release and subsequently "
+            "removed after that.",
+            category=PendingDeprecationWarning,
+            stacklevel=2,
+        )
         if isinstance(quantum_instance, Backend):
             self._quantum_instance = QuantumInstance(quantum_instance)
         else:
@@ -162,11 +195,16 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
             The result of the optimizer applied to the problem.
 
         Raises:
+            ValueError: If a quantum instance or a sampler has not been provided.
+            ValueError: If both a quantum instance and sampler are set.
             AttributeError: If the quantum instance has not been set.
             QiskitOptimizationError: If the problem is incompatible with the optimizer.
         """
-        if self.quantum_instance is None:
-            raise AttributeError("The quantum instance or backend has not been set.")
+        if self._sampler is None and self._quantum_instance is None:
+            raise ValueError("A quantum instance or sampler must be provided.")
+
+        if self._quantum_instance is not None and self._sampler is not None:
+            raise ValueError("Only one of quantum_instance or sampler can be passed, not both!")
 
         self._verify_compatibility(problem)
 
@@ -199,7 +237,7 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
         # Initialize oracle helper object.
         qr_key_value = QuantumRegister(self._num_key_qubits + self._num_value_qubits)
         orig_constant = problem_.objective.constant
-        measurement = not self.quantum_instance.is_statevector
+        measurement = self._quantum_instance is None or not self._quantum_instance.is_statevector
         oracle, is_good_state = self._get_oracle(qr_key_value)
 
         while not optimum_found:
@@ -246,15 +284,19 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
                     threshold = optimum_value
 
                     # trace out work qubits and store samples
-                    if self._quantum_instance.is_statevector:
-                        indices = list(range(n_key, len(outcome)))
-                        rho = partial_trace(self._circuit_results, indices)
-                        self._circuit_results = cast(Dict, np.diag(rho.data) ** 0.5)
-                    else:
+                    if self._sampler is not None:
                         self._circuit_results = {
                             i[-1 * n_key :]: v for i, v in self._circuit_results.items()
                         }
-
+                    else:
+                        if self._quantum_instance.is_statevector:
+                            indices = list(range(n_key, len(outcome)))
+                            rho = partial_trace(self._circuit_results, indices)
+                            self._circuit_results = cast(Dict, np.diag(rho.data) ** 0.5)
+                        else:
+                            self._circuit_results = {
+                                i[-1 * n_key :]: v for i, v in self._circuit_results.items()
+                            }
                     raw_samples = self._eigenvector_to_solutions(
                         self._circuit_results, problem_init
                     )
@@ -312,33 +354,52 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
 
     def _measure(self, circuit: QuantumCircuit) -> str:
         """Get probabilities from the given backend, and picks a random outcome."""
-        probs = self._get_probs(circuit)
+        probs = self._get_prob_dist(circuit)
         logger.info("Frequencies: %s", probs)
         # Pick a random outcome.
         return algorithm_globals.random.choice(list(probs.keys()), 1, p=list(probs.values()))[0]
 
-    def _get_probs(self, qc: QuantumCircuit) -> Dict[str, float]:
+    def _get_prob_dist(self, qc: QuantumCircuit) -> Dict[str, float]:
         """Gets probabilities from a given backend."""
         # Execute job and filter results.
-        result = self.quantum_instance.execute(qc)
-        if self.quantum_instance.is_statevector:
-            state = result.get_statevector(qc)
-            if not isinstance(state, np.ndarray):
-                state = state.data
-            keys = [
-                bin(i)[2::].rjust(int(np.log2(len(state))), "0")[::-1] for i in range(0, len(state))
-            ]
-            probs = [abs(a) ** 2 for a in state]
-            total = math.fsum(probs)
-            probs = [p / total for p in probs]
-            hist = {key: prob for key, prob in zip(keys, probs) if prob > 0}
-            self._circuit_results = state
+        if self._sampler is not None:
+            job = self._sampler.run([qc])
+
+            try:
+                result = job.result()
+            except Exception as exc:
+                raise QiskitOptimizationError("Sampler job failed.") from exc
+            quasi_dist = result.quasi_dists[0]
+            bit_length = (len(quasi_dist) - 1).bit_length()
+            prob_dist = {f"{i:0{bit_length}b}"[::-1]: v for i, v in quasi_dist.items()}
+            self._circuit_results = {
+                f"{i:0{bit_length}b}": v**0.5
+                for i, v in quasi_dist.items()
+                if not np.isclose(v, 0)
+            }
         else:
-            state = result.get_counts(qc)
-            shots = self.quantum_instance.run_config.shots
-            hist = {key[::-1]: val / shots for key, val in sorted(state.items()) if val > 0}
-            self._circuit_results = {b: (v / shots) ** 0.5 for (b, v) in state.items()}
-        return hist
+            result = self._quantum_instance.execute(qc)
+            if self._quantum_instance.is_statevector:
+                state = result.get_statevector(qc)
+                if not isinstance(state, np.ndarray):
+                    state = state.data
+                keys = [
+                    bin(i)[2::].rjust(int(np.log2(len(state))), "0")[::-1]
+                    for i in range(0, len(state))
+                ]
+                probs = [abs(a) ** 2 for a in state]
+                total = math.fsum(probs)
+                probs = [p / total for p in probs]
+                prob_dist = {key: prob for key, prob in zip(keys, probs) if prob > 0}
+                self._circuit_results = state
+            else:
+                state = result.get_counts(qc)
+                shots = self._quantum_instance.run_config.shots
+                prob_dist = {
+                    key[::-1]: val / shots for key, val in sorted(state.items()) if val > 0
+                }
+                self._circuit_results = {b: (v / shots) ** 0.5 for (b, v) in state.items()}
+        return prob_dist
 
     @staticmethod
     def _bin_to_int(v: str, num_value_bits: int) -> int:

qiskit_optimization/algorithms/minimum_eigen_optimizer.py

@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2020, 2021.
+# (C) Copyright IBM 2020, 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
@@ -11,24 +11,37 @@
 # that they have been altered from the originals.
 
 """A wrapper for minimum eigen solvers to be used within the optimization module."""
-from typing import Optional, Union, List, cast
+from typing import List, Optional, Union, cast
 
 import numpy as np
+from qiskit.algorithms.minimum_eigen_solvers import MinimumEigensolver as LegacyMinimumEigensolver
+from qiskit.algorithms.minimum_eigen_solvers import (
+    MinimumEigensolverResult as LegacyMinimumEigensolverResult,
+)
+from qiskit.algorithms.minimum_eigensolvers import (
+    NumPyMinimumEigensolver,
+    NumPyMinimumEigensolverResult,
+    SamplingMinimumEigensolver,
+    SamplingMinimumEigensolverResult,
+)
+from qiskit.opflow import OperatorBase, PauliOp, PauliSumOp
 
-from qiskit.algorithms import MinimumEigensolver, MinimumEigensolverResult
-from qiskit.opflow import OperatorBase
+from ..converters.quadratic_program_to_qubo import QuadraticProgramConverter, QuadraticProgramToQubo
+from ..exceptions import QiskitOptimizationError
+from ..problems.quadratic_program import QuadraticProgram, Variable
 from .optimization_algorithm import (
-    OptimizationResultStatus,
     OptimizationAlgorithm,
     OptimizationResult,
+    OptimizationResultStatus,
     SolutionSample,
 )
-from ..exceptions import QiskitOptimizationError
-from ..converters.quadratic_program_to_qubo import (
-    QuadraticProgramToQubo,
-    QuadraticProgramConverter,
-)
-from ..problems.quadratic_program import QuadraticProgram, Variable
+
+MinimumEigensolver = Union[
+    SamplingMinimumEigensolver, NumPyMinimumEigensolver, LegacyMinimumEigensolver
+]
+MinimumEigensolverResult = Union[
+    SamplingMinimumEigensolverResult, NumPyMinimumEigensolverResult, LegacyMinimumEigensolverResult
+]
 
 
 class MinimumEigenOptimizationResult(OptimizationResult):
@@ -101,7 +114,7 @@ class MinimumEigenOptimizer(OptimizationAlgorithm):
 
     .. code-block::
 
-        from qiskit.algorithms import QAOA
+        from qiskit.algorithms.minimum_eigensolver import QAOA
         from qiskit_optimization.problems import QuadraticProgram
         from qiskit_optimization.algorithms import MinimumEigenOptimizer
         problem = QuadraticProgram()
@@ -141,7 +154,7 @@ def __init__(
         if not min_eigen_solver.supports_aux_operators():
             raise QiskitOptimizationError(
                 "Given MinimumEigensolver does not return the eigenstate "
-                + "and is not supported by the MinimumEigenOptimizer."
+                "and is not supported by the MinimumEigenOptimizer."
             )
         self._min_eigen_solver = min_eigen_solver
         self._penalty = penalty
@@ -206,6 +219,9 @@ def _solve_internal(
         # only try to solve non-empty Ising Hamiltonians
         eigen_result: Optional[MinimumEigensolverResult] = None
         if operator.num_qubits > 0:
+            # NumPyEigensolver does not accept PauliOp but PauliSumOp
+            if isinstance(operator, PauliOp):
+                operator = PauliSumOp.from_list([(operator.primitive.to_label(), operator.coeff)])
             # approximate ground state of operator using min eigen solver
             eigen_result = self._min_eigen_solver.compute_minimum_eigenvalue(operator)
             # analyze results