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

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

רקע

מדריך זה מדגים כיצד להשתמש בנוסחת מכפלת-רב (MPF) כדי להשיג שגיאת Trotter נמוכה יותר על האובזרבל שלנו בהשוואה לזו שנגרמת על ידי מעגל ה-Trotter העמוק ביותר שבאמת נריץ. MPF מפחיתים את שגיאת ה-Trotter של דינמיקת המילטוניאן דרך שילוב משוקלל של מספר ביצועי מעגל. שקול את המשימה של מציאת ערכי ציפייה של אובזרבל עבור המצב הקוונטי $\rho(t)=e^{-i H t} \rho(0) e^{i H t}$ עם ההמילטוניאן $H$. אפשר להשתמש בנוסחאות מכפלה (PF) כדי לקרב את אבולוציית הזמן $e^{-i H t}$ על ידי ביצוע הדברים הבאים:

• כתוב את ההמילטוניאן $H$ כ-$H=\sum_{a=1}^d F_a,$ כאשר $F_a$ הם אופרטורים הרמיטיים כך שכל יוניטרי מתאים יכול להיות מיושם ביעילות על מכשיר קוונטי.
• קרב את האיברים $F_a$ שאינם מתחלפים זה עם זה.

אז, ה-PF מסדר ראשון (נוסחת Lie-Trotter) הוא:

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

שיש לו איבר שגיאה ריבועי $S_1(t)=e^{-i H t}+\mathcal{O}\left(t^{2}\right)$. אפשר גם להשתמש ב-PF מסדרים גבוהים יותר (נוסחאות Lie-Trotter-Suzuki), שמתכנסים מהר יותר, ומוגדרים באופן רקורסיבי כ:

$S_2(t):=\prod_{a=1}^d e^{-i F_a t/2}\prod_{a=1}^d e^{-i F_a t/2}$

$S_{2 \chi}(t):= S_{2 \chi -2}(s_{\chi}t)^2 S_{2 \chi -2}((1-4s_{\chi})t)S_{2 \chi -2}(s_{\chi}t)^2,$

כאשר $\chi$ הוא הסדר של ה-PF הסימטרי ו-$s_p = \left( 4 - 4^{1/(2p-1)} \right)^{-1}$. עבור אבולוציות זמן ארוכות, אפשר לפצל את מרווח הזמן $t$ ל-$k$ מרווחים, הנקראים צעדי Trotter, בעלי משך $t/k$ ולקרב את אבולוציית הזמן בכל מרווח עם נוסחת מכפלה מסדר $\chi$ $S_{\chi}$. לפיכך, ה-PF מסדר $\chi$ עבור אופרטור אבולוציית זמן על פני $k$ צעדי Trotter הוא:

$S_{\chi}^{k}(t) = \left[ S_{\chi} \left( \frac{t}{k} \right)\right]^k = e^{-i H t}+O\left(t \left( \frac{t}{k} \right)^{\chi} \right)$

כאשר איבר השגיאה פוחת עם מספר צעדי ה-Trotter $k$ והסדר $\chi$ של ה-PF.

בהינתן מספר שלם $k \geq 1$ ונוסחת מכפלה $S_{\chi}(t)$, המצב המתפתח בזמן המקורב $\rho_k(t)$ ניתן לקבלה מ-$\rho_0$ על ידי הפעלת $k$ איטרציות של נוסחת המכפלה $S_{\chi}\left(\frac{t}{k}\right)$.

$$\rho_k(t)=S_{\chi}\left(\frac{t}{k}\right)^k \rho_0 S_{\chi}\left(\frac{t}{k}\right)^{-k}$$

$\rho_k(t)$ הוא קירוב ל-$\rho(t)$ עם שגיאת קירוב Trotter ||$\rho_k(t)-\rho(t) ||$. אם אנחנו שוקלים שילוב לינארי של קירובי Trotter של $\rho(t)$:

$$\mu(t) = \sum_{j}^{l} 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}_{\chi}$, הכוללת $k_j$ צעדי Trotter, ו-$j \in {1, ..., l}$ מאנדקס את מספר ה-PF שמרכיבים את ה-MPF. כל האיברים ב-$\mu(t)$ משתמשים באותה נוסחת מכפלה $S_{\chi}(t)$ כבסיס שלהם. המטרה היא לשפר את ||$\rho_k(t)-\rho(t) \|$ על ידי מציאת $\mu(t)$ עם $\|\mu(t)-\rho(t)\|$ אפילו נמוך יותר.

• $\mu(t)$ אינו חייב להיות מצב פיזיקלי מכיוון ש-$x_i$ אינם חייבים להיות חיוביים. המטרה כאן היא למזער את השגיאה בערך הצפייה של האובזרבלים ולא למצוא תחליף פיזיקלי ל-$\rho(t)$.
• $k_j$ קובע גם את עומק המעגל וגם את רמת קירוב Trotter. ערכים קטנים יותר של $k_j$ מובילים למעגלים קצרים יותר, שגורמים לפחות שגיאות מעגל אך יהיו קירוב פחות מדויק למצב הרצוי.

המפתח כאן הוא ששגיאת ה-Trotter הנותרת שניתנת על ידי $\mu(t)$ קטנה יותר משגיאת ה-Trotter שהיינו מקבלים פשוט על ידי שימוש בערך ה-$k_j$ הגדול ביותר.

