주 콘텐츠로 건너뛰기

Trotter 오차를 줄이기 위한 다중곱 공식

사용 시간 추정치: Heron r2 프로세서에서 4분 (참고: 이는 추정치일 뿐입니다. 실제 실행 시간은 다를 수 있습니다.)

학습 목표​

  • 다중곱 공식(MPF)이 여러 얕은 회로의 기댓값을 결합하여 해밀토니안 시뮬레이션의 Trotter 오차를 줄이는 방법

  • MPF가 표준 곱 공식보다 유리한 경우와 적절한 도구가 아닌 경우

  • qiskit_addon_mpf 패키지를 사용하여 정적 및 동적 MPF 계수를 계산하는 방법

  • 트랜스파일, 오류 완화, 후처리를 포함하여 IBM Quantum® 하드웨어에서 MPF 워크플로를 처음부터 끝까지 실행하는 방법

사전 준비 사항​

배경​

다중곱 공식이란 무엇인가?​

양자 컴퓨터에서 양자계를 시뮬레이션할 때 핵심 과제는 해밀토니안 HH에 대한 시간 전개 연산자 e−iHte^{-iHt}를 근사하는 것입니다. 표준 접근법은 Trotter-Suzuki 분해로도 알려진 곱 공식(PF)을 사용합니다. 이는 H=∑a=1dFaH = \sum_{a=1}^d F_a를 개별 유니터리 e−iFate^{-iF_a t}를 효율적으로 구현할 수 있는 항들로 분해한 다음, 이러한 더 단순한 유니터리들의 순서 있는 곱으로 전체 전개를 근사합니다.

1차 곱 공식(Lie-Trotter)은 다음과 같습니다:

S1(t):=∏a=1de−iFat,S_1(t) := \prod_{a=1}^d e^{-i F_a t},

이는 이차 오차를 발생시킵니다: S1(t)=e−iHt+O(t2)S_1(t) = e^{-iHt} + \mathcal{O}(t^2). 대칭 곱 공식의 차수를 나타내는 χ\chi를 갖는 고차 대칭 공식 S2χ(t)S_{2\chi}(t)(참고 문헌 [1] 참조)는 e−iHt+O(t2χ+1)e^{-iHt} + \mathcal{O}(t^{2\chi+1})로 더 빠르게 수렴하지만, 스텝당 더 깊은 회로가 필요합니다.

고정된 차수 χ\chi에서 오차를 줄이기 위해 일반적으로 전체 전개 시간 tt를 kk개의 더 작은 Trotter 스텝으로 분할합니다. 각 스텝은 곱 공식으로 e−iHt/ke^{-iHt/k}를 근사하며 스텝들이 연결됩니다:

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

2χ2\chi차 대칭 공식의 경우 잔여 Trotter 오차는 O ⁣(t2χ+1/k2χ)\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)로 스케일링됩니다. 따라서 kk를 증가시키면 Trotter 오차가 빠르게 억제됩니다 — 하지만 회로가 선형적으로 깊어지기도 하며, 잡음이 있는 하드웨어에서는 이것이 더 많은 누적 게이트 잡음을 의미합니다. **Trotter 오차(더 큰 kk를 선호)**와 하드웨어 잡음(더 작은 kk를 선호) 사이의 이러한 긴장 관계가 바로 다중곱 공식이 해결하도록 설계된 문제입니다. MPF는 고정된 차수 χ\chi에서 서로 다른 kk 선택의 결과를 결합하는 것이며, 기본 곱 공식의 차수를 바꾸지 않는다는 점에 유의하세요.

다중곱 공식(MPF) [1]은 서로 다른 수의 Trotter 스텝 k1,k2,…,krk_1, k_2, \ldots, k_r (rr개의 스텝 수 집합)을 각각 사용하는 여러 개의 더 얕은 Trotter 회로에서 얻은 기댓값들의 가중 선형 결합을 구성합니다:

⟨A⟩MPF(t)=∑j=1rxj ⟨A⟩kj(t),\langle A \rangle_{\text{MPF}}(t) = \sum_{j=1}^r x_j \, \langle A \rangle_{k_j}(t),

