Ekstrapolacja do zerowego szumu z probabilistycznym wzmacnianiem błędów

Ekstrapolacja do zerowego szumu (ZNE) to technika łagodzenia szumu kwantowego w obliczeniach wartości oczekiwanych obserwabli. Polega na wykonaniu obwodu kwantowego przy rosnących skalach szumu, obliczeniu wartości oczekiwanych względem tych skal szumu i ekstrapolacji w celu oszacowania wartości oczekiwanej przy zerowym szumie. Zarówno etap wzmacniania szumu, jak i etap ekstrapolacji można zaimplementować na różne sposoby. W tym przewodniku wzmocnisz szum za pomocą probabilistycznego wzmacniania błędów (PEA), zgodnie z propozycją w "Evidence for the utility of quantum computing before fault tolerance". Aby wzmacniać za pomocą PEA, musisz najpierw zamodelować szum wpływający na bramki splątujące obwodu. Wzmacnianie można też zaimplementować za pomocą cyfrowego zagęszczania bramek, które nie wymaga modelu szumu i jest omówione w przewodniku o zagęszczaniu bramek. Nauczysz się szumu za pomocą qiskit-noise-learning i wyekstrapolujesz za pomocą wbudowanej wykładniczej strategii dopasowania.

Pamiętaj, że choć ZNE może skutecznie łagodzić szum kwantowy przy niższym koszcie próbkowania niż inne metody, nie daje żadnych ograniczeń na obciążenie (bias).

Przepływ pracy wygląda następująco:

1. Określ obwód i obserwablę
2. Oblicz dokładne wartości oczekiwane za pomocą symulatora wektora stanu Qiskit
3. Oblicz zaszumione, niełagodzone wartości oczekiwane
4. Naucz się szumu warstw splątujących za pomocą qiskit-noise-learning
5. Przygotuj program ZNE z probabilistycznym wzmacnianiem błędów i próbkuj obwody o różnych skalach szumu
6. Wyekstrapoluj do granicy zerowego szumu
7. Porównaj złagodzone wartości oczekiwane

Określ obwód i obserwablę

W tym przewodniku obliczysz średnią po węzłach magnetyzację $\langle Z \rangle_{\text{avg}} = \frac{1}{N} \sum_{i=1}^{N} \langle Z_i \rangle$ trotteryzowanego modelu Isinga z polem poprzecznym (TFIM) na ośmiu kubitach w topologii pierścienia. Przeprowadzisz przemiatanie po kilku wartościach natężenia pola poprzecznego i narysujesz magnetyzację jako funkcję natężenia pola. W topologii pierścienia na parzystej liczbie kubitów warstwy bramek dwukubitowych można podzielić na dwie unikalne warstwy, co pokazują czerwony i zielony kolor na poniższej mapie połączeń. Projektowanie obwodów tak, aby miały minimalny zestaw unikalnych warstw splątujących, jest ważne, jeśli chcesz uczyć się modelu szumu, ponieważ unikalne warstwy muszą być uczone niezależnie. Jeśli w obwodzie jest zbyt wiele unikalnych warstw, uczenie szumu może stać się niewykonalne z powodu kosztu próbkowania przy uczeniu wielu warstw.

# 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

Oblicz dokładne wartości oczekiwane za pomocą symulatora wektora stanu Qiskit

Użyj klasy Statevector z Qiskit, aby obliczyć dokładną wartość oczekiwaną $Z_i$ w każdym węźle dla każdej wartości $\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()

Oblicz zaszumione, niełagodzone wartości oczekiwane

QuantumProgram zawiera obwód kwantowy i samplex w formacie odpowiednim dla prymitywu Executor. MitigationTask.prepare buduje QuantumProgram z obwodu, obserwabli i wartości parametrów oraz otacza obwód pudełkami zgodnie z podanymi opcjami boxingu. Ustaw opcję boxingu measure_annotations na "change_basis", która obraca pomiary do wymaganych baz, i dodaj add_tags: "unique_box", aby symulator mógł odwzorować model szumu na pudełka obwodu. Program jest domyślnie poddawany Pauli twirlingowi, co na sprzęcie zamienia szum koherentny w szum Pauliego, którego zakładają metody łagodzenia. num_randomizations określa, ile zatwirlowanych instancji każdego obwodu jest próbkowanych, a shots_per_randomization określa, ile strzałów otrzymuje każda instancja, więc ich iloczyn to budżet strzałów na obwód.

