Formules multi-produits pour réduire l'erreur de Trotter

Utilisation estimée : quatre minutes sur un processeur Heron r2 (REMARQUE : il s'agit d'une estimation uniquement. Votre temps d'exécution peut varier.)

Résultats d'apprentissage

• Comment les formules multi-produits (MPF) réduisent l'erreur de Trotter dans la simulation hamiltonienne en combinant les valeurs moyennes de plusieurs circuits peu profonds

• Quand les MPF sont avantageuses par rapport aux formules produits standard et quand elles ne sont pas le bon outil

• Comment calculer les coefficients MPF statiques et dynamiques à l'aide du package qiskit_addon_mpf

• Comment exécuter un flux de travail MPF de bout en bout sur du matériel IBM Quantum®, y compris la transpilation, l'atténuation des erreurs et le post-traitement

Prérequis

• Méthodes de compilation pour les circuits de simulation hamiltonienne — présente les circuits de Trotter (formule produit) dans Qiskit.

• Les formules produits dans Qiskit, en particulier les classes de synthèse SuzukiTrotter et LieTrotter.

• Les primitives Qiskit et l'interface Estimator.

Contexte

Que sont les formules multi-produits ?

Lors de la simulation de systèmes quantiques sur un ordinateur quantique, une tâche centrale consiste à approximer l'opérateur d'évolution temporelle $e^{-iHt}$ pour un hamiltonien $H$. L'approche standard utilise des formules produits (PF), également appelées décompositions de Trotter-Suzuki. Celles-ci décomposent $H = \sum_{a=1}^d F_a$ en termes dont les unitaires individuels $e^{-iF_a t}$ sont efficaces à mettre en œuvre, puis approximent l'évolution complète comme un produit ordonné de ces unitaires plus simples.

La formule produit du premier ordre (Lie-Trotter) est :

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

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 produit symétrique (voir Réf. [1]), convergent plus rapidement sous la forme $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$, mais au prix de circuits plus profonds par pas.

Pour réduire l'erreur à un ordre $\chi$ fixe, on divise généralement le temps d'évolution total $t$ en $k$ pas de Trotter plus petits. Chaque pas approxime $e^{-iHt/k}$ par une formule produit et les pas sont concaténés :

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

Pour une formule symétrique d'ordre $2\chi$, l'erreur de Trotter résiduelle évolue alors comme $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$. Ainsi, augmenter $k$ supprime rapidement l'erreur de Trotter — mais cela approfondit également linéairement le circuit, et sur du matériel bruité, cela signifie plus de bruit de porte accumulé. Cette tension entre l'erreur de Trotter (favorise un $k$ plus grand) et le bruit matériel (favorise un $k$ plus petit) est précisément ce que les formules multi-produits sont conçues pour résoudre. Notez que les MPF consistent à combiner les résultats de différents choix de $k$ à un ordre $\chi$ fixe — elles ne changent pas l'ordre de la formule produit sous-jacente.

Les formules multi-produits (MPF) [1] construisent une combinaison linéaire pondérée de valeurs moyennes obtenues à partir de plusieurs circuits de Trotter moins profonds, chacun utilisant un nombre différent de pas de Trotter $k_1, k_2, \ldots, k_r$ (un ensemble de $r$ nombres de pas) :

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

où $\langle A \rangle_{k_j}(t)$ est la valeur moyenne d'une observable $A$ au temps $t$ estimée à partir d'un circuit de Trotter avec $k_j$ pas, et les coefficients $\{x_j\}_{j=1}^r$ sont choisis de sorte que les termes d'erreur de Trotter dominants dans la combinaison s'annulent. Nous reviendrons sur cette expression à l'Étape 4, où nous l'évaluons explicitement pour combiner nos résultats de Trotter. Le point pratique clé est que le circuit le plus profond de la MPF n'a besoin que de $k_{\max}$ pas, ce qui est bien inférieur au $k$ unique qui serait nécessaire pour atteindre directement la même erreur de Trotter effective. Les circuits moins profonds rendent l'approche MPF mieux adaptée au matériel bruité.

Comment les coefficients sont-ils déterminés ?

Il existe deux familles de coefficients MPF :

Les coefficients statiques sont indépendants du hamiltonien, de l'état initial et du temps d'évolution. Ils sont trouvés en résolvant un système linéaire $Ax = b$ qui impose l'annulation des termes d'erreur de Trotter dominants. Pour un ensemble de pas de Trotter $\{k_j\}_{j=1}^r$ utilisé avec une formule produit symétrique d'ordre $2\chi$, développer l'erreur de Trotter en puissances inverses de $k_j$ conduit à des équations de contrainte de la forme :

$$\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),$$

