การทแยงเมตริกซ์เชิงควอนตัมแบบ Krylov ของ Lattice Hamiltonians

ประมาณการเวลาใช้งาน: 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

พื้นหลัง

บทช่วยสอนนี้จะสาธิตวิธีการนำอัลกอริทึม Krylov Quantum Diagonalization (KQD) ไปใช้ภายในบริบทของ Qiskit patterns ก่อนอื่นจะได้เรียนรู้ทฤษฎีเบื้องหลังอัลกอริทึม จากนั้นจะเห็นการสาธิตการรันจริงบน QPU

ในหลากหลายสาขาวิชา เราสนใจการหาคุณสมบัติของสถานะพื้น (ground state) ของระบบควอนตัม ตัวอย่างเช่น การทำความเข้าใจธรรมชาติพื้นฐานของอนุภาคและแรง การทำนายและเข้าใจพฤติกรรมของวัสดุที่ซับซ้อน และการทำความเข้าใจอันตรกิริยาและปฏิกิริยาทางชีวเคมี เนื่องจากการเติบโตแบบเอกซ์โปเนนเชียลของ Hilbert space และความสัมพันธ์ที่เกิดขึ้นในระบบที่พัน (entangled) กัน อัลกอริทึมแบบคลาสสิกจึงต้องเผชิญความยากลำบากในการแก้ปัญหานี้สำหรับระบบควอนตัมที่มีขนาดใหญ่ขึ้นเรื่อยๆ ในด้านหนึ่งมีแนวทางที่ใช้ประโยชน์จากฮาร์ดแวร์ควอนตัมโดยเน้นที่วิธีการเชิงแปรผัน (variational quantum methods) (เช่น variational quantum eigensolver) เทคนิคเหล่านี้เผชิญความท้าทายกับอุปกรณ์ในปัจจุบันเนื่องจากต้องเรียกใช้ฟังก์ชันจำนวนมากในกระบวนการหาค่าที่เหมาะสม ซึ่งเพิ่มภาระทรัพยากรอย่างมากเมื่อนำเทคนิคการลดข้อผิดพลาดขั้นสูงมาใช้ ทำให้ประสิทธิภาพของมันจำกัดอยู่กับระบบขนาดเล็ก อีกด้านหนึ่งมีวิธีการควอนตัมที่รองรับข้อผิดพลาด (fault-tolerant) ซึ่งมีการรับประกันประสิทธิภาพ (เช่น quantum phase estimation) ที่ต้องการ Circuit ลึกซึ่งรันได้เฉพาะบนอุปกรณ์ fault-tolerant เท่านั้น ด้วยเหตุผลเหล่านี้ เราจึงนำเสนออัลกอริทึมควอนตัมที่อิงบนวิธีการ subspace (ตามที่อธิบายใน review paper นี้) ซึ่งก็คืออัลกอริทึม Krylov quantum diagonalization (KQD) อัลกอริทึมนี้ทำงานได้ดีในระดับขนาดใหญ่ [1] บนฮาร์ดแวร์ควอนตัมในปัจจุบัน มีการรับประกันประสิทธิภาพที่คล้ายคลึงกับ phase estimation รองรับเทคนิคการลดข้อผิดพลาดขั้นสูง และอาจให้ผลลัพธ์ที่ไม่สามารถเข้าถึงได้ด้วยการคำนวณแบบคลาสสิก

ข้อกำหนดเบื้องต้น

ก่อนเริ่มบทช่วยสอนนี้ ตรวจสอบให้แน่ใจว่าได้ติดตั้งสิ่งต่อไปนี้แล้ว:

• Qiskit SDK v2.0 หรือใหม่กว่า พร้อมรองรับ visualization
• Qiskit Runtime v0.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: แปลงอินพุตแบบคลาสสิกเป็นปัญหาควอนตัม

Krylov space

