נוסחאות מכפלת-רב להפחתת שגיאת Trotter

הערכת שימוש: ארבע דקות על מעבד Heron r2 (הערה: זוהי הערכה בלבד. זמן הריצה שלך עשוי להשתנות.)

תוצרי למידה

• כיצד נוסחאות מכפלה מרובות (Multi-Product Formulas - MPFs) מפחיתות את שגיאת Trotter בסימולציית המילטוניאן על ידי שילוב ערכי תוחלת מכמה Circuit-ים רדודים

• מתי MPFs מועילות מעל נוסחאות מוצר סטנדרטיות ומתי הן לא הכלי הנכון

• כיצד לחשב מקדמי MPF סטטיים ודינמיים באמצעות החבילה qiskit_addon_mpf

• כיצד להריץ תהליך עבודה של MPF מקצה לקצה על חומרת IBM Quantum®, כולל טרנספילציה, הקלת שגיאות ועיבוד לאחר מכן

דרישות קדם

• שיטות קומפילציה עבור מעגלי סימולציית המילטוניאן — מציג Circuit-ים של Trotter (נוסחת מוצר) ב-Qiskit.

• נוסחאות מוצר ב-Qiskit, בפרט מחלקות הסינתזה SuzukiTrotter ו-LieTrotter.

• פרימיטיבים של Qiskit וממשק ה-Estimator.

רקע

מהן נוסחאות מכפלה מרובות?

בעת סימולציית מערכות קוונטיות על מחשב קוונטי, משימה מרכזית היא לקרב את אופרטור התפתחות הזמן $e^{-iHt}$ עבור המילטוניאן $H$. הגישה הסטנדרטית משתמשת בנוסחאות מוצר (product formulas - PFs), הידועות גם כפירוקי Trotter-Suzuki. אלה מפרקות $H = \sum_{a=1}^d F_a$ לאיברים שהאוניטרים הבודדים שלהם $e^{-iF_a t}$ יעילים ליישום, ולאחר מכן מקרבות את ההתפתחות המלאה כמכפלה מסודרת של האוניטרים הפשוטים יותר הללו.

נוסחת המוצר מהסדר הראשון (Lie-Trotter) היא:

$$S_1(t) := \prod_{a=1}^d e^{-i F_a t},$$

אשר גורמת לשגיאה ריבועית: $S_1(t) = e^{-iHt} + \mathcal{O}(t^2)$. נוסחאות סימטריות מסדר גבוה יותר $S_{2\chi}(t)$, כאשר $\chi$ מציין את הסדר של נוסחת המוצר הסימטרית (ראה מקור [1]), מתכנסות מהר יותר כ-$e^{-iHt} + \mathcal{O}(t^{2\chi+1})$, אך במחיר של Circuit-ים עמוקים יותר לכל צעד.

כדי להפחית את השגיאה בסדר $\chi$ קבוע, בדרך כלל מחלקים את זמן ההתפתחות הכולל $t$ ל-$k$ צעדי Trotter קטנים יותר. כל צעד מקרב את $e^{-iHt/k}$ על ידי נוסחת מוצר והצעדים משורשרים יחד:

$$e^{-iHt} \approx \left[S_{2\chi}(t/k)\right]^k.$$

עבור נוסחה סימטרית מסדר $2\chi$, שגיאת Trotter השיורית מתנהגת אז כ-$\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$. כך שהגדלת $k$ מדכאת במהירות את שגיאת Trotter — אך היא גם מעמיקה את ה-Circuit באופן ליניארי, ועל חומרה רועשת זה אומר יותר רעש Gate מצטבר. המתח הזה בין שגיאת Trotter (מעדיפה $k$ גדול יותר) לבין רעש חומרה (מעדיף $k$ קטן יותר) הוא בדיוק מה שנוסחאות מכפלה מרובות נועדו לפתור. שים לב ש-MPFs עוסקות בשילוב תוצאות מבחירות שונות של $k$ בסדר $\chi$ קבוע — הן לא משנות את הסדר של נוסחת המכפלה הבסיסית.

נוסחאות מכפלה מרובות (MPFs) [1] בונות צירוף לינארי משוקלל של ערכי תוחלת המתקבלים מכמה Circuit-ים רדודים יותר של Trotter, כל אחד משתמש במספר שונה של צעדי Trotter $k_1, k_2, \ldots, k_r$ (קבוצה של $r$ מספרי צעדים):

$$\langle A \rangle_{\text{MPF}}(t) = \sum_{j=1}^r x_j \, \langle A \rangle_{k_j}(t),$$

כאשר $\langle A \rangle_{k_j}(t)$ הוא ערך התוחלת של תצפית $A$ בזמן $t$ המוערך מ-Circuit של Trotter עם $k_j$ צעדים, והמקדמים $\{x_j\}_{j=1}^r$ נבחרים כך שאיברי שגיאת Trotter המובילים בצירוף מתבטלים. נחזור לביטוי זה בשלב 4, שם נעריך אותו במפורש כדי לשלב את תוצאות ה-Trotter שלנו. הנקודה המעשית המרכזית היא של-Circuit העמוק ביותר ב-MPF דרוש רק $k_{\max}$ צעדים, שהוא הרבה יותר קטן מה-$k$ היחיד שהיה נדרש כדי להגיע לאותה שגיאת Trotter אפקטיבית באופן ישיר. ה-Circuit-ים הרדודים יותר הופכים את גישת ה-MPF למתאימה יותר לחומרה רועשת.

כיצד נקבעים המקדמים?

קיימות שתי משפחות של מקדמי MPF:

מקדמים סטטיים הם בלתי תלויים בהמילטוניאן, במצב ההתחלתי ובזמן ההתפתחות. הם נמצאים על ידי פתרון מערכת ליניארית $Ax = b$ שאוכפת ביטול של איברי שגיאת Trotter המובילים. עבור קבוצה של צעדי Trotter $\{k_j\}_{j=1}^r$ בשימוש עם נוסחת מוצר סימטרית מסדר $2\chi$, פיתוח שגיאת Trotter בחזקות הופכיות של $k_j$ מוביל למשוואות אילוץ בצורה:

$$\sum_{j=1}^r x_j = 1, \quad \sum_{j=1}^r \frac{x_j}{k_j^{\eta_n}} = 0 \quad (n = 0, \ldots, r-2),$$

