Zero-Noise-Extrapolation mit probabilistischer Fehlerverstärkung

Zero-Noise-Extrapolation (ZNE) ist eine Technik zur Minderung von Quantenrauschen bei der Berechnung von Erwartungswerten von Observablen. Sie funktioniert, indem ein Quantencircuit bei zunehmenden Rauschskalen ausgeführt wird, Erwartungswerte bezüglich dieser Rauschskalen bestimmt werden und daraus auf einen Erwartungswert ohne Rauschen extrapoliert wird. Sowohl die Rauschverstärkung als auch die Extrapolation lassen sich auf unterschiedliche Weise umsetzen. In dieser Anleitung verstärkst du das Rauschen mit probabilistischer Fehlerverstärkung (PEA), wie in "Evidence for the utility of quantum computing before fault tolerance" vorgeschlagen. Um mit PEA zu verstärken, musst du zuerst das Rauschen modellieren, das die verschränkenden Gates des Circuits beeinflusst. Du kannst die Verstärkung auch mit digitalem Gate Folding umsetzen, das kein Rauschmodell benötigt und in der Anleitung zu Gate Folding behandelt wird. Du lernst das Rauschen mit qiskit-noise-learning und extrapolierst mit der eingebauten exponentiellen Fit-Strategie.

Beachte: ZNE kann Quantenrauschen zwar bei geringeren Sampling-Kosten als andere Methoden wirksam mindern, bietet aber keine Schranken für den Bias.

Der Workflow sieht so aus:

1. Einen Circuit und eine Observable festlegen
2. Exakte Erwartungswerte mit dem Statevector-Simulator von Qiskit berechnen
3. Verrauschte, ungeminderte Erwartungswerte berechnen
4. Das Rauschen der verschränkenden Schichten mit qiskit-noise-learning lernen
5. Ein ZNE-Programm mit probabilistischer Fehlerverstärkung vorbereiten und Circuits mit unterschiedlichen Rauschskalen samplen
6. Auf den Zero-Noise-Grenzwert extrapolieren
7. Die geminderten Erwartungswerte vergleichen

Einen Circuit und eine Observable festlegen

In dieser Anleitung berechnest du die standortgemittelte Magnetisierung $\langle Z \rangle_{\text{avg}} = \frac{1}{N} \sum_{i=1}^{N} \langle Z_i \rangle$ eines trotterisierten Ising-Modells mit transversalem Feld (TFIM) auf acht Qubits in einer Ringtopologie. Du durchläufst mehrere Werte der Stärke des transversalen Felds und stellst die Magnetisierung als Funktion der Feldstärke dar. In einer Ringtopologie mit gerader Qubit-Anzahl lassen sich die Zwei-Qubit-Gate-Schichten in zwei eindeutige Schichten aufteilen, wie die roten und grünen Farben in der Kopplungskarte unten zeigen. Circuits so zu entwerfen, dass sie eine minimale Menge eindeutiger verschränkender Schichten haben, ist wichtig, wenn du ein Rauschmodell lernen möchtest, weil die eindeutigen Schichten unabhängig voneinander gelernt werden müssen. Wenn dein Circuit zu viele eindeutige Schichten hat, kann das Rauschlernen aufgrund der Sampling-Kosten für das Lernen vieler Schichten undurchführbar werden.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"IPython": "IPython", "matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "qiskit_mitigation": "qiskit-mitigation", "qiskit_noise_learning": "qiskit-noise-learning", "samplomatic": "samplomatic", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import warnings

import IPython.display
import numpy as np
from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.visualization import plot_coupling_map

warnings.filterwarnings("ignore")

def trotter_circuit(num_qubits, theta, layers, num_steps):
    """Trotterized transverse-field Ising evolution: RX on every site, then a CZ layer per coupling layer."""
    circuit = QuantumCircuit(num_qubits)
    for _ in range(num_steps):
        circuit.rx(theta, range(num_qubits))
        for layer in layers:
            circuit.barrier()
            for edge in layer:
                circuit.cz(*edge)
                circuit.sdg(edge)
            circuit.barrier()
    return circuit