여기서 ⟨A⟩kj(t)\langle A \rangle_{k_j}(t)는 kjk_j개의 스텝을 갖는 Trotter 회로로부터 추정된 시간 tt에서의 관측량 AA의 기댓값이며, 계수 {xj}j=1r\{x_j\}_{j=1}^r는 결합에서 선행 Trotter 오차 항이 상쇄되도록 선택됩니다. 이 식은 4단계에서 다시 다루며, 그곳에서 Trotter 결과를 결합하기 위해 명시적으로 평가합니다. 핵심적인 실용적 포인트는 MPF에서 가장 깊은 회로가 kmax⁡k_{\max} 스텝만 필요로 한다는 것이며, 이는 동일한 유효 Trotter 오차에 직접 도달하기 위해 필요한 단일 kk보다 훨씬 작습니다. 더 얕은 회로들 덕분에 MPF 접근법은 잡음이 있는 하드웨어에 더 적합합니다.

계수는 어떻게 결정되는가?​

MPF 계수에는 두 가지 계열이 있습니다:

정적 계수는 해밀토니안, 초기 상태, 전개 시간과 무관합니다. 이는 선행 Trotter 오차 항의 상쇄를 강제하는 선형계 Ax=bAx = b를 풀어서 구합니다. 2χ2\chi차 대칭 곱 공식과 함께 사용되는 Trotter 스텝 집합 {kj}j=1r\{k_j\}_{j=1}^r에 대해, kjk_j의 역거듭제곱 급수로 Trotter 오차를 전개하면 다음 형태의 제약 방정식이 도출됩니다:

∑j=1rxj=1,∑j=1rxjkjηn=0(n=0,…,r−2),\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),

여기서 정수 지수 {ηn}\{\eta_n\}은 선택된 곱 공식에 대한 연속적인 Trotter 오차 항들의 차수입니다. 대칭 2χ2\chi차 PF의 경우 [S2χ(t/k)]k\left[S_{2\chi}(t/k)\right]^k의 선행 오차는 1/k2χ1/k^{2\chi}로 스케일링되며, 이후 보정 항은 1/k2χ+2,1/k2χ+4,…1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots에서 나타납니다 — 따라서 지수는 ηn=2χ+2n\eta_n = 2\chi + 2n입니다. 비대칭 PF의 경우 홀수와 짝수 거듭제곱이 모두 기여하며 ηn=2χ+n\eta_n = 2\chi + n입니다. 전체 유도는 참고 문헌 [1]을 참조하세요. 위 시스템의 첫 번째 방정식은 비편향성(MPF가 kj→∞k_j \to \infty 극한에서 정확한 기댓값을 재현함)을 보장하며, 나머지 r−1r-1개의 방정식은 처음 r−1r-1개의 Trotter 오차 항을 차례로 상쇄합니다. 결과로 나온 L1L_1-노름 ∥x∥1\|x\|_1이 너무 커서(이는 샘플링 잡음을 증폭시킵니다) 문제가 되는 경우, ∥Ax−b∥\|Ax - b\|를 최소화하면서 ∥x∥1\|x\|_1을 제한하는 근사 최적화를 대신 풀 수 있습니다.

동적 계수 [2], [3]는 추가로 해밀토니안, 초기 상태, 전개 시간 tt에 의존합니다. 이는 실제 시간 전개된 상태와 MPF 근사 사이의 프로베니우스 노름 거리를 최소화합니다:

∥ρ(t)−μD(t)∥F2=1+∑i,jMij(t) xi(t) xj(t)−2∑iLi(t) xi(t),\|\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),

여기서 Mij(t)=Tr[ρki(t) ρkj(t)]M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]는 서로 다른 스텝 수 ki,kjk_i, k_j에 대해 Trotter 전개된 상태들 사이의 겹침의 그람 행렬이며, Li(t)=Tr[ρ(t) ρki(t)]L_i(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]는 (근사) 정확한 상태와의 겹침을 측정합니다. 이 튜토리얼에서는 이러한 양을 텐서 네트워크 방법, 특히 qiskit_addon_mpf의 TeNPy 기반 백엔드를 사용하여 효율적으로 계산합니다.

MPF를 사용해야 할 때​

MPF는 다음의 경우에 가장 유용합니다:

  • 회로 깊이가 병목인 경우. 하드웨어 잡음이 실행 가능한 깊이를 제한하는 경우, MPF를 사용하여 더 얕은 회로로부터 더 높은 유효 Trotter 정확도를 얻을 수 있습니다.

  • 완전한 상태 준비가 아닌 정확한 기댓값이 필요한 경우. MPF는 기댓값 수준에서 작동합니다 — 양자 상태가 아니라 고전적인 숫자를 결합합니다. 따라서 Estimator 프리미티브를 사용한 관측량 추정에 이상적입니다.

  • 적당한 수의 Trotter 스텝 수를 결합하는 경우. 일반적으로 r=3r = 3–55개의 서로 다른 스텝 수 kjk_j를 결합하면 ∥x∥1\|x\|_1을 관리 가능하게 유지하면서 여러 선행 Trotter 오차 항을 상쇄하기에 충분합니다.

