Réduire l'erreur de Trotter de la dynamique hamiltonienne avec des formules multi-produits

Dans ce notebook, tu vas apprendre à utiliser une formule multi-produits (MPF) pour obtenir une erreur de Trotter plus faible sur notre observable par rapport à celle engendrée par le circuit de Trotter le plus profond que nous allons réellement exécuter. Tu vas le faire en suivant les étapes d'un motif Qiskit :

• Étape 1 : Mapper vers le problème quantique
  • Initialiser le Hamiltonien de notre problème
  • Utiliser une MPF pour générer les circuits d'évolution temporelle de Trotter
• Étape 2 : Optimiser le problème
  • Ici, nous transpirons nos circuits pour un GenericBackendV2
• Étape 3 : Exécuter les expériences
  • Utiliser un StatevectorEstimator pour simplifier les choses dans ce notebook
• Étape 4 : Reconstruire les résultats
  • Calculer la valeur d'espérance de la MPF

Étape 1 : Mapper vers le problème quantique

1a : Mise en place de notre Hamiltonien

Nous utilisons le modèle d'Ising sur une chaîne de 10 sites :

$$\hat{\mathcal{H}}_{\text{Ising}} = \sum_{i=1}^{9} J_{i,(i+1)} Z_i Z_{(i+1)} + \sum_{i=1}^{10} h_i X_i \, ,$$

où $J$ est la force de couplage entre deux sites et $h$ est le champ magnétique externe. Le paquet qiskit_addon_utils fournit des fonctionnalités réutilisables pour diverses applications.

Son module qiskit_addon_utils.problem_generators fournit des fonctions pour générer des Hamiltoniens de type Heisenberg sur un graphe de connectivité donné. Ce graphe peut être soit un rustworkx.PyGraph, soit un CouplingMap, ce qui le rend facile à utiliser dans des flux de travail centrés sur Qiskit.

Dans ce qui suit, nous créons une simple chaîne de 10 qubits en utilisant la méthode CouplingMap.from_line.

# Added by doQumentation — required packages for this notebook
!pip install -q numpy qiskit qiskit-addon-mpf qiskit-addon-utils rustworkx scipy

from qiskit.transpiler import CouplingMap

# Generate some coupling map to use for this example
coupling_map = CouplingMap.from_line(10, bidirectional=False)

from rustworkx.visualization import graphviz_draw

graphviz_draw(coupling_map.graph, method="circo")

Ensuite, nous générons le SparsePauliOp sur la connectivité fournie avec les constantes désirées.

from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

# Get a qubit operator describing the Ising field model
hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(0.0, 0.0, 1.0),
    ext_magnetic_field=(0.4, 0.0, 0.0),
)
print(hamiltonian)

SparsePauliOp(['IIIIIIIZZI', 'IIIIIZZIII', 'IIIZZIIIII', 'IZZIIIIIII', 'IIIIIIIIZZ', 'IIIIIIZZII', 'IIIIZZIIII', 'IIZZIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIXI', 'IIIIIIIXII', 'IIIIIIXIII', 'IIIIIXIIII', 'IIIIXIIIII', 'IIIXIIIIII', 'IIXIIIIIII', 'IXIIIIIIII', 'XIIIIIIIII'],
              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, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j, 0.4+0.j,
 0.4+0.j, 0.4+0.j, 0.4+0.j])

L'observable que nous allons mesurer est la magnétisation totale, que nous pouvons simplement construire comme indiqué ci-dessous :

from qiskit.quantum_info import SparsePauliOp

L = coupling_map.size()
observable = SparsePauliOp.from_sparse_list([("Z", [i], 1 / L / 2) for i in range(L)], num_qubits=L)
print(observable)

SparsePauliOp(['IIIIIIIIIZ', 'IIIIIIIIZI', 'IIIIIIIZII', 'IIIIIIZIII', 'IIIIIZIIII', 'IIIIZIIIII', 'IIIZIIIIII', 'IIZIIIIIII', 'IZIIIIIIII', 'ZIIIIIIIII'],
              coeffs=[0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j, 0.05+0.j,
 0.05+0.j, 0.05+0.j, 0.05+0.j])

1b : Formules multi-produits

Les MPFs réduisent l'erreur de Trotter de la dynamique hamiltonienne grâce à une combinaison pondérée de plusieurs exécutions de circuits.

Pour rendre cela plus concret, nous définissons une MPF comme :

$$\mu(t) = \sum_{j} x_j \rho^{k_j}_{j}\left(\frac{t}{k_j}\right) + \text{some remaining Trotter error} \, ,$$

