{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "454a9dfd",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Mehrproduktformeln zur Reduzierung von Trotter-Fehlern\"\n",
        "description: \"Verwenden Sie bei der Schätzung von Observablen Mehrproduktformeln, um den Trotter-Fehler zu reduzieren, oder implementieren Sie die Zeitentwicklung mit festem Trotter-Fehler bei geringerer Tiefe.\"\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",
        "# Mehrproduktformeln zur Reduzierung von Trotter-Fehlern\n",
        "\n",
        "*Geschätzte Laufzeit: Vier Minuten auf einem Heron- r2 -Prozessor (HINWEIS: Dies ist lediglich eine Schätzung. (Die Laufzeit kann variieren.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c4d0b2f2",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Lernziele\n",
        "\n",
        "* Wie Multi-Product-Formeln (MPFs) den Trotter-Fehler bei der Hamilton-Simulation durch die Kombination von Erwartungswerten aus mehreren flachen Schaltkreisen reduzieren\n",
        "* Wann MPFs gegenüber Standardproduktformeln Vorteile bieten und wann sie nicht das richtige Mittel der Wahl sind\n",
        "* So berechnen Sie statische und dynamische MPF-Koeffizienten mithilfe des `qiskit_addon_mpf` Pakets\n",
        "* So führen Sie einen MPF-Workflow von Anfang bis Ende auf einer „ IBM Quantum® “-Hardware durch, einschließlich Transpilation, Fehlerbehebung und Nachbearbeitung\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f5dfd316",
      "metadata": {},
      "source": [
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Voraussetzungen\n",
        "\n",
        "* [Kompilierungsmethoden für Hamilton-Simulationsschaltungen](/docs/tutorials/compilation-methods-for-hamiltonian-simulation-circuits) – Einführung in Trotter-Schaltungen (Produktformel) in Qiskit.\n",
        "* Produktformeln in Qiskit, insbesondere die [`SuzukiTrotter`](/docs/api/qiskit/qiskit.synthesis.SuzukiTrotter) Synthese-Klassen und [`LieTrotter`](/docs/api/qiskit/qiskit.synthesis.LieTrotter) .\n",
        "* [Qiskit primitives und die Estimator-Schnittstelle](/docs/guides/primitives).\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "07273b26",
      "metadata": {},
      "source": [
        "<span id=\"background\" />\n",
        "\n",
        "## Hintergrund\n",
        "\n",
        "<span id=\"what-are-multi-product-formulas\" />\n",
        "\n",
        "### Was sind Mehrproduktformeln?\n",
        "\n",
        "Bei der Simulation von Quantensystemen auf einem Quantencomputer besteht eine zentrale Aufgabe darin, den Zeitentwicklungsoperator $e^{-iHt}$ für einen Hamilton-Operator $H$ zu approximieren. Der Standardansatz nutzt *Produktformeln* (PFs), die auch als Trotter-Suzuki-Zerlegungen bekannt sind. Diese zerlegen $H = \\sum_{a=1}^d F_a$ in Terme, deren einzelne Unitarien $e^{-iF_a t}$ effizient zu implementieren sind, und approximieren dann die vollständige Entwicklung als geordnetes Produkt dieser einfacheren Unitarien.\n",
        "\n",
        "Die Produktformel erster Ordnung (Lie-Trotter) lautet:\n",
        "\n",
        "$$\n",
        "S_1(t) := \\prod_{a=1}^d e^{-i F_a t},\n",
        "$$\n",
        "\n",
        "was zu einem quadratischen Fehler führt: $S_1(t) = e^{-iHt} + \\mathcal{O}(t^2)$. Symmetrische Formeln höherer Ordnung $S_{2\\chi}(t)$, wobei $\\chi$ die Ordnung der symmetrischen Produktformel bezeichnet (siehe Ref. [\\[1\\]](#references) ), konvergieren schneller als $e^{-iHt} + \\mathcal{O}(t^{2\\chi+1})$, allerdings auf Kosten tieferer Schaltkreise pro Schritt.\n",
        "\n",
        "Um den Fehler bei einem $\\chi$ *fester* Ordnung zu verringern, wird die gesamte Evolutionszeit $t$ üblicherweise in $k$ kleinere Trotter-Schritte unterteilt. Jeder Schritt approximiert $e^{-iHt/k}$ durch eine Produktformel, und die Schritte werden aneinandergereiht:\n",
        "\n",
        "$$\n",
        "e^{-iHt} \\approx \\left[S_{2\\chi}(t/k)\\right]^k.\n",
        "$$\n",
        "\n",
        "Bei einer symmetrischen Formel d $2\\chi$ -ter Ordnung skaliert der verbleibende Trotter-Fehler dann wie folgt: $\\mathcal{O}\\!\\left(t^{2\\chi+1} / k^{2\\chi}\\right)$. Eine Erhöhung von $k$ führt also zu einer raschen Unterdrückung des Trotter-Fehlers – allerdings wird dadurch auch die Schaltung linear tiefer, was auf rauschbehafteter Hardware zu mehr kumuliertem Gatterrauschen führt. Genau diese Spannung zwischen **dem Trotter-Fehler (begünstigt größere $k$ )** und **dem Hardware-Rauschen (begünstigt kleinere $k$ )** soll durch Multiproduktformeln aufgelöst werden. Beachten Sie, dass es bei MPFs darum geht, Ergebnisse aus *verschiedenen Auswahlen von „ $k$ “* in einer festgelegten Reihenfolge $\\chi$ zu kombinieren – sie ändern die Reihenfolge der zugrunde liegenden Produktformel nicht.\n",
        "\n",
        "**Multi-Produkt-Formeln (MPFs)** [\\[1\\]](#references) bilden eine *gewichtete lineare Kombination* von Erwartungswerten, die aus mehreren flacheren Trotter-Schaltkreisen gewonnen werden, wobei jeder eine unterschiedliche Anzahl von Trotter-Schritten verwendet $k_1, k_2, \\ldots, k_r$ (eine Reihe von $r$ Schrittzahlen):\n",
        "\n",
        "$$\n",
        "\\langle A \\rangle_{\\text{MPF}}(t) = \\sum_{j=1}^r x_j \\, \\langle A \\rangle_{k_j}(t),\n",
        "$$\n",
        "\n",
        "wobei $\\langle A \\rangle_{k_j}(t)$ der Erwartungswert einer Observablen $A$ zum Zeitpunkt $t$ ist, der anhand eines Trotter-Kreises mit $k_j$ Schritten geschätzt wurde, und die Koeffizienten $\\{x_j\\}_{j=1}^r$ so gewählt sind, dass sich die führenden Trotter-Fehlerterme in der Kombination aufheben. Wir werden in [Schritt 4](#small-scale-step-4) auf diesen Ausdruck zurückkommen, wo wir ihn explizit auswerten, um unsere Trotter-Ergebnisse zusammenzufassen. Der entscheidende praktische Aspekt ist, dass die tiefste Schaltung im MPF nur $k_{\\max}$ Schritte benötigt, was deutlich weniger ist als der einzelne Schritt $k$, der erforderlich wäre, um denselben effektiven Trotter-Fehler direkt zu erreichen. Durch die flacheren Schaltkreise eignet sich der MPF-Ansatz besser für rauschbehaftete Hardware.\n",
        "\n",
        "<span id=\"how-are-the-coefficients-determined\" />\n",
        "\n",
        "### Wie werden die Koeffizienten ermittelt?\n",
        "\n",
        "Es gibt zwei Familien von MPF-Koeffizienten:\n",
        "\n",
        "**Statische Koeffizienten** sind unabhängig vom Hamilton-Operator, vom Anfangszustand und von der Entwicklungszeit. Sie werden durch die Lösung eines linearen Gleichungssystems $Ax = b$ ermittelt, das die Aufhebung der führenden Trotter-Fehlerterme gewährleistet. Für eine Reihe von Trotter-Schritten $\\{k_j\\}_{j=1}^r$ in Verbindung mit einer symmetrischen Produktformel der Ordnung $2\\chi$ führt die Entwicklung des Trotter-Fehlers in umgekehrten Potenzen von $k_j$ zu Nebenbedingungen der Form:\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",
        "wobei die ganzzahligen Exponenten $\\{\\eta_n\\}$ die Ordnungen der aufeinanderfolgenden Trotter-Fehlerterme für die gewählte Produktformel sind. Für eine *symmetrische* PF der Ordnung $2\\chi$ gilt, dass der Hauptfehler in $\\left[S_{2\\chi}(t/k)\\right]^k$ wie folgt skaliert: $1/k^{2\\chi}$, mit nachfolgenden Korrekturen gemäß $1/k^{2\\chi+2}, 1/k^{2\\chi+4}, \\ldots$ – die Exponenten lauten also $\\eta_n = 2\\chi + 2n$. Bei nicht-symmetrischen PFs tragen sowohl ungerade als auch gerade Potenzen bei, und es gilt $\\eta_n = 2\\chi + n$. Die vollständige Herleitung findet sich in Ref. [\\[1\\]](#references). Die erste Gleichung im obigen System gewährleistet die Unverzerrtheit (die MPF gibt den exakten Erwartungswert im Grenzfall „ $k_j \\to \\infty$ “ wieder), und die übrigen $r-1$ Gleichungen heben nacheinander die ersten $r-1$ Trotter-Fehlerterme auf. Wenn die resultierende $L_1$ -Norm $\\|x\\|_1$ zu groß ist (was das Abtastrauschen verstärkt), können Sie stattdessen eine approximative Optimierung lösen, die $\\|x\\|_1$ begrenzt und gleichzeitig $\\|Ax - b\\|$ minimiert.\n",
        "\n",
        "**Die dynamischen Koeffizienten** [\\[2\\]](#references), [\\[3\\]](#references) hängen zudem vom Hamilton-Operator, vom Anfangszustand und von der Entwicklungszeit $t$ ab. Sie minimieren den Abstand in der Frobenius-Norm zwischen dem tatsächlichen, zeitlich entwickelten Zustand und der MPF-Näherung:\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",
        "wobei $M_{ij}(t) = \\mathrm{Tr}[\\rho_{k_i}(t)\\,\\rho_{k_j}(t)]$ die Gram-Matrix der Überlappungen zwischen den durch Trotter-Evolution erzeugten Zuständen für verschiedene Schrittzahlen $k_i, k_j$ ist und $L_i(t) = \\mathrm{Tr}[\\rho(t)\\,\\rho_{k_i}(t)]$ die Überlappung mit dem (annähernd) exakten Zustand misst. In diesem Tutorial werden diese Größen mithilfe von Tensor-Netzwerk-Methoden effizient berechnet, insbesondere mithilfe der Backends „ TeNPy-based “ in `qiskit_addon_mpf`.\n",
        "\n",
        "<span id=\"when-to-use-mpfs\" />\n",
        "\n",
        "### Wann sollten MPFs eingesetzt werden?\n",
        "\n",
        "MPFs sind am vorteilhaftesten, wenn:\n",
        "\n",
        "* **Die Schaltungstiefe ist der Engpass.** Wenn die durch die Hardware bedingten Rauschgrenzen die mögliche Tiefe einschränken, verwenden Sie MPFs, um bei flacheren Schaltungen eine höhere effektive Trotter-Genauigkeit zu erzielen.\n",
        "* **Sie benötigen genaue Erwartungswerte, keine vollständige Zustandsvorbereitung.** MPFs arbeiten auf der Ebene der Erwartungswerte – sie kombinieren klassische Zahlen, keine Quantenzustände. Sie eignen sich daher ideal für die Schätzung anhand von Beobachtungsdaten bei Verwendung der Estimator-Primitive.\n",
        "* **Du kombinierst eine überschaubare Anzahl von Trotter-Schritten.** In der Regel reicht es aus, $r = 3$ – $5$ e unterschiedliche Schrittzahlen $k_j$ zu kombinieren, um mehrere führende Trotter-Fehlerterme zu eliminieren und gleichzeitig $\\|x\\|_1$ überschaubar zu halten.\n",
        "\n",
        "<span id=\"when-mpfs-might-not-help\" />\n",
        "\n",
        "### Wann MPFs möglicherweise nicht helfen\n",
        "\n",
        "* **Sehr kurze Entwicklungszeiten.** Wenn $t$ so klein ist, dass bereits eine einzige Trotter-Formel niedriger Ordnung ausreichend genau ist, entfällt der Mehraufwand für die Ausführung mehrerer Durchläufe.\n",
        "* **Aufgaben zur Vorbereitung auf die Prüfung.** MPFs erzeugen einen korrigierten *Erwartungswert*, keinen korrigierten Quantenzustand. Wenn Sie den aktuellen, zeitlich entwickelten Zustand benötigen (beispielsweise als Eingabe für eine andere Quanten-Subroutine), kommen MPFs nicht in Frage.\n",
        "* **Trotter-Schrittzahlen, die gegen das Konvergenzregime verstoßen.** Bei der Herleitung des statischen Koeffizienten wird jeder einzelne Term $\\left[S_{2\\chi}(t/k_j)\\right]^{k_j}$ als Reihe in $t/k_j$ entwickelt; diese Entwicklung konvergiert nur dann gut, wenn $t/k_{\\min} \\lesssim 1$ gilt. Wird $k_{\\min}$ für die gegebenen $t$ zu klein gewählt, liegt der flachste Schaltkreis weit außerhalb des perturbativen Bereichs, die Fehlerterme höherer Ordnung, die durch die MPF nicht aufgehoben werden, werden groß, und die Aufhebung kann große Koeffizienten erfordern. Die „ $L_1$ “-Norm $\\|x\\|_1$ ist ein praktisches Diagnosekriterium: Wenn $\\|x\\|_1 \\gg 1$ gilt, könnte der Overhead durch die Abtastung $\\propto \\|x\\|_1^2$ die Verringerung des Trotter-Fehlers überwiegen. Weitere Informationen finden Sie im [Leitfaden zur Auswahl von Trotter-Stufen](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html).\n",
        "\n",
        "<span id=\"what-this-tutorial-covers\" />\n",
        "\n",
        "### Inhalt dieses Tutorials\n",
        "\n",
        "Dieses Tutorial führt Sie in zwei Schritten durch einen MPF-Workflow von Anfang bis Ende. Zunächst zeigt **ein Beispiel mit einem Simulator im kleinen Maßstab** (10-Qubit-Heisenberg-Kette), wie man das Problem aufstellt, statische und dynamische MPF-Koeffizienten berechnet und die resultierenden Erwartungswerte mit der exakten Diagonalisierung vergleicht. Anschließend zeigt ein **groß angelegtes Hardware-Beispiel** (50-Qubit-XXZ-Kette), wie man den Code transpiliert, auf einer „ IBM Quantum “-Hardware mit Fehlerminderung ausführt und die Ergebnisse mithilfe der MPF-Koeffizienten nachbearbeitet. Wir verwenden das `qiskit_addon_mpf` Paket durchgehend zusammen mit den Standard-Tools von Qiskit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b5d478ce",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## Anforderungen\n",
        "\n",
        "Vergewissern Sie sich vor Beginn dieses Tutorials, dass Sie die folgenden Komponenten installiert haben:\n",
        "\n",
        "* Qiskit SDK v2.0 oder höher, mit [Visual](/docs/api/qiskit/visualization) isierungsunterstützung\n",
        "* Qiskit Runtime v0.22 oder später (`pip install qiskit-ibm-runtime`)\n",
        "* Qiskit Aer-Simulator (`pip install qiskit-aer`)\n",
        "* MPF-Qiskit-Add-on mit dem „ TeNPy “-Backend (`pip install \"qiskit-addon-mpf[tenpy]\"`)\n",
        "* Qiskit-Add-on-Dienstprogramme (`pip install qiskit-addon-utils`)\n",
        "* SciPy (`pip install scipy`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2584c37e",
      "metadata": {},
      "source": [
        "<span id=\"setup\" />\n",
        "\n",
        "## Einrichtung\n",
        "\n",
        "Im Folgenden fassen wir *alle* in diesem Tutorial verwendeten Paketimporte in einer einzigen Zelle zusammen. `XXPlusYYGate`Außerdem definieren wir einen `CollectAndCollapse` Transpiler-Durchlauf, der benachbarte `rxx` und `ryy` Rotationen zu einer einzigen zusammenfasst. Dieser Durchlauf wird sowohl beim Aufbau der Schaltung in Schritt 1 (um die Anzahl der Gatter gering zu halten) als auch indirekt bei der Extraktion der Schichtstruktur für das dynamische MPF in Schritt 4 angewendet (die Funktion „ TeNPy “ erwartet Zwei-Qubit-Gatter und keine Paare von nicht fusionierten Rotationen).\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",
        "## Beispiel für einen Simulator im kleinen Maßstab\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "378e82ba",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### Schritt 1: Klassische Eingaben einem Quantenproblem zuordnen\n",
        "\n",
        "Wir beginnen mit einem 10-Qubit-Heisenberg-Modell auf einer Linie und verwenden dabei die Néel-Zustands $\\vert 0101\\ldots01 \\rangle$ als Anfangszustand. Der Hamilton-Operator lautet:\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",
        "wobei $J$ die Kopplungsstärke zwischen den nächsten Nachbarn ist. Wir messen den ZZ-Korrelator $Z_{L/2-1} Z_{L/2}$ an einem Qubit-Paar in der Mitte der Kette und verwenden Trotter-Schritte $k_j = [1, 2, 4]$ mit einer Produktformel zweiter Ordnung.\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",
        "#### Trotter-Schaltungen aufbauen\n",
        "\n",
        "Wir erstellen die Schaltungen, die die approximativen Trotter-Zeitentwicklungen für jeden Zeitpunkt und jede Anzahl von Trotter-Schritten umsetzen. Der `CollectAndCollapse` im Abschnitt „Setup“ definierte Durchlauf fasst XX- und YY-Rotationen zu einzelnen XX+YY-Gates zusammen, um später eine effizientere Tensor-Netzwerk-Simulation zu ermöglichen.\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",
        "### Schritt 2: Optimierung des Problems für die Ausführung auf Quantenhardware\n",
        "\n",
        "Für das Beispiel im kleinen Maßstab nehmen wir den Aer-Simulator als Ziel. Bevor die Schaltkreise zur Ausführung bereit sind, finden zwei Umwandlungen statt:\n",
        "\n",
        "1. **Gate-Erfassung auf der Ebene der Hamilton-Simulation.** `XXPlusYYGate`In der Setup-Zelle haben wir einen `CollectAndCollapse` Durchgang erstellt, `ryy` der benachbarte `rxx` Rotationen zu einer einzigen zusammenführt. Diesen Durchlauf haben wir bereits bei der Erstellung der Trotter-Schaltungen in Schritt 1 (dem `pm.run(...)` Aufruf) angewendet. Dadurch wird sowohl die Anzahl der Zwei-Qubit-Gatter reduziert als auch eine Struktur geschaffen, die sich besser für die Tensor-Netzwerk-Simulation zur späteren Berechnung der dynamischen Koeffizienten eignet.\n",
        "\n",
        "2. **Anpassung an die ISA des Simulators.** Im Folgenden führen wir den Qiskit-Preset-Pass-Manager aus, `optimization_level=3` um jede Trotter-Schaltung an die Befehlssatzarchitektur (ISA) des Simulators anzupassen.\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",
        "### Schritt 3: Mit „ Qiskit primitives “ ausführen\n",
        "\n",
        "Für das Beispiel in kleinem Maßstab lassen wir die ISA-optimierten Trotter-Schaltkreise durch die `EstimatorV2` von Aer unterstützte Primitive laufen. Auf diese Weise erhalten wir für jedes $(k_j, t)$ -Paar einen *rauschfreien* Referenzwert – dies sind die $\\langle A \\rangle_{k_j}(t)$ -Werte, die der MPF in Schritt 4 kombinieren wird. Wir geben einen Überblick über die Entwicklungszeiträume, damit wir später die vollständige Zeitreihenkurve jeder einzelnen Produktformel und der MPF darstellen können.\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",
        "### Schritt 4: Nachbearbeitung und Rückgabe des Ergebnisses im gewünschten klassischen Format\n",
        "\n",
        "In Schritt 4 wird das MPF tatsächlich erstellt. Auch wenn die Koeffizienten $x_j$ hier *berechnet* werden (und diese Berechnung bei der dynamischen Variante sehr rechenintensiv sein kann), stellen sie konzeptionell ein klassisches Verfahren dar, um die Quantenmessungen aus Schritt 3 zu einem einzigen korrigierten Erwartungswert zu kombinieren – daher betrachten wir den gesamten Arbeitsablauf rund um die Koeffizienten und die Kombination als Nachbearbeitung.\n",
        "\n",
        "Um zu beurteilen, wie gut das MPF-Modell die tatsächliche Dynamik abbildet, berechnen wir zunächst die exakten, zeitabhängigen Erwartungswerte, indem wir den Hamilton-Operator direkt potenzieren. Dies ist nur deshalb machbar, weil $L = 10$ gilt; im folgenden Beispiel für groß angelegte Hardware müssen wir uns stattdessen auf Schätzungen mittels Tensor-Netzwerken stützen.\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",
        "#### Statische MPF-Koeffizienten\n",
        "\n",
        "Statische MPFs verwenden Koeffizienten $x_j$, die unabhängig von der Evolutionszeit, dem Hamilton-Operator und dem Anfangszustand sind. Wir stellen das im Hintergrund beschriebene lineare System $Ax = b$ auf und bestimmen die Koeffizienten. Die Matrix $A$ wird durch die Trotter-Schrittzahlen $k_j$, die Ordnung $\\chi$ der Produktformel und die Frage bestimmt, ob die Formel symmetrisch ist (was die Exponenten $\\eta_n$ bestimmt).\n",
        "\n",
        "Für unser einfaches Beispiel verwenden wir $k_j = [1, 2, 4]$ mit einer asymmetrischen Suzuki-Trotter-Formel der Ordnung $2\\chi=2$ (also $\\chi=1$ und $\\eta_n = 2 + n$, was $\\eta_0 = 2,\\, \\eta_1 = 3$ ergibt). Das System lautet dann:\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",
        "Die erste Zeile gewährleistet die Unvoreingenommenheit ( $\\sum_j x_j = 1$ ); die zweite und dritte Zeile heben jeweils die führenden $1/k^2$ und die Trotter-Fehlerterme nächster Ordnung $1/k^3$ auf.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4f2ca1e2",
      "metadata": {},
      "source": [
        "<span id=\"set-up-the-lse\" />\n",
        "\n",
        "##### Richten Sie die LSE ein\n",
        "\n",
        "Wir verwenden `setup_static_lse` von `qiskit_addon_mpf.static` bis, um die Matrix $A$ und den Vektor auf der rechten Seite $b$, wie oben beschrieben, zu bilden. Die Matrix $A$ hängt nicht nur von $k_j$ ab, sondern auch von der Wahl der Produktformel – insbesondere von deren *Ordnung* $\\chi$ und davon, ob sie *symmetrisch* ist. Das `symmetric` Flag steuert das Exponentenmuster $\\eta_n$ (symmetrische Formeln erzeugen ausschließlich Trotter-Fehlerterme mit gerader Potenz; siehe Ref. [\\[1\\]](#references) ). Beachten Sie, dass, wie in [\\[2\\]](#references) gezeigt, die Festlegung `symmetric=True` von nicht unbedingt erforderlich ist, selbst wenn die zugrunde liegende PF symmetrisch ist – die nicht-symmetrische LSE bleibt gültig (sie erzwingt zusätzliche, unnötige Einschränkungen).\n",
        "\n",
        "Für unser Beispiel haben wir und `symmetric = False` bereits in Schritt 1 festgelegt `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": [
        "Überprüfen Sie die erstellte Matrix $A$ und den Vektor $b$, um sicherzustellen, dass sie mit dem oben angegebenen System übereinstimmen.\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": [
        "Mit der LSE in der Hand bestimmen wir die statischen Koeffizienten $x_j$ über `lse.solve()` (dies ist die direkte Lösung nach $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",
        "##### Optimieren Sie für „ $x$ “ mithilfe eines exakten Modells\n",
        "\n",
        "Alternativ zur Berechnung von $x = A^{-1}b$ können Sie [`setup\\_exact\\_model`](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_exact_model.html) verwenden, um eine [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem)-Instanz zu erstellen, die die LSE als Nebenbedingungen verwendet und deren optimale Lösung $x$ ergibt.\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",
        "##### Optimieren Sie für „ $x$ “ mithilfe eines Näherungsmodells\n",
        "\n",
        "Es kann vorkommen, dass die Norm „ $L_1$ “ für die gewählte Menge von „ $k_j$ “-Werten als zu hoch angesehen wird. Wenn dies der Fall ist und Sie keine anderen Werte für „ $k_j$ “ auswählen können, können Sie eine Näherungslösung verwenden, die die „ $L_1$ “-Norm auf einen gewählten Schwellenwert begrenzt und gleichzeitig „ $\\|Ax - b\\|$ “ minimiert. Lesen Sie dazu die Anleitung [zur](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html) Verwendung des Näherungsmodells.\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",
        "#### Dynamische MPF-Koeffizienten\n",
        "\n",
        "Die statische MPF eliminiert Trotter-Fehlerterme auf eine vom Hamilton-Operator und vom Zustand unabhängige Weise, sodass sie für einen gegebenen Hamilton-Operator und einen gegebenen Anfangszustand nicht zwangsläufig den kleinstmöglichen Approximationsfehler erzeugt. Das dynamische MPF-Verfahren (Lit. [\\[2\\]](#references), [\\[3\\]](#references) ) ermittelt stattdessen zeitabhängige Koeffizienten $x_i(t)$, die den Abstand in der Frobenius-Norm $\\|\\rho(t) - \\mu^D(t)\\|_F^2$ zu jedem Zeitpunkt $t$ minimieren. Wie im Abschnitt „Hintergrund“ dargelegt, erfordert dies die Überlappungsmatrix $M_{ij}(t)$ zwischen den nach Trotter entwickelten Zuständen und die Überlappung $L_i(t)$ mit dem exakten Zustand – beides schätzen wir mithilfe von Tensor-Network-Backends ( TeNPy ) in `qiskit_addon_mpf`.\n",
        "\n",
        "Um die dynamische LSE einzurichten, benötigen wir drei Komponenten:\n",
        "\n",
        "1. **Eine** „Approximate Evolver Factory“, die das Add-on für jeden „ $k_j$ “ ausführt, um „ $\\rho_{k_j}(t)$ “ als MPS/MPO zu erzeugen. Wir konstruieren ihn ausgehend von der Schichtstruktur der Ordnung $2$ Trotter-Schaltung (eine Schicht pro `slice_by_depth`), eingekapselt als mit den `LayerwiseEvolver` Trunkierungsparametern TeNPy.\n",
        "2. Eine **exakte Evolver-Fabrik**, die ein hochgenaues Referenz $\\rho(t)$ erzeugt. Wir verwenden eine Suzuki-Trotter-Schaltung vierter Ordnung mit kleinem Zeitschritt (`dt=0.1`, `order=4`) als Näherung für die exakte Evolution.\n",
        "3. Eine **Identitätsfabrik** und ein **MPS im Anfangszustand,** die als Startdaten für die Simulation „ TeNPy “ dienen.\n",
        "\n",
        "Die folgende Zelle erstellt die Fabrik für die approximativen Evolver.\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",
        "  Die Optionen von `LayerwiseEvolver` , die die Details der Tensornetzsimulation bestimmen, müssen sorgfältig ausgewählt werden, um zu vermeiden, dass ein schlecht definiertes Optimierungsproblem entsteht.\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1c702f67",
      "metadata": {},
      "source": [
        "`dt=0.1`Wir approximieren den exakten zeitabhängigen Zustand mit einer Suzuki-Trotter-Formel vierter Ordnung unter Verwendung eines kleinen Zeitschritts. Die Trunkierungsparameter von „ TeNPy “ können die Genauigkeit beeinflussen, daher ist es wichtig, verschiedene Werte zu testen.\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": [
        "Abschließend definieren wir ein `identity_factory` , das den anfänglichen MPO-Zustand liefert, und bereiten den Néel-Ausgangszustand als MPS vor, der mit dem vom geschichteten Trotter-Modell verwendeten Gitter übereinstimmt.\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": [
        "Nachdem die Fabriken eingerichtet sind, berechnen wir nun die dynamischen Koeffizienten zu jedem Evolutionszeitpunkt. Für jeden $t$ erstellt das Programm die entsprechenden Überlappungsmatrizen über TeNPy, und `setup_frobenius_problem` gibt einen `cvxpy.Problem` Wert zurück, `setup_dynamic_lse` der die Kosten der Frobenius-Norm minimiert. Der Solver gibt auf diesen Zeitpunkt zugeschnittene Koeffizienten $x_j(t)$ zurück; wir sammeln sie in `mpf_dynamic_coeffs_list`. Wenn der Solver bei einem bestimmten „ $t$ “ fehlschlägt, greifen wir auf Nullkoeffizienten zurück, sodass die Schleife fortgesetzt wird.\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",
        "#### Kombinieren Sie die Erwartungswerte nach Trotter mit den MPF-Koeffizienten\n",
        "\n",
        "Nun berechnen wir für jeden Satz von Koeffizienten (static-exact, static-approximate und dynamic) den Wert von „ $\\langle A \\rangle_{\\text{MPF}}(t) = \\sum_j x_j \\, \\langle A \\rangle_{k_j}(t)$ “, übertragen die Standardfehler pro Schaltung und stellen die resultierenden Zeitreihen gegen die Kurve der exakten Diagonalisierung grafisch dar.\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": [
        "Die obige Grafik veranschaulicht das Zusammenspiel zwischen Trotter-Fehler und Stichprobenfehler.\n",
        "\n",
        "* **Trotter-Fehler.** Die einzelnen Produktformeln (graue Markierungen) weichen mit zunehmender Zeit immer stärker von der exakten Kurve ab. Der Regelkreis „ $k=1$ “ weist die größte Abweichung auf und ist am flachsten, befindet sich jedoch bereits in dem Bereich, in dem $t/k \\gtrsim 1$ gilt, sodass der führende Fehlerterm $1/k^{2}$ groß ist. Die MPF-Kombinationen (farbige Markierungen) heben mehrere dieser führenden Trotter-Fehlerterme auf, sodass sie der exakten Kurve viel genauer folgen als jede einzelne „ $k_j$ “-Schaltung. Die verbleibende Lücke spiegelt die Trotter-Terme höherer Ordnung wider, die durch die MPF *nicht* eliminiert werden: Ein statisches MPF der Ordnung $2$, $r=3$ eliminiert lediglich die ersten beiden Fehlerordnungen, und bei großen $t/k_{\\min}$ dominiert schließlich der nicht eliminierte Rest – daher garantiert die MPF nicht, dass sehr flache Schaltkreise zu beliebigen Zeitpunkten genau bleiben.\n",
        "\n",
        "* **Stichprobenfehler.** Die breiteren Fehlerbalken bei den MPF-Kurven sind eine direkte Folge der linearen Kombination: Die Übertragung unabhängiger Standardfehler pro Schaltung $\\sigma_{k_j}$ ergibt eine Gesamtvarianz $\\sigma_{\\text{MPF}}^2 = \\sum_j x_j^2 \\, \\sigma_{k_j}^2$. Je größer also der „ $\\|x\\|_2$ “ (und in der Praxis der „ $\\|x\\|_1$ “, den wir steuern) ist, desto mehr Messungen sind erforderlich, um eine bestimmte Zielunsicherheit zu erreichen. Dies ist der Kompromiss, der hinter der Option „approximate-solver“ im Hintergrund steht: Wir begrenzen „ $\\|x\\|_1$ “, um diesen Mehraufwand überschaubar zu halten. Entscheidend ist, dass der Stichprobenfehler – anders als der Trotter-Fehler – mit $1/\\sqrt{N_{\\text{shots}}}$ abnimmt, sodass er durch den Einsatz weiterer Schüsse stets verringert werden kann.\n",
        "\n",
        "Im folgenden Beispiel für eine groß angelegte Hardware-Implementierung tritt Hardware-Rauschen als zusätzliche Fehlerquelle bei jedem „ $\\langle A \\rangle_{k_j}$ “ auf, das ebenfalls durch die MPF-Koeffizienten verstärkt wird. In diesem Abschnitt werden wir sehen, wie die Fehlerminderung mit MPFs zusammenwirkt.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6fa763ff",
      "metadata": {},
      "source": [
        "<span id=\"large-scale-hardware-example\" />\n",
        "\n",
        "## Beispiel für groß angelegte Hardware\n",
        "\n",
        "In diesem Abschnitt erweitern wir das Problem über das hinaus, was sich exakt simulieren lässt. Wir reproduzieren einige der in Ref. [\\[3\\]](#references) dargestellten Ergebnisse unter Verwendung einer 50-Qubit-XXZ-Kette zum Zeit $t = 3$. Wir folgen demselben vierstufigen Arbeitsablauf wie im Beispiel im kleinen Maßstab, richten uns nun jedoch auf echte Quantenhardware mit Fehlerminderung aus. Wie in der Vorlage ist jeder Schritt direkt im Code gekennzeichnet, und ein einzelner Schritt kann sich über mehrere Zellen erstrecken, wenn seine Zwischenergebnisse eine genauere Betrachtung wert sind.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "481fea70",
      "metadata": {},
      "source": [
        "Die Abbildung entspricht dem Beispiel im kleinen Maßstab: Definieren Sie einen Hamilton-Operator, wählen Sie Trotter-Parameter aus, berechnen Sie die MPF-Koeffizienten (statisch und dynamisch) und bauen Sie Schaltungen auf. Die wichtigsten Unterschiede sind:\n",
        "\n",
        "* Ein **XXZ-Hamilton-Operator** an 50 Knoten mit zufälligen Kopplungen, entnommen aus $\\mathcal{U}(0.5, 1.5)$ (Ref. [\\[3\\]](#references) ).\n",
        "* Eine **symmetrische** Trotter-Formel zweiter Ordnung mit $k_j = [3, 4, 6]$ (also $\\chi=1$, `symmetric=True`).\n",
        "* Ein einziger fester Evolutionszeit $t = 3$. Mit $k_{\\min}=3$ ergibt sich daraus $t/k_{\\min}=1$, wodurch die flachen Anteile im Trotter-Konvergenzbereich verbleiben, in dem das vom MPF verwendete Modell der führenden Fehlerquelle gültig ist.\n",
        "* Ein zusätzlicher **Einzelkreisvergleichslauf mit $k = 10$ Trotter-Schritten**, der als Basislinie verwendet wird. Wir haben uns für „ $k = 10$ “ entschieden, da dessen 2-Qubit-Tiefe auf der Hardware größer ist als die des tiefsten MPF-Bestandteils ( $k_{\\max}=6$ ) zuzüglich des Overheads für die Ausführung mehrerer MPF-Schaltungen – tief genug, um rauschbegrenzt zu sein. In diesem Bereich wird erwartet, dass die MPF-Kombination die Basislinie mit einer einzelnen Schaltung übertrifft. Es handelt sich um einen Vergleich mit einem „einzelnen tiefen Schaltkreis“ im Vergleich zur MPF-Kombination, nicht um einen Schaltkreis, der auf den effektiven Trotter-Fehler des MPF abzielt (was wesentlich mehr Schritte erfordern würde).\n",
        "\n",
        "Beachten Sie, dass wir hier zwar noch in Schritt 1 sind (Abbildung und Schaltungsaufbau), in dieser Zelle jedoch neben den statischen auch die dynamischen Koeffizienten vorab berechnen. Die dynamischen Koeffizienten hängen von $H$ und $t$ ab, jedoch nicht von den Quantenmessungen; daher können sie jederzeit vor Schritt 4 berechnet werden. Wir machen das jetzt, um alle MPF-spezifischen Einstellungen an einem Ort zu bündeln.\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": [
        "Nun optimieren wir die Schaltungen für das ausgewählte Backend. Wir verwenden den Qiskit-Preset-Pass-Manager unter `optimization_level=3`, der automatisch eine geeignete Auswahl physikalischer Qubits auswählt und jede Schaltung auf die Gerätetopologie abbildet.\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": [
        "Um komplexere Schaltkreise auf echter Hardware auszuführen, sind umfassende Maßnahmen zur Fehlerbehebung erforderlich. Wir ermöglichen dynamische Entkopplung, Gate- und Mess-Twirling, Messfehlerminderung sowie rauschfreie Extrapolation (ZNE). Beachten Sie, dass die hier verwendeten ZNE-Rauschfaktoren (`1, 1.2, 1.4`) kleiner sind als in einem Szenario mit flachen Schaltkreisen, da die tiefer liegenden MPF-Bestandteile bereits nahe an der Rauschschwelle liegen und starke Rauschverstärkungen sie über den Punkt hinausdrücken würden, ab dem eine ZNE-Extrapolation noch zuverlässig ist.\n",
        "\n",
        "Wir reichen alle vier Schaltungen (drei MPF-Komponenten unter $k_j = [3, 4, 6]$ sowie die Basislinie „ $k = 10$ “) in einem einzigen Estimator-Auftrag ein.\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": [
        "Wir entnehmen dem Job-Ergebnis die Erwartungswerte und Standardabweichungen pro Schaltkreis und kombinieren diese dann mit den jeweiligen MPF-Koeffizienten genau wie im Beispiel für den kleinen Maßstab: $\\langle A \\rangle_{\\text{MPF}} = \\sum_j x_j \\, \\langle A \\rangle_{k_j}$, mit fortgeführter Varianz $\\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": [
        "Einige Anmerkungen zu den oben aufgeführten Hardware-Ergebnissen:\n",
        "\n",
        "* **Auf der Hardware ist eine tiefergehende Analyse nicht kostenlos.** Die Basislinien der Einzelschaltungen sprechen eine deutliche Sprache: Die Schaltung „ $k = 6$ “ ist im Wesentlichen exakt ( $-0.256$ im Vergleich zur Referenz $-0.244$ ), während die Basislinie der tieferen Schaltung „ $k = 10$\\* “ schlechter\\* ist ( $-0.061$, Abweichung von $\\sim 0.18$ ) und nicht besser. Sobald der Trotter-Fehler bereits gering ist, führt das Hinzufügen weiterer Schritte meist zu einer Vertiefung der Schaltung und zu einer Zunahme von Gate-Rauschen und Dekohärenz. Genau für diesen Anwendungsfall sind MPFs konzipiert: die Genauigkeit einer tiefen Schaltung allein mit flachen Bausteinen zu erreichen.\n",
        "\n",
        "* **Ein MPF mit kleiner Norm schlägt den Deep Single Circuit.** Der annähernd statische MPF (begrenzt auf $\\|x\\|_1 \\approx 2$ ) liegt bei $-0.259$, nur $\\sim 0.015$ vom Referenzwert entfernt und damit weit näher als die $k = 10$ -Basislinie. Auch der dynamische MPF ( $-0.127$ ) übertrifft diese Basislinie deutlich. Beide kombinieren zwar nur die flachen „ $k_j = [3, 4, 6]$ “-Schaltkreise, liefern jedoch eine Antwort, die der tiefe Einzelschaltkreis nicht liefern konnte.\n",
        "\n",
        "* **Die Koeffizienten-Norm ist wichtiger als die mathematische Optimalität.** Das exakt-statische MPF weist eine Norm von $\\|x\\|_1 = 4.66$ auf und ist der *schlechteste* Schätzer von allen ( $-0.567$, mit einer Abweichung von mehr als $0.3$ ): Die große Koeffizienten-Norm verstärkt das Rest-Gate-Rauschen, die Dekohärenz und den ZNE-Fehler bei jedem $\\langle A \\rangle_{k_j}$ um ungefähr denselben Faktor, wodurch die dadurch erzielte Aufhebung des Trotter-Fehlers zunichte gemacht wird. Durch die Begrenzung der Norm (des approximativ-statischen Lösers, $\\|x\\|_1 \\approx 2$ ) wird diese Überlastung beseitigt und die beste Schätzung erzielt – auch wenn sich dessen Koeffizienten nicht mehr exakt mit dem führenden Trotter-Fehler aufheben.\n",
        "\n",
        "* **Auch einzelne Flachwasser-Wettkämpfe können nach wie vor wettbewerbsfähig sein.** Die einzige Komponente „ $k = 6$ “ ( $-0.256$ ) ist hier im Grunde genommen exakt – bei diesem Durchlauf liegt sie sogar geringfügig näher als der annähernd statische MPF. Der Haken daran ist, dass man im Voraus nicht weiß\\*, welcher\\* einzelne $k$ genau den „Sweet Spot“ zwischen „konvergiert, aber noch nicht rauschbegrenzt“ trifft, und dass die auf den ersten Blick sichere Entscheidung, einfach tiefer zu gehen ( $k = 10$ ), um die Trotter-Konvergenz zu gewährleisten, genau die ist, die fehlschlägt. Das MPF bietet eine prinzipielle Kombination aus flachen Schaltkreisen, bei der die richtige Tiefe nicht erraten werden muss.\n",
        "\n",
        "Die praktische Schlussfolgerung lautet, dass MPFs auf der Hardware mit einer starken Fehlerminderung für jeden einzelnen „ $\\langle A \\rangle_{k_j}$ “ kombiniert werden sollten, die Koeffizienten- $L_1$ -Norm moderat gehalten werden sollte (Verwendung des approximativen Lösers oder des dynamischen MPF) und die Trotter-Schritte $k_j$ so gewählt werden sollten, dass $t/k_{\\min} \\lesssim 1$ — hier liefert $k_{\\min} = 3$ unter $t = 3$ das Ergebnis $t/k_{\\min} = 1$, wodurch die Komponenten im konvergenten Bereich bleiben, in dem das Modell für den führenden Fehler, auf das sich der statische MPF stützt, gültig ist. Mit diesen Einstellungen erreichen die MPFs mit kleiner Norm hier die Leistung eines konvergierten Einzelkreis-Modells, während dies bei der naiven Basislinie „einfach tiefer gehen“ nicht der Fall ist, wodurch der in Ref. [\\[3\\]](#references) gezeigte Vorteil von Tiefe gegenüber Genauigkeit wiederhergestellt wird. Beachten Sie außerdem, dass einzelne Durchläufe mit Rauschen behaftet sind – bei einer anderen Einreichung desselben Jobs (oder auf einem anderen Backend) kann sich die genaue Reihenfolge verschieben; die robusten Trends zeigen, dass kleine MPFs ( $\\|x\\|_1$ ) gut abschneiden, der große, exakt-statische MPF ( $\\|x\\|_1$ ) durch Hardware-Rauschen verstärkt wird und der übertiefe Einzelkreis durch Rauschen begrenzt ist.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5ac2f8a8",
      "metadata": {},
      "source": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## Nächste Schritte\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Empfehlungen\">\n",
        "  Wenn Sie diese Arbeit interessant fanden, könnten Sie auch an folgenden Materialien interessiert sein:\n",
        "\n",
        "  * [So wählt man die Trotter-Schritte für ein MPF aus](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html) – praktische Hinweise zur Auswahl der Werte für „ $k_j$ “, um Instabilitäten zu vermeiden\n",
        "  * [Anwendung des Näherungsmodells](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html) – Anpassung der $L_1$ -Norm-Beschränkung und der Solver-Optionen für das statische MPF-Näherungsmodell\n",
        "  * [`qiskit-addon-mpf` API-Referenz](https://qiskit.github.io/qiskit-addon-mpf/) – vollständige Dokumentation zu statischen, dynamischen und Backend-Modulen\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "70be41e1",
      "metadata": {},
      "source": [
        "<span id=\"references\" />\n",
        "\n",
        "## Referenzen\n",
        "\n",
        "\\[1] Vázquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. Gut konditionierte Mehrproduktformeln für eine hardwarefreundliche Hamilton-Simulation. [Quantum, 7, 1067 (2023)](https://quantum-journal.org/papers/q-2023-07-25-1067/)\n",
        "\n",
        "\\[2] Zhuk, S., Robertson, N. F., & Bravyi, S. Trotter-Fehlergrenzen und dynamische Mehrproduktformeln für die Hamilton-Simulation. [Physical Review Research, 6(3), 033309 (2024)](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309)\n",
        "\n",
        "\\[3] Robertson, N. F., et al. Durch Tensor-Netzwerke erweiterte dynamische Mehrproduktformeln. [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
}