Add BackendEstimatorV2

qiskit/primitives/__init__.py

@@ -387,14 +387,15 @@
 Primitives API
 ==============
 
-Primitives V2
+Estimator V2
 -------------
 
 .. autosummary::
    :toctree: ../stubs/
 
    BaseEstimatorV2
    StatevectorEstimator
+   BackendEstimatorV2
 
 Sampler V2
 ----------
@@ -473,3 +474,4 @@
 from .sampler import Sampler
 from .statevector_estimator import StatevectorEstimator
 from .statevector_sampler import StatevectorSampler
+from .backend_estimator_v2 import BackendEstimatorV2

qiskit/primitives/backend_estimator.py

@@ -413,14 +413,12 @@ def _paulis2inds(paulis: PauliList) -> list[int]:
     # Treat Z, X, Y the same
     nonid = paulis.z | paulis.x
 
-    inds = [0] * paulis.size
     # bits are packed into uint8 in little endian
     # e.g., i-th bit corresponds to coefficient 2^i
     packed_vals = np.packbits(nonid, axis=1, bitorder="little")
-    for i, vals in enumerate(packed_vals):
-        for j, val in enumerate(vals):
-            inds[i] += val.item() * (1 << (8 * j))
-    return inds
+    power_uint8 = 1 << (8 * np.arange(packed_vals.shape[1], dtype=object))
+    inds = packed_vals @ power_uint8
+    return inds.tolist()
 
 
 def _parity(integer: int) -> int:

qiskit/primitives/backend_estimator_v2.py