Użytkownik z dostępem do QPU wykonałby program za pomocą Executor zbudowanego z backendu QiskitRuntimeService. Tutaj przekażesz do Executor symulator AerSimulator, który uruchamia program lokalnie, i użyjesz jego opcji symulatora layer_noise_model, aby wstrzyknąć losowy model szumu Pauliego-Lindblada do każdej warstwy bramek splątujących, wygenerowany za pomocą random_pauli_lindblad_noise z qiskit-mitigation z błędem na warstwę bramek (EPLG) równym 0.003, typowym dla obecnych QPU Heron.

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

Poznaj szum warstw splątujących za pomocą qiskit-noise-learning

PEA zaczyna od nauczenia się modelu szumu po twirlingu dla każdej warstwy bramek splątujących w obwodzie. Zakładamy markowowski model szumu, dlatego należy nauczyć się tylko modelu szumu każdej unikalnej warstwy bramek splątujących, ponieważ zakłada się, że głębokość obwodu nie wpływa na kanał szumu danej warstwy splątującej. Im dokładniejszy nauczony model szumu, tym dokładniejsze będzie wzmocnienie.

Do uczenia się szumu używamy qiskit-noise-learning. Learner powtarza każdą warstwę fragment_depths razy w swoich obwodach uczących, wykonuje twirling każdego obwodu num_randomizations razy i dopasowuje zanik zmierzonych wierności do modelu Pauliego-Lindblada. Poniższy kod rysuje nauczony model szumu, z dwiema warstwami połączonymi w jedną mapę. Każdy kubit jest wykresem kołowym swoich współczynników X, Y i Z o wadze 1, zsumowanych po obu warstwach, a każda krawędź pokazuje nauczone współczynniki o wadze 2; współczynnik 3e-4 oznacza, że błąd ten występuje średnio przy około trzech na każde 10 000 użyć warstwy.

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

Przygotuj program ZNE z probabilistycznym wzmacnianiem błędów i próbkuj obwody o różnych skalach szumu

Tutaj przygotujesz program ZNE zawierający obwody o rosnących skalach szumu. W PEA każdy współczynnik szumu oznacza krotność nauczonego szumu, której doświadcza obwód: współczynnik szumu 2.0 wstrzykuje jedną dodatkową kopię nauczonego kanału szumu każdej warstwy po tej warstwie, oprócz szumu, który warstwa już ma. Ponieważ PEA skaluje nauczone współczynniki szumu bezpośrednio, zamiast powtarzać bramki, możliwy jest dowolny współczynnik szumu, także niecałkowity, a sam obwód nie jest modyfikowany: elementy programu współdzielą jeden obwód szablonowy i różnią się tylko szumem wstrzykiwanym przez samplex. Wstrzykiwany szum jest próbkowany raz na randomizację, więc błąd statystyczny PEA zależy od liczby randomizacji, a nie od liczby strzałów na randomizację; dlatego używamy wielu randomizacji z kilkoma strzałami każda, utrzymując łączną liczbę strzałów taką samą jak w programie bez mitygacji. Ważne, aby nie używać zbyt małej liczby randomizacji, ponieważ w przeciwnym razie możesz napotkać obciążenie próbkowania we wstrzykiwaniu szumu.

Całkowity narzut próbkowania rośnie liniowo wraz z liczbą współczynników szumu.

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

Na koniec uruchom program na zaszumionym symulatorze. Jego wyniki zasilają ekstrapolację w następnym kroku.

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

Ekstrapolacja do granicy zerowego szumu

Teraz, gdy masz dane dla każdej skali szumu, możesz dopasować funkcję do danych i wyekstrapolować z powrotem do wartości przy zerowym szumie, aby oszacować wartość oczekiwaną wolnej od szumu obserwabli. Przy szumie Pauliego obserwabla zanika wykładniczo wraz ze współczynnikiem szumu, więc jako pierwszy wybór wybierzesz ekstrapolator wykładniczy, a jeśli dopasowanie wykładnicze się nie powiedzie, wrócisz do liniowego. Aby zwizualizować ekstrapolację, narysuj uśrednioną po węzłach magnetyzację obliczoną dla każdego współczynnika szumu, dopasowanie wykładnicze oraz średnią po węzłach wyekstrapolowanych wartości oczekiwanych w punkcie zerowego szumu. Widzimy, że dopasowanie dla tej konkretnej obserwabli nieco zaniża wynik.

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

Porównanie zmitygowanych wartości oczekiwanych

Narysuj uśrednioną po węzłach magnetyzację każdej metody w całym zakresie kątów obrotu Isinga. PEA odtwarza dokładną magnetyzację w granicach swoich słupków błędów, podczas gdy wynik bez mitygacji jest obciążony i leży znacznie poniżej niej.

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