الخوارزميات الكمية: الخوارزميات الكمية التبايُنية

ملاحظة

Takashi Imamichi (24 May 2024)

تحميل ملف PDF للمحاضرة الأصلية. لاحظ أن بعض مقتطفات الشيفرة قد تصبح قديمة لأنها صور ثابتة.

الوقت التقريبي لوحدة المعالجة الكمية لتشغيل هذه التجربة هو 9 دقائق (تم الاختبار على معالج Eagle).

(قد لا يتم تقييم هذا الدفتر في الوقت المسموح به في الخطة المفتوحة. يرجى استخدام موارد الحوسبة الكمية بحكمة.)

1. مقدمة

يقدم هذا الدرس نظرة عامة على خوارزمية هجينة كمية-كلاسيكية، مع التركيز بشكل خاص على حلّال القيمة الذاتية الكمي التبايُني (VQE) وخوارزمية التحسين التقريبي الكمي (QAOA). الهدف الأساسي لهذه الخوارزميات هو معالجة مسائل التحسين من خلال استخدام دوائر كمية ذات بوابات كمية ذات معاملات.

على الرغم من التقدم في الحوسبة الكمية، فإن وجود الضوضاء في الأجهزة الكمية الحالية يجعل من الصعب استخلاص نتائج ذات معنى من الدوائر الكمية العميقة. للتغلب على هذا التحدي، يتبنى VQE و QAOA نهجًا هجينًا كميًا-كلاسيكيًا، يتضمن تنفيذ دوائر كمية قصيرة نسبيًا بشكل تكراري باستخدام الحوسبة الكمية وتحسين معاملات الدوائر الكمية ذات المعاملات المستهدفة باستخدام الحوسبة الكلاسيكية.

تتمتع QAOA بإمكانية تقديم الحلول المثلى للمسائل المستهدفة على نطاق المنفعة، بفضل تطبيق تقنيات مختلفة لتخفيف الأخطاء وقمعها. لدى VQE تطبيقات عديدة (مثل الكيمياء الكمية) حيث يكون أقل قابلية للتوسع. لكن ظهر عدد من المناهج المتعلقة بالقيم الذاتية لتكمل وتعزز VQE، بما في ذلك قطرنة فضاء Krylov الفرعي والقطرنة الكمومية القائم على أخذ العينات (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.]

للحصول على القيم الذاتية باستخدام خوارزمية كمية تبايُنية، نبني دائرة ذات بوابات تأخذ معاملات تبايُنية:

# 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

نجحت هذه الدائرة، لكن قيم المعاملات المختارة لم تقابل حالة ذات طاقة منخفضة جدًا (أو قيمة ذاتية منخفضة). القيمة الذاتية المحسوبة أعلى بكثير من الحد الأدنى. النتيجة مشابهة عند استخدام estimator.

لاحظ أن Estimator يأخذ دوائر كمية بدون قياسات.

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

للحصول على القيم الذاتية باستخدام خوارزمية كمية تبايُنية، نبني دائرة ذات بوابات تأخذ معاملات تبايُنية:

# 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

نجحت هذه الدائرة، لكن قيم المعاملات المختارة لم تقابل حالة ذات طاقة منخفضة جدًا (أو قيمة ذاتية منخفضة). القيمة الذاتية المحسوبة أعلى بكثير من الحد الأدنى. النتيجة مشابهة عند استخدام 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-hard) مع عدد من التطبيقات المختلفة في التجميع وعلم الشبكات والفيزياء الإحصائية. يعتبر هذا الدرس رسمًا بيانيًا من عُقد متصلة بحواف ويهدف إلى تقسيم العُقد إلى مجموعتين عن طريق "قطع" الحواف، بحيث يتم تعظيم عدد الحواف المقطوعة. لتوفير بعض السياق قبل تعيين هذه المسألة إلى خوارزمية كمية، يمكنك فهم كيف تصبح مسألة 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. تعيين المُدخلات الكلاسيكية إلى مسألة كمية

الخطوة الأولى من النمط هي تعيين المسألة الكلاسيكية (الرسم البياني) إلى دوائر و مؤثرات كمية. للقيام بذلك، هناك ثلاث خطوات رئيسية يجب اتخاذها:

1. استخدام سلسلة من إعادات الصياغة الرياضية، لتمثيل هذه المسألة باستخدام تدوين مسائل التحسين الثنائي غير المقيد التربيعي (QUBO).
2. إعادة كتابة مسألة التحسين كهاملتوني تقابل حالته الأساسية الحل الذي يصغّر دالة التكلفة.
3. إنشاء دائرة كمية ستحضّر الحالة الأساسية لهذا الهاملتوني عبر عملية مشابهة للتلدين الكمي.

