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

Estimation de la durée d'exécution : 37 minutes sur un processeur Heron (REMARQUE : il s'agit uniquement d'une estimation. (La durée d'exécution peut varier.)

Résultats d'apprentissage

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

• Comment modéliser une chaîne de spins de Heisenberg sous forme d'hamiltonien quantique à l'aide de 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 le matériel quantique d’ IBM® à l’aide des primitives et des sessions d’ IBM Quantum

Prérequis

Nous vous recommandons de vous familiariser avec les sujets suivants :

• Notions de base sur l’information quantique
• Introduction aux modèles Qiskit
• Conception d'algorithmes variationnels

Arrière-plan

La chaîne de spins 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. Il décrit un réseau unidimensionnel de spins quantiques en interaction, dans lequel les spins voisins sont couplés par des interactions d'échange. L'hamiltonien du modèle de Heisenberg isotrope soumis à 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$ porte sur les paires de voisins immédiats, $J_x = J_y = J_z = 0.5$ sont les constantes de couplage d'échange (isotropes 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 de manière aléatoire sur l'intervalle $[-1, 1]$. Il convient de noter 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 premiers qubits de l' $N$, qui peuvent ne pas former une chaîne linéaire stricte selon la topologie du dispositif.

Comprendre l'énergie de l'état fondamental de cet hamiltonien revêt une importance fondamentale en physique. L'état fondamental contient des informations sur les transitions de phase quantiques, la structure d'intrication et l'ordre magnétique. En règle générale, le calcul de l'énergie exacte de l'état fondamental devient impossible à réaliser lorsque le nombre de spins augmente, car la dimension de l'espace de Hilbert croît de manière exponentielle selon la formule $2^N$ pour $N$ spins. Cela en fait un candidat tout désigné 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. Le procédé consiste à préparer un état quantique paramétré $|\psi(\theta)\rangle$ (appelé « ansatz ») sur un ordinateur quantique, puis à mesurer la valeur d'espérance $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Un optimiseur classique ajuste ensuite les paramètres de manière itérative $\theta$ afin de minimiser cette énergie, en s'appuyant sur le principe variationnel qui garantit que l'énergie mesurée constitue toujours une borne supérieure de l'énergie réelle 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 de qubits uniques et de portes d'intrication. L'optimisation est réalisée à l'aide de l'algorithme SPSA (Simultaneous Perturbation Stochastic Approximation), particulièrement adapté au matériel quantique sujet au bruit, car il estime les gradients en n'utilisant que deux évaluations de fonction par itération, quel que soit le nombre de paramètres.

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.44 ou version ultérieure (pip install qiskit-ibm-runtime)

Configuration

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 passons en revue chaque étape du modèle Qiskit à petite échelle, en expliquant les éléments clés au fur et à mesure que nous mettons en place le flux de travail.

Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique

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

Élaborez un ansatz et un hamiltonien permettant de modéliser une chaîne de Heisenberg à 10 spins. Au cours de cette étape, nous allons construire un hamiltonien de Heisenberg à 10 spins sur la carte de couplage du backend le moins sollicité 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")

Output:

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

• Entrée : Circuit abstrait, observable
• Résultat : Circuit cible et observable, optimisés pour la QPU sélectionnée

Utilisez la fonction generate_preset_pass_manager de Qiskit pour générer automatiquement une routine d'optimisation pour notre circuit par rapport à la QPU sélectionnée. Nous choisissons optimization_level=3, qui offre le plus haut niveau d'optimisation des gestionnaires de passage prédéfinis. Nous incluons également 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)

Output:

Étape 3 : Exécutez à l'aide d' Qiskit primitives

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

Réduire au minimum l'énergie estimée de l'état fondamental du système en optimisant les paramètres du circuit. Utilisez la primitive Estimator « IBM Quantum » pour évaluer la fonction de coût pendant l'optimisation.

