Fix betas calculation on VQD and async. cost evaluation (#9245)

* Fix betas

* Add async

* Fix seed test

* Fix lint

* Add reno

* Add enumerate

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

* Apply suggestions Julien

* Fix lint

Co-authored-by: Julien Gacon <[email protected]>
Co-authored-by: mergify[bot] <37929162+mergify[bot]@users.noreply.github.com>
(cherry picked from commit f339f71)

qiskit/algorithms/eigensolvers/vqd.py

@@ -24,11 +24,12 @@
 
 import numpy as np
 
+from qiskit.algorithms.state_fidelities import BaseStateFidelity
 from qiskit.circuit import QuantumCircuit
-from qiskit.quantum_info.operators.base_operator import BaseOperator
 from qiskit.opflow import PauliSumOp
 from qiskit.primitives import BaseEstimator
-from qiskit.algorithms.state_fidelities import BaseStateFidelity
+from qiskit.quantum_info.operators.base_operator import BaseOperator
+from qiskit.quantum_info import SparsePauliOp
 
 from ..list_or_dict import ListOrDict
 from ..optimizers import Optimizer, Minimizer, OptimizerResult
@@ -176,12 +177,12 @@ def _check_operator_ansatz(self, operator: BaseOperator | PauliSumOp):
                 # try to set the number of qubits on the ansatz, if possible
                 try:
                     self.ansatz.num_qubits = operator.num_qubits
-                except AttributeError as ex:
+                except AttributeError as exc:
                     raise AlgorithmError(
                         "The number of qubits of the ansatz does not match the "
                         "operator, and the ansatz does not allow setting the "
                         "number of qubits using `num_qubits`."
-                    ) from ex
+                    ) from exc
 
     @classmethod
     def supports_aux_operators(cls) -> bool:
@@ -205,7 +206,7 @@ def compute_eigenvalues(
 
         # We need to handle the array entries being zero or Optional i.e. having value None
         if aux_operators:
-            zero_op = PauliSumOp.from_list([("I" * self.ansatz.num_qubits, 0)])
+            zero_op = SparsePauliOp.from_list([("I" * self.ansatz.num_qubits, 0)])
 
             # Convert the None and zero values when aux_operators is a list.
             # Drop None and convert zero values when aux_operators is a dict.
@@ -225,17 +226,21 @@ def compute_eigenvalues(
             aux_operators = None
 
         if self.betas is None:
+
             if isinstance(operator, PauliSumOp):
-                upper_bound = abs(operator.coeff) * sum(
-                    abs(operation.coeff) for operation in operator
-                )
-                betas = [upper_bound * 10] * (self.k)
-                logger.info("beta autoevaluated to %s", betas[0])
-            else:
+                operator = operator.coeff * operator.primitive
+
+            try:
+                upper_bound = sum(np.abs(operator.coeffs))
+
+            except Exception as exc:
                 raise NotImplementedError(
-                    r"Beta autoevaluation is only supported for operators"
-                    f"of type PauliSumOp, found {type(operator)}."
-                )
+                    r"Beta autoevaluation is not supported for operators"
+                    f"of type {type(operator)}."
+                ) from exc
+
+            betas = [upper_bound * 10] * (self.k)
+            logger.info("beta autoevaluated to %s", betas[0])
         else:
             betas = self.betas
 
@@ -334,6 +339,7 @@ def _get_evaluate_energy(
         Args:
             step: level of energy being calculated. 0 for ground, 1 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.
 
         Returns:
@@ -360,27 +366,31 @@ def _get_evaluate_energy(
 
         def evaluate_energy(parameters: np.ndarray) -> np.ndarray | float:
 
-            try:
-                estimator_job = self.estimator.run(
-                    circuits=[self.ansatz], observables=[operator], parameter_values=[parameters]
-                )
-                estimator_result = estimator_job.result()
-                values = estimator_result.values
-
-            except Exception as exc:
-                raise AlgorithmError("The primitive job to evaluate the energy failed!") from exc
+            estimator_job = self.estimator.run(
+                circuits=[self.ansatz], observables=[operator], parameter_values=[parameters]
+            )
 
+            total_cost = 0
             if step > 1:
-                # Compute overlap cost
+                # compute overlap cost
                 fidelity_job = self.fidelity.run(
                     [self.ansatz] * (step - 1),
                     prev_states,
                     [parameters] * (step - 1),
                 )
+
                 costs = fidelity_job.result().fidelities
 
-                for (state, cost) in zip(range(step - 1), costs):
-                    values += np.real(betas[state] * cost)
+                for state, cost in enumerate(costs):
+                    total_cost += np.real(betas[state] * cost)
+
+            try:
+                estimator_result = estimator_job.result()
+
+            except Exception as exc:
+                raise AlgorithmError("The primitive job to evaluate the energy failed!") from exc
+
+            values = estimator_result.values + total_cost
 
             if self.callback is not None:
                 metadata = estimator_result.metadata

releasenotes/notes/fix-vqd-betas-async-df99ab6e26e9da1e.yaml

@@ -0,0 +1,10 @@
+---
+fixes:
+  - |
+    Fixed existing betas autoevaluation code on class
+    :class:`~qiskit.algorithms.eigensolvers.VQD`, added support for
+    :class:`~qiskit.quantum_info.SparsePauliOp` inputs, and fixed
+    the energy evaluation function to leverage the async execution
+    of primitives, by only retrieving the job results after both
+    jobs have been submitted.
+

test/python/algorithms/eigensolvers/test_vqd.py

@@ -26,28 +26,27 @@
     L_BFGS_B,
     SLSQP,
 )
-
+from qiskit.algorithms.state_fidelities import ComputeUncompute
 from qiskit.circuit.library import TwoLocal, RealAmplitudes
 from qiskit.opflow import PauliSumOp
 from qiskit.primitives import Sampler, Estimator
-from qiskit.algorithms.state_fidelities import ComputeUncompute
-from qiskit.utils import algorithm_globals
+from qiskit.quantum_info import SparsePauliOp
 from qiskit.quantum_info.operators import Operator
+from qiskit.utils import algorithm_globals
 
 
-I = PauliSumOp.from_list([("I", 1)])  # pylint: disable=invalid-name
-X = PauliSumOp.from_list([("X", 1)])  # pylint: disable=invalid-name
-Z = PauliSumOp.from_list([("Z", 1)])  # pylint: disable=invalid-name
-
-H2_PAULI = (
-    -1.052373245772859 * (I ^ I)
-    + 0.39793742484318045 * (I ^ Z)
-    - 0.39793742484318045 * (Z ^ I)
-    - 0.01128010425623538 * (Z ^ Z)
-    + 0.18093119978423156 * (X ^ X)
+H2_SPARSE_PAULI = SparsePauliOp.from_list(
+    [
+        ("II", -1.052373245772859),
+        ("IZ", 0.39793742484318045),
+        ("ZI", -0.39793742484318045),
+        ("ZZ", -0.01128010425623538),
+        ("XX", 0.18093119978423156),
+    ]
 )
+H2_OP = Operator(H2_SPARSE_PAULI.to_matrix())
 
-H2_OP = Operator(H2_PAULI.to_matrix())
+H2_PAULI = PauliSumOp(H2_SPARSE_PAULI)
 
 
 @ddt
@@ -68,11 +67,11 @@ def setUp(self):
         self.ry_wavefunction = TwoLocal(rotation_blocks="ry", entanglement_blocks="cz")
 
         self.estimator = Estimator()
-        self.estimator_shots = Estimator(options={"shots": 2048, "seed": self.seed})
+        self.estimator_shots = Estimator(options={"shots": 1024, "seed": self.seed})
         self.fidelity = ComputeUncompute(Sampler())
         self.betas = [50, 50]
 
-    @data(H2_PAULI, H2_OP)
+    @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_basic_operator(self, op):
         """Test the VQD without aux_operators."""
         wavefunction = self.ryrz_wavefunction
@@ -116,7 +115,21 @@ def test_full_spectrum(self):
             result.eigenvalues.real, self.h2_energy_excited, decimal=2
         )
 
-    @data(H2_PAULI, H2_OP)
+    @data(H2_PAULI, H2_SPARSE_PAULI)
+    def test_beta_autoeval(self, op):
+        """Test beta autoevaluation for different operator types."""
+
+        with self.assertLogs(level="INFO") as logs:
+            vqd = VQD(
+                self.estimator_shots, self.fidelity, self.ryrz_wavefunction, optimizer=L_BFGS_B()
+            )
+            _ = vqd.compute_eigenvalues(op)
+
+        # the first log message shows the value of beta[0]
+        beta = float(logs.output[0].split()[-1])
+        self.assertAlmostEqual(beta, 20.40459399499687, 4)
+
+    @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_mismatching_num_qubits(self, op):
         """Ensuring circuit and operator mismatch is caught"""
         wavefunction = QuantumCircuit(1)
@@ -132,7 +145,7 @@ def test_mismatching_num_qubits(self, op):
         with self.assertRaises(AlgorithmError):
             _ = vqd.compute_eigenvalues(operator=op)
 
-    @data(H2_PAULI, H2_OP)
+    @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_missing_varform_params(self, op):
         """Test specifying a variational form with no parameters raises an error."""
         circuit = QuantumCircuit(op.num_qubits)
@@ -147,7 +160,7 @@ def test_missing_varform_params(self, op):
         with self.assertRaises(AlgorithmError):
             vqd.compute_eigenvalues(operator=op)
 
-    @data(H2_PAULI, H2_OP)
+    @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_callback(self, op):
         """Test the callback on VQD."""
         history = {"eval_count": [], "parameters": [], "mean": [], "metadata": [], "step": []}
@@ -163,7 +176,7 @@ def store_intermediate_result(eval_count, parameters, mean, metadata, step):
         wavefunction = self.ry_wavefunction
 
         vqd = VQD(
-            estimator=self.estimator,
+            estimator=self.estimator_shots,
             fidelity=self.fidelity,
             ansatz=wavefunction,
             optimizer=optimizer,
@@ -191,7 +204,7 @@ def store_intermediate_result(eval_count, parameters, mean, metadata, step):
         np.testing.assert_array_almost_equal(history["mean"], ref_mean, decimal=2)
         np.testing.assert_array_almost_equal(history["step"], ref_step, decimal=0)
 
-    @data(H2_PAULI, H2_OP)
+    @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(
@@ -218,7 +231,7 @@ def run_check():
             vqd.optimizer = L_BFGS_B()
             run_check()
 
-    @data(H2_PAULI, H2_OP)
+    @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_aux_operators_list(self, op):
         """Test list-based aux_operators."""
         wavefunction = self.ry_wavefunction
@@ -239,8 +252,8 @@ def test_aux_operators_list(self, op):
         self.assertIsNone(result.aux_operators_evaluated)
 
         # Go again with two auxiliary operators
-        aux_op1 = PauliSumOp.from_list([("II", 2.0)])
-        aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
+        aux_op1 = SparsePauliOp.from_list([("II", 2.0)])
+        aux_op2 = SparsePauliOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
         aux_ops = [aux_op1, aux_op2]
         result = vqd.compute_eigenvalues(op, aux_operators=aux_ops)
         np.testing.assert_array_almost_equal(
@@ -271,7 +284,7 @@ def test_aux_operators_list(self, op):
         self.assertIsInstance(result.aux_operators_evaluated[0][1][1], dict)
         self.assertIsInstance(result.aux_operators_evaluated[0][3][1], dict)
 
-    @data(H2_PAULI, H2_OP)
+    @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_aux_operators_dict(self, op):
         """Test dictionary compatibility of aux_operators"""
         wavefunction = self.ry_wavefunction
@@ -291,8 +304,8 @@ def test_aux_operators_dict(self, op):
         self.assertIsNone(result.aux_operators_evaluated)
 
         # Go again with two auxiliary operators
-        aux_op1 = PauliSumOp.from_list([("II", 2.0)])
-        aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
+        aux_op1 = SparsePauliOp.from_list([("II", 2.0)])
+        aux_op2 = SparsePauliOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
         aux_ops = {"aux_op1": aux_op1, "aux_op2": aux_op2}
         result = vqd.compute_eigenvalues(op, aux_operators=aux_ops)
         self.assertEqual(len(result.eigenvalues), 2)
@@ -325,7 +338,7 @@ def test_aux_operators_dict(self, op):
         self.assertIsInstance(result.aux_operators_evaluated[0]["aux_op2"][1], dict)
         self.assertIsInstance(result.aux_operators_evaluated[0]["zero_operator"][1], dict)
 
-    @data(H2_PAULI, H2_OP)
+    @data(H2_PAULI, H2_OP, H2_SPARSE_PAULI)
     def test_aux_operator_std_dev(self, op):
         """Test non-zero standard deviations of aux operators."""
         wavefunction = self.ry_wavefunction
@@ -348,8 +361,8 @@ def test_aux_operator_std_dev(self, op):
         )
 
         # Go again with two auxiliary operators
-        aux_op1 = PauliSumOp.from_list([("II", 2.0)])
-        aux_op2 = PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
+        aux_op1 = SparsePauliOp.from_list([("II", 2.0)])
+        aux_op2 = SparsePauliOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)])
         aux_ops = [aux_op1, aux_op2]
         result = vqd.compute_eigenvalues(op, aux_operators=aux_ops)
         self.assertEqual(len(result.aux_operators_evaluated), 2)