Stima dell'energia dello stato fondamentale della catena di Heisenberg con VQE

Stima dell'utilizzo: 37 minuti su un processore Heron (NOTA: Questa è solo una stima. Il tuo tempo di esecuzione potrebbe variare.)

Obiettivi di apprendimento

Dopo aver completato questo tutorial, potrai capire le seguenti informazioni:

• Come modellare una catena di spin di Heisenberg come Hamiltoniana quantistica usando Qiskit
• Come usare l'ottimizzatore SPSA per stimare l'energia dello stato fondamentale di un sistema quantistico
• Come eseguire workflow variazionali su hardware quantistico IBM® usando le primitive e le sessioni di Qiskit Runtime

Prerequisiti

È consigliato familiarizzare con questi argomenti:

• Fondamenti di informazione quantistica
• Introduzione ai pattern Qiskit
• Progettazione di algoritmi variazionali

Contesto

La catena di spin di Heisenberg è uno dei modelli più studiati nella fisica della materia condensata e nel magnetismo quantistico. Descrive un reticolo unidimensionale di spin quantistici interagenti, dove spin vicini sono accoppiati tramite interazioni di scambio. L'Hamiltoniana per il modello di Heisenberg isotropo con un campo magnetico esterno è data da:

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

dove $X_i$, $Y_i$ e $Z_i$ sono gli operatori di Pauli che agiscono sul sito $i$, la somma $\langle i,j \rangle$ è sui vicini più prossimi, $J_x = J_y = J_z = 0.5$ sono le costanti di accoppiamento di scambio (isotrope in questo tutorial) e $h_i$ rappresenta un campo magnetico esterno dipendente dal sito. In questo tutorial, i valori del campo magnetico sono campionati casualmente dall'intervallo $[-1, 1]$. Nota che nell'implementazione seguente, l'insieme delle coppie "vicini più prossimi" è determinato dall'accoppiamento nativo del backend hardware tra i primi $N$ qubit, che potrebbe non formare una catena lineare strettamente lineare a seconda della topologia del dispositivo.

Comprendere l'energia dello stato fondamentale di questa Hamiltoniana è di fondamentale importanza in fisica. Lo stato fondamentale codifica informazioni sulle transizioni di fase quantistiche, la struttura di entanglement e l'ordine magnetico. Classicamente, calcolare l'esatta energia dello stato fondamentale diventa intrattabile con l'aumentare del numero di spin, poiché la dimensione dello spazio di Hilbert scala esponenzialmente come $2^N$ per $N$ spin. Questo lo rende un candidato naturale per la simulazione quantistica.

Il Variational Quantum Eigensolver (VQE) è un algoritmo ibrido quantistico-classico progettato per stimare l'energia dello stato fondamentale di un'Hamiltoniana. Funziona preparando uno stato quantistico parametrizzato $|\psi(\theta)\rangle$ (chiamato ansatz) su un computer quantistico e misurando il valore di aspettazione $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Un ottimizzatore classico regola quindi iterativamente i parametri $\theta$ per minimizzare questa energia, sfruttando il principio variazionale che garantisce che l'energia misurata sia sempre un limite superiore alla vera energia dello stato fondamentale.

In questo tutorial, usiamo l'ansatz efficient_su2 dalla libreria di circuiti di Qiskit, che costruisce strati di rotazioni a qubit singolo e porte di entanglement. L'ottimizzazione viene eseguita usando l'algoritmo SPSA (Simultaneous Perturbation Stochastic Approximation), ben adatto per hardware quantistico rumoroso perché stima i gradienti usando solo due valutazioni della funzione per iterazione, indipendentemente dal numero di parametri.

Requisiti

Prima di iniziare questo tutorial, assicurati di avere installato quanto segue:

• Qiskit SDK v2.0 o successivo, con supporto per la visualizzazione
• Qiskit Runtime v0.44 o successivo (pip install qiskit-ibm-runtime)

Configurazione

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

Esempio su piccola scala

In questa sezione, percorriamo ogni passo del pattern Qiskit su piccola scala, spiegando i componenti chiave mentre costruiamo il workflow.

Passaggio 1: Mappare input classici a un problema quantistico

• Input: Numero di spin
• Output: Ansatz e Hamiltoniana che modellano la catena di Heisenberg

Costruisci un ansatz e un'Hamiltoniana che modellano una catena di Heisenberg a 10 spin. In questo passaggio, costruiremo un'Hamiltoniana di Heisenberg a 10 spin sulla mappa di accoppiamento del backend meno occupato e prepareremo 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")

Passaggio 2: Ottimizzare il problema per l'esecuzione su hardware quantistico

• Input: Circuito astratto, osservabile
• Output: Circuito target e osservabile, ottimizzati per la QPU selezionata

Utilizza la funzione generate_preset_pass_manager di Qiskit per generare automaticamente una routine di ottimizzazione per il nostro circuito rispetto alla QPU selezionata. Scegliamo optimization_level=3, che fornisce il livello più alto di ottimizzazione dei pass manager preimpostati. Includiamo anche i passaggi di scheduling ALAPScheduleAnalysis e PadDynamicalDecoupling per sopprimere gli errori di decoerenza.

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)

Passaggio 3: Eseguire utilizzando le primitive Qiskit

• Input: Circuito target e osservabile
• Output: Risultati dell'ottimizzazione

Minimizza l'energia stimata dello stato fondamentale del sistema ottimizzando i parametri del circuito. Utilizza la primitiva Estimator di Qiskit Runtime per valutare la funzione di costo durante l'ottimizzazione.

