트로터 오류를 줄이기 위한 다중 곱 공식

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

배경

이 튜토리얼에서는 다중 곱 공식(MPF, Multi-Product Formula)을 사용하여 실제로 실행할 가장 깊은 트로터 Circuit으로 인한 오류보다 관측값에 대한 트로터 오류를 낮추는 방법을 설명합니다. MPF는 여러 Circuit 실행의 가중 조합을 통해 해밀토니안 동역학의 트로터 오류를 줄입니다. 해밀토니안 $H$에 대해 양자 상태 $\rho(t)=e^{-i H t} \rho(0) e^{i H t}$의 관측값 기댓값을 구하는 작업을 생각해 보겠습니다. 곱 공식(PF, Product Formula)을 사용하면 다음과 같이 시간 발전 $e^{-i H t}$를 근사할 수 있습니다:

• 해밀토니안 $H$를 $H=\sum_{a=1}^d F_a$로 작성합니다. 여기서 $F_a$는 에르미트 연산자로, 각 대응 유니터리가 양자 디바이스에서 효율적으로 구현될 수 있습니다.
• 서로 가환하지 않는 항 $F_a$를 근사합니다.

그러면 1차 PF(Lie-Trotter 공식)는 다음과 같습니다:

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

이는 2차 오류 항 $S_1(t)=e^{-i H t}+\mathcal{O}\left(t^{2}\right)$을 가집니다. 더 빠르게 수렴하는 고차 PF(Lie-Trotter-Suzuki 공식)도 사용할 수 있으며, 이는 다음과 같이 재귀적으로 정의됩니다:

$S_2(t):=\prod_{a=1}^d e^{-i F_a t/2}\prod_{a=1}^d e^{-i F_a t/2}$

$S_{2 \chi}(t):= S_{2 \chi -2}(s_{\chi}t)^2 S_{2 \chi -2}((1-4s_{\chi})t)S_{2 \chi -2}(s_{\chi}t)^2,$

여기서 $\chi$는 대칭 PF의 차수이고 $s_p = \left( 4 - 4^{1/(2p-1)} \right)^{-1}$입니다. 장시간 발전의 경우 시간 구간 $t$를 $t/k$ 길이의 $k$개 구간(트로터 스텝)으로 나누고, 각 구간에서 $\chi$차 곱 공식 $S_{\chi}$로 시간 발전을 근사할 수 있습니다. 따라서 $k$개의 트로터 스텝에 걸친 시간 발전 연산자에 대한 $\chi$차 PF는 다음과 같습니다:

$S_{\chi}^{k}(t) = \left[ S_{\chi} \left( \frac{t}{k} \right)\right]^k = e^{-i H t}+O\left(t \left( \frac{t}{k} \right)^{\chi} \right)$

여기서 오류 항은 트로터 스텝 수 $k$와 PF의 차수 $\chi$가 증가할수록 감소합니다.

정수 $k \geq 1$과 곱 공식 $S_{\chi}(t)$가 주어지면, 근사 시간 발전 상태 $\rho_k(t)$는 $\rho_0$에 곱 공식 $S_{\chi}\left(\frac{t}{k}\right)$를 $k$번 적용하여 얻을 수 있습니다.

$$\rho_k(t)=S_{\chi}\left(\frac{t}{k}\right)^k \rho_0 S_{\chi}\left(\frac{t}{k}\right)^{-k}$$

$\rho_k(t)$는 트로터 근사 오류 ||$\rho_k(t)-\rho(t) ||$를 가지는 $\rho(t)$의 근사값입니다. $\rho(t)$의 트로터 근사들의 선형 조합을 고려하면:

$$\mu(t) = \sum_{j}^{l} x_j \rho^{k_j}_{j}\left(\frac{t}{k_j}\right) + \text{some remaining Trotter error} \, ,$$

여기서 $x_j$는 가중 계수, $\rho^{k_j}_j$는 곱 공식 $S^{k_j}_{\chi}$로 초기 상태를 $k_j$개의 트로터 스텝만큼 발전시켜 얻은 순수 상태에 해당하는 밀도 행렬이며, $j \in {1, ..., l}$은 MPF를 구성하는 PF의 수에 대한 인덱스입니다. $\mu(t)$의 모든 항은 동일한 곱 공식 $S_{\chi}(t)$를 기반으로 합니다. 목표는 $\|\mu(t)-\rho(t)\|$가 더 낮은 $\mu(t)$를 찾아 ||$\rho_k(t)-\rho(t) \|$를 개선하는 것입니다.

• $x_i$가 양수일 필요가 없으므로 $\mu(t)$는 물리적 상태일 필요가 없습니다. 여기서의 목표는 관측값의 기댓값 오류를 최소화하는 것이며, $\rho(t)$의 물리적 대체물을 찾는 것이 아닙니다.
• $k_j$는 Circuit 깊이와 트로터 근사 수준 모두를 결정합니다. $k_j$ 값이 작을수록 Circuit이 짧아져 Circuit 오류는 줄지만 원하는 상태에 대한 근사는 덜 정확해집니다.

핵심은 $\mu(t)$로 주어지는 나머지 트로터 오류가 가장 큰 $k_j$ 값을 단순히 사용했을 때 얻는 트로터 오류보다 작다는 점입니다.

이 방법의 유용성은 두 가지 관점에서 볼 수 있습니다:

1. 실행 가능한 고정된 트로터 스텝 예산에 대해, 총 트로터 오류가 더 작은 결과를 얻을 수 있습니다.
2. 실행하기에 너무 큰 목표 트로터 스텝 수가 주어졌을 때, MPF를 사용하여 유사한 트로터 오류를 달성하는 더 낮은 깊이의 Circuit 컬렉션을 실행할 수 있습니다.

요구 사항

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

• Qiskit SDK v1.0 이상 (시각화 지원 포함)
• Qiskit Runtime v0.22 이상 (pip install qiskit-ibm-runtime)
• MPF Qiskit 애드온 (pip install qiskit_addon_mpf)
• Qiskit 애드온 유틸리티 (pip install qiskit_addon_utils)
• Quimb 라이브러리 (pip install quimb)
• Qiskit Quimb 라이브러리 (pip install qiskit-quimb)
• 패키지 간 호환성을 위한 Numpy v0.21 (pip install numpy==0.21)