כאשר המעריכים השלמים $\{\eta_n\}$ הם הסדרים של איברי שגיאת Trotter העוקבים עבור נוסחת המוצר הנבחרת. עבור PF סימטרית מסדר $2\chi$, השגיאה המובילה ב-$\left[S_{2\chi}(t/k)\right]^k$ מתקנתת כ-$1/k^{2\chi}$, עם תיקונים עוקבים ב-$1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ — כך שהמעריכים הם $\eta_n = 2\chi + 2n$. עבור PF-ים לא-סימטריות, גם חזקות אי-זוגיות וגם זוגיות תורמות ו-$\eta_n = 2\chi + n$. ראה מקור [1] לגזירה המלאה. המשוואה הראשונה במערכת לעיל מבטיחה חוסר הטיה (ה-MPF משחזר את ערך התוחלת המדויק בגבול $k_j \to \infty$), ו-$r-1$ המשוואות הנותרות מבטלות ברצף את $r-1$ איברי שגיאת ה-Trotter הראשונים. כאשר נורמת ה-$L_1$ המתקבלת $\|x\|_1$ גדולה מדי (מה שמגביר רעש דגימה), ניתן במקום זאת לפתור אופטימיזציה מקורבת שמגבילה את $\|x\|_1$ תוך מזעור $\|Ax - b\|$.

מקדמים דינמיים [2], [3] תלויים גם בהמילטוניאן, במצב ההתחלתי ובזמן ההתפתחות $t$. הם ממזערים את מרחק נורמת פרובניוס בין המצב המתפתח בזמן האמיתי לקירוב ה-MPF:

$$\|\rho(t) - \mu^D(t)\|_F^2 = 1 + \sum_{i,j} M_{ij}(t)\, x_i(t)\, x_j(t) - 2\sum_i L_i(t)\, x_i(t),$$

כאשר $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ היא מטריצת גראם של חפיפות בין מצבים שהתפתחו באמצעות Trotter עבור מספרי צעדים שונים $k_i, k_j$, ו-$L_i(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$ מודדת חפיפה עם המצב המדויק (המקורב). במדריך זה כמויות אלה מחושבות ביעילות באמצעות שיטות רשת טנזורית, ובפרט ה-backends מבוססי TeNPy ב-qiskit_addon_mpf.

מתי להשתמש ב-MPFs

MPFs מועילות ביותר כאשר:

• עומק ה-Circuit הוא צוואר הבקבוק. אם רעש החומרה מגביל את העומק שניתן להריץ, השתמש ב-MPFs כדי להשיג דיוק Trotter אפקטיבי גבוה יותר מ-Circuit-ים רדודים יותר.

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

• אתה משלב מספר צנוע של מספרי צעדי Trotter. בדרך כלל שילוב של $r = 3$–$5$ מספרי צעדים שונים $k_j$ מספיק כדי לבטל כמה איברי שגיאת Trotter מובילים תוך שמירה על $\|x\|_1$ בר-ניהול.

מתי MPFs עשויות לא לעזור

• זמני התפתחות קצרים מאוד. כאשר $t$ קטן מספיק כך שנוסחת Trotter מסדר נמוך יחיד כבר מדויקת, התקורה של הרצת Circuit-ים מרובים אינה נחוצה.

• משימות הכנת מצב. MPFs מייצרות ערך תוחלת מתוקן, לא מצב קוונטי מתוקן. אם אתה זקוק למצב שהתפתח בזמן האמיתי (למשל, כקלט לתת-שגרה קוונטית אחרת), MPFs אינן ישימות.

• מספרי צעדי Trotter המפרים את משטר ההתכנסות. הגזירה של המקדם הסטטי מפתחת כל $\left[S_{2\chi}(t/k_j)\right]^{k_j}$ בודד כטור ב-$t/k_j$; פיתוח זה מתכנס היטב רק כאשר $t/k_{\min} \lesssim 1$. אם $k_{\min}$ נבחר קטן מדי עבור ה-$t$ הנתון, ה-Circuit הרדוד ביותר נמצא הרחק מחוץ למשטר הפרטורבטיבי, איברי השגיאה מסדר גבוה יותר שה-MPF משאיר לא מבוטלים הופכים גדולים, והביטול עשוי לדרוש מקדמים גדולים. נורמת ה-$L_1$ $\|x\|_1$ היא האבחון המעשי: כאשר $\|x\|_1 \gg 1$, תקורת הדגימה $\propto \|x\|_1^2$ עלולה לגבור על הפחתת שגיאת Trotter. ראה המדריך לבחירת צעדי Trotter לפרטים.

מה מדריך זה מכסה

מדריך זה עובר על תהליך עבודה של MPF מקצה לקצה בשני שלבים. ראשית, דוגמה בקנה מידה קטן על סימולטור (שרשרת הייזנברג בת 10 Qubit-ים) מדגימה כיצד להגדיר את הבעיה, לחשב מקדמי MPF סטטיים ודינמיים, ולהשוות את ערכי התוחלת המתקבלים כנגד דיאגונליזציה מדויקת. לאחר מכן, דוגמה בקנה מידה גדול על חומרה (שרשרת XXZ בת 50 Qubit-ים) מראה כיצד לבצע טרנספילציה, להריץ על חומרת IBM Quantum עם הקלת שגיאות, ולעבד תוצאות לאחר מכן באמצעות מקדמי ה-MPF. לאורך הדרך, אנו משתמשים בחבילה qiskit_addon_mpf לצד כלי Qiskit סטנדרטיים.

דרישות

לפני שמתחילים עם מדריך זה, ודא שהדברים הבאים מותקנים:

• Qiskit SDK v2.0 ומעלה, עם תמיכה ב-visualization

• Qiskit Runtime v0.22 ומעלה (pip install qiskit-ibm-runtime)

• סימולטור Qiskit Aer (pip install qiskit-aer)

• תוסף MPF Qiskit עם ה-backend של TeNPy (pip install "qiskit-addon-mpf[tenpy]")

• כלי השירות (utilities) של תוספי Qiskit (pip install qiskit-addon-utils)

• SciPy (pip install scipy)

הגדרה

בהמשך אנו אוספים את כל ייבואי החבילות המשמשים לאורך מדריך זה בתא אחד. אנו גם מגדירים מעבר Transpiler בשם CollectAndCollapse שמאחד רוטציות rxx ו-ryy סמוכות ל-XXPlusYYGate יחיד. מעבר זה מיושם הן במהלך בניית ה-Circuit בשלב 1 (כדי לשמור על מספר Gate-ים נמוך) והן בעקיפין כאשר אנו מחלצים את המבנה השכבתי עבור ה-MPF הדינמי בשלב 4 (TeNPy מצפה ל-Gate-ים דו-Qubit, לא לזוגות רוטציות לא-מאוחדות).

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

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_addon_mpf": "qiskit-addon-mpf", "qiskit_addon_utils": "qiskit-addon-utils", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy", "tenpy": "physics-tenpy"}
_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
import matplotlib.pyplot as plt
from functools import partial
from copy import deepcopy

