מקרים ורחבות

פרק זה יסקור מספר אלגוריתמים קוונטיים וריאציוניים, כולל

• Variational Quantum Eigensolver (VQE)
• Subspace Search VQE (SSVQE)
• Variational Quantum Deflation (VQD)
• Quantum Sampling Regression (QSR)

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

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

Variational Quantum Eigensolver (VQE)

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

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

המבנה של VQE הוא פשוט:

• הכן אופרטורי ייחוס $U_R$
  • אנחנו מתחילים מהמצב $|0\rangle$ ועוברים למצב הייחוס $|\rho\rangle$
• הפעל את הצורה הוריאציונית $U_V(\vec\theta_{i,j})$ כדי ליצור ansatz $U_A(\vec\theta_{i,j})$
  • אנחנו עוברים מהמצב $|\rho\rangle$ אל $U_V(\vec\theta_{i,j})|\rho\rangle = |\psi(\vec\theta_{i,j})\rangle$
• בצע אתחול בנקודה $i=0$ אם יש לנו בעיה דומה (בדרך כלל נמצאת באמצעות סימולציה קלאסית או דגימה)
  • כל מייעל יאותחל באופן שונה, ויוצר קבוצה ראשונית של וקטורי פרמטרים $\Theta_0 := \\{ {\vec\theta_{0,j} | j \in \mathcal{J}_\text{opt}^0} \\}$ (למשל, מנקודת התחלה $\vec\theta_0$).
• חשב את פונקציית העלות $C(\vec\theta_{i,j}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ עבור כל המצבים המוכנים במחשב קוונטי.
• השתמש במייעל קלאסי לבחירת קבוצת הפרמטרים הבאה $\Theta_{i+1}$.
• חזור על התהליך עד להשגת התכנסות.

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

הנה הדוגמה עבור האוברוול הבא:

$$\hat{O}_1 = 2 II - 2 XX + 3 YY - 3 ZZ,$$

מבנה תיאורטי

יישום

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

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.library import TwoLocal
import numpy as np

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
observable = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
    2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)

ansatz = reference_circuit.compose(variational_form)

ansatz.decompose().draw("mpl")

