多積公式によるトロッター誤差の低減

使用量の目安：Heron r2プロセッサで約4分（注意：これは推定値です。実際の実行時間は異なる場合があります。）

学習の成果

• 複数の浅い回路からの期待値を組み合わせることで、多積公式（MPF）がハミルトニアン・シミュレーションにおけるトロッター誤差をどのように削減するか

• 標準的な積公式に対して MPF が有利な場合と、適切なツールではない場合

• qiskit_addon_mpf パッケージを使用して、静的および動的な MPF 係数を計算する方法

• トランスパイル、エラー緩和（エラー軽減）、後処理を含め、MPF ワークフローを IBM Quantum® ハードウェア上でエンドツーエンドに実行する方法

前提条件

• ハミルトニアン・シミュレーション回路のコンパイル手法 — Qiskit におけるトロッター（積公式）回路を紹介します。

• Qiskit における積公式、特に SuzukiTrotter および LieTrotter 合成クラス。

• Qiskit プリミティブと Estimator インターフェース。

背景

多積公式とは？

量子コンピュータ上で量子系をシミュレートする際、中心となる課題は、ハミルトニアン $H$ に対する時間発展演算子 $e^{-iHt}$ を近似することです。標準的なアプローチでは、トロッター・鈴木分解としても知られる積公式（PF）を使用します。これらは $H = \sum_{a=1}^d F_a$ を、個々のユニタリ $e^{-iF_a t}$ が効率的に実装できる項に分解し、それらの単純なユニタリの順序付き積として完全な発展を近似します。

1次積公式（Lie-Trotter）は以下のようになります：

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

これは2次の誤差を伴います：$S_1(t) = e^{-iHt} + \mathcal{O}(t^2)$。より高次の対称公式 $S_{2\chi}(t)$（$\chi$ は対称積公式の次数を表す。参考文献 [1] を参照）は、$e^{-iHt} + \mathcal{O}(t^{2\chi+1})$ としてより速く収束しますが、ステップあたりの回路がより深くなるという代償があります。

固定された次数 $\chi$ で誤差を減らすために、通常は全発展時間 $t$ を $k$ 個のより小さなトロッター・ステップに分割します。各ステップは積公式によって $e^{-iHt/k}$ を近似し、それらのステップを連結します：

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

$2\chi$次の対称公式の場合、残留トロッター誤差は $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$ としてスケールします。したがって、$k$ を増やすとトロッター誤差は急速に抑制されます — しかし同時に回路の深さも線形に増加し、ノイズのあるハードウェアでは蓄積されるゲートノイズも増えることを意味します。**トロッター誤差（大きな $k$ を好む）とハードウェアノイズ（小さな $k$ を好む）**の間のこの緊張関係こそ、多積公式がまさに解決しようとしているものです。MPF は、固定された次数 $\chi$ における異なる $k$ の選択からの結果を組み合わせるものであり、基礎となる積公式の次数を変えるものではないことに注意してください。

多積公式（MPF） [1] は、それぞれ異なるトロッター・ステップ数 $k_1, k_2, \ldots, k_r$（$r$ 個のステップ数の集合）を使用する、複数のより浅いトロッター回路から得られた期待値の重み付き線形結合を構築します：

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

ここで、$\langle A \rangle_{k_j}(t)$ は、$k_j$ ステップのトロッター回路から推定された時刻 $t$ におけるオブザーバブル $A$ の期待値であり、係数 $\{x_j\}_{j=1}^r$ は、組み合わせにおける主要なトロッター誤差項が打ち消し合うように選ばれます。この表現には ステップ4 で再び触れ、トロッター結果を組み合わせるために明示的に評価します。実用上の重要な点は、MPF の中で最も深い回路が必要とするのは $k_{\max}$ ステップのみであり、これは同じ有効トロッター誤差に直接到達するために必要となる単一の $k$ よりもはるかに小さいということです。より浅い回路であることにより、MPF アプローチはノイズのあるハードウェアにより適したものになります。

係数はどのように決定されるか？

MPF 係数には2つのファミリーがあります：

