diff --git a/cirq-core/cirq/experiments/readout_confusion_matrix.py b/cirq-core/cirq/experiments/readout_confusion_matrix.py index d99925aad14..a11efee1f03 100644 --- a/cirq-core/cirq/experiments/readout_confusion_matrix.py +++ b/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(" diff --git a/cirq-core/cirq/experiments/readout_confusion_matrix_test.py b/cirq-core/cirq/experiments/readout_confusion_matrix_test.py index bfd13e30cb6..48c0ce889d1 100644 --- a/cirq-core/cirq/experiments/readout_confusion_matrix_test.py +++ b/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 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 + - The statistical uncertainty of the previous output + - The unmitigated expectation value of + - 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)