סימולציה של מערכות קוונטיות רועשות באמצעות Pauli propagation

במדריך הזה נשתמש בחבילה pauli-prop כדי לסמלץ בצורה קלאסית את הדינמיקה בזמן של מודל Ising עם שדה רוחבי (TFIM) רועש על תשעה qubit ברשת ריבועית של 3x3. אנחנו משתמשים בהוראות PauliLindbladError כדי להגדיר ערוץ רעש, $\Lambda$, שפועל על קבוצה של שכבות סבך, $\mathcal{U}$. לאחר מכן אנחנו מקדמים את האובזרבבל, $O$, אחורה דרך המעגל הרועש ומעריכים ערכי תצפית עבור מגוון מודלי רעש, וכן עבור המקרה נטול הרעש.

כשהאובזרבבל מקודם אחורה דרך המעגל, כל ערוץ רעש, $\Lambda_k$, המשויך לשכבת סבך, $\mathcal{U}_k$, מדכא את איברי Pauli ב-$O$ שמתחלפים בסימן הפוך (anti-commute) עם יוצרי Pauli-Lindblad שלו. ליתר דיוק, אם $G_{k,i}$ הוא יוצר Pauli של $\Lambda_k$ עם קצב, $\gamma_{k,i}$, אז איבר Pauli, $P$, ב-$O$ הופך כך: $c_P \mapsto c_P e^{-2\gamma_{k,i}} \quad \text{if } \{P, G_{k,i}\}=0$, כאשר $c_P$ הוא המקדם של $P$. ברגע ש-$O$ קודם אל תחילת המעגל, אפשר לחשב בקלות את ערך התצפית ביחס למצב האפס, $|0\rangle^{\otimes N}$, על ידי סכימת המקדמים של כל איבר אלכסוני ב-$O$ (איברים שמכילים $Z$ או $I$ בכל ה-qubit).

זרימת עבודה:

• מגדירים את רשת ה-TFIM, ומשתמשים בצביעת קשתות כדי לזהות קבוצה מינימלית של שכבות סבך
• יוצרים מודלי רעש סינתטיים, $\Lambda_k$, לכל שכבת סבך ייחודית, $U_k$
  • יוצרים מודלי רעש בסקאלות שונות כדי לחקור את ההשפעה של רעש השערים על המערכת
• יוצרים מעגלים קוונטיים נטולי רעש ורועשים עבור העומקים וסקאלות הרעש שמעניינים אותנו
  • במעגלים הרועשים מוכנסות הוראות PauliLindbladError לפני כל שכבת סבך
• משתמשים ב-Pauli propagation כדי לסמלץ ערכי תצפית מדויקים של המערכת בעומקים שונים
  • עבור תשעה qubit, זה נעשה בכך שנותנים ל-$O$ לגדול עד $4^9$ איברים, שמכסים את כל מרחב Pauli
• משתמשים ב-Pauli propagation כדי לסמלץ ערכי תצפית רועשים
• מתבוננים איך הגדלת רעש השערים פוגעת בדיוק של המודל הקוונטי

יצירת רשת ריבועית של 3x3 ומציאת צביעה ב-4 צבעים של הקשתות

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

זיהוי קבוצה מינימלית של שכבות סבך ייחודיות חשוב לרוב ליישום יעיל של פרוטוקולי למידת רעש, כי צריך ללמוד את הרעש של כל שכבה בנפרד. ככל שצריך ללמוד יותר שכבות, כך צריך לקחת יותר shots מה-QPU. בהדגמה הזאת אנחנו משתמשים במידע על השכבות כדי לבנות מעגלים רועשים ולהזריק הוראות PauliLindbladError מתוך qiskit-aer לפני כל שכבת סבך, כדי לדמות את רעש השערים של ה-QPU.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "pauli_prop": "pauli-prop", "qiskit": "qiskit", "qiskit_addon_utils": "qiskit-addon-utils", "qiskit_aer": "qiskit-aer", "rustworkx": "rustworkx"}
_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")

from collections import defaultdict

import numpy as np
from qiskit.transpiler import CouplingMap
from qiskit_addon_utils.coloring import auto_color_edges