静的係数は、ハミルトニアン、初期状態、発展時間に依存しません。これらは、主要なトロッター誤差項の打ち消しを強制する線形系 $Ax = b$ を解くことで求められます。$2\chi$次の対称積公式で使用されるトロッター・ステップの集合 $\{k_j\}_{j=1}^r$ について、トロッター誤差を $k_j$ の逆べきで展開すると、次の形の制約方程式が得られます：

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

ここで、整数指数 $\{\eta_n\}$ は、選択した積公式に対する連続するトロッター誤差項の次数です。対称的な $2\chi$次の PF の場合、$\left[S_{2\chi}(t/k)\right]^k$ における主要な誤差は $1/k^{2\chi}$ としてスケールし、以降の補正は $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ となります — したがって指数は $\eta_n = 2\chi + 2n$ です。非対称な PF の場合、奇数べきと偶数べきの両方が寄与し、$\eta_n = 2\chi + n$ となります。完全な導出については参考文献 [1] を参照してください。上記の系の最初の方程式は不偏性を保証し（MPF は $k_j \to \infty$ の極限で厳密な期待値を再現します）、残りの $r-1$ 個の方程式は最初の $r-1$ 個のトロッター誤差項を順次打ち消します。結果として得られる $L_1$ ノルム $\|x\|_1$ が大きすぎる場合（サンプリングノイズを増幅させる）、代わりに $\|x\|_1$ に上限を設けながら $\|Ax - b\|$ を最小化する近似最適化を解くことができます。

動的係数 [2]、[3] は、さらにハミルトニアン、初期状態、発展時間 $t$ に依存します。これらは、真の時間発展状態と MPF 近似との間のフロベニウス・ノルム距離を最小化します：

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

ここで、$M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ は、異なるステップ数 $k_i, k_j$ に対するトロッター発展状態間の重なりのグラム行列であり、$L_i(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$ は（近似的な）厳密状態との重なりを測定します。このチュートリアルでは、これらの量はテンソルネットワーク法、具体的には qiskit_addon_mpf における TeNPy ベースのバックエンドを使用して効率的に計算されます。

MPF をいつ使用するか

MPFが特に有益なのは、次のような場合です。

• 回路の深さがボトルネックとなる場合。 ハードウェアノイズによって実行可能な深さが制限されている場合、MPFを使用することで、より浅い回路からより高い実効トロッター精度を達成できます。

• フル状態準備ではなく、正確な期待値が必要な場合。 MPFは期待値のレベルで動作します — 量子状態ではなく、古典的な数値を組み合わせます。そのため、Estimatorプリミティブを使用したオブザーバブル推定に理想的です。

• 中程度の数のトロッターステップ数を組み合わせる場合。 通常、$r = 3$–$5$個の異なるステップ数$k_j$を組み合わせれば、$\|x\|_1$を管理可能な範囲に保ちながら、いくつかの主要なトロッター誤差項を打ち消すのに十分です。

MPFが役立たない場合

• 非常に短い時間発展の場合。 $t$が十分に小さく、単一の低次トロッター公式ですでに十分な精度がある場合、複数の回路を実行するオーバーヘッドは不要です。

• 状態準備タスクの場合。 MPFは補正された期待値を生成しますが、補正された量子状態は生成しません。実際の時間発展した状態が必要な場合（たとえば、別の量子サブルーチンへの入力として）、MPFは適用できません。

• 収束領域に違反するトロッターステップ数の場合。 静的係数の導出では、各$\left[S_{2\chi}(t/k_j)\right]^{k_j}$を$t/k_j$の級数として展開します。この展開は$t/k_{\min} \lesssim 1$の場合にのみうまく収束します。$k_{\min}$が与えられた$t$に対して小さすぎる場合、最も浅い回路は摂動領域から大きく外れ、MPFが打ち消せない高次の誤差項が大きくなり、打ち消しには大きな係数が必要になる場合があります。$L_1$ノルム$\|x\|_1$は実用的な診断指標です。$\|x\|_1 \gg 1$の場合、サンプリングオーバーヘッド$\propto \|x\|_1^2$がトロッター誤差の削減効果を上回る可能性があります。詳細はトロッターステップの選択に関するガイドを参照してください。

このチュートリアルの内容

このチュートリアルでは、2つの段階でエンドツーエンドの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）

• TeNPyバックエンドを備えたMPF Qiskitアドオン（pip install "qiskit-addon-mpf[tenpy]"）

• Qiskit アドオン・ユーティリティ（pip install qiskit-addon-utils）

• SciPy (pip install scipy)

セットアップ

以下では、このチュートリアル全体で使用するすべてのパッケージのインポートを1つのセルにまとめています。また、隣接するrxx回転とryy回転を1つのXXPlusYYGateに融合するCollectAndCollapseトランスパイラ・パスも定義します。このパスは、ステップ1での回路構築時（ゲート数を低く抑えるため）と、ステップ4で動的MPF用に層構造を抽出する際（TeNPyは融合されていない回転のペアではなく、2量子ビットゲートを想定しているため）に間接的に適用されます。

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_addon_mpf": "qiskit-addon-mpf", "qiskit_addon_utils": "qiskit-addon-utils", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy", "tenpy": "physics-tenpy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

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: 古典的な入力を量子問題にマッピングする

まず、初期状態としてネール状態$\vert 0101\ldots01 \rangle$を用いた、直線上の10量子ビットのハイゼンベルグ模型から始めます。ハミルトニアンは次のとおりです。

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

ここで、$J$は最近接結合強度です。鎖の中央にある1組の量子ビットでZZ相関子$Z_{L/2-1} Z_{L/2}$を測定し、2次の積公式を用いてトロッターステップ$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)

