صيغ متعددة الحدود لتقليل خطأ Trotter

الاستخدام التقديري لوحدة المعالجة الكمومية (QPU): أربع دقائق على معالج Heron r2 (ملاحظة: هذا تقدير فحسب. قد يختلف وقت التشغيل الفعلي.)

الخلفية

يستعرض هذا البرنامج التعليمي كيفية استخدام صيغة متعددة الحدود (MPF) لتحقيق خطأ Trotter أقل في الرصدة مقارنةً بالخطأ الناجم عن أعمق دائرة Trotter التي سننفذها فعلياً. تُقلِّل صيغ MPF من خطأ Trotter في ديناميكيات هاملتوني من خلال تركيبة موزونة من عدة عمليات تنفيذ للدوائر. نتأمل مهمة إيجاد القيم المتوقعة للرصدات للحالة الكمومية $\rho(t)=e^{-i H t} \rho(0) e^{i H t}$ مع الهاملتوني $H$. يمكن استخدام صيغ الضرب (PFs) لتقريب عملية التطور الزمني $e^{-i H t}$ بالطريقة الآتية:

• كتابة الهاملتوني $H$ على الصورة $H=\sum_{a=1}^d F_a,$ حيث $F_a$ مؤثرات Hermitian بحيث يمكن تنفيذ كل وحدة موحدة مقابلة بكفاءة على الجهاز الكمومي.
• تقريب الحدود $F_a$ التي لا تتبادل التبديل مع بعضها.

وعليه، تكون صيغة الضرب من الرتبة الأولى (صيغة Lie-Trotter) على النحو:

$S_1(t):=\prod_{a=1}^d e^{-i F_a t},$

وهي تمتلك حد خطأ تربيعي $S_1(t)=e^{-i H t}+\mathcal{O}\left(t^{2}\right)$. يمكن أيضاً استخدام صيغ الضرب ذات الرتب الأعلى (صيغ Lie-Trotter-Suzuki)، التي تتقارب بصورة أسرع، وتُعرَّف بصورة تكرارية على النحو:

$S_2(t):=\prod_{a=1}^d e^{-i F_a t/2}\prod_{a=1}^d e^{-i F_a t/2}$

$S_{2 \chi}(t):= S_{2 \chi -2}(s_{\chi}t)^2 S_{2 \chi -2}((1-4s_{\chi})t)S_{2 \chi -2}(s_{\chi}t)^2,$

حيث $\chi$ هو رتبة صيغة الضرب المتماثلة و $s_p = \left( 4 - 4^{1/(2p-1)} \right)^{-1}$. في حالة التطورات الزمنية الطويلة، يمكن تقسيم الفترة الزمنية $t$ إلى $k$ فترات تُسمى خطوات Trotter، مدة كل منها $t/k$، وتقريب التطور الزمني في كل فترة بصيغة ضرب من الرتبة $\chi$ وهي $S_{\chi}$. وبذلك، تكون صيغة الضرب من الرتبة $\chi$ لمؤثر التطور الزمني عبر $k$ خطوة Trotter:

$S_{\chi}^{k}(t) = \left[ S_{\chi} \left( \frac{t}{k} \right)\right]^k = e^{-i H t}+O\left(t \left( \frac{t}{k} \right)^{\chi} \right)$

حيث يتناقص حد الخطأ مع زيادة عدد خطوات Trotter $k$ ورتبة $\chi$ لصيغة الضرب.

بالنظر إلى عدد صحيح $k \geq 1$ وصيغة ضرب $S_{\chi}(t)$، يمكن الحصول على الحالة المتطورة زمنياً التقريبية $\rho_k(t)$ من $\rho_0$ بتطبيق $k$ تكراراً من صيغة الضرب $S_{\chi}\left(\frac{t}{k}\right)$.

$$\rho_k(t)=S_{\chi}\left(\frac{t}{k}\right)^k \rho_0 S_{\chi}\left(\frac{t}{k}\right)^{-k}$$

$\rho_k(t)$ هو تقريب لـ $\rho(t)$ بخطأ تقريب Trotter ||$\rho_k(t)-\rho(t) ||$. إذا اعتبرنا تركيبة خطية من تقريبات Trotter لـ $\rho(t)$:

$$\mu(t) = \sum_{j}^{l} x_j \rho^{k_j}_{j}\left(\frac{t}{k_j}\right) + \text{some remaining Trotter error} \, ,$$

حيث $x_j$ هي معاملات الترجيح لدينا، و $\rho^{k_j}_j$ هي مصفوفة الكثافة المقابلة للحالة النقية المتحصل عليها بتطوير الحالة الابتدائية بصيغة الضرب $S^{k_j}_{\chi}$ المتضمنة $k_j$ خطوة Trotter، و $j \in {1, ..., l}$ يُفهرس عدد صيغ الضرب التي تتكون منها صيغة MPF. تستخدم جميع الحدود في $\mu(t)$ صيغة الضرب ذاتها $S_{\chi}(t)$ أساساً لها. الهدف هو التحسين على ||$\rho_k(t)-\rho(t) \|$ بإيجاد $\mu(t)$ ذات $\|\mu(t)-\rho(t)\|$ أصغر.

• لا يشترط أن يكون $\mu(t)$ حالة فيزيائية لأن $x_i$ لا يلزم أن تكون موجبة. الهدف هنا هو تقليل الخطأ في القيمة المتوقعة للرصدات وليس إيجاد بديل فيزيائي لـ $\rho(t)$.
• يحدد $k_j$ كلاً من عمق الدائرة ومستوى تقريب Trotter. تؤدي القيم الأصغر لـ $k_j$ إلى دوائر أقصر، مما يُفضي إلى أخطاء أقل في الدوائر لكنه يجعلها تقريباً أقل دقة للحالة المطلوبة.

المحور الأساسي هنا هو أن خطأ Trotter المتبقي الذي تعطيه $\mu(t)$ أصغر من خطأ Trotter الذي يمكن الحصول عليه باستخدام أكبر قيمة لـ $k_j$ مباشرةً.

يمكن النظر إلى فائدة هذا النهج من منظورين:

1. بالنسبة لميزانية ثابتة من خطوات Trotter القابلة للتنفيذ، يمكن الحصول على نتائج بخطأ Trotter إجمالي أصغر.
2. بالنظر إلى عدد مستهدف من خطوات Trotter يفوق الطاقة التنفيذية، يمكن استخدام صيغة MPF للعثور على مجموعة من الدوائر ذات العمق الأقل لتشغيلها بما يؤدي إلى خطأ Trotter مماثل.

المتطلبات

قبل البدء في هذا البرنامج التعليمي، تأكد من تثبيت ما يلي:

• Qiskit SDK الإصدار v1.0 أو أحدث، مع دعم التصور
• Qiskit Runtime الإصدار v0.22 أو أحدث (pip install qiskit-ibm-runtime)
• إضافات MPF لـ Qiskit (pip install qiskit_addon_mpf)
• أدوات إضافات Qiskit (pip install qiskit_addon_utils)
• مكتبة Quimb (pip install quimb)
• مكتبة Qiskit Quimb (pip install qiskit-quimb)
• Numpy الإصدار v0.21 لضمان التوافق بين الحزم (pip install numpy==0.21)

الجزء الأول. مثال على نطاق صغير

استكشاف استقرار صيغة MPF

لا يوجد قيد واضح على اختيار عدد خطوات Trotter $k_j$ التي تشكّل حالة MPF وهي $\mu(t)$. غير أنه يجب اختيارها بعناية لتجنب عدم الاستقرار في قيم التوقع المحسوبة من $\mu(t)$. قاعدة عامة جيدة هي تعيين أصغر خطوة Trotter $k_{\text{min}}$ بحيث يكون $t/k_{\text{min}} \lt 1$. إذا أردت معرفة المزيد عن هذا الموضوع وكيفية اختيار قيم $k_j$ الأخرى، راجع دليل كيفية اختيار خطوات Trotter لصيغة MPF.

