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

Tiempo estimado de ejecución: 37 minutos en un procesador Heron (NOTA: Se trata únicamente de una estimación). (El tiempo de ejecución puede variar.)

Resultados del aprendizaje

Una vez completado este tutorial, habrás adquirido los siguientes conocimientos:

• Cómo modelar una cadena de espín de Heisenberg como un hamiltoniano cuántico utilizando Qiskit
• Cómo utilizar el optimizador SPSA para calcular la energía del estado fundamental de un sistema cuántico
• Cómo ejecutar flujos de trabajo variacionales en el hardware cuántico de IBM® utilizando primitivas y sesiones de Qiskit Runtime

Requisitos previos

Se recomienda que te familiarices con estos temas:

• Fundamentos de la información cuántica
• Introducción a los patrones de Qiskit
• Diseño de algoritmos variacionales

En segundo plano

La cadena de espín de Heisenberg es uno de los modelos más estudiados en la física de la materia condensada y el magnetismo cuántico. Describe una red unidimensional de espines cuánticos que interactúan entre sí, en la que los espines más cercanos están acoplados mediante interacciones de intercambio. El hamiltoniano del modelo de Heisenberg isotrópico con un campo magnético externo viene dado por:

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

donde $X_i$, $Y_i$ y $Z_i$ son los operadores de Pauli que actúan en el sitio $i$, la suma $\langle i,j \rangle$ se realiza sobre los pares de vecinos más cercanos, $J_x = J_y = J_z = 0.5$ son las constantes de acoplamiento de intercambio (isotrópicas en este tutorial) y $h_i$ representa un campo magnético externo dependiente del sitio. En este tutorial, los valores del campo magnético se obtienen mediante muestreo aleatorio del intervalo $[-1, 1]$. Cabe señalar que, en la implementación que se muestra a continuación, el conjunto de pares de «vecinos más cercanos» viene determinado por el acoplamiento nativo del backend de hardware entre los primeros $N$ qubits, lo que puede que no forme una cadena lineal estricta dependiendo de la topología del dispositivo.

Comprender la energía del estado fundamental de este hamiltoniano reviste una importancia fundamental en física. El estado fundamental contiene información sobre las transiciones de fase cuánticas, la estructura del entrelazamiento y el orden magnético. Por lo general, calcular la energía exacta del estado fundamental resulta inabarcable a medida que aumenta el número de espines, ya que la dimensión del espacio de Hilbert crece exponencialmente como $2^N$ para $N$ espines. Esto lo convierte en un candidato ideal para la simulación cuántica.

El Variational Quantum Eigensolver (VQE) es un algoritmo híbrido cuántico-clásico diseñado para estimar la energía del estado fundamental de un hamiltoniano. El método consiste en preparar un estado cuántico parametrizado $|\psi(\theta)\rangle$ (denominado «ansatz») en un ordenador cuántico y medir el valor esperado $\langle \psi(\theta) | H | \psi(\theta) \rangle$. A continuación, un optimizador clásico ajusta de forma iterativa los parámetros $\theta$ para minimizar esta energía, aprovechando el principio variacional, que garantiza que la energía medida sea siempre un límite superior de la energía real del estado fundamental.

En este tutorial utilizamos el efficient_su2 enfoque de la biblioteca de circuitos de Qiskit, que construye capas de rotaciones de un solo qubit y puertas de entrelazamiento. La optimización se lleva a cabo mediante el algoritmo de aproximación estocástica por perturbación simultánea (SPSA), que resulta muy adecuado para hardware cuántico con ruido, ya que estima los gradientes utilizando solo dos evaluaciones de la función por iteración, independientemente del número de parámetros.

Requisitos

Antes de empezar este tutorial, asegúrate de que tienes instalado lo siguiente:

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

Configuración

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

Ejemplo a pequeña escala

En esta sección, repasamos cada paso del patrón de Qiskit a pequeña escala, explicando los componentes clave a medida que vamos construyendo el flujo de trabajo.

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

• Entrada: Número de giros
• Salida: Ansatz y Hamiltoniano modelando la cadena de Heisenberg

Elabora un enfoque y un hamiltoniano que modelen una cadena de Heisenberg de 10 espines. En este paso, construiremos un hamiltoniano de Heisenberg de 10 espines sobre el mapa de acoplamiento del backend menos ocupado y prepararemos el efficient_su2 ansatz.

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:

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

• Entrada: Circuito abstracto, observable
• Salida: Circuito objetivo y observable, optimizado 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 mayor nivel de optimización de los gestores de pases preestablecidos. También incluimos ALAPScheduleAnalysis y PadDynamicalDecoupling pases de programación para suprimir los 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)

Output:

Paso 3: Ejecutar utilizando Qiskit primitives

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

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