אתה יכול לראות את התועלת של זה משתי נקודות מבט:

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

דרישות

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

• Qiskit SDK v1.0 ומעלה, עם תמיכה ב-visualization
• Qiskit Runtime v0.22 ומעלה (pip install qiskit-ibm-runtime)
• תוספי MPF Qiskit (pip install qiskit_addon_mpf)
• תוספי utils של Qiskit (pip install qiskit_addon_utils)
• ספריית Quimb (pip install quimb)
• ספריית Qiskit Quimb (pip install qiskit-quimb)
• Numpy v0.21 עבור תאימות בין חבילות (pip install numpy==0.21)

חלק I. דוגמה בקנה מידה קטן

חקור את היציבות של MPF

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

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

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

import numpy as np

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)

עבור דוגמה זו נשתמש במצב Neel כמצב ההתחלתי $\vert \text{Neel} \rangle = \vert 0101...01 \rangle$ ובמודל Heisenberg על קו של 10 אתרים עבור ההמילטוניאן השולט באבולוציית הזמן

$$\hat{\mathcal{H}}_{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$ הוא חוזק הצימוד עבור קצוות שכן-קרוב ביותר.

from qiskit.transpiler import CouplingMap
from rustworkx.visualization import graphviz_draw
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian
import numpy as np

L = 10

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

# Get a qubit operator describing the Heisenberg field model
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])

האובזרבל שנמדוד הוא מגנטיזציה על זוג קיוביטים באמצע השרשרת.

from qiskit.quantum_info import SparsePauliOp

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

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

אנחנו מגדירים מעבר טרנספילר לאסוף את סיבובי ה-XX וה-YY במעגל כשער XX+YY יחיד. זה יאפשר לנו למנף את תכונות שימור הספין של TeNPy במהלך חישוב ה-MPO, להאיץ משמעותית את החישוב.

from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler import PassManager
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)
from functools import partial

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

אז אנחנו יוצרים את המעגלים המיישמים את אבולוציות הזמן המקורבות של Trotter.

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

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

לאחר מכן, אנחנו מחשבים את ערכי הצפייה המתפתחים בזמן ממעגלי ה-Trotter.

from copy import deepcopy
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

aer_sim = AerSimulator()
estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for t, mpf_circuits in zip(trotter_times, all_circs):
    mpf_expvals = []
    circuits = [deepcopy(circuit) for circuit in mpf_circuits]
    pm_sim = generate_preset_pass_manager(
        backend=aer_sim, optimization_level=3
    )
    isa_circuits = pm_sim.run(circuits)
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals = [res.data.evs for res in result]
    mpf_stds = [res.data.stds for res in result]
    mpf_expvals_all_times.append(mpf_expvals)
    mpf_stds_all_times.append(mpf_stds)

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

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

exact_expvals = []
for t in exact_evolution_times:
    # Exact expectation values
    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 סטטיים

MPF סטטיים הם אלה שבהם ערכי ה-$x_j$ אינם תלויים בזמן האבולוציה, $t$. בואו נשקול את ה-PF מסדר $\chi = 1$ עם $k_j$ צעדי Trotter, זה יכול להיכתב כ:

$S_1^{k_j}\left( \frac{t}{k_j} \right)=e^{-i H t}+ \sum_{n=1}^{\infty} A_n \frac{t^{n+1}}{k_j^n}$

כאשר $A_n$ הן מטריצות שתלויות בקומוטטורים של איברי $F_a$ בפירוק של ההמילטוניאן. חשוב לציין ש-$A_n$ עצמם אינם תלויים בזמן ובמספר צעדי ה-Trotter $k_j$. לכן, אפשר לבטל איברי שגיאה מסדר נמוך יותר התורמים ל-$\mu(t)$ עם בחירה קפדנית של המשקלים $x_j$ של השילוב הלינארי. כדי לבטל את שגיאת ה-Trotter עבור $l-1$ האיברים הראשונים (אלה יתנו את התרומות הגדולות ביותר מכיוון שהם מתאימים למספר הנמוך יותר של צעדי Trotter) בביטוי עבור $\mu(t)$, המקדמים $x_j$ חייבים לספק את המשוואות הבאות:

$\sum_{j=1}^l x_j = 1$ $\sum_{j=1}^{l-1} \frac{x_j}{k_j^{n}} = 0$

עם $n=0, ... l-2$. המשוואה הראשונה מבטיחה שאין הטיה במצב הבנוי $\mu(t)$, בעוד שהמשוואה השנייה מבטיחה את הביטול של שגיאות ה-Trotter. עבור PF מסדר גבוה יותר, המשוואה השנייה הופכת ל-$\sum_{j=1}^{l-1} \frac{x_j}{k_j^{\eta}} = 0$ כאשר $\eta = \chi + 2n$ עבור PF סימטריים ו-$\eta = \chi + n$ אחרת, עם $n=0, ..., l-2$. השגיאה המתקבלת (ראה [1],[2]) היא אז

$\epsilon = \mathcal{O} \left( \frac{t^{l+1}}{k_1^l} \right).$

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

$A_{0,j} = 1$ $A_{i>0,j} = k_{j}^{-(\chi + s(i-1))}$ $b_0 = 1$ $b_{i>0} = 0$

