Diagonalización cuántica de Krylov basada en muestras (SKQD)

Esta lección sobre Diagonalización cuántica de Krylov basada en muestras (SKQD) combina métodos explicados en métodos anteriores. Consiste en un único ejemplo que aprovecha el marco de patrones Qiskit:

• Paso 1: Asignar el problema a circuitos y operadores cuánticos
• Paso 2: Optimización para el hardware de destino
• Paso 3: Ejecutar utilizando primitivas de « IBM Quantum »
• Paso 4: Tratamiento posterior

Un paso importante en el método de diagonalización cuántica basado en muestras es generar vectores de calidad para el subespacio. En la lección anterior, utilizamos el ansatz LUCJ para generar vectores subespaciales para un Hamiltoniano de química. En esta lección, utilizaremos los estados cuánticos de Krylov [1] tal y como se discutió en la lección 2. En primer lugar, repasaremos cómo crear el espacio de Krylov en un ordenador cuántico utilizando operaciones de evolución temporal. A continuación, tomaremos una muestra. Proyectaremos el Hamiltoniano del sistema en el subespacio muestreado y lo diagonalizaremos para estimar la energía del estado base. El algoritmo converge de forma demostrable y eficiente al estado fundamental, bajo los supuestos descritos en la lección 2.

0. El espacio de Krylov

Recordemos que un espacio de Krylov $\mathcal{K}^r$ de orden $r$ es el espacio abarcado por los vectores obtenidos multiplicando las potencias superiores de una matriz $A$, hasta $r-1$, por un vector de referencia $\vert v \rangle$.

Si la matriz $A$ es el hamiltoniano $H$, el espacio correspondiente se denomina espacio de Krylov de potencia $\mathcal{K}_P$. Si $A$ es el operador de evolución temporal generado por el hamiltoniano $U=e^{-iH(dt)}$, el espacio se denomina espacio de Krylov unitario $\mathcal{K}_U$. El subespacio de Krylov de potencia no puede generarse directamente en una computadora cuántica, ya que $H$ no es un operador unitario. En su lugar, podemos utilizar el operador de evolución temporal $U = e^{-iH(dt)}$, que puede demostrarse que ofrece garantías de convergencia similares a las del espacio de Krylov de potencia. Las potencias de $U$ se convierten entonces en diferentes pasos temporales $U^k = e^{-iH(k dt)}$ donde $k = 0, 1, 2, ..., (r-1)$.

1. Asignar el problema a circuitos y operadores cuánticos

En esta lección, consideramos el Hamiltoniano para la cadena antiferromagnética XX-Z spin-1/2 con $L = 22$ sitios con la condición de frontera periódica:

from qiskit.transpiler import CouplingMap
from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian

num_spins = 22
coupling_map = CouplingMap.from_ring(num_spins)
H_op = generate_xyz_hamiltonian(coupling_map, coupling_constants=(0.3, 0.3, 1.0))

Para construir el espacio de Krylov, necesitamos tres ingredientes principales:

1. Una elección de la dimensión de Krylov ( $r$ ) y del paso temporal ( $dt$ ).
2. Un estado inicial (de referencia) (vector $\vert v \rangle$ anterior) con solapamiento polinómico con el estado objetivo (base), donde el estado objetivo es disperso. Este requisito de solapamiento polinómico es el mismo que en el algoritmo de estimación de fase cuántica.
3. Operadores de evolución temporal $U^{k}=e^{-iH(k * dt)}$ ( $k = 0, 1, 2, ..., r-1$ ).

Para un valor elegido de $r$ (y, $dt$ ), crearemos $r$ circuitos cuánticos separados y tomaremos muestras de ellos. Cada circuito cuántico se crea uniendo la representación del circuito cuántico del estado de referencia y el operador de evolución temporal para un valor $k$.

Una dimensión de Krylov mayor mejora la convergencia de la energía estimada. En esta lección fijamos la dimensión en $5$ para ilustrar la tendencia a la convergencia.

La referencia [2] demostró que un paso de tiempo suficientemente pequeño para KQD es $\pi / \vert \vert H \vert \vert$, y que es preferible subestimar este valor que sobreestimarlo. Por otro lado, elegir $dt$ demasiado pequeño conduce a un peor acondicionamiento del subespacio de Krylov, ya que los vectores base de Krylov difieren menos de un paso temporal a otro. Además, aunque se ha demostrado que esta elección de $dt$ es adecuada para la convergencia de SKQD, en este contexto basado en el muestreo la elección óptima de $dt$ en la práctica es un tema de estudio en curso. En esta lección, configuramos $dt = 0.15$.

Además de la dimensión de Krylov y el paso temporal, necesitamos establecer el número de pasos de Trotter para la evolución temporal. Utilizar muy pocos pasos conduce a mayores errores de trotterización, mientras que demasiados pasos conducen a circuitos más profundos. En esta lección, fijamos el número de pasos de Trotter en $6$.

# Set parameters for quantum Krylov algorithm
krylov_dim = 5  # size of krylov subspace
dt = 0.15
num_trotter_steps = 6

