אומדן אנרגיית מצב-היסוד של שרשרת היזנברג עם VQE

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

תוצאות למידה

לאחר השלמת מדריך זה, תוכלו לצפות להבין את המידע הבא:

• כיצד לדגם שרשרת ספין של היזנברג כהמילטוניאן קוונטי באמצעות Qiskit
• כיצד להשתמש ב-optimizer SPSA לאומדן אנרגיית מצב-היסוד של מערכת קוונטית
• כיצד להריץ תהליכי עבודה וריאציוניים על חומרת קוונטום של IBM® באמצעות Qiskit Runtime primitives וסשנים

דרישות מוקדמות

מומלץ להתוודע לנושאים הבאים:

• יסודות מידע קוונטי
• מבוא לתבניות Qiskit
• עיצוב אלגוריתם וריאציוני

רקע

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

$H = \sum_{\langle i,j \rangle} \left( J_x X_i X_j + J_y Y_i Y_j + J_z Z_i Z_j \right) + \sum_{i} h_i Z_i,$

כאשר $X_i$, $Y_i$, ו-$Z_i$ הם אופרטורי פאולי הפועלים על אתר $i$, הסכום $\langle i,j \rangle$ רץ על פני זוגות שכנים קרובים, $J_x = J_y = J_z = 0.5$ הם קבועי הצימוד החילופין (איזוטרופי במדריך זה), ו-$h_i$ מייצג שדה מגנטי חיצוני תלוי-אתר. במדריך זה, ערכי השדה המגנטי נדגמים אקראית מהקטע $[-1, 1]$. שימו לב שבמימוש שלהלן, קבוצת זוגות "שכנים קרובים" נקבעת על פי הצימוד הטבעי של ה-backend לחומרה בין $N$ ה-qubits הראשונים, שאינם בהכרח מהווים שרשרת לינארית קפדנית בהתאם לטופולוגיית המכשיר.

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

ה-Variational Quantum Eigensolver (VQE) הוא אלגוריתם היברידי קוונטי-קלאסי המיועד לאומדן אנרגיית מצב-היסוד של המילטוניאן. הוא פועל על ידי הכנת מצב קוונטי פרמטרי $|\psi(\theta)\rangle$ (הנקרא ansatz) על מחשב קוונטי ומדידת ערך הציפייה $\langle \psi(\theta) | H | \psi(\theta) \rangle$. optimizer קלאסי מכוונן אז איטרטיבית את הפרמטרים $\theta$ כדי למזעֵר אנרגיה זו, מנצל את עקרון הוריאציה המבטיח שהאנרגיה הנמדדת היא תמיד חסם עליון לאנרגיית מצב-היסוד האמיתית.

במדריך זה, אנו משתמשים ב-ansatz ‏efficient_su2 מספריית המעגלים של Qiskit, אשר בונה שכבות של סיבובי qubit יחיד ושערי שזירה. האופטימיזציה מתבצעת באמצעות אלגוריתם Simultaneous Perturbation Stochastic Approximation (SPSA), המתאים במיוחד לחומרה קוונטית רועשת מפני שהוא מאמד גרדיאנטים תוך שימוש בשתי הערכות פונקציה בלבד לכל איטרציה ללא תלות במספר הפרמטרים.

דרישות

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

• Qiskit SDK v2.0 או מאוחר יותר, עם תמיכה בהדמייה
• Qiskit Runtime v0.44 או מאוחר יותר (pip install qiskit-ibm-runtime)

הגדרה

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-ibm-runtime

import numpy as np
import matplotlib.pyplot as plt
from typing import Sequence

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import BaseEstimatorV2
from qiskit.circuit.library import XGate
from qiskit.circuit.library import efficient_su2
from qiskit.transpiler import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes.scheduling import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
)
from qiskit_ibm_runtime import QiskitRuntimeService, Session, EstimatorV2

def visualize_results(results):
    plt.plot(results["cost_history"], lw=2)
    plt.xlabel("Number of function evaluations")
    plt.ylabel("Energy")
    plt.show()

דוגמה בקנה מידה קטן

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

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

• קלט: מספר הספינים
• פלט: Ansatz והמילטוניאן המדגימים את שרשרת היזנברג

נבנה ansatz והמילטוניאן המדגימים שרשרת היזנברג בעלת 10 ספינים. בשלב זה, נבנה המילטוניאן של היזנברג בעל 10 ספינים על מפת הצימוד של ה-backend הפחות עמוס ונכין את ה-ansatz ‏efficient_su2.

num_spins = 10
ansatz = efficient_su2(num_qubits=num_spins, reps=2)

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=num_spins, simulator=False
)

coupling = backend.target.build_coupling_map()
reduced_coupling = coupling.reduce(list(range(num_spins)))

edge_list = reduced_coupling.graph.edge_list()
ham_list = []

for edge in edge_list:
    ham_list.append(("ZZ", edge, 0.5))
    ham_list.append(("YY", edge, 0.5))
    ham_list.append(("XX", edge, 0.5))

for qubit in reduced_coupling.physical_qubits:
    ham_list.append(("Z", [qubit], np.random.random() * 2 - 1))