トロッター回路の構築

各時間点と各トロッターステップ数について、近似的なトロッター時間発展を実装する回路を作成します。セットアップセクションで定義した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)

ステップ 2: 量子ハードウェア実行のための問題の最適化

小規模な例では、Aerシミュレーターを対象とします。回路が実行可能になるまでに、2つの変換が行われます。

1. ハミルトニアン・シミュレーションレベルでのゲート収集。 セットアップセルで、隣接するrxx回転とryy回転を1つのXXPlusYYGateに融合するCollectAndCollapseパスを構築しました。このパスは、ステップ1でトロッター回路を構築する際（pm.run(...)呼び出し）にすでに適用しています。これにより、2量子ビットゲート数が削減されるとともに、後の動的係数計算のためのテンソルネットワークシミュレーションに、より適した構造が生成されます。

2. シミュレーターのISAへの変換。 以下では、Qiskitのプリセット・パスマネージャーをoptimization_level=3で実行し、各トロッター回路をシミュレーターの命令セット・アーキテクチャ（ISA）に変換します。

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

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

ステップ 3: Qiskitプリミティブを使用した実行

小規模な例では、ISAに変換されたトロッター回路を、Aerを利用したEstimatorV2プリミティブで実行します。これにより、各$(k_j, t)$の組み合わせに対するノイズのない参照値が得られます — これらは、ステップ4でMPFが組み合わせる$\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])

ステップ 4: 後処理と所望の古典形式での結果の返却

ステップ4では、実際にMPFを構築します。係数$x_j$はここで計算されます（動的変種では、この計算は負荷が高くなる場合があります）が、概念的にはこれらは、ステップ3の量子測定結果を単一の補正済み期待値に組み合わせるための古典的なレシピです — そのため、係数と組み合わせのワークフロー全体を後処理として扱います。

MPFが真のダイナミクスをどれだけ正確に追跡できているかを評価するために、まずハミルトニアンを直接指数化して、厳密な時間発展した期待値を計算します。これは$L = 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は、時間発展の時間、ハミルトニアン、初期状態に依存しない係数$x_j$を使用します。背景セクションで説明した線形システム$Ax = b$を設定し、係数を解きます。行列$A$は、トロッターステップ数$k_j$、プロダクト公式の次数$\chi$、および公式が対称かどうか（指数$\eta_n$を決定します）によって決まります。

この小規模な例では、$k_j = [1, 2, 4]$と非対称な次数$2\chi=2$の鈴木・トロッター公式を使用します（つまり$\chi=1$、$\eta_n = 2 + n$となり、$\eta_0 = 2,\, \eta_1 = 3$）。システムは次のようになります。

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

1行目は不偏性を保証します（$\sum_j x_j = 1$）。2行目と3行目は、それぞれ主要な$1/k^2$項と次の$1/k^3$トロッター誤差項を打ち消します。

LSEの設定

