{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "454a9dfd",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"トロッター誤差を低減する多製品式\"\n",
        "description: \"観測可能推定において多製品式を使用し、トロッター誤差を低減するか、あるいはより浅い深さにおいてトロッター誤差を一定に保ちながら時間発展を実装する。\"\n",
        "---\n",
        "\n",
        "{/* cspell:ignore ncol circo Layerwise markersize unbiasedness infty ndash lesssim propto tenpy unfused Néel correlator Neel exponentiating gtrsim */}\n",
        "\n",
        "<span id=\"multi-product-formulas-to-reduce-trotter-error\" />\n",
        "\n",
        "# トロッター誤差を低減する多製品式\n",
        "\n",
        "*推定所要時間：Heron r2 プロセッサで 4 分（注：これはあくまで推定値です。 （実行時間は状況によって異なる場合があります。）*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c4d0b2f2",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## 学習成果\n",
        "\n",
        "このチュートリアルを修了すると、以下の内容を理解できるようになります：\n",
        "\n",
        "* マルチプロダクト式（MPF）が、複数の浅い回路からの期待値を組み合わせることで、ハミルトニアンシミュレーションにおけるトロッター誤差をどのように低減するか\n",
        "* MPFが標準的な製品処方よりも優れている場合と、適切な手段ではない場合\n",
        "* この `qiskit_addon_mpf` パッケージを使用して、静的および動的なMPF係数を計算する方法\n",
        "* IBM Quantum® ハードウェア上で、トランスパイル、エラーの軽減、後処理を含むMPFワークフローをエンドツーエンドで実行する方法\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f5dfd316",
      "metadata": {},
      "source": [
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## 前提条件\n",
        "\n",
        "このチュートリアルを進める前に、以下のトピックについてあらかじめ理解しておいていただくことをお勧めします：\n",
        "\n",
        "* [ハミルトニアンシミュレーション回路の構成手法](/docs/tutorials/compilation-methods-for-hamiltonian-simulation-circuits) — Qiskitにおけるトロッター（積の公式）回路について紹介します。\n",
        "* Qiskit のプロダクト式、特に および [`LieTrotter`](/docs/api/qiskit/qiskit.synthesis.LieTrotter) の [`SuzukiTrotter`](/docs/api/qiskit/qiskit.synthesis.SuzukiTrotter) 合成クラス。\n",
        "* [Qiskit primitives およびEstimatorインターフェース](/docs/guides/primitives)。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "07273b26",
      "metadata": {},
      "source": [
        "<span id=\"background\" />\n",
        "\n",
        "## 背景\n",
        "\n",
        "<span id=\"what-are-multi-product-formulas\" />\n",
        "\n",
        "### マルチプロダクト・フォーミュラとは何ですか？\n",
        "\n",
        "量子コンピュータ上で量子系をシミュレーションする際、中心的な課題は、ハミルトニアン $H$ に対する時間発展演算子 $e^{-iHt}$ を近似することである。標準的なアプローチでは、トロッター・鈴木分解としても知られる*積の公式* （PF）が用いられる。 これらは、 $H = \\sum_{a=1}^d F_a$ を、個々のユニタリー演算 $e^{-iF_a t}$ が効率的に実装できる項に分解し、その後、これらのより単純なユニタリー演算の順序付き積として、完全な進化を近似する。\n",
        "\n",
        "一階積の公式（リー・トロッターの公式）は次のとおりである：\n",
        "\n",
        "$$\n",
        "S_1(t) := \\prod_{a=1}^d e^{-i F_a t},\n",
        "$$\n",
        "\n",
        "これにより、二次誤差が生じます： $S_1(t) = e^{-iHt} + \\mathcal{O}(t^2)$。高次の対称式 $S_{2\\chi}(t)$ （ここで、 $\\chi$ は対称積の式の次数を表します［参考文献 [\\[1\\]](#references) ］）は、 $e^{-iHt} + \\mathcal{O}(t^{2\\chi+1})$ としてより速く収束しますが、その代償として、1ステップあたりの回路の深さが増加します。\n",
        "\n",
        "$\\chi$ *の固定*次数における誤差を低減するために、通常、総進化時間 $t$ を $k$ の小さなトロッターステップに分割する。 各ステップでは、 $e^{-iHt/k}$ を積の公式を用いて近似し、これらのステップを連結します：\n",
        "\n",
        "$$\n",
        "e^{-iHt} \\approx \\left[S_{2\\chi}(t/k)\\right]^k.\n",
        "$$\n",
        "\n",
        "$2\\chi$ 次の対称式の場合、残留トロッター誤差は $\\mathcal{O}\\!\\left(t^{2\\chi+1} / k^{2\\chi}\\right)$ の割合で増加する。したがって、 $k$ を増加させるとトロッター誤差は急速に抑制されるが、同時に回路の深さも線形に増大し、ノイズの多いハードウェアでは、ゲートノイズの累積が増加することになる。 **トロッター誤差（ $k$ の値が大きくなる傾向）**\\*\\* とハードウェアノイズ（ $k$ の値が小さくなる傾向）\\*\\* との間のこの緊張関係こそが、多製品式が解決するために設計された問題そのものである。 なお、MPFは、固定された順序 $\\chi$ で、 *$k$ の異なる選択肢*からの結果を組み合わせるものであり、基礎となる積の公式の順序を変えるものではないことに注意してください。\n",
        "\n",
        "**マルチプロダクト式（MPF）**[ \\[1\\]](#references) は、それぞれ異なる数のトロッターステップ $k_1, k_2, \\ldots, k_r$ （ $r$ のステップ数の集合）を用いる、いくつかの浅いトロッター回路から得られた期待値の*重み付き線形結合*を構成する：\n",
        "\n",
        "$$\n",
        "\\langle A \\rangle_{\\text{MPF}}(t) = \\sum_{j=1}^r x_j \\, \\langle A \\rangle_{k_j}(t),\n",
        "$$\n",
        "\n",
        "ここで、 $\\langle A \\rangle_{k_j}(t)$ は、時刻 $t$ における観測量 $A$ の期待値であり、 $k_j$ ステップのトロッター回路から推定されたものである。また、係数 $\\{x_j\\}_{j=1}^r$ は、組み合わせにおける主要なトロッター誤差項が相殺されるように選ばれている。 この式については[ステップ4](#small-scale-step-4) で改めて取り上げ、そこで明示的に評価して、トロッターの定理の結果と組み合わせます。 実用上の重要なポイントは、MPFの最深層回路に必要なステップ数が $k_{\\max}$ に過ぎないという点であり、これは、同じ有効トロッター誤差に直接到達するために必要となる単一の $k$ よりもはるかに少ない。 回路の深さが浅いことで、MPFアプローチはノイズの多いハードウェアにより適したものとなります。\n",
        "\n",
        "<span id=\"how-are-the-coefficients-determined\" />\n",
        "\n",
        "### 係数はどのように決定されるのですか？\n",
        "\n",
        "MPF係数には、次の2つの系統があります：\n",
        "\n",
        "**静的係数は**、ハミルトニアン、初期状態、および進化時間とは無関係である。 これらは、先行するトロッター誤差項の消去を強制する線形連立方程式 $Ax = b$ を解くことによって求められます。 $2\\chi$ 次の対称積の公式と組み合わせて用いられる一連のトロッター段階（ $\\{k_j\\}_{j=1}^r$ ）について、 $k_j$ の逆数で表されるトロッター誤差を展開すると、次のような形式の制約方程式が得られる：\n",
        "\n",
        "$$\n",
        "\\sum_{j=1}^r x_j = 1, \\quad \\sum_{j=1}^r \\frac{x_j}{k_j^{\\eta_n}} = 0 \\quad (n = 0, \\ldots, r-2),\n",
        "$$\n",
        "\n",
        "ここで、整数の指数 $\\{\\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\\]](#references) を参照のこと。 上記の方程式系の最初の方程式は、不偏性を保証する（ $k_j \\to \\infty$ の極限において、MPFは正確な期待値を再現する）ものであり、残りの $r-1$ の方程式は、最初の $r-1$ のトロッター誤差項を順次打ち消していく。 結果として得られる $L_1$ -ノルム $\\|x\\|_1$ が大きすぎる場合（これによりサンプリングノイズが増幅される）、代わりに、 $\\|x\\|_1$ を上限としつつ、 $\\|Ax - b\\|$ を最小化する近似最適化問題を解くことができます。\n",
        "\n",
        "**動的係数** [\\[2\\]](#references)、 [\\[3\\]](#references) は、さらにハミルトニアン、初期状態、および進化時間 $t$ にも依存する。これらは、真の時間発展状態とMPF近似との間のフロベニウスノルム距離を最小化する：\n",
        "\n",
        "$$\n",
        "\\|\\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),\n",
        "$$\n",
        "\n",
        "ここで、 $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-based のバックエンドを用いて効率的に計算します。\n",
        "\n",
        "<span id=\"when-to-use-mpfs\" />\n",
        "\n",
        "### MPFをいつ利用すべきか\n",
        "\n",
        "MPFは、次のような場合に最も有益です：\n",
        "\n",
        "* **回路の深さがボトルネックとなっている。** ハードウェアノイズによって実行可能な深さに制限がある場合は、MPF を使用することで、より浅い回路からより高い実効トロッター精度を実現できます。\n",
        "* **必要なのは正確な期待値であり、完全な状態の準備ではない。** MPFは期待値のレベルで動作する――つまり、量子状態ではなく、古典的な数を組み合わせるものである。 したがって、これらはEstimatorプリミティブを使用する場合の観測可能推定に最適です。\n",
        "* **トロッターの歩数を、それほど多くない数だけ組み合わせます。** 通常、 $r = 3$ と $5$ の異なるステップ数を組み合わせる $k_j$ ことで、 $\\|x\\|_1$ を扱いやすい範囲に保ちつつ、先行するいくつかのトロッター誤差項を相殺するのに十分である。\n",
        "\n",
        "<span id=\"when-mpfs-might-not-help\" />\n",
        "\n",
        "### MPFが役に立たない場合\n",
        "\n",
        "* **進化の時間が非常に短い。** $t$ が十分に小さく、単一の下位順トロッター公式だけで十分な精度が得られる場合、複数の回路を実行するオーバーヘッドは不要となる。\n",
        "* **試験対策の課題。** MPFは、補正された量子状態ではなく、補正された*期待*値を生成する。 （例えば、別の量子サブルーチンの入力として）時間発展後の実際の状態が必要な場合、MPFは適用されません。\n",
        "* **収束条件に違反するトロッター法による歩数。** 静的係数の導出では、個々の $\\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$ の場合、サンプリングによるオーバーヘッド $\\propto \\|x\\|_1^2$ が、トロッター誤差の低減効果を上回る可能性がある。 詳細については[、「トロッターステップの選び方」ガイド](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html)をご覧ください。\n",
        "\n",
        "<span id=\"what-this-tutorial-covers\" />\n",
        "\n",
        "### このチュートリアルの内容\n",
        "\n",
        "このチュートリアルでは、MPFのワークフロー全体を2つの段階に分けて解説します。 まず、 **小規模な**シミュレータの例（10キュービットのハイゼンベルグ鎖）を用いて、問題の設定方法、静的および動的なMPF係数の計算方法、そして得られた期待値を厳密な対角化結果と比較する方法を示します。 続いて、 **大規模なハードウェア例** （50キュービットのXXZチェーン）を用いて、トランスパイルの方法、 IBM Quantum ハードウェア上でのエラー緩和機能付きの実行方法、およびMPF係数を用いた結果の後処理について解説します。 本稿では、標準のQiskitツールと併せて、この `qiskit_addon_mpf` パッケージを使用しています。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b5d478ce",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## 要件\n",
        "\n",
        "このチュートリアルを始める前に、以下のものがインストールされていることを確認してください：\n",
        "\n",
        "* Qiskit SDK v2.0 またはそれ以降で、 [可視化](/docs/api/qiskit/visualization)機能をサポートしているもの\n",
        "* Qiskit Runtime v0.22 またはそれ以降 (`pip install qiskit-ibm-runtime`)\n",
        "* Qiskit Aer シミュレータ (`pip install qiskit-aer`)\n",
        "* TeNPy バックエンドを備えたMPF Qiskitアドオン (`pip install \"qiskit-addon-mpf[tenpy]\"`)\n",
        "* Qiskit アドオンユーティリティ (`pip install qiskit-addon-utils`)\n",
        "* SciPy (`pip install scipy`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2584c37e",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## セットアップ\n",
        "\n",
        "以下では、このチュートリアル全体で使用されている*すべての*パッケージのインポートを、1つのセルにまとめています。 `XXPlusYYGate`また、隣接する `rxx` および `ryy` の回転を単一の に融合させるトランスパイラー・パスを `CollectAndCollapse` 定義する。 この処理は、ステップ1の回路構築時（ゲート数を少なく抑えるため）と、ステップ4で動的MPFの層構造を抽出する際（ TeNPy は、融合されていない回転のペアではなく、2量子ビットゲートを期待するため）の両方で間接的に適用されます。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "bf79f9e7",
      "metadata": {},
      "outputs": [],
      "source": [
        "import warnings\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from functools import partial\n",
        "from copy import deepcopy\n",
        "\n",
        "from qiskit import QuantumCircuit\n",
        "from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector\n",
        "from qiskit.synthesis import SuzukiTrotter\n",
        "from qiskit.transpiler import CouplingMap, PassManager\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "from qiskit.circuit.library import XXPlusYYGate\n",
        "from qiskit.transpiler.passes.optimization.collect_and_collapse import (\n",
        "    CollectAndCollapse,\n",
        "    collect_using_filter_function,\n",
        "    collapse_to_operation,\n",
        ")\n",
        "\n",
        "from qiskit_aer import AerSimulator\n",
        "from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService\n",
        "\n",
        "from qiskit_addon_utils.problem_generators import (\n",
        "    generate_xyz_hamiltonian,\n",
        "    generate_time_evolution_circuit,\n",
        ")\n",
        "from qiskit_addon_utils.slicing import slice_by_depth\n",
        "from qiskit_addon_mpf.static import setup_static_lse\n",
        "from qiskit_addon_mpf.dynamic import setup_dynamic_lse\n",
        "from qiskit_addon_mpf.costs import (\n",
        "    setup_exact_problem,\n",
        "    setup_sum_of_squares_problem,\n",
        "    setup_frobenius_problem,\n",
        ")\n",
        "from qiskit_addon_mpf.backends.tenpy_layers import (\n",
        "    LayerModel,\n",
        "    LayerwiseEvolver,\n",
        ")\n",
        "from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state\n",
        "\n",
        "from scipy.linalg import expm\n",
        "\n",
        "# Suppress TeNPy's `unit_cell_width` future-API warning. The default\n",
        "# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what\n",
        "# `CouplingMap.from_line(...)` produces here, so the warning is informational.\n",
        "warnings.filterwarnings(\n",
        "    \"ignore\",\n",
        "    message=r\".*unit_cell_width.*\",\n",
        "    category=UserWarning,\n",
        ")\n",
        "\n",
        "\n",
        "# --- Helper: collect XX + YY rotations into a single gate ---\n",
        "def filter_function(node):\n",
        "    return node.op.name in {\"rxx\", \"ryy\"}\n",
        "\n",
        "\n",
        "collect_function = partial(\n",
        "    collect_using_filter_function,\n",
        "    filter_function=filter_function,\n",
        "    split_blocks=True,\n",
        "    min_block_size=1,\n",
        ")\n",
        "\n",
        "\n",
        "def collapse_to_xx_plus_yy(block):\n",
        "    param = 0.0\n",
        "    for node in block.data:\n",
        "        param += node.operation.params[0]\n",
        "    return XXPlusYYGate(param)\n",
        "\n",
        "\n",
        "collapse_function = partial(\n",
        "    collapse_to_operation,\n",
        "    collapse_function=collapse_to_xx_plus_yy,\n",
        ")\n",
        "\n",
        "pm = PassManager()\n",
        "pm.append(CollectAndCollapse(collect_function, collapse_function))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "24f08467",
      "metadata": {},
      "source": [
        "<span id=\"small-scale-simulator-example\" />\n",
        "\n",
        "## 小規模シミュレータの例\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "378e82ba",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### ステップ1：古典的な入力を量子問題にマッピングする\n",
        "\n",
        "まず、直線上の10キュービットのハイゼンベルクモデルについて、初期状態としてネール状態 $\\vert 0101\\ldots01 \\rangle$ を用いる。 ハミルトニアンは次のとおりである：\n",
        "\n",
        "$$\n",
        "\\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),\n",
        "$$\n",
        "\n",
        "ここで、 $J$ は最近傍結合強度である。 チェーンの中央にある1組の量子ビットについて、ZZ相関関数 $Z_{L/2-1} Z_{L/2}$ を測定し、2次積公式を用いたトロッター法 $k_j = [1, 2, 4]$ を適用する。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "bdd0d4fc",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "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'],\n",
            "              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,\n",
            " 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,\n",
            " 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])\n"
          ]
        }
      ],
      "source": [
        "L = 10\n",
        "\n",
        "# Generate coupling map and Hamiltonian\n",
        "coupling_map = CouplingMap.from_line(L, bidirectional=False)\n",
        "\n",
        "hamiltonian = generate_xyz_hamiltonian(\n",
        "    coupling_map,\n",
        "    coupling_constants=(1.0, 1.0, 1.0),\n",
        "    ext_magnetic_field=(0.0, 0.0, 0.0),\n",
        ")\n",
        "print(hamiltonian)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "fd3dc9c8",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "SparsePauliOp(['IIIIZZIIII'],\n",
            "              coeffs=[1.+0.j])\n"
          ]
        }
      ],
      "source": [
        "# Observable: ZZ on the middle pair of qubits\n",
        "observable = SparsePauliOp.from_sparse_list(\n",
        "    [(\"ZZ\", (L // 2 - 1, L // 2), 1.0)], num_qubits=L\n",
        ")\n",
        "print(observable)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "398c33b2",
      "metadata": {},
      "outputs": [],
      "source": [
        "# MPF parameters\n",
        "mpf_trotter_steps = [1, 2, 4]\n",
        "order = 2\n",
        "symmetric = False\n",
        "\n",
        "trotter_times = np.arange(0.5, 1.55, 0.1)\n",
        "exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1512ba0e",
      "metadata": {},
      "source": [
        "<span id=\"build-trotter-circuits\" />\n",
        "\n",
        "#### トロッター回路を構築する\n",
        "\n",
        "各時点および各トロッターステップ数について、近似トロッター時間発展を実装する回路を作成します。 「セットアップ」セクションで定義されたパスは `CollectAndCollapse` 、XXおよびYYの回転を単一のXX+YYゲートにまとめ、後のテンソルネットワークのシミュレーションをより効率的に行うための準備を行います。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "1c194d2b",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Initial Neel state preparation\n",
        "initial_state_circ = QuantumCircuit(L)\n",
        "initial_state_circ.x([i for i in range(L) if i % 2 != 0])\n",
        "\n",
        "\n",
        "all_circs = []\n",
        "for total_time in trotter_times:\n",
        "    mpf_trotter_circs = [\n",
        "        generate_time_evolution_circuit(\n",
        "            hamiltonian,\n",
        "            time=total_time,\n",
        "            synthesis=SuzukiTrotter(reps=num_steps, order=order),\n",
        "        )\n",
        "        for num_steps in mpf_trotter_steps\n",
        "    ]\n",
        "\n",
        "    mpf_trotter_circs = pm.run(\n",
        "        mpf_trotter_circs\n",
        "    )  # Collect XX and YY into XX + YY\n",
        "\n",
        "    mpf_circuits = [\n",
        "        initial_state_circ.compose(circuit) for circuit in mpf_trotter_circs\n",
        "    ]\n",
        "    all_circs.append(mpf_circuits)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "c7ee61e7",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/multi-product-formula/extracted-outputs/c7ee61e7-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "mpf_circuits[-1].draw(\"mpl\", fold=-1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4cd6c782",
      "metadata": {},
      "source": [
        "<span id=\"step-2-optimize-problem-for-quantum-hardware-execution\" />\n",
        "\n",
        "### ステップ2：量子ハードウェア実行に向けた問題の最適化\n",
        "\n",
        "小規模な例として、Aerシミュレータを対象とします。 回路が実行可能になる前に、2つの変換が行われます：\n",
        "\n",
        "1. **ハミルトニアンシミュレーションレベルにおけるゲート集積。** `XXPlusYYGate``CollectAndCollapse` 「セットアップ」セルでは、隣接する `rxx` と `ryy` の回転を1つの回転に融合させるパスを作成しました。 この処理は、ステップ1でトロッター回路を構築した際（その `pm.run(...)` 呼び出し）に、すでに適用済みです。 これにより、2量子ビットゲートの数が減少するだけでなく、後の動的係数の計算において、テンソルネットワークによるシミュレーションに適した構造が得られる。\n",
        "\n",
        "2. **シミュレータのISAに合わせて下位互換性を確保する。** 以下では、Qiskitのプリセット・パス・マネージャーを実行し、 `optimization_level=3` 各トロッター回路をシミュレータの命令セットアーキテクチャ（ISA）に展開します。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "03590a05",
      "metadata": {},
      "outputs": [],
      "source": [
        "aer_sim = AerSimulator()\n",
        "pm_sim = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)\n",
        "\n",
        "isa_circs_all_times = [\n",
        "    pm_sim.run([deepcopy(c) for c in mpf_circuits])\n",
        "    for mpf_circuits in all_circs\n",
        "]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ab6588d3",
      "metadata": {},
      "source": [
        "<span id=\"step-3-execute-using-qiskit-primitives\" />\n",
        "\n",
        "### ステップ3: `Qiskit primitives`を使用して実行する\n",
        "\n",
        "小規模な例として、ISA低減されたトロッター回路を、Aerがバックエンドとなるプリミティブ `EstimatorV2` を通じて実行します。 これにより、各 $(k_j, t)$ のペアに対して*ノイズのない*基準値が得られます。これらは、ステップ4でMPFが組み合わせる $\\langle A \\rangle_{k_j}(t)$ の値となります。 後で各製品処方の時系列曲線全体およびMPFの時系列曲線全体をプロットできるように、進化の経過時間を概観します。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "7225d782",
      "metadata": {},
      "outputs": [],
      "source": [
        "estimator = Estimator(mode=aer_sim)\n",
        "\n",
        "mpf_expvals_all_times, mpf_stds_all_times = [], []\n",
        "for isa_circuits in isa_circs_all_times:\n",
        "    result = estimator.run(\n",
        "        [(circuit, observable) for circuit in isa_circuits], precision=0.005\n",
        "    ).result()\n",
        "    mpf_expvals_all_times.append([res.data.evs for res in result])\n",
        "    mpf_stds_all_times.append([res.data.stds for res in result])"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a384f017",
      "metadata": {},
      "source": [
        "<span id=\"small-scale-step-4\" />\n",
        "\n",
        "<span id=\"step-4-post-process-and-return-result-in-desired-classical-format\" />\n",
        "\n",
        "### ステップ4：後処理を行い、結果を希望の古典形式で返す\n",
        "\n",
        "ステップ4では、MPFが実際に構築されます。 ここでは係数 $x_j$\\* が計算\\*されますが（動的バリエーションの場合、この計算は負荷がかかることがあります）、概念的には、これらはステップ3の量子測定結果を単一の補正済み期待値に組み合わせるための古典的な手法であるため、係数の算出と組み合わせのワークフロー全体を後処理として扱います。\n",
        "\n",
        "MPFが実際のダイナミクスをどの程度正確に追跡しているかを評価するために、まずハミルトニアンを直接指数関数化することで、時間発展を経た厳密な期待値を計算する。 これが処理可能であるのは、 $L = 10$ であるからに過ぎない。以下の大規模ハードウェアの例では、代わりにテンソルネットワークによる推定に頼らざるを得ない。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "a223a360",
      "metadata": {},
      "outputs": [],
      "source": [
        "exact_expvals = []\n",
        "for t in exact_evolution_times:\n",
        "    exp_H = expm(-1j * t * hamiltonian.to_matrix())\n",
        "    initial_state = Statevector(initial_state_circ).data\n",
        "    time_evolved_state = exp_H @ initial_state\n",
        "\n",
        "    exact_obs = (\n",
        "        time_evolved_state.conj()\n",
        "        @ observable.to_matrix()\n",
        "        @ time_evolved_state\n",
        "    ).real\n",
        "    exact_expvals.append(exact_obs)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "361e726e",
      "metadata": {},
      "source": [
        "<span id=\"static-mpf-coefficients\" />\n",
        "\n",
        "#### 静的MPF係数\n",
        "\n",
        "静的MPFでは、進化時間、ハミルトニアン、および初期状態に依存しない係数 $x_j$ が用いられる。 「背景」で説明した線形システム $Ax = b$ を設定し、係数を求めます。 行列 $A$ は、トロッターのステップ数 $k_j$、積の公式の次数 $\\chi$、およびその公式が対称であるかどうか（これにより指数 $\\eta_n$ が決まる）によって決定される。\n",
        "\n",
        "この小規模な例では、 $k_j = [1, 2, 4]$ を用い、非対称な $2\\chi=2$ のスズキ・トロッターの公式（したがって、 $\\chi=1$ および $\\eta_n = 2 + n$ となり、 $\\eta_0 = 2,\\, \\eta_1 = 3$ となる）を採用する。これにより、系は次のようになる：\n",
        "\n",
        "$$\n",
        "A =\n",
        "\\begin{bmatrix}\n",
        "1 & 1 & 1\\\\\n",
        "1 & \\frac{1}{2^2} & \\frac{1}{4^2}  \\\\\n",
        "1 & \\frac{1}{2^3} & \\frac{1}{4^3}  \\\\\n",
        "\\end{bmatrix}, \\quad\n",
        "b =\n",
        "\\begin{bmatrix}\n",
        "1 \\\\\n",
        "0 \\\\\n",
        "0\n",
        "\\end{bmatrix}.\n",
        "$$\n",
        "\n",
        "1行目は不偏性を保証する（ $\\sum_j x_j = 1$ ）。2行目と3行目は、それぞれ、先行する $1/k^2$ および次次の $1/k^3$ のトロッター誤差項を相殺する。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4f2ca1e2",
      "metadata": {},
      "source": [
        "<span id=\"set-up-the-lse\" />\n",
        "\n",
        "##### LSEを設定する\n",
        "\n",
        "上述した行列 $A$ および右辺ベクトル $b$ を構成するために、から `qiskit_addon_mpf.static` を使用 `setup_static_lse` する。 行列 $A$ は、 $k_j$ だけでなく、積の公式の選び方――特にその*次数* $\\chi$ や、それが*対称*であるかどうか――にも依存する。 この `symmetric` フラグは、 $\\eta_n$ の指数パターンを制御します（対称式では、トロッター誤差項は偶数次のものしか生成されません。参考文献 [\\[1\\]](#references) を参照）。 なお、参考文献 [\\[2\\]](#references) に示されているように、基礎となるPFが対称である場合でも、 を設定 `symmetric=True` することは厳密には必須ではない。非対称のLSEは依然として有効である（ただし、不要な追加の制約が課されることになる）。\n",
        "\n",
        "この例では、ステップ1ですでに と `symmetric = False` を設定 `order = 2` しています。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "827b0b42",
      "metadata": {},
      "outputs": [],
      "source": [
        "lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "003e4bdd",
      "metadata": {},
      "source": [
        "構築された行列 $A$ およびベクトル $b$ を確認し、それらが上記で記述したシステムと一致していることを確認してください。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "6f879978",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[1.      , 1.      , 1.      ],\n",
              "       [1.      , 0.25    , 0.0625  ],\n",
              "       [1.      , 0.125   , 0.015625]])"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "lse.A"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "64fe7db9",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([1., 0., 0.])"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "lse.b"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "16ab9f79",
      "metadata": {},
      "source": [
        "LSEを適用して、静的係数 $x_j$ を次式を用いて `lse.solve()` 求める（これが $x = A^{-1}b$ の直接解法である）。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "7b69192a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The static coefficients associated with the ansatze are: [ 0.04761905 -0.57142857  1.52380952]\n"
          ]
        }
      ],
      "source": [
        "mpf_coeffs = lse.solve()\n",
        "print(\n",
        "    f\"The static coefficients associated with the ansatze are: {mpf_coeffs}\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1c238a59",
      "metadata": {},
      "source": [
        "<span id=\"optimize-for-$x$-using-an-exact-model\" />\n",
        "\n",
        "##### $x$ 向けに正確なモデルを用いて最適化\n",
        "\n",
        "$x = A^{-1}b$ を計算する代わりに、 [setup\\_exact\\_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_exact_model.html) を使用して、LSE を制約条件とする [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) インスタンスを構築することもできます。このインスタンスの最適解は、 $x$ となります。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "993465e9",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[ 0.04761905 -0.57142857  1.52380952]\n"
          ]
        }
      ],
      "source": [
        "model_exact, coeffs_exact = setup_exact_problem(lse)\n",
        "model_exact.solve()\n",
        "print(coeffs_exact.value)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "eb61ea70",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "L1 norm of the exact coefficients: 2.1428571428556378\n"
          ]
        }
      ],
      "source": [
        "print(\n",
        "    \"L1 norm of the exact coefficients:\",\n",
        "    np.linalg.norm(coeffs_exact.value, ord=1),\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "dac0472b",
      "metadata": {},
      "source": [
        "<span id=\"optimize-for-$x$-using-an-approximate-model\" />\n",
        "\n",
        "##### 近似モデルを用いた $x$ の最適化\n",
        "\n",
        "選択した $k_j$ の値の集合に対する $L_1$ ノルムが、高すぎるとみなされる場合がある。 その場合、 $k_j$ の値を別の組み合わせに選択できないときは、 $L_1$ ノルムを所定の閾値以下に制限しつつ、 $\\|Ax - b\\|$ を最小化する近似解を使用することができます。 [「近似モデルの使用方法」](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html) に関するガイドをご覧ください。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "0cd7dea4",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[-1.10294118e-03 -2.48897059e-01  1.25000000e+00]\n",
            "L1 norm of the approximate coefficients: 1.5\n"
          ]
        }
      ],
      "source": [
        "model_approx, coeffs_approx = setup_sum_of_squares_problem(\n",
        "    lse, max_l1_norm=1.5\n",
        ")\n",
        "model_approx.solve()\n",
        "print(coeffs_approx.value)\n",
        "print(\n",
        "    \"L1 norm of the approximate coefficients:\",\n",
        "    np.linalg.norm(coeffs_approx.value, ord=1),\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "10fd05eb",
      "metadata": {},
      "source": [
        "<span id=\"dynamic-mpf-coefficients\" />\n",
        "\n",
        "#### 動的MPF係数\n",
        "\n",
        "静的MPFは、ハミルトニアンや状態に依存しない方法でトロッター誤差項を打ち消すため、与えられたハミルトニアンと初期状態に対して、必ずしも最小の近似誤差をもたらすとは限らない。 `qiskit_addon_mpf`一方、動的MPF（参考文献 [\\[2\\]](#references)、 [\\[3\\]](#references) ）は、各時刻 $t$ において、フロベニウスノルムの距離 $\\|\\rho(t) - \\mu^D(t)\\|_F^2$ を最小化する時間依存係数 $x_i(t)$ を求める。背景の項で示したように、これには、トロッター進化状態間の重なり行列 $M_{ij}(t)$ と、正確な状態との重なり $L_i(t)$ が必要となる。これらはいずれも、本論文においてテンソルネットワーク（ TeNPy ）バックエンドを用いて推定する。\n",
        "\n",
        "動的なLSEを設定するには、次の3つの要素が必要です：\n",
        "\n",
        "1. このアドオンが、各 $k_j$ に対して実行し、 $\\rho_{k_j}(t)$ をMPS/MPOとして生成する**近似エヴォルバーファクトリ**。 `slice_by_depth`これは、 $2$ の順序を持つトロッター回路の層状構造（1層につき1つ）から構成され、 TeNPy の切り捨てパラメータを持つ形でラップされています。 `LayerwiseEvolver`\n",
        "2. 高精度な基準 $\\rho(t)$ を生成する、 **厳密な進化ファクトリ**。厳密な進化の近似として、微小時間ステップの4次スズキ・トロッター回路（`dt=0.1`, `order=4`）を用いる。\n",
        "3. TeNPy シミュレーションの初期化に用いる、 **ID** ファクトリと**初期状態**のMPS。\n",
        "\n",
        "以下のセルは、近似エヴォルバーファクトリを構築します。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "52d78403",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Create approximate time-evolution circuits\n",
        "single_2nd_order_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)\n",
        ")\n",
        "single_2nd_order_circ = pm.run(single_2nd_order_circ)  # collect XX and YY\n",
        "\n",
        "# Find layers in the circuit\n",
        "layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)\n",
        "\n",
        "# Create tensor network models\n",
        "models = [\n",
        "    LayerModel.from_quantum_circuit(layer, conserve=\"Sz\") for layer in layers\n",
        "]\n",
        "\n",
        "# Create the time-evolution object\n",
        "approx_factory = partial(\n",
        "    LayerwiseEvolver,\n",
        "    layers=models,\n",
        "    options={\n",
        "        \"preserve_norm\": False,\n",
        "        \"trunc_params\": {\n",
        "            \"chi_max\": 64,\n",
        "            \"svd_min\": 1e-8,\n",
        "            \"trunc_cut\": None,\n",
        "        },\n",
        "        \"max_delta_t\": 2,\n",
        "    },\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "91d2783e",
      "metadata": {},
      "source": [
        "<Admonition type=\"warning\">\n",
        "  テンソルネットワークシミュレーションの詳細を決定する `LayerwiseEvolver` のオプションは、定義されていない最適化問題を設定しないように注意深く選択する必要がある。\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1c702f67",
      "metadata": {},
      "source": [
        "`dt=0.1`小さな時間ステップを用いて、4次の鈴木・トロッターの公式により、時間発展を経た正確な状態を近似する。 TeNPy の切り捨てパラメータは精度に影響を与える可能性があるため、さまざまな値を試してみることが重要です。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "abab8bfc",
      "metadata": {},
      "outputs": [],
      "source": [
        "single_4th_order_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)\n",
        ")\n",
        "single_4th_order_circ = pm.run(single_4th_order_circ)\n",
        "exact_model_layers = [\n",
        "    LayerModel.from_quantum_circuit(layer, conserve=\"Sz\")\n",
        "    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)\n",
        "]\n",
        "\n",
        "exact_factory = partial(\n",
        "    LayerwiseEvolver,\n",
        "    layers=exact_model_layers,\n",
        "    dt=0.1,\n",
        "    options={\n",
        "        \"preserve_norm\": False,\n",
        "        \"trunc_params\": {\n",
        "            \"chi_max\": 64,\n",
        "            \"svd_min\": 1e-8,\n",
        "            \"trunc_cut\": None,\n",
        "        },\n",
        "        \"max_delta_t\": 2,\n",
        "    },\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2486594c",
      "metadata": {},
      "source": [
        "最後に、初期のMPO状態をもたらす を `identity_factory` 定義し、層状トロッターモデルで使用される格子と一致するMPSとして、ネール初期状態を準備する。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "1e216575",
      "metadata": {},
      "outputs": [],
      "source": [
        "def identity_factory():\n",
        "    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)\n",
        "\n",
        "\n",
        "mps_initial_state = MPS_neel_state(models[0].lat)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5658315c",
      "metadata": {},
      "source": [
        "工場の配置が決まったので、各進化時点における動的係数を計算する。 各 $t$ について、 `setup_dynamic_lse`TeNPy, を用いて関連するオーバーラップ行列を構築し、 `setup_frobenius_problem` フロベニウスノルムコストを最小化するa `cvxpy.Problem` を返す。 `mpf_dynamic_coeffs_list`ソルバーは、その時刻に合わせた係数 $x_j(t)$ を返す。これらを.に収集する。 特定の $t$ に対してソルバーが失敗した場合、係数をゼロに設定してループを継続させます。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "b05dc012",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Computing dynamic coefficients for time=0.5\n",
            "\n",
            "Computing dynamic coefficients for time=0.6\n",
            "\n",
            "Computing dynamic coefficients for time=0.7\n",
            "\n",
            "Computing dynamic coefficients for time=0.7999999999999999\n",
            "\n",
            "Computing dynamic coefficients for time=0.8999999999999999\n",
            "\n",
            "Computing dynamic coefficients for time=0.9999999999999999\n",
            "\n",
            "Computing dynamic coefficients for time=1.0999999999999999\n",
            "\n",
            "Computing dynamic coefficients for time=1.1999999999999997\n",
            "\n",
            "Computing dynamic coefficients for time=1.2999999999999998\n",
            "\n",
            "Computing dynamic coefficients for time=1.4\n",
            "\n",
            "Computing dynamic coefficients for time=1.4999999999999998\n",
            "\n"
          ]
        }
      ],
      "source": [
        "mpf_dynamic_coeffs_list = []\n",
        "for t in trotter_times:\n",
        "    print(f\"Computing dynamic coefficients for time={t}\")\n",
        "    lse = setup_dynamic_lse(\n",
        "        mpf_trotter_steps,\n",
        "        t,\n",
        "        identity_factory,\n",
        "        exact_factory,\n",
        "        approx_factory,\n",
        "        mps_initial_state,\n",
        "    )\n",
        "    problem, coeffs = setup_frobenius_problem(lse)\n",
        "    try:\n",
        "        problem.solve()\n",
        "        mpf_dynamic_coeffs_list.append(coeffs.value)\n",
        "    except Exception as error:\n",
        "        mpf_dynamic_coeffs_list.append(np.zeros(len(mpf_trotter_steps)))\n",
        "        print(error, \"Calculation Failed for time\", t)\n",
        "    print(\"\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4f8814e3",
      "metadata": {},
      "source": [
        "<span id=\"combine-trotter-expectation-values-with-the-mpf-coefficients\" />\n",
        "\n",
        "#### トロッターの期待値とMPF係数を組み合わせる\n",
        "\n",
        "ここで、各係数セット（static-exact、static-approximate、および dynamic）について、 $\\langle A \\rangle_{\\text{MPF}}(t) = \\sum_j x_j \\, \\langle A \\rangle_{k_j}(t)$ を評価し、回路ごとの標準誤差を伝播させ、その結果得られた時系列を、厳密対角化曲線と対比してプロットする。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "id": "35042576",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/multi-product-formula/extracted-outputs/35042576-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "sym = {1: \"^\", 2: \"s\", 4: \"p\"}\n",
        "# Get expectation values at all times for each Trotter step\n",
        "for k, step in enumerate(mpf_trotter_steps):\n",
        "    trotter_curve, trotter_curve_error = [], []\n",
        "    for trotter_expvals, trotter_stds in zip(\n",
        "        mpf_expvals_all_times, mpf_stds_all_times\n",
        "    ):\n",
        "        trotter_curve.append(trotter_expvals[k])\n",
        "        trotter_curve_error.append(trotter_stds[k])\n",
        "\n",
        "    plt.errorbar(\n",
        "        trotter_times,\n",
        "        trotter_curve,\n",
        "        yerr=trotter_curve_error,\n",
        "        alpha=0.5,\n",
        "        markersize=4,\n",
        "        marker=sym[step],\n",
        "        color=\"grey\",\n",
        "        label=f\"{mpf_trotter_steps[k]} Trotter steps\",\n",
        "    )\n",
        "\n",
        "# Get expectation values at all times for the static MPF with exact coeffs\n",
        "exact_mpf_curve, exact_mpf_curve_error = [], []\n",
        "for trotter_expvals, trotter_stds in zip(\n",
        "    mpf_expvals_all_times, mpf_stds_all_times\n",
        "):\n",
        "    mpf_std = np.sqrt(\n",
        "        sum(\n",
        "            [\n",
        "                (coeff**2) * (std**2)\n",
        "                for coeff, std in zip(coeffs_exact.value, trotter_stds)\n",
        "            ]\n",
        "        )\n",
        "    )\n",
        "    exact_mpf_curve_error.append(mpf_std)\n",
        "    exact_mpf_curve.append(trotter_expvals @ coeffs_exact.value)\n",
        "\n",
        "plt.errorbar(\n",
        "    trotter_times,\n",
        "    exact_mpf_curve,\n",
        "    yerr=exact_mpf_curve_error,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    label=\"Static MPF - Exact\",\n",
        "    color=\"purple\",\n",
        ")\n",
        "\n",
        "\n",
        "# Get expectation values at all times for the static MPF with approximate coeffs\n",
        "approx_mpf_curve, approx_mpf_curve_error = [], []\n",
        "for trotter_expvals, trotter_stds in zip(\n",
        "    mpf_expvals_all_times, mpf_stds_all_times\n",
        "):\n",
        "    mpf_std = np.sqrt(\n",
        "        sum(\n",
        "            [\n",
        "                (coeff**2) * (std**2)\n",
        "                for coeff, std in zip(coeffs_approx.value, trotter_stds)\n",
        "            ]\n",
        "        )\n",
        "    )\n",
        "    approx_mpf_curve_error.append(mpf_std)\n",
        "    approx_mpf_curve.append(trotter_expvals @ coeffs_approx.value)\n",
        "\n",
        "plt.errorbar(\n",
        "    trotter_times,\n",
        "    approx_mpf_curve,\n",
        "    yerr=approx_mpf_curve_error,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    label=\"Static MPF - Approx\",\n",
        "    color=\"orange\",\n",
        ")\n",
        "\n",
        "\n",
        "# Get expectation values at all times for the dynamic MPF\n",
        "dynamic_mpf_curve, dynamic_mpf_curve_error = [], []\n",
        "for trotter_expvals, trotter_stds, dynamic_coeffs in zip(\n",
        "    mpf_expvals_all_times, mpf_stds_all_times, mpf_dynamic_coeffs_list\n",
        "):\n",
        "    mpf_std = np.sqrt(\n",
        "        sum(\n",
        "            [\n",
        "                (coeff**2) * (std**2)\n",
        "                for coeff, std in zip(dynamic_coeffs, trotter_stds)\n",
        "            ]\n",
        "        )\n",
        "    )\n",
        "    dynamic_mpf_curve_error.append(mpf_std)\n",
        "    dynamic_mpf_curve.append(trotter_expvals @ dynamic_coeffs)\n",
        "\n",
        "plt.errorbar(\n",
        "    trotter_times,\n",
        "    dynamic_mpf_curve,\n",
        "    yerr=dynamic_mpf_curve_error,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    label=\"Dynamic MPF\",\n",
        "    color=\"pink\",\n",
        ")\n",
        "\n",
        "\n",
        "# Exact expectation values\n",
        "plt.plot(\n",
        "    exact_evolution_times,\n",
        "    exact_expvals,\n",
        "    color=\"red\",\n",
        "    linestyle=\"--\",\n",
        "    label=\"Exact time-evolution\",\n",
        ")\n",
        "\n",
        "plt.title(f\"$\\\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\\\rangle$ vs time\")\n",
        "plt.xlabel(\"Time\")\n",
        "plt.ylabel(\"Expectation Value\")\n",
        "plt.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.2), ncol=2)\n",
        "plt.grid(alpha=0.1)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "34923748",
      "metadata": {},
      "source": [
        "上の図は、トロッター誤差とサンプリング誤差の相互作用を示しています。\n",
        "\n",
        "* **トロッターのエラー。** 個々の製品の配合（灰色のマーカー）は、時間が経つにつれて正確な曲線からますます乖離していく。 $k=1$ の回路は偏差が最も大きく、深さも最も浅いですが、すでに $t/k \\gtrsim 1$ という領域にあるため、主項である $1/k^{2}$ の誤差項が大きくなっています。 MPFの組み合わせ（色付きのマーカー）は、これらの先行するトロッター誤差項のいくつかを打ち消すため、単一の $k_j$ 回路よりもはるかに正確に曲線に追従します。 残りの誤差は、MPFでは*打ち*消されない高次のトロッター項を反映している。 $2$、 $r=3$ の静的MPFは最初の2次の誤差しか除去できず、 $t/k_{\\min}$ が十分に大きい場合、打ち消されなかった尾部項が最終的に支配的となる。したがって、MPFでは、非常に浅い回路が任意の時点で正確さを維持することが保証されない。\n",
        "\n",
        "* **標本誤差。** MPF曲線の誤差バーが広くなっているのは、線形結合の直接的な結果である。回路ごとの独立した標準誤差 $\\sigma_{k_j}$ を伝播させると、総分散 $\\sigma_{\\text{MPF}}^2 = \\sum_j x_j^2 \\, \\sigma_{k_j}^2$ が得られる。したがって、 $\\|x\\|_2$ （実際には、我々が制御対象としている $\\|x\\|_1$ ）が大きければ大きいほど、所定の目標不確実性に到達するために必要な測定回数も多くなる。 これが「Background」の近似ソルバーオプションの背後にあるトレードオフです。このオーバーヘッドを許容範囲内に抑えるため、 $\\|x\\|_1$ に上限を設けています。 重要な点として、トロッター誤差とは異なり、サンプリング誤差は $1/\\sqrt{N_{\\text{shots}}}$ に比例して小さくなるため、ショット数を増やすことで常に低減させることができる。\n",
        "\n",
        "以下の大規模なハードウェアの例では、各 $\\langle A \\rangle_{k_j}$ にハードウェアノイズが追加の誤差源として混入し、これも同様にMPF係数によって増幅されます。 そのセクションでは、エラー緩和がMPFとどのように相互作用するかについて見ていきます。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6fa763ff",
      "metadata": {},
      "source": [
        "<span id=\"large-scale-hardware-example\" />\n",
        "\n",
        "## 大規模なハードウェアの例\n",
        "\n",
        "このセクションでは、問題を、正確にシミュレーションすることが不可能な規模まで拡大します。 参考文献 [\\[3\\]](#references) に示された結果の一部を、 $t = 3$ 時点における50量子ビットのXXZ鎖を用いて再現する。小規模な例と同様の4段階のワークフローに従い、今回はエラー緩和機能を備えた実際の量子ハードウェアを対象とする。 テンプレートと同様に、各ステップはコード内にインラインでマークされており、中間出力を確認する価値がある場合は、1つのステップが複数のセルにまたがることもあります。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "481fea70",
      "metadata": {},
      "source": [
        "このマッピングは、小規模な例と同様の手順を踏んでいます。すなわち、ハミルトニアンを定義し、トロッターパラメータを選択し、（静的および動的な）MPF係数を計算し、回路を構築します。 主な違いは以下の通りです：\n",
        "\n",
        "* $\\mathcal{U}(0.5, 1.5)$ （参考文献 [\\[3\\]](#references) ）から抽出したランダムな結合を持つ、 **5** 0サイトからなるXXZハミルトニアン。\n",
        "* `symmetric=True`$k_j = [3, 4, 6]$ （したがって、 $\\chi=1$、）を**満た**す対称な2次トロッターの公式。\n",
        "* 単一の固定進化時間 $t = 3$。 $k_{\\min}=3$ を用いると、 $t/k_{\\min}=1$ となり、浅い構成要素をトロッター収束領域内に収めることができる。この領域では、MPFが依存する主誤差モデルが有効である。\n",
        "* **$k = 10$ のトロッター法を用いた単一回路の比較実験**を1回追加で実施し、これを基準とした。 $k = 10$ を選んだ理由は、そのハードウェア上の2量子ビットの深さが、最も深いMPF構成要素（ $k_{\\max}=6$ ）に、複数のMPF回路を実行するためのオーバーヘッドを加えた値よりも深いからである。この深さはノイズ制限領域に達しており、この領域では、MPFの組み合わせが単一回路のベースラインよりも優れた性能を発揮すると期待される。 これは、MPFの組み合わせに対する「単一の深層回路」による比較であり、MPFの実効トロッター誤差をターゲットとした回路ではありません（後者の場合、はるかに多くのステップが必要となります）。\n",
        "\n",
        "なお、ここではまだステップ1（マッピングと回路の構築）の段階ですが、このセルでは静的係数とともに動的係数も事前に計算しています。 動的係数は、 $H$ および $t$ に依存するが、量子測定には依存しないため、ステップ4以前の任意の時点で計算することができる。 MPF固有の設定をすべて一か所にまとめるために、今これを行っています。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "id": "a019ac32",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Static coefficients: [ 0.42857143 -1.82857143  2.4       ]\n",
            "L1 norm: 4.65714285714286\n",
            "Approximate coefficients: [-0.4942491   0.40206845  1.09218065]\n",
            "L1 norm (approx): 1.9884981979026675\n",
            "Computing dynamic coefficients for time=3\n"
          ]
        }
      ],
      "source": [
        "# -------------------------Step 1-------------------------\n",
        "L = 50\n",
        "coupling_map = CouplingMap.from_line(L, bidirectional=False)\n",
        "\n",
        "# XXZ Hamiltonian with random couplings (Ref. [3])\n",
        "np.random.seed(0)\n",
        "even_edges = list(coupling_map.get_edges())[::2]\n",
        "odd_edges = list(coupling_map.get_edges())[1::2]\n",
        "\n",
        "Js = np.random.uniform(0.5, 1.5, size=L)\n",
        "hamiltonian = SparsePauliOp(Pauli(\"I\" * L))\n",
        "for i, edge in enumerate(even_edges + odd_edges):\n",
        "    hamiltonian += SparsePauliOp.from_sparse_list(\n",
        "        [\n",
        "            (\"XX\", (edge), 2 * Js[i]),\n",
        "            (\"YY\", (edge), 2 * Js[i]),\n",
        "            (\"ZZ\", (edge), 4 * Js[i]),\n",
        "        ],\n",
        "        num_qubits=L,\n",
        "    )\n",
        "\n",
        "observable = SparsePauliOp.from_sparse_list(\n",
        "    [(\"ZZ\", (L // 2 - 1, L // 2), 1.0)], num_qubits=L\n",
        ")\n",
        "\n",
        "total_time = 3\n",
        "mpf_trotter_steps = [3, 4, 6]\n",
        "order = 2\n",
        "symmetric = True\n",
        "\n",
        "# Static coefficients\n",
        "lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)\n",
        "mpf_coeffs = lse.solve()\n",
        "print(f\"Static coefficients: {mpf_coeffs}\")\n",
        "print(f\"L1 norm: {np.linalg.norm(mpf_coeffs, ord=1)}\")\n",
        "\n",
        "model_approx, coeffs_approx = setup_sum_of_squares_problem(\n",
        "    lse, max_l1_norm=2.0\n",
        ")\n",
        "model_approx.solve()\n",
        "print(f\"Approximate coefficients: {coeffs_approx.value}\")\n",
        "print(f\"L1 norm (approx): {np.linalg.norm(coeffs_approx.value, ord=1)}\")\n",
        "\n",
        "# -------------------------Dynamic coefficients-------------------------\n",
        "single_2nd_order_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)\n",
        ")\n",
        "single_2nd_order_circ = pm.run(single_2nd_order_circ)\n",
        "\n",
        "layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)\n",
        "models = [\n",
        "    LayerModel.from_quantum_circuit(layer, conserve=\"Sz\") for layer in layers\n",
        "]\n",
        "\n",
        "approx_factory = partial(\n",
        "    LayerwiseEvolver,\n",
        "    layers=models,\n",
        "    options={\n",
        "        \"preserve_norm\": False,\n",
        "        \"trunc_params\": {\"chi_max\": 64, \"svd_min\": 1e-8, \"trunc_cut\": None},\n",
        "        \"max_delta_t\": 4,\n",
        "    },\n",
        ")\n",
        "\n",
        "single_4th_order_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)\n",
        ")\n",
        "single_4th_order_circ = pm.run(single_4th_order_circ)\n",
        "exact_model_layers = [\n",
        "    LayerModel.from_quantum_circuit(layer, conserve=\"Sz\")\n",
        "    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)\n",
        "]\n",
        "\n",
        "exact_factory = partial(\n",
        "    LayerwiseEvolver,\n",
        "    layers=exact_model_layers,\n",
        "    dt=0.1,\n",
        "    options={\n",
        "        \"preserve_norm\": False,\n",
        "        \"trunc_params\": {\"chi_max\": 64, \"svd_min\": 1e-8, \"trunc_cut\": None},\n",
        "        \"max_delta_t\": 3,\n",
        "    },\n",
        ")\n",
        "\n",
        "\n",
        "def identity_factory():\n",
        "    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)\n",
        "\n",
        "\n",
        "mps_initial_state = MPS_neel_state(models[0].lat)\n",
        "\n",
        "print(f\"Computing dynamic coefficients for time={total_time}\")\n",
        "lse_dyn = setup_dynamic_lse(\n",
        "    mpf_trotter_steps,\n",
        "    total_time,\n",
        "    identity_factory,\n",
        "    exact_factory,\n",
        "    approx_factory,\n",
        "    mps_initial_state,\n",
        ")\n",
        "problem, coeffs_dyn = setup_frobenius_problem(lse_dyn)\n",
        "try:\n",
        "    problem.solve()\n",
        "    mpf_dynamic_coeffs = coeffs_dyn.value\n",
        "except Exception as error:\n",
        "    mpf_dynamic_coeffs = np.zeros(len(mpf_trotter_steps))\n",
        "    print(error, \"Calculation Failed\")\n",
        "\n",
        "# -------------------------Step 1 (cont): Build circuits-------------------------\n",
        "mpf_circuits = []\n",
        "for k in mpf_trotter_steps:\n",
        "    circuit = QuantumCircuit(L)\n",
        "    circuit.x([i for i in range(L) if i % 2])\n",
        "    trotter_circ = generate_time_evolution_circuit(\n",
        "        hamiltonian,\n",
        "        synthesis=SuzukiTrotter(reps=k, order=order),\n",
        "        time=total_time,\n",
        "    )\n",
        "    circuit.compose(trotter_circ, qubits=range(L), inplace=True)\n",
        "    mpf_circuits.append(circuit)\n",
        "\n",
        "# Baseline \"single deep circuit\" comparison run with k=10 Trotter steps.\n",
        "# Its two-qubit depth is deeper than the deepest MPF constituent (k_max=6) plus\n",
        "# the overhead of running multiple circuits, pushing it into the noise-limited\n",
        "# regime where MPF is expected to outperform. It does NOT target the MPF's effective\n",
        "# Trotter error (which would require many more steps).\n",
        "comp_circuit = QuantumCircuit(L)\n",
        "comp_circuit.x([i for i in range(L) if i % 2])\n",
        "trotter_circ = generate_time_evolution_circuit(\n",
        "    hamiltonian,\n",
        "    synthesis=SuzukiTrotter(reps=10, order=order),\n",
        "    time=total_time,\n",
        ")\n",
        "comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)\n",
        "mpf_circuits.append(comp_circuit)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b5d388b1",
      "metadata": {},
      "source": [
        "ここで、選択したバックエンドに合わせて回路を最適化します。 `optimization_level=3`ここでは、Qiskitのプリセット・パス・マネージャーを使用しており、これにより適切な物理量子ビットのセットが自動的に選択され、各回路がデバイスのトポロジー上にルーティングされます。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "id": "05bad997",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "<IBMBackend('ibm_fez')>\n"
          ]
        }
      ],
      "source": [
        "# -------------------------Step 2-------------------------\n",
        "service = QiskitRuntimeService()\n",
        "# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=L)\n",
        "backend = service.backend(\"ibm_fez\")\n",
        "print(backend)\n",
        "\n",
        "transpiler = generate_preset_pass_manager(\n",
        "    optimization_level=3, backend=backend\n",
        ")\n",
        "transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]\n",
        "\n",
        "isa_observables = [\n",
        "    observable.apply_layout(circ.layout) for circ in transpiled_circuits\n",
        "]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b578c285",
      "metadata": {},
      "source": [
        "実際のハードウェア上でより複雑な回路を実行するには、徹底的なエラー緩和対策が必要となる。 当手法では、動的デカップリング、ゲートおよび測定のツイリング、測定誤差の低減、およびゼロノイズ外挿（ZNE）を実現します。 なお、ここで用いているZNEノイズ係数（`1, 1.2, 1.4`）は、浅い回路のシナリオの場合よりも小さいことに留意されたい。これは、より深いMPF構成要素がすでにノイズ閾値に近い位置にあり、ノイズの増幅が大きくなると、ZNEによる外挿が信頼できる範囲を超えてしまうためである。\n",
        "\n",
        "4つの回路すべて（ $k_j = [3, 4, 6]$ 上の3つのMPF構成要素と、 $k = 10$ のベースライン）を、1つのEstimatorジョブとして提出します。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "id": "e2722b61",
      "metadata": {},
      "outputs": [],
      "source": [
        "# -------------------------Step 3-------------------------\n",
        "estimator = Estimator(mode=backend)\n",
        "estimator.options.default_shots = 30000\n",
        "\n",
        "# Error suppression/mitigation\n",
        "estimator.options.dynamical_decoupling.enable = True\n",
        "estimator.options.twirling.enable_gates = True\n",
        "estimator.options.twirling.enable_measure = True\n",
        "estimator.options.twirling.num_randomizations = \"auto\"\n",
        "estimator.options.twirling.strategy = \"active-accum\"\n",
        "estimator.options.resilience.measure_mitigation = True\n",
        "estimator.options.experimental.execution_path = \"gen3-turbo\"\n",
        "\n",
        "estimator.options.resilience.zne_mitigation = True\n",
        "estimator.options.resilience.zne.noise_factors = (1, 1.2, 1.4)\n",
        "estimator.options.resilience.zne.extrapolator = \"linear\"\n",
        "\n",
        "estimator.options.environment.job_tags = [\"TUT_MPF\"]\n",
        "\n",
        "job_50 = estimator.run(\n",
        "    [\n",
        "        (circ, observable)\n",
        "        for circ, observable in zip(transpiled_circuits, isa_observables)\n",
        "    ]\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "afc0c029",
      "metadata": {},
      "source": [
        "ジョブ結果から回路ごとの期待値と標準偏差を抽出し、それらを小規模な例とまったく同じように各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$ となります。\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "id": "a924d79c",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[array(-0.07916195), array(-0.04479681), array(-0.2560756), array(-0.06045848)]\n",
            "[array(0.04605538), array(0.10056336), array(0.14426151), array(0.04059092)]\n"
          ]
        }
      ],
      "source": [
        "# -------------------------Step 4-------------------------\n",
        "result = job_50.result()\n",
        "evs = [res.data.evs for res in result]\n",
        "std = [res.data.stds for res in result]\n",
        "\n",
        "print(evs)\n",
        "print(std)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "id": "1071de0d",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Exact static MPF expectation value:  -0.5665938395816946 +- 0.3925273058119915\n",
            "Approximate static MPF expectation value:  -0.25856647611537903 +- 0.164249927266166\n",
            "Dynamic MPF expectation value:  -0.12667812062949296 +- 0.06059471006973169\n"
          ]
        }
      ],
      "source": [
        "exact_mpf_std = np.sqrt(\n",
        "    sum([(coeff**2) * (std**2) for coeff, std in zip(mpf_coeffs, std[:3])])\n",
        ")\n",
        "print(\n",
        "    \"Exact static MPF expectation value: \",\n",
        "    evs[:3] @ mpf_coeffs,\n",
        "    \"+-\",\n",
        "    exact_mpf_std,\n",
        ")\n",
        "approx_mpf_std = np.sqrt(\n",
        "    sum(\n",
        "        [\n",
        "            (coeff**2) * (std**2)\n",
        "            for coeff, std in zip(coeffs_approx.value, std[:3])\n",
        "        ]\n",
        "    )\n",
        ")\n",
        "print(\n",
        "    \"Approximate static MPF expectation value: \",\n",
        "    evs[:3] @ coeffs_approx.value,\n",
        "    \"+-\",\n",
        "    approx_mpf_std,\n",
        ")\n",
        "dynamic_mpf_std = np.sqrt(\n",
        "    sum(\n",
        "        [\n",
        "            (coeff**2) * (std**2)\n",
        "            for coeff, std in zip(mpf_dynamic_coeffs, std[:3])\n",
        "        ]\n",
        "    )\n",
        ")\n",
        "print(\n",
        "    \"Dynamic MPF expectation value: \",\n",
        "    evs[:3] @ mpf_dynamic_coeffs,\n",
        "    \"+-\",\n",
        "    dynamic_mpf_std,\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "id": "64360d85",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/multi-product-formula/extracted-outputs/64360d85-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "sym = {3: \"^\", 4: \"s\", 6: \"p\"}\n",
        "for k, step in enumerate(mpf_trotter_steps):\n",
        "    plt.errorbar(\n",
        "        k,\n",
        "        evs[k],\n",
        "        yerr=std[k],\n",
        "        alpha=0.5,\n",
        "        markersize=4,\n",
        "        marker=sym[step],\n",
        "        color=\"grey\",\n",
        "        label=f\"{mpf_trotter_steps[k]} Trotter steps\",\n",
        "    )\n",
        "\n",
        "plt.errorbar(\n",
        "    3,\n",
        "    evs[-1],\n",
        "    yerr=std[-1],\n",
        "    alpha=0.5,\n",
        "    markersize=8,\n",
        "    marker=\"x\",\n",
        "    color=\"blue\",\n",
        "    label=\"10 Trotter steps\",\n",
        ")\n",
        "\n",
        "plt.errorbar(\n",
        "    4,\n",
        "    evs[:3] @ mpf_coeffs,\n",
        "    yerr=exact_mpf_std,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    color=\"purple\",\n",
        "    label=\"Static MPF\",\n",
        ")\n",
        "\n",
        "plt.errorbar(\n",
        "    5,\n",
        "    evs[:3] @ coeffs_approx.value,\n",
        "    yerr=approx_mpf_std,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    color=\"orange\",\n",
        "    label=\"Approximate static MPF\",\n",
        ")\n",
        "\n",
        "plt.errorbar(\n",
        "    6,\n",
        "    evs[:3] @ mpf_dynamic_coeffs,\n",
        "    yerr=dynamic_mpf_std,\n",
        "    markersize=4,\n",
        "    marker=\"o\",\n",
        "    color=\"pink\",\n",
        "    label=\"Dynamic MPF\",\n",
        ")\n",
        "\n",
        "exact_obs = -0.24384471447172074  # Calculated via Tensor Network calculation\n",
        "plt.axhline(\n",
        "    y=exact_obs, linestyle=\"--\", color=\"red\", label=\"Exact time-evolution\"\n",
        ")\n",
        "\n",
        "plt.title(\n",
        "    f\"$\\\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\\\rangle$ at time {total_time} for the different methods\"\n",
        ")\n",
        "plt.xlabel(\"Method\")\n",
        "plt.ylabel(\"Expectation Value\")\n",
        "plt.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.2), ncol=2)\n",
        "plt.grid(alpha=0.1)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "31effe5f",
      "metadata": {},
      "source": [
        "上記のハードウェアに関する結果について、いくつか指摘しておきます：\n",
        "\n",
        "* **ハードウェアにおいて、より深く掘り下げるにはコストがかかります。** シングル回路のベースラインは、その状況を如実に物語っています。 $k = 6$ 回路は実質的に正確ですが（ $-0.256$ 対 参照値 $-0.244$ ）、一方、より深い $k = 10$ のベースラインは、改善されているどころか、 *むしろ悪化*しています（ $-0.061$、誤差 $\\sim 0.18$ ）。 トロッター誤差がすでに小さい場合、ステップを追加しても、主に回路が深くなるだけで、ゲートノイズやデコヒーレンスがさらに蓄積されることになる。 これこそが、MPFが構築された本来の目的、すなわち、浅い構成要素のみを用いて、深い回路と同等の精度を達成することなのです。\n",
        "\n",
        "* **小ノルムMPFは、ディープ・シングル・サーキットよりも優れた性能を発揮する。** 近似静的MPF（上限 $\\|x\\|_1 \\approx 2$ ）は $-0.259$ となり、基準値から $\\sim 0.015$ の範囲内に収まっており、 $k = 10$ のベースラインよりもはるかに近い値となっている。 ダイナミックMPF（ $-0.127$ ）も、その基準値を余裕で上回っています。 どちらも、浅い $k_j = [3, 4, 6]$ 回路のみを組み合わせているにもかかわらず、深い単一回路では導き出せなかった答えを導き出している。\n",
        "\n",
        "* **係数のノルムは、数学的な最適性よりも重要である。** 正確静的MPFは $\\|x\\|_1 = 4.66$ となり、すべての推定法の中で*最悪の*性能を示す（ $-0.567$、誤差は $0.3$ 以上）。係数のノルムが大きいため、各 $\\langle A \\rangle_{k_j}$ における残留ゲートノイズ、デコヒーレンス、およびZNE誤差がほぼ同じ倍率で増幅され、それによって得られるトロッター誤差の相殺効果を圧倒してしまう。 ノルムに上限を設ける（近似静的ソルバー、 $\\|x\\|_1 \\approx 2$ ）ことで、この過大な影響が解消され、最良の推定値が得られる――たとえその係数が、主たるトロッター誤差を正確に相殺しなくなっても。\n",
        "\n",
        "* **個々の浅いコースであっても、依然として競争力を持つことは可能です。** 唯一の $k = 6$ 構成要素（ $-0.256$ ）は、ここではそれ自体が実質的に正確である――今回の実行では、近似静的MPFよりもわずかに近い結果となっている。 問題は、 *どの*$k$ が「収束しているが、まだノイズ制限を受けていない」という最適な範囲にあるのか、事前にわからないという点にある。また、トロッター収束を保証するために単に深さを増す（ $k = 10$ ）という、一見安全に見える選択こそが、まさに失敗に終わる選択なのである。 MPFは、適切な深さを推測する必要のない、浅い回路の原理に基づいた組み合わせを提供します。\n",
        "\n",
        "実用上のポイントとしては、ハードウェア上では、MPFを個々の $\\langle A \\rangle_{k_j}$ に対する強力な誤差緩和策と組み合わせるべきであり、係数 $L_1$ のノルムは適度な値に保つべきである（近似ソルバーまたは動的MPFを使用する）。また、Trotterステップ $k_j$ は、 $t/k_{\\min} \\lesssim 1$ となるように選択すべきである。ここで、 $t = 3$ の $k_{\\min} = 3$ によると、 $t/k_{\\min} = 1$ となり、静的MPFが依存する先行誤差モデルが有効である収束領域内に構成要素を保つことができる。 これらの選択により、ここでの小ノルムMPFは収束した単一回路と同等の性能を示すのに対し、単純な「深さを増すだけ」というベースラインはそうならず、参考文献 [\\[3\\]](#references) で示された「深さ対精度」の利点が再現される。 また、個々の実行結果にはノイズが含まれる点にも留意してください。同じジョブを別のサブミッション（または別のバックエンド）で実行すると、正確な順位が変動する可能性があります。堅牢な傾向としては、small- $\\|x\\|_1$ のMPFは良好な結果を示し、large- $\\|x\\|_1$ のexact-static MPFはハードウェアノイズの影響を強く受け、over-deepの単一回路はノイズによって性能が制限されることが挙げられます。\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5ac2f8a8",
      "metadata": {},
      "source": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## 次のステップ\n",
        "\n",
        "<Admonition type=\"tip\" title=\"推奨事項\">\n",
        "  この作品に興味を持たれた方は、以下の資料もご参考になるかもしれません：\n",
        "\n",
        "  * [MPF用のトロッターステップの選び方](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html) — 不安定性を回避するための $k_j$ 値の選定に関する実践的な指針\n",
        "  * [近似モデルの使用方法](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html) — 近似静的MPFにおける $L_1$ -ノルム制約およびソルバーオプションの調整\n",
        "  * [`qiskit-addon-mpf` APIリファレンス](https://qiskit.github.io/qiskit-addon-mpf/) — 静的モジュール、動的モジュール、およびバックエンドモジュールに関する完全なドキュメント\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "70be41e1",
      "metadata": {},
      "source": [
        "<span id=\"references\" />\n",
        "\n",
        "## 参照\n",
        "\n",
        "\\[1] バスケス, A. C., エガー, D. J., Ochsner, D., & Woerner, S. ハードウェアに適したハミルトニアンシミュレーションのための、条件の整った多製品公式。 [『Quantum』第7巻、1067頁（2023年）](https://quantum-journal.org/papers/q-2023-07-25-1067/)\n",
        "\n",
        "\\[2] ジュク, S., Robertson, N. F., & Bravyi, S. ハミルトニアンシミュレーションのためのトロッター誤差の上界と動的マルチプロダクト公式。 [『Physical Review Research』, 6(3), 033309 (2024)](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309)\n",
        "\n",
        "\\[3] ロバートソン, N. F., et al. テンソルネットワークによる動的マルチプロダクト公式の拡張。 [arXiv:2407.17405 (2024)](https://arxiv.org/abs/2407.17405)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    },
    "hours": 1.5,
    "qpuSeconds": 240
  },
  "nbformat": 4,
  "nbformat_minor": 5
}