num_qubits = 8
trotter_steps = 15
ring_coupling = [(i, (i + 1) % num_qubits) for i in range(num_qubits)]
layers = [
    ring_coupling[0::2],
    ring_coupling[1::2],
]  # two disjoint layers of CZ gates
theta = Parameter("theta")
circuit = trotter_circuit(num_qubits, theta, layers, trotter_steps)

# Parameter values: theta in [0, 0.25 pi] - values for the transverse-field sweep
data_points = 6
parameter_values = np.linspace(0, 0.25 * np.pi, data_points).reshape(
    (data_points, 1)
)

# Define observables - Z on each site
observables = [
    SparsePauliOp("I" * (num_qubits - q - 1) + "Z" + "I" * q)
    for q in range(num_qubits)
]

print("Edges of the same color may execute their gates simultaneously.")
IPython.display.display(
    plot_coupling_map(
        num_qubits=num_qubits,
        qubit_coordinates=None,
        coupling_map=ring_coupling,
        figsize=(4, 4),
        qubit_size=150,
        line_width=5,
        line_color=[
            "red" if edge in layers[0] else "green" for edge in ring_coupling
        ],
        font_size=90,
    )
)
print("Visualization of 1 Trotter step")
trotter_circuit(num_qubits, theta, layers, 1).draw("mpl", scale=0.7)

Edges of the same color may execute their gates simultaneously.

Visualization of 1 Trotter step

Exakte Erwartungswerte mit dem Statevector-Simulator von Qiskit berechnen

Verwende die Klasse Statevector von Qiskit, um den exakten Erwartungswert von $Z_i$ an jedem Gitterplatz für jeden Wert von $\theta$ zu berechnen.

from qiskit.quantum_info import Statevector
import matplotlib.pyplot as plt

exact_evs = np.array(
    [
        [
            Statevector(circuit.assign_parameters([value]))
            .expectation_value(obs)
            .real
            for value in parameter_values.ravel()
        ]
        for obs in observables
    ]
)  # shape (observables, parameter values)

x_vals = parameter_values.ravel() / np.pi
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, exact_evs.mean(axis=0), "ko-", label="Exact")
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel(r"$\langle Z \rangle_{\text{avg}}$")
ax.legend()
plt.tight_layout()
plt.show()

Verrauschte, ungeminderte Erwartungswerte berechnen

Ein QuantumProgram enthält einen Quantencircuit und ein Samplex in einem Format, das für das Executor-Primitive geeignet ist. MitigationTask.prepare erstellt aus dem Circuit, den Observablen und den Parameterwerten ein QuantumProgram und umschließt den Circuit mit den angegebenen Boxing-Optionen. Setze die Boxing-Option measure_annotations auf "change_basis", wodurch die Messungen in die benötigten Basen gedreht werden, und füge add_tags: "unique_box" hinzu, damit der Simulator ein Rauschmodell auf die Boxen des Circuits abbilden kann. Das Programm wird standardmäßig Pauli-getwirlt, was auf Hardware kohärentes Rauschen in das Pauli-Rauschen umwandelt, das die Mitigationsmethoden voraussetzen. num_randomizations legt fest, wie viele getwirlte Instanzen jedes Circuits gesampelt werden, und shots_per_randomization legt fest, wie viele Shots jede Instanz erhält, sodass ihr Produkt das Shot-Budget pro Circuit ist.

Ein Nutzer mit QPU-Zugang würde das Programm mit einem Executor ausführen, der aus einem Backend des QiskitRuntimeService erstellt wird. Hier übergibst du dem Executor einen AerSimulator, der das Programm lokal ausführt, und nutzt dessen Simulatoroption layer_noise_model, um auf jeder Schicht verschränkender Gates ein zufälliges Pauli-Lindblad-Rauschmodell einzuspeisen. Es wird mit random_pauli_lindblad_noise aus qiskit-mitigation erzeugt, mit einem Fehler pro Gate-Schicht (EPLG) von 0.003, typisch für aktuelle Heron-QPUs.

from qiskit.circuit import CircuitInstruction
from qiskit.quantum_info import PauliLindbladMap, pauli_basis
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import Executor
from qiskit_mitigation import MitigationTask, find_combined_unique_layers