上記の行列$A$と右辺のベクトル$b$を組み立てるために、qiskit_addon_mpf.staticのsetup_static_lseを使用します。行列$A$は$k_j$だけでなく、選択したプロダクト公式にも依存します — 特にその次数$\chi$と対称かどうかです。symmetricフラグは、指数パターン$\eta_n$を制御します（対称な公式は偶数次のトロッター誤差項のみを生成します。参考文献[1]を参照）。参考文献[2]に示されているように、基礎となるPFが対称であっても、symmetric=Trueに設定することは厳密には必要ではないことに注意してください — 非対称のLSEは依然として有効です（余分な不要な制約を課すことになります）。

この例では、ステップ1ですでにorder = 2とsymmetric = Falseを設定しています。

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

構築された行列$A$とベクトル$b$を確認し、上記のシステムと一致していることを確認します。

lse.A

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

lse.b

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

LSEが得られたら、lse.solve()を介して静的係数$x_j$を解きます（これは直接的な$x = A^{-1}b$の解です）。

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

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

厳密モデルを使用した$x$の最適化

$x = A^{-1}b$を計算する代わりに、setup_exact_modelを使用して、LSEを制約として使用し、その最適解が$x$を与える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

近似モデルを使用した$x$の最適化

選択した$k_j$値の集合に対して$L_1$ノルムが高すぎると判断される場合があります。その場合、異なる$k_j$値の集合を選択できないのであれば、$\|Ax - b\|$を最小化しながら$L_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は、ハミルトニアンや状態に依存しない方法でトロッター誤差項を打ち消すため、特定のハミルトニアンと初期状態に対して必ずしも最小の近似誤差を生成するわけではありません。動的MPF（参考文献[2]、[3]）は代わりに、各時刻$t$においてフロベニウスノルム距離$\|\rho(t) - \mu^D(t)\|_F^2$を最小化する時間依存係数$x_i(t)$を求めます。背景セクションで示したように、これにはトロッター発展した状態間のオーバーラップ行列$M_{ij}(t)$と、厳密な状態とのオーバーラップ$L_i(t)$が必要です — どちらもqiskit_addon_mpf内のテンソルネットワーク（TeNPy）バックエンドを使用して推定します。

動的LSEを設定するには、3つの要素が必要です。

1. アドオンが各$k_j$に対して実行し、$\rho_{k_j}(t)$をMPS/MPOとして生成する近似発展器ファクトリ。これは、2次トロッター回路の層構造から構築します（slice_by_depthごとに1層）。TeNPyの切り捨てパラメーターを持つLayerwiseEvolverとしてラップします。

2. 高精度な参照$\rho(t)$を生成する厳密発展因子。厳密な発展の代理として、小さな時間ステップの4次鈴木・トロッター回路（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の切り捨てパラメーターは精度に影響を与える可能性があるため、値の範囲を検討することが重要です。

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を定義し、層状トロッターモデルで使用される格子に一致するMPSとしてネール初期状態を準備します。

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がTeNPyを介して関連するオーバーラップ行列を構築し、setup_frobenius_problemがフロベニウスノルムコストを最小化するcvxpy.Problemを返します。ソルバーはその時刻に合わせた係数$x_j(t)$を返し、それをmpf_dynamic_coeffs_listに収集します。特定の$t$でソルバーが失敗した場合は、ゼロ係数にフォールバックしてループを続行します。

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

トロッター期待値とMPF係数の組み合わせ

次に、各係数の組み合わせ（静的厳密、静的近似、動的）について$\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()

上記のプロットは、トロッター誤差とサンプリング誤差の相互作用を示しています。

• トロッター誤差。 個々の積公式（灰色のマーカー）は、時間が経つにつれて厳密な曲線からますます逸脱していきます。$k=1$回路は最も浅い一方で逸脱も最も大きく、これはすでに$t/k \gtrsim 1$の領域にあって主要な$1/k^{2}$誤差項が大きいためです。MPFの組み合わせ（色付きマーカー）は、これらの主要なトロッター誤差項のいくつかを打ち消すため、単一の$k_j$回路よりもはるかに厳密な曲線に近づきます。残るギャップは、MPFが打ち消さない高次のトロッター項を反映しています。次数$2$、$r=3$の静的MPFは最初の2つの誤差次数しか除去せず、大きな$t/k_{\min}$では、打ち消されなかったテール項がやがて支配的になります — そのため、MPFは非常に浅い回路が任意の時刻で正確であり続けることを保証するものではありません。

