אלכסון קוונטי של קריילוב מבוסס-דגימה (SKQD)

שיעור זה על אלכסון קוונטי של קריילוב מבוסס-דגימה (SKQD) משלב שיטות שהוסברו בשיעורים קודמים. הוא מורכב מדוגמה אחת המשתמשת במסגרת התבניות של Qiskit:

• שלב 1: מיפוי הבעיה ל-Circuit קוונטיים ואופרטורים

• שלב 2: אופטימיזציה עבור חומרה יעד

• שלב 3: הרצה באמצעות IBM Quantum primitives

• שלב 4: עיבוד התוצאות

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

0. מרחב קריילוב

נזכור שמרחב קריילוב $\mathcal{K}^r$ מסדר $r$ הוא המרחב הנפרש על-ידי וקטורים המתקבלים מהכפלת חזקות גבוהות יותר של מטריצה $A$, עד $r-1$, בוקטור ייחוס $\vert v \rangle$.

$$\mathcal{K}^r = \left\{ \vert v \rangle, A \vert v \rangle, A^2 \vert v \rangle, ..., A^{r-1} \vert v \rangle \right\}$$

אם המטריצה $A$ היא ההמילטוניאן $H$, המרחב המתאים נקרא מרחב קריילוב של חזקות $\mathcal{K}_P$. במקרה שבו $A$ הוא אופרטור האבולוציה בזמן הנוצר על-ידי ההמילטוניאן $U=e^{-iH(dt)}$, המרחב מכונה מרחב קריילוב האוניטרי $\mathcal{K}_U$. תת-מרחב קריילוב של חזקות אינו יכול להיווצר ישירות במחשב קוונטי, מאחר ש-$H$ אינו אופרטור אוניטרי. במקום זאת, נוכל להשתמש באופרטור האבולוציה בזמן $U = e^{-iH(dt)}$ שניתן להראות כי הוא מספק ערבויות התכנסות דומות למרחב קריילוב של חזקות. חזקות של $U$ הופכות אז לצעדי זמן שונים $U^k = e^{-iH(k dt)}$ כאשר $k = 0, 1, 2, ..., (r-1)$.

$$\mathcal{K}_U^r = \left\{ \vert \psi \rangle, U \vert \psi \rangle, U^2 \vert \psi \rangle, ..., U^{r-1} \vert \psi \rangle \right\}$$

1. מיפוי הבעיה ל-Circuit קוונטיים ואופרטורים

בשיעור זה אנחנו בוחנים את ההמילטוניאן לשרשרת ספין-$1/2$ של XX-Z אנטיפרומגנטית עם $L = 22$ אתרים ותנאי שפה מחזוריים:

$$H = \sum_{i, j}^{N} J_{xy} (X_{i} X_{j} + Y_{i} Y_{j}) + Z_{i} Z_{j}$$

# 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_sqd": "qiskit-addon-sqd", "qiskit_addon_utils": "qiskit-addon-utils", "qiskit_ibm_runtime": "qiskit-ibm-runtime"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

from qiskit.transpiler import CouplingMap
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

num_spins = 22
coupling_map = CouplingMap.from_ring(num_spins)
H_op = generate_xyz_hamiltonian(coupling_map, coupling_constants=(0.3, 0.3, 1.0))

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

1. בחירת מימד קריילוב ($r$) וצעד זמן ($dt$).

2. מצב ראשוני (ייחוס) (הוקטור $\vert v \rangle$ לעיל) עם חפיפה פולינומית עם מצב היעד (מצב היסוד), כאשר מצב היעד הוא דליל. דרישת החפיפה הפולינומית זהה לזו שבאלגוריתם אמדן הפאזה הקוונטית.

3. אופרטורי אבולוציה בזמן $U^{k}=e^{-iH(k * dt)}$ ($k = 0, 1, 2, ..., r-1$).

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

מימד קריילוב גדול יותר משפר את ההתכנסות של האנרגיה המוערכת. אנחנו קובעים את המימד ל-$5$ בשיעור זה כדי להמחיש את מגמת ההתכנסות.