def cost_func_vqe(parameters, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """

    estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])
    estimator_result = estimator_job.result()[0]

    cost = estimator_result.data.evs[0]
    return cost

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()

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

# SciPy minimizer routine
from scipy.optimize import minimize

x0 = np.ones(8)

result = minimize(
    cost_func_vqe, x0, args=(ansatz, observable, estimator), method="COBYLA"
)

result

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -5.999999982445723
       x: [ 1.741e+00  9.606e-01  1.571e+00  2.115e-05  1.899e+00
            1.243e+00  6.063e-01  6.063e-01]
    nfev: 136
   maxcv: 0.0

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

נבחר את ה-Backend הפחות עמוס, ונייבא את הרכיבים הנדרשים מ-qiskit_ibm_runtime.

from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import Session, EstimatorOptions
from qiskit_ibm_runtime import QiskitRuntimeService

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

<IBMBackend('ibm_brisbane')>

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

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = observable.apply_layout(layout=isa_ansatz.layout)

שלב 3: ביצוע באמצעות Qiskit Runtime primitives

אנחנו מוכנים כעת להריץ את החישוב שלנו על חומרת IBM Quantum®. מכיוון שמיזור פונקציית העלות הוא תהליך איטרטיבי מאוד, נפתח Session. כך נצטרך לחכות בתור פעם אחת בלבד. לאחר שהמשימה תחל לרוץ, כל איטרציה עם עדכוני פרמטרים תרוץ מיידית.

x0 = np.ones(8)

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    result = minimize(
        cost_func_vqe,
        x0,
        args=(isa_ansatz, isa_observable, estimator),
        method="COBYLA",
        options={"maxiter": 200, "disp": True},
    )
session.close()
print(result)

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

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

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

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

Subspace Search VQE (SSVQE)

SSVQE הוא גרסה של VQE המאפשרת קבלת $k$ ערכים עצמיים ראשונים של אוברוול $\hat{H}$ עם ערכים עצמיים $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, כאשר $N\geq k$. ללא אובדן כלליות, אנחנו מניחים ש-$\lambda_0<\lambda_1<...<\lambda_{N-1}$. SSVQE מציגה רעיון חדש על ידי הוספת משקלים כדי לעזור לתת עדיפות לאופטימיזציה של המונה עם המשקל הגבוה ביותר.

כדי לממש אלגוריתם זה, אנחנו זקוקים ל-$k$ מצבי ייחוס אורתוגונליים זה לזה $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, כלומר $\langle \rho_j | \rho_l \rangle = \delta_{jl}$ עבור $j,l<k$. מצבים אלה ניתן לבנות באמצעות אופרטורי Pauli. פונקציית העלות של אלגוריתם זה היא אז:

$$\begin{aligned} C(\vec{\theta}) & := \sum_{j=0}^{k-1} w_j \langle \rho_j | U_{V}^{\dagger}(\vec{\theta})\hat{H} U_{V}(\vec{\theta})|\rho_j \rangle \\[1mm] & := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle \\[1mm] \end{aligned}$$

כאשר $w_j$ הוא מספר חיובי שרירותי כך שאם $j<l<k$ אז $w_j>w_l$, ו-$U_V(\vec{\theta})$ הוא הצורה הוריאציונית המוגדרת על ידי המשתמש.

אלגוריתם SSVQE מסתמך על העובדה שמצבים עצמיים התואמים לערכים עצמיים שונים הם אורתוגונליים זה לזה. באופן ספציפי, המכפלה הפנימית של $U_V(\vec{\theta})|\rho_j\rangle$ ו-$U_V(\vec{\theta})|\rho_l\rangle$ ניתנת לביטוי כ:

$$\begin{aligned} \langle \rho_j | U_{V}^{\dagger}(\vec{\theta})U_{V}(\vec{\theta})|\rho_l \rangle & = \langle \rho_j | I |\rho_l \rangle \\[1mm] & = \langle \rho_j | \rho_l \rangle \\[1mm] & = \delta_{jl} \end{aligned}$$

השוויון הראשון מתקיים מכיוון ש-$U_{V}(\vec{\theta})$ הוא אופרטור קוונטי ולכן הוא יוניטרי. השוויון האחרון מתקיים בשל האורתוגונליות של מצבי הייחוס $|\rho_j\rangle$. העובדה שאורתוגונליות נשמרת דרך טרנספורמציות יוניטריות קשורה עמוקות לעיקרון שימור המידע, כפי שהוא מבוטא במדע המידע הקוונטי. לפי תפיסה זו, טרנספורמציות לא-יוניטריות מייצגות תהליכים בהם מידע אובד או מוזרק.

משקלים $w_j$ עוזרים להבטיח שכל המצבים הם מצבים עצמיים. אם המשקלים שונים מספיק, למונה עם המשקל הגבוה ביותר (כלומר, $w_0$) תינתן עדיפות במהלך האופטימיזציה על פני האחרים. כתוצאה מכך, המצב המתקבל $U_{V}(\vec{\theta})|\rho_0 \rangle$ יהפוך למצב העצמי התואם ל-$\lambda_0$. מכיוון ש-$\{ U_{V}(\vec{\theta})|\rho_j\rangle \}_{j=0}^{k-1}$ אורתוגונליים זה לזה, שאר המצבים יהיו אורתוגונליים לו ולכן יהיו כלולים במרחב-החלק התואם לערכים העצמיים $\{\lambda_1,...,\lambda_{N-1}\}$.

אם נחיל את אותו הטיעון על שאר המונות, העדיפות הבאה תהיה למונה עם המשקל $w_1$, כך ש-$U_{V}(\vec{\theta})|\rho_1 \rangle$ יהיה המצב העצמי התואם ל-$\lambda_1$, והמונות האחרים יהיו כלולים במרחב העצמי של $\{\lambda_2,...,\lambda_{N-1}\}$.

על ידי הסקה אינדוקטיבית, אנחנו מסיקים ש-$U_{V}(\vec{\theta})|\rho_j \rangle$ יהיה מצב עצמי קירובי של $\lambda_j$ עבור $0\leq j < k.$

מבנה תיאורטי

ניתן לסכם את SSVQE כך:

• הכן מספר מצבי ייחוס על ידי הפעלת יוניטרי U_R על k מצבי בסיס חישוביים שונים
  • אלגוריתם זה דורש שימוש ב-$k$ מצבי ייחוס אורתוגונליים זה לזה $\{ |\rho_j\rangle \}_{j=0}^{k-1}$, כך ש-$\langle \rho_j | \rho_l \rangle = \delta_{jl}$ עבור $j,l<k$.
• הפעל את הצורה הוריאציונית $U_V(\vec\theta_{i,j})$ על כל מצב ייחוס, ויוצר את ה-ansatz הבא $U_A(\vec\theta_{i,j})$.
• בצע אתחול בנקודה $i=0$ אם בעיה דומה זמינה (בדרך כלל נמצאת באמצעות סימולציה קלאסית או דגימה).
• חשב את פונקציית העלות $C(\vec\theta_{i,j}) := \sum_{j=0}^{k-1} w_j \langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ עבור כל המצבים המוכנים במחשב קוונטי.
  • ניתן לפצל זאת לחישוב ערך הציפייה עבור אוברוול $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ ולהכפיל את התוצאה ב-$w_j$.
  • לאחר מכן, פונקציית העלות מחזירה את סכום כל ערכי הציפייה המשוקללים.
• השתמש במייעל קלאסי לקביעת קבוצת הפרמטרים הבאה $\Theta_{i+1}$.
• חזור על הצעדים הנ"ל עד להשגת התכנסות.

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

import numpy as np

def cost_func_ssvqe(
    params, initialized_anastz_list, weights, ansatz, hamiltonian, estimator
):
    # """Return estimate of energy from estimator

    # Parameters:
    #     params (ndarray): Array of ansatz parameters
    #     initialized_anastz_list (list QuantumCircuit): Array of initialised ansatz with reference
    #     weights (list): List of weights
    #     ansatz (QuantumCircuit): Parameterized ansatz circuit
    #     hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
    #     estimator (Estimator): Estimator primitive instance

    # Returns:
    #     float: Weighted energy estimate
    # """

    energies = []

    # Define SSVQE

    weighted_energy_sum = np.dot(energies, weights)
    return weighted_energy_sum

Variational Quantum Deflation (VQD)

VQD היא שיטה איטרטיבית שמרחיבה את VQE כדי לקבל את $k$ ערכי העצמי הראשונים של אופרטור $\hat{H}$ עם ערכי עצמי $\{\lambda_0, \lambda_1,...,\lambda_{N-1}\}$, כאשר $N\geq k$, ולא רק את הראשון. לשאר החלק הזה נניח, ללא אבדן כלליות, ש-$\lambda_0\leq\lambda_1\leq...\leq\lambda_{N-1}$. VQD מכניסה את הרעיון של עלות עונשין כדי לכוון את תהליך האופטימיזציה.

VQD מכניסה איבר עונשין, המסומן ב-$\beta$, כדי לאזן את תרומת כל איבר חפיפה לעלות. איבר עונשין זה נועד להעניש את תהליך האופטימיזציה אם לא מושגת אורתוגונליות. אנחנו מטילים אילוץ זה מכיוון שמצבי העצמי של אופרטור אבחוני, או אופרטור הרמיטי, המתאימים לערכי עצמי שונים הם תמיד אורתוגונליים זה לזה, או שניתן לעשותם כאלה במקרה של ניוון או ערכי עצמי חוזרים. לכן, על ידי אכיפת אורתוגונליות עם מצב העצמי המתאים ל-$\lambda_0$, אנחנו בעצם מבצעים אופטימיזציה על תת-המרחב המתאים לשאר ערכי העצמי $\{\lambda_1, \lambda_2,..., \lambda_{N-1}\}$. כאן, $\lambda_1$ הוא ערך העצמי הנמוך ביותר מבין שאר ערכי העצמי ולכן הפתרון האופטימלי של הבעיה החדשה ניתן להשגה באמצעות משפט הווריאציה.

הרעיון הכללי מאחורי VQD הוא להשתמש ב-VQE כרגיל כדי לקבל את ערך העצמי הנמוך ביותר $\lambda_0 := C_0(\vec\theta^0) \equiv C_\text{VQE}(\vec\theta^0)$ יחד עם מצב העצמי (המקורב) המתאים $|\psi(\vec{\theta^0})\rangle$ עבור וקטור פרמטרים אופטימלי כלשהו $\vec{\theta^0}$. לאחר מכן, כדי לקבל את ערך העצמי הבא $\lambda_1 > \lambda_0$, במקום למזמן את פונקציית העלות $C_0(\vec{\theta}) := \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$, אנחנו ממטבים:

$$C_1(\vec{\theta}) := C_0(\vec{\theta})+ \beta_0 |\langle \psi(\vec{\theta})| \psi(\vec{\theta^0})\rangle |^2$$

הערך החיובי $\beta_0$ צריך באופן אידיאלי להיות גדול מ-$\lambda_1-\lambda_0$.

זה מכניס פונקציית עלות חדשה שניתן להתייחס אליה כבעיה עם אילוצים, כאשר אנחנו ממזמנים את $C_\text{VQE}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ בכפוף לאילוץ שהמצב חייב להיות אורתוגונלי ל-$|\psi(\vec{\theta^0})\rangle$ שהתקבל קודם, כאשר $\beta_0$ פועל כאיבר עונשין אם האילוץ לא מתקיים.

לחלופין, ניתן לפרש את הבעיה החדשה הזו כהרצת VQE על האופרטור החדש:

$$\hat{H_1} := \hat{H} + \beta_0 |\psi(\vec{\theta^0})\rangle \langle \psi(\vec{\theta^0})| \quad \Rightarrow \quad C_1(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H_1} | \psi(\vec{\theta})\rangle,$$

בהנחה שהפתרון לבעיה החדשה הוא $|\psi(\vec{\theta^1})\rangle$, הערך הצפוי של $\hat{H}$ (לא $\hat{H_1}$) צריך להיות $\langle \psi(\vec{\theta^1}) | \hat{H} | \psi(\vec{\theta^1})\rangle = \lambda_1$. כדי לקבל את ערך העצמי השלישי $\lambda_2$, פונקציית העלות לאופטימיזציה היא:

$$C_2(\vec{\theta}) := C_1(\vec{\theta}) + \beta_1 |\langle \psi(\vec{\theta})| \psi(\vec{\theta^1})\rangle |^2$$

כאשר $\beta_1$ הוא קבוע חיובי גדול מספיק כדי לאכוף אורתוגונליות של מצב הפתרון הן ל-$|\psi(\vec{\theta^0})\rangle$ והן ל-$|\psi(\vec{\theta^1})\rangle$. זה מעניש מצבים במרחב החיפוש שאינם עומדים בדרישה זו, ובכך מגביל למעשה את מרחב החיפוש. לכן, הפתרון האופטימלי של הבעיה החדשה צריך להיות מצב העצמי המתאים ל-$\lambda_2$.

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

$$\hat{H_2} := \hat{H_1} + \beta_1 |\psi(\vec{\theta^1})\rangle \langle \psi(\vec{\theta^1})| \quad \Rightarrow \quad C_2(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H_2} | \psi(\vec{\theta})\rangle.$$

אם הפתרון לבעיה החדשה הזו הוא $|\psi(\vec{\theta^2})\rangle$, הערך הצפוי של $\hat{H}$ (לא $\hat{H_2}$) צריך להיות $\langle \psi(\vec{\theta^2}) | \hat{H} | \psi(\vec{\theta^2})\rangle = \lambda_2$. באנלוגיה, כדי לקבל את ערך העצמי ה-$k$, $\lambda_{k-1}$, תמזמן את פונקציית העלות:

$$C_{k-1}(\vec{\theta}) := C_{k-2}(\vec{\theta}) + \beta_{k-2} |\langle \psi(\vec{\theta})| \psi(\vec{\theta^{k-2}})\rangle |^2,$$

זכור שהגדרנו את $\vec{\theta^j}$ כך ש-$\langle \psi(\vec{\theta^j}) | \hat{H} | \psi(\vec{\theta^j})\rangle = \lambda_j, \forall j<k$. בעיה זו שקולה למזמון $C(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle$ אך עם האילוץ שהמצב חייב להיות אורתוגונלי ל-$|\psi(\vec{\theta^j})\rangle ; \forall j \in {0, \cdots, k-1}$, ובכך מגביל את מרחב החיפוש לתת-המרחב המתאים לערכי העצמי $\{\lambda_{k-1},\cdots,\lambda_{N-1}\}$.

בעיה זו שקולה ל-VQE עם האופרטור:

$$\hat{H}_{k-1} := \hat{H}_{k-2} + \beta_{k-2} |\psi(\vec{\theta^{k-2}})\rangle \langle \psi(\vec{\theta^{k-2}})| \quad \Rightarrow \quad C_{k-1}(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H}_{k-1} | \psi(\vec{\theta})\rangle,$$

כפי שאפשר לראות מהתהליך, כדי לקבל את ערך העצמי ה-$k$, אתה צריך את מצבי העצמי (המקורבים) של $k-1$ ערכי העצמי הקודמים, ולכן תצטרך להריץ VQE סך הכל $k$ פעמים. לכן, פונקציית העלות של VQD היא כדלקמן:

$$C_k(\vec{\theta}) = \langle \psi(\vec{\theta}) | \hat{H} | \psi(\vec{\theta})\rangle + \sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$$

מבנה תיאורטי

ניתן לסכם את המבנה של VQD כדלקמן:

• הכן אופרטור ייחוס $U_R$
• הפעל את הצורה הווריאציונית $U_V(\vec\theta_{i,j})$ על מצב הייחוס, ויצור את ה-ansatz הבא $U_A(\vec\theta_{i,j})$
• אתחל ב-$i=0$ אם יש לנו בעיה דומה (נמצאת בדרך כלל על ידי סימולציה קלאסית או דגימה).
• חשב את פונקציית העלות $C_k(\vec{\theta})$, הכוללת חישוב $k$ מצבים מעוררים ומערך של $\beta$-ים המגדירים את עונשין החפיפה לכל איבר חפיפה.
  • חשב את ערך הציפייה של אופרטור $\langle \psi_{j}(\vec{\theta}) | \hat{H} | \psi_{j}(\vec{\theta}) \rangle$ עבור כל $k$
  • חשב את העונשין $\sum_{j=0}^{k-1}\beta_j |\langle \psi(\vec{\theta})| \psi(\vec{\theta^j})\rangle |^2$.
  • פונקציית העלות צריכה להחזיר את סכום שני האיברים הללו
• השתמש ב-optimizer קלאסי לבחירת הקבוצה הבאה של פרמטרים $\Theta_{i+1}$.
• חזור על התהליך עד להתכנסות.

מימוש

עבור המימוש הזה, ניצור פונקציה לעונשין חפיפה. עונשין זה ישמש בפונקציית העלות בכל איטרציה. תהליך זה יחזור על עצמו עבור כל מצב מעורר

from qiskit.circuit.library import TwoLocal

ansatz = TwoLocal(2, rotation_blocks=["ry", "rz"], entanglement_blocks="cz", reps=1)

ansatz.decompose().draw("mpl")

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

import numpy as np

def calculate_overlaps(ansatz, prev_circuits, parameters, sampler):
    def create_fidelity_circuit(circuit_1, circuit_2):
        """
        Constructs the list of fidelity circuits to be evaluated.
        These circuits represent the state overlap between pairs of input circuits,
        and their construction depends on the fidelity method implementations.
        """

        if len(circuit_1.clbits) > 0:
            circuit_1.remove_final_measurements()
        if len(circuit_2.clbits) > 0:
            circuit_2.remove_final_measurements()

        circuit = circuit_1.compose(circuit_2.inverse())
        circuit.measure_all()
        return circuit

    overlaps = []

    for prev_circuit in prev_circuits:
        fidelity_circuit = create_fidelity_circuit(ansatz, prev_circuit)
        sampler_job = sampler.run([(fidelity_circuit, parameters)])
        meas_data = sampler_job.result()[0].data.meas

        counts_0 = meas_data.get_int_counts().get(0, 0)
        shots = meas_data.num_shots
        overlap = counts_0 / shots
        overlaps.append(overlap)

    return np.array(overlaps)

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

def cost_func_vqd(
    parameters, ansatz, prev_states, step, betas, estimator, sampler, hamiltonian
):
    estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])

    total_cost = 0

    if step > 1:
        overlaps = calculate_overlaps(ansatz, prev_states, parameters, sampler)
        total_cost = np.sum(
            [np.real(betas[state] * overlap) for state, overlap in enumerate(overlaps)]
        )

    estimator_result = estimator_job.result()[0]

    value = estimator_result.data.evs[0] + total_cost

    return value