def generate_random_pauli_lindblad_noise(
    layers,
    num_qubits: int,
    eplg: float = 3e-3,
    seed: int | None = None,
) -> list[tuple[CircuitInstruction, PauliLindbladMap]]:
    """Generate random Pauli-Lindblad noise over the full Pauli basis."""
    layer_noise_model = []
    for layer in layers:
        layer_instructions = layer.operation.params[0]
        contain_measurement = any(
            inst.name.startswith("measure") for inst in layer_instructions
        )
        if contain_measurement:  # do not add measurement noise
            layer_noise = PauliLindbladMap.identity(num_qubits)
        else:
            edges = [
                (inst.qubits[0]._index, inst.qubits[1]._index)
                for inst in layer_instructions
            ]  # extract the edges of the layer

            basis_paulis = [p for p in pauli_basis(2) if np.sum(p.x + p.z)]
            rng = np.random.default_rng(seed=seed)
            num_generators = len(basis_paulis)
            rates = rng.random(len(basis_paulis)) * eplg / num_generators

            layer_noise = PauliLindbladMap.from_sparse_list(
                [
                    (pauli.to_label(), edge, rate)
                    for pauli, rate in zip(basis_paulis, rates, strict=True)
                    for edge in edges
                ],
                num_qubits=num_qubits,
            )
        layer_noise_model.append((layer, layer_noise))
    return layer_noise_model

num_randomizations = 50
shots_per_randomization = 64

boxing_options = {
    "measure_annotations": "change_basis",
    "add_tags": "unique_box",
}

# find the unique layers to generate the layered noise for
box_layers = find_combined_unique_layers(
    [circuit], custom_boxing_options=boxing_options
)
layer_noise_model = generate_random_pauli_lindblad_noise(
    box_layers, num_qubits, seed=42
)

# a local noisy simulator in place of a QPU backend
noisy_executor = Executor(
    AerSimulator(),
    options={"simulator": {"layer_noise_model": layer_noise_model}},
)

task = MitigationTask()
program = task.prepare(
    circuit=circuit,
    observables=observables,
    parameters=parameter_values,
    num_randomizations=num_randomizations,
    shots_per_randomization=shots_per_randomization,
    custom_boxing_options=boxing_options,
)
program._semantic_role = "simulation"  # flag the program for the simulator
unmitigated_res = task.postprocess(noisy_executor.run(program).result())

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, exact_evs.mean(axis=0), "k-", label="Exact")
ax.errorbar(
    x_vals,
    np.mean(unmitigated_res.data.evs, axis=0),
    np.mean(unmitigated_res.data.stds, axis=0),
    fmt="o-",
    label="Unmitigated",
)
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel(r"$\langle Z \rangle_{\text{avg}}$")
ax.legend()
plt.tight_layout()
plt.show()

Lerne das Rauschen der verschränkenden Schichten mit qiskit-noise-learning

PEA beginnt damit, das getwirlte Rauschmodell jeder Schicht verschränkender Gates im Circuit zu lernen. Wir nehmen ein Markov-Rauschmodell an, daher muss nur das Rauschmodell jeder einzigartigen Schicht verschränkender Gates gelernt werden, da angenommen wird, dass die Tiefe des Circuits den Rauschkanal einer gegebenen verschränkenden Schicht nicht beeinflusst. Je genauer das gelernte Rauschmodell ist, desto genauer wird die Verstärkung sein.

Wir verwenden qiskit-noise-learning, um das Rauschen zu lernen. Der Lerner wiederholt jede Schicht fragment_depths-mal in seinen Lern-Circuits, twirlt jeden Circuit num_randomizations-mal und passt den Abfall der gemessenen Fidelitäten an ein Pauli-Lindblad-Modell an. Der folgende Code zeichnet das gelernte Rauschmodell, wobei die beiden Schichten in einer Karte zusammengeführt sind. Jedes Qubit ist ein Kreisdiagramm aus seinen über die beiden Schichten summierten Gewicht-1-Raten für X, Y und Z, und jede Kante zeigt die gelernten Gewicht-2-Raten; eine Rate von 3e-4 bedeutet, dass dieser Fehler bei etwa drei von 10.000 Anwendungen der Schicht auftritt.

