Add uncorrelated readout error mitigation for arbitrarily long Pauli …

…strings (quantumlib#6654)

* Add uncorrelated readout error mitigation for arbitrarily long pauli strings

cirq-core/cirq/experiments/readout_confusion_matrix.py

@@ -26,12 +26,26 @@
     import cirq
 
 
+def _mitigate_single_bitstring(z_rinv_all: list[np.ndarray], bitstring: np.ndarray) -> float:
+    """Return the mitigated Pauli expectation value for a single observed bitstring.
+
+    Args:
+        z_rinv_all: A list of single-qubit pauli-Z (as a vector) contracted with the single-qubit
+            inverse response matrix, for each qubit.
+        bitstring: The measured bitstring.
+
+    Returns:
+        The corrected expectation value of ZZZZZ... given the single measured bitstring.
+    """
+    return np.prod([z_rinv[bit] for z_rinv, bit in zip(z_rinv_all, bitstring)])
+
+
 class TensoredConfusionMatrices:
     """Store and use confusion matrices for readout error mitigation on sets of qubits.
 
     The confusion matrix (CM) for one qubit is:
 
-        [ Pr(0|0) Pr(1|0) ]
+        [ Pr(0|0) Pr(0|1) ]
         [ Pr(1|0) Pr(1|1) ]
 
     where Pr(i | j) = Probability of observing state "i" given state "j" was prepared.
@@ -71,7 +85,7 @@ def __init__(
             confusion_matrices: Sequence of confusion matrices, computed for qubit patterns present
                                 in `measure_qubits`. A single confusion matrix is also accepted.
             measure_qubits: Sequence of smaller qubit patterns, for which the confusion matrices
-                            were computed. A single qubit pattern is also accepted. Note that the
+                            were computed. A single qubit pattern is also accepted. Note that
                             each qubit pattern is a sequence of qubits used to label the axes of
                             the corresponding confusion matrix.
             repetitions:    The number of repetitions that were used to estimate the confusion
@@ -298,6 +312,72 @@ def func(x):
             )  # pragma: no cover
         return res.x
 
+    def readout_mitigation_pauli_uncorrelated(
+        self, qubits: Sequence['cirq.Qid'], measured_bitstrings: np.ndarray
+    ) -> tuple[float, float]:
+        r"""Uncorrelated readout error mitigation for a multi-qubit Pauli operator.
+
+        This function scalably performs readout error mitigation on an arbitrary-length Pauli
+        operator. It is a reimplementation of https://github.com/eliottrosenberg/correlated_SPAM
+        but specialized to the case in which readout is uncorrelated. We require that the confusion
+        matrix is a tensor product of single-qubit confusion matrices. We then invert the confusion
+        matrix by inverting each of the $C^{(q)}$ Then, in a bit-by-bit fashion, we apply the
+        inverses of the single-site confusion matrices to the bits of the measured bitstring,
+        contract them with the single-site Pauli operator, and take the product over all of the
+        bits. This could be generalized to tensor product spaces that are larger than single qubits,
+        but the essential simplification is that each tensor product space is small, so that none of
+        the response matrices is exponentially large.
+
+        This can result in mitigated Pauli operators that are not in the range [-1, 1], but if
+        the readout error is indeed uncorrelated and well-characterized, then it should converge
+        to being within this range. Results are improved both by a more precise characterization
+        of the response matrices (whose statistical uncertainty is not accounted for in the error
+        propagation here) and by increasing the number of measured bitstrings.
+
+        Args:
+            qubits: The qubits on which the Pauli operator acts.
+            measured_bitstrings: The experimentally measured bitstrings in the eigenbasis of the
+                Pauli operator. measured_bitstrings[i,j] is the ith bitstring, qubit j.
+
+        Returns:
+            The error-mitigated expectation value of the Pauli operator and its statistical
+            uncertainty (not including the uncertainty in the confusion matrices for now).
+
+        Raises:
+            NotImplementedError: If the confusion matrix is not a tensor product of single-qubit
+                                 confusion matrices for all of `qubits`.
+        """
+
+        # first, get all of the confusion matrices
+        cm_all = []
+        for qubit in qubits:
+            try:
+                idx = self.measure_qubits.index((qubit,))
+            except:  # pragma: no cover
+                raise NotImplementedError(  # pragma: no cover
+                    "The response matrix must be a tensor product of single-qu"  # pragma: no cover
+                    + f"bit response matrices, including that of qubit {qubit}."  # pragma: no cover
+                )  # pragma: no cover
+            cm_all.append(self.confusion_matrices[idx])
+
+        # get the correction matrices, assuming uncorrelated readout:
+        cminv_all = [np.linalg.inv(cm) for cm in cm_all]
+
+        # next, contract them with the single-qubit Pauli operators:
+        z = np.array([1, -1])
+        z_cminv_all = [z @ cminv for cminv in cminv_all]
+
+        # finally, mitigate each bitstring:
+        z_mit_all_shots = np.array(
+            [
+                _mitigate_single_bitstring(z_cminv_all, bitstring)
+                for bitstring in measured_bitstrings
+            ]
+        )
+
+        # return mean and statistical uncertainty:
+        return np.mean(z_mit_all_shots), np.std(z_mit_all_shots) / np.sqrt(len(measured_bitstrings))
+
     def __repr__(self) -> str:
         return (
             f"cirq.TensoredConfusionMatrices("

cirq-core/cirq/experiments/readout_confusion_matrix_test.py

@@ -17,6 +17,34 @@
 import pytest
 
 from cirq.experiments.single_qubit_readout_calibration_test import NoisySingleQubitReadoutSampler
+from cirq.experiments.readout_confusion_matrix import TensoredConfusionMatrices
+
+
+def add_readout_error(
+    measurements: np.ndarray,
+    zero_errors: np.ndarray,
+    one_errors: np.ndarray,
+    rng: np.random.Generator,
+) -> np.ndarray:
+    """Add readout errors to measured (or simulated) bitstrings.
+
+    Args:
+        measurements: The bitstrings to which we will add readout errors. measurements[i,j] is the
+                      ith bitstring, qubit j.
+        zero_errors: zero_errors[i] is the probability of a 0->1 readout error on qubit i.
+        one_errors: one_errors[i] is the probability of a 1->0 readout error on qubit i.
+        rng: The pseudorandom number generator to use.
+
+    Returns:
+        New measurements but with readout errors added.
+    """
+    num_bitstrs, n = measurements.shape
+    assert len(zero_errors) == len(one_errors) == n
+    # compute the probability that each bit is 1 after adding readout errors:
+    p1 = measurements * (1 - one_errors) + (1 - measurements) * zero_errors
+    r = rng.random((num_bitstrs, n))
+    noisy_measurements = r < p1
+    return noisy_measurements.astype(int)
 
 
 def get_expected_cm(num_qubits: int, p0: float, p1: float):
@@ -159,3 +187,78 @@ def test_readout_confusion_matrix_repr_and_equality():
     eq.add_equality_group(a, a)
     eq.add_equality_group(b, b)
     eq.add_equality_group(c, c)
+
+
+def _sample_ghz(n: int, repetitions: int, rng: np.random.Generator) -> np.ndarray:
+    """Sample a GHZ state in the z basis.
+    Args:
+        n: The number of qubits.
+        repetitions: The number of repetitions.
+        rng: The pseudorandom number generator to use.
+    Returns:
+        An array of the measurement outcomes.
+    """
+    return np.tile(rng.integers(0, 2, size=repetitions), (n, 1)).T
+
+
+def _add_noise_and_mitigate_ghz(
+    n: int,
+    repetitions: int,
+    zero_errors: np.ndarray,
+    one_errors: np.ndarray,
+    rng: np.random.Generator | None = None,
+) -> tuple[float, float, float, float]:
+    """Add readout error to GHZ-like bitstrings and measure <ZZZ...> with and
+    without readout error mitigation.
+    Args:
+        n: The number of qubits.
+        repetitions: The number of repetitions.
+        zero_errors: zero_errors[i] is the probability of a 0->1 readout error on qubit i.
+        one_errors: one_errors[i] is the probability of a 1->0 readout error on qubit i.
+        rng: The pseudorandom number generator to use.
+    Returns:
+        A tuple of:
+        - The mitigated expectation value of <ZZZ...>
+        - The statistical uncertainty of the previous output
+        - The unmitigated expectation value of <ZZZ...>
+        - The statstical uncertainty of the previous output
+    """
+    if rng is None:
+        rng = np.random.default_rng(0)
+    confusion_matrices = [
+        np.array([[1 - e0, e1], [e0, 1 - e1]]) for e0, e1 in zip(zero_errors, one_errors)
+    ]
+    qubits = cirq.LineQubit.range(n)
+    tcm = TensoredConfusionMatrices(
+        confusion_matrices, [[q] for q in qubits], repetitions=0, timestamp=0.0
+    )
+
+    measurements = _sample_ghz(n, repetitions, rng)
+    noisy_measurements = add_readout_error(measurements, zero_errors, one_errors, rng)
+    # unmitigated:
+    p1 = np.mean(np.sum(noisy_measurements, axis=1) % 2)
+    z = 1 - 2 * np.mean(p1)
+    dz = 2 * np.sqrt(p1 * (1 - p1) / repetitions)
+    # return mitigated and unmitigated:
+    return (*tcm.readout_mitigation_pauli_uncorrelated(qubits, noisy_measurements), z, dz)
+
+
+def test_uncorrelated_readout_mitigation_pauli():
+    n_all = np.arange(2, 35)
+    z_all_mit = []
+    dz_all_mit = []
+    z_all_raw = []
+    dz_all_raw = []
+    repetitions = 10_000
+    for n in n_all:
+        e0 = np.ones(n) * 0.005
+        e1 = np.ones(n) * 0.03
+        z_mit, dz_mit, z_raw, dz_raw = _add_noise_and_mitigate_ghz(n, repetitions, e0, e1)
+        z_all_mit.append(z_mit)
+        dz_all_mit.append(dz_mit)
+        z_all_raw.append(z_raw)
+        dz_all_raw.append(dz_raw)
+
+    for n, z, dz in zip(n_all, z_all_mit, dz_all_mit):
+        ideal = 1.0 if n % 2 == 0 else 0.0
+        assert np.isclose(z, ideal, atol=4 * dz)