from qiskit import QuantumCircuit
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector
from qiskit.synthesis import SuzukiTrotter
from qiskit.transpiler import CouplingMap, PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)

from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService

from qiskit_addon_utils.problem_generators import (
    generate_xyz_hamiltonian,
    generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.static import setup_static_lse
from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import (
    setup_exact_problem,
    setup_sum_of_squares_problem,
    setup_frobenius_problem,
)
from qiskit_addon_mpf.backends.tenpy_layers import (
    LayerModel,
    LayerwiseEvolver,
)
from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state

from scipy.linalg import expm

# Suppress TeNPy's `unit_cell_width` future-API warning. The default
# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what
# `CouplingMap.from_line(...)` produces here, so the warning is informational.
warnings.filterwarnings(
    "ignore",
    message=r".*unit_cell_width.*",
    category=UserWarning,
)

# --- Helper: collect XX + YY rotations into a single gate ---
def filter_function(node):
    return node.op.name in {"rxx", "ryy"}

collect_function = partial(
    collect_using_filter_function,
    filter_function=filter_function,
    split_blocks=True,
    min_block_size=1,
)

def collapse_to_xx_plus_yy(block):
    param = 0.0
    for node in block.data:
        param += node.operation.params[0]
    return XXPlusYYGate(param)

collapse_function = partial(
    collapse_to_operation,
    collapse_function=collapse_to_xx_plus_yy,
)

pm = PassManager()
pm.append(CollectAndCollapse(collect_function, collapse_function))

דוגמה בקנה מידה קטן על סימולטור

שלב 1: מיפוי קלטים קלאסיים לבעיה קוונטית

אנו מתחילים עם מודל הייזנברג בן 10 Qubit-ים על קו, תוך שימוש במצב נל $\vert 0101\ldots01 \rangle$ כמצב ההתחלתי. ההמילטוניאן הוא:

$$\hat{\mathcal{H}}_{\text{Heis}} = J \sum_{i=1}^{L-1} \left(X_i X_{i+1} + Y_i Y_{i+1} + Z_i Z_{i+1}\right),$$

כאשר $J$ הוא עוצמת הצימוד בין שכנים קרובים. אנו מודדים את מתאם ה-ZZ $Z_{L/2-1} Z_{L/2}$ על זוג Qubit-ים באמצע השרשרת, ומשתמשים בצעדי Trotter $k_j = [1, 2, 4]$ עם נוסחת מוצר מסדר שני.

L = 10

# Generate coupling map and Hamiltonian
coupling_map = CouplingMap.from_line(L, bidirectional=False)

hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(1.0, 1.0, 1.0),
    ext_magnetic_field=(0.0, 0.0, 0.0),
)
print(hamiltonian)

SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

# Observable: ZZ on the middle pair of qubits
observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)

SparsePauliOp(['IIIIZZIIII'],
              coeffs=[1.+0.j])

# MPF parameters
mpf_trotter_steps = [1, 2, 4]
order = 2
symmetric = False

trotter_times = np.arange(0.5, 1.55, 0.1)
exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)

בניית Circuit-ים של Trotter

אנו יוצרים את ה-Circuit-ים המיישמים את התפתחויות הזמן המקורבות של Trotter עבור כל נקודת זמן וכל מספר צעדי Trotter. המעבר CollectAndCollapse המוגדר בסעיף ההגדרה אוסף רוטציות XX ו-YY ל-Gate-ים בודדים מסוג XX+YY, כדי להתכונן לסימולציית רשת טנזורית יעילה יותר בהמשך.

# Initial Neel state preparation
initial_state_circ = QuantumCircuit(L)
initial_state_circ.x([i for i in range(L) if i % 2 != 0])

all_circs = []
for total_time in trotter_times:
    mpf_trotter_circs = [
        generate_time_evolution_circuit(
            hamiltonian,
            time=total_time,
            synthesis=SuzukiTrotter(reps=num_steps, order=order),
        )
        for num_steps in mpf_trotter_steps
    ]

    mpf_trotter_circs = pm.run(
        mpf_trotter_circs
    )  # Collect XX and YY into XX + YY

    mpf_circuits = [
        initial_state_circ.compose(circuit) for circuit in mpf_trotter_circs
    ]
    all_circs.append(mpf_circuits)

mpf_circuits[-1].draw("mpl", fold=-1)

שלב 2: ייעול הבעיה לביצוע על חומרה קוונטית

עבור הדוגמה בקנה מידה קטן אנו מכוונים לסימולטור Aer. שני שינויים מתרחשים לפני שה-Circuit-ים מוכנים להרצה:

1. איסוף Gate-ים ברמת סימולציית ההמילטוניאן. בתא ההגדרה בנינו מעבר CollectAndCollapse שמאחד רוטציות rxx ו-ryy סמוכות ל-XXPlusYYGate יחיד. כבר יישמנו מעבר זה כאשר בנינו את ה-Circuit-ים של Trotter בשלב 1 (הקריאה ל-pm.run(...)). זה גם מפחית את מספר ה-Gate-ים הדו-Qubit-יים וגם מייצר מבנה שנוח יותר לסימולציית רשת טנזורית עבור חישוב המקדם הדינמי בהמשך.

2. הורדה ל-ISA של הסימולטור. בהמשך אנו מריצים את ה-pass manager המוגדר מראש של Qiskit ברמת optimization_level=3 כדי להוריד כל Circuit של Trotter לארכיטקטורת קבוצת ההוראות (ISA) של הסימולטור.

aer_sim = AerSimulator()
pm_sim = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)

isa_circs_all_times = [
    pm_sim.run([deepcopy(c) for c in mpf_circuits])
    for mpf_circuits in all_circs
]

שלב 3: ביצוע באמצעות פרימיטיבים של Qiskit

עבור הדוגמה בקנה מידה קטן אנו מריצים את ה-Circuit-ים של Trotter שהורדו ל-ISA דרך הפרימיטיב EstimatorV2 הנתמך על ידי Aer. עשיית זאת נותנת לנו ערך ייחוס נטול רעש עבור כל זוג $(k_j, t)$ — אלה הם ערכי $\langle A \rangle_{k_j}(t)$ שה-MPF ישלב בשלב 4. אנו סורקים על פני זמני התפתחות כך שנוכל בהמשך לצייר את עקומת סדרת הזמן המלאה של כל נוסחת מוצר בודדת ושל ה-MPF.

estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for isa_circuits in isa_circs_all_times:
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals_all_times.append([res.data.evs for res in result])
    mpf_stds_all_times.append([res.data.stds for res in result])

