マルチプロダクト公式によるトロッター誤差の低減

マルチプロダクト公式によるトロッター誤差の低減#

観測量推定においてマルチプロダクト公式を用いてトロッター誤差を低減する、あるいは一定のトロッター誤差のもとでより浅い深さの時間発展を実装します。

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

学習成果#

このチュートリアルを完了すると、次の内容が理解できるようになります。

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

  • 標準的なプロダクト公式に対して MPF が有利な場面と、MPF が適切な手法ではない場面

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

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

前提知識#

このチュートリアルに進む前に、次のトピックについて理解しておくことをお勧めします。

背景#

マルチプロダクト公式とは#

量子コンピューター上で量子系をシミュレートする際の中心的な課題は、ハミルトニアン \(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})\) としてより速く収束しますが、その代わり 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] を参照してください。上記の方程式系の 1 番目の式は不偏性（\(k_j \to \infty\) の極限で MPF が厳密な期待値を再現すること）を保証し、残りの \(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\) が十分小さく、低次のプロダクト公式 1 つでも既に十分な精度が得られるときは、複数の回路を実行するオーバーヘッドは不要です。

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

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

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

このチュートリアルでは、MPF のエンドツーエンドのワークフローを 2 段階に分けて見ていきます。まず、小規模なシミュレーターでの例（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）

セットアップ#

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

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

# TeNPy の `unit_cell_width` に関する将来 API の警告を抑制する。デフォルト値
# (`unit_cell_width=len(sites)`) は Chain 格子に対して正しく、ここで
# `CouplingMap.from_line(...)` が生成するのはまさにそれなので、この警告は情報提供にすぎない。
warnings.filterwarnings(
    "ignore",
    message=r".*unit_cell_width.*",
    category=UserWarning,
)


# --- 補助: XX と YY の回転を 1 つのゲートにまとめる ---
def filter_function(node):
    return node.op.name in {"rxx", "ryy"}


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


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


collapse_function = partial(
    collapse_to_operation,
    collapse_function=collapse_to_xx_plus_yy,
)

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

小規模なシミュレーターでの例#

ステップ 1: 古典的な入力を量子問題へマッピングする#

まず、直線上の 10 量子ビットのハイゼンベルクモデルから始めます。初期状態にはネール状態 \(\vert 0101\ldots01 \rangle\) を使います。ハミルトニアンは次のとおりです。

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

# カップリングマップとハミルトニアンを生成する
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])

# 観測量: 中央の 2 量子ビットに対する ZZ
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 のパラメーター
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 の回転をまとめて 1 つの XX+YY ゲートにし、後段でのテンソルネットワークシミュレーションをより効率的に行えるように準備します。

# 初期状態としてネール状態を準備する
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)

<Image src="/docs/images/tutorials/multi-product-formula/extracted-outputs/c7ee61e7-0.avif" alt="Output of the previous code cell" />

ステップ 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 の量子測定結果を 1 つの補正された期待値へ組み合わせるための古典的な手順書です。そのため、係数の計算と組み合わせのワークフロー全体を後処理として扱います。

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\) を決めます）によって定まります。

今回の小規模な例では、非対称な \(2\chi=2\) 次の鈴木・トロッター公式とともに \(k_j = [1, 2, 4]\) を使います（したがって \(\chi=1\)、\(\eta_n = 2 + n\) であり、\(\eta_0 = 2,\, \eta_1 = 3\) となります）。方程式系は次のようになります。

\[\begin{split} 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}. \end{split}\]

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

LSE を設定する#

