สูตร Multi-product เพื่อลดความผิดพลาด Trotter

ประมาณเวลาใช้งาน: สี่นาทีบนโปรเซสเซอร์ Heron r2 (หมายเหตุ: นี่เป็นเพียงการประมาณเท่านั้น เวลาจริงอาจแตกต่างกัน)

ผลลัพธ์การเรียนรู้

• multi-product formula (MPF) ช่วยลด Trotter error ในการจำลอง Hamiltonian ได้อย่างไรโดยการรวมค่าคาดหวังจาก Circuit ที่ตื้นหลายตัว

• เมื่อใด MPF จึงมีประโยชน์เหนือกว่า product formula มาตรฐาน และเมื่อใดที่ไม่ใช่เครื่องมือที่เหมาะสม

• วิธีคำนวณสัมประสิทธิ์ MPF แบบ static และ dynamic โดยใช้แพ็กเกจ qiskit_addon_mpf

• วิธีรัน workflow ของ MPF แบบครบวงจรบนฮาร์ดแวร์ IBM Quantum® รวมถึงการ transpile การลด error และการประมวลผลภายหลัง

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

• วิธี Compilation สำหรับ Circuit การจำลอง Hamiltonian — แนะนำ Circuit แบบ Trotter (product formula) ใน Qiskit

• Product formula ใน Qiskit โดยเฉพาะคลาส synthesis SuzukiTrotter และ LieTrotter

• Qiskit primitive และอินเทอร์เฟซ Estimator

พื้นหลัง

multi-product formula คืออะไร

เมื่อจำลองระบบควอนตัมบนคอมพิวเตอร์ควอนตัม งานสำคัญคือการประมาณตัวดำเนินการวิวัฒนาการตามเวลา $e^{-iHt}$ สำหรับ Hamiltonian $H$ วิธีมาตรฐานใช้ product formula (PF) หรือที่รู้จักกันในชื่อ Trotter-Suzuki decomposition วิธีนี้แยก $H = \sum_{a=1}^d F_a$ ออกเป็นพจน์ที่ unitary แต่ละตัว $e^{-iF_a t}$ นำไปใช้ได้อย่างมีประสิทธิภาพ แล้วประมาณการวิวัฒนาการทั้งหมดเป็นผลคูณที่เรียงลำดับของ unitary ที่ง่ายกว่าเหล่านี้

product formula อันดับที่หนึ่ง (Lie-Trotter) คือ:

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

ซึ่งก่อให้เกิด error แบบกำลังสอง: $S_1(t) = e^{-iHt} + \mathcal{O}(t^2)$ formula สมมาตรอันดับสูง $S_{2\chi}(t)$ ซึ่ง $\chi$ ระบุอันดับของ symmetric product formula (ดู Ref. [1]) จะลู่เข้าเร็วขึ้นเป็น $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$ แต่แลกมาด้วย Circuit ที่ลึกขึ้นต่อขั้นตอน

เพื่อลด error ที่อันดับ $\chi$ คงที่ โดยทั่วไปจะแบ่งเวลาการวิวัฒนาการทั้งหมด $t$ ออกเป็นขั้นตอน Trotter ที่เล็กลง $k$ ขั้นตอน แต่ละขั้นตอนประมาณ $e^{-iHt/k}$ ด้วย product formula และขั้นตอนต่าง ๆ ถูกนำมาต่อกัน:

$$e^{-iHt} \approx \left[S_{2\chi}(t/k)\right]^k.$$

สำหรับ symmetric formula อันดับ $2\chi$ Trotter error ที่เหลือจะปรับตัวตาม $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$ ดังนั้นการเพิ่ม $k$ จะลด Trotter error ได้อย่างรวดเร็ว — แต่ก็ทำให้ Circuit ลึกขึ้นเป็นเชิงเส้นเช่นกัน และบนฮาร์ดแวร์ที่มี noise นั่นหมายถึง gate noise ที่สะสมมากขึ้น ความตึงเครียดระหว่าง Trotter error (สนับสนุน $k$ ที่ใหญ่กว่า) และ noise ของฮาร์ดแวร์ (สนับสนุน $k$ ที่เล็กกว่า) นี้เองคือสิ่งที่ multi-product formula ถูกออกแบบมาเพื่อแก้ไข โปรดทราบว่า MPF เกี่ยวกับการรวมผลลัพธ์จาก ตัวเลือกที่ต่างกันของ $k$ ที่อันดับ $\chi$ คงที่ — มันไม่ได้เปลี่ยนอันดับของ product formula ที่เป็นพื้นฐาน

Multi-product formula (MPF) [1] สร้าง combination เชิงเส้นแบบถ่วงน้ำหนัก ของค่าคาดหวังที่ได้จาก Trotter circuit ที่ตื้นกว่าหลายตัว แต่ละตัวใช้จำนวนขั้นตอน Trotter ที่ต่างกัน $k_1, k_2, \ldots, k_r$ (ชุดของจำนวนขั้นตอน $r$ ชุด):

$$\langle A \rangle_{\text{MPF}}(t) = \sum_{j=1}^r x_j \, \langle A \rangle_{k_j}(t),$$

โดยที่ $\langle A \rangle_{k_j}(t)$ คือค่าคาดหวังของ observable $A$ ที่เวลา $t$ ซึ่งประมาณจาก Trotter circuit ที่มีขั้นตอน $k_j$ และสัมประสิทธิ์ $\{x_j\}_{j=1}^r$ ถูกเลือกเพื่อให้พจน์ Trotter error หลักในการรวมยกเลิกกัน เราจะกลับมาดูนิพจน์นี้อีกครั้งในขั้นตอนที่ 4 ซึ่งเราจะประเมินมันโดยชัดแจ้งเพื่อรวมผลลัพธ์ Trotter ของเรา ประเด็นสำคัญในทางปฏิบัติคือ Circuit ที่ลึกที่สุดใน MPF ต้องการเพียง $k_{\max}$ ขั้นตอน ซึ่งเล็กกว่า $k$ เดี่ยวที่จำเป็นในการเข้าถึง Trotter error ที่มีผลเท่ากันโดยตรงมาก Circuit ที่ตื้นกว่าทำให้วิธี MPF เหมาะกับฮาร์ดแวร์ที่มี noise มากกว่า

สัมประสิทธิ์ถูกกำหนดอย่างไร

มีสัมประสิทธิ์ MPF สองตระกูล:

สัมประสิทธิ์แบบ Static ไม่ขึ้นกับ Hamiltonian, สถานะเริ่มต้น และเวลาการวิวัฒนาการ ถูกหาโดยการแก้ระบบเชิงเส้น $Ax = b$ ที่บังคับให้พจน์ Trotter error หลักยกเลิกกัน สำหรับชุดขั้นตอน Trotter $\{k_j\}_{j=1}^r$ ที่ใช้กับ symmetric product formula อันดับ $2\chi$ การขยาย Trotter error ในกำลังผกผันของ $k_j$ นำไปสู่สมการข้อจำกัดในรูปแบบ:

$$\sum_{j=1}^r x_j = 1, \quad \sum_{j=1}^r \frac{x_j}{k_j^{\eta_n}} = 0 \quad (n = 0, \ldots, r-2),$$

โดยที่เลขชี้กำลังจำนวนเต็ม $\{\eta_n\}$ คืออันดับของพจน์ Trotter error ที่ต่อเนื่องกันสำหรับ product formula ที่เลือก สำหรับ PF อันดับ $2\chi$ ที่ สมมาตร error หลักใน $\left[S_{2\chi}(t/k)\right]^k$ ปรับตัวตาม $1/k^{2\chi}$ โดยมีการแก้ไขถัดมาที่ $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ — ดังนั้นเลขชี้กำลังคือ $\eta_n = 2\chi + 2n$ สำหรับ PF ที่ไม่สมมาตร ทั้งกำลังคี่และคู่มีส่วนร่วม และ $\eta_n = 2\chi + n$ ดู Ref. [1] สำหรับการอนุมานแบบเต็ม สมการแรกในระบบข้างต้นรับประกันความไม่ลำเอียง (MPF สร้างค่าคาดหวังที่แม่นยำในลิมิต $k_j \to \infty$) และสมการที่เหลือ $r-1$ สมการยกเลิกพจน์ Trotter error แรก $r-1$ พจน์ตามลำดับ เมื่อค่า $L_1$-norm $\|x\|_1$ ที่ได้มีค่ามากเกินไป (ซึ่งขยาย sampling noise) คุณสามารถแก้ปัญหาการปรับให้เหมาะสมโดยประมาณแทน โดยจำกัด $\|x\|_1$ ในขณะที่ลด $\|Ax - b\|$ ให้น้อยที่สุด

สัมประสิทธิ์แบบ Dynamic [2], [3] ยังขึ้นกับ Hamiltonian, สถานะเริ่มต้น และเวลาการวิวัฒนาการ $t$ ด้วย มันลดระยะห่างแบบ Frobenius-norm ระหว่างสถานะที่วิวัฒนาการตามเวลาจริงกับการประมาณ MPF:

$$\|\rho(t) - \mu^D(t)\|_F^2 = 1 + \sum_{i,j} M_{ij}(t)\, x_i(t)\, x_j(t) - 2\sum_i L_i(t)\, x_i(t),$$

โดยที่ $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ คือ Gram matrix ของการซ้อนทับระหว่างสถานะที่วิวัฒนาการแบบ Trotter สำหรับจำนวนขั้นตอนต่าง ๆ $k_i, k_j$ และ $L_i(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$ วัดการซ้อนทับกับสถานะที่แม่นยำ (โดยประมาณ) ในแบบฝึกหัดนี้ ปริมาณเหล่านี้ถูกคำนวณอย่างมีประสิทธิภาพโดยใช้วิธี tensor-network โดยเฉพาะ backend ที่อิงจาก TeNPy ใน qiskit_addon_mpf

เมื่อใดควรใช้ MPF

MPF มีประโยชน์มากที่สุดเมื่อ:

• ความลึกของ Circuit เป็นคอขวด หากความลึกที่คุณสามารถรันถูกจำกัดด้วย noise ของฮาร์ดแวร์ ให้ใช้ MPF เพื่อให้ได้ความแม่นยำ Trotter ที่มีประสิทธิผลสูงขึ้นจาก Circuit ที่ตื้นกว่า

• คุณต้องการค่าคาดหวังที่แม่นยำ ไม่ใช่การเตรียมสถานะแบบเต็มรูปแบบ MPF ทำงานที่ระดับค่าคาดหวัง — มันรวมตัวเลขแบบคลาสสิก ไม่ใช่สถานะควอนตัม ดังนั้นจึงเหมาะสำหรับการประมาณ observable เมื่อใช้ Estimator primitive

• คุณรวมจำนวนขั้นตอน Trotter ที่พอประมาณ โดยทั่วไปการรวมจำนวนขั้นตอน $k_j$ ที่ต่างกัน $r = 3$–$5$ ค่าก็เพียงพอที่จะยกเลิกพจน์ Trotter error หลักได้หลายพจน์ ในขณะที่รักษา $\|x\|_1$ ให้จัดการได้

เมื่อใดที่ MPF อาจไม่ช่วย

• เวลาการวิวัฒนาการที่สั้นมาก เมื่อ $t$ เล็กพอที่ Trotter formula อันดับต่ำเดียวก็แม่นยำอยู่แล้ว ภาระในการรัน Circuit หลายตัวก็ไม่จำเป็น

• งานเตรียมสถานะ MPF สร้างค่า ค่าคาดหวัง ที่แก้ไขแล้ว ไม่ใช่สถานะควอนตัมที่แก้ไขแล้ว หากคุณต้องการสถานะที่วิวัฒนาการตามเวลาจริง (ตัวอย่างเช่น เป็นข้อมูลนำเข้าสำหรับ subroutine ควอนตัมอื่น) MPF ใช้ไม่ได้

• จำนวนขั้นตอน Trotter ที่ละเมิด regime การลู่เข้า การอนุมานสัมประสิทธิ์ static ขยาย $\left[S_{2\chi}(t/k_j)\right]^{k_j}$ แต่ละตัวเป็นอนุกรมใน $t/k_j$; การขยายนี้ลู่เข้าได้ดีเฉพาะเมื่อ $t/k_{\min} \lesssim 1$ หาก $k_{\min}$ ถูกเลือกให้เล็กเกินไปสำหรับ $t$ ที่กำหนด Circuit ที่ตื้นที่สุดจะอยู่นอก regime การรบกวนอย่างมาก พจน์ error อันดับสูงที่ MPF ปล่อยให้ไม่ถูกยกเลิกจะมีค่ามาก และการยกเลิกอาจต้องการสัมประสิทธิ์ที่มาก $L_1$-norm $\|x\|_1$ เป็นตัววินิจฉัยในทางปฏิบัติ: เมื่อ $\|x\|_1 \gg 1$ ภาระ sampling $\propto \|x\|_1^2$ อาจมีค่ามากกว่าการลด Trotter error ดูแนวทางการเลือกขั้นตอน Trotter สำหรับรายละเอียด

สิ่งที่แบบฝึกหัดนี้ครอบคลุม

แบบฝึกหัดนี้พาไปดู workflow ของ MPF แบบครบวงจรในสองขั้นตอน อันดับแรกคือ ตัวอย่าง simulator ขนาดเล็ก (10-qubit Heisenberg chain) สาธิตวิธีตั้งค่าปัญหา คำนวณสัมประสิทธิ์ MPF แบบ static และ dynamic และเปรียบเทียบค่าคาดหวังที่ได้กับการทำ diagonalization ที่แม่นยำ จากนั้น ตัวอย่างฮาร์ดแวร์ขนาดใหญ่ (50-qubit XXZ chain) แสดงวิธี transpile รันบนฮาร์ดแวร์ IBM Quantum พร้อมการลด error และประมวลผลผลลัพธ์ภายหลังโดยใช้สัมประสิทธิ์ MPF ตลอดทั้งบทเรียน เราใช้แพ็กเกจ qiskit_addon_mpf ควบคู่ไปกับเครื่องมือมาตรฐานของ Qiskit

ข้อกำหนด

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

• Qiskit SDK v2.0 ขึ้นไป พร้อมรองรับ visualization

• Qiskit Runtime v0.22 ขึ้นไป (pip install qiskit-ibm-runtime)

• Qiskit Aer simulator (pip install qiskit-aer)

• MPF Qiskit addon พร้อม backend TeNPy (pip install "qiskit-addon-mpf[tenpy]")

• ยูทิลิตีของส่วนเสริม Qiskit (pip install qiskit-addon-utils)

• SciPy (pip install scipy)

การตั้งค่า

ด้านล่างนี้เรารวบรวมการ import แพ็กเกจ ทั้งหมด ที่ใช้ตลอดแบบฝึกหัดนี้ไว้ในเซลล์เดียว เรายังกำหนด Transpiler pass ชื่อ CollectAndCollapse ที่รวม rxx และ ryy rotation ที่อยู่ติดกันเข้าเป็น XXPlusYYGate เดียว pass นี้ถูกใช้ทั้งระหว่างการสร้าง Circuit ในขั้นตอนที่ 1 (เพื่อรักษาจำนวน gate ให้ต่ำ) และโดยอ้อมเมื่อเราแยกโครงสร้างแบบเลเยอร์ออกมาสำหรับ MPF แบบ dynamic ในขั้นตอนที่ 4 (TeNPy ต้องการ two-qubit gate ไม่ใช่คู่ของ rotation ที่ไม่ได้รวมกัน)

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

import warnings

import numpy as np
import matplotlib.pyplot as plt
from functools import partial
from copy import deepcopy

from qiskit import QuantumCircuit
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector
from qiskit.synthesis import SuzukiTrotter
from qiskit.transpiler import CouplingMap, PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)

