Fix VQD with SPSA optimizer (Qiskit#9538)

* Add parameter batching

* Add spsa to tests

* Add reno

* Remove batch size from evals

* Apply suggestions from code review

* Apply suggestion from code review

* Update tests

* Delete Untitled.ipynb

* Fix batch size SPSA

* Small fixes

* fix batching

* update test

* Update reno

* Change SPSA test

* Apply reno suggestion

Co-authored-by: Julien Gacon <[email protected]>

---------

Co-authored-by: Julien Gacon <[email protected]>
Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com>

qiskit/algorithms/eigensolvers/vqd.py

@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2022.
+# (C) Copyright IBM 2022, 2023.
 #
 # 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
@@ -193,7 +193,6 @@ def compute_eigenvalues(
         operator: BaseOperator | PauliSumOp,
         aux_operators: ListOrDict[BaseOperator | PauliSumOp] | None = None,
     ) -> VQDResult:
-
         super().compute_eigenvalues(operator, aux_operators)
 
         # this sets the size of the ansatz, so it must be called before the initial point
@@ -226,7 +225,6 @@ def compute_eigenvalues(
             aux_operators = None
 
         if self.betas is None:
-
             if isinstance(operator, PauliSumOp):
                 operator = operator.coeff * operator.primitive
 
@@ -254,7 +252,6 @@ def compute_eigenvalues(
         prev_states = []
 
         for step in range(1, self.k + 1):
-
             # update list of optimal circuits
             if step > 1:
                 prev_states.append(self.ansatz.bind_parameters(result.optimal_points[-1]))
@@ -365,22 +362,28 @@ def _get_evaluate_energy(
         self._check_operator_ansatz(operator)
 
         def evaluate_energy(parameters: np.ndarray) -> np.ndarray | float:
+            # handle broadcasting: ensure parameters is of shape [array, array, ...]
+            if len(parameters.shape) == 1:
+                parameters = np.reshape(parameters, (-1, num_parameters))
+            batch_size = len(parameters)
 
             estimator_job = self.estimator.run(
-                circuits=[self.ansatz], observables=[operator], parameter_values=[parameters]
+                batch_size * [self.ansatz], batch_size * [operator], parameters
             )
 
-            total_cost = 0
+            total_cost = np.zeros(batch_size)
+
             if step > 1:
                 # compute overlap cost
+                batched_prev_states = [state for state in prev_states for _ in range(batch_size)]
                 fidelity_job = self.fidelity.run(
-                    [self.ansatz] * (step - 1),
-                    prev_states,
-                    [parameters] * (step - 1),
+                    batch_size * [self.ansatz] * (step - 1),
+                    batched_prev_states,
+                    np.tile(parameters, (step - 1, 1)),
                 )
-
                 costs = fidelity_job.result().fidelities
 
+                costs = np.reshape(costs, (step - 1, -1))
                 for state, cost in enumerate(costs):
                     total_cost += np.real(betas[state] * cost)
 
@@ -394,7 +397,7 @@ def evaluate_energy(parameters: np.ndarray) -> np.ndarray | float:
 
             if self.callback is not None:
                 metadata = estimator_result.metadata
-                for params, value, meta in zip([parameters], values, metadata):
+                for params, value, meta in zip(parameters, values, metadata):
                     self._eval_count += 1
                     self.callback(self._eval_count, params, value, meta, step)
             else:
@@ -406,7 +409,6 @@ def evaluate_energy(parameters: np.ndarray) -> np.ndarray | float:
 
     @staticmethod
     def _build_vqd_result() -> VQDResult:
-
         result = VQDResult()
         result.optimal_points = []
         result.optimal_parameters = []
@@ -420,7 +422,6 @@ def _build_vqd_result() -> VQDResult:
 
     @staticmethod
     def _update_vqd_result(result, opt_result, eval_time, ansatz) -> VQDResult:
-
         result.optimal_points.append(opt_result.x)
         result.optimal_parameters.append(dict(zip(ansatz.parameters, opt_result.x)))
         result.optimal_values.append(opt_result.fun)

test/python/algorithms/eigensolvers/test_vqd.py

@@ -19,13 +19,9 @@
 from ddt import data, ddt
 
 from qiskit import QuantumCircuit
-from qiskit.algorithms.eigensolvers import VQD
+from qiskit.algorithms.eigensolvers import VQD, VQDResult
 from qiskit.algorithms import AlgorithmError
-from qiskit.algorithms.optimizers import (
-    COBYLA,
-    L_BFGS_B,
-    SLSQP,
-)
+from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B, SLSQP, SPSA
 from qiskit.algorithms.state_fidelities import ComputeUncompute
 from qiskit.circuit.library import TwoLocal, RealAmplitudes
 from qiskit.opflow import PauliSumOp
@@ -207,6 +203,7 @@ def store_intermediate_result(eval_count, parameters, mean, metadata, step):
     @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_vqd_optimizer(self, op):
         """Test running same VQD twice to re-use optimizer, then switch optimizer"""
+
         vqd = VQD(
             estimator=self.estimator,
             fidelity=self.fidelity,
@@ -231,6 +228,17 @@ def run_check():
             vqd.optimizer = L_BFGS_B()
             run_check()
 
+        with self.subTest("Batched optimizer replace"):
+            vqd.optimizer = SLSQP(maxiter=60, max_evals_grouped=10)
+            run_check()
+
+        with self.subTest("SPSA replace"):
+            # SPSA takes too long to converge, so we will
+            # only check that it runs with no errors.
+            vqd.optimizer = SPSA(maxiter=5, learning_rate=0.01, perturbation=0.01)
+            result = vqd.compute_eigenvalues(operator=op)
+            self.assertIsInstance(result, VQDResult)
+
     @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_aux_operators_list(self, op):
         """Test list-based aux_operators."""