MPF가 도움이 되지 않을 수 있는 경우​

  • 매우 짧은 전개 시간. tt가 충분히 작아서 단일 저차 Trotter 공식만으로도 이미 정확한 경우, 여러 회로를 실행하는 오버헤드는 불필요합니다.

  • 상태 준비 작업. MPF는 보정된 양자 상태가 아니라 보정된 기댓값을 생성합니다. 실제 시간 전개된 상태가 필요한 경우(예: 다른 양자 서브루틴의 입력으로), MPF는 적용되지 않습니다.

  • 수렴 영역을 벗어나는 Trotter 스텝 수. 정적 계수 유도는 각 개별 [S2χ(t/kj)]kj\left[S_{2\chi}(t/k_j)\right]^{k_j}를 t/kjt/k_j의 급수로 전개합니다. 이 전개는 t/kmin⁡≲1t/k_{\min} \lesssim 1일 때만 잘 수렴합니다. 주어진 tt에 대해 kmin⁡k_{\min}이 너무 작게 선택되면 가장 얕은 회로가 섭동 영역을 크게 벗어나게 되고, MPF가 상쇄하지 못하고 남긴 고차 오차 항이 커지며, 상쇄를 위해 큰 계수가 필요할 수 있습니다. L1L_1-노름 ∥x∥1\|x\|_1이 실용적인 진단 지표입니다: ∥x∥1≫1\|x\|_1 \gg 1인 경우, 샘플링 오버헤드 ∝∥x∥12\propto \|x\|_1^2가 Trotter 오차 감소보다 클 수 있습니다. 자세한 내용은 Trotter 스텝 선택에 대한 가이드를 참조하세요.

이 튜토리얼에서 다루는 내용​

이 튜토리얼은 두 단계로 처음부터 끝까지의 MPF 워크플로를 안내합니다. 먼저 소규모 시뮬레이터 예제(10큐비트 하이젠베르크 체인)를 통해 문제를 설정하고, 정적 및 동적 MPF 계수를 계산하며, 결과 기댓값을 정확한 대각화와 비교하는 방법을 시연합니다. 그다음 대규모 하드웨어 예제(50큐비트 XXZ 체인)를 통해 트랜스파일하고, 오류 완화를 사용하여 IBM Quantum 하드웨어에서 실행하고, MPF 계수를 사용하여 결과를 후처리하는 방법을 보여줍니다. 전체적으로 표준 Qiskit 도구와 함께 qiskit_addon_mpf 패키지를 사용합니다.

요구 사항​

이 튜토리얼을 시작하기 전에 다음 패키지가 설치되어 있는지 확인하세요:

  • Qiskit SDK v2.0 이상 (시각화 지원 포함)

  • Qiskit Runtime v0.22 이상 (pip install qiskit-ibm-runtime)

  • Qiskit Aer 시뮬레이터 (pip install qiskit-aer)

  • MPF Qiskit 애드온(TeNPy 백엔드 포함) (pip install "qiskit-addon-mpf[tenpy]")

  • Qiskit 애드온 유틸리티 (pip install qiskit-addon-utils)

  • SciPy (pip install scipy)

설정​

아래에서 이 튜토리얼 전체에서 사용되는 모든 패키지 가져오기를 하나의 셀에 모읍니다. 또한 인접한 rxx 및 ryy 회전을 단일 XXPlusYYGate로 융합하는 CollectAndCollapse 트랜스파일러 패스를 정의합니다. 이 패스는 1단계에서 회로를 구성하는 동안(게이트 수를 낮게 유지하기 위해) 적용되며, 4단계에서 동적 MPF를 위한 계층화된 구조를 추출할 때도 간접적으로 적용됩니다(TeNPy는 융합되지 않은 회전 쌍이 아니라 2큐비트 게이트를 필요로 합니다).

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

소규모 시뮬레이터 예제​

1단계: 고전 입력을 양자 문제로 매핑​

선 위의 10큐비트 하이젠베르크 모델로 시작하며, 초기 상태로 Néel 상태 ∣0101…01⟩\vert 0101\ldots01 \rangle을 사용합니다. 해밀토니안은 다음과 같습니다:

H^Heis=J∑i=1L−1(XiXi+1+YiYi+1+ZiZi+1),\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),

여기서 JJ는 최근접 이웃 결합 세기입니다. 체인 중간에 있는 한 쌍의 큐비트에서 ZZ 상관자 ZL/2−1ZL/2Z_{L/2-1} Z_{L/2}를 측정하고, 2차 곱 공식과 함께 Trotter 스텝 kj=[1,2,4]k_j = [1, 2, 4]를 사용합니다.

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 회로 구축하기​

