Use primitives for solving the relaxed problem (#52)

* WIP: use primitives for solving the relaxed problem

* Use QAOAAnsatz with VQE

The new QAOA implementation in qiskit-terra assumes the problem
Hamiltonian is diagonal so that it can benefit from certain
additional features, like applying CVaR or other aggregation
functions.  However, for us the problem Hamiltonian is not
diagonal.  So, we implement QAOA using the `QAOAAnsatz` together
with `VQE`.

* Use the default mixer in QAOA, rather than the problem Hamiltonian

* The optimizer now works with both the old and new eigen[ ]solvers

* Update optimizer to accept new eigensolvers

Parts of this follow qiskit-community/qiskit-optimization#436

* Restore tests for the other legacy eigen_solvers

(and suppress the `PendingDeprecationWarning`s using pytest)

* Make the code path for `LegacyNumPyMinimumEigensolver` more explicit

This way, it will be more clear when this branch/block can be
removed in the future

* Fix a typo

qrao/encoding.py

@@ -39,7 +39,6 @@
     X,
     Y,
     Z,
-    PauliOp,
     PauliSumOp,
     PrimitiveOp,
     CircuitOp,
@@ -48,6 +47,7 @@
     StateFn,
     CircuitStateFn,
 )
+from qiskit.quantum_info import SparsePauliOp
 
 from qiskit_optimization.problems.quadratic_program import QuadraticProgram
 
@@ -278,7 +278,7 @@ def __init__(self, max_vars_per_qubit: int = 3):
             raise ValueError("max_vars_per_qubit must be 1, 2, or 3")
         self._ops = self.OPERATORS[max_vars_per_qubit - 1]
 
-        self._qubit_op: Optional[Union[PauliOp, PauliSumOp]] = None
+        self._qubit_op: Optional[PauliSumOp] = None
         self._offset: Optional[float] = None
         self._problem: Optional[QuadraticProgram] = None
         self._var2op: Dict[int, Tuple[int, PrimitiveOp]] = {}
@@ -329,7 +329,7 @@ def minimum_recovery_probability(self) -> float:
         return (1 + 1 / np.sqrt(n)) / 2
 
     @property
