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

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

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

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

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

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

لنعتبر هاميلتونية هايزنبرغ لـ $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}$ عبر تقريب لي-تروتر.

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>

اختبار هادامار

$$\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

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

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: تحسين المشكلة لتنفيذها على العتاد الكمي

اختبار هادامار الفعال

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

نفترض أنه يمكننا حساب $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$ وليس تطورات زمنية متحكم بها. إعادة صياغة حساباتنا كما ورد أعلاه ستسمح لنا بتقليل عمق الدوائر الناتجة بشكل كبير.

تفكيك مؤثر التطور الزمني بتفكيك تروتر

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

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

# 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

لقد قللنا بشكل كبير من عمق اختبار هادامار بمزيج من تقريب تروتر والعمليات الوحدوية غير المتحكم بها

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

إنشاء مثيل للخلفية وتعيين معاملات وقت التشغيل

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)

التحويل إلى وحدة معالجة كمية

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

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

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

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()

الملحق: فضاء كريلوف الفرعي من التطورات الزمنية الحقيقية

يُعرَّف فضاء كريلوف الوحدوي كالتالي

$$\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$. لاحظ أنه عندما نسقط الهاميلتونية في فضاء كريلوف أعلاه، فإنه لا يمكن تمييزه عن فضاء كريلوف

$$\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$$

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

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

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

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

لإكمال الحجة، نلاحظ أولاً أن ما سبق هو فقط خطأ الطاقة للحالة المحددة $f(H)|\psi\rangle$، وليس خطأ الطاقة لأدنى حالة طاقة في فضاء كريلوف. ومع ذلك، بموجب مبدأ (رايلي-ريتز) التباينية، فإن خطأ الطاقة لأدنى حالة طاقة في فضاء كريلوف محدود من الأعلى بخطأ الطاقة لأي حالة في فضاء كريلوف، لذا فإن ما سبق هو أيضاً حد أعلى لخطأ الطاقة لأدنى حالة طاقة، أي ناتج خوارزمية التقطير الكمي بطريقة كريلوف.

يمكن إجراء تحليل مماثل لما سبق يأخذ في الاعتبار أيضاً الضوضاء وإجراء العتبة الموصوف في الدفتر. انظر [2] و [4] لهذا التحليل.

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

ما يلي مشتق في الغالب من [3]، المبرهنة 3.1: ليكن $0 < a < b$ وليكن $\Pi^*_d$ فضاء كثيرات الحدود المتبقية (كثيرات الحدود التي قيمتها عند 0 تساوي 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).$$

نريد تحويل هذا إلى دالة يمكن التعبير عنها بشكل طبيعي بدلالة الأسيات المعقدة، لأنها التطورات الزمنية الحقيقية التي تولّد فضاء كريلوف الكمي. للقيام بذلك، من الملائم إدخال التحويل التالي للطاقات ضمن النطاق الطيفي للهاميلتونية إلى أرقام في المجال $[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 وb وd المعيّنة كـ $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).

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

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

رابط الاستبيان