Implements high-level functions for QSE (#1902)

* Adds execute_with_qse, qse_decorator and mitigate_executor

* fix type checking list -> List

* adds Executor as a possible type for qse

---------

Co-authored-by: Sahmoud <sahmoua@amazon.com>

mitiq/qse/__init__.py

@@ -3,4 +3,8 @@
 # This source code is licensed under the GPL license (v3) found in the
 # LICENSE file in the root directory of this source tree.
 
-from mitiq.qse.qse_utils import get_projector
+from mitiq.qse.qse import execute_with_qse, mitigate_executor, qse_decorator
+from mitiq.qse.qse_utils import (
+    get_projector,
+    get_expectation_value_for_observable,
+)

mitiq/qse/qse.py

@@ -0,0 +1,161 @@
+# Copyright (C) Unitary Fund
+#
+# This source code is licensed under the GPL license (v3) found in the
+# LICENSE file in the root directory of this source tree.
+
+
+"""High-level Quantum Susbapce Expansion tools."""
+
+from functools import wraps
+from typing import Callable, Sequence, Dict, List, Union
+
+from mitiq import Executor, Observable, QPROGRAM, QuantumResult, PauliString
+from mitiq.qse.qse_utils import (
+    get_projector,
+    get_expectation_value_for_observable,
+)
+
+
+def execute_with_qse(
+    circuit: QPROGRAM,
+    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
+    check_operators: Sequence[PauliString],
+    code_hamiltonian: Observable,
+    observable: Observable,
+    pauli_string_to_expectation_cache: Dict[PauliString, complex] = {},
+) -> float:
+    """Function for the calculation of an observable from some circuit of
+        interest to be mitigated with Quantum Subspace Expansion
+    Args:
+        circuit: Quantum program to execute with error mitigation.
+        executor: Executes a circuit and returns a `QuantumResult`.
+        check_operators: List of check operators that define the
+        stabilizer code space.
+        code_hamiltonian: Hamiltonian of the code space.
+        observable: Observable to compute the mitigated expectation value of.
+        pauli_string_to_expectation_cache: Cache for expectation values of
+        Pauli strings used to compute the projector and the observable.
+    Returns:
+        The expectation value estimated with QSE.
+    """
+    projector = get_projector(
+        circuit,
+        executor,
+        check_operators,
+        code_hamiltonian,
+        pauli_string_to_expectation_cache,
+    )
+    # Compute the expectation value of the observable: <P O P>
+    pop = get_expectation_value_for_observable(
+        circuit,
+        executor,
+        projector * observable * projector,
+        pauli_string_to_expectation_cache,
+    )
+    # Compute the normalization factor: <P P>
+    pp = get_expectation_value_for_observable(
+        circuit,
+        executor,
+        projector * projector,
+        pauli_string_to_expectation_cache,
+    )
+    return pop / pp
+
+
+def mitigate_executor(
+    executor: Callable[[QPROGRAM], QuantumResult],
+    check_operators: Sequence[PauliString],
+    code_hamiltonian: Observable,
+    observable: Observable,
+    pauli_string_to_expectation_cache: Dict[PauliString, complex] = {},
+) -> Callable[[QPROGRAM], float]:
+    """Returns a modified version of the input 'executor' which is
+    error-mitigated with zero-noise extrapolation (ZNE).
+
+    Args:
+        circuit: Quantum program to execute with error mitigation.
+        executor: Executes a circuit and returns a `QuantumResult`.
+        check_operators: List of check operators that define the
+        stabilizer code space.
+        code_hamiltonian: Hamiltonian of the code space.
+        observable: Observable to compute the mitigated expectation value for.
+        pauli_string_to_expectation_cache: Cache for expectation values of
+        Pauli strings used to compute the projector and the observable.
+        share_cache: Only applicable for batched executors. If True, the
+        cache is shared between the all circuits in the batch.
+    Returns:
+        The error-mitigated version of the input executor.
+    """
+    executor_obj = Executor(executor)
+    if not executor_obj.can_batch:
+
+        @wraps(executor)
+        def new_executor(circuit: QPROGRAM) -> float:
+            return execute_with_qse(
+                circuit,
+                executor,
+                check_operators,
+                code_hamiltonian,
+                observable,
+                pauli_string_to_expectation_cache,
+            )
+
+    else:
+
+        @wraps(executor)
+        def new_executor(circuits: List[QPROGRAM]) -> List[float]:
+            return [
+                execute_with_qse(
+                    circuit,
+                    executor,
+                    check_operators,
+                    code_hamiltonian,
+                    observable,
+                    pauli_string_to_expectation_cache,
+                )
+                for circuit in circuits
+            ]
+
+    return new_executor
+
+
+def qse_decorator(
+    check_operators: Sequence[PauliString],
+    code_hamiltonian: Observable,
+    observable: Observable,
+    pauli_string_to_expectation_cache: Dict[PauliString, complex] = {},
+) -> Callable[
+    [Callable[[QPROGRAM], QuantumResult]], Callable[[QPROGRAM], float]
+]:
+    """Decorator which adds an error-mitigation layer based on quantum
+    subspace expansion (QSE) to an executor function, i.e., a function which
+    executes a quantum circuit with an arbitrary backend and returns a
+    ``QuantumResult``.
+
+    Args:
+        check_operators: List of check operators that define the
+        stabilizer code space.
+        code_hamiltonian: Hamiltonian of the code space.
+        observable: Observable to compute the mitigated expectation value of.
+        pauli_string_to_expectation_cache: Cache for expectation values of
+        Pauli strings used to compute the projector and the observable.
+        share_cache: Only applicable for batched executors. If True, the
+        cache is shared between the all circuits in the batch.
+
+    Returns:
+        The error-mitigating decorator to be applied to an executor function.
+    """
+
+    def decorator(
+        executor: Callable[[QPROGRAM], QuantumResult],
+    ) -> Callable[[QPROGRAM], float]:
+        val = mitigate_executor(
+            executor,
+            check_operators,
+            code_hamiltonian,
+            observable,
+            pauli_string_to_expectation_cache,
+        )
+        return val
+
+    return decorator

mitiq/qse/qse_utils.py

@@ -5,9 +5,8 @@
 
 """Functions for computing the projector for subspace expansion."""
 
-from typing import Callable, Sequence, Union
-from mitiq import Observable, QPROGRAM, QuantumResult, PauliString
-from typing import Dict
+from typing import Callable, Sequence, Union, Dict, List
+from mitiq import Observable, QPROGRAM, QuantumResult, PauliString, Executor
 import numpy as np
 import numpy.typing as npt
 from scipy.linalg import eigh
@@ -16,16 +15,16 @@
 
 def get_projector(
     circuit: QPROGRAM,
-    executor: Callable[[QPROGRAM], QuantumResult],
+    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
     check_operators: Sequence[PauliString],
     code_hamiltonian: Observable,
+    pauli_string_to_expectation_cache: Dict[PauliString, complex] = {},
 ) -> Observable:
     """Computes the projector onto the code space defined by the
     check_operators provided that minimizes the code_hamiltonian.
 
     Returns: Projector as an Observable.
     """
-    pauli_string_to_expectation_cache: Dict[PauliString, complex] = {}
     S = _compute_overlap_matrix(
         circuit, executor, check_operators, pauli_string_to_expectation_cache
     )
@@ -48,24 +47,33 @@ def get_projector(
     return projector
 
 
-def _get_expectation_value_for_observable(
+def get_expectation_value_for_observable(
     circuit: QPROGRAM,
-    executor: Callable[[QPROGRAM], QuantumResult],
+    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
     observable: Union[PauliString, Observable],
     pauli_string_to_expectation_cache: Dict[PauliString, complex] = {},
 ) -> float:
+    """Provide pauli_string_to_expectation_cache if you want to take advantage
+    of caching.
+
+    This function modifies pauli_string_to_expectation_cache in place.
+    """
+
     def get_expectation_value_for_one_pauli(
         pauli_string: PauliString,
     ) -> float:
         cache_key = pauli_string.with_coeff(1)
-        pauli_string_to_expectation_cache[cache_key] = Observable(
-            cache_key
-        ).expectation(circuit, executor)
+        pauli_string_to_expectation_cache[cache_key] = final_executor.evaluate(
+            circuit, Observable(cache_key)
+        )[0]
         return (
             pauli_string_to_expectation_cache[cache_key] * pauli_string.coeff
         ).real
 
     total_expectation_value_for_observable = 0.0
+    final_executor: Executor = (
+        executor if isinstance(executor, Executor) else Executor(executor)
+    )
 
     if isinstance(observable, PauliString):
         pauli_string = observable
@@ -82,16 +90,16 @@ def get_expectation_value_for_one_pauli(
 
 def _compute_overlap_matrix(
     circuit: QPROGRAM,
-    executor: Callable[[QPROGRAM], QuantumResult],
+    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
     check_operators: Sequence[PauliString],
     pauli_string_to_expectation_cache: Dict[PauliString, complex] = {},
 ) -> npt.NDArray[np.float64]:
-    S: list[list[float]] = []
+    S: List[List[float]] = []
     # S_ij = <Ψ|Mi† Mj|Ψ>
     for i in range(len(check_operators)):
         S.append([])
         for j in range(len(check_operators)):
-            sij = _get_expectation_value_for_observable(
+            sij = get_expectation_value_for_observable(
                 circuit,
                 executor,
                 check_operators[i] * check_operators[j],
@@ -103,18 +111,18 @@ def _compute_overlap_matrix(
 
 def _compute_hamiltonian_overlap_matrix(
     circuit: QPROGRAM,
-    executor: Callable[[QPROGRAM], QuantumResult],
+    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
     check_operators: Sequence[PauliString],
     code_hamiltonian: Observable,
     pauli_string_to_expectation_cache: Dict[PauliString, complex] = {},
 ) -> npt.NDArray[np.float64]:
-    H: list[list[float]] = []
+    H: List[List[float]] = []
     # H_ij = <Ψ|Mi† H Mj|Ψ>
     for i in range(len(check_operators)):
         H.append([])
         for j in range(len(check_operators)):
             H[-1].append(
-                _get_expectation_value_for_observable(
+                get_expectation_value_for_observable(
                     circuit,
                     executor,
                     check_operators[i] * code_hamiltonian * check_operators[j],