تقطير كريلوف الكمي لهاميلتونيات الشبكة

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

# This cell is hidden from users – it disables some lint rules
# ruff: noqa: E402 E722 F601

الخلفية النظرية

يوضّح هذا البرنامج التعليمي كيفية تنفيذ خوارزمية تقطير كريلوف الكمي (KQD) في سياق أنماط Qiskit. ستتعرّف أولاً على النظرية التي تقوم عليها الخوارزمية، ثم ترى عرضاً توضيحياً لتنفيذها على وحدة معالجة كمية (QPU).

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

المتطلبات

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

• Qiskit SDK الإصدار 2.0 أو أحدث، مع دعم التصوير البياني
• Qiskit Runtime الإصدار 0.22 أو أحدث ( pip install qiskit-ibm-runtime )

الإعداد

import numpy as np
import scipy as sp
import matplotlib.pylab as plt
from typing import Union, List
import itertools as it
import copy
from sympy import Matrix
import warnings

warnings.filterwarnings("ignore")

from qiskit.quantum_info import SparsePauliOp, Pauli, StabilizerState
from qiskit.circuit import Parameter, IfElseOp
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter
from qiskit.transpiler import Target, CouplingMap
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    EstimatorV2 as Estimator,
)

def solve_regularized_gen_eig(
    h: np.ndarray,
    s: np.ndarray,
    threshold: float,
    k: int = 1,
    return_dimn: bool = False,
) -> Union[float, List[float]]:
    """
    Method for solving the generalized eigenvalue problem with regularization

    Args:
        h (numpy.ndarray):
            The effective representation of the matrix in the Krylov subspace
        s (numpy.ndarray):
            The matrix of overlaps between vectors of the Krylov subspace
        threshold (float):
            Cut-off value for the eigenvalue of s
        k (int):
            Number of eigenvalues to return
        return_dimn (bool):
            Whether to return the size of the regularized subspace

    Returns:
        lowest k-eigenvalue(s) that are the solution of the regularized generalized eigenvalue problem

    """
    s_vals, s_vecs = sp.linalg.eigh(s)
    s_vecs = s_vecs.T
    good_vecs = np.array(
        [vec for val, vec in zip(s_vals, s_vecs) if val > threshold]
    )
    h_reg = good_vecs.conj() @ h @ good_vecs.T
    s_reg = good_vecs.conj() @ s @ good_vecs.T
    if k == 1:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][0], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][0]
    else:
        if return_dimn:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k], len(good_vecs)
        else:
            return sp.linalg.eigh(h_reg, s_reg)[0][:k]

def single_particle_gs(H_op, n_qubits):
    """
    Find the ground state of the single particle(excitation) sector
    """
    H_x = []
    for p, coeff in H_op.to_list():
        H_x.append(set([i for i, v in enumerate(Pauli(p).x) if v]))

    H_z = []
    for p, coeff in H_op.to_list():
        H_z.append(set([i for i, v in enumerate(Pauli(p).z) if v]))

    H_c = H_op.coeffs

    print("n_sys_qubits", n_qubits)

    n_exc = 1
    sub_dimn = int(sp.special.comb(n_qubits + 1, n_exc))
    print("n_exc", n_exc, ", subspace dimension", sub_dimn)

    few_particle_H = np.zeros((sub_dimn, sub_dimn), dtype=complex)

    sparse_vecs = [
        set(vec) for vec in it.combinations(range(n_qubits + 1), r=n_exc)
    ]  # list all of the possible sets of n_exc indices of 1s in n_exc-particle states

    m = 0
    for i, i_set in enumerate(sparse_vecs):
        for j, j_set in enumerate(sparse_vecs):
            m += 1

            if len(i_set.symmetric_difference(j_set)) <= 2:
                for p_x, p_z, coeff in zip(H_x, H_z, H_c):
                    if i_set.symmetric_difference(j_set) == p_x:
                        sgn = ((-1j) ** len(p_x.intersection(p_z))) * (
                            (-1) ** len(i_set.intersection(p_z))
                        )
                    else:
                        sgn = 0

                    few_particle_H[i, j] += sgn * coeff

    gs_en = min(np.linalg.eigvalsh(few_particle_H))
    print("single particle ground state energy: ", gs_en)
    return gs_en

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

فضاء كريلوف

فضاء كريلوف $\mathcal{K}^r$ من الرتبة $r$ هو الفضاء الممتد بالمتجهات الناتجة عن ضرب القوى المتصاعدة لمصفوفة $A$، حتى $r-1$، في متجه مرجعي $\vert v \rangle$.

$$\mathcal{K}^r = \left\{ \vert v \rangle, A \vert v \rangle, A^2 \vert v \rangle, ..., A^{r-1} \vert v \rangle \right\}$$

إذا كانت المصفوفة $A$ هي هاميلتوني $H$، فسنشير إلى الفضاء المقابل باسم فضاء كريلوف القوّي $\mathcal{K}_P$. أما في حالة كون $A$ هي عامل التطور الزمني المولَّد بواسطة الهاميلتوني $U=e^{-iHt}$، فسنشير إلى الفضاء باسم فضاء كريلوف الأحادي $\mathcal{K}_U$. لا يمكن توليد الفضاء الجزئي لكريلوف القوّي الذي نستخدمه كلاسيكياً بصورة مباشرة على حاسوب كمي، إذ إن $H$ ليست عاملاً أحادياً. غير أنه يمكننا استخدام عامل التطور الزمني $U = e^{-iHt}$ الذي يمكن إثبات أنه يُعطي ضمانات تقارب مماثلة للطريقة القوّية. وتصبح قوى $U$ حينئذٍ خطوات زمنية مختلفة $U^k = e^{-iH(kt)}$.

$$\mathcal{K}_U^r = \left\{ \vert \psi \rangle, U \vert \psi \rangle, U^2 \vert \psi \rangle, ..., U^{r-1} \vert \psi \rangle \right\}$$

راجع الملحق للاطلاع على اشتقاق تفصيلي لكيفية تمكين فضاء كريلوف الأحادي من تمثيل الحالات الذاتية منخفضة الطاقة بدقة.

خوارزمية تقطير كريلوف الكمي