A continuación, tenemos que elegir un estado de referencia $\vert \psi \rangle$ que tenga cierto solapamiento con el estado de tierra. Para este Hamiltoniano, utilizamos el estado Neel con alternancia de 1s y 0s $\vert ...101...010...101 \rangle$ como nuestro estado de referencia.

# Prep `Neel` state as the reference state for evolution
from qiskit import QuantumCircuit

qc_state_prep = QuantumCircuit(num_spins)
for i in range(num_spins):
    if i % 2 == 0:
        qc_state_prep.x(i)

Por último, necesitamos asignar el operador de evolución temporal a un circuito cuántico. Esto se hizo en la lección 2, pero aquí vamos a aprovechar los métodos en Qiskit, específicamente un método llamado síntesis. Existen diferentes métodos para sintetizar operadores matemáticos en circuitos cuánticos con puertas cuánticas. Muchas de estas técnicas están disponibles en el módulo de síntesis de Qiskit. Utilizaremos el LieTrotter enfoque de síntesis [3] [4].

from qiskit.circuit import QuantumRegister
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.synthesis import LieTrotter

evol_gate = PauliEvolutionGate(
    H_op, time=(dt / num_trotter_steps), synthesis=LieTrotter(reps=num_trotter_steps)
)  # `U` operator

qr = QuantumRegister(num_spins)
qc_evol = QuantumCircuit(qr)
qc_evol.append(evol_gate, qargs=qr)

circuits = []
for rep in range(krylov_dim):
    circ = qc_state_prep.copy()

    # Repeating the `U` operator to implement U^0, U^1, U^2, and so on, for power Krylov space
    for _ in range(rep):
        circ.compose(other=qc_evol, inplace=True)

    circ.measure_all()
    circuits.append(circ)

circuits[1].decompose().draw("mpl", fold=-1)

Output:

circuits[2].decompose().draw("mpl", fold=-1)

Output:

2. Optimizar para el hardware de destino

Ahora que hemos creado los circuitos, podemos optimizarlos para un hardware de destino. Elegimos una QPU a escala de servicio público.

import warnings

from qiskit import generate_preset_pass_manager
from qiskit_ibm_runtime import QiskitRuntimeService

warnings.filterwarnings("ignore")

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

Ahora, transpilamos los circuitos al backend de destino utilizando un gestor de pases preestablecido.

pm = generate_preset_pass_manager(backend=backend, optimization_level=3)
isa_circuits = pm.run(circuits=circuits)

3. Ejecutar en el hardware de destino

Una vez optimizados los circuitos para su ejecución en hardware, estamos listos para ejecutarlos en el hardware de destino y recoger muestras para la estimación de la energía del estado fundamental.

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
job = sampler.run(isa_circuits, shots=100_000)  # Takes approximately 2m 58s of QPU time

counts_all = [job.result()[k].data.meas.get_counts() for k in range(krylov_dim)]

4. Resultados posteriores al procesamiento

A continuación, agregamos los recuentos para dimensiones de Krylov crecientes de forma acumulativa. Utilizando los recuentos acumulativos, abarcaremos subespacios de dimensión de Krylov creciente y analizaremos el comportamiento de convergencia.

from collections import Counter

counts_cumulative = []
for i in range(krylov_dim):
    counter = Counter()
    for d in counts_all[: i + 1]:
        counter.update(d)

    counts = dict(counter)
    counts_cumulative.append(counts)

Para proyectar y diagonalizar el Hamiltoniano, utilizamos las capacidades de qiskit-addon-sqd. El complemento ofrece funcionalidades para proyectar hamiltonianos basados en cadenas de Pauli en un subespacio y resuelve los valores propios utilizando SciPy.

from qiskit_addon_sqd.counts import counts_to_arrays
from qiskit_addon_sqd.qubit import solve_qubit

En principio, podemos filtrar las cadenas de bits con un patrón incorrecto antes de abarcar el subespacio. Por ejemplo, el estado base para el Hamiltoniano antiferromagnético de esta lección suele tener el mismo número de espines "arriba" y "abajo", es decir, el número de "1 "s en la cadena de bits debe ser exactamente la mitad del número total de bits (espines) en el sistema. La siguiente función filtra de los recuentos las cadenas de bits con un número incorrecto de "1".

# Filters out bitstrings that do not have specified number (`num_ones`) of `1` bits.
def postselect_counts(counts, num_ones):
    filtered_counts = {}
    for bitstring, freq in counts.items():
        if bitstring.count("1") == num_ones:
            filtered_counts[bitstring] = freq

    return filtered_counts

Utilizando cadenas de bits con el número correcto de electrones arriba/abajo, abarcamos subespacios y calculamos valores propios para una dimensión de Krylov creciente. Dependiendo del tamaño del problema y de los recursos clásicos disponibles, puede que tengamos que adoptar el submuestreo (similar a la lección sobre SQD ) para mantener la dimensión del subespacio bajo control. Además, podemos aplicar la noción de recuperación de la configuración de forma similar a la Lección 4. Podemos calcular la ocupación de electrones por sitio a partir de los estados propios reconstruidos y utilizar la información para corregir cadenas de bits con un número incorrecto de electrones arriba/abajo. Dejamos esto como ejercicio para los lectores interesados.

