Fix VQD callback not being used (#8576)

* Added support for non-hermitian operators in AerPauliExpectation

Fixes #415

QEOM creates a dictionary of operators to evaluate on the groundstate.
When using the noisy simulators (qasm_simulator or aer_simulator), one
could either use PauliExpectation (with noise) or AerPauliExpectation
(without noise). PauliExpectation works with non-hermitian operators
but internal methods of AerPauliExpectation raised an Error.
This is a workaround to this limitation.
Note that using include_custom=True on qasm allows the VQE to use a
local AerPauliExpectation without using the "expectation" input.
This does not apply to QEOM and one should explicitly define the
"expectation" input of the VQE for it to apply globally.

* Added support for non-hermitian operators in AerPauliExpectation

Fixes #415

QEOM creates a dictionary of operators to evaluate on the groundstate.
When using the noisy simulators (qasm_simulator or aer_simulator), one
could either use PauliExpectation (with noise) or AerPauliExpectation
(without noise). PauliExpectation works with non-hermitian operators
but internal methods of AerPauliExpectation raised an Error.
This is a workaround to this limitation.
Note that using include_custom=True on qasm allows the VQE to use a
local AerPauliExpectation without using the "expectation" input.
This does not apply to QEOM and one should explicitly define the
"expectation" input of the VQE for it to apply globally.

* Add a test case for non-hermitian operators.

* Add a test case for non-hermitian operators.

* Add a test case for non-hermitian operators.

* Update test/python/opflow/test_aer_pauli_expectation.py

Co-authored-by: Julien Gacon <gaconju@gmail.com>

* Update aer_pauli_expectation.py

Use a generator instead of list

* Update qiskit/opflow/expectations/aer_pauli_expectation.py

Co-authored-by: Julien Gacon <gaconju@gmail.com>

* Update releasenotes/notes/add-support-non-hermitian-op-aerpauliexpectation-653d8e16de4eca07.yaml

Co-authored-by: Julien Gacon <gaconju@gmail.com>

* Add a test case for PauliOp

* Change the test cases from using ~StateFn() to using StateFn(, is_measurement=True)

* Fix the formatting

* Working point for QEOM

* Small changes + Release note

* Undesired change

* Undesired change

* Indentation error

* Parenthesis

* Minor changes in docstring

* Minor changes in docstring

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

qiskit/algorithms/eigen_solvers/vqd.py

@@ -108,10 +108,11 @@ class VQD(VariationalAlgorithm, Eigensolver):
         Args:
             ansatz: A parameterized circuit used as ansatz for the wave function.
             k: the number of eigenvalues to return. Returns the lowest k eigenvalues.
-            betas: beta parameter in the VQD paper. Should have size k -1, the number of excited states.
+            betas: beta parameters in the VQD paper.
-                It is a hyperparameter that balances the contribution of the overlap
+                Should have length k - 1, with k the number of excited states.
-                term to the cost function and has a default value computed as
+                These hyperparameters balance the contribution of each overlap term to the cost
-                mean square sum of coefficients of observable.
+                function and have a default value computed as the mean square sum of the
+                coefficients of the observable.
             optimizer: A classical optimizer. Can either be a Qiskit optimizer or a callable
                 that takes an array as input and returns a Qiskit or SciPy optimization result.
             initial_point: An optional initial point (i.e. initial parameter values)
@@ -141,8 +142,8 @@ class VQD(VariationalAlgorithm, Eigensolver):
             callback: 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 for its current set of parameters as it works towards the minimum.
-                These are: the evaluation count, the optimizer parameters for the
+                These are: the evaluation count, the optimizer parameters for the ansatz, the
-                ansatz, the evaluated mean and the evaluated standard deviation.`
+                evaluated mean, the evaluated standard deviation, and the current step.
             quantum_instance: Quantum Instance or Backend
 
         """
@@ -264,12 +265,12 @@ class VQD(VariationalAlgorithm, Eigensolver):
             self.expectation = None
 
     @property
-    def callback(self) -> Optional[Callable[[int, np.ndarray, float, float], None]]:
+    def callback(self) -> Optional[Callable[[int, np.ndarray, float, float, int], None]]:
         """Returns callback"""
         return self._callback
 
     @callback.setter
-    def callback(self, callback: Optional[Callable[[int, np.ndarray, float, float], None]]):
+    def callback(self, callback: Optional[Callable[[int, np.ndarray, float, float, int], None]]):
         """Sets callback"""
         self._callback = callback
 
@@ -475,7 +476,8 @@ class VQD(VariationalAlgorithm, Eigensolver):
         aux_op_results = zip(aux_op_means, std_devs)
 
         # Return None eigenvalues for None operators if aux_operators is a list.
-        # None operators are already dropped in compute_minimum_eigenvalue if aux_operators is a dict.
+        # None operators are already dropped in compute_minimum_eigenvalue if aux_operators is a
+        # dict.
         if isinstance(aux_operators, list):
             aux_operator_eigenvalues = [None] * len(aux_operators)
             key_value_iterator = enumerate(aux_op_results)
@@ -608,7 +610,8 @@ class VQD(VariationalAlgorithm, Eigensolver):
             if step == 1:
 
                 logger.info(
-                    "Ground state optimization complete in %s seconds.\nFound opt_params %s in %s evals",
+                    "Ground state optimization complete in %s seconds.\n"
+                    "Found opt_params %s in %s evals",
                     eval_time,
                     result.optimal_point,
                     self._eval_count,
@@ -616,7 +619,8 @@ class VQD(VariationalAlgorithm, Eigensolver):
             else:
                 logger.info(
                     (
-                        "%s excited state optimization complete in %s s.\nFound opt_parms %s in %s evals"
+                        "%s excited state optimization complete in %s s.\n"
+                        "Found opt_params %s in %s evals"
                     ),
                     str(step - 1),
                     eval_time,
@@ -624,7 +628,7 @@ class VQD(VariationalAlgorithm, Eigensolver):
                     self._eval_count,
                 )
 
-        # To match the siignature of NumpyEigenSolver Result
+        # To match the signature of NumpyEigenSolver Result
         result.eigenstates = ListOp([StateFn(vec) for vec in result.eigenstates])
         result.eigenvalues = np.array(result.eigenvalues)
         result.optimal_point = np.array(result.optimal_point)
@@ -649,7 +653,7 @@ class VQD(VariationalAlgorithm, Eigensolver):
         This return value is the objective function to be passed to the optimizer for evaluation.
 
         Args:
-            step: level of enegy being calculated. 0 for ground, 1 for first excited state and so on.
+            step: level of energy being calculated. 0 for ground, 1 for first excited state...
             operator: The operator whose energy to evaluate.
             return_expectation: If True, return the ``ExpectationBase`` expectation converter used
                 in the construction of the expectation value. Useful e.g. to evaluate other
@@ -695,22 +699,29 @@ class VQD(VariationalAlgorithm, Eigensolver):
 
         def energy_evaluation(parameters):
             parameter_sets = np.reshape(parameters, (-1, num_parameters))
-            # Create dict associating each parameter with the lists of parameterization values for it
+            # Dict associating each parameter with the lists of parameterization values for it
             param_bindings = dict(zip(ansatz_params, parameter_sets.transpose().tolist()))
 
             sampled_expect_op = self._circuit_sampler.convert(expect_op, params=param_bindings)
-            mean = np.real(sampled_expect_op.eval())
+            means = np.real(sampled_expect_op.eval())
 
             for state in range(step - 1):
                 sampled_final_op = self._circuit_sampler.convert(
                     overlap_op[state], params=param_bindings
                 )
                 cost = sampled_final_op.eval()
-                mean += np.real(self.betas[state] * np.conj(cost) * cost)
+                means += np.real(self.betas[state] * np.conj(cost) * cost)
 
-            self._eval_count += len(mean)
+            if self._callback is not None:
+                variance = np.real(expectation.compute_variance(sampled_expect_op))
+                estimator_error = np.sqrt(variance / self.quantum_instance.run_config.shots)
+                for i, param_set in enumerate(parameter_sets):
+                    self._eval_count += 1
+                    self._callback(self._eval_count, param_set, means[i], estimator_error[i], step)
+            else:
+                self._eval_count += len(means)
 
-            return mean if len(mean) > 1 else mean[0]
+            return means if len(means) > 1 else means[0]
 
         if return_expectation:
             return energy_evaluation, expectation

releasenotes/notes/make-use-of-callback-in-vqd-99e3c85f03181298.yaml

@@ -0,0 +1,9 @@
+features:
+  - |
+    Add calls of the callback function during the :meth:`energy_evaluation` to track the progress of
+    the algorithm. Also adds a ``step`` argument to the callback to track which eigenstates of the 
+    Hamiltonian is currently being optimized.
+issues:
+  - |
+    The callback function in the :class:`VQD` was defined but never used. 
+

test/python/algorithms/test_vqd.py

@@ -257,13 +257,14 @@ class TestVQD(QiskitAlgorithmsTestCase):
 
     def test_callback(self):
         """Test the callback on VQD."""
-        history = {"eval_count": [], "parameters": [], "mean": [], "std": []}
+        history = {"eval_count": [], "parameters": [], "mean": [], "std": [], "step": []}
 
-        def store_intermediate_result(eval_count, parameters, mean, std):
+        def store_intermediate_result(eval_count, parameters, mean, std, step):
             history["eval_count"].append(eval_count)
             history["parameters"].append(parameters)
             history["mean"].append(mean)
             history["std"].append(std)
+            history["step"].append(step)
 
         optimizer = COBYLA(maxiter=3)
         wavefunction = self.ry_wavefunction
@@ -279,9 +280,20 @@ class TestVQD(QiskitAlgorithmsTestCase):
         self.assertTrue(all(isinstance(count, int) for count in history["eval_count"]))
         self.assertTrue(all(isinstance(mean, float) for mean in history["mean"]))
         self.assertTrue(all(isinstance(std, float) for std in history["std"]))
+        self.assertTrue(all(isinstance(count, int) for count in history["step"]))
         for params in history["parameters"]:
             self.assertTrue(all(isinstance(param, float) for param in params))
 
+        ref_eval_count = [1, 2, 3, 1, 2, 3]
+        ref_mean = [-1.063, -1.457, -1.360, 37.340, 48.543, 28.586]
+        ref_std = [0.011, 0.010, 0.014, 0.011, 0.010, 0.015]
+        ref_step = [1, 1, 1, 2, 2, 2]
+
+        np.testing.assert_array_almost_equal(history["eval_count"], ref_eval_count, decimal=0)
+        np.testing.assert_array_almost_equal(history["mean"], ref_mean, decimal=2)
+        np.testing.assert_array_almost_equal(history["std"], ref_std, decimal=2)
+        np.testing.assert_array_almost_equal(history["step"], ref_step, decimal=0)
+
     def test_reuse(self):
         """Test re-using a VQD algorithm instance."""
         vqd = VQD(k=1)