بالنظر إلى هاميلتوني $H$ المراد قطره، نبدأ بالنظر في فضاء كريلوف الأحادي المقابل $\mathcal{K}_U$. الهدف هو إيجاد تمثيل مضغوط للهاميلتوني في $\mathcal{K}_U$، وهو ما سنشير إليه بـ $\tilde{H}$. يمكن حساب عناصر المصفوفة لـ $\tilde{H}$، أي إسقاط الهاميلتوني في فضاء كريلوف، عن طريق حساب القيم المتوقعة التالية

$$\tilde{H}_{mn} = \langle \psi_m \vert H \vert \psi_n \rangle =$$ $$= \langle \psi \vert e^{i H t_m} H e^{-i H t_n} \vert \psi \rangle$$ $$= \langle \psi \vert e^{i H m dt} H e^{-i H n dt} \vert \psi \rangle$$

حيث $\vert \psi_n \rangle = e^{-i H t_n} \vert \psi \rangle$ هي متجهات فضاء كريلوف الأحادي و $t_n = n dt$ هي مضاعفات خطوة الزمن $dt$ المختارة. على الحاسوب الكمي، يمكن حساب كل عنصر من عناصر المصفوفة بأي خوارزمية تتيح الحصول على التداخل بين الحالات الكمية. يركز هذا البرنامج التعليمي على اختبار Hadamard. وبالنظر إلى أن $\mathcal{K}_U$ ذات بُعد $r$، سيكون للهاميلتوني المُسقَط في الفضاء الجزئي أبعاد $r \times r$. مع صغر $r$ بما فيه الكفاية (يكفي عموماً $r<<100$ للحصول على تقارب في تقديرات طاقات القيم الذاتية)، يمكننا حينئذٍ قطر الهاميلتوني المُسقَط $\tilde{H}$ بسهولة. غير أننا لا نستطيع قطر $\tilde{H}$ مباشرةً بسبب عدم تعامد متجهات فضاء كريلوف. سيتوجب علينا قياس تداخلاتها وبناء مصفوفة $\tilde{S}$

$$\tilde{S}_{mn} = \langle \psi_m \vert \psi_n \rangle$$

يتيح لنا ذلك حل مسألة القيمة الذاتية في فضاء غير متعامد (تُعرف أيضاً بمسألة القيمة الذاتية المعممة)

$$\tilde{H} \ \vec{c} = E \ \tilde{S} \ \vec{c}$$

يمكن حينئذٍ الحصول على تقديرات للقيم الذاتية والحالات الذاتية لـ $H$ بالنظر في مقابلاتها في $\tilde{H}$. على سبيل المثال، يُحصَل على تقدير طاقة الحالة الأساسية بأخذ أصغر قيمة ذاتية $c$ والحالة الأساسية من المتجه الذاتي المقابل $\vec{c}$. تحدد المعاملات في $\vec{c}$ مساهمة المتجهات المختلفة التي تمتد إلى $\mathcal{K}_U$.

يُظهر الشكل تمثيلاً دائرياً لاختبار Hadamard المعدَّل، وهو أسلوب يُستخدم لحساب التداخل بين الحالات الكمية المختلفة. لكل عنصر من عناصر المصفوفة $\tilde{H}_{i,j}$، يُجرى اختبار Hadamard بين الحالتين $\vert \psi_i \rangle$ و $\vert \psi_j \rangle$. يُبرز ذلك مخطط الألوان لعناصر المصفوفة والعمليتين $\text{Prep} \; \psi_i$ و $\text{Prep} \; \psi_j$ المقابلتين. وبذلك، يلزم إجراء مجموعة من اختبارات Hadamard لجميع التوليفات الممكنة لمتجهات فضاء كريلوف لحساب جميع عناصر مصفوفة الهاميلتوني المُسقَطة $\tilde{H}$. السلك العلوي في دائرة اختبار Hadamard هو كيوبت مساعد (ancilla) يُقاس إما في الأساس X أو الأساس Y، وتحدد قيمته المتوقعة قيمة التداخل بين الحالتين. يمثل السلك السفلي جميع كيوبتات هاميلتوني النظام. تُعدّ عملية $\text{Prep} \; \psi_i$ كيوبت النظام في الحالة $\vert \psi_i \rangle$ بشكل مشروط بحالة الكيوبت المساعد (وبالمثل لعملية $\text{Prep} \; \psi_j$)، فيما تمثل العملية $P$ تحليل باولي لهاميلتوني النظام $H = \sum_i P_i$. يُعطى اشتقاق أكثر تفصيلاً للعمليات التي يحسبها اختبار Hadamard أدناه.

تعريف الهاميلتوني

لنأخذ هاميلتوني هايزنبرغ لـ $N$ كيوبتاً على سلسلة خطية: $H= \sum_{i,j}^N X_i X_j + Y_i Y_j - J Z_i Z_j$

# Define problem Hamiltonian.
n_qubits = 30
J = 1  # coupling strength for ZZ interaction

# Define the Hamiltonian:
H_int = [["I"] * n_qubits for _ in range(3 * (n_qubits - 1))]
for i in range(n_qubits - 1):
    H_int[i][i] = "Z"
    H_int[i][i + 1] = "Z"
for i in range(n_qubits - 1):
    H_int[n_qubits - 1 + i][i] = "X"
    H_int[n_qubits - 1 + i][i + 1] = "X"
for i in range(n_qubits - 1):
    H_int[2 * (n_qubits - 1) + i][i] = "Y"
    H_int[2 * (n_qubits - 1) + i][i + 1] = "Y"
H_int = ["".join(term) for term in H_int]
H_tot = [(term, J) if term.count("Z") == 2 else (term, 1) for term in H_int]

# Get operator
H_op = SparsePauliOp.from_list(H_tot)
print(H_tot)

[('ZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 1), ('XXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 1), ('YYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 1), ('IIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 1)]

ضبط معاملات الخوارزمية

نختار قيمة تجريبية لخطوة الزمن dt (استناداً إلى حدود علوية لمعيار الهاميلتوني). أثبت المرجع [2] أن خطوة الزمن الصغيرة بما يكفي هي $\pi/\vert \vert H \vert \vert$، وأنه يُفضَّل إلى حدٍّ ما تقليل تقدير هذه القيمة بدلاً من تضخيمه، إذ يمكن لتضخيم القيمة أن يسمح لمساهمات الحالات عالية الطاقة بإفساد حتى الحالة المثلى في فضاء كريلوف. من ناحية أخرى، اختيار $dt$ صغيراً جداً يؤدي إلى تكييف أسوأ للفضاء الجزئي لكريلوف، إذ تتباين متجهات أساس كريلوف بدرجة أقل من خطوة زمنية إلى أخرى.