שים לב במיוחד שפונקציית העלות לעיל מתייחסת לפונקציה calculate_overlaps, שבפועל יוצרת Circuit קוונטי חדש. אם אנחנו רוצים להריץ על חומרה אמיתית, ה-Circuit החדש הזה חייב גם הוא לעבור Transpile, באופן אופטימלי ככל האפשר, כדי לרוץ על ה-Backend שנבחר. שים לב שהTranspile לא שולב בפונקציות calculate_overlaps או cost_func_vqd. תרגיש חופשי לנסות לשנות את הקוד בעצמך כדי לשלב Transpile נוסף (ומותנה) זה — אך זה ייעשה עבורך גם בשיעור הבא.

בשיעור זה, נריץ את אלגוריתם VQD באמצעות Statevector Sampler ו-Statevector Estimator:

from qiskit.primitives import StatevectorEstimator as Estimator

sampler = Sampler()
estimator = Estimator()

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

from qiskit.quantum_info import SparsePauliOp

observable = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])

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

k = 3
betas = [33, 33, 33]
x0 = np.zeros(8)

כעת נוכל להריץ את החישוב:

from scipy.optimize import minimize

prev_states = []
prev_opt_parameters = []
eigenvalues = []

for step in range(1, k + 1):
    if step > 1:
        prev_states.append(ansatz.assign_parameters(prev_opt_parameters))

    result = minimize(
        cost_func_vqd,
        x0,
        args=(ansatz, prev_states, step, betas, estimator, sampler, observable),
        method="COBYLA",
        options={
            "maxiter": 200,
        },
    )
    print(result)

    prev_opt_parameters = result.x
    eigenvalues.append(result.fun)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -5.999999979545955
       x: [-5.150e-01 -5.452e-02 -1.571e+00 -2.853e-05  2.671e-01
           -2.672e-01 -8.509e-01 -8.510e-01]
    nfev: 131
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: 4.024550284767612
       x: [-3.745e-01  1.041e+00  8.637e-01  1.202e+00 -8.847e-02
            1.181e-02  7.611e-01 -3.006e-01]
    nfev: 110
   maxcv: 0.0
 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: 5.608925562838559
       x: [-2.670e-01  1.280e+00  1.070e+00 -8.031e-01 -1.524e-01
           -6.956e-02  7.018e-01  1.514e+00]
    nfev: 90
   maxcv: 0.0

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

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

