Combine error mitigation options with the Estimator primitive

Usage estimate: 7 minutes on a Heron r2 processor (NOTE: This is an estimate only. Your runtime might vary.)

Learning outcomes

We suggest that users are familiar with the following topics before going through this tutorial:

• The basics of dynamical decoupling, measurement error mitigation, gate twirling, and zero-noise extrapolation, as described in this guide.

Prerequisites

After going through this tutorial, users should understand:

• How the aforementioned error mitigation techniques are selectively implemented on hardware.
• How they compare in terms of their ability to mitigate hardware noise.

Background

This tutorial explores the error suppression and error mitigation options available with the Estimator primitive from Qiskit Runtime. This tutorial shows how to implement each of the follow methods individually:

• Dynamical decoupling
• Measurement error mitigation
• Gate twirling
• Zero-noise extrapolation (ZNE)

Note that an alternative to implementing these techniques individually is to implement them using a resilience level, whereby resilience_level takes values 0, 1, 2:

• 0 : No mitigation is implemented.
• 1 : Measurement error mitigation is implemented.
• 2 : Gate twirling, measurement error mitigation, and ZNE are implemented.

In this tutorial, you will construct a circuit and observable and submit jobs using the Estimator primitive using different combinations of error mitigation settings. Then, you will plot the results to observe the effects of the various settings. Most of the tutorial uses a 10-qubit circuit to make visualizations easier, and at the end, you will scale up the workflow to 50 qubits.

Requirements

Before starting this walkthrough, ensure that you have the following installed:

• Qiskit SDK v2.1 or later, with visualization support
• Qiskit Runtime v0.40 or later (pip install qiskit-ibm-runtime)

Setup

import matplotlib.pyplot as plt
import numpy as np

from qiskit.circuit.library import efficient_su2, unitary_overlap
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Batch, EstimatorV2 as Estimator

Small-scale simulator example

We will forgo this step since runtime error mitigation is not supported on simulators.

Hardware example

Step 1: Map classical inputs to a quantum problem

This walkthrough assumes that the classical problem has already been mapped to quantum. Begin by constructing a circuit and observable to measure. While the techniques used here apply to many different kinds of circuits, for simplicity this walkthrough uses the efficient_su2 circuit included in the Qiskit circuit library.

efficient_su2 is a parameterized quantum circuit designed to be efficiently executable on quantum hardware with limited qubit connectivity, while still being expressive enough to solve problems in application domains like optimization and chemistry. It is built by alternating layers of parameterized single-qubit gates with a layer containing a fixed pattern of two-qubit gates, for a chosen number of repetitions. The pattern of two-qubit gates can be specified by the user. Here you can use the built-in pairwise pattern because it minimizes the circuit depth by packing the two-qubit gates as densely as possible. This pattern can be executed using only linear qubit connectivity.

n_qubits = 10
reps = 1

circuit = efficient_su2(n_qubits, entanglement="pairwise", reps=reps)

circuit.decompose().draw("mpl", scale=0.7)

Output:

As our observable, let's take the Pauli $Z$ operator acting on the last qubit, $Z I \cdots I$. Note that the fact that the last qubit corresponds to the first element of this string is due to Qiskit's use of little-endian notation.

# Z on the last qubit (index -1) with coefficient 1.0
observable = SparsePauliOp.from_sparse_list(
    [("Z", [-1], 1.0)], num_qubits=n_qubits
)

At this point, you could proceed to run your circuit and measure the observable. However, you also want to compare the output of the quantum device with the correct answer — that is, the theoretical value of the observable, if the circuit had been executed without error. For small quantum circuits you can calculate this value by simulating the circuit on a classical computer, but this is not possible for larger, utility-scale circuits. You can work around this issue with the "mirror circuit" technique (also known as "compute-uncompute"), which is useful for benchmarking the performance of quantum devices.

Mirror circuit

In the mirror circuit technique, you concatenate the circuit with its inverse circuit, which is formed by inverting each gate of the circuit in reverse order. The resulting circuit implements the identity operator, which can trivially be simulated. Because the structure of the original circuit is preserved in the mirror circuit, executing the mirror circuit still gives an idea of how the quantum device would perform on the original circuit.

The following code cell assigns random parameters to your circuit, and then constructs the mirror circuit using the unitary_overlap class. Before mirroring the circuit, append a barrier instruction to it to prevent the transpiler from merging the two parts of the circuit on either side of the barrier and resulting in a transpiled circuit without any gates.

# Generate random parameters
rng = np.random.default_rng(1234)
params = rng.uniform(-np.pi, np.pi, size=circuit.num_parameters)

# Assign the parameters to the circuit
assigned_circuit = circuit.assign_parameters(params)

# Add a barrier to prevent circuit optimization of mirrored operators
assigned_circuit.barrier()

# Construct mirror circuit
mirror_circuit = unitary_overlap(assigned_circuit, assigned_circuit)

mirror_circuit.decompose().draw("mpl", scale=0.7)