כאשר $\chi$ הוא ה-order, $s$ הוא $2$ אם symmetric הוא True ו-$1$ אחרת, $k_{j}$ הם ה-trotter_steps, ו-$x$ הם המשתנים לפתור עבורם. האינדקסים $i$ ו-$j$ מתחילים ב-$0$. אנחנו יכולים גם לדמיין זאת בצורת מטריצה:

$$A = \begin{bmatrix} A_{0,0} & A_{0,1} & A_{0,2} & ... \\ A_{1,0} & A_{1,1} & A_{1,2} & ... \\ A_{2,0} & A_{2,1} & A_{2,2} & ... \\ ... & ... & ... & ... \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & ... \\ k_{0}^{-(\chi + s(1-1))} & k_{1}^{-(\chi + s(1-1))} & k_{2}^{-(\chi + s(1-1))} & ... \\ k_{0}^{-(\chi + s(2-1))} & k_{1}^{-(\chi + s(2-1))} & k_{2}^{-(\chi + s(2-1))} & ... \\ ... & ... & ... & ... \end{bmatrix}$$

ו

$$b = \begin{bmatrix} b_{0} \\ b_{1} \\ b_{2} \\ ... \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ ... \end{bmatrix}$$

לפרטים נוספים, ראה בתיעוד של מערכת המשוואות הלינאריות (LSE).

אנחנו יכולים למצוא פתרון עבור $x$ באופן אנליטי כ-$x = A^{-1}b$; ראה לדוגמה [1] או [2]. עם זאת, הפתרון המדויק הזה יכול להיות "לא-מותנה היטב", מה שגורם לנורמות L1 גדולות מאוד של המקדמים שלנו, $x$, שיכולות להוביל לביצועים גרועים של ה-MPF. במקום זאת, אפשר גם לקבל פתרון מקורב שממזער את נורמת ה-L1 של $x$ כדי לנסות לייעל את התנהגות ה-MPF.

הגדר את ה-LSE

כעת לאחר שבחרנו את ערכי ה-$k_j$ שלנו, עלינו קודם לבנות את ה-LSE, $Ax=b$ כפי שהוסבר לעיל. המטריצה $A$ תלויה לא רק ב-$k_j$ אלא גם בבחירת ה-PF שלנו, במיוחד ה-order שלו. בנוסף, אתה עשוי לקחת בחשבון האם ה-PF סימטרי או לא (ראה [1]) על ידי הגדרת symmetric=True/False. עם זאת, זה לא נדרש, כפי שמוצג על ידי [2].

from qiskit_addon_mpf.static import setup_static_lse

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

בואו נעבור דרך הערכים שנבחרו למעלה כדי לבנות את מטריצת ה-$A$ ואת וקטור ה-$b$. עם $j=0,1, 2$ צעדי Trotter $k_j = [1, 2, 4]$, סדר $\chi = 2$ ובחירה של צעדי Trotter לא-סימטריים ($s=1$), יש לנו שאלמנטי המטריצה של $A$ מתחת לשורה הראשונה נקבעים על ידי הביטוי $A_{i>0,j} = k_{j}^{-(2 + 1(i-1))}$, באופן ספציפי:

$A_{0,0} = A_{0,1} = A_{0,2} = 1$ $A_{1,j} = k_{j}^{-1} \rightarrow A_{1,0} = \frac{1}{1^2}, \;, A_{1,1} = \frac{1}{2^2}, \;, A_{1,2} = \frac{1}{4^2}$ $A_{2,j} = k_{j}^{-2} \rightarrow A_{2,0} = \frac{1}{1^3}, \;, A_{2,1} = \frac{1}{2^3}, \;, A_{2,2} = \frac{1}{4^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}$$

זה אפשרי לראות על ידי בדיקת האובייקט lse:

lse.A

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

בעוד שלוקטור האילוצים $b$ יש את האלמנטים הבאים: $b_{0} = 1$ $b_1 = b_2 = 0$

לפיכך,

$$b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$$

ובאופן דומה ב-lse:

lse.b

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

לאובייקט lse יש מתודות למציאת המקדמים הסטטיים $x_j$ המספקים את מערכת המשוואות.

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$.

בסעיף הבא, יהיה ברור מדוע הממשק הזה קיים.

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.04761905 -0.57142857  1.52380952]

כאינדיקטור האם MPF שנבנה עם מקדמים אלה ייתן תוצאות טובות, אנחנו יכולים להשתמש בנורמת L1 (ראה גם [1]).

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

L1 norm of the exact coefficients: 2.1428571428556378

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

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

כדי לעשות זאת, פשוט השתמש ב-setup_approximate_model כדי לבנות מופע 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=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

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

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

בסעיף הקודם, הצגנו MPF סטטי שמשפר על קירוב ה-Trotter הסטנדרטי. עם זאת, הגרסה הסטטית הזו אינה בהכרח ממזערת את שגיאת הקירוב. באופן קונקרטי, ה-MPF הסטטי, המסומן $\mu^S(t)$, אינו ההיטל האופטימלי של $\rho(t)$ על תת-המרחב הנפרש על ידי מצבי נוסחת-המכפלה $\{\rho_{k_i}(t)\}_{i=1}^r$.

