Error mitigation

This lesson provides background on various error mitigation techniques to help users develop intuition for how and when they are best used. Importantly, modern applications of these techniques will leverage more sophisticated tools, like Samplomatic, Shaded Lightcones, and Executor. Those workflows will be addressed in later lessons.

This video gives an overview of error mitigation techniques. The text below adds more detail and explores these methods using Qiskit.

Overview of error mitigation techniques

Error mitigation refers to any process that attempts to reduce the effects of errors in quantum computations after they have occurred. It is different from suppression, which attempts to prevent errors, and from error correction, which encodes information across many physical qubits and checks the consistency in that information to preserve information in a single logical qubit. Throughout this lesson, error mitigation techniques can best be understood as methods for constructing better estimators of ideal expectation values or better samplers of probability distributions, not as methods for correcting quantum states during execution. We will examine the error mitigation methods below, all of which are available in Qiskit and its addons, and which are promising tools in the pursuit of quantum advantage.

Twirled readout error extinction (TREX) randomizes the measurement process by applying randomly selected Pauli operations before measurement and then classically undoing their known effects on the observed bitstrings. Under suitable measurement-noise assumptions, this converts the readout error channel into a simpler and more symmetric form that is easier to calibrate and mitigate.

Matrix-free measurement mitigation (M3) is a scalable measurement-error mitigation technique that avoids explicitly constructing the assignment matrix or its inverse.

Zero noise extrapolation (ZNE) approximates the value an observable would have in the absence of noise by artificially scaling noise by a known factor and extrapolating back to the zero-noise case. This requires that noise scale in a predictable way. That predictability might require other methods discussed in this lesson (pulse stretching, gate folding, and probabilistic error amplification (PEA)).

Probabilistic error cancellation (PEC) is among the most powerful, but also the most computationally expensive, error-mitigation techniques. It works by learning the noise behavior of the system and attempting to undo the noise by statistically constructing an unbiased estimate of the ideal computation. PEC outputs statistically corrected observables, not error-corrected results.

We will now dive into each of these methods.

M3

Measurement errors are one of the most persistent and hardware-specific sources of noise in near-term quantum devices. Even when a quantum state is prepared and evolved correctly, the final readout step can misreport the true computational basis outcome. These errors arise from imperfect discrimination between the physical states corresponding to $|0\rangle$ and $|1\rangle$, relaxation during measurement, and imperfections in measurement electronics. As a result, the observed distribution of bitstrings is often a distorted version of the true (ideal) distribution, with probability incorrectly reassigned between outcomes.

Measurement mitigation addresses this problem by explicitly modeling how measurement outcomes are misassigned. In the standard formulation, this behavior is captured by an assignment matrix $A$, whose elements $A_{ij}$ represent the probability of recording outcome $i$ when the true basis state was $j$. If the measurement noise is classical and stationary, the observed probability vector $p_{\text{obs}}$ is related to the ideal distribution $p_{\text{ideal}}$ by a simple linear relation,

Note that the assignment matrix $A$ is a classical stochastic map, not a quantum channel. It acts on classical readouts, not on the underlying quantum states that were sampled to obtain the readouts.

Matrix-free measurement mitigation (M3) uses this structure to construct a mitigated estimate of $p_{\text{ideal}}$ without ever forming or inverting the full matrix $A$, which would be infeasible for large numbers of qubits. Instead, M3 exploits the fact that many near-term quantum workloads produce highly peaked or low-entropy distributions, so that mitigation focuses only on the subset of outcomes that significantly contribute to the observed distribution.

To build the assignment model, additional calibration experiments must be run on the quantum hardware. These typically involve preparing known computational basis states and recording how they are reported by the measurement system. The calibration overhead scales with the number of measured qubits, while the computational cost of mitigation depends on the structure of the observed distribution. This cost is far smaller than what would be required to characterize the full $2^n \times 2^n$ assignment matrix explicitly. While this calibration overhead is real, it is lightweight compared to circuit-level techniques like zero-noise extrapolation or probabilistic error cancellation, making M3 a practical and widely used tool for mitigating readout errors in near-term quantum experiments.

Note that calibration is only useful to the extent that the measurement process remains sufficiently stable between the calibration runs and the target experiment.

We will first carry out a simple case by hand to illustrate how it works, then use tools built into Qiskit to show the use of M3 at a more realistic scale.

M3 small example

We start by building a simple Bell state.

from qiskit import QuantumCircuit
import numpy as np

qc = QuantumCircuit(2, 2)
qc.h(0)
qc.cx(0, 1)
qc.measure([0, 1], [0, 1])

Output:

<qiskit.circuit.instructionset.InstructionSet at 0x714a5fa4d630>

We can use a simulator with a noise model because we are restricting this example to a very small circuit. We explicitly configure our noise model to have readout error of the form:

if true=0: [95% chance of reading 0, 5% chance of reading 1]

if true=1: [8% chance of reading 0, 92% chance of reading 1]

These readout error rates are much higher than what one sees in real IBM® quantum computers. The high values have been chosen for pedagogical reasons, so that we can clearly see the effect of M3.

from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, ReadoutError

ro_error = ReadoutError([[0.95, 0.05], [0.08, 0.92]])

noise_model = NoiseModel()
noise_model.add_readout_error(ro_error, [0])
noise_model.add_readout_error(ro_error, [1])

backend_noisy = AerSimulator(noise_model=noise_model)

We can run the circuit on the noisy simulator to observe the effects of readout error.

result = backend_noisy.run(qc, shots=5000, seed_simulator=123).result()
raw_counts = result.get_counts()
raw_counts

Output:

{'01': 290, '11': 2165, '10': 304, '00': 2241}

In our toy example, we defined the readout error for our noise model. But when using a real quantum computer, the readout error is not known a priori. Instead, it must be estimated experimentally. For this, we create circuits to prepare a qubit in a known computational basis state $|0\rangle$ or $|1\rangle$ and measure it. In the absence of readout error we would measure exactly the prepared state. By collecting statistics from these measurements, we can build an assignment matrix describing the readout error. In this toy example, it should look similar to the one we defined.

For a more realistic system, statistics would need to be accumulated for each physical qubit. This is sometimes non-trivial, because transpilation maps an abstract circuit to hardware and the qubit number changes. We will revisit that in the next example. In this simulator example, we assigned the same readout-error model to every qubit, so a single characterization is sufficient.

# Define a function to make a state and measure it.
def measure_state(state):
    qc = QuantumCircuit(1, 1)
    if state == 1:
        qc.x(0)
    qc.measure(0, 0)
    return backend_noisy.run(qc, shots=5000, seed_simulator=123).result().get_counts()

# Use the function to get statistics on readout errors
cal_counts_0 = measure_state(0)  # True state = 0
cal_counts_1 = measure_state(1)  # True state = 1

cal_counts_0, cal_counts_1

Output:

({'1': 244, '0': 4756}, {'0': 420, '1': 4580})

# Create the normalized readout error matrix

# def normalize(d):
#    s = sum(d.values())
#    return {k: v/s for k,v in d.items()}
# p0 = normalize(cal_counts_0)
# p1 = normalize(cal_counts_1)

p0 = {k: v / sum(cal_counts_0.values()) for k, v in cal_counts_0.items()}
p1 = {k: v / sum(cal_counts_1.values()) for k, v in cal_counts_1.items()}


A_1q = np.array(
    [
        [p0["0"], p1["0"]],  # row: measured 0, columns true 0/1
        [p0["1"], p1["1"]],
    ]
)
A_1q

Output:

array([[0.9512, 0.084 ],
       [0.0488, 0.916 ]])

Note that this is very close to the readout error matrix used to construct the simulator.

Again, on a real quantum computer we would have statistics like these for multiple physical qubits. For this simulator, we assume the two qubits have identical and independent readout errors, allowing us to construct the two-qubit assignment matrix A_2q as a tensor product of the corresponding one-qubit matrices. Also note that the point of M3 is not to construct the entire assignment matrix. Instead, M3 focuses on the subset of outcomes that significantly contribute to the observed distribution. We will build the entire matrix for our very small example to see how it works. The implementation of M3 in qiskit-ibm-runtime is more sophisticated than this toy example because it must handle many qubits with distinct readout characteristics while still avoiding construction of the full assignment matrix.

A_2q = np.kron(A_1q, A_1q)
print(A_2q)

Output:

[[0.89794576 0.07107    0.07107    0.005625  ]
 [0.04965424 0.87653    0.00393    0.069375  ]
 [0.04965424 0.00393    0.87653    0.069375  ]
 [0.00274576 0.04847    0.04847    0.855625  ]]

For this example we assume that readout errors on different qubits are independent, allowing the two-qubit assignment matrix to be constructed as a tensor product of the corresponding one-qubit matrices.

We obtain the raw and mitigated results for comparison.

bitstrings = ["00", "01", "10", "11"]
# Normalize raw counts
raw_probs = np.array([raw_counts.get(b, 0) for b in bitstrings], dtype=float)
raw_probs /= raw_probs.sum()
raw_probs

Output:

array([0.4482, 0.0592, 0.0598, 0.4328])

At this point, we have experimentally characterized how the measurement process distorts ideal outcome probabilities. Our goal is now to undo this distortion and apply the resulting correction to the raw probabilities.