הפניה [2] הראתה שצעד זמן קטן מספיק עבור KQD הוא $\pi / \vert \vert H \vert \vert$, ושעדיף להמעיט בהערכת ערך זה מאשר להגזים בו. מצד שני, בחירת $dt$ קטן מדי מובילה לקונדיציה גרועה יותר של תת-מרחב קריילוב, שכן וקטורי הבסיס של קריילוב שונים פחות מצעד זמן לצעד זמן. בנוסף, בעוד שבחירה זו של $dt$ מוכחת כמספקת להתכנסות של SKQD, בהקשר של דגימה זה, הבחירה האופטימלית של $dt$ בפועל היא נושא לחקר מתמשך. בשיעור זה אנחנו קובעים $dt = 0.15$.

מעבר למימד קריילוב ולצעד הזמן, עלינו לקבוע את מספר צעדי Trotter לאבולוציה בזמן. שימוש בצעדים מעטים מדי מוביל לשגיאות Trotterization גדולות יותר, בעוד שצעדים רבים מדי מובילים ל-Circuit עמוקים יותר. בשיעור זה אנחנו קובעים את מספר צעדי Trotter ל-$6$.

# Set parameters for quantum Krylov algorithm
krylov_dim = 5  # size of krylov subspace
dt = 0.15
num_trotter_steps = 6

לאחר מכן, עלינו לבחור מצב ייחוס $\vert \psi \rangle$ שיש לו חפיפה כלשהי עם מצב היסוד. עבור ההמילטוניאן הזה, אנחנו משתמשים במצב Neel עם "1"ים ו-"0"ים לסירוגין $\vert ...101...010...101 \rangle$ כמצב הייחוס שלנו.

# Prep `Neel` state as the reference state for evolution
from qiskit import QuantumCircuit

qc_state_prep = QuantumCircuit(num_spins)
for i in range(num_spins):
    if i % 2 == 0:
        qc_state_prep.x(i)

לבסוף, עלינו למפות את אופרטור האבולוציה בזמן ל-Circuit קוונטי. זה נעשה בשיעור 2, אך כאן נמנף שיטות ב-Qiskit, ספציפית שיטה בשם synthesis. קיימות שיטות שונות לסינתזה של אופרטורים מתמטיים ל-Circuit קוונטיים עם Gate קוונטיים. טכניקות רבות כאלה זמינות במודול הסינתזה של Qiskit. נשתמש בגישת LieTrotter לסינתזה [3] [4].

from qiskit.circuit import QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

evol_gate = PauliEvolutionGate(
    H_op, time=(dt / num_trotter_steps), synthesis=LieTrotter(reps=num_trotter_steps)
)  # `U` operator

qr = QuantumRegister(num_spins)
qc_evol = QuantumCircuit(qr)
qc_evol.append(evol_gate, qargs=qr)

circuits = []
for rep in range(krylov_dim):
    circ = qc_state_prep.copy()

    # Repeating the `U` operator to implement U^0, U^1, U^2, and so on, for power Krylov space
    for _ in range(rep):
        circ.compose(other=qc_evol, inplace=True)

    circ.measure_all()
    circuits.append(circ)

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

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

2. אופטימיזציה עבור חומרה יעד

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

import warnings

from qiskit import generate_preset_pass_manager
from qiskit_ibm_runtime import QiskitRuntimeService

warnings.filterwarnings("ignore")

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

כעת, אנחנו מבצעים transpilation של ה-Circuit ל-Backend היעד באמצעות מנהל מעברים מוגדר מראש.

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_circuits = pm.run(circuits=circuits)

3. הרצה על חומרה יעד

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

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
job = sampler.run(isa_circuits, shots=100_000)  # Takes approximately 2m 58s of QPU time

counts_all = [job.result()[k].data.meas.get_counts() for k in range(krylov_dim)]

4. עיבוד התוצאות

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

from collections import Counter

counts_cumulative = []
for i in range(krylov_dim):
    counter = Counter()
    for d in counts_all[: i + 1]:
        counter.update(d)

    counts = dict(counter)
    counts_cumulative.append(counts)

כדי להטיל ולאלכסן את ההמילטוניאן, אנחנו משתמשים ביכולות מ-qiskit-addon-sqd. ה-addon מציע פונקציונליות להטלת המילטוניאנים מבוססי מחרוזות פאולי על תת-מרחב ופותר ערכים עצמיים באמצעות SciPy.

from qiskit_addon_sqd.counts import counts_to_arrays
from qiskit_addon_sqd.qubit import solve_qubit

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

