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

Contesto

Questo tutorial mostra come costruire, implementare ed eseguire un workflow di sviluppo chiamato pattern Qiskit per simulare una catena di Heisenberg e stimare la sua energia dello stato fondamentale usando l'ottimizzatore SPSA.

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

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. Per prima cosa, importiamo alcuni pacchetti generici e creiamo alcune funzioni di supporto.

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

service = QiskitRuntimeService(
    channel="ibm_cloud",
    token="<YOUR_API_TOKEN>",  # Replace with your actual API token
    instance="<YOUR_INSTANCE_NAME>",  # Replace with your instance name if needed
)
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 = a/(A + i + 1)^\alpha$ e $c_i = 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
        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: -2.19621]

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: [-2.19621239]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

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:

• Trova l'energia dello stato fondamentale di un'Hamiltoniana sparsa seguendo il tutorial SQD
• Segui il corso Progettazione di algoritmi variazionali su IBM Learning