Екстраполяція до нульового шуму з ймовірнісним підсиленням помилок

Екстраполяція до нульового шуму (ZNE) — це техніка пом'якшення квантового шуму в обчисленнях очікуваних значень спостережуваних. Вона працює шляхом виконання квантової схеми при зростаючих масштабах шуму, взяття очікуваних значень для цих масштабів і екстраполяції для оцінки очікуваного значення при нульовому шумі. Як підсилення шуму, так і етапи екстраполяції можна реалізувати по-різному. У цьому посібнику ти підсиліш шум за допомогою ймовірнісного підсилення помилок (PEA), як запропоновано в "Evidence for the utility of quantum computing before fault tolerance". Щоб підсилювати за допомогою PEA, спочатку потрібно змоделювати шум, що впливає на заплутуючі вентилі схеми. Підсилення також можна реалізувати за допомогою цифрового згортання вентилів, яке не потребує моделі шуму і розглянуте в посібнику зі згортання вентилів. Ти навчиш шум за допомогою qiskit-noise-learning і виконаєш екстраполяцію вбудованою експоненціальною стратегією апроксимації.

Пам'ятай: хоча ZNE може ефективно пом'якшувати квантовий шум за нижчої вартості вибірки, ніж інші методи, вона не дає жодних меж для зміщення.

Робочий процес такий:

1. Задати схему та спостережувану
2. Обчислити точні очікувані значення за допомогою симулятора вектора стану Qiskit
3. Обчислити зашумлені непом'якшені очікувані значення
4. Навчити шум заплутуючих шарів за допомогою qiskit-noise-learning
5. Підготувати програму ZNE з ймовірнісним підсиленням помилок та вибіркою схем різних масштабів шуму
6. Екстраполювати до межі нульового шуму
7. Порівняти пом'якшені очікувані значення

Задай схему та спостережувану

У цьому посібнику ти обчислиш усереднену за вузлами намагніченість $\langle Z \rangle_{\text{avg}} = \frac{1}{N} \sum_{i=1}^{N} \langle Z_i \rangle$ тротеризованої моделі Ізинга з поперечним полем (TFIM) на восьми кубітах у кільцевій топології. Ти переберешеш кілька значень напруженості поперечного поля та побудуєш графік намагніченості як функції напруженості поля. У кільцевій топології з парним числом кубітів шари двокубітних вентилів можна розбити на два унікальні шари, як показують червоний і зелений кольори на карті зв'язності нижче. Проєктувати схеми з мінімальним набором унікальних заплутуючих шарів важливо, якщо ти хочеш навчити модель шуму, адже унікальні шари потрібно навчати незалежно. Якщо у твоїй схемі забагато унікальних шарів, навчання шуму може стати нездійсненним через вартість вибірки для навчання багатьох шарів.

# 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

Обчисли точні очікувані значення за допомогою симулятора вектора стану Qiskit

Використай клас Qiskit Statevector, щоб обчислити точне очікуване значення $Z_i$ на кожному вузлі для кожного значення $\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()

Обчисли зашумлені непом'якшені очікувані значення

QuantumProgram містить квантову схему та samplex у форматі, придатному для примітива Executor. MitigationTask.prepare будує QuantumProgram зі схеми, спостережуваних і значень параметрів та огортає схему в бокси із заданими параметрами боксування. Задай параметр боксування measure_annotations рівним "change_basis", що повертає вимірювання в потрібні базиси, і додай add_tags: "unique_box", щоб симулятор міг накласти модель шуму на бокси схеми. Програма за замовчуванням піддається паулівському твірлінгу, що на апаратному забезпеченні перетворює когерентний шум на паулівський шум, який передбачають методи пом'якшення. num_randomizations задає, скільки твірлованих екземплярів кожної схеми береться у вибірку, а shots_per_randomization задає, скільки шотів отримує кожен екземпляр, тож їх добуток — це бюджет шотів на схему.

