Genera circuitos de evolución temporal de Krylov (SKQD)

Importante

Los conceptos de esta guía solo están disponibles actualmente en la API de Python. La funcionalidad equivalente estará disponible en la API de C en una versión futura.

La diagonalización cuántica de Krylov basada en muestras (SKQD) es una variante cuántico-céntrica de la diagonalización por subespacios. En lugar de muestrear un único ansatz variacional (como en la diagonalización cuántica basada en muestras simple (SQD)) o un conjunto de circuitos aleatorizados de evolución temporal (como en SqDRIFT), muestrea una base de Krylov: la familia de estados

$$\lvert \psi_k \rangle = e^{-i k \, \Delta t \, H} \lvert \psi_0 \rangle, \qquad k = 0, 1, \dots, D - 1,$$

obtenidos al evolucionar un estado de referencia $\lvert \psi_0 \rangle$ durante múltiplos crecientes de un paso de tiempo fijo $\Delta t$. Las cadenas de bits muestreadas de estos $D$ circuitos pueblan el subespacio de Krylov en el que posteriormente se diagonaliza el hamiltoniano de forma clásica.

Esta guía reproduce el paso de construcción de circuitos del algoritmo SKQD para el modelo de Anderson de una sola impureza (SIAM), una impureza magnética acoplada a un baño sin interacciones, siguiendo la publicación de SKQD.

Nota

El paso de muestreo de SKQD solo funciona si el estado fundamental es disperso en la base computacional: la diagonalización clásica debe poder reconstruirlo a partir de un número manejable de cadenas de bits muestreadas. Como se muestra más abajo, esa dispersión no es automática; es una propiedad de la base de partícula única en la que se expresa el problema. Elegir bien la base es la clave de esta construcción.

1. Construye el hamiltoniano del SIAM como operador fermiónico

El modelo de Anderson de una sola impureza coloca un orbital de impureza con interacción (índice $0$) en contacto con un baño sin interacciones de $L$ sitios dispuestos en una cadena. Su hamiltoniano se divide en una parte de un cuerpo $H_1$ y una parte de dos cuerpos $H_2$,

$$H = H_1 + H_2,$$

donde la parte de un cuerpo reúne la energía local de la impureza (el potencial químico $\mu$), la hibridación impureza-baño $V$ y el salto entre vecinos más cercanos del baño $t$,

$$H_1 = \mu \sum_\sigma a^\dagger_{0\sigma} a_{0\sigma} - V \sum_\sigma \left( a^\dagger_{0\sigma} a_{1\sigma} + \text{h.c.} \right) - t \sum_{\sigma} \sum_{i=1}^{L-1} \left( a^\dagger_{i\sigma} a_{i+1,\sigma} + \text{h.c.} \right),$$

y la parte de dos cuerpos aporta la repulsión de Coulomb local $U$ solo en la impureza,

$$H_2 = U \, a^\dagger_{0\uparrow} a_{0\uparrow} a^\dagger_{0\downarrow} a_{0\downarrow}.$$

Esta guía trabaja en el punto de simetría partícula-hueco $\mu = -U/2$. Estas dos partes se corresponden directamente con las integrales de un cuerpo h1e (salto del baño, hibridación de la impureza, energía local de la impureza) y la única integral de dos cuerpos h2e (el término de Coulomb de la impureza). Primero ensámblalas en la base de posición, en la que el modelo se define de forma natural.

