Quantum Error Mitigation: ZNE vs PEC vs Clifford Data Regression

Error mitigation is not error correction. Error correction encodes logical qubits in physical qubits and actively detects and corrects errors; it requires hundreds to thousands of physical qubits per logical qubit and hardware we do not yet have at scale. Error mitigation is a different approach for NISQ-era hardware: run noisy circuits, then use classical post-processing to recover a less-noisy expectation value. No extra qubits required.

The three main mitigation methods (zero-noise extrapolation, ZNE; probabilistic error cancellation, PEC; and Clifford data regression, CDR) have different assumptions, costs, and failure modes. Understanding when each works (and when it doesn’t) is essential for getting useful results from real hardware today.

Setup

pip install mitiq qiskit qiskit-aer

import numpy as np
import mitiq
from mitiq import zne, pec, cdr
from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error

The test circuit and noisy executor

All three methods need a function that takes a circuit and returns an expectation value. We’ll use a simple H-CX circuit with a depolarizing noise model as a controlled test environment.

from qiskit import QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, depolarizing_error

def make_bell_circuit():
    qc = QuantumCircuit(2)
    qc.h(0)
    qc.cx(0, 1)
    return qc

def make_noise_model(error_rate: float):
    noise = NoiseModel()
    noise.add_all_qubit_quantum_error(depolarizing_error(error_rate, 1), ['h'])
    noise.add_all_qubit_quantum_error(depolarizing_error(error_rate, 2), ['cx'])
    return noise

# The noisy executor: runs circuit, returns <ZZ> expectation value
def noisy_executor(circuit, error_rate=0.02):
    """Execute circuit on noisy Aer simulator, return <Z0 Z1>."""
    from qiskit import transpile
    from qiskit.quantum_info import SparsePauliOp
    
    sim = AerSimulator(noise_model=make_noise_model(error_rate))
    qc = circuit.copy()
    qc.save_statevector()
    
    # For Mitiq compatibility, work with shot-based sampling
    qc_meas = circuit.copy()
    qc_meas.measure_all()
    result = sim.run(transpile(qc_meas, sim), shots=4096).result()
    counts = result.get_counts()
    
    # <ZZ> = P(00) + P(11) - P(01) - P(10)
    total = sum(counts.values())
    zz = 0.0
    for bitstring, count in counts.items():
        # Qiskit bit order: rightmost = qubit 0
        q0 = int(bitstring[-1])
        q1 = int(bitstring[-2]) if len(bitstring) > 1 else 0
        sign = 1 if (q0 == q1) else -1
        zz += sign * count / total
    return zz

# Ideal value for Bell state: <ZZ> = 1.0
# (Both qubits always correlated: 00 or 11)
ideal = 1.0
noisy = noisy_executor(make_bell_circuit(), error_rate=0.02)
print(f"Ideal <ZZ>: {ideal:.4f}")
print(f"Noisy <ZZ>: {noisy:.4f}")
print(f"Error: {abs(ideal - noisy):.4f}")
# Example output:
# Ideal <ZZ>: 1.0000
# Noisy <ZZ>: ~0.97  (varies with shots)
# Error: ~0.03

Method 1: Zero-Noise Extrapolation (ZNE)

Idea. Run the circuit at several artificially higher noise levels, then extrapolate back to zero noise. The intuition: if you know how the result degrades as noise increases, you can reverse-engineer the noiseless answer.

Noise scaling. The standard approach is gate-level folding: replace each gate G with G * G† * G (which is equivalent to G but adds noise). Folding by factor 3 triples the gate count and roughly triples the noise. Mitiq handles this automatically.

Extrapolation. Fit a polynomial (usually linear or Richardson extrapolation) to the (noise_level, expectation_value) pairs and evaluate at zero noise.

from mitiq import zne

# Define the executor compatible with Mitiq (takes a mitiq.Circuit)
def mitiq_executor(circuit):
    """Mitiq-compatible executor: converts Mitiq circuit to Qiskit and runs."""
    # Mitiq uses Cirq circuits internally; use Qiskit frontend
    from mitiq.interface.mitiq_qiskit import qiskit_utils
    # For simplicity, directly use the Qiskit circuit version
    return noisy_executor(make_bell_circuit())

