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

Resultados de aprendizaje

Después de completar este tutorial, podrás entender la siguiente información:

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

Prerrequisitos

Se recomienda que te familiarices con estos temas:

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

Contexto

La cadena de espines de Heisenberg es uno de los modelos más ampliamente estudiados en física de la materia condensada y magnetismo cuántico. Describe una red unidimensional de espines cuánticos en interacción, donde los espines vecinos más cercanos están acoplados mediante interacciones de intercambio. El hamiltoniano para el modelo de Heisenberg isotrópico con un campo magnético externo es:

$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 sobre el sitio $i$, la suma $\langle i,j \rangle$ recorre 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 muestrean aleatoriamente del intervalo $[-1, 1]$. Ten en cuenta que en la implementación siguiente, el conjunto de pares "vecinos más cercanos" está determinado por el acoplamiento nativo del backend de hardware entre los primeros $N$ qubits, que puede no formar una cadena lineal estricta dependiendo de la topología del dispositivo.

Entender la energía del estado fundamental de este hamiltoniano es de importancia fundamental en física. El estado fundamental codifica información sobre transiciones de fase cuánticas, estructura de entrelazamiento y ordenamiento magnético. De manera clásica, calcular la energía exacta del estado fundamental se vuelve intratable a medida que crece el número de espines, ya que la dimensión del espacio de Hilbert escala exponencialmente como $2^N$ para $N$ espines. Esto lo convierte en un candidato natural 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. Funciona preparando un estado cuántico parametrizado $|\psi(\theta)\rangle$ (llamado ansatz) en un ordenador cuántico y midiendo el valor esperado $\langle \psi(\theta) | H | \psi(\theta) \rangle$. Un optimizador clásico ajusta iterativamente los parámetros $\theta$ para minimizar esta energía, aprovechando el principio variacional que garantiza que la energía medida es siempre una cota superior de la energía real del estado fundamental.

En este tutorial, usamos el ansatz efficient_su2 de la biblioteca de circuitos de Qiskit, que construye capas de rotaciones de un solo qubit y puertas de entrelazamiento. La optimización se realiza usando el algoritmo Simultaneous Perturbation Stochastic Approximation (SPSA), que es adecuado para hardware cuántico ruidoso porque estima gradientes usando solo dos evaluaciones de función por iteración independientemente del número de parámetros.

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

Ejemplo a pequeña escala

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

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. En este paso, construiremos un hamiltoniano de Heisenberg de 10 espines sobre el mapa de acoplamiento del backend menos ocupado y prepararemos el 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")

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 = \frac{a}{(A + i + 1)^\alpha} \quad \text{and} \quad c_i = \frac{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
        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í 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: -3.03843]

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

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

Ejemplo a gran escala en hardware

Un ejemplo a gran escala en hardware no está incluido en este tutorial. A medida que aumenta el número de qubits, el VQE enfrenta desafíos significativos debido al fenómeno de las mesetas áridas: el gradiente de la función de costo desaparece exponencialmente con el tamaño del sistema, haciendo que la optimización sea prácticamente inviable para circuitos grandes. Combinado con el ruido del hardware, esto significa que escalar VQE a cadenas de espines más grandes no produce resultados fiablemente reproducibles. Para enfoques que superan estas limitaciones, consulta la sección de Próximos pasos más abajo.

Desafío

Ahora que tienes una implementación funcional de VQE para la cadena de Heisenberg, intenta lo siguiente:

1. Experimenta con la profundidad del ansatz: Modifica el parámetro reps en efficient_su2 (por ejemplo, prueba reps=1 y reps=3). ¿Cómo afecta la profundidad del ansatz a la energía estimada del estado fundamental y la velocidad de convergencia? ¿En qué punto observas rendimientos decrecientes o inestabilidad?
2. Ajusta los hiperparámetros de SPSA: Ajusta los parámetros del programa de tasa de aprendizaje (a, c, alpha, gamma, A) y observa cómo impactan la convergencia. ¿Puedes encontrar una configuración que converja más rápido que los valores predeterminados usados aquí?
3. Compara topologías de acoplamiento: En lugar de usar el mapa de acoplamiento nativo del backend, intenta construir una cadena lineal de vecinos más cercanos simple y compara los resultados. ¿Cómo afecta la conectividad del hardware físico a la profundidad del circuito transpilado y la estimación final de la energía?

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:

• Prueba la Diagonalización Cuántica Basada en Muestras (SQD): Como se demostró en este tutorial, VQE enfrenta desafíos a escala debido a las mesetas áridas y el alto costo de medición. IBM ha desarrollado la Diagonalización Cuántica Basada en Muestras (SQD) como una alternativa más escalable. A diferencia de VQE, SQD evita completamente la optimización variacional; en cambio, un ordenador cuántico genera muestras y un ordenador clásico proyecta el hamiltoniano sobre un subespacio generado por esas muestras y lo diagonaliza. Esto proporciona una cota superior de la energía del estado fundamental con significativamente menos mediciones y sin susceptibilidad a las mesetas áridas. Sigue el tutorial de SQD para ver este enfoque en acción.
• Explora el curso de Algoritmos de Diagonalización Cuántica: Profundiza tu comprensión de VQE y SQD, incluyendo sus ventajas e inconvenientes, en el curso de Algoritmos de diagonalización cuántica en IBM Quantum Learning.