في المثال أدناه، نستكشف استقرار حل MPF بحساب القيمة المتوقعة للمغنطة لنطاق زمني باستخدام حالات متطورة زمنياً مختلفة. تحديداً، نقارن القيم المتوقعة المحسوبة من كل تطور زمني تقريبي مُنفَّذ بخطوات Trotter المقابلة له، ومن نماذج MPF المختلفة (المعاملات الثابتة والديناميكية)، بالقيم الدقيقة للرصدة المتطورة زمنياً. أولاً، لنحدد معلمات صيغ Trotter وأوقات التطور

# Added by doQumentation — installs packages not in the Binder environment
%pip install -q qiskit-addon-mpf

import numpy as np

mpf_trotter_steps = [1, 2, 4]
order = 2
symmetric = False

trotter_times = np.arange(0.5, 1.55, 0.1)
exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)

في هذا المثال سنستخدم حالة Neel كالحالة الابتدائية $\vert \text{Neel} \rangle = \vert 0101...01 \rangle$ ونموذج Heisenberg على سلسلة مكوّنة من 10 مواقع باعتباره الهاملتوني المحكوم للتطور الزمني

$$\hat{\mathcal{H}}_{Heis} = J \sum_{i=1}^{L-1} \left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ Z_i Z_{(i+1)} \right) \, ,$$

حيث $J$ هو ثابت الاقتران للحواف المجاورة.

from qiskit.transpiler import CouplingMap
from rustworkx.visualization import graphviz_draw
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian
import numpy as np

L = 10

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

# Get a qubit operator describing the Heisenberg field model
hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(1.0, 1.0, 1.0),
    ext_magnetic_field=(0.0, 0.0, 0.0),
)

print(hamiltonian)

SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

الرصدة التي سنقيسها هي المغنطة على زوج من الكيوبتات في منتصف السلسلة.

from qiskit.quantum_info import SparsePauliOp

observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)

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

نُعرِّف تمريرة transpiler لجمع دورات XX وYY في الدائرة في بوابة XX+YY واحدة. سيتيح لنا هذا الاستفادة من خصائص حفظ الدوران في TeNPy أثناء حساب MPO، مما يُسرِّع الحساب تسريعاً ملحوظاً.

from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler import PassManager
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)
from functools import partial

def filter_function(node):
    return node.op.name in {"rxx", "ryy"}

collect_function = partial(
    collect_using_filter_function,
    filter_function=filter_function,
    split_blocks=True,
    min_block_size=1,
)

def collapse_to_xx_plus_yy(block):
    param = 0.0
    for node in block.data:
        param += node.operation.params[0]
    return XXPlusYYGate(param)

collapse_function = partial(
    collapse_to_operation,
    collapse_function=collapse_to_xx_plus_yy,
)

pm = PassManager()
pm.append(CollectAndCollapse(collect_function, collapse_function))

ثم ننشئ الدوائر التي تُنفِّذ التطورات الزمنية التقريبية بصيغة Trotter.

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

# Initial Neel state preparation
initial_state_circ = QuantumCircuit(L)
initial_state_circ.x([i for i in range(L) if i % 2 != 0])

all_circs = []
for total_time in trotter_times:
    mpf_trotter_circs = [
        generate_time_evolution_circuit(
            hamiltonian,
            time=total_time,
            synthesis=SuzukiTrotter(reps=num_steps, order=order),
        )
        for num_steps in mpf_trotter_steps
    ]

    mpf_trotter_circs = pm.run(
        mpf_trotter_circs
    )  # Collect XX and YY into XX + YY

    mpf_circuits = [
        initial_state_circ.compose(circuit) for circuit in mpf_trotter_circs
    ]
    all_circs.append(mpf_circuits)

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

بعد ذلك، نحسب قيم التوقع المتطورة زمنياً من دوائر Trotter.

from copy import deepcopy
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

aer_sim = AerSimulator()
estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for t, mpf_circuits in zip(trotter_times, all_circs):
    mpf_expvals = []
    circuits = [deepcopy(circuit) for circuit in mpf_circuits]
    pm_sim = generate_preset_pass_manager(
        backend=aer_sim, optimization_level=3
    )
    isa_circuits = pm_sim.run(circuits)
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals = [res.data.evs for res in result]
    mpf_stds = [res.data.stds for res in result]
    mpf_expvals_all_times.append(mpf_expvals)
    mpf_stds_all_times.append(mpf_stds)

نحسب أيضاً القيم المتوقعة الدقيقة للمقارنة.

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

exact_expvals = []
for t in exact_evolution_times:
    # Exact expectation values
    exp_H = expm(-1j * t * hamiltonian.to_matrix())
    initial_state = Statevector(initial_state_circ).data
    time_evolved_state = exp_H @ initial_state

    exact_obs = (
        time_evolved_state.conj()
        @ observable.to_matrix()
        @ time_evolved_state
    ).real
    exact_expvals.append(exact_obs)

معاملات MPF الثابتة

صيغ MPF الثابتة هي تلك التي لا تعتمد فيها قيم $x_j$ على وقت التطور $t$. لنأخذ صيغة الضرب من الرتبة $\chi = 1$ مع $k_j$ خطوة Trotter، والتي يمكن كتابتها على الصورة:

$S_1^{k_j}\left( \frac{t}{k_j} \right)=e^{-i H t}+ \sum_{n=1}^{\infty} A_n \frac{t^{n+1}}{k_j^n}$

حيث $A_n$ مصفوفات تعتمد على المبدِّلات (commutators) بين حدود $F_a$ في تحليل الهاملتوني. من المهم الإشارة إلى أن $A_n$ نفسها مستقلة عن الوقت وعدد خطوات Trotter $k_j$. لذلك، يمكن إلغاء حدود الخطأ ذات الرتب الدنيا المساهمة في $\mu(t)$ باختيار دقيق للأوزان $x_j$ للتركيبة الخطية. لإلغاء خطأ Trotter لأول $l-1$ حد (وهي ذات الإسهام الأكبر لأنها تقابل عدداً أقل من خطوات Trotter) في تعبير $\mu(t)$، يجب أن تحقق المعاملات $x_j$ المعادلات الآتية:

$\sum_{j=1}^l x_j = 1$ $\sum_{j=1}^{l-1} \frac{x_j}{k_j^{n}} = 0$

حيث $n=0, ... l-2$. تضمن المعادلة الأولى انعدام التحيز في الحالة المُنشأة $\mu(t)$، بينما تكفل المعادلة الثانية إلغاء أخطاء Trotter. في حالة صيغ الضرب ذات الرتب الأعلى، تصبح المعادلة الثانية $\sum_{j=1}^{l-1} \frac{x_j}{k_j^{\eta}} = 0$ حيث $\eta = \chi + 2n$ لصيغ الضرب المتماثلة و $\eta = \chi + n$ لغيرها، مع $n=0, ..., l-2$. ويكون الخطأ الناتج (المراجع [1],[2]) حينئذٍ:

$\epsilon = \mathcal{O} \left( \frac{t^{l+1}}{k_1^l} \right).$

يتمثل تحديد معاملات MPF الثابتة لمجموعة من قيم $k_j$ في حل النظام الخطي من المعادلات المعرَّف بالمعادلتين أعلاه بالنسبة للمتغيرات $x_j$: $Ax=b$. حيث $x$ هي معاملاتنا المطلوبة، و $A$ مصفوفة تعتمد على $k_j$ ونوع صيغة الضرب المستخدمة ($S$)، و $b$ شعاع من القيود. تحديداً:

$A_{0,j} = 1$ $A_{i>0,j} = k_{j}^{-(\chi + s(i-1))}$ $b_0 = 1$ $b_{i>0} = 0$