שלב 4: עיבוד לאחר והחזרת תוצאה בפורמט קלאסי רצוי

שלב 4 הוא המקום שבו ה-MPF בפועל נבנה. למרות שהמקדמים $x_j$ מחושבים כאן (ועבור הגרסה הדינמית, חישוב זה יכול להיות אינטנסיבי), באופן מושגי הם מתכון קלאסי לשילוב המדידות הקוונטיות משלב 3 לערך תוחלת מתוקן יחיד — לכן אנו מתייחסים לכל תהליך העבודה של המקדם והשילוב כעיבוד לאחר מכן.

כדי להעריך עד כמה ה-MPF עוקב אחר הדינמיקה האמיתית, אנו קודם מחשבים את ערכי התוחלת המדויקים שהתפתחו בזמן על ידי החזקה ישירה של ההמילטוניאן. זה ניתן ליישום רק בגלל ש-$L = 10$; בדוגמה בקנה מידה גדול על חומרה בהמשך נצטרך להסתמך על הערכות רשת טנזורית במקום זאת.

exact_expvals = []
for t in exact_evolution_times:
    exp_H = expm(-1j * t * hamiltonian.to_matrix())
    initial_state = Statevector(initial_state_circ).data
    time_evolved_state = exp_H @ initial_state

    exact_obs = (
        time_evolved_state.conj()
        @ observable.to_matrix()
        @ time_evolved_state
    ).real
    exact_expvals.append(exact_obs)

מקדמי MPF סטטיים

MPFs סטטיים משתמשים במקדמים $x_j$ שהם בלתי תלויים בזמן ההתפתחות, בהמילטוניאן ובמצב ההתחלתי. אנו מגדירים את המערכת הליניארית $Ax = b$ המתוארת ברקע ופותרים עבור המקדמים. המטריצה $A$ נקבעת על ידי מספרי צעדי Trotter $k_j$, הסדר $\chi$ של נוסחת המוצר, והאם הנוסחה סימטרית (מה ששולט במעריכים $\eta_n$).

עבור הדוגמה בקנה מידה קטן שלנו אנו משתמשים ב-$k_j = [1, 2, 4]$ עם נוסחת Suzuki-Trotter לא-סימטרית מסדר $2\chi=2$ (כך ש-$\chi=1$ ו-$\eta_n = 2 + n$, נותן $\eta_0 = 2,\, \eta_1 = 3$). המערכת הופכת ל:

$$A = \begin{bmatrix} 1 & 1 & 1\\ 1 & \frac{1}{2^2} & \frac{1}{4^2} \\ 1 & \frac{1}{2^3} & \frac{1}{4^3} \\ \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}.$$

השורה הראשונה אוכפת חוסר הטיה ($\sum_j x_j = 1$); השורה השנייה והשלישית מבטלות את איברי שגיאת ה-Trotter המובילים $1/k^2$ ו-$1/k^3$ מהסדר הבא, בהתאמה.

הגדר את ה-LSE

אנו משתמשים ב-setup_static_lse מ-qiskit_addon_mpf.static כדי להרכיב את המטריצה $A$ ואת הווקטור בצד ימין $b$ המתוארים לעיל. המטריצה $A$ תלויה לא רק ב-$k_j$ אלא גם בבחירת נוסחת המוצר שלנו — בפרט הסדר שלה $\chi$ והאם היא סימטרית. הדגל symmetric שולט בתבנית המעריך $\eta_n$ (נוסחאות סימטריות מייצרות רק איברי שגיאת Trotter בחזקות זוגיות; ראה מקור [1]). שים לב שכפי שמוצג במקור [2], הגדרת symmetric=True אינה הכרחית באופן מוחלט אפילו כאשר ה-PF הבסיסית סימטרית — ה-LSE הלא-סימטרי נשאר תקף (הוא אוכף אילוצים נוספים שאינם נחוצים).

עבור הדוגמה שלנו כבר הגדרנו order = 2 ו-symmetric = False בשלב 1.

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)

בדוק את המטריצה $A$ ואת הווקטור $b$ שנבנו כדי לוודא שהם תואמים למערכת שנכתבה לעיל.

lse.A

array([[1.      , 1.      , 1.      ],
       [1.      , 0.25    , 0.0625  ],
       [1.      , 0.125   , 0.015625]])

lse.b

array([1., 0., 0.])

לאחר שה-LSE בידינו, אנו פותרים עבור המקדמים הסטטיים $x_j$ באמצעות lse.solve() (זהו הפתרון הישיר $x = A^{-1}b$).

mpf_coeffs = lse.solve()
print(
    f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)

The static coefficients associated with the ansatze are: [ 0.04761905 -0.57142857  1.52380952]

ייעל עבור $x$ באמצעות מודל מדויק

כחלופה לחישוב $x = A^{-1}b$, אתה יכול להשתמש ב-setup_exact_model כדי לבנות מופע cvxpy.Problem שמשתמש ב-LSE כאילוצים ושהפתרון האופטימלי שלו ייתן $x$.

model_exact, coeffs_exact = setup_exact_problem(lse)
model_exact.solve()
print(coeffs_exact.value)

[ 0.04761905 -0.57142857  1.52380952]

print(
    "L1 norm of the exact coefficients:",
    np.linalg.norm(coeffs_exact.value, ord=1),
)

L1 norm of the exact coefficients: 2.1428571428556378

ייעל עבור $x$ באמצעות מודל מקורב

ייתכן שנורמת ה-$L_1$ עבור קבוצת ערכי $k_j$ שנבחרה תיחשב גבוהה מדי. אם זה המקרה ואתה לא יכול לבחור קבוצה שונה של ערכי $k_j$, אתה יכול להשתמש בפתרון מקורב שמגביל את נורמת ה-$L_1$ לסף נבחר תוך מזעור $\|Ax - b\|$. עיין במדריך כיצד להשתמש במודל המקורב.

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=1.5
)
model_approx.solve()
print(coeffs_approx.value)
print(
    "L1 norm of the approximate coefficients:",
    np.linalg.norm(coeffs_approx.value, ord=1),
)

[-1.10294118e-03 -2.48897059e-01  1.25000000e+00]
L1 norm of the approximate coefficients: 1.5

מקדמי MPF דינמיים