# ZNE with linear extrapolation
# scale_factors: run at noise 1x, 2x, 3x then extrapolate to 0x
factory = zne.LinearFactory(scale_factors=[1.0, 2.0, 3.0])

# Manual demonstration (without Mitiq's full interface for clarity):
scale_factors = [1.0, 2.0, 3.0]
results = []
for scale in scale_factors:
    # Approximate higher noise by increasing error rate
    scaled_noise = 0.02 * scale
    result = noisy_executor(make_bell_circuit(), error_rate=scaled_noise)
    results.append(result)
    print(f"Scale {scale:.1f}x | noise {scaled_noise:.3f} | <ZZ> = {result:.4f}")

# Linear extrapolation to zero noise
coeffs = np.polyfit(scale_factors, results, deg=1)
zne_estimate = np.polyval(coeffs, 0.0)  # evaluate at scale = 0
print(f"\nZNE estimate: {zne_estimate:.4f} (ideal: {ideal:.4f})")
print(f"ZNE error:    {abs(ideal - zne_estimate):.4f}")
print(f"Raw error:    {abs(ideal - results[0]):.4f}")

When ZNE works well:

• Circuit depth is small to moderate (fewer gates = cleaner scaling behavior)
• Noise is approximately gate-level depolarizing (scales predictably with gate count)
• You can run 3-5x more shots than a single circuit run

When ZNE fails:

• Very deep circuits where noise is dominated by T1/T2 (amplitude damping), not gate errors. T1/T2 decay is not linear in gate count.
• When the true expectation value is near 0 and the extrapolated value has high uncertainty (extrapolation amplifies variance)
• When noise is non-Markovian (correlated over time)

Method 2: Probabilistic Error Cancellation (PEC)

Idea. Represent each noisy gate as a linear combination of implementable noisy operations. Then sample the circuit from this representation and take a quasi-probability-weighted average. This produces an unbiased estimate of the ideal circuit output.

Cost. PEC is unbiased but exponentially expensive in sampling: the number of shots needed scales as e^(2 * gamma * circuit_depth) where gamma is the one-norm of the quasi-probability representation. For practical circuits, PEC is typically only feasible for shallow circuits with well-characterized noise.

The sampling overhead. For a circuit with 10 two-qubit gates at 1% error rate each, the PEC overhead is approximately e^(2 * 0.01 * 10) ≈ 1.2x. For 100 gates at 1%, it’s e^2 ≈ 7x. For 100 gates at 5%, it’s e^10 ≈ 22,000x: completely impractical.

# PEC requires knowing the noise model in detail to construct quasi-probability representations.
# This simplified example shows the structure without full noise tomography.

from mitiq import pec

# PEC with Mitiq requires OperationRepresentations -- structured descriptions of how
# each gate behaves under noise and what implementable operations it decomposes into.
# For demonstration, we estimate the overhead:

def pec_sampling_overhead(n_two_qubit_gates: int, two_qubit_error_rate: float) -> float:
    """Estimate PEC shot overhead for a circuit."""
    # gamma per gate ~ error_rate for depolarizing noise
    gamma_total = n_two_qubit_gates * two_qubit_error_rate
    # PEC overhead ~ exp(2 * gamma_total)
    return np.exp(2 * gamma_total)

print("PEC sampling overhead estimates:")
configs = [
    (5, 0.01, "5 gates, 1% error"),
    (10, 0.01, "10 gates, 1% error"),
    (10, 0.05, "10 gates, 5% error"),
    (50, 0.01, "50 gates, 1% error"),
]
for n_gates, err, label in configs:
    overhead = pec_sampling_overhead(n_gates, err)
    print(f"  {label}: {overhead:.1f}x more shots")
# Output:
#   5 gates, 1% error: 1.1x more shots
#   10 gates, 1% error: 1.2x more shots
#   10 gates, 5% error: 2.7x more shots
#   50 gates, 1% error: 2.7x more shots

