Estimation de l'énergie de l'état fondamental de la chaîne de Heisenberg avec VQE

Estimation d'utilisation : 37 minutes sur un processeur Heron (REMARQUE : il s'agit d'une estimation uniquement. Ton temps d'exécution peut varier.)

Objectifs d'apprentissage

Après avoir terminé ce tutoriel, tu pourras t'attendre à comprendre les informations suivantes :

• Comment modéliser une chaîne de spin de Heisenberg en tant que hamiltonien quantique avec Qiskit
• Comment utiliser l'optimiseur SPSA pour estimer l'énergie de l'état fondamental d'un système quantique
• Comment exécuter des workflows variationnels sur du matériel quantique IBM® à l'aide des primitives Qiskit Runtime et des sessions

Prérequis

Il est recommandé de se familiariser avec ces sujets :

• Bases de l'information quantique
• Introduction aux patrons Qiskit
• Conception d'algorithmes variationnels

Contexte

La chaîne de spin de Heisenberg est l'un des modèles les plus étudiés en physique de la matière condensée et en magnétisme quantique. Elle décrit un réseau unidimensionnel de spins quantiques en interaction, où les spins voisins sont couplés par des interactions d'échange. Le hamiltonien pour le modèle de Heisenberg isotrope avec un champ magnétique externe est donné par :

$H = \sum_{\langle i,j \rangle} \left( J_x X_i X_j + J_y Y_i Y_j + J_z Z_i Z_j \right) + \sum_{i} h_i Z_i,$

où $X_i$, $Y_i$ et $Z_i$ sont les opérateurs de Pauli agissant sur le site $i$, la somme $\langle i,j \rangle$ parcourt les paires de voisins les plus proches, $J_x = J_y = J_z = 0.5$ sont les constantes de couplage d'échange (isotrope dans ce tutoriel), et $h_i$ représente un champ magnétique externe dépendant du site. Dans ce tutoriel, les valeurs du champ magnétique sont échantillonnées aléatoirement dans l'intervalle $[-1, 1]$. Note que dans l'implémentation ci-dessous, l'ensemble des paires de « voisins les plus proches » est déterminé par le couplage natif du backend matériel entre les $N$ premiers qubits, ce qui peut ne pas former une chaîne linéaire stricte selon la topologie du dispositif.

Comprendre l'énergie de l'état fondamental de ce hamiltonien revêt une importance fondamentale en physique. L'état fondamental encode des informations sur les transitions de phase quantiques, la structure d'intrication et l'ordre magnétique. Classiquement, calculer l'énergie exacte de l'état fondamental devient insoluble à mesure que le nombre de spins augmente, car la dimension de l'espace de Hilbert croît exponentiellement comme $2^N$ pour $N$ spins. Cela en fait un candidat naturel pour la simulation quantique.

Le Variational Quantum Eigensolver (VQE) est un algorithme hybride quantique-classique conçu pour estimer l'énergie de l'état fondamental d'un hamiltonien. Il fonctionne en préparant un état quantique paramétré $|\psi(\theta)\rangle$ (appelé ansatz) sur un ordinateur quantique et en mesurant la valeur d'expectation $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Un optimiseur classique ajuste ensuite itérativement les paramètres $\theta$ pour minimiser cette énergie, en exploitant le principe variationnel qui garantit que l'énergie mesurée est toujours une borne supérieure à la vraie énergie de l'état fondamental.

Dans ce tutoriel, nous utilisons l'ansatz efficient_su2 de la bibliothèque de circuits de Qiskit, qui construit des couches de rotations à qubit unique et de portes d'intrication. L'optimisation est réalisée à l'aide de l'algorithme d'approximation stochastique par perturbation simultanée (SPSA), qui est bien adapté au matériel quantique bruité car il estime les gradients en utilisant seulement deux évaluations de fonction par itération, quel que soit le nombre de paramètres.

Exigences

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

• Qiskit SDK v2.0 ou version ultérieure, avec le support de visualisation
• Qiskit Runtime v0.44 ou version ultérieure (pip install qiskit-ibm-runtime)

Configuration

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-ibm-runtime

