Extrapolation à bruit nul avec amplification probabiliste d'erreur

L'extrapolation à bruit nul (ZNE) est une technique d'atténuation du bruit quantique dans les calculs de valeurs d'espérance d'observables. Elle consiste à exécuter un circuit quantique à des échelles de bruit croissantes, à calculer les valeurs d'espérance pour ces échelles de bruit, puis à extrapoler pour estimer une valeur d'espérance à bruit nul. Les étapes d'amplification du bruit et d'extrapolation peuvent être implémentées de différentes manières. Dans ce guide, tu vas amplifier le bruit par amplification probabiliste d'erreur (PEA), comme suggéré dans "Evidence for the utility of quantum computing before fault tolerance". Pour amplifier avec la PEA, tu dois d'abord modéliser le bruit affectant les portes d'intrication du circuit. Tu peux aussi implémenter l'amplification par repliement numérique de portes, qui ne nécessite pas de modèle de bruit et est traité dans le guide sur le repliement de portes. Tu vas apprendre le bruit avec qiskit-noise-learning et extrapoler avec la stratégie d'ajustement exponentielle intégrée.

Garde à l'esprit que, si la ZNE peut être efficace pour atténuer le bruit quantique avec un coût d'échantillonnage inférieur à celui d'autres méthodes, elle ne donne aucune borne sur le biais.

Le flux de travail est le suivant :

1. Spécifier un circuit et une observable
2. Calculer les valeurs d'espérance exactes avec le simulateur de vecteur d'état de Qiskit
3. Calculer les valeurs d'espérance bruitées, non atténuées
4. Apprendre le bruit des couches d'intrication avec qiskit-noise-learning
5. Préparer un programme ZNE avec amplification probabiliste d'erreur et échantillonner des circuits d'échelles de bruit variées
6. Extrapoler vers la limite à bruit nul
7. Comparer les valeurs d'espérance atténuées

Spécifier un circuit et une observable

Dans ce guide, tu vas calculer l'aimantation moyenne par site $\langle Z \rangle_{\text{avg}} = \frac{1}{N} \sum_{i=1}^{N} \langle Z_i \rangle$ d'un modèle d'Ising à champ transverse (TFIM) trotterisé sur huit qubits en topologie en anneau. Tu vas balayer plusieurs valeurs de l'intensité du champ transverse et tracer l'aimantation en fonction de l'intensité du champ. Dans une topologie en anneau sur un nombre pair de qubits, les couches de portes à deux qubits peuvent être scindées en deux couches distinctes, comme le montrent les couleurs rouge et verte dans la carte de couplage ci-dessous. Concevoir des circuits avec un ensemble minimal de couches d'intrication distinctes est important si tu veux apprendre un modèle de bruit, car les couches distinctes doivent être apprises indépendamment. Si ton circuit contient trop de couches distinctes, l'apprentissage du bruit peut devenir infaisable en raison du coût d'échantillonnage nécessaire pour apprendre de nombreuses couches.

# 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

Calculer les valeurs d'espérance exactes avec le simulateur de vecteur d'état de Qiskit

Utilise la classe Statevector de Qiskit pour calculer la valeur d'espérance exacte de $Z_i$ sur chaque site pour chaque valeur de $\theta$.

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

Calculer les valeurs d'espérance bruitées, non atténuées

Un QuantumProgram contient un circuit quantique et un samplex dans un format adapté à la primitive Executor. MitigationTask.prepare construit un QuantumProgram à partir du circuit, des observables et des valeurs de paramètres, et encadre le circuit dans des boîtes avec les options de boxing données. Définis l'option de boxing measure_annotations sur "change_basis", qui fait pivoter les mesures vers les bases requises, et ajoute add_tags: "unique_box" pour que le simulateur puisse appliquer un modèle de bruit aux boîtes du circuit. Le programme est par défaut soumis à un twirling de Pauli, qui sur du matériel transforme le bruit cohérent en bruit de Pauli, celui que supposent les méthodes d'atténuation. num_randomizations définit combien d'instances twirlées de chaque circuit sont échantillonnées, et shots_per_randomization définit combien de shots reçoit chaque instance ; leur produit est donc le budget de shots par circuit.