# Define rectangular square-lattice on 20 qubits
num_rows = 3
num_cols = 3
num_qubits = num_rows * num_cols

coupling_map = CouplingMap.from_grid(
    num_rows=num_rows, num_columns=num_cols, bidirectional=False
)

# Create mapping from color to edge list
coloring = auto_color_edges(coupling_map.get_edges())
color_to_edge = defaultdict(list)
for edge, color in coloring.items():
    color_to_edge[color].append(edge)

from rustworkx import PyDiGraph
from rustworkx.visualization import graphviz_draw

# Inspect graph coupling and unique entangling layers
print(
    f"The circuit will have {num_qubits} qubits and {len(color_to_edge)} unique entangling layers."
)
sq_lattice = PyDiGraph()
sq_lattice.extend_from_weighted_edge_list(
    [
        (source, target, color)
        for ((source, target), color) in coloring.items()
    ]
)

def color_edge_4color(edge):
    color_dict = {0: "red", 1: "green", 2: "blue", 3: "orange"}
    return {"color": color_dict[edge]}

graphviz_draw(sq_lattice, edge_attr_fn=color_edge_4color, method="neato")

The circuit will have 9 qubits and 4 unique entangling layers.

יצירת מודלי רעש סינתטיים

לפני שניצור את המעגלים הקוונטיים, אנחנו יוצרים מודל רעש (מופע של PauliLindbladError) לכל אחת משכבות הסבך. נשלב אותם מאוחר יותר כהוראות במעגלים הקוונטיים שלנו. לכל שכבה אנחנו יוצרים ערוצי רעש בסקאלות שונות. ליתר דיוק, אנחנו יוצרים מודלי רעש עם Error Per Layered Gate (EPLG) של בערך .0004, .0008, .0012, .0016, ו-.002.

from qiskit.quantum_info import SparsePauliOp, pauli_basis
from qiskit_aer.noise import PauliLindbladError

# Pauli-Lindblad noise parameters
seed = 1764
target_EPLGs = [0.0004, 0.0008, 0.0012, 0.0016, 0.002]

def generate_random_pauli_lindblad_noise(
    edges,
    num_qubits: int | None = None,
    noise_scale: float = 1e-3,
    seed: int | None = None,
) -> PauliLindbladError:
    """Generate random Pauli-Lindblad noise over the full Pauli basis."""
    if num_qubits is None:
        num_qubits = np.max(edges)

    basis_paulis = [p for p in pauli_basis(2) if np.sum(p.x + p.z)]
    basis_paulis = SparsePauliOp.from_sparse_list(
        [
            (pauli.to_label(), edge, 1)
            for pauli in basis_paulis
            for edge in edges
        ],
        num_qubits=num_qubits,
    )
    basis_paulis = basis_paulis.simplify()
    basis_paulis = basis_paulis.paulis

    rng = np.random.default_rng(seed=seed)
    rates = rng.random(len(basis_paulis)) * noise_scale

    return PauliLindbladError(generators=basis_paulis, rates=rates)

num_generators = (
    (num_rows * num_cols)
    + (num_rows - 1) * num_cols
    + num_rows * (num_cols - 1)
)
noise_scales = [
    EPLG * (num_rows * num_cols) / num_generators for EPLG in target_EPLGs
]
noise_models_per_EPLG = [
    [
        generate_random_pauli_lindblad_noise(
            color_to_edge[color],
            num_qubits=num_qubits,
            noise_scale=noise_scale,
            seed=seed,
        )
        for color in range(len(color_to_edge))
    ]
    for noise_scale in noise_scales
]

יצירת המעגלים הקוונטיים

בהדגמה הזאת אנחנו מסמלצים את הדינמיקה בזמן של מודל Ising עם שדה רוחבי (TFIM) עבור מספר גדל והולך של צעדי Trotter (1-10 צעדים). לכל אחד מ-10 עומקי המעגל, אנחנו מסמלצים את ההשפעה של רעש השערים, בהינתן מודלי רעש בסקאלות שונות (EPLGs = .0004, .0008, .0012, .0016, .002). הרעש מוכנס אל ה-QuantumCircuit כהוראת PauliLindbladError מתוך Qiskit Aer. ההמילטוניאן שבו נתמקד הוא:

$H = -J\sum\limits_{\langle i,j \rangle} Z_iZ_j + h\sum\limits_iX_i$

כאשר $J>0$ מתאר את הצימוד בין ספינים שכנים, $i<j$, ו-$h$ הוא השדה הרוחבי הגלובלי.

כאן אנחנו מממשים את ההמילטוניאן המתפתח בזמן על פני סקאלות שונות של זמן ורעש. אנחנו יוצרים 60 מעגלים בסך הכול: 10 מעגלים נטולי רעש שמשתנים בעומק Trotter, ו-50 מעגלים רועשים עבור 10 עומקי Trotter על פני חמש סקאלות רעש. בהינתן גרף קישוריות, המודל מוגדר על ידי כמה משתנים:

• num_steps: מספר צעדי Trotter
• J: עוצמת הצימוד של אתרים מחוברים
• h: עוצמת השדה המגנטי החיצוני
• dt: השינוי בזמן בצעד Trotter אחד
• initial_state_angle: עירור התחלתי, $R_y(\theta)$, שמופעל באופן אחיד על כל ה-qubit

from typing import Any

from qiskit import QuantumCircuit

# Ising model parameters
num_steps = 10
J = -1.0
dt = 0.25 / abs(J)
h = 2.0 * abs(J)
initial_state_angle = np.pi / 18.0
rx_angle = 2.0 * h * dt
rzz_angle = 2.0 * J * dt

def generate_ising_circuit(
    num_qubits: int,
    num_steps: int,
    rx_angle: float,
    rzz_angle: float,
    coloring: dict[Any, list[tuple[int, int]]],
    layer_noise_models: list[PauliLindbladError] | None = None,
    initial_state_angle: float | None = None,
) -> QuantumCircuit:
    """Generate a quantum circuit implementing a transverse-field Ising model"""
    qc = QuantumCircuit(num_qubits)
    if initial_state_angle:
        qc.ry(initial_state_angle, range(num_qubits))
    qc.rx(rx_angle / 2, range(num_qubits))
    for i in range(num_steps):
        for j, layer in enumerate(coloring):
            edges = coloring[layer]
            if layer_noise_models:
                qc.append(layer_noise_models[j], qargs=range(num_qubits))
            for edge in edges:
                qc.rzz(rzz_angle, *edge)
        if i == num_steps - 1:
            qc.rx(rx_angle / 2, range(num_qubits))
        else:
            qc.rx(rx_angle, range(num_qubits))
    return qc

# Create the noiseless and noisy circuits
noiseless_circs = []
noisy_circs = []
for steps in range(1, num_steps + 1):
    noiseless_circs.append(
        generate_ising_circuit(
            num_qubits,
            steps,
            rx_angle,
            rzz_angle,
            color_to_edge,
            initial_state_angle=initial_state_angle,
        )
    )
    noisy_circs_per_step = []
    for noise_models in noise_models_per_EPLG:
        noisy_circs_per_step.append(
            generate_ising_circuit(
                num_qubits,
                steps,
                rx_angle,
                rzz_angle,
                color_to_edge,
                layer_noise_models=noise_models,
                initial_state_angle=initial_state_angle,
            )
        )
    noisy_circs.append(noisy_circs_per_step)
print(
    f"{num_steps} noiseless and {num_steps * len(target_EPLGs)} noisy Trotter circuits generated. {num_steps} different depths across {len(target_EPLGs)} different noise models"
)
print("\nBelow: Initial state and one noisy Trotter step.")
noisy_circs[0][0].draw("mpl", fold=-1)

10 noiseless and 50 noisy Trotter circuits generated. 10 different depths across 5 different noise models

Below: Initial state and one noisy Trotter step.

הגדרת האובזרבבל והרצת הסימולציות

בהדגמה הזאת אנחנו מסמלצים ערכי תצפית של מתאם שני האתרים הממוצע:

$\langle O \rangle = \langle Z_{tot}^2(s) \rangle = \frac{1}{N^2}\sum \langle \Psi(\theta)|(\mathscr{U}^{\dagger})^sZ_jZ_k(\mathscr{U})^s|\Psi(\theta) \rangle$