hamiltonian = SparsePauliOp.from_sparse_list(ham_list, num_qubits=num_spins)

ansatz.draw("mpl", style="iqp")

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

• קלט: מעגל מופשט, observable
• פלט: מעגל יעד ו-observable, מותאמים עבור ה-QPU שנבחר

השתמשו בפונקציה generate_preset_pass_manager מ-Qiskit כדי ליצור אוטומטית שגרת אופטימיזציה עבור המעגל שלנו ביחס ל-QPU שנבחר. אנחנו בוחרים optimization_level=3, המספק את הרמה הגבוהה ביותר של אופטימיזציה מבין מנהלי ההעברות הקבועים מראש. אנחנו גם כוללים העברות תזמון ALAPScheduleAnalysis ו-PadDynamicalDecoupling כדי לדכא שגיאות דה-קוהרנטיות.

target = backend.target
pm = generate_preset_pass_manager(optimization_level=3, target=target)
pm.scheduling = PassManager(
    [
        ALAPScheduleAnalysis(durations=target.durations()),
        PadDynamicalDecoupling(
            durations=target.durations(),
            dd_sequence=[XGate(), XGate()],
            pulse_alignment=target.pulse_alignment,
        ),
    ]
)
isa_ansatz = pm.run(ansatz)
isa_observable = hamiltonian.apply_layout(isa_ansatz.layout)
isa_ansatz.draw("mpl", scale=0.6, style="iqp", fold=-1, idle_wires=False)

שלב 3: הרצה באמצעות primitives של Qiskit

• קלט: מעגל יעד ו-observable
• פלט: תוצאות האופטימיזציה

מזעור אנרגיית מצב-היסוד המוערכת של המערכת על ידי אופטימיזציה של פרמטרי המעגל. השתמשו ב-primitive Estimator מ-Qiskit Runtime כדי להעריך את פונקציית העלות במהלך האופטימיזציה.

מאחר שאופטמנו את המעגל עבור ה-backend בשלב 2, ניתן להימנע מטרנספילציה בשרת Runtime על ידי הגדרת skip_transpilation=True והעברת המעגל המאופטם. להדגמה זו, נריץ על QPU באמצעות primitives של qiskit-ibm-runtime. כדי להריץ עם primitives מבוססי statevector של qiskit, החליפו את בלוק הקוד המשתמש ב-primitives של Qiskit Runtime בבלוק המוערת.

במדריך זה אנו משתמשים ב-Simultaneous Perturbation Stochastic Approximation (SPSA), שהוא optimizer מבוסס גרדיאנט. להלן הצגה קצרה שלו, וקוד לממש SPSA עם Qiskit v2.0.

הצגת SPSA

Simultaneous Perturbation Stochastic Approximation (SPSA) [1] הוא אלגוריתם אופטימיזציה המקרב את וקטור הגרדיאנט כולו בשימוש בשתי קריאות פונקציה בלבד בכל איטרציה. תהי $f:\mathbb{R}^p\rightarrow \mathbb{R}$ פונקציית העלות עם $p$ פרמטרים לאופטימיזציה, ו-$x_i\in \mathbb{R}^p$ וקטור הפרמטרים בשלב ה-$i^{th}$ של האיטרציה. לחישוב הגרדיאנט, נוצר וקטור אקראי $\Delta_i$ בגודל $p$, כאשר כל איבר $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$, נדגם אחיד מ-$\{-1, 1\}$. לאחר מכן, כל איבר של הוקטור האקראי $\Delta_i$ מוכפל בערך קטן $c_i$ כדי ליצור הפרעה אקראית. הגרדיאנט מוערך אז כ-

$[\nabla f(x_i)]_j \approx \frac{f(x_i + c_i \Delta_i) - f(x_i - c_i \Delta_i)}{2c_i\Delta_{ij}}.$

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

במדריך זה, ה-hyperparameters עבור האיטרציה ה-$i^{th}$, $a_i$ ו-$c_i$, מחושבים כ-

$a_i = \frac{a}{(A + i + 1)^\alpha} \quad \text{and} \quad c_i = \frac{c}{(i+1)^\gamma},$

כאשר הערכים הקבועים הם $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$, ו-$\gamma = 0.4$. ערכים אלה נבחרו מ-[2]. כיוון נאות של hyperparameters הכרחי להפקת ביצועים טובים מ-SPSA.

def spsa(
    fun, x0, args=(), A=30, alpha=0.9, a=0.3, c=0.1, gamma=0.4, maxiter=100
):
    nparams = len(x0)
    x = np.copy(x0)

    for i in range(maxiter):
        a_i = a / (A + i + 1) ** alpha
        c_i = c / (i + 1) ** gamma
        delta_i = np.random.choice([-1, 1], nparams)

        # two hardware calls
        eval_1 = fun(x + c_i * delta_i, *args)
        eval_2 = fun(x - c_i * delta_i, *args)

        # compute the gradient and update the parameters
        grad = (eval_1 - eval_2) / (2 * c_i) * np.reciprocal(delta_i)
        x = x - a_i * grad

    return x

