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

Stima dell'utilizzo: 37 minuti su un processore Heron (NOTA: si tratta solo di una stima. (La durata potrebbe variare.)

Risultati di apprendimento

Una volta completato questo tutorial, avrai acquisito le seguenti conoscenze:

• Come modellare una catena di spin di Heisenberg come hamiltoniano quantistico utilizzando Qiskit
• Come utilizzare l'ottimizzatore SPSA per stimare l'energia dello stato fondamentale di un sistema quantistico
• Come eseguire flussi di lavoro variazionali su hardware quantistico IBM® utilizzando primitive e sessioni IBM Quantum

Prerequisiti

Si consiglia di approfondire i seguenti argomenti:

• Fondamenti dell'informazione quantistica
• Introduzione ai modelli di Qiskit
• Progettazione di algoritmi variazionali

Sfondo

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, in cui gli spin più vicini sono accoppiati tramite interazioni di scambio. L'hamiltoniano del modello di Heisenberg isotropo con un campo magnetico esterno è dato 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$ si estende sulle coppie di vicini più prossimi, $J_x = J_y = J_z = 0.5$ sono le costanti di accoppiamento di scambio (isotropiche in questo tutorial) e $h_i$ rappresenta un campo magnetico esterno dipendente dal sito. In questo tutorial, i valori del campo magnetico vengono campionati in modo casuale nell'intervallo $[-1, 1]$. Si noti che nell'implementazione riportata di seguito, l'insieme delle coppie "di vicini più prossimi" è determinato dall'accoppiamento nativo del backend hardware tra i primi $N$ qubit, che potrebbero non formare una catena lineare rigorosa a seconda della topologia del dispositivo.

Comprendere l'energia dello stato fondamentale di questo hamiltoniano riveste un'importanza fondamentale in fisica. Lo stato fondamentale codifica informazioni relative alle transizioni di fase quantistiche, alla struttura dell'entanglement e all'ordinamento magnetico. In linea di massima, il calcolo dell’energia esatta dello stato fondamentale diventa impraticabile all’aumentare del numero di spin, poiché la dimensione dello spazio di Hilbert cresce in modo esponenziale secondo la formula $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 hamiltoniano. Il metodo consiste nel preparare uno stato quantistico parametrizzato $|\psi(\theta)\rangle$ (detto "ansatz") su un computer quantistico e nel misurare il valore atteso $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Un ottimizzatore classico regola poi iterativamente i parametri $\theta$ per minimizzare tale energia, sfruttando il principio variazionale che garantisce che l'energia misurata sia sempre un limite superiore dell'energia effettiva dello stato fondamentale.

In questo tutorial utilizziamo l'ansatz efficient_su2 della libreria circuit di Qiskit, che costruisce livelli di rotazioni di singoli qubit e porte di entanglement. L'ottimizzazione viene eseguita utilizzando l'algoritmo SPSA (Simultaneous Perturbation Stochastic Approximation), particolarmente adatto all'hardware quantistico soggetto a rumore poiché stima i gradienti utilizzando solo due valutazioni della funzione per ogni iterazione, indipendentemente dal numero di parametri.

Requisiti

Prima di iniziare questa esercitazione, assicuratevi di aver installato quanto segue:

• Qiskit SDK v2.0 o versioni successive, con supporto alla visualizzazione
• Qiskit Runtime v0.44 o versioni successive (pip install qiskit-ibm-runtime)

Configura

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, esamineremo passo dopo passo il modello Qiskit su piccola scala, illustrando gli elementi chiave man mano che sviluppiamo il flusso di lavoro.

Fase 1: mappare gli input classici su un problema quantistico

• Ingresso: Numero di rotazioni
• Uscita: Ansatz e Hamiltoniana che modella la catena di Heisenberg

Costruire un ansatz e un hamiltoniano che descrivano una catena di Heisenberg a 10 spin. In questa fase, costruiremo un hamiltoniano di Heisenberg a 10 spin sulla mappa di accoppiamento del backend meno carico 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")

Output:

Fase 2: Ottimizzazione del problema per l'esecuzione su hardware quantistico

• Ingresso: Circuito astratto, osservabile
• Uscita: Circuito di destinazione e osservabile, ottimizzato per la QPU selezionata

Utilizzare la funzione generate_preset_pass_manager di Qiskit per generare automaticamente una routine di ottimizzazione per il nostro circuito rispetto alla QPU selezionata. Abbiamo scelto optimization_level=3, che offre il massimo livello di ottimizzazione dei gestori di passaggi preimpostati. Includiamo anche passaggi di programmazione 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)

Output:

Passaggio 3: eseguire utilizzando Qiskit primitives

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

Ridurre al minimo l'energia stimata dello stato fondamentale del sistema ottimizzando i parametri del circuito. Utilizza la primitiva Estimator IBM Quantum per calcolare la funzione di costo durante l'ottimizzazione.