כדי לטפל בכך, אנחנו שוקלים MPF דינמי (שהוצג ב-[2] והוכח ניסיונית ב-[3]) שכן ממזער את שגיאת הקירוב בנורמת Frobenius. באופן פורמלי, אנחנו מתמקדים במזעור

$$\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr],$$

ביחס למקדמים מסוימים $x_i(t)$ בכל זמן $t$. המקרין האופטימלי בנורמת Frobenius הוא אז $\mu^D(t) \;=\; \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t)$, ואנו קוראים ל-$\mu^D(t)$ ה-MPF הדינמי. הצבת ההגדרות למעלה:

$$\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr] \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t) \right) \left( \rho(t) - \sum_{j=1}^r x_j(t)\,\rho_{k_j}(t) \right) \bigr] \;=\; \\ = 1 \;+\; \sum_{i,j=1}^r M_{i,j}(t)\,x_i(t)\,x_j(t) \;-\; 2 \sum_{i=1}^r L_i^{\mathrm{exact}}(t)\,x_i(t),$$

כאשר $M_{i,j}(t)$ היא מטריצת Gram, המוגדרת על ידי

$$M_{i,j}(t) \;=\; \mathrm{Tr}\bigl[\rho_{k_i}(t)\,\rho_{k_j}(t)\bigr] \;=\; \bigl|\langle \psi_{\mathrm{in}} \!\mid S\bigl(t/k_i\bigr)^{-k_i}\,S\bigl(t/k_j\bigr)^{k_j} \!\mid \psi_{\mathrm{in}} \rangle \bigr|^2.$$

ו

$$L_i^{\mathrm{exact}}(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$$

מייצג את החפיפה בין המצב המדויק $\rho(t)$ וכל קירוב נוסחת-מכפלה $\rho_{k_i}(t)$. בתרחישים מעשיים, חפיפות אלו עשויות להימדד רק באופן מקורב בגלל רעש או גישה חלקית ל-$\rho(t)$.

כאן, $\lvert\psi_{\mathrm{in}}\rangle$ הוא המצב ההתחלתי, ו-$S(\cdot)$ היא הפעולה המופעלת בנוסחת המכפלה. על ידי בחירת המקדמים $x_i(t)$ שממזערים את הביטוי הזה (וטיפול בנתוני חפיפה מקורבים כאשר $\rho(t)$ אינו ידוע במלואו), אנחנו מקבלים את הקירוב הדינמי "הטוב ביותר" (במובן נורמת Frobenius) של $\rho(t)$ בתוך תת-מרחב ה-MPF. הכמויות $L_i(t)$ ו-$M_{i,j}(t)$ ניתנות לחישוב ביעילות באמצעות שיטות רשת טנזורית [3]. תוסף ה-MPF של Qiskit מספק מספר "backends" לביצוע החישוב. הדוגמה למטה מראה את הדרך הגמישה ביותר לעשות זאת, והתיעוד של TeNPy layer-based backend גם מסביר בפירוט רב. כדי להשתמש בשיטה זו, התחל מהמעגל המיישם את אבולוציית הזמן הרצויה וצור מודלים המייצגים את הפעולות האלה מהשכבות של המעגל המתאים. לבסוף, נוצר אובייקט Evolver שניתן להשתמש בו כדי לייצר את הכמויות המתפתחות בזמן $M_{i,j}(t)$ ו-$L_i(t)$. אנחנו מתחילים ביצירת האובייקט Evolver המתאים לאבולוציית הזמן המקורבת (ApproxEvolverFactory) המיושמת על ידי המעגלים. במיוחד, שים תשומת לב מיוחדת למשתנה order כך שהם יתאימו. שים לב שביצירת המעגלים המתאימים לאבולוציית הזמן המקורבת, אנחנו משתמשים בערכי placeholder עבור time = 1.0 ומספר צעדי Trotter (reps=1). המעגלים המקרבים הנכונים מיוצרים אז על ידי פותר הבעיה הדינמית ב-setup_dynamic_lse.

from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.backends.tenpy_layers import LayerModel
from qiskit_addon_mpf.backends.tenpy_layers import LayerwiseEvolver
from functools import partial

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

לאחר מכן אנחנו מגדירים את ה-evolver המדויק (לדוגמה, ExactEvolverFactory), שמחזיר אובייקט Evolver המחשב את אבולוציית הזמן האמיתית או "ייחוס". באופן ריאלי, היינו מקרבים את האבולוציה המדויקת על ידי שימוש בנוסחת Suzuki-Trotter מסדר גבוה יותר או שיטה אמינה אחרת עם צעד זמן קטן. למטה, אנחנו מקרבים את המצב המדויק שהתפתח בזמן עם נוסחת Suzuki-Trotter מסדר רביעי באמצעות צעד זמן קטן dt=0.1, מה שאומר שמספר צעדי ה-Trotter המשמשים בזמן $t$ הוא $k=t/dt$. אנחנו גם מציינים כמה אופציות חיתוך ספציפיות ל-TeNPy כדי לתחום את מימד הקשר המקסימלי של רשת הטנזורים הבסיסית, כמו גם את הערכים הסינגולריים המינימליים של קשרי רשת הטנזורים המפוצלים. פרמטרים אלה יכולים להשפיע על הדיוק של ערך הצפייה המחושב עם מקדמי ה-MPF הדינמיים, לכן חשוב לחקור טווח של ערכים כדי למצוא את האיזון האופטימלי בין זמן חישוב ודיוק. שים לב שחישוב מקדמי ה-MPF אינו מסתמך על ערך הצפייה של ה-PF שהתקבל מביצוע חומרה, ולכן ניתן לכוונן אותו בעיבוד-לאחר.

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,
    },
)

