אלגוריתמים קוונטיים: אלגוריתמים קוונטיים וריאציוניים

הערה

Takashi Imamichi (24 May 2024)

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

זמן QPU משוער להרצת הניסוי הזה הוא 9 דקות (נבדק על מעבד Eagle).

(מחברת זו עשויה שלא להסתיים בזמן המוקצב בתוכנית Open Plan. אנא השתמש במשאבי המחשוב הקוונטי בתבונה.)

1. מבוא

מדריך זה מספק סקירה של אלגוריתם היברידי קוונטי-קלאסי, תוך התמקדות ספציפית ב-variational quantum eigensolver (VQE) וב-quantum approximate optimization algorithm (QAOA). המטרה העיקרית של אלגוריתמים אלו היא להתמודד עם בעיות אופטימיזציה תוך שימוש ב-Circuit קוונטי עם Gate קוונטי בעל פרמטרים.

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

ל-QAOA יש פוטנציאל לספק פתרונות אופטימליים לבעיות הנבחרות בסקאלה של שימושיות, הודות ליישום טכניקות שונות להפחתת ודיכוי שגיאות. ל-VQE יישומים רבים (כמו כימיה קוונטית) שבהם הוא פחות ניתן לסקאלה. אך צמחו מספר גישות הקשורות לערכים עצמיים שמשלימות ומרחיבות את VQE, כולל Krylov subspace diagonalization ו-sampling-based quantum diagonalization (SQD). הבנת VQE היא צעד ראשון חשוב בהבנת מגוון הרחב של אלגוריתמים היברידיים קלאסיים-קוונטיים שצמחו.

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

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

# % pip install 'qiskit[visualization]' qiskit-ibm-runtime

2. חישוב הערך העצמי המינימלי של המילטוניאן פשוט

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

import numpy as np
from qiskit.circuit import ParameterVector, QuantumCircuit
from qiskit.primitives import StatevectorEstimator, StatevectorSampler
from qiskit.quantum_info import SparsePauliOp
from scipy.optimize import minimize

כעת נגדיר את האופרטור הנבחר ונצפה בו בצורת מטריצה.

op = SparsePauliOp("Z")
op.to_matrix()

array([[ 1.+0.j,  0.+0.j],
       [ 0.+0.j, -1.+0.j]])

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

# compute eigenvalues with numpy
result = np.linalg.eigh(op.to_matrix())
print("Eigenvalues:", result.eigenvalues)

Eigenvalues: [-1.  1.]

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

# define a variational form
param = ParameterVector("a", 3)
qc = QuantumCircuit(1, 1)
qc.u(param[0], param[1], param[2], 0)
qc_estimator = qc.copy()
qc.measure(0, 0)
qc.draw("mpl")

אם רוצים להעריך את ערך הציפייה של אופרטור (כמו $Z$), כדאי להשתמש ב-Estimator. אם רוצים לבחון את המצבים של המערכת, משתמשים ב-Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

ניתן לחשב ספירות של מחרוזות ביט 0 ו-1 עם ערכי פרמטרים אקראיים [1, 2, 3] באמצעות Sampler.

# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3])]).result()
counts = result[0].data.c.get_counts()
counts

{'0': 783, '1': 241}

אנו יודעים שניתן לחשב את ערך הציפייה של Z על ידי $\langle Z \rangle = p_0 - p_1$ עם ההסתברויות $\{0: p_0, 1: p_1\}$.

# compute the expectation value of Z based on the counts
(counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())

0.529296875

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

שים לב ש-Estimator מקבל Circuit קוונטי ללא מדידות.

result = estimator.run([(qc_estimator, op, [1, 2, 3])]).result()
result[0].data.evs

array(0.54030231)

נצטרך לחפש בין הפרמטרים ולמצוא את אלה שמניבים את הערך העצמי הנמוך ביותר. נגדיר פונקציה שמקבלת ערכי פרמטרים של הצורה הוריאציונית ומחזירה את ערך הציפייה $\langle Z \rangle$.

# define a cost function to look for the minimum eigenvalue of Z
def cost(x):
    result = sampler.run([(qc, x)]).result()
    counts = result[0].data.c.get_counts()
    expval = (counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())
    # the following line shows the trajectory of the optimization
    print(expval, counts)
    return expval