לבסוף, נציין שכל שלושת המינימיזציות התכנסו לתוך הסובלנות ברירת המחדל של ה-optimizer הקלאסי (כאן COBYLA). הן דרשו 131, 110 ו-90 הערכות פונקציה, בהתאמה.

Quantum Sampling Regression (QSR)

אחת הבעיות העיקריות עם VQE היא הקריאות המרובות למחשב קוונטי הנדרשות לקבלת הפרמטרים עבור כל שלב, למשל, $k$, $k-1$ וכן הלאה. זה בעייתי במיוחד כאשר הגישה למכשירים קוונטיים מצויה בתור. בעוד ש-Session ניתן לשימוש לקיבוץ מספר קריאות איטרטיביות, גישה חלופית היא להשתמש בדגימה. על ידי ניצול יותר משאבים קלאסיים, ניתן להשלים את תהליך האופטימיזציה המלא בקריאה בודדת. כאן נכנסת Quantum Sampling Regression לתמונה. מכיוון שהגישה למחשבים קוונטיים היא עדיין מצרך עם היצע נמוך וביקוש גבוה, אנחנו מוצאים שהפשרה הזו אפשרית ונוחה למחקרים רבים בהווה. גישה זו מנצלת את כל היכולות הקלאסיות הזמינות תוך לכידת רבים מהעקרונות הפנימיים והתכונות המהותיות של חישובים קוונטיים שאינם מופיעים בסימולציה.