import numpy as np

num_batches = 10
rand_seed = 0
scipy_kwargs = {"k": 2, "which": "SA"}

ground_state_energies = []
for idx, counts in enumerate(counts_cumulative):
    counts = postselect_counts(counts, num_ones=num_spins // 2)
    bitstring_matrix, probs = counts_to_arrays(counts=counts)

    eigenvals, eigenstates = solve_qubit(
        bitstring_matrix, H_op, verbose=False, **scipy_kwargs
    )
    gs_en = np.min(eigenvals)
    ground_state_energies.append(gs_en)

A continuación, trazamos la energía calculada en función de la dimensión de Krylov y la comparamos con la energía exacta. La energía exacta se calcula por separado utilizando un método clásico de fuerza bruta. Podemos ver que la energía estimada del estado base converge con el aumento de la dimensión del espacio de Krylov. Aunque la dimensión de Krylov de $5$ es limitante, los resultados siguen mostrando una convergencia impresionante, que se espera mejore con una dimensión de Krylov mayor [1].

import matplotlib.pyplot as plt

exact_gs_en = -23.934184
plt.plot(
    range(1, krylov_dim + 1),
    ground_state_energies,
    color="blue",
    linestyle="-.",
    label="estimate",
)
plt.plot(
    range(1, krylov_dim + 1),
    [exact_gs_en] * krylov_dim,
    color="red",
    linestyle="-",
    label="exact",
)
plt.xticks(range(1, krylov_dim + 1), range(1, krylov_dim + 1))
plt.legend()
plt.xlabel("Krylov space dimension")
plt.ylabel("Energy")
plt.ylim([-24, -22.50])
plt.title(
    "Estimating Ground state energy with Sample-based Krylov Quantum Diagonalization"
)
plt.show()

Output:

Comprueba tu comprensión

Lee las preguntas siguientes, piensa tus respuestas y haz clic en los triángulos para descubrir las soluciones.

• ¿Qué se podría hacer para mejorar la convergencia en el gráfico anterior?

  Respuesta:

  Aumentar la dimensión de Krylov. En general, también se podría aumentar el número de disparos, pero éste ya es bastante elevado en el cálculo anterior.

• ¿Cuáles son las principales ventajas de SKQD sobre (a) SQD y (b) KQD?

  Respuesta:

  Puede haber otras respuestas válidas, pero las respuestas completas deben incluir lo siguiente:

  (a) SKQD ofrece garantías de convergencia que SQD no ofrece. En SQD, o bien tienes que hacer una suposición muy buena para tu ansatz que tenga un solapamiento excelente con el soporte del estado básico en la base computacional, o bien tienes que introducir un componente variacional en el cálculo para muestrear una familia de ansatz.

  (b) SKQD requiere mucho menos tiempo de QPU, porque evita el costoso cálculo de los elementos de la matriz mediante la prueba Hadamard.

5. Resumen

• La estimación de la energía del estado básico mediante el muestreo de estados base de Krylov es muy adecuada para los modelos de celosía, incluidos los sistemas de espín, los problemas de materia condensada y las teorías gauge de celosía. Este enfoque se escala mucho mejor que VQE, porque no requiere la optimización sobre muchos parámetros en un ansatz variacional como en VQE, o en SQD basado en ansatz heurístico (por ejemplo, el problema de química de la lección anterior).
  • Para mantener baja la profundidad de los circuitos, conviene abordar los problemas de entramado que son susceptibles de hardware tolerante a fallos previos.
• SKQD no incurre en un problema de medición cuántica como en VQE. No hay grupos de operadores de Pauli conmutativos que estimar.
• SKQD es resistente a las muestras ruidosas, ya que se puede utilizar una rutina de poselección específica para cada problema (por ejemplo, filtrar las cadenas de bits que no se ajustan a los patrones propios del problema) o asumir la sobrecarga que supone la diagonalización clásica (es decir, diagonalizar en un subespacio más amplio) para eliminar de forma eficaz el efecto del ruido.

Referencias

[1] Jeffery Yu et al., "Algoritmo cuántico-céntrico para la diagonalización de Krylov basada en muestras" (2025). arxiv:quant-ph/2501.09702.

[2] Ethan N. Epperly, Lin Lin y Yuji Nakatsukasa. "Una teoría de diagonalización de subespacios cuánticos". SIAM Journal on Matrix Analysis and Applications 43, 1263-1290 (2022).

[2] N. Hatano y M. Suzuki, "Finding Exponential Product Formulas of Higher Orders" (2005). arXiv:math-ph/0506007.

[4] D. Berry, G. Ahokas, R. Cleve y B. Sanders, "Efficient quantum algorithms for simulating sparse Hamiltonians" (2006). arXiv:quant-ph/0508139.