ملاحظة: في منهجية QAOA، تريد في النهاية أن يكون لديك مؤثر (هاملتوني) يمثل دالة التكلفة لخوارزميتنا الهجينة، بالإضافة إلى دائرة ذات معاملات (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$ إلى مصفوفة Pauli $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

الهاملتوني → الدائرة الكمية

يحتوي الهاملتوني $H_C$ على التعريف الكمي لمسألتك. يمكنك الآن إنشاء دائرة كمية ستساعد في أخذ عينات من حلول جيدة من الحاسوب الكمي. تستوحي QAOA من التلدين الكمي وتطبق طبقات متناوبة من المؤثرات في الدائرة الكمية.

الفكرة العامة هي البدء من الحالة الأساسية لنظام معروف، $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~$.

الدائرة الكمية التي تنشئها ذات معاملات بـ $\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. تحسين الدوائر للتنفيذ على الأجهزة الكمية

تحتوي الدائرة أعلاه على سلسلة من التجريدات المفيدة للتفكير في الخوارزميات الكمية، لكنها غير قابلة للتشغيل على الأجهزة. لتتمكن من التشغيل على وحدة المعالجة الكمية، تحتاج الدائرة إلى سلسلة من العمليات التي تشكل خطوة التحويل أو تحسين الدائرة من النمط.

تقدم مكتبة Qiskit سلسلة من تمريرات التحويل التي تلبي نطاقًا واسعًا من تحويلات الدوائر. تحتاج إلى التأكد من أن دائرتك محسّنة لغرضك.

قد يتضمن التحويل عدة خطوات، مثل:

• التعيين الأولي للكيوبتات في الدائرة (مثل متغيرات القرار) إلى الكيوبتات الفيزيائية على الجهاز.
• التفكيك للتعليمات في الدائرة الكمية إلى التعليمات الأصلية للأجهزة التي تفهمها الخلفية.
• التوجيه لأي كيوبتات في الدائرة تتفاعل مع كيوبتات فيزيائية متجاورة مع بعضها البعض.
• قمع الأخطاء بإضافة بوابات أحادية الكيوبت لقمع الضوضاء بالفصل الديناميكي.

المزيد من المعلومات حول التحويل متاح في الوثائق.

تحوّل الشيفرة التالية الدائرة المجردة وتحسّنها إلى صيغة جاهزة للتنفيذ على أحد الأجهزة المتاحة عبر السحابة باستخدام خدمة 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 المثلى في حلقة تحسين تكرارية، تنفّذ سلسلة من تقييمات الدوائر وتستخدم محسّنًا كلاسيكيًا لإيجاد المعاملات المثلى $\beta_k$ و $\gamma_k$. يتم تنفيذ حلقة التنفيذ هذه عبر الخطوات التالية:

1. تحديد المعاملات الأولية
2. إنشاء Session جديدة تحتوي على حلقة التحسين والبدائية المستخدمة لأخذ عينات من الدائرة
3. بمجرد إيجاد مجموعة مثلى من المعاملات، تنفيذ الدائرة مرة أخيرة للحصول على توزيع نهائي سيُستخدم في خطوة المعالجة اللاحقة.

تعريف الدائرة بالمعاملات الأولية

نبدأ بمعاملات مختارة عشوائيًا.

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

تعريف الخلفية وبدائية التنفيذ

استخدم بدائيات Qiskit Runtime للتفاعل مع خلفيات IBM®. البدائيتان هما Sampler و Estimator، واختيار البدائية يعتمد على نوع القياس الذي تريد تشغيله على الحاسوب الكمي. لتصغير $H_C$، استخدم Estimator لأن قياس دالة التكلفة هو ببساطة القيمة المتوقعة لـ $\langle H_C \rangle$.

التشغيل

تقدم البدائيات مجموعة متنوعة من أوضاع التنفيذ لجدولة أعباء العمل على الأجهزة الكمية، ويعمل سير عمل QAOA بشكل تكراري في جلسة. يمكنك توصيل دالة التكلفة المبنية على 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

تمكّن المحسّن من تقليل التكلفة وإيجاد معاملات أفضل للدائرة.

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

بمجرد إيجاد المعاملات المثلى للدائرة، يمكنك تعيين هذه المعاملات وأخذ عينات من التوزيع النهائي المحصول عليه بالمعاملات المحسّنة. هنا يجب استخدام بدائية 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: توسيع النطاق!" من درس خوارزمية التحسين التقريبي الكمي.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'