from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService

from qiskit_addon_utils.problem_generators import (
    generate_xyz_hamiltonian,
    generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.static import setup_static_lse
from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import (
    setup_exact_problem,
    setup_sum_of_squares_problem,
    setup_frobenius_problem,
)
from qiskit_addon_mpf.backends.tenpy_layers import (
    LayerModel,
    LayerwiseEvolver,
)
from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state

from scipy.linalg import expm

# Suppress TeNPy's `unit_cell_width` future-API warning. The default
# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what
# `CouplingMap.from_line(...)` produces here, so the warning is informational.
warnings.filterwarnings(
    "ignore",
    message=r".*unit_cell_width.*",
    category=UserWarning,
)

# --- Helper: collect XX + YY rotations into a single gate ---
def filter_function(node):
    return node.op.name in {"rxx", "ryy"}

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

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

collapse_function = partial(
    collapse_to_operation,
    collapse_function=collapse_to_xx_plus_yy,
)

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

ตัวอย่าง simulator ขนาดเล็ก

ขั้นตอนที่ 1: แมปข้อมูลนำเข้าแบบคลาสสิกไปยังปัญหาควอนตัม

เราเริ่มด้วยแบบจำลอง Heisenberg แบบ 10-qubit บนเส้นตรง โดยใช้สถานะ Néel $\vert 0101\ldots01 \rangle$ เป็นสถานะเริ่มต้น Hamiltonian คือ:

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

โดยที่ $J$ คือความแรงของการมีปฏิสัมพันธ์ระหว่างเพื่อนบ้านใกล้เคียง เราวัด ZZ correlator $Z_{L/2-1} Z_{L/2}$ บนคู่ qubit ที่กึ่งกลางของ chain และใช้ขั้นตอน Trotter $k_j = [1, 2, 4]$ กับ product formula อันดับที่สอง

L = 10

# Generate coupling map and Hamiltonian
coupling_map = CouplingMap.from_line(L, bidirectional=False)

hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(1.0, 1.0, 1.0),
    ext_magnetic_field=(0.0, 0.0, 0.0),
)
print(hamiltonian)

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

# Observable: ZZ on the middle pair of qubits
observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)

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

# MPF parameters
mpf_trotter_steps = [1, 2, 4]
order = 2
symmetric = False

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

สร้าง Trotter circuit

เราสร้าง Circuit ที่นำ Trotter time-evolution แบบประมาณไปใช้สำหรับแต่ละจุดเวลาและแต่ละจำนวนขั้นตอน Trotter pass CollectAndCollapse ที่กำหนดไว้ในส่วนการตั้งค่ารวบรวม XX และ YY rotation เข้าเป็น XX+YY gate เดียว เพื่อเตรียมสำหรับการจำลอง tensor-network ที่มีประสิทธิภาพมากขึ้นในภายหลัง

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

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

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

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

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

ขั้นตอนที่ 2: ปรับแต่งปัญหาสำหรับการรันบนฮาร์ดแวร์ควอนตัม

สำหรับตัวอย่างขนาดเล็ก เรามุ่งเป้าไปที่ Aer simulator การแปลงสองอย่างเกิดขึ้นก่อนที่ Circuit จะพร้อมสำหรับการรัน:

1. การรวบรวม gate ที่ระดับการจำลอง Hamiltonian ในเซลล์การตั้งค่า เราสร้าง pass CollectAndCollapse ที่รวม rxx และ ryy rotation ที่อยู่ติดกันเข้าเป็น XXPlusYYGate เดียว เราใช้ pass นี้ไปแล้วเมื่อเราสร้าง Trotter circuit ในขั้นตอนที่ 1 (การเรียก pm.run(...)) ทำให้ทั้งลดจำนวน two-qubit gate และสร้างโครงสร้างที่เหมาะกับการจำลอง tensor-network สำหรับการคำนวณสัมประสิทธิ์แบบ dynamic ในภายหลังมากขึ้น

