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

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.)

Acquis d'apprentissage

À l'issue de ce tutoriel, vous devriez être en mesure de comprendre les éléments suivants :

• 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
• 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é?
• Comment calculer les coefficients MPF statiques et dynamiques à l'aide du qiskit_addon_mpf package
• 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

Prérequis

Nous recommandons aux utilisateurs de se familiariser avec les sujets suivants avant de suivre ce tutoriel :

• Méthodes de compilation pour les circuits de simulation hamiltoniens — présentation des circuits de Trotter (formule du produit) dans Qiskit.
• Les formules de produits dans Qiskit, en particulier les SuzukiTrotter classes de synthèse et LieTrotter .
• Qiskit primitives et l'interface Estimator.

Arrière-plan

Qu'est-ce qu'une formule multiproduit?

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.

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

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] ), convergent plus rapidement selon $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$, mais au prix de circuits plus complexes à chaque étape.

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 :

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.

Les formules multiproduits (MPF) [1] 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$ ) :

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, 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.

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

Il existe deux familles de coefficients MPF :

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 :

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] 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\|$.

Les coefficients dynamiques [2], [3] 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 :

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.

Quand utiliser les MPF?

Les fonds de prévoyance (MPF) sont particulièrement avantageux lorsque :

• 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.
• 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 ».
• 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.

Quand les fonds de prévoyance (MPF) ne sont pas forcément la solution

• 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.
• 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.
• 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.

Contenu de ce tutoriel

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.

Exigences

Avant de commencer ce tutoriel, assurez-vous que les éléments suivants sont installés :

• Qiskit SDK v2.0 ou version ultérieure, avec prise en charge de la visualisation
• Qiskit Runtime v0.22 ou plus tard (pip install qiskit-ibm-runtime)
• Simulateur Qiskit Aer (pip install qiskit-aer)
• Module complémentaire MPF Qiskit avec le backend « TeNPy » (pip install "qiskit-addon-mpf[tenpy]")
• Utilitaires complémentaires de Qiskit (pip install qiskit-addon-utils)
• SciPy (pip install scipy)

Configuration

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).

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 : Mettre en correspondance les entrées classiques avec un problème quantique

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 :

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.

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)

Output:

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)

Output:

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 des circuits Trotter

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.

# 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)

Output:

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

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 :

1. Collecte de données au niveau de la simulation hamiltonienne. XXPlusYYGateDans 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.

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.

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écutez à l'aide d' Qiskit primitives

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.

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-traitement et restitution du résultat dans le format classique souhaité

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.

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.

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$ 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$ ).

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 :

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$.

Configurer le LSE

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] ). Il convient de noter que, comme le montre la référence [2], 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).

Dans 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)

Vérifiez que la matrice $A$ et le vecteur $b$ correspondent bien au système présenté ci-dessus.

lse.A

Output:

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

lse.b

Output:

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

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$ ).

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

Output:

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

Optimiser l' $x$ s à l'aide d'un modèle exact

Au lieu de calculer $x = A^{-1}b$, vous pouvez 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)

Output:

[ 0.04761905 -0.57142857  1.52380952]

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

Output:

L1 norm of the exact coefficients: 2.1428571428556378

Optimiser l' $x$ s à l'aide d'un modèle approximatif

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 ».

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

Output:

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

Coefficients MPF dynamiques

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], [3] ) 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.

Pour configurer le LSE dynamique, il nous faut trois éléments :

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.
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.
3. Une « identity factory » et un MPS d'état initial qui servent de base à la simulation « TeNPy ».

La cellule ci-dessous permet de créer l'usine d'évoluteurs approximatifs.

# 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,
    },
)

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.

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 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.

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


mps_initial_state = MPS_neel_state(models[0].lat)

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.

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("")

Output:

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 attendues de Trotter avec les coefficients MPF

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.

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()

Output:

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

• 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.

• 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.

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.

Exemple de matériel à grande échelle

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], 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.

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 :

• Un hamiltonien XXZ à 50 sites avec des couplages aléatoires tirés de $\mathcal{U}(0.5, 1.5)$ (réf. [3] ).
• Une formule de Trotter symétrique du second ordre avec l' $k_j = [3, 4, 6]$ e (donc $\chi=1$, symmetric=True).
• 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.
• 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).

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.

# -------------------------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)

Output:

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 à présent les circuits pour le backend choisi. optimization_level=3Nous 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.

# -------------------------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
]

Output:

<IBMBackend('ibm_fez')>

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.

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.

# -------------------------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 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$.

# -------------------------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)

Output:

[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,
)

Output:

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()

Output:

Quelques remarques concernant les résultats matériels ci-dessus :

• 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.

• 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.

• 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.

• 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.

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]. 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.

Etapes suivantes

Références

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

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

[3] Robertson, N. F., et al. Formules dynamiques multiproduits améliorées par un réseau de tenseurs. arXiv:2407.17405 (2024)