נשתמש בפונקציית minimize של SciPy כדי למצוא את הערך העצמי המינימלי של Z.

# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0], method="COBYLA", tol=1e-8)
min_result

1.0 {'0': 1024}
0.494140625 {'0': 765, '1': 259}
0.466796875 {'0': 751, '1': 273}
0.564453125 {'0': 801, '1': 223}
-0.4296875 {'1': 732, '0': 292}
-0.984375 {'1': 1016, '0': 8}
-0.8984375 {'1': 972, '0': 52}
-0.990234375 {'1': 1019, '0': 5}
-0.892578125 {'1': 969, '0': 55}
-0.986328125 {'1': 1017, '0': 7}
-0.861328125 {'1': 953, '0': 71}
-1.0 {'1': 1024}
-0.982421875 {'1': 1015, '0': 9}
-0.99609375 {'1': 1022, '0': 2}
-0.986328125 {'1': 1017, '0': 7}
-1.0 {'1': 1024}
-0.990234375 {'1': 1019, '0': 5}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.994140625 {'1': 1021, '0': 3}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -1.0
       x: [ 3.182e+00  1.338e+00  1.664e-01]
    nfev: 63
   maxcv: 0.0

# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()

{'0': 1, '1': 1023}

2.1 תרגיל

חשב את ערך הייגן המינימלי של $Z \otimes Z$ עם VQE.

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())

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

# compute eigenvalues with numpy

# define a variational form
# qc = ...

# compute counts of bitstrings with a random parameter values by Sampler
# result = sampler.run(...)
# result

# compute the expectation value of ZZ based on the counts

# verify the expectation value of ZZ with Estimator

# define a cost function to look for the minimum eigenvalue of ZZ
# def cost(x):
#    expval = ...
#    return expval

# minimize the cost function with scipy's minimize
# min_result = minimize(cost, [...], method="COBYLA", tol=1e-8)
# min_result

# check counts of bitstrings with the optimal parameter values
# result = sampler.run(qc, min_result.x).result()
# result

פתרונות התרגיל

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

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())

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

כדי לקבל ערכי ייגן באמצעות אלגוריתם קוונטי וריאציוני, אנחנו בונים Circuit עם Gate-ים שמקבלים פרמטרים וריאציוניים:

# define a variational form
param = ParameterVector("a", 6)
qc = QuantumCircuit(2, 2)
qc.u(param[0], param[1], param[2], 0)
qc.u(param[3], param[4], param[5], 1)
qc_estimator = qc.copy()
qc.measure([0, 1], [0, 1])
qc.draw("mpl")

אם אנחנו רוצים להעריך את ערך הציפייה של אופרטור (כמו $Z \otimes Z$), נשתמש ב-Estimator. אם אנחנו רוצים לבחון את המצבים של המערכת, נשתמש ב-Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3, 4, 5, 6])]).result()
counts = result[0].data.c.get_counts()
counts

{'10': 661, '11': 203, '01': 47, '00': 113}

# compute the expectation value of ZZ based on the counts
(
    counts.get("00", 0)
    - counts.get("01", 0)
    - counts.get("10", 0)
    + counts.get("11", 0)
) / sum(counts.values())

-0.3828125

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

# verify the expectation value of ZZ with Estimator
result = estimator.run([(qc_estimator, z2, [1, 2, 3, 4, 5, 6])]).result()
result[0].data.evs

array(-0.35316516)

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

# define a cost function to look for the minimum eigenvalue of ZZ
def cost(x):
    result = sampler.run([(qc, x)]).result()
    counts = result[0].data.c.get_counts()
    expval = (
        counts.get("00", 0)
        - counts.get("01", 0)
        - counts.get("10", 0)
        + counts.get("11", 0)
    ) / sum(counts.values())
    print(expval, counts)
    return expval

# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0, 0, 0, 0], method="COBYLA", tol=1e-8)
min_result