각 시간 지점과 각 Trotter 스텝 수에 대해 근사 Trotter 시간 전개를 구현하는 회로를 생성합니다. 설정 섹션에서 정의한 CollectAndCollapse 패스는 나중에 더 효율적인 텐서 네트워크 시뮬레이션을 준비하기 위해 XX 및 YY 회전을 단일 XX+YY 게이트로 모읍니다.

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

Output of the previous code cell

Step 2: 양자 하드웨어 실행을 위한 문제 최적화​

소규모 예제에서는 Aer 시뮬레이터를 대상으로 합니다. 회로가 실행 준비가 되기 전에 두 가지 변환이 일어납니다:

  1. 해밀토니안 시뮬레이션 수준에서의 게이트 수집. 설정 셀에서 인접한 rxx 및 ryy 회전을 단일 XXPlusYYGate로 융합하는 CollectAndCollapse 패스를 만들었습니다. 1단계에서 Trotter 회로를 구축할 때(pm.run(...) 호출) 이미 이 패스를 적용했습니다. 이는 2큐비트 게이트 수를 줄이는 동시에 나중에 동적 계수 계산을 위한 텐서 네트워크 시뮬레이션에 더 적합한 구조를 만들어냅니다.

  2. 시뮬레이터의 ISA로 낮추기. 아래에서는 각 Trotter 회로를 시뮬레이터의 명령어 집합 아키텍처(ISA)로 낮추기 위해 optimization_level=3으로 Qiskit 사전 설정 패스 매니저를 실행합니다.

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
]

Step 3: Qiskit 프리미티브를 사용한 실행​

소규모 예제에서는 ISA로 낮춰진 Trotter 회로를 Aer가 뒷받침하는 EstimatorV2 프리미티브를 통해 실행합니다. 이렇게 하면 각 (kj,t)(k_j, t) 쌍에 대해 노이즈 없는 참조 값을 얻을 수 있습니다 — 이것이 4단계에서 MPF가 결합할 ⟨A⟩kj(t)\langle A \rangle_{k_j}(t) 값들입니다. 나중에 각 개별 곱 공식과 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])

Step 4: 결과를 원하는 고전적 형식으로 후처리하고 반환하기​

4단계는 MPF가 실제로 구성되는 곳입니다. 계수 xjx_j가 여기서 계산되지만(동적 변형의 경우 이 계산은 집약적일 수 있습니다), 개념적으로 이들은 3단계의 양자 측정을 단일 보정된 기댓값으로 결합하기 위한 고전적인 레시피입니다 — 따라서 전체 계수 및 결합 워크플로를 후처리로 취급합니다.

MPF가 실제 동역학을 얼마나 잘 추적하는지 평가하기 위해, 먼저 해밀토니안을 직접 지수화하여 정확한 시간 전개 기댓값을 계산합니다. 이는 L=10L = 10이기 때문에만 가능합니다. 아래의 대규모 하드웨어 예제에서는 대신 텐서 네트워크 추정치에 의존해야 합니다.

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는 전개 시간, 해밀토니안, 초기 상태와 무관한 계수 xjx_j를 사용합니다. 배경에서 설명한 선형계 Ax=bAx = b를 설정하고 계수를 구합니다. 행렬 AA는 Trotter 스텝 수 kjk_j, 곱 공식의 차수 χ\chi, 그리고 공식이 대칭인지 여부(지수 ηn\eta_n을 제어함)에 의해 결정됩니다.

소규모 예제에서는 비대칭 2χ=22\chi=2차 Suzuki-Trotter 공식과 함께 kj=[1,2,4]k_j = [1, 2, 4]를 사용합니다(따라서 χ=1\chi=1이고 ηn=2+n\eta_n = 2 + n이며, η0=2, η1=3\eta_0 = 2,\, \eta_1 = 3이 됩니다). 시스템은 다음과 같습니다:

A=[11111221421123143],b=[100].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}.

첫 번째 행은 비편향성(∑jxj=1\sum_j x_j = 1)을 강제하며, 두 번째와 세 번째 행은 각각 선행 1/k21/k^2와 다음 차수 1/k31/k^3 Trotter 오차 항을 상쇄합니다.

LSE 설정하기​