# Get Hamiltonian restricted to single-particle states
single_particle_H = np.zeros((n_qubits, n_qubits))
for i in range(n_qubits):
    for j in range(i + 1):
        for p, coeff in H_op.to_list():
            p_x = Pauli(p).x
            p_z = Pauli(p).z
            if all(
                p_x[k] == ((i == k) + (j == k)) % 2 for k in range(n_qubits)
            ):
                sgn = (
                    (-1j) ** sum(p_z[k] and p_x[k] for k in range(n_qubits))
                ) * ((-1) ** p_z[i])
            else:
                sgn = 0
            single_particle_H[i, j] += sgn * coeff
for i in range(n_qubits):
    for j in range(i + 1, n_qubits):
        single_particle_H[i, j] = np.conj(single_particle_H[j, i])

# Set dt according to spectral norm
dt = np.pi / np.linalg.norm(single_particle_H, ord=2)
dt

np.float64(0.10833078115826875)

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

# Set parameters for quantum Krylov algorithm
krylov_dim = 5  # size of Krylov subspace
num_trotter_steps = 6
dt_circ = dt / num_trotter_steps

تحضير الحالة

اختر حالة مرجعية $\vert \psi \rangle$ تمتلك بعض التداخل مع الحالة الأساسية. لهذا الهاميلتوني، نستخدم حالة ذات إثارة في الكيوبت الأوسط $\vert 00..010...00 \rangle$ كحالتنا المرجعية.

qc_state_prep = QuantumCircuit(n_qubits)
qc_state_prep.x(int(n_qubits / 2) + 1)
qc_state_prep.draw("mpl", scale=0.5)

التطور الزمني

يمكننا تحقيق عامل التطور الزمني المولَّد بواسطة هاميلتوني معطى: $U=e^{-iHt}$ عبر تقريب Lie-Trotter.

t = Parameter("t")

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=t, synthesis=LieTrotter(reps=num_trotter_steps)
)

qr = QuantumRegister(n_qubits)
qc_evol = QuantumCircuit(qr)
qc_evol.append(evol_gate, qargs=qr)

<qiskit.circuit.instructionset.InstructionSet at 0x11eef9be0>

اختبار Hadamard

$$\begin{equation*} |0\rangle|0\rangle^N \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle + |1\rangle \Big)|0\rangle^N \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle|0\rangle^N+|1\rangle |\psi_i\rangle\Big) \quad\longrightarrow\quad \frac{1}{\sqrt{2}}\Big(|0\rangle |0\rangle^N+|1\rangle P |\psi_i\rangle\Big) \quad\longrightarrow\quad\frac{1}{\sqrt{2}}\Big(|0\rangle |\psi_j\rangle+|1\rangle P|\psi_i\rangle\Big) \end{equation*}$$

حيث $P$ هو أحد حدود تحليل الهاميلتوني $H=\sum P$ و $\text{Prep} \; \psi_i$ و $\text{Prep} \; \psi_j$ هي عمليات مشروطة تُعدّ متجهات $|\psi_i\rangle$ و $|\psi_j\rangle$ من فضاء كريلوف الأحادي، مع $|\psi_k\rangle = e^{-i H k dt } \vert \psi \rangle = e^{-i H k dt } U_{\psi} \vert 0 \rangle^N$. لقياس $X$، نطبّق أولاً $H$...

$$\begin{equation*} \longrightarrow\quad\frac{1}{2}|0\rangle\Big( |\psi_j\rangle + P|\psi_i\rangle\Big) + \frac{1}{2}|1\rangle\Big(|\psi_j\rangle - P|\psi_i\rangle\Big) \end{equation*}$$

... ثم القياس:

$$\begin{equation*} \begin{split} \Rightarrow\quad\langle X\rangle &= \frac{1}{4}\Bigg(\Big\|| \psi_j\rangle + P|\psi_i\rangle \Big\|^2-\Big\||\psi_j\rangle - P|\psi_i\rangle\Big\|^2\Bigg) \\ &= \text{Re}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{split} \end{equation*}$$

من المتطابقة $|a + b\|^2 = \langle a + b | a + b \rangle = \|a\|^2 + \|b\|^2 + 2\text{Re}\langle a | b \rangle$. وبالمثل، يُعطي قياس $Y$

$$\begin{equation*} \langle Y\rangle = \text{Im}\Big[\langle\psi_j| P|\psi_i\rangle\Big]. \end{equation*}$$

## Create the time-evo op circuit
evol_gate = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)

## Create the time-evo op dagger circuit
evol_gate_d = PauliEvolutionGate(
    H_op, time=dt, synthesis=LieTrotter(reps=num_trotter_steps)
)
evol_gate_d = evol_gate_d.inverse()

# Put pieces together
qc_reg = QuantumRegister(n_qubits)
qc_temp = QuantumCircuit(qc_reg)
qc_temp.compose(qc_state_prep, inplace=True)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate, qargs=qc_reg)
for _ in range(num_trotter_steps):
    qc_temp.append(evol_gate_d, qargs=qc_reg)
qc_temp.compose(qc_state_prep.inverse(), inplace=True)

# Create controlled version of the circuit
controlled_U = qc_temp.to_gate().control(1)

# Create hadamard test circuit for real part
qr = QuantumRegister(n_qubits + 1)
qc_real = QuantumCircuit(qr)
qc_real.h(0)
qc_real.append(controlled_U, list(range(n_qubits + 1)))
qc_real.h(0)

print(
    "Circuit for calculating the real part of the overlap in S via Hadamard test"
)
qc_real.draw("mpl", fold=-1, scale=0.5)

Circuit for calculating the real part of the overlap in S via Hadamard test

يمكن أن تكون دائرة اختبار Hadamard دائرة عميقة بمجرد تحليلها إلى البوابات الأصلية (وهو ما سيزداد أكثر عند مراعاة طولوجيا الجهاز)

print(
    "Number of layers of 2Q operations",
    qc_real.decompose(reps=2).depth(lambda x: x[0].num_qubits == 2),
)

Number of layers of 2Q operations 112753

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

اختبار Hadamard الفعّال