où les exposants entiers $\{\eta_n\}$ sont les ordres des termes d'erreur de Trotter successifs pour la formule produit choisie. Pour une PF symétrique d'ordre $2\chi$, l'erreur dominante dans $\left[S_{2\chi}(t/k)\right]^k$ évolue comme $1/k^{2\chi}$, avec des corrections suivantes à $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ — donc les exposants sont $\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 Réf. [1] pour la dérivation complète. La première équation du système ci-dessus garantit l'absence de biais (la MPF reproduit la valeur d'expectation exacte dans la limite $k_j \to \infty$), et les $r-1$ équations restantes annulent successivement les $r-1$ premiers termes d'erreur de Trotter. Lorsque la norme $L_1$ résultante $\|x\|_1$ est trop grande (ce qui amplifie le bruit d'échantillonnage), tu peux à la place résoudre une optimisation approximative qui plafonne $\|x\|_1$ tout en minimisant $\|Ax - b\|$.

Les coefficients dynamiques [2], [3] dépendent en plus du hamiltonien, de l'état initial et du temps d'évolution $t$. Ils minimisent la distance de norme de Frobenius entre l'état réellement évolué dans le temps et l'approximation MPF :

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

où $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ est la matrice de Gram des chevauchements entre les états évolués par Trotter pour différents nombres de pas $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 quantités sont calculées efficacement à l'aide de méthodes de réseaux de tenseurs, en particulier les backends basés sur TeNPy dans qiskit_addon_mpf.

Quand utiliser les MPF

Les MPF sont les plus avantageuses lorsque :

• La profondeur du circuit est le goulot d'étranglement. Si le bruit matériel limite la profondeur que tu peux exécuter, utilise les MPF pour obtenir une précision de Trotter effective plus élevée à partir de circuits moins profonds.

• Tu as besoin de valeurs moyennes précises, pas d'une préparation d'état complète. Les MPF opèrent au niveau des valeurs moyennes — elles combinent des nombres classiques, pas des états quantiques. Elles sont donc idéales pour l'estimation d'observables lors de l'utilisation de la primitive Estimator.

• Tu combines un nombre modeste de nombres de pas de Trotter. En général, combiner $r = 3$–$5$ nombres de pas différents $k_j$ suffit pour annuler plusieurs termes d'erreur de Trotter dominants tout en gardant $\|x\|_1$ gérable.

Quand les MPF pourraient ne pas aider

• Temps d'évolution très courts. Lorsque $t$ est suffisamment petit pour qu'une seule formule de Trotter d'ordre faible soit déjà précise, le surcoût lié à l'exécution de plusieurs circuits est inutile.

• Tâches de préparation d'état. Les MPF produisent une valeur moyenne corrigée, pas un état quantique corrigé. Si tu as besoin de l'état réellement évolué dans le temps (par exemple, comme entrée d'une autre sous-routine quantique), les MPF ne s'appliquent pas.

• Nombres de pas de Trotter qui violent le régime de convergence. La dérivation des coefficients statiques développe chaque $\left[S_{2\chi}(t/k_j)\right]^{k_j}$ individuel comme une série en $t/k_j$ ; ce développement ne converge bien que lorsque $t/k_{\min} \lesssim 1$. Si $k_{\min}$ est choisi trop petit pour le $t$ donné, le circuit le moins profond se trouve bien en dehors du régime perturbatif, les termes d'erreur d'ordre supérieur que la MPF laisse non annulés deviennent grands, et l'annulation peut nécessiter de grands coefficients. La norme $L_1$ $\|x\|_1$ est le diagnostic pratique : lorsque $\|x\|_1 \gg 1$, le surcoût d'échantillonnage $\propto \|x\|_1^2$ pourrait l'emporter sur la réduction de l'erreur de Trotter. Voir le guide sur le choix des pas de Trotter pour plus de détails.

Ce que couvre ce tutoriel

Ce tutoriel présente un flux de travail MPF complet, de bout en bout, en deux temps. D'abord, un exemple à petite échelle sur simulateur (chaîne de Heisenberg à 10 qubits) montre comment configurer le problème, calculer les coefficients MPF statiques et dynamiques, et comparer les valeurs moyennes résultantes à la diagonalisation exacte. Ensuite, un exemple à grande échelle sur du matériel (chaîne XXZ à 50 qubits) montre comment transpiler, exécuter sur du matériel IBM Quantum avec atténuation des erreurs, et post-traiter les résultats à l'aide des coefficients MPF. Tout au long, nous utilisons le package qiskit_addon_mpf avec les outils Qiskit standard.

Prérequis

Avant de commencer ce tutoriel, assure-toi d'avoir les éléments suivants installés :

• Qiskit SDK v2.0 ou version ultérieure, avec le support de visualisation

• Qiskit Runtime v0.22 ou version ultérieure (pip install qiskit-ibm-runtime)

• Simulateur Qiskit Aer (pip install qiskit-aer)

• Addon Qiskit MPF avec le backend TeNPy (pip install "qiskit-addon-mpf[tenpy]")

• Addons Qiskit utilities (pip install qiskit-addon-utils)

• SciPy (pip install scipy)

Configuration

Ci-dessous, nous rassemblons tous les imports de packages utilisés tout au long de ce tutoriel dans une seule cellule. Nous définissons également une passe de transpileur CollectAndCollapse qui fusionne les rotations rxx et ryy adjacentes en une seule XXPlusYYGate. Cette passe est appliquée à la fois lors de la construction du circuit à l'Étape 1 (pour maintenir un faible nombre de portes) et indirectement lorsque nous extrayons la structure en couches pour la MPF dynamique à l'Étape 4 (TeNPy attend des portes à deux qubits, pas des paires de rotations non fusionnées).

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

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

import warnings

import numpy as np
import matplotlib.pyplot as plt
from functools import partial
from copy import deepcopy

from qiskit import QuantumCircuit
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector
from qiskit.synthesis import SuzukiTrotter
from qiskit.transpiler import CouplingMap, PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)

from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService

from qiskit_addon_utils.problem_generators import (
    generate_xyz_hamiltonian,
    generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.static import setup_static_lse
from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import (
    setup_exact_problem,
    setup_sum_of_squares_problem,
    setup_frobenius_problem,
)
from qiskit_addon_mpf.backends.tenpy_layers import (
    LayerModel,
    LayerwiseEvolver,
)
from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state

from scipy.linalg import expm

# Suppress TeNPy's `unit_cell_width` future-API warning. The default
# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what
# `CouplingMap.from_line(...)` produces here, so the warning is informational.
warnings.filterwarnings(
    "ignore",
    message=r".*unit_cell_width.*",
    category=UserWarning,
)

# --- Helper: collect XX + YY rotations into a single gate ---
def filter_function(node):
    return node.op.name in {"rxx", "ryy"}

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

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

collapse_function = partial(
    collapse_to_operation,
    collapse_function=collapse_to_xx_plus_yy,
)

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

Exemple de simulateur à petite échelle

Étape 1 : Convertir les entrées classiques en un problème quantique

Nous commençons avec un modèle de Heisenberg à 10 qubits sur une ligne, en utilisant l'état de Néel $\vert 0101\ldots01 \rangle$ comme état initial. L'hamiltonien est :

$$\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),$$

où $J$ est la force de couplage entre plus proches voisins. Nous mesurons le corrélateur ZZ $Z_{L/2-1} Z_{L/2}$ sur une paire de qubits au milieu de la chaîne, et utilisons des pas de Trotter $k_j = [1, 2, 4]$ avec une formule produit du second ordre.

L = 10

# Generate coupling map and Hamiltonian
coupling_map = CouplingMap.from_line(L, bidirectional=False)

hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(1.0, 1.0, 1.0),
    ext_magnetic_field=(0.0, 0.0, 0.0),
)
print(hamiltonian)

SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

# Observable: ZZ on the middle pair of qubits
observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)