Because we are performing a linear inversion rather than a physical correction, the resulting "probabilities" need not always be positive. Very small negative quasi-probabilities are expected and reflect the fact that measurement mitigation produces an estimator rather than a physical probability distribution.

# Explicit inversion is feasible only for very small systems
# This is why M3 avoids creating the entire A matrix.

A_inv = np.linalg.inv(A_2q)

# Note: mitigated_probs might contain negative values.
# This is expected and reflects estimator, not physical, probabilities.

mitigated_probs = A_inv @ raw_probs
print(mitigated_probs)

Output:

[ 4.96000451e-01 -4.69852308e-04  2.17747967e-04  5.04251654e-01]

We convert these to counts for plotting.

shots = sum(raw_counts.values())  # e.g., 5000
ideal_probs = [0.50, 0.00, 0.00, 0.50]  # ideal Bell distribution

# Convert probabilities to integer counts for plotting only.
# Not physically meaningful when probabilities are negative


def probs_to_counts(bitstrings, probs, shots):
    """
    Convert a probability vector into a Qiskit-style counts dictionary.
    Ensures integer shot counts and correct normalization.
    """
    counts = {}
    # First convert to raw integer counts
    int_counts = [round(p * shots) for p in probs]
    # Fix rounding error by adjusting the last entry
    diff = shots - sum(int_counts)
    int_counts[-1] += diff
    # Build dict
    for b, c in zip(bitstrings, int_counts):
        if c > 0:
            counts[b] = c
    return counts


mitigated_counts = probs_to_counts(bitstrings, mitigated_probs, shots)
ideal_counts = probs_to_counts(bitstrings, ideal_probs, shots)

Finally we look at the raw, mitigated, and ideal distributions. The mitigated "counts" below are only a visualization aid obtained by rescaling the mitigated quasi-probabilities.

from qiskit.visualization import plot_histogram

plot_histogram(
    [raw_counts, mitigated_counts, ideal_counts],
    legend=["Raw", "Mitigated", "Ideal"],
    sort="asc",
    title="Raw vs Mitigated vs Ideal",
)

Output:

Clearly, the M3-mitigated results are much closer to the ideal distributions, and we can see why.

Applying M3 to more realistic workflows is much more complex. Consider:

• Real quantum computers contain many physical qubits, each with its own readout-error characteristics. These errors must be modeled and combined appropriately.
• The full assignment matrix $A$ grows exponentially with the number of qubits and quickly becomes impractical to construct explicitly. This is why M3 focuses on the subset of outcomes that contribute significantly to the observed distribution rather than constructing the entire assignment matrix.
• For real problems of interest, the ideal distribution might not be known.

For all of these reasons, the full mitigation workflow is encapsulated in the mthree Qiskit addon.

Larger M3 example

We now apply M3 to a larger circuit, one for which determining the relevant subspace of $A$ by hand would be cumbersome, but which is still simple enough that we know the ideal distribution: a GHZ state. We start importing the mthree package and defining a GHZ quantum circuit.

import numpy as np
from qiskit import QuantumCircuit

# from qiskit_ibm_runtime.fake_provider import FakeCasablancaV2
import mthree

num_qubits = 10

qc = QuantumCircuit(num_qubits)
qc.reset(range(num_qubits))
qc.h(0)
for i in range(num_qubits - 1):
    qc.cx(i, i + 1)
qc.measure_all()
qc.draw("mpl")

Output:

To run this on a real quantum computer, we import qiskit_ibm_runtime and other necessary packages. We use SamplerV2 to obtain a distribution of outputs.

# Import packages needed for using real quantum computers

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime.executor_sampler import SamplerV2 as Sampler
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager


# --------- Configuration ----------
service = QiskitRuntimeService(channel="ibm_quantum_platform", instance="your crn here")

# backend_real = service.backend("ibm_fez")     # adjust if needed
backend_real = service.least_busy(simulator=False, operational=True)

Now we transpile the circuit to our backend. We are setting the optimization level to 0 here in order to focus on the effects of M3, without adding in additional tools. This means the final fidelity will not be representative of what is achievable on IBM Quantum® hardware using the full suite of available tools, but it allows us to focus specifically on the usefulness of M3.

# Exclude any optimization to focus only on M3 effects
opt_level = 0
pm = generate_preset_pass_manager(optimization_level=opt_level, backend=backend_real)
qc_isa = pm.run(qc)

We now configure the Sampler, again disabling additional error-suppression techniques such as dynamical decoupling and twirling.

shots = 8192

# No twirling
sampler_no_twirl = Sampler(mode=backend_real)
sampler_no_twirl.options.twirling.enable_gates = False
sampler_no_twirl.options.twirling.enable_measure = False
sampler_no_twirl.options.default_shots = shots
sampler_no_twirl.options.dynamical_decoupling.enable = False

Finally, we execute our job.

# This job took 7 seconds on ibm_fez. Your actual execution time can vary.
job = sampler_no_twirl.run([qc_isa])
res_real = job.result()
job_id = job.job_id()
print("job number for real  ", job_id)

Output:

job number for real   d7ark4e5nvhs73a2fu6g

We can extract the raw counts from the circuit above.

counts_real = []
for i in range(len(res_real)):
    counts_real.append(res_real[i].data.meas.get_counts())

Now we run mthree on our backend to characterize the assignment model used for measurement mitigation. We rescale the raw distribution using that information. Note that this can sometimes result in very small negative probabilities (hence the term quasi-probabilities).

# Added 6 sec on ibm_fez, on top of original 7 sec to run circuits. Your time might differ.
mit_real = mthree.M3Mitigation(backend_real)
mappings_real = mthree.utils.final_measurement_mapping(qc_isa)
mit_real.cals_from_system(mappings_real)
quasis_real = mit_real.apply_correction(counts_real, mappings_real)

The mthree correction produces quasi-probabilities, rather than counts. We can rescale to counts by hand. The helper function below converts to counts and rounds to integer values. Also, for large circuits, there will be many erroneous strings with very few counts, but which clutter a histogram. The second helper function below trims away any strings with fewer counts than some threshold we choose, to give us a clear view of the dominant terms in the distribution.

def quasidistribution_to_counts(quasi_list, shots):
    return [
        {bitstr: int(round(float(prob) * shots)) for bitstr, prob in quasi.items()}
        for quasi in quasi_list
    ]


def filter_counts_by_threshold(counts_list, t):
    filtered = []

    for counts in counts_list:
        new_counts = {}
        for bitstr, val in counts.items():
            c = int(round(float(val)))  # FORCE conversion
            if abs(c) >= t:
                new_counts[bitstr] = c
        filtered.append(new_counts)

    return filtered

quasi_counts = quasidistribution_to_counts(quasis_real, shots)
trim_threshold = 30
trimmed_quasi = filter_counts_by_threshold(quasi_counts, trim_threshold)
trimmed_raw = filter_counts_by_threshold(counts_real, trim_threshold)

We can also easily produce the theoretical, ideal distribution in which we have equal probabilities for the all-zero and all-one states, and zero probability for all other outcomes, for comparison.

ideal = [
    {
        "0" * num_qubits: int(round(float(shots) / 2)),
        "1" * num_qubits: int(round(float(shots) / 2)),
    }
]

Finally, we can visualize our results.

from qiskit.visualization import plot_histogram

plot_histogram(
    [trimmed_raw[0], trimmed_quasi[0], ideal[0]],
    legend=["Raw", "Mitigated", "Ideal"],
    sort="asc",
    title="Raw vs Mitigated vs Ideal",
)

Output:

The mitigated results are closer to the ideal distribution than the unmitigated results. We can quantify this better using counts for the ideal bitstrings.

def extract_all_zero_and_all_one(counts):
    """
    Extract counts for the all-0 and all-1 bitstrings.

    Parameters
    ----------
    counts : dict[str, int]
        Counts dictionary.

    Returns
    -------
    dict[str, int]
        Dictionary containing only 00...0 and 11...1 entries if present.
    """
    n = len(next(iter(counts)))  # infer bitstring length
    zero_state = "0" * n
    one_state = "1" * n

    return {k: counts[k] for k in (zero_state, one_state) if k in counts}


correct_counts_raw = extract_all_zero_and_all_one(trimmed_raw[0])
correct_counts_mitigated = extract_all_zero_and_all_one(trimmed_quasi[0])
keys = ["0" * num_qubits, "1" * num_qubits]
fidelity_raw = (
    correct_counts_raw.get(keys[0]) + correct_counts_raw.get(keys[1])
) / shots
fidelity_mit = (
    correct_counts_mitigated.get(keys[0]) + correct_counts_mitigated.get(keys[1])
) / shots
print("Percent correct raw counts: ", fidelity_raw)
print("Percent correct mitigated counts: ", fidelity_mit)

Output:

Percent correct raw counts:  0.6654052734375
Percent correct mitigated counts:  0.838134765625

This simple tool improved our fidelity by 26% relative to the unmitigated fidelity, with fairly low overhead. Clearly M3 can be very valuable. But the GHZ circuit was selected as a good use case, because it is highly peaked on two output states. Some problems result in more evenly-distributed counts (high entropy states). For very deep circuits with many qubits, gate errors and noise unrelated to measurement errors will often result in a high-entropy state. In such cases, M3 operates on a state space that is quite large. There are three possible approaches to this:

(1) Accept that this is a limitation of M3.