يمكننا تحسين الدوائر العميقة لاختبار Hadamard الذي حصلنا عليه من خلال إدخال بعض التقريبات والاعتماد على بعض الافتراضات حول هاميلتوني النموذج. على سبيل المثال، انظر الدائرة التالية لاختبار Hadamard:

افترض أنه يمكننا حساب $E_0$ كلاسيكياً، وهو القيمة الذاتية لـ $|0\rangle^N$ تحت هاملتوني $H$. يتحقق هذا الشرط عندما يحافظ الهاملتوني على تماثل U(1). وعلى الرغم من أن هذا قد يبدو افتراضاً قوياً، إلا أن هناك حالات كثيرة يمكن فيها افتراض وجود حالة فراغ (تُرسَم على الحالة $|0\rangle^N$ في هذه الحالة) لا يؤثر فيها فعل الهاملتوني. ينطبق هذا على سبيل المثال على هاملتونيات الكيمياء التي تصف جزيءً مستقراً (حيث يكون عدد الإلكترونات محفوظاً). نظراً لأن البوابة $\text{Prep} \; \psi$ تُهيئ حالة المرجع المطلوبة $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$، فمثلاً لتهيئة حالة HF في الكيمياء تكون $\text{Prep} \; \psi$ حاصل ضرب بوابات NOT أحادية الكيوبت، وبالتالي فإن $\text{Prep} \; \psi$ المتحكَّم بها هو مجرد حاصل ضرب بوابات CNOT. وعندئذ تُطبّق الدارة الكهربية أعلاه الحالة التالية قبل القياس:

$$\begin{equation} \begin{split} \ket{0} \ket{0}^N\xrightarrow{H}&\frac{1}{\sqrt{2}} \left( \ket{0}\ket{0}^N+ \ket{1} \ket{0}^N \right)\\ \xrightarrow{\text{1-ctrl-init}}&\frac{1}{\sqrt{2}}\left(|0\rangle|0\rangle^N+|1\rangle|\psi\rangle\right)\\ \xrightarrow{U}&\frac{1}{\sqrt{2}}\left(e^{i\phi}\ket{0}\ket{0}^N+\ket{1} U\ket{\psi}\right)\\ \xrightarrow{\text{0-ctrl-init}}&\frac{1}{\sqrt{2}} \left( e^{i\phi}\ket{0} \ket{\psi} +\ket{1} U\ket{\psi} \right)\\ =&\frac{1}{2} \left( \ket{+}\left(e^{i\phi}\ket{\psi}+U\ket{\psi}\right) +\ket{-}\left(e^{i\phi}\ket{\psi}-U\ket{\psi}\right) \right)\\ =&\frac{1}{2} \left( \ket{+i}\left(e^{i\phi}\ket{\psi}-iU\ket{\psi}\right) +\ket{-i}\left(e^{i\phi}\ket{\psi}+iU\ket{\psi}\right) \right) \end{split} \end{equation}$$

حيث استخدمنا إزاحة الطور القابلة للمحاكاة كلاسيكياً $U\ket{0}^N = e^{i\phi}\ket{0}^N$ في السطر الثالث. وبالتالي تُحسب قيم التوقع على النحو التالي:

$$\begin{equation} \begin{split} \langle X\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi}+\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}+U\ket{\psi}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi}-\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}-U\ket{\psi}\right) \Big)\\ &=\text{Re}\left[e^{-i\phi}\bra{\psi}PU\ket{\psi}\right], \end{split} \end{equation}$$ $$\begin{equation} \begin{split} \langle Y\otimes P\rangle&=\frac{1}{4} \Big( \left(e^{-i\phi}\bra{\psi}+i\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}-iU\ket{\psi}\right) \\ &\qquad-\left(e^{-i\phi}\bra{\psi}-i\bra{\psi}U^\dagger\right)P\left(e^{i\phi}\ket{\psi}+iU\ket{\psi}\right) \Big)\\ &=\text{Im}\left[e^{-i\phi}\bra{\psi}PU\ket{\psi}\right]. \end{split} \end{equation}$$

باستخدام هذه الافتراضات، تمكّنّا من كتابة قيم توقع المؤثرات المعنية باستخدام عدد أقل من العمليات المتحكَّم بها. في الواقع، نحتاج فقط إلى تطبيق تهيئة الحالة المتحكَّم بها $\text{Prep} \; \psi$ دون الحاجة إلى تطورات زمنية متحكَّم بها. إن إعادة صياغة حساباتنا على النحو أعلاه ستُمكّننا من تقليص عمق الدارات الناتجة بشكل كبير.

تحليل مؤثر التطور الزمني باستخدام تحليل Trotter

بدلاً من تطبيق مؤثر التطور الزمني بصورة دقيقة، يمكننا استخدام تحليل Trotter لتطبيق تقريب له. إن تكرار تحليل Trotter بترتيب معين عدة مرات يُتيح لنا تقليصاً إضافياً للخطأ الناجم عن التقريب. فيما يلي، نبني تطبيق Trotter مباشرةً بالطريقة الأكثر كفاءةً لبيان تفاعل الهاملتوني الذي ندرسه (تفاعلات الجيران الأقرب فحسب). عملياً، نُدرج دورانات Pauli وهي $R_{xx}$، $R_{yy}$، $R_{zz}$ بزاوية مُعامَلة $t$ تقابل التطبيق التقريبي لـ $e^{-i (XX + YY + ZZ) t}$. نظراً للاختلاف في تعريف دورانات Pauli والتطور الزمني الذي نسعى إلى تطبيقه، سيتعين علينا استخدام المعامل $2*dt$ لتحقيق تطور زمني مقداره $dt$. علاوة على ذلك، نعكس ترتيب العمليات عند الأعداد الفردية من خطوات Trotter، وهو ما يكافئ وظيفياً الترتيب الأصلي لكنه يُتيح تجميع العمليات المتجاورة في وحدة $SU(2)$ موحدة. ينتج عن ذلك دارة أقل عمقاً بكثير مما يُنتجه استخدام الوظيفة العامة PauliEvolutionGate().

t = Parameter("t")

# Create instruction for rotation about XX+YY-ZZ:
Rxyz_circ = QuantumCircuit(2)
Rxyz_circ.rxx(t, 0, 1)
Rxyz_circ.ryy(t, 0, 1)
Rxyz_circ.rzz(t, 0, 1)
Rxyz_instr = Rxyz_circ.to_instruction(label="RXX+YY+ZZ")