הרעיון מאחורי QSR הוא שפונקציית העלות $C(\theta) := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle$ ניתן לבטאה כסדרת פורייה באופן הבא:

$$\begin{aligned} C(\vec{\theta}) & := \langle \psi(\theta) | \hat{H} | \psi(\theta)\rangle \\[1mm] & := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)] \\[1mm] \end{aligned}$$

בהתאם לתקופתיות וברוחב הפס של הפונקציה המקורית, הקבוצה $S$ עשויה להיות סופית או אינסופית. לצורך דיון זה, נניח שהיא אינסופית. השלב הבא הוא לדגום את פונקציית העלות $C(\theta)$ מספר פעמים כדי לקבל את מקדמי פורייה $\{a_0, a_k, b_k\}_{k=1}^S$. ספציפית, מכיוון שיש לנו $2S+1$ נעלמים, נצטרך לדגום את פונקציית העלות $2S+1$ פעמים.

אם לאחר מכן נדגום את פונקציית העלות עבור $2S+1$ ערכי פרמטר $\{\theta_1,...,\theta_{2S+1}\}$, נוכל לקבל את המערכת הבאה:

$$\begin{pmatrix} 1 & \cos(\theta_1) & \sin(\theta_1) & \cos(2\theta_1) & ... & \sin(S\theta_1) \\ 1 & \cos(\theta_2) & \sin(\theta_2) & \cos(2\theta_2) & \cdots & \sin(S\theta_2)\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & \cos(\theta_{2S+1}) & \sin(\theta_{2S+1}) & \cos(2\theta_{2S+1}) & \cdots & \sin(S\theta_{2S+1}) \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ b_1 \\ a_2 \\ \vdots \\ b_S \end{pmatrix} = \begin{pmatrix} C(\theta_1) \\ C(\theta_2) \\ \vdots \\ C(\theta_{2S+1}) \end{pmatrix},$$