1.0 {'00': 1024}
0.578125 {'00': 808, '01': 216}
0.5234375 {'00': 780, '01': 244}
0.548828125 {'00': 793, '01': 231}
0.3515625 {'00': 637, '10': 164, '11': 55, '01': 168}
0.3359375 {'00': 638, '11': 46, '10': 174, '01': 166}
0.283203125 {'00': 602, '10': 181, '01': 186, '11': 55}
-0.087890625 {'01': 414, '00': 184, '10': 143, '11': 283}
0.236328125 {'10': 27, '11': 623, '01': 364, '00': 10}
-0.0625 {'11': 261, '01': 403, '00': 219, '10': 141}
0.248046875 {'01': 366, '11': 628, '00': 11, '10': 19}
-0.0625 {'10': 145, '11': 254, '01': 399, '00': 226}
0.228515625 {'01': 373, '11': 609, '00': 20, '10': 22}
0.0546875 {'11': 376, '10': 273, '01': 211, '00': 164}
-0.447265625 {'01': 731, '10': 10, '11': 267, '00': 16}
-0.71484375 {'01': 871, '11': 99, '00': 47, '10': 7}
-0.46484375 {'01': 741, '00': 253, '10': 9, '11': 21}
-0.87890625 {'01': 962, '00': 39, '11': 23}
-0.640625 {'00': 176, '01': 837, '11': 8, '10': 3}
-0.88671875 {'01': 966, '00': 41, '11': 17}
-0.994140625 {'01': 1021, '11': 3}
-0.91796875 {'01': 982, '11': 35, '00': 7}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.939453125 {'01': 993, '00': 31}
-0.990234375 {'01': 1019, '11': 5}
-0.90234375 {'01': 974, '00': 21, '11': 29}
-0.98046875 {'01': 1014, '11': 10}
-0.994140625 {'01': 1021, '00': 3}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.98828125 {'01': 1018, '11': 6}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '00': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '00': 1, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '11': 1, '00': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.994140625 {'01': 1021, '00': 3}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.998046875
       x: [ 3.167e+00  6.940e-01  1.033e+00 -2.894e-02  8.933e-01
            1.885e+00]
    nfev: 128
   maxcv: 0.0

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.99609375
       x: [ 3.098e+00 -5.402e-01  1.091e+00 -1.004e-02  3.615e-01
            6.913e-01]
    nfev: 115
   maxcv: 0.0

קיבלנו ערך ייגן קרוב מאוד למינימום שנתן לנו numpy.

# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()

{'01': 1024}

3. אופטימיזציה קוונטית עם תבניות Qiskit

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

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

3.1 תבנית Qiskit (בקנה מידה קטן) לאופטימיזציה

חלק זה ישתמש בבעיית Max-Cut בקנה מידה קטן כדי להמחיש את השלבים הנדרשים לפתרון בעיית אופטימיזציה באמצעות מחשב קוונטי. בעיית Max-Cut היא בעיית אופטימיזציה שקשה לפתור (ביתר דיוק, זוהי בעיה NP-קשה) עם מגוון יישומים בקיבוץ, במדע הרשתות ובפיזיקה סטטיסטית. מדריך זה שוקל גרף של צמתים המחוברים בקשתות ומטרתו לחלק את הצמתים לשתי קבוצות על ידי "חיתוך" קשתות, כך שמספר הקשתות הנחתכות יהיה מרבי. כדי לתת קצת הקשר לפני מיפוי הבעיה לאלגוריתם קוונטי, ניתן להבין טוב יותר כיצד בעיית Max-Cut הופכת לבעיית אופטימיזציה קומבינטורית קלאסית על ידי התחשבות תחילה במינימיזציה של פונקציה $f(x)$

$$\min_{x\in \{0, 1\}^n}f(x),$$

כאשר הקלט $x$ הוא וקטור שרכיביו מתאימים לכל צומת של גרף. לאחר מכן, מגבילים כל אחד מהרכיבים הללו להיות $0$ או $1$ (המייצגים הכללה או אי-הכללה בחיתוך). דוגמת קנה המידה הקטן זה משתמשת בגרף עם $n=5$ צמתים.

ניתן לכתוב פונקציה של זוג צמתים $i,j$ המציינת האם הקשת המתאימה $(i,j)$ נמצאת בחיתוך. לדוגמה, הפונקציה $x_i + x_j - 2 x_i x_j$ שווה ל-1 רק אם אחד מ-$x_i$ או $x_j$ שווה ל-1 (כלומר הקשת נמצאת בחיתוך) ואפס אחרת. בעיית המקסימום של הקשתות בחיתוך ניתנת לניסוח כ:

$$\max_{x\in \{0, 1\}^n} \sum_{(i,j)} x_i + x_j - 2 x_i x_j,$$

אשר ניתן לשכתב כמינימיזציה בצורה

$$\min_{x\in \{0, 1\}^n} \sum_{(i,j)} 2 x_i x_j - x_i - x_j.$$

המינימום של $f(x)$ במקרה זה הוא כאשר מספר הקשתות שחוצה החיתוך הוא מרבי. כפי שניתן לראות, עדיין אין שום קשר למחשוב קוונטי. יש לנסח מחדש את הבעיה למשהו שמחשב קוונטי יוכל להבין. אתחלו את הבעיה על ידי יצירת גרף עם $n=5$ צמתים.

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 5

graph = rx.PyGraph()
graph.add_nodes_from(range(1, n + 1))
edge_list = [
    (0, 1, 1.0),
    (0, 2, 1.0),
    (1, 2, 1.0),
    (1, 3, 1.0),
    (2, 4, 1.0),
    (3, 4, 1.0),
]
graph.add_edges_from(edge_list)
pos = rx.spring_layout(graph, seed=2)
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str)

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

השלב הראשון של התבנית הוא למפות את הבעיה הקלאסית (גרף) ל-Circuitים ו-אופרטורים קוונטיים. לשם כך, ישנם שלושה שלבים עיקריים:

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

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

גרף → בעיית אופטימיזציה

השלב הראשון של המיפוי הוא שינוי סימון. הביטוי הבא מביע את הבעיה בסימון QUBO:

$$\min_{x\in \{0, 1\}^n}x^T Q x,$$

כאשר $Q$ הוא מטריצה של $n\times n$ מספרים ממשיים, $n$ מתאים למספר הצמתים בגרף, $x$ הוא וקטור המשתנים הבינאריים שהוצג לעיל, ו-$x^T$ מציין את הטרנספוז של הוקטור $x$.

Problem name: maxcut

Minimize
  2*x_1*x_2 + 2*x_1*x_3 + 2*x_2*x_3 + 2*x_2*x_4 + 2*x_3*x_5 + 2*x_4*x_5 - 2*x_1
  - 3*x_2 - 3*x_3 - 2*x_4 - 2*x_5

Subject to
  No constraints

  Binary variables (5)
    x_1 x_2 x_3 x_4 x_5

בעיית אופטימיזציה → המילטוניאן

לאחר מכן ניתן לנסח מחדש את בעיית QUBO כ-המילטוניאן (כאן, מטריצה המייצגת את האנרגיה של מערכת):

$$H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.$$

שלבי הניסוח מחדש מבעיית QAOA להמילטוניאן

כדי להדגים כיצד ניתן לשכתב את בעיית QAOA בדרך זו, ראשית יש להחליף את המשתנים הבינאריים $x_i$ בקבוצה חדשה של משתנים $z_i\in\{-1, 1\}$ באמצעות

$$x_i = \frac{1-z_i}{2}.$$

כאן ניתן לראות שאם $x_i$ הוא $0$, אז $z_i$ חייב להיות $1$. כאשר ה-$x_i$-ים מוחלפים ב-$z_i$-ים בבעיית האופטימיזציה ($x^TQx$), ניתן לקבל ניסוח שקול.

$$x^TQx=\sum_{ij}Q_{ij}x_ix_j \\ =\frac{1}{4}\sum_{ij}Q_{ij}(1-z_i)(1-z_j) \\=\frac{1}{4}\sum_{ij}Q_{ij}z_iz_j-\frac{1}{4}\sum_{ij}(Q_{ij}+Q_{ji})z_i + \frac{n^2}{4}.$$

כעת אם נגדיר $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$, נסיר את הגורם הקדמי ואת האיבר הקבוע $n^2$, נגיע לשני ניסוחים שקולים של אותה בעיית אופטימיזציה.

$$min_{x\in\{0,1\}^n} x^TQx\Longleftrightarrow \min_{z\in\{-1,1\}^n}z^TQz + b^Tz$$

כאן, $b$ תלוי ב-$Q$. שימו לב שכדי לקבל את $z^TQz + b^Tz$ השמטנו את גורם ה-1/4 ואת ההיסט הקבוע $n^2$ שאינם משחקים תפקיד באופטימיזציה.

