אקסטרפולציה לאפס רעש עם הגברת שגיאות הסתברותית

אקסטרפולציה לאפס רעש (ZNE) היא טכניקה לצמצום רעש קוונטי בחישובי ערכי תצפית של observable. היא פועלת על ידי הרצת מעגל קוונטי בסקאלות רעש הולכות וגדלות, חישוב ערכי תצפית ביחס לסקאלות הרעש האלה, ואקסטרפולציה כדי להעריך ערך תצפית באפס רעש. אפשר לממש בדרכים שונות גם את שלב הגברת הרעש וגם את שלב האקסטרפולציה. במדריך הזה תגברו את הרעש באמצעות הגברת שגיאות הסתברותית (PEA), כפי שהוצע ב-"Evidence for the utility of quantum computing before fault tolerance". כדי להגביר עם PEA, צריך קודם למדל את הרעש שמשפיע על השערים המסבכים של המעגל. אפשר לממש את ההגברה גם באמצעות קיפול שערים דיגיטלי, שלא דורש מודל רעש ומכוסה במדריך קיפול השערים. תלמדו את הרעש באמצעות qiskit-noise-learning ותבצעו אקסטרפולציה עם אסטרטגיית ההתאמה האקספוננציאלית המובנית.

זכרו שאף ש-ZNE יכולה להיות יעילה בצמצום רעש קוונטי בעלות דגימה נמוכה יותר משיטות אחרות, אין לה שום חסם על ההטיה.

זרימת העבודה היא כדלקמן:

1. מגדירים מעגל ו-observable
2. מחשבים ערכי תצפית מדויקים באמצעות סימולטור statevector של Qiskit
3. מחשבים ערכי תצפית רועשים, ללא צמצום רעש
4. לומדים את הרעש של שכבות השערים המסבכים באמצעות qiskit-noise-learning
5. מכינים תוכנית ZNE עם הגברת שגיאות הסתברותית ודוגמים מעגלים בסקאלות רעש שונות
6. מבצעים אקסטרפולציה לגבול של אפס רעש
7. משווים את ערכי התצפית אחרי צמצום הרעש

הגדרת מעגל ו-observable

במדריך הזה תחשבו את המגנטיזציה הממוצעת על פני האתרים $\langle Z \rangle_{\text{avg}} = \frac{1}{N} \sum_{i=1}^{N} \langle Z_i \rangle$ של מודל איזינג בשדה רוחבי (TFIM) בגרסת Trotter על שמונה קיוביטים בטופולוגיית טבעת. תסרקו כמה ערכים של עוצמת השדה הרוחבי ותשרטטו את המגנטיזציה כפונקציה של עוצמת השדה. בטופולוגיית טבעת עם מספר זוגי של קיוביטים, אפשר לפצל את שכבות השערים של שני קיוביטים לשתי שכבות ייחודיות, כפי שהצבעים האדום והירוק מראים במפת הצימוד שלמטה. חשוב לתכנן מעגלים עם קבוצה מינימלית של שכבות מסבכות ייחודיות אם רוצים ללמוד מודל רעש, כי את השכבות הייחודיות צריך ללמוד בנפרד. אם יש יותר מדי שכבות ייחודיות במעגל, למידת הרעש עלולה להפוך לבלתי אפשרית בגלל עלות הדגימה של למידת שכבות רבות.

# 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

חישוב ערכי תצפית מדויקים באמצעות סימולטור statevector של Qiskit

השתמשו במחלקה Statevector של Qiskit כדי לחשב את ערך התצפית המדויק של $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 מהמעגל, מה-observables ומערכי הפרמטרים, ועוטפת את המעגל בתיבות (boxes) לפי אפשרויות ה-boxing שניתנו. הגדירו את אפשרות ה-boxing ‏measure_annotations לערך "change_basis", שמסובב את המדידות לבסיסים הנדרשים, והוסיפו add_tags: "unique_box" כדי שהסימולטור יוכל למפות מודל רעש על התיבות של המעגל. התוכנית עוברת Pauli twirling כברירת מחדל, מה שעל חומרה הופך רעש קוהרנטי לרעש Pauli שאליו שיטות צמצום הרעש מתייחסות. num_randomizations קובע כמה מופעים עם twirling של כל מעגל נדגמים, ו-shots_per_randomization קובע כמה shots מקבל כל מופע, כך שהמכפלה שלהם היא תקציב ה-shots לכל מעגל.