Output:

Step 2: Optimize problem for quantum hardware execution

You must optimize your circuit before running it on hardware. This process involves a few steps:

• Pick a qubit layout that maps the virtual qubits of your circuit to physical qubits on the hardware.
• Insert swap gates as needed to route interactions between qubits that are not connected.
• Translate the gates in your circuit to Instruction Set Architecture (ISA) instructions that can directly be executed on the hardware.
• Perform circuit optimizations to minimize the circuit depth and gate count.

The transpiler built into Qiskit can perform all of these steps for you. Because this example uses a hardware-efficient circuit, the transpiler should be able to pick a qubit layout that does not require any swap gates to be inserted for routing interactions.

You need to choose the hardware device to use before you optimize your circuit. The following code cell requests the least-busy device with at least 127 qubits.

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=127
)

print(backend)

Output:

<IBMBackend('ibm_fez')>

You can transpile your circuit to your chosen backend by creating a pass manager and then running the pass manager on the circuit. An easy way to create a pass manager is to use the generate_preset_pass_manager function. See Transpile with pass managers for a more detailed explanation of transpiling with pass managers.

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend, seed_transpiler=1234
)
isa_circuit = pass_manager.run(mirror_circuit)

isa_circuit.draw("mpl", idle_wires=False, scale=0.7, fold=-1)

Output:

The transpiled circuit now contains only ISA instructions. All gates have been decomposed in terms of $\sqrt{X}$ gates and $R_z$ rotations, and CZ gates.

The transpilation process has mapped the virtual qubits of the circuit to physical qubits on the hardware. The information about the qubit layout is stored in the layout attribute of the transpiled circuit. The observable was also defined in terms of the virtual qubits, so you need to apply this layout to the observable, which you can do with the apply_layout method of SparsePauliOp.

isa_observable = observable.apply_layout(isa_circuit.layout)

print("Original observable:")
print(observable)
print()
print("Observable with layout applied:")
print(isa_observable)

Output:

Original observable:
SparsePauliOp(['ZIIIIIIIII'],
              coeffs=[1.+0.j])

Observable with layout applied:
SparsePauliOp(['IIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII'],
              coeffs=[1.+0.j])

Step 3: Execute using Qiskit primitives

You are now ready to run your circuit using the Estimator primitive.

Here you will submit five separate jobs, starting with no error suppression or mitigation, and successively enabling various error suppression and mitigation options available in Qiskit Runtime. For information about the options, refer to the following pages:

• Overview of all options
• Dynamical decoupling
• Resilience, including measurement error mitigation and zero-noise extrapolation (ZNE)
• Twirling

Because these jobs can run independently of each other, you can use batch mode to allow Qiskit Runtime to optimize the timing of their execution.

pub = (isa_circuit, isa_observable)

jobs = []

with Batch(backend=backend) as batch:
    estimator = Estimator(mode=batch)
    estimator.options.environment.job_tags = [
        "TUT_CEM_SS"
    ]  # add tag for this small scale job
    # Set number of shots
    estimator.options.default_shots = 100_000
    # Disable runtime compilation and error mitigation
    estimator.options.resilience_level = 0

    # Run job with no error mitigation
    job0 = estimator.run([pub])
    jobs.append(job0)

    # Add dynamical decoupling (DD)
    estimator.options.dynamical_decoupling.enable = True
    estimator.options.dynamical_decoupling.sequence_type = "XpXm"
    job1 = estimator.run([pub])
    jobs.append(job1)

    # Add readout error mitigation (DD + TREX)
    estimator.options.resilience.measure_mitigation = True
    job2 = estimator.run([pub])
    jobs.append(job2)

    # Add gate twirling (DD + TREX + Gate Twirling)
    estimator.options.twirling.enable_gates = True
    estimator.options.twirling.num_randomizations = "auto"
    job3 = estimator.run([pub])
    jobs.append(job3)

    # Add zero-noise extrapolation (DD + TREX + Gate Twirling + ZNE)
    estimator.options.resilience.zne_mitigation = True
    estimator.options.resilience.zne.noise_factors = (1, 3, 5)
    estimator.options.resilience.zne.extrapolator = ("exponential", "linear")
    job4 = estimator.run([pub])
    jobs.append(job4)

Step 4: Post-process and return result in desired classical format

Finally, you can analyze the data. Here you will retrieve the job results, extract the measured expectation values from them, and plot the values, including error bars of one standard deviation.

# Retrieve the job results
results = [job.result() for job in jobs]

# Unpack the PUB results (there's only one PUB result in each job result)
pub_results = [result[0] for result in results]

# Unpack the expectation values and standard errors
expectation_vals = np.array(
    [float(pub_result.data.evs) for pub_result in pub_results]
)
standard_errors = np.array(
    [float(pub_result.data.stds) for pub_result in pub_results]
)