כעת, כדי לקבל ניסוח קוונטי של הבעיה, יש לקדם את משתני $z_i$ למטריצת פאולי $Z$, כגון מטריצה $2\times 2$ בצורה

$$Z_i = \begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}.$$

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

$$H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.$$

כמו כן, זכרו שמטריצות ה-$Z$ מוטבעות במרחב החישובי של המחשב הקוונטי, כלומר מרחב הילברט בגודל $2^n\times 2^n$. לפיכך, יש להבין מונחים כמו $Z_iZ_j$ כמכפלת הטנסור $Z_i\otimes Z_j$ המוטבעת במרחב הילברט $2^n\times 2^n$. לדוגמה, בבעיה עם חמישה משתני החלטה המונח $Z_1Z_3$ מובן כ-$I\otimes Z_3\otimes I\otimes Z_1\otimes I$ כאשר $I$ היא מטריצת הזהות $2\times 2$.

המילטוניאן זה נקרא המילטוניאן פונקציית העלות. יש לו את התכונה שהמצב הבסיסי שלו מתאים לפתרון הממזער את פונקציית העלות $f(x)$. לפיכך, כדי לפתור את בעיית האופטימיזציה עכשיו צריך להכין את המצב הבסיסי של $H_C$ (או מצב עם חפיפה גבוהה אליו) במחשב הקוונטי. לאחר מכן, דגימה ממצב זה תניב, בהסתברות גבוהה, את הפתרון ל-$\min~f(x)$.

def build_max_cut_operator(graph: rx.PyGraph) -> tuple[SparsePauliOp, float]:
    sp_list = []
    constant = 0
    for s, t in graph.edge_list():
        w = graph.get_edge_data(s, t)
        sp_list.append(("ZZ", [s, t], w / 2))
        constant -= 1 / 2
    return SparsePauliOp.from_sparse_list(
        sp_list, num_qubits=graph.num_nodes()
    ), constant

cost_hamiltonian, constant = build_max_cut_operator(graph)
print("Cost Function Hamiltonian:", cost_hamiltonian)
print("Constant:", constant)

Cost Function Hamiltonian: SparsePauliOp(['IIIZZ', 'IIZIZ', 'IIZZI', 'IZIZI', 'ZIZII', 'ZZIII'],
              coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])
Constant: -3.0

המילטוניאן → Circuit קוונטי

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

הרעיון הכללי הוא להתחיל במצב הבסיסי של מערכת ידועה, $H^{\otimes n}|0\rangle$ לעיל, ולאחר מכן להוביל את המערכת אל המצב הבסיסי של אופרטור העלות שמתעניינים בו. זה נעשה על ידי החלת האופרטורים $\exp\{-i\gamma_k H_C\}$ ו-$\exp\{-i\beta_k H_m\}$ עם זוויות $\gamma_1,...,\gamma_p$ ו-$\beta_1,...,\beta_p~$.

ה-Circuit הקוונטי שיוצרים הוא פרמטרי ב-$\gamma_i$ ו-$\beta_i$, כך שניתן לנסות ערכים שונים של $\gamma_i$ ו-$\beta_i$ ולדגום מהמצב המתקבל. במקרה זה ננסה דוגמה עם שכבת QAOA אחת המכילה שני פרמטרים: $\gamma_1$ ו-$\beta_1$.

from qiskit.circuit.library import QAOAAnsatz

circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=1)
circuit.measure_all()
circuit.draw("mpl")

circuit.decompose(reps=3).draw("mpl", fold=-1)

circuit.parameters

ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(γ[0])])

3.3 שלב 2. אופטימיזציה של Circuit-ים לביצוע על חומרה קוונטית

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

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

ה-Transpilation עשויה לכלול מספר שלבים, כגון:

• מיפוי ראשוני של ה-Qubit-ים ב-Circuit (כגון משתני החלטה) ל-Qubit-ים פיזיים על המכשיר.
• פריסה של ההוראות ב-Circuit הקוונטי להוראות ילידיות לחומרה שה-Backend מבין.
• ניתוב של כל Qubit-ים ב-Circuit המתקשרים ל-Qubit-ים פיזיים הסמוכים זה לזה.
• דיכוי שגיאות על ידי הוספת Gate-ים של Qubit בודד לדיכוי רעש עם פענוח דינמי.

