Zero-Noise-Extrapolation mit Gate Folding und benutzerdefinierter Extrapolation

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. Diese Anleitung verstärkt das Rauschen mit "digitalem Gate Folding", bei dem Zwei-Qubit-Gates durch äquivalente Folgen aus dem Gate und seinem Inversen ersetzt werden. Ersetzt man zum Beispiel ein unitäres $U$ durch $U U^† U$, ergibt sich ein Rauschverstärkungsfaktor von 3. Du kannst die Verstärkung auch mit probabilistischer Fehlerverstärkung (PEA) umsetzen, die einen zusätzlichen Schritt zum Lernen eines Rauschmodells erfordert und in der Anleitung zur probabilistischen Fehlerverstärkung behandelt wird. Diese Anleitung extrapoliert mit einigen eingebauten Fit-Strategien (linear und exponentiell) und zeigt außerdem, wie du deine eigene benutzerdefinierte Fit-Funktion implementierst.

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. Ein ZNE-Programm mit digitalem Gate Folding vorbereiten und Circuits mit unterschiedlichen Rauschskalen samplen
5. Mit eingebauten und benutzerdefinierten Fit-Funktionen auf den Zero-Noise-Grenzwert extrapolieren
6. Die geminderten Erwartungswerte vergleichen

Einen Circuit und eine Observable festlegen

Berechne 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.

# 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", "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 numpy as np
from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

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("Visualization of 1 Trotter step")
trotter_circuit(num_qubits, theta, layers, 1).draw("mpl", scale=0.7)

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.

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

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. Du setzt die Boxing-Option measure_annotations auf "change_basis", wodurch die Messungen in die benötigten Basen gedreht werden, und fügst 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-Zugriff 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, das mit random_pauli_lindblad_noise aus qiskit-mitigation mit einem Fehler pro geschichtetem Gate (EPLG) von 0.003 erzeugt wird, 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()

Ein ZNE-Programm mit digitalem Gate Folding vorbereiten und Circuits mit unterschiedlichen Rauschskalen samplen

Bereite ein ZNE-Programm vor, das Circuits mit zunehmenden Rauschskalen enthält. Jeder Rauschfaktor steht für das Vielfache der verschränkenden Gates, die der Circuit enthält (ein Rauschfaktor von 3.0 bedeutet dreimal so viele verschränkende Gates wie im ursprünglichen Circuit). Der Circuit jedes Rauschfaktors wird mit derselben Anzahl an Randomisierungen und Shots getwirlt wie das ungeminderte Programm.

from qiskit_mitigation import GateFolding

task_gf = GateFolding()
program_zne = task_gf.prepare(
    circuit=circuit,
    observables=observables,
    parameters=parameter_values,
    num_randomizations=num_randomizations,
    shots_per_randomization=shots_per_randomization,
    custom_boxing_options=boxing_options,
    noise_factors=[1, 3, 5],  # ZNE option
)

Zeichne den ersten Trotter-Schritt der getwirlten Vorlage mit Rauschfaktor 1 und mit Rauschfaktor 3, um die Verstärkung zu veranschaulichen. Bei Rauschfaktor 3.0 wird jede Schicht des Circuits durch die Schicht, ihr Inverses und erneut die Schicht ersetzt.

import IPython.display

def first_trotter_step(template, noise_factor):
    """The template up to the closing barrier of its first Trotter step (two layers, each folded ``noise_factor`` times)."""
    toy = template.copy_empty_like()
    for instruction in template.data:
        toy.append(instruction)
        if (
            instruction.operation.name == "barrier"
            and instruction.operation.label.startswith(
                f"R{2 * noise_factor - 1}@"
            )
        ):
            break
    return toy

IPython.display.display(
    first_trotter_step(program_zne.items[0].circuit, 1).draw(
        "mpl", scale=0.7, plot_barriers=False
    )
)
first_trotter_step(program_zne.items[1].circuit, 3).draw(
    "mpl", fold=-1, plot_barriers=False
)

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

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

Mit eingebauten und benutzerdefinierten Fit-Funktionen auf den Zero-Noise-Grenzwert extrapolieren