interaction_list = [
    [[i, i + 1] for i in range(0, n_qubits - 1, 2)],
    [[i, i + 1] for i in range(1, n_qubits - 1, 2)],
]  # linear chain

qr = QuantumRegister(n_qubits)
trotter_step_circ = QuantumCircuit(qr)
for i, color in enumerate(interaction_list):
    for interaction in color:
        trotter_step_circ.append(Rxyz_instr, interaction)
    if i < len(interaction_list) - 1:
        trotter_step_circ.barrier()
reverse_trotter_step_circ = trotter_step_circ.reverse_ops()

qc_evol = QuantumCircuit(qr)
for step in range(num_trotter_steps):
    if step % 2 == 0:
        qc_evol = qc_evol.compose(trotter_step_circ)
    else:
        qc_evol = qc_evol.compose(reverse_trotter_step_circ)

qc_evol.decompose().draw("mpl", fold=-1, scale=0.5)

استخدام دارة مُحسَّنة لتهيئة الحالة

control = 0
excitation = int(n_qubits / 2) + 1
controlled_state_prep = QuantumCircuit(n_qubits + 1)
controlled_state_prep.cx(control, excitation)
controlled_state_prep.draw("mpl", fold=-1, scale=0.5)

دارات النماذج لحساب عناصر المصفوفة لـ $\tilde{S}$ و $\tilde{H}$ عبر اختبار Hadamard

الفارق الوحيد بين الدارات المستخدمة في اختبار Hadamard هو الطور في مؤثر التطور الزمني والمتغيرات التي يجري قياسها. لذلك يمكننا إعداد دارة نموذج تمثّل الدارة العامة لاختبار Hadamard، مع عناصر نائبة للبوابات التي تعتمد على مؤثر التطور الزمني.

# Parameters for the template circuits
parameters = []
for idx in range(1, krylov_dim):
    parameters.append(2 * dt_circ * (idx))

# Create modified hadamard test circuit
qr = QuantumRegister(n_qubits + 1)
qc = QuantumCircuit(qr)
qc.h(0)
qc.compose(controlled_state_prep, list(range(n_qubits + 1)), inplace=True)
qc.barrier()
qc.compose(qc_evol, list(range(1, n_qubits + 1)), inplace=True)
qc.barrier()
qc.x(0)
qc.compose(
    controlled_state_prep.inverse(), list(range(n_qubits + 1)), inplace=True
)
qc.x(0)

qc.decompose().draw("mpl", fold=-1)

print(
    "The optimized circuit has 2Q gates depth: ",
    qc.decompose().decompose().depth(lambda x: x[0].num_qubits == 2),
)

The optimized circuit has 2Q gates depth:  74

لقد قلّصنا عمق اختبار Hadamard تقليصاً ملحوظاً من خلال الجمع بين تقريب Trotter والوحدات غير المتحكَّم بها.

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

إنشاء نسخة الـ backend وضبط معاملات وقت التشغيل

service = QiskitRuntimeService()
backend = service.least_busy(operational=True, simulator=False)
if (
    "if_else" not in backend.target.operation_names
):  # Needed as "op_name" could be "if_else"
    backend.target.add_instruction(IfElseOp, name="if_else")
print(backend.name)

التحويل إلى QPU

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

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

cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    t2 = target.qubit_properties[q].t2 * 1e6
    if meas_err > 0.02 or t2 < 100:
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except:
                    continue

for q in cmap_list:
    op_name = list(target.operation_names_for_qargs(q))[0]
    twoq_gate_err = target[f"{op_name}"][q].error
    if twoq_gate_err > 0.005:
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q)
                except:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)
cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates,
    coupling_map=cust_cmap,
)

ثم نُحوَّل الدارة الافتراضية إلى أفضل تخطيط فيزيائي في هذا الهدف الجديد

basis_gates = list(target.operation_names)
pm = generate_preset_pass_manager(
    optimization_level=3,
    target=cust_target,
    basis_gates=basis_gates,
)

qc_trans = pm.run(qc)

print("depth", qc_trans.depth(lambda x: x[0].num_qubits == 2))
print("num 2q ops", qc_trans.count_ops())
print(
    "physical qubits",
    sorted(
        [
            idx
            for idx, qb in qc_trans.layout.initial_layout.get_physical_bits().items()
            if qb._register.name != "ancilla"
        ]
    ),
)

depth 52
num 2q ops OrderedDict([('rz', 2058), ('sx', 1703), ('cz', 728), ('x', 84), ('barrier', 8)])
physical qubits [91, 92, 93, 94, 95, 98, 99, 108, 109, 110, 111, 113, 114, 115, 119, 127, 132, 133, 134, 135, 137, 139, 147, 148, 149, 150, 151, 152, 153, 154, 155]

إنشاء PUBs للتنفيذ مع Estimator

# Define observables to measure for S
observable_S_real = "I" * (n_qubits) + "X"
observable_S_imag = "I" * (n_qubits) + "Y"

observable_op_real = SparsePauliOp(
    observable_S_real
)  # define a sparse pauli operator for the observable
observable_op_imag = SparsePauliOp(observable_S_imag)

layout = qc_trans.layout  # get layout of transpiled circuit
observable_op_real = observable_op_real.apply_layout(
    layout
)  # apply physical layout to the observable
observable_op_imag = observable_op_imag.apply_layout(layout)
observable_S_real = (
    observable_op_real.paulis.to_labels()
)  # get the label of the physical observable
observable_S_imag = observable_op_imag.paulis.to_labels()

observables_S = [[observable_S_real], [observable_S_imag]]

# Define observables to measure for H
# Hamiltonian terms to measure
observable_list = []
for pauli, coeff in zip(H_op.paulis, H_op.coeffs):
    # print(pauli)
    observable_H_real = pauli[::-1].to_label() + "X"
    observable_H_imag = pauli[::-1].to_label() + "Y"
    observable_list.append([observable_H_real])
    observable_list.append([observable_H_imag])

layout = qc_trans.layout

observable_trans_list = []
for observable in observable_list:
    observable_op = SparsePauliOp(observable)
    observable_op = observable_op.apply_layout(layout)
    observable_trans_list.append([observable_op.paulis.to_labels()])

observables_H = observable_trans_list

# Define a sweep over parameter values
params = np.vstack(parameters).T