מידע נוסף על Transpilation זמין בתיעוד.

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

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

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Use a quantum device
service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
# backend = service.backend("ibm_kingston")

# You can test your programs locally with a fake backend (local testing mode)
# backend = FakeBrisbane()

print(backend)

# Create pass manager for transpilation
pm = generate_preset_pass_manager(optimization_level=3, backend=backend)

candidate_circuit = pm.run(circuit)
candidate_circuit.draw("mpl", fold=False, idle_wires=False)

service = QiskitRuntimeService(channel="ibm_quantum_platform")

<IBMBackend('ibm_strasbourg')>

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

בזרימת העבודה של QAOA, הפרמטרים האופטימליים של QAOA נמצאים בלולאת אופטימיזציה איטרטיבית, שמריצה סדרה של הערכות Circuit ומשתמשת במייעל קלאסי למציאת הפרמטרים האופטימליים $\beta_k$ ו-$\gamma_k$. לולאת ביצוע זו מתבצעת דרך השלבים הבאים:

1. הגדרת הפרמטרים הראשוניים
2. יצירת Session חדשה המכילה את לולאת האופטימיזציה ואת הפרימיטיב המשמש לדגימת ה-Circuit
3. לאחר מציאת קבוצה אופטימלית של פרמטרים, הרצת ה-Circuit פעם אחרונה לקבלת התפלגות סופית שתשמש בשלב עיבוד ה-Post.

הגדרת Circuit עם פרמטרים ראשוניים

מתחילים עם פרמטרים שנבחרו באופן שרירותי.

initial_gamma = np.pi
initial_beta = np.pi / 2
init_params = [initial_gamma, initial_beta]

הגדרת Backend ופרימיטיב הביצוע

השתמשו ב-פרימיטיבים של Qiskit Runtime לאינטראקציה עם Backend-ים של IBM®. שני הפרימיטיבים הם Sampler ו-Estimator, והבחירה בפרימיטיב תלויה בסוג המדידה שרוצים להריץ על המחשב הקוונטי. למינימיזציה של $H_C$, השתמשו ב-Estimator מאחר שמדידת פונקציית העלות היא פשוט ערך הציפייה של $\langle H_C \rangle$.

הרצה

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

def cost_func_estimator(params, ansatz, hamiltonian, estimator):
    # transform the observable defined on virtual qubits to
    # an observable defined on all physical qubits
    isa_hamiltonian = hamiltonian.apply_layout(ansatz.layout)

    pub = (ansatz, isa_hamiltonian, params)
    job = estimator.run([pub])

    results = job.result()[0]
    cost = results.data.evs

    objective_func_vals.append(cost)

    return cost

from qiskit_ibm_runtime import Session, EstimatorV2
from scipy.optimize import minimize

objective_func_vals = []  # Global variable
with Session(backend=backend) as session:
    # If using qiskit-ibm-runtime<0.24.0, change `mode=` to `session=`
    estimator = EstimatorV2(mode=session)
    estimator.options.default_shots = 1000

    # Set simple error suppression/mitigation options
    estimator.options.dynamical_decoupling.enable = True
    estimator.options.dynamical_decoupling.sequence_type = "XY4"
    estimator.options.twirling.enable_gates = True
    estimator.options.twirling.num_randomizations = "auto"

    result = minimize(
        cost_func_estimator,
        init_params,
        args=(candidate_circuit, cost_hamiltonian, estimator),
        method="COBYLA",
        tol=1e-2,
    )
    print(result)

message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.6557925874481715
       x: [ 2.873e+00  9.414e-01]
    nfev: 21
   maxcv: 0.0

המייעל הצליח להפחית את העלות ולמצוא פרמטרים טובים יותר עבור ה-Circuit.

plt.figure(figsize=(12, 6))
plt.plot(objective_func_vals)
plt.xlabel("Iteration")
plt.ylabel("Cost")
plt.show()

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

הערה: משמעות הדבר היא הכנת מצב קוונטי $\psi$ במחשב ולאחר מכן מדידתו. מדידה תכווץ את המצב למצב בסיס חישובי בודד - לדוגמה, 010101110000... - המתאים לפתרון מועמד $x$ לבעיית האופטימיזציה הראשונית שלנו ($\max f(x)$ או $\min f(x)$ בהתאם למשימה).