לאחר מכן, צור את המצב ההתחלתי של המערכת שלך בפורמט תואם ל-TeNPy (לדוגמה, MPS_neel_state=$\vert 0101...01 \rangle$). זה מגדיר את פונקציית הגל של גוף-רב שתתפתח בזמן $\lvert\psi_{\mathrm{in}}\rangle$ כטנזור.

from qiskit_addon_mpf.backends.tenpy_tebd import MPOState
from qiskit_addon_mpf.backends.tenpy_tebd import MPS_neel_state

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. האובייקט המתאים מכיל את המידע על בעיית ה-MPF הדינמית: lse.A נותן את מטריצת Gram $M$ בעוד lse.b נותן את החפיפה $L$. אנחנו יכולים אז לפתור את ה-LSE (כאשר לא מוגדר לא-היטב) כדי למצוא את המקדמים הדינמיים באמצעות setup_frobenius_problem. חשוב לציין את ההבדל עם המקדמים הסטטיים, שתלויים רק בפרטים של נוסחת המכפלה המשמשת ואינם תלויים בפרטים של אבולוציית הזמן (המילטוניאן ומצב התחלתי).

from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import setup_frobenius_problem

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

לבסוף, תרשם את ערכי הצפייה האלה לאורך זמן האבולוציה.

import matplotlib.pyplot as plt

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
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",
    color="orange",
    label="Static MPF - Approximate",
)

# # 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",
    color="pink",
    label="Dynamic MPF",
)

plt.plot(
    exact_evolution_times,
    exact_expvals,
    lw=3,
    color="red",
    label="Exact time-evolution",
)

plt.title(
    f"Expectation values for (ZZ,{(L//2-1, L//2)}) as a function of time"
)
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend()
plt.grid()

במקרים כמו הדוגמה למעלה, שבהם ה-PF $k=1$ מתנהג בצורה גרועה בכל הזמנים, איכות תוצאות ה-MPF הדינמיות גם מושפעת מאוד. במצבים כאלה, שימושי לחקור את האפשרות של שימוש ב-PF בודדים עם מספר גבוה יותר של צעדי Trotter כדי לשפר את האיכות הכוללת של התוצאות. בסימולציות אלו, אנחנו רואים את השילוב של סוגים שונים של שגיאות: שגיאה מדגימה סופית, ושגיאת Trotter מנוסחאות המכפלה. MPF עוזר להפחית את שגיאת ה-Trotter בגלל נוסחאות המכפלה אך גורם לשגיאת דגימה גבוהה יותר בהשוואה לנוסחאות המכפלה. זה יכול להיות יתרון, מכיוון שנוסחאות מכפלה יכולות להפחית את שגיאת הדגימה עם דגימה מוגברת, אך השגיאה השיטתית בגלל קירוב ה-Trotter נשארת ללא נגיעה.

התנהגות מעניינת נוספת שאנו יכולים לראות מהתרשים היא שערך הצפייה עבור ה-PF עבור $k=1$ מתחיל להתנהג בצורה לא סדירה (על גבי היותו לא קירוב טוב למצב המדויק) בזמנים שבהם $t/k > 1$, כפי שהוסבר ב-מדריך על איך לבחור את מספר צעדי ה-Trotter.

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

בואו נשקול כעת זמן בודד $t=1.0$ ונחשב את ערך הצפייה של המגנטיזציה עם השיטות השונות תוך שימוש ב-QPU אחד. הבחירה הספציפית של $t$ נעשתה כדי למקסם את ההבדל בין השיטות השונות ולהתבונן ביעילות היחסית שלהן. כדי לקבוע את חלון הזמן שבו MPF דינמי מובטח לייצר אובזרבלים עם שגיאה נמוכה יותר מכל אחד מנוסחאות ה-Trotter הבודדות בתוך המכפלת-הרב, אנחנו יכולים ליישם את "בדיקת MPF" - ראה משוואה (17) והטקסט הסובב אותה ב-[3].

הגדר את מעגלי Trotter

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

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

total_time = 1.0
mpf_circuits = []
for k in mpf_trotter_steps:
    # Initial Neel state preparation
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2 != 0])

    trotter_circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(order=order, reps=k),
        time=total_time,
    )

    circuit.compose(trotter_circ, inplace=True)

    mpf_circuits.append(pm.run(circuit))

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

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

בואו נחזור לחישוב ערך הצפייה עבור נקודת זמן בודדת. נבחר backend לביצוע הניסוי על חומרה.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
print(backend)

qubits = list(range(backend.num_qubits))

אז אנחנו מסירים קיוביטים חריגים ממפת הצימוד כדי להבטיח ששלב ה-layout של הטרנספילר לא יכלול אותם. למטה אנחנו משתמשים במאפייני ה-backend המדווחים המאוחסנים באובייקט target ומסירים קיוביטים שיש להם שגיאת מדידה או שער דו-קיוביטי מעל סף מסוים (max_meas_err, max_twoq_err) או זמן $T_2$ (שקובע את אובדן הקוהרנטיות) מתחת לסף מסוים (min_t2).

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