Un utilisateur ayant accès à un QPU exécuterait le programme avec un Executor construit à partir d'un Backend du QiskitRuntimeService. Ici, tu vas passer un AerSimulator à l'Executor, qui exécute le programme localement, et utiliser son option de simulateur layer_noise_model pour injecter un modèle de bruit de Pauli-Lindblad aléatoire sur chaque couche de Gates d'intrication, généré avec random_pauli_lindblad_noise de qiskit-mitigation avec une erreur par Gate en couche (EPLG) de 0.003, typique des QPU Heron actuels.

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

Apprendre le bruit des couches d'intrication avec qiskit-noise-learning

La PEA commence par apprendre le modèle de bruit twirlé de chaque couche de Gates d'intrication du circuit. Nous supposons un modèle de bruit markovien ; seul le modèle de bruit de chaque couche unique de Gates d'intrication doit donc être appris, car la profondeur du circuit est supposée ne pas affecter le canal de bruit d'une couche d'intrication donnée. Plus le modèle de bruit appris est précis, plus l'amplification sera précise.

Nous utilisons qiskit-noise-learning pour apprendre le bruit. Le module d'apprentissage répète chaque couche fragment_depths fois dans ses circuits d'apprentissage, applique le twirling à chaque circuit num_randomizations fois et ajuste la décroissance des fidélités mesurées à un modèle de Pauli-Lindblad. Le code suivant dessine le modèle de bruit appris, avec les deux couches fusionnées en une seule carte. Chaque qubit est un diagramme circulaire de ses taux X, Y et Z de poids 1 additionnés sur les deux couches, et chaque arête indique les taux de poids 2 appris ; un taux de 3e-4 signifie que l'erreur se produit environ trois fois sur 10 000 utilisations de la couche.

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

Préparer un programme ZNE avec amplification probabiliste des erreurs et échantillonner des circuits à différentes échelles de bruit

Ici, tu vas préparer un programme ZNE contenant des circuits d'échelles de bruit croissantes. Avec la PEA, chaque facteur de bruit représente le multiple du bruit appris que subit le circuit : un facteur de bruit de 2.0 injecte une copie supplémentaire du canal de bruit appris de chaque couche après cette couche, en plus du bruit déjà présent dans la couche. Comme la PEA met à l'échelle directement les taux de bruit appris au lieu de répéter des Gates, tout facteur de bruit est possible, y compris non entier, et le circuit lui-même n'est pas modifié : les éléments du programme partagent un même circuit modèle et ne diffèrent que par le bruit injecté par le samplex. Le bruit injecté est échantillonné une fois par randomisation ; l'erreur statistique de la PEA est donc contrôlée par le nombre de randomisations plutôt que par les shots par randomisation ; nous utilisons donc beaucoup de randomisations avec quelques shots chacune, en gardant le nombre total de shots identique à celui du programme non atténué. Il est important de ne pas utiliser trop peu de randomisations, sinon tu pourrais rencontrer un biais d'échantillonnage dans l'injection de bruit.

Le surcoût total d'échantillonnage augmente linéairement avec le nombre de facteurs de bruit.

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

Enfin, exécute le programme sur le simulateur bruité. Ses résultats alimentent l'extrapolation de l'étape suivante.

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

Extrapoler vers la limite de bruit nul

Maintenant que tu disposes de données pour chaque échelle de bruit, tu peux ajuster une fonction aux données et extrapoler jusqu'à la valeur à bruit nul pour estimer une valeur moyenne de l'observable sans bruit. Sous bruit de Pauli, une observable décroît exponentiellement avec le facteur de bruit ; tu choisiras donc l'extrapolateur exponentiel en premier choix, et tu te rabattras sur l'extrapolateur linéaire si l'ajustement exponentiel échoue. Pour visualiser l'extrapolation, trace l'aimantation moyennée sur les sites calculée à chaque facteur de bruit, l'ajustement exponentiel et la moyenne sur les sites des valeurs moyennes extrapolées au point de bruit nul. Nous voyons que l'ajustement sous-estime légèrement pour cette observable particulière.

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

Comparer les valeurs moyennes atténuées

Trace l'aimantation moyennée sur les sites de chaque méthode sur l'ensemble des angles de rotation d'Ising. La PEA retrouve l'aimantation exacte dans ses barres d'erreur, tandis que le résultat non atténué est biaisé bien en dessous.

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