위에서 설명한 행렬 AA와 우변 벡터 bb를 조립하기 위해 qiskit_addon_mpf.static의 setup_static_lse를 사용합니다. 행렬 AA는 kjk_j뿐만 아니라 우리가 선택한 곱 공식, 특히 그 차수 χ\chi와 대칭인지 여부에도 의존합니다. symmetric 플래그는 지수 패턴 ηn\eta_n을 제어합니다(대칭 공식은 짝수 거듭제곱의 Trotter 오차 항만 생성합니다. 참고 문헌 [1] 참조). 참고 문헌 [2]에서 보인 것처럼, 기본 PF가 대칭이더라도 symmetric=True로 설정하는 것이 반드시 필요한 것은 아닙니다 — 비대칭 LSE는 여전히 유효합니다(불필요한 추가 제약을 강제할 뿐입니다).

우리 예제에서는 1단계에서 이미 order = 2와 symmetric = False를 설정했습니다.

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

구성된 행렬 AA와 벡터 bb를 검사하여 위에 작성된 시스템과 일치하는지 확인합니다.

lse.A
array([[1.      , 1.      , 1.      ],
       [1.      , 0.25    , 0.0625  ],
       [1.      , 0.125   , 0.015625]])
lse.b
array([1., 0., 0.])

LSE를 확보한 상태에서 lse.solve()를 통해 정적 계수 xjx_j를 구합니다(이는 직접적인 x=A−1bx = 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]
정확한 모델을 이용한 xx 최적화​

x=A−1bx = A^{-1}b를 계산하는 대신, setup_exact_model을 사용하여 LSE를 제약 조건으로 사용하고 최적해가 xx를 산출하는 cvxpy.Problem 인스턴스를 구성할 수 있습니다.

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
근사 모델을 이용한 xx 최적화​

선택한 kjk_j 값 집합에 대한 L1L_1 노름이 너무 높다고 판단될 수 있습니다. 그런 경우이면서 다른 kjk_j 값 집합을 선택할 수 없다면, ∥Ax−b∥\|Ax - b\|를 최소화하면서 L1L_1-노름을 선택한 임계값으로 제한하는 근사 해를 사용할 수 있습니다. 근사 모델 사용 방법 가이드를 확인하세요.

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=1.5
)
model_approx.solve()
print(coeffs_approx.value)
print(
    "L1 norm of the approximate coefficients:",
    np.linalg.norm(coeffs_approx.value, ord=1),
)
[-1.10294118e-03 -2.48897059e-01  1.25000000e+00]
L1 norm of the approximate coefficients: 1.5

동적 MPF 계수​

정적 MPF는 해밀토니안과 상태에 무관한 방식으로 Trotter 오차 항을 상쇄하므로, 주어진 해밀토니안과 초기 상태에 대해 반드시 가능한 가장 작은 근사 오차를 만들어내지는 않습니다. 대신 동적 MPF(참고문헌 [2], [3])는 각 시간 tt에서 Frobenius-노름 거리 ∥ρ(t)−μD(t)∥F2\|\rho(t) - \mu^D(t)\|_F^2를 최소화하는 시간 의존 계수 xi(t)x_i(t)를 찾습니다. 배경에서 보였듯이, 이를 위해서는 Trotter로 시간 전개된 상태 간의 겹침 행렬 Mij(t)M_{ij}(t)와 정확한 상태와의 겹침 Li(t)L_i(t)가 필요합니다 — 둘 다 qiskit_addon_mpf에서 텐서 네트워크(TeNPy) 백엔드를 사용해 추정합니다.

동적 LSE를 설정하려면 세 가지 요소가 필요합니다.

  1. 각 kjk_j에 대해 애드온이 실행하여 ρkj(t)\rho_{k_j}(t)를 MPS/MPO로 생성하는 근사 전개기 팩토리입니다. 이는 2차 Trotter 회로의 레이어 구조(slice_by_depth당 하나의 레이어)로부터 구축되며, TeNPy 절단 매개변수를 사용해 LayerwiseEvolver로 래핑됩니다.

  2. 고정밀 기준 ρ(t)\rho(t)를 생성하는 정확한 전개기 팩토리입니다. 정확한 전개의 대리값으로 작은 시간 단계의 4차 Suzuki-Trotter 회로(dt=0.1, order=4)를 사용합니다.

  3. TeNPy 시뮬레이션의 시드가 되는 항등원 팩토리와 초기 상태 MPS입니다.

아래 셀은 근사 전개기 팩토리를 구성합니다.

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

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

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

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

텐서 네트워크 시뮬레이션의 세부 사항을 결정하는 LayerwiseEvolver의 옵션은 잘못 정의된 최적화 문제가 설정되지 않도록 신중하게 선택해야 합니다.

