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

אקסטרפולציה לאפס רעש (ZNE) היא טכניקה לצמצום רעש קוונטי בחישובי ערכי תצפית של observable. היא פועלת על ידי הרצת מעגל קוונטי בסקאלות רעש הולכות וגדלות, חישוב ערכי תצפית ביחס לסקאלות הרעש האלה, ואקסטרפולציה כדי להעריך ערך תצפית באפס רעש. אפשר לממש בדרכים שונות גם את שלב הגברת הרעש וגם את שלב האקסטרפולציה. המדריך הזה מגביר את הרעש באמצעות "קיפול שערים דיגיטלי", כלומר שערים של שני קיוביטים מוחלפים ברצפים שקולים של השער וההופכי שלו. לדוגמה, החלפת אופרטור אוניטרי $U$ ב-$U U^† U$ תניב גורם הגברת רעש של 3. אפשר לממש את ההגברה גם באמצעות הגברת שגיאות הסתברותית (PEA), שדורשת שלב נוסף של למידת מודל רעש ומכוסה במדריך הגברת שגיאות הסתברותית. המדריך הזה מבצע אקסטרפולציה בעזרת כמה אסטרטגיות התאמה מובנות (לינארית ואקספוננציאלית), וגם מדגים איך לממש פונקציית התאמה משלכם.

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

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

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

הגדרת מעגל ו-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", "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

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

השתמשו במחלקה Statevector של Qiskit כדי לחשב את ערך התצפית המדויק של $Z_i$ בכל אתר, עבור כל ערך של $\theta$.

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

חישוב ערכי תצפית רועשים, ללא צמצום רעש

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 אקראי על כל שכבה של שערים מסבכים, שנוצר באמצעות random_pauli_lindblad_noise מ-qiskit-mitigation עם שגיאה לשער בשכבה (EPLG) של 0.003, אופיינית ל-QPUs של 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()

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

הכינו תוכנית ZNE שמכילה מעגלים בסקאלות רעש הולכות וגדלות. כל גורם רעש מייצג את הכפולה של השערים המסבכים שהמעגל מכיל (גורם רעש 3.0 מחזיק פי שלושה שערים מסבכים מהמעגל המקורי). המעגל של כל גורם רעש עובר twirling עם אותו מספר randomizations ו-shots כמו התוכנית ללא צמצום רעש.

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
)

שרטטו את צעד Trotter הראשון של התבנית אחרי twirling עם גורם רעש 1 ועם גורם רעש 3 כדי להמחיש את ההגברה. עבור גורם רעש 3.0, כל שכבה של המעגל מוחלפת בשכבה, בהופכית שלה, ושוב בשכבה.

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
)

לבסוף, הריצו את התוכנית על הסימולטור הרועש. התוצאות שלה מזינות את האקסטרפולציה בשלב הבא.

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

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

עכשיו, כשיש לכם נתונים לכל סקאלת רעש, אפשר להתאים פונקציה לנתונים ולבצע אקסטרפולציה חזרה לערך באפס רעש כדי להעריך ערך תצפית של ה-observable נטול הרעש. יש דרכים רבות להתאים נתונים, והבחירה הנכונה של פונקציה תלויה באופי הנתונים שנדגמו. GateFolding.postprocess מספקת כמה פונקציות אקסטרפולציה מובנות, בהן linear, exponential ו-double_exponential. המתודה postprocess מקבלת רשימה של שיטות אקסטרפולציה, מעריכה אותן לפי הסדר, ומחזירה את הערך הראשון שהוא תקף פיזיקלית ובעל שונות סופית. כלומר, הסדר שבו אתם מציינים את המאקסטרפולטורים משפיע מאוד על הבחירה. העברת fallback כמאקסטרפולטור אומרת שאם לא נמצאו תשובות תקפות, ייבחר ויוחזר ערך התצפית של גורם הרעש הקטן ביותר. אם fallback לא הועבר ושום מאקסטרפולטור לא מניב ערך תקף פיזיקלית, מוחזרת התוצאה עם ה-std הקטן ביותר. תחת רעש Pauli observable דועך אקספוננציאלית עם גורם הרעש, ולכן עם 30 שכבות (15 צעדי Trotter, 2 שכבות לצעד) של שערים מסבכים, המאקסטרפולטור האקספוננציאלי הוא הבחירה הטבעית; אפשר גם לבצע אקסטרפולציה לינארית, בחירה נפוצה כשהדעיכה קטנה.

כדי להשתמש בפונקציית התאמה מותאמת אישית, העבירו את המאקסטרפולטור "custom" יחד עם tuple בשם custom_fit שמכיל את פונקציית ההתאמה, רשימה עם ניחוש התחלתי לכל אחד מהפרמטרים שלה, וגבולות. אלה מועברים ל-scipy.optimize.curve_fit. אפשר לחזור על העיבוד המשלים עם מאקסטרפולטורים שונים בלי להריץ מחדש את התוכנית. כאן משתמשים בפונקציית התאמה מותאמת אישית של שורש ריבועי. כדי להמחיש את האקסטרפולציה, שרטטו את המגנטיזציה הממוצעת על פני האתרים שחושבה בכל גורם רעש עבור $\theta = 0.15\pi$, את שלוש ההתאמות, ואת הממוצע על פני האתרים של ערכי התצפית שעברו אקסטרפולציה בנקודת אפס הרעש. שלוש נקודות הנתונים נראות כמעט לינאריות, ועדיין האקסטרפולציה הלינארית נוחתת באופן ברור מתחת לערך המדויק. כפי שחזינו, ההתאמה האקספוננציאלית מעריכה בצורה היעילה ביותר את ערך התצפית נטול הרעש.

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

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

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

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