Krylov space $\mathcal{K}^r$ ลำดับที่ $r$ คือปริภูมิที่ถูก span โดยเวกเตอร์ที่ได้จากการคูณกำลังสูงกว่าของเมตริกซ์ $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$ คือ Hamiltonian $H$ เราจะเรียกปริภูมิที่สอดคล้องกันว่า power Krylov space $\mathcal{K}_P$ ในกรณีที่ $A$ คือตัวดำเนินการวิวัฒนาการตามเวลาที่สร้างขึ้นโดย Hamiltonian $U=e^{-iHt}$ เราจะเรียกปริภูมินั้นว่า unitary Krylov space $\mathcal{K}_U$ power Krylov subspace ที่เราใช้ในการคำนวณแบบคลาสสิกไม่สามารถสร้างขึ้นโดยตรงบนคอมพิวเตอร์ควอนตัม เนื่องจาก $H$ ไม่ใช่ตัวดำเนินการยูนิทารี แต่เราสามารถใช้ตัวดำเนินการวิวัฒนาการตามเวลา $U = e^{-iHt}$ ซึ่งสามารถแสดงให้เห็นว่าให้การรับประกันการลู่เข้าที่คล้ายคลึงกับ power method กำลังของ $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\}$$

ดูภาคผนวกสำหรับการอนุมานโดยละเอียดว่า unitary Krylov space ช่วยให้แสดง eigenstates พลังงานต่ำได้อย่างแม่นยำได้อย่างไร

อัลกอริทึม Krylov quantum diagonalization

เมื่อกำหนด Hamiltonian $H$ ที่เราต้องการทแยงเมตริกซ์ อันดับแรกเราพิจารณา unitary Krylov space $\mathcal{K}_U$ ที่สอดคล้องกัน เป้าหมายคือการหาการแทนค่าแบบกระชับของ Hamiltonian ใน $\mathcal{K}_U$ ซึ่งเราจะเรียกว่า $\tilde{H}$ องค์ประกอบเมตริกซ์ของ $\tilde{H}$ ซึ่งเป็น projection ของ Hamiltonian ลงใน Krylov space สามารถคำนวณได้โดยการคำนวณค่าความคาดหวังต่อไปนี้

$$\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$ คือเวกเตอร์ของ unitary Krylov space และ $t_n = n dt$ คือผลคูณของขั้นเวลา $dt$ ที่เลือก บนคอมพิวเตอร์ควอนตัม การคำนวณองค์ประกอบเมตริกซ์แต่ละตัวสามารถทำได้ด้วยอัลกอริทึมใดก็ตามที่ช่วยให้ได้ค่า overlap ระหว่างสถานะควอนตัม บทช่วยสอนนี้เน้นที่ Hadamard test เนื่องจาก $\mathcal{K}_U$ มีมิติ $r$ Hamiltonian ที่ projected ลงใน subspace จะมีมิติ $r \times r$ ด้วยค่า $r$ ที่เล็กพอ (โดยทั่วไป $r<<100$ เพียงพอที่จะได้การลู่เข้าของการประมาณค่า eigenenergies) เราสามารถทแยงเมตริกซ์ projected Hamiltonian $\tilde{H}$ ได้อย่างง่ายดาย อย่างไรก็ตาม เราไม่สามารถทแยงเมตริกซ์ $\tilde{H}$ โดยตรงได้เนื่องจากความไม่ตั้งฉากกันของเวกเตอร์ Krylov space เราต้องวัดค่า overlap และสร้างเมตริกซ์ $\tilde{S}$

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

สิ่งนี้ช่วยให้เราแก้ปัญหา eigenvalue ในปริภูมิที่ไม่ตั้งฉาก (เรียกอีกอย่างว่า generalized eigenvalue problem)

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

จากนั้นเราสามารถหาค่าประมาณของ eigenvalues และ eigenstates ของ $H$ ได้โดยดูที่ค่าต่างๆ ของ $\tilde{H}$ ตัวอย่างเช่น ค่าประมาณพลังงานสถานะพื้นได้มาจากการนำ eigenvalue ที่เล็กที่สุด $c$ และสถานะพื้นจาก eigenvector $\vec{c}$ ที่สอดคล้องกัน สัมประสิทธิ์ใน $\vec{c}$ กำหนดการมีส่วนร่วมของเวกเตอร์ต่างๆ ที่ span $\mathcal{K}_U$