ה-MPF הסטטי מבטל איברי שגיאת Trotter באופן בלתי תלוי בהמילטוניאן ובמצב, כך שהוא לא בהכרח מייצר את שגיאת הקירוב הקטנה ביותר האפשרית עבור המילטוניאן ומצב התחלתי נתונים. ה-MPF הדינמי (מקורות [2], [3]) מוצא במקום זאת מקדמים תלויי-זמן $x_i(t)$ שממזערים את מרחק נורמת פרובניוס $\|\rho(t) - \mu^D(t)\|_F^2$ בכל זמן $t$. כפי שהוצג ברקע, זה דורש את מטריצת החפיפה $M_{ij}(t)$ בין מצבים שהתפתחו באמצעות Trotter ואת החפיפה $L_i(t)$ עם המצב המדויק — שניהם אנו מעריכים באמצעות backends של רשת טנזורית (TeNPy) ב-qiskit_addon_mpf.

כדי להגדיר את ה-LSE הדינמי אנו זקוקים לשלושה מרכיבים:

1. פונקציית factory ליוצר ההתפתחות (evolver) המקורב שהתוסף יריץ עבור כל $k_j$ כדי לייצר $\rho_{k_j}(t)$ כ-MPS/MPO. אנו בונים אותו מהמבנה השכבתי של ה-Circuit של Trotter מסדר $2$ (שכבה אחת עבור כל slice_by_depth), עטוף כ-LayerwiseEvolver עם פרמטרי קיצוץ של TeNPy.

2. מפעל מפתח מדויק שמייצר ייחוס $\rho(t)$ בדיוק גבוה. אנו משתמשים ב-Circuit Suzuki-Trotter מסדר רביעי עם צעד זמן קטן (dt=0.1, order=4) כפרוקסי להתפתחות מדויקת.

3. מפעל זהות ו-MPS מצב התחלתי שמזרעים את סימולציית TeNPy.

התא שלהלן בונה את מפעל המפתח המקורב.

# Create approximate time-evolution circuits
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)  # collect XX and YY

# Find layers in the circuit
layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)

# Create tensor network models
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

# Create the time-evolution object
approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

אזהרה

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

אנו מקרבים את המצב המדויק המתפתח בזמן באמצעות נוסחת Suzuki-Trotter מסדר רביעי, תוך שימוש בצעד זמן קטן dt=0.1. פרמטרי הקיטום (truncation) של TeNPy יכולים להשפיע על הדיוק, ולכן חשוב לבחון טווח של ערכים.

single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

לבסוף, אנו מגדירים identity_factory שמניב את מצב ה-MPO ההתחלתי ומכינים את מצב נל ההתחלתי כ-MPS שתואם לסריג המשמש את מודל Trotter השכבתי.

def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)

mps_initial_state = MPS_neel_state(models[0].lat)

עם המפעלים במקום אנו כעת מחשבים את המקדמים הדינמיים בכל זמן התפתחות. עבור כל $t$, setup_dynamic_lse בונה את מטריצות החפיפה הרלוונטיות באמצעות TeNPy, ו-setup_frobenius_problem מחזיר cvxpy.Problem שממזער את עלות נורמת פרובניוס. הפותר מחזיר מקדמים $x_j(t)$ מותאמים לאותו זמן; אנו אוספים אותם ב-mpf_dynamic_coeffs_list. אם הפותר נכשל עבור $t$ נתון, אנו חוזרים למקדמי אפס כך שהלולאה ממשיכה.

mpf_dynamic_coeffs_list = []
for t in trotter_times:
    print(f"Computing dynamic coefficients for time={t}")
    lse = setup_dynamic_lse(
        mpf_trotter_steps,
        t,
        identity_factory,
        exact_factory,
        approx_factory,
        mps_initial_state,
    )
    problem, coeffs = setup_frobenius_problem(lse)
    try:
        problem.solve()
        mpf_dynamic_coeffs_list.append(coeffs.value)
    except Exception as error:
        mpf_dynamic_coeffs_list.append(np.zeros(len(mpf_trotter_steps)))
        print(error, "Calculation Failed for time", t)
    print("")

Computing dynamic coefficients for time=0.5

Computing dynamic coefficients for time=0.6

Computing dynamic coefficients for time=0.7

Computing dynamic coefficients for time=0.7999999999999999

Computing dynamic coefficients for time=0.8999999999999999

Computing dynamic coefficients for time=0.9999999999999999

Computing dynamic coefficients for time=1.0999999999999999

Computing dynamic coefficients for time=1.1999999999999997

Computing dynamic coefficients for time=1.2999999999999998

Computing dynamic coefficients for time=1.4

Computing dynamic coefficients for time=1.4999999999999998

שילוב ערכי תוחלת של Trotter עם מקדמי ה-MPF

כעת אנו מעריכים $\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)$ עבור כל קבוצת מקדמים (סטטי-מדויק, סטטי-מקורב, ודינמי), מפיצים את השגיאות התקניות לכל Circuit, ומציירים את סדרת הזמן המתקבלת כנגד עקומת הדיאגונליזציה המדויקת.

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    trotter_curve, trotter_curve_error = [], []
    for trotter_expvals, trotter_stds in zip(
        mpf_expvals_all_times, mpf_stds_all_times
    ):
        trotter_curve.append(trotter_expvals[k])
        trotter_curve_error.append(trotter_stds[k])

    plt.errorbar(
        trotter_times,
        trotter_curve,
        yerr=trotter_curve_error,
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

# Get expectation values at all times for the static MPF with exact coeffs
exact_mpf_curve, exact_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_exact.value, trotter_stds)
            ]
        )
    )
    exact_mpf_curve_error.append(mpf_std)
    exact_mpf_curve.append(trotter_expvals @ coeffs_exact.value)

plt.errorbar(
    trotter_times,
    exact_mpf_curve,
    yerr=exact_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Exact",
    color="purple",
)

# Get expectation values at all times for the static MPF with approximate coeffs
approx_mpf_curve, approx_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_approx.value, trotter_stds)
            ]
        )
    )
    approx_mpf_curve_error.append(mpf_std)
    approx_mpf_curve.append(trotter_expvals @ coeffs_approx.value)

plt.errorbar(
    trotter_times,
    approx_mpf_curve,
    yerr=approx_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Approx",
    color="orange",
)

# Get expectation values at all times for the dynamic MPF
dynamic_mpf_curve, dynamic_mpf_curve_error = [], []
for trotter_expvals, trotter_stds, dynamic_coeffs in zip(
    mpf_expvals_all_times, mpf_stds_all_times, mpf_dynamic_coeffs_list
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(dynamic_coeffs, trotter_stds)
            ]
        )
    )
    dynamic_mpf_curve_error.append(mpf_std)
    dynamic_mpf_curve.append(trotter_expvals @ dynamic_coeffs)

plt.errorbar(
    trotter_times,
    dynamic_mpf_curve,
    yerr=dynamic_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Dynamic MPF",
    color="pink",
)