SparsePauliOp(['IIIIZZIIII'],
              coeffs=[1.+0.j])

# MPF parameters
mpf_trotter_steps = [1, 2, 4]
order = 2
symmetric = False

trotter_times = np.arange(0.5, 1.55, 0.1)
exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)

Construire les circuits de Trotter

Nous créons les circuits implémentant les évolutions temporelles de Trotter approximatives pour chaque point temporel et chaque nombre de pas de Trotter. La passe CollectAndCollapse définie dans la section Configuration regroupe les rotations XX et YY en portes XX+YY uniques, afin de préparer une simulation par réseau de tenseurs plus efficace par la suite.

# Initial Neel state preparation
initial_state_circ = QuantumCircuit(L)
initial_state_circ.x([i for i in range(L) if i % 2 != 0])

all_circs = []
for total_time in trotter_times:
    mpf_trotter_circs = [
        generate_time_evolution_circuit(
            hamiltonian,
            time=total_time,
            synthesis=SuzukiTrotter(reps=num_steps, order=order),
        )
        for num_steps in mpf_trotter_steps
    ]

    mpf_trotter_circs = pm.run(
        mpf_trotter_circs
    )  # Collect XX and YY into XX + YY

    mpf_circuits = [
        initial_state_circ.compose(circuit) for circuit in mpf_trotter_circs
    ]
    all_circs.append(mpf_circuits)

mpf_circuits[-1].draw("mpl", fold=-1)

Étape 2 : Optimiser le problème pour l'exécution sur matériel quantique

Pour l'exemple à petite échelle, nous ciblons le simulateur Aer. Deux transformations ont lieu avant que les circuits soient prêts à être exécutés :

