Quantum noise and error mitigation
Toshinari Itoko (28 June 2024)
Download the pdf of the original lecture. Note that some code snippets might become deprecated since these are static images.
Approximate QPU time to run this experiment is 1 m 40 s.
1. Introduction
Throughout this lesson, we will examine noise and how it can be mitigated on quantum computers. We will begin by looking at the effects of noise using a simulator that can simulate noise in a few ways, including using noise profiles from real quantum computers. Then we will move on to real quantum computers, in which noise is inherent. We will look at the effects of error mitigation, including combinations of things like zero-noise extrapolation (ZNE) and gate-twirling.
We will start by loading some packages.
# !pip install qiskit qiskit_aer qiskit_ibm_runtime
# !pip install jupyter
# !pip install matplotlib pylatexenc
import qiskit
qiskit.__version__
Output:
'2.0.2'
import qiskit_aer
qiskit_aer.__version__
Output:
'0.17.1'
import qiskit_ibm_runtime
qiskit_ibm_runtime.__version__
Output:
'0.40.1'
2. Noisy simulation without error mitigation
Qiskit Aer is a classical simulator for quantum computing. It can simulate not only ideal execution but also noisy execution of quantum circuits. This notebook demonstrates how to run noisy simulation using Qiskit Aer:
- Build a noise model
- Build a noisy sampler (simulator) with the noise model
- Run a quantum circuit on the noisy sampler
noise_model = NoiseModel()
...
noisy_sampler = Sampler(options={"backend_options": {"noise_model": noise_model}})
job = noisy_sampler.run([circuit])
2.1 Build a test circuit
We consider toy 1-qubit circuits which just repeat X gates d
times (d
=0 ... 100) and measure the Z
observable.
from qiskit.circuit import QuantumCircuit
MAX_DEPTH = 100
circuits = []
for d in range(MAX_DEPTH + 1):
circ = QuantumCircuit(1)
for _ in range(d):
circ.x(0)
circ.barrier(0)
circ.measure_all()
circuits.append(circ)
display(circuits[3].draw(output="mpl"))
Output:

from qiskit.quantum_info import SparsePauliOp
obs = SparsePauliOp.from_list([("Z", 1.0)])
obs
Output:
SparsePauliOp(['Z'],
coeffs=[1.+0.j])
2.2 Build a noise model
To do noisy simulation, we need to specify NoiseModel
. We show how to build NoiseModel
in this section.
We first need to define quantum (or readout) errors to add to a noise model.
from qiskit_aer.noise.errors import (
coherent_unitary_error,
amplitude_damping_error,
ReadoutError,
)
from qiskit.circuit.library import RXGate
# Coherent (unitary) error: Over X-rotation error
# https://qiskit.github.io/qiskit-aer/stubs/qiskit_aer.noise.coherent_unitary_error.html#qiskit_aer.noise.coherent_unitary_error
OVER_ROTATION_ANGLE = 0.05
coherent_error = coherent_unitary_error(RXGate(OVER_ROTATION_ANGLE).to_matrix())
# Incoherent error: Amplitude dumping error
# https://qiskit.github.io/qiskit-aer/stubs/qiskit_aer.noise.amplitude_damping_error.html#qiskit_aer.noise.amplitude_damping_error
AMPLITUDE_DAMPING_PARAM = 0.02 # in [0, 1] (0: no error)
incoherent_error = amplitude_damping_error(AMPLITUDE_DAMPING_PARAM)
# Readout (measurement) error: Readout error
# https://qiskit.github.io/qiskit-aer/stubs/qiskit_aer.noise.ReadoutError.html#qiskit_aer.noise.ReadoutError
PREP0_MEAS1 = 0.03 # P(1|0): Probability of preparing 0 and measuring 1
PREP1_MEAS0 = 0.08 # P(0|1): Probability of preparing 1 and measuring 0
readout_error = ReadoutError(
[[1 - PREP0_MEAS1, PREP0_MEAS1], [PREP1_MEAS0, 1 - PREP1_MEAS0]]
)
from qiskit_aer.noise import NoiseModel
noise_model = NoiseModel()
noise_model.add_quantum_error(coherent_error.compose(incoherent_error), "x", (0,))
noise_model.add_readout_error(readout_error, (0,))
2.3 Build a noisy sampler with the noise model
from qiskit_aer.primitives import SamplerV2 as Sampler
noisy_sampler = Sampler(options={"backend_options": {"noise_model": noise_model}})
2.4 Run quantum circuits on the noisy sampler
job = noisy_sampler.run(circuits, shots=400)
result = job.result()
result[0].data.meas.get_counts()
Output:
{'0': 389, '1': 11}
2.5 Plot results
import matplotlib.pyplot as plt
plt.title("Noisy simulation")
ds = list(range(MAX_DEPTH + 1))
plt.plot(
ds,
[result[d].data.meas.expectation_values(["Z"]) for d in ds],
color="gray",
linestyle="-",
)
plt.scatter(ds, [result[d].data.meas.expectation_values(["Z"]) for d in ds], marker="o")
plt.hlines(0, xmin=0, xmax=MAX_DEPTH, colors="black")
plt.ylim(-1, 1)
plt.xlabel("Circuit depth")
plt.ylabel("Measured <Z>")
plt.show()
2.6 Ideal simulation
ideal_sampler = Sampler()
job_ideal = ideal_sampler.run(circuits)
result_ideal = job_ideal.result()
plt.title("Ideal simulation")
ds = list(range(MAX_DEPTH + 1))
plt.plot(
ds,
[result_ideal[d].data.meas.expectation_values(["Z"]) for d in ds],
color="gray",
linestyle="-",
)
plt.scatter(
ds, [result_ideal[d].data.meas.expectation_values(["Z"]) for d in ds], marker="o"
)
plt.hlines(0, xmin=0, xmax=MAX_DEPTH, colors="black")
plt.xlabel("Circuit depth")
plt.ylabel("Measured <Z>")
plt.show()
Output:

2.7 Exercise
By tweaking the code below,
- Try 25x number of shots (= 10_000 shots) and ensure that a smoother plot is obtained
- Change noise parameters (OVER_ROTATION_ANGLE, AMPLITUDE_DAMPING_PARAM, PREP0_MEAS1, or PREP1_MEAS0) and see how the plot changes
OVER_ROTATION_ANGLE = 0.05
coherent_error = coherent_unitary_error(RXGate(OVER_ROTATION_ANGLE).to_matrix())
AMPLITUDE_DAMPING_PARAM = 0.02 # in [0, 1] (0: no error)
incoherent_error = amplitude_damping_error(AMPLITUDE_DAMPING_PARAM)
PREP0_MEAS1 = 0.1 # P(1|0): Probability of preparing 0 and measuring 1
PREP1_MEAS0 = 0.05 # P(0|1): Probability of preparing 1 and measuring 0
readout_error = ReadoutError(
[[1 - PREP0_MEAS1, PREP0_MEAS1], [PREP1_MEAS0, 1 - PREP1_MEAS0]]
)
noise_model = NoiseModel()
noise_model.add_quantum_error(coherent_error.compose(incoherent_error), "x", (0,))
noise_model.add_readout_error(readout_error, (0,))
options = {
"backend_options": {"noise_model": noise_model},
}
noisy_sampler = Sampler(options=options)
job = noisy_sampler.run(circuits, shots=400)
result = job.result()
plt.title("Noisy simulation")
ds = list(range(MAX_DEPTH + 1))
plt.plot(
ds,
[result[d].data.meas.expectation_values(["Z"]) for d in ds],
marker="o",
linestyle="-",
)
plt.hlines(0, xmin=0, xmax=MAX_DEPTH, colors="black")
plt.ylim(-1, 1)
plt.xlabel("Depth")
plt.ylabel("Measured <Z>")
plt.show()
Output:

2.8 More realistic noisy simulation
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import SamplerV2 as Sampler, QiskitRuntimeService
service = QiskitRuntimeService()
real_backend = service.least_busy(
operational=True, simulator=False, min_num_qubits=127
) # Eagle
Output:
<IBMBackend('ibm_strasbourg')>
aer = AerSimulator.from_backend(real_backend)
noisy_sampler = Sampler(mode=aer)
job = noisy_sampler.run(circuits)
result = job.result()
plt.title("Noisy simulation with noise model from real backend")
ds = list(range(MAX_DEPTH + 1))
plt.plot(
ds,
[result[d].data.meas.expectation_values(["Z"]) for d in ds],
marker="o",
linestyle="-",
)
plt.hlines(0, xmin=0, xmax=MAX_DEPTH, colors="black")
plt.ylim(-1, 1)
plt.xlabel("Depth")
plt.ylabel("Measured <Z>")
plt.show()
Output:

3. Real quantum computation with error mitigation
In this part, we demonstrate how to obtain error mitigated results (expectation values) using Qiskit Estimator. We consider 6-qubit Trotterized circuits for simulating the time evolution of one dimensional Ising model and see how the error scales with respect to the number of time steps.
backend = service.least_busy(
operational=True, simulator=False, min_num_qubits=127
) # Eagle
backend
Output:
<IBMBackend('ibm_strasbourg')>
NUM_QUBITS = 6
NUM_TIME_STEPS = list(range(8))
RX_ANGLE = 0.1
RZZ_ANGLE = 0.1
3.1 Build circuits
# Build circuits with different number of time steps
circuits = []
for n_steps in NUM_TIME_STEPS:
circ = QuantumCircuit(NUM_QUBITS)
for i in range(n_steps):
# rx layer
for q in range(NUM_QUBITS):
circ.rx(RX_ANGLE, q)
# 1st rzz layer
for q in range(1, NUM_QUBITS - 1, 2):
circ.rzz(RZZ_ANGLE, q, q + 1)
# 2nd rzz layer
for q in range(0, NUM_QUBITS - 1, 2):
circ.rzz(RZZ_ANGLE, q, q + 1)
circ.barrier() # need not to optimize the circuit
# Uncompute stage
for i in range(n_steps):
for q in range(0, NUM_QUBITS - 1, 2):
circ.rzz(-RZZ_ANGLE, q, q + 1)
for q in range(1, NUM_QUBITS - 1, 2):
circ.rzz(-RZZ_ANGLE, q, q + 1)
for q in range(NUM_QUBITS):
circ.rx(-RX_ANGLE, q)
circuits.append(circ)
To know the ideal output in advance, we use compute-uncompute circuits that consist of a first stage where the original circuit is applied, and a second stage where it is reversed . Note that the ideal outcome of such circuits will trivially be the input state , which has the trivial expectation values for any Pauli observables, e.g. .
# Print the circuit with 2 time steps
circuits[2].draw(output="mpl")
Output:

Note: As shown above, the circuit with time steps will have two-qubit gate layers.
obs = SparsePauliOp.from_sparse_list([("Z", [0], 1.0)], num_qubits=NUM_QUBITS)
obs
Output:
SparsePauliOp(['IIIIIZ'],
coeffs=[1.+0.j])
3.2 Transpile the circuits
We transpile the circuits for the backend with optimization (optimization_level=1
).
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
pm = generate_preset_pass_manager(optimization_level=1, backend=backend)
isa_circuits = pm.run(circuits)
display(isa_circuits[2].draw("mpl", idle_wires=False, fold=-1))
Output:

3.3 Execute using Estimator (with different resilience levels)
Setting the resilienece level (estimator.options.resilience_level
) is the easiest way to apply error mitigation when using Qiskit Estimator. Estimator supports the following resilience levels (as of 2024/06/28). See more details in the Configure error mitigation guide.

from qiskit_ibm_runtime import Batch
from qiskit_ibm_runtime import EstimatorV2 as Estimator
jobs = []
job_ids = []
with Batch(backend=backend):
for resilience_level in [0, 1, 2]:
estimator = Estimator()
estimator.options.resilience_level = resilience_level
job = estimator.run(
[(circ, obs.apply_layout(circ.layout)) for circ in isa_circuits]
)
job_ids.append(job.job_id())
print(f"Job ID (rl={resilience_level}): {job.job_id()}")
jobs.append(job)
Output:
Job ID (rl=0): d146vcnmya70008emprg
Job ID (rl=1): d146vdnqf56g0081sva0
Job ID (rl=2): d146ven5z6q00087c61g
# check job status
for job in jobs:
print(job.status())
Output:
DONE
DONE
DONE
# REPLACE WITH YOUR OWN JOB IDS
jobs = [service.job(job_id) for job_id in job_ids]
# Get results
results = [job.result() for job in jobs]
3.4 Plot results
plt.title("Error mitigation with different resilience levels")
labels = ["0 (No mitigation)", "1 (TREX)", "2 (ZNE + Gate twirling)"]
steps = NUM_TIME_STEPS
for result, label in zip(results, labels):
plt.errorbar(
x=steps,
y=[result[s].data.evs for s in steps],
yerr=[result[s].data.stds for s in steps],
marker="o",
linestyle="-",
capsize=4,
label=label,
)
plt.hlines(
1.0, min(steps), max(steps), linestyle="dashed", label="Ideal", colors="black"
)
plt.xlabel("Time steps")
plt.ylabel("Mitigated <IIIIIZ>")
plt.legend()
plt.show()
Output:

4. (Optional) Customize error mitigation options
We can customize the application of error mitigation techniques via options as shown below.
# TREX
estimator.options.twirling.enable_measure = True
estimator.options.twirling.num_randomizations = "auto"
estimator.options.twirling.shots_per_randomization = "auto"
# Gate twirling
estimator.options.twirling.enable_gates = True
# ZNE
estimator.options.resilience.zne_mitigation = True
estimator.options.resilience.zne.noise_factors = [1, 3, 5]
estimator.options.resilience.zne.extrapolator = ("exponential", "linear")
# Dynamical decoupling
estimator.options.dynamical_decoupling.enable = True # Default: False
estimator.options.dynamical_decoupling.sequence_type = "XX"
# Other options
estimator.options.default_shots = 10_000
See the following guides and API reference for the details of error mitigation options.