{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "454a9dfd",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Formules multi-produits pour réduire l'erreur Trotter\"\n",
        "description: \"Utiliser des formules multiproduits dans l'estimation observable afin de réduire l'erreur de Trotter, ou mettre en œuvre une évolution temporelle avec une erreur de Trotter fixe à une profondeur moindre.\"\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",
        "# Formules multi-produits pour réduire l'erreur Trotter\n",
        "\n",
        "*Estimation de la durée d'exécution : quatre minutes sur un processeur Heron r2 (REMARQUE : il s'agit uniquement d'une estimation. (La durée d'exécution peut varier.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c4d0b2f2",
      "metadata": {},
      "source": [
        "<span id=\"learning-outcomes\" />\n",
        "\n",
        "## Acquis d'apprentissage\n",
        "\n",
        "À l'issue de ce tutoriel, vous devriez être en mesure de comprendre les éléments suivants :\n",
        "\n",
        "* Comment les formules multiproduits (MPF) réduisent l'erreur de Trotter dans la simulation hamiltonienne en combinant les valeurs d'espérance issues de plusieurs circuits peu profonds\n",
        "* Dans quels cas les MPF présentent-ils un avantage par rapport aux formules de produits standard, et dans quels cas ne constituent-ils pas l'outil approprié?\n",
        "* Comment calculer les coefficients MPF statiques et dynamiques à l'aide du `qiskit_addon_mpf` package\n",
        "* Comment exécuter un workflow MPF de bout en bout sur du matériel d' IBM Quantum®, y compris la transpilation, la correction des erreurs et le post-traitement\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f5dfd316",
      "metadata": {},
      "source": [
        "<span id=\"prerequisites\" />\n",
        "\n",
        "## Prérequis\n",
        "\n",
        "Nous recommandons aux utilisateurs de se familiariser avec les sujets suivants avant de suivre ce tutoriel :\n",
        "\n",
        "* [Méthodes de compilation pour les circuits de simulation hamiltoniens](/docs/tutorials/compilation-methods-for-hamiltonian-simulation-circuits) — présentation des circuits de Trotter (formule du produit) dans Qiskit.\n",
        "* Les formules de produits dans Qiskit, en particulier les [`SuzukiTrotter`](/docs/api/qiskit/qiskit.synthesis.SuzukiTrotter) classes de synthèse et [`LieTrotter`](/docs/api/qiskit/qiskit.synthesis.LieTrotter) .\n",
        "* [Qiskit primitives et l'interface Estimator](/docs/guides/primitives).\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "07273b26",
      "metadata": {},
      "source": [
        "<span id=\"background\" />\n",
        "\n",
        "## Arrière-plan\n",
        "\n",
        "<span id=\"what-are-multi-product-formulas\" />\n",
        "\n",
        "### Qu'est-ce qu'une formule multiproduit?\n",
        "\n",
        "Lorsqu'on simule des systèmes quantiques sur un ordinateur quantique, l'une des tâches principales consiste à approximer l'opérateur d'évolution temporelle $e^{-iHt}$ pour un hamiltonien $H$. L'approche standard utilise *les formules de produit* (PF), également appelées décompositions de Trotter-Suzuki. Celles-ci décomposent $H = \\sum_{a=1}^d F_a$ en termes dont les opérateurs unitaires individuels $e^{-iF_a t}$ sont faciles à mettre en œuvre, puis approximent l'évolution complète sous la forme d'un produit ordonné de ces opérateurs unitaires plus simples.\n",
        "\n",
        "La formule du produit du premier ordre (Lie-Trotter) est la suivante :\n",
        "\n",
        "$$\n",
        "S_1(t) := \\prod_{a=1}^d e^{-i F_a t},\n",
        "$$\n",
        "\n",
        "ce qui entraîne une erreur quadratique : $S_1(t) = e^{-iHt} + \\mathcal{O}(t^2)$. Les formules symétriques d’ordre supérieur $S_{2\\chi}(t)$, où $\\chi$ désigne l’ordre de la formule du produit symétrique (voir réf. [\\[1\\]](#references) ), convergent plus rapidement selon $e^{-iHt} + \\mathcal{O}(t^{2\\chi+1})$, mais au prix de circuits plus complexes à chaque étape.\n",
        "\n",
        "Pour réduire l'erreur à un ordre *fixe* $\\chi$, on divise généralement le temps d'évolution total $t$ en $k$ petites étapes de Trotter. Chaque étape fournit une approximation de $e^{-iHt/k}$ à l'aide d'une formule de produit, et les étapes sont enchaînées :\n",
        "\n",
        "$$\n",
        "e^{-iHt} \\approx \\left[S_{2\\chi}(t/k)\\right]^k.\n",
        "$$\n",
        "\n",
        "Pour une formule symétrique d’ $2\\chi$ e, l’erreur résiduelle de Trotter évolue alors selon la loi $\\mathcal{O}\\!\\left(t^{2\\chi+1} / k^{2\\chi}\\right)$. Ainsi, l’augmentation de l’ $k$ réduit rapidement l’erreur de Trotter, mais elle accroît également de manière linéaire la profondeur du circuit, ce qui, sur un matériel sujet au bruit, se traduit par une accumulation plus importante de bruit de porte. C'est précisément cette tension entre **l'erreur de Trotter (qui favorise les valeurs de l' $k$ plus élevées)** et **le bruit matériel (qui favorise les valeurs de l' $k$ plus faibles)** que les formules multiproduits sont censées résoudre. Notez que les MPF consistent à combiner les résultats de *différents choix de $k$* dans un ordre fixe $\\chi$ — ils ne modifient pas l'ordre de la formule du produit sous-jacent.\n",
        "\n",
        "**Les formules multiproduits (MPF)** [\\[1\\]](#references) constituent une *combinaison linéaire pondérée* des valeurs attendues obtenues à partir de plusieurs circuits de Trotter moins profonds, chacun utilisant un nombre différent d’étapes de Trotter $k_1, k_2, \\ldots, k_r$ (un ensemble de nombres d’étapes $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",
        "où $\\langle A \\rangle_{k_j}(t)$ est la valeur attendue d’une observable $A$ à l’instant $t$, estimée à partir d’un circuit de Trotter comportant $k_j$ étapes, et où les coefficients $\\{x_j\\}_{j=1}^r$ sont choisis de manière à ce que les termes dominants de l’erreur de Trotter dans la combinaison s’annulent. Nous reviendrons sur cette expression à [l'étape 4](#small-scale-step-4), où nous l'évaluerons explicitement afin d'intégrer nos résultats de Trotter. Le point essentiel d'un point de vue pratique est que le circuit le plus profond du MPF ne nécessite qu' $k_{\\max}$ s d'étapes, ce qui est bien inférieur à l' $k$ unique qui serait nécessaire pour atteindre directement la même erreur de Trotter effective. Grâce à des circuits moins profonds, l'approche MPF est mieux adaptée aux matériels bruyants.\n",
        "\n",
        "<span id=\"how-are-the-coefficients-determined\" />\n",
        "\n",
        "### Comment les coefficients sont-ils déterminés?\n",
        "\n",
        "Il existe deux familles de coefficients MPF :\n",
        "\n",
        "**Les coefficients statiques** sont indépendants de l'hamiltonien, de l'état initial et du temps d'évolution. On les obtient en résolvant un système linéaire $Ax = b$ qui garantit l'annulation des termes d'erreur de Trotter dominants. Pour un ensemble d'étapes de Trotter $\\{k_j\\}_{j=1}^r$ utilisé avec une formule de produit symétrique d'ordre $2\\chi$, le développement de l'erreur de Trotter en puissances inverses de $k_j$ conduit à des équations de contrainte de la forme :\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",
        "où les exposants entiers $\\{\\eta_n\\}$ correspondent aux ordres des termes d'erreur de Trotter successifs pour la formule de produit choisie. Pour un PF *symétrique* d’ $2\\chi$ e d’ordre n, l’erreur dominante dans $\\left[S_{2\\chi}(t/k)\\right]^k$ est proportionnelle à $1/k^{2\\chi}$, avec des corrections successives de l’ordre de $1/k^{2\\chi+2}, 1/k^{2\\chi+4}, \\ldots$ — les exposants sont donc $\\eta_n = 2\\chi + 2n$. Pour les PF non symétriques, les puissances impaires et paires contribuent toutes deux, et $\\eta_n = 2\\chi + n$. Voir la référence [\\[1\\]](#references) pour la démonstration complète. La première équation du système ci-dessus garantit l'absence de biais (le MPF reproduit la valeur exacte de l'espérance dans la limite d' $k_j \\to \\infty$), et les autres équations d' $r-1$ annulent successivement les premiers termes d'erreur de Trotter de type $r-1$. Lorsque la norme d' $L_1$ - $\\|x\\|_1$ qui en résulte est trop élevée (ce qui amplifie le bruit d'échantillonnage), vous pouvez plutôt résoudre un problème d'optimisation approximative qui plafonne $\\|x\\|_1$ tout en minimisant $\\|Ax - b\\|$.\n",
        "\n",
        "**Les coefficients dynamiques** [\\[2\\]](#references), [\\[3\\]](#references) dépendent en outre de l'hamiltonien, de l'état initial et de la durée d'évolution $t$. Ils minimisent la distance, mesurée par la norme de Frobenius, entre l'état réel issu de l'évolution temporelle et l'approximation 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",
        "où $M_{ij}(t) = \\mathrm{Tr}[\\rho_{k_i}(t)\\,\\rho_{k_j}(t)]$ est la matrice de Gram des chevauchements entre les états issus de l'évolution de Trotter pour différents nombres d'étapes $k_i, k_j$, et $L_i(t) = \\mathrm{Tr}[\\rho(t)\\,\\rho_{k_i}(t)]$ mesure le chevauchement avec l'état exact (approximatif). Dans ce tutoriel, ces grandeurs sont calculées efficacement à l'aide de méthodes de réseaux de tenseurs, et plus précisément des backends « TeNPy-based » dans `qiskit_addon_mpf`.\n",
        "\n",
        "<span id=\"when-to-use-mpfs\" />\n",
        "\n",
        "### Quand utiliser les MPF?\n",
        "\n",
        "Les fonds de prévoyance (MPF) sont particulièrement avantageux lorsque :\n",
        "\n",
        "* **La profondeur du circuit constitue le goulot d'étranglement.** Si le bruit généré par le matériel limite la profondeur à laquelle vous pouvez travailler, utilisez des MPF pour obtenir une meilleure précision effective de Trotter à partir de circuits moins profonds.\n",
        "* **Vous avez besoin de valeurs d'espérance précises, et non d'une préparation complète de l'état.** Les MPF opèrent au niveau des valeurs attendues : elles combinent des nombres classiques, et non des états quantiques. Ils sont donc parfaits pour l'estimation à partir d'observables lorsqu'on utilise la primitive « Estimator ».\n",
        "* **Vous enchaînez un nombre modeste de pas de trotteur.** En général, la combinaison d’ $r = 3$ – $5$, avec différents nombres de pas $k_j$, suffit à annuler plusieurs termes d’erreur de Trotter de premier ordre tout en conservant $\\|x\\|_1$ à un niveau raisonnable.\n",
        "\n",
        "<span id=\"when-mpfs-might-not-help\" />\n",
        "\n",
        "### Quand les fonds de prévoyance (MPF) ne sont pas forcément la solution\n",
        "\n",
        "* **Des temps d'évolution très courts.** Lorsque $t$ est suffisamment petit pour qu’une seule formule de Trotter d’ordre inférieur soit déjà précise, la charge supplémentaire liée à l’exécution de plusieurs circuits n’est pas nécessaire.\n",
        "* **Exercices de préparation aux examens d'État.** Les MPF produisent une *valeur attendue* corrigée, et non un état quantique corrigé. Si vous avez besoin de l'état réel évolué dans le temps (par exemple, pour l'utiliser comme entrée d'une autre sous-routine quantique), les MPF ne s'appliquent pas.\n",
        "* **Les nombres d'étapes de trot qui ne respectent pas le régime de convergence.** Le calcul du coefficient statique consiste à développer chaque « $\\left[S_{2\\chi}(t/k_j)\\right]^{k_j}$ » en une série en $t/k_j$; ce développement ne converge correctement que lorsque $t/k_{\\min} \\lesssim 1$. Si $k_{\\min}$ est choisi trop petit pour l’ $t$ donnée, le circuit le moins profond se trouve bien en dehors du régime perturbatif, les termes d’erreur d’ordre supérieur que le MPF ne parvient pas à annuler deviennent importants, et l’annulation peut nécessiter des coefficients élevés. La norme d' $L_1$ $\\|x\\|_1$ constitue un critère de diagnostic pratique : lorsque $\\|x\\|_1 \\gg 1$, la surcharge liée à l'échantillonnage $\\propto \\|x\\|_1^2$ peut l'emporter sur la réduction de l'erreur de Trotter. Pour plus de détails, consultez le [guide sur le choix des marches Trotter](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html).\n",
        "\n",
        "<span id=\"what-this-tutorial-covers\" />\n",
        "\n",
        "### Contenu de ce tutoriel\n",
        "\n",
        "Ce tutoriel présente, en deux étapes, un workflow MPF complet. Tout d'abord, un **exemple de simulation à petite échelle** (chaîne de Heisenberg à 10 qubits) montre comment formuler le problème, calculer les coefficients MPF statiques et dynamiques, et comparer les valeurs attendues obtenues à celles issues de la diagonalisation exacte. Ensuite, **un exemple de calcul sur matériel à grande échelle** (chaîne XXZ de 50 qubits) montre comment transcompiler, exécuter sur un matériel de type « IBM Quantum » avec atténuation des erreurs, puis traiter les résultats à l'aide des coefficients MPF. Tout au long de ce document, nous utilisons ce `qiskit_addon_mpf` package en complément des outils standard de Qiskit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b5d478ce",
      "metadata": {},
      "source": [
        "<span id=\"requirements\" />\n",
        "\n",
        "## Exigences\n",
        "\n",
        "Avant de commencer ce tutoriel, assurez-vous que les éléments suivants sont installés :\n",
        "\n",
        "* Qiskit SDK v2.0 ou version ultérieure, avec prise en charge [de la visualisation](/docs/api/qiskit/visualization)\n",
        "* Qiskit Runtime v0.22 ou plus tard (`pip install qiskit-ibm-runtime`)\n",
        "* Simulateur Qiskit Aer (`pip install qiskit-aer`)\n",
        "* Module complémentaire MPF Qiskit avec le backend « TeNPy » (`pip install \"qiskit-addon-mpf[tenpy]\"`)\n",
        "* Utilitaires complémentaires de 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",
        "## Configuration\n",
        "\n",
        "Nous avons regroupé ci-dessous, dans une seule cellule\\*, toutes\\* les importations de paquets utilisées tout au long de ce tutoriel. `XXPlusYYGate` Nous définissons également un `CollectAndCollapse` passage de transpileur qui fusionne les rotations adjacentes `rxx` et `ryy` en une seule. Cette étape est appliquée à la fois lors de la construction du circuit à l’étape 1 (pour limiter le nombre de portes) et indirectement lorsque nous extrayons la structure en couches pour le MPF dynamique à l’étape 4 (la fonction « TeNPy » attend des portes à deux qubits, et non des paires de rotations non fusionnées).\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",
        "## Exemple de simulateur à petite échelle\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "378e82ba",
      "metadata": {},
      "source": [
        "<span id=\"step-1-map-classical-inputs-to-a-quantum-problem\" />\n",
        "\n",
        "### Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique\n",
        "\n",
        "Nous commençons par un modèle de Heisenberg à 10 qubits sur une ligne, en prenant comme état initial l' $\\vert 0101\\ldots01 \\rangle$ de l'état de Néel. L'hamiltonien est le suivant :\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",
        "où $J$ correspond à l'intensité du couplage entre voisins les plus proches. Nous mesurons le corrélateur ZZ $Z_{L/2-1} Z_{L/2}$ sur une paire de qubits située au milieu de la chaîne, et utilisons les pas de Trotter $k_j = [1, 2, 4]$ avec une formule de produit du second ordre.\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",
        "#### Construire des circuits Trotter\n",
        "\n",
        "Nous créons les circuits mettant en œuvre les évolutions temporelles approximatives de Trotter pour chaque instant et chaque nombre d'étapes de Trotter. Le `CollectAndCollapse` passage défini dans la section « Configuration » regroupe les rotations XX et YY en portes uniques XX+YY, afin de permettre ultérieurement une simulation plus efficace du réseau de tenseurs.\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",
        "### Étape 2 : Optimiser le problème pour l'exécution sur du matériel quantique\n",
        "\n",
        "Pour cet exemple à petite échelle, nous nous concentrons sur le simulateur Aer. Deux transformations ont lieu avant que les circuits ne soient prêts à s'exécuter :\n",
        "\n",
        "1. **Collecte de données au niveau de la simulation hamiltonienne.** `XXPlusYYGate`Dans la cellule « Setup », nous avons créé un `CollectAndCollapse` enchaînement qui fusionne les rotations adjacentes `rxx` et `ryy` en une seule. Nous avons déjà appliqué cette étape lors de la création des circuits Trotter à l'étape 1 (l'appel `pm.run(...)` ). Cela permet à la fois de réduire le nombre de portes à deux qubits et d'obtenir une structure qui se prête mieux à la simulation par réseau de tenseurs pour le calcul ultérieur des coefficients dynamiques.\n",
        "\n",
        "2. **Adaptation à l'ISA du simulateur.** Ci-dessous, nous exécutons le gestionnaire de passes prédéfini de Qiskit afin `optimization_level=3` de transposer chaque circuit de Trotter vers l'architecture du jeu d'instructions (ISA) du simulateur.\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",
        "### Étape 3 : Exécutez à l'aide d' Qiskit primitives\n",
        "\n",
        "Pour l'exemple à petite échelle, nous appliquons les circuits de Trotter convertis en ISA à la `EstimatorV2` primitive implémentée par Aer. Cela nous permet d'obtenir une valeur de référence *sans bruit* pour chaque paire « $(k_j, t)$ » — ce sont ces valeurs « $\\langle A \\rangle_{k_j}(t)$ » que le MPF combinera à l'étape 4. Nous passons rapidement en revue les périodes d'évolution afin de pouvoir tracer ultérieurement la courbe complète de la série chronologique de chaque formule de produit et de la 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",
        "### Étape 4 : Post-traitement et restitution du résultat dans le format classique souhaité\n",
        "\n",
        "C'est à l'étape 4 que le MPF est effectivement élaboré. Même si les coefficients $x_j$ sont *calculés* ici (et que, pour la variante dynamique, ce calcul peut s'avérer très gourmand en ressources), ils constituent, d'un point de vue conceptuel, une méthode classique permettant de combiner les mesures quantiques de l'étape 3 en une seule valeur d'espérance corrigée — nous considérons donc l'ensemble du processus de calcul des coefficients et de combinaison comme un post-traitement.\n",
        "\n",
        "Pour évaluer dans quelle mesure le MPF reflète fidèlement la dynamique réelle, nous calculons tout d'abord les valeurs attendues exactes en fonction du temps en élevant directement l'hamiltonien à la puissance. Cela n'est possible que parce que $L = 10$; dans l'exemple de matériel à grande échelle ci-dessous, nous devrons plutôt nous appuyer sur des estimations issues de réseaux de tenseurs.\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",
        "#### Coefficients MPF statiques\n",
        "\n",
        "Les MPF statiques utilisent des coefficients $x_j$ qui sont indépendants du temps d'évolution, de l'hamiltonien et de l'état initial. Nous établissons le système linéaire $Ax = b$ décrit dans la section « Contexte » et déterminons les coefficients. La matrice $A$ est déterminée par le nombre d'étapes de Trotter $k_j$, l'ordre $\\chi$ de la formule du produit, et le fait que la formule soit symétrique ou non (ce qui détermine les exposants $\\eta_n$ ).\n",
        "\n",
        "Pour notre exemple à petite échelle, nous utilisons $k_j = [1, 2, 4]$ avec une formule de Suzuki-Trotter non symétrique d’ordre $2\\chi=2$ (d’où $\\chi=1$ et $\\eta_n = 2 + n$, ce qui donne $\\eta_0 = 2,\\, \\eta_1 = 3$ ). Le système se présente alors comme suit :\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",
        "La première ligne garantit l'absence de biais ($\\sum_j x_j = 1$); les deuxième et troisième lignes annulent respectivement les termes d'erreur de Trotter de premier ordre $1/k^2$ et de second ordre $1/k^3$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4f2ca1e2",
      "metadata": {},
      "source": [
        "<span id=\"set-up-the-lse\" />\n",
        "\n",
        "##### Configurer le LSE\n",
        "\n",
        "Nous utilisons `setup_static_lse` de `qiskit_addon_mpf.static` pour construire la matrice $A$ et le vecteur du côté droit $b$ décrits ci-dessus. La matrice $A$ dépend non seulement de $k_j$, mais aussi du choix de la formule du produit — en particulier de son *ordre* $\\chi$ et du fait qu’elle soit *symétrique* ou non. Le `symmetric` drapeau contrôle le schéma de l'exposant $\\eta_n$ (les formules symétriques ne produisent que des termes d'erreur de Trotter de puissance paire; voir la réf. [\\[1\\]](#references) ). Il convient de noter que, comme le montre la référence [\\[2\\]](#references), il n’est pas strictement nécessaire de définir `symmetric=True` même lorsque la fonction de performance sous-jacente est symétrique — l’estimation LSE non symétrique reste valable (elle impose toutefois des contraintes supplémentaires inutiles).\n",
        "\n",
        "Dans notre exemple, nous avons déjà défini `order = 2` et `symmetric = False` à l'étape 1.\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": [
        "Vérifiez que la matrice $A$ et le vecteur $b$ correspondent bien au système présenté ci-dessus.\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": [
        "Une fois l'équation de LSE établie, on détermine les coefficients statiques $x_j$ à l'aide de `lse.solve()` (il s'agit de la solution directe de l' $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",
        "##### Optimiser l' $x$ s à l'aide d'un modèle exact\n",
        "\n",
        "Au lieu de calculer $x = A^{-1}b$, vous pouvez utiliser [setup\\_exact\\_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_exact_model.html) pour construire une instance de [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) qui utilise le LSE comme contraintes et dont la solution optimale donnera $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",
        "##### Optimiser l' $x$ s à l'aide d'un modèle approximatif\n",
        "\n",
        "Il peut arriver que la norme d’ $L_1$, pour l’ensemble choisi de valeurs d’ $k_j$, soit jugée trop élevée. Si tel est le cas et que vous ne pouvez pas choisir un autre ensemble de valeurs pour « $k_j$ », vous pouvez utiliser une solution approximative qui limite la norme d’ $L_1$ à un seuil donné tout en minimisant « $\\|Ax - b\\|$ ». Consultez le guide intitulé «[ Comment utiliser le modèle approximatif](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",
        "#### Coefficients MPF dynamiques\n",
        "\n",
        "La méthode MPF statique annule les termes d'erreur de Trotter d'une manière indépendante de l'hamiltonien et de l'état; elle ne produit donc pas nécessairement l'erreur d'approximation la plus faible possible pour un hamiltonien et un état initial donnés. La méthode MPF dynamique (réf. [\\[2\\]](#references), [\\[3\\]](#references) ) détermine quant à elle des coefficients dépendants du temps $x_i(t)$ qui minimisent la distance de la norme de Frobenius $\\|\\rho(t) - \\mu^D(t)\\|_F^2$ à chaque instant $t$. Comme indiqué dans la section « Contexte », cela nécessite la matrice de chevauchement $M_{ij}(t)$ entre les états évolués selon la méthode de Trotter et le chevauchement $L_i(t)$ avec l’état exact — que nous estimons tous deux à l’aide de backends de réseaux de tenseurs ( TeNPy ) dans `qiskit_addon_mpf`.\n",
        "\n",
        "Pour configurer le LSE dynamique, il nous faut trois éléments :\n",
        "\n",
        "1. **Une fabrique d'évoluteurs approximatifs** que l'extension exécutera pour chaque $k_j$ afin de générer $\\rho_{k_j}(t)$ sous forme de MPS/MPO. Nous le construisons à partir de la structure en couches du circuit de Trotter d’ordre $2$ (une couche par `slice_by_depth`), enveloppé sous la forme d’un `LayerwiseEvolver` avec des paramètres de troncature de type TeNPy.\n",
        "2. **Une fonction d'évolution exacte** qui produit une référence de haute précision $\\rho(t)$. Nous utilisons un circuit de Suzuki-Trotter du quatrième ordre à petit pas de temps (`dt=0.1`, `order=4`) comme approximation de l'évolution exacte.\n",
        "3. Une **« identity factory »** et un **MPS d'état initial** qui servent de base à la simulation « TeNPy ».\n",
        "\n",
        "La cellule ci-dessous permet de créer l'usine d'évoluteurs approximatifs.\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",
        "  Les options de `LayerwiseEvolver` qui déterminent les détails de la simulation du réseau tensoriel doivent être choisies avec soin pour éviter de créer un problème d'optimisation mal défini.\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1c702f67",
      "metadata": {},
      "source": [
        "`dt=0.1` Nous approximons l'état exact évolué dans le temps à l'aide d'une formule de Suzuki-Trotter du quatrième ordre, en utilisant un petit pas de temps. Les paramètres de troncature de l' TeNPy e peuvent avoir une incidence sur la précision; il est donc important d'explorer toute une gamme de valeurs.\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": [
        "Enfin, nous définissons un `identity_factory` qui donne l'état MPO initial et préparons l'état initial de Néel sous la forme d'un MPS correspondant au réseau utilisé par le modèle de Trotter en couches.\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": [
        "Une fois les modèles mis en place, nous calculons désormais les coefficients dynamiques à chaque instant d'évolution. Pour chaque $t$, `setup_dynamic_lse` le programme construit les matrices de chevauchement correspondantes à l'aide de TeNPy, et `setup_frobenius_problem` renvoie une solution `cvxpy.Problem` qui minimise le coût de la norme de Frobenius. Le solveur renvoie des coefficients $x_j(t)$ adaptés à cette période; nous les rassemblons dans `mpf_dynamic_coeffs_list`. Si le solveur échoue pour une « $t$ » donnée, nous revenons à des coefficients nuls afin que la boucle se poursuive.\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",
        "#### Combiner les valeurs attendues de Trotter avec les coefficients MPF\n",
        "\n",
        "Nous évaluons à présent la valeur de « $\\langle A \\rangle_{\\text{MPF}}(t) = \\sum_j x_j \\, \\langle A \\rangle_{k_j}(t)$ » pour chaque ensemble de coefficients (static-exact, static-approximate et dynamic), nous propageons les erreurs-types par circuit, puis nous représentons graphiquement les séries chronologiques obtenues par rapport à la courbe de diagonalisation exacte.\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": [
        "Le graphique ci-dessus illustre l'interaction entre l'erreur de Trotter et l'erreur d'échantillonnage.\n",
        "\n",
        "* **Erreur de Trotter.** Les courbes de chaque produit (marqueurs gris) s'écartent de plus en plus de la courbe exacte à mesure que le temps passe. Le circuit « $k=1$ » présente l'écart le plus important et est le moins profond, mais il se situe déjà dans le régime où « $t/k \\gtrsim 1$ », de sorte que le terme d'erreur principal « $1/k^{2}$ » est important. Les combinaisons MPF (marqueurs colorés) annulent plusieurs de ces termes d'erreur de Trotter dominants, ce qui leur permet de suivre la courbe exacte de bien plus près que n'importe quel circuit d' $k_j$ ation pris isolément. L'écart résiduel reflète les termes de Trotter d'ordre supérieur que la MPF *n* 'annule pas : d'ordre $2$, $r=3$ une MPF statique n'élimine que les deux premiers ordres d'erreur, et à de grandes $t/k_{\\min}$, la partie résiduelle non annulée finit par dominer — la MPF ne garantit donc pas que les circuits très peu profonds restent précis à des instants arbitraires.\n",
        "\n",
        "* **Erreur d'échantillonnage.** Les barres d’erreur plus larges sur les courbes MPF sont une conséquence directe de la combinaison linéaire : la propagation des erreurs-types indépendantes par circuit $\\sigma_{k_j}$ donne une variance totale $\\sigma_{\\text{MPF}}^2 = \\sum_j x_j^2 \\, \\sigma_{k_j}^2$. Par conséquent, plus l’ $\\|x\\|_2$ (et, en pratique, l’ $\\|x\\|_1$, que nous contrôlons) est grande, plus il faut de mesures pour atteindre une incertitude cible donnée. C'est là le compromis qui sous-tend l'option « approximate-solver » dans la section « Contexte » : nous limitons la valeur de $\\|x\\|_1$ afin que cette surcharge reste gérable. Il est essentiel de noter que, contrairement à l'erreur de Trotter, l'erreur d'échantillonnage diminue avec l' $1/\\sqrt{N_{\\text{shots}}}$; elle peut donc toujours être réduite en effectuant davantage de tirs.\n",
        "\n",
        "Dans l'exemple de matériel à grande échelle ci-dessous, le bruit matériel constitue une source d'erreur supplémentaire à chaque é $\\langle A \\rangle_{k_j}$ e, qui est elle aussi amplifiée par les coefficients du MPF. Nous verrons dans cette section comment l'atténuation des erreurs interagit avec les MPF.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6fa763ff",
      "metadata": {},
      "source": [
        "<span id=\"large-scale-hardware-example\" />\n",
        "\n",
        "## Exemple de matériel à grande échelle\n",
        "\n",
        "Dans cette section, nous élargissons le problème au-delà de ce qu'il est possible de simuler avec exactitude. Nous reproduisons certains des résultats présentés dans la référence [\\[3\\]](#references), en utilisant une chaîne XXZ de 50 qubits aux instan $t = 3$. Nous suivons la même procédure en quatre étapes que dans l'exemple à petite échelle, en ciblant cette fois-ci du matériel quantique réel doté d'un système d'atténuation des erreurs. Comme dans le modèle, chaque étape est indiquée directement dans le code, et une même étape peut s'étendre sur plusieurs cellules lorsque ses résultats intermédiaires méritent d'être examinés.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "481fea70",
      "metadata": {},
      "source": [
        "La procédure suit le même principe que l'exemple à petite échelle : définir un hamiltonien, choisir les paramètres de Trotter, calculer les coefficients MPF (statiques et dynamiques), puis construire les circuits. Les principales différences sont les suivantes :\n",
        "\n",
        "* Un **hamiltonien XXZ** à 50 sites avec des couplages aléatoires tirés de $\\mathcal{U}(0.5, 1.5)$ (réf. [\\[3\\]](#references) ).\n",
        "* Une formule de Trotter **symétrique** du second ordre avec l' $k_j = [3, 4, 6]$ e (donc $\\chi=1$, `symmetric=True`).\n",
        "* Un temps d'évolution fixe unique $t = 3$. Avec $k_{\\min}=3$, on obtient $t/k_{\\min}=1$, ce qui maintient les constituants de faible profondeur dans le régime de convergence de Trotter, où le modèle d'erreur dominante sur lequel repose le MPF est valide.\n",
        "* Une **comparaison supplémentaire sur un seul circuit avec des étapes Trotter $k = 10$**, utilisée comme référence. Nous avons choisi « $k = 10$ » car sa profondeur de deux qubits sur le matériel est supérieure à celle du constituant MPF le plus profond ( $k_{\\max}=6$ ), à laquelle s’ajoute la surcharge liée à l’exécution de plusieurs circuits MPF — une profondeur suffisante pour que le système soit limité par le bruit, régime dans lequel la combinaison MPF devrait surpasser la référence à circuit unique. Il s'agit d'une comparaison de « circuit profond unique » par rapport à la combinaison MPF, et non d'un circuit visant à corriger l'erreur de Trotter effective du MPF (ce qui nécessiterait bien plus d'étapes).\n",
        "\n",
        "Notez que, même si nous en sommes encore à l'étape 1 (mappage et construction du circuit), nous précalculons également les coefficients dynamiques en même temps que les coefficients statiques dans cette cellule. Les coefficients dynamiques dépendent de $H$ et $t$, mais pas des mesures quantiques; ils peuvent donc être calculés à tout moment avant l'étape 4. Nous procédons ainsi afin de regrouper tous les paramètres spécifiques au MPF en un seul endroit.\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": [
        "Nous optimisons à présent les circuits pour le backend choisi. `optimization_level=3`Nous utilisons le gestionnaire de passes prédéfini de Qiskit, qui sélectionne automatiquement un ensemble optimal de qubits physiques et achemine chaque circuit vers la topologie du dispositif.\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": [
        "L'exécution de circuits plus complexes sur du matériel réel nécessite des mesures énergiques d'atténuation des erreurs. Nous permettons le découplage dynamique, la rotation des portes et des mesures, l'atténuation des erreurs de mesure et l'extrapolation sans bruit (ZNE). Il convient de noter que les facteurs de bruit ZNE que nous utilisons ici (`1, 1.2, 1.4`) sont inférieurs à ceux d'un scénario de circuit peu profond, car les constituants MPF situés plus en profondeur sont déjà proches du seuil de bruit et des amplifications importantes du bruit les feraient dépasser le seuil à partir duquel l'extrapolation ZNE n'est plus fiable.\n",
        "\n",
        "Nous soumettons les quatre circuits (les trois composants MPF disponibles à l'adresse $k_j = [3, 4, 6]$ ainsi que la configuration de référence de l' $k = 10$ ) dans une seule tâche 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": [
        "Nous extrayons les valeurs attendues et les écarts-types par circuit à partir des résultats du calcul, puis nous les combinons avec chaque ensemble de coefficients MPF exactement comme dans l'exemple à petite échelle : $\\langle A \\rangle_{\\text{MPF}} = \\sum_j x_j \\, \\langle A \\rangle_{k_j}$, avec la variance propagée $\\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": [
        "Quelques remarques concernant les résultats matériels ci-dessus :\n",
        "\n",
        "* **Aller plus loin n'est pas sans conséquence sur le matériel.** Les courbes de référence à circuit unique parlent d'elles-mêmes : le circuit $k = 6$ est pratiquement exact ( $-0.256$ par rapport à la référence $-0.244$ ), tandis que la courbe de référence $k = 10$, plus approfondie, est *moins bonne* ( $-0.061$, avec un écart de $\\sim 0.18$ ), et non meilleure. Une fois que l'erreur de Trotter est déjà faible, l'ajout d'étapes ne fait pour l'essentiel qu'alourdir le circuit et accroître le bruit de porte ainsi que la décohérence. C'est précisément le type de configuration pour lequel les MPF ont été conçus : atteindre la précision d'un circuit profond en utilisant uniquement des composants de faible profondeur.\n",
        "\n",
        "* **Un MPF à petite norme surpasse le circuit unique profond.** Le MPF « approximativement statique » (plafonné à $\\|x\\|_1 \\approx 2$ ) s'établit à $-0.259$, soit à $\\sim 0.015$ de la valeur de référence, et bien plus proche que la valeur de référence de $k = 10$. Le MPF dynamique ( $-0.127$ ) dépasse lui aussi largement ce seuil de référence. Les deux ne combinent que les circuits d’ $k_j = [3, 4, 6]$ s peu profonds, mais parviennent néanmoins à trouver une réponse que le circuit unique profond n’aurait pas pu trouver.\n",
        "\n",
        "* **La norme des coefficients est plus importante que l'optimalité mathématique.** Le MPF statique exact présente une norme de coefficient de $\\|x\\|_1 = 4.66$ et constitue le *pire* estimateur de tous ( $-0.567$, avec un écart supérieur à $0.3$ ) : cette norme élevée amplifie le bruit résiduel de la porte, la décohérence et l'erreur ZNE sur chaque $\\langle A \\rangle_{k_j}$ d'un facteur à peu près identique, annulant ainsi l'effet de l'annulation de l'erreur de Trotter qu'il permet d'obtenir. Le plafonnement de la norme (le solveur quasi-statique, $\\|x\\|_1 \\approx 2$ ) élimine ce problème de surcharge et fournit la meilleure estimation — même si ses coefficients n'annulent plus exactement l'erreur de Trotter dominante.\n",
        "\n",
        "* **Les circuits individuels peu profonds peuvent tout de même être compétitifs.** Le seul composant $k = 6$ ($-0.256$) est lui-même pratiquement exact ici — lors de cette exécution, il est même légèrement plus précis que le MPF « approximate-static ». Le problème, c'est qu'on ne sait pas à l'avance *quelle* $k$ se situe exactement dans cette zone idéale où « la convergence est atteinte mais sans que le bruit ne soit encore limitant », et que le choix qui semble le plus sûr, à savoir simplement aller plus loin ( $k = 10$ ) pour garantir la convergence de Trotter, est précisément celui qui échoue. Le MPF propose une combinaison raisonnée de circuits peu profonds qui ne nécessite pas de deviner la profondeur adéquate.\n",
        "\n",
        "Concrètement, cela signifie qu’au niveau matériel, les MPF doivent être associés à une atténuation d’erreur efficace sur chaque é $\\langle A \\rangle_{k_j}$ individuelle, que la norme du coefficient $L_1$ doit rester modérée (utiliser le solveur approximatif ou le MPF dynamique), et que le nombre d’étapes de Trotter $k_j$ doit être choisi de manière à ce que $t/k_{\\min} \\lesssim 1$ — ici $k_{\\min} = 3$ sur $t = 3$ donne $t/k_{\\min} = 1$, en maintenant les composantes dans le régime de convergence où le modèle d’erreur dominante sur lequel repose le MPF statique est valide. Avec ces choix, les MPF à petite norme présentés ici égalent les performances d’un circuit unique convergent, contrairement à la méthode de référence naïve consistant simplement à « aller plus en profondeur », ce qui confirme l’avantage « profondeur contre précision » mis en évidence dans la référence [\\[3\\]](#references). Il convient également de noter que les exécutions individuelles sont sujettes au bruit : lors d’une soumission différente du même travail (ou sur un backend différent), l’ordre exact peut varier; les tendances stables montrent que les MPF de type « small- $\\|x\\|_1$ » donnent de bons résultats, que le MPF de type « large- $\\|x\\|_1$ » (exact-static) est amplifié par le bruit matériel, et que le circuit unique « over-deep » est limité par le bruit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5ac2f8a8",
      "metadata": {},
      "source": [
        "<span id=\"next-steps\" />\n",
        "\n",
        "## Etapes suivantes\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Recommandations\">\n",
        "  Si ce travail vous a paru intéressant, les documents suivants pourraient vous intéresser :\n",
        "\n",
        "  * [Comment choisir les pas de Trotter pour un MPF](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html) — conseils pratiques pour sélectionner les valeurs d’ $k_j$ s afin d’éviter les instabilités\n",
        "  * [Comment utiliser le modèle approximatif](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html) — réglage de la contrainte de norme « $L_1$ » et des options du solveur pour le MPF statique approximatif\n",
        "  * [`qiskit-addon-mpf` Référence API](https://qiskit.github.io/qiskit-addon-mpf/) — documentation complète sur les modules statiques, dynamiques et backend\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "70be41e1",
      "metadata": {},
      "source": [
        "<span id=\"references\" />\n",
        "\n",
        "## Références\n",
        "\n",
        "\\[1] Vazquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. Formules multiproduits bien conditionnées pour une simulation hamiltonienne adaptée au matériel. [Quantum, vol. 7, p. 1067 (2023)](https://quantum-journal.org/papers/q-2023-07-25-1067/)\n",
        "\n",
        "\\[2] Zhuk, S., Robertson, N. F., & Bravyi, S. Limites d'erreur de Trotter et formules dynamiques multiproduits pour la simulation hamiltonienne. [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. Formules dynamiques multiproduits améliorées par un réseau de tenseurs. [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
}