حيث $\chi$ هو order، و $s$ يساوي $2$ إذا كان symmetric هو True وإلا فيساوي $1$، و $k_{j}$ هي trotter_steps، و $x$ هي المتغيرات المراد إيجادها. تبدأ المؤشرات $i$ و $j$ من $0$. يمكننا أيضاً تمثيل هذا في الصورة المصفوفية:

$$A = \begin{bmatrix} A_{0,0} & A_{0,1} & A_{0,2} & ... \\ A_{1,0} & A_{1,1} & A_{1,2} & ... \\ A_{2,0} & A_{2,1} & A_{2,2} & ... \\ ... & ... & ... & ... \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & ... \\ k_{0}^{-(\chi + s(1-1))} & k_{1}^{-(\chi + s(1-1))} & k_{2}^{-(\chi + s(1-1))} & ... \\ k_{0}^{-(\chi + s(2-1))} & k_{1}^{-(\chi + s(2-1))} & k_{2}^{-(\chi + s(2-1))} & ... \\ ... & ... & ... & ... \end{bmatrix}$$

و

$$b = \begin{bmatrix} b_{0} \\ b_{1} \\ b_{2} \\ ... \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ ... \end{bmatrix}$$

لمزيد من التفاصيل، راجع توثيق نظام المعادلات الخطية (LSE).

يمكن إيجاد حل تحليلي لـ $x$ بالصورة $x = A^{-1}b$؛ انظر مثلاً المراجع [1] أو [2]. غير أن هذا الحل الدقيق قد يكون "سيئ التهيئة"، مما يؤدي إلى قيم كبيرة جداً لـ L1-norm لمعاملاتنا $x$، وهو ما قد يُفضي إلى أداء رديء لصيغة MPF. بدلاً من ذلك، يمكن الحصول على حل تقريبي يُصغِّر L1-norm لـ $x$ في محاولة لتحسين سلوك صيغة MPF.

إعداد نظام المعادلات الخطية (LSE)

بعد أن اخترنا قيم $k_j$، يجب أولاً بناء نظام المعادلات الخطية $Ax=b$ كما شُرح أعلاه. تعتمد المصفوفة $A$ ليس فقط على $k_j$ بل أيضاً على اختيارنا لصيغة الضرب، ولا سيما رتبتها. بالإضافة إلى ذلك، يمكنك مراعاة ما إذا كانت صيغة الضرب متماثلة أم لا (انظر [1]) بضبط symmetric=True/False. إلا أن ذلك ليس شرطاً، كما أوضح المرجع [2].

from qiskit_addon_mpf.static import setup_static_lse

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)

لنستعرض القيم المختارة أعلاه لبناء المصفوفة $A$ وشعاع $b$. مع $j=0,1, 2$ وخطوات Trotter $k_j = [1, 2, 4]$، والرتبة $\chi = 2$ وخيار خطوات Trotter غير المتماثلة ($s=1$)، نجد أن عناصر المصفوفة $A$ أسفل الصف الأول تُحدَّد بالتعبير $A_{i>0,j} = k_{j}^{-(2 + 1(i-1))}$، وبالتحديد:

$A_{0,0} = A_{0,1} = A_{0,2} = 1$ $A_{1,j} = k_{j}^{-1} \rightarrow A_{1,0} = \frac{1}{1^2}, \;, A_{1,1} = \frac{1}{2^2}, \;, A_{1,2} = \frac{1}{4^2}$ $A_{2,j} = k_{j}^{-2} \rightarrow A_{2,0} = \frac{1}{1^3}, \;, A_{2,1} = \frac{1}{2^3}, \;, A_{2,2} = \frac{1}{4^3}$

أو في الصورة المصفوفية:

$$A = \begin{bmatrix} 1 & 1 & 1\\ 1 & \frac{1}{2^2} & \frac{1}{4^2} \\ 1 & \frac{1}{2^3} & \frac{1}{4^3} \\ \end{bmatrix}$$

يمكن التحقق من هذا بفحص الكائن lse:

lse.A

array([[1.      , 1.      , 1.      ],
       [1.      , 0.25    , 0.0625  ],
       [1.      , 0.125   , 0.015625]])

أما شعاع القيود $b$ فعناصره على النحو الآتي: $b_{0} = 1$ $b_1 = b_2 = 0$

وبذلك،

$$b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$$

وبالمثل في lse:

```python
lse.b

array([1., 0., 0.])

يحتوي الكائن lse على دوال لإيجاد المعاملات الساكنة $x_j$ التي تحقق نظام المعادلات.

mpf_coeffs = lse.solve()
print(
    f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)

The static coefficients associated with the ansatze are: [ 0.04761905 -0.57142857  1.52380952]

تحسين $x$ باستخدام نموذج دقيق

كبديل لحساب $x=A^{-1}b$، يمكنك أيضاً استخدام setup_exact_model لبناء نسخة من cvxpy.Problem تستخدم نظام المعادلات الخطية كقيود وتكون نتيجتها المثلى هي $x$.

سيتضح في القسم التالي سبب وجود هذه الواجهة.

from qiskit_addon_mpf.costs import setup_exact_problem

model_exact, coeffs_exact = setup_exact_problem(lse)
model_exact.solve()
print(coeffs_exact.value)

[ 0.04761905 -0.57142857  1.52380952]

كمؤشر على ما إذا كان MPF المبني بهذه المعاملات سيعطي نتائج جيدة، يمكننا استخدام معيار L1 (انظر أيضاً المرجع [1]).

print(
    "L1 norm of the exact coefficients:",
    np.linalg.norm(coeffs_exact.value, ord=1),
)  # ord specifies the norm. ord=1 is for L1

L1 norm of the exact coefficients: 2.1428571428556378

تحسين $x$ باستخدام نموذج تقريبي

قد يحدث أن يُعدّ معيار L1 لمجموعة قيم $k_j$ المختارة مرتفعاً جداً. إذا كان الأمر كذلك ولم تتمكن من اختيار مجموعة مختلفة من قيم $k_j$، يمكنك استخدام حل تقريبي لنظام المعادلات الخطية بدلاً من الحل الدقيق.

للقيام بذلك، استخدم ببساطة setup_approximate_model لبناء نسخة مختلفة من cvxpy.Problem، والتي تقيد معيار L1 بحد أقصى محدد مع تصغير الفارق بين $Ax$ و $b$.

from qiskit_addon_mpf.costs import setup_sum_of_squares_problem

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=1.5
)
model_approx.solve()
print(coeffs_approx.value)
print(
    "L1 norm of the approximate coefficients:",
    np.linalg.norm(coeffs_approx.value, ord=1),
)

[-1.10294118e-03 -2.48897059e-01  1.25000000e+00]
L1 norm of the approximate coefficients: 1.5

لاحظ أن لديك حرية كاملة في طريقة حل مسألة التحسين هذه، مما يعني أنك تستطيع تغيير محلّ التحسين وعتبات التقارب الخاصة به وما إلى ذلك. راجع الدليل المختص بـ كيفية استخدام النموذج التقريبي.

معاملات MPF الديناميكية

في القسم السابق، قدمنا MPF ساكناً يتحسن على تقريب Trotter القياسي. غير أن هذا الإصدار الساكن لا يقلل بالضرورة خطأ التقريب. بشكل محدد، إن MPF الساكن، المُعبَّر عنه بـ $\mu^S(t)$، ليس الإسقاط الأمثل لـ $\rho(t)$ على الفضاء الجزئي الممتد بحالات صيغة الجداء $\{\rho_{k_i}(t)\}_{i=1}^r$.

لمعالجة ذلك، نأخذ في الاعتبار MPF ديناميكياً (مُقدَّم في المرجع [2] وموضَّح تجريبياً في المرجع [3]) يُقلل خطأ التقريب بمعيار Frobenius. رسمياً، نركز على تصغير

$$\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr],$$

بالنسبة لبعض المعاملات $x_i(t)$ عند كل زمن $t$. إن المُسقِط الأمثل بمعيار Frobenius هو $\mu^D(t) \;=\; \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t)$، ونسمي $\mu^D(t)$ بـ MPF الديناميكي. وبالتعويض بالتعريفات أعلاه:

$$\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr] \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t) \right) \left( \rho(t) - \sum_{j=1}^r x_j(t)\,\rho_{k_j}(t) \right) \bigr] \;=\; \\ = 1 \;+\; \sum_{i,j=1}^r M_{i,j}(t)\,x_i(t)\,x_j(t) \;-\; 2 \sum_{i=1}^r L_i^{\mathrm{exact}}(t)\,x_i(t),$$

حيث $M_{i,j}(t)$ هي مصفوفة Gram، المعرَّفة بـ

$$M_{i,j}(t) \;=\; \mathrm{Tr}\bigl[\rho_{k_i}(t)\,\rho_{k_j}(t)\bigr] \;=\; \bigl|\langle \psi_{\mathrm{in}} \!\mid S\bigl(t/k_i\bigr)^{-k_i}\,S\bigl(t/k_j\bigr)^{k_j} \!\mid \psi_{\mathrm{in}} \rangle \bigr|^2.$$

و

$$L_i^{\mathrm{exact}}(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$$

تمثل التداخل بين الحالة الدقيقة $\rho(t)$ وكل تقريب لصيغة الجداء $\rho_{k_i}(t)$. في السيناريوهات العملية، قد لا تُقاس هذه التداخلات إلا بصورة تقريبية بسبب الضوضاء أو محدودية الوصول إلى $\rho(t)$.

هنا، $\lvert\psi_{\mathrm{in}}\rangle$ هي الحالة الابتدائية، و $S(\cdot)$ هي العملية المطبقة في صيغة الجداء. باختيار المعاملات $x_i(t)$ التي تقلل هذا التعبير (مع التعامل مع بيانات التداخل التقريبية عندما لا تكون $\rho(t)$ معروفة كلياً)، نحصل على أفضل تقريب ديناميكي لـ $\rho(t)$ داخل الفضاء الجزئي لـ MPF (بمعنى معيار Frobenius). يمكن حساب الكميات $L_i(t)$ و $M_{i,j}(t)$ بكفاءة باستخدام طرق شبكة الموترات [3]. تُوفر إضافة MPF لـ Qiskit عدة "backends" لتنفيذ هذا الحساب. يوضح المثال أدناه الطريقة الأكثر مرونة للقيام بذلك، كما تشرح وثائق TeNPy layer-based backend تفاصيل إضافية. لاستخدام هذه الطريقة، ابدأ من الدائرة التي تُطبق التطور الزمني المطلوب وأنشئ نماذج تمثل هذه العمليات من طبقات الدائرة المقابلة. وأخيراً، يُنشأ كائن Evolver يمكن استخدامه لتوليد الكميات المتطورة زمنياً $M_{i,j}(t)$ و $L_i(t)$. نبدأ بإنشاء كائن Evolver المقابل للتطور الزمني التقريبي (ApproxEvolverFactory) الذي تُنفذه الدوائر. بصفة خاصة، انتبه جيداً لمتغير order للتأكد من تطابقه. لاحظ أنه عند توليد الدوائر المقابلة للتطور الزمني التقريبي، نستخدم قيماً وهمية للوسيط time = 1.0 وعدد خطوات Trotter (reps=1). ثم يُنتج حلّال المسألة الديناميكية في setup_dynamic_lse الدوائر التقريبية الصحيحة.

from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.backends.tenpy_layers import LayerModel
from qiskit_addon_mpf.backends.tenpy_layers import LayerwiseEvolver
from functools import partial

# Create approximate time-evolution circuits
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)  # collect XX and YY

# Find layers in the circuit
layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)

# Create tensor network models
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

# Create the time-evolution object
approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

تحذير

يجب اختيار خيارات LayerwiseEvolver التي تحدد تفاصيل محاكاة شبكة الموترات بعناية لتفادي إعداد مسألة تحسين غير محددة بشكل صحيح.

ثم نُعدّ المطوِّر الدقيق (على سبيل المثال، ExactEvolverFactory)، الذي يُعيد كائن Evolver يحسب التطور الزمني الحقيقي أو "المرجعي". من الناحية العملية، سنقرّب التطور الدقيق باستخدام صيغة Suzuki-Trotter من رتبة أعلى أو طريقة موثوقة أخرى بخطوة زمنية صغيرة. نُقرّب أدناه الحالة المتطورة زمنياً الدقيقة بصيغة Suzuki-Trotter من الرتبة الرابعة باستخدام خطوة زمنية صغيرة dt=0.1، مما يعني أن عدد خطوات Trotter المستخدمة عند الزمن $t$ هو $k=t/dt$. نحدد أيضاً بعض خيارات الاقتطاع الخاصة بـ TeNPy لتحديد أقصى أبعاد الرابطة للشبكة الموترية الأساسية، فضلاً عن أدنى القيم الفردية لرابطات الشبكة الموترية المشقوقة. يمكن أن تؤثر هذه المعاملات على دقة قيمة التوقع المحسوبة بمعاملات MPF الديناميكية، لذا من المهم استكشاف نطاق من القيم لإيجاد التوازن الأمثل بين وقت الحساب والدقة. لاحظ أن حساب معاملات MPF لا يعتمد على قيمة توقع صيغة الجداء المحصّلة من تنفيذ الأجهزة، وبالتالي يمكن ضبطه في مرحلة المعالجة اللاحقة.

single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

بعد ذلك، أنشئ الحالة الابتدائية لنظامك بصيغة متوافقة مع TeNPy (على سبيل المثال، MPS_neel_state=$\vert 0101...01 \rangle$). يُعدّ هذا الدالة الموجية للجسم المتعدد التي ستُطوَّر زمنياً $\lvert\psi_{\mathrm{in}}\rangle$ بوصفها موتراً.

from qiskit_addon_mpf.backends.tenpy_tebd import MPOState
from qiskit_addon_mpf.backends.tenpy_tebd import MPS_neel_state

def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)

mps_initial_state = MPS_neel_state(models[0].lat)

لكل خطوة زمنية $t$ نُعدّ نظام المعادلات الخطية الديناميكي باستخدام طريقة setup_dynamic_lse. يحتوي الكائن المقابل على معلومات حول مسألة MPF الديناميكية: يُعطي lse.A مصفوفة Gram المتمثلة في $M$ بينما يُعطي lse.b التداخل $L$. يمكننا بعد ذلك حل نظام المعادلات الخطية (عندما لا يكون غير محدد) لإيجاد المعاملات الديناميكية باستخدام setup_frobenius_problem. من المهم ملاحظة الفارق مع المعاملات الساكنة، التي تعتمد فقط على تفاصيل صيغة الجداء المستخدمة وهي مستقلة عن تفاصيل التطور الزمني (الهاميلتوني والحالة الابتدائية).

from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import setup_frobenius_problem

mpf_dynamic_coeffs_list = []
for t in trotter_times:
    print(f"Computing dynamic coefficients for time={t}")
    lse = setup_dynamic_lse(
        mpf_trotter_steps,
        t,
        identity_factory,
        exact_factory,
        approx_factory,
        mps_initial_state,
    )
    problem, coeffs = setup_frobenius_problem(lse)
    try:
        problem.solve()
        mpf_dynamic_coeffs_list.append(coeffs.value)
    except Exception as error:
        mpf_dynamic_coeffs_list.append(np.zeros(len(mpf_trotter_steps)))
        print(error, "Calculation Failed for time", t)
    print("")

Computing dynamic coefficients for time=0.5

Computing dynamic coefficients for time=0.6

Computing dynamic coefficients for time=0.7

Computing dynamic coefficients for time=0.7999999999999999

Computing dynamic coefficients for time=0.8999999999999999

Computing dynamic coefficients for time=0.9999999999999999

Computing dynamic coefficients for time=1.0999999999999999

Computing dynamic coefficients for time=1.1999999999999997

Computing dynamic coefficients for time=1.2999999999999998

Computing dynamic coefficients for time=1.4

Computing dynamic coefficients for time=1.4999999999999998

أخيراً، ارسم قيم التوقع هذه عبر زمن التطور.

import matplotlib.pyplot as plt

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    trotter_curve, trotter_curve_error = [], []
    for trotter_expvals, trotter_stds in zip(
        mpf_expvals_all_times, mpf_stds_all_times
    ):
        trotter_curve.append(trotter_expvals[k])
        trotter_curve_error.append(trotter_stds[k])

    plt.errorbar(
        trotter_times,
        trotter_curve,
        yerr=trotter_curve_error,
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )  # , , )

# Get expectation values at all times for the static MPF with exact coeffs
exact_mpf_curve, exact_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_exact.value, trotter_stds)
            ]
        )
    )
    exact_mpf_curve_error.append(mpf_std)
    exact_mpf_curve.append(trotter_expvals @ coeffs_exact.value)

plt.errorbar(
    trotter_times,
    exact_mpf_curve,
    yerr=exact_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Exact",
    color="purple",
)

# Get expectation values at all times for the static MPF with approximate
approx_mpf_curve, approx_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_approx.value, trotter_stds)
            ]
        )
    )
    approx_mpf_curve_error.append(mpf_std)
    approx_mpf_curve.append(trotter_expvals @ coeffs_approx.value)

plt.errorbar(
    trotter_times,
    approx_mpf_curve,
    yerr=approx_mpf_curve_error,
    markersize=4,
    marker="o",
    color="orange",
    label="Static MPF - Approximate",
)

# # Get expectation values at all times for the dynamic MPF
dynamic_mpf_curve, dynamic_mpf_curve_error = [], []
for trotter_expvals, trotter_stds, dynamic_coeffs in zip(
    mpf_expvals_all_times, mpf_stds_all_times, mpf_dynamic_coeffs_list
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(dynamic_coeffs, trotter_stds)
            ]
        )
    )
    dynamic_mpf_curve_error.append(mpf_std)
    dynamic_mpf_curve.append(trotter_expvals @ dynamic_coeffs)

plt.errorbar(
    trotter_times,
    dynamic_mpf_curve,
    yerr=dynamic_mpf_curve_error,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

plt.plot(
    exact_evolution_times,
    exact_expvals,
    lw=3,
    color="red",
    label="Exact time-evolution",
)

plt.title(
    f"Expectation values for (ZZ,{(L//2-1, L//2)}) as a function of time"
)
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend()
plt.grid()

في حالات كالمثال أعلاه، حيث تؤدي صيغة الجداء $k=1$ أداءً سيئاً في جميع الأزمنة، تتأثر جودة نتائج MPF الديناميكي بشكل كبير أيضاً. في مثل هذه الحالات، يُفيد التحقيق في إمكانية استخدام صيغ جداء منفردة بعدد أكبر من خطوات Trotter لتحسين الجودة الإجمالية للنتائج. في هذه المحاكاة، نلاحظ التفاعل بين أنواع مختلفة من الأخطاء: خطأ ناتج عن أخذ عينات محدودة، وخطأ Trotter الناتج عن صيغ الجداء. يُساعد MPF في تقليل خطأ Trotter الناتج عن صيغ الجداء، لكنه يُفضي إلى خطأ أخذ عينات أعلى مقارنةً بصيغ الجداء. قد يكون هذا مفيداً، إذ يمكن لصيغ الجداء تقليل خطأ أخذ العينات بزيادة حجم العينة، غير أن الخطأ المنهجي الناتج عن تقريب Trotter يظل دون تغيير.

سلوك مثير للاهتمام آخر يمكن ملاحظته من الرسم هو أن قيمة التوقع لصيغة الجداء عند $k=1$ تبدأ في التصرف بشكل عشوائي (فضلاً عن كونها تقريباً سيئاً للقيمة الدقيقة) عند الأزمنة التي يكون فيها $t/k > 1$، كما هو موضح في الدليل حول كيفية اختيار عدد خطوات Trotter.

الخطوة 1: تحويل المدخلات الكلاسيكية إلى مسألة كمية

لنأخذ الآن زمناً واحداً $t=1.0$ ونحسب قيمة توقع المغنطة بالطرق المختلفة باستخدام وحدة معالجة كمية واحدة (QPU). تم اختيار $t$ المحدد هذا لتعظيم الفارق بين الطرق المختلفة ومراقبة كفاءتها النسبية. لتحديد النافذة الزمنية التي يُضمن فيها أن ينتج MPF الديناميكي مقادير بخطأ أقل من أي من صيغ Trotter الفردية ضمن الجداء المتعدد، يمكننا تطبيق "اختبار MPF" - انظر المعادلة (17) والنص المحيط بها في المرجع [3].

إعداد دوائر Trotter

عند هذه النقطة، وجدنا معاملات التوسع $x$، ولم يبقَ سوى توليد الدوائر الكمية المُطَوترة. مرة أخرى، يأتي وحدة qiskit_addon_utils.problem_generators للإنقاذ بدالة مفيدة للقيام بذلك:

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

total_time = 1.0
mpf_circuits = []
for k in mpf_trotter_steps:
    # Initial Neel state preparation
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2 != 0])

    trotter_circ = generate_time_evolution_circuit(

hamiltonian, synthesis=SuzukiTrotter(order=order, reps=k), time=total_time, )

circuit.compose(trotter_circ, inplace=True)

mpf_circuits.append(pm.run(circuit))

```python
mpf_circuits[-1].draw("mpl", fold=-1, scale=0.4)