שנכתוב מחדש בתור

$$Fa=c.$$

בפועל, מערכת זו בדרך כלל אינה עקבית מכיוון שערכי פונקציית העלות $c$ אינם מדויקים. לכן, בדרך כלל טוב לנרמל אותם על ידי הכפלתם ב-$F^\dagger$ משמאל, אשר מניב:

$$F^\dagger Fa = F^\dagger c.$$

המערכת החדשה הזו תמיד עקבית, ופתרונה הוא פתרון הריבועים הפחותים של הבעיה המקורית. אם יש לנו $k$ פרמטרים במקום אחד בלבד, וכל פרמטר $\theta^i$ יש לו $S_i$ משלו עבור $i \in {1,...,k}$, אז המספר הכולל של הדגימות הנדרשות הוא:

$$T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n,$$

כאשר $S_{\max} = \max_i(S_i)$. יתר על כן, התאמת $S_{\max}$ כפרמטר ניתן לכיוון (במקום לאמוד אותו) פותחת אפשרויות חדשות, כגון:

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

מבנה תיאורטי

ניתן לסכם את המבנה של QSR כדלקמן:

• הכן אופרטורי ייחוס $U_R$.
  • נעבור מהמצב $|0\rangle$ למצב הייחוס $|\rho\rangle$
• הפעל את הצורה הווריאציונית $U_V(\vec\theta_{i,j})$ ליצירת ansatz $U_A(\vec\theta_{i,j})$.
  • קבע את רוחב הפס הקשור לכל פרמטר ב-ansatz. חסם עליון מספיק.
• אתחל ב-$i=0$ אם יש לנו בעיה דומה (נמצאת בדרך כלל על ידי סימולציה קלאסית או דגימה).
• דגום את פונקציית העלות $C(\vec\theta) := a_0 + \sum_{k=1}^S[a_k\cos(k\theta)+ b_k\sin(k\theta)]$ לפחות $T$ פעמים.
  • $T=\prod_{i=1}^k(2S_i+1)\leq \prod_{i=1}^k(2S_{max}+1) = (2S_{max}+1)^n$
  • החלט אם לדגום יתר/דגום חסר כדי לאזן בין מהירות לדיוק על ידי התאמת $T$.
• חשב את מקדמי פורייה מהדגימות (כלומר, פתור את מערכת המשוואות הליניארית המנורמלת).
• פתור עבור המינימום הגלובלי של פונקציית הרגרסיה המתקבלת במחשב קלאסי.

סיכום

עם שיעור זה, למדת על מספר מופעים וריאציוניים זמינים:

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