import numpy as np
import matplotlib.pyplot as plt
from typing import Sequence

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import BaseEstimatorV2
from qiskit.circuit.library import XGate
from qiskit.circuit.library import efficient_su2
from qiskit.transpiler import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes.scheduling import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
)
from qiskit_ibm_runtime import QiskitRuntimeService, Session, EstimatorV2

def visualize_results(results):
    plt.plot(results["cost_history"], lw=2)
    plt.xlabel("Number of function evaluations")
    plt.ylabel("Energy")
    plt.show()

Exemple à petite échelle

Dans cette section, nous parcourons chaque étape du patron Qiskit à petite échelle, en expliquant les composants clés au fur et à mesure que nous construisons le workflow.

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

• Entrée : Nombre de spins
• Sortie : Ansatz et hamiltonien modélisant la chaîne de Heisenberg

Construis un ansatz et un hamiltonien qui modélisent une chaîne de Heisenberg à 10 spins. Dans cette étape, nous allons construire un hamiltonien de Heisenberg à 10 spins sur la carte de couplage du backend le moins occupé et préparer l'ansatz efficient_su2.

num_spins = 10
ansatz = efficient_su2(num_qubits=num_spins, reps=2)

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=num_spins, simulator=False
)

coupling = backend.target.build_coupling_map()
reduced_coupling = coupling.reduce(list(range(num_spins)))

edge_list = reduced_coupling.graph.edge_list()
ham_list = []

for edge in edge_list:
    ham_list.append(("ZZ", edge, 0.5))
    ham_list.append(("YY", edge, 0.5))
    ham_list.append(("XX", edge, 0.5))

for qubit in reduced_coupling.physical_qubits:
    ham_list.append(("Z", [qubit], np.random.random() * 2 - 1))

hamiltonian = SparsePauliOp.from_sparse_list(ham_list, num_qubits=num_spins)

ansatz.draw("mpl", style="iqp")

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

• Entrée : Circuit abstrait, observable
• Sortie : Circuit et observable cibles, optimisés pour le QPU sélectionné

Utilise la fonction generate_preset_pass_manager de Qiskit pour générer automatiquement une routine d'optimisation pour notre circuit par rapport au QPU sélectionné. Nous choisissons optimization_level=3, qui fournit le niveau d'optimisation le plus élevé des gestionnaires de passes prédéfinis. Nous incluons également les passes de planification ALAPScheduleAnalysis et PadDynamicalDecoupling pour supprimer les erreurs de décohérence.

target = backend.target
pm = generate_preset_pass_manager(optimization_level=3, target=target)
pm.scheduling = PassManager(
    [
        ALAPScheduleAnalysis(durations=target.durations()),
        PadDynamicalDecoupling(
            durations=target.durations(),
            dd_sequence=[XGate(), XGate()],
            pulse_alignment=target.pulse_alignment,
        ),
    ]
)
isa_ansatz = pm.run(ansatz)
isa_observable = hamiltonian.apply_layout(isa_ansatz.layout)
isa_ansatz.draw("mpl", scale=0.6, style="iqp", fold=-1, idle_wires=False)

Étape 3 : Exécuter à l'aide des primitives Qiskit

• Entrée : Circuit et observable cibles
• Sortie : Résultats de l'optimisation

Minimise l'énergie estimée de l'état fondamental du système en optimisant les paramètres du circuit. Utilise la primitive Estimator de Qiskit Runtime pour évaluer la fonction de coût pendant l'optimisation.

Puisque nous avons optimisé le circuit pour le backend à l'étape 2, nous pouvons éviter la transpilation sur le serveur Runtime en définissant skip_transpilation=True et en passant le circuit optimisé. Pour cette démonstration, nous exécuterons sur un QPU en utilisant les primitives qiskit-ibm-runtime. Pour exécuter avec les primitives basées sur le vecteur d'état de qiskit, remplace le bloc de code utilisant les primitives Qiskit Runtime par le bloc commenté.

Dans ce tutoriel, nous utilisons l'approximation stochastique par perturbation simultanée (SPSA), qui est un optimiseur basé sur le gradient. Ci-après, nous donnons une brève introduction et fournissons le code pour implémenter SPSA avec Qiskit v2.0.

Présentation de SPSA