(2) Note that the upper bound on the space of the matrix formed by M3 is num_shots $\times$ num_shots, rather than $2^N \times 2^N$, with the former being much smaller for large number of qubits $N$.

(3) Only with great caution, remove very low-probability bitstrings from the raw counts fed into M3. This might introduce bias, and the limitation of the support of your states can interfere with convergence in approaches like the Krylov method. This should only be done when you know exactly why you are doing it and can anticipate all its effects. That is why such trimming is not a built-in feature of M3.

For more on the mthree Qiskit addon see examples, installation instructions, background, and API references in its documentation.

Twirled readout error extinction (TREX)

Twirled Readout Error Extinction (TREX) is a measurement-error mitigation technique used on modern quantum computers. Its goal is to remove bias introduced by classical readout errors without requiring full characterization or inversion of the measurement assignment matrix.

TREX is not error correction and does not improve state fidelity. Instead, it ensures that expectation value estimators are unbiased with respect to readout noise, up to statistical error, by averaging over randomized measurement bases and rescaling. On IBM Quantum hardware, TREX is the default mitigation method for expectation-value estimation and corresponds to resilience_level = 1 in the IBM Quantum Estimator primitive.

Nature of readout errors on superconducting hardware

On IBM superconducting qubits, readout errors are dominated by:

• Asymmetric bit-flip probabilities (for example, $|1\rangle \rightarrow |0\rangle$ more frequently than the reverse)
• State-dependent relaxation during measurement itself
• Mild but non-negligible crosstalk between neighboring readout resonators

Importantly, these errors can often be modeled as a classical stochastic process acting on measurement outcomes. TREX exploits a standard idea from noise tailoring: randomization (twirling) of the measurement procedure. TREX typically proceeds in two stages:

1. Twirling and classical "untwirling"
2. Rescaling of expectation values using calibration circuits

Suppose an expectation value of observable $O$ is related to noisy estimates, like the following:

where $a$ is a fixed factor that scales the estimate and $b$ is a bias from sources described above. We are using $\approx$ here as a reminder that we are only addressing readout error here; there is still gate error to consider.

Just before measurement, each qubit is randomly conjugated with a Pauli operator (like I or X). Measurement outcomes are classically untwirled afterward.

Averaging over these randomizations symmetrizes the readout channel. As a result, systematic asymmetries in the readout process are converted into a more symmetric error channel whose remaining effect can be treated through calibration and rescaling. In simple terms, if for example there was a bias toward the $|0\rangle$ state, we would be missing some measurements of $|1\rangle$. By twirling around the measurement, we randomly interchange which logical outcome is associated with the physical readout process, then undo that interchange classically afterward. Averaging over many such randomized measurements removes the systematic preference for one reported outcome over the other. This removes the dominant measurement bias under the assumptions of the TREX model, so that we have:

However, expectation values might still need to be rescaled. The second stage is using calibration circuits (circuits with known outcomes) to estimate the scaling factor $a$. This can then be used to obtain a greatly enhanced estimate.

Expectation values of Pauli observables become unbiased estimators of their ideal values. No explicit noise model is learned or inverted, which is why TREX is numerically stable and inexpensive. Again, gate error has not been addressed here.

In the estimation of an expectation value $\langle O\rangle$ of observable $O$, there can be variance due to finite sampling, there can be many types of incoherent and coherent errors, and there can also be a bias, $\delta$, due to readout error. This bias arises from sources like the asymmetric transitions we've discussed. In practice, this bias is typically well-modeled as a classical stochastic channel acting after the quantum evolution. TREX transforms the effective measurement channel so that the estimator of $\langle O\rangle$ is unbiased with respect to readout error. This is relatively easy to see, using the simplest example of twirling gates: X. If there is asymmetry in $|1\rangle \rightarrow |0\rangle$ compared to $|0\rangle \rightarrow |1\rangle$, but we randomly switch which state is being read (flipping back later, classically), then this obviously removes the systematic asymmetry (though finite sampling still means there might be slight stochastic asymmetry in measurement). In practice, twirling gates for TREX are not limited to I and X, but this suffices for an explanation.

TREX allows us to obtain an estimate of an expectation value that is free of bias due to readout error. There could still be a bias from some other source (like coherent over-rotation); TREX targets only readout. Because TREX replaces systematic bias with randomized symmetrization, it generally converts bias into additional stochastic variance, increasing the variance of the estimator at fixed sampling cost.

Check your understanding

Suppose we use estimator on a circuit that has prepared the state $|\psi\rangle = \frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ in the estimation of $\langle Z\rangle$. Suppose further that measurement errors cause outcomes corresponding to $|1\rangle$ to be reported as $0$ 10% of the time, while outcomes corresponding to $|0\rangle$ are reported as $1$ only 1% of the time. The asymmetry between these reporting errors introduces a measurement bias.

(a) In the absence of any such measurement errors, what would be the ideal expectation value of the observable $\langle Z \rangle$?

(b) With the asymmetric measurement or readout errors, what is the expectation value of the observable $\langle Z \rangle$?

(c) Consider a limited set of Pauli twirling gates: I and X. Assume I and X are applied randomly with 50% probability in multiple shots. The same measurement error rates still apply. Then the Pauli gate (I or X) is undone classically. Make a table of all possible outcomes and their probabilities due to initial superposition, asymmetric measurement errors, and randomly assigned Pauli gates.

(d) Use the table to determine the expectation value $\langle Z \rangle$ using the twirling strategy.

(e) Compare the expectation values: ideal, without twirling, and with twirling.

(f) What has been omitted from this treatment that is found in a full TREX application?

• Answer

  (a) You can simply use intuition and symmetry to say $\langle Z \rangle = 0$. To show explicitly:

  $$\begin{aligned} \langle Z \rangle & = \frac{1}{\sqrt{2}}(\langle0|+\langle1|) Z \frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)\\ & = \frac{1}{2}(\langle0|Z|0\rangle+\langle0|Z|1\rangle+\langle1|Z|0\rangle+\langle1|Z|1\rangle)\\ & = \frac{1}{2}((1)\langle0|0\rangle+(-1)\langle0|1\rangle+(1)\langle1|0\rangle+(-1)\langle1|1\rangle)\\ & = \frac{1}{2}(1+0+0+-1)\\ & = 0 \end{aligned}$$

  (b) The state could collapse to either $|0\rangle$ or $|1\rangle$ when measured, each with 50% chance. But of the 50% that collapse to the $|0\rangle$ state, 1% of those (or 0.5% of all cases) will be mis-reported as $|1\rangle$. Of the 50% that collapse to $|1\rangle$, 10% (or 5% of all cases) will be mis-reported as $|0\rangle$. This leaves us with:

  • $|0\rangle$ 54.5%
  • $|1\rangle$ 45.5% Ignoring all other types of possible errors, we would obtain an estimate of

  $$\langle Z \rangle = (0.545)(1)+(0.455)(-1) = 0.09$$

  (c)

  Pauli	Ideal measurement outcome	Outcome after Pauli	Readout error	Undo Pauli	Final reported outcome	Probability
  I	0	0	None	I	0	50% $\times$ 50% $\times$ 99% = 24.75%
  I	0	0	$0 \rightarrow 1$	I	1	50% $\times$ 50% $\times$ 1% = 0.25%
  I	1	1	None	I	1	50% $\times$ 50% $\times$ 90% = 22.5%
  I	1	1	$1 \rightarrow 0$	I	0	50% $\times$ 50% $\times$ 10% = 2.5%
  X	0	1	None	X	0	50% $\times$ 50% $\times$ 90% = 22.5%
  X	0	1	$1 \rightarrow 0$	X	1	50% $\times$ 50% $\times$ 10% = 2.5%
  X	1	0	None	X	1	50% $\times$ 50% $\times$ 99% = 24.75%
  X	1	0	$0 \rightarrow 1$	X	0	50% $\times$ 50% $\times$ 1% = 0.25%

  (d) From previous parts, we know the expectation values of Z in each of the final states. We need only take the weighted average:

  $$\begin{aligned} \langle Z \rangle & = (0.2475)(1)+(0.0025)(-1)+(0.225)(-1)+(0.025)(1)\\ & + (0.225)(1)+(0.025)(-1)+(0.2475)(-1)+(0.0025)(1)\\ & = 0 \end{aligned}$$

  (e) In this very simple case our application of twirling recovered exactly the ideal expectation value. That is not typical. What is typical is that the twirling undoes biases like those introduced through asymmetric transitions, and returns an expectation value closer to the ideal.

  (f) This is a toy model application. There could be several things missing. Notably, we have not run calibration circuits here to estimate a rescaling factor $a$. This is typically necessary.

Implementation details

TREX is automatically applied when using the EstimatorV2 primitive with resilience_level = 1 (the default). Note that it is not applied automatically to the Sampler primitive, because Sampler returns full probability distributions rather than expectation values.

When TREX is enabled, the runtime generates measurement-twirled circuit variants by inserting random single-qubit Paulis immediately before measurement. Classical post-processing reverses these Paulis. These circuit variants are random across different shots and also across different qubits within a shot. All circuit variants are then combined into a single estimator.

TREX is very scalable. It requires additional calibration circuits, such that the sampling overhead is typically a small constant factor (often 2-4x). But the mitigation is observable-specific, not distribution-preserving.

TREX is designed to integrate cleanly with many error suppression and error mitigation methods, including the following:

• Dynamical decoupling (DD) - acts during idle times, orthogonal to readout mitigation.
• Gate or Pauli twirling - often complementary; TREX already relies on twirling but targets only readout noise.
• Zero-Noise Extrapolation (ZNE) - ZNE targets gate noise, while TREX handles measurement bias. This combination is used at resilience level 2.
• Probabilistic error amplification and probabilistic error cancellation in Estimator workflows.

However, some methods cannot (or should not) be combined with TREX:

• Sampler-level probability correction: TREX does not return corrected bitstring probabilities.
• Post-hoc mitigation on already-mitigated TREX data: double correction can reintroduce bias.

Combining TREX with assignment-matrix-based readout mitigation like M3 is tricky. M3 assumes a fixed measurement channel. TREX deliberately randomizes measurement, invalidating assumption used by M3. In that sense, you cannot simply combine these techniques. However, if you apply measurement twirling during the calibration used in M3, you might be able to leverage both an assignment matrix-like approach and a twirling mitigation approach simultaneously (though this twirling approach is not quite the same as TREX). This is a recurring theme in combining mitigation and suppression techniques: blindly adding techniques fails; adding them carefully such that they are "aware" of each other can work and can be highly beneficial.

TREX only guarantees expectation values that are free of bias due to readout errors; it does not reduce variance, nor does it protect against other types of bias. It certainly does not correct all errors. However, it is often a substantial improvement, and it is achievable with much lower overhead than some other methods that require learning full noise models. TREX is often the lowest-cost and most robust readout mitigation available, making it an excellent default in experiments.

TREX with Qiskit

Here we create a simple two-qubit circuit, and use Estimator to estimate the expectation value of a very simple observable, ZX. The ideal expectation value is 1. We will compare the expectation values with and without TREX.

# Import packages needed for using real quantum computers

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2, Session
import numpy as np


# --------- Configuration ----------

backend = service.backend("ibm_fez")  # adjust if needed

Create the simple quantum circuit.

qc = QuantumCircuit(2)
qc.x(0)
qc.h(0)
qc.cx(0, 1)
qc.h(0)
qc.barrier()
qc.draw("mpl")

Output:

Transpile the circuit using a pass manager.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pass_manager = generate_preset_pass_manager(optimization_level=1, backend=backend)

qc_t = pass_manager.run(qc)

Apply the same layout from transpilation to the observable.

# Observable: X on qubit 0, Z on qubit 1, ideal expectation value is 1.0.
observable = SparsePauliOp("ZX")
observable_t = observable.apply_layout(qc_t.layout)

Run the experiment on a real quantum computer. First we proceed with reslience_level = 1; this is the default, but we have made it explicit for distinguishing the two cases.

# First, we run this on an ideal simulator to confirm the ideal expectation value.

from qiskit.primitives import StatevectorEstimator as Estimator

estimator = Estimator()

job = estimator.run(pubs=[(qc, observable)])

res = job.result()
val = res[0].data.evs


print("⟨ZX⟩ ideal:", val)

Output:

⟨Z₀⟩ ideal: 0.9999999999999996

# This took 27 sec on ibm_fez. Your runtime might differ.

with Session(backend=backend) as session:
    estimator = EstimatorV2(
        mode=session,
        options={
            "resilience_level": 1,  # TREX
            "default_precision": 0.005,
        },
    )

    job = estimator.run(pubs=[(qc_t, observable_t)])

    res = job.result()
    val = res[0].data.evs

print("⟨ZX⟩ with TREX:", val)

Output:

⟨Z₀⟩ with TREX: 0.9520029143897996

Now we set the resilience_level to 0 to run without TREX.

# This took 12 sec on ibm_fez. Your runtime might differ.
with Session(backend=backend) as session:
    estimator = EstimatorV2(
        mode=session,
        options={
            "resilience_level": 0,  # TREX off
            "default_precision": 0.005,
        },
    )

    job = estimator.run(pubs=[(qc_t, observable_t)])

    res = job.result()
    val = res[0].data.evs

print("⟨ZX⟩ without TREX:", val)

Output:

⟨Z₀⟩ without TREX: 0.8651

The expectation value is substantially closer to the ideal (1) with TREX, as expected.

Of course, this was a single expectation value on a single qubit. We now prepare a circuit using non-Clifford gates and many qubits to test whether the TREX improvement remains robust in a different setting.

import numpy as np

num_qubits = 20
qc = QuantumCircuit(num_qubits)
for i in range(num_qubits):
    qc.rx(np.pi * 7 / 8, i)
qc.barrier()

Output:

CircuitInstruction(operation=Instruction(name='barrier', num_qubits=20, num_clbits=0, params=[]), qubits=(<Qubit register=(20, "q"), index=0>, <Qubit register=(20, "q"), index=1>, <Qubit register=(20, "q"), index=2>, <Qubit register=(20, "q"), index=3>, <Qubit register=(20, "q"), index=4>, <Qubit register=(20, "q"), index=5>, <Qubit register=(20, "q"), index=6>, <Qubit register=(20, "q"), index=7>, <Qubit register=(20, "q"), index=8>, <Qubit register=(20, "q"), index=9>, <Qubit register=(20, "q"), index=10>, <Qubit register=(20, "q"), index=11>, <Qubit register=(20, "q"), index=12>, <Qubit register=(20, "q"), index=13>, <Qubit register=(20, "q"), index=14>, <Qubit register=(20, "q"), index=15>, <Qubit register=(20, "q"), index=16>, <Qubit register=(20, "q"), index=17>, <Qubit register=(20, "q"), index=18>, <Qubit register=(20, "q"), index=19>), clbits=())

Here, there is nothing special about the rotation angle $7\pi/8$. We simply chose any unequal superposition of $|0\rangle$ and $|1\rangle$ in order to probe effects of asymmetric erroneous transitions and how they are treated by TREX. Given this state, the ideal expectation value of Z on every qubit, in the absence of any errors (readout or otherwise) would be:

We therefore hope to obtain expectation values close to -0.924 on all qubits. We will see if TREX improves each expectation value individually and also whether it improves the average over all qubits.

Check your understanding

Verify the expectation value above through direct calculation.

• Answer
  $$\begin{aligned} \langle Z \rangle & = \langle \psi |Z|\psi \rangle\\ & = \langle 0 |RX^\dagger (7\pi/8) \quad Z \quad RX(7\pi/8)|0 \rangle\\ & = \begin{pmatrix}1 & 0\end{pmatrix}\begin{pmatrix}\cos\left(\frac{7\pi}{16}\right) & -i\sin\left(\frac{7\pi}{16}\right) \\ -i\sin\left(\frac{7\pi}{16}\right) & \cos\left(\frac{7\pi}{16}\right)\end{pmatrix} \begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix} \begin{pmatrix}\cos\left(\frac{7\pi}{16}\right) & -i\sin\left(\frac{7\pi}{16}\right) \\ -i\sin\left(\frac{7\pi}{16}\right) & \cos\left(\frac{7\pi}{16}\right)\end{pmatrix}\begin{pmatrix}1 \\ 0\end{pmatrix}\\ & = \begin{pmatrix}1 & 0\end{pmatrix}\begin{pmatrix}\cos\left(\frac{7\pi}{16}\right) & -i\sin\left(\frac{7\pi}{16}\right) \\ -i\sin\left(\frac{7\pi}{16}\right) & \cos\left(\frac{7\pi}{16}\right)\end{pmatrix} \begin{pmatrix}\cos\left(\frac{7\pi}{16}\right) \\ -i\sin\left(\frac{7\pi}{16}\right)\end{pmatrix}\\ & = \begin{pmatrix}1 & 0\end{pmatrix}\begin{pmatrix}\cos^2\left(\frac{7\pi}{16}\right)-\sin^2\left(\frac{7\pi}{16}\right) \\ -2\sin\left(\frac{7\pi}{16}\right)\cos\left(\frac{7\pi}{16}\right)\end{pmatrix}\\ & = \cos^2\left(\frac{7\pi}{16}\right)-\sin^2\left(\frac{7\pi}{16}\right)\\ & = -0.9239... \end{aligned}$$

We transpile using a pass manager.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pass_manager = generate_preset_pass_manager(optimization_level=1, backend=backend_real)

qc_t = pass_manager.run(qc)

Now we construct one observable per qubit so that we can estimate all twenty $\langle Z_i\rangle$ expectation values in a single workflow. That is, for the $i^\text{th}$ qubit, we want an observable with all identities $I$, except at the $i^\text{th}$ position, where we want a $Z$. Remember to apply your transpiled layout to your observables as well.

# Observables: Z on qubit i for all i.
observables = []
for i in range(num_qubits):
    obs_string = "I" * i + "Z" + "I" * (num_qubits - i - 1)
    obs = SparsePauliOp(obs_string)
    observables.append(obs)
observables_t = [obs.apply_layout(qc_t.layout) for obs in observables]

Now we can execute. The code cell below uses session mode; if you do not have access to sessions, you can replace it with batch mode to run the same workflow but with possible queue times between circuits.

# This took 20 sec on ibm_fez. Your runtime might differ.

with Session(backend=backend_real) as session:
    estimator = EstimatorV2(
        mode=session,
        options={
            "resilience_level": 0,  # TREX turned off
            "default_precision": 0.005,
        },
    )

    job = estimator.run(pubs=[(qc_t, observables_t)])

    res_no_trex = job.result()
    val_no_trex = res_no_trex[0].data.evs