# Exact expectation values
plt.plot(
    exact_evolution_times,
    exact_expvals,
    color="red",
    linestyle="--",
    label="Exact time-evolution",
)

plt.title(f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ vs time")
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

התרשים לעיל ממחיש את יחסי הגומלין בין שגיאת Trotter לשגיאת דגימה.

• שגיאת Trotter. נוסחאות המוצר הבודדות (סמנים אפורים) סוטות מהעקומה המדויקת יותר ויותר ככל שהזמן גדל. ל-Circuit $k=1$ יש את הסטייה הגדולה ביותר והוא הרדוד ביותר, אך הוא גם כבר במשטר שבו $t/k \gtrsim 1$, כך שאיבר השגיאה המוביל $1/k^{2}$ גדול. שילובי ה-MPF (סמנים צבעוניים) מבטלים כמה מאיברי שגיאת ה-Trotter המובילים הללו, כך שהם עוקבים אחר העקומה המדויקת הרבה יותר מקרוב מכל Circuit $k_j$ בודד. הפער הנותר משקף את איברי ה-Trotter מסדר גבוה יותר שה-MPF אינו מבטל: MPF סטטי מסדר $2$, $r=3$ הורג רק את שני סדרי השגיאה הראשונים, ובערכי $t/k_{\min}$ גדולים הזנב הלא-מבוטל שולט בסופו של דבר — כך שה-MPF לא מבטיח ש-Circuit-ים רדודים מאוד יישארו מדויקים בזמנים שרירותיים.

• שגיאת דגימה. מוטות השגיאה הרחבים יותר בעקומות ה-MPF הם תוצאה ישירה של הצירוף הליניארי: הפצת שגיאות תקניות עצמאיות לכל Circuit $\sigma_{k_j}$ נותנת שונות כוללת $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$. לכן, ככל ש-$\|x\|_2$ גדול יותר (ובפועל $\|x\|_1$, שזה מה שאנו שולטים בו), כך יידרשו יותר shots כדי להגיע לאי-ודאות יעד נתונה. זהו הפשרה שמאחורי אפשרות הפותר המקורב ברקע: אנו מגבילים את $\|x\|_1$ כדי לשמור על תקורה זו בת-ניהול. באופן מכריע, בניגוד לשגיאת Trotter, שגיאת דגימה מצטמצמת עם $1/\sqrt{N_{\text{shots}}}$, כך שתמיד ניתן להפחית אותה על ידי הוצאת יותר shots.

בדוגמה בקנה מידה גדול על חומרה בהמשך, רעש חומרה נכנס כמקור שגיאה נוסף בכל $\langle A \rangle_{k_j}$, אשר באופן דומה מוגבר על ידי מקדמי ה-MPF. נראה כיצד הקלת שגיאות פועלת בגומלין עם MPFs באותו סעיף.

דוגמה בקנה מידה גדול על חומרה

בסעיף זה אנו מגדילים את הבעיה מעבר למה שניתן לסמלץ באופן מדויק. אנו משחזרים חלק מהתוצאות המוצגות במקור [3], תוך שימוש בשרשרת XXZ בת 50 Qubit-ים בזמן $t = 3$. אנו עוקבים אחר אותו תהליך עבודה בן ארבעה שלבים כמו הדוגמה בקנה מידה קטן, כעת מכוונים לחומרה קוונטית אמיתית עם הקלת שגיאות. כמו בתבנית, כל שלב מסומן בתוך הקוד, וצעד יחיד יכול להשתרע על פני מספר תאים כאשר התוצרים הביניים שלו שווים בדיקה. המיפוי משקף את הדוגמה בקנה מידה קטן: הגדרת המילטוניאן, בחירת פרמטרי Trotter, חישוב מקדמי MPF (סטטיים ודינמיים), ובניית Circuit-ים. ההבדלים המרכזיים הם:

• המילטוניאן XXZ על 50 אתרים עם צימודים אקראיים שנשלפו מ-$\mathcal{U}(0.5, 1.5)$ (מקור [3]).

• נוסחת Trotter סימטרית מסדר שני עם $k_j = [3, 4, 6]$ (כך ש-$\chi=1$, symmetric=True).

• זמן התפתחות קבוע יחיד $t = 3$. עם $k_{\min}=3$ זה נותן $t/k_{\min}=1$, ושומר על המרכיבים הרדודים בתוך משטר ההתכנסות של Trotter שבו מודל השגיאה המובילה שה-MPF מסתמך עליו תקף.

• הרצת השוואה נוספת של Circuit יחיד עם $k = 10$ צעדי Trotter, המשמשת כבסיס. בחרנו $k = 10$ כיוון שהעומק הדו-Qubit-י שלו על חומרה עמוק יותר מהמרכיב העמוק ביותר של MPF ($k_{\max}=6$) בתוספת התקורה של הרצת מספר Circuit-ים של MPF — עמוק מספיק כדי להיות מוגבל-רעש, שזה המשטר שבו שילוב ה-MPF צפוי להתעלות על הבסיס של Circuit יחיד. זוהי השוואת "Circuit עמוק יחיד" מול שילוב ה-MPF, לא Circuit שמכוון לשגיאת Trotter האפקטיבית של ה-MPF (שיידרש הרבה יותר צעדים).

שים לב שלמרות שאנו עדיין בשלב 1 כאן (מיפוי ובניית Circuit), אנו גם מחשבים מראש את המקדמים הדינמיים לצד הסטטיים בתא זה. מקדמים דינמיים תלויים ב-$H$ וב-$t$ אך לא במדידות הקוונטיות, כך שניתן לחשב אותם בכל עת לפני שלב 4. אנו עושים זאת כעת כדי לשמור על כל הגדרת ה-MPF במקום אחד.

# -------------------------Step 1-------------------------
L = 50
coupling_map = CouplingMap.from_line(L, bidirectional=False)

# XXZ Hamiltonian with random couplings (Ref. [3])
np.random.seed(0)
even_edges = list(coupling_map.get_edges())[::2]
odd_edges = list(coupling_map.get_edges())[1::2]

Js = np.random.uniform(0.5, 1.5, size=L)
hamiltonian = SparsePauliOp(Pauli("I" * L))
for i, edge in enumerate(even_edges + odd_edges):
    hamiltonian += SparsePauliOp.from_sparse_list(
        [
            ("XX", (edge), 2 * Js[i]),
            ("YY", (edge), 2 * Js[i]),
            ("ZZ", (edge), 4 * Js[i]),
        ],
        num_qubits=L,
    )

observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)

total_time = 3
mpf_trotter_steps = [3, 4, 6]
order = 2
symmetric = True