2. การลดระดับให้เป็น ISA ของ simulator ด้านล่างนี้เรารัน Qiskit preset pass manager ที่ optimization_level=3 เพื่อลดระดับ Trotter circuit แต่ละตัวให้เป็นสถาปัตยกรรมชุดคำสั่ง (ISA) ของ simulator

aer_sim = AerSimulator()
pm_sim = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)

isa_circs_all_times = [
    pm_sim.run([deepcopy(c) for c in mpf_circuits])
    for mpf_circuits in all_circs
]

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

สำหรับตัวอย่างขนาดเล็ก เรารัน Trotter circuit ที่ลดระดับเป็น ISA แล้วผ่าน primitive EstimatorV2 ที่ใช้ Aer เป็นฐาน การทำเช่นนี้ให้ค่าอ้างอิงที่ ไม่มี noise สำหรับแต่ละคู่ $(k_j, t)$ — ค่าเหล่านี้คือค่า $\langle A \rangle_{k_j}(t)$ ที่ MPF จะรวมในขั้นตอนที่ 4 เรากวาดผ่านเวลาการวิวัฒนาการเพื่อให้ภายหลังเราสามารถพล็อตกราฟอนุกรมเวลาแบบเต็มของ product formula แต่ละตัวและของ MPF

estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for isa_circuits in isa_circs_all_times:
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals_all_times.append([res.data.evs for res in result])
    mpf_stds_all_times.append([res.data.stds for res in result])

ขั้นตอนที่ 4: Post-process และแสดงผลลัพธ์ในรูปแบบคลาสสิกที่ต้องการ

ขั้นตอนที่ 4 คือจุดที่ MPF ถูกสร้างขึ้นจริง แม้ว่าสัมประสิทธิ์ $x_j$ จะถูก คำนวณ ที่นี่ (และสำหรับตัวแปร dynamic การคำนวณนี้อาจใช้ทรัพยากรมาก) แต่ในเชิงแนวคิด สัมประสิทธิ์เหล่านี้เป็นสูตรคลาสสิกสำหรับการรวมการวัดควอนตัมจากขั้นตอนที่ 3 เข้าเป็นค่าคาดหวังที่แก้ไขแล้วค่าเดียว — ดังนั้นเราจึงถือว่า workflow ของสัมประสิทธิ์และการรวมกันทั้งหมดเป็นการประมวลผลภายหลัง

เพื่อประเมินว่า MPF ติดตามพลศาสตร์จริงได้ดีเพียงใด ก่อนอื่นเราคำนวณค่าคาดหวังที่วิวัฒนาการตามเวลาจริงโดยการทำ exponentiate ของ Hamiltonian โดยตรง วิธีนี้ทำได้ในทางปฏิบัติเพราะ $L = 10$; ในตัวอย่างฮาร์ดแวร์ขนาดใหญ่ด้านล่าง เราจะต้องพึ่งพาการประมาณด้วย tensor-network แทน

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

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

สัมประสิทธิ์ MPF แบบคงที่

MPF แบบ static ใช้สัมประสิทธิ์ $x_j$ ที่ไม่ขึ้นกับเวลาการวิวัฒนาการ, Hamiltonian และสถานะเริ่มต้น เราตั้งค่าระบบเชิงเส้น $Ax = b$ ที่อธิบายไว้ในพื้นหลังและแก้หาสัมประสิทธิ์ matrix $A$ ถูกกำหนดโดยจำนวนขั้นตอน Trotter $k_j$, อันดับ $\chi$ ของ product formula และว่า formula นั้นสมมาตรหรือไม่ (ซึ่งควบคุมเลขชี้กำลัง $\eta_n$)

สำหรับตัวอย่างขนาดเล็กของเรา เราใช้ $k_j = [1, 2, 4]$ กับ Suzuki-Trotter formula อันดับ $2\chi=2$ ที่ไม่สมมาตร (ดังนั้น $\chi=1$ และ $\eta_n = 2 + n$ ทำให้ $\eta_0 = 2,\, \eta_1 = 3$) ระบบกลายเป็น:

$$A = \begin{bmatrix} 1 & 1 & 1\\ 1 & \frac{1}{2^2} & \frac{1}{4^2} \\ 1 & \frac{1}{2^3} & \frac{1}{4^3} \\ \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}.$$

แถวแรกบังคับความไม่ลำเอียง ($\sum_j x_j = 1$); แถวที่สองและสามยกเลิกพจน์ Trotter error หลัก $1/k^2$ และอันดับถัดไป $1/k^3$ ตามลำดับ

ตั้งค่า LSE

เราใช้ setup_static_lse จาก qiskit_addon_mpf.static เพื่อประกอบ matrix $A$ และเวกเตอร์ด้านขวา $b$ ที่อธิบายไว้ข้างต้น matrix $A$ ไม่ได้ขึ้นกับ $k_j$ เท่านั้น แต่ยังขึ้นกับตัวเลือก product formula ของเราด้วย — โดยเฉพาะ อันดับ $\chi$ และว่ามัน สมมาตร หรือไม่ แฟล็ก symmetric ควบคุมรูปแบบเลขชี้กำลัง $\eta_n$ (formula ที่สมมาตรสร้างพจน์ Trotter error แบบกำลังคู่เท่านั้น; ดู Ref. [1]) โปรดทราบว่า ดังที่แสดงใน Ref. [2] การตั้งค่า symmetric=True ไม่จำเป็นอย่างเคร่งครัดแม้ว่า PF พื้นฐานจะสมมาตร — LSE ที่ไม่สมมาตรยังคงใช้ได้ (มันบังคับข้อจำกัดเพิ่มเติมที่ไม่จำเป็น)

สำหรับตัวอย่างของเรา เราได้ตั้งค่า order = 2 และ symmetric = False ไว้แล้วในขั้นตอนที่ 1

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

ตรวจสอบ matrix $A$ และเวกเตอร์ $b$ ที่สร้างขึ้นเพื่อยืนยันว่าตรงกับระบบที่เขียนไว้ข้างต้น

lse.A

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

lse.b

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

เมื่อมี LSE พร้อมแล้ว เราแก้หาสัมประสิทธิ์ static $x_j$ ผ่าน lse.solve() (นี่คือคำตอบโดยตรง $x = A^{-1}b$)

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

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