Comme nous avons optimisé le circuit pour le backend à l'étape 2, nous pouvons éviter la transpilation sur le serveur Quantum Compute en définissant et skip_transpilation=True en transmettant le circuit optimisé. Pour cette démonstration, nous allons effectuer l'exécution sur un QPU à l'aide de primitives « IBM Quantum ». Pour utiliser des primitives basées qiskit sur des vecteurs d'état, remplacez le bloc de code utilisant les primitives de type « IBM Quantum » 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. Nous allons maintenant vous en présenter brièvement les principes et vous fournir le code permettant de mettre en œuvre l'algorithme SPSA à l'aide de Qiskit v2.0.

Présentation de SPSA

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

$[\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, étant donné qu'une perturbation aléatoire est appliquée lors de l'estimation du gradient, on peut s'attendre à ce que de légers écarts par rapport aux valeurs exactes d' $f$, dus au bruit, puissent être tolérés et pris en compte. En effet, l'algorithme SPSA est particulièrement réputé pour sa résistance au bruit et ne nécessite que deux appels matériels par itération. Il s'agit donc de l'un des optimiseurs les plus prisés pour la mise en œuvre d'algorithmes variationnels.

Dans ce tutoriel, les hyperparamètres pour l'itération de l' $i^{th}$, $a_i$ et $c_i$, sont calculés comme suit :

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

où les valeurs des constantes sont respectivement de l'ordre de $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 adéquat des hyperparamètres est nécessaire pour obtenir de bonnes performances avec le 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)

C'est ici que nous définissons le maxiter = 50. Il convient de noter que, comme chaque itération nécessite deux appels à la fonction pour calculer le gradient, le nombre total d'appels à la fonction sera de $2 \times \text{maxiter}$. La valeur de maxiter peut être augmentée à une valeur supérieure pour obtenir une meilleure estimation de l'énergie.

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

Output:

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

Étape 4 : Post-traitement et restitution du résultat dans le format classique souhaité

• Données d'entrée : estimations de l'énergie de l'état fondamental pendant l'optimisation
• Résultat : Énergie estimée de l'état fondamental

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

Output:

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Output:

Exemple de matériel à grande échelle

Ce tutoriel ne comprend pas d'exemple concret à grande échelle. À mesure que le nombre de qubits augmente, l'optimisation VQE se heurte à des difficultés majeures en raison du phénomène de plateau stérile : le gradient de la fonction de coût s'annule de manière exponentielle avec la taille du système, rendant l'optimisation pratiquement impossible pour les grands circuits. Si l'on ajoute à cela le bruit matériel, cela signifie que l'extension de la VQE à des chaînes de spin plus longues ne permet pas d'obtenir des résultats reproductibles de manière fiable. Pour découvrir des approches permettant de surmonter ces limites, consultez la section « Étapes suivantes » ci-dessous.

Le défi

Maintenant que vous disposez d'une implémentation VQE fonctionnelle pour la chaîne de Heisenberg, essayez ce qui suit :

1. Jouez sur la profondeur de l'approximation : modifiez le reps paramètre dans efficient_su2 (par exemple, essayez reps=1 et reps=3). Dans quelle mesure la profondeur de l'approximation influe-t-elle sur l'estimation de l'énergie de l'état fondamental et sur la vitesse de convergence? À partir de quel moment constatez-vous une baisse de rendement ou une instabilité?
2. Régler les hyperparamètres SPSA : ajuster les paramètres du programme de taux d'apprentissage (a, c, alpha, gamma, A) et observer leur incidence sur la convergence. Pouvez-vous trouver une configuration qui converge plus rapidement que les paramètres par défaut utilisés ici?
3. Comparez les topologies de couplage : au lieu d'utiliser la carte de couplage native du backend, essayez de construire une simple chaîne linéaire de voisins les plus proches et comparez les résultats. En quoi la connectivité du matériel physique influe-t-elle sur la profondeur du circuit transpilé et sur l'estimation finale de la consommation énergétique?

Références

[1] Spall, J. C. (2002). Mise en œuvre de l'algorithme de perturbation simultanée pour l'optimisation stochastique. IEEE Transactions on Aerospace and Electronic Systems, 34(3), p. 817-823.

[2] Sahin, M. Emre, et al. (2025). Machine Learning Qiskit : une bibliothèque open source destinée aux tâches d'apprentissage automatique quantique à grande échelle sur du matériel quantique et des simulateurs classiques. arXiv:2505.17756.

Etapes suivantes