optimized_circuit = candidate_circuit.assign_parameters(result.x)
optimized_circuit.draw("mpl", fold=False, idle_wires=False)

from qiskit_ibm_runtime import SamplerV2

# If using qiskit-ibm-runtime<0.24.0, change `mode=` to `backend=`
sampler = SamplerV2(mode=backend)

# Set simple error suppression/mitigation options
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XY4"
sampler.options.twirling.enable_gates = True
sampler.options.twirling.num_randomizations = "auto"

pub = (optimized_circuit,)
job = sampler.run([pub], shots=int(1e4))
counts_int = job.result()[0].data.meas.get_int_counts()
counts_bin = job.result()[0].data.meas.get_counts()
shots = sum(counts_int.values())
final_distribution_int = {key: val / shots for key, val in counts_int.items()}
final_distribution_bin = {key: val / shots for key, val in counts_bin.items()}
print(final_distribution_int)

{12: 0.0652, 31: 0.0089, 4: 0.0085, 13: 0.0731, 26: 0.0256, 28: 0.0246, 17: 0.0405, 25: 0.0591, 20: 0.031, 15: 0.0221, 8: 0.017, 21: 0.0371, 14: 0.0461, 16: 0.0229, 19: 0.0723, 23: 0.0199, 22: 0.0478, 18: 0.0708, 24: 0.0165, 6: 0.0525, 7: 0.0155, 5: 0.0245, 3: 0.0231, 29: 0.0121, 30: 0.0062, 10: 0.0363, 1: 0.0097, 9: 0.042, 27: 0.0094, 11: 0.0349, 0: 0.0129, 2: 0.0119}

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

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

# auxiliary functions to sample most likely bitstring
def to_bitstring(integer, num_bits):
    result = np.binary_repr(integer, width=num_bits)
    return [int(digit) for digit in result]

keys = list(final_distribution_int.keys())
values = list(final_distribution_int.values())
most_likely = keys[np.argmax(np.abs(values))]
most_likely_bitstring = to_bitstring(most_likely, len(graph))
most_likely_bitstring.reverse()

print("Result bitstring:", most_likely_bitstring)

Result bitstring: [1, 0, 1, 1, 0]

import matplotlib.pyplot as plt

matplotlib.rcParams.update({"font.size": 10})
final_bits = final_distribution_bin
values = np.abs(list(final_bits.values()))
top_4_values = sorted(values, reverse=True)[:4]
positions = []
for value in top_4_values:
    positions.append(np.where(values == value)[0])
fig = plt.figure(figsize=(11, 6))
ax = fig.add_subplot(1, 1, 1)
plt.xticks(rotation=45)
plt.title("Result Distribution")
plt.xlabel("Bitstrings (reversed)")
plt.ylabel("Probability")
ax.bar(list(final_bits.keys()), list(final_bits.values()), color="tab:grey")
for p in positions:
    ax.get_children()[p[0].item()].set_color("tab:purple")
plt.show()

ויזואליזציה של החיתוך הטוב ביותר

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

colors = ["tab:grey" if i == 0 else "tab:purple" for i in most_likely_bitstring]
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str, node_color=colors)

ולאחר מכן לחשב את ערך החיתוך. הפתרון אינו אופטימלי בשל רעש (ערך החיתוך של הפתרון האופטימלי הוא 5).

from typing import Sequence

def evaluate_sample(x: Sequence[int], graph: rx.PyGraph) -> float:
    assert len(x) == len(
        list(graph.nodes())
    ), "The length of x must coincide with the number of nodes in the graph."
    return sum(
        x[u] * (1 - x[v]) + x[v] * (1 - x[u]) for u, v in list(graph.edge_list())
    )

cut_value = evaluate_sample(most_likely_bitstring, graph)
print("The value of the cut is:", cut_value)

The value of the cut is: 5

זהו סיום מדריך QAOA בקנה מידה קטן. תלמד כיצד להתאים QAOA בקנה מידה שימושי ב"חלק 2: להגדיל את הקנה מידה!" של המדריך Quantum approximate optimization algorithm.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'