Nachdem du nun Daten für jede Rauschskala hast, kannst du eine Funktion an die Daten fitten und zurück auf den Zero-Noise-Wert extrapolieren, um einen Erwartungswert der rauschfreien Observable abzuschätzen. Es gibt viele Möglichkeiten, Daten zu fitten, und die richtige Wahl der Funktion hängt von der Art der gesampelten Daten ab. GateFolding.postprocess bietet mehrere eingebaute Extrapolationsfunktionen, darunter linear, exponential und double_exponential. Die Methode postprocess nimmt eine Liste von Extrapolationsmethoden entgegen, wertet sie der Reihe nach aus und gibt den ersten Wert mit physikalisch gültigen Werten und endlicher Varianz zurück. Das bedeutet, dass die Reihenfolge, in der du die Extrapolatoren angibst, großen Einfluss darauf hat, welcher Extrapolator gewählt wird. Wenn du fallback als Extrapolator angibst, wird der Erwartungswert des kleinsten Rauschfaktors gewählt und zurückgegeben, falls keine gültigen Antworten gefunden wurden. Wird fallback nicht angegeben und liefert kein Extrapolator einen physikalisch gültigen Wert, wird das Ergebnis mit dem kleinsten std zurückgegeben. Unter Pauli-Rauschen zerfällt eine Observable exponentiell mit dem Rauschfaktor, daher ist bei 30 Schichten (15 Trotter-Schritte, 2 Schichten pro Schritt) verschränkender Gates der exponentielle Extrapolator die naheliegende Wahl; du kannst auch linear extrapolieren, was eine gängige Wahl ist, wenn der Zerfall gering ist.

Um eine benutzerdefinierte Fit-Funktion zu verwenden, übergib den Extrapolator "custom" zusammen mit einem custom_fit-Tupel, das die Fit-Funktion, eine Liste mit einem Startwert für jeden ihrer Parameter und Schranken enthält. Diese werden an scipy.optimize.curve_fit weitergereicht. Die Nachbearbeitung kann mit anderen Extrapolatoren wiederholt werden, ohne das Programm erneut auszuführen. Hier verwendest du eine benutzerdefinierte Quadratwurzel-Fit-Funktion. Um die Extrapolation zu visualisieren, stelle die standortgemittelte Magnetisierung dar, die bei jedem Rauschfaktor für $\theta = 0.15\pi$ berechnet wurde, dazu die drei Fits und den Standortdurchschnitt der extrapolierten Erwartungswerte am Zero-Noise-Punkt. Die drei Datenpunkte wirken nahezu linear, doch die lineare Extrapolation landet sichtbar unter dem exakten Wert. Wie vorhergesagt schätzt der exponentielle Fit den rauschfreien Erwartungswert am besten.

from scipy.optimize import curve_fit

zne_mitigated_exp = task_gf.postprocess(
    job_zne_results, extrapolator=["exponential", "linear"]
)
zne_mitigated_lin = task_gf.postprocess(
    job_zne_results, extrapolator=["linear"]
)

def sqrt_fit(x, a, b):
    return a * np.sqrt(x) + b

custom_fit = (sqrt_fit, [1.0, 1.0], (-np.inf, np.inf))
zne_mitigated_custom = task_gf.postprocess(
    job_zne_results,
    extrapolator=["custom", "exponential", "linear"],
    custom_fit=custom_fit,
)

# The result holds the expectation value at each noise factor and the extrapolated value, but
# not the fit parameters, so we refit the same models here to draw the curves.
noise_factors = np.array([1, 3, 5])
angle = 3  # theta = 0.15 pi, the angle with the largest unmitigated error
amplifications = zne_mitigated_exp.data.evs_noise_factors[:, angle, :].mean(
    axis=0
)
x_fit = np.linspace(0, 5, 51)
fits = {
    "exponential": (
        lambda x, a, b: a * np.exp(b * x),
        zne_mitigated_exp,
        "black",
    ),
    "linear": (lambda x, a, b: a * x + b, zne_mitigated_lin, "blue"),
    "square-root": (sqrt_fit, zne_mitigated_custom, "green"),
}

plt.plot(
    noise_factors, amplifications, "o", color="red", label="amplifications"
)
for name, (fit, res, color) in fits.items():
    params, _ = curve_fit(fit, noise_factors, amplifications, p0=[1.0, -0.1])
    plt.plot(
        x_fit, fit(x_fit, *params), "--", color=color, label=f"{name} fit"
    )
    plt.plot(
        [0],
        res.data.evs[:, angle].mean(),
        "*",
        color=color,
        markersize=15,
        label=f"{name} 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()

Die geminderten Erwartungswerte vergleichen

Stelle die standortgemittelte Magnetisierung jeder Methode über den gesamten Durchlauf der Ising-Drehwinkel dar. Die exponentielle Extrapolation stellt die exakte Magnetisierung innerhalb ihrer Fehlerbalken wieder her, während die lineare Extrapolation einen Großteil des Bias beseitigt, aber durchgängig zu niedrig liegt, weil der wahre Zerfall exponentiell ist. Der Quadratwurzel-Fit überschießt deutlich.

results = {
    "Unmitigated": unmitigated_res,
    "ZNE exponential": zne_mitigated_exp,
    "ZNE linear": zne_mitigated_lin,
    "ZNE square-root": zne_mitigated_custom,
}

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()