Adding a convergence threshold to VQD to filter non-sensible values

docs/tutorials/04_vqd.ipynb

@@ -15,7 +15,7 @@
    "source": [
     "## Introduction\n",
     "\n",
-    "VQD is a quantum algorithm that uses a variational technique to find the *k* eigenvalues of the Hamiltonian *H* of a given system.\n",
+    "VQD is a quantum algorithm that uses a variational technique to find the *k* lowest eigenvalues of the Hamiltonian *H* of a given system.\n",
     "\n",
     "The algorithm computes excited state energies of generalized hamiltonians by optimizing over a modified cost function. Each successive eigenvalue is calculated iteratively by introducing an overlap term with all the previously computed eigenstates that must be minimized. This ensures that higher energy eigenstates are found."
    ]
@@ -83,11 +83,11 @@
    ],
    "source": [
     "from qiskit.circuit.library import TwoLocal\n",
-    "from qiskit_algorithms.optimizers import SLSQP\n",
+    "from qiskit_algorithms.optimizers import COBYLA\n",
     "\n",
     "ansatz = TwoLocal(2, rotation_blocks=[\"ry\", \"rz\"], entanglement_blocks=\"cz\", reps=1)\n",
     "\n",
-    "optimizer = SLSQP()\n",
+    "optimizer = COBYLA()\n",
     "ansatz.decompose().draw(\"mpl\")"
    ]
   },
@@ -128,7 +128,7 @@
    "outputs": [],
    "source": [
     "k = 3\n",
-    "betas = [33, 33, 33]"
+    "betas = [3, 3, 3]"
    ]
   },
   {
@@ -360,7 +360,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.0"
+   "version": "3.12.6"
   }
  },
  "nbformat": 4,

qiskit_algorithms/eigensolvers/vqd.py

@@ -50,7 +50,7 @@ class VQD(VariationalAlgorithm, Eigensolver):
 
     `VQD <https://arxiv.org/abs/1805.08138>`__ is a quantum algorithm that uses a
     variational technique to find
-    the k eigenvalues of the Hamiltonian :math:`H` of a given system.
+    the k lowest eigenvalues of the Hamiltonian :math:`H` of a given system.
 
     The algorithm computes excited state energies of generalised hamiltonians
     by optimizing over a modified cost function where each successive eigenvalue
@@ -108,6 +108,9 @@ class VQD(VariationalAlgorithm, Eigensolver):
                 follows during each evaluation by the optimizer: the evaluation count,
                 the optimizer parameters for the ansatz, the estimated value, the estimation
                 metadata, and the current step.
+            convergence_threshold: A threshold under which the algorithm is considered to have
+                converged. It corresponds to the maximal average fidelity an eigenstate is allowed
+                to have with the previous eigenstates. If set to None, no check is performed.
     """
 
     def __init__(
@@ -121,6 +124,7 @@ def __init__(
         betas: np.ndarray | None = None,
         initial_point: np.ndarray | list[np.ndarray] | None = None,
         callback: Callable[[int, np.ndarray, float, dict[str, Any], int], None] | None = None,
+        convergence_threshold: float | None = None,
     ) -> None:
         """
 