1. Regroupement de portes au niveau de la simulation hamiltonienne. Dans la cellule de Configuration, nous avons construit une passe CollectAndCollapse qui fusionne les rotations rxx et ryy adjacentes en une seule XXPlusYYGate. Nous avons déjà appliqué cette passe lors de la construction des circuits de Trotter à l'étape 1 (l'appel pm.run(...)). Cela réduit à la fois le nombre de portes à deux qubits et produit une structure plus propice à la simulation par réseau de tenseurs pour le calcul des coefficients dynamiques par la suite.

2. Abaissement vers l'ISA du simulateur. Ci-dessous, nous exécutons le gestionnaire de passes prédéfini de Qiskit avec optimization_level=3 pour abaisser chaque circuit de Trotter vers l'architecture de jeu d'instructions (ISA) du simulateur.

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

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

Étape 3 : Exécuter en utilisant les primitives Qiskit

Pour l'exemple à petite échelle, nous exécutons les circuits de Trotter abaissés vers l'ISA via la primitive EstimatorV2 reposant sur Aer. Cela nous donne une valeur de référence sans bruit pour chaque paire $(k_j, t)$ — ce sont les valeurs $\langle A \rangle_{k_j}(t)$ que le MPF combinera à l'étape 4. Nous balayons les temps d'évolution afin de pouvoir ensuite tracer la courbe complète de la série temporelle de chaque formule produit individuelle et du MPF.

estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for isa_circuits in isa_circs_all_times:
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals_all_times.append([res.data.evs for res in result])
    mpf_stds_all_times.append([res.data.stds for res in result])

Étape 4 : Post-traiter et renvoyer le résultat dans le format classique souhaité

L'étape 4 est celle où le MPF est effectivement construit. Même si les coefficients $x_j$ sont calculés ici (et, pour la variante dynamique, ce calcul peut être intensif), sur le plan conceptuel ils constituent une recette classique pour combiner les mesures quantiques de l'étape 3 en une seule valeur d'espérance corrigée — nous traitons donc l'ensemble du flux de calcul et de combinaison des coefficients comme du post-traitement.

Pour évaluer dans quelle mesure le MPF suit la dynamique réelle, nous calculons d'abord les valeurs d'espérance exactes évoluées dans le temps en exponentiant directement l'hamiltonien. Cela n'est faisable que parce que $L = 10$ ; dans l'exemple matériel à grande échelle ci-dessous, nous devrons plutôt nous appuyer sur des estimations par réseau de tenseurs.

exact_expvals = []
for t in exact_evolution_times:
    exp_H = expm(-1j * t * hamiltonian.to_matrix())
    initial_state = Statevector(initial_state_circ).data
    time_evolved_state = exp_H @ initial_state

    exact_obs = (
        time_evolved_state.conj()
        @ observable.to_matrix()
        @ time_evolved_state
    ).real
    exact_expvals.append(exact_obs)

Coefficients MPF statiques

Les MPF statiques utilisent des coefficients $x_j$ 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 le Contexte et résolvons pour les coefficients. La matrice $A$ est déterminée par les nombres de pas de Trotter $k_j$, l'ordre $\chi$ de la formule produit, et le fait que la formule soit symétrique ou non (ce qui contrôle les exposants $\eta_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$ (donc $\chi=1$ et $\eta_n = 2 + n$, ce qui donne $\eta_0 = 2,\, \eta_1 = 3$). Le système devient :

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

La première ligne impose l'absence de biais ($\sum_j x_j = 1$) ; la deuxième et la troisième ligne annulent respectivement les termes d'erreur de Trotter dominants en $1/k^2$ et en $1/k^3$ de l'ordre suivant.

Configurer le LSE