หาค่า $x$ ที่เหมาะสมโดยใช้ exact model

นอกจากการคำนวณ $x=A^{-1}b$ แล้ว ยังสามารถใช้ setup_exact_model เพื่อสร้าง instance ของ cvxpy.Problem ที่ใช้ LSE เป็น constraints และ optimal solution ของมันจะให้ค่า $x$

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

[ 0.04761905 -0.57142857  1.52380952]

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

L1 norm of the exact coefficients: 2.1428571428556378

หาค่า $x$ ที่เหมาะสมโดยใช้ approximate model

อาจเกิดขึ้นได้ที่ $L_1$ norm สำหรับชุดค่า $k_j$ ที่เลือกถูกพิจารณาว่าสูงเกินไป หากเป็นเช่นนั้นและไม่สามารถเลือกชุดค่า $k_j$ ที่ต่างออกไปได้ คุณสามารถใช้คำตอบโดยประมาณที่จำกัด $L_1$-norm ให้อยู่ที่ threshold ที่เลือกไว้ ในขณะที่ลด $\|Ax - b\|$ ให้น้อยที่สุด ดูแนวทางวิธีใช้โมเดลโดยประมาณ

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

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

สัมประสิทธิ์ MPF แบบไดนามิก

MPF แบบ static ยกเลิกพจน์ Trotter error ในลักษณะที่ไม่ขึ้นกับ Hamiltonian และสถานะ ดังนั้นจึงอาจไม่ให้ error การประมาณที่เล็กที่สุดเท่าที่เป็นไปได้เสมอไปสำหรับ Hamiltonian และสถานะเริ่มต้นที่กำหนด MPF แบบ dynamic (Ref. [2], [3]) หาสัมประสิทธิ์ที่ขึ้นกับเวลา $x_i(t)$ ที่ลดระยะห่างแบบ Frobenius-norm $\|\rho(t) - \mu^D(t)\|_F^2$ ที่แต่ละเวลา $t$ แทน ดังที่แสดงในพื้นหลัง สิ่งนี้ต้องการ matrix ของการซ้อนทับ $M_{ij}(t)$ ระหว่างสถานะที่วิวัฒนาการแบบ Trotter และการซ้อนทับ $L_i(t)$ กับสถานะที่แม่นยำ — ทั้งสองอย่างเราประมาณโดยใช้ backend tensor-network (TeNPy) ใน qiskit_addon_mpf

ในการตั้งค่า LSE แบบ dynamic เราต้องการส่วนประกอบสามอย่าง:

1. Factory ของ evolver โดยประมาณ ที่ addon จะรันสำหรับแต่ละ $k_j$ เพื่อสร้าง $\rho_{k_j}(t)$ เป็น MPS/MPO เราสร้างมันจากโครงสร้างแบบเลเยอร์ของ Trotter circuit อันดับ 2 (หนึ่งเลเยอร์ต่อ slice_by_depth) ห่อหุ้มเป็น LayerwiseEvolver พร้อมพารามิเตอร์การตัด TeNPy

2. Factory ของ evolver ที่แม่นยำ ที่สร้างค่าอ้างอิงความแม่นยำสูง $\rho(t)$ เราใช้ Suzuki-Trotter circuit อันดับที่สี่ที่มีขั้นตอนเวลาเล็ก (dt=0.1, order=4) เป็นตัวแทนของการวิวัฒนาการที่แม่นยำ

3. Factory ของ identity และ MPS ของสถานะเริ่มต้น ที่ใช้เริ่มการจำลอง TeNPy

เซลล์ด้านล่างสร้าง factory ของ evolver โดยประมาณ

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

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

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

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

คำเตือน

ต้องเลือก options ของ LayerwiseEvolver ที่กำหนดรายละเอียดของการจำลอง tensor network อย่างระมัดระวัง เพื่อหลีกเลี่ยงการตั้งค่าปัญหา optimization ที่ไม่ชัดเจน

เราประมาณสถานะที่วิวัฒนาการตามเวลาอย่างแม่นยำด้วยสูตร fourth-order Suzuki-Trotter โดยใช้ time step dt=0.1 ที่มีค่าน้อย พารามิเตอร์ truncation ของ TeNPy สามารถส่งผลต่อความแม่นยำได้ ดังนั้นจึงเป็นสิ่งสำคัญที่จะสำรวจช่วงของค่าต่าง ๆ

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

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

สุดท้าย เรากำหนด identity_factory ที่ให้สถานะ MPO เริ่มต้น และเตรียมสถานะเริ่มต้นแบบ Néel เป็น MPS ที่ตรงกับ lattice ที่ใช้ในแบบจำลอง Trotter แบบเลเยอร์

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

mps_initial_state = MPS_neel_state(models[0].lat)

เมื่อ factory พร้อมแล้ว ตอนนี้เราคำนวณสัมประสิทธิ์ dynamic ที่แต่ละเวลาการวิวัฒนาการ สำหรับแต่ละ $t$ setup_dynamic_lse สร้าง matrix ของการซ้อนทับที่เกี่ยวข้องผ่าน TeNPy และ setup_frobenius_problem ส่งคืน cvxpy.Problem ที่ลดต้นทุน Frobenius-norm ตัวแก้ปัญหาส่งคืนสัมประสิทธิ์ $x_j(t)$ ที่ปรับให้เข้ากับเวลานั้น เราเก็บสัมประสิทธิ์เหล่านี้ไว้ใน mpf_dynamic_coeffs_list หากตัวแก้ปัญหาล้มเหลวสำหรับ $t$ ที่กำหนด เราจะกลับไปใช้สัมประสิทธิ์ศูนย์เพื่อให้ loop ดำเนินต่อไป

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

Computing dynamic coefficients for time=0.5

Computing dynamic coefficients for time=0.6

Computing dynamic coefficients for time=0.7

Computing dynamic coefficients for time=0.7999999999999999

Computing dynamic coefficients for time=0.8999999999999999

Computing dynamic coefficients for time=0.9999999999999999

Computing dynamic coefficients for time=1.0999999999999999

Computing dynamic coefficients for time=1.1999999999999997

Computing dynamic coefficients for time=1.2999999999999998

Computing dynamic coefficients for time=1.4

Computing dynamic coefficients for time=1.4999999999999998

รวมค่าคาดหวัง Trotter เข้ากับสัมประสิทธิ์ MPF

