Améliorer l'estimation d'énergie d'un modèle de réseau fermionique avec SQD

Dans ce tutoriel, nous implémentons un motif Qiskit montrant comment post-traiter des échantillons quantiques bruités pour trouver une approximation de l'état fondamental d'un Hamiltonien de réseau fermionique connu sous le nom de modèle d'Anderson à une impureté. Nous suivrons une approche de diagonalisation quantique basée sur des échantillons pour traiter des échantillons tirés d'un ensemble d'états de base Krylov à 16 qubits sur des intervalles de temps croissants. Ces états sont réalisés sur le dispositif quantique à l'aide de la trotterisation de l'évolution temporelle. Afin de tenir compte de l'effet du bruit quantique, la technique de récupération de configuration est utilisée. En supposant un bon état initial et la parcimonie de l'état fondamental, cette approche est prouvée converger efficacement.

Le motif peut être décrit en quatre étapes :

1. Étape 1 : Mapper le problème quantique

   • Générer un ensemble d'états de base Krylov (c'est-à-dire, des circuits d'évolution temporelle trotterisés) sur des intervalles de temps croissants pour estimer l'état fondamental

2. Étape 2 : Optimiser le problème

   • Transpiler les circuits pour le Backend

3. Étape 3 : Exécuter les expériences

   • Tirer des échantillons des circuits en utilisant la primitive Sampler

4. Étape 4 : Post-traiter les résultats

   • Boucle de récupération de configuration auto-cohérente
     • Post-traiter l'ensemble complet des échantillons de chaînes de bits, en utilisant la connaissance préalable du nombre de particules et l'occupation orbitale moyenne calculée lors de la dernière itération
     • Créer de façon probabiliste des lots de sous-échantillons à partir de chaînes de bits récupérées
     • Projeter et diagonaliser le Hamiltonien de réseau fermionique sur chaque sous-espace échantillonné
     • Sauvegarder l'énergie minimale de l'état fondamental trouvée dans tous les lots et mettre à jour l'occupation orbitale moyenne

Étape 1 : Mapper le problème vers un Circuit quantique

Tout d'abord, nous allons créer les Hamiltoniens à un et deux corps décrivant le modèle d'Anderson à une impureté unidimensionnel (SIAM) avec 7 sites de bain (8 électrons dans 8 orbitales). Ce modèle est utilisé pour décrire des impuretés magnétiques incorporées dans des métaux.

Ensuite, nous allons créer les circuits de Trotter à 16 qubits utilisés pour générer le sous-espace de Krylov quantique.

# 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

Ensuite, nous allons générer le sous-espace de Krylov quantique avec un ensemble de circuits quantiques trotterisés. Ici, nous créons des fonctions auxiliaires pour générer l'état initial (de référence) ainsi que l'évolution temporelle pour les parties à un et deux corps du Hamiltonien. Pour une description plus détaillée de ce modèle et de la conception des circuits, consulte l'article.

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],
        )

Génère d états soumis à évolution temporelle qui spécifient le sous-espace de Krylov quantique. Ici, nous avons choisi d=8. L'erreur due à l'échantillonnage des états de base Krylov converge avec l'augmentation de d. Note que les particularités de cette instance de problème permettent une compilation efficace de l'évolution à un corps en utilisant OrbitalRotationJW ; cependant, en général, on pourrait utiliser la PauliEvolutionGate de Qiskit pour implémenter l'évolution temporelle trotterisée du Hamiltonien complet.

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

Étape 2 : Optimiser le problème

Après avoir créé les circuits trotterisés, nous allons les optimiser pour un matériel cible. Nous devons choisir le dispositif matériel à utiliser avant l'optimisation. Nous utiliserons un faux Backend de 127 qubits de qiskit_ibm_runtime pour émuler un vrai dispositif. Pour exécuter sur un vrai QPU, il suffit de remplacer le faux Backend par un vrai Backend. Consulte la documentation Qiskit IBM Runtime pour plus d'informations.

from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke

backend = FakeSherbrooke()

Ensuite, nous allons transpiler les circuits vers le Backend cible en utilisant 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)

Étape 3 : Exécuter les expériences

Après avoir optimisé les circuits pour l'exécution matérielle, nous sommes prêts à les exécuter sur le matériel cible et à collecter des échantillons pour l'estimation de l'énergie de l'état fondamental. Ici, nous utilisons SamplerV2 de qiskit-ibm-runtime pour simuler des échantillons bruités tirés du Backend ibm_sherbrooke. Nous combinons ensuite les comptes de chacun des états de base Krylov en un seul dictionnaire de comptes et traçons les 20 chaînes de bits les plus fréquemment échantillonnées.

Remarque : Qiskit Aer est requis pour simuler des échantillons à partir de circuits transpilés.

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)

Étape 4 : Post-traiter les résultats

Maintenant, nous exécutons l'algorithme SQD en utilisant la fonction diagonalize_fermionic_hamiltonian. Consulte la documentation de l'API pour les explications des arguments de cette fonction.

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

Maintenant, nous traçons les résultats. Le premier graphe montre qu'après quelques itérations, nous obtenons l'énergie exacte de l'état fondamental.

Cet exemple est suffisamment petit pour que nous puissions explorer l'espace de Hilbert complet, comme le montrent les instructions d'affichage ci-dessus. Rappelle-toi que l'espace de Hilbert complet est de dimension (num_orbitals choose nelec_a) * (num_orbitals choose nelec_b). Donc pour ce problème : (8 choose 4)**2 = 4900. Les dimensions du sous-espace augmentent grâce à une récupération de configuration améliorée, ainsi qu'au fait que l'algorithme SQD transporte les configurations importantes d'une itération à l'autre. À la dernière itération, nous diagonalisons dans l'espace de Hilbert entier.

Le deuxième graphe montre l'occupation moyenne de chaque orbitale spatiale à travers toutes les solutions des lots. Nous voyons qu'avec une forte probabilité, chaque orbitale contient un électron.

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