الخطوة 2: تحسين المسألة لتنفيذها على العتاد الكمي

لنعُد إلى حساب قيمة التوقع لنقطة زمنية واحدة. سنختار backend لتنفيذ التجربة على العتاد الفعلي.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
print(backend)

qubits = list(range(backend.num_qubits))

بعد ذلك نُزيل البِتات الكمية الشاذة من خريطة الترابط لضمان ألا تشملها مرحلة التخطيط في المُصرِّف. نستخدم فيما يلي خصائص الـ backend المُسجَّلة المُخزَّنة في الكائن target، ونُزيل البِتات الكمية التي يتجاوز خطأ قياسها أو خطأ بوابتها ثنائية البِت حدًا معيَّنًا (max_meas_err، max_twoq_err)، أو يقل زمن $T_2$ لديها (الذي يحدد فقدان التماسك) عن حد أدنى معيَّن (min_t2).

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

max_meas_err = 0.012
min_t2 = 40
max_twoq_err = 0.005

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    if target.qubit_properties[q].t2 is not None:
        t2 = target.qubit_properties[q].t2 * 1e6
    else:
        t2 = 0
    if meas_err > max_meas_err or t2 < min_t2:
        # print(q)
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
    twoq_gate_err = target[instruction_2q][q].error
    if twoq_gate_err > max_twoq_err:
        # print(q)
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates
    + ["measure"],  # or whatever new set of gates
    coupling_map=cust_cmap,
)

sorted_components = sorted(
    [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
    reverse=True,
)
print("size of largest component", len(sorted_components[0]))

size of largest component 10

نريد ضبط max_meas_err وmin_t2 وmax_twoq_err بحيث نجد مجموعة فرعية كبيرة بما يكفي من البِتات الكمية تدعم تشغيل الدائرة. يكفينا في حالتنا إيجاد سلسلة أحادية البعد مؤلفة من 10 بِتات كمية.

cust_cmap.draw()

يمكننا حينئذٍ تعيين الدائرة والمُراقِبة على البِتات الكمية الفيزيائية للجهاز.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
    optimization_level=3, target=cust_target
)

transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