Користувач з доступом до QPU виконав би програму за допомогою Executor, створеного з бекенда QiskitRuntimeService. Тут ти передаси AerSimulator до Executor, який виконує програму локально, і за допомогою його опції симулятора layer_noise_model додаси випадкову модель шуму Паулі-Ліндбладу на кожен шар заплутуючих гейтів. Модель згенерована за допомогою random_pauli_lindblad_noise з qiskit-mitigation із помилкою на шаровий гейт (EPLG) 0.003, типовою для сучасних 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()

Вивчи шум заплутуючих шарів за допомогою qiskit-noise-learning

PEA починається з вивчення твірльованої моделі шуму кожного шару заплутуючих гейтів у схемі. Ми припускаємо марковську модель шуму, тому слід вивчати лише модель шуму кожного унікального шару заплутуючих гейтів, адже вважається, що глибина схеми не впливає на канал шуму заданого заплутуючого шару. Що точніше вивчена модель шуму, то точнішим буде підсилення.

Для вивчення шуму ми використовуємо qiskit-noise-learning. Вивчальник повторює кожен шар fragment_depths разів у своїх навчальних схемах, твірлить кожну схему num_randomizations разів і підганяє спад виміряних точностей до моделі Паулі-Ліндбладу. Наведений нижче код малює вивчену модель шуму, де два шари об'єднано в одну карту. Кожен кубіт зображено як діаграму-круг із його ваговими-1 швидкостями X, Y і Z, підсумованими за двома шарами, а кожне ребро показує вивчені вагові-2 швидкості; швидкість 3e-4 означає, що помилка виникає приблизно три рази на кожні 10 000 використань шару.

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

Підготуй програму ZNE з імовірнісним підсиленням помилок і семпли схем із різними масштабами шуму

Тут ти підготуєш програму ZNE, що містить схеми зі зростаючими масштабами шуму. У PEA кожен коефіцієнт шуму означає кратність вивченого шуму, який відчуває схема: коефіцієнт шуму 2.0 додає одну додаткову копію вивченого каналу шуму кожного шару після цього шару, понад шум, який шар уже має. Оскільки PEA масштабує вивчені швидкості шуму безпосередньо, а не повторює гейти, можливий будь-який коефіцієнт шуму, зокрема нецілий, а сама схема не змінюється: елементи програми спільно використовують один шаблонний circuit і відрізняються лише шумом, який додає samplex. Доданий шум семплується один раз на рандомізацію, тому статистична похибка PEA контролюється кількістю рандомізацій, а не кількістю шотів на рандомізацію; тому ми використовуємо багато рандомізацій з кількома шотами кожна, зберігаючи загальну кількість шотів такою ж, як для програми без пом'якшення. Важливо не використовувати надто мало рандомізацій, інакше можна зіткнутися зі зміщенням семплування під час додавання шуму.

Загальні витрати на семплування зростають лінійно з кількістю коефіцієнтів шуму.

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

Нарешті, запусти програму на шумному симуляторі. Його результати використовуються для екстраполяції на наступному кроці.

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

Екстраполюй до межі нульового шуму

Тепер, коли ти маєш дані для кожного масштабу шуму, можна підігнати функцію до даних і екстраполювати назад до значення при нульовому шумі, щоб оцінити очікуване значення безшумної observable. За шуму Паулі observable спадає експоненційно з коефіцієнтом шуму, тому першим вибором буде експоненційний екстраполятор, а якщо експоненційна підгонка не вдається, слід повернутися до лінійного. Щоб візуалізувати екстраполяцію, побудуй графік усередненої за вузлами намагніченості, обчисленої при кожному коефіцієнті шуму, експоненційної підгонки та усередненого за вузлами екстрапольованого очікуваного значення в точці нульового шуму. Ми бачимо, що для цієї конкретної observable підгонка трохи занижує значення.

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

Порівняй пом'якшені очікувані значення

Побудуй графік усередненої за вузлами намагніченості кожного методу для набору кутів обертання Ізінга. PEA відновлює точну намагніченість у межах своїх смуг похибок, тоді як результат без пом'якшення значно зміщений униз від неї.

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