max_meas_err = 0.012
min_t2 = 40
max_twoq_err = 0.005

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    if target.qubit_properties[q].t2 is not None:
        t2 = target.qubit_properties[q].t2 * 1e6
    else:
        t2 = 0
    if meas_err > max_meas_err or t2 < min_t2:
        # print(q)
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
    twoq_gate_err = target[instruction_2q][q].error
    if twoq_gate_err > max_twoq_err:
        # print(q)
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates
    + ["measure"],  # or whatever new set of gates
    coupling_map=cust_cmap,
)

sorted_components = sorted(
    [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
    reverse=True,
)
print("size of largest component", len(sorted_components[0]))

size of largest component 10

אנחנו רוצים להגדיר את max_meas_err, min_t2 ו-max_twoq_err כך שנמצא תת-קבוצה גדולה מספיק של קיוביטים התומכת במעגל להרצה. במקרה שלנו מספיק למצוא שרשרת 1D של 10 קיוביטים.

cust_cmap.draw()

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

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
    optimization_level=3, target=cust_target
)

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

qubits_layouts = [
    [
        idx
        for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
        if qb._register.name != "ancilla"
    ]
    for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
    transpiler = generate_preset_pass_manager(
        optimization_level=3, backend=backend, initial_layout=layout
    )
    transpiled_circuit = transpiler.run(circuit)
    transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

print(transpiled_circuits[-1].depth(lambda x: x.operation.num_qubits == 2))
print(transpiled_circuits[-1].count_ops())
transpiled_circuits[-1].draw("mpl", idle_wires=False, fold=False)

51
OrderedDict([('sx', 310), ('rz', 232), ('cz', 132), ('x', 19)])

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

עם הפרימיטיב Estimator אנחנו יכולים לקבל את ההערכה של ערך הצפייה מה-QPU. אנחנו מבצעים את מעגלי ה-AQC המותאמים עם טכניקות נוספות של הפחתת שגיאות ודיכוי.

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
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, 3, 5)
estimator.options.resilience.zne.extrapolator = ("exponential", "linear")

estimator.options.environment.job_tags = ["mpf small"]

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

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

שלב העיבוד-לאחר היחיד הוא לשלב את ערך הצפייה שהתקבל מהפרימיטיבים של Qiskit Runtime בצעדי Trotter שונים תוך שימוש במקדמי ה-MPF המתאימים. עבור אובזרבל $A$ יש לנו:

$\langle A \rangle_{\text{mpf}} = \text{Tr} [A \mu(t)] = \sum_{j} x_j \text{Tr} [A \rho_{k_j}] = \sum_{j} x_j \langle A \rangle_j$

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

result_exp = job.result()
evs_exp = [res.data.evs for res in result_exp]
evs_std = [res.data.stds for res in result_exp]

print(evs_exp)

[array(-0.06361607), array(-0.23820448), array(-0.50271805)]

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

exact_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_exact.value, evs_std)
        ]
    )
)
print(
    "Exact static MPF expectation value: ",
    evs_exp @ coeffs_exact.value,
    "+-",
    exact_mpf_std,
)
approx_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_approx.value, evs_std)
        ]
    )
)
print(
    "Approximate static MPF expectation value: ",
    evs_exp @ coeffs_approx.value,
    "+-",
    approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(mpf_dynamic_coeffs_list[7], evs_std)
        ]
    )
)
print(
    "Dynamic MPF expectation value: ",
    evs_exp @ mpf_dynamic_coeffs_list[7],
    "+-",
    dynamic_mpf_std,
)

Exact static MPF expectation value:  -0.6329590442738475 +- 0.012798249760406036
Approximate static MPF expectation value:  -0.5690390035339492 +- 0.010459559917168473
Dynamic MPF expectation value:  -0.4655579758795695 +- 0.007639139186720507

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

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

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

initial_state_circuit = QuantumCircuit(L)
initial_state_circuit.x([i for i in range(L) if i % 2 != 0])
initial_state = Statevector(initial_state_circuit).data

time_evolved_state = exp_H @ initial_state

exact_obs = (
    time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
)
print("Exact expectation value ", exact_obs.real)

Exact expectation value  -0.39909900734489434

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs_exp[k],
        yerr=evs_std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )  # , , )

plt.errorbar(
    3,
    evs_exp @ coeffs_exact.value,
    yerr=exact_mpf_std,
    markersize=4,
    marker="o",
    color="purple",
    label="Static MPF",
)

plt.errorbar(
    4,
    evs_exp @ coeffs_approx.value,
    yerr=approx_mpf_std,
    markersize=4,
    marker="o",
    color="orange",
    label="Approximate static MPF",
)