>>> import numpy as np
>>>
>>> norb = 8
>>> nelec = (norb // 2, norb // 2)  # half filling
>>> hopping, onsite, hybridization = 1.0, 10.0, 1.0
>>> chemical_potential = -0.5 * onsite
>>>
>>> # one-body integrals: bath hopping, impurity hybridization, and impurity on-site energy
>>> h1e = np.zeros((norb, norb))
>>> np.fill_diagonal(h1e[:, 1:], -hopping)
>>> np.fill_diagonal(h1e[1:, :], -hopping)
>>> h1e[0, 1] = h1e[1, 0] = -hybridization
>>> h1e[0, 0] = chemical_potential
>>>
>>> # two-body integrals: on-site Coulomb repulsion on the impurity orbital
>>> h2e = np.zeros((norb,) * 4)
>>> h2e[0, 0, 0, 0] = onsite

2. Cambia a la base de momento

La parte del baño del SIAM es una cadena de fermiones libres, por lo que se diagonaliza mediante un cambio a la base de momento. Este es el paso decisivo para SKQD. Como el baño es invariante bajo traslaciones, el estado fundamental no es disperso en la base de posición, pero sí lo es en la base de momento, donde el hamiltoniano del baño es diagonal. La rotación orbital de partícula única que realiza este cambio es álgebra lineal específica del problema. Diagonaliza la matriz de salto del baño sin tocar el orbital de la impureza y después reubica la impureza en un sitio central:

>>> def momentum_basis(norb):
...     """Orbital rotation from the position to the momentum basis."""
...     n_bath = norb - 1
...     hopping_matrix = np.zeros((n_bath, n_bath))
...     np.fill_diagonal(hopping_matrix[:, 1:], -1)
...     np.fill_diagonal(hopping_matrix[1:, :], -1)
...     _, vecs = np.linalg.eigh(hopping_matrix)
...     orbital_rotation = np.zeros((norb, norb))
...     orbital_rotation[0, 0] = 1.0
...     orbital_rotation[1:, 1:] = vecs
...     new_index = n_bath // 2
...     perm = np.r_[1:(new_index + 1), 0, (new_index + 1):norb]
...     return orbital_rotation[:, perm]
...
>>> orbital_rotation = momentum_basis(norb)

A diferencia de un flujo de trabajo variacional, aquí se rotan las integrales del hamiltoniano a la base de momento, de modo que todo el circuito (tanto el hamiltoniano como el estado) vive en la base en la que el estado fundamental es disperso y, por tanto, las cadenas de bits muestreadas son informativas. Rotar las integrales de uno y dos cuerpos es un cambio de base estándar de los tensores de estructura electrónica:

>>> rotation = orbital_rotation.T.conj()
>>> h1e_momentum = np.einsum("ab,Aa,Bb->AB", h1e, rotation, rotation.conj(), optimize="greedy")
>>> h2e_momentum = np.einsum(
...     "abcd,Aa,Bb,Cc,Dd->ABCD",
...     h2e,
...     rotation,
...     rotation.conj(),
...     rotation,
...     rotation.conj(),
...     optimize="greedy",
... )

Como la rotación deja fijo el orbital de la impureza (actúa solo sobre el baño), la interacción local no cambia. h2e_momentum sigue siendo un único término número-número, igual que h2e. El único cambio es que la impureza se ha reubicado en un índice central. Esto se ve directamente en los tensores de integrales, que tienen cada uno una única entrada distinta de cero:

>>> print(f"h2e nonzero at:          {np.argwhere(np.abs(h2e) > 1e-12).tolist()}")
h2e nonzero at:          [[0, 0, 0, 0]]
>>> print(f"h2e_momentum nonzero at: {np.argwhere(np.abs(h2e_momentum) > 1e-12).tolist()}")
h2e_momentum nonzero at: [[3, 3, 3, 3]]

Esta localidad es lo que mantiene barata la síntesis de la interacción más adelante.

Ensambla las integrales de un cuerpo en un FermionOperator con el constructor de integrales electrónicas from_1body_tril_spin_sym(), que espera que estén empaquetadas en orden triangular inferior. La parte de dos cuerpos es un único término número-número, así que, en lugar de pasarla por la maquinaria de integrales de dos cuerpos, constrúyela directamente con cre()/ann(): el operador de Coulomb local $U \, n_{p\uparrow} n_{p\downarrow}$ para el modo de impureza $p$ (su modo $\uparrow$ en el índice p, su modo $\downarrow$ en el índice p + norb en la disposición de espín por bloques). Ambos constructores duplican el espín de los norb orbitales espaciales en 2 * norb modos fermiónicos.

Mantén la parte de un cuerpo $H_1$ y la parte de dos cuerpos $H_2$ como operadores separados en lugar de sumarlas de inmediato, de acuerdo con los términos $H_1$ y $H_2$ de las ecuaciones de la sección 1. Construir cada pieza una sola vez (en ambas bases) permite reutilizarlas más abajo. Sumadas, dan el hamiltoniano completo para la diagonalización exacta en cualquiera de las bases. Por separado, $H_2$ (base de momento) es el operador que el paso de Trotter evoluciona en el paso 4.

>>> from qiskit_fermions.operators import FermionOperator, ann, cre
>>>
>>> # U n_imp_up n_imp_down: the impurity's up/down modes sit at `imp` and `imp + norb` in the
>>> # block-spin layout (the impurity is orbital 0 in the position basis, orbital 3 in the
>>> # momentum basis -- the single relocated index found above)
>>> imp_pos, imp_mom = 0, 3
>>>
>>> one_body_position = FermionOperator.from_1body_tril_spin_sym(
...     h1e[np.tril_indices(norb)], norb
... )
>>> two_body_position = FermionOperator.from_dict(
...     {(cre(imp_pos), ann(imp_pos), cre(imp_pos + norb), ann(imp_pos + norb)): onsite}
... )
>>> one_body_momentum = FermionOperator.from_1body_tril_spin_sym(
...     h1e_momentum[np.tril_indices(norb)], norb
... )
>>> two_body_momentum = FermionOperator.from_dict(
...     {(cre(imp_mom), ann(imp_mom), cre(imp_mom + norb), ann(imp_mom + norb)): onsite}
... )
>>>
>>> hamiltonian_position = (
...     (one_body_position + two_body_position).normal_ordered().simplify(atol=1e-14)
... )
>>> hamiltonian_momentum = (
...     (one_body_momentum + two_body_momentum).normal_ordered().simplify(atol=1e-14)
... )

El cambio de base es una transformación de semejanza unitaria, así que deja intacto el espectro. Confirma que el operador es correcto comparando su valor propio más bajo en el sector $(\text{norb}, \text{nelec})$ con una diagonalización exacta. El FermionOperator expone un LinearOperator de SciPy (respaldado por un núcleo nativo de producto matriz-vector de interacción de configuraciones completa (FCI)), por lo que puede pasarse directamente a scipy.sparse.linalg.eigsh().

>>> import ffsim
>>> import scipy.sparse.linalg
>>>
>>> linop = ffsim.linear_operator(hamiltonian_momentum, norb, nelec)
>>> reference_energy, ground_state = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
>>> reference_energy = reference_energy[0]
>>> ground_state = ground_state[:, 0]
>>> print(f"exact ground-state energy: {reference_energy:.6f}")
exact ground-state energy: -13.422492

La ventaja del cambio de base es la dispersión. Para verlo de forma concreta, cuenta cuántos determinantes de la base computacional se necesitan para capturar el 99 % del peso del estado fundamental y compara la base de momento con la base de posición. La energía es independiente de la base, pero el número de determinantes no:

>>> def determinants_for(vec, fraction=0.99):
...     weights = np.sort(np.abs(vec) ** 2)[::-1]
...     return int(np.searchsorted(np.cumsum(weights), fraction) + 1)
...
>>> # the same Hamiltonian in the position basis has an identical spectrum, but a dense ground state
>>> position_linop = ffsim.linear_operator(hamiltonian_position, norb, nelec)
>>> position_energy, position_ground = scipy.sparse.linalg.eigsh(position_linop, k=1, which="SA")
>>>
>>> print(f"position-basis energy: {position_energy[0]:.6f}")
position-basis energy: -13.422492
>>> print(f"determinants for 99% weight: {determinants_for(position_ground[:, 0])} (position)")
determinants for 99% weight: 2563 (position)
>>> print(f"determinants for 99% weight: {determinants_for(ground_state)} (momentum)")
determinants for 99% weight: 21 (momentum)

Los dos operadores comparten la misma energía, pero el estado fundamental en la base de posición se reparte entre miles de determinantes, mientras que la base de momento lo concentra en unos 20, y por eso el muestreo de SKQD funciona en la base de momento y no en la de posición.

3. Prepara el estado de referencia

La referencia de SKQD $\lvert \psi_0 \rangle$ es la superposición de todas las excitaciones de los electrones más cercanos al nivel de Fermi hacia los modos de momento vacíos más próximos. La publicación de SKQD la construye con una red escrita a mano de XXPlusYYGate. Parece fuertemente correlacionada (tiene cientos de amplitudes distintas de cero en la base computacional), pero una red de rotaciones de dos modos que conservan el número es una operación de partícula única (gaussiana fermiónica), por lo que el estado que produce es un único determinante de Slater. Eso significa que puedes prepararlo con una sola compuerta PrepareSlaterDeterminant, dada la rotación de partícula única que implementa la red.

Como cada XXPlusYYGate actúa como una rotación de Givens sobre un par de modos, esa rotación puede construirse directamente en la imagen de partícula única (una matriz de Givens por compuerta) en lugar de emitir a mano las compuertas de dos qubits. Las compuertas forman un patrón de ladrillos de rotaciones entre vecinos más cercanos que se abre desde el nivel de Fermi en capas cada vez más amplias, extendiendo cada electrón cercano al nivel de Fermi por los modos vacíos situados justo por encima:

>>> def fermi_level_rotation(norb, nocc, n_excited=3):
...     """Orbital rotation exciting the near-Fermi electrons into the nearest empty modes."""
...     cos, sin = np.cos(np.pi / 4), np.sin(np.pi / 4)
...     fermi = nocc - 1  # highest occupied mode; the Fermi level sits between it and mode `nocc`
...     # Nearest-neighbor Givens rotations fanning out from the Fermi level in expanding layers:
...     # layer `depth` couples the mode `depth` steps below the level to the one above it, mixing
...     # the `n_excited` highest-occupied modes with the empty modes just above them. A final
...     # rotation carries amplitude up into the topmost empty mode of the window.
...     network = [
...         (lower, lower + 1)
...         for depth in range(n_excited)
...         for lower in range(fermi - depth, fermi + depth + 1, 2)
...     ]
...     network.append((fermi + n_excited, fermi + n_excited + 1))
...     unitary = np.eye(norb, dtype=complex)
...     for i, j in network:
...         givens = np.eye(norb, dtype=complex)
...         givens[i, i] = givens[j, j] = cos
...         givens[i, j], givens[j, i] = sin, -sin
...         unitary = givens @ unitary
...     return unitary
...
>>> reference_rotation = fermi_level_rotation(norb, nelec[0])

Nota

La rotación anterior se expresa en la base de momento, a juego con el hamiltoniano. Los electrones cercanos al nivel de Fermi son los modos de momento ocupados de mayor energía, y la red los excita hacia los modos vacíos de menor energía. Por eso la referencia debe construirse después del cambio de base del paso anterior, no antes.

La preparación del determinante de Slater se hace por sector de espín, ya que cada sector tiene su propia forma (electrones, orbitales). Una sola compuerta PrepareSlaterDeterminant (una ocupación, con los nocc modos más bajos llenos, junto con la rotación de partícula única) prepara un sector; la misma compuerta se aplica una vez sobre los modos alfa 0..norb y una vez sobre los modos beta norb..2*norb, ya que aquí ambos sectores de espín comparten la misma ocupación de referencia y la misma rotación:

>>> from qiskit_fermions.circuit.library import PrepareSlaterDeterminant
>>>
>>> occupation = [i < nelec[0] for i in range(norb)]
>>> reference_preparation = PrepareSlaterDeterminant(occupation, reference_rotation)

Esta compuerta declarativa sustituye a la red de XXPlusYYGate escrita a mano; el transpilador genera la síntesis de rotaciones de Givens.

4. Ensambla los circuitos de Krylov

Cada circuito de Krylov es la preparación del estado de referencia seguida de la evolución bajo $e^{-i k \, \Delta t \, H}$. Divide el hamiltoniano en sus partes de uno y de dos cuerpos, $H = H_1 + H_2$, y aplica Trotter a la evolución sobre esa división. Esta es la elección decisiva para la profundidad del circuito. La evolución de un cuerpo $e^{-i t H_1}$ es una operación de fermiones libres (gaussiana fermiónica), por lo que es una rotación de base de partícula única (una OrbitalRotation con unitaria $e^{-i t \, \mathbf{h}_1}$) y no algo que haya que trotterizar con Pauli término a término. La evolución de dos cuerpos $e^{-i t H_2}$ actúa solo sobre el único término de interacción local. Construir el paso de este modo permite que el transpilador emita un patrón de ladrillos poco profundo de rotaciones de Givens más una única compuerta de interacción de dos qubits, la estructura poco profunda que la publicación de SKQD construye a mano.

Ya tenemos la parte de dos cuerpos. Es el operador two_body_momentum construido en el paso 2 (el único término local). La evolución de un cuerpo $e^{-i t H_1}$ es independiente del espín, por lo que es una única OrbitalRotation por sector con unitaria $e^{-i t \, \mathbf{h}_1}$ aplicada una vez a cada sector de espín, la misma estructura paralela que la preparación del estado de referencia.

Cada circuito de Krylov evoluciona durante un tiempo total $k \, \Delta t$. Esto se realiza con un producto de Trotter de segundo orden de $k$ pasos de tamaño $\Delta t$, de modo que el error por paso se mantiene fijo a medida que crece la dimensión de Krylov. Cada paso intercala una rotación completa de un cuerpo (una por sector de espín) entre dos medios pasos de la evolución de dos cuerpos (barata, de un solo término). La dimensión de Krylov $D$ es el número de estos circuitos, con dim que va de 0 (solo el estado de referencia) a $D - 1$.

>>> import scipy.linalg
>>>
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution, OrbitalRotation
>>>
>>> num_modes = 2 * norb
>>> time_step = 0.2
>>>
>>> half_exp_h2 = Evolution(num_modes, two_body_momentum, time_step / 2)
>>> full_exp_h1 = OrbitalRotation(scipy.linalg.expm(-1j * time_step * h1e_momentum))
>>>
>>> def krylov_circuit(dim):
...     # closes over the fermionic gates defined in the outer scope:
...     # `reference_preparation`, `half_exp_h2` and `full_exp_h1`
...     circuit = FermionicCircuit(num_modes)
...     circuit.append(reference_preparation, range(norb))
...     circuit.append(reference_preparation, range(norb, num_modes))
...     for _ in range(dim):  # second-order Trotter product of `dim` steps
...         circuit.append(half_exp_h2, circuit.modes)
...         circuit.append(full_exp_h1, range(norb))
...         circuit.append(full_exp_h1, range(norb, num_modes))
...         circuit.append(half_exp_h2, circuit.modes)
...     return circuit

Nota

Los circuitos construidos aquí todavía no incluyen mediciones, porque los mismos circuitos fermiónicos también se usan en la simulación exacta del paso 6, y una medición no es una operación unitaria. El paso 5 las añade con FermionicCircuit.measure_all(), que lee el modo i en el bit clásico i.

Es preferible medir el circuito fermiónico en lugar del transpilado. Una medición fermiónica nombra el modo que lee, de modo que la disposición y el enrutamiento la llevan junto con el qubit en el que terminó ese modo, y la cadena de bits muestreada queda indexada por modo. Si en cambio mides el circuito transpilado, te toca deshacer final_index_layout() por tu cuenta, y se devuelven resultados permutados sin avisar si lo olvidas. Esto se apoya en la codificación de Jordan-Wigner en la base de ocupación que usa el preset, donde cada modo se representa con un qubit.

5. Transpila a circuitos de qubits

Transpilar un FermionicCircuit mapea sus compuertas fermiónicas a compuertas de qubits. Para ello se usa el mapeo de Jordan-Wigner jordan_wigner(). Las compuertas OrbitalRotation se sintetizan en un patrón de ladrillos de rotaciones de Givens, y la Evolution de un solo término se reduce a una única compuerta de interacción de dos qubits.

Como cada compuerta fermiónica aquí tiene una síntesis predeterminada, puedes usar directamente el generate_preset_jw_pass_manager() ya preparado, en lugar de ensamblar las etapas a mano. Cualquier argumento de palabra clave se reenvía a la etapa de qubits; pasa optimization_level=0 para obtener una imagen fiel y sin optimizar de la profundidad sintetizada. Generar toda la familia de circuitos de Krylov es entonces un bucle sobre la dimensión de Krylov. Conserva los FermionicCircuit sin transpilar junto con los circuitos de qubits transpilados. Los fermiónicos se usan en la simulación exacta (de vector de estado) del paso 6, por lo que deben permanecer sin mediciones; measure_all(inplace=False) devuelve una copia con mediciones para transpilar y deja intacto el original, y esa copia es la que ejecutaría un backend.

>>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
>>>
>>> pass_manager = generate_preset_jw_pass_manager(optimization_level=0)
>>>
>>> def krylov_circuits(krylov_dim):
...     fermionic_circuits, circuits = [], []
...     for dim in range(krylov_dim):
...         fermionic_circuit = krylov_circuit(dim)
...         fermionic_circuits.append(fermionic_circuit)
...         measured = fermionic_circuit.measure_all(inplace=False)
...         circuits.append(pass_manager.run(measured))
...     return fermionic_circuits, circuits
...
>>> krylov_dim = 5
>>> fermionic_circuits, circuits = krylov_circuits(krylov_dim)

La dimensión de Krylov se refleja directamente en la profundidad de dos qubits de los circuitos. La preparación del estado de referencia es un costo fijo, y cada potencia de Krylov adicional añade otro paso de Trotter de segundo orden del operador de evolución temporal.

>>> two_qubit = lambda instr: instr.operation.num_qubits == 2
>>> [circuit.depth(two_qubit) for circuit in circuits]
[4, 12, 21, 30, 39]

Nota

Dividir el hamiltoniano mantiene estas profundidades moderadas. La evolución de un cuerpo es una auténtica rotación orbital, por lo que se sintetiza en un patrón de ladrillos de XXPlusYYGate sobre qubits adyacentes. Cada medio paso de la interacción local es una única RZZGate que acopla los dos modos de espín de la impureza; el término barato que el producto de segundo orden intercala alrededor de la rotación. Entregar todo el hamiltoniano a una sola compuerta Evolution haría, en cambio, que la parte de un cuerpo se trotterizara con Pauli término a término, propagando largas cadenas Z de Jordan-Wigner a través de los acoplamientos impureza-baño (de largo alcance) e inflando la profundidad en un orden de magnitud.

El primer circuito prepara solo el determinante de referencia, mientras que el último contiene la mayor cantidad de pasos de Trotter (cuatro en este caso):

>>> circuits[0].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

>>> circuits[-1].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

6. Muestrea cadenas de bits de los circuitos de Krylov

En hardware, ejecutarías los circuitos transpilados y medirías. Esta es la versión sin ruido (de vector de estado) de ese paso, que muestra las cadenas de bits que recopilaría SKQD. Simula cada FermionicCircuit sin transpilar con ffsim.apply_unitary() y muestrea el vector de estado resultante con ffsim.sample_state_vector().

Una compuerta PrepareSlaterDeterminant funciona en simulación validando primero y rotando después (comprueba que el estado de entrada ocupa su determinante de referencia y luego aplica la rotación). Por tanto, la entrada es el determinante de ocupación simple: los nocc modos más bajos llenos por espín, sin rotación. La compuerta aplica reference_rotation por sí misma.

>>> from collections import Counter
>>>
>>> reference = ffsim.slater_determinant(
...     norb, (list(range(nelec[0])), list(range(nelec[1])))
... )
>>>
>>> shots = 1000
>>> counts = Counter()
>>> for dim, fermionic_circuit in enumerate(fermionic_circuits):
...     statevector = ffsim.apply_unitary(reference, fermionic_circuit, norb=norb, nelec=nelec)
...     samples = ffsim.sample_state_vector(
...         statevector, norb=norb, nelec=nelec, shots=shots, seed=dim
...     )
...     counts += Counter(samples)

Al reunir las muestras de los cinco circuitos de Krylov, el soporte es diminuto: unos pocos cientos de cadenas de bits distintas de entre los $\binom{8}{4}^2 = 4900$ determinantes del sector (4, 4). Esta concentración, consecuencia directa de la dispersión en la base de momento establecida en el paso 2, hace tratable la posterior diagonalización clásica.

Graficar los conteos con qiskit.visualization.plot_histogram() hace visible la dispersión. Unas pocas configuraciones dominan, encabezadas por el determinante de referencia.

>>> from qiskit.visualization import plot_histogram
>>>
>>> plot_histogram(dict(counts), number_to_keep=25)
<Figure size ... with 1 Axes>

7. Diagonalizar en el subespacio muestreado

Las cadenas de bits muestreadas son la entrada de la mitad clásica de SKQD. El hamiltoniano se proyecta sobre el subespacio que generan las configuraciones muestreadas y se diagonaliza allí. Esto se pasa al paquete qiskit-addon-sqd, cuya función diagonalize_fermionic_hamiltonian() ejecuta el bucle completo de diagonalización cuántica basada en muestras. Construye un subespacio a partir de lotes de las configuraciones muestreadas, diagonaliza el hamiltoniano en él y refina el subespacio de forma iterativa mediante la recuperación de configuraciones: se invierten ocupaciones para reparar las configuraciones que violan la simetría conocida del número de partículas. La recuperación de configuraciones está pensada para deshacer el daño del ruido del hardware, que corrompe las cadenas de bits muestreadas hacia un sector de número de partículas incorrecto; aquí el muestreo no tiene ruido, así que no hay nada que reparar en ese aspecto.

El solucionador consume los conteos muestreados como un BitArray, construido directamente a partir de los counts acumulados del paso 6. Pasa num_bits de forma explícita: sin él, el ancho se infiere de la cadena de bits muestreada más ancha, que es más estrecha que 2 * norb siempre que ninguna muestra ocupe el último modo, y entonces el solucionador rechaza el arreglo. El solucionador diagonaliza en la misma base de momentos que usa el resto de la construcción, por lo que toma directamente las integrales h1e_momentum y h2e_momentum del paso 2, sin reempaquetar. La energía recuperada coincide con la referencia exacta con una diferencia muy inferior a un mili-Hartree, a partir de solo unos cientos de configuraciones muestreadas. Como la energía de SQD es una cota superior variacional del estado fundamental real, siempre se sitúa en reference_energy o por encima (todavía disponible desde el paso 2):

>>> from qiskit.primitives import BitArray
>>> from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
>>>
>>> bit_array = BitArray.from_counts(dict(counts), num_bits=2 * norb)
>>>
>>> # record each iteration's batch results so we can watch the energy converge
>>> result_history = []
>>> result = diagonalize_fermionic_hamiltonian(
...     h1e_momentum,
...     h2e_momentum,
...     bit_array,
...     samples_per_batch=100,
...     norb=norb,
...     nelec=nelec,
...     num_batches=3,
...     max_iterations=3,
...     symmetrize_spin=True,
...     callback=result_history.append,
...     seed=np.random.default_rng(24),
... )
>>> print(f"SQD estimate: {result.energy:.4f} (exact: {reference_energy:.4f})")
SQD estimate: -13.4224 (exact: -13.4225)

Visualiza la ejecución: la convergencia de la energía a lo largo de las iteraciones y la ocupación promedio de cada orbital espacial en el estado fundamental recuperado. La energía por iteración es la más baja entre los lotes de esa iteración, y la ocupación suma las orbital_occupancies de los dos sectores de espín. Incluso sin ruido que corregir, la energía disminuye de una iteración a otra. Cada iteración diagonaliza dentro de num_batches lotes de solo samples_per_batch configuraciones submuestreadas del conjunto acumulado, y el refinamiento iterativo debido al arrastre de cadenas de bits mejora de forma constante qué configuraciones caen en esos lotes, de modo que la energía recuperada mejora:

>>> import matplotlib.pyplot as plt
>>> from matplotlib.ticker import MaxNLocator
>>>
>>> min_energies = [min(batch, key=lambda r: r.energy).energy for batch in result_history]
>>> occupancies = np.sum(result.orbital_occupancies, axis=0)
>>>
>>> fig, axs = plt.subplots(1, 2, figsize=(12, 5))
>>> _ = axs[0].plot(range(len(min_energies)), min_energies, marker="o", label="SQD energy")
>>> _ = axs[0].axhline(reference_energy, color="gray", linestyle="--", label="exact energy")
>>> _ = axs[0].set_xlabel("iteration")
>>> _ = axs[0].set_ylabel("energy")
>>> _ = axs[0].set_title("Energy vs SQD iteration")
>>> _ = axs[0].xaxis.set_major_locator(MaxNLocator(integer=True))
>>> _ = axs[0].legend()
>>> _ = axs[1].bar(range(norb), occupancies)
>>> _ = axs[1].set_xlabel("spatial orbital")
>>> _ = axs[1].set_ylabel("average occupancy")
>>> _ = axs[1].set_title("Occupancy per spatial orbital")
>>> fig
<Figure size ... with 2 Axes>

Esto cierra el bucle de SKQD. Los circuitos en base de momentos de los pasos 1–5 producen una muestra dispersa (paso 6), y la diagonalización basada en muestras anterior convierte esa muestra de nuevo en una energía del estado fundamental.

Próximos pasos

Esta guía ejecutó todo el flujo de SKQD en un simulador de vector de estado sin ruido. En hardware, los circuits transpilados se ejecutarían y medirían en lugar del muestreo del vector de estado del paso 6, y los conteos resultantes alimentarían la diagonalización del paso 7 sin cambios. Consulta la documentación de Qiskit para obtener ayuda con la ejecución de circuitos, y los tutoriales del addon SQD para información sobre el posprocesamiento de diagonalización de subespacios, incluido su comportamiento con muestras ruidosas, donde la recuperación de configuraciones hace la mayor parte del trabajo.

Lee la guía de la relación con ffsim para entender por qué ffsim.apply_unitary() y ffsim.sample_state_vector() funcionan de forma nativa con los circuitos fermiónicos de este paquete, y por qué el sector de número de partículas fijo usado en el paso 6 hace que la dispersión en base de momentos establecida en el paso 2 sea eficiente en muestras.