כאשר $\Psi(\theta)$ מתאים לסיבוב אחיד $R_y(\theta)$ על כל ה-qubit, $\mathscr{U}^s$ מתאר $s$ שכבות Trotter, ו-$(j,k)$ מאנדקסים את כל זוגות הקודקודים המחוברים ברשת.

לבסוף, אנחנו משתמשים ב-pauli_prop כדי לסמלץ ערכי תצפית של האובזרבבל עבור כל אחד מהמעגלים. בהדגמה הזאת עם תשעה qubit אנחנו מבצעים את כל הסימולציות במדויק. לא תתבצע שום קטיעה ב-Pauli propagation; לכן אפשר לייחס את ההבדלים בערכי התצפית בין מודלי הרעש השונים כולם לשגיאת השערים. תהליך הסימולציה מתבצע בארבעה צעדים:

• מפתחים את שערי ה-Clifford במעגל אל תחילת המעגל באמצעות pauli_prop.evolve_through_cliffords
• מקדמים את האובזרבבל דרך החלק של המעגל שאינו Clifford באמצעות pauli_prop.propagate_through_circuit
  • אנחנו מבצעים סימולציות מדויקות בכך שנותנים לאובזרבבל לגדול עד גודל מרחב Pauli המלא, $4^9$
• מקדמים את האובזרבבל המפותח דרך חלק ה-Clifford של המעגל באמצעות SparsePauliOp.evolve של Qiskit
• מעריכים את ערך התצפית ביחס למצב האפס, $|0\rangle^{\otimes N}$, על ידי סכימת המקדמים של כל איבר אלכסוני ב-$O$ (איברים שמכילים $Z$ או $I$ בכל ה-qubit)

import time

from pauli_prop import evolve_through_cliffords, propagate_through_circuit
from qiskit.quantum_info import Pauli

# Average ZZ-correlator observable
id_pauli = Pauli("I" * num_qubits)
observable = 2 * SparsePauliOp(
    [
        id_pauli.dot(Pauli("ZZ"), [i, j])
        for i in range(num_qubits)
        for j in range(i + 1, num_qubits)
    ]
)
observable /= num_qubits**2

# Pauli propagation parameters
max_terms = 4**num_qubits  # Exact propagation
atol = 1e-12

# Run simulations
exact_evs = []
noisy_evs = [[] for _ in range(len(target_EPLGs))]
st = time.perf_counter()
for i, noiseless_circ in enumerate(noiseless_circs):
    cliff, non_cliff = evolve_through_cliffords(noiseless_circ)
    evolved_obs = propagate_through_circuit(
        observable, non_cliff, max_terms=max_terms, atol=atol, frame="h"
    )[0]
    evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame="h")
    exact_evs.append(
        float(evolved_obs.coeffs[~evolved_obs.paulis.x.any(axis=1)].sum())
    )
    for j in range(len(target_EPLGs)):
        noisy_circ = noisy_circs[i][j]
        cliff, non_cliff = evolve_through_cliffords(noisy_circ)
        evolved_obs = propagate_through_circuit(
            observable, non_cliff, max_terms=max_terms, atol=1e-12, frame="h"
        )[0]
        evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame="h")
        noisy_evs[j].append(
            float(evolved_obs.coeffs[~evolved_obs.paulis.x.any(axis=1)].sum())
        )
print(
    f"Ran {len(noiseless_circs)} noiseless and {len(target_EPLGs) * num_steps} noisy simulations in {int(time.perf_counter() - st)}s."
)

Ran 10 noiseless and 50 noisy simulations in 103s.

התבוננות בהשפעה של שגיאת שערים על המודל

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

import matplotlib.pyplot as plt

xs = range(1, num_steps + 1)
plt.plot(xs, exact_evs, label="Noiseless", color="black", marker="o")
colors = [".3", ".4", ".5", ".6", ".7"]
for i, evs in enumerate(noisy_evs):
    plt.plot(
        xs,
        evs,
        label=f"{target_EPLGs[i]} EPLG",
        linestyle="--",
        color=colors[i],
        marker="o",
    )
plt.xlabel("# Trotter steps")
plt.ylabel(r"$\langle Z_{tot}^2 \rangle$")
plt.legend()
plt.show()