ภาพแสดง Circuit representation ของ modified Hadamard test ซึ่งเป็นวิธีการที่ใช้คำนวณค่า overlap ระหว่างสถานะควอนตัมต่างๆ สำหรับองค์ประกอบเมตริกซ์แต่ละตัว $\tilde{H}_{i,j}$ จะมีการทำ Hadamard test ระหว่างสถานะ $\vert \psi_i \rangle$ และ $\vert \psi_j \rangle$ ซึ่งเน้นให้เห็นในภาพด้วยรูปแบบสีของส่วนประกอบเมตริกซ์และการดำเนินการ $\text{Prep} \; \psi_i$ และ $\text{Prep} \; \psi_j$ ที่สอดคล้องกัน ดังนั้นจึงต้องมีชุด Hadamard tests สำหรับทุกคู่ที่เป็นไปได้ของเวกเตอร์ Krylov space เพื่อคำนวณองค์ประกอบเมตริกซ์ทั้งหมดของ projected Hamiltonian $\tilde{H}$ สายบนใน Hadamard test circuit คือ ancilla Qubit ที่ถูกวัดในฐาน X หรือ Y ค่าความคาดหวังของมันกำหนดค่า overlap ระหว่างสถานะต่างๆ สายล่างแทน Qubit ทั้งหมดของระบบ Hamiltonian การดำเนินการ $\text{Prep} \; \psi_i$ เตรียม system qubit ในสถานะ $\vert \psi_i \rangle$ โดยควบคุมด้วยสถานะของ ancilla qubit (เช่นเดียวกันสำหรับ $\text{Prep} \; \psi_j$) และการดำเนินการ $P$ แทน Pauli decomposition ของระบบ Hamiltonian $H = \sum_i P_i$ การอนุมานโดยละเอียดเพิ่มเติมของการดำเนินการที่คำนวณโดย Hadamard test ให้ไว้ด้านล่าง

กำหนด Hamiltonian

มาพิจารณา Heisenberg Hamiltonian สำหรับ $N$ Qubits บน linear chain: $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 สำหรับขั้นเวลาด้วยวิธีการฮิวริสติก (อิงจากขอบเขตบนของ norm ของ Hamiltonian) Ref [2] แสดงว่าขั้นเวลาที่เล็กพอคือ $\pi/\vert \vert H \vert \vert$ และว่าการประมาณค่าต่ำกว่าจริงนั้นดีกว่าการประมาณค่าสูงกว่าจริงในระดับหนึ่ง เนื่องจากการประมาณสูงเกินไปอาจทำให้มีส่วนร่วมจากสถานะพลังงานสูงมาทำลายแม้แต่สถานะที่เหมาะสมที่สุดใน Krylov space ในทางกลับกัน การเลือก $dt$ ที่เล็กเกินไปทำให้ conditioning ของ Krylov subspace แย่ลง เนื่องจากเวกเตอร์ฐาน Krylov แตกต่างกันน้อยลงจากขั้นเวลาหนึ่งไปสู่อีกขั้นเวลาหนึ่ง

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

และกำหนดพารามิเตอร์อื่นๆ ของอัลกอริทึม เพื่อวัตถุประสงค์ของบทช่วยสอนนี้ เราจะจำกัดตัวเองให้ใช้ Krylov space ที่มีเพียงห้ามิติ ซึ่งค่อนข้างจำกัด

# 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$ ที่มีค่า overlap บางส่วนกับสถานะพื้น สำหรับ Hamiltonian นี้ เราใช้สถานะที่มีการกระตุ้นที่ Qubit กลาง $\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)

วิวัฒนาการตามเวลา

