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

במחברת זו תלמד/י כיצד להשתמש ב-נוסחת רב-מכפלה (MPF) כדי להשיג שגיאת טרוטר נמוכה יותר באובזרבל שלנו בהשוואה לשגיאה שמייצר המעגל העמוק ביותר של טרוטר שנריץ בפועל. תעשה/י זאת על ידי עבודה לפי השלבים של תבנית Qiskit:

• שלב 1: מיפוי לבעיה קוונטית
  • אתחול ההמילטוניאן של הבעיה שלנו
  • שימוש ב-MPF ליצירת מעגלי האבולוציה בזמן בשיטת טרוטר
• שלב 2: אופטימיזציה של הבעיה
  • כאן נבצע Transpile למעגלים שלנו עבור GenericBackendV2
• שלב 3: ביצוע ניסויים
  • שימוש ב-StatevectorEstimator לשם הפשטות במחברת זו
• שלב 4: שחזור תוצאות
  • חישוב ערך הציפייה של ה-MPF

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

1א: הגדרת ההמילטוניאן שלנו

אנחנו משתמשים במודל איזינג על שרשרת של 10 אתרים:

$$\hat{\mathcal{H}}_{\text{Ising}} = \sum_{i=1}^{9} J_{i,(i+1)} Z_i Z_{(i+1)} + \sum_{i=1}^{10} h_i X_i \, ,$$

כאשר $J$ הוא עוצמת הצימוד בין שני אתרים ו-$h$ הוא השדה המגנטי החיצוני. החבילה qiskit_addon_utils מספקת פונקציות שימושיות לשימוש חוזר למטרות שונות.

המודול qiskit_addon_utils.problem_generators מספק פונקציות ליצירת המילטוניאנים מסוג היזנברג על גרף קישוריות נתון. גרף זה יכול להיות rustworkx.PyGraph או CouplingMap, מה שמקל על השימוש בו בתהליכי עבודה ממוקדי Qiskit.

בהמשך, ניצור שרשרת פשוטה של 10 Qubit באמצעות המתודה CouplingMap.from_line.

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-addon-mpf qiskit-addon-utils rustworkx scipy

from qiskit.transpiler import CouplingMap

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(10, bidirectional=False)

from rustworkx.visualization import graphviz_draw

graphviz_draw(coupling_map.graph, method="circo")

לאחר מכן, אנחנו יוצרים את SparsePauliOp על הקישוריות הנתונה עם הקבועים הרצויים.

from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

# Get a qubit operator describing the Ising field model
hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(0.0, 0.0, 1.0),
    ext_magnetic_field=(0.4, 0.0, 0.0),
)
print(hamiltonian)

SparsePauliOp(['IIIIIIIZZI', 'IIIIIZZIII', 'IIIZZIIIII', 'IZZIIIIIII', 'IIIIIIIIZZ', 'IIIIIIZZII', 'IIIIZZIIII', 'IIZZIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIXI', 'IIIIIIIXII', 'IIIIIIXIII', 'IIIIIXIIII', 'IIIIXIIIII', 'IIIXIIIIII', 'IIXIIIIIII', 'IXIIIIIIII', 'XIIIIIIIII'],
              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, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j,
 0.4+0.j, 0.4+0.j, 0.4+0.j])

האובזרבל שנמדוד הוא המגנטיזציה הכוללת, אותה ניתן לבנות בפשטות כמוצג להלן:

from qiskit.quantum_info import SparsePauliOp

L = coupling_map.size()
observable = SparsePauliOp.from_sparse_list([("Z", [i], 1 / L / 2) for i in range(L)], num_qubits=L)
print(observable)

SparsePauliOp(['IIIIIIIIIZ', 'IIIIIIIIZI', 'IIIIIIIZII', 'IIIIIIZIII', 'IIIIIZIIII', 'IIIIZIIIII', 'IIIZIIIIII', 'IIZIIIIIII', 'IZIIIIIIII', 'ZIIIIIIIII'],
              coeffs=[0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j,
 0.05+0.j, 0.05+0.j, 0.05+0.j])

1ב: נוסחאות רב-מכפלה

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

כדי להמחיש זאת, נגדיר MPF כך:

$$\mu(t) = \sum_{j} x_j \rho^{k_j}_{j}\left(\frac{t}{k_j}\right) + \text{some remaining Trotter error} \, ,$$