작은 시간 간격 dt=0.1을 사용하는 4차 Suzuki-Trotter 공식으로 정확한 시간 전개 상태를 근사합니다. TeNPy 절단(truncation) 매개변수는 정확도에 영향을 줄 수 있으므로 다양한 값을 탐색하는 것이 중요합니다.

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

마지막으로, 초기 MPO 상태를 산출하는 identity_factory를 정의하고, 레이어드 Trotter 모델에 사용된 격자와 일치하는 MPS로 Néel 초기 상태를 준비합니다.

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

mps_initial_state = MPS_neel_state(models[0].lat)

팩토리가 준비되면 이제 각 전개 시간에서 동적 계수를 계산합니다. 각 tt에 대해 setup_dynamic_lse는 TeNPy를 통해 관련 겹침 행렬을 구축하고, setup_frobenius_problem은 Frobenius-노름 비용을 최소화하는 cvxpy.Problem을 반환합니다. 솔버는 해당 시간에 맞춘 계수 xj(t)x_j(t)를 반환하며, 이를 mpf_dynamic_coeffs_list에 수집합니다. 특정 tt에서 솔버가 실패하면 루프가 계속되도록 0 계수로 대체합니다.

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 계수 결합​

이제 각 계수 집합(정적-정확, 정적-근사, 동적)에 대해 ⟨A⟩MPF(t)=∑jxj ⟨A⟩kj(t)\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)를 평가하고, 회로별 표준 오차를 전파한 다음, 결과 시계열을 정확한 대각화 곡선과 비교하여 그립니다.

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

Output of the previous code cell

위 그래프는 Trotter 오차와 샘플링 오차 간의 상호작용을 보여줍니다.

  • Trotter 오차. 개별 곱 공식(회색 마커)은 시간이 지남에 따라 정확한 곡선에서 점점 더 벗어납니다. k=1k=1 회로가 가장 얕지만 가장 큰 편차를 보이는데, 이는 이미 t/k≳1t/k \gtrsim 1인 영역에 있어 선행 1/k21/k^{2} 오차 항이 크기 때문입니다. MPF 조합(색상 마커)은 이러한 선행 Trotter 오차 항 중 여러 개를 상쇄하므로, 어떤 단일 kjk_j 회로보다 정확한 곡선을 훨씬 더 밀접하게 따라갑니다. 남은 간극은 MPF가 상쇄하지 않는 고차 Trotter 항을 반영합니다. 2차, r=3r=3의 정적 MPF는 처음 두 오차 차수만 제거하며, 큰 t/kmin⁡t/k_{\min}에서는 상쇄되지 않은 나머지가 결국 지배하게 됩니다 — 따라서 MPF는 매우 얕은 회로가 임의의 시간에서 정확성을 유지한다는 것을 보장하지 않습니다.

  • 샘플링 오차. MPF 곡선의 더 넓은 오차 막대는 선형 결합의 직접적인 결과입니다. 독립적인 회로별 표준 오차 σkj\sigma_{k_j}를 전파하면 총 분산은 σMPF2=∑jxj2 σkj2\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2가 됩니다. 따라서 ∥x∥2\|x\|_2(그리고 실제로 우리가 제어하는 ∥x∥1\|x\|_1)가 클수록 주어진 목표 불확실성에 도달하기 위해 더 많은 샷이 필요합니다. 이것이 배경에서 다룬 근사 솔버 옵션 뒤에 있는 트레이드오프입니다. 즉, 이 오버헤드를 관리 가능하게 유지하기 위해 ∥x∥1\|x\|_1을 제한합니다. 결정적으로, Trotter 오차와 달리 샘플링 오차는 1/Nshots1/\sqrt{N_{\text{shots}}}로 줄어들므로, 더 많은 샷을 사용하여 항상 줄일 수 있습니다.

아래의 대규모 하드웨어 예제에서는 하드웨어 노이즈가 각 ⟨A⟩kj\langle A \rangle_{k_j}에 추가 오차 원천으로 유입되며, 이 역시 MPF 계수에 의해 증폭됩니다. 해당 섹션에서 오차 완화가 MPF와 어떻게 상호작용하는지 살펴보겠습니다.

대규모 하드웨어 예제​