# Estimate the expectation value for all combinations of
# observables and parameter values, where the pub result will have
# shape (# observables, # parameter values).
pub = (qc_trans, observables_S + observables_H, params)

تشغيل الدارات

دارات $t=0$ قابلة للحساب كلاسيكياً

qc_cliff = qc.assign_parameters({t: 0})

# Get expectation values from experiment
S_expval_real = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "X")
)
S_expval_imag = StabilizerState(qc_cliff).expectation_value(
    Pauli("I" * (n_qubits) + "Y")
)

# Get expectation values
S_expval = S_expval_real + 1j * S_expval_imag

H_expval = 0
for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Get expectation values from experiment
    expval_real = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "X")
    )
    expval_imag = StabilizerState(qc_cliff).expectation_value(
        Pauli(pauli[::-1].to_label() + "Y")
    )
    expval = expval_real + 1j * expval_imag

    # Fill-in matrix elements
    H_expval += coeff * expval

print(H_expval)

(25+0j)

تنفيذ دارات $S$ و $\tilde{H}$ باستخدام Estimator

# Experiment options
num_randomizations = 300
num_randomizations_learning = 30
shots_per_randomization = 100
noise_factors = [1, 1.2, 1.4]
learning_pair_depths = [0, 4, 24, 48]

experimental_opts = {}
experimental_opts["resilience"] = {
    "measure_mitigation": True,
    "measure_noise_learning": {
        "num_randomizations": num_randomizations_learning,
        "shots_per_randomization": shots_per_randomization,
    },
    "zne_mitigation": True,
    "zne": {"noise_factors": noise_factors},
    "layer_noise_learning": {
        "max_layers_to_learn": 10,
        "layer_pair_depths": learning_pair_depths,
        "shots_per_randomization": shots_per_randomization,
        "num_randomizations": num_randomizations_learning,
    },
    "zne": {
        "amplifier": "pea",
        "extrapolated_noise_factors": [0] + noise_factors,
    },
}
experimental_opts["twirling"] = {
    "num_randomizations": num_randomizations,
    "shots_per_randomization": shots_per_randomization,
    "strategy": "all",
}

estimator = Estimator(mode=backend, options=experimental_opts)

job = estimator.run([pub])

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

results = job.result()[0]

حساب مصفوفة الهاملتوني الفعّال ومصفوفة التراكب

أولاً، نحسب الطور المتراكم لدى الحالة $\vert 0 \rangle$ خلال التطور الزمني غير المتحكَّم به

prefactors = [
    np.exp(-1j * sum([c for p, c in H_op.to_list() if "Z" in p]) * i * dt)
    for i in range(1, krylov_dim)
]

بمجرد الحصول على نتائج تنفيذ الدارات، يمكننا المعالجة اللاحقة للبيانات لحساب عناصر مصفوفة $S$

# Assemble S, the overlap matrix of dimension D:
S_first_row = np.zeros(krylov_dim, dtype=complex)
S_first_row[0] = 1 + 0j

# Add in ancilla-only measurements:
for i in range(krylov_dim - 1):
    # Get expectation values from experiment