Nous utilisons setup_static_lse de qiskit_addon_mpf.static pour assembler la matrice $A$ et le vecteur $b$ du membre de droite décrits ci-dessus. La matrice $A$ dépend non seulement de $k_j$, mais aussi de notre choix de formule produit — en particulier son ordre $\chi$ et le fait qu'elle soit symétrique. L'indicateur symmetric contrôle le motif des exposants $\eta_n$ (les formules symétriques ne produisent que des termes d'erreur de Trotter de puissance paire ; voir Réf. [1]). Notez que, comme le montre la Réf. [2], définir symmetric=True n'est pas strictement nécessaire même lorsque la PF sous-jacente est symétrique — le LSE non symétrique reste valide (il impose des contraintes supplémentaires non nécessaires).

Pour notre exemple, nous avons déjà défini order = 2 et symmetric = False à l'étape 1.

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

Inspectez la matrice construite $A$ et le vecteur $b$ pour confirmer qu'ils correspondent au système écrit ci-dessus.

lse.A

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

lse.b

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

Une fois le LSE en main, nous résolvons pour les coefficients statiques $x_j$ via lse.solve() (il s'agit de la solution directe $x = A^{-1}b$).

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

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

Optimiser pour $x$ en utilisant un modèle exact

Alternativement au calcul de $x = A^{-1}b$, tu peux utiliser setup_exact_model pour construire une instance de cvxpy.Problem qui utilise le LSE comme contraintes et dont la solution optimale donnera $x$.

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

[ 0.04761905 -0.57142857  1.52380952]

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

L1 norm of the exact coefficients: 2.1428571428556378

Optimiser pour $x$ en utilisant un modèle approximatif

Il peut arriver que la norme $L_1$ pour l'ensemble choisi de valeurs $k_j$ soit jugée trop élevée. Si c'est le cas et que tu ne peux pas choisir un ensemble différent de valeurs $k_j$, tu peux utiliser une solution approximative qui contraint la norme $L_1$ à un seuil choisi tout en minimisant $\|Ax - b\|$. Consulte le guide Comment utiliser le modèle approximatif.

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

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

Coefficients MPF dynamiques

Le MPF statique annule les termes d'erreur de Trotter d'une manière indépendante de l'hamiltonien et de l'état, si bien qu'il ne produit pas nécessairement la plus petite erreur d'approximation possible pour un hamiltonien et un état initial donnés. Le MPF dynamique (Réf. [2], [3]) trouve à la place des coefficients dépendants du temps $x_i(t)$ qui minimisent la distance en norme de Frobenius $\|\rho(t) - \mu^D(t)\|_F^2$ à chaque instant $t$. Comme le montre le Contexte, cela nécessite la matrice de recouvrement $M_{ij}(t)$ entre états évolués par Trotter et le recouvrement $L_i(t)$ avec l'état exact — que nous estimons tous deux à l'aide de moteurs de réseau de tenseurs (TeNPy) dans qiskit_addon_mpf.

Pour établir le LSE dynamique, nous avons besoin de trois ingrédients :

1. Une fabrique d'évolueur approximatif que l'addon exécutera pour chaque $k_j$ afin de produire $\rho_{k_j}(t)$ sous forme de MPS/MPO. Nous la construisons à partir de la structure en couches du circuit de Trotter d'ordre $2$ (une couche par slice_by_depth), enveloppée sous forme de LayerwiseEvolver avec les paramètres de troncature de TeNPy.

2. Une fabrique d'évolueur exact qui produit une référence de haute précision $\rho(t)$. Nous utilisons un circuit de Suzuki-Trotter d'ordre quatre à petit pas de temps (dt=0.1, order=4) comme approximation de l'évolution exacte.

3. Une fabrique d'identité et un MPS d'état initial qui amorcent la simulation TeNPy.

La cellule ci-dessous construit la fabrique d'évolueur approximatif.

# Create approximate time-evolution circuits
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)  # collect XX and YY

# Find layers in the circuit
layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)

# Create tensor network models
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

# Create the time-evolution object
approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

avertissement

Les options de LayerwiseEvolver qui déterminent les détails de la simulation par réseau de tenseurs doivent être choisies avec soin afin d'éviter de mettre en place un problème d'optimisation mal défini.

Nous approximons l'état évolué dans le temps exact avec une formule de Suzuki-Trotter du quatrième ordre en utilisant un petit pas de temps dt=0.1. Les paramètres de troncature de TeNPy peuvent affecter la précision, il est donc important d'explorer une gamme de valeurs.

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

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

Enfin, nous définissons une identity_factory qui produit l'état MPO initial et préparons l'état initial de Néel sous forme de MPS correspondant au réseau utilisé par le modèle de Trotter en couches.

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

mps_initial_state = MPS_neel_state(models[0].lat)

Les fabriques étant en place, nous calculons maintenant les coefficients dynamiques à chaque temps d'évolution. Pour chaque $t$, setup_dynamic_lse construit les matrices de recouvrement pertinentes via TeNPy, et setup_frobenius_problem renvoie un cvxpy.Problem qui minimise le coût en norme de Frobenius. Le solveur renvoie des coefficients $x_j(t)$ adaptés à cet instant ; nous les rassemblons dans mpf_dynamic_coeffs_list. Si le solveur échoue pour un $t$ donné, nous revenons à des coefficients nuls afin que la boucle se poursuive.

mpf_dynamic_coeffs_list = []
for t in trotter_times:
    print(f"Computing dynamic coefficients for time={t}")
    lse = setup_dynamic_lse(
        mpf_trotter_steps,
        t,
        identity_factory,
        exact_factory,
        approx_factory,
        mps_initial_state,
    )
    problem, coeffs = setup_frobenius_problem(lse)
    try:
        problem.solve()
        mpf_dynamic_coeffs_list.append(coeffs.value)
    except Exception as error:
        mpf_dynamic_coeffs_list.append(np.zeros(len(mpf_trotter_steps)))
        print(error, "Calculation Failed for time", t)
    print("")

Computing dynamic coefficients for time=0.5

Computing dynamic coefficients for time=0.6