이 섹션에서는 정확하게 시뮬레이션할 수 있는 범위를 넘어 문제를 확장합니다. 시간 t=3t = 3에서 50-큐비트 XXZ 사슬을 사용하여 참고문헌 [3]에 나온 일부 결과를 재현합니다. 소규모 예제와 동일한 4단계 워크플로를 따르되, 이번에는 오차 완화를 적용하여 실제 양자 하드웨어를 대상으로 합니다. 템플릿에서와 마찬가지로 각 단계는 코드에 인라인으로 표시되어 있으며, 중간 출력을 검토할 가치가 있을 때는 한 단계가 여러 셀에 걸쳐 있을 수 있습니다. 매핑은 소규모 예제를 그대로 따릅니다. 즉, 해밀토니안을 정의하고, Trotter 매개변수를 선택하고, MPF 계수(정적 및 동적)를 계산하고, 회로를 구축합니다. 주요 차이점은 다음과 같습니다.

  • U(0.5,1.5)\mathcal{U}(0.5, 1.5)에서 무작위로 추출한 결합을 갖는 50개 사이트의 XXZ 해밀토니안입니다(참고문헌 [3]).

  • kj=[3,4,6]k_j = [3, 4, 6]을 갖는 대칭 2차 Trotter 공식입니다(따라서 χ=1\chi=1, symmetric=True).

  • 단일 고정 전개 시간 t=3t = 3입니다. kmin⁡=3k_{\min}=3이므로 t/kmin⁡=1t/k_{\min}=1이 되어, 얕은 구성 요소들이 MPF가 의존하는 선행 오차 모델이 유효한 Trotter 수렴 영역 내에 유지됩니다.

  • 기준선으로 사용되는 k=10k = 10 Trotter 단계를 갖는 추가 단일 회로 비교 실행입니다. k=10k = 10을 선택한 이유는 하드웨어에서의 2-큐비트 깊이가 가장 깊은 MPF 구성 요소(kmax⁡=6k_{\max}=6)에 여러 MPF 회로를 실행하는 오버헤드를 더한 것보다 더 깊기 때문입니다 — 이는 노이즈로 제한되는 영역에 충분히 깊어서, MPF 조합이 단일 회로 기준선보다 우수한 성능을 보일 것으로 예상되는 영역입니다. 이는 MPF의 유효 Trotter 오차를 목표로 하는 회로(훨씬 더 많은 단계가 필요함)가 아니라, MPF 조합에 대한 "단일 깊은 회로" 비교입니다.

여기서는 아직 1단계(매핑과 회로 구성)에 있지만, 이 셀에서 정적 계수와 함께 동적 계수도 미리 계산한다는 점에 유의하세요. 동적 계수는 HH와 tt에 의존하지만 양자 측정에는 의존하지 않으므로, 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

이제 선택한 백엔드에 맞게 회로를 최적화합니다. optimization_level=3의 Qiskit 사전 설정 패스 매니저를 사용하며, 이는 자동으로 좋은 물리적 큐비트 집합을 선택하고 각 회로를 디바이스 토폴로지에 라우팅합니다.

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

실제 하드웨어에서 더 깊은 회로를 실행하려면 공격적인 오차 완화가 필요합니다. 동적 디커플링, 게이트 및 측정 트와이링, 측정 오차 완화, 그리고 제로-노이즈 외삽(ZNE)을 활성화합니다. 여기서 사용하는 ZNE 노이즈 계수(1, 1.2, 1.4)는 얕은 회로 시나리오보다 작다는 점에 유의하세요. 더 깊은 MPF 구성 요소가 이미 노이즈 임계값에 가까워, 큰 노이즈 증폭이 ZNE 외삽이 신뢰할 수 있는 지점을 넘어서게 만들 수 있기 때문입니다.

네 개의 회로 모두(kj=[3,4,6]k_j = [3, 4, 6]의 세 MPF 구성 요소에 k=10k = 10 기준선을 더한 것)를 단일 Estimator 작업으로 제출합니다.

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

작업 결과에서 회로별 기댓값과 표준 편차를 가져온 다음, 소규모 예제에서와 정확히 동일하게 각 MPF 계수 집합과 결합합니다: ⟨A⟩MPF=∑jxj ⟨A⟩kj\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}, 전파된 분산은 σ2=∑jxj2σkj2\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()

Output of the previous code cell