plt.errorbar(
    5,
    evs_exp @ mpf_dynamic_coeffs_list[7],
    yerr=dynamic_mpf_std,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

plt.axhline(
    y=exact_obs.real,
    linestyle="--",
    color="red",
    label="Exact time-evolution",
)

plt.title(
    f"Expectation values for (ZZ,{(L//2-1, L//2)}) 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()

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

def relative_error(ev, exact_ev):
    return abs(ev - exact_ev)

relative_error_k = [relative_error(ev, exact_obs.real) for ev in evs_exp]
relative_error_mpf = relative_error(evs_exp @ mpf_coeffs, exact_obs.real)
relative_error_approx_mpf = relative_error(
    evs_exp @ coeffs_approx.value, exact_obs.real
)
relative_error_dynamic_mpf = relative_error(
    evs_exp @ mpf_dynamic_coeffs_list[7], exact_obs.real
)

print("relative error for each trotter steps", relative_error_k)
print("relative error with MPF exact coeffs", relative_error_mpf)
print("relative error with MPF approx coeffs", relative_error_approx_mpf)
print("relative error with MPF dynamic coeffs", relative_error_dynamic_mpf)

relative error for each trotter steps [0.33548293650112293, 0.16089452939226306, 0.10361904247828346]
relative error with MPF exact coeffs 0.2338600369291003
relative error with MPF approx coeffs 0.16993999618905486
relative error with MPF dynamic coeffs 0.06645896853467514

חלק II: הגדל את קנה המידה

בואו נגדיל את הבעיה מעבר למה שאפשר לדמות באופן מדויק. בסעיף זה נתמקד בשחזור חלק מהתוצאות המוצגות ב-[3].

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

המילטוניאן

עבור הדוגמה בקנה מידה גדול, אנחנו משתמשים במודל XXZ על קו של 50 אתרים:

$$\hat{\mathcal{H}}_{XXZ} = \sum_{i=1}^{L-1} J_{i,(i+1)}\left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ 2\cdot Z_i Z_{(i+1)} \right) \, ,$$

כאשר $J_{i,(i+1)}$ הוא מקדם אקראי המתאים לקצה $(i, i+1)$. זה ההמילטוניאן שנשקל בהדגמה המוצגת ב-[3].

L = 50
# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

import numpy as np
from qiskit.quantum_info import SparsePauliOp, Pauli

# Generate random coefficients for our XXZ Hamiltonian
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,
    )

print(hamiltonian)

SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'XXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'YYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'ZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII'],
              coeffs=[1.        +0.j, 2.09762701+0.j, 2.09762701+0.j, 4.19525402+0.j,
 2.43037873+0.j, 2.43037873+0.j, 4.86075747+0.j, 2.20552675+0.j,
 2.20552675+0.j, 4.4110535 +0.j, 2.08976637+0.j, 2.08976637+0.j,
 4.17953273+0.j, 1.8473096 +0.j, 1.8473096 +0.j, 3.6946192 +0.j,
 2.29178823+0.j, 2.29178823+0.j, 4.58357645+0.j, 1.87517442+0.j,
 1.87517442+0.j, 3.75034885+0.j, 2.783546  +0.j, 2.783546  +0.j,
 5.567092  +0.j, 2.92732552+0.j, 2.92732552+0.j, 5.85465104+0.j,
 1.76688304+0.j, 1.76688304+0.j, 3.53376608+0.j, 2.58345008+0.j,
 2.58345008+0.j, 5.16690015+0.j, 2.05778984+0.j, 2.05778984+0.j,
 4.11557968+0.j, 2.13608912+0.j, 2.13608912+0.j, 4.27217824+0.j,
 2.85119328+0.j, 2.85119328+0.j, 5.70238655+0.j, 1.14207212+0.j,
 1.14207212+0.j, 2.28414423+0.j, 1.1742586 +0.j, 1.1742586 +0.j,
 2.3485172 +0.j, 1.04043679+0.j, 1.04043679+0.j, 2.08087359+0.j,
 2.66523969+0.j, 2.66523969+0.j, 5.33047938+0.j, 2.5563135 +0.j,
 2.5563135 +0.j, 5.112627  +0.j, 2.7400243 +0.j, 2.7400243 +0.j,
 5.48004859+0.j, 2.95723668+0.j, 2.95723668+0.j, 5.91447337+0.j,
 2.59831713+0.j, 2.59831713+0.j, 5.19663426+0.j, 1.92295872+0.j,
 1.92295872+0.j, 3.84591745+0.j, 2.56105835+0.j, 2.56105835+0.j,
 5.12211671+0.j, 1.23654885+0.j, 1.23654885+0.j, 2.4730977 +0.j,
 2.27984204+0.j, 2.27984204+0.j, 4.55968409+0.j, 1.28670657+0.j,
 1.28670657+0.j, 2.57341315+0.j, 2.88933783+0.j, 2.88933783+0.j,
 5.77867567+0.j, 2.04369664+0.j, 2.04369664+0.j, 4.08739329+0.j,
 1.82932388+0.j, 1.82932388+0.j, 3.65864776+0.j, 1.52911122+0.j,
 1.52911122+0.j, 3.05822245+0.j, 2.54846738+0.j, 2.54846738+0.j,
 5.09693476+0.j, 1.91230066+0.j, 1.91230066+0.j, 3.82460133+0.j,
 2.1368679 +0.j, 2.1368679 +0.j, 4.2737358 +0.j, 1.0375796 +0.j,
 1.0375796 +0.j, 2.0751592 +0.j, 2.23527099+0.j, 2.23527099+0.j,
 4.47054199+0.j, 2.22419145+0.j, 2.22419145+0.j, 4.44838289+0.j,
 2.23386799+0.j, 2.23386799+0.j, 4.46773599+0.j, 2.88749616+0.j,
 2.88749616+0.j, 5.77499231+0.j, 2.3636406 +0.j, 2.3636406 +0.j,
 4.7272812 +0.j, 1.7190158 +0.j, 1.7190158 +0.j, 3.4380316 +0.j,
 1.87406391+0.j, 1.87406391+0.j, 3.74812782+0.j, 2.39526239+0.j,
 2.39526239+0.j, 4.79052478+0.j, 1.12045094+0.j, 1.12045094+0.j,
 2.24090189+0.j, 2.33353343+0.j, 2.33353343+0.j, 4.66706686+0.j,
 2.34127574+0.j, 2.34127574+0.j, 4.68255148+0.j, 1.42076512+0.j,
 1.42076512+0.j, 2.84153024+0.j, 1.2578526 +0.j, 1.2578526 +0.j,
 2.51570519+0.j, 1.6308567 +0.j, 1.6308567 +0.j, 3.2617134 +0.j])