# Static coefficients
lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(f"Static coefficients: {mpf_coeffs}")
print(f"L1 norm: {np.linalg.norm(mpf_coeffs, ord=1)}")

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=2.0
)
model_approx.solve()
print(f"Approximate coefficients: {coeffs_approx.value}")
print(f"L1 norm (approx): {np.linalg.norm(coeffs_approx.value, ord=1)}")

# -------------------------Dynamic coefficients-------------------------
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)

layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {"chi_max": 64, "svd_min": 1e-8, "trunc_cut": None},
        "max_delta_t": 4,
    },
)

single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {"chi_max": 64, "svd_min": 1e-8, "trunc_cut": None},
        "max_delta_t": 3,
    },
)

def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)

mps_initial_state = MPS_neel_state(models[0].lat)

print(f"Computing dynamic coefficients for time={total_time}")
lse_dyn = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs_dyn = setup_frobenius_problem(lse_dyn)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs_dyn.value
except Exception as error:
    mpf_dynamic_coeffs = np.zeros(len(mpf_trotter_steps))
    print(error, "Calculation Failed")

# -------------------------Step 1 (cont): Build circuits-------------------------
mpf_circuits = []
for k in mpf_trotter_steps:
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2])
    trotter_circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(reps=k, order=order),
        time=total_time,
    )
    circuit.compose(trotter_circ, qubits=range(L), inplace=True)
    mpf_circuits.append(circuit)

# Baseline "single deep circuit" comparison run with k=10 Trotter steps.
# Its two-qubit depth is deeper than the deepest MPF constituent (k_max=6) plus
# the overhead of running multiple circuits, pushing it into the noise-limited
# regime where MPF is expected to outperform. It does NOT target the MPF's effective
# Trotter error (which would require many more steps).
comp_circuit = QuantumCircuit(L)
comp_circuit.x([i for i in range(L) if i % 2])
trotter_circ = generate_time_evolution_circuit(
    hamiltonian,
    synthesis=SuzukiTrotter(reps=10, order=order),
    time=total_time,
)
comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)
mpf_circuits.append(comp_circuit)

Static coefficients: [ 0.42857143 -1.82857143  2.4       ]
L1 norm: 4.65714285714286
Approximate coefficients: [-0.4942491   0.40206845  1.09218065]
L1 norm (approx): 1.9884981979026675
Computing dynamic coefficients for time=3

כעת אנו מבצעים אופטימיזציה של ה-Circuit-ים עבור ה-Backend הנבחר. אנו משתמשים ב-pass manager המוגדר מראש של Qiskit ברמת optimization_level=3, אשר בוחר באופן אוטומטי קבוצה טובה של Qubit-ים פיזיים ומנתב כל Circuit לטופולוגיית ההתקן.

# -------------------------Step 2-------------------------
service = QiskitRuntimeService()
# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=L)
backend = service.backend("ibm_fez")
print(backend)

transpiler = generate_preset_pass_manager(
    optimization_level=3, backend=backend
)
transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

<IBMBackend('ibm_fez')>

הרצת Circuit-ים עמוקים יותר על חומרה אמיתית דורשת הקלת שגיאות אגרסיבית. אנו מפעילים ניתוק דינמי, טוויירלינג של Gate ומדידה, הקלת שגיאות מדידה, ואקסטרפולציית רעש אפס (ZNE). שים לב שגורמי הרעש של ZNE שבהם אנו משתמשים כאן (1, 1.2, 1.4) קטנים יותר מאשר בתרחיש Circuit רדוד, מכיוון שהמרכיבים העמוקים יותר של ה-MPF כבר קרובים לסף הרעש והגברות רעש גדולות היו דוחפות אותם מעבר לנקודה שבה אקסטרפולציית ZNE אמינה.

אנו שולחים את כל ארבעת ה-Circuit-ים (שלושה מרכיבי MPF ב-$k_j = [3, 4, 6]$ בתוספת הבסיס $k = 10$) בעבודת Estimator יחידה.

# -------------------------Step 3-------------------------
estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Error suppression/mitigation
estimator.options.dynamical_decoupling.enable = True
estimator.options.twirling.enable_gates = True
estimator.options.twirling.enable_measure = True
estimator.options.twirling.num_randomizations = "auto"
estimator.options.twirling.strategy = "active-accum"
estimator.options.resilience.measure_mitigation = True
estimator.options.experimental.execution_path = "gen3-turbo"

estimator.options.resilience.zne_mitigation = True
estimator.options.resilience.zne.noise_factors = (1, 1.2, 1.4)
estimator.options.resilience.zne.extrapolator = "linear"

estimator.options.environment.job_tags = ["TUT_MPF"]

job_50 = estimator.run(
    [
        (circ, observable)
        for circ, observable in zip(transpiled_circuits, isa_observables)
    ]
)

אנו שולפים את ערכי התוחלת לכל Circuit ואת סטיות התקן מתוצאת העבודה, ולאחר מכן משלבים אותם עם כל קבוצת מקדמי MPF בדיוק כמו בדוגמה בקנה מידה קטן: $\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}$, עם שונות מופצת $\sigma^2 = \sum_j x_j^2 \sigma_{k_j}^2$.

# -------------------------Step 4-------------------------
result = job_50.result()
evs = [res.data.evs for res in result]
std = [res.data.stds for res in result]

print(evs)
print(std)

[array(-0.07916195), array(-0.04479681), array(-0.2560756), array(-0.06045848)]
[array(0.04605538), array(0.10056336), array(0.14426151), array(0.04059092)]

exact_mpf_std = np.sqrt(
    sum([(coeff**2) * (std**2) for coeff, std in zip(mpf_coeffs, std[:3])])
)
print(
    "Exact static MPF expectation value: ",
    evs[:3] @ mpf_coeffs,
    "+-",
    exact_mpf_std,
)
approx_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_approx.value, std[:3])
        ]
    )
)
print(
    "Approximate static MPF expectation value: ",
    evs[:3] @ coeffs_approx.value,
    "+-",
    approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(mpf_dynamic_coeffs, std[:3])
        ]
    )
)
print(
    "Dynamic MPF expectation value: ",
    evs[:3] @ mpf_dynamic_coeffs,
    "+-",
    dynamic_mpf_std,
)

Exact static MPF expectation value:  -0.5665938395816946 +- 0.3925273058119915
Approximate static MPF expectation value:  -0.25856647611537903 +- 0.164249927266166
Dynamic MPF expectation value:  -0.12667812062949296 +- 0.06059471006973169

sym = {3: "^", 4: "s", 6: "p"}
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs[k],
        yerr=std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