ตอนนี้เราประเมิน $\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)$ สำหรับชุดสัมประสิทธิ์แต่ละชุด (static-exact, static-approximate และ dynamic) แพร่ความคลาดเคลื่อนมาตรฐานต่อ Circuit และพล็อตอนุกรมเวลาที่ได้เทียบกับกราฟ exact-diagonalization

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

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

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

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

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

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

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

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

# Exact expectation values
plt.plot(
    exact_evolution_times,
    exact_expvals,
    color="red",
    linestyle="--",
    label="Exact time-evolution",
)

plt.title(f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ vs time")
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

กราฟด้านบนแสดงให้เห็นปฏิสัมพันธ์ระหว่าง Trotter error และ sampling error

• Trotter error product formula แต่ละตัว (จุดสีเทา) เบี่ยงเบนจากกราฟที่แม่นยำมากขึ้นเรื่อย ๆ เมื่อเวลาเพิ่มขึ้น Circuit ที่ $k=1$ มีการเบี่ยงเบนมากที่สุดและตื้นที่สุด แต่มันก็อยู่ใน regime ที่ $t/k \gtrsim 1$ แล้ว ดังนั้นพจน์ error หลัก $1/k^{2}$ จึงมีค่ามาก การรวม MPF (จุดสี) ยกเลิกพจน์ Trotter error หลักเหล่านี้หลายพจน์ ดังนั้นจึงติดตามกราฟที่แม่นยำได้ใกล้กว่า Circuit $k_j$ เดี่ยวใด ๆ มาก ช่องว่างที่เหลืออยู่สะท้อนถึงพจน์ Trotter อันดับสูงที่ MPF ไม่ ยกเลิก: MPF แบบ static อันดับ $2$, $r=3$ ยกเลิกเพียงอันดับ error สองอันดับแรกเท่านั้น และที่ $t/k_{\min}$ มาก หางที่ไม่ถูกยกเลิกในที่สุดก็มีอิทธิพลเหนือกว่า — ดังนั้น MPF จึงไม่รับประกันว่า Circuit ที่ตื้นมากจะยังคงแม่นยำที่เวลาใดก็ได้

• Sampling error ค่าความคลาดเคลื่อนที่กว้างกว่าบนกราฟ MPF เป็นผลโดยตรงจากการรวมเชิงเส้น: การแพร่ความคลาดเคลื่อนมาตรฐานต่อ Circuit ที่เป็นอิสระ $\sigma_{k_j}$ ให้ความแปรปรวนรวม $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$ ดังนั้น ยิ่ง $\|x\|_2$ ใหญ่ขึ้น (และในทางปฏิบัติคือ $\|x\|_1$ ซึ่งเป็นสิ่งที่เราควบคุม) ยิ่งต้องการ shot มากขึ้นเพื่อให้ได้ความไม่แน่นอนเป้าหมายที่กำหนด นี่คือข้อแลกเปลี่ยนที่อยู่เบื้องหลังตัวเลือก approximate-solver ในพื้นหลัง: เราจำกัด $\|x\|_1$ เพื่อรักษาภาระนี้ให้จัดการได้ ที่สำคัญคือ ต่างจาก Trotter error, sampling error หดตัวลงตาม $1/\sqrt{N_{\text{shots}}}$ ดังนั้นจึงสามารถลดลงได้เสมอโดยการใช้ shot มากขึ้น

ในตัวอย่างฮาร์ดแวร์ขนาดใหญ่ด้านล่าง noise ของฮาร์ดแวร์เข้ามาเป็นแหล่ง error เพิ่มเติมในแต่ละ $\langle A \rangle_{k_j}$ ซึ่งถูกขยายโดยสัมประสิทธิ์ MPF ในลักษณะเดียวกัน เราจะดูว่าการลด error มีปฏิสัมพันธ์กับ MPF อย่างไรในส่วนนั้น

ตัวอย่างฮาร์ดแวร์ขนาดใหญ่

ในส่วนนี้เราขยายปัญหาให้ใหญ่ขึ้นเกินกว่าที่จะจำลองได้อย่างแม่นยำ เราสร้างผลลัพธ์บางส่วนที่แสดงใน Ref. [3] ขึ้นมาใหม่ โดยใช้ 50-qubit XXZ chain ที่เวลา $t = 3$ เราทำตาม workflow สี่ขั้นตอนเดียวกับตัวอย่างขนาดเล็ก โดยตอนนี้มุ่งเป้าไปที่ฮาร์ดแวร์ควอนตัมจริงพร้อมการลด error เช่นเดียวกับในเทมเพลต แต่ละขั้นตอนถูกทำเครื่องหมายไว้แบบ inline ในโค้ด และขั้นตอนเดียวสามารถครอบคลุมหลายเซลล์ได้เมื่อผลลัพธ์ระหว่างทางของมันควรค่าแก่การตรวจสอบ การแมปสะท้อนตัวอย่างขนาดเล็ก: กำหนด Hamiltonian, เลือกพารามิเตอร์ Trotter, คำนวณสัมประสิทธิ์ MPF (static และ dynamic) และสร้าง Circuit ความแตกต่างหลักคือ:

• Hamiltonian แบบ XXZ บน 50 ตำแหน่งพร้อม coupling แบบสุ่มที่ดึงมาจาก $\mathcal{U}(0.5, 1.5)$ (Ref. [3])

• Trotter formula อันดับสองแบบ สมมาตร ด้วย $k_j = [3, 4, 6]$ (ดังนั้น $\chi=1$, symmetric=True)

• เวลาการวิวัฒนาการคงที่เดียว $t = 3$ ด้วย $k_{\min}=3$ นี้ให้ $t/k_{\min}=1$ ซึ่งรักษาส่วนประกอบที่ตื้นให้อยู่ใน regime การลู่เข้าของ Trotter ที่โมเดล error หลักที่ MPF พึ่งพานั้นใช้ได้

• การรันเปรียบเทียบแบบ Circuit เดียวเพิ่มเติมด้วยขั้นตอน Trotter $k = 10$ ใช้เป็น baseline เราเลือก $k = 10$ เพราะความลึกแบบ two-qubit บนฮาร์ดแวร์ของมันลึกกว่าส่วนประกอบ MPF ที่ลึกที่สุด ($k_{\max}=6$) บวกกับภาระในการรัน Circuit MPF หลายตัว — ลึกพอที่จะถูกจำกัดด้วย noise ซึ่งเป็น regime ที่การรวม MPF ถูกคาดหวังให้มีประสิทธิภาพเหนือกว่า baseline แบบ Circuit เดียว มันคือการเปรียบเทียบ "Circuit เดียวที่ลึก" กับการรวม MPF ไม่ใช่ Circuit ที่มุ่งเป้าไปที่ Trotter error ที่มีประสิทธิผลของ MPF (ซึ่งจะต้องการขั้นตอนมากกว่านี้อีกมาก)

โปรดทราบว่าแม้เรายังอยู่ในขั้นตอนที่ 1 ที่นี่ (การแมปและการสร้าง Circuit) เรายังคำนวณสัมประสิทธิ์ dynamic ล่วงหน้าควบคู่ไปกับสัมประสิทธิ์ static ในเซลล์นี้ด้วย สัมประสิทธิ์ dynamic ขึ้นกับ $H$ และ $t$ แต่ไม่ขึ้นกับการวัดควอนตัม ดังนั้นจึงสามารถคำนวณได้ทุกเมื่อก่อนขั้นตอนที่ 4 เราทำมันตอนนี้เพื่อให้การตั้งค่าเฉพาะของ MPF ทั้งหมดอยู่ในที่เดียว

# -------------------------Step 1-------------------------
L = 50
coupling_map = CouplingMap.from_line(L, bidirectional=False)

# XXZ Hamiltonian with random couplings (Ref. [3])
np.random.seed(0)
even_edges = list(coupling_map.get_edges())[::2]
odd_edges = list(coupling_map.get_edges())[1::2]

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

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

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

# Static coefficients
lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(f"Static coefficients: {mpf_coeffs}")
print(f"L1 norm: {np.linalg.norm(mpf_coeffs, ord=1)}")

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=2.0
)
model_approx.solve()
print(f"Approximate coefficients: {coeffs_approx.value}")
print(f"L1 norm (approx): {np.linalg.norm(coeffs_approx.value, ord=1)}")