เราสามารถทำให้ตัวดำเนินการวิวัฒนาการตามเวลาที่สร้างโดย Hamiltonian ที่กำหนด $U=e^{-iHt}$ เป็นจริงได้ผ่านทาง Lie-Trotter approximation

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 test

$$\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$ คือหนึ่งในพจน์ใน decomposition ของ Hamiltonian $H=\sum P$ และ $\text{Prep} \; \psi_i$, $\text{Prep} \; \psi_j$ คือการดำเนินการแบบ controlled ที่เตรียม $|\psi_i\rangle$, $|\psi_j\rangle$ เวกเตอร์ของ unitary Krylov space โดยที่ $|\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

Circuit ของ Hadamard test อาจเป็น Circuit ที่ลึกมากเมื่อ decompose ไปเป็น native gates (ซึ่งจะเพิ่มขึ้นอีกหากคำนึงถึง topology ของอุปกรณ์)

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: ปรับปรุงปัญหาให้เหมาะสมสำหรับการรันบน Quantum hardware

Hadamard test แบบมีประสิทธิภาพ

เราสามารถปรับปรุง Circuit ที่ลึกมากสำหรับ Hadamard test ที่ได้มาได้ โดยการนำการประมาณค่าบางอย่างมาใช้และอาศัยสมมติฐานบางประการเกี่ยวกับ Hamiltonian ของโมเดล ตัวอย่างเช่น ลองพิจารณา Circuit ต่อไปนี้สำหรับ Hadamard test:

สมมติว่าเราสามารถคำนวณ $E_0$ ซึ่งเป็น eigenvalue ของ $|0\rangle^N$ ภายใต้ Hamiltonian $H$ ได้แบบคลาสสิก เงื่อนไขนี้จะเป็นจริงเมื่อ Hamiltonian รักษา U(1) symmetry แม้ว่านี่อาจดูเหมือนเป็นสมมติฐานที่เข้มแข็ง แต่มีหลายกรณีที่ปลอดภัยที่จะสมมติว่ามี vacuum state (ในกรณีนี้จะ map ไปยัง state $|0\rangle^N$) ซึ่งไม่ได้รับผลกระทบจาก Hamiltonian ตัวอย่างเช่น Hamiltonian ทางเคมีที่อธิบายโมเลกุลที่เสถียร (ซึ่งจำนวนอิเล็กตรอนได้รับการอนุรักษ์) เป็นไปตามเงื่อนไขนี้ เมื่อกำหนดว่า Gate $\text{Prep} \; \psi$ เตรียม reference state ที่ต้องการ $\ket{psi} = \text{Prep} \; \psi \ket{0} = e^{-i H 0 dt} U_{\psi} \ket{0}$ เช่น เพื่อเตรียม HF state สำหรับเคมี $\text{Prep} \; \psi$ จะเป็นผลคูณของ single-qubit NOT ดังนั้น controlled-$\text{Prep} \; \psi$ จะเป็นแค่ผลคูณของ CNOT จากนั้น Circuit ด้านบนจะ implement state ต่อไปนี้ก่อนการวัด:

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

โดยเราได้ใช้ phase shift ที่จำลองแบบคลาสสิกได้ $U\ket{0}^N = e^{i\phi}\ket{0}^N$ ในบรรทัดที่สาม ดังนั้น expectation value จึงได้รับเป็น

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

ด้วยสมมติฐานเหล่านี้ เราสามารถเขียน expectation value ของ operator ที่สนใจโดยใช้ controlled operations น้อยลง ที่จริงแล้ว เราต้องการ implement เฉพาะ controlled state preparation $\text{Prep} \; \psi$ เท่านั้น ไม่ใช่ controlled time evolution การกำหนดกรอบการคำนวณใหม่ดังกล่าวจะช่วยให้เราลด depth ของ Circuit ที่ได้ออกมาได้อย่างมาก

แยกสลาย time-evolution operator ด้วย Trotter decomposition