Dado que en el paso 2 optimizamos el circuito para el backend, podemos evitar la transpilación en el servidor de tiempo de ejecución configurando skip_transpilation=True y pasando el circuito optimizado. Para esta demostración, la ejecutaremos en una QPU utilizando qiskit-ibm-runtime primitivas. Para ejecutar el código con primitivas qiskit basadas en vectores de estado, sustituye el bloque de código que utiliza primitivas de tipo « Qiskit Runtime » por el bloque comentado.

En este tutorial utilizamos la aproximación estocástica por perturbación simultánea (SPSA), que es un optimizador basado en el gradiente. A continuación ofrecemos una breve introducción al tema y proporcionamos el código para implementar el SPSA utilizando Qiskit v2.0.

Presentamos SPSA

La aproximación estocástica por perturbación simultánea (SPSA) [1] es un algoritmo de optimización que aproxima el vector gradiente completo utilizando solo dos llamadas a la función en cada iteración. Sea $f:\mathbb{R}^p\rightarrow \mathbb{R}$ la función de coste, con $p$ los parámetros que se van a optimizar, y $x_i\in \mathbb{R}^p$ el vector de parámetros en el paso $i^{th}$ de la iteración. Para calcular el gradiente, se crea un vector aleatorio $\Delta_i$ de tamaño $p$, en el que cada elemento $\Delta_{ij}$, $\forall$ $j\in \{1,2,...,p\}$ se obtiene mediante muestreo uniforme de $\{-1, 1\}$. A continuación, cada elemento del vector aleatorio $\Delta_i$ se multiplica por un valor pequeño $c_i$ para generar una perturbación aleatoria. A continuación, el gradiente se calcula 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 durante la estimación del gradiente se aplica una perturbación aleatoria, cabe esperar que se puedan tolerar y tener en cuenta las pequeñas desviaciones en los valores exactos de l $f$, debidas al ruido. De hecho, el SPSA destaca especialmente por su resistencia al ruido y solo requiere dos llamadas de hardware por cada iteración. Por lo tanto, es uno de los optimizadores más utilizados para implementar algoritmos variacionales.

En este tutorial, los hiperparámetros para la iteración « $i^{th}$ », « $a_i$ » y « $c_i$ », se calculan como

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

donde los valores constantes se toman como $A = 30$, $\alpha = 0.9$, $a = 0.3$, $c = 0.1$ y $\gamma = 0.4$. Estos valores se han seleccionado de [2]. Es necesario ajustar adecuadamente los hiperparámetros para obtener un buen rendimiento del 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)

Aquí configuramos el 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á de $2 \times \text{maxiter}$. El valor de maxiter puede aumentarse a cualquier valor superior para obtener una mejor estimación de la energía.

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

Output:

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

Paso 4: Procesamiento posterior y devolución del resultado en el formato clásico deseado

• Datos de entrada: estimaciones de la energía del estado fundamental durante la optimización
• Resultado: Energía estimada del estado fundamental

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

Output:

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Output:

Ejemplo de hardware a gran escala

En este tutorial no se incluye ningún ejemplo de hardware a gran escala. A medida que aumenta el número de qubits, la optimización del VQE se enfrenta a importantes retos debido al fenómeno de la meseta estéril: el gradiente de la función de coste desaparece exponencialmente con el tamaño del sistema, lo que hace que la optimización resulte prácticamente inviable para circuitos de gran tamaño. Si a esto le sumamos el ruido del hardware, esto significa que ampliar el VQE a cadenas de espín más largas no produce resultados fiables y reproducibles. Para conocer los enfoques que superan estas limitaciones, consulta la sección «Próximos pasos» que figura a continuación.

Reto

Ahora que ya tienes una implementación de VQE que funciona para la cadena de Heisenberg, prueba lo siguiente:

1. Prueba con diferentes profundidades de aproximación: modifica el reps parámetro en efficient_su2 (por ejemplo, prueba con reps=1 y reps=3). ¿Cómo influye la profundidad del enfoque en la estimación de la energía del estado fundamental y en la velocidad de convergencia? ¿En qué momento se observan rendimientos decrecientes o inestabilidad?
2. Ajustar los hiperparámetros de SPSA: Modifica los parámetros del programa de tasa de aprendizaje (a, c, alpha, gamma, A) y observa cómo influyen en la convergencia. ¿Puedes encontrar una configuración que converja más rápido que los valores predeterminados que se utilizan aquí?
3. Compara diferentes topologías de acoplamiento: en lugar de utilizar el mapa de acoplamiento nativo del backend, intenta construir una cadena lineal simple de «vecinos más cercanos» y compara los resultados. ¿Cómo afecta la conectividad del hardware físico a la profundidad del circuito transpilado y a la estimación final de la energía?

Referencias

[1] Spall, J. C. (2002). Aplicación del algoritmo de perturbación simultánea para la optimización estocástica. IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817-823.

[2] Sahin, M. Emre, et al. (2025). Qiskit Machine Learning : una biblioteca de código abierto para tareas de aprendizaje automático cuántico a gran escala en hardware cuántico y simuladores clásicos. arXiv:2505.17756.

Próximos pasos