qiskit_addon_mpf.static の setup_static_lse を使い、上で述べた行列 \(A\) と右辺ベクトル \(b\) を組み立てます。行列 \(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 を制約として用いる cvxpy.Problem のインスタンスを構成し、その最適解として \(x\) を得ることもできます。

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

[ 0.04761905 -0.57142857  1.52380952]

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

L1 norm of the exact coefficients: 2.1428571428556378
近似モデルを使って \(x\) を最適化する#

選んだ \(k_j\) の集合に対する \(L_1\) ノルムが高すぎると判断される場合もあります。その場合で、かつ別の \(k_j\) の集合を選べないときは、\(L_1\) ノルムを指定した閾値に制約しつつ \(\|Ax - b\|\) を最小化する近似解を使えます。近似モデルの使い方 のガイドをご覧ください。

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

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

動的 MPF 係数#

静的 MPF はハミルトニアンや状態に依存しない形でトロッター誤差項を打ち消すため、与えられたハミルトニアンと初期状態に対して必ずしも最小の近似誤差をもたらすわけではありません。動的 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。

下のセルでは、近似発展演算子のファクトリーを構成します。

# 近似的な時間発展回路を作成する
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

# 回路内の層を見つける
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": 2,
    },
)
テンソルネットワークシミュレーションの詳細を決める `LayerwiseEvolver` のオプションは、不適切な最適化問題を設定してしまわないよう慎重に選ぶ必要があります。

厳密な時間発展状態は、小さな時間刻み dt=0.1 を用いた 4 次の鈴木・トロッター公式で近似します。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"}
# 各トロッターステップについて、全時刻の期待値を取得する
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",
    )

# 厳密係数による静的 MPF について、全時刻の期待値を取得する
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",
)


# 近似係数による静的 MPF について、全時刻の期待値を取得する
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",
)


# 動的 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",
)


# 厳密な期待値
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()

<Image src="/docs/images/tutorials/multi-product-formula/extracted-outputs/35042576-0.avif" alt="Output of the previous code cell" />

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

  • トロッター誤差。 個々のプロダクト公式（灰色のマーカー）は、時間が経つにつれて厳密な曲線からますます外れていきます。\(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 係数（静的・動的）を計算し、回路を構築します。主な違いは次の点です。

  • 50 サイト上の XXZ ハミルトニアン。結合は \(\mathcal{U}(0.5, 1.5)\) からランダムに引きます（文献 [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 の実効的なトロッター誤差を狙った回路ではありません（それにはもっと多くのステップが必要になります）。

なお、ここはまだステップ 1（マッピングと回路の構成）ですが、このセルでは静的係数とあわせて動的係数も事前に計算しています。動的係数は \(H\) と \(t\) に依存しますが量子測定には依存しないため、ステップ 4 より前の任意の時点で計算できます。MPF 固有の設定を 1 か所にまとめるため、ここで計算しておきます。

# -------------------------ステップ 1-------------------------
L = 50
coupling_map = CouplingMap.from_line(L, bidirectional=False)

# ランダムな結合をもつ XXZ ハミルトニアン（文献 [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

# 静的係数
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)}")

# -------------------------動的係数-------------------------
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")

# -------------------------ステップ 1（続き）: 回路を構築する-------------------------
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)

# k=10 トロッターステップによる「単一の深い回路」とのベースライン比較実行。
# その 2 量子ビット深さは、MPF の最も深い構成要素（k_max=6）に複数回路を実行する
# オーバーヘッドを加えたものより深く、MPF が上回ると期待されるノイズ律速の領域に
# 入る。これは MPF の実効的なトロッター誤差を狙ったものではない
# （それにはもっと多くのステップが必要になる）。
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 で使用します。これにより、適切な物理量子ビットの集合が自動的に選ばれ、各回路がデバイスのトポロジー上にルーティングされます。

# -------------------------ステップ 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]\) の MPF 構成要素 3 つと \(k = 10\) のベースライン）を、1 つの Estimator ジョブとして投入します。

# -------------------------ステップ 3-------------------------
estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# 誤差抑制・誤差低減
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\) として伝播させます。

# -------------------------ステップ 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()

<Image src="/docs/images/tutorials/multi-product-formula/extracted-outputs/64360d85-0.avif" alt="Output of the previous code cell" />

上のハードウェア結果について、いくつか観察できることがあります。

  • ハードウェアでは深くすることに代償が伴う。 単一回路のベースラインがそれを直接物語っています。\(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 がハードウェアノイズによって増幅されること、そして深すぎる単一回路がノイズ律速になることです。

次のステップ#

この内容に興味を持たれた方には、次の資料もお勧めします。

参考文献#

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

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

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

© IBM Corp., 2017-2026