plt.errorbar(
    3,
    evs[-1],
    yerr=std[-1],
    alpha=0.5,
    markersize=8,
    marker="x",
    color="blue",
    label="10 Trotter steps",
)

plt.errorbar(
    4,
    evs[:3] @ mpf_coeffs,
    yerr=exact_mpf_std,
    markersize=4,
    marker="o",
    color="purple",
    label="Static MPF",
)

plt.errorbar(
    5,
    evs[:3] @ coeffs_approx.value,
    yerr=approx_mpf_std,
    markersize=4,
    marker="o",
    color="orange",
    label="Approximate static MPF",
)

plt.errorbar(
    6,
    evs[:3] @ mpf_dynamic_coeffs,
    yerr=dynamic_mpf_std,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

exact_obs = -0.24384471447172074  # Calculated via Tensor Network calculation
plt.axhline(
    y=exact_obs, linestyle="--", color="red", label="Exact time-evolution"
)

plt.title(
    f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ at time {total_time} for the different methods"
)
plt.xlabel("Method")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

כמה תצפיות לגבי תוצאות החומרה לעיל:

• להעמיק אינו חינם על חומרה. הבסיסים של Circuit יחיד מספרים את הסיפור באופן ישיר: ה-Circuit $k = 6$ הוא כמעט מדויק לחלוטין ($-0.256$ לעומת הייחוס $-0.244$), ובכל זאת הבסיס העמוק יותר $k = 10$ הוא גרוע יותר ($-0.061$, סטייה של כ-$0.18$), לא טוב יותר. ברגע ששגיאת Trotter כבר קטנה, הוספת צעדים בעיקר מעמיקה את ה-Circuit וצוברת יותר רעש Gate ודקוהרנטיות. זהו בדיוק המשטר שעבורו MPFs נבנו: להגיע לדיוק של Circuit עמוק תוך שימוש רק במרכיבים רדודים.

• MPF בעל נורמה קטנה מנצח את ה-Circuit העמוק היחיד. ה-MPF הסטטי-המקורב (מוגבל ל-$\|x\|_1 \approx 2$) נוחת על $-0.259$, בטווח של כ-$0.015$ מהייחוס והרבה יותר קרוב מהבסיס $k = 10$. ה-MPF הדינמי ($-0.127$) גם הוא מנצח בנוחות את הבסיס ההוא. שניהם משלבים רק את ה-Circuit-ים הרדודים $k_j = [3, 4, 6]$, ובכל זאת משחזרים תשובה שה-Circuit העמוק היחיד לא הצליח להשיג.

• נורמת המקדם חשובה יותר מאופטימליות מתמטית. ל-MPF הסטטי-המדויק יש $\|x\|_1 = 4.66$ והוא הגרוע ביותר מבין כל האומדנים ($-0.567$, סטייה של יותר מ-$0.3$): נורמת המקדם הגדולה מגבירה את רעש ה-Gate השיורי, הדקוהרנטיות, ושגיאת ה-ZNE על כל $\langle A \rangle_{k_j}$ בערך באותו גורם, ומכריעה על ביטול שגיאת ה-Trotter שהיא קונה. הגבלת הנורמה (הפותר הסטטי-המקורב, $\|x\|_1 \approx 2$) מסירה את ההכרעה הזו ונותנת את ההערכה הטובה ביותר — אף על פי שהמקדמים שלה כבר לא מבטלים את שגיאת ה-Trotter המובילה באופן מדויק.

• Circuit-ים רדודים בודדים עדיין יכולים להיות תחרותיים. המרכיב הבודד $k = 6$ ($-0.256$) הוא בעצמו כמעט מדויק לחלוטין כאן — בהרצה זו הוא אפילו קרוב יותר במעט מה-MPF הסטטי-המקורב. הבעיה היא שאתה לא יודע מראש איזה $k$ בודד נמצא בנקודה המתוקה של "מתכנס אך לא עדיין מוגבל-רעש", והבחירה הנראית בטוחה של פשוט להעמיק יותר ($k = 10$) כדי להבטיח התכנסות Trotter היא בדיוק זו שנכשלת. ה-MPF נותן שילוב מבוסס-עקרונות של Circuit-ים רדודים שלא דורש ניחוש העומק הנכון.

המסקנה המעשית היא שעל חומרה, MPFs צריכים להיות מזווגים עם הקלת שגיאות חזקה על כל $\langle A \rangle_{k_j}$ בודד, נורמת ה-$L_1$ של המקדם צריכה להישאר צנועה (השתמש בפותר המקורב, או ב-MPF הדינמי), וצעדי Trotter $k_j$ צריכים להיבחר כך ש-$t/k_{\min} \lesssim 1$ — כאן $k_{\min} = 3$ ב-$t = 3$ נותן $t/k_{\min} = 1$, ושומר על המרכיבים בתוך המשטר המתכנס שבו מודל השגיאה המובילה שה-MPF הסטטי מסתמך עליו תקף. עם בחירות אלה, ה-MPFs בעלי הנורמה הקטנה כאן תואמים Circuit יחיד מתכנס בעוד שהבסיס הנאיבי "פשוט להעמיק יותר" לא, ומשחזר את יתרון העומק-מול-דיוק המוצג במקור [3]. שים לב גם שהרצות בודדות רועשות — בשליחה שונה של אותה עבודה (או Backend שונה), הסדר המדויק יכול להשתנות; המגמות החזקות הן ש-MPFs בעלי $\|x\|_1$ קטן מתפקדים היטב, ה-MPF הסטטי-המדויק בעל $\|x\|_1$ גדול מוגבר על ידי רעש חומרה, וה-Circuit היחיד העמוק מדי מוגבל-רעש.

צעדים הבאים

המלצות

אם עבודה זו עניינה אותך, ייתכן שתתעניין בחומר הבא:

• כיצד לבחור את צעדי Trotter עבור MPF — הדרכה מעשית לבחירת ערכי $k_j$ כדי להימנע מאי-יציבות
• כיצד להשתמש במודל המקורב — כוונון אילוץ נורמת ה-$L_1$ ואפשרויות הפותר עבור ה-MPF הסטטי המקורב
• תיעוד API של qiskit-addon-mpf — תיעוד מלא עבור מודולי סטטי, דינמי, ו-backend

אסמכתאות

[1] Vazquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation. Quantum, 7, 1067 (2023)

[2] Zhuk, S., Robertson, N. F., & Bravyi, S. Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation. Physical Review Research, 6(3), 033309 (2024)

[3] Robertson, N. F., et al. Tensor network enhanced dynamic multiproduct formulas. arXiv:2407.17405 (2024)