Migliorare la stima dell'energia di un modello di reticolo Fermionico con SQD

In questo tutorial implementiamo un pattern Qiskit che mostra come post-elaborare campioni quantistici rumorosi per trovare un'approssimazione allo stato fondamentale di un Hamiltoniano di reticolo Fermionico noto come modello di Anderson a singola impurità. Seguiremo un approccio di diagonalizzazione quantistica basata su campioni per elaborare campioni prelevati da un insieme di 16 stati base di Krylov su qubit con intervalli di tempo crescenti. Questi stati vengono realizzati sul dispositivo quantistico tramite la Trotterizzazione dell'evoluzione temporale. Per tenere conto dell'effetto del rumore quantistico, viene utilizzata la tecnica di recupero della configurazione. Assumendo un buon stato iniziale e la sparsità dello stato fondamentale, questo approccio è dimostrato convergere in modo efficiente.

Il pattern può essere descritto in quattro passi:

1. Passo 1: Mappa al problema quantistico

   • Genera un insieme di stati base di Krylov (ovvero circuiti di evoluzione temporale Trotterizzata) su intervalli di tempo crescenti per stimare lo stato fondamentale

2. Passo 2: Ottimizza il problema

   • Transpila i circuiti per il backend

3. Passo 3: Esegui gli esperimenti

   • Estrai campioni dai circuiti usando la primitiva Sampler

4. Passo 4: Post-elabora i risultati

   • Loop di recupero della configurazione auto-consistente
     • Post-elabora l'intero insieme di campioni di bitstring, usando la conoscenza a priori del numero di particelle e l'occupazione orbitale media calcolata nell'iterazione più recente
     • Crea in modo probabilistico batch di sottocampioni dalle bitstring recuperate
     • Proietta e diagonalizza l'Hamiltoniano di reticolo Fermionico su ogni sottospazio campionato
     • Salva l'energia dello stato fondamentale minima trovata in tutti i batch e aggiorna l'occupazione orbitale media

Passo 1: Mappa il problema a un circuito quantistico

Per prima cosa creeremo gli Hamiltoniani a uno e due corpi che descrivono il modello di Anderson a singola impurità (SIAM) unidimensionale con 7 siti di bagno (8 elettroni in 8 orbitali). Questo modello viene utilizzato per descrivere impurità magnetiche incorporate nei metalli.

Poi creeremo i circuiti di Trotter a 16 qubit usati per generare il sottospazio di Krylov quantistico.

# 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

Successivamente genereremo il sottospazio di Krylov quantistico con un insieme di circuiti quantistici Trotterizzati. Qui creiamo degli helper per generare lo stato iniziale (di riferimento) e l'evoluzione temporale per le parti a uno e due corpi dell'Hamiltoniano. Per una descrizione più dettagliata di questo modello e di come i circuiti sono progettati, fare riferimento a l'articolo.

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 stati con evoluzione temporale che specificano il sottospazio di Krylov quantistico. Qui abbiamo scelto d=8. L'errore derivante dal campionamento degli stati base di Krylov converge all'aumentare di d. Si noti che le particolarità di questa istanza del problema consentono una compilazione efficiente dell'evoluzione a un corpo usando OrbitalRotationJW; tuttavia, in generale, si potrebbe usare il PauliEvolutionGate di Qiskit per implementare l'evoluzione temporale Trotterizzata dell'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)

Passo 2: Ottimizza il problema

Dopo aver creato i circuiti Trotterizzati, li ottimizzeremo per un hardware di destinazione. Dobbiamo scegliere il dispositivo hardware da usare prima dell'ottimizzazione. Useremo un backend fittizio a 127 qubit da qiskit_ibm_runtime per emulare un dispositivo reale. Per eseguire su una QPU reale, basta sostituire il backend fittizio con un backend reale. Consulta la documentazione di Qiskit IBM Runtime per maggiori informazioni.

from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke

backend = FakeSherbrooke()

Successivamente transpileremo i circuiti verso il backend di destinazione 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)

Passo 3: Esegui gli esperimenti

Dopo aver ottimizzato i circuiti per l'esecuzione hardware, siamo pronti a eseguirli sull'hardware di destinazione e raccogliere campioni per la stima dell'energia dello stato fondamentale. Qui usiamo SamplerV2 da qiskit-ibm-runtime per simulare campioni rumorosi prelevati dal backend ibm_sherbrooke. Poi combiniamo i conteggi da ciascuno degli stati base di Krylov in un unico dizionario di conteggi e tracciamo le 20 bitstring campionate più comunemente.

Nota: Qiskit Aer è richiesto per simulare campioni da circuiti transpilati.

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)

Passo 4: Post-elabora i risultati

Ora eseguiamo l'algoritmo SQD usando la funzione diagonalize_fermionic_hamiltonian. Consulta la documentazione API per le spiegazioni degli argomenti di questa funzione.

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

Ora tracciamo i risultati. Il primo grafico mostra che dopo alcune iterazioni otteniamo l'energia esatta dello stato fondamentale.

Questo esempio è abbastanza piccolo da permetterci di esplorare l'intero spazio di Hilbert, come si vede nelle istruzioni di stampa sopra. Ricorda, l'intero spazio di Hilbert ha dimensione (num_orbitals choose nelec_a) * (num_orbitals choose nelec_b). Quindi per questo problema: (8 choose 4)**2 = 4900. Le dimensioni del sottospazio aumentano grazie al recupero della configurazione migliorato, e anche al fatto che l'algoritmo SQD trasporta le configurazioni importanti da un'iterazione alla successiva. Nell'ultima iterazione, stiamo diagonalizzando nell'intero spazio di Hilbert.

Il secondo grafico mostra l'occupazione media di ogni orbitale spaziale in tutte le soluzioni dei batch. Vediamo che con alta probabilità ogni orbitale contiene un elettrone.

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