L'approximation stochastique par perturbation simultanée (SPSA) [1] est un algorithme d'optimisation qui approxime l'intégralité du vecteur gradient en utilisant seulement deux appels de fonction à chaque itération. Soit $f:\mathbb{R}^p\rightarrow \mathbb{R}$ la fonction de coût avec $p$ paramètres à optimiser, et $x_i\in \mathbb{R}^p$ le vecteur de paramètres à la $i^{\text{ième}}$ étape de l'itération. Pour calculer le gradient, un vecteur aléatoire $\Delta_i$ de taille $p$ est créé, où chaque élément $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$, est échantillonné uniformément dans $\{-1, 1\}$. Ensuite, chaque élément du vecteur aléatoire $\Delta_i$ est multiplié par une petite valeur $c_i$ pour créer une perturbation aléatoire. Le gradient est alors estimé comme

$[\nabla f(x_i)]_j \approx \frac{f(x_i + c_i \Delta_i) - f(x_i - c_i \Delta_i)}{2c_i\Delta_{ij}}.$

Intuitivement, puisqu'une perturbation aléatoire est appliquée lors de l'estimation du gradient, on s'attend à ce que de petites déviations dans les valeurs exactes de $f$ dues au bruit puissent être tolérées et prises en compte. En fait, SPSA est particulièrement connu pour sa robustesse au bruit, et ne nécessite que deux appels au matériel par itération. Il est donc l'un des optimiseurs très préférés pour l'implémentation d'algorithmes variationnels.

Dans ce tutoriel, les hyperparamètres pour la $i^{\text{ième}}$ itération, $a_i$ et $c_i$, sont calculés comme

$a_i = \frac{a}{(A + i + 1)^\alpha} \quad \text{and} \quad c_i = \frac{c}{(i+1)^\gamma},$

où les valeurs constantes sont $A = 30$, $\alpha = 0,9$, $a = 0,3$, $c = 0,1$, et $\gamma = 0,4$. Ces valeurs sont tirées de [2]. Un réglage approprié des hyperparamètres est nécessaire pour obtenir de bonnes performances avec SPSA.

def spsa(
    fun, x0, args=(), A=30, alpha=0.9, a=0.3, c=0.1, gamma=0.4, maxiter=100
):
    nparams = len(x0)
    x = np.copy(x0)

    for i in range(maxiter):
        a_i = a / (A + i + 1) ** alpha
        c_i = c / (i + 1) ** gamma
        delta_i = np.random.choice([-1, 1], nparams)

        # two hardware calls
        eval_1 = fun(x + c_i * delta_i, *args)
        eval_2 = fun(x - c_i * delta_i, *args)

        # compute the gradient and update the parameters
        grad = (eval_1 - eval_2) / (2 * c_i) * np.reciprocal(delta_i)
        x = x - a_i * grad

    return x

def cost_func(
    params: Sequence,
    ansatz: QuantumCircuit,
    hamiltonian: SparsePauliOp,
    estimator: BaseEstimatorV2,
    cost_history_dict: dict,
) -> float:
    """Ground state energy evaluation."""
    energy = (
        estimator.run([(ansatz, hamiltonian, [params])]).result()[0].data.evs
    )

    cost_history_dict["iters"] += 1
    cost_history_dict["prev_vector"] = list(params)
    cost_history_dict["cost_history"].append(float(energy[0]))

    print(
        f"Fx Iters. done: {cost_history_dict['iters']} [Current cost: {round(energy[0], 5)}]",
        end="\r",
    )

    return energy

def solve(x0, isa_ansatz, isa_observable, maxiter=150):
    cost_history_dict = {
        "prev_vector": None,
        "iters": 0,
        "cost_history": [],
        "y_min": None,
    }

    # Evaluate the problem using a QPU via Qiskit IBM Runtime
    with Session(backend=backend) as session:
        estimator = EstimatorV2(mode=session)
        estimator.skip_transpilation = True
        estimator.options.environment.job_tags = ["TUT_HSVQE"]
        x_opt = spsa(
            cost_func,
            x0=x0,
            args=(isa_ansatz, isa_observable, estimator, cost_history_dict),
            maxiter=maxiter,
        )

        y_min = cost_func(
            x_opt, isa_ansatz, isa_observable, estimator, cost_history_dict
        )

    return y_min, cost_history_dict