@@ -0,0 +1,262 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2024.
+#
+# 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.
+
+"""Estimator V2 implementation for an arbitrary BackendV2 object."""
+
+from __future__ import annotations
+
+from collections import defaultdict
+from collections.abc import Iterable
+
+import numpy as np
+
+from qiskit.circuit import ClassicalRegister, QuantumCircuit, QuantumRegister
+from qiskit.exceptions import QiskitError
+from qiskit.providers import BackendV2
+from qiskit.quantum_info import Pauli, PauliList
+from qiskit.transpiler import PassManager, PassManagerConfig
+from qiskit.transpiler.passes import Optimize1qGatesDecomposition
+
+from .backend_estimator import _pauli_expval_with_variance, _prepare_counts, _run_circuits
+from .base import BaseEstimatorV2
+from .containers import EstimatorPubLike, PrimitiveResult, PubResult
+from .containers.bindings_array import BindingsArray
+from .containers.estimator_pub import EstimatorPub
+from .primitive_job import PrimitiveJob
+
+
+class BackendEstimatorV2(BaseEstimatorV2):
+    """Evaluates expectation value using Pauli rotation gates.
+
+    The :class:`~.BackendEstimator` class is a generic implementation of the
+    :class:`~.BaseEstimator` interface that is used to wrap a :class:`~.BackendV2`
+    (or :class:`~.BackendV1`) object in the :class:`~.BaseEstimator` API. It
+    facilitates using backends that do not provide a native
+    :class:`~.BaseEstimator` implementation in places that work with
+    :class:`~.BaseEstimator`. However,
+    if you're using a provider that has a native implementation of
+    :class:`~.BaseEstimator`, it is a better choice to leverage that native
+    implementation as it will likely include additional optimizations and be
+    a more efficient implementation. The generic nature of this class
+    precludes doing any provider- or backend-specific optimizations.
+
+    Implementation of :class:`BaseEstimatorV2` using a backend.
+
+    This class provides a EstimatorV2 interface from any :class:`~.BackendV2` backend
+    and doesn't do any measurement mitigation, it just computes the expectation values.
+    At present, this implementation is only compatible with Pauli-based observables.
+
+    Each tuple of ``(circuit, observables, <optional> parameter values, <optional> precision)``,
+    called an estimator primitive unified bloc (PUB), produces its own array-based result. The
+    :meth:`~.EstimatorV2.run` method can be given a sequence of pubs to run in one call.
+    """
+
+    _VARIANCE_UPPER_BOUND: float = 4.0
+
+    def __init__(
+        self,
+        *,
+        backend: BackendV2,
+        default_precision: float = 0.0,
+        abelian_grouping: bool = True,
+    ):
+        """
+        Args:
+            backend: Required: the backend to run the primitive on
+            default_precision: Default precision
+            abelian_grouping: Whether the observable should be grouped into commuting
+        """
+        super().__init__()
+        self._backend = backend
+        basis = PassManagerConfig.from_backend(backend).basis_gates
+        self._passmanager = PassManager(
+            [Optimize1qGatesDecomposition(basis=basis, target=backend.target)]
+        )
+        self._default_precision = default_precision
+        self._abelian_grouping = abelian_grouping
+
+    @property
+    def default_precision(self) -> float:
+        """Return the default precision"""
+        return self._default_precision
+
+    @property
+    def ablian_grouping(self) -> bool:
+        """Return whether Abelian grouping is used or not."""
+        return self._abelian_grouping
+
+    @property
+    def backend(self) -> BackendV2:
+        """Returns the backend which this sampler object based on."""
+        return self._backend
+
+    def run(
+        self, pubs: Iterable[EstimatorPubLike], *, precision: float | None = None
+    ) -> PrimitiveJob[PrimitiveResult[PubResult]]:
+        if precision is None:
+            precision = self._default_precision
+        coerced_pubs = [EstimatorPub.coerce(pub, precision) for pub in pubs]
+        self._validate_pubs(coerced_pubs)
+        job = PrimitiveJob(self._run, coerced_pubs)
+        job._submit()
+        return job
+
+    def _validate_pubs(self, pubs: list[EstimatorPub]):
+        for i, pub in enumerate(pubs):
+            if pub.precision == 0.0:
+                raise ValueError(
+                    f"The {i}-th pub has precision 0. But precision should be larger than 0."
+                )
+
+    def _run(self, pubs: list[EstimatorPub]) -> PrimitiveResult[PubResult]:
+        return PrimitiveResult([self._run_pub(pub) for pub in pubs])
+
+    def _run_pub(self, pub: EstimatorPub) -> PubResult:
+        shots = int(np.ceil(self._VARIANCE_UPPER_BOUND / pub.precision**2))
+        circuit = pub.circuit
+        observables = pub.observables
+        parameter_values = pub.parameter_values
+
+        # calculate broadcasting of parameters and observables
+        param_shape = parameter_values.shape
+        param_indices = np.fromiter(np.ndindex(param_shape), dtype=object).reshape(param_shape)
+        bc_param_ind, bc_obs = np.broadcast_arrays(param_indices, observables)
+
+        # calculate expectation values for each pair of parameter value set and pauli
+        param_obs_map = defaultdict(set)
+        for index in np.ndindex(*bc_param_ind.shape):
+            param_index = bc_param_ind[index]
+            param_obs_map[param_index].update(bc_obs[index].keys())
+        expval_map = self._calc_expval_paulis(circuit, parameter_values, param_obs_map, shots)
+
+        # calculate expectation values and standard deviations
+        evs = np.zeros_like(bc_param_ind, dtype=float)
+        variances = np.zeros_like(bc_param_ind, dtype=float)
+        for index in np.ndindex(*bc_param_ind.shape):
+            param_index = bc_param_ind[index]
+            for pauli, coeff in bc_obs[index].items():
+                expval, variance = expval_map[param_index, pauli]
+                evs[index] += expval * coeff
+                variances[index] += variance * coeff**2
+        stds = np.sqrt(variances / shots)
+        data_bin_cls = self._make_data_bin(pub)
+        data_bin = data_bin_cls(evs=evs, stds=stds)
+        return PubResult(data_bin, metadata={"precision": pub.precision})
+
+    def _calc_expval_paulis(
+        self,
+        circuit: QuantumCircuit,
+        parameter_values: BindingsArray,
+        param_obs_map: dict[tuple[int, ...], set[str]],
+        shots: int,
+    ) -> dict[tuple[tuple[int, ...], str], tuple[float, float]]:
+        # generate circuits
+        circuits = []
+        for param_index, pauli_strings in param_obs_map.items():
+            bound_circuit = parameter_values.bind(circuit, param_index)
+            meas_paulis = PauliList(pauli_strings)
+            new_circuits = self._preprocessing(bound_circuit, meas_paulis, param_index)
+            circuits.extend(new_circuits)
+
+        # run circuits
+        result, metadata = _run_circuits(circuits, self._backend, shots=shots)
+
+        # postprocessing results
+        expval_map: dict[tuple[tuple[int, ...], str], tuple[float, float]] = {}
+        counts = _prepare_counts(result)
+        for count, meta in zip(counts, metadata):
+            orig_paulis = meta["orig_paulis"]
+            meas_paulis = meta["meas_paulis"]
+            param_index = meta["param_index"]
+            expvals, variances = _pauli_expval_with_variance(count, meas_paulis)
+            for pauli, expval, variance in zip(orig_paulis, expvals, variances):
+                expval_map[param_index, pauli.to_label()] = (expval, variance)
+        return expval_map
+
+    def _preprocessing(
+        self, circuit: QuantumCircuit, observable: PauliList, param_index: tuple[int, ...]
+    ) -> list[QuantumCircuit]:
+        # generate measurement circuits with metadata
+        meas_circuits: list[QuantumCircuit] = []
+        if self._abelian_grouping:
+            for obs in observable.group_commuting(qubit_wise=True):
+                basis = Pauli((np.logical_or.reduce(obs.z), np.logical_or.reduce(obs.x)))
+                meas_circuit, indices = _measurement_circuit(circuit.num_qubits, basis)
+                paulis = PauliList.from_symplectic(
+                    obs.z[:, indices],
+                    obs.x[:, indices],
+                    obs.phase,
+                )
+                meas_circuit.metadata = {
+                    "orig_paulis": obs,
+                    "meas_paulis": paulis,
+                    "param_index": param_index,
+                }
+                meas_circuits.append(meas_circuit)
+        else:
+            for basis in observable:
+                meas_circuit, indices = _measurement_circuit(circuit.num_qubits, basis)
+                obs = PauliList(basis)
+                paulis = PauliList.from_symplectic(
+                    obs.z[:, indices],
+                    obs.x[:, indices],
+                    obs.phase,
+                )
+                meas_circuit.metadata = {
+                    "orig_paulis": obs,
+                    "meas_paulis": paulis,
+                    "param_index": param_index,
+                }
+                meas_circuits.append(meas_circuit)
+
+        # unroll basis gates
+        meas_circuits = self._passmanager.run(meas_circuits)
+
+        # combine measurement circuits
+        preprocessed_circuits = []
+        for meas_circuit in meas_circuits:
+            circuit_copy = circuit.copy()
+            # meas_circuit is supposed to have a classical register whose name is different from
+            # those of the transpiled_circuit
+            clbits = meas_circuit.cregs[0]
+            for creg in circuit_copy.cregs:
+                if clbits.name == creg.name:
+                    raise QiskitError(
+                        "Classical register for measurements conflict with those of the input "
+                        f"circuit: {clbits}. "
+                        "Recommended to avoid register names starting with '__'."
+                    )
+            circuit_copy.add_register(clbits)
+            circuit_copy.compose(meas_circuit, clbits=clbits, inplace=True)
+            circuit_copy.metadata = meas_circuit.metadata
+            preprocessed_circuits.append(circuit_copy)
+        return preprocessed_circuits
+
+
+def _measurement_circuit(num_qubits: int, pauli: Pauli):
+    # Note: if pauli is I for all qubits, this function generates a circuit to measure only
+    # the first qubit.
+    # Although such an operator can be optimized out by interpreting it as a constant (1),
+    # this optimization requires changes in various methods. So it is left as future work.
+    qubit_indices = np.arange(pauli.num_qubits)[pauli.z | pauli.x]
+    if not np.any(qubit_indices):
+        qubit_indices = [0]
+    meas_circuit = QuantumCircuit(
+        QuantumRegister(num_qubits, "q"), ClassicalRegister(len(qubit_indices), f"__c_{pauli}")
+    )
+    for clbit, i in enumerate(qubit_indices):
+        if pauli.x[i]:
+            if pauli.z[i]:
+                meas_circuit.sdg(i)
+            meas_circuit.h(i)
+        meas_circuit.measure(i, clbit)
+    return meas_circuit, qubit_indices

releasenotes/notes/add-backend-estimator-v2-26cf14a3612bb81a.yaml

@@ -0,0 +1,4 @@
+---
+features:
+  - |
+    The implementation :class:`~.BackendEstimatorV2` of :class:`~.BaseEstimatorV2` was added.

test/python/primitives/test_backend_estimator.py

@@ -101,7 +101,7 @@ def test_estimator_run(self, backend):
         # Note that passing objects has an overhead
         # since the corresponding indices need to be searched.
         # User can append a circuit and observable.
-        # calculate [ <psi2(theta2)|H2|psi2(theta2)> ]
+        # calculate [ <psi2(theta2)|H1|psi2(theta2)> ]
         result2 = estimator.run([psi2], [hamiltonian1], [theta2]).result()
         np.testing.assert_allclose(result2.values, [2.97797666], rtol=0.5, atol=0.2)