# Plot the expectation values
fig, ax = plt.subplots()
labels = ["No mitigation", "+ DD", "+ TREX", "+ Twirling", "+ ZNE"]
ax.bar(
    range(len(labels)),
    expectation_vals,
    yerr=standard_errors,
    label="experiment",
)
ax.axhline(y=1.0, color="gray", linestyle="--", label="ideal")
ax.set_xticks(range(len(labels)))
ax.set_xticklabels(labels)
ax.set_ylabel("Expectation value")
ax.legend(loc="upper left")

plt.show()

Output:

At this small scale, it is difficult to see the effect of most of the error mitigation techniques, but zero-noise extrapolation does give a noticeable improvement. However, note that this improvement does not come for free, because the ZNE result also has a larger error bar.

Large-scale hardware example

When developing an experiment, it's useful to start with a small circuit to make visualizations and simulations easier. Now that you've developed and tested our workflow on a 10-qubit circuit, you can scale it up to 50 qubits. The following code cell repeats all of the steps in this walkthrough, but now applies them to a 50-qubit circuit.

n_qubits = 50
reps = 1

# Construct circuit and observable
circuit = efficient_su2(n_qubits, entanglement="pairwise", reps=reps)
observable = SparsePauliOp.from_sparse_list(
    [("Z", [-1], 1.0)], num_qubits=n_qubits
)

# Assign parameters to circuit
params = rng.uniform(-np.pi, np.pi, size=circuit.num_parameters)
assigned_circuit = circuit.assign_parameters(params)
assigned_circuit.barrier()

# Construct mirror circuit
mirror_circuit = unitary_overlap(assigned_circuit, assigned_circuit)

# Transpile circuit and observable
isa_circuit = pass_manager.run(mirror_circuit)
isa_observable = observable.apply_layout(isa_circuit.layout)

# Run jobs
pub = (isa_circuit, isa_observable)

jobs = []

with Batch(backend=backend) as batch:
    estimator = Estimator(mode=batch)
    estimator.options.environment.job_tags = [
        "TUT_CEM_LS"
    ]  # add tag for this large scale job
    # Set number of shots
    estimator.options.default_shots = 100_000
    # Disable runtime compilation and error mitigation
    estimator.options.resilience_level = 0

    # Run job with no error mitigation
    job0 = estimator.run([pub])
    jobs.append(job0)

    # Add dynamical decoupling (DD)
    estimator.options.dynamical_decoupling.enable = True
    estimator.options.dynamical_decoupling.sequence_type = "XpXm"
    job1 = estimator.run([pub])
    jobs.append(job1)

    # Add readout error mitigation (DD + TREX)
    estimator.options.resilience.measure_mitigation = True
    job2 = estimator.run([pub])
    jobs.append(job2)

    # Add gate twirling (DD + TREX + Gate Twirling)
    estimator.options.twirling.enable_gates = True
    estimator.options.twirling.num_randomizations = "auto"
    job3 = estimator.run([pub])
    jobs.append(job3)

    # Add zero-noise extrapolation (DD + TREX + Gate Twirling + ZNE)
    estimator.options.resilience.zne_mitigation = True
    estimator.options.resilience.zne.noise_factors = (1, 3, 5)
    estimator.options.resilience.zne.extrapolator = ("exponential", "linear")
    job4 = estimator.run([pub])
    jobs.append(job4)

# Retrieve the job results
results = [job.result() for job in jobs]

# Unpack the PUB results (there's only one PUB result in each job result)
pub_results = [result[0] for result in results]

# Unpack the expectation values and standard errors
expectation_vals = np.array(
    [float(pub_result.data.evs) for pub_result in pub_results]
)
standard_errors = np.array(
    [float(pub_result.data.stds) for pub_result in pub_results]
)

# Plot the expectation values
fig, ax = plt.subplots()
labels = ["No mitigation", "+ DD", "+ TREX", "+ Twirling", "+ ZNE"]
ax.bar(
    range(len(labels)),
    expectation_vals,
    yerr=standard_errors,
    label="experiment",
)
ax.axhline(y=1.0, color="gray", linestyle="--", label="ideal")
ax.set_xticks(range(len(labels)))
ax.set_xticklabels(labels)
ax.set_ylabel("Expectation value")
ax.legend(loc="upper left")

plt.show()

Output:

When you compare the 50-qubit results with the 10-qubit results from earlier, you might note the following (your results might differ across runs):

• All experiments give results closer to the ideal value and all the error bars are smaller.
• The addition of dynamical decoupling might have worsened performance compared to the no-mitigation case. This is not surprising, because the circuit is very dense. Dynamical decoupling is primarily useful when there are large gaps in the circuit, during which qubits sit idle without gates being applied to them. When these gaps are not present, dynamical decoupling is not effective, and can actually worsen performance due to errors in the dynamical decoupling pulses themselves. The 10-qubit circuit may have been too small for us to observe this effect.
• With zero-noise extrapolation, the result is very close to the ideal value. This demonstrates the power of ZNE.

Next steps