Estimación de la energía del estado fundamental de la cadena de Heisenberg con VQE

Estimación de uso: 37 minutos en un procesador Heron (NOTA: Esto es solo una estimación. Tu tiempo de ejecución puede variar.)

Contexto

Este tutorial muestra cómo construir, desplegar y ejecutar un flujo de trabajo de desarrollo llamado patrón de Qiskit para simular una cadena de Heisenberg y estimar su energía del estado fundamental utilizando el optimizador SPSA.

Requisitos

Antes de comenzar este tutorial, asegúrate de tener instalado lo siguiente:

• Qiskit SDK v2.0 o posterior, con soporte de visualización
• Qiskit Runtime v0.44 o posterior (pip install qiskit-ibm-runtime)

Configuración

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

Paso 1: Mapear las entradas clásicas a un problema cuántico

• Entrada: Número de espines
• Salida: Ansatz y hamiltoniano que modelan la cadena de Heisenberg

Construye un ansatz y un hamiltoniano que modelen una cadena de Heisenberg de 10 espines. Primero, importamos algunos paquetes genéricos y creamos algunas funciones auxiliares.

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

Paso 2: Optimizar el problema para la ejecución en hardware cuántico

• Entrada: Circuito abstracto, observable
• Salida: Circuito y observable objetivo, optimizados para la QPU seleccionada

Utiliza la función generate_preset_pass_manager de Qiskit para generar automáticamente una rutina de optimización para nuestro circuito con respecto a la QPU seleccionada. Elegimos optimization_level=3, que proporciona el nivel más alto de optimización de los gestores de pases predefinidos. También incluimos los pases de planificación ALAPScheduleAnalysis y PadDynamicalDecoupling para suprimir errores de decoherencia.

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)

Paso 3: Ejecutar utilizando primitivas de Qiskit

• Entrada: Circuito y observable objetivo
• Salida: Resultados de la optimización

Minimiza la energía estimada del estado fundamental del sistema optimizando los parámetros del circuito. Utiliza la primitiva Estimator de Qiskit Runtime para evaluar la función de costo durante la optimización.

Como ya optimizamos el circuito para el backend en el Paso 2, podemos evitar la transpilación en el servidor de Runtime estableciendo skip_transpilation=True y pasando el circuito optimizado. Para esta demostración, ejecutaremos en una QPU utilizando las primitivas de Qiskit Runtime. Para ejecutar con las primitivas basadas en statevector de qiskit, reemplaza el bloque de código que utiliza las primitivas de Qiskit Runtime con el bloque comentado.

En este tutorial usamos Simultaneous Perturbation Stochastic Approximation (SPSA), que es un optimizador basado en gradiente. A continuación damos una breve introducción y el código para implementar SPSA con Qiskit v2.0.

Introduciendo SPSA

Simultaneous Perturbation Stochastic Approximation (SPSA) [1] es un algoritmo de optimización que aproxima el vector de gradiente completo usando solo dos llamadas a la función en cada iteración. Sea $f:\mathbb{R}^p\rightarrow \mathbb{R}$ la función de costo con $p$ parámetros a optimizar, y $x_i\in \mathbb{R}^p$ el vector de parámetros en el $i^{th}$ paso de la iteración. Para calcular el gradiente, se crea un vector aleatorio $\Delta_i$ de tamaño $p$, donde cada elemento $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$, se muestrea uniformemente de $\{-1, 1\}$. A continuación, cada elemento del vector aleatorio $\Delta_i$ se multiplica por un valor pequeño $c_i$ para crear una perturbación aleatoria. El gradiente se estima entonces como

$[\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, dado que se aplica una perturbación aleatoria durante la estimación del gradiente, se espera que las pequeñas desviaciones en los valores exactos de $f$ provenientes del ruido puedan tolerarse y tenerse en cuenta. De hecho, SPSA es especialmente conocido por ser robusto frente al ruido y requiere solo dos llamadas al hardware por iteración. Por ello, es uno de los optimizadores más preferidos para implementar algoritmos variacionales.

En este tutorial, los hiperparámetros para la $i^{th}$ iteración, $a_i$ y $c_i$, se calculan como $a_i = a/(A + i + 1)^\alpha$ y $c_i = c / (i+1)^\gamma$, donde los valores constantes son $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ y $\gamma = 0.4$. Estos valores se toman de [2]. El ajuste apropiado de los hiperparámetros es necesario para obtener un buen rendimiento de 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)

Aquí establecemos maxiter = 50. Ten en cuenta que, dado que cada iteración requiere dos llamadas a la función para calcular el gradiente, el número total de llamadas a la función será $2 \times \text{maxiter}$. El valor de maxiter puede aumentarse a cualquier valor mayor para una mejor estimación de la energía.

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

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

Paso 4: Post-procesar y devolver el resultado en el formato clásico deseado

• Entrada: Estimaciones de la energía del estado fundamental durante la optimización
• Salida: Energía estimada del estado fundamental

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

Estimated ground state energy: [-2.19621239]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Referencias

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

Próximos pasos

Recomendaciones

Si este trabajo te resultó interesante, puede que te interese el siguiente material:

• Encuentra la energía del estado fundamental de un hamiltoniano disperso siguiendo el tutorial de SQD
• Realiza el curso de Diseño de algoritmos variacionales en IBM Learning