משתמש עם גישה ל-QPU יריץ את התוכנית באמצעות Executor שנבנה מ-Backend של QiskitRuntimeService. כאן תעביר AerSimulator ל-Executor, שמריץ את התוכנית באופן מקומי, ותשתמש באפשרות הסימולטור layer_noise_model שלו כדי להזריק מודל רעש Pauli-Lindblad אקראי בכל שכבה של שערי details שזירה (entangling). המודל נוצר באמצעות 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 מתחיל בלימוד מודל הרעש המסובב (twirled) של כל שכבה של שערי שזירה במעגל. אנחנו מניחים מודל רעש מרקובי, ולכן צריך ללמוד רק את מודל הרעש של כל שכבה ייחודית של שערי שזירה, כי מניחים שעומק המעגל לא משפיע על ערוץ הרעש של שכבת שזירה נתונה. ככל שמודל הרעש שנלמד מדויק יותר, כך ההגברה תהיה מדויקת יותר.

אנחנו משתמשים ב-qiskit-noise-learning כדי ללמוד את הרעש. הלומד חוזר על כל שכבה fragment_depths פעמים במעגלי הלמידה שלו, מסובב (twirl) כל מעגל num_randomizations פעמים, ומתאים את דעיכת הנאמנויות (fidelities) הנמדדות למודל Pauli-Lindblad. הקוד הבא משרטט את מודל הרעש שנלמד, כשתי השכבות ממוזגות למפה אחת. כל Qubit הוא עוגה של שיעורי X, Y ו-Z ממשקל 1 שסוכמו על פני שתי השכבות, וכל צלע מראה את השיעורים ממשקל 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 מכפיל ישירות את שיעורי הרעש שנלמדו במקום לחזור על שערים, אפשר להשתמש בכל גורם רעש, כולל כאלה שאינם שלמים, והמעגל עצמו לא משתנה: פריטי התוכנית חולקים מעגל תבנית אחד ומשתנים רק ברעש שה-samplex מזריק. הרעש המוזרק נדגם פעם אחת לכל randomization, ולכן השגיאה הסטטיסטית של PEA נקבעת לפי מספר ה-randomizations ולא לפי מספר ה-shots לכל randomization; לכן אנחנו משתמשים בהרבה randomizations עם מעט shots בכל אחד, ומשאירים את המספר הכולל של ה-shots זהה לזה של התוכנית ללא מיתון. חשוב לא להשתמש במעט מדי randomizations, אחרת אתה עלול להיתקל בהטיית דגימה בהזרקת הרעש.

תקורת הדגימה הכוללת גדלה באופן ליניארי עם מספר גורמי הרעש.

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

אקסטרפולציה לגבול של אפס רעש

עכשיו כשיש לך נתונים לכל קנה מידה של רעש, אפשר להתאים פונקציה לנתונים ולבצע אקסטרפולציה חזרה לערך של אפס רעש כדי להעריך ערך תצפית של האובזרבבל נטול הרעש. תחת רעש Pauli, אובזרבבל דועך אקספוננציאלית עם גורם הרעש, ולכן תבחר באקסטרפולטור האקספוננציאלי כבחירה ראשונה, ותחזור לליניארי אם ההתאמה האקספוננציאלית נכשלת. כדי להמחיש את האקסטרפולציה, שרטט את המגנטיזציה הממוצעת על פני האתרים שחושבה בכל גורם רעש, את ההתאמה האקספוננציאלית, ואת הממוצע על פני האתרים של ערכי התצפית שעברו אקסטרפולציה בנקודת אפס הרעש. אנחנו רואים שההתאמה נמוכה מעט מדי עבור האובזרבבל הספציפי הזה.

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

השוואת ערכי התצפית אחרי מיתון

שרטט את המגנטיזציה הממוצעת על פני האתרים של כל שיטה, על פני סריקת זוויות הסיבוב של Ising. 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()