Allow list of optimizers+initial points in VQD (Qiskit#9151)

* test commit

* added optional lists of init
ial point and optimizer

* fixed initial points and add more compact expressions

* fixed bug initial_point input

* fix initial points

* formatting

* fix none initial point

* Added releasenote

* Ensure initial_point sampled once only if None

* moved release note

* Update typehints

* Update typehints

* Update typehints

* Update typehints

* Update typehints

* Apply suggestions from code review

* Fix conflict, add test

* Fix black

* Update releasenotes/notes/vqd-list-initial-points-list-optimizers-033d7439f86bbb71.yaml

---------

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

qiskit/algorithms/eigensolvers/vqd.py

@@ -18,9 +18,9 @@
 from __future__ import annotations
 
 from collections.abc import Callable, Sequence
+from typing import Any
 import logging
 from time import time
-from typing import Any
 
 import numpy as np
 
@@ -61,7 +61,7 @@ class VQD(VariationalAlgorithm, Eigensolver):
     An instance of VQD requires defining three algorithmic sub-components:
     an integer k denoting the number of eigenstates to calculate, a trial
     state (a.k.a. ansatz) which is a :class:`QuantumCircuit`,
-    and one of the classical :mod:`~qiskit.algorithms.optimizers`.
+    and one instance (or list of) classical :mod:`~qiskit.algorithms.optimizers`.
     The optimizer varies the circuit parameters
     The trial state :math:`|\psi(\vec\theta)\rangle` is varied by the optimizer,
     which modifies the set of ansatz parameters :math:`\vec\theta`
@@ -79,6 +79,7 @@ class VQD(VariationalAlgorithm, Eigensolver):
     of ``None``, then VQD will look to the ansatz for a preferred value, based on its
     given initial state. If the ansatz returns ``None``,
     then a random point will be generated within the parameter bounds set, as per above.
+    It is also possible to give a list of initial points, one for every kth eigenvalue.
     If the ansatz provides ``None`` as the lower bound, then VQD
     will default it to :math:`-2\pi`; similarly, if the ansatz returns ``None``
     as the upper bound, the default value will be :math:`2\pi`.
@@ -92,18 +93,21 @@ class VQD(VariationalAlgorithm, Eigensolver):
             fidelity (BaseStateFidelity): The fidelity class instance used to compute the
                 overlap estimation as indicated in the VQD paper.
             ansatz (QuantumCircuit): A parameterized circuit used as ansatz for the wave function.
-            optimizer(Optimizer): A classical optimizer. Can either be a Qiskit optimizer or a callable
+            optimizer(Optimizer | Sequence[Optimizer]): A classical optimizer or a list of optimizers,
+            one for every k-th eigenvalue. Can either be a Qiskit optimizer or a callable
                 that takes an array as input and returns a Qiskit or SciPy optimization result.
             k (int): the number of eigenvalues to return. Returns the lowest k eigenvalues.
             betas (list[float]): Beta parameters in the VQD paper.
                 Should have length k - 1, with k the number of excited states.
                 These hyper-parameters balance the contribution of each overlap term to the cost
                 function and have a default value computed as the mean square sum of the
                 coefficients of the observable.
-            initial point (list[float]): An optional initial point (i.e. initial parameter values)
-                for the optimizer. If ``None`` then VQD will look to the ansatz for a preferred
-                point and if not will simply compute a random one.
-            callback (Callable[[int, np.ndarray, float, dict[str, Any], int], None] | None):
+            initial point (Sequence[float] | Sequence[Sequence[float]] | None): An optional initial
+                point (i.e. initial parameter values) or a list of initial points
+                (one for every k-th eigenvalue) for the optimizer.
+                If ``None`` then VQD will look to the ansatz for a
+                preferred point and if not will simply compute a random one.
+            callback (Callable[[int, np.ndarray, float, dict[str, Any]], None] | None):
                 A callback that can access the intermediate data
                 during the optimization. Four parameter values are passed to the callback as
                 follows during each evaluation by the optimizer: the evaluation count,
@@ -116,20 +120,21 @@ def __init__(
         estimator: BaseEstimator,
         fidelity: BaseStateFidelity,
         ansatz: QuantumCircuit,
-        optimizer: Optimizer | Minimizer,
+        optimizer: Optimizer | Minimizer | Sequence[Optimizer | Minimizer],
         *,
         k: int = 2,
         betas: Sequence[float] | None = None,
-        initial_point: Sequence[float] | None = None,
-        callback: Callable[[int, np.ndarray, float, dict[str, Any], int], None] | None = None,
+        initial_point: Sequence[float] | Sequence[Sequence[float]] | None = None,
+        callback: Callable[[int, np.ndarray, float, dict[str, Any]], None] | None = None,
     ) -> None:
         """
 
         Args:
             estimator: The estimator primitive.
             fidelity: The fidelity class using primitives.
             ansatz: A parameterized circuit used as ansatz for the wave function.
-            optimizer: A classical optimizer. Can either be a Qiskit optimizer or a callable
+            optimizer: A classical optimizer or a list of optimizers, one for every k-th eigenvalue.
+            Can either be a Qiskit optimizer or a callable
                 that takes an array as input and returns a Qiskit or SciPy optimization result.
             k: The number of eigenvalues to return. Returns the lowest k eigenvalues.
             betas: Beta parameters in the VQD paper.
@@ -138,7 +143,9 @@ def __init__(
                 function and have a default value computed as the mean square sum of the
                 coefficients of the observable.
             initial_point: An optional initial point (i.e. initial parameter values)
-                for the optimizer. If ``None`` then VQD will look to the ansatz for a preferred
+                or a list of initial points (one for every k-th eigenvalue)
+                for the optimizer.
+                If ``None`` then VQD will look to the ansatz for a preferred
                 point and if not will simply compute a random one.
             callback: A callback that can access the intermediate data
                 during the optimization. Four parameter values are passed to the callback as
@@ -161,12 +168,12 @@ def __init__(
         self._eval_count = 0
 
     @property
-    def initial_point(self) -> Sequence[float] | None:
+    def initial_point(self) -> Sequence[float] | Sequence[Sequence[float]] | None:
         """Returns initial point."""
         return self._initial_point
 
     @initial_point.setter
-    def initial_point(self, initial_point: Sequence[float]):
+    def initial_point(self, initial_point: Sequence[float] | Sequence[Sequence[float]] | None):
         """Sets initial point"""
         self._initial_point = initial_point
 
@@ -199,8 +206,6 @@ def compute_eigenvalues(
         # validation
         self._check_operator_ansatz(operator)
 
-        initial_point = validate_initial_point(self.initial_point, self.ansatz)
-
         bounds = validate_bounds(self.ansatz)
 
         # We need to handle the array entries being zero or Optional i.e. having value None
@@ -251,8 +256,20 @@ def compute_eigenvalues(
         # the same parameters to the ansatz if we do multiple steps
         prev_states = []
 
+        num_initial_points = 0
+        if self.initial_point is not None:
+            initial_points = np.reshape(self.initial_point, (-1, self.ansatz.num_parameters))
+            num_initial_points = len(initial_points)
+
+        # 0 just means the initial point is ``None`` and ``validate_initial_point``
+        # will select a random point
+        if num_initial_points <= 1:
+            initial_point = validate_initial_point(self.initial_point, self.ansatz)
+
         for step in range(1, self.k + 1):
-            # update list of optimal circuits
+            if num_initial_points > 1:
+                initial_point = validate_initial_point(initial_points[step - 1], self.ansatz)
+
             if step > 1:
                 prev_states.append(self.ansatz.bind_parameters(result.optimal_points[-1]))
 
@@ -264,20 +281,27 @@ def compute_eigenvalues(
             start_time = time()
 
             # TODO: add gradient support after FidelityGradients are implemented
-            if callable(self.optimizer):
-                opt_result = self.optimizer(fun=energy_evaluation, x0=initial_point, bounds=bounds)
+            if isinstance(self.optimizer, Sequence):
+                optimizer = self.optimizer[step - 1]
+            else:
+                optimizer = self.optimizer  # fall back to single optimizer if not list
+
+            if callable(optimizer):
+                opt_result = optimizer(  # pylint: disable=not-callable
+                    fun=energy_evaluation, x0=initial_point, bounds=bounds
+                )
             else:
                 # we always want to submit as many estimations per job as possible for minimal
                 # overhead on the hardware
-                was_updated = _set_default_batchsize(self.optimizer)
+                was_updated = _set_default_batchsize(optimizer)
 
-                opt_result = self.optimizer.minimize(
+                opt_result = optimizer.minimize(
                     fun=energy_evaluation, x0=initial_point, bounds=bounds
                 )
 
                 # reset to original value
                 if was_updated:
-                    self.optimizer.set_max_evals_grouped(None)
+                    optimizer.set_max_evals_grouped(None)
 
             eval_time = time() - start_time
 

releasenotes/notes/vqd-list-initial-points-list-optimizers-033d7439f86bbb71.yaml

@@ -0,0 +1,10 @@
+---
+features:
+  - |
+    Added extensions to the :class:`~.eigensolvers.VQD` algorithm, which allow
+    to pass a list of optimizers and initial points for the different
+    minimization runs. For example, the ``k``-th initial point and 
+    ``k``-th optimizer will be used for the optimization of the
+    ``k-1``-th exicted state.
+
+

test/python/algorithms/eigensolvers/test_vqd.py

@@ -239,6 +239,46 @@ def run_check():
             result = vqd.compute_eigenvalues(operator=op)
             self.assertIsInstance(result, VQDResult)
 
+    @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
+    def test_optimizer_list(self, op):
+        """Test sending an optimizer list"""
+
+        optimizers = [SLSQP(), L_BFGS_B()]
+        initial_point_1 = [
+            1.70256666,
+            -5.34843975,
+            -0.39542903,
+            5.99477786,
+            -2.74374986,
+            -4.85284669,
+            0.2442925,
+            -1.51638917,
+        ]
+        initial_point_2 = [
+            0.5,
+            0.5,
+            0.5,
+            0.5,
+            0.5,
+            0.5,
+            0.5,
+            0.5,
+        ]
+        vqd = VQD(
+            estimator=self.estimator,
+            fidelity=self.fidelity,
+            ansatz=RealAmplitudes(),
+            optimizer=optimizers,
+            initial_point=[initial_point_1, initial_point_2],
+            k=2,
+            betas=self.betas,
+        )
+
+        result = vqd.compute_eigenvalues(operator=op)
+        np.testing.assert_array_almost_equal(
+            result.eigenvalues.real, self.h2_energy_excited[:2], decimal=3
+        )
+
     @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_aux_operators_list(self, op):
         """Test list-based aux_operators."""