print("⟨Z_i⟩ without TREX:", val_no_trex)

Output:

⟨Z_i⟩ without TREX: [-0.90645 -0.89755 -0.89125 -0.9165  -0.87655 -0.91035 -0.8757  -0.87905
 -0.66035 -0.8876  -0.87125 -0.89975 -0.83745 -0.89755 -0.82645 -0.90325
 -0.84385 -0.90195 -0.85865 -0.8884 ]

# This took 36 sec on ibm_fez. Your runtime might differ.

with Session(backend=backend_real) as session:
    estimator = EstimatorV2(
        mode=session,
        options={
            "resilience_level": 1,  # TREX
            "default_precision": 0.005,
        },
    )

    job = estimator.run(pubs=[(qc_t, observables_t)])

    res_trex = job.result()
    val_trex = res_trex[0].data.evs

print("⟨Z_i⟩ with TREX:", val_trex)

Output:

⟨Z_i⟩ with TREX: [-0.92086303 -0.90253698 -0.91688854 -0.93308327 -0.91273953 -0.9221892
 -0.903782   -0.93046052 -0.90348934 -0.90284174 -0.89499677 -0.91630959
 -0.89508481 -0.92233292 -0.90497141 -0.91867138 -0.91421506 -0.91646851
 -0.9268364  -0.91669818]

Here we compare the mean absolute errors with and without TREX.

ideal = np.cos(7 * np.pi / 8)

mean_abs_error_no_trex = np.mean(np.abs(val_no_trex - ideal))
mean_abs_error_trex = np.mean(np.abs(val_trex - ideal))

print("The mean absolute error without TREX is: ", mean_abs_error_no_trex)
print("The mean absolute error with TREX is: ", mean_abs_error_trex)

Output:

The mean absolute error without TREX is:  0.05238453251128674
The mean absolute error with TREX is:  0.011980732757900719

zx = list(range(num_qubits))

import matplotlib.pyplot as plt

fig, ax = plt.subplots()

# Add values with no DD
ax.scatter(zx, val_no_trex, c="blue", linestyle="-", label="No TREX", alpha=0.7)

## Add values with DD
ax.scatter(zx, val_trex, c="red", linestyle="-", label="With TREX", alpha=0.7)


# Add ideal reference line
ax.axhline(y=-0.924, color="black", linestyle="--", linewidth=2, label="Ideal ⟨Z⟩")


# Add labels and plot
ax.set_xlabel("Qubits")
ax.set_ylabel("<Z>")
ax.legend()
ax.set_title("<Z> in RX(7pi/8)|0> state with and without TREX")
ax.grid(True)

plt.show()

Output:

TREX certainly appears to have helped. The expectation value on every qubit has improved, with the sole exception of qubit 16. Such an outlier is to be expected. We can quantify the overall application of TREX to the circuit by taking the average values.

mean_no_trex = sum(val_no_trex) / num_qubits
mean_trex = sum(val_trex) / num_qubits

print("Mean without TREX", mean_no_trex)
print("Mean with TREX", mean_trex)

Output:

Mean without TREX -0.8714950000000001
Mean with TREX -0.9137729590315635

On average, TREX moved our results about 2.9% closer to the ideal results, and left us with a result that is, on average, different from the theoretical ideal by just 0.4%. This is great, and this low-overhead improvement is exactly why TREX is the default error mitigation on Estimator. But we must once more emphasize that TREX is removing a bias in readout error, not correcting all our errors. We chose an initial state that had a good chance of being affected by asymmetry in readout. The effect would have been different if we had prepared each qubit in a different state. For a trivial argument why this must be the case, consider what would happen if we left each qubit in the ground state $|0\rangle$. There would be no opportunity to observe (or correct) an asymmetry in readout error between the states $|0\rangle$ and $|1\rangle$.

Zero-noise extrapolation (ZNE)

Zero-noise extrapolation (ZNE) is one of the most widely used error-mitigation techniques. The general strategy is also easy to understand, though the details of its implementation are complex. ZNE works by:

• Running a target circuit at multiple "noise levels" — modifications to the circuit that cause its noise to scale in predictable ways
• Estimating the observable of interest for each noise scale.
• Extrapolating back to zero noise using a fitting model (linear, polynomial, exponential, and so on).

The result is an estimate of what the observable would have been in the zero-noise limit implied by the extrapolation model.

ZNE does not guarantee that noise increases smoothly or linearly as we modify circuits. Instead, it assumes that noise scaling is sufficiently well-behaved to allow extrapolation. The complexity in ZNE is in how one modifies circuits to satisfy this assumption as well as possible. We will examine some approaches in this lesson.

ZNE resource scaling

In ZNE, the overhead is mostly quantum overhead. That is, to obtain expectation values at $k$ noise levels, you have to carry out your experiment on $k$ circuits, each with many shots. To maintain a fixed statistical precision, the total number of samples generally increases as additional noise-scaled circuits are introduced.

The number of noise levels you examine, $k$, will depend on the method used to scale the noise. For one simple noise-scaling technique (called "gate-folding", discussed below), $k$ is commonly between 3 and 5. A different method (probabilistic noise amplification (PNA)) allows for non-integer noise factors, and typically results in lower noise levels than gate-folding. The total overhead scales roughly as:

where $S$ is the baseline shot count. Another thing to keep in mind is that if noise scaling involves repeated gates, the circuit depth might increase by a factor proportional to the noise-scaling value. The scaling of noise that comes with increased depth is the point, but keep the total depth in mind.

ZNE is well-suited to problems in which the following points apply:

• You need a general-purpose mitigation method effective on many error types.
• Circuit depth is not too large; ZNE fails when the scaled circuit becomes too noisy to extrapolate.
• You can tolerate moderate runtime overhead for improved accuracy.

In practice this makes ZNE ideal for many variational circuits, expectation-value estimation, and other near-term quantum algorithms. It is implemented in Estimator using resilience_level = 2. Note that ZNE extrapolates in a space that has a single outcome as a function of noise throughout the circuit. That means it is amenable to estimation problems, but not to sampling problems. It is not designed to work on distributions of bitstrings.

Noise factors

How exactly can we increase the level of noise in a circuit? There are two main ways, although many variants exist.

1. Replace a gate or block of gates $G$ with $GG^\dagger G$. Because gates on quantum computers are unitary, we know that $G^\dagger = G^{-1}$ so that logically: $GG^\dagger G = G$. However, every repetition of $G$ or $G^\dagger$ adds similar or identical error, allowing us to scale up noise/error levels. Note that the replacement shown above takes us from the original noise of gate $G$ (noise level 1) to three times that noise level in $GG^\dagger G$ (noise level 3). In practice, we can repeat the last two gates $\ell$ times: $G(G^\dagger G)^\ell$ to obtain a noise level $2\ell+1$. This is called gate folding, and it is a very common noise scaling technique.

2. There are other methods for increasing the noise level, such as probabilistic error amplification (PEA). In PEA, we first learn a suitable noise model for the circuit, typically expressed as a Pauli channel. Noise is then amplified by sampling from a rescaled version of this learned noise channel. By adjusting the sampling, you can continuously scale the effective noise strength, including to non-integer values.

In practice, gate folding can very reliably scale errors, but it can also lead to very deep circuits. If you only need three noise levels (1, 3, and 5), it means increasing your circuit depth by a factor of 5 and still obtaining reasonable results. Thus, if your target circuit has a transpiled, two-qubit depth of 60 or more, PEA is likely a better fit than gate-folding.

Gates can be folded either locally (meaning you repeat $GG^\dagger G$ for each two-qubit gate as soon as it occurs), or they can be grouped into a block of gates and the block can be repeated (like gates $ABC\rightarrow ABC(ABC)^\dagger ABC$). Often local folding is preferable for assumptions made by ZNE, but both are possible.

A lot of nuance goes into the statistical treatment of noise-scaled results in ZNE. For example, at five times the noise level your results are certain to have higher variance and be less reliable. Should that be treated with the same weight as lower-noise data points? No. However, for the purposes of seeing ZNE in action, we will ignore some of this complexity for a toy model experiment, in which we will build folded circuits by hand, and do a naive extrapolation with equal weightings to illustrate the key steps. But we will follow this with an experiment using the ZNE functionality built into Qiskit, and we will see that it does a substantially better job of handling the complexity we are skipping.

ZNE by hand

We start by selecting a quantum circuit and an observable for which we can calculate the ideal expectation value in the absence of noise and errors. We choose a GHZ state $|\psi\rangle = \frac{1}{\sqrt{2}}\left(|0\rangle^{\otimes N}+|1\rangle^{\otimes N}\right)$, with an even number of qubits. Our observable is $Z^{\otimes N}$. We choose this for a number of reasons:

• The expectation value $\langle Z^{\otimes N}\rangle = 1$ is trivial to calculate in the ideal case. Because the number of qubits is even, we either get a 1 from every $Z$, multiplied together, or we get a -1 from every $Z$, multiplied together an even number of times.
• This will be fairly sensitive to errors and noise, since a single error on one qubit can flip the overall outcome from that shot from +1 to -1.
• The expectation value is far from 0. ZNE can encounter difficulties when the ideal expectation value is very near zero, because high noise levels will tend to randomize outcomes, making them positive as often as negative. The result is that high noise values will end up yielding values close to the ideal; this can make extrapolation very tricky.

