Mejora de la estimación de energía de un modelo de red fermiónico con SQD

En este tutorial implementamos un patrón de Qiskit que muestra cómo posprocesar muestras cuánticas ruidosas para encontrar una aproximación al estado base de un Hamiltoniano de red fermiónico conocido como el modelo de Anderson de una sola impureza. Seguiremos un enfoque de diagonalización cuántica basado en muestras para procesar muestras tomadas de un conjunto de estados de la base de Krylov de 16 qubits en intervalos de tiempo crecientes. Estos estados se realizan en el dispositivo cuántico mediante la Troterización de la evolución temporal. Para tener en cuenta el efecto del ruido cuántico, se utiliza la técnica de recuperación de configuración. Asumiendo un buen estado inicial y la dispersión del estado base, se ha demostrado que este enfoque converge eficientemente.

El patrón se puede describir en cuatro pasos:

1. Paso 1: Mapear al problema cuántico

   • Generar un conjunto de estados de la base de Krylov (es decir, circuitos de evolución temporal troterizada) en intervalos de tiempo crecientes para estimar el estado base

2. Paso 2: Optimizar el problema

   • Transpilar los circuitos para el backend

3. Paso 3: Ejecutar experimentos

   • Extraer muestras de los circuitos usando la primitiva Sampler

4. Paso 4: Posprocesar los resultados

   • Bucle de recuperación de configuración autoconsonante
     • Posprocesar el conjunto completo de muestras de cadenas de bits, usando el conocimiento previo del número de partículas y la ocupación orbital promedio calculada en la iteración más reciente
     • Crear probabilísticamente lotes de submuestras a partir de cadenas de bits recuperadas
     • Proyectar y diagonalizar el Hamiltoniano de red fermiónico sobre cada subespacio muestreado
     • Guardar la energía mínima del estado base encontrada en todos los lotes y actualizar la ocupación orbital promedio

Paso 1: Mapear el problema a un Circuit cuántico

Primero, crearemos los Hamiltonianos de uno y dos cuerpos que describen el modelo de Anderson de una sola impureza (SIAM, por sus siglas en inglés) unidimensional con 7 sitios de baño (8 electrones en 8 orbitales). Este modelo se utiliza para describir impurezas magnéticas embebidas en metales.

Luego crearemos los circuitos de Trotter de 16 qubits utilizados para generar el subespacio de Krylov cuántico.

# Added by doQumentation — required packages for this notebook
!pip install -q ffsim matplotlib numpy qiskit qiskit-addon-sqd qiskit-ibm-runtime scipy

import numpy as np

n_bath = 7  # number of bath sites

V = 1  # hybridization amplitude
t = 1  # bath hopping amplitude
U = 10  # Impurity onsite repulsion
eps = -U / 2  # Chemical potential for the impurity

# Place the impurity on the first qubit
impurity_index = 0

# One body matrix elements in the "position" basis
h1e = -t * np.diag(np.ones(n_bath), k=1) - t * np.diag(np.ones(n_bath), k=-1)
h1e[impurity_index, impurity_index + 1] = -V
h1e[impurity_index + 1, impurity_index] = -V
h1e[impurity_index, impurity_index] = eps

# Two body matrix elements in the "position" basis
h2e = np.zeros((n_bath + 1, n_bath + 1, n_bath + 1, n_bath + 1))
h2e[impurity_index, impurity_index, impurity_index, impurity_index] = U

A continuación, generaremos el subespacio de Krylov cuántico con un conjunto de circuitos cuánticos troterizados. Aquí creamos funciones auxiliares para generar el estado inicial (de referencia) así como la evolución temporal para las partes de uno y dos cuerpos del Hamiltoniano. Para una descripción más detallada de este modelo y cómo se diseñan los circuitos, consulta el artículo.

import ffsim
import scipy
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate

n_modes = n_bath + 1
nelec = (n_modes // 2, n_modes // 2)

dt = 0.2
Utar = scipy.linalg.expm(-1j * dt * h1e)

# The reference state
def initial_state(q_circuit, norb, nocc):
    """Prepare an initial state."""
    for i in range(nocc):
        q_circuit.append(XGate(), [i])
        q_circuit.append(XGate(), [norb + i])
    rot = XXPlusYYGate(np.pi / 2, -np.pi / 2)

    for i in range(3):
        for j in range(nocc - i - 1, nocc + i, 2):
            q_circuit.append(rot, [j, j + 1])
            q_circuit.append(rot, [norb + j, norb + j + 1])
    q_circuit.append(rot, [j + 1, j + 2])
    q_circuit.append(rot, [norb + j + 1, norb + j + 2])

# The one-body time evolution
free_fermion_evolution = ffsim.qiskit.OrbitalRotationJW(n_modes, Utar)

# The two-body time evolution
def append_diagonal_evolution(dt, U, impurity_qubit, num_orb, q_circuit):
    """Append two-body time evolution to a quantum circuit."""
    if U != 0:
        q_circuit.append(
            CPhaseGate(-dt / 2 * U),
            [impurity_qubit, impurity_qubit + num_orb],
        )

Genera d estados con evolución temporal que especifican el subespacio de Krylov cuántico. Aquí hemos elegido d=8. El error al muestrear los estados de la base de Krylov converge con d creciente. Ten en cuenta que las particularidades de esta instancia del problema permiten una compilación eficiente de la evolución de un cuerpo usando OrbitalRotationJW; sin embargo, en general, se podría usar la PauliEvolutionGate de Qiskit para implementar la evolución temporal troterizada del Hamiltoniano completo.

# Generate the initial state
qubits = QuantumRegister(2 * n_modes, name="q")
init_state = QuantumCircuit(qubits)
initial_state(init_state, n_modes, n_modes // 2)
init_state.draw("mpl", scale=0.4, fold=-1)

d = 8  # Number of Krylov basis states
circuits = []
for i in range(d):
    circ = init_state.copy()
    circuits.append(circ)
    for _ in range(i):
        append_diagonal_evolution(dt, U, impurity_index, n_modes, circ)
        circ.append(free_fermion_evolution, qubits)
        append_diagonal_evolution(dt, U, impurity_index, n_modes, circ)
    circ.measure_all()

circuits[0].draw("mpl", scale=0.4, fold=-1)

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

Paso 2: Optimizar el problema

Una vez que hemos creado los circuitos troterizados, los optimizaremos para un hardware de destino. Necesitamos elegir el dispositivo de hardware a usar antes de la optimización. Usaremos un backend falso de 127 qubits de qiskit_ibm_runtime para emular un dispositivo real. Para ejecutar en una QPU real, simplemente reemplaza el backend falso por uno real. Consulta la documentación de Qiskit IBM Runtime para más información.

from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke

backend = FakeSherbrooke()

A continuación, transpilaremos los circuitos al backend de destino usando Qiskit.

from qiskit.transpiler import generate_preset_pass_manager

# The circuit needs to be transpiled to the AerSimulator target
pass_manager = generate_preset_pass_manager(optimization_level=3, backend=backend)
isa_circuits = pass_manager.run(circuits)

Paso 3: Ejecutar experimentos

Después de optimizar los circuitos para la ejecución en hardware, estamos listos para ejecutarlos en el hardware de destino y recopilar muestras para la estimación de la energía del estado base. Aquí usamos SamplerV2 de qiskit-ibm-runtime para simular muestras ruidosas tomadas del backend ibm_sherbrooke. Luego combinamos los conteos de cada uno de los estados de la base de Krylov en un único diccionario de conteos y graficamos las 20 cadenas de bits más frecuentemente muestreadas.

Nota: Se requiere Qiskit Aer para simular muestras de circuitos transpilados.

from qiskit.primitives import BitArray
from qiskit.visualization import plot_histogram
from qiskit_ibm_runtime import SamplerV2 as Sampler

# Sample from the circuits
noisy_sampler = Sampler(backend, options={"simulator": {"seed_simulator": 24}})
job = noisy_sampler.run(isa_circuits, shots=500)

# Combine the counts from the individual Trotter circuits
bit_array = BitArray.concatenate_shots([result.data.meas for result in job.result()])

plot_histogram(bit_array.get_counts(), number_to_keep=20)

Paso 4: Posprocesar los resultados

Ahora ejecutamos el algoritmo SQD usando la función diagonalize_fermionic_hamiltonian. Consulta la documentación de la API para obtener explicaciones de los argumentos de esta función.

from qiskit_addon_sqd.fermion import SCIResult, diagonalize_fermionic_hamiltonian

# List to capture intermediate results
result_history = []

def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy}")
        print(f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}")

rng = np.random.default_rng(24)
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=300,
    norb=n_modes,
    nelec=nelec,
    num_batches=3,
    max_iterations=10,
    symmetrize_spin=True,
    callback=callback,
    seed=rng,
)