파트 I. 소규모 예제

MPF의 안정성 탐색

MPF 상태 $\mu(t)$를 구성하는 트로터 스텝 수 $k_j$의 선택에는 명백한 제한이 없습니다. 하지만 $\mu(t)$로부터 계산되는 기댓값의 불안정성을 피하기 위해 신중하게 선택해야 합니다. 일반적으로 좋은 규칙은 최소 트로터 스텝 $k_{\text{min}}$을 $t/k_{\text{min}} \lt 1$이 되도록 설정하는 것입니다. 이에 대한 자세한 내용과 다른 $k_j$ 값 선택 방법은 MPF에 대한 트로터 스텝 선택 방법 가이드를 참고하세요.

아래 예제에서는 다양한 시간 범위에 걸쳐 서로 다른 시간 발전 상태를 사용하여 자화의 기댓값을 계산함으로써 MPF 솔루션의 안정성을 탐색합니다. 구체적으로, 각 트로터 스텝으로 구현된 근사 시간 발전에서 계산된 기댓값과 다양한 MPF 모델(정적 및 동적 계수)의 기댓값을 시간 발전된 관측값의 정확한 값과 비교합니다. 먼저 트로터 공식의 매개변수와 발전 시간을 정의해 보겠습니다.

# 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 rustworkx scipy

import numpy as np

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)

이 예제에서는 초기 상태로 Neel 상태 $\vert \text{Neel} \rangle = \vert 0101...01 \rangle$를 사용하고, 10개 사이트로 이루어진 선형 Heisenberg 모델을 시간 발전을 지배하는 해밀토니안으로 사용합니다.

$$\hat{\mathcal{H}}_{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$는 최근접 이웃 간의 결합 강도입니다.

from qiskit.transpiler import CouplingMap
from rustworkx.visualization import graphviz_draw
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian
import numpy as np

L = 10

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

# Get a qubit operator describing the Heisenberg field model
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])

측정할 관측값은 체인 중간에 있는 Qubit 쌍의 자화입니다.

from qiskit.quantum_info import SparsePauliOp

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

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

Circuit에서 XX 및 YY 회전을 단일 XX+YY Gate로 수집하는 Transpiler 패스를 정의합니다. 이를 통해 MPO 계산 중 TeNPy의 스핀 보존 특성을 활용하여 계산 속도를 크게 높일 수 있습니다.

from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler import PassManager
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)
from functools import partial

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

그런 다음 근사 트로터 시간 발전을 구현하는 Circuit을 생성합니다.

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

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

다음으로 트로터 Circuit에서 시간 발전된 기댓값을 계산합니다.

from copy import deepcopy
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

aer_sim = AerSimulator()
estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for t, mpf_circuits in zip(trotter_times, all_circs):
    mpf_expvals = []
    circuits = [deepcopy(circuit) for circuit in mpf_circuits]
    pm_sim = generate_preset_pass_manager(
        backend=aer_sim, optimization_level=3
    )
    isa_circuits = pm_sim.run(circuits)
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals = [res.data.evs for res in result]
    mpf_stds = [res.data.stds for res in result]
    mpf_expvals_all_times.append(mpf_expvals)
    mpf_stds_all_times.append(mpf_stds)

비교를 위해 정확한 기댓값도 계산합니다.

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

exact_expvals = []
for t in exact_evolution_times:
    # Exact expectation values
    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는 $x_j$ 값이 진화 시간 $t$에 의존하지 않는 경우입니다. 차수 $\chi = 1$인 PF에서 $k_j$ 트로터 스텝을 사용하는 경우를 고려하면 다음과 같이 쓸 수 있습니다:

$S_1^{k_j}\left( \frac{t}{k_j} \right)=e^{-i H t}+ \sum_{n=1}^{\infty} A_n \frac{t^{n+1}}{k_j^n}$

여기서 $A_n$은 해밀토니안 분해에서 $F_a$ 항들의 교환자(commutator)에 의존하는 행렬입니다. $A_n$ 자체는 시간과 트로터 스텝 수 $k_j$에 독립적이라는 점이 중요합니다. 따라서 가중치 $x_j$를 신중하게 선택하면 $\mu(t)$에 기여하는 낮은 차수의 오류 항을 상쇄할 수 있습니다. $\mu(t)$ 표현식에서 처음 $l-1$개 항(트로터 스텝 수가 적어 가장 큰 기여를 하는 항)의 트로터 오류를 상쇄하려면, 계수 $x_j$가 다음 방정식을 만족해야 합니다:

$\sum_{j=1}^l x_j = 1$ $\sum_{j=1}^{l-1} \frac{x_j}{k_j^{n}} = 0$

여기서 $n=0, ... l-2$입니다. 첫 번째 방정식은 구성된 상태 $\mu(t)$에 편향이 없음을 보장하며, 두 번째 방정식은 트로터 오류의 상쇄를 보장합니다. 고차 PF의 경우, 두 번째 방정식은 $\sum_{j=1}^{l-1} \frac{x_j}{k_j^{\eta}} = 0$이 되며, 여기서 $\eta = \chi + 2n$ (대칭 PF의 경우) 또는 $\eta = \chi + n$ (그 외의 경우), $n=0, ..., l-2$입니다. 결과적인 오류 (참고문헌 [1], [2])는 다음과 같습니다:

$\epsilon = \mathcal{O} \left( \frac{t^{l+1}}{k_1^l} \right).$

주어진 $k_j$ 값 집합에 대한 정적 MPF 계수를 결정하는 것은, 변수 $x_j$에 대해 위 두 방정식으로 정의된 선형 연립방정식 $Ax=b$를 푸는 것과 같습니다. 여기서 $x$는 구하고자 하는 계수, $A$는 $k_j$와 사용하는 PF 유형($S$)에 의존하는 행렬, $b$는 제약 조건 벡터입니다. 구체적으로:

$A_{0,j} = 1$ $A_{i>0,j} = k_{j}^{-(\chi + s(i-1))}$ $b_0 = 1$ $b_{i>0} = 0$

여기서 $\chi$는 order(차수), $s$는 symmetric이 True이면 $2$, 아니면 $1$, $k_{j}$는 trotter_steps, $x$는 풀어야 할 변수입니다. 인덱스 $i$와 $j$는 $0$에서 시작합니다. 이를 행렬 형태로도 시각화할 수 있습니다:

$$A = \begin{bmatrix} A_{0,0} & A_{0,1} & A_{0,2} & ... \\ A_{1,0} & A_{1,1} & A_{1,2} & ... \\ A_{2,0} & A_{2,1} & A_{2,2} & ... \\ ... & ... & ... & ... \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & ... \\ k_{0}^{-(\chi + s(1-1))} & k_{1}^{-(\chi + s(1-1))} & k_{2}^{-(\chi + s(1-1))} & ... \\ k_{0}^{-(\chi + s(2-1))} & k_{1}^{-(\chi + s(2-1))} & k_{2}^{-(\chi + s(2-1))} & ... \\ ... & ... & ... & ... \end{bmatrix}$$

그리고

$$b = \begin{bmatrix} b_{0} \\ b_{1} \\ b_{2} \\ ... \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ ... \end{bmatrix}$$

자세한 내용은 선형 연립방정식(LSE) 문서를 참조하세요.

$x = A^{-1}b$로 분석적 해를 구할 수 있습니다. 예를 들어 참고문헌 [1] 또는 [2]를 참고하세요. 단, 이 정확한 해는 "불량 조건(ill-conditioned)"일 수 있어 계수 $x$의 L1-노름이 매우 커질 수 있으며, 이는 MPF의 성능 저하로 이어질 수 있습니다. 대신, L1-노름을 최소화하는 근사 해를 구하여 MPF 동작을 최적화하려 시도할 수도 있습니다.

LSE 설정

$k_j$ 값을 선택한 후에는 위에서 설명한 대로 먼저 LSE인 $Ax=b$를 구성해야 합니다. 행렬 $A$는 $k_j$뿐만 아니라 PF의 선택, 특히 그 _차수(order)_에도 의존합니다. 또한 symmetric=True/False를 설정하여 PF의 대칭 여부를 고려할 수 있습니다(참고문헌 [1] 참조). 단, 참고문헌 [2]에서 보여주듯이 이는 필수는 아닙니다.

from qiskit_addon_mpf.static import setup_static_lse

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

위에서 선택한 값들로 $A$ 행렬과 $b$ 벡터를 구성해 보겠습니다. $j=0,1,2$인 트로터 스텝 $k_j = [1, 2, 4]$, 차수 $\chi = 2$, 비대칭 트로터 스텝 선택($s=1$)에서, 첫 번째 행 아래의 $A$ 행렬 원소는 $A_{i>0,j} = k_{j}^{-(2 + 1(i-1))}$ 표현식으로 결정됩니다:

$A_{0,0} = A_{0,1} = A_{0,2} = 1$ $A_{1,j} = k_{j}^{-1} \rightarrow A_{1,0} = \frac{1}{1^2}, \;, A_{1,1} = \frac{1}{2^2}, \;, A_{1,2} = \frac{1}{4^2}$ $A_{2,j} = k_{j}^{-2} \rightarrow A_{2,0} = \frac{1}{1^3}, \;, A_{2,1} = \frac{1}{2^3}, \;, A_{2,2} = \frac{1}{4^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}$$

이는 lse 객체를 검사하여 확인할 수 있습니다:

lse.A

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

제약 조건 벡터 $b$의 원소는 다음과 같습니다: $b_{0} = 1$ $b_1 = b_2 = 0$

따라서,

$$b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$$

마찬가지로 lse에서도 확인할 수 있습니다:

lse.b

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

lse 객체에는 연립방정식을 만족하는 정적 계수 $x_j$를 찾는 메서드가 있습니다.

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$ 최적화

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

다음 섹션에서 이 인터페이스가 왜 존재하는지 명확해질 것입니다.

from qiskit_addon_mpf.costs import setup_exact_problem

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

[ 0.04761905 -0.57142857  1.52380952]

이 계수들로 구성된 MPF가 좋은 결과를 낼지 여부의 지표로 L1-노름을 사용할 수 있습니다(참고문헌 [1] 참조).

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

L1 norm of the exact coefficients: 2.1428571428556378

근사 모델을 이용한 $x$ 최적화

선택한 $k_j$ 값 집합에 대한 L1 노름이 너무 높다고 판단될 수 있습니다. 그런 경우이면서 다른 $k_j$ 값 집합을 선택할 수 없다면, 정확한 해 대신 근사 해를 사용할 수 있습니다.

이를 위해 setup_approximate_model을 사용하여 L1-노름을 선택한 임계값으로 제한하면서 $Ax$와 $b$의 차이를 최소화하는 다른 cvxpy.Problem 인스턴스를 구성합니다.

from qiskit_addon_mpf.costs import setup_sum_of_squares_problem

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

이 최적화 문제를 해결하는 방법에는 완전한 자유가 있으므로, 최적화 솔버, 수렴 임계값 등을 변경할 수 있습니다. 근사 모델 사용 방법 가이드를 확인하세요.

Dynamic MPF coefficients

이전 섹션에서는 표준 Trotter 근사를 개선하는 정적 MPF를 소개했습니다. 그러나 이 정적 버전이 반드시 근사 오류를 최소화하는 것은 아닙니다. 구체적으로, $\mu^S(t)$로 표기되는 정적 MPF는 곱 공식 상태 $\{\rho_{k_i}(t)\}_{i=1}^r$로 생성되는 부분 공간으로의 $\rho(t)$의 최적 사영이 아닐 수 있습니다.

이를 해결하기 위해, Frobenius 노름 의미에서 근사 오류를 최소화하는 동적 MPF(참고문헌 [2]에서 도입되었고 참고문헌 [3]에서 실험적으로 검증됨)를 고려합니다. 형식적으로, 다음을 최소화하는 데 집중합니다:

$$\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr],$$

각 시간 $t$에서의 계수 $x_i(t)$에 대해. Frobenius 노름에서의 최적 사영자는 $\mu^D(t) \;=\; \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t)$이며, 이를 동적 MPF라고 합니다. 위의 정의를 대입하면:

$$\|\rho(t) - \mu^D(t)\|_F^2 \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \mu^D(t)\right)^2 \bigr] \;=\; \\ = \mathrm{Tr}\bigl[ \left( \rho(t) - \sum_{i=1}^r x_i(t)\,\rho_{k_i}(t) \right) \left( \rho(t) - \sum_{j=1}^r x_j(t)\,\rho_{k_j}(t) \right) \bigr] \;=\; \\ = 1 \;+\; \sum_{i,j=1}^r M_{i,j}(t)\,x_i(t)\,x_j(t) \;-\; 2 \sum_{i=1}^r L_i^{\mathrm{exact}}(t)\,x_i(t),$$

여기서 $M_{i,j}(t)$는 *그람 행렬(Gram matrix)*로, 다음과 같이 정의됩니다:

$$M_{i,j}(t) \;=\; \mathrm{Tr}\bigl[\rho_{k_i}(t)\,\rho_{k_j}(t)\bigr] \;=\; \bigl|\langle \psi_{\mathrm{in}} \!\mid S\bigl(t/k_i\bigr)^{-k_i}\,S\bigl(t/k_j\bigr)^{k_j} \!\mid \psi_{\mathrm{in}} \rangle \bigr|^2.$$

그리고

$$L_i^{\mathrm{exact}}(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$$

는 정확한 상태 $\rho(t)$와 각 곱 공식 근사 $\rho_{k_i}(t)$ 사이의 중첩(overlap)을 나타냅니다. 실제 상황에서 이 중첩들은 노이즈나 $\rho(t)$에 대한 부분적 접근으로 인해 근사적으로만 측정될 수 있습니다.

여기서 $\lvert\psi_{\mathrm{in}}\rangle$은 초기 상태이고, $S(\cdot)$는 곱 공식에서 적용되는 연산입니다. 이 식을 최소화하는 계수 $x_i(t)$를 선택하고 ($\rho(t)$가 완전히 알려지지 않은 경우 근사적인 중첩 데이터를 처리하여), MPF 부분 공간 내에서 $\rho(t)$의 (Frobenius 노름 의미의) "최상의" 동적 근사를 얻습니다. $L_i(t)$와 $M_{i,j}(t)$의 양은 텐서 네트워크 방법으로 효율적으로 계산할 수 있습니다 [3]. MPF Qiskit 애드온은 이 계산을 수행하기 위한 여러 "백엔드"를 제공합니다. 아래 예시는 가장 유연한 방법을 보여주며, TeNPy 레이어 기반 백엔드 문서에도 자세히 설명되어 있습니다. 이 방법을 사용하려면, 원하는 시간 진화를 구현하는 Circuit에서 시작하여 해당 Circuit의 레이어에서 이러한 연산을 나타내는 모델을 생성합니다. 마지막으로, 시간 진화된 양 $M_{i,j}(t)$와 $L_i(t)$를 생성하는 데 사용할 수 있는 Evolver 객체를 생성합니다. 먼저 Circuit으로 구현된 근사 시간 진화(ApproxEvolverFactory)에 해당하는 Evolver 객체를 생성합니다. 특히, order 변수가 일치하는지 각별히 주의하세요. 근사 시간 진화에 해당하는 Circuit을 생성할 때, time = 1.0 및 Trotter 단계 수(reps=1)에 대한 플레이스홀더 값을 사용합니다. 올바른 근사 Circuit은 이후 setup_dynamic_lse의 동적 문제 해결자에 의해 생성됩니다.

from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.backends.tenpy_layers import LayerModel
from qiskit_addon_mpf.backends.tenpy_layers import LayerwiseEvolver
from functools import partial

# 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의 옵션은 잘못 정의된 최적화 문제를 설정하지 않도록 신중하게 선택해야 합니다.

그런 다음 정확한 진화자(예: ExactEvolverFactory)를 설정합니다. 이는 실제 또는 "참조" 시간 진화를 계산하는 Evolver 객체를 반환합니다. 현실적으로는, 더 높은 차수의 Suzuki–Trotter 공식 또는 작은 시간 단계를 사용하는 다른 신뢰할 수 있는 방법으로 정확한 진화를 근사할 것입니다. 아래에서는 작은 시간 단계 dt=0.1을 사용하는 4차 Suzuki-Trotter 공식으로 정확한 시간 진화 상태를 근사하는데, 이는 시간 $t$에서 사용되는 Trotter 단계 수가 $k=t/dt$임을 의미합니다. 또한 기저 텐서 네트워크의 최대 본드 차원과 분할된 텐서 네트워크 본드의 최소 특이값을 제한하기 위해 일부 TeNPy 전용 절단 옵션을 지정합니다. 이 매개변수들은 동적 MPF 계수로 계산된 기댓값의 정확도에 영향을 미칠 수 있으므로, 계산 시간과 정확도 사이의 최적 균형을 찾기 위해 다양한 값 범위를 탐색하는 것이 중요합니다. MPF 계수의 계산은 하드웨어 실행으로 얻은 PF의 기댓값에 의존하지 않으므로, 후처리에서 조정할 수 있습니다.

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

다음으로, TeNPy와 호환되는 형식(예: MPS_neel_state=$\vert 0101...01 \rangle$)으로 시스템의 초기 상태를 생성합니다. 이를 통해 시간에 따라 진화시킬 다체 파동 함수 $\lvert\psi_{\mathrm{in}}\rangle$을 텐서로 설정합니다.

from qiskit_addon_mpf.backends.tenpy_tebd import MPOState
from qiskit_addon_mpf.backends.tenpy_tebd import MPS_neel_state

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

mps_initial_state = MPS_neel_state(models[0].lat)

각 시간 단계 $t$에 대해 setup_dynamic_lse 방법으로 동적 연립 방정식을 설정합니다. 해당 객체에는 동적 MPF 문제에 관한 정보가 포함되어 있습니다: lse.A는 그람 행렬 $M$을 제공하고, lse.b는 중첩 $L$을 제공합니다. 그런 다음 setup_frobenius_problem을 사용하여 LSE를 풀어(비정상적이지 않을 때) 동적 계수를 구할 수 있습니다. 곱 공식의 세부 사항에만 의존하고 시간 진화의 세부 사항(해밀토니안과 초기 상태)과 독립적인 정적 계수와의 차이에 주목하는 것이 중요합니다.

from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import setup_frobenius_problem

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

마지막으로, 진화 시간 전체에 걸쳐 이 기댓값들을 플롯합니다.

import matplotlib.pyplot as plt

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
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",
    color="orange",
    label="Static MPF - Approximate",
)

