Moved and refactored `vqe._eval_aux_ops` method. (#7758)

* Extracted and refactored aux_ops_evaluator.py

* Code refactoring

* Implemented unit test.

* Extended unit test.

* Date fixed.

* Reno added.

* Switched to bound ansatz.

* Fix reno.

* Added docs.

* Refactored unit test.

* Lint fixed.

* Code review edits.

* Added unit test cases for dicts.

* Fixed reno reference.

* Improved unit test.

* Fixed quantum_state types and conversion.

* Fixed quantum_state types and conversion.

* use list of expectations instead of exp. of a list

* Added unit tests for Statevector input.

* Fixed test_vqe.py

* Black fix.

* Code review changes.

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

qiskit/algorithms/__init__.py

@@ -182,6 +182,7 @@ Algorithms that estimate the phases of eigenstates of a unitary.
    PhaseEstimationResult
    IterativePhaseEstimation
 
+
 Exceptions
 ==========
 
@@ -189,6 +190,17 @@ Exceptions
    :toctree: ../stubs/
 
    AlgorithmError
+
+
+Utility methods
+---------------
+
+Utility methods used by algorithms.
+
+.. autosummary::
+   :toctree: ../stubs/
+
+   eval_observables
 """
 
 from .algorithm_result import AlgorithmResult
@@ -230,6 +242,7 @@ from .phase_estimators import (
     IterativePhaseEstimation,
 )
 from .exceptions import AlgorithmError
+from .aux_ops_evaluator import eval_observables
 
 __all__ = [
     "AlgorithmResult",
@@ -276,4 +289,5 @@ __all__ = [
     "PhaseEstimationResult",
     "IterativePhaseEstimation",
     "AlgorithmError",
+    "eval_observables",
 ]

qiskit/algorithms/aux_ops_evaluator.py

@@ -0,0 +1,183 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# 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
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""Evaluator of auxiliary operators for algorithms."""
+
+from typing import Tuple, Union, List
+
+import numpy as np
+
+from qiskit import QuantumCircuit
+from qiskit.algorithms.minimum_eigen_solvers.minimum_eigen_solver import ListOrDict
+from qiskit.opflow import (
+    CircuitSampler,
+    ListOp,
+    StateFn,
+    OperatorBase,
+    ExpectationBase,
+)
+from qiskit.providers import BaseBackend, Backend
+from qiskit.quantum_info import Statevector
+from qiskit.utils import QuantumInstance
+
+
+def eval_observables(
+    quantum_instance: Union[QuantumInstance, BaseBackend, Backend],
+    quantum_state: Union[
+        Statevector,
+        QuantumCircuit,
+        OperatorBase,
+    ],
+    observables: ListOrDict[OperatorBase],
+    expectation: ExpectationBase,
+    threshold: float = 1e-12,
+) -> ListOrDict[Tuple[complex, complex]]:
+    """
+    Accepts a list or a dictionary of operators and calculates their expectation values - means
+    and standard deviations. They are calculated with respect to a quantum state provided. A user
+    can optionally provide a threshold value which filters mean values falling below the threshold.
+
+    Args:
+        quantum_instance: A quantum instance used for calculations.
+        quantum_state: An unparametrized quantum circuit representing a quantum state that
+            expectation values are computed against.
+        observables: A list or a dictionary of operators whose expectation values are to be
+            calculated.
+        expectation: An instance of ExpectationBase which defines a method for calculating
+            expectation values.
+        threshold: A threshold value that defines which mean values should be neglected (helpful for
+            ignoring numerical instabilities close to 0).
+
+    Returns:
+        A list or a dictionary of tuples (mean, standard deviation).
+
+    Raises:
+        ValueError: If a ``quantum_state`` with free parameters is provided.
+    """
+
+    if (
+        isinstance(
+            quantum_state, (QuantumCircuit, OperatorBase)
+        )  # Statevector cannot be parametrized
+        and len(quantum_state.parameters) > 0
+    ):
+        raise ValueError(
+            "A parametrized representation of a quantum_state was provided. It is not "
+            "allowed - it cannot have free parameters."
+        )
+
+    # Create new CircuitSampler to avoid breaking existing one's caches.
+    sampler = CircuitSampler(quantum_instance)
+
+    list_op = _prepare_list_op(quantum_state, observables)
+    observables_expect = expectation.convert(list_op)
+    observables_expect_sampled = sampler.convert(observables_expect)
+
+    # compute means
+    values = np.real(observables_expect_sampled.eval())
+
+    # compute standard deviations
+    std_devs = _compute_std_devs(
+        observables_expect_sampled, observables, expectation, quantum_instance
+    )
+
+    # Discard values below threshold
+    observables_means = values * (np.abs(values) > threshold)
+    # zip means and standard deviations into tuples
+    observables_results = list(zip(observables_means, std_devs))
+
+    # Return None eigenvalues for None operators if observables is a list.
+    # None operators are already dropped in compute_minimum_eigenvalue if observables is a dict.
+
+    return _prepare_result(observables_results, observables)
+
+
+def _prepare_list_op(
+    quantum_state: Union[
+        Statevector,
+        QuantumCircuit,
+        OperatorBase,
+    ],
+    observables: ListOrDict[OperatorBase],
+) -> ListOp:
+    """
+    Accepts a list or a dictionary of operators and converts them to a ``ListOp``.
+
+    Args:
+        quantum_state: An unparametrized quantum circuit representing a quantum state that
+            expectation values are computed against.
+        observables: A list or a dictionary of operators.
+
+    Returns:
+        A ``ListOp`` that includes all provided observables.
+    """
+    if isinstance(observables, dict):
+        observables = list(observables.values())
+
+    state = StateFn(quantum_state)
+    return ListOp([StateFn(obs, is_measurement=True).compose(state) for obs in observables])
+
+
+def _prepare_result(
+    observables_results: List[Tuple[complex, complex]],
+    observables: ListOrDict[OperatorBase],
+) -> ListOrDict[Tuple[complex, complex]]:
+    """
+    Prepares a list or a dictionary of eigenvalues from ``observables_results`` and
+    ``observables``.
+
+    Args:
+        observables_results: A list of of tuples (mean, standard deviation).
+        observables: A list or a dictionary of operators whose expectation values are to be
+            calculated.
+
+    Returns:
+        A list or a dictionary of tuples (mean, standard deviation).
+    """
+    if isinstance(observables, list):
+        observables_eigenvalues = [None] * len(observables)
+        key_value_iterator = enumerate(observables_results)
+    else:
+        observables_eigenvalues = {}
+        key_value_iterator = zip(observables.keys(), observables_results)
+    for key, value in key_value_iterator:
+        if observables[key] is not None:
+            observables_eigenvalues[key] = value
+    return observables_eigenvalues
+
+
+def _compute_std_devs(
+    observables_expect_sampled: OperatorBase,
+    observables: ListOrDict[OperatorBase],
+    expectation: ExpectationBase,
+    quantum_instance: Union[QuantumInstance, BaseBackend, Backend],
+) -> List[complex]:
+    """
+    Calculates a list of standard deviations from expectation values of observables provided.
+
+    Args:
+        observables_expect_sampled: Expected values of observables.
+        observables: A list or a dictionary of operators whose expectation values are to be
+            calculated.
+        expectation: An instance of ExpectationBase which defines a method for calculating
+            expectation values.
+        quantum_instance: A quantum instance used for calculations.
+
+    Returns:
+        A list of standard deviations.
+    """
+    variances = np.real(expectation.compute_variance(observables_expect_sampled))
+    if not isinstance(variances, np.ndarray) and variances == 0.0:
+        # when `variances` is a single value equal to 0., our expectation value is exact and we
+        # manually ensure the variances to be a list of the correct length
+        variances = np.zeros(len(observables), dtype=float)
+    std_devs = np.sqrt(variances / quantum_instance.run_config.shots)
+    return std_devs

qiskit/algorithms/minimum_eigen_solvers/vqe.py

@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2018, 2020.
+# (C) Copyright IBM 2018, 2022.
 #
 # 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
@@ -45,6 +45,7 @@ from ..optimizers import Optimizer, SLSQP, OptimizerResult
 from ..variational_algorithm import VariationalAlgorithm, VariationalResult
 from .minimum_eigen_solver import MinimumEigensolver, MinimumEigensolverResult
 from ..exceptions import AlgorithmError
+from ..aux_ops_evaluator import eval_observables
 
 logger = logging.getLogger(__name__)
 
@@ -480,59 +481,6 @@ class VQE(VariationalAlgorithm, MinimumEigensolver):
     def supports_aux_operators(cls) -> bool:
         return True
 
-    def _eval_aux_ops(
-        self,
-        parameters: np.ndarray,
-        aux_operators: ListOrDict[OperatorBase],
-        expectation: ExpectationBase,
-        threshold: float = 1e-12,
-    ) -> ListOrDict[Tuple[complex, complex]]:
-        # Create new CircuitSampler to avoid breaking existing one's caches.
-        sampler = CircuitSampler(self.quantum_instance)
-
-        if isinstance(aux_operators, dict):
-            list_op = ListOp(list(aux_operators.values()))
-        else:
-            list_op = ListOp(aux_operators)
-
-        aux_op_expect = expectation.convert(
-            StateFn(list_op, is_measurement=True).compose(
-                CircuitStateFn(self.ansatz.bind_parameters(parameters))
-            )
-        )
-        aux_op_expect_sampled = sampler.convert(aux_op_expect)
-
-        # compute means
-        values = np.real(aux_op_expect_sampled.eval())
-
-        # compute standard deviations
-        variances = np.real(expectation.compute_variance(aux_op_expect_sampled))
-        if not isinstance(variances, np.ndarray) and variances == 0.0:
-            # when `variances` is a single value equal to 0., our expectation value is exact and we
-            # manually ensure the variances to be a list of the correct length
-            variances = np.zeros(len(aux_operators), dtype=float)
-        std_devs = np.sqrt(variances / self.quantum_instance.run_config.shots)
-
-        # Discard values below threshold
-        aux_op_means = values * (np.abs(values) > threshold)
-        # zip means and standard deviations into tuples
-        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.
-        if isinstance(aux_operators, list):
-            aux_operator_eigenvalues = [None] * len(aux_operators)
-            key_value_iterator = enumerate(aux_op_results)
-        else:
-            aux_operator_eigenvalues = {}
-            key_value_iterator = zip(aux_operators.keys(), aux_op_results)
-
-        for key, value in key_value_iterator:
-            if aux_operators[key] is not None:
-                aux_operator_eigenvalues[key] = value
-
-        return aux_operator_eigenvalues
-
     def compute_minimum_eigenvalue(
         self, operator: OperatorBase, aux_operators: Optional[ListOrDict[OperatorBase]] = None
     ) -> MinimumEigensolverResult:
@@ -639,7 +587,11 @@ class VQE(VariationalAlgorithm, MinimumEigensolver):
         self._ret = result
 
         if aux_operators is not None:
-            aux_values = self._eval_aux_ops(opt_result.x, aux_operators, expectation=expectation)
+            bound_ansatz = self.ansatz.bind_parameters(result.optimal_point)
+
+            aux_values = eval_observables(
+                self.quantum_instance, bound_ansatz, aux_operators, expectation=expectation
+            )
             result.aux_operator_eigenvalues = aux_values
 
         return result

releasenotes/notes/refactor-aux-operators-79d790f8a693a7c0.yaml

@@ -0,0 +1,8 @@
+---
+other:
+  - |
+    A new convenience method :func:`qiskit.algorithms.eval_observables` to evaluate observables
+    given a bound circuit is added. It originates from a method previously residing in
+    :class:`qiskit.algorithms.VQE`. The method is general enough so that it can be used
+    in other algorithms, for example time evolution algorithms. The method is also refactored to be
+    more modular.

test/python/algorithms/test_aux_ops_evaluator.py

@@ -0,0 +1,160 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# 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
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""Tests evaluator of auxiliary operators for algorithms."""
+
+import unittest
+from typing import Tuple, Union
+
+from test.python.algorithms import QiskitAlgorithmsTestCase
+import numpy as np
+from ddt import ddt, data
+
+from qiskit.algorithms.list_or_dict import ListOrDict
+from qiskit.providers import BaseBackend, Backend
+from qiskit.quantum_info import Statevector
+from qiskit.algorithms import eval_observables
+from qiskit import BasicAer, QuantumCircuit
+from qiskit.circuit.library import EfficientSU2
+from qiskit.opflow import PauliSumOp, X, Z, I, ExpectationFactory, OperatorBase, ExpectationBase
+from qiskit.utils import QuantumInstance, algorithm_globals
+
+
+@ddt
+class TestAuxOpsEvaluator(QiskitAlgorithmsTestCase):
+    """Tests evaluator of auxiliary operators for algorithms."""
+
+    def setUp(self):
+        super().setUp()
+        self.seed = 50
+        algorithm_globals.random_seed = self.seed
+        self.h2_op = (
+            -1.052373245772859 * (I ^ I)
+            + 0.39793742484318045 * (I ^ Z)
+            - 0.39793742484318045 * (Z ^ I)
+            - 0.01128010425623538 * (Z ^ Z)
+            + 0.18093119978423156 * (X ^ X)
+        )
+
+        self.threshold = 1e-8
+        self.backend_names = ["statevector_simulator", "qasm_simulator"]
+
+    def get_exact_expectation(self, ansatz: QuantumCircuit, observables: ListOrDict[OperatorBase]):
+        """
+        Calculates the exact expectation to be used as an expected result for unit tests.
+        """
+        if isinstance(observables, dict):
+            observables_list = list(observables.values())
+        else:
+            observables_list = observables
+
+        # the exact value is a list of (mean, variance) where we expect 0 variance
+        exact = [
+            (Statevector(ansatz).expectation_value(observable), 0)
+            for observable in observables_list
+        ]
+
+        if isinstance(observables, dict):
+            return dict(zip(observables.keys(), exact))
+
+        return exact
+
+    def _run_test(
+        self,
+        expected_result: ListOrDict[Tuple[complex, complex]],
+        quantum_state: Union[QuantumCircuit, Statevector],
+        decimal: int,
+        expectation: ExpectationBase,
+        observables: ListOrDict[OperatorBase],
+        quantum_instance: Union[QuantumInstance, BaseBackend, Backend],
+    ):
+        result = eval_observables(
+            quantum_instance, quantum_state, observables, expectation, self.threshold
+        )
+
+        if isinstance(observables, dict):
+            np.testing.assert_equal(list(result.keys()), list(expected_result.keys()))
+            np.testing.assert_array_almost_equal(
+                list(result.values()), list(expected_result.values()), decimal=decimal
+            )
+        else:
+            np.testing.assert_array_almost_equal(result, expected_result, decimal=decimal)
+
+    @data(
+        [
+            PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)]),
+            PauliSumOp.from_list([("II", 2.0)]),
+        ],
+        [
+            PauliSumOp.from_list([("ZZ", 2.0)]),
+        ],
+        {
+            "op1": PauliSumOp.from_list([("II", 2.0)]),
+            "op2": PauliSumOp.from_list([("II", 0.5), ("ZZ", 0.5), ("YY", 0.5), ("XX", -0.5)]),
+        },
+        {
+            "op1": PauliSumOp.from_list([("ZZ", 2.0)]),
+        },
+    )
+    def test_eval_observables(self, observables: ListOrDict[OperatorBase]):
+        """Tests evaluator of auxiliary operators for algorithms."""
+
+        ansatz = EfficientSU2(2)
+        parameters = np.array(
+            [1.2, 4.2, 1.4, 2.0, 1.2, 4.2, 1.4, 2.0, 1.2, 4.2, 1.4, 2.0, 1.2, 4.2, 1.4, 2.0],
+            dtype=float,
+        )
+
+        bound_ansatz = ansatz.bind_parameters(parameters)
+        expected_result = self.get_exact_expectation(bound_ansatz, observables)
+
+        for backend_name in self.backend_names:
+            shots = 4096 if backend_name == "qasm_simulator" else 1
+            decimal = (
+                1 if backend_name == "qasm_simulator" else 6
+            )  # to accommodate for qasm being imperfect
+            with self.subTest(msg=f"Test {backend_name} backend."):
+                backend = BasicAer.get_backend(backend_name)
+                quantum_instance = QuantumInstance(
+                    backend=backend,
+                    shots=shots,
+                    seed_simulator=self.seed,
+                    seed_transpiler=self.seed,
+                )
+                expectation = ExpectationFactory.build(
+                    operator=self.h2_op,
+                    backend=quantum_instance,
+                )
+
+                with self.subTest(msg="Test QuantumCircuit."):
+                    self._run_test(
+                        expected_result,
+                        bound_ansatz,
+                        decimal,
+                        expectation,
+                        observables,
+                        quantum_instance,
+                    )
+
+                with self.subTest(msg="Test Statevector."):
+                    statevector = Statevector(bound_ansatz)
+                    self._run_test(
+                        expected_result,
+                        statevector,
+                        decimal,
+                        expectation,
+                        observables,
+                        quantum_instance,
+                    )
+
+
+if __name__ == "__main__":
+    unittest.main()

test/python/algorithms/test_vqe.py

@@ -15,6 +15,7 @@
 import logging
 import unittest
 from test.python.algorithms import QiskitAlgorithmsTestCase
+from test.python.transpiler._dummy_passes import DummyAP
 
 from functools import partial
 import numpy as np
@@ -52,7 +53,6 @@ from qiskit.quantum_info import Statevector
 from qiskit.transpiler import PassManager, PassManagerConfig
 from qiskit.transpiler.preset_passmanagers import level_1_pass_manager
 from qiskit.utils import QuantumInstance, algorithm_globals, has_aer
-from ..transpiler._dummy_passes import DummyAP
 
 if has_aer():
     from qiskit import Aer
@@ -616,12 +616,10 @@ class TestVQE(QiskitAlgorithmsTestCase):
         self.assertEqual(len(result.aux_operator_eigenvalues), 2)
         # expectation values
         self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0], 2.0, places=6)
-        self.assertAlmostEqual(result.aux_operator_eigenvalues[1][0], 0.5784419552370315, places=6)
+        self.assertAlmostEqual(result.aux_operator_eigenvalues[1][0], 0.6796875, places=6)
         # standard deviations
         self.assertAlmostEqual(result.aux_operator_eigenvalues[0][1], 0.0)
-        self.assertAlmostEqual(
+        self.assertAlmostEqual(result.aux_operator_eigenvalues[1][1], 0.02534712219145965, places=6)
-            result.aux_operator_eigenvalues[1][1], 0.015183867579396111, places=6
-        )
 
         # Go again with additional None and zero operators
         aux_ops = [*aux_ops, None, 0]
@@ -629,12 +627,14 @@ class TestVQE(QiskitAlgorithmsTestCase):
         self.assertEqual(len(result.aux_operator_eigenvalues), 4)
         # expectation values
         self.assertAlmostEqual(result.aux_operator_eigenvalues[0][0], 2.0, places=6)
-        self.assertAlmostEqual(result.aux_operator_eigenvalues[1][0], 0.56640625, places=6)
+        self.assertAlmostEqual(result.aux_operator_eigenvalues[1][0], 0.57421875, places=6)
         self.assertEqual(result.aux_operator_eigenvalues[2][0], 0.0)
         self.assertEqual(result.aux_operator_eigenvalues[3][0], 0.0)
         # # standard deviations
         self.assertAlmostEqual(result.aux_operator_eigenvalues[0][1], 0.0)
-        self.assertAlmostEqual(result.aux_operator_eigenvalues[1][1], 0.01548658094658011, places=6)
+        self.assertAlmostEqual(
+            result.aux_operator_eigenvalues[1][1], 0.026562146577166837, places=6
+        )
         self.assertAlmostEqual(result.aux_operator_eigenvalues[2][1], 0.0)
         self.assertAlmostEqual(result.aux_operator_eigenvalues[3][1], 0.0)