כאשר $x_j$ הם מקדמי המשקל שלנו, $\rho^{k_j}_j$ היא מטריצת הצפיפות המתאימה למצב הטהור המתקבל על ידי אבולוציה של המצב ההתחלתי עם נוסחת המכפלה $S^{k_j}$ הכוללת $k_j$ צעדי טרוטר, ו-$j$ הוא אינדקס מספר נוסחאות המכפלה המרכיבות את ה-MPF.

הנקודה המרכזית היא ששגיאת הטרוטר הנותרת קטנה יותר מהשגיאה שהיינו מקבלים באמצעות שימוש פשוט בערך $k_j$ הגדול ביותר!

ניתן להסתכל על התועלת של MPF משתי זוויות:

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

מבוא ל-MPF סטטיות

MPF סטטיות הן אלה שבהן ערכי $x_j$ אינם תלויים בזמן האבולוציה, $t$.

קביעת מקדמי ה-MPF הסטטית עבור קבוצה נתונה של ערכי $k_j$ שקולה לפתרון מערכת משוואות לינארית: $Ax=b$, כאשר $x$ הם המקדמים הנדרשים, $A$ היא מטריצה התלויה ב-$k_j$ ובסוג ה-PF שבו אנחנו משתמשים ($S$), ו-$b$ הוא וקטור אילוצים. לשם קיצור, לא נרחיב כאן ובמקום זאת נפנה אותך לתיעוד של LSE.

ניתן למצוא פתרון ל-$x$ אנליטית כ-$x = A^{-1}b$, ראה למשל Carrera Vazquez et al., 2023 או Zhuk et al., 2023. אולם פתרון מדויק זה עלול להיות "בעייתי מבחינת קונדיציה" ולהוביל לנורמות L1 גדולות מאוד של המקדמים $x$, מה שעלול לגרום לביצועים גרועים של ה-MPF. חלופה לכך היא קבלת פתרון מקורב המינימיזציה את נורמת ה-L1 של $x$ במטרה לייעל את התנהגות ה-MPF.

בהמשך תלמד/י כיצד לעשות את כל זאת.

בחירת $k_j$

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

כאן פשוט נבחר כמה ערכים קבועים ל-$k_j$. הערך הקטן ביותר מונחה על ידי זמן האבולוציה המטרה $t=8.0$, שבדרך כלל מחייב לעמוד בתנאי $t/k_{\text{min}} \lt 1$, אך אמפירית אנחנו יודעים שהגדרתו שווה ל-$1$ בדרך כלל עובדת אף היא. אם ברצונך ללמוד עוד על כך וכיצד לבחור את שאר ערכי $k_j$, ראה במדריך הרלוונטי: How to choose the Trotter steps for an MPF.

time = 8.0
trotter_steps = (8, 12, 19)

הגדרת ה-LSE

כעת שבחרנו את ערכי $k_j$, עלינו תחילה לבנות את ה-LSE, $Ax=b$ כפי שהוסבר לעיל. המטריצה $A$ תלויה לא רק ב-$k_j$ אלא גם בבחירת נוסחת המכפלה (PF) שלנו -- ובפרט ב_סדר_ שלה. בנוסף, ניתן לקחת בחשבון האם ה-PF סימטרי או לא (ראה Carrera Vazquez et al., 2023), על ידי הגדרת symmetric=True. אולם זה אינו חובה כפי שמוצג ב-Zhuk et al., 2023.

כאן נשתמש בנוסחת סוזוקי-טרוטר מסדר שני המניבה order=2 ונגדיר symmetric=True.

from qiskit_addon_mpf.static import setup_static_lse

lse = setup_static_lse(trotter_steps, order=2, symmetric=True)
print(lse)

LSE(A=array([[1.00000000e+00, 1.00000000e+00, 1.00000000e+00],
       [1.56250000e-02, 6.94444444e-03, 2.77008310e-03],
       [2.44140625e-04, 4.82253086e-05, 7.67336039e-06]]), b=array([1., 0., 0.]))

פתרון אנליטי של $x$

כפי שהוזכר קודם, ניתן למצוא את $x$ אנליטית:

import numpy as np

coeffs_analytical = lse.solve()
print(coeffs_analytical)

[ 0.17239057 -1.19447005  2.02207947]

אופטימיזציה של $x$ באמצעות מודל מדויק

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

בחלק הבא יתברר מדוע ממשק זה קיים.

from qiskit_addon_mpf.costs import setup_exact_problem

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

[ 0.17239057 -1.19447005  2.02207947]