@@ -147,6 +151,9 @@ def __init__(
                 follows during each evaluation by the optimizer: the evaluation count,
                 the optimizer parameters for the ansatz, the estimated value,
                 the estimation metadata, and the current step.
+            convergence_threshold: A threshold under which the algorithm is considered to have
+                converged. It corresponds to the maximal average fidelity an eigenstate is allowed
+                to have with the previous eigenstates. If set to None, no check is performed.
         """
         super().__init__()
 
@@ -159,6 +166,7 @@ def __init__(
         # this has to go via getters and setters due to the VariationalAlgorithm interface
         self.initial_point = initial_point
         self.callback = callback
+        self.convergence_threshold = convergence_threshold
 
         self._eval_count = 0
 
@@ -267,16 +275,24 @@ def compute_eigenvalues(
                 self.initial_point, self.ansatz  # type: ignore[arg-type]
             )
 
+        current_optimal_point: dict[str, Any] = {"optimal_value": float("inf")}
+
         for step in range(1, self.k + 1):
+            current_optimal_point["optimal_value"] = float("inf")
+
             if num_initial_points > 1:
                 initial_point = validate_initial_point(initial_points[step - 1], self.ansatz)
 
             if step > 1:
-                prev_states.append(self.ansatz.assign_parameters(result.optimal_points[-1]))
+                prev_states.append(self.ansatz.assign_parameters(current_optimal_point["x"]))
 
             self._eval_count = 0
             energy_evaluation = self._get_evaluate_energy(
-                step, operator, betas, prev_states=prev_states
+                step,
+                operator,
+                betas,
+                prev_states=prev_states,
+                current_optimal_point=current_optimal_point,
             )
 
             start_time = time()
@@ -309,11 +325,13 @@ def compute_eigenvalues(
 
             eval_time = time() - start_time
 
-            self._update_vqd_result(result, opt_result, eval_time, self.ansatz.copy())
+            self._update_vqd_result(
+                result, opt_result, eval_time, self.ansatz.copy(), current_optimal_point
+            )
 
             if aux_operators is not None:
                 aux_value = estimate_observables(
-                    self.estimator, self.ansatz, aux_operators, result.optimal_points[-1]
+                    self.estimator, self.ansatz, aux_operators, current_optimal_point["x"]
                 )
                 aux_values.append(aux_value)
 
@@ -326,6 +344,29 @@ def compute_eigenvalues(
                     self._eval_count,
                 )
             else:
+                average_fidelity = current_optimal_point["total_fidelity"][0] / (step - 1)
+
+                if (
+                    self.convergence_threshold is not None
+                    and average_fidelity > self.convergence_threshold
+                ):
+                    last_digit = step % 10
+
+                    if last_digit == 1 and step % 100 != 11:
+                        suffix = "st"
+                    elif last_digit == 2:
+                        suffix = "nd"
+                    elif last_digit == 3:
+                        suffix = "rd"
+                    else:
+                        suffix = "th"
+
+                    raise AlgorithmError(
+                        f"Convergence threshold is set to {self.convergence_threshold} but an "
+                        f"average fidelity of {average_fidelity:.5f} with the previous eigenstates"
+                        f"has been observed during the evaluation of the {step}{suffix} lowest"
+                        f"eigenvalue."
+                    )
                 logger.info(
                     (
                         "%s excited state optimization complete in %s s.\n"
@@ -345,21 +386,24 @@ def compute_eigenvalues(
 
         return result
 
-    def _get_evaluate_energy(
+    def _get_evaluate_energy(  # pylint: disable=too-many-positional-arguments
         self,
         step: int,
         operator: BaseOperator,
         betas: np.ndarray,
+        current_optimal_point: dict["str", Any],
         prev_states: list[QuantumCircuit] | None = None,
     ) -> Callable[[np.ndarray], float | np.ndarray]:
         """Returns a function handle to evaluate the ansatz's energy for any given parameters.
             This is the objective function to be passed to the optimizer that is used for evaluation.
 
         Args:
-            step: level of energy being calculated. 0 for ground, 1 for first excited state...
+            step: level of energy being calculated. 1 for ground, 2 for first excited state...
             operator: The operator whose energy to evaluate.
             betas: Beta parameters in the VQD paper.
             prev_states: List of optimal circuits from previous rounds of optimization.
+            current_optimal_point: A dict to keep track of the current optimal point, which is used
+                to check the algorithm's convergence.
 
         Returns:
             A callable that computes and returns the energy of the hamiltonian
@@ -425,6 +469,17 @@ def evaluate_energy(parameters: np.ndarray) -> float | np.ndarray:
             else:
                 self._eval_count += len(values)
 
+            for param, value in zip(parameters, values):
+                if value < current_optimal_point["optimal_value"]:
+                    current_optimal_point["optimal_value"] = value
+                    current_optimal_point["x"] = param
+
+                    if step > 1:
+                        current_optimal_point["total_fidelity"] = total_cost
+                        current_optimal_point["eigenvalue"] = (value - total_cost)[0]
+                    else:
+                        current_optimal_point["eigenvalue"] = value
+
             return values if len(values) > 1 else values[0]
 
         return evaluate_energy
@@ -444,20 +499,22 @@ def _build_vqd_result() -> VQDResult:
 
     @staticmethod
     def _update_vqd_result(
-        result: VQDResult, opt_result: OptimizerResult, eval_time, ansatz
+        result: VQDResult, opt_result: OptimizerResult, eval_time, ansatz, optimal_point
     ) -> VQDResult:
         result.optimal_points = (
-            np.concatenate([result.optimal_points, [opt_result.x]])
+            np.concatenate([result.optimal_points, [optimal_point["x"]]])
             if len(result.optimal_points) > 0
-            else np.array([opt_result.x])
+            else np.array([optimal_point["x"]])
         )
         result.optimal_parameters.append(
-            dict(zip(ansatz.parameters, cast(np.ndarray, opt_result.x)))
+            dict(zip(ansatz.parameters, cast(np.ndarray, optimal_point["x"])))
+        )
+        result.optimal_values = np.concatenate(
+            [result.optimal_values, [optimal_point["optimal_value"]]]
         )
-        result.optimal_values = np.concatenate([result.optimal_values, [opt_result.fun]])
         result.cost_function_evals = np.concatenate([result.cost_function_evals, [opt_result.nfev]])
         result.optimizer_times = np.concatenate([result.optimizer_times, [eval_time]])
-        result.eigenvalues.append(opt_result.fun + 0j)  # type: ignore[attr-defined]
+        result.eigenvalues.append(optimal_point["eigenvalue"] + 0j)  # type: ignore[attr-defined]
         result.optimizer_results.append(opt_result)
         result.optimal_circuits.append(ansatz)
         return result

test/eigensolvers/test_vqd.py

@@ -59,8 +59,8 @@ def setUp(self):
 
         self.estimator = Estimator()
         self.estimator_shots = Estimator(options={"shots": 1024, "seed": self.seed})
-        self.fidelity = ComputeUncompute(Sampler())
-        self.betas = [50, 50]
+        self.fidelity = ComputeUncompute(Sampler(options={"shots": 100_000, "seed": self.seed}))
+        self.betas = [3]
 
     @data(H2_SPARSE_PAULI)
     def test_basic_operator(self, op):
@@ -105,7 +105,14 @@ def test_basic_operator(self, op):
 
     def test_full_spectrum(self):
         """Test obtaining all eigenvalues."""
-        vqd = VQD(self.estimator, self.fidelity, self.ryrz_wavefunction, optimizer=L_BFGS_B(), k=4)
+        vqd = VQD(
+            self.estimator,
+            self.fidelity,
+            self.ryrz_wavefunction,
+            optimizer=COBYLA(),
+            k=4,
+            betas=[3, 3, 3],
+        )
         result = vqd.compute_eigenvalues(H2_SPARSE_PAULI)
         np.testing.assert_array_almost_equal(
             result.eigenvalues.real, self.h2_energy_excited, decimal=2
@@ -190,9 +197,8 @@ def store_intermediate_result(eval_count, parameters, mean, metadata, step):
             self.assertTrue(all(isinstance(param, float) for param in params))
 
         ref_eval_count = [1, 2, 3, 1, 2, 3]
-        ref_mean = [-1.07, -1.45, -1.37, 37.43, 48.55, 28.94]
-        # new ref_mean for statevector simulator. The old unit test was on qasm
-        # and the ref_mean values were slightly different.
+        ref_mean = [-1.07, -1.45, -1.36, 1.24, 1.55, 1.07]
+        # new ref_mean since the betas were changed
 
         ref_step = [1, 1, 1, 2, 2, 2]
 
@@ -208,15 +214,15 @@ def test_vqd_optimizer(self, op):
             estimator=self.estimator,
             fidelity=self.fidelity,
             ansatz=RealAmplitudes(),
-            optimizer=SLSQP(),
+            optimizer=COBYLA(),
             k=2,
             betas=self.betas,
         )
 
         def run_check():
             result = vqd.compute_eigenvalues(operator=op)
             np.testing.assert_array_almost_equal(
-                result.eigenvalues.real, self.h2_energy_excited[:2], decimal=3
+                result.eigenvalues.real, self.h2_energy_excited[:2], decimal=2
             )
 
         run_check()
@@ -225,11 +231,11 @@ def run_check():
             run_check()
 
         with self.subTest("Optimizer replace"):
-            vqd.optimizer = L_BFGS_B()
+            vqd.optimizer = SPSA()
             run_check()
 
         with self.subTest("Batched optimizer replace"):
-            vqd.optimizer = SLSQP(maxiter=60, max_evals_grouped=10)
+            vqd.optimizer = COBYLA(maxiter=60, max_evals_grouped=10)
             run_check()
 
         with self.subTest("SPSA replace"):
@@ -243,7 +249,7 @@ def run_check():
     def test_optimizer_list(self, op):
         """Test sending an optimizer list"""
 
-        optimizers = [SLSQP(), L_BFGS_B()]
+        optimizers = [COBYLA(), SPSA()]
         initial_point_1 = [
             1.70256666,
             -5.34843975,
@@ -287,7 +293,7 @@ def test_aux_operators_list(self, op):
             estimator=self.estimator,
             fidelity=self.fidelity,
             ansatz=wavefunction,
-            optimizer=SLSQP(),
+            optimizer=COBYLA(),
             k=2,
             betas=self.betas,
         )
@@ -340,7 +346,7 @@ def test_aux_operators_dict(self, op):
             estimator=self.estimator,
             fidelity=self.fidelity,
             ansatz=wavefunction,
-            optimizer=SLSQP(),
+            optimizer=COBYLA(),
             betas=self.betas,
         )
 
@@ -440,6 +446,28 @@ def test_aux_operator_std_dev(self, op):
         self.assertIsInstance(result.aux_operators_evaluated[0][2][1], dict)
         self.assertIsInstance(result.aux_operators_evaluated[0][3][1], dict)
 
+    def test_convergence_threshold(self):
+        """Test the convergence threshold raises an error if and only if too high"""
+        vqd = VQD(
+            self.estimator,
+            self.fidelity,
+            RealAmplitudes(),
+            SLSQP(),
+            k=2,
+            betas=self.betas,
+            convergence_threshold=1e-3,
+        )
+        with self.subTest("Failed convergence"):
+            with self.assertRaises(AlgorithmError):
+                vqd.compute_eigenvalues(operator=H2_SPARSE_PAULI)
+
+        with self.subTest("Convergence accepted"):
+            vqd.convergence_threshold = 1e-1
+            result = vqd.compute_eigenvalues(operator=H2_SPARSE_PAULI)
+            np.testing.assert_array_almost_equal(
+                result.eigenvalues.real, self.h2_energy_excited[:2], decimal=1
+            )
+
 
 if __name__ == "__main__":
     unittest.main()