où $x_j$ sont nos coefficients de pondération, $\rho^{k_j}_j$ est la matrice densité correspondant à l'état pur obtenu en faisant évoluer l'état initial avec la formule produit, $S^{k_j}$, impliquant $k_j$ pas de Trotter, et $j$ indexe le nombre de formules produits qui constituent la MPF.

L'idée clé est que l'erreur de Trotter restante est plus petite que celle qu'on obtiendrait en utilisant simplement la plus grande valeur $k_j$ !

Tu peux voir l'utilité de la MPF sous deux angles :

1. Pour un budget fixe de pas de Trotter que tu es capable d'exécuter, tu peux obtenir des résultats avec une erreur de Trotter globalement plus petite.
2. Pour un nombre de pas de Trotter qui aboutit à des circuits profonds, tu peux utiliser la MPF pour trouver plusieurs circuits moins profonds à exécuter qui donnent une erreur de Trotter similaire.

Introduction aux MPFs statiques

Les MPFs statiques sont celles où les valeurs $x_j$ ne dépendent PAS du temps d'évolution, $t$.

Déterminer les coefficients MPF statiques pour un ensemble donné de valeurs $k_j$ revient à résoudre un système linéaire d'équations : $Ax=b$, où $x$ sont nos coefficients d'intérêt, $A$ est une matrice dépendant de $k_j$ et du type de PF que nous utilisons ($S$), et $b$ est un vecteur de contraintes. Par souci de concision, nous n'allons pas entrer dans les détails ici et te renvoyons plutôt à la documentation de LSE.

Nous pouvons trouver une solution pour $x$ analytiquement comme $x = A^{-1}b$, voir par exemple Carrera Vazquez et al., 2023 ou Zhuk et al., 2023. Cependant, cette solution exacte peut être "mal conditionnée", ce qui entraîne de très grandes normes L1 de nos coefficients $x$, pouvant mener à de mauvaises performances de la MPF. À la place, on peut aussi obtenir une solution approchée qui minimise la norme L1 de $x$ afin de tenter d'optimiser le comportement de la MPF.

Dans ce qui suit, tu vas apprendre à faire tout cela.

Choisir $k_j$

Le choix de $k_j$ est laissé à l'utilisateur final. En principe, n'importe quelles valeurs peuvent être choisies, mais certains $k_j$ entraîneront une amplification du bruit plus importante sur les appareils réels que d'autres choix. Il est donc important d'essayer de trouver de "bonnes" valeurs de $k_j$.

Ici, nous allons simplement choisir des valeurs fixes pour $k_j$. La valeur minimale est motivée par le temps d'évolution cible $t=8.0$ qui nous dit normalement de satisfaire $t/k_{\text{min}} \lt 1$, mais empiriquement nous savons que le fixer égal à $1$ fonctionne aussi généralement. Si tu veux en savoir plus à ce sujet et comment choisir tes autres valeurs $k_j$, réfère-toi au guide correspondant : How to choose the Trotter steps for an MPF.

time = 8.0
trotter_steps = (8, 12, 19)

Mise en place du LSE

Maintenant que nous avons choisi nos $k_j$, nous devons d'abord construire le LSE, $Ax=b$ comme expliqué ci-dessus. La matrice $A$ dépend non seulement de $k_j$ mais aussi de notre choix de formule produit (PF) -- en particulier de son ordre. De plus, on peut prendre en compte si la PF est symétrique ou non (voir Carrera Vazquez et al., 2023), en définissant symmetric=True. Cependant, cela n'est pas obligatoire comme le montrent Zhuk et al., 2023.

Ici, nous allons utiliser une formule de Suzuki-Trotter du second ordre donnant order=2 et nous allons définir symmetric=True.

from qiskit_addon_mpf.static import setup_static_lse

lse = setup_static_lse(trotter_steps, order=2, symmetric=True)
print(lse)

LSE(A=array([[1.00000000e+00, 1.00000000e+00, 1.00000000e+00],
       [1.56250000e-02, 6.94444444e-03, 2.77008310e-03],
       [2.44140625e-04, 4.82253086e-05, 7.67336039e-06]]), b=array([1., 0., 0.]))

Résoudre analytiquement pour $x$

Comme mentionné précédemment, nous pouvons trouver $x$ analytiquement :

import numpy as np

coeffs_analytical = lse.solve()
print(coeffs_analytical)

[ 0.17239057 -1.19447005  2.02207947]

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

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

Dans la section suivante, il sera clair pourquoi cette interface existe.

from qiskit_addon_mpf.costs import setup_exact_problem

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

[ 0.17239057 -1.19447005  2.02207947]

Comme indicateur pour savoir si une MPF construite avec ces coefficients donnera de bons résultats, nous pouvons utiliser la norme L1 (voir aussi Carrera Vazquez et al., 2023).