כמדד לשאלה האם MPF שנבנית עם מקדמים אלו תניב תוצאות טובות, ניתן להשתמש בנורמת ה-L1 (ראה גם Carrera Vazquez et al., 2023).

print(np.linalg.norm(coeffs_exact.value, ord=1))

3.3889400921655914

אופטימיזציה של $x$ באמצעות מודל מקורב

ייתכן שנורמת ה-L1 עבור קבוצת ערכי $k_j$ הנבחרת תיחשב גבוהה מדי. אם כך הוא המצב ולא ניתן לבחור קבוצת ערכי $k_j$ שונה, ניתן להשתמש בפתרון מקורב ל-LSE במקום בפתרון מדויק.

לשם כך, פשוט השתמש/י ב-setup_sum_of_squares_problem לבניית מופע cvxpy.Problem שונה המגביל את נורמת ה-L1 לסף נבחר תוך מינימיזציה של ההפרש בין $Ax$ ל-$b$.

from qiskit_addon_mpf.costs import setup_sum_of_squares_problem

model_approx, coeffs_approx = setup_sum_of_squares_problem(lse, max_l1_norm=3.0)
model_approx.solve()
print(coeffs_approx.value)
print(np.linalg.norm(coeffs_approx.value, ord=1))

[-0.40454257  0.57553173  0.8290123 ]
1.8090865903790838

שים/י לב שיש לך חופש מלא על אופן פתרון בעיית האופטימיזציה, כלומר ניתן לשנות את ה-solver, סף ההתכנסות שלו, וכן הלאה. ראה/י במדריך הרלוונטי How to use the approximate model.

1ג: הגדרת מעגלי הטרוטר

בשלב זה מצאנו את מקדמי הפיתוח $x$, וכל שנותר הוא לייצר את מעגלי הקוונטום בשיטת טרוטר. שוב, המודול qiskit_addon_utils.problem_generators בא לעזרתנו בדיוק לשם כך:

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit

circuits = []
for k in trotter_steps:
    circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(order=2, reps=k),
        time=time,
    )
    circuits.append(circ)

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

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

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

שלב 2: אופטימיזציה של הבעיה

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

from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.transpiler import generate_preset_pass_manager

backend = GenericBackendV2(num_qubits=10)
transpiler = generate_preset_pass_manager(optimization_level=2, backend=backend)

transpiled_circuits = [transpiler.run(circ) for circ in circuits]

שלב 3: ביצוע ניסויים קוונטיים

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

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()
job = estimator.run([(circ, observable) for circ in transpiled_circuits])
result = job.result()

שלב 4: שחזור תוצאות

תחילה, נחלץ את ערכי הציפייה הבודדים שהתקבלו עבור כל אחד ממעגלי הטרוטר:

evs = [res.data.evs for res in result]
print(evs)

[array(0.23799162), array(0.35754312), array(0.38649906)]

לאחר מכן, פשוט נשלב אותם עם מקדמי ה-MPF שלנו כדי לקבל את ערכי הציפייה הכוללים של ה-MPF. להלן נעשה זאת עבור כל אחת מהדרכים השונות שבהן חישבנו את $x$.

print("Analytical    solution:", evs @ coeffs_analytical)
print("Exact model   solution:", evs @ coeffs_exact.value)
print("Approx. model solution:", evs @ coeffs_approx.value)

Analytical    solution: 0.3954847855980006
Exact model   solution: 0.39548478559800204
Approx. model solution: 0.42991214253489807

לבסוף, עבור בעיה קטנה זו ניתן לחשב את ערך הייחוס המדויק באמצעות scipy.linalg.expm כדלקמן:

from scipy.linalg import expm

exp_H = expm(-1j * time * hamiltonian.to_matrix())

initial_state = np.zeros(exp_H.shape[0])
initial_state[0] = 1.0

time_evolved_state = exp_H @ initial_state

exact_obs = time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
print(exact_obs.real)

0.40060242487899755

ניתן לראות בבירור שה-MPF הצליח לצמצם את שגיאת הטרוטר בהשוואה לזו שהתקבלה עם ה-PF הבודד העמוק ביותר עם $k_j=19$. אולם אנחנו רואים גם שהמודל המקורב אינו ללא פגמים, שכן הוא הניב ערך ציפייה גרוע יותר מהפתרון המדויק. זה מדגיש את החשיבות של שימוש בקריטריוני התכנסות הדוקים במודל המקורב, כפי שתלמד/י במדריך How to use the approximate model.