# Filters out bitstrings that do not have specified number (`num_ones`) of `1` bits.
def postselect_counts(counts, num_ones):
    filtered_counts = {}
    for bitstring, freq in counts.items():
        if bitstring.count("1") == num_ones:
            filtered_counts[bitstring] = freq

    return filtered_counts

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

import numpy as np

num_batches = 10
rand_seed = 0
scipy_kwargs = {"k": 2, "which": "SA"}

ground_state_energies = []
for idx, counts in enumerate(counts_cumulative):
    counts = postselect_counts(counts, num_ones=num_spins // 2)
    bitstring_matrix, probs = counts_to_arrays(counts=counts)

    eigenvals, eigenstates = solve_qubit(
        bitstring_matrix, H_op, verbose=False, **scipy_kwargs
    )
    gs_en = np.min(eigenvals)
    ground_state_energies.append(gs_en)

לאחר מכן, אנחנו מציגים את האנרגיה המחושבת כפונקציה של מימד קריילוב ומשווים עם האנרגיה המדויקת. האנרגיה המדויקת מחושבת בנפרד באמצעות שיטת כוח גס קלאסית. ניתן לראות שאנרגיית מצב היסוד המוערכת מתכנסת עם גידול מימד מרחב קריילוב. למרות שמימד קריילוב של $5$ הוא מגביל, התוצאות עדיין מראות התכנסות מרשימה, אשר צפויה להשתפר עם מימד קריילוב גדול יותר [1].

import matplotlib.pyplot as plt

exact_gs_en = -23.934184
plt.plot(
    range(1, krylov_dim + 1),
    ground_state_energies,
    color="blue",
    linestyle="-.",
    label="estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [exact_gs_en] * krylov_dim,
    color="red",
    linestyle="-",
    label="exact",
)
plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))
plt.legend()
plt.xlabel("Krylov space dimension")
plt.ylabel("Energy")
plt.ylim([-24, -22.50])
plt.title(
    "Estimating Ground state energy with Sample-based Krylov Quantum Diagonalization"
)
plt.show()

בדוק את הבנתך

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

מה אפשר לעשות כדי לשפר את ההתכנסות בגרף שלמעלה?

Answer:

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

מהם היתרונות העיקריים של SKQD על פני (א) SQD, ו-(ב) KQD?

Answer:

ייתכנו תשובות תקפות נוספות, אך תשובות מלאות צריכות לכלול את הבאים:

(א) ל-SKQD יש ערבויות התכנסות שאין ל-SQD. ב-SQD אתה צריך לנחש בצורה טובה מאוד את ה-ansatz שלך שיש לו חפיפה מצוינת עם התמיכה של מצב היסוד בבסיס החישובי, או שאתה צריך להכניס רכיב וריאציוני לחישוב כדי לדגום משפחה של ansaetze.

(ב) SKQD דורש הרבה פחות זמן QPU, מפני שהוא נמנע מהחישוב היקר של איברי מטריצה באמצעות מבחן הדמר.

5. סיכום

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

  • כדי לשמור על עומקי Circuit נמוכים, כדאי לעסוק בבעיות סריג המתאימות לחומרה טרום-סובלת-שגיאות.

• SKQD אינו סובל מבעיית מדידה קוונטית כמו ב-VQE. אין קבוצות של אופרטורי פאולי מתחלפים שיש להעריך.

• SKQD חסין לדגימות רועשות שכן ניתן להשתמש בשגרת בחירה לאחר מדידה ספציפית לבעיה (לדוגמה, לסנן מחרוזות סיביות שאינן עומדות בתבניות ספציפיות לבעיה) או לשאת בעלות אלכסון קלאסי (כלומר, לאלכסן בתת-מרחב גדול יותר) כדי להסיר למעשה את השפעת הרעש.

מקורות

[1] Jeffery Yu et al., "Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization" (2025). arxiv:quant-ph/2501.09702.

[2] Ethan N. Epperly, Lin Lin, and Yuji Nakatsukasa. "A theory of quantum subspace diagonalization". SIAM Journal on Matrix Analysis and Applications 43, 1263–1290 (2022).

[2] N. Hatano and M. Suzuki, "Finding Exponential Product Formulas of Higher Orders" (2005). arXiv:math-ph/0506007.

[4] D. Berry, G. Ahokas, R. Cleve and B. Sanders, "Efficient quantum algorithms for simulating sparse Hamiltonians" (2006). arXiv:quant-ph/0508139.