print(np.linalg.norm(coeffs_exact.value, ord=1))

3.3889400921655914

Optimiser pour $x$ à l'aide d'un modèle approché

Il peut arriver que la norme L1 pour l'ensemble de valeurs $k_j$ choisi soit jugée trop élevée. Dans ce cas, si tu ne peux pas choisir un autre ensemble de valeurs $k_j$, tu peux utiliser une solution approchée du LSE plutôt qu'une solution exacte.

Pour ce faire, utilise simplement setup_sum_of_squares_problem pour construire une instance cvxpy.Problem différente qui contraint la norme L1 à un seuil choisi tout en minimisant la différence entre $Ax$ et $b$.

from qiskit_addon_mpf.costs import setup_sum_of_squares_problem

model_approx, coeffs_approx = setup_sum_of_squares_problem(lse, max_l1_norm=3.0)
model_approx.solve()
print(coeffs_approx.value)
print(np.linalg.norm(coeffs_approx.value, ord=1))

[-0.40454257  0.57553173  0.8290123 ]
1.8090865903790838

Note que tu as une liberté totale sur la façon de résoudre ce problème d'optimisation, ce qui signifie que tu peux changer le solveur d'optimisation, ses seuils de convergence, et ainsi de suite. Consulte le guide correspondant sur How to use the approximate model.

1c : Mise en place des circuits de Trotter

À ce stade, nous avons trouvé nos coefficients d'expansion, $x$, et il ne reste plus qu'à générer les circuits quantiques de Trotter. Une fois de plus, le module qiskit_addon_utils.problem_generators vient à la rescousse pour faire exactement cela :

from qiskit.synthesis import SuzukiTrotter
from qiskit_addon_utils.problem_generators import generate_time_evolution_circuit

circuits = []
for k in trotter_steps:
    circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(order=2, reps=k),
        time=time,
    )
    circuits.append(circ)

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

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

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

Étape 2 : Optimiser le problème

Normalement, c'est l'étape du motif pendant laquelle tu optimises tes circuits pour l'exécution sur du matériel. Ici, comme nous n'utilisons qu'un simulateur sans bruit, nous transpirons simplement notre circuit pour un GenericBackendV2.

from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.transpiler import generate_preset_pass_manager

backend = GenericBackendV2(num_qubits=10)
transpiler = generate_preset_pass_manager(optimization_level=2, backend=backend)

transpiled_circuits = [transpiler.run(circ) for circ in circuits]

Étape 3 : Exécuter les expériences quantiques

Comme expliqué au tout début, nous allons sauter l'étape d'optimisation 2 parce que nous allons simplement calculer les valeurs d'espérance de notre observable cible à l'aide d'un simulateur sans bruit, à savoir le StatevectorEstimator.

from qiskit.primitives import StatevectorEstimator

estimator = StatevectorEstimator()
job = estimator.run([(circ, observable) for circ in transpiled_circuits])
result = job.result()

Étape 4 : Reconstruire les résultats

D'abord, nous extrayons les valeurs d'espérance individuelles obtenues pour chacun des circuits de Trotter :

evs = [res.data.evs for res in result]
print(evs)

[array(0.23799162), array(0.35754312), array(0.38649906)]

Ensuite, nous les recombinons simplement avec nos coefficients MPF pour obtenir les valeurs d'espérance totales de la MPF. Ci-dessous, nous le faisons pour chacune des différentes façons dont nous avons calculé $x$.

print("Analytical    solution:", evs @ coeffs_analytical)
print("Exact model   solution:", evs @ coeffs_exact.value)
print("Approx. model solution:", evs @ coeffs_approx.value)

Analytical    solution: 0.3954847855980006
Exact model   solution: 0.39548478559800204
Approx. model solution: 0.42991214253489807

Enfin, pour ce petit problème, nous pouvons calculer la valeur de référence exacte en utilisant scipy.linalg.expm comme suit :

from scipy.linalg import expm

exp_H = expm(-1j * time * hamiltonian.to_matrix())

initial_state = np.zeros(exp_H.shape[0])
initial_state[0] = 1.0

time_evolved_state = exp_H @ initial_state

exact_obs = time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state
print(exact_obs.real)

0.40060242487899755

Nous pouvons clairement voir que la MPF a réduit l'erreur de Trotter par rapport à celle obtenue avec la PF individuelle la plus profonde avec $k_j=19$. Cependant, nous constatons également que le modèle approché n'est pas parfait, car il a en réalité abouti à une valeur d'espérance moins bonne que la solution exacte. Cela montre l'importance d'utiliser des critères de convergence stricts sur le modèle approché, comme tu l'apprendras dans le guide How to use the approximate model.