qubits_layouts = [
    [
        idx
        for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
        if qb._register.name != "ancilla"
    ]
    for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
    transpiler = generate_preset_pass_manager(
        optimization_level=3, backend=backend, initial_layout=layout
    )
    transpiled_circuit = transpiler.run(circuit)
    transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
 isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

print(transpiled_circuits[-1].depth(lambda x: x.operation.num_qubits == 2))
print(transpiled_circuits[-1].count_ops())
transpiled_circuits[-1].draw("mpl", idle_wires=False, fold=False)

51
OrderedDict([('sx', 310), ('rz', 232), ('cz', 132), ('x', 19)])

الخطوة 3: التنفيذ باستخدام Qiskit primitives

باستخدام الأداة الأولية Estimator يمكننا الحصول على تقدير قيمة التوقع من وحدة المعالجة الكمية QPU. ننفِّذ دوائر AQC المُحسَّنة مع تقنيات إضافية للتخفيف من الأخطاء وقمعها.

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
estimator.options.dynamical_decoupling.enable = True
estimator.options.twirling.enable_gates = True
estimator.options.twirling.enable_measure = True
estimator.options.twirling.num_randomizations = "auto"
estimator.options.twirling.strategy = "active-accum"
estimator.options.resilience.measure_mitigation = True
estimator.options.experimental.execution_path = "gen3-turbo"

estimator.options.resilience.zne_mitigation = True
estimator.options.resilience.zne.noise_factors = (1, 3, 5)
estimator.options.resilience.zne.extrapolator = ("exponential", "linear")

estimator.options.environment.job_tags = ["mpf small"]

job = estimator.run(
    [
        (circ, observable)
        for circ, observable in zip(transpiled_circuits, isa_observables)
    ]
)

الخطوة 4: المعالجة اللاحقة وإعادة النتيجة بالصيغة الكلاسيكية المطلوبة

خطوة المعالجة اللاحقة الوحيدة هي دمج قيمة التوقع المُستخرَجة من الأدوات الأولية لـ Qiskit Runtime عند خطوات Trotter المختلفة باستخدام معاملات MPF المقابلة. للمُراقِبة $A$ يكون لدينا:

$\langle A \rangle_{\text{mpf}} = \text{Tr} [A \mu(t)] = \sum_{j} x_j \text{Tr} [A \rho_{k_j}] = \sum_{j} x_j \langle A \rangle_j$

أولاً، نستخرج قيم التوقع الفردية المُتحصَّل عليها لكل دائرة من دوائر Trotter:

result_exp = job.result()
evs_exp = [res.data.evs for res in result_exp]
evs_std = [res.data.stds for res in result_exp]

print(evs_exp)

[array(-0.06361607), array(-0.23820448), array(-0.50271805)]

بعد ذلك، نُعيد دمجها ببساطة مع معاملات MPF لدينا للحصول على إجمالي قيم التوقع لـ MPF. فيما يلي، نفعل ذلك لكل طريقة من الطرق المختلفة التي حسبنا بها $x$.

exact_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_exact.value, evs_std)
        ]
    )
)
print(
    "Exact static MPF expectation value: ",
    evs_exp @ coeffs_exact.value,
    "+-",
    exact_mpf_std,
)
approx_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_approx.value, evs_std)
        ]
    )
)
print(
    "Approximate static MPF expectation value: ",
    evs_exp @ coeffs_approx.value,
    "+-",
    approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(mpf_dynamic_coeffs_list[7], evs_std)
        ]
    )
)
print(
    "Dynamic MPF expectation value: ",
    evs_exp @ mpf_dynamic_coeffs_list[7],
    "+-",
    dynamic_mpf_std,
)

Exact static MPF expectation value:  -0.6329590442738475 +- 0.012798249760406036
Approximate static MPF expectation value:  -0.5690390035339492 +- 0.010459559917168473
Dynamic MPF expectation value:  -0.4655579758795695 +- 0.007639139186720507

أخيرًا، يمكننا لهذه المسألة الصغيرة حساب القيمة المرجعية الدقيقة باستخدام scipy.linalg.expm على النحو التالي:

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

exp_H = expm(-1j * total_time * hamiltonian.to_matrix())

initial_state_circuit = QuantumCircuit(L)
initial_state_circuit.x([i for i in range(L) if i % 2 != 0])
initial_state = Statevector(initial_state_circuit).data

time_evolved_state = exp_H @ initial_state

exact_obs = (
    time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
)
print("Exact expectation value ", exact_obs.real)

Exact expectation value  -0.39909900734489434

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs_exp[k],
        yerr=evs_std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )  # , , )

plt.errorbar(
    3,
    evs_exp @ coeffs_exact.value,
    yerr=exact_mpf_std,
    markersize=4,
    marker="o",
    color="purple",
    label="Static MPF",
)

plt.errorbar(
    4,
    evs_exp @ coeffs_approx.value,
    yerr=approx_mpf_std,
    markersize=4,
    marker="o",
    color="orange",
    label="Approximate static MPF",
)

plt.errorbar(
    5,
    evs_exp @ mpf_dynamic_coeffs_list[7],
    yerr=dynamic_mpf_std,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

plt.axhline(
    y=exact_obs.real,
    linestyle="--",
    color="red",
    label="Exact time-evolution",
)

plt.title(
    f"Expectation values for (ZZ,{(L//2-1, L//2)}) at time {total_time} for the different methods "
)
plt.xlabel("Method")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

في المثال أعلاه، تُحقِّق طريقة Dynamic MPF أفضل أداء من حيث قيمة التوقع، متفوقةً على ما كنا سنحصل عليه باستخدام أعلى عدد من خطوات Trotter وحده. وعلى الرغم من أن تقنيات MPF المختلفة لا تُحقِّق دائمًا قيمةَ توقع أفضل مقارنةً بأعلى عدد من خطوات Trotter (كالنموذجين الدقيق والتقريبي في الرسم البياني أعلاه)، فإن الانحراف المعياري لهذه القيم يعكس بدقة التباين المتزايد الناجم عن استخدام تقنية MPF. يُبرز هذا عدم اليقين المحيط بقيمة التوقع المُتحصَّل عليها، التي تشتمل دائمًا على قيمة التوقع المتوقعة من التطور الزمني الدقيق للنظام. في المقابل، تفشل قيم التوقع المحسوبة بعدد أقل من خطوات Trotter في احتواء قيمة التوقع الدقيقة ضمن نطاق عدم يقينها، مما يعني أنها تُعيد نتيجة خاطئة بثقة.

def relative_error(ev, exact_ev):
    return abs(ev - exact_ev)

relative_error_k = [relative_error(ev, exact_obs.real) for ev in evs_exp]
relative_error_mpf = relative_error(evs_exp @ mpf_coeffs, exact_obs.real)
relative_error_approx_mpf = relative_error(
    evs_exp @ coeffs_approx.value, exact_obs.real
)
relative_error_dynamic_mpf = relative_error(
    evs_exp @ mpf_dynamic_coeffs_list[7], exact_obs.real
)

print("relative error for each trotter steps", relative_error_k)
print("relative error with MPF exact coeffs", relative_error_mpf)
print("relative error with MPF approx coeffs", relative_error_approx_mpf)
print("relative error with MPF dynamic coeffs", relative_error_dynamic_mpf)

relative error for each trotter steps [0.33548293650112293, 0.16089452939226306, 0.10361904247828346]
relative error with MPF exact coeffs 0.2338600369291003
relative error with MPF approx coeffs 0.16993999618905486
relative error with MPF dynamic coeffs 0.06645896853467514

الجزء الثاني: توسيع النطاق

لنُوسِّع نطاق المسألة إلى ما هو أبعد مما يمكن محاكاته بدقة. في هذا القسم سنركز على إعادة إنتاج بعض النتائج الواردة في المرجع [3].

الخطوة 1: تعيين المدخلات الكلاسيكية إلى مسألة كمية

Hamiltonian

في المثال واسع النطاق، نستخدم نموذج XXZ على سلسلة مكوَّنة من 50 موقعًا:

$$\hat{\mathcal{H}}_{XXZ} = \sum_{i=1}^{L-1} J_{i,(i+1)}\left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ 2\cdot Z_i Z_{(i+1)} \right) \, ,$$