def cost_func(
    params: Sequence,
    ansatz: QuantumCircuit,
    hamiltonian: SparsePauliOp,
    estimator: BaseEstimatorV2,
    cost_history_dict: dict,
) -> float:
    """Ground state energy evaluation."""
    energy = (
        estimator.run([(ansatz, hamiltonian, [params])]).result()[0].data.evs
    )

    cost_history_dict["iters"] += 1
    cost_history_dict["prev_vector"] = list(params)
    cost_history_dict["cost_history"].append(float(energy[0]))

    print(
        f"Fx Iters. done: {cost_history_dict['iters']} [Current cost: {round(energy[0], 5)}]",
        end="\r",
    )

    return energy

def solve(x0, isa_ansatz, isa_observable, maxiter=150):
    cost_history_dict = {
        "prev_vector": None,
        "iters": 0,
        "cost_history": [],
        "y_min": None,
    }

    # Evaluate the problem using a QPU via Qiskit IBM Runtime
    with Session(backend=backend) as session:
        estimator = EstimatorV2(mode=session)
        estimator.skip_transpilation = True
        estimator.options.environment.job_tags = ["TUT_HSVQE"]
        x_opt = spsa(
            cost_func,
            x0=x0,
            args=(isa_ansatz, isa_observable, estimator, cost_history_dict),
            maxiter=maxiter,
        )

        y_min = cost_func(
            x_opt, isa_ansatz, isa_observable, estimator, cost_history_dict
        )

    return y_min, cost_history_dict

np.random.seed(42)
num_params = ansatz.num_parameters
params = 2 * np.pi * np.random.random(num_params)

כאן אנו מגדירים את maxiter = 50. שימו לב שמאחר שכל איטרציה מצריכה שתי קריאות לפונקציה לחישוב הגרדיאנט, סך מספר קריאות הפונקציה יהיה $2 \times \text{maxiter}$. ניתן להגדיל את maxiter לכל ערך גבוה יותר לאומדן אנרגיה טוב יותר.

maxiter = 50
spsa_min, spsa_history = solve(
    params, isa_ansatz, isa_observable, maxiter=maxiter
)

Fx Iters. done: 101 [Current cost: -3.03843]

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

• קלט: אומדני אנרגיית מצב-היסוד במהלך האופטימיזציה
• פלט: אנרגיית מצב-יסוד מוערכת

print(f"Estimated ground state energy: {spsa_min}")

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

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

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

אתגר

עכשיו שיש לכם מימוש VQE עובד עבור שרשרת היזנברג, נסו את הדברים הבאים:

1. ניסוי עם עומק ה-ansatz: שנו את פרמטר reps ב-efficient_su2 (לדוגמה, נסו reps=1 ו-reps=3). כיצד עומק ה-ansatz משפיע על אנרגיית מצב-היסוד המוערכת ומהירות ההתכנסות? באיזו נקודה אתם מבחינים בתשואות פוחתות או חוסר יציבות?
2. כוונון hyperparameters של SPSA: כוונו את פרמטרי לוח הזמנים של קצב הלמידה (a, c, alpha, gamma, A) והתבוננו כיצד הם משפיעים על ההתכנסות. האם תוכלו למצוא הגדרות שמתכנסות מהר יותר מהברירות המחדל שבשימוש כאן?
3. השוואת טופולוגיות צימוד: במקום להשתמש במפת הצימוד הטבעית של ה-backend, נסו לבנות שרשרת לינארית פשוטה של שכנים קרובים והשוו את התוצאות. כיצד הקישוריות של החומרה הפיזית משפיעה על עומק המעגל המוטרנספל ואומדן האנרגיה הסופי?

מקורות

[1] Spall, J. C. (2002). Implementation of the simultaneous perturbation algorithm for stochastic optimization. IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817-823.

[2] Sahin, M. Emre, et al. (2025). Qiskit Machine Learning: an open-source library for quantum machine learning tasks at scale on quantum hardware and classical simulators. arXiv:2505.17756.

צעדים הבאים

המלצות

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

• נסו Sample-based Quantum Diagonalization (SQD): כפי שהודגם במדריך זה, VQE מתמודד עם אתגרים בקנה מידה גדול בשל מישורים בורים ועומס מדידות גבוה. IBM פיתחה Sample-based Quantum Diagonalization (SQD) כחלופה ניתנת להרחבה יותר. בשונה מ-VQE, SQD נמנעת כלל מאופטימיזציה וריאציונית; במקום זאת, מחשב קוונטי מייצר דגימות ומחשב קלאסי מקרין את ההמילטוניאן על תת-מרחב שנפרש על ידי הדגימות הללו ומבצע עליו דיאגונליזציה. זה מספק חסם עליון לאנרגיית מצב-היסוד עם מדידות פחות משמעותית וללא רגישות למישורים בורים. עקבו אחר מדריך SQD כדי לראות גישה זו בפעולה.
• חקרו את קורס אלגוריתמי דיאגונליזציה קוונטית: העמיקו את ההבנה שלכם הן ב-VQE והן ב-SQD, כולל פשרותיהם, בקורס אלגוריתמי דיאגונליזציה קוונטית ב-IBM Quantum Learning.