แทนที่จะ implement time-evolution operator แบบแม่นยำ เราสามารถใช้ Trotter decomposition เพื่อ implement การประมาณค่าของมันได้ การทำซ้ำ Trotter decomposition ลำดับหนึ่งหลายครั้งจะช่วยลด error ที่เกิดจากการประมาณค่าได้อีกด้วย ต่อไปนี้ เราจะสร้าง Trotter implementation โดยตรงในวิธีที่มีประสิทธิภาพสูงสุดสำหรับ interaction graph ของ Hamiltonian ที่เรากำลังพิจารณา (มีเฉพาะ nearest neighbor interactions) ในทางปฏิบัติ เราแทรก Pauli rotations $R_{xx}$, $R_{yy}$, $R_{zz}$ ด้วย parametrized angle $t$ ซึ่งสอดคล้องกับการ implement โดยประมาณของ $e^{-i (XX + YY + ZZ) t}$ เนื่องจากความแตกต่างในนิยามของ Pauli rotations และ time-evolution ที่เราพยายาม implement เราจึงต้องใช้พารามิเตอร์ $2*dt$ เพื่อให้ได้ time-evolution ของ $dt$ นอกจากนี้ เราสลับลำดับของ operations สำหรับ Trotter step ที่เป็นเลขคี่ ซึ่งให้ผลลัพธ์เทียบเท่าแต่ช่วยให้สังเคราะห์ operations ที่อยู่ติดกันเป็น unitary เดียวของ $SU(2)$ ได้ ส่งผลให้ได้ Circuit ที่ตื้นกว่ามากเมื่อเทียบกับการใช้ฟังก์ชัน 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)

ใช้ Circuit ที่ปรับปรุงแล้วสำหรับการเตรียม state

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)

Template circuits สำหรับคำนวณ matrix elements ของ $\tilde{S}$ และ $\tilde{H}$ ผ่าน Hadamard test

ความแตกต่างเพียงอย่างเดียวระหว่าง Circuit ที่ใช้ใน Hadamard test คือ phase ใน time-evolution operator และ observable ที่วัด ดังนั้นเราจึงสามารถเตรียม template circuit ซึ่งแทน Circuit ทั่วไปสำหรับ Hadamard test โดยมี placeholder สำหรับ Gate ที่ขึ้นอยู่กับ time-evolution operator

# 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

เราลด depth ของ Hadamard test ได้อย่างมาก ด้วยการรวม Trotter approximation และ uncontrolled unitaries เข้าด้วยกัน

ขั้นตอนที่ 3: รันด้วย Qiskit primitives

สร้าง instance ของ Backend และกำหนดพารามิเตอร์ runtime

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)

Transpile ไปยัง QPU

ก่อนอื่น เราเลือก subset ของ coupling map ที่มี Qubit "คุณภาพดี" (โดยที่ "ดี" ที่นี่ค่อนข้างยืดหยุ่น เราต้องการหลีกเลี่ยง Qubit ที่ทำงานได้แย่มากเป็นหลัก) และสร้าง target ใหม่สำหรับการ Transpile

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

จากนั้น Transpile virtual Circuit ไปยัง physical layout ที่ดีที่สุดใน target ใหม่นี้

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]

สร้าง PUB สำหรับรันด้วย 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)

รัน Circuit

Circuit สำหรับ $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)

รัน Circuit สำหรับ $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]

คำนวณเมทริกซ์ Hamiltonian ที่มีประสิทธิภาพและเมทริกซ์ Overlap

ขั้นแรกให้คำนวณเฟสที่สะสมโดยสถานะ $\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)
]

เมื่อได้ผลลัพธ์จากการรัน Circuit แล้ว เราสามารถประมวลผลข้อมูลเพื่อคำนวณองค์ประกอบเมทริกซ์ของ $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]$$

สุดท้าย เราสามารถแก้ปัญหา eigenvalue แบบ generalized สำหรับ $\tilde{H}$:

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

และได้ค่าประมาณของพลังงาน ground state $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