With this motivation in mind, let's build our circuits.

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler import generate_preset_pass_manager

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Estimator,
    EstimatorOptions,
)

num_qubits = 10  # must be even

# # Noise amplification factors used for folding (2* integer +1 for gate folding)
noise_steps = [1, 3, 5]
ghz_amp = []

for j in noise_steps:
    qc = QuantumCircuit(num_qubits)
    qc.h(0)
    # Since CX^2 = I, odd repetitions implement folded versions of the gate.
    for i in range(num_qubits - 1):
        for gates in range(j):
            qc.cx(i, i + 1)
    qc.barrier()
    ghz_amp.append(qc)

Note the subtlety that here $\text{CNOT}^{-1} = \text{CNOT}^\dagger = \text{CNOT}$. Even though gate folding was advertised as consisting of operators of the form $GG^\dagger G$, here that same structure looks like CNOT CNOT CNOT.

Let's make sure that our $0^\text{th}$ circuit is indeed the normal GHZ circuit:

ghz_amp[0].draw("mpl")

Output:

You can look at the other circuits to verify they have the expected structure.

Now we build our observables.

pauli_string = "Z" * num_qubits
observable = SparsePauliOp(pauli_string)
observables = []
for j in range(len(noise_steps)):
    observables.append(observable)

We build a pass manager and apply it to our circuits. We then take the layouts and apply them to our observables.

pm = generate_preset_pass_manager(
    backend=backend_real,
    optimization_level=0,
)

qcs_transpiled = []
obs_transpiled = []

for j in range(len(ghz_amp)):
    # Transpile the circuit
    qc_pm = pm.run(ghz_amp[j])
    qcs_transpiled.append(qc_pm)
    # Extract the final layout
    layout = qc_pm.layout
    # Apply the layout to the observable
    obs_transpiled.append(observable.apply_layout(layout))

Format the pubs.

pubs = [(qc, obs) for qc, obs in zip(qcs_transpiled, obs_transpiled)]

Now we run Estimator. We would like to focus on the role of ZNE, and so we might want to set resilience_level = 0 to avoid conflating effects of twirling with the effect of ZNE. However, we also want to compare our "by-hand" ZNE to the ZNE settings built into Qiskit, which include twirling. We thus set resilience_level = 1 here. Setting resilience_level = 2 invokes the Qiskit ZNE settings, which we will use shortly.

# This took 21 sec on ibm_fez. Your actual runtime might vary.
options = EstimatorOptions(resilience_level=1)

estimator = Estimator(
    mode=backend_real,
    options=options,
)

job = estimator.run(pubs)

res = job.result()

Now we print out the noise levels and expectation values.

# We save this for later comparisons
value_no_zne = res[0].data.evs

# Print all values used for extrapolation
print("Noise level 1: ", res[0].data.evs)
print("Noise level 3: ", res[1].data.evs)
print("Noise level 5: ", res[2].data.evs)

Output:

Noise level 1:  0.8327970939283861
Noise level 3:  0.5612869745718734
Noise level 5:  0.3130254281266217

This is the part where we are being the most cavalier: we do a linear fit to these data without any attention being paid to the variance of each value we obtained, and without trying to theoretically predict whether the fit needs to be linear. Again, the built-in ZNE does a much better job than what we're doing by hand.

# Data
noise_factors = np.array([1, 3, 5])
values = np.array([res[0].data.evs, res[1].data.evs, res[2].data.evs])

# Linear fit: E(lambda) = m * lambda + b
coeffs = np.polyfit(noise_factors, values, deg=1)
slope, intercept = coeffs

print(f"Best-fit slope: {slope}")
print(f"Zero-noise extrapolated value: {intercept}")

Output:

Best-fit slope: -0.12994291645044112
Zero-noise extrapolated value: 0.9588652482269501

We go ahead and plot this prior to interpretation.

import matplotlib.pyplot as plt

x_fit = np.linspace(0, 5.5, 100)
y_fit = slope * x_fit + intercept
y_ideal = 1.0

plt.figure()
plt.plot(noise_factors, values, "o", label="Measured data")
plt.plot(x_fit, y_fit, "-", label="Linear fit")
plt.axhline(y_ideal, linestyle="--", label="Ideal")
plt.axvline(0, linestyle="--", color="gray")
plt.scatter(0, intercept, color="red", zorder=3, label="Zero-noise estimate")

plt.xlabel("Noise amplification factor")
plt.ylabel(r"$\langle Z^{\otimes N} \rangle$")
plt.legend()
plt.show()

Output:

At the making of this lesson, from the extrapolation to the zero-noise case, we get an expectation value of 0.959. Your values might differ. This is substantially better than the raw circuit's expectation value of 0.83, around 13% closer to the ideal value. But you can also see that we are using one data point where the expectation value was down at 0.31. Obviously this data point is of poorer quality. If that data point had been given lower weighting in our linear approximation, and the line were fitted more to the first two data points, we can see the zero-extrapolated estimate would have been higher (nearer the ideal value). This sort of consideration is handled automatically in the ZNE built into Estimator.

ZNE built into Estimator

Now let's see what happens when we use the ZNE functionality built into Estimator. We start by defining a simple GHZ circuit, without any gate folding.

num_qubits = 10  # must be even

qc = QuantumCircuit(num_qubits)
qc.h(0)

for i in range(num_qubits - 1):
    qc.cx(i, i + 1)

qc.barrier()

Output:

CircuitInstruction(operation=Instruction(name='barrier', num_qubits=10, num_clbits=0, params=[]), qubits=(<Qubit register=(10, "q"), index=0>, <Qubit register=(10, "q"), index=1>, <Qubit register=(10, "q"), index=2>, <Qubit register=(10, "q"), index=3>, <Qubit register=(10, "q"), index=4>, <Qubit register=(10, "q"), index=5>, <Qubit register=(10, "q"), index=6>, <Qubit register=(10, "q"), index=7>, <Qubit register=(10, "q"), index=8>, <Qubit register=(10, "q"), index=9>), clbits=())

We define the same observable as before.

pauli_string = "Z" * num_qubits
observable = SparsePauliOp(pauli_string)

We once again define a pass manager, using optimization_level = 1. We will return to this choice below.

pm = generate_preset_pass_manager(
    backend=backend_real,
    optimization_level=1,
)

Transpile the circuit and apply the layout to the observables.

# Transpile the circuit
qc_transpiled = pm.run(qc)

# Extract the final layout
layout = qc_transpiled.layout

# Apply the layout to the observable
obs_transpiled = observable.apply_layout(layout)

We run the Estimator using ZNE and twirling, and we compare the outcomes to our ZNE by-hand implementation. In particular, we compare them to the result without ZNE (value_no_zne).

# This used 18 sec on ibm_fez. Your actual QPU time might vary.

options_zne = {
    "resilience_level": 2,
    "resilience": {
        "zne": {
            "noise_factors": [1, 3, 5],
        }
    },
    "twirling": {
        "enable_gates": True,
        "enable_measure": True,
        "num_randomizations": "auto",
    },
}

estimator_zne = Estimator(
    mode=backend_real,
    options=options_zne,
)

job_zne = estimator_zne.run([(qc_transpiled, obs_transpiled)])

value_zne = job_zne.result()[0].data.evs

We can retrieve the values at each noise factor using the syntax job_zne.result()[0].data.evs_noise_factors. As before, we could naively use these results to get a linear fit and extrapolate to the zero-noise case, but the estimator with resilience_level = 2 has already done a better job, taking variance and other considerations into account. Below we show what a linear extrapolation would yield, and we show what built-in ZNE actually yielded.

# Data
noise_factors = np.array([1, 3, 5])
values_qiskit_zne = np.array(job_zne.result()[0].data.evs_noise_factors)

# Linear fit: E(lambda) = m * lambda + b
coeffs = np.polyfit(noise_factors, values_qiskit_zne, deg=1)
slope, intercept = coeffs

print(f"Best-fit slope: {slope}")
print(f"Zero-noise extrapolated value: {intercept}")

Output:

Best-fit slope: -0.11745451044574956
Zero-noise extrapolated value: 0.9045922788100509

print("Exp val ideal: 1.0")
print("Exp val with no ZNE: ", value_no_zne)
print("Exp val with ZNE: ", value_zne)
print("Fractional improvement: ", (value_zne - value_no_zne) * 100 / 1.0, "%")

Output:

Exp val ideal: 1.0
Exp val with no ZNE:  0.8327970939283861
Exp val with ZNE:  0.9877508802932561
Fractional improvement:  15.495378636487 %

import matplotlib.pyplot as plt

x_fit = np.linspace(0, 5.5, 100)
y_fit = slope * x_fit + intercept
y_ideal = 1.0

plt.figure()
plt.plot(noise_factors, values, "o", label="Measured data")
plt.plot(x_fit, y_fit, "-", label="Linear fit")
plt.axhline(y_ideal, linestyle="--", label="Ideal")
plt.axvline(0, linestyle="--", color="gray")
plt.scatter(0, intercept, color="red", zorder=3, label="Linear")
plt.scatter(0, value_zne, color="black", zorder=3, label="Proper ZNE")