Poiché abbiamo ottimizzato il circuito per il backend nel Passo 2, possiamo evitare di eseguire la transpilazione sul server Quantum Compute impostando e skip_transpilation=True passando il circuito ottimizzato. Per questa demo, eseguiremo il codice su una QPU utilizzando le primitive " IBM Quantum ". Per eseguire il codice con primitive basate qiskit su vettori di stato, sostituire il blocco di codice che utilizza primitive di tipo “ IBM Quantum ” con il blocco commentato.

In questo tutorial utilizziamo l'SPSA (Simultaneous Perturbation Stochastic Approximation), un ottimizzatore basato sul gradiente. Di seguito ne forniremo una breve introduzione e il codice per implementare l'SPSA utilizzando Qiskit v2.0.

Vi presentiamo SPSA

L'approssimazione stocastica a perturbazione simultanea (SPSA) [1] è un algoritmo di ottimizzazione che approssima l'intero vettore del gradiente utilizzando solo due chiamate di funzione ad ogni iterazione. Sia $f:\mathbb{R}^p\rightarrow \mathbb{R}$ la funzione di costo con $p$ i parametri da ottimizzare, e sia $x_i\in \mathbb{R}^p$ il vettore dei parametri al $i^{th}$ passo dell'iterazione. Per calcolare il gradiente, viene creato un vettore casuale $\Delta_i$ di dimensione $p$, in cui ciascun elemento $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$ è campionato uniformemente da $\{-1, 1\}$. Successivamente, ciascun elemento del vettore casuale $\Delta_i$ viene moltiplicato per un valore piccolo $c_i$ per generare una perturbazione casuale. Il gradiente viene quindi 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, dato che durante la stima del gradiente viene applicata una perturbazione casuale, ci si aspetta che piccole deviazioni nei valori esatti dell' $f$ e, dovute al rumore, possano essere tollerate e prese in considerazione. Infatti, l'SPSA è noto soprattutto per la sua resistenza al rumore e richiede solo due operazioni hardware per ogni iterazione. Si tratta, quindi, di uno degli ottimizzatori più apprezzati per l'implementazione di algoritmi variazionali.

In questo tutorial, gli iperparametri per l'iterazione dell' $i^{th}$, $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 fissati rispettivamente a $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ e $\gamma = 0.4$. Questi valori sono stati ricavati da [2]. Per ottenere buone prestazioni dall'SPSA è necessario un adeguato aggiustamento degli iperparametri.

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 il maxiter = 50. Si noti che, poiché ogni iterazione richiede due chiamate alla funzione per calcolare il gradiente, il numero totale di chiamate alla funzione sarà pari a $2 \times \text{maxiter}$. Il valore maxiter può essere aumentato a un valore qualsiasi superiore per ottenere una stima più accurata dell'energia.

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

Output:

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

Fase 4: Post-elaborazione e restituzione del risultato nel formato classico desiderato

• Dati di input: stime dell'energia dello stato fondamentale durante l'ottimizzazione
• Risultato: Energia stimata dello stato fondamentale

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

Output:

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Output:

Esempio di hardware su larga scala

In questo tutorial non è incluso alcun esempio di hardware su larga scala. Con l'aumentare del numero di qubit, il VQE deve affrontare sfide significative a causa del fenomeno del "plateau sterile ": il gradiente della funzione di costo svanisce in modo esponenziale con l'aumentare delle dimensioni del sistema, rendendo l'ottimizzazione praticamente impossibile per i circuiti di grandi dimensioni. Se a ciò si aggiunge il rumore dell'hardware, ciò significa che l'estensione del VQE a catene di spin più lunghe non produce risultati affidabili e riproducibili. Per gli approcci che superano tali limiti, consultare la sezione "Prossimi passi" qui di seguito.

Sfida

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

1. Prova a variare la profondità dell'ansatz: modifica il reps parametro in efficient_su2 (ad esempio, prova reps=1 con e reps=3). In che modo la profondità dell'ansatz influisce sull'energia stimata dello stato fondamentale e sulla velocità di convergenza? A che punto si osservano rendimenti decrescenti o instabilità?
2. Ottimizzare gli iperparametri di SPSA: regolare i parametri della curva di apprendimento (a, c, alpha, gamma, A) e osservare come influiscono sulla convergenza. Riesci a trovare una configurazione che converga più rapidamente rispetto alle impostazioni predefinite utilizzate qui?
3. Confronta le topologie di accoppiamento: invece di utilizzare la mappa di accoppiamento nativa del backend, prova a costruire una semplice catena lineare dei vicini più prossimi e confronta i risultati. In che modo la connettività dell'hardware fisico influisce sulla profondità del circuito transpilato e sulla stima finale del consumo energetico?

Riferimenti

[1] Spall, J. C. (2002). Implementazione dell'algoritmo di perturbazione simultanea per l'ottimizzazione stocastica. IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817-823.

[2] Sahin, M. Emre, et al. (2025). Qiskit Machine Learning : una libreria open source per attività di machine learning quantistico su larga scala su hardware quantistico e simulatori classici. arXiv:2505.17756.

Passi successivi