حيث $J_{i,(i+1)}$ معامل عشوائي يقابل الحافة $(i, i+1)$. هذا هو الـ Hamiltonian المُستخدَم في العرض التوضيحي المقدَّم في المرجع [3].

L = 50
# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

import numpy as np
from qiskit.quantum_info import SparsePauliOp, Pauli

# Generate random coefficients for our XXZ Hamiltonian
np.random.seed(0)

even_edges = list(coupling_map.get_edges())[::2]
odd_edges = list(coupling_map.get_edges())[1::2]

Js = np.random.uniform(0.5, 1.5, size=L)
hamiltonian = SparsePauliOp(Pauli("I" * L))
for i, edge in enumerate(even_edges + odd_edges):
    hamiltonian += SparsePauliOp.from_sparse_list(
        [
            ("XX", (edge), 2 * Js[i]),
            ("YY", (edge), 2 * Js[i]),
            ("ZZ", (edge), 4 * Js[i]),
        ],
        num_qubits=L,
    )

print(hamiltonian)

SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'XXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'YYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'ZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII'],
              coeffs=[1.        +0.j, 2.09762701+0.j, 2.09762701+0.j, 4.19525402+0.j,
 2.43037873+0.j, 2.43037873+0.j, 4.86075747+0.j, 2.20552675+0.j,
 2.20552675+0.j, 4.4110535 +0.j, 2.08976637+0.j, 2.08976637+0.j,
 4.17953273+0.j, 1.8473096 +0.j, 1.8473096 +0.j, 3.6946192 +0.j,
 2.29178823+0.j, 2.29178823+0.j, 4.58357645+0.j, 1.87517442+0.j,
 1.87517442+0.j, 3.75034885+0.j, 2.783546  +0.j, 2.783546  +0.j,
 5.567092  +0.j, 2.92732552+0.j, 2.92732552+0.j, 5.85465104+0.j,
 1.76688304+0.j, 1.76688304+0.j, 3.53376608+0.j, 2.58345008+0.j,
 2.58345008+0.j, 5.16690015+0.j, 2.05778984+0.j, 2.05778984+0.j,
 4.11557968+0.j, 2.13608912+0.j, 2.13608912+0.j, 4.27217824+0.j,
 2.85119328+0.j, 2.85119328+0.j, 5.70238655+0.j, 1.14207212+0.j,
 1.14207212+0.j, 2.28414423+0.j, 1.1742586 +0.j, 1.1742586 +0.j,
 2.3485172 +0.j, 1.04043679+0.j, 1.04043679+0.j, 2.08087359+0.j,
 2.66523969+0.j, 2.66523969+0.j, 5.33047938+0.j, 2.5563135 +0.j,
 2.5563135 +0.j, 5.112627  +0.j, 2.7400243 +0.j, 2.7400243 +0.j,
 5.48004859+0.j, 2.95723668+0.j, 2.95723668+0.j, 5.91447337+0.j,
 2.59831713+0.j, 2.59831713+0.j, 5.19663426+0.j, 1.92295872+0.j,
 1.92295872+0.j, 3.84591745+0.j, 2.56105835+0.j, 2.56105835+0.j,
 5.12211671+0.j, 1.23654885+0.j, 1.23654885+0.j, 2.4730977 +0.j,
 2.27984204+0.j, 2.27984204+0.j, 4.55968409+0.j, 1.28670657+0.j,
 1.28670657+0.j, 2.57341315+0.j, 2.88933783+0.j, 2.88933783+0.j,
 5.77867567+0.j, 2.04369664+0.j, 2.04369664+0.j, 4.08739329+0.j,
 1.82932388+0.j, 1.82932388+0.j, 3.65864776+0.j, 1.52911122+0.j,
 1.52911122+0.j, 3.05822245+0.j, 2.54846738+0.j, 2.54846738+0.j,
 5.09693476+0.j, 1.91230066+0.j, 1.91230066+0.j, 3.82460133+0.j,
 2.1368679 +0.j, 2.1368679 +0.j, 4.2737358 +0.j, 1.0375796 +0.j,
 1.0375796 +0.j, 2.0751592 +0.j, 2.23527099+0.j, 2.23527099+0.j,
 4.47054199+0.j, 2.22419145+0.j, 2.22419145+0.j, 4.44838289+0.j,
 2.23386799+0.j, 2.23386799+0.j, 4.46773599+0.j, 2.88749616+0.j,
 2.88749616+0.j, 5.77499231+0.j, 2.3636406 +0.j, 2.3636406 +0.j,
 4.7272812 +0.j, 1.7190158 +0.j, 1.7190158 +0.j, 3.4380316 +0.j,
 1.87406391+0.j, 1.87406391+0.j, 3.74812782+0.j, 2.39526239+0.j,
 2.39526239+0.j, 4.79052478+0.j, 1.12045094+0.j, 1.12045094+0.j,
 2.24090189+0.j, 2.33353343+0.j, 2.33353343+0.j, 4.66706686+0.j,
 2.34127574+0.j, 2.34127574+0.j, 4.68255148+0.j, 1.42076512+0.j,
 1.42076512+0.j, 2.84153024+0.j, 1.2578526 +0.j, 1.2578526 +0.j,
 2.51570519+0.j, 1.6308567 +0.j, 1.6308567 +0.j, 3.2617134 +0.j])

بالنسبة للمؤثر (observable) نختار $Z_{24}Z_{25}$، كما هو موضح في اللوحة السفلية من الشكل 5 في المرجع [3].

observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)

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

اختيار خطوات Trotter

تستخدم التجربة الموضحة في الشكل 4 من المرجع [3] خطوات Trotter المتماثلة $k_j = [2, 3, 4]$ من الرتبة $2$. نركّز على النتائج عند الزمن $t=3$، حيث يكون لكل من MPF وصيغة Trotter بعدد أكبر من الخطوات (6 في هذه الحالة) خطأ Trotter متساوٍ. غير أن قيمة توقع MPF تُحسب من الدوائر المقابلة لعدد أقل من خطوات Trotter، وبالتالي تكون هذه الدوائر أكثر ضحالة (shallower). من الناحية العملية، حتى لو كان لكل من MPF ودائرة خطوات Trotter الأعمق نفس خطأ Trotter، فإننا نتوقع أن تكون قيمة التوقع التجريبية المحسوبة من دوائر MPF أقرب إلى القيمة النظرية، نظرًا لأنها تنطوي على تشغيل دوائر ضحلة أقل تعرضًا لضوضاء الأجهزة مقارنةً بالدائرة المقابلة لصيغة Trotter ذات عدد الخطوات الأعلى.

total_time = 3
mpf_trotter_steps = [2, 3, 4]
order = 2
symmetric = True

إعداد نظام المعادلات الخطية (LSE)

هنا نستعرض معاملات MPF الثابتة (static MPF coefficients) لهذه المسألة.

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(
    f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)
print("L1 norm:", np.linalg.norm(mpf_coeffs, ord=1))

The static coefficients associated with the ansatze are: [ 0.26666667 -2.31428571  3.04761905]
L1 norm: 5.628571428571431

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=2.0
)
model_approx.solve()
print(coeffs_approx.value)
print(
    "L1 norm of the approximate coefficients:",
    np.linalg.norm(coeffs_approx.value, ord=1),
)

[-0.24255546 -0.25744454  1.5       ]
L1 norm of the approximate coefficients: 2.0

المعاملات الديناميكية (Dynamic coefficients)

# Create approximate time-evolution circuits
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)  # collect XX and YY

# Find layers in the circuit
layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)

# Create tensor network models
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