• サンプリング誤差。 MPF曲線でエラーバーが広くなるのは、線形結合の直接的な帰結です。独立した回路ごとの標準誤差$\sigma_{k_j}$を伝播させると、合計分散は$\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$になります。したがって、$\|x\|_2$（実際には制御対象である$\|x\|_1$）が大きいほど、目標とする不確かさに到達するために必要なショット数が多くなります。これが背景セクションで説明した近似ソルバーの選択肢の背後にあるトレードオフです。私たちは、このオーバーヘッドを管理可能な範囲に保つために$\|x\|_1$を制限します。重要なのは、トロッター誤差とは異なり、サンプリング誤差は$1/\sqrt{N_{\text{shots}}}$で縮小するため、ショット数を増やせば常に減らせるという点です。

以下の大規模なハードウェアの例では、ハードウェアノイズが各$\langle A \rangle_{k_j}$に追加の誤差源として加わり、同様にMPF係数によって増幅されます。そのセクションでは、エラー緩和がMPFとどのように相互作用するかを見ていきます。

大規模なハードウェアの例

このセクションでは、厳密にシミュレーションできる範囲を超えて問題を拡大します。参考文献[3]に示された結果のいくつかを、時刻$t = 3$の50量子ビットXXZ鎖を用いて再現します。小規模な例と同じ4ステップのワークフローに従いますが、今度はエラー緩和を用いて実際の量子ハードウェアを対象とします。テンプレートと同様に、各ステップはコード内にインラインで示されており、中間出力を確認する価値がある場合、1つのステップが複数のセルにまたがることがあります。 このマッピングは、小規模な例を反映しています。ハミルトニアンを定義し、トロッターパラメーターを選択し、MPF係数（静的および動的）を計算し、回路を構築します。主な違いは以下のとおりです。

• $\mathcal{U}(0.5, 1.5)$からランダムに抽出された結合を持つ、50サイト上のXXZハミルトニアン（参考文献[3]）。

• $k_j = [3, 4, 6]$を持つ対称な2次トロッター公式（つまり$\chi=1$、symmetric=True）。

• 単一の固定された時間発展の時刻$t = 3$。$k_{\min}=3$では$t/k_{\min}=1$となり、浅い構成要素がトロッター収束領域内に留まり、MPFが依存する主要誤差モデルが有効であることが保証されます。

• ベースラインとして使用される、追加の単一回路比較実行（$k = 10$トロッターステップ）。$k = 10$を選んだのは、ハードウェア上での2量子ビット深度が、最も深いMPF構成要素（$k_{\max}=6$）に複数のMPF回路を実行するオーバーヘッドを加えたものよりも深いためです。これは、MPFの組み合わせが単一回路のベースラインを上回ることが期待される領域である、ノイズが制限要因となるほど十分に深いものです。これは、MPFの実効トロッター誤差を目標とする回路（そのためにはさらに多くのステップが必要になります）ではなく、MPFの組み合わせと比較する「単一の深い回路」の比較です。

ここではまだステップ1（マッピングと回路構築）の段階ですが、このセルで静的係数とともに動的係数も事前に計算していることに注意してください。動的係数は$H$と$t$に依存しますが、量子測定には依存しないため、ステップ4より前であればいつでも計算できます。MPF固有のセットアップを一箇所にまとめるために、ここで計算を行います。

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

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

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

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

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

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

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

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

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

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

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

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

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

mps_initial_state = MPS_neel_state(models[0].lat)

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

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

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

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

次に、選択したバックエンドに合わせて回路を最適化します。Qiskitのプリセット・パスマネージャーをoptimization_level=3で使用します。これにより、良好な物理量子ビットの組が自動的に選択され、各回路がデバイスのトポロジーにルーティングされます。

# -------------------------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外挿が信頼できる範囲を超えてしまうためです。

4つの回路すべて（$k_j = [3, 4, 6]$の3つのMPF構成要素と$k = 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係数の集合と組み合わせます。$\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}$、伝播された分散は$\sigma^2 = \sum_j x_j^2 \sigma_{k_j}^2$です。

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