위 하드웨어 결과에 대한 몇 가지 관찰 사항입니다.

  • 더 깊이 들어가는 것이 하드웨어에서는 공짜가 아닙니다. 단일 회로 기준선이 그 이야기를 직접 보여줍니다. k=6k = 6 회로는 본질적으로 정확한 반면(−0.256-0.256 대 기준값 −0.244-0.244), 더 깊은 k=10k = 10 기준선은 더 좋아지는 것이 아니라 더 나빠집니다(−0.061-0.061, ∼0.18\sim 0.18 벗어남). Trotter 오차가 이미 작을 때 단계를 추가하는 것은 대부분 회로를 깊게만 만들고 더 많은 게이트 노이즈와 결맞음 붕괴를 누적시킵니다. 이것이 바로 MPF가 만들어진 영역입니다. 즉, 얕은 구성 요소만을 사용하여 깊은 회로의 정확도에 도달하는 것입니다.

  • 작은 노름의 MPF가 깊은 단일 회로를 능가합니다. 근사-정적 MPF(∥x∥1≈2\|x\|_1 \approx 2로 제한)는 −0.259-0.259에 도달하며, 기준값에서 ∼0.015\sim 0.015 이내이고 k=10k = 10 기준선보다 훨씬 가깝습니다. 동적 MPF(−0.127-0.127)도 그 기준선을 편안하게 능가합니다. 둘 다 얕은 kj=[3,4,6]k_j = [3, 4, 6] 회로만 결합했음에도, 깊은 단일 회로가 얻지 못한 답을 회복합니다.

  • 계수 노름이 수학적 최적성보다 더 중요합니다. 정확-정적 MPF는 ∥x∥1=4.66\|x\|_1 = 4.66이며 모든 것 중 가장 나쁜 추정치입니다(−0.567-0.567, 0.30.3 이상 벗어남). 큰 계수 노름은 각 ⟨A⟩kj\langle A \rangle_{k_j}의 잔여 게이트 노이즈, 결맞음 붕괴, ZNE 오차를 거의 같은 배율로 증폭시켜, 그 대가로 얻는 Trotter 오차 상쇄 효과를 압도합니다. 노름을 제한하면(근사-정적 솔버, ∥x∥1≈2\|x\|_1 \approx 2) 이 압도 효과가 사라지고 최선의 추정치를 제공합니다 — 비록 그 계수가 더 이상 선행 Trotter 오차를 정확히 상쇄하지 않더라도 말입니다.

  • 개별 얕은 회로도 여전히 경쟁력이 있을 수 있습니다. 단독 k=6k = 6 구성 요소(−0.256-0.256)는 여기서도 본질적으로 정확하며, 이번 실행에서는 근사-정적 MPF보다 근소하게 더 가깝기까지 합니다. 문제는 "수렴했지만 아직 노이즈로 제한되지 않은" 최적 지점에 어떤 단일 kk가 위치하는지 미리 알 수 없다는 것이며, Trotter 수렴을 보장하기 위해 단순히 더 깊게 가는(k=10k = 10) 안전해 보이는 선택이 바로 실패하는 선택이라는 것입니다. MPF는 올바른 깊이를 추측할 필요 없이 얕은 회로들의 원칙적인 조합을 제공합니다.

실용적인 결론은, 하드웨어에서 MPF는 각 개별 ⟨A⟩kj\langle A \rangle_{k_j}에 대한 강력한 오차 완화와 함께 사용해야 하며, 계수 L1L_1-노름은 적당하게 유지해야 하고(근사 솔버 또는 동적 MPF 사용), Trotter 단계 kjk_j는 t/kmin⁡≲1t/k_{\min} \lesssim 1이 되도록 선택해야 한다는 것입니다 — 여기서는 t=3t = 3에서 kmin⁡=3k_{\min} = 3이 t/kmin⁡=1t/k_{\min} = 1을 주어, 정적 MPF가 의존하는 선행 오차 모델이 유효한 수렴 영역 내에 구성 요소들을 유지합니다. 이러한 선택을 하면, 여기서 작은 노름의 MPF는 수렴된 단일 회로와 일치하는 반면 순진하게 "그냥 더 깊게" 가는 기준선은 그렇지 않으며, 참고문헌 [3]에 나온 깊이 대 정확도 이점을 회복합니다. 또한 개별 실행에는 노이즈가 있다는 점에 유의하세요 — 동일 작업의 다른 제출(또는 다른 백엔드)에서는 정확한 순서가 바뀔 수 있습니다. 견고한 경향은 작은 ∥x∥1\|x\|_1의 MPF가 잘 작동하고, 큰 ∥x∥1\|x\|_1의 정확-정적 MPF는 하드웨어 노이즈에 의해 증폭되며, 지나치게 깊은 단일 회로는 노이즈로 제한된다는 것입니다.

다음 단계​

권장 사항

이 작업이 흥미로웠다면, 다음 자료에도 관심이 있을 수 있습니다:

참고 문헌​

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

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

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

출처: IBM Quantum docs — 2026년 9월 9일 업데이트
doQumentation의 영어 버전 — 2026년 9월 9일 업데이트
사이트 또는 번역 문제?콘텐츠 문제? IBM Quantum docs에서 편집