Poiché abbiamo ottimizzato il circuito per il backend nel Passaggio 2, possiamo evitare la transpilazione sul server Runtime impostando skip_transpilation=True e passando il circuito ottimizzato. Per questa demo, eseguiremo su una QPU utilizzando le primitive qiskit-ibm-runtime. Per eseguire con le primitive basate su statevector di qiskit, sostituite il blocco di codice che utilizza le primitive Qiskit Runtime con il blocco commentato.

In questo tutorial utilizziamo l'Approssimazione Stocastica per Perturbazione Simultanea (SPSA), che è un ottimizzatore basato sul gradiente. Di seguito ne diamo una breve introduzione e forniamo il codice per implementare SPSA usando Qiskit v2.0.

Introduzione a SPSA

L'Approssimazione Stocastica per Perturbazione Simultanea (SPSA) [1] è un algoritmo di ottimizzazione che approssima l'intero vettore gradiente usando solo due chiamate alla funzione a ogni iterazione. Sia $f:\mathbb{R}^p\rightarrow \mathbb{R}$ la funzione di costo con $p$ parametri da ottimizzare, e $x_i\in \mathbb{R}^p$ il vettore dei parametri all'$i^{mo}$ passo dell'iterazione. Per calcolare il gradiente, viene creato un vettore casuale $\Delta_i$ di dimensione $p$, dove ogni elemento $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$, è campionato uniformemente da $\{-1, 1\}$. Successivamente, ogni elemento del vettore casuale $\Delta_i$ viene moltiplicato per un piccolo valore $c_i$ per creare una perturbazione casuale. Il gradiente viene poi stimato come

$[\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}}.$

Intuitivamente, poiché durante la stima del gradiente viene applicata una perturbazione casuale, ci si aspetta che piccole deviazioni nei valori esatti di $f$ dovute al rumore possano essere tollerate e compensate. In effetti, SPSA è particolarmente noto per essere robusto al rumore e richiede solo due chiamate hardware per ogni iterazione. È quindi uno degli ottimizzatori più preferiti per l'implementazione di algoritmi variazionali.

In questo tutorial, gli iperparametri per l'$i^{mo}$ passo, $a_i$ e $c_i$, vengono calcolati come

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

dove i valori costanti sono $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ e $\gamma = 0.4$. Questi valori sono tratti da [2]. Una corretta regolazione degli iperparametri è necessaria per ottenere buone prestazioni da 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)

Qui impostiamo maxiter = 50. Nota che poiché ogni iterazione richiede due chiamate alla funzione per calcolare il gradiente, il numero totale di chiamate alla funzione sarà $2 \times \text{maxiter}$. Il maxiter può essere aumentato a qualsiasi valore più alto per una migliore stima dell'energia.

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

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

Passaggio 4: Post-elaborare e restituire il risultato nel formato classico desiderato

• Input: Stime dell'energia dello stato fondamentale durante l'ottimizzazione
• Output: Energia stimata dello stato fondamentale

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

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Esempio su hardware su larga scala

Un esempio su hardware su larga scala non è incluso in questo tutorial. Con l'aumentare del numero di qubit, il VQE incontra sfide significative a causa del fenomeno del barren plateau: il gradiente della funzione di costo svanisce esponenzialmente con la dimensione del sistema, rendendo l'ottimizzazione praticamente impossibile per circuiti grandi. Combinato con il rumore hardware, ciò significa che scalare il VQE a catene di spin più grandi non produce risultati affidabilmente riproducibili. Per approcci che superano queste limitazioni, consulta la sezione Prossimi passi di seguito.

Sfida

Ora che hai un'implementazione VQE funzionante per la catena di Heisenberg, prova quanto segue:

1. Sperimenta con la profondità dell'ansatz: Modifica il parametro reps in efficient_su2 (ad esempio, prova reps=1 e reps=3). Come influisce la profondità dell'ansatz sull'energia stimata dello stato fondamentale e sulla velocità di convergenza? A che punto osservi rendimenti decrescenti o instabilità?
2. Regola gli iperparametri SPSA: Modifica i parametri del programma di tasso di apprendimento (a, c, alpha, gamma, A) e osserva come influiscono sulla convergenza. Riesci a trovare una configurazione che converge più velocemente dei valori predefiniti usati qui?
3. Confronta le topologie di accoppiamento: Invece di usare la mappa di accoppiamento nativa del backend, prova a costruire una semplice catena lineare nearest-neighbor e confronta i risultati. In che modo la connettività dell'hardware fisico influisce sulla profondità del circuito transpilato e sulla stima finale dell'energia?

Riferimenti

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

Prossimi passi

Suggerimenti

Se hai trovato interessante questo lavoro, potresti essere interessato al seguente materiale:

• Prova la Sample-based Quantum Diagonalization (SQD): Come dimostrato in questo tutorial, il VQE incontra sfide su scala a causa dei barren plateau e dell'elevato overhead di misurazione. IBM ha sviluppato la Sample-based Quantum Diagonalization (SQD) come alternativa più scalabile. A differenza del VQE, SQD evita completamente l'ottimizzazione variazionale; invece, un computer quantistico genera campioni e un computer classico proietta l'Hamiltoniana su un sottospazio generato da quei campioni e la diagonalizza. Questo fornisce un limite superiore all'energia dello stato fondamentale con significativamente meno misurazioni e senza suscettibilità ai barren plateau. Segui il tutorial SQD per vedere questo approccio in azione.
• Esplora il corso Quantum Diagonalization Algorithms: Approfondisci la tua comprensione sia del VQE che di SQD, inclusi i loro compromessi, nel corso Algoritmi di diagonalizzazione quantistica su IBM Quantum Learning.