# -------------------------Dynamic coefficients-------------------------
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)

layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

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

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

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

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

mps_initial_state = MPS_neel_state(models[0].lat)

print(f"Computing dynamic coefficients for time={total_time}")
lse_dyn = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs_dyn = setup_frobenius_problem(lse_dyn)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs_dyn.value
except Exception as error:
    mpf_dynamic_coeffs = np.zeros(len(mpf_trotter_steps))
    print(error, "Calculation Failed")

# -------------------------Step 1 (cont): Build circuits-------------------------
mpf_circuits = []
for k in mpf_trotter_steps:
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2])
    trotter_circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(reps=k, order=order),
        time=total_time,
    )
    circuit.compose(trotter_circ, qubits=range(L), inplace=True)
    mpf_circuits.append(circuit)

# Baseline "single deep circuit" comparison run with k=10 Trotter steps.
# Its two-qubit depth is deeper than the deepest MPF constituent (k_max=6) plus
# the overhead of running multiple circuits, pushing it into the noise-limited
# regime where MPF is expected to outperform. It does NOT target the MPF's effective
# Trotter error (which would require many more steps).
comp_circuit = QuantumCircuit(L)
comp_circuit.x([i for i in range(L) if i % 2])
trotter_circ = generate_time_evolution_circuit(
    hamiltonian,
    synthesis=SuzukiTrotter(reps=10, order=order),
    time=total_time,
)
comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)
mpf_circuits.append(comp_circuit)

Static coefficients: [ 0.42857143 -1.82857143  2.4       ]
L1 norm: 4.65714285714286
Approximate coefficients: [-0.4942491   0.40206845  1.09218065]
L1 norm (approx): 1.9884981979026675
Computing dynamic coefficients for time=3

ตอนนี้เราปรับ Circuit ให้เหมาะกับ Backend ที่เลือกไว้ เราใช้ Qiskit preset pass manager ที่ optimization_level=3 ซึ่งเลือกชุด physical qubit ที่ดีโดยอัตโนมัติและกำหนดเส้นทาง Circuit แต่ละตัวไปยัง topology ของอุปกรณ์

# -------------------------Step 2-------------------------
service = QiskitRuntimeService()
# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=L)
backend = service.backend("ibm_fez")
print(backend)

transpiler = generate_preset_pass_manager(
    optimization_level=3, backend=backend
)
transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

<IBMBackend('ibm_fez')>

การรัน Circuit ที่ลึกกว่าบนฮาร์ดแวร์จริงต้องการการลด error ที่เข้มข้น เราเปิดใช้ dynamical decoupling, gate และ measurement twirling, การลด error ในการวัด และ zero-noise extrapolation (ZNE) โปรดทราบว่าตัวประกอบ noise ของ ZNE ที่เราใช้ที่นี่ (1, 1.2, 1.4) เล็กกว่าในสถานการณ์ Circuit ตื้น เนื่องจากส่วนประกอบ MPF ที่ลึกกว่านั้นใกล้เคียงกับ threshold ของ noise อยู่แล้ว และการขยาย noise ขนาดใหญ่จะผลักดันให้เกินจุดที่การประมาณ ZNE เชื่อถือได้

เราส่ง Circuit ทั้งสี่ตัว (ส่วนประกอบ MPF สามตัวที่ $k_j = [3, 4, 6]$ บวก baseline $k = 10$) ใน Estimator job เดียว

# -------------------------Step 3-------------------------
estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

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

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

estimator.options.environment.job_tags = ["TUT_MPF"]

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

เราดึงค่าคาดหวังและค่าเบี่ยงเบนมาตรฐานต่อ Circuit จากผลลัพธ์ job แล้วรวมเข้ากับสัมประสิทธิ์ MPF แต่ละชุดเช่นเดียวกับในตัวอย่างขนาดเล็ก: $\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}$ พร้อมความแปรปรวนที่แพร่แล้ว $\sigma^2 = \sum_j x_j^2 \sigma_{k_j}^2$

# -------------------------Step 4-------------------------
result = job_50.result()
evs = [res.data.evs for res in result]
std = [res.data.stds for res in result]

print(evs)
print(std)

[array(-0.07916195), array(-0.04479681), array(-0.2560756), array(-0.06045848)]
[array(0.04605538), array(0.10056336), array(0.14426151), array(0.04059092)]

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

Exact static MPF expectation value:  -0.5665938395816946 +- 0.3925273058119915
Approximate static MPF expectation value:  -0.25856647611537903 +- 0.164249927266166
Dynamic MPF expectation value:  -0.12667812062949296 +- 0.06059471006973169

sym = {3: "^", 4: "s", 6: "p"}
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs[k],
        yerr=std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

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

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

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

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

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

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

ข้อสังเกตบางประการเกี่ยวกับผลลัพธ์ฮาร์ดแวร์ข้างต้น:

• การลึกขึ้นไม่ใช่ของฟรีบนฮาร์ดแวร์ baseline แบบ Circuit เดียวบอกเล่าเรื่องราวโดยตรง: Circuit ที่ $k = 6$ แม่นยำเกือบสมบูรณ์ ($-0.256$ เทียบกับค่าอ้างอิง $-0.244$) แต่ baseline ที่ $k = 10$ ที่ลึกกว่ากลับ แย่กว่า ($-0.061$ คลาดเคลื่อนไปประมาณ $0.18$) ไม่ได้ดีกว่า เมื่อ Trotter error มีค่าน้อยอยู่แล้ว การเพิ่มขั้นตอนส่วนใหญ่จะทำให้ Circuit ลึกขึ้นและสะสม gate noise กับ decoherence มากขึ้น นี่คือ regime ที่ MPF ถูกสร้างขึ้นมาโดยเฉพาะ: เข้าถึงความแม่นยำของ Circuit ที่ลึกโดยใช้เพียงส่วนประกอบที่ตื้น

• MPF ที่มี norm เล็กเอาชนะ Circuit เดี่ยวที่ลึก MPF แบบ approximate-static (จำกัดที่ $\|x\|_1 \approx 2$) ได้ค่า $-0.259$ ซึ่งอยู่ในระยะประมาณ $0.015$ จากค่าอ้างอิง และใกล้เคียงกว่า baseline $k = 10$ มาก MPF แบบ dynamic ($-0.127$) ก็เอาชนะ baseline นั้นได้อย่างสบาย ๆ เช่นกัน ทั้งสองรวมเฉพาะ Circuit ตื้น $k_j = [3, 4, 6]$ เท่านั้น แต่กู้คืนคำตอบที่ Circuit เดี่ยวที่ลึกทำไม่ได้

• Norm ของสัมประสิทธิ์สำคัญกว่า mathematical optimality exact-static MPF มี $\|x\|_1 = 4.66$ และเป็นตัวประมาณค่าที่ แย่ที่สุด ในบรรดาทั้งหมด ($-0.567$ ผิดพลาดไปมากกว่า $0.3$): norm ของสัมประสิทธิ์ที่มีขนาดใหญ่จะขยาย residual gate noise, decoherence, และ ZNE error บน $\langle A \rangle_{k_j}$ แต่ละตัวด้วยปัจจัยที่ใกล้เคียงกัน จนกลบการหักล้างของ Trotter-error ที่มันควรได้มา การจำกัดขนาด norm (approximate-static solver, $\|x\|_1 \approx 2$) ช่วยขจัดผลกระทบนี้และให้ค่าประมาณที่ดีที่สุด — แม้ว่า สัมประสิทธิ์ ของมันจะไม่หักล้าง leading Trotter error ได้อย่างแม่นยำอีกต่อไปก็ตาม

• Circuit ตื้นเดี่ยว ๆ ก็ยังสามารถแข่งขันได้ ส่วนประกอบเดี่ยวที่ $k = 6$ ($-0.256$) นั้นแม่นยำเกือบสมบูรณ์ในกรณีนี้ — ในการรันนี้มันยังใกล้เคียงกว่า approximate-static MPF เล็กน้อยด้วยซ้ำ ปัญหาคือคุณไม่รู้ล่วงหน้าว่า $k$ ตัวไหนที่อยู่ในจุดที่เหมาะสมของ "converged แต่ยังไม่ถูกจำกัดด้วย noise" และทางเลือกที่ดูปลอดภัยอย่างการเพิ่มความลึก ($k = 10$) เพื่อรับประกัน Trotter convergence นั้นกลับเป็นทางเลือกที่ล้มเหลวพอดี MPF ให้การผสมผสาน Circuit ตื้นอย่างมีหลักการโดยไม่ต้องเดาความลึกที่ถูกต้อง

ข้อสรุปเชิงปฏิบัติคือ บนฮาร์ดแวร์จริง MPF ควรใช้คู่กับการลด error ที่แข็งแกร่งบน $\langle A \rangle_{k_j}$ แต่ละตัว, coefficient $L_1$-norm ควรถูกควบคุมให้อยู่ในระดับปานกลาง (ใช้ approximate solver หรือ dynamic MPF) และ Trotter step $k_j$ ควรถูกเลือกให้ $t/k_{\min} \lesssim 1$ — ในที่นี้ $k_{\min} = 3$ ที่ $t = 3$ ให้ $t/k_{\min} = 1$ ซึ่งทำให้องค์ประกอบต่าง ๆ อยู่ในช่วงที่ converge ซึ่งเป็นช่วงที่ leading-error model ที่ static MPF พึ่งพาใช้ได้จริง ด้วยตัวเลือกเหล่านี้ MPF ที่มี norm เล็กในที่นี้ให้ผลตรงกับวงจรเดี่ยวที่ converged ในขณะที่แนวทาง "เพิ่มความลึกเรื่อย ๆ" แบบไร้เดียงสาทำไม่ได้ ซึ่งกู้คืนข้อได้เปรียบด้าน depth-versus-accuracy ที่แสดงไว้ใน Ref. [3] นอกจากนี้โปรดสังเกตว่าการรันแต่ละครั้งมี noise — ในการส่งงานเดียวกันครั้งอื่น (หรือ backend อื่น) ลำดับที่แน่นอนอาจเปลี่ยนแปลงได้ แนวโน้มที่มั่นคงคือ MPF ที่มี $\|x\|_1$ เล็กให้ผลดี, exact-static MPF ที่มี $\|x\|_1$ ใหญ่ถูกขยายด้วย noise ของฮาร์ดแวร์ และวงจรเดี่ยวที่ลึกเกินไปถูกจำกัดด้วย noise

ขั้นตอนถัดไป

คำแนะนำ

หากคุณพบว่างานนี้น่าสนใจ คุณอาจสนใจเนื้อหาต่อไปนี้:

• วิธีเลือก Trotter step สำหรับ MPF — คำแนะนำเชิงปฏิบัติเกี่ยวกับการเลือกค่า $k_j$ เพื่อหลีกเลี่ยงความไม่เสถียร
• วิธีใช้ approximate model — การปรับ constraint ของ $L_1$-norm และตัวเลือกของ solver สำหรับ approximate static MPF
• เอกสารอ้างอิง API ของ qiskit-addon-mpf — เอกสารฉบับเต็มสำหรับโมดูล static, dynamic และ backend

อ้างอิง

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

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

[3] Robertson, N. F., et al. (2024). Tensor network enhanced dynamic multiproduct formulas. arXiv preprint arXiv:2407.17405.