Computing dynamic coefficients for time=0.7

Computing dynamic coefficients for time=0.7999999999999999

Computing dynamic coefficients for time=0.8999999999999999

Computing dynamic coefficients for time=0.9999999999999999

Computing dynamic coefficients for time=1.0999999999999999

Computing dynamic coefficients for time=1.1999999999999997

Computing dynamic coefficients for time=1.2999999999999998

Computing dynamic coefficients for time=1.4

Computing dynamic coefficients for time=1.4999999999999998

Combiner les valeurs d'espérance de Trotter avec les coefficients du MPF

Nous évaluons maintenant $\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)$ pour chaque jeu de coefficients (statique-exact, statique-approximatif et dynamique), propageons les erreurs standard par circuit, et traçons la série temporelle résultante par rapport à la courbe de diagonalisation exacte.

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    trotter_curve, trotter_curve_error = [], []
    for trotter_expvals, trotter_stds in zip(
        mpf_expvals_all_times, mpf_stds_all_times
    ):
        trotter_curve.append(trotter_expvals[k])
        trotter_curve_error.append(trotter_stds[k])

    plt.errorbar(
        trotter_times,
        trotter_curve,
        yerr=trotter_curve_error,
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

# Get expectation values at all times for the static MPF with exact coeffs
exact_mpf_curve, exact_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_exact.value, trotter_stds)
            ]
        )
    )
    exact_mpf_curve_error.append(mpf_std)
    exact_mpf_curve.append(trotter_expvals @ coeffs_exact.value)

plt.errorbar(
    trotter_times,
    exact_mpf_curve,
    yerr=exact_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Exact",
    color="purple",
)

# Get expectation values at all times for the static MPF with approximate coeffs
approx_mpf_curve, approx_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_approx.value, trotter_stds)
            ]
        )
    )
    approx_mpf_curve_error.append(mpf_std)
    approx_mpf_curve.append(trotter_expvals @ coeffs_approx.value)

plt.errorbar(
    trotter_times,
    approx_mpf_curve,
    yerr=approx_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Approx",
    color="orange",
)

# Get expectation values at all times for the dynamic MPF
dynamic_mpf_curve, dynamic_mpf_curve_error = [], []
for trotter_expvals, trotter_stds, dynamic_coeffs in zip(
    mpf_expvals_all_times, mpf_stds_all_times, mpf_dynamic_coeffs_list
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(dynamic_coeffs, trotter_stds)
            ]
        )
    )
    dynamic_mpf_curve_error.append(mpf_std)
    dynamic_mpf_curve.append(trotter_expvals @ dynamic_coeffs)

plt.errorbar(
    trotter_times,
    dynamic_mpf_curve,
    yerr=dynamic_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Dynamic MPF",
    color="pink",
)

# Exact expectation values
plt.plot(
    exact_evolution_times,
    exact_expvals,
    color="red",
    linestyle="--",
    label="Exact time-evolution",
)

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

Le graphique ci-dessus illustre l'interaction entre l'erreur de Trotter et l'erreur d'échantillonnage.

• Erreur de Trotter. Les formules produits individuelles (marqueurs gris) s'écartent de plus en plus de la courbe exacte à mesure que le temps augmente. Le circuit $k=1$ présente le plus grand écart et est le moins profond, mais il se trouve déjà dans le régime où $t/k \gtrsim 1$, si bien que le terme d'erreur dominant en $1/k^{2}$ est important. Les combinaisons MPF (marqueurs colorés) annulent plusieurs de ces termes d'erreur de Trotter dominants, si bien qu'elles suivent la courbe exacte de beaucoup plus près qu'un seul circuit $k_j$. L'écart restant reflète les termes de Trotter d'ordre supérieur que le MPF n'annule pas : un MPF statique d'ordre $2$, $r=3$ n'élimine que les deux premiers ordres d'erreur, et pour des $t/k_{\min}$ élevés, la queue non annulée finit par dominer — le MPF ne garantit donc pas que des circuits très peu profonds restent précis à des instants arbitraires.

• 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 d'erreurs standard indépendantes par circuit $\sigma_{k_j}$ donne une variance totale $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$. Ainsi, plus $\|x\|_2$ est grand (et en pratique $\|x\|_1$, que nous contrôlons), plus il faut de shots pour atteindre une incertitude cible donnée. C'est le compromis derrière l'option de solveur approximatif dans le Contexte : nous plafonnons $\|x\|_1$ pour garder ce surcoût gérable. Fait crucial, contrairement à l'erreur de Trotter, l'erreur d'échantillonnage diminue en $1/\sqrt{N_{\text{shots}}}$, si bien qu'elle peut toujours être réduite en dépensant plus de shots.