-    def qubit_op(self) -> Union[PauliOp, PauliSumOp]:
+    def qubit_op(self) -> PauliSumOp:
         """Relaxed Hamiltonian operator"""
         if self._qubit_op is None:
             raise AttributeError(
@@ -400,7 +400,7 @@ def _add_term(self, w: float, *variables: int) -> None:
         # Pauli operator.  To generalize, we replace the `d` in that equation
         # with `d_prime`, defined as follows:
         d_prime = np.sqrt(self.max_vars_per_qubit) ** len(variables)
-        op = w * d_prime * self.term2op(*variables)
+        op = self.term2op(*variables).mul(w * d_prime)
         # We perform the following short-circuit *after* calling term2op so at
         # least we have confirmed that the user provided a valid variables list.
         if w == 0.0:
@@ -410,8 +410,8 @@ def _add_term(self, w: float, *variables: int) -> None:
         else:
             self._qubit_op += op
 
-    def term2op(self, *variables: int) -> PauliOp:
-        """Construct a ``PauliOp`` that is a product of encoded decision ``variable``\\(s).
+    def term2op(self, *variables: int) -> PauliSumOp:
+        """Construct a ``PauliSumOp`` that is a product of encoded decision ``variable``\\(s).
 
         The decision variables provided must all be encoded on different qubits.
         """
@@ -423,7 +423,9 @@ def term2op(self, *variables: int) -> PauliOp:
                 raise RuntimeError(f"Collision of variables: {variables}")
             ops[pos] = op
             done.add(pos)
-        return reduce(lambda x, y: x ^ y, ops)
+        pauli_op = reduce(lambda x, y: x ^ y, ops)
+        # Convert from PauliOp to PauliSumOp
+        return PauliSumOp(SparsePauliOp(pauli_op.primitive, coeffs=[pauli_op.coeff]))
 
     @staticmethod
     def _generate_ising_terms(

qrao/quantum_random_access_optimizer.py

@@ -18,11 +18,16 @@
 import numpy as np
 
 from qiskit import QuantumCircuit
-from qiskit.algorithms import (
+from qiskit.algorithms.minimum_eigensolvers import (
     MinimumEigensolver,
     MinimumEigensolverResult,
     NumPyMinimumEigensolver,
 )
+from qiskit.algorithms.minimum_eigen_solvers import (
+    MinimumEigensolver as LegacyMinimumEigensolver,
+    MinimumEigensolverResult as LegacyMinimumEigensolverResult,
+    NumPyMinimumEigensolver as LegacyNumPyMinimumEigensolver,
+)
 
 from qiskit_optimization.algorithms import (
     OptimizationAlgorithm,
@@ -37,6 +42,16 @@
 from .semideterministic_rounding import SemideterministicRounding
 
 
+def _get_aux_operators_evaluated(relaxed_results):
+    try:
+        # Must be using the new "minimum_eigensolvers"
+        # https://github.com/Qiskit/qiskit-terra/blob/main/releasenotes/notes/0.22/add-eigensolvers-with-primitives-8b3a9f55f5fd285f.yaml
+        return relaxed_results.aux_operators_evaluated
+    except AttributeError:
+        # Must be using the old (deprecated) "minimum_eigen_solvers"
+        return relaxed_results.aux_operator_eigenvalues
+
+
 class QuantumRandomAccessOptimizationResult(OptimizationResult):
     """Result of Quantum Random Access Optimization procedure."""
 
@@ -48,7 +63,9 @@ def __init__(
         variables: List[Variable],
         status: OptimizationResultStatus,
         samples: Optional[List[SolutionSample]],
-        relaxed_results: MinimumEigensolverResult,
+        relaxed_results: Union[
+            MinimumEigensolverResult, LegacyMinimumEigensolverResult
+        ],
         rounding_results: RoundingResult,
         relaxed_results_offset: float,
         sense: int,
@@ -78,7 +95,9 @@ def __init__(
         self._sense = sense
 
     @property
-    def relaxed_results(self) -> MinimumEigensolverResult:
+    def relaxed_results(
+        self,
+    ) -> Union[MinimumEigensolverResult, LegacyMinimumEigensolverResult]:
         """Variationally obtained ground state of the relaxed Hamiltonian"""
         return self._relaxed_results
 
@@ -90,7 +109,9 @@ def rounding_results(self) -> RoundingResult:
     @property
     def trace_values(self):
         """List of expectation values, one corresponding to each decision variable"""
-        trace_values = [v[0] for v in self._relaxed_results.aux_operator_eigenvalues]
+        trace_values = [
+            v[0] for v in _get_aux_operators_evaluated(self._relaxed_results)
+        ]
         return trace_values
 
     @property
@@ -120,7 +141,7 @@ class QuantumRandomAccessOptimizer(OptimizationAlgorithm):
 
     def __init__(
         self,
-        min_eigen_solver: MinimumEigensolver,
+        min_eigen_solver: Union[MinimumEigensolver, LegacyMinimumEigensolver],
         encoding: QuantumRandomAccessEncoding,
         rounding_scheme: Optional[RoundingScheme] = None,
     ):
@@ -144,12 +165,14 @@ def __init__(
         self.rounding_scheme = rounding_scheme
 
     @property
-    def min_eigen_solver(self) -> MinimumEigensolver:
+    def min_eigen_solver(self) -> Union[MinimumEigensolver, LegacyMinimumEigensolver]:
         """The minimum eigensolver."""
         return self._min_eigen_solver
 
     @min_eigen_solver.setter
-    def min_eigen_solver(self, min_eigen_solver: MinimumEigensolver) -> None:
+    def min_eigen_solver(
+        self, min_eigen_solver: Union[MinimumEigensolver, LegacyMinimumEigensolver]
+    ) -> None:
         """Set the minimum eigensolver."""
         if not min_eigen_solver.supports_aux_operators():
             raise TypeError(
@@ -184,7 +207,11 @@ def get_compatibility_msg(self, problem: QuadraticProgram) -> str:
             )
         return ""
 
-    def solve_relaxed(self) -> Tuple[MinimumEigensolverResult, RoundingContext]:
+    def solve_relaxed(
+        self,
+    ) -> Tuple[
+        Union[MinimumEigensolverResult, LegacyMinimumEigensolverResult], RoundingContext
+    ]:
         """Solve the relaxed Hamiltonian given the ``encoding`` provided to the constructor."""
         # Get the ordered list of operators that correspond to each decision
         # variable.  This line assumes the variables are numbered consecutively
@@ -202,7 +229,7 @@ def solve_relaxed(self) -> Tuple[MinimumEigensolverResult, RoundingContext]:
         stop_time_relaxed = time.time()
         relaxed_results.time_taken = stop_time_relaxed - start_time_relaxed
 
-        trace_values = [v[0] for v in relaxed_results.aux_operator_eigenvalues]
+        trace_values = [v[0] for v in _get_aux_operators_evaluated(relaxed_results)]
 
         # Collect inputs for rounding
         # double check later that there's no funny business with the
@@ -218,10 +245,17 @@ def solve_relaxed(self) -> Tuple[MinimumEigensolverResult, RoundingContext]:
             circuit = self.min_eigen_solver.ansatz.bind_parameters(
                 relaxed_results.optimal_point
             )
-        elif isinstance(self.min_eigen_solver, NumPyMinimumEigensolver):
+        elif isinstance(
+            self.min_eigen_solver,
+            (NumPyMinimumEigensolver, LegacyNumPyMinimumEigensolver),
+        ):
             statevector = relaxed_results.eigenstate
+            if isinstance(self.min_eigen_solver, LegacyNumPyMinimumEigensolver):
+                # statevector is a StateFn in this case, so we must convert it
+                # to a Statevector
+                statevector = statevector.primitive
             circuit = QuantumCircuit(self.encoding.num_qubits)
-            circuit.initialize(statevector.primitive)
+            circuit.initialize(statevector)
         else:
             circuit = None
 

requirements.txt

@@ -1,6 +1,6 @@
 certifi>=2021.5.30
 importlib_metadata>=4.8.1,!=5.0.0
-qiskit-terra>=0.21.0
+qiskit-terra>=0.22.0
 qiskit-aer>=0.11.0
 qiskit-ibm-provider>=0.1.0
 qiskit-optimization>=0.4.0

tests/test_quantum_random_access_optimizer.py

@@ -13,16 +13,23 @@
 """Tests for QuantumRandomAccessOptimizer."""
 
 from unittest import TestCase
+import pytest
 
 from qiskit.utils import QuantumInstance
-from qiskit.algorithms.minimum_eigen_solvers import (
+from qiskit.algorithms.minimum_eigensolvers import (
     VQE,
-    QAOA,
     NumPyMinimumEigensolver,
     MinimumEigensolverResult,
 )
-from qiskit.circuit.library import RealAmplitudes
+from qiskit.algorithms.minimum_eigen_solvers import (
+    VQE as LegacyVQE,
+    QAOA as LegacyQAOA,
+    NumPyMinimumEigensolver as LegacyNumPyMinimumEigensolver,
+    MinimumEigensolverResult as LegacyMinimumEigensolverResult,
+)
+from qiskit.circuit.library import RealAmplitudes, QAOAAnsatz
 from qiskit.algorithms.optimizers import SPSA
+from qiskit.primitives import Estimator
 from qiskit_aer import Aer
 
 from qiskit_optimization.translators import from_docplex_mp
@@ -50,11 +57,28 @@ def setUp(self):
         self.ansatz = RealAmplitudes(self.encoding.num_qubits)  # for VQE
 
     def test_solve_relaxed_vqe(self):
-        """Test QuantumRandomAccessOptimizer."""
+        """Test QuantumRandomAccessOptimizer with VQE."""
+        estimator = Estimator(options={"shots": 100})
+        vqe = VQE(
+            ansatz=self.ansatz,
+            optimizer=SPSA(maxiter=1, learning_rate=0.01, perturbation=0.1),
+            estimator=estimator,
+        )
+
+        qrao = QuantumRandomAccessOptimizer(
+            encoding=self.encoding, min_eigen_solver=vqe
+        )
+        relaxed_results, rounding_context = qrao.solve_relaxed()
+        self.assertIsInstance(relaxed_results, MinimumEigensolverResult)
+        self.assertIsInstance(rounding_context, RoundingContext)
+
+    @pytest.mark.filterwarnings("ignore::PendingDeprecationWarning")
+    def test_solve_relaxed_legacy_vqe(self):
+        """Test QuantumRandomAccessOptimizer with legacy VQE."""
         relaxed_qi = QuantumInstance(
             backend=Aer.get_backend("aer_simulator"), shots=100
         )
-        vqe = VQE(
+        vqe = LegacyVQE(
             ansatz=self.ansatz,
             optimizer=SPSA(maxiter=1, learning_rate=0.01, perturbation=0.1),
             quantum_instance=relaxed_qi,
@@ -64,15 +88,34 @@ def test_solve_relaxed_vqe(self):
             encoding=self.encoding, min_eigen_solver=vqe
         )
         relaxed_results, rounding_context = qrao.solve_relaxed()
-        self.assertIsInstance(relaxed_results, MinimumEigensolverResult)
+        self.assertIsInstance(relaxed_results, LegacyMinimumEigensolverResult)
         self.assertIsInstance(rounding_context, RoundingContext)
 
     def test_solve_relaxed_qaoa(self):
-        """Test QuantumRandomAccessOptimizer."""
+        """Test QuantumRandomAccessOptimizer with QAOA."""
+        estimator = Estimator(options={"shots": 100})
+        qaoa_ansatz = QAOAAnsatz(
+            cost_operator=self.encoding.qubit_op,
+        )
+        qaoa = VQE(
+            ansatz=qaoa_ansatz,
+            optimizer=SPSA(maxiter=1, learning_rate=0.01, perturbation=0.1),
+            estimator=estimator,
+        )
+        qrao = QuantumRandomAccessOptimizer(
+            encoding=self.encoding, min_eigen_solver=qaoa
+        )
+        relaxed_results, rounding_context = qrao.solve_relaxed()
+        self.assertIsInstance(relaxed_results, MinimumEigensolverResult)
+        self.assertIsInstance(rounding_context, RoundingContext)
+
+    @pytest.mark.filterwarnings("ignore::PendingDeprecationWarning")
+    def test_solve_relaxed_legacy_qaoa(self):
+        """Test QuantumRandomAccessOptimizer with legacy QAOA."""
         relaxed_qi = QuantumInstance(
             backend=Aer.get_backend("aer_simulator"), shots=100
         )
-        qaoa = QAOA(
+        qaoa = LegacyQAOA(
             optimizer=SPSA(maxiter=1, learning_rate=0.01, perturbation=0.1),
             quantum_instance=relaxed_qi,
             mixer=self.encoding.qubit_op,
@@ -81,11 +124,11 @@ def test_solve_relaxed_qaoa(self):
             encoding=self.encoding, min_eigen_solver=qaoa
         )
         relaxed_results, rounding_context = qrao.solve_relaxed()
-        self.assertIsInstance(relaxed_results, MinimumEigensolverResult)
+        self.assertIsInstance(relaxed_results, LegacyMinimumEigensolverResult)
         self.assertIsInstance(rounding_context, RoundingContext)
 
     def test_solve_relaxed_numpy(self):
-        """Test QuantumRandomAccessOptimizer."""
+        """Test QuantumRandomAccessOptimizer with NumPyMinimumEigensolver."""
         np_solver = NumPyMinimumEigensolver()
         qrao = QuantumRandomAccessOptimizer(
             encoding=self.encoding, min_eigen_solver=np_solver
@@ -94,6 +137,17 @@ def test_solve_relaxed_numpy(self):
         self.assertIsInstance(relaxed_results, MinimumEigensolverResult)
         self.assertIsInstance(rounding_context, RoundingContext)
 
+    @pytest.mark.filterwarnings("ignore::PendingDeprecationWarning")
+    def test_solve_relaxed_legacy_numpy(self):
+        """Test QuantumRandomAccessOptimizer with legacy NumPyMinimumEigensolver."""
+        np_solver = LegacyNumPyMinimumEigensolver()
+        qrao = QuantumRandomAccessOptimizer(
+            encoding=self.encoding, min_eigen_solver=np_solver
+        )
+        relaxed_results, rounding_context = qrao.solve_relaxed()
+        self.assertIsInstance(relaxed_results, LegacyMinimumEigensolverResult)
+        self.assertIsInstance(rounding_context, RoundingContext)
+
     def test_different_problem(self):
         """Test passing a different problem to solve() than the encoding."""
         np_solver = NumPyMinimumEigensolver()
@@ -123,30 +177,26 @@ class ModifiedVQE(VQE):
             def supports_aux_operators(cls) -> bool:
                 return False
 
-        relaxed_qi = QuantumInstance(
-            backend=Aer.get_backend("aer_simulator"), shots=100
-        )
+        estimator = Estimator(options={"shots": 100})
         vqe = ModifiedVQE(
             ansatz=self.ansatz,
             optimizer=SPSA(maxiter=1, learning_rate=0.01, perturbation=0.1),
-            quantum_instance=relaxed_qi,
+            estimator=estimator,
         )
         with self.assertRaises(TypeError):
             QuantumRandomAccessOptimizer(encoding=self.encoding, min_eigen_solver=vqe)
 
     def test_solve_without_args(self):
         """Test solve() without arguments"""
-        relaxed_qi = QuantumInstance(
-            backend=Aer.get_backend("aer_simulator_statevector"), shots=100
-        )
+        estimator = Estimator(options={"shots": 100})
         rounding_qi = QuantumInstance(
             backend=Aer.get_backend("aer_simulator"), shots=100
         )
 
         vqe = VQE(
             ansatz=self.ansatz,
             optimizer=SPSA(maxiter=1, learning_rate=0.01, perturbation=0.1),
-            quantum_instance=relaxed_qi,
+            estimator=estimator,
         )
         rounding_scheme = MagicRounding(quantum_instance=rounding_qi)
 
@@ -164,13 +214,11 @@ def test_solve_without_args(self):
 
     def test_empty_encoding(self):
         """Test that an exception is raised if the encoding has no qubits"""
-        relaxed_qi = QuantumInstance(
-            backend=Aer.get_backend("aer_simulator_statevector"), shots=100
-        )
+        estimator = Estimator(options={"shots": 100})
         vqe = VQE(
             ansatz=self.ansatz,
             optimizer=SPSA(maxiter=1, learning_rate=0.01, perturbation=0.1),
-            quantum_instance=relaxed_qi,
+            estimator=estimator,
         )
         encoding = QuantumRandomAccessEncoding(3)
         with self.assertRaises(ValueError):