np.random.seed(42)
num_params = ansatz.num_parameters
params = 2 * np.pi * np.random.random(num_params)

Ici nous définissons maxiter = 50. Remarque que chaque itération nécessite deux appels à la fonction pour calculer le gradient, donc le nombre total d'appels de fonction sera $2 \times \text{maxiter}$. Le paramètre maxiter peut être augmenté à n'importe quelle valeur supérieure pour une meilleure estimation de l'énergie.

maxiter = 50
spsa_min, spsa_history = solve(
    params, isa_ansatz, isa_observable, maxiter=maxiter
)

Fx Iters. done: 101 [Current cost: -3.03843]

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

• Entrée : Estimations de l'énergie de l'état fondamental pendant l'optimisation
• Sortie : Énergie estimée de l'état fondamental

print(f"Estimated ground state energy: {spsa_min}")

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Exemple matériel à grande échelle

Un exemple matériel à grande échelle n'est pas inclus dans ce tutoriel. À mesure que le nombre de qubits augmente, VQE rencontre des défis importants en raison du phénomène de plateau stérile : le gradient de la fonction de coût disparaît exponentiellement avec la taille du système, rendant l'optimisation pratiquement irréalisable pour les grands circuits. Combiné au bruit matériel, cela signifie que mettre à l'échelle VQE pour des chaînes de spin plus grandes ne produit pas de résultats fiablement reproductibles. Pour des approches qui surmontent ces limitations, voir la section Prochaines étapes ci-dessous.

Défi

Maintenant que tu as une implémentation VQE fonctionnelle pour la chaîne de Heisenberg, essaie les points suivants :

1. Expérimenter avec la profondeur de l'ansatz : Modifie le paramètre reps dans efficient_su2 (par exemple, essaie reps=1 et reps=3). Comment la profondeur de l'ansatz affecte-t-elle l'énergie estimée de l'état fondamental et la vitesse de convergence ? À quel moment observes-tu des rendements décroissants ou une instabilité ?
2. Régler les hyperparamètres SPSA : Ajuste les paramètres du calendrier de taux d'apprentissage (a, c, alpha, gamma, A) et observe leur impact sur la convergence. Peux-tu trouver une configuration qui converge plus rapidement que les valeurs par défaut utilisées ici ?
3. Comparer les topologies de couplage : Au lieu d'utiliser la carte de couplage native du backend, essaie de construire une simple chaîne linéaire de voisins les plus proches et compare les résultats. Comment la connectivité du matériel physique affecte-t-elle la profondeur du circuit transpilé et l'estimation finale de l'énergie ?

Références

[1] Spall, J. C. (2002). Implementation of the simultaneous perturbation algorithm for stochastic optimization. IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817-823.

[2] Sahin, M. Emre, et al. (2025). Qiskit Machine Learning: an open-source library for quantum machine learning tasks at scale on quantum hardware and classical simulators. arXiv:2505.17756.

Prochaines étapes

Recommandations

Si ce travail t'a intéressé, tu pourrais être intéressé par le matériel suivant :

• Essayer la diagonalisation quantique basée sur des échantillons (SQD) : Comme démontré dans ce tutoriel, VQE fait face à des défis à l'échelle en raison des plateaux stériles et d'un overhead de mesure élevé. IBM a développé la diagonalisation quantique basée sur des échantillons (SQD) comme alternative plus évolutive. Contrairement à VQE, SQD évite entièrement l'optimisation variationnelle ; à la place, un ordinateur quantique génère des échantillons et un ordinateur classique projette le hamiltonien sur un sous-espace délimité par ces échantillons et le diagonalise. Cela fournit une borne supérieure à l'énergie de l'état fondamental avec beaucoup moins de mesures et sans susceptibilité aux plateaux stériles. Suis le tutoriel SQD pour voir cette approche en action.
• Explorer le cours sur les algorithmes de diagonalisation quantique : Approfondis ta compréhension de VQE et de SQD, y compris leurs compromis, dans le cours Algorithmes de diagonalisation quantique sur IBM Quantum Learning.