Iteration 1
	Subsample 0
		Energy: -13.257128325607997
		Subspace dimension: 3969
	Subsample 1
		Energy: -13.257128325607997
		Subspace dimension: 3969
	Subsample 2
		Energy: -13.257128325607997
		Subspace dimension: 3969
Iteration 2
	Subsample 0
		Energy: -13.319666127542039
		Subspace dimension: 4096
	Subsample 1
		Energy: -13.420534292304595
		Subspace dimension: 4624
	Subsample 2
		Energy: -9.136171594591085
		Subspace dimension: 4624
Iteration 3
	Subsample 0
		Energy: -13.422491814612833
		Subspace dimension: 4900
	Subsample 1
		Energy: -13.422491814612833
		Subspace dimension: 4900
	Subsample 2
		Energy: -13.422491814612833
		Subspace dimension: 4900
Iteration 4
	Subsample 0
		Energy: -13.422491814612833
		Subspace dimension: 4900
	Subsample 1
		Energy: -13.422491814612833
		Subspace dimension: 4900
	Subsample 2
		Energy: -13.422491814612833
		Subspace dimension: 4900

Ahora graficamos los resultados. El primer gráfico muestra que después de algunas iteraciones, obtenemos la energía exacta del estado base.

Este ejemplo es lo suficientemente pequeño como para explorar el espacio de Hilbert completo, como se ve en las declaraciones de impresión anteriores. Recuerda, el espacio de Hilbert completo tiene dimensión (num_orbitals choose nelec_a) * (num_orbitals choose nelec_b). Entonces para este problema: (8 choose 4)**2 = 4900. Las dimensiones del subespacio aumentan debido a la recuperación de configuración mejorada, y también al hecho de que el algoritmo SQD lleva las configuraciones importantes de una iteración a la siguiente. En la última iteración, estamos diagonalizando en el espacio de Hilbert completo.

El segundo gráfico muestra la ocupación promedio de cada orbital espacial en todas las soluciones de los lotes. Vemos que con alta probabilidad cada orbital contiene un electrón.

import matplotlib.pyplot as plt

exact_energy = -13.422491814605827
min_es = [min(result, key=lambda res: res.energy).energy for result in result_history]
min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])

# Data for energies plot
x1 = range(len(result_history))
yt1 = list(np.arange(-13.5, -13.1, 0.1))
ytl = [f"{i:.1f}" for i in yt1]

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, min_es, label="SQD energy", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(ytl)
axs[0].axhline(y=exact_energy, color="#BF5700", linestyle="--", label="Exact energy")
axs[0].set_title("Approximated Ground State Energy vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

print(f"Exact energy: {exact_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - exact_energy):.5f}")
plt.tight_layout()
plt.show()

Exact energy: -13.42249
SQD energy: -13.42249
Absolute error: 0.00000