Adds observable evaluator with primitives.

qiskit/algorithms/observables_evaluator.py

@@ -0,0 +1,162 @@
+# 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 __future__ import annotations
+
+from typing import Tuple, List
+
+import numpy as np
+
+from qiskit import QuantumCircuit
+from qiskit.opflow import (
+    PauliSumOp,
+)
+from . import AlgorithmError
+from .list_or_dict import ListOrDict
+from ..primitives import EstimatorResult, BaseEstimator
+from ..quantum_info.operators.base_operator import BaseOperator
+
+
+def eval_observables(
+    estimator: BaseEstimator,
+    quantum_state: QuantumCircuit,
+    observables: ListOrDict[BaseOperator | PauliSumOp],
+    threshold: float = 1e-12,
+) -> ListOrDict[Tuple[complex, complex]]:
+    """
+    Accepts a sequence 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:
+        estimator: An estimator primitive 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.
+        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.
+        AlgorithmError: If a primitive job is not successful.
+    """
+
+    if (
+        isinstance(quantum_state, QuantumCircuit)  # State 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."
+        )
+    if isinstance(observables, dict):
+        observables_list = list(observables.values())
+    else:
+        observables_list = observables
+
+    observables_list = _handle_zero_ops(observables_list)
+    quantum_state = [quantum_state] * len(observables)
+    try:
+        estimator_job = estimator.run(quantum_state, observables_list)
+        expectation_values = estimator_job.result().values
+    except Exception as exc:
+        raise AlgorithmError("The primitive job failed!") from exc
+
+    std_devs = _compute_std_devs(estimator_job, len(expectation_values))
+
+    # Discard values below threshold
+    observables_means = expectation_values * (np.abs(expectation_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.
+    return _prepare_result(observables_results, observables)
+
+
+def _handle_zero_ops(
+    observables_list: list[BaseOperator | PauliSumOp],
+) -> list[BaseOperator | PauliSumOp]:
+    """Replaces all occurrence of operators equal to 0 in the list with an equivalent ``PauliSumOp``
+    operator."""
+    zero_op = PauliSumOp.from_list([("I" * observables_list[0].num_qubits, 0)])
+    for ind, observable in enumerate(observables_list):
+        if observable == 0:
+            observables_list[ind] = zero_op
+    return observables_list
+
+
+def _prepare_result(
+    observables_results: List[Tuple[complex, complex]],
+    observables: ListOrDict[BaseOperator | PauliSumOp],
+) -> ListOrDict[Tuple[complex, complex]]:
+    """
+    Prepares a list of eigenvalues and standard deviations 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(
+    estimator_result: EstimatorResult,
+    results_length: int,
+) -> List[complex | None]:
+    """
+    Calculates a list of standard deviations from expectation values of observables provided. If
+    there is no variance data available from a primitive, the standard deviation values will be set
+    to ``None``.
+
+    Args:
+        estimator_result: An estimator result.
+        results_length: Number of expectation values calculated.
+
+    Returns:
+        A list of standard deviations.
+    """
+    if not estimator_result.metadata:
+        return [0] * results_length
+
+    std_devs = []
+    for metadata in estimator_result.metadata:
+        if metadata and "variance" in metadata.keys() and "shots" in metadata.keys():
+            variance = metadata["variance"]
+            shots = metadata["shots"]
+            if variance is None or shots is None:
+                std_devs.append(None)
+            else:
+                std_devs.append(np.sqrt(variance / shots))
+        else:
+            std_devs.append(0)
+
+    return std_devs

releasenotes/notes/observable-eval-primitives-e1fd989e15c7760c.yaml

@@ -0,0 +1,9 @@
+---
+features:
+  - |
+    Added `~qiskit.algorithms.observables_evaluator` with
+    :class:`~qiskit.primitives.BaseEstimator` as ``init`` parameter. It will soon replace
+    `~qiskit.algorithms.aux_ops_evaluator`.
+deprecations:
+  - |
+    Using `~qiskit.algorithms.aux_ops_evaluator` will now issue a ``PendingDeprecationWarning``.

test/python/algorithms/test_observables_evaluator.py

@@ -0,0 +1,142 @@
+# 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."""
+from __future__ import annotations
+import unittest
+from typing import Tuple
+
+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.quantum_info.operators.base_operator import BaseOperator
+from qiskit.algorithms.observables_evaluator import eval_observables
+from qiskit.primitives import Estimator
+from qiskit.quantum_info import Statevector
+from qiskit import QuantumCircuit
+from qiskit.circuit.library import EfficientSU2
+from qiskit.opflow import (
+    PauliSumOp,
+    X,
+    Y,
+)
+from qiskit.utils import algorithm_globals
+
+
+@ddt
+class TestObservablesEvaluator(QiskitAlgorithmsTestCase):
+    """Tests evaluator of auxiliary operators for algorithms."""
+
+    def setUp(self):
+        super().setUp()
+        self.seed = 50
+        algorithm_globals.random_seed = self.seed
+
+        self.threshold = 1e-8
+
+    def get_exact_expectation(
+        self, ansatz: QuantumCircuit, observables: ListOrDict[BaseOperator | PauliSumOp]
+    ):
+        """
+        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: QuantumCircuit,
+        decimal: int,
+        observables: ListOrDict[BaseOperator | PauliSumOp],
+        estimator: Estimator,
+    ):
+        result = eval_observables(estimator, quantum_state, observables, 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[BaseOperator | PauliSumOp]):
+        """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)
+        states = bound_ansatz
+        expected_result = self.get_exact_expectation(bound_ansatz, observables)
+        estimator = Estimator()
+        decimal = 6
+        self._run_test(
+            expected_result,
+            states,
+            decimal,
+            observables,
+            estimator,
+        )
+
+    def test_eval_observables_zero_op(self):
+        """Tests if a zero operator is handled correctly."""
+        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)
+        state = bound_ansatz
+        estimator = Estimator()
+        observables = [(X ^ X) + (Y ^ Y), 0]
+        result = eval_observables(estimator, state, observables, self.threshold)
+        expected_result = [(0.015607318055509564, 0), (0.0, 0)]
+        np.testing.assert_array_almost_equal(result, expected_result)
+
+
+if __name__ == "__main__":
+    unittest.main()