Dans l'exemple matériel à grande échelle ci-dessous, le bruit matériel apparaît comme une source d'erreur supplémentaire sur chaque $\langle A \rangle_{k_j}$, qui est de même amplifiée par les coefficients du MPF. Nous verrons comment l'atténuation des erreurs interagit avec les MPF dans cette section.

Exemple matériel à grande échelle

Dans cette section, nous faisons passer le problème à une échelle qui dépasse ce qui peut être simulé exactement. Nous reproduisons certains des résultats présentés dans la Réf. [3], en utilisant une chaîne XXZ à 50 qubits au temps $t = 3$. Nous suivons le même flux de travail en quatre étapes que l'exemple à petite échelle, en ciblant maintenant du matériel quantique réel avec atténuation d'erreurs. Comme dans le modèle, chaque étape est indiquée en ligne dans le code, et une seule étape peut s'étendre sur plusieurs cellules lorsque ses sorties intermédiaires méritent d'être examinées. La correspondance reflète l'exemple à petite échelle : définir un hamiltonien, choisir les paramètres de Trotter, calculer les coefficients du MPF (statiques et dynamiques), et construire les circuits. Les différences clés sont :

• Un hamiltonien XXZ sur 50 sites avec des couplages aléatoires tirés de $\mathcal{U}(0.5, 1.5)$ (Réf. [3]).

• Une formule de Trotter du second ordre symétrique avec $k_j = [3, 4, 6]$ (donc $\chi=1$, symmetric=True).

• Un unique temps d'évolution fixe $t = 3$. Avec $k_{\min}=3$, cela donne $t/k_{\min}=1$, ce qui maintient les constituants peu profonds dans le régime de convergence de Trotter où le modèle d'erreur dominante sur lequel le MPF s'appuie est valide.

• Une exécution de comparaison à circuit unique supplémentaire avec $k = 10$ pas de Trotter, utilisée comme référence. Nous avons choisi $k = 10$ car sa profondeur à deux qubits sur le matériel est plus grande que celle du constituant MPF le plus profond ($k_{\max}=6$) plus le surcoût lié à l'exécution de plusieurs circuits MPF — suffisamment profond pour être limité par le bruit, ce qui est le régime dans lequel la combinaison MPF est censée surpasser la référence à circuit unique. Il s'agit d'une comparaison de type « circuit unique profond » contre la combinaison MPF, et non d'un circuit ciblant l'erreur de Trotter effective du MPF (ce qui nécessiterait beaucoup plus de pas).

Notez que même si nous en sommes encore à l'étape 1 ici (correspondance 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 de $t$ mais pas des mesures quantiques, donc ils peuvent être calculés à tout moment avant l'étape 4. Nous le faisons maintenant afin de regrouper toute la configuration spécifique au MPF au même endroit.

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

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

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

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

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

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

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

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

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

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

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

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

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

mps_initial_state = MPS_neel_state(models[0].lat)

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

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

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

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

Nous optimisons maintenant les circuits pour le backend choisi. Nous utilisons le gestionnaire de passes prédéfini de Qiskit avec optimization_level=3, qui sélectionne automatiquement un bon ensemble de qubits physiques et achemine chaque circuit sur la topologie de l'appareil.

# -------------------------Step 2-------------------------
service = QiskitRuntimeService()
# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=L)
backend = service.backend("ibm_fez")
print(backend)

transpiler = generate_preset_pass_manager(
    optimization_level=3, backend=backend
)
transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

<IBMBackend('ibm_fez')>

Exécuter des circuits plus profonds sur du matériel réel nécessite une atténuation d'erreurs agressive. Nous activons le découplage dynamique, le twirling de portes et de mesures, l'atténuation des erreurs de mesure, et l'extrapolation à bruit nul (ZNE). Notez que les facteurs de bruit ZNE que nous utilisons ici (1, 1.2, 1.4) sont plus petits que dans un scénario de circuit peu profond, car les constituants MPF plus profonds sont déjà proches du seuil de bruit et de fortes amplifications de bruit les feraient dépasser le point où l'extrapolation ZNE reste fiable.

Nous soumettons les quatre circuits (trois constituants MPF à $k_j = [3, 4, 6]$ plus la référence $k = 10$) dans un seul job Estimator.

# -------------------------Step 3-------------------------
estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Error suppression/mitigation
estimator.options.dynamical_decoupling.enable = True
estimator.options.twirling.enable_gates = True
estimator.options.twirling.enable_measure = True
estimator.options.twirling.num_randomizations = "auto"
estimator.options.twirling.strategy = "active-accum"
estimator.options.resilience.measure_mitigation = True
estimator.options.experimental.execution_path = "gen3-turbo"

estimator.options.resilience.zne_mitigation = True
estimator.options.resilience.zne.noise_factors = (1, 1.2, 1.4)
estimator.options.resilience.zne.extrapolator = "linear"

estimator.options.environment.job_tags = ["TUT_MPF"]

