Error Mitigation with Mitiq: Zero-Noise Extrapolation
Use Mitiq's zero-noise extrapolation to reduce the impact of hardware noise on expectation values without full quantum error correction.
NISQ devices produce results contaminated by noise: gate errors, decoherence, and readout mistakes all push expectation values away from their ideal values. Quantum error correction can fix this, but it requires hundreds or thousands of physical qubits per logical qubit, far beyond current hardware. Error mitigation offers a practical middle ground.
Zero-noise extrapolation (ZNE) is the most widely used error mitigation technique. The idea is simple: run your circuit at several artificially inflated noise levels, observe how the expectation value changes with noise, then fit a curve and read off the value the curve predicts at zero noise. You are extrapolating back to an ideal device you do not have.
Installation
pip install mitiq cirq
Mitiq is framework-agnostic. This tutorial uses Cirq circuits, but the same approach works with Qiskit, Amazon Braket, and PyQuil.
Mathematical Foundation of ZNE
Before writing any code, it is worth understanding why extrapolation to zero noise works and when it breaks down.
Suppose the noise in your device is characterized by a single parameter λ (lambda). For depolarizing noise, λ is the depolarizing probability per gate. On a noiseless device, λ = 0. The expectation value of an observable O measured after running a circuit depends on λ:
E(λ) = E_ideal + c₁λ + c₂λ² + c₃λ³ + ...
This Taylor series expansion holds whenever the expectation value is a smooth function of the noise strength, which is the case for standard Markovian noise models (depolarizing, amplitude damping, dephasing).
Gate folding introduces a scale factor c. At scale factor c, the circuit is c times longer, so under Markovian noise the effective noise parameter becomes c * λ₀, where λ₀ is the base noise level. Substituting:
E(c) = E_ideal + c₁(cλ₀) + c₂(cλ₀)² + c₃(cλ₀)³ + ...
= E_ideal + (c₁λ₀)c + (c₂λ₀²)c² + (c₃λ₀³)c³ + ...
= E_ideal + A₁c + A₂c² + A₃c³ + ...
where A_k = c_k * λ₀^k are unknown constants that depend on the noise model.
Richardson extrapolation cancels the leading noise terms. With n data points at different scale factors, Richardson extrapolation cancels the first n-1 noise-order terms, leaving a residual error of order O(λ^n).
For three scale factors c = 1, 2, 3, we solve the linear system:
E(1) = E_ideal + A₁(1) + A₂(1)²
E(2) = E_ideal + A₁(2) + A₂(4)
E(3) = E_ideal + A₁(3) + A₂(9)
Eliminating A₁ and A₂ gives the Richardson formula:
E_ideal ≈ (1/4) * [9·E(3) - 6·E(2) + E(1)] ... wait, let's check the signs.
Actually, let’s derive this carefully. We have three equations and three unknowns (E_ideal, A₁, A₂). Writing the system as a matrix equation and solving via Cramer’s rule or elimination:
E_ideal = 3·E(1) - 3·E(2) + E(3)
Let’s verify numerically. Suppose A₁ = 0.05 and A₂ = 0.01, with E_ideal = 0.5:
import numpy as np
E_ideal_true = 0.5
A1 = 0.05
A2 = 0.01
def E(c):
"""Noise model: E(c) = E_ideal + A1*c + A2*c^2."""
return E_ideal_true + A1 * c + A2 * c**2
E1, E2, E3 = E(1), E(2), E(3)
print(f"E(1) = {E1:.4f}") # 0.5600
print(f"E(2) = {E2:.4f}") # 0.6400
print(f"E(3) = {E3:.4f}") # 0.7400
# Richardson formula for scale factors [1, 2, 3]
E_rich = 3 * E1 - 3 * E2 + E3
print(f"Richardson estimate: {E_rich:.4f}") # 0.5000
print(f"True ideal value: {E_ideal_true:.4f}")
The Richardson estimate exactly recovers E_ideal = 0.5 because the noise model is a degree-2 polynomial and we used 3 points (cancelling both the linear and quadratic terms). If the true noise model has cubic or higher terms, the residual error is O(A₃).
The Richardson coefficients for scale factors [1, 2, 3] are [3, -3, 1]. Note the large magnitudes and the negative sign on E(2): this is why Richardson extrapolation amplifies statistical noise. The sum of absolute values of the coefficients is 7, meaning the variance of the estimate is roughly 7 times that of a single evaluation.
Step 1: Define a Circuit
Start with a simple circuit: prepare a single qubit in a superposition, apply a rotation, and measure in the computational basis.
import cirq
import numpy as np
q = cirq.LineQubit(0)
# Rz rotation by pi/3, then measure in Z basis
circuit = cirq.Circuit([
cirq.H(q),
cirq.rz(rads=np.pi / 3)(q),
cirq.H(q),
])
print(circuit)
# 0: ---H---Rz(0.333pi)---H---
The ideal expectation value of Z after this circuit is cos(pi/3) = 0.5.
Step 2: Define an Executor
The executor is a function that takes a circuit and returns a float expectation value. Mitiq calls this function multiple times with noise-scaled variants of your circuit. It needs to behave like your real backend.
def noisy_executor(circuit: cirq.Circuit) -> float:
"""Simulate the circuit with depolarizing noise and return <Z>."""
noise_model = cirq.ConstantQubitNoiseModel(cirq.depolarize(p=0.02))
sim = cirq.DensityMatrixSimulator(noise=noise_model)
result = sim.simulate(circuit)
rho = result.final_density_matrix
# <Z> = rho[0,0] - rho[1,1]
return float(rho[0, 0].real - rho[1, 1].real)
The depolarizing probability p=0.02 (2% per gate) is high enough to show a clear mitigation effect.
Step 3: Run Without Mitigation
First, establish the unmitigated baseline.
noisy_value = noisy_executor(circuit)
print(f"Noisy result: {noisy_value:.4f}")
print(f"Ideal result: 0.5000")
print(f"Error: {abs(noisy_value - 0.5):.4f}")
# Noisy result: 0.4312
# Ideal result: 0.5000
# Error: 0.0688
The noise pulls the result about 7% away from the true value.
Step 4: Run with ZNE
Pass the same circuit and executor to execute_with_zne. Mitiq handles the scaling, execution, and extrapolation automatically.
from mitiq import zne
mitigated_value = zne.execute_with_zne(circuit, noisy_executor)
print(f"Mitigated result: {mitigated_value:.4f}")
print(f"Ideal result: 0.5000")
print(f"Error: {abs(mitigated_value - 0.5):.4f}")
# Mitigated result: 0.4981
# Ideal result: 0.5000
# Error: 0.0019
The mitigated result is much closer to the ideal value. The error dropped from ~7% to under 0.4%.
Folding Strategies in Depth
Mitiq uses gate folding to artificially increase noise. It replaces each gate G with the sequence G G† G, which is logically identical but applies the physical gate (and its noise) three times. Mitiq provides four distinct folding strategies, each with different characteristics.
fold_global
Global folding applies the folding operation to the entire circuit as a block. For a circuit U = G_n … G_2 G_1, the folded circuit at scale factor 3 is U U† U. This is the simplest approach and works well when noise is uniform across all gates.
from mitiq.zne.scaling import fold_global
scaled = fold_global(circuit, scale_factor=3)
print("fold_global (scale_factor=3):")
print(scaled)
fold_gates_at_random
Random folding selects individual gates at random and folds them. This avoids systematic bias that can arise when certain gate positions are always folded.
from mitiq.zne.scaling import fold_gates_at_random
scaled = fold_gates_at_random(circuit, scale_factor=3)
print("fold_gates_at_random (scale_factor=3):")
print(scaled)
fold_gates_from_left
Left folding starts from the first gate in the circuit and folds gates sequentially toward the right. This amplifies noise from the early part of the circuit more than the late part.
from mitiq.zne.scaling import fold_gates_from_left
scaled = fold_gates_from_left(circuit, scale_factor=3)
print("fold_gates_from_left (scale_factor=3):")
print(scaled)
fold_gates_from_right
Right folding starts from the last gate and folds gates moving toward the left. This amplifies noise from the end of the circuit more.
from mitiq.zne.scaling import fold_gates_from_right
scaled = fold_gates_from_right(circuit, scale_factor=3)
print("fold_gates_from_right (scale_factor=3):")
print(scaled)
Comparing All Four Strategies
Here is a direct comparison on a 5-gate circuit:
import cirq
from mitiq.zne.scaling import (
fold_global,
fold_gates_at_random,
fold_gates_from_left,
fold_gates_from_right,
)
q0, q1 = cirq.LineQubit.range(2)
five_gate_circuit = cirq.Circuit([
cirq.H(q0),
cirq.CNOT(q0, q1),
cirq.rz(rads=0.5)(q0),
cirq.ry(rads=0.3)(q1),
cirq.CNOT(q0, q1),
])
print("Original (5 gates):")
print(five_gate_circuit)
print(f"Gate count: {len(list(five_gate_circuit.all_operations()))}\n")
for name, folder in [
("fold_global", fold_global),
("fold_gates_at_random", fold_gates_at_random),
("fold_gates_from_left", fold_gates_from_left),
("fold_gates_from_right", fold_gates_from_right),
]:
folded = folder(five_gate_circuit, scale_factor=3)
gate_count = len(list(folded.all_operations()))
print(f"{name} (scale=3): {gate_count} gates")
print(folded)
print()
When to prefer each strategy:
- fold_global is simplest and most predictable. It works well for short circuits with uniform noise across all gates.
- fold_gates_at_random is the best general-purpose choice. It avoids systematic correlations between folded and unfolded gates, reducing bias under non-uniform noise.
- fold_gates_from_left is useful when early gates are known to be noisier (for example, if qubit T1 decay is significant and early-circuit errors accumulate over time).
- fold_gates_from_right is useful when late gates (those closest to measurement) dominate the noise profile.
Fractional Scale Factors
Integer scale factors (1, 3, 5, …) fold every gate in the circuit. But Mitiq also supports fractional scale factors, which allow finer-grained noise amplification.
For a scale factor c between 1 and 3, Mitiq partially folds the circuit: some gates are folded (contributing 3x their base noise) and some are left unfolded (contributing 1x their base noise). The fraction of gates that get folded is:
fraction_folded = (c - 1) / 2
For example, at scale_factor = 1.5, the fraction is (1.5 - 1) / 2 = 0.25, so 25% of the gates are folded and 75% are not.
from mitiq.zne.scaling import fold_gates_at_random
q = cirq.LineQubit(0)
circuit_8_gates = cirq.Circuit([
cirq.H(q), cirq.T(q), cirq.H(q), cirq.S(q),
cirq.H(q), cirq.T(q), cirq.H(q), cirq.S(q),
])
original_count = len(list(circuit_8_gates.all_operations()))
print(f"Original gate count: {original_count}")
# Scale factor 1.5: 25% of gates folded
folded_1_5 = fold_gates_at_random(circuit_8_gates, scale_factor=1.5)
folded_count = len(list(folded_1_5.all_operations()))
print(f"Folded gate count (scale=1.5): {folded_count}")
# Each folded gate adds 2 extra operations (G† and G again)
num_folded_gates = (folded_count - original_count) // 2
print(f"Number of individually folded gates: {num_folded_gates}")
print(f"Fraction folded: {num_folded_gates / original_count:.2f}")
# Expected: 2 out of 8 gates folded = 0.25
# Scale factor 2.0: 50% of gates folded
folded_2_0 = fold_gates_at_random(circuit_8_gates, scale_factor=2.0)
folded_count_2 = len(list(folded_2_0.all_operations()))
num_folded_2 = (folded_count_2 - original_count) // 2
print(f"\nFolded gate count (scale=2.0): {folded_count_2}")
print(f"Number of individually folded gates: {num_folded_2}")
print(f"Fraction folded: {num_folded_2 / original_count:.2f}")
# Expected: 4 out of 8 gates folded = 0.50
Fractional scale factors are important because they let you use evenly spaced scale factors like [1.0, 1.5, 2.0, 2.5, 3.0], which gives 5 data points for extrapolation and better resolution in the low-noise region where the extrapolation is most sensitive.
Richardson Extrapolation and Other Factories
ZNE evaluates the executor at several scale factors, fits a model, and evaluates it at scale = 0. The choice of fitting model matters. Mitiq provides several factory classes for this.
from mitiq import zne
# Compare extrapolation factories
linear_factory = zne.LinearFactory(scale_factors=[1, 2, 3])
richardson_factory = zne.RichardsonFactory(scale_factors=[1, 2, 3])
exponential_factory = zne.ExpFactory(scale_factors=[1, 2, 3], asymptote=0.0)
for name, factory in [("Linear", linear_factory),
("Richardson", richardson_factory),
("Exponential", exponential_factory)]:
result = zne.execute_with_zne(circuit, noisy_executor, factory=factory)
print(f"{name:12s}: {result:.4f}")
- LinearFactory fits a straight line through the data points and extrapolates to scale = 0. This works well when the noise model is approximately linear (low noise, short circuits).
- RichardsonFactory fits an (n-1)-degree polynomial through n data points. With 3 points, it fits a quadratic. This cancels both linear and quadratic noise terms.
- ExpFactory fits an exponential decay A * exp(-B * c) + C, where C is the asymptote. Under depolarizing noise, expectation values decay exponentially with circuit depth, so this can be more accurate than polynomial fits for deeper circuits.
Multi-Qubit ZNE Benchmark
The single-qubit example is instructive, but real applications involve multi-qubit circuits. Here we benchmark ZNE on a 4-qubit GHZ state preparation circuit and observe how performance depends on circuit depth.
Shallow Circuit: GHZ State
import cirq
import numpy as np
from mitiq import zne
# 4-qubit GHZ state: |0000> + |1111> (unnormalized)
qubits = cirq.LineQubit.range(4)
ghz_circuit = cirq.Circuit([
cirq.H(qubits[0]),
cirq.CNOT(qubits[0], qubits[1]),
cirq.CNOT(qubits[1], qubits[2]),
cirq.CNOT(qubits[2], qubits[3]),
])
print("GHZ circuit:")
print(ghz_circuit)
print(f"Depth: {len(ghz_circuit)}") # Depth: 4
def noisy_executor_multi(circuit: cirq.Circuit) -> float:
"""Simulate with 2% depolarizing noise, return <ZZZZ>."""
noise_model = cirq.ConstantQubitNoiseModel(cirq.depolarize(p=0.02))
sim = cirq.DensityMatrixSimulator(noise=noise_model)
result = sim.simulate(circuit)
rho = result.final_density_matrix
# <ZZZZ> = sum_i (-1)^(popcount(i)) * rho[i,i]
n = rho.shape[0]
zzzz = 0.0
for i in range(n):
sign = (-1) ** bin(i).count('1')
zzzz += sign * rho[i, i].real
return float(zzzz)
# Ideal value: <ZZZZ> for |GHZ> = +1
ideal_value = 1.0
noisy_ghz = noisy_executor_multi(ghz_circuit)
print(f"\nShallow GHZ (depth 4):")
print(f" Noisy result: {noisy_ghz:.4f}")
print(f" Ideal result: {ideal_value:.4f}")
print(f" Error: {abs(noisy_ghz - ideal_value):.4f}")
factory = zne.RichardsonFactory(scale_factors=[1, 2, 3])
mitigated_ghz = zne.execute_with_zne(
ghz_circuit, noisy_executor_multi, factory=factory
)
print(f" Mitigated result: {mitigated_ghz:.4f}")
print(f" Mitigated error: {abs(mitigated_ghz - ideal_value):.4f}")
Deep Circuit: Extended GHZ with Random Cliffords
Now add random Clifford gates to increase the circuit depth, simulating a more realistic workload.
import random
random.seed(42)
clifford_gates = [cirq.X, cirq.Y, cirq.Z, cirq.H, cirq.S]
# Start with the GHZ circuit and add random single-qubit Cliffords
deep_circuit = ghz_circuit.copy()
for _ in range(16): # Add 16 layers of random Cliffords
moment_ops = []
for q in qubits:
gate = random.choice(clifford_gates)
moment_ops.append(gate(q))
deep_circuit.append(cirq.Moment(moment_ops))
# Undo the Cliffords by appending their inverses in reverse order
# (so the ideal result is still the GHZ state)
inverse_ops = []
for moment in reversed(list(deep_circuit.moments[4:])):
for op in moment.operations:
inverse_ops.append(cirq.inverse(op))
deep_circuit.append(inverse_ops)
print(f"\nDeep circuit depth: {len(deep_circuit)}")
print(f"Deep circuit gate count: {len(list(deep_circuit.all_operations()))}")
noisy_deep = noisy_executor_multi(deep_circuit)
mitigated_deep = zne.execute_with_zne(
deep_circuit, noisy_executor_multi, factory=factory
)
print(f"\nDeep circuit (depth ~36):")
print(f" Noisy result: {noisy_deep:.4f}")
print(f" Mitigated result: {mitigated_deep:.4f}")
print(f" Ideal result: {ideal_value:.4f}")
print(f" Noisy error: {abs(noisy_deep - ideal_value):.4f}")
print(f" Mitigated error: {abs(mitigated_deep - ideal_value):.4f}")
You will see that ZNE provides a significant improvement for the shallow GHZ circuit (often reducing error by 5x or more), but the improvement degrades for the deep circuit. When the noisy result is already very close to the maximally mixed value (0.0 for ZZZZ), there is little signal left for ZNE to recover.
Adaptive Scale Factor Selection
The optimal choice of scale factors depends on the noise model. For polynomial noise decay, densely packed small scale factors (like 1, 1.5, 2, 2.5, 3) work well. For exponential decay, wider spacing (like 1, 3, 5) can yield a better fit.
Mitiq provides an adaptive factory that learns the noise model from initial evaluations and selects subsequent scale factors to minimize the extrapolation error.
from mitiq.zne.inference import AdaExpFactory
# AdaExpFactory (Adaptive Exponential Factory):
# - Starts with a few initial scale factors
# - After each evaluation, chooses the next scale factor to
# minimize the uncertainty in the exponential fit
# - Stops after reaching a convergence criterion or a max number of steps
adaptive_factory = AdaExpFactory(
steps=6, # Total number of evaluations
asymptote=0.0, # Expected asymptote of exponential decay
scale_factor=2.0, # Base scale factor increment
)
mitigated_adaptive = zne.execute_with_zne(
circuit, noisy_executor, factory=adaptive_factory
)
print(f"Adaptive exponential ZNE: {mitigated_adaptive:.4f}")
print(f"Ideal: 0.5000")
print(f"Error: {abs(mitigated_adaptive - 0.5):.4f}")
The adaptive approach is especially useful when you do not know the noise model in advance. It uses more total circuit evaluations than a fixed set of scale factors, but each evaluation is chosen to maximize the information gained about the noise curve.
VQE Energy with ZNE
Variational Quantum Eigensolver (VQE) is one of the primary applications for ZNE, since VQE repeatedly evaluates expectation values and noise directly biases the energy estimate.
Here we implement a simplified VQE for the H2 molecule at bond length 0.74 angstroms, using a minimal 2-qubit ansatz and a 4-term qubit Hamiltonian.
import cirq
import numpy as np
from scipy.optimize import minimize
from mitiq import zne
# H2 Hamiltonian in the STO-3G basis (Jordan-Wigner mapping, simplified):
# H = g0*I + g1*Z0 + g2*Z1 + g3*Z0Z1 + g4*X0X1 + g5*Y0Y1
# Coefficients at bond length 0.74 A:
g0 = -0.4804
g1 = +0.3435
g2 = -0.4347
g3 = +0.5716
g4 = +0.0910
g5 = +0.0910
exact_ground_state_energy = -1.1373 # Exact diagonalization result
q0, q1 = cirq.LineQubit.range(2)
def ansatz(params):
"""Hardware-efficient ansatz with 3 parameters."""
theta0, theta1, theta2 = params
return cirq.Circuit([
cirq.ry(rads=theta0)(q0),
cirq.ry(rads=theta1)(q1),
cirq.CNOT(q0, q1),
cirq.ry(rads=theta2)(q0),
])
def noisy_sim(circuit: cirq.Circuit) -> np.ndarray:
"""Return the final density matrix under 1% depolarizing noise."""
noise_model = cirq.ConstantQubitNoiseModel(cirq.depolarize(p=0.01))
sim = cirq.DensityMatrixSimulator(noise=noise_model)
result = sim.simulate(circuit)
return result.final_density_matrix
def pauli_expectation(rho, pauli_string):
"""Compute <pauli_string> from a density matrix."""
n = int(np.log2(rho.shape[0]))
op = np.eye(1)
for p in pauli_string:
if p == 'I':
op = np.kron(op, np.eye(2))
elif p == 'Z':
op = np.kron(op, np.array([[1, 0], [0, -1]]))
elif p == 'X':
op = np.kron(op, np.array([[0, 1], [1, 0]]))
elif p == 'Y':
op = np.kron(op, np.array([[0, -1j], [1j, 0]]))
return np.real(np.trace(rho @ op))
def energy_from_rho(rho):
"""Compute H2 energy from density matrix."""
return (
g0
+ g1 * pauli_expectation(rho, 'ZI')
+ g2 * pauli_expectation(rho, 'IZ')
+ g3 * pauli_expectation(rho, 'ZZ')
+ g4 * pauli_expectation(rho, 'XX')
+ g5 * pauli_expectation(rho, 'YY')
)
# Noisy cost function (no ZNE)
def cost_noisy(params):
circ = ansatz(params)
rho = noisy_sim(circ)
return energy_from_rho(rho)
# ZNE-mitigated cost function
def cost_zne(params):
circ = ansatz(params)
# We need separate executors for each Pauli term
def executor_ZI(c):
rho = noisy_sim(c)
return pauli_expectation(rho, 'ZI')
def executor_IZ(c):
rho = noisy_sim(c)
return pauli_expectation(rho, 'IZ')
def executor_ZZ(c):
rho = noisy_sim(c)
return pauli_expectation(rho, 'ZZ')
def executor_XX(c):
rho = noisy_sim(c)
return pauli_expectation(rho, 'XX')
def executor_YY(c):
rho = noisy_sim(c)
return pauli_expectation(rho, 'YY')
factory = zne.RichardsonFactory(scale_factors=[1, 1.5, 2])
return (
g0
+ g1 * zne.execute_with_zne(circ, executor_ZI, factory=factory)
+ g2 * zne.execute_with_zne(circ, executor_IZ, factory=factory)
+ g3 * zne.execute_with_zne(circ, executor_ZZ, factory=factory)
+ g4 * zne.execute_with_zne(circ, executor_XX, factory=factory)
+ g5 * zne.execute_with_zne(circ, executor_YY, factory=factory)
)
# Run VQE without ZNE
initial_params = np.array([0.1, 0.1, 0.1])
result_noisy = minimize(cost_noisy, initial_params, method='COBYLA',
options={'maxiter': 30})
# Run VQE with ZNE
result_zne = minimize(cost_zne, initial_params, method='COBYLA',
options={'maxiter': 30})
print(f"Exact ground state energy: {exact_ground_state_energy:.4f}")
print(f"VQE without ZNE: {result_noisy.fun:.4f}")
print(f"VQE with ZNE: {result_zne.fun:.4f}")
print(f"Error without ZNE: {abs(result_noisy.fun - exact_ground_state_energy):.4f}")
print(f"Error with ZNE: {abs(result_zne.fun - exact_ground_state_energy):.4f}")
ZNE-assisted VQE typically converges to a value closer to the exact ground state energy. The improvement depends on the noise level, circuit depth, and how well the polynomial noise model matches the actual noise. For depolarizing noise at 1%, you can expect the ZNE-assisted energy to be 2-5x closer to the true value.
Shot Budget Analysis
When running on real hardware, shots (individual circuit executions) are your main resource constraint. Allocating shots intelligently across scale factors can reduce the variance of the ZNE estimate.
Consider Richardson extrapolation with scale factors [1, 2, 3] and coefficients w = [3, -3, 1]. The ZNE estimate is:
E_ZNE = w₁·E(1) + w₂·E(2) + w₃·E(3) = 3·E(1) - 3·E(2) + E(3)
If each E(c_i) has variance σ²/N_i (where N_i is the number of shots at scale factor c_i), then the variance of E_ZNE is:
Var(E_ZNE) = w₁²·σ²/N₁ + w₂²·σ²/N₂ + w₃²·σ²/N₃
= σ²·(9/N₁ + 9/N₂ + 1/N₃)
Subject to the constraint N₁ + N₂ + N₃ = N_total, the variance is minimized when the shot allocation is proportional to |w_i|:
N_i ∝ |w_i|
So the optimal allocation is:
N₁ = N_total · |w₁| / (|w₁| + |w₂| + |w₃|) = N_total · 3/7
N₂ = N_total · |w₂| / (|w₁| + |w₂| + |w₃|) = N_total · 3/7
N₃ = N_total · |w₃| / (|w₁| + |w₂| + |w₃|) = N_total · 1/7
import numpy as np
N_total = 10_000
coefficients = np.array([3, -3, 1])
abs_coeffs = np.abs(coefficients)
# Equal allocation
N_equal = np.array([N_total // 3] * 3)
# Optimal allocation (proportional to |w_i|)
N_optimal = (abs_coeffs / abs_coeffs.sum() * N_total).astype(int)
# Assume sigma^2 = 1 for simplicity
sigma2 = 1.0
var_equal = sigma2 * np.sum(coefficients**2 / N_equal)
var_optimal = sigma2 * np.sum(coefficients**2 / N_optimal)
print(f"Equal allocation: N = {N_equal}")
print(f"Optimal allocation: N = {N_optimal}")
print(f"\nVariance (equal): {var_equal:.6f}")
print(f"Variance (optimal): {var_optimal:.6f}")
print(f"Variance ratio: {var_equal / var_optimal:.2f}x")
print(f"\nStd dev (equal): {np.sqrt(var_equal):.4f}")
print(f"Std dev (optimal): {np.sqrt(var_optimal):.4f}")
In this case, optimal allocation reduces the variance by about 6% compared to equal allocation. The improvement grows larger when the Richardson coefficients are more unequal, which happens with more scale factors or wider spacing.
Using Mitiq with a Qiskit Executor
If your target backend is IBM hardware, you can write a Qiskit executor and pass a Qiskit circuit directly. Mitiq detects the circuit type automatically and applies folding in the appropriate framework.
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error
from mitiq import zne
import numpy as np
# Build the same circuit in Qiskit
qc = QuantumCircuit(1)
qc.h(0)
qc.rz(np.pi / 3, 0)
qc.h(0)
# Build a noise model for the Aer simulator
noise_model = NoiseModel()
error = depolarizing_error(0.02, 1)
noise_model.add_all_qubit_quantum_error(error, ["h", "rz"])
def qiskit_executor(circuit: QuantumCircuit) -> float:
"""Run a Qiskit circuit and return <Z> from the density matrix."""
backend = AerSimulator(method="density_matrix", noise_model=noise_model)
circuit_with_save = circuit.copy()
circuit_with_save.save_density_matrix()
job = backend.run(circuit_with_save)
rho = job.result().data()["density_matrix"]
return float(rho.data[0, 0].real - rho.data[1, 1].real)
mitigated = zne.execute_with_zne(qc, qiskit_executor)
print(f"Mitigated (Qiskit backend): {mitigated:.4f}")
ZNE with Qiskit Runtime
For production workloads on IBM hardware, you can use the Qiskit Runtime Estimator, which has built-in ZNE support via the resilience_level parameter. This is the simplest way to apply ZNE without managing the folding and extrapolation yourself.
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import EstimatorV2, QiskitRuntimeService, Options
from qiskit_ibm_runtime.fake_provider import FakeAlmadenV2
# Set up a fake backend for demonstration
backend = FakeAlmadenV2()
# Transpile the circuit for the target backend
pm = generate_preset_pass_manager(optimization_level=1, backend=backend)
qc = QuantumCircuit(1)
qc.h(0)
qc.rz(3.14159 / 3, 0)
qc.h(0)
isa_circuit = pm.run(qc)
# Define the observable
observable = SparsePauliOp("Z")
# Run WITHOUT error mitigation (resilience_level=0)
estimator_raw = EstimatorV2(mode=backend)
estimator_raw.options.resilience_level = 0
job_raw = estimator_raw.run([(isa_circuit, observable)])
raw_result = job_raw.result()[0].data.evs
# Run WITH ZNE (resilience_level=2)
estimator_zne = EstimatorV2(mode=backend)
estimator_zne.options.resilience_level = 2
job_zne = estimator_zne.run([(isa_circuit, observable)])
zne_result = job_zne.result()[0].data.evs
print(f"Raw result (resilience=0): {raw_result:.4f}")
print(f"ZNE result (resilience=2): {zne_result:.4f}")
print(f"Ideal result: 0.5000")
Note: resilience_level=1 applies measurement error mitigation only. resilience_level=2 adds ZNE on top of that. The Qiskit Runtime implementation uses digital ZNE (gate twirling and noise amplification) under the hood.
ZNE with Noise Characterization
If you have noise characterization data from randomized benchmarking (RB), you can use it to predict ZNE performance and validate your results.
Randomized benchmarking measures the average error rate per Clifford gate, typically denoted r. Under depolarizing noise, the expectation value of an n-qubit Pauli observable after d gates decays as:
E(d) ≈ E_ideal * (1 - 2r)^d
For a circuit with d gates at scale factor c, the effective depth is c * d, so:
E(c) ≈ E_ideal * (1 - 2r)^(c*d)
This is an exponential function of c. After ZNE with polynomial Richardson extrapolation of order n, the residual bias scales as the (n+1)th derivative of the noise function evaluated at c = 0, which for exponential decay is:
bias_ZNE ≈ E_ideal * [d * ln(1 - 2r)]^(n+1) / (n+1)!
import numpy as np
# Hypothetical RB data
rb_error_rate = 0.005 # 0.5% error per Clifford gate
d = 10 # Circuit depth (gate count)
E_ideal = 0.5 # Ideal expectation value
# Predicted noisy result (no mitigation)
noisy_prediction = E_ideal * (1 - 2 * rb_error_rate) ** d
print(f"Predicted noisy value: {noisy_prediction:.4f}")
# Predicted ZNE residual bias (Richardson order n=2, using 3 scale factors)
n = 2
log_factor = d * np.log(1 - 2 * rb_error_rate)
bias_zne = E_ideal * abs(log_factor) ** (n + 1) / np.math.factorial(n + 1)
print(f"Predicted ZNE bias: {bias_zne:.6f}")
print(f"Predicted ZNE result: {E_ideal - bias_zne:.4f}")
# Verify against actual simulation
import cirq
from mitiq import zne
q = cirq.LineQubit(0)
circuit = cirq.Circuit([cirq.H(q), cirq.rz(rads=np.pi / 3)(q), cirq.H(q)])
def executor_rb(circuit: cirq.Circuit) -> float:
noise_model = cirq.ConstantQubitNoiseModel(cirq.depolarize(p=rb_error_rate))
sim = cirq.DensityMatrixSimulator(noise=noise_model)
result = sim.simulate(circuit)
rho = result.final_density_matrix
return float(rho[0, 0].real - rho[1, 1].real)
factory = zne.RichardsonFactory(scale_factors=[1, 2, 3])
actual_zne = zne.execute_with_zne(circuit, executor_rb, factory=factory)
actual_noisy = executor_rb(circuit)
print(f"\nActual noisy value: {actual_noisy:.4f}")
print(f"Actual ZNE value: {actual_zne:.4f}")
print(f"Actual ZNE error: {abs(actual_zne - E_ideal):.6f}")
When the RB-predicted bias and the actual ZNE error are in the same ballpark, you can be confident that the noise is well-characterized and ZNE is working as expected. A large discrepancy indicates that the noise model is more complex than simple depolarizing (for example, coherent errors or correlated noise).
PEC Comparison: When to Switch from ZNE
Probabilistic Error Cancellation (PEC) is another error mitigation technique available in Mitiq. While ZNE extrapolates to zero noise, PEC directly inverts the noise channel by sampling from a quasi-probability distribution of circuits. Here is a side-by-side comparison.
import cirq
import numpy as np
from mitiq import zne, pec
from mitiq.pec.representations import represent_operations_in_circuit_with_local_depolarizing_noise
q = cirq.LineQubit(0)
circuit = cirq.Circuit([
cirq.H(q),
cirq.rz(rads=np.pi / 3)(q),
cirq.H(q),
])
ideal_value = 0.5
noise_level = 0.02
def noisy_exec(circuit: cirq.Circuit) -> float:
noise_model = cirq.ConstantQubitNoiseModel(cirq.depolarize(p=noise_level))
sim = cirq.DensityMatrixSimulator(noise=noise_model)
result = sim.simulate(circuit)
rho = result.final_density_matrix
return float(rho[0, 0].real - rho[1, 1].real)
# ZNE
factory = zne.RichardsonFactory(scale_factors=[1, 2, 3])
zne_result = zne.execute_with_zne(circuit, noisy_exec, factory=factory)
# PEC requires a noise representation
reps = represent_operations_in_circuit_with_local_depolarizing_noise(
circuit, noise_level=noise_level
)
pec_result = pec.execute_with_pec(
circuit, noisy_exec,
representations=reps,
num_samples=1000,
random_state=42,
)
noisy_result = noisy_exec(circuit)
print(f"{'Method':<12} {'Result':>8} {'Error':>8} {'Evals':>8}")
print(f"{'Noisy':<12} {noisy_result:>8.4f} {abs(noisy_result - ideal_value):>8.4f} {'1':>8}")
print(f"{'ZNE':<12} {zne_result:>8.4f} {abs(zne_result - ideal_value):>8.4f} {'3':>8}")
print(f"{'PEC':<12} {pec_result:>8.4f} {abs(pec_result - ideal_value):>8.4f} {'1000':>8}")
The practical tradeoffs:
- ZNE requires ~3x the circuit evaluations of an unmitigated run (one per scale factor). The overhead is fixed regardless of circuit size. The accuracy is limited by the polynomial noise model assumption.
- PEC requires gamma^2 * N samples, where gamma is the one-norm of the quasi-probability representation and grows exponentially with the number of noisy gates and the noise rate. For a circuit with n two-qubit gates at error rate p, gamma ≈ (1 + 4p/3)^n. At n = 10 gates and p = 0.01, gamma ≈ 1.14, so the overhead is about 1.3x. At n = 50 gates, gamma ≈ 1.94, and the overhead is 3.8x. At n = 100 gates, gamma ≈ 3.77, and the overhead is 14.2x.
Practical guideline: for circuits with fewer than 20 two-qubit gates and error rates below 1%, PEC provides more accurate results and the sampling overhead is manageable. For deeper circuits or higher noise rates, ZNE is the more practical choice because its overhead does not scale with circuit size.
ZNE Failure Modes
ZNE is not universally reliable. Understanding when it fails is as important as knowing how to use it.
Failure Mode 1: Coherent (Systematic) Errors
ZNE assumes that noise behaves monotonically with scale factor: more noise should push the expectation value smoothly toward the maximally mixed value. Coherent errors (systematic over-rotations or under-rotations) violate this assumption because they cause oscillatory behavior.
import cirq
import numpy as np
from mitiq import zne
q = cirq.LineQubit(0)
# Circuit with a systematic over-rotation error
# Instead of Rz(pi/3), the hardware applies Rz(pi/3 + 0.1)
epsilon = 0.1 # Systematic over-rotation
circuit_with_coherent_error = cirq.Circuit([
cirq.H(q),
cirq.rz(rads=np.pi / 3 + epsilon)(q),
cirq.H(q),
])
def coherent_executor(circuit: cirq.Circuit) -> float:
"""Simulate with ONLY coherent errors (no stochastic noise)."""
# The coherent error is baked into the circuit above.
# We simulate without additional noise.
sim = cirq.Simulator()
result = sim.simulate(circuit)
state = result.final_state_vector
return float((abs(state[0])**2 - abs(state[1])**2).real)
ideal_value = np.cos(np.pi / 3) # 0.5
# Under gate folding, the coherent error also gets folded:
# G G† G applies the over-rotation 3 times, giving 3*epsilon
# This means the error at scale c is cos(pi/3 + c*epsilon), not
# a polynomial function of c.
print("Scale factor -> Expectation value:")
for scale in [1, 2, 3, 4, 5]:
# At scale c, the over-rotation is c * epsilon
value = np.cos(np.pi / 3 + scale * epsilon)
print(f" c={scale}: {value:.4f}")
# ZNE extrapolation
factory = zne.RichardsonFactory(scale_factors=[1, 2, 3])
# For this to work with mitiq, we need a proper executor
def coherent_exec_for_mitiq(circuit: cirq.Circuit) -> float:
sim = cirq.Simulator()
result = sim.simulate(circuit)
state = result.final_state_vector
return float((abs(state[0])**2 - abs(state[1])**2).real)
zne_result = zne.execute_with_zne(
circuit_with_coherent_error, coherent_exec_for_mitiq, factory=factory
)
print(f"\nIdeal value: {ideal_value:.4f}")
print(f"Noisy value: {np.cos(np.pi / 3 + epsilon):.4f}")
print(f"ZNE result: {zne_result:.4f}")
print(f"ZNE error: {abs(zne_result - ideal_value):.4f}")
Because cos(pi/3 + c * epsilon) is not a low-degree polynomial in c, the Richardson extrapolation does not converge to the correct value. The extrapolation can even overshoot, giving a result further from the ideal than the raw noisy value.
Failure Mode 2: Noise Floor Below Shot Noise
ZNE relies on detecting a signal in how the expectation value changes across scale factors. If the noise contribution per gate is smaller than the shot noise, the differences between scale factors are indistinguishable from random fluctuations, and the extrapolation becomes unreliable.
import numpy as np
# Scenario: 50-gate circuit, 1% depolarizing noise per gate
d = 50
p = 0.01
n_shots = 1000
# Noisy expectation value (exponential decay)
E_ideal = 1.0
E_noisy = E_ideal * (1 - 2*p)**d
print(f"Noisy expectation value: {E_noisy:.4f}")
# Shot noise (standard error of the mean)
# For a ±1 observable, variance per shot ≈ 1
shot_noise = 1.0 / np.sqrt(n_shots)
print(f"Shot noise (1-sigma): {shot_noise:.4f}")
# Signal available for ZNE: difference between scale 1 and scale 3
E_scale3 = E_ideal * (1 - 2*p)**(3*d)
noise_signal = abs(E_noisy - E_scale3)
print(f"Signal (E(1) - E(3)): {noise_signal:.4f}")
# ZNE reliability criterion
print(f"\nReliability check:")
print(f" noise_signal / shot_noise = {noise_signal / shot_noise:.2f}")
if noise_signal > 2 * shot_noise:
print(f" ZNE is likely reliable (ratio > 2)")
else:
print(f" ZNE is UNRELIABLE (ratio < 2)")
print(f" Need at least {int((2 / (noise_signal / 1.0))**2)} shots per scale point")
For this 50-gate circuit, the noisy value is very close to zero (approximately 0.36), and at scale factor 3 it drops even further (approximately 0.005). The difference is detectable with 1000 shots. But if the circuit were deeper (say, 100 gates), the noisy value would be ~0.13 and the scale-3 value would be ~0.0002, making the signal-to-noise ratio too small for reliable extrapolation.
Failure Mode 3: Extrapolation Amplifies Shot Noise
Even when the noise signal is detectable, Richardson extrapolation amplifies the shot noise by the sum of absolute values of its coefficients. With coefficients [3, -3, 1], the amplification factor is |3| + |-3| + |1| = 7. This means you need roughly 7x more shots per scale factor to achieve the same precision as an unmitigated estimate.
For very deep circuits where the noisy expectation value is already close to the maximally mixed state, the ZNE “correction” can be dominated by amplified shot noise, making the mitigated result worse than the raw result.
import numpy as np
# Compare ZNE standard error to raw standard error
n_shots_per_scale = 1000
n_scale_factors = 3
richardson_coeffs = np.array([3, -3, 1])
# Raw standard error (single evaluation)
sigma_raw = 1.0 / np.sqrt(n_shots_per_scale)
# ZNE standard error
sigma_zne = np.sqrt(np.sum(richardson_coeffs**2)) / np.sqrt(n_shots_per_scale)
print(f"Raw std error: {sigma_raw:.4f}")
print(f"ZNE std error: {sigma_zne:.4f}")
print(f"Amplification: {sigma_zne / sigma_raw:.1f}x")
print(f"Total shots used: {n_shots_per_scale * n_scale_factors}")
Common Mistakes
Mistake 1: ZNE with Mid-Circuit Measurements
ZNE assumes the circuit is a single coherent unitary block. The noise scaling argument (replacing G with G G† G preserves the logical operation while tripling the noise) relies on the circuit being a sequence of unitary gates. Mid-circuit measurements collapse the quantum state, breaking the noise scaling assumption: the post-measurement state does not scale with lambda in the same way as a unitary channel, and the folded circuit produces a logically different operation.
If your circuit has mid-circuit measurements, either: (a) defer all measurements to the end, (b) split the circuit into coherent blocks and apply ZNE to each block separately, or (c) use a different mitigation technique like PEC.
Mistake 2: Insufficient Shots per Scale Level
The standard error of the Richardson estimate is:
sigma_ZNE = sqrt(sum_i c_i^2 * sigma_i^2 / N_i)
For scale factors [1, 2, 3] with Richardson coefficients [3, -3, 1], the sum of absolute coefficients is |3| + |-3| + |1| = 7. With equal shot allocation, the combined amplification factor is sqrt(9 + 9 + 1) ≈ 4.4. This means you need about 19x more total shots (4.4^2 * 1, spread across 3 scale factors) to achieve the same standard error as a single unmitigated evaluation.
In practice, 1000-8000 shots per scale point is typical for 1-5 qubit circuits. For larger circuits or higher precision requirements, plan for 10,000+ shots per scale point.
Mistake 3: Confusing Noise Parameter with Scale Factor
The scale factor c means the circuit is c times longer. Under Markovian noise, the effective noise parameter lambda becomes c * lambda_0 (linear in c, not c^2 or exponential). The polynomial noise model is:
E(c) = E_ideal + A₁c + A₂c² + ...
Here, c is the scale factor (input variable), not the noise parameter itself. The coefficients A_k encode the noise model. A common error is to write E(c) = E_ideal * exp(-c * lambda), which conflates the scale factor with the noise parameter. The correct exponential model is E(c) = E_ideal * exp(-c * d * lambda), where d is the circuit depth and lambda is the per-gate noise rate.
Mistake 4: Applying ZNE to Sampling Tasks
ZNE only works for expectation values (scalar quantities). It cannot improve full output probability distributions or individual bitstring probabilities. If your algorithm requires sampling specific bitstrings (like Shor’s algorithm or random circuit sampling), ZNE does not apply.
For sampling tasks, consider Probabilistic Error Cancellation (PEC), which directly corrects the output distribution, or measurement error mitigation, which corrects the classical readout step.
Mistake 5: Trusting a Single ZNE Estimate
ZNE estimates have statistical variance. A single run might happen to land close to the ideal value (or far from it) by chance. Always run multiple independent ZNE estimates and compute their spread.
import cirq
import numpy as np
from mitiq import zne
q = cirq.LineQubit(0)
circuit = cirq.Circuit([
cirq.H(q),
cirq.rz(rads=np.pi / 3)(q),
cirq.H(q),
])
def shot_based_executor(circuit: cirq.Circuit) -> float:
"""Executor with finite shot noise (1000 shots)."""
noise_model = cirq.ConstantQubitNoiseModel(cirq.depolarize(p=0.02))
sim = cirq.DensityMatrixSimulator(noise=noise_model)
# Use sampling instead of exact density matrix
circuit_with_measure = circuit + cirq.measure(circuit.all_qubits(), key='m')
result = sim.run(circuit_with_measure, repetitions=1000)
counts = result.histogram(key='m')
total = sum(counts.values())
p0 = counts.get(0, 0) / total
p1 = counts.get(1, 0) / total
return p0 - p1 # <Z>
# Run 5 independent ZNE estimates
factory = zne.RichardsonFactory(scale_factors=[1, 2, 3])
estimates = []
for i in range(5):
val = zne.execute_with_zne(circuit, shot_based_executor, factory=factory)
estimates.append(val)
print(f" Run {i+1}: {val:.4f}")
mean_val = np.mean(estimates)
std_val = np.std(estimates, ddof=1)
print(f"\nMean: {mean_val:.4f}")
print(f"Std dev: {std_val:.4f}")
print(f"Ideal: 0.5000")
# Check if the spread is smaller than the expected mitigation improvement
expected_improvement = 0.07 # From the unmitigated example
if std_val > expected_improvement:
print("WARNING: ZNE variance exceeds the expected improvement.")
print("Increase your shot count before trusting these results.")
else:
print("ZNE variance is within acceptable bounds.")
If the standard deviation across runs is larger than the expected mitigation improvement (~0.07 in our running example), the shot count is insufficient and ZNE is not providing reliable results.
What to Try Next
- Apply ZNE to a multi-qubit circuit and inspect how mitigation quality changes with circuit depth
- Try PEC (
mitiq.pec.execute_with_pec) when you have a calibrated noise model - Experiment with different folding strategies on real hardware to find which works best for your device’s noise profile
- Combine ZNE with measurement error mitigation for additional improvement
- See the Mitiq Reference for a full list of supported techniques and backends
Was this tutorial helpful?