עבור אובזרבל אנחנו בוחרים $Z_{24}Z_{25}$, כפי שנראה בפאנל התחתון של איור 5 ב-[3].

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

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

בחר צעדי Trotter

הניסוי המוצג באיור 4 של [3] משתמש ב-$k_j = [2, 3, 4]$ צעדי Trotter סימטריים מסדר $2$. אנחנו מתמקדים בתוצאות עבור זמן $t=3$, שבו ל-MPF ול-PF עם מספר גבוה יותר של צעדי Trotter (6 במקרה זה) יש אותה שגיאת Trotter. עם זאת, ערך הצפייה של ה-MPF מחושב ממעגלים המתאימים למספר הנמוך יותר של צעדי Trotter ולכן רדודים יותר. בפועל, גם אם ל-MPF ולמעגל צעדי ה-Trotter העמוקים יותר יש אותה שגיאת Trotter, אנחנו מצפים שערך הצפייה הניסיוני המחושב ממעגלי ה-MPF יהיה קרוב יותר לתאורטי, מכיוון שהוא כרוך בהרצת מעגלים רדודים יותר שפחות חשופים לרעש חומרה בהשוואה למעגל המתאים ל-PF של צעד Trotter גבוה יותר.

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

הגדר את ה-LSE

כאן אנחנו בוחנים את מקדמי ה-MPF הסטטיים עבור בעיה זו.

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(
    f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)
print("L1 norm:", np.linalg.norm(mpf_coeffs, ord=1))

The static coefficients associated with the ansatze are: [ 0.26666667 -2.31428571  3.04761905]
L1 norm: 5.628571428571431

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

[-0.24255546 -0.25744454  1.5       ]
L1 norm of the approximate coefficients: 2.0

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

# 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": 4,
    },
)

# Create exact time-evolution circuits
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)
]

# Create the time-evolution object
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)

lse = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs = setup_frobenius_problem(lse)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs.value
except Exception as error:
    print(error, "Calculation Failed for time", total_time)
print("")

בנה כל אחד ממעגלי ה-Trotter בפירוק ה-MPF שלנו

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

mpf_circuits = []
for k in mpf_trotter_steps:
    # Initial state preparation |1010..>
    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)

בנה מעגל Trotter עם שגיאת Trotter דומה ל-MPF

k = 6

# Initial state preparation |1010..>
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=k, order=order),
    time=total_time,
)

comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)

mpf_circuits.append(comp_circuit)

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

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

max_meas_err = 0.055
min_t2 = 30
max_twoq_err = 0.01

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    if target.qubit_properties[q].t2 is not None:
        t2 = target.qubit_properties[q].t2 * 1e6
    else:
        t2 = 0
    if meas_err > max_meas_err or t2 < min_t2:
        # print(q)
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
    twoq_gate_err = target[instruction_2q][q].error
    if twoq_gate_err > max_twoq_err:
        # print(q)
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates
    + ["measure"],  # or whatever new set of gates
    coupling_map=cust_cmap,
)

sorted_components = sorted(
    [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
    reverse=True,
)
print("size of largest component", len(sorted_components[0]))

size of largest component 73

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
    optimization_level=3, target=cust_target
)

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

qubits_layouts = [
    [
        idx
        for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
        if qb._register.name != "ancilla"
    ]
    for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
    transpiler = generate_preset_pass_manager(
        optimization_level=3, backend=backend, initial_layout=layout
    )
    transpiled_circuit = transpiler.run(circuit)
    transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

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

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
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 = ["mpf large"]

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

שלב 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.08034071), array(-0.00605026), array(-0.15345759), array(-0.18127293)]
[array(0.04482517), array(0.03438413), array(0.21540776), array(0.21520829)]

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.47510243192011536 +- 0.6613940032465087
Approximate static MPF expectation value:  -0.20914170384216998 +- 0.32341567460419135
Dynamic MPF expectation value:  -0.07994951978722761 +- 0.07423091963310202

sym = {2: "^", 3: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
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="6 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"Expectation values for (ZZ,{(L//2-1, L//2)}) 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()

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

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

חשוב גם לציין שבדוגמה שמשחזרת את הניסוי ב-[3], נקודת הזמן $t=3$ היא מעבר לגבול שבו מצופה שה-PF עם $k=2$ יתנהג היטב, שהוא $t/k>1$ כפי שנדון ב-מדריך זה.

אסמכתאות

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

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

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