# Conclusion: PEC is practical only for shallow circuits with low error rates

When PEC works well:

• Shallow circuits (< 20 two-qubit gates)
• Well-characterized noise (you’ve done gate set tomography)
• Low error rates (< 1% two-qubit gate error)
• You need an unbiased estimate with controlled variance

When PEC fails:

• Deep circuits (overhead becomes astronomical)
• Unknown or complex noise (quasi-probability decomposition requires accurate noise model)
• When the overhead makes it more expensive than just running more shots

Method 3: Clifford Data Regression (CDR)

Idea. Clifford circuits (circuits using only H, CNOT, S gates) can be classically simulated efficiently. Near-Clifford circuits are circuits close to Clifford: they contain mostly Clifford gates with a few non-Clifford rotations.

CDR works by:

1. Generating a set of “training circuits” by replacing non-Clifford gates in your circuit with Clifford approximations
2. Running both the ideal (classically simulated) and noisy (hardware) versions of training circuits
3. Fitting a linear regression model from noisy to ideal outputs using training data
4. Applying the learned mapping to the noisy output of your actual circuit

# CDR demonstration: replace T-gate rotations with Clifford approximations
import numpy as np
from qiskit import QuantumCircuit

def make_near_clifford_circuit():
    """A circuit with some non-Clifford (RZ) rotations."""
    qc = QuantumCircuit(3)
    qc.h(0)
    qc.cx(0, 1)
    qc.rz(0.3, 0)   # non-Clifford: small rotation
    qc.cx(1, 2)
    qc.rz(0.15, 2)  # non-Clifford: small rotation
    return qc

def make_clifford_approximation(circuit):
    """Replace RZ rotations with nearby Clifford gates for training circuits."""
    # For demonstration, replace RZ with S gate (pi/2 rotation)
    qc_clifford = QuantumCircuit(circuit.num_qubits)
    for instruction in circuit.data:
        gate = instruction.operation
        qubits = instruction.qubits
        if gate.name == 'rz':
            qc_clifford.s(qubits[0])  # Replace RZ with S (Clifford)
        else:
            qc_clifford.append(gate, qubits)
    return qc_clifford

# Generate training pairs
def cdr_training(n_training=8, error_rate=0.02):
    """Collect (ideal, noisy) pairs from Clifford approximations."""
    base_circuit = make_near_clifford_circuit()
    pairs = []
    
    for i in range(n_training):
        training_circuit = make_clifford_approximation(base_circuit)
        
        # Ideal: simulate classically (error-free)
        from qiskit_aer import AerSimulator
        ideal_sim = AerSimulator()
        qc_meas = training_circuit.copy()
        qc_meas.measure_all()
        from qiskit import transpile
        ideal_result = ideal_sim.run(transpile(qc_meas, ideal_sim), shots=8192).result()
        ideal_counts = ideal_result.get_counts()
        # Compute <Z0> = P(q0=0) - P(q0=1)
        total = sum(ideal_counts.values())
        ideal_z0 = sum((1 if b[-1] == '0' else -1) * c / total for b, c in ideal_counts.items())
        
        # Noisy: run on noisy simulator
        noisy_sim = AerSimulator(noise_model=make_noise_model(error_rate))
        noisy_result = noisy_sim.run(transpile(qc_meas, noisy_sim), shots=8192).result()
        noisy_counts = noisy_result.get_counts()
        noisy_total = sum(noisy_counts.values())
        noisy_z0 = sum((1 if b[-1] == '0' else -1) * c / noisy_total for b, c in noisy_counts.items())
        
        pairs.append((noisy_z0, ideal_z0))
    
    return pairs

# Fit linear regression
pairs = cdr_training(n_training=10, error_rate=0.02)
noisy_vals = np.array([p[0] for p in pairs])
ideal_vals = np.array([p[1] for p in pairs])

# Linear fit: ideal = a * noisy + b
a, b = np.polyfit(noisy_vals, ideal_vals, deg=1)

# Apply to our actual circuit
actual_circuit = make_near_clifford_circuit()
qc_meas = actual_circuit.copy()
qc_meas.measure_all()
from qiskit import transpile
from qiskit_aer import AerSimulator

