أمثلة وتطبيقات
في هذا الدرس، سنستعرض بعض أمثلة الخوارزميات التغايرية وكيفية تطبيقها:
- كيفية كتابة خوارزمية تغايرية مخصصة
- كيفية تطبيق خوارزمية تغايرية لإيجاد أدنى القيم الذاتية
- كيفية توظيف الخوارزميات التغايرية لحل حالات استخدام تطبيقية
لاحظ أن إطار عمل أنماط Qiskit يمكن تطبيقه على جميع المسائل التي نقدمها هنا. غير أنه تجنبًا للتكرار، سنُشير صراحةً إلى خطوات الإطار في مثال واحد فقط، يُشغَّل على عتاد حقيقي.
تعريف المسائل
تخيّل أننا نريد استخدام خوارزمية تغايرية لإيجاد القيمة الذاتية للمرصود التالي:
يمتلك هذا المرصود القيم الذاتية التالية:
والحالات الذاتية:
# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-ibm-runtime rustworkx scipy
from qiskit.quantum_info import SparsePauliOp
observable_1 = SparsePauliOp.from_list([("II", 2), ("XX", -2), ("YY", 3), ("ZZ", -3)])
VQE مخصص
سنستكشف أولًا كيفية بناء نموذج VQE يدويًا لإيجاد أدنى قيمة ذاتية لـ . سيستوعب هذا مجموعة من التقنيات التي تناولناها على مدار هذه الدورة.
def cost_func_vqe(params, 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
"""
pub = (ansatz, hamiltonian, params)
cost = estimator.run([pub]).result()[0].data.evs
return cost
from qiskit.circuit.library.n_local import n_local
from qiskit import QuantumCircuit
import numpy as np
reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)
variational_form = n_local(
num_qubits=2,
rotation_blocks=["rz", "ry"],
entanglement_blocks="cx",
entanglement="linear",
reps=1,
)
raw_ansatz = reference_circuit.compose(variational_form)
raw_ansatz.decompose().draw("mpl")

سنبدأ بتنقيح الكود على محاكيات محلية.
from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler
estimator = Estimator()
sampler = Sampler()
الآن نضبط مجموعة أولية من المعاملات:
x0 = np.ones(raw_ansatz.num_parameters)
print(x0)
[1. 1. 1. 1. 1. 1. 1. 1.]
يمكننا تصغير دالة التكلفة هذه لحساب المعاملات المثلى
# SciPy minimizer routine
from scipy.optimize import minimize
import time
start_time = time.time()
result = minimize(
cost_func_vqe,
x0,
args=(raw_ansatz, observable_1, estimator),
method="COBYLA",
options={"maxiter": 1000, "disp": True},
)
end_time = time.time()
execution_time = end_time - start_time
Return from COBYLA because the trust region radius reaches its lower bound.
Number of function values = 103 Least value of F = -5.999999998357189
The corresponding X is:
[2.27483579e+00 8.37593091e-01 1.57080508e+00 5.82932911e-06
2.49973063e+00 6.41884255e-01 6.33686904e-01 6.33688223e-01]
result
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -5.999999998357189
x: [ 2.275e+00 8.376e-01 1.571e+00 5.829e-06 2.500e+00
6.419e-01 6.337e-01 6.337e-01]
nfev: 103
maxcv: 0.0
بما أن هذه المسألة التجريبية تستخدم كيوبتين فقط، يمكننا التحقق من ذلك باستخدام حلّال القيم الذاتية في الجبر الخطي بـ NumPy.
from numpy.linalg import eigvalsh
solution_eigenvalue = min(eigvalsh(observable_1.to_matrix()))
print(f"""Number of iterations: {result.nfev}""")
print(f"""Time (s): {execution_time}""")
print(
f"Percent error: {100*abs((result.fun - solution_eigenvalue)/solution_eigenvalue):.2e}"
)
Number of iterations: 103
Time (s): 0.4394676685333252
Percent error: 2.74e-08
كما ترى، النتيجة قريبة جدًا من القيمة المثالية.
التجريب لتحسين السرعة والدقة
إضافة حالة مرجعية
في المثال السابق لم نستخدم أي معامل مرجعي . الآن دعنا نفكر في كيفية الحصول على الحالة الذاتية المثالية . تأمل الدائرة التالية.
from qiskit import QuantumCircuit
ideal_qc = QuantumCircuit(2)
ideal_qc.h(0)
ideal_qc.cx(0, 1)
ideal_qc.draw("mpl")
يمكننا التحقق بسرعة من أن هذه الدائرة تعطينا الحالة المطلوبة.
from qiskit.quantum_info import Statevector
Statevector(ideal_qc)
Statevector([0.70710678+0.j, 0. +0.j, 0. +0.j,
0.70710678+0.j],
dims=(2, 2))
الآن بعد أن رأينا كيف تبدو الدائرة التي تُحضّر حالة الحل، يبدو من المنطقي استخدام بوابة Hadamard كدائرة مرجعية، بحيث يصبح ansatz الكامل:
reference = QuantumCircuit(2)
reference.h(0)
reference.cx(0, 1)
# Include barrier to separate reference from variational form
reference.barrier()
ref_ansatz = variational_form.decompose().compose(reference, front=True)
ref_ansatz.draw("mpl")

بالنسبة لهذه الدائرة الجديدة، يمكن بلوغ الحل المثالي حين تكون جميع المعاملات مساوية لـ ، وهذا يؤكد أن اختيار الدائرة المرجعية كان معقولاً.
الآن دعنا نقارن عدد تقييمات دالة التكلفة، وعدد تكرارات المُحسِّن، والوقت المستغرق، بما كانت عليه في المحاولة السابقة.
import time
start_time = time.time()
ref_result = minimize(
cost_func_vqe, x0, args=(ref_ansatz, observable_1, estimator), method="COBYLA"
)
end_time = time.time()
execution_time = end_time - start_time
باستخدام المعاملات المثلى لحساب القيمة الذاتية الدنيا:
experimental_min_eigenvalue_ref = cost_func_vqe(
ref_result.x, ref_ansatz, observable_1, estimator
)
print(experimental_min_eigenvalue_ref)
-5.999999996759607
print("ADDED REFERENCE STATE:")
print(f"""Number of iterations: {ref_result.nfev}""")
print(f"""Time (s): {execution_time}""")
print(
f"Percent error: {100*abs((experimental_min_eigenvalue_ref - solution_eigenvalue)/solution_eigenvalue):.2e}"
)
ADDED REFERENCE STATE:
Number of iterations: 127
Time (s): 0.5620882511138916
Percent error: 5.40e-08
قد يُحقق هذا تحسينًا في السرعة أو الدقة أو لا، وذلك اعتمادًا على نظامك المحدد في هذا المثال الصغير جدًا. الفكرة الأساسية هي أن البدء بحالات مرجعية مستوحاة فيزيائيًا يكتسب أهمية متزايدة في تحسين السرعة والدقة كلما تضخّمت المسائل.
تغيير نقطة البداية
الآن بعد أن رأينا تأثير إضافة الحالة المرجعية، سنستعرض ما يحدث حين نختار نقاط بداية مختلفة . سنستخدم تحديدًا و .
تذكّر أنه كما نوقش عند تقديم الحالة المرجعية، يُعثر على الحل المثالي حين تكون جميع المعاملات ، لذا ينبغي لنقطة البداية الأولى أن تستلزم تقييمات أقل.
import time
start_time = time.time()
x0 = [0, 0, 0, 0, 6, 0, 0, 0]
x0_1_result = minimize(
cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="COBYLA"
)
end_time = time.time()
execution_time = end_time - start_time
print("INITIAL POINT 1:")
print(f"""Number of iterations: {x0_1_result.nfev}""")
print(f"""Time (s): {execution_time}""")
INITIAL POINT 1:
Number of iterations: 108
Time (s): 0.4492197036743164
ضبط نقطة البداية إلى :
import time
start_time = time.time()
x0 = 6 * np.ones(raw_ansatz.num_parameters)
x0_2_result = minimize(
cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="COBYLA"
)
end_time = time.time()
execution_time = end_time - start_time
print("INITIAL POINT 2:")
print(f"""Number of iterations: {x0_2_result.nfev}""")
print(f"""Time (s): {execution_time}""")
INITIAL POINT 2:
Number of iterations: 107
Time (s): 0.40889453887939453
من خلال التجريب بنقاط بداية مختلفة، قد تتمكن من تحقيق التقارب بشكل أسرع ومع تقييمات أقل للدوال.
التجريب بمُحسِّنات مختلفة
يمكننا ضبط المُحسِّن باستخدام وسيطة method الخاصة بـ minimize في SciPy، ويمكن الاطلاع على مزيد من الخيارات هنا. استخدمنا في الأصل مُحسِّنًا مقيَّدًا (COBYLA). في هذا المثال، سنستكشف استخدام مُحسِّن غير مقيَّد (BFGS) بدلاً منه.
import time
start_time = time.time()
result = minimize(
cost_func_vqe, x0, args=(raw_ansatz, observable_1, estimator), method="BFGS"
)
end_time = time.time()
execution_time = end_time - start_time
print("CHANGED TO BFGS OPTIMIZER:")
print(f"""Number of iterations: {result.nfev}""")
print(f"""Time (s): {execution_time}""")
CHANGED TO BFGS OPTIMIZER:
Number of iterations: 117
Time (s): 0.31656408309936523
مثال VQD
هنا نُطبّق إطار عمل Qiskit patterns بشكل صريح.
الخطوة 1: تحويل المدخلات الكلاسيكية إلى مسألة كمومية
بدلًا من البحث عن أصغر قيمة ذاتية فقط للمرصودات، سنبحث هنا عن جميع القيم الـ (حيث ).
تذكّر أن دوال تكلفة VQD هي:
هذا مهم بشكل خاص لأنه يجب تمرير المتجه (في هذه الحالة ) كوسيطة عند تعريف كائن VQD.
علاوةً على ذلك، في تنفيذ Qiskit لـ VQD، وبدلًا من اعتبار المؤثرات الفعالة المُوضَّحة في الدرس السابق، تُحسَب الإخلاصات مباشرةً عبر خوارزمية ComputeUncompute، التي تستعين بأداة Sampler الأولية لأخذ عينات من احتمال الحصول على للدائرة
. ويعمل هذا بالضبط لأن هذا الاحتمال هو
ansatz = n_local(
num_qubits=2,
rotation_blocks=["ry", "rz"],
entanglement_blocks="cz",
# entanglement="linear",
reps=1,
)
ansatz.decompose().draw("mpl")

لنبدأ بفحص المرصود التالي:
لهذا المرصود القيم الذاتية التالية:
والحالات الذاتية:
from qiskit.quantum_info import SparsePauliOp
observable_2 = SparsePauliOp.from_list([("II", 2), ("XX", -3), ("YY", 2), ("ZZ", -4)])
سنستخدم الدالة التالية لحساب عقوبة التداخل. لاحظ أن هذا لا يزال جزءًا من تحويل المسألة إلى دوائر كمومية. غير أنه، كما نوقش في الدرس السابق، تحسب هذه الدالة التداخل بين الدائرة المتغيرة الحالية والدائرة المُحسَّنة من حالة سابقة ذات طاقة/تكلفة أدنى. كما يجب ترجمة (transpile) الدائرة الجديدة المُولَّدة للتشغيل على عتاد حقيقي. رأينا هذه الدالة من قبل عند استخدام محاكي. هنا، يجب أن نأخذ بعين الاعتبار مسبقًا عملية الترجمة والتحسين المرتبط بها عند استخدام backend حقيقي، ومن هنا جاءت الأسطر المتعلقة بـ if realbackend == 1. هذا يمزج قليلًا بين الخطوتين، لكننا سنُشير إلى الخطوة 2 صراحةً لاحقًا.
import numpy as np
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
def calculate_overlaps(
ansatz, prev_circuits, parameters, sampler, realbackend, backend
):
def create_fidelity_circuit(circuit_1, circuit_2):
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)
if realbackend == 1:
pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
fidelity_circuit = pm.run(fidelity_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. لاحظ أنه مقارنةً بالدرس السابق، لدينا الآن وسيطتان إضافيتان (realbackend وbackend) لمساعدتنا في الترجمة عند استخدام backends حقيقية.
def cost_func_vqd(
parameters,
ansatz,
prev_states,
step,
betas,
estimator,
sampler,
hamiltonian,
realbackend,
backend,
):
estimator_job = estimator.run([(ansatz, hamiltonian, [parameters])])
total_cost = 0
if step > 1:
overlaps = calculate_overlaps(
ansatz, prev_states, parameters, sampler, realbackend, backend
)
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
مجددًا، سنستخدم المحاكيات للتصحيح أولًا، ثم ننتقل إلى العتاد الحقيقي.
from qiskit.primitives import StatevectorSampler
from qiskit.primitives import StatevectorEstimator
sampler = StatevectorSampler(default_shots=4092)
estimator = StatevectorEstimator()
هنا نُحدد عدد الحالات التي نرغب في حسابها، والعقوبات، ومجموعة المعاملات الأولية x0.
from qiskit.quantum_info import SparsePauliOp
k = 4
betas = [50, 60, 40]
x0 = np.ones(8)
سنختبر الآن الخوارزمية باستخدام المحاكيات:
from scipy.optimize import minimize
prev_states = []
prev_opt_parameters = []
eigenvalues = []
realbackend = 0
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_2,
realbackend,
None,
),
method="COBYLA",
options={"maxiter": 200, "tol": 0.000001},
)
print(result)
prev_opt_parameters = result.x
eigenvalues.append(result.fun)
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -6.9999999999996
x: [ 1.571e+00 1.571e+00 2.519e+00 2.100e+00 1.242e+00
6.935e-01 2.298e+00 1.991e+00]
nfev: 151
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: 3.698974255258432
x: [ 1.269e+00 1.109e+00 1.080e+00 1.200e+00 1.094e+00
1.163e+00 9.752e-01 9.519e-01]
nfev: 103
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: 4.731320121938101
x: [ 1.533e+00 2.451e+00 2.526e+00 2.406e+00 1.968e+00
2.105e+00 8.537e-01 8.442e-01]
nfev: 110
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: 7.008239313655201
x: [ 4.150e+00 2.120e+00 3.495e+00 7.262e-01 1.953e+00
-1.982e-01 3.263e-01 2.563e+00]
nfev: 126
maxcv: 0.0
eigenvalues
[np.float64(-6.9999999999996),
np.float64(3.698974255258432),
np.float64(4.731320121938101),
np.float64(7.008239313655201)]
هذه النتائج قريبة إلى حدٍّ ما من النتائج المتوقعة، باستثناء خطأ التقريب والطور العام. يمكننا ضبط مستوى التسامح في المحسِّن الكلاسيكي و/أو عقوبات تداخل متجهات الحالة للحصول على قيم أكثر دقة.
solution_eigenvalues = [-7, 3, 5, 7]
for index, experimental_eigenvalue in enumerate(eigenvalues):
solution_eigenvalue = solution_eigenvalues[index]
print(
f"Percent error: {abs((experimental_eigenvalue - solution_eigenvalue)/solution_eigenvalue):.2e}"
)
Percent error: 5.71e-14
Percent error: 2.33e-01
Percent error: 5.37e-02
Percent error: 1.18e-03
تغيير بيتا
كما ذُكر في الدرس السابق، يجب أن تكون قيم أكبر من الفرق بين القيم الذاتية. لنرَ ماذا يحدث حين لا تستوفي هذا الشرط مع
بقيم ذاتية
from qiskit.quantum_info import SparsePauliOp
k = 4
betas = np.ones(3)
x0 = np.zeros(8)
from scipy.optimize import minimize
prev_states = []
prev_opt_parameters = []
eigenvalues = []
realbackend = 0
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_2,
realbackend,
None,
),
method="COBYLA",
options={"tol": 0.01, "maxiter": 200},
)
print(result)
prev_opt_parameters = result.x
eigenvalues.append(result.fun)
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -6.999916534745094
x: [ 1.568e+00 -1.569e+00 1.385e-01 1.398e-01 -7.972e-01
7.835e-01 -2.375e-01 4.539e-02]
nfev: 125
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -1.515139929812874
x: [-5.317e-04 -2.514e-03 1.016e+00 9.998e-01 3.890e-04
1.772e-04 1.568e-04 8.497e-04]
nfev: 35
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: -0.509948114293115
x: [-3.796e-03 8.853e-03 3.015e-04 9.997e-01 6.271e-04
-2.554e-03 1.017e-04 2.766e-04]
nfev: 37
maxcv: 0.0
message: Return from COBYLA because the trust region radius reaches its lower bound.
success: True
status: 0
fun: 0.4914672235935682
x: [-7.178e-03 -8.652e-03 1.125e+00 -5.428e-02 -1.586e-03
2.031e-03 -3.462e-03 5.734e-03]
nfev: 35
maxcv: 0.0
solution_eigenvalues = [-7, 3, 5, 7]
for index, experimental_eigenvalue in enumerate(eigenvalues):
solution_eigenvalue = solution_eigenvalues[index]
print(
f"Percent error: {abs((experimental_eigenvalue - solution_eigenvalue)/solution_eigenvalue):.2e}"
)
Percent error: 1.19e-05
Percent error: 1.51e+00
Percent error: 1.10e+00
Percent error: 9.30e-01
هذه المرة، يُعيد المحسِّن نفس الحالة باعتبارها الحل المقترح لجميع الحالات الذاتية: وهذا خطأ واضح. يحدث هذا لأن قيم betas كانت صغيرة جدًا بحيث لم تُعاقب الحالة الذاتية الدنيا في دوال التكلفة المتعاقبة. لذلك، لم تُستبعد من فضاء البحث الفعال في التكرارات اللاحقة للخوارزمية، وظلت تُختار دائمًا كأفضل حل ممكن.
نوصي بتجربة قيم مختلفة لـ ، والتأكد من أنها أكبر من الفرق بين القيم الذاتية.
الخطوة 2: تحسين المسألة للتنفيذ الكمومي
لتشغيل هذا على عتاد حقيقي، يجب علينا تحسين الدوائر الكمومية لحاسوبنا الكمومي المختار. لأغراضنا هنا، سنستخدم ببساطة الـ backend الأقل انشغالًا.
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)
# Or use a specific backend
# backend = service.backend("ibm_brisbane")
print(backend)
<IBMBackend('ibm_brisbane')>
سنُترجم دائرتنا باستخدام مُدير مسار جاهز ومستوى تحسين 3.
pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_ansatz = pm.run(ansatz)
isa_observable = observable_2.apply_layout(layout=isa_ansatz.layout)
الخطوة 3: التنفيذ باستخدام Qiskit primitives
بعد التأكد من إعادة ضبط قيم betas لتكون مرتفعة بما يكفي، يمكننا الآن تشغيل حسابنا على عتاد كمومي حقيقي.
# Estimated compute resource usage: 25 minutes.
# Benchmarked at 24 min, 30 sec on an Eagle r3 processor on 5-30-24
k = 2
betas = [30, 50, 80]
x0 = np.zeros(8)
real_prev_states = []
real_prev_opt_parameters = []
real_eigenvalues = []
realbackend = 1
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)
with Session(backend=backend) as session:
estimator = Estimator(mode=session, options=estimator_options)
sampler = Sampler(mode=session)
for step in range(1, k + 1):
if step > 1:
real_prev_states.append(isa_ansatz.assign_parameters(prev_opt_parameters))
result = minimize(
cost_func_vqd,
x0,
args=(
isa_ansatz,
real_prev_states,
step,
betas,
estimator,
sampler,
isa_observable,
realbackend,
backend,
),
method="COBYLA",
options={"maxiter": 200},
)
print(result)
real_prev_opt_parameters = result.x
real_eigenvalues.append(result.fun)
session.close()
print(real_eigenvalues)
الخطوة 4: المعالجة اللاحقة، وإعادة النتيجة بصيغة كلاسيكية
مخرجاتنا مشابهة هيكليًا لما نوقش في الدروس والأمثلة السابقة. لكن ثمة إشكالية في النتائج أعلاه يمكننا استخلاص منها درس تحذيري في سياق الحالات المُثارة. لتقليل وقت الحساب المستخدم في هذا المثال التعليمي، حددنا حدًا أقصى لعدد تكرارات المحسِّن الكلاسيكي كان منخفضًا بشكل محتمل: 200 تكرار. حساب سابق على محاكٍ لم يتقارب في 200 تكرار. هنا، تقارب الحساب لدينا... لكن إلى أي مستوى من الدقة؟ لم نُحدد مستوى تسامح لـ COBYLA ليعتبر نفسه قد "تقارب". نظرة سريعة على قيمة الدالة ومقارنتها بعمليات تشغيل سابقة تُخبرنا بأن COBYLA لم يكن قريبًا من التقارب إلى الدقة المطلوبة.
ثمة إشكالية أخرى: تبدو طاقة الحالة المُثارة الأولى أقل من طاقة حالة الأساس! انظر إن كنت تستطيع تفسير كيف يمكن حدوث ذلك. تلميح: يرتبط بنقطة التقارب التي تناولناها للتو. يُشرح هذا السلوك بالتفصيل أدناه بعد تطبيق VQD على جزيء H2.
الكيمياء الكمية: حل حالة الطاقة الأساسية والمثارة
هدفنا هو تصغير قيمة التوقع للمرصود الذي يمثّل الطاقة (هاميلتوني ):
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.library import efficient_su2
H2_op = SparsePauliOp.from_list(
[
("II", -1.052373245772859),
("IZ", 0.39793742484318045),
("ZI", -0.39793742484318045),
("ZZ", -0.01128010425623538),
("XX", 0.18093119978423156),
]
)
chem_ansatz = efficient_su2(H2_op.num_qubits)
chem_ansatz.decompose().draw("mpl")

from qiskit import QuantumCircuit
def cost_func_vqe(params, 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
"""
pub = (ansatz, hamiltonian, params)
cost = estimator.run([pub]).result()[0].data.evs
# cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
return cost
نضع الآن مجموعة أولية من المعاملات:
import numpy as np
x0 = np.ones(chem_ansatz.num_parameters)
يمكننا تصغير دالة التكلفة هذه للحصول على المعاملات المثلى، ويمكننا التحقق من الكود أولاً باستخدام محاكٍ محلي.
from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler
estimator = Estimator()
sampler = Sampler()
# SciPy minimizer routine
from scipy.optimize import minimize
import time
start_time = time.time()
result = minimize(
cost_func_vqe, x0, args=(chem_ansatz, H2_op, estimator), method="COBYLA"
)
end_time = time.time()
execution_time = end_time - start_time
result
message: Optimization terminated successfully.
success: True
status: 1
fun: -1.857275029048451
x: [ 7.326e-01 1.354e+00 ... 1.040e+00 1.508e+00]
nfev: 242
maxcv: 0.0
القيمة الدنيا لدالة التكلفة (-1.857...) هي طاقة الحالة الأساسية لجزيء H2، بوحدة هارتري.
الحالات المثارة
يمكننا أيضاً الاستفادة من VQD لحل حالة إجمالية (الحالة الأساسية والحالة المثارة الأولى).
from qiskit.quantum_info import SparsePauliOp
import numpy as np
k = 2
betas = [33, 33]
# x0 = np.zeros(ansatz.num_parameters)
x0 = [
1.164e00,
-2.438e-01,
9.358e-04,
6.745e-02,
1.990e00,
9.810e-02,
6.154e-01,
5.454e-01,
]
سنضيف حساب التداخل الخاص بنا:
from scipy.optimize import minimize
prev_states = []
prev_opt_parameters = []
eigenvalues = []
realbackend = 0
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,
H2_op,
realbackend,
None,
),
method="COBYLA",
options={"tol": 0.001, "maxiter": 2000},
)
print(result)
prev_opt_parameters = result.x
eigenvalues.append(result.fun)
message: Optimization terminated successfully.
success: True
status: 1
fun: -1.8572671093941977
x: [ 1.164e+00 -2.437e-01 2.118e-03 6.448e-02 1.990e+00
9.870e-02 6.167e-01 5.476e-01]
nfev: 58
maxcv: 0.0
message: Optimization terminated successfully.
success: True
status: 1
fun: -1.0322873777662176
x: [ 3.205e+00 1.502e+00 1.699e+00 -1.107e-02 3.086e+00
1.530e+00 4.445e-02 7.013e-02]
nfev: 99
maxcv: 0.0
eigenvalues
[-1.8572671093941977, -1.0322873777662176]
العتاد الحقيقي وتحذير ختامي
لتشغيل هذا على عتاد حقيقي، يجب علينا تحسين الدوائر الكمية لجهاز الحوسبة الكمية الذي اخترناه. لأغراضنا هنا، سنستخدم ببساطة الـ backend الأقل انشغالًا.
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)
سنستخدم مدير تمريرات جاهزاً للترجمة، وسنُحسّن دائرتنا إلى أقصى حد باستخدام مستوى التحسين 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 = H2_op.apply_layout(layout=isa_ansatz.layout)
نظراً لأن VQD يعتمد اعتماداً كبيراً على التكرار، سنُنجز جميع الخطوات داخل جلسة Runtime، بحيث تُضاف وظائفنا إلى قائمة الانتظار في البداية فقط، وليس بين كل تحديث للمعاملات. لا يتغير شيء آخر في صياغة دالة التكلفة أو المُقدِّر.
x0 = [
1.306e00,
-2.284e-01,
6.913e-02,
-2.530e-02,
1.849e00,
7.433e-02,
6.366e-01,
5.600e-01,
]
# Estimated hardware usage: 20 min benchmarked on an Eagle r3 processor on 5-30-24
real_prev_states = []
real_prev_opt_parameters = []
real_eigenvalues = []
realbackend = 1
estimator_options = EstimatorOptions(resilience_level=1, default_shots=4096)
with Session(backend=backend) as session:
estimator = Estimator(mode=session)
sampler = Sampler(mode=session)
for step in range(1, k + 1):
if step > 1:
real_prev_states.append(
isa_ansatz.assign_parameters(real_prev_opt_parameters)
)
result = minimize(
cost_func_vqd,
x0,
args=(
isa_ansatz,
real_prev_states,
step,
betas,
estimator,
sampler,
isa_observable,
realbackend,
backend,
),
method="COBYLA",
options={"tol": 0.001, "maxiter": 300},
)
print(result)
real_prev_opt_parameters = result.x
real_eigenvalues.append(result.fun)
session.close()
print(real_eigenvalues)
طاقة الحالة الأساسية المحصّل عليها (-1.83 هارتري) ليست بعيدة كثيراً عن القيمة الصحيحة (-1.85 هارتري). غير أن طاقة الحالة المثارة بعيدة جداً عن القيمة الصحيحة. وهذا يشبه السلوك الخاطئ الذي رأيناه سابقاً في هذا الدرس. فالطاقة المُبلَّغ عنها للحالة المثارة تكاد تطابق طاقة الحالة الأساسية. وفي الحالة السابقة، رأينا حتى طاقة حالة مثارة كانت أقل من طاقة الحالة الأساسية المُبلَّغ عنها.
لا يمكن للحساب التغايري أن يُعطي طاقةً أقل من طاقة الحالة الأساسية الحقيقية. وفي الحالة السابقة، لم تكن طاقة الحالة الأساسية التي حصلنا عليها قريبة جداً من الطاقة الأساسية الحقيقية. وبما أننا لم نحصل على طاقة الحالة الأساسية الحقيقية في تلك الحالة، فلا يوجد تناقض. أما في الحالة الراهنة، فقد كانت طاقة الحالة الأساسية قريبةً إلى حدٍّ ما من القيمة الصحيحة، ومع ذلك يبدو أن طاقة الحالة المثارة تقترب بشكل غريب من تلك القيمة ذاتها.
لكي نفهم كيف حدث هذا، لنتذكّر أن طريقة إيجاد الحالة المثارة تقوم على اشتراط أن تكون الحالة التغايرية متعامدة مع الحالة الأساسية (باستخدام دوائر التداخل وعبارات العقوبة). فإذا فشلنا في الحصول على طاقة حالة أساسية دقيقة (أو أخطأنا بنسبة مئوية بسيطة)، فإننا نفشل أيضاً في الحصول على متجه الحالة الأساسية الدقيق! فعندما نشترط أن تكون الحالة المثارة متعامدة مع الحالة الأولى التي وجدناها، لم نكن نفرض التعامد مع الحالة الأساسية الحقيقية، بل مع تقريب لها (وأحياناً تقريب رديء). وهكذا، لم تُضطر الحالة المثارة إلى أن تكون متعامدة مع الحالة الأساسية الحقيقية، وكانت تقديرات طاقتنا للحالات المثارة قريبةً فعلاً من طاقة الحالة الأساسية.
سيظل هذا مصدر قلق دائم في VQD. لكن يمكن من حيث المبدأ تصحيح ذلك بزيادة الحد الأقصى لعدد التكرارات في المُحسِّن الكلاسيكي، وفرض تفاوت أدنى للمُحسِّن الكلاسيكي، وربما تجربة ansatz مختلف إذا كنا نخطئ باستمرار في الوصول إلى الحالة الأساسية الحقيقية. وكما رأينا، قد يلزم أيضاً تعديل عقوبات التداخل (betas). لكن هذه مسألة مختلفة في الواقع. لن تُبعدك أي عقوبة على التداخل عن الحالة الأساسية الحقيقية، إذا لم تجد تقديراً جيداً للحالة الأساسية الحقيقية في دائرة التداخل.
التحسين: أقصى قطع (Max-Cut)
مسألة الحد الأقصى للقطع (Max-Cut) هي مسألة تحسين تركيبي تتمثل في تقسيم رؤوس رسم بياني إلى مجموعتين منفصلتين بحيث يكون عدد الحواف بين المجموعتين أقصى ما يمكن. بصياغة أكثر دقة، بافتراض رسم بياني غير موجَّه ، حيث هي مجموعة الرؤوس و هي مجموعة الحواف، تطلب منا مسألة Max-Cut تقسيم الرؤوس إلى مجموعتين فرعيتين منفصلتين، و، بحيث يتعظَّم عدد الحواف التي طرفها الأول في والآخر في .
يمكن تطبيق Max-Cut لحل مسائل متنوعة كالتجميع (clustering)، وتصميم الشبكات، وانتقالات الطور. سنبدأ بإنشاء رسم بياني للمسألة:
import rustworkx as rx
from rustworkx.visualization import mpl_draw
n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)
mpl_draw(
G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)
يمكن التعبير عن هذه المسألة بوصفها مسألة تحسين ثنائي. لكل عقدة ، حيث هو عدد عقد الرسم البياني (في هذه الحالة )، نأخذ المتغير الثنائي . قيمة هذا المتغير هي إذا كانت العقدة في المجموعة التي نسميها ، و إذا كانت في المجموعة الأخرى التي نسميها . كذلك نُعرِّف (العنصر من مصفوفة التجاور ) بوصفه وزن الحافة من العقدة إلى العقدة . وبما أن الرسم البياني غير موجَّه، فإن . يمكننا حينئذٍ صياغة مسألتنا على أنها تعظيم دالة التكلفة التالية:
لحل هذه المسألة بحاسوب كمي، سنُعبِّر عن دالة التكلفة على شكل القيمة المتوقعة للمرصود. غير أن المرصودات التي يقبلها Qiskit بشكل أصلي تتكون من مؤثرات باولي، التي تمتلك القيم الذاتية و بدلاً من و. لهذا السبب سنجري تغيير المتغير التالي:
حيث . يمكننا استخدام مصفوفة التجاور للوصول المريح إلى أوزان جميع الحواف، وسيُستخدَم ذلك للحصول على دالة التكلفة:
وهذا يعني:
إذن دالة التكلفة الجديدة التي نريد تعظيمها هي:
علاوة على ذلك، الحاسوب الكمي بطبيعته يجد الحدود الدنيا (عادةً أدنى طاقة) لا القصوى، لذا بدلاً من تعظيم سنعمل على تصغير:
الآن بعد أن حصلنا على دالة تكلفة للتصغير تأخذ متغيراتها القيمتين و، يمكننا إجراء القياس التالي مع مؤثر باولي :
بعبارة أخرى، المتغير سيكون مكافئاً لبوابة تعمل على الكيوبت . وعلاوة على ذلك:
وبالتالي المرصود الذي سنأخذه بالاعتبار هو:
ويجب إضافة الحد المستقل لاحقاً:
المؤثر هو مجموعة خطية من الحدود التي تحوي مؤثرات Z على العقد المرتبطة بحافة (تذكر أن الكيوبت رقم 0 هو الأبعد إلى اليمين): . بعد بناء المؤثر، يمكن بناء دارة التهيئة (ansatz) لخوارزمية QAOA بسهولة باستخدام دارة QAOAAnsatz من مكتبة دارات Qiskit.
from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp
max_hamiltonian = SparsePauliOp.from_list(
[("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)
max_ansatz = QAOAAnsatz(max_hamiltonian, reps=2)
# Draw
max_ansatz.decompose(reps=3).draw("mpl")
# Sum the weights, and divide by 2
offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")
Offset: -2.5
def cost_func(params, 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
"""
pub = (ansatz, hamiltonian, params)
cost = estimator.run([pub]).result()[0].data.evs
# cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
return cost
from qiskit.primitives import StatevectorEstimator as Estimator
from qiskit.primitives import StatevectorSampler as Sampler
estimator = Estimator()
sampler = Sampler()
نُحدِّد الآن مجموعة أولية من المعاملات العشوائية:
import numpy as np
x0 = 2 * np.pi * np.random.rand(max_ansatz.num_parameters)
print(x0)
[6.0252949 0.58448176 2.15785731 1.13646074]
يمكن استخدام أي مُحسِّن كلاسيكي لتصغير دالة التكلفة. على نظام كمي حقيقي، عادةً ما يُعطي المُحسِّن المصمَّم للتعامل مع مشاهد دالة التكلفة غير الناعمة نتائج أفضل. هنا نستخدم روتين COBYLA من SciPy عبر دالة minimize.
بما أننا ننفِّذ استدعاءات متكررة كثيرة إلى Runtime، نستخدم جلسة (session) لتنفيذ جميع الاستدعاءات ضمن كتلة واحدة. علاوة على ذلك، في QAOA يكون الحل مُشفَّراً في توزيع المخرجات لدارة التهيئة المرتبطة بالمعاملات المثلى الناتجة عن التصغير. لذلك سنحتاج إلى أداة Sampler الأولية، ونُنشئ مثيلاً منها مع الجلسة نفسها، ثم نُشغِّل روتين التصغير:
result = minimize(
cost_func, x0, args=(max_ansatz, max_hamiltonian, estimator), method="COBYLA"
)
print(result)
message: Optimization terminated successfully.
success: True
status: 1
fun: -2.585287311689236
x: [ 7.332e+00 3.904e-01 2.045e+00 1.028e+00]
nfev: 80
maxcv: 0.0
متجه الحل من زوايا المعاملات (x)، عند تعويضه في دارة التهيئة، يعطينا تقسيم الرسم البياني الذي كنا نبحث عنه.
eigenvalue = cost_func(result.x, max_ansatz, max_hamiltonian, estimator)
print(f"""Eigenvalue: {eigenvalue}""")
print(f"""Max-Cut Objective: {eigenvalue + offset}""")
Eigenvalue: -2.585287311689236
Max-Cut Objective: -5.085287311689235
from qiskit.result import QuasiDistribution
from qiskit.primitives import StatevectorSampler
sampler = StatevectorSampler()
# Assign solution parameters to ansatz
qc = max_ansatz.assign_parameters(result.x)
# Add measurements to our circuit
qc.measure_all()
# Sample ansatz at optimal parameters
# samp_dist = sampler.run(qc).result().quasi_dists[0]
shots = 1024
job = sampler.run([qc], shots=shots)
qc.decompose().draw("mpl")
data_pub = job.result()[0].data
bitstrings = data_pub.meas.get_bitstrings()
counts = data_pub.meas.get_counts()
quasi_dist = QuasiDistribution(
{outcome: freq / shots for outcome, freq in counts.items()}
)
probabilities = quasi_dist
# Close the session since we are now done with it
# session.close()
from qiskit.visualization import plot_distribution
plot_distribution(counts)
binary_string = max(counts.items(), key=lambda kv: kv[1])[0]
x = np.asarray([int(y) for y in reversed(list(binary_string))])
colors = ["r" if x[i] == 0 else "c" for i in range(n)]
mpl_draw(
G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color=colors
)
الملخص
من خلال هذا الدرس، تعلمت:
- كيفية كتابة خوارزمية تغايرية مخصصة
- كيفية تطبيق خوارزمية تغايرية لإيجاد أدنى القيم الذاتية
- كيفية الاستفادة من الخوارزميات التغايرية لحل حالات تطبيقية عملية