สำหรับ sector อนุภาคเดี่ยว เราสามารถคำนวณ ground state ของ sector นี้ของ 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 subspace จากการวิวัฒนาการเวลาจริง

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

สำหรับ timestep $dt$ บางค่าที่เราจะกำหนดในภายหลัง ให้สมมติชั่วคราวว่า $r$ เป็นเลขคู่: แล้วกำหนด $d=r/2$ สังเกตว่าเมื่อเรา project Hamiltonian ลงใน Krylov space ข้างต้น มันแยกแยะไม่ออกจาก Krylov space

$$\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$ timesteps สาเหตุที่แยกแยะไม่ออกก็เพราะองค์ประกอบเมทริกซ์

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

ไม่เปลี่ยนแปลงภายใต้การเลื่อนโดยรวมของเวลาวิวัฒนาการ เนื่องจากการวิวัฒนาการเวลา commute กับ Hamiltonian สำหรับ $r$ คี่ เราสามารถใช้การวิเคราะห์สำหรับ $r-1$

เราต้องการแสดงว่าที่ใดที่หนึ่งใน Krylov space นี้รับประกันว่าต้องมีสถานะพลังงานต่ำอยู่ เราทำโดยผ่านผลลัพธ์ต่อไปนี้ ซึ่งได้มาจาก Theorem 3.1 ใน [3]:

ข้อเรียกร้องที่ 1: มีฟังก์ชัน $f$ อยู่ ซึ่งสำหรับพลังงาน $E$ ในช่วง spectral ของ Hamiltonian (นั่นคือ ระหว่างพลังงาน ground state และพลังงานสูงสุด)...

1. $f(E_0)=1$
2. $|f(E)|\le2\left(1 + \delta\right)^{-d}$ สำหรับค่า $E$ ทั้งหมดที่ห่างจาก $E_0$ ไม่น้อยกว่า $\delta$ กล่าวคือ มันถูก suppress แบบ exponential
3. $f(E)$ เป็น linear combination ของ $e^{ijE\,dt}$ สำหรับ $j=-d,-d+1,...,d-1,d$

เราจะให้หลักฐานด้านล่าง แต่สามารถข้ามได้อย่างปลอดภัยหากไม่ต้องการเข้าใจการโต้แย้งเชิงเข้มงวดอย่างครบถ้วน ตอนนี้เรามุ่งเน้นไปที่ผลกระทบของข้อเรียกร้องข้างต้น จากคุณสมบัติ 3 ด้านบน เราเห็นได้ว่า Krylov space ที่เลื่อนแล้วข้างต้นบรรจุสถานะ $f(H)|\psi\rangle$ ไว้ นี่คือสถานะพลังงานต่ำของเรา เพื่อดูว่าทำไม ให้เขียน $|\psi\rangle$ ใน energy eigenbasis:

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

โดยที่ $|E_k\rangle$ คือ energy eigenstate ลำดับที่ k และ $\gamma_k$ คือ amplitude ของมันในสถานะเริ่มต้น $|\psi\rangle$ เมื่อแสดงในรูปนี้ $f(H)|\psi\rangle$ จะเป็น

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

โดยใช้ข้อเท็จจริงที่ว่าเราสามารถแทน $H$ ด้วย $E_k$ เมื่อมันกระทำบน eigenstate $|E_k\rangle$ ดังนั้น energy error ของสถานะนี้จึงเป็น

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

เพื่อแปลงสิ่งนี้ให้เป็น upper bound ที่เข้าใจได้ง่ายขึ้น เราแยกผลรวมในตัวเศษออกเป็นพจน์ที่ $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}.$$