from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.quantum_info import PauliList
from qiskit_ibm_runtime.results import LayerError, PauliLindbladError
from qiskit_ibm_runtime.visualization.embeddings import Embedding
from qiskit_mitigation import PEA
from qiskit_noise_learning.aer_executor import AerExecutor
from qiskit_noise_learning.noise_learner import LearningOptions, NoiseLearner
from samplomatic import InjectNoise, Tag
from samplomatic.utils import get_annotation

backend = GenericBackendV2(
    num_qubits=num_qubits,
    basis_gates=["id", "rz", "sx", "x", "cz"],
    coupling_map=ring_coupling,
)

task_pea = PEA()
# learn only the layers PEA will inject noise into, i.e. not the measurement layer
pea_layers = [
    layer
    for layer in task_pea.find_unique_layers(circuit, boxing_options)
    if get_annotation(layer.operation, InjectNoise)
]

learning_options = LearningOptions(
    fragment_depths=[2, 4, 8, 16],
    num_randomizations=50,
    shots_per_randomizations=64,
)
# qiskit-noise-learning is not compatible with the Executor in local mode,
# and requires the simulator present in the package for running simulation
noise_dict = {
    get_annotation(layer.operation, Tag).ref: noise_layer
    for layer, noise_layer in layer_noise_model
}
noise_learner_simulator = AerExecutor(
    qasm_simulator=AerSimulator(), noise_dict=noise_dict
)
learner = NoiseLearner(
    backend, options=learning_options, executor=noise_learner_simulator
)
learned_noise = (
    learner.run(pea_layers).result().to_dict()
)  # keyed by InjectNoise reference

# Draw the learned model on the ring, laid out as a square with the qubits in ring order.
# Composing the two layers' maps sums their rates, giving one map over all the ring's edges.
embedding = Embedding(
    [(0, 0), (0, 2), (0, 4), (2, 4), (4, 4), (4, 2), (4, 0), (2, 0)],
    ring_coupling,
)
first, second = learned_noise.values()
learned_model = first.compose(second).simplify()
error = PauliLindbladError(
    PauliList(
        learned_model.get_qubit_sparse_pauli_list_copy().to_pauli_list()
    ),
    learned_model.rates,
)
fig = LayerError(circuit, list(range(num_qubits)), error).draw_map(
    embedding, radius=0.4, height=460, width=620
)
fig.update_traces(
    marker={
        "cmin": 0,
        "cmax": max(learned_model.rates),
        "colorbar": {"tickformat": ".1e", "title": "rate"},
    },
    selector={"mode": "markers"},
)
fig.show()

Bereite ein ZNE-Programm mit probabilistischer Fehlerverstärkung vor und sample Circuits mit unterschiedlichen Rauschskalen

Hier bereitest du ein ZNE-Programm vor, das Circuits mit zunehmenden Rauschskalen enthält. Bei PEA steht jeder Rauschfaktor für das Vielfache des gelernten Rauschens, das der Circuit erfährt: Ein Rauschfaktor von 2.0 speist nach jeder Schicht eine zusätzliche Kopie des gelernten Rauschkanals dieser Schicht ein, zusätzlich zum Rauschen, das die Schicht bereits hat. Da PEA die gelernten Rauschraten direkt skaliert, statt Gates zu wiederholen, ist jeder Rauschfaktor möglich, auch nicht ganzzahlige, und der Circuit selbst wird nicht verändert: Die Programmelemente teilen sich einen Vorlagen-Circuit und unterscheiden sich nur im Rauschen, das vom Samplex eingespeist wird. Das eingespeiste Rauschen wird einmal pro Randomisierung gesampelt, daher wird der statistische Fehler von PEA durch die Anzahl der Randomisierungen statt durch die Shots pro Randomisierung bestimmt; wir verwenden deshalb viele Randomisierungen mit wenigen Shots und halten die Gesamtzahl der Shots gleich wie beim nicht mitigierten Programm. Es ist wichtig, nicht zu wenige Randomisierungen zu verwenden, da sonst ein Sampling-Bias bei der Rauscheinspeisung auftreten kann.