noisy_sim = AerSimulator(noise_model=make_noise_model(0.02))
result = noisy_sim.run(transpile(qc_meas, noisy_sim), shots=4096).result()
counts = result.get_counts()
total = sum(counts.values())
noisy_z0 = sum((1 if b[-1] == '0' else -1) * c / total for b, c in counts.items())

cdr_estimate = a * noisy_z0 + b

ideal_sim = AerSimulator()
result = ideal_sim.run(transpile(qc_meas, ideal_sim), shots=4096).result()
counts = result.get_counts()
total = sum(counts.values())
ideal_z0 = sum((1 if b[-1] == '0' else -1) * c / total for b, c in counts.items())

print(f"Ideal <Z0>:        {ideal_z0:.4f}")
print(f"Noisy <Z0>:        {noisy_z0:.4f}")
print(f"CDR estimate:      {cdr_estimate:.4f}")
print(f"Noisy error:       {abs(ideal_z0 - noisy_z0):.4f}")
print(f"CDR error:         {abs(ideal_z0 - cdr_estimate):.4f}")

When CDR works well:

• Your circuit is close to Clifford (most gates are Clifford, few non-Clifford rotations)
• You can generate many meaningful training circuits
• The noise model is consistent between training and target circuits (same hardware, same session)

When CDR fails:

• Circuits far from Clifford (many large-angle rotations; deep VQE circuits, for example)
• When training circuits are not structurally similar enough to the target to transfer the learned correction
• When noise varies significantly across the chip or across time

Comparison Summary

Method	Assumptions	Shot overhead	Best circuit type	Bias
ZNE	Noise scales predictably with gate count	3-5x	Moderate depth, depolarizing noise	Small
PEC	Accurate noise model known	Exponential in depth	Shallow, well-characterized	None (unbiased)
CDR	Circuit is near-Clifford	~10x (training)	Near-Clifford, small rotations	Model-dependent

Combining Methods

ZNE and CDR can be combined: use CDR to build the training data at multiple noise levels, then extrapolate. Mitiq calls this “variable noise CDR.” This hybrid can outperform either method alone when the circuit is near-Clifford and the noise scales somewhat predictably.

# Using Mitiq's built-in ZNE (works with Cirq or Qiskit circuits natively)
# Requires: pip install mitiq
# See: mitiq.zne.execute_with_zne for the full interface

# The simplest Mitiq ZNE usage with a Cirq circuit:
import cirq
import mitiq

# Create a simple Bell circuit in Cirq
q0, q1 = cirq.LineQubit.range(2)
bell = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1)])

def cirq_executor(circuit: cirq.Circuit) -> float:
    """Noisy Cirq executor returning <ZZ>."""
    result = cirq.DensityMatrixSimulator(noise=cirq.depolarize(0.02)).simulate(circuit)
    dm = result.final_density_matrix
    zz_op = cirq.tensor_product(cirq.unitary(cirq.Z), cirq.unitary(cirq.Z))
    return float(np.real(np.trace(dm @ zz_op)))

# ZNE with Richardson extrapolation
mitigated = mitiq.zne.execute_with_zne(
    bell,
    cirq_executor,
    factory=mitiq.zne.RichardsonFactory(scale_factors=[1, 2, 3]),
)
print(f"Mitigated <ZZ>: {mitigated:.4f} (ideal: 1.0)")
# Expected output: Mitigated <ZZ>: ~1.0000 (ideal: 1.0)

Practical Recommendation

Start with ZNE. It requires the fewest assumptions, works with existing noise models, and has a predictable implementation path. Move to CDR if your circuit is near-Clifford and you can generate good training data. Use PEC only if you have well-characterized noise and circuits shallow enough that the exponential overhead is manageable, typically fewer than 20 two-qubit gates at sub-1% error rates.

None of these methods replaces hardware improvement. They reduce the effective noise, but they cannot exceed the information-theoretic limits set by the shot count and the noise level. For precision work, treat mitigation as a way to extend the useful depth of near-term devices, not as a path to fault tolerance.