# Create the time-evolution object
approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 4,
    },
)

# Create exact time-evolution circuits
single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

# Create the time-evolution object
exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 3,
    },
)

def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)

mps_initial_state = MPS_neel_state(models[0].lat)

lse = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs = setup_frobenius_problem(lse)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs.value
except Exception as error:
    print(error, "Calculation Failed for time", total_time)
print("")

إنشاء كل دائرة Trotter في تحليل MPF الخاص بنا

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

mpf_circuits = []
for k in mpf_trotter_steps:
    # Initial state preparation |1010..>
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2])

    trotter_circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(reps=k, order=order),
        time=total_time,
    )

    circuit.compose(trotter_circ, qubits=range(L), inplace=True)

    mpf_circuits.append(circuit)

إنشاء دائرة Trotter ذات خطأ Trotter مماثل لـ MPF

k = 6

# Initial state preparation |1010..>
comp_circuit = QuantumCircuit(L)
comp_circuit.x([i for i in range(L) if i % 2])

trotter_circ = generate_time_evolution_circuit(
    hamiltonian,
    synthesis=SuzukiTrotter(reps=k, order=order),
    time=total_time,
)

comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)

mpf_circuits.append(comp_circuit)

الخطوة 2: تحسين المسألة لتنفيذها على الأجهزة الكمومية

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

cmap = target.build_coupling_map(filter_idle_qubits=True)
cmap_list = list(cmap.get_edges())

max_meas_err = 0.055
min_t2 = 30
max_twoq_err = 0.01

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    if target.qubit_properties[q].t2 is not None:
        t2 = target.qubit_properties[q].t2 * 1e6
    else:
        t2 = 0
    if meas_err > max_meas_err or t2 < min_t2:
        # print(q)
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
    twoq_gate_err = target[instruction_2q][q].error
    if twoq_gate_err > max_twoq_err:
        # print(q)
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates
    + ["measure"],  # or whatever new set of gates
    coupling_map=cust_cmap,
)

sorted_components = sorted(
    [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
    reverse=True,
)
print("size of largest component", len(sorted_components[0]))

size of largest component 73

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
    optimization_level=3, target=cust_target
)

transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

qubits_layouts = [
    [
        idx
        for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
        if qb._register.name != "ancilla"
    ]
    for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
    transpiler = generate_preset_pass_manager(
        optimization_level=3, backend=backend, initial_layout=layout
    )
    transpiled_circuit = transpiler.run(circuit)
    transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

الخطوة 3: التنفيذ باستخدام Qiskit primitives

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
estimator.options.dynamical_decoupling.enable = True
estimator.options.twirling.enable_gates = True
estimator.options.twirling.enable_measure = True
estimator.options.twirling.num_randomizations = "auto"
estimator.options.twirling.strategy = "active-accum"
estimator.options.resilience.measure_mitigation = True
estimator.options.experimental.execution_path = "gen3-turbo"

estimator.options.resilience.zne_mitigation = True
estimator.options.resilience.zne.noise_factors = (1, 1.2, 1.4)
estimator.options.resilience.zne.extrapolator = "linear"

estimator.options.environment.job_tags = ["mpf large"]

job_50 = estimator.run(
    [
        (circ, observable)
        for circ, observable in zip(transpiled_circuits, isa_observables)
    ]
)

الخطوة 4: المعالجة اللاحقة وإرجاع النتيجة بالتنسيق الكلاسيكي المطلوب

result = job_50.result()
evs = [res.data.evs for res in result]
std = [res.data.stds for res in result]

print(evs)
print(std)

[array(-0.08034071), array(-0.00605026), array(-0.15345759), array(-0.18127293)]
[array(0.04482517), array(0.03438413), array(0.21540776), array(0.21520829)]

exact_mpf_std = np.sqrt(
    sum([(coeff**2) * (std**2) for coeff, std in zip(mpf_coeffs, std[:3])])
)
print(
    "Exact static MPF expectation value: ",
    evs[:3] @ mpf_coeffs,
    "+-",
    exact_mpf_std,
)
approx_mpf_std = np.sqrt(
    sum(
        [

            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_approx.value, std[:3])
        ]
    )
)
print(
    "Approximate static MPF expectation value: ",
    evs[:3] @ coeffs_approx.value,
    "+-",
    approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(mpf_dynamic_coeffs, std[:3])
        ]
    )
)
print(
    "Dynamic MPF expectation value: ",
    evs[:3] @ mpf_dynamic_coeffs,
    "+-",
    dynamic_mpf_std,
)

Exact static MPF expectation value:  -0.47510243192011536 +- 0.6613940032465087
Approximate static MPF expectation value:  -0.20914170384216998 +- 0.32341567460419135
Dynamic MPF expectation value:  -0.07994951978722761 +- 0.07423091963310202

sym = {2: "^", 3: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs[k],
        yerr=std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

plt.errorbar(
    3,
    evs[-1],
    yerr=std[-1],
    alpha=0.5,
    markersize=8,
    marker="x",
    color="blue",
    label="6 Trotter steps",
)

plt.errorbar(
    4,
    evs[:3] @ mpf_coeffs,
    yerr=exact_mpf_std,
    markersize=4,
    marker="o",
    color="purple",
    label="Static MPF",
)

plt.errorbar(
    5,
    evs[:3] @ coeffs_approx.value,
    yerr=approx_mpf_std,
    markersize=4,
    marker="o",
    color="orange",
    label="Approximate static MPF",
)

plt.errorbar(
    6,
    evs[:3] @ mpf_dynamic_coeffs,
    yerr=dynamic_mpf_std,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

exact_obs = -0.24384471447172074  # Calculated via Tensor Network calculation
plt.axhline(
    y=exact_obs, linestyle="--", color="red", label="Exact time-evolution"
)

plt.title(
    f"Expectation values for (ZZ,{(L//2-1, L//2)}) at time {total_time} for the different methods "
)
plt.xlabel("Method")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

عند تنفيذ الدوائر الكمومية على الأجهزة الفعلية، قد نواجه تحديات إضافية في الحصول على قيم توقع دقيقة بسبب وجود الضوضاء في الأجهزة. لا يأخذ فورمالية MPF هذا العامل بعين الاعتبار، وقد يؤدي ذلك إلى نتائج تعمل ضد حل MPF. على سبيل المثال، قد يكون هذا هو السبب وراء فشل المعاملات الديناميكية في تقديم تقدير أفضل لقيمة التوقع مقارنةً بالمعامل الساكن التقريبي في الرسم البياني. أي أن المطوّر التقريبي، الذي يحاكي الدائرة التقريبية، لا يعكس بدقة النتائج التي يتم الحصول عليها عند تنفيذ الدوائر التقريبية في وجود ضوضاء الأجهزة. لهذه الأسباب، يُوصى بدمج تقنيات مختلفة للتخفيف من الأخطاء للحصول على نتائج أقرب ما يمكن إلى القيم المثالية لكل من صيغ الضرب. وسيُظهر ذلك فوائد متسقة من نهج MPF.

بشكل عام، لا تزال المعاملات الساكنة التقريبية تُعطي حلًا أكثر دقة من صيغة الضرب ذات العدد الأعلى من خطوات Trotter مع نفس مقدار خطأ Trotter في الإعداد الخالي من الضوضاء.

من المهم أيضًا الإشارة إلى أن المثال الذي يُعيد إنتاج التجربة الواردة في المرجع [3]، تقع النقطة الزمنية $t=3$ خارج الحد الذي يُتوقع عنده أن يؤدي PF مع $k=2$ بشكل جيد، وهو $t/k>1$ كما هو مناقش في هذا الدليل.

المراجع

[1] Vazquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. (2023). Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation. Quantum, 7, 1067.

[2] Zhuk, S., Robertson, N. F., & Bravyi, S. (2024). Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation. Physical Review Research, 6(3), 033309.