Der gesamte Sampling-Overhead skaliert linear mit der Anzahl der Rauschfaktoren.

noise_factors = [1, 2, 3]
program_pea = task_pea.prepare(
    circuit=circuit,
    observables=observables,
    parameters=parameter_values,
    num_randomizations=200,
    shots_per_randomization=16,
    custom_boxing_options=boxing_options,
    noise_maps=learned_noise,
    noise_factors=noise_factors,
)

Führe das Programm abschließend auf dem verrauschten Simulator aus. Seine Ergebnisse speisen die Extrapolation im nächsten Schritt.

program_pea._semantic_role = (
    "simulation"  # flag the program for the simulator
)
job_pea_results = noisy_executor.run(program_pea).result()

Extrapoliere auf den Grenzwert ohne Rauschen

Nachdem du nun Daten für jede Rauschskala hast, kannst du eine Funktion an die Daten anpassen und auf den Wert ohne Rauschen zurückextrapolieren, um einen Erwartungswert der rauschfreien Observablen abzuschätzen. Unter Pauli-Rauschen fällt eine Observable exponentiell mit dem Rauschfaktor ab, daher wählst du als erste Wahl den exponentiellen Extrapolator und weichst auf den linearen aus, falls die exponentielle Anpassung fehlschlägt. Um die Extrapolation zu visualisieren, zeichne die über die Plätze gemittelte Magnetisierung bei jedem Rauschfaktor, die exponentielle Anpassung und den Mittelwert über die Plätze der extrapolierten Erwartungswerte am Punkt ohne Rauschen. Wir sehen, dass die Anpassung bei dieser speziellen Observablen leicht zu niedrig liegt.

from scipy.optimize import curve_fit

pea_mitigated_res = task_pea.postprocess(
    job_pea_results, extrapolator=["exponential", "linear"]
)

# The result holds the expectation value at each noise factor and the extrapolated value for
# every observable; we average them over the sites and fit the same model to draw the curve.
angle = 3  # theta = 0.15 pi, the angle with the largest unmitigated error
amplifications = pea_mitigated_res.data.evs_noise_factors[:, angle, :].mean(
    axis=0
)
x_fit = np.linspace(0, max(noise_factors), 51)
params, _ = curve_fit(
    lambda x, a, b: a * np.exp(b * x),
    noise_factors,
    amplifications,
    p0=[1.0, -0.1],
)

plt.plot(
    noise_factors, amplifications, "o", color="red", label="amplifications"
)
plt.plot(
    x_fit,
    params[0] * np.exp(params[1] * x_fit),
    "--",
    color="black",
    label="exponential fit",
)
plt.plot(
    [0],
    pea_mitigated_res.data.evs[:, angle].mean(),
    "*",
    color="black",
    markersize=15,
    label="extrapolation",
)
plt.plot(
    [0],
    exact_evs[:, angle].mean(),
    "_",
    color="red",
    markersize=15,
    markeredgewidth=3,
    label="exact",
)
plt.xlabel("noise factor")
plt.ylabel(r"$\langle Z \rangle_{\text{avg}}$")
plt.legend()
plt.show()

Vergleiche die mitigierten Erwartungswerte

Zeichne die über die Plätze gemittelte Magnetisierung jeder Methode über den Bereich der Ising-Rotationswinkel. PEA stellt die exakte Magnetisierung innerhalb seiner Fehlerbalken wieder her, während das nicht mitigierte Ergebnis deutlich darunter verzerrt ist.

results = {
    "Unmitigated": unmitigated_res,
    "PEA": pea_mitigated_res,
}

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, exact_evs.mean(axis=0), "k-", label="Exact")
for label, res in results.items():
    ax.errorbar(
        x_vals,
        np.mean(res.data.evs, axis=0),
        np.mean(res.data.stds, axis=0),
        fmt="o-",
        label=label,
    )
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel(r"$\langle Z \rangle_{\text{avg}}$")
ax.legend()
ax.set_title(
    f"Error-Mitigated Magnetization [{num_qubits}-qubit, {trotter_steps} Trotter steps]"
)
plt.tight_layout()
plt.show()