# # 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",
    color="pink",
    label="Dynamic MPF",
)

plt.plot(
    exact_evolution_times,
    exact_expvals,
    lw=3,
    color="red",
    label="Exact time-evolution",
)

plt.title(
    f"Expectation values for (ZZ,{(L//2-1, L//2)}) as a function of time"
)
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend()
plt.grid()

위의 예시처럼 $k=1$ PF가 모든 시간에서 제대로 동작하지 않는 경우, 동적 MPF 결과의 품질도 크게 영향을 받습니다. 이러한 상황에서는 더 많은 Trotter 단계를 가진 개별 PF를 사용하여 전체 결과 품질을 개선할 가능성을 조사하는 것이 유용합니다. 이 시뮬레이션에서, 우리는 서로 다른 유형의 오류들 사이의 상호 작용을 확인할 수 있습니다: 유한 샘플링으로 인한 오류와 곱 공식으로 인한 Trotter 오류. MPF는 곱 공식에서 발생하는 Trotter 오류를 줄이는 데 도움이 되지만, 곱 공식에 비해 더 높은 샘플링 오류가 발생합니다. 이는 유리하게 작용할 수 있는데, 곱 공식은 샘플링을 늘려 샘플링 오류를 줄일 수 있지만, Trotter 근사로 인한 체계적 오류는 그대로 남아 있기 때문입니다.

플롯에서 관찰할 수 있는 또 다른 흥미로운 동작은, Trotter 단계 수 선택 방법에 관한 가이드에서 설명된 것처럼 $t/k > 1$인 시간부터 $k=1$의 PF 기댓값이 (정확한 값에 대한 좋은 근사가 아닌 것에 더해) 불규칙하게 동작하기 시작한다는 점입니다.

Step 1: Map classical inputs to a quantum problem

이제 단일 시간 $t=1.0$을 고려하고 하나의 QPU를 사용하여 다양한 방법으로 자화(magnetization)의 기댓값을 계산해 보겠습니다. $t$의 특정 선택은 다양한 방법들 간의 차이를 최대화하고 상대적인 효과를 관찰하기 위해 이루어졌습니다. 동적 MPF가 다중 곱 공식 내의 개별 Trotter 공식보다 더 낮은 오류를 가진 관측값을 생성하도록 보장되는 시간 범위를 결정하기 위해 "MPF 검정(MPF test)"을 구현할 수 있습니다 — 참고문헌 [3]의 식 (17)과 주변 텍스트를 참조하세요.

Set up the Trotter circuits

이 시점에서, 우리는 확장 계수 $x$를 찾았으며, 이제 남은 것은 Trotterized 양자 Circuit을 생성하는 것입니다. 다시 한번, qiskit_addon_utils.problem_generators 모듈이 이를 위한 유용한 함수를 제공합니다:

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

total_time = 1.0
mpf_circuits = []
for k in mpf_trotter_steps:
    # Initial Neel state preparation
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2 != 0])

    trotter_circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(order=order, reps=k),
        time=total_time,
    )

    circuit.compose(trotter_circ, inplace=True)

    mpf_circuits.append(pm.run(circuit))

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

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

단일 시간 지점에 대한 기댓값 계산으로 돌아가 보겠습니다. 하드웨어에서 실험을 실행하기 위한 Backend를 선택합니다.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
print(backend)

qubits = list(range(backend.num_qubits))

그런 다음, Transpiler의 레이아웃 단계에서 이상값 큐비트가 포함되지 않도록 커플링 맵에서 해당 큐비트를 제거합니다. 아래에서는 target 객체에 저장된 백엔드 속성 정보를 사용하여, 측정 오류 또는 2큐비트 게이트 오류가 특정 임계값(max_meas_err, max_twoq_err)을 초과하거나 $T_2$ 시간(결맞음 손실을 결정)이 특정 임계값(min_t2) 미만인 큐비트를 제거합니다.

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

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

max_meas_err = 0.012
min_t2 = 40
max_twoq_err = 0.005

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    if target.qubit_properties[q].t2 is not None:
        t2 = target.qubit_properties[q].t2 * 1e6
    else:
        t2 = 0
    if meas_err > max_meas_err or t2 < min_t2:
        # print(q)
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
    twoq_gate_err = target[instruction_2q][q].error
    if twoq_gate_err > max_twoq_err:
        # print(q)
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates
    + ["measure"],  # or whatever new set of gates
    coupling_map=cust_cmap,
)

sorted_components = sorted(
    [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
    reverse=True,
)
print("size of largest component", len(sorted_components[0]))

size of largest component 10

회로를 실행하기에 충분히 큰 큐비트 부분집합을 찾을 수 있도록 max_meas_err, min_t2, max_twoq_err 값을 설정해야 합니다. 이 경우에는 10큐비트 1차원 체인을 찾는 것으로 충분합니다.

cust_cmap.draw()

그런 다음 회로와 옵저버블을 장치의 물리적 큐비트에 매핑할 수 있습니다.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
    optimization_level=3, target=cust_target
)

transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