plt.xlabel("Noise amplification factor")
plt.ylabel(r"$\langle Z^{\otimes N} \rangle$")
plt.legend()
plt.show()

Output:

We see that the built-in ZNE improved our estimate by over 15%, and worked quite a bit better than our attempt by hand.

Before moving on, note that throughout this ZNE exploration, we have set optimization_level = 1 in the pass manager. Higher levels of optimization reduce error, and sometimes the error reduction can make the errors scale less predictably. As is often the case, there is a balance: reduce errors, except when you are intentionally increasing them by a known factor.

Combining methods

ZNE relies on the following:

• A well-defined noise-scaling parameter $\lambda$
• Sufficiently smooth and predictable scaling of expectation values with $\lambda$
• The same effective noise channel, at different strengths
• No bias introduced that does not scale with $\lambda$

ZNE combines cleanly with techniques that remove non-scaling biases (like M3) or make noise more regular (like twirling), but conflicts with techniques that either attempt to cancel noise entirely (PEC) or reshape it dynamically (DD).

In particular, M3 addresses only measurement/readout error, which in practice can be treated as separate from gate errors (which are the focus of error scaling processes in ZNE). M3 is largely independent of gate folding; it can be cleanly combined with ZNE. TREX is very similar, and it can be combined with ZNE.

Pauli twirling does focus on gates, unlike M3 and TREX. However, recall that Pauli twirling does not directly mitigate error/noise; rather, it changes it to make it more predictable. Pauli twirling might actually be beneficial to the ZNE process in making gate errors scale more linearly and thus make the extrapolation more accurate.

Dynamical decoupling (DD) is more nuanced. DD dynamically alters noise behavior by partially cancelling certain errors through the addition of gate sequences. If it is possible to construct circuits in such a way that DD behaves the same even when gates or blocks are folded, then in principle one could combine ZNE and DD, but in practice DD sometimes interferes with the scaling. If you want to combine them, things to consider include the following:

• Is DD applied consistently at all noise levels?
• Does folding preserve the DD structure?
• Does the resulting noise scaling still behave monotonically?

The last error mitigation method we will discuss as a standalone approach is not compatible with ZNE: probabilistic error cancellation.

Probabilistic Error Cancellation (PEC)

Unlike ZNE or measurement mitigation, PEC attempts to cancel noise at the gate level by reconstructing the ideal channel in expectation. Probabilistic Error Cancellation (PEC) is the most powerful and most expensive error-mitigation technique covered in this lesson. PEC is built around explicit inversion of the noise channel, carried out statistically rather than physically. PEC works by:

• Characterizing the noise model in terms of a basis of "error operations."
• Expressing the inverse of the noise channel as a quasi-probability distribution over physical operations.
• Sampling from this distribution during circuit execution to statistically construct an unbiased estimate of the ideal computation.

Conceptually, PEC is like virtual error correction — it "undoes" noise by combining noisy operations in carefully weighted random patterns. Procedurally, however, PEC is not closely related to error correction or detection code. PEC's unbiased estimation of the ideal, noiseless expectation value comes at the cost of a very large sampling overhead and a strong dependence on accurate noise characterization.

On IBM quantum computers, PEC is best viewed as a precision tool rather than a general-purpose mitigation layer: it is most appropriate for small circuits or small portions of larger circuits, benchmarking tasks, and high-accuracy studies where high quantum overhead is acceptable.

Characterizing the noise channel

In PEC, each ideal gate in a quantum circuit is modeled as an ideal unitary channel followed by a noise process. Let $\rho$ denote the density matrix of the quantum state. The ideal gate is represented by a unitary matrix $U$, which induces a quantum channel:

On real hardware, this ideal evolution is corrupted by noise. We model the physical implementation of the gate as an effective noisy channel:

where $\mathcal{N}$ is a noise channel capturing the stochastic errors associated with the gate. Throughout this discussion, $\mathcal{U}$, $\mathcal{N}$, and $\mathcal{E}$ are linear maps acting on density operators. The ideal channel $\mathcal{U}$ is unitary and reversible, while the noise channel $\mathcal{N}$ (and hence $\mathcal{E}$) is a completely positive, trace-preserving (CPTP) map.

PEC assumes that the noise channel $\mathcal{N}$ is known or approximately known. On IBM Quantum hardware, this information is typically obtained from device calibration data, error rates inferred from benchmarking experiments, and simplified noise models such as Pauli or depolarizing approximations for small gate sets. While PEC does not require a perfectly accurate noise model, its performance degrades rapidly as the assumed noise channel diverges from the true hardware behavior.

Quasi-probability decomposition of the inverse noise channel

The central idea of PEC is to cancel the effect of noise in expectation by constructing a stochastic approximation to the inverse of the noise channel. Formally, the goal is to apply a linear map that approximately undoes $\mathcal{N}$, such that:

in expectation over many circuit executions.

The inverse $\mathcal{N}^{-1}$ is defined here as the algebraic inverse of $\mathcal{N}$ treated as a linear superoperator acting on operators. This inverse exists only if $\mathcal{N}$ is full-rank; non-invertible noise processes cannot be exactly mitigated using PEC.

When inversion is possible, PEC expresses the inverse noise channel as a linear combination of physically implementable operations:

Here, each $\mathcal{O}_i$ is a noisy operation that can be executed on the hardware (typically native gates or short gate sequences), and the real coefficients $\eta_i$ form a quasi-probability distribution. Unlike ordinary probabilities, the coefficients $\eta_i$ might be negative or exceed 1.

The quasi-probability norm:

quantifies the cost of the noise inversion. For any nontrivial noise channel, $\gamma > 1$. On IBM Quantum hardware, this decomposition is typically performed:

• per gate or per small gate set,
• assuming restricted noise models (often Pauli-diagonal approximations),
• using backend-specific native operations and scheduling constraints.

Monte Carlo reconstruction of the ideal circuit

To estimate the expectation value of an observable $O$, PEC uses Monte Carlo sampling based on the decomposition above. For each occurrence of a noisy gate in the circuit, the following procedure is applied:

• Randomly select an operation $\mathcal{O}_i$ with probability: $$p_i = \frac{|\eta_i|}{\gamma}.$$
• Assign a classical weight: $$w_i = \mathrm{sign}(\eta_i).$$
• Execute the resulting randomized noisy circuit on the quantum backend.
• Multiply the measured outcome by the product of all accumulated weights.
• Average the weighted results over many circuit executions.

The resulting estimator $\hat{O}_{\mathrm{PEC}}$ satisfies

provided the noise model $\mathcal{N}$ accurately describes the hardware. PEC therefore introduces no systematic bias in expectation, but replaces bias with large statistical variance. Each run produces a noisy but weighted estimate; only after averaging many such runs does the ideal expectation value emerge.

Conceptual interpretation

PEC is sometimes described as virtual error correction. While this phrasing is not entirely incorrect, it should not be confused with true quantum error correction. PEC involves no logical encoding, no syndrome extraction, and no real-time feedback. Errors are not prevented or corrected during circuit execution.

Instead, PEC cancels noise statistically, by combining many noisy circuit realizations with signed classical weights so that errors cancel only after classical post-processing. Information is not protected during execution; ideal results are reconstructed after the fact.

Resource scaling and cost

The dominant cost of PEC stems from variance amplification. The variance of the PEC estimator scales as:

implying:

This scaling is fundamental and cannot be improved through classical post-processing.

For circuits containing many noisy gates, each gate contributes a factor $\gamma_k > 1$, causing the total quasi-probability norm to grow approximately multiplicatively:

where $d$ is the circuit depth. On IBM Quantum hardware, this typically leads to exponential growth in sampling cost with circuit depth, with particularly strong sensitivity to two-qubit gate errors. This scaling places practical limits on the applicability of PEC.

PEC is currently practical only when one or all of the following conditions are met:

• Circuits are very small or shallow (often $\leq 10-20$ gates),
• Noise models are simple and well-validated (for example, Pauli noise),
• High statistical cost is acceptable,
• The goal is accuracy rather than throughput.

For this reason, PEC is often applied to:

• Isolated subcircuits within larger workflows,
• Short-time evolution steps,
• State preparation or measurement routines,
• Benchmarking and validation experiments.

The IBM software stack (IBM Quantum Compute and error-mitigation tooling) reflects this by treating PEC as a specialized option, not a default mitigation layer.

Typical PEC use cases include:

• Precision benchmarking of quantum hardware,
• Verification of small quantum algorithms,
• High-accuracy quantum chemistry on reduced systems,
• Cross-checking ZNE or other mitigation methods.

PEC's relationship to other mitigation methods

Many other mitigation techniques trade accuracy for scalability; PEC does just the opposite, focusing on accuracy above all else.

Because PEC might only be used on a subcircuit, it is especially important here to ask what other mitigation methods can be combined with it. To reason clearly about compatibility with PEC, it helps to categorize techniques by what they do to the noise. Consider these general categories:

1. Noise cancellation - attempt to invert error channels
2. Noise scaling - deliberately increase noise to extrapolate
3. Noise reshaping - change the structure or spectrum of noise
4. Bias removal - correct errors that do not scale with gate noise

PEC sits squarely in category (1). Any method that changes or rescales the noise channel must be examined carefully.

Summary table (high-level)

PEC and ZNE address the same error sources (gate errors) via mutually incompatible assumptions or workstreams, or both. For example, in ZNE, we assume error can be expanded in a parameter $\lambda$:

and that the error can be effectively scaled to model this expansion. PEC would undo the very error we rely on for the scaling. The parameter $\lambda$ would become meaningless. Another way of thinking about it is:

• Applying ZNE before PEC changes the noise channel you are trying to invert.
• Applying PEC before ZNE removes the very noise ZNE relies on to extrapolate.

PEC and ZNE should be considered alternative error-mitigation paradigms, not composable layers. Comparing them side-by-side is instructive, while combining them is not.

PEC and Pauli twirling are generally compatible and often beneficial. Twirling might help PEC in that twirling converts coherent errors into stochastic Pauli errors. This simplifies the noise channel, making it easier to represent and invert. For PEC, this is extremely valuable, for the following reasons:

• Quasi-probability decompositions are much simpler for Pauli channels
• The quasi-probability norm $\gamma$ is often significantly smaller
• Sampling variance is reduced

Many PEC demonstrations assume or build in Pauli-twirled noise. It's important to clarify that Pauli twirling does not mitigate errors by itself; it reshapes noise into a form that probabilistic error cancellation can handle more efficiently.

PEC can and should be combined with TREX. These approaches are fully compatible and play orthogonal roles. TREX operates only at measurement, removes readout bias via randomized twirls of the measurement channel, and does not affect gate noise or circuit dynamics. PEC operates on gate noise. Because measurement error is additive rather than compositional, it is independent of gate sequences. It is safe and often desirable to apply TREX alongside PEC.

The ordering of these methods is what you might expect:

1. Execute PEC-randomized circuit
2. Apply TREX to mitigate readout error
3. Apply PEC weighting and averaging

PEC and DD are generally incompatible unless the noise characterization used for PEC is performed with DD active. PEC assumes that $N\text{hardware noise} = \mathcal N$ is stationary, repeatable, and accurately characterizable. DD explicitly reshapes the noise channel in time by modulating system-environment interactions and introducing non-trivial pulse-timing effects. As a result, the effective noise with DD is not the same channel as without DD. This is made even more complex by the fact that the noise channel can depend sensitively on circuit scheduling. The end result is that noise inversion becomes ill-defined or invalid.

It is theoretically possible to characterize noise with DD active and build a PEC model accordingly, but this is highly advanced and not currently a recommended workflow.

Applying PEC by hand

PEC is complex and it is not recommended that you implement it by hand. Still, we want to work through the math of a basic application on a single qubit to understand what steps the process entails. We choose a mixed but nontrivial state so that we can see numerical effects clearly.

Take a qubit that should be in the pure state

The ideal density matrix is:

We've chosen this density matrix because it is not diagonal and is only for a single qubit so that it is tractable to examine "by hand".

We now apply a very simple, but realistic noise model: a single-qubit depolarizing channel, with small strength $p$.

To keep the arithmetic friendly, take $p = 0.15$. This channel is completely positive, trace-preserving (CPTP), as it must be. It is a realistic noise channel for real hardware, or at least realistic for describing some of the noise. Further, this noise channel is analytically invertible, as a linear map. In other words, there exists an inverse matrix, but that doesn't mean that the matrix corresponds to a CPTP physical channel, and it certainly doesn't correspond to a unitary that could be implemented on hardware.

Now we apply our depolarizing noise channel to our single-qubit density matrix and obtain:

With $p = 0.15$: $\frac{4p}{3} = 0.20$ and $\frac{2p}{3} = 0.10$.

So:

Applying this to our state:

Compute explicitly:

This is still a physical density matrix, and it still has trace one — but it is more mixed now.

Now we compute the ideal inverse channel, but being very prepared for the possibility that it might not be physical. We will check whether it is. Mathematically, the inverse linear map is:

If this inverse action were itself a valid physical channel, then we would be tempted to try to implement it exactly. Of course, we would still need to figure out how to implement it using unitary operations and classical steps, but implementing this would be the goal. However, it turns out this is not CPTP, and there is generally no guarantee that such inverse mappings will be CPTP. That is critical to our current discussion, because it is the reason PEC exists.

You might be able to see that it is not CPTP just by looking at it. A sure-fire check is to apply it to a density matrix and see if it returns physical results. This is left as an exercise to check your understanding.

Check your understanding

Check whether the inverse mapping from above $\mathcal{E}^{-1}$ is CPTP by applying it to a new density matrix:

Interpret the results as being physical or non-physical.

• Answer
  $$\begin{aligned} \mathcal{E}^{-1}(\rho_{\text{test}}) & = \frac{1}{0.8} \left(\rho_{\text{test}} - 0.1 I\right)\\ & = \frac{1}{0.8} \left(\begin{pmatrix} 0.95 & 0 \\ 0 & 0.05 \end{pmatrix} - \begin{pmatrix} 0.1 & 0 \\ 0 & 0.1 \end{pmatrix}\right)\\ & = \begin{pmatrix} 1.0625 & 0 \\ 0 & -0.0625 \end{pmatrix} \end{aligned}$$

  As a physical density matrix, this corresponds to a classical probability of greater than 1 that the system is in the state $|0\rangle$ and a negative probability of being in the state $|1\rangle$. Clearly this mapping can yield nonsense density matrices. Further, by inspection we can see that the eigenvalues of the resulting density matrix are 1.0625 and -0.0625, meaning not all eigenvalues are positive. Thus the mapping is not "completely positive", the "CP" in CPTP.

Check your understanding

Check by explicit calculation that the inverse operation above $\mathcal{E}^{-1}$ undoes the noise channel and returns the noise-free state.

• Answer

  Apply it to the noisy state:

  $$\begin{aligned} \rho_{\text{recovered}} & = \mathcal{E}^{-1}(\rho_\text{noisy})\\ & =\frac{1}{0.8} \left( \begin{pmatrix} 0.7 & 0.346 \\ 0.346 & 0.3 \end{pmatrix} - \begin{pmatrix} 0.1 & 0 \\ 0 & 0.1 \end{pmatrix} \right)\\ & = \frac{1}{0.8} \begin{pmatrix} 0.6 & 0.346 \\ 0.346 & 0.2 \end{pmatrix}\\ & \approx \begin{pmatrix} 0.75 & 0.433 \\ 0.433 & 0.25 \end{pmatrix} \end{aligned}$$

  This is exactly the original state (up to rounding).

We restate to be absolutely clear: we cannot invert the noise channel exactly; the inverse is not CPTP, nor is it a valid physical channel. However, we can effectively invert the channel as a linear combination of physical channels.

For the depolarizing channel:

where: $\mathcal{X}(\rho) = X\rho X$ etc.

The inverse can be written as:

with:

• $\mathcal{P}_j$ = physically implementable Clifford channels
• $\eta_j$ = quasi-probabilities
• Some $\eta_j < 0$

This is the core idea behind PEC:

• Sample channels according to $|\eta_j|$∣
• Attach signs to measurement outcomes
• Expectation values converge to the noise-free value

We apply PEC "by hand" to the estimation of $\langle Z \rangle$ by using the density matrix $\rho_\text{noisy}$ and see if we can recover the results for $\rho_\text{ideal}$. For $\rho_\text{ideal}$ we expect:

Note also that for the noisy density matrix, we obtain:

It's the right order of magnitude but significantly far from the ideal case.

Now we step through how, with the PEC process, we can reconstruct the unbiased expectation value from the noisy density matrix. In real workflows, this is probabilistic and will be automated. In this "by-hand" workflow, we will choose probabilistic values that are realistic, but perhaps convenient for our presentation. This should not be taken as a claim that PEC is perfect. This is simply to show you a plausible scenario that could happen inside the PEC packages.

Recall that our noise channel is of the form:

The inverse can be written (exactly):

based on symmetry of the transformation and the completeness of the Pauli operators. Inverting the noise channel somehow in the estimation process should yield an estimator that is unbiased and normalized. To do this, we must at least determine the coefficents $\alpha$ and $\beta$ and then think about how to implement the inverse transformation. When done by a noise learner, the inversion of this noise channel would be done at the channel level, independent of a specific case like our density matrix. The coefficients $\alpha$ and $\beta$ would come from this inversion.

It is simpler for this introduction to focus on the channel's effect on our density matrix. This allows us to work with familiar objects like density matrices and expectation values. To that end, consider what must be true about the effects of this inverse operation on our estimation process: the estimator produced by PEC must be normalized and unbiased. This corresponds to two conditions:

1. Normalization: To ensure that this operation does not change the scale of the observable, we require that $\alpha+3\beta = 1$. To see why this is the case, we could either consider the action of the inverse on a density matrix, or consider its action on the expectation value of something trivial like $\mathcal{I}$ itself.

2. Calibration: In a real workflow, we use calibration circuits to learn noise behavior at a gate level. We then use the gate-level knowledge of noise behavior, and the linearity of expectation values to apply the learned noise model to the larger circuit. In our case, we started with artificial knowledge that the ideal expectation value $\langle Z \rangle = 0.5$. But it would be very easy to come up with a state for which we know theoretically that this would be the ideal value, and then use actual measurements for calibration. A key point is that for PEC, one typically does not learn a full noise model of the entire circuit. Instead, one learns local noise models for individual gates or circuit layers and uses those to construct the quasi-probability decomposition.

We use these conditions to determine $\alpha$ and $\beta$. We consider the ansatz