```python
    expval_real = results.data.evs[0][0][
        i
    ]  # automatic extrapolated evs if ZNE is used
    expval_imag = results.data.evs[1][0][
        i
    ]  # automatic extrapolated evs if ZNE is used

    # Get expectation values
    expval = expval_real + 1j * expval_imag
    S_first_row[i + 1] += prefactors[i] * expval

S_first_row_list = S_first_row.tolist()  # for saving purposes

S_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        S_circ[j, i] = S_first_row[i - j]
    else:
        S_circ[j, i] = np.conj(S_first_row[j - i])

Matrix(S_circ)

$$\displaystyle \left[\begin{matrix}1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i & 0.0012070853532697 + 0.312052218182462 i\\-0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i & -0.180546747798251 - 0.492624093654174 i\\0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i & 0.467051960502366 + 0.516197865254034 i\\-0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0 & -0.723052998582984 - 0.345085413575966 i\\0.0012070853532697 - 0.312052218182462 i & -0.180546747798251 + 0.492624093654174 i & 0.467051960502366 - 0.516197865254034 i & -0.723052998582984 + 0.345085413575966 i & 1.0\end{matrix}\right]$$

وعناصر المصفوفة $\tilde{H}$

# Assemble S, the overlap matrix of dimension D:
H_first_row = np.zeros(krylov_dim, dtype=complex)
H_first_row[0] = H_expval

for obs_idx, (pauli, coeff) in enumerate(zip(H_op.paulis, H_op.coeffs)):
    # Add in ancilla-only measurements:
    for i in range(krylov_dim - 1):
        # Get expectation values from experiment
        expval_real = results.data.evs[2 + 2 * obs_idx][0][
            i
        ]  # automatic extrapolated evs if ZNE is used
        expval_imag = results.data.evs[2 + 2 * obs_idx + 1][0][
            i
        ]  # automatic extrapolated evs if ZNE is used

        # Get expectation values
        expval = expval_real + 1j * expval_imag
        H_first_row[i + 1] += prefactors[i] * coeff * expval

H_first_row_list = H_first_row.tolist()

H_eff_circ = np.zeros((krylov_dim, krylov_dim), dtype=complex)

# Distribute entries from first row across matrix:
for i, j in it.product(range(krylov_dim), repeat=2):
    if i >= j:
        H_eff_circ[j, i] = H_first_row[i - j]
    else:
        H_eff_circ[j, i] = np.conj(H_first_row[j - i])

Matrix(H_eff_circ)

$$\displaystyle \left[\begin{matrix}25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i & 1.98818301405581 + 5.8897614762563 i\\-14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i & -5.15587257589417 - 8.88280836036843 i\\10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i & 10.2857217968584 + 9.0431912203186 i\\-5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0 & -14.2437089383409 - 6.50486277982165 i\\1.98818301405581 - 5.8897614762563 i & -5.15587257589417 + 8.88280836036843 i & 10.2857217968584 - 9.0431912203186 i & -14.2437089383409 + 6.50486277982165 i & 25.0\end{matrix}\right]$$

أخيرًا، يمكننا حل مسألة القيم الذاتية المعممة لـ $\tilde{H}$:

$\tilde{H} \vec{c} = c S \vec{c}$

والحصول على تقدير لطاقة الحالة الأساسية $c_{min}$

gnd_en_circ_est_list = []
for d in range(1, krylov_dim + 1):
    # Solve generalized eigenvalue problem for different size of the Krylov space
    gnd_en_circ_est = solve_regularized_gen_eig(
        H_eff_circ[:d, :d], S_circ[:d, :d], threshold=9e-1
    )
    gnd_en_circ_est_list.append(gnd_en_circ_est)
    print("The estimated ground state energy is: ", gnd_en_circ_est)

The estimated ground state energy is:  25.0
The estimated ground state energy is:  22.572154819954875
The estimated ground state energy is:  21.691509219286587
The estimated ground state energy is:  21.23882298756386
The estimated ground state energy is:  20.965499325470294

بالنسبة لقطاع الجسيم الواحد، يمكننا حساب الحالة الأساسية لهذا القطاع من Hamiltonian بكفاءة باستخدام الأساليب الكلاسيكية

gs_en = single_particle_gs(H_op, n_qubits)

n_sys_qubits 30
n_exc 1 , subspace dimension 31
single particle ground state energy:  21.021912418526906

plt.plot(
    range(1, krylov_dim + 1),
    gnd_en_circ_est_list,
    color="blue",
    linestyle="-.",
    label="KQD estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [gs_en] * krylov_dim,
    color="red",
    linestyle="-",
    label="exact",
)
plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))
plt.legend()
plt.xlabel("Krylov space dimension")
plt.ylabel("Energy")
plt.title(
    "Estimating Ground state energy with Krylov Quantum Diagonalization"
)
plt.show()

ملحق: الفضاء الجزئي لـ Krylov من التطورات الزمنية الحقيقية

يُعرَّف فضاء Krylov الأحادي (unitary Krylov space) بالصيغة التالية:

$$\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ |\psi\rangle, e^{-iH\,dt} |\psi\rangle, \dots, e^{-irH\,dt} |\psi\rangle \right\}$$

لخطوة زمنية $dt$ سنحددها لاحقًا. لنفترض مؤقتًا أن $r$ عددٌ زوجي، ونعرّف $d=r/2$. لاحظ أنه عند إسقاط Hamiltonian على فضاء Krylov أعلاه، فإنه لا يمكن تمييزه عن فضاء Krylov

$$\mathcal{K}_U(H, |\psi\rangle) = \text{span}\left\{ e^{i\,d\,H\,dt}|\psi\rangle, e^{i(d-1)H\,dt} |\psi\rangle, \dots, e^{-i(d-1)H\,dt} |\psi\rangle, e^{-i\,d\,H\,dt} |\psi\rangle \right\},$$

أي حيث تُزاح جميع التطورات الزمنية للخلف بمقدار $d$ خطوة زمنية. والسبب في عدم إمكانية التمييز هو أن عناصر المصفوفة

$$\tilde{H}_{j,k} = \langle\psi|e^{i\,j\,H\,dt}He^{-i\,k\,H\,dt}|\psi\rangle=\langle\psi|He^{i(j-k)H\,dt}|\psi\rangle$$

تبقى ثابتة تحت الإزاحات الكلية لزمن التطور، إذ إن التطورات الزمنية تتبادل مع Hamiltonian. وبالنسبة للأعداد الفردية $r$، يمكننا استخدام تحليل $r-1$.

نريد أن نُثبت أنه في مكان ما داخل فضاء Krylov هذا، تُضمَن دائمًا وجود حالة منخفضة الطاقة. ونفعل ذلك عن طريق النتيجة التالية، المشتقة من النظرية 3.1 في [3]:

الادعاء 1: توجد دالة $f$ بحيث أنه للطاقات $E$ ضمن النطاق الطيفي لـ Hamiltonian (أي بين طاقة الحالة الأساسية والطاقة القصوى)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ لجميع قيم $E$ التي تبعد مسافة $\ge\delta$ عن $E_0$، أي أنها تتناقص أسيًا
3. $f(E)$ هي تركيبة خطية من $e^{ijE\,dt}$ بحيث $j=-d,-d+1,...,d-1,d$

نُقدّم إثباتًا لذلك أدناه، غير أنه يمكن تخطيه بأمان ما لم يُرد المرء فهم الحجة الكاملة الصارمة. نركّز الآن على مضامين الادعاء أعلاه. بموجب الخاصية 3 أعلاه، يتضح أن فضاء Krylov المُزاح أعلاه يحتوي على الحالة $f(H)|\psi\rangle$. وهذه هي حالتنا منخفضة الطاقة. لمعرفة السبب، نكتب $|\psi\rangle$ في أساس القيم الذاتية للطاقة:

$$|\psi\rangle = \sum_{k=0}^{N}\gamma_k|E_k\rangle,$$

حيث $|E_k\rangle$ هي الحالة الذاتية للطاقة رقم $k$، و$\gamma_k$ هي سعتها في الحالة الابتدائية $|\psi\rangle$. بالتعبير عن ذلك، تُعطَى $f(H)|\psi\rangle$ بالصيغة

$$f(H)|\psi\rangle = \sum_{k=0}^{N}\gamma_kf(E_k)|E_k\rangle,$$

مستعملين حقيقة أنه يمكن استبدال $H$ بـ $E_k$ حين تؤثر على الحالة الذاتية $|E_k\rangle$. وبالتالي، يكون خطأ الطاقة لهذه الحالة هو

$$\text{energy error} = \frac{\langle\psi|f(H)(H-E_0)f(H)|\psi\rangle}{\langle\psi|f(H)^2|\psi\rangle}$$ $$= \frac{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.$$

لتحويل هذا إلى حدٍّ أعلى أسهل في الفهم، نفصل أولًا مجموع البسط إلى حدود تحقق $E_k-E_0\le\delta$ وحدود تحقق $E_k-E_0>\delta$:

$$\text{energy error} = \frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} + \frac{\sum_{E_k> E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2}.$$

يمكننا تقييد الحد الأول بحد أعلى قيمته $\delta$،

$$\frac{\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0)}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} < \frac{\delta\sum_{E_k\le E_0+\delta}|\gamma_k|^2f(E_k)^2}{\sum_{k=0}^{N}|\gamma_k|^2f(E_k)^2} \le \delta,$$

حيث تنبثق الخطوة الأولى من كون $E_k-E_0\le\delta$ لكل $E_k$ في المجموع، وتنبثق الخطوة الثانية من كون المجموع في البسط جزءًا من المجموع في المقام. أما الحد الثاني، فنُقيّد المقام بحد أدنى هو $|\gamma_0|^2$، إذ إن $f(E_0)^2=1$: وبجمع كل شيء معًا، يعطينا ذلك

$$\text{energy error} \le \delta + \frac{1}{|\gamma_0|^2}\sum_{E_k>E_0+\delta}|\gamma_k|^2f(E_k)^2(E_k-E_0).$$

لتبسيط ما تبقى، نلاحظ أنه لجميع هذه القيم $E_k$، يستتبع من تعريف $f$ أن $f(E_k)^2 \le 4\left(1 + \delta\right)^{-2d}$. وبتقييد $E_k-E_0<2\|H\|$ بحد أعلى، وتقييد $\sum_{E_k>E_0+\delta}|\gamma_k|^2<1$ بحد أعلى، نحصل على

$$\text{energy error} \le \delta + \frac{8}{|\gamma_0|^2}\|H\|\left(1 + \delta\right)^{-2d}.$$

يسري هذا لأي $\delta>0$، لذا إذا ضبطنا $\delta$ مساويًا لهدف الخطأ المنشود، فإن الحد الأعلى للخطأ يتقارب نحوه أسيًا مع بُعد فضاء Krylov $2d=r$. ولاحظ أيضًا أنه إذا كان $\delta<E_1-E_0$ فإن حد $\delta$ يختفي كليًا من الحد الأعلى أعلاه.

لإتمام الحجة، نشير أولًا إلى أن ما سبق ليس سوى خطأ الطاقة للحالة المحددة $f(H)|\psi\rangle$، لا خطأ الطاقة لأدنى حالة طاقة في فضاء Krylov. غير أنه بموجب مبدأ التنويع (Rayleigh-Ritz)، يكون خطأ الطاقة لأدنى حالة طاقة في فضاء Krylov مقيدًا بحد أعلى هو خطأ الطاقة لأي حالة في فضاء Krylov، وبالتالي فإن ما سبق يُعدّ أيضًا حدًا أعلى لخطأ الطاقة لأدنى حالة طاقة، أي ناتج خوارزمية Krylov quantum diagonalization.

يمكن إجراء تحليل مشابه لما سبق يأخذ في الحسبان أيضًا الضوضاء وإجراء العتبة (thresholding) المُناقَش في هذا الكتيب. راجع [2] و[4] لهذا التحليل.

ملحق: إثبات الادعاء 1

يُشتق ما يلي في معظمه من [3]، النظرية 3.1: ليكن $0 < a < b$ وليكن $\Pi^*_d$ فضاء المتعددات الحدية المتبقية (polynomials whose value at 0 is 1) ذات الدرجة الأكثر $d$. الحل لـ

$$\beta(a, b, d) = \min_{p \in \Pi^*_d} \max_{x \in [a, b]} |p(x)| \quad$$

هو

$$p^*(x) = \frac{T_d\left(\frac{b + a - 2x}{b - a}\right)}{T_d\left(\frac{b + a}{b - a}\right)}, \quad$$

والقيمة الدنيا المقابلة هي

$$\beta(a, b, d) = T_d^{-1}\left(\frac{b + a}{b - a}\right).$$

نريد تحويل هذا إلى دالة يمكن التعبير عنها بصورة طبيعية من خلال الأسيات المركبة، لأن تلك هي التطورات الزمنية الحقيقية التي تولّد فضاء Krylov الكمي. لتحقيق ذلك، من المناسب تقديم التحويل التالي للطاقات ضمن النطاق الطيفي لـ Hamiltonian إلى أعداد في المدى $[0,1]$: نعرّف

$$g(E) = \frac{1-\cos\big((E-E_0)dt\big)}{2},$$

حيث $dt$ خطوة زمنية بحيث $-\pi < E_0dt < E_\text{max}dt < \pi$. لاحظ أن $g(E_0)=0$ وأن $g(E)$ تتزايد كلما ابتعدت $E$ عن $E_0$.

الآن باستخدام المتعددة الحدية $p^*(x)$ مع الضبط $a = g(E_0 + \delta)$، $b = 1$، و d = int(r/2)، نعرّف الدالة:

$$f(E) = p^* \left( g(E) \right) = \frac{T_d\left(1 + 2\frac{\cos\big((E-E_0)dt\big) - \cos\big(\delta\,dt\big)}{1 +\cos\big(\delta\,dt\big)}\right)}{T_d\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)}$$

حيث $E_0$ هي طاقة الحالة الأساسية. يمكننا أن نرى بإدراج $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ أن $f(E)$ هي متعددة حدية مثلثية من الدرجة $d$، أي تركيبة خطية من $e^{ijE\,dt}$ بحيث $j=-d,-d+1,...,d-1,d$. فضلًا عن ذلك، من تعريف $p^*(x)$ أعلاه لدينا أن $f(E_0)=p(0)=1$، وأنه لأي $E$ في النطاق الطيفي بحيث $\vert E-E_0 \vert > \delta$ يكون

$$|f(E)| \le \beta(a, b, d) = T_d^{-1}\left(1 + 2\frac{1-\cos\big(\delta\,dt\big)}{1 + \cos\big(\delta\,dt\big)}\right)$$ $$\leq 2\left(1 + \delta\right)^{-d} = 2\left(1 + \delta\right)^{-\lfloor k/2\rfloor}.$$

المراجع

[1] N. Yoshioka, M. Amico, W. Kirby et al. "Diagonalization of large many-body Hamiltonians on a quantum processor". arXiv:2407.14431

[2] Ethan N. Epperly, Lin Lin, and Yuji Nakatsukasa. "A theory of quantum subspace diagonalization". SIAM Journal on Matrix Analysis and Applications 43, 1263–1290 (2022).

[3] Å. Björck. "Numerical methods in matrix computations". Texts in Applied Mathematics. Springer International Publishing. (2014).

[4] William Kirby. "Analysis of quantum Krylov algorithms with errors". Quantum 8, 1457 (2024).

استبيان البرنامج التعليمي

يُرجى إجراء هذا الاستبيان القصير لتقديم ملاحظاتك حول هذا البرنامج التعليمي. ستساعدنا آراؤك في تحسين محتوانا وتجربة المستخدمين.