qubits_layouts = [
    [
        idx
        for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
        if qb._register.name != "ancilla"
    ]
    for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
    transpiler = generate_preset_pass_manager(
        optimization_level=3, backend=backend, initial_layout=layout
    )
    transpiled_circuit = transpiler.run(circuit)
    transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

print(transpiled_circuits[-1].depth(lambda x: x.operation.num_qubits == 2))
print(transpiled_circuits[-1].count_ops())
transpiled_circuits[-1].draw("mpl", idle_wires=False, fold=False)

51
OrderedDict([('sx', 310), ('rz', 232), ('cz', 132), ('x', 19)])

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

Estimator 프리미티브를 사용하면 QPU에서 기댓값 추정치를 구할 수 있습니다. 추가적인 오류 완화 및 억제 기법을 적용하여 최적화된 AQC 회로를 실행합니다.

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
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, 3, 5)
estimator.options.resilience.zne.extrapolator = ("exponential", "linear")

estimator.options.environment.job_tags = ["mpf small"]

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

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

유일한 후처리 단계는 각각의 MPF 계수를 사용하여 서로 다른 Trotter 단계에서 Qiskit Runtime 프리미티브로부터 얻은 기댓값을 결합하는 것입니다. 관측량 $A$에 대해 다음이 성립합니다:

$\langle A \rangle_{\text{mpf}} = \text{Tr} [A \mu(t)] = \sum_{j} x_j \text{Tr} [A \rho_{k_j}] = \sum_{j} x_j \langle A \rangle_j$

먼저, 각 Trotter Circuit에 대해 얻은 개별 기댓값을 추출합니다:

result_exp = job.result()
evs_exp = [res.data.evs for res in result_exp]
evs_std = [res.data.stds for res in result_exp]

print(evs_exp)

[array(-0.06361607), array(-0.23820448), array(-0.50271805)]

다음으로, MPF 계수와 결합하여 MPF의 전체 기댓값을 구합니다. 아래에서는 $x$를 계산한 각각의 방법에 대해 이를 수행합니다.

exact_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_exact.value, evs_std)
        ]
    )
)
print(
    "Exact static MPF expectation value: ",
    evs_exp @ coeffs_exact.value,
    "+-",
    exact_mpf_std,
)
approx_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_approx.value, evs_std)
        ]
    )
)
print(
    "Approximate static MPF expectation value: ",
    evs_exp @ coeffs_approx.value,
    "+-",
    approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(mpf_dynamic_coeffs_list[7], evs_std)
        ]
    )
)
print(
    "Dynamic MPF expectation value: ",
    evs_exp @ mpf_dynamic_coeffs_list[7],
    "+-",
    dynamic_mpf_std,
)

Exact static MPF expectation value:  -0.6329590442738475 +- 0.012798249760406036
Approximate static MPF expectation value:  -0.5690390035339492 +- 0.010459559917168473
Dynamic MPF expectation value:  -0.4655579758795695 +- 0.007639139186720507

마지막으로, 이 소규모 문제에 대해서는 scipy.linalg.expm을 사용하여 정확한 참조값을 다음과 같이 계산할 수 있습니다:

from scipy.linalg import expm
from qiskit.quantum_info import Statevector

exp_H = expm(-1j * total_time * hamiltonian.to_matrix())

initial_state_circuit = QuantumCircuit(L)
initial_state_circuit.x([i for i in range(L) if i % 2 != 0])
initial_state = Statevector(initial_state_circuit).data

time_evolved_state = exp_H @ initial_state

exact_obs = (
    time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
)
print("Exact expectation value ", exact_obs.real)

Exact expectation value  -0.39909900734489434

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs_exp[k],
        yerr=evs_std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )  # , , )

plt.errorbar(
    3,
    evs_exp @ coeffs_exact.value,
    yerr=exact_mpf_std,
    markersize=4,
    marker="o",
    color="purple",
    label="Static MPF",
)

plt.errorbar(
    4,
    evs_exp @ coeffs_approx.value,
    yerr=approx_mpf_std,
    markersize=4,
    marker="o",
    color="orange",
    label="Approximate static MPF",
)

plt.errorbar(
    5,
    evs_exp @ mpf_dynamic_coeffs_list[7],
    yerr=dynamic_mpf_std,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

plt.axhline(
    y=exact_obs.real,
    linestyle="--",
    color="red",
    label="Exact time-evolution",
)

plt.title(
    f"Expectation values for (ZZ,{(L//2-1, L//2)}) 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()

위 예시에서 동적 MPF 방법은 기댓값 측면에서 가장 좋은 성능을 보이며, 가장 높은 Trotter 단계 수만을 사용했을 때보다 개선된 결과를 제공합니다. 다양한 MPF 기법이 항상 가장 높은 Trotter 단계 수 대비 향상된 기댓값을 달성하지는 않더라도(위 플롯의 정확 및 근사 모델의 경우처럼), 이 값들의 표준 편차는 MPF 기법 사용 시 발생하는 분산 증가를 잘 포착합니다. 이는 얻어진 기댓값 주변의 불확실성을 강조하며, 항상 시스템의 정확한 시간 진화에서 기대되는 기댓값을 포함합니다. 반면, 더 낮은 Trotter 단계 수로 계산된 기댓값은 불확실성 범위 내에서 정확한 기댓값을 포착하지 못하여, 잘못된 결과를 자신 있게 반환하게 됩니다.

def relative_error(ev, exact_ev):
    return abs(ev - exact_ev)

relative_error_k = [relative_error(ev, exact_obs.real) for ev in evs_exp]
relative_error_mpf = relative_error(evs_exp @ mpf_coeffs, exact_obs.real)
relative_error_approx_mpf = relative_error(
    evs_exp @ coeffs_approx.value, exact_obs.real
)
relative_error_dynamic_mpf = relative_error(
    evs_exp @ mpf_dynamic_coeffs_list[7], exact_obs.real
)

print("relative error for each trotter steps", relative_error_k)
print("relative error with MPF exact coeffs", relative_error_mpf)
print("relative error with MPF approx coeffs", relative_error_approx_mpf)
print("relative error with MPF dynamic coeffs", relative_error_dynamic_mpf)

relative error for each trotter steps [0.33548293650112293, 0.16089452939226306, 0.10361904247828346]
relative error with MPF exact coeffs 0.2338600369291003
relative error with MPF approx coeffs 0.16993999618905486
relative error with MPF dynamic coeffs 0.06645896853467514

파트 II: 규모 확장

정확한 시뮬레이션이 가능한 범위를 넘어 문제의 규모를 확장해 보겠습니다. 이 섹션에서는 참고문헌 [3]에 제시된 결과 중 일부를 재현하는 데 집중합니다.

Step 1: 고전적 입력을 양자 문제로 매핑하기

해밀토니안

대규모 예시에서는 50개 사이트로 구성된 선상의 XXZ 모델을 사용합니다:

$$\hat{\mathcal{H}}_{XXZ} = \sum_{i=1}^{L-1} J_{i,(i+1)}\left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ 2\cdot Z_i Z_{(i+1)} \right) \, ,$$

여기서 $J_{i,(i+1)}$는 에지 $(i, i+1)$에 대응하는 무작위 계수입니다. 이것이 참고문헌 [3]의 데모에서 사용된 해밀토니안입니다.

L = 50
# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(L, bidirectional=False)
graphviz_draw(coupling_map.graph, method="circo")

import numpy as np
from qiskit.quantum_info import SparsePauliOp, Pauli

# Generate random coefficients for our XXZ Hamiltonian
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,
    )

print(hamiltonian)

SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'XXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'YYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'ZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII'],
              coeffs=[1.        +0.j, 2.09762701+0.j, 2.09762701+0.j, 4.19525402+0.j,
 2.43037873+0.j, 2.43037873+0.j, 4.86075747+0.j, 2.20552675+0.j,
 2.20552675+0.j, 4.4110535 +0.j, 2.08976637+0.j, 2.08976637+0.j,
 4.17953273+0.j, 1.8473096 +0.j, 1.8473096 +0.j, 3.6946192 +0.j,
 2.29178823+0.j, 2.29178823+0.j, 4.58357645+0.j, 1.87517442+0.j,
 1.87517442+0.j, 3.75034885+0.j, 2.783546  +0.j, 2.783546  +0.j,
 5.567092  +0.j, 2.92732552+0.j, 2.92732552+0.j, 5.85465104+0.j,
 1.76688304+0.j, 1.76688304+0.j, 3.53376608+0.j, 2.58345008+0.j,
 2.58345008+0.j, 5.16690015+0.j, 2.05778984+0.j, 2.05778984+0.j,
 4.11557968+0.j, 2.13608912+0.j, 2.13608912+0.j, 4.27217824+0.j,
 2.85119328+0.j, 2.85119328+0.j, 5.70238655+0.j, 1.14207212+0.j,
 1.14207212+0.j, 2.28414423+0.j, 1.1742586 +0.j, 1.1742586 +0.j,
 2.3485172 +0.j, 1.04043679+0.j, 1.04043679+0.j, 2.08087359+0.j,
 2.66523969+0.j, 2.66523969+0.j, 5.33047938+0.j, 2.5563135 +0.j,
 2.5563135 +0.j, 5.112627  +0.j, 2.7400243 +0.j, 2.7400243 +0.j,
 5.48004859+0.j, 2.95723668+0.j, 2.95723668+0.j, 5.91447337+0.j,
 2.59831713+0.j, 2.59831713+0.j, 5.19663426+0.j, 1.92295872+0.j,
 1.92295872+0.j, 3.84591745+0.j, 2.56105835+0.j, 2.56105835+0.j,
 5.12211671+0.j, 1.23654885+0.j, 1.23654885+0.j, 2.4730977 +0.j,
 2.27984204+0.j, 2.27984204+0.j, 4.55968409+0.j, 1.28670657+0.j,
 1.28670657+0.j, 2.57341315+0.j, 2.88933783+0.j, 2.88933783+0.j,
 5.77867567+0.j, 2.04369664+0.j, 2.04369664+0.j, 4.08739329+0.j,
 1.82932388+0.j, 1.82932388+0.j, 3.65864776+0.j, 1.52911122+0.j,
 1.52911122+0.j, 3.05822245+0.j, 2.54846738+0.j, 2.54846738+0.j,
 5.09693476+0.j, 1.91230066+0.j, 1.91230066+0.j, 3.82460133+0.j,
 2.1368679 +0.j, 2.1368679 +0.j, 4.2737358 +0.j, 1.0375796 +0.j,
 1.0375796 +0.j, 2.0751592 +0.j, 2.23527099+0.j, 2.23527099+0.j,
 4.47054199+0.j, 2.22419145+0.j, 2.22419145+0.j, 4.44838289+0.j,
 2.23386799+0.j, 2.23386799+0.j, 4.46773599+0.j, 2.88749616+0.j,
 2.88749616+0.j, 5.77499231+0.j, 2.3636406 +0.j, 2.3636406 +0.j,
 4.7272812 +0.j, 1.7190158 +0.j, 1.7190158 +0.j, 3.4380316 +0.j,
 1.87406391+0.j, 1.87406391+0.j, 3.74812782+0.j, 2.39526239+0.j,
 2.39526239+0.j, 4.79052478+0.j, 1.12045094+0.j, 1.12045094+0.j,
 2.24090189+0.j, 2.33353343+0.j, 2.33353343+0.j, 4.66706686+0.j,
 2.34127574+0.j, 2.34127574+0.j, 4.68255148+0.j, 1.42076512+0.j,
 1.42076512+0.j, 2.84153024+0.j, 1.2578526 +0.j, 1.2578526 +0.j,
 2.51570519+0.j, 1.6308567 +0.j, 1.6308567 +0.j, 3.2617134 +0.j])

관측량으로는 $Z_{24}Z_{25}$를 선택합니다. 참고문헌 [3]의 Fig. 5 하단 패널에서 확인할 수 있습니다.

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

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

Trotter 단계 선택

참고문헌 [3]의 Fig. 4에 제시된 실험에서는 2차(order $2$)의 대칭 Trotter 단계 $k_j = [2, 3, 4]$를 사용합니다. 여기서는 MPF와 더 많은 Trotter 단계 수(이 경우 6)를 사용한 PF가 동일한 Trotter 오차를 갖는 시간 $t=3$의 결과에 집중합니다. 그러나 MPF 기댓값은 더 적은 Trotter 단계 수에 해당하는 Circuit, 즉 더 얕은 Circuit에서 계산됩니다. 실제로, MPF와 더 깊은 Trotter 단계 Circuit이 동일한 Trotter 오차를 갖더라도, MPF Circuit에서 계산된 실험적 기댓값이 이론값에 더 가까울 것으로 기대됩니다. 이는 더 높은 Trotter 단계 PF에 해당하는 Circuit에 비해 하드웨어 노이즈에 덜 노출된 더 얕은 Circuit을 실행하기 때문입니다.

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

LSE 설정

여기서는 이 문제에 대한 정적 MPF 계수를 살펴봅니다.

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(
    f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)
print("L1 norm:", np.linalg.norm(mpf_coeffs, ord=1))

The static coefficients associated with the ansatze are: [ 0.26666667 -2.31428571  3.04761905]
L1 norm: 5.628571428571431

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

[-0.24255546 -0.25744454  1.5       ]
L1 norm of the approximate coefficients: 2.0

동적 계수

# 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": 4,
    },
)

# Create exact time-evolution circuits
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)
]

# Create the time-evolution object
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)

lse = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs = setup_frobenius_problem(lse)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs.value
except Exception as error:
    print(error, "Calculation Failed for time", total_time)
print("")

MPF 분해에서 각 Trotter Circuit 구성

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit import QuantumCircuit

mpf_circuits = []
for k in mpf_trotter_steps:
    # Initial state preparation |1010..>
    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)

MPF와 비슷한 Trotter 오류를 가진 Trotter Circuit 구성

k = 6

# Initial state preparation |1010..>
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=k, order=order),
    time=total_time,
)

comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)

mpf_circuits.append(comp_circuit)

2단계: 양자 하드웨어 실행을 위해 문제 최적화

import copy
from qiskit.transpiler import Target, CouplingMap

target = backend.target
instruction_2q = "cz"

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

max_meas_err = 0.055
min_t2 = 30
max_twoq_err = 0.01

# Remove qubits with bad measurement or t2
cust_cmap_list = copy.deepcopy(cmap_list)
for q in range(target.num_qubits):
    meas_err = target["measure"][(q,)].error
    if target.qubit_properties[q].t2 is not None:
        t2 = target.qubit_properties[q].t2 * 1e6
    else:
        t2 = 0
    if meas_err > max_meas_err or t2 < min_t2:
        # print(q)
        for q_pair in cmap_list:
            if q in q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

# Remove qubits with bad 2q gate or t2
for q in cmap_list:
    twoq_gate_err = target[instruction_2q][q].error
    if twoq_gate_err > max_twoq_err:
        # print(q)
        for q_pair in cmap_list:
            if q == q_pair:
                try:
                    cust_cmap_list.remove(q_pair)
                except ValueError:
                    continue

cust_cmap = CouplingMap(cust_cmap_list)

cust_target = Target.from_configuration(
    basis_gates=backend.configuration().basis_gates
    + ["measure"],  # or whatever new set of gates
    coupling_map=cust_cmap,
)

sorted_components = sorted(
    [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],
    reverse=True,
)
print("size of largest component", len(sorted_components[0]))

size of largest component 73

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

transpiler = generate_preset_pass_manager(
    optimization_level=3, target=cust_target
)

transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

qubits_layouts = [
    [
        idx
        for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()
        if qb._register.name != "ancilla"
    ]
    for circuit in transpiled_circuits
]

transpiled_circuits = []
for circuit, layout in zip(mpf_circuits, qubits_layouts):
    transpiler = generate_preset_pass_manager(
        optimization_level=3, backend=backend, initial_layout=layout
    )
    transpiled_circuit = transpiler.run(circuit)
    transpiled_circuits.append(transpiled_circuit)

# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

3단계: Qiskit 프리미티브를 사용하여 실행

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Set simple error suppression/mitigation options
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 = ["mpf large"]

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

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.08034071), array(-0.00605026), array(-0.15345759), array(-0.18127293)]
[array(0.04482517), array(0.03438413), array(0.21540776), array(0.21520829)]

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.47510243192011536 +- 0.6613940032465087
Approximate static MPF expectation value:  -0.20914170384216998 +- 0.32341567460419135
Dynamic MPF expectation value:  -0.07994951978722761 +- 0.07423091963310202

sym = {2: "^", 3: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
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="6 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"Expectation values for (ZZ,{(L//2-1, L//2)}) 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()

하드웨어에서 Circuit을 실행할 때는 하드웨어 노이즈로 인해 정확한 기댓값을 얻는 데 추가적인 어려움이 생길 수 있습니다. 이는 MPF 형식 체계에서 고려되지 않으며, MPF 솔루션에 불리하게 작용할 수 있습니다. 예를 들어, 플롯에서 동적 계수가 근사 정적 계수에 비해 기댓값의 더 나은 추정치를 제공하지 못한 이유가 여기에 있을 수 있습니다. 즉, 근사 Circuit을 시뮬레이션하는 근사 에볼버가 하드웨어 노이즈 환경에서 근사 Circuit을 실행하여 얻은 결과를 정확하게 반영하지 못할 수 있습니다. 이러한 이유로, 각 곱 공식에 대해 이상적인 값에 최대한 가까운 결과를 얻으려면 다양한 오류 경감 기법을 함께 사용하는 것이 권장됩니다. 이를 통해 MPF 접근법의 일관된 이점을 확인할 수 있을 것입니다.

전반적으로, 근사 정적 계수는 노이즈가 없는 환경에서 동일한 Trotter 오류를 가지면서도 Trotter 스텝 수가 더 많은 곱 공식보다 여전히 더 정확한 해를 제공합니다.

또한 참고 문헌 [3]의 실험을 재현하는 예제에서, 시간점 $t=3$은 $k=2$인 PF가 제대로 동작할 것으로 기대되는 한계인 $t/k>1$을 초과한다는 점에 유의하는 것이 중요합니다. 이는 이 가이드에서 논의된 내용입니다.

참고 문헌

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