เราสามารถ upper bound พจน์แรกด้วย $\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$ ในผลรวม และขั้นตอนที่สองตามมาเพราะผลรวมในตัวเศษเป็น subset ของผลรวมในตัวส่วน สำหรับพจน์ที่สอง ขั้นแรกเรา lower bound ตัวส่วนด้วย $|\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}$ นอกจากนี้ การ upper bound $E_k-E_0<2\|H\|$ และ upper bound $\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$ ให้เท่ากับ error เป้าหมายของเรา ค่า error bound ข้างต้นจะ converge ไปสู่ค่านั้นแบบ exponential ตาม Krylov dimension $2d=r$ นอกจากนี้ โปรดสังเกตว่าหาก $\delta<E_1-E_0$ แล้วพจน์ $\delta$ จะหายไปโดยสิ้นเชิงใน bound ข้างต้น

เพื่อให้การโต้แย้งสมบูรณ์ เราสังเกตก่อนว่าข้างต้นเป็นเพียง energy error ของสถานะเฉพาะ $f(H)|\psi\rangle$ ไม่ใช่ energy error ของสถานะพลังงานต่ำสุดใน Krylov space อย่างไรก็ตาม ตาม variational principle (Rayleigh-Ritz) energy error ของสถานะพลังงานต่ำสุดใน Krylov space ถูก upper bound โดย energy error ของสถานะใดๆ ใน Krylov space ดังนั้นข้างต้นยังเป็น upper bound ของ energy error ของสถานะพลังงานต่ำสุดด้วย นั่นคือ output ของ Krylov quantum diagonalization algorithm

การวิเคราะห์ที่คล้ายกันกับข้างต้นสามารถดำเนินการได้ซึ่งยังคำนึงถึง noise และขั้นตอน thresholding ที่กล่าวถึงใน notebook อีกด้วย ดู [2] และ [4] สำหรับการวิเคราะห์นี้

ภาคผนวก: หลักฐานของข้อเรียกร้องที่ 1

สิ่งต่อไปนี้ส่วนใหญ่ได้มาจาก [3], Theorem 3.1: ให้ $0 < a < b$ และให้ $\Pi^*_d$ เป็น space ของ residual polynomials (polynomials ที่มีค่าที่ 0 เท่ากับ 1) ที่มี degree ไม่เกิน $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).$$

เราต้องการแปลงสิ่งนี้ให้เป็นฟังก์ชันที่สามารถแสดงออกตามธรรมชาติในรูปของ complex exponentials เพราะนั่นคือการวิวัฒนาการเวลาจริงที่สร้าง quantum Krylov space ในการทำเช่นนั้น เป็นประโยชน์ที่จะแนะนำ transformation ต่อไปนี้ของพลังงานภายใน spectral range ของ Hamiltonian ไปสู่ตัวเลขในช่วง $[0,1]$: กำหนด

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

โดยที่ $dt$ คือ timestep ที่ $-\pi < E_0dt < E_\text{max}dt < \pi$ สังเกตว่า $g(E_0)=0$ และ $g(E)$ เพิ่มขึ้นเมื่อ $E$ ห่างออกไปจาก $E_0$

ตอนนี้โดยใช้ polynomial $p^*(x)$ กับ parameters 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$ คือพลังงาน ground state เราเห็นได้โดยการแทน $\cos(x)=\frac{e^{ix}+e^{-ix}}{2}$ ว่า $f(E)$ เป็น trigonometric polynomial ของ degree $d$ นั่นคือ linear combination ของ $e^{ijE\,dt}$ สำหรับ $j=-d,-d+1,...,d-1,d$ นอกจากนี้ จากนิยามของ $p^*(x)$ ข้างต้น เรามีว่า $f(E_0)=p(0)=1$ และสำหรับ $E$ ใดๆ ใน spectral range ที่ $\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).

แบบสำรวจบทเรียน

Link to survey

Note: This survey is provided by IBM Quantum and relates to the original English content. To give feedback on doQumentation's website, translations, or code execution, please open a GitHub issue.

กรุณาทำแบบสำรวจสั้นนี้เพื่อให้ feedback เกี่ยวกับบทเรียนนี้ ความคิดเห็นของคุณจะช่วยให้เราปรับปรุงเนื้อหาและประสบการณ์การใช้งาน