print(evs)
print(std)

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

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

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

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

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

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

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

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

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

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

上記のハードウェア結果について、いくつかの観察事項があります。

• ハードウェア上では、より深くすることは無料ではありません。 単一回路のベースラインが直接それを物語っています。$k = 6$回路は本質的に厳密です（参照値$-0.244$に対して$-0.256$）が、より深い$k = 10$のベースラインはより良いのではなく悪くなっています（$-0.061$で、約$0.18$の誤差）。トロッター誤差がすでに小さい場合、ステップを追加することは主に回路を深くし、より多くのゲートノイズとデコヒーレンスを蓄積するだけです。これはまさに、MPFが構築された領域です。浅い構成要素のみを使用して、深い回路の精度を達成することです。

• ノルムの小さいMPFは、深い単一回路を上回ります。 近似静的MPF（$\|x\|_1 \approx 2$に制限）は$-0.259$に到達し、参照値からの誤差は約$0.015$と、$k = 10$のベースラインよりもはるかに近い値です。動的MPF（$-0.127$）も、そのベースラインを楽に上回ります。どちらも浅い$k_j = [3, 4, 6]$回路のみを組み合わせているにもかかわらず、深い単一回路では得られなかった答えを回復しています。

• 係数ノルムは数学的最適性よりも重要です。 厳密静的MPFは$\|x\|_1 = 4.66$であり、すべての推定量の中で最悪です（$-0.567$で、誤差は$0.3$を超えます）。大きな係数ノルムは、各$\langle A \rangle_{k_j}$に残る残留ゲートノイズ、デコヒーレンス、ZNE誤差をほぼ同じ係数で増幅し、それによって得られるトロッター誤差の打ち消し効果を圧倒してしまいます。ノルムを制限すること（近似静的ソルバー、$\|x\|_1 \approx 2$）で、この圧倒的な影響が取り除かれ、最良の推定値が得られます — ただし、その係数はもはや主要なトロッター誤差を厳密には打ち消しません。

• 個々の浅い回路も依然として競争力を持つ場合があります。 単独の$k = 6$構成要素（$-0.256$）自体も、ここでは本質的に厳密です — この実行では、近似静的MPFよりもわずかに近いくらいです。問題は、「収束しているがまだノイズに制限されていない」という最適な領域に、どの単一の$k$が位置するかを事前に知ることができないという点であり、トロッター収束を保証するために単に深くする（$k = 10$）という一見安全に見える選択こそが、まさに失敗する選択なのです。MPFは、正しい深さを推測する必要のない、浅い回路の原理に基づいた組み合わせを提供します。

実用的な教訓は、ハードウェア上ではMPFを各個別の$\langle A \rangle_{k_j}$に対する強力なエラー緩和と組み合わせるべきであり、係数の$L_1$ノルムは控えめに保つべきであり（近似ソルバーまたは動的MPFを使用）、トロッターステップ$k_j$は$t/k_{\min} \lesssim 1$となるように選択すべきだということです — ここでは、$t = 3$における$k_{\min} = 3$により$t/k_{\min} = 1$となり、構成要素が、静的MPFが依存する主要誤差モデルが有効な収束領域内に留まっています。これらの選択により、ここでのノルムの小さいMPFは収束した単一回路に匹敵する一方で、単純に「もっと深くする」という素朴なベースラインはそうならず、参考文献[3]に示された深さと精度のトレードオフの利点が再現されています。また、個々の実行にはノイズがあることにも注意してください — 同じジョブの別の送信（または別のバックエンド）では、正確な順序が変わる可能性があります。堅牢な傾向としては、小さな$\|x\|_1$のMPFはうまく機能し、大きな$\|x\|_1$の厳密静的MPFはハードウェアノイズによって増幅され、深すぎる単一回路はノイズに制限されるという点です。

次のステップ

推奨事項

この作業に興味を持たれた方は、以下の資料にも興味があるかもしれません。

• MPFのためのトロッターステップの選び方 — 不安定性を回避するための$k_j$値の選択に関する実用的なガイダンス
• 近似モデルの使用方法 — 近似静的MPFの$L_1$ノルム制約とソルバーオプションの調整
• qiskit-addon-mpf APIリファレンス — 静的、動的、バックエンドの各モジュールの完全なドキュメント

参考文献

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