job_50 = estimator.run(
    [
        (circ, observable)
        for circ, observable in zip(transpiled_circuits, isa_observables)
    ]
)

Nous récupérons les valeurs d'espérance et les écarts types par circuit à partir du résultat du job, puis les combinons avec chaque jeu 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 une variance propagée $\sigma^2 = \sum_j x_j^2 \sigma_{k_j}^2$.

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

print(evs)
print(std)

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

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

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

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

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

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

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

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

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

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

Quelques observations sur les résultats matériels ci-dessus :

• Aller plus profond n'est pas gratuit sur le matériel. Les références à circuit unique le montrent directement : le circuit $k = 6$ est essentiellement exact ($-0.256$ contre la référence $-0.244$), pourtant la référence plus profonde $k = 10$ est pire ($-0.061$, un écart d'environ $0.18$), et non meilleure. Une fois que l'erreur de Trotter est déjà faible, ajouter des pas ne fait surtout qu'approfondir le circuit et accumuler davantage de bruit de porte et de décohérence. C'est précisément le régime pour lequel les MPF sont conçus : atteindre la précision d'un circuit profond en n'utilisant que des constituants peu profonds.

• Un MPF à norme faible surpasse le circuit unique profond. Le MPF statique-approximatif (plafonné à $\|x\|_1 \approx 2$) atterrit à $-0.259$, à environ $0.015$ de la référence et bien plus proche que la référence $k = 10$. Le MPF dynamique ($-0.127$) surpasse également confortablement cette référence. Les deux ne combinent que les circuits peu profonds $k_j = [3, 4, 6]$, et pourtant retrouvent une réponse que le circuit unique profond n'a pas pu obtenir.

• La norme des coefficients importe plus que l'optimalité mathématique. Le MPF statique-exact a $\|x\|_1 = 4.66$ et est le pire estimateur de tous ($-0.567$, un écart de plus de $0.3$) : la grande norme des coefficients amplifie le bruit de porte résiduel, la décohérence et l'erreur ZNE sur chaque $\langle A \rangle_{k_j}$ à peu près du même facteur, ce qui l'emporte sur l'annulation de l'erreur de Trotter qu'il apporte. Plafonner la norme (le solveur statique-approximatif, $\|x\|_1 \approx 2$) supprime ce problème et donne la meilleure estimation — même si ses coefficients n'annulent plus exactement l'erreur de Trotter dominante.

• Des circuits peu profonds individuels peuvent encore être compétitifs. Le seul constituant $k = 6$ ($-0.256$) est lui-même essentiellement exact ici — sur cette exécution, il est même légèrement plus proche que le MPF statique-approximatif. Le hic, c'est que l'on ne sait pas à l'avance quel $k$ unique se situe dans la zone idéale de « convergé mais pas encore limité par le bruit », et le choix apparemment sûr consistant simplement à aller plus profond ($k = 10$) pour garantir la convergence de Trotter est précisément celui qui échoue. Le MPF offre une combinaison rigoureuse de circuits peu profonds qui ne nécessite pas de deviner la bonne profondeur.

En pratique, la leçon à retenir est que sur du matériel, les MPF doivent être associés à une atténuation d'erreurs solide sur chaque $\langle A \rangle_{k_j}$ individuel, la norme $L_1$ des coefficients doit rester modeste (utiliser le solveur approximatif, ou le MPF dynamique), et les pas de Trotter $k_j$ doivent être choisis de sorte que $t/k_{\min} \lesssim 1$ — ici $k_{\min} = 3$ à $t = 3$ donne $t/k_{\min} = 1$, maintenant les constituants dans le régime convergent où le modèle d'erreur dominante sur lequel s'appuie le MPF statique est valide. Avec ces choix, les MPF à norme faible ici égalent un circuit unique convergé alors que la référence naïve « aller simplement plus profond » n'y parvient pas, retrouvant l'avantage profondeur-contre-précision montré dans la Réf. [3]. Notez également que les exécutions individuelles sont bruitées — sur une soumission différente du même job (ou un backend différent), l'ordre exact peut changer ; les tendances robustes sont que les MPF à faible $\|x\|_1$ se comportent bien, le MPF statique-exact à grand $\|x\|_1$ est amplifié par le bruit matériel, et le circuit unique trop profond est limité par le bruit.

Prochaines étapes

Recommandations

Si ce travail t'a intéressé, tu pourrais être intéressé par les ressources suivantes :

• Comment choisir les pas de Trotter pour un MPF — conseils pratiques pour sélectionner des valeurs de $k_j$ afin d'éviter les instabilités
• Comment utiliser le modèle approximatif — ajuster la contrainte de norme $L_1$ et les options du solveur pour le MPF statique approximatif
• Référence de l'API qiskit-addon-mpf — documentation complète pour les modules statique, dynamique et backend

Références

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

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

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