Îmbunătățirea estimării energiei unui model de rețea fermionică cu SQD

În acest tutorial implementăm un tipar Qiskit care arată cum să post-procesezi eșantioane cuantice zgomotoase pentru a găsi o aproximare a stării fundamentale a unui Hamiltonian de rețea fermionică cunoscut sub numele de modelul Anderson cu o singură impuritate. Vom urma o abordare de diagonalizare cuantică bazată pe eșantionare pentru a procesa eșantioanele preluate dintr-un set de stări din baza Krylov cu 16 Qubit-uri pe intervale de timp din ce în ce mai mari. Aceste stări sunt realizate pe dispozitivul cuantic folosind Trotterizarea evoluției temporale. Pentru a ține cont de efectul zgomotului cuantic, se folosește tehnica de recuperare a configurației. Presupunând o stare inițială bună și sparsitatea stării fundamentale, această abordare s-a dovedit a converge eficient.

Tiparul poate fi descris în patru pași:

1. Pasul 1: Maparea la problema cuantică
   • Generează un set de stări din baza Krylov (adică Circuit-uri de evoluție temporală Trotterizată) pe intervale de timp din ce în ce mai mari pentru estimarea stării fundamentale
2. Pasul 2: Optimizarea problemei
   • Transpilează Circuit-urile pentru Backend
3. Pasul 3: Execuția experimentelor
   • Extrage eșantioane din Circuit-uri folosind primitivul Sampler
4. Pasul 4: Post-procesarea rezultatelor
   • Buclă de recuperare a configurației auto-consistentă
     • Post-procesează întregul set de eșantioane de biți, folosind cunoașterea prealabilă a numărului de particule și ocupanța orbitală medie calculată la cea mai recentă iterație
     • Creează în mod probabilistic loturi de sub-eșantioane din biții recuperați
     • Proiectează și diagonalizează Hamiltonianul de rețea fermionică în fiecare subspațiu eșantionat
     • Salvează energia minimă a stării fundamentale găsită în toate loturile și actualizează ocupanța orbitală medie

Pasul 1: Maparea problemei la un Circuit cuantic

Mai întâi, vom crea Hamiltonienele cu unul și două corpuri care descriu modelul Anderson unidimensional cu o singură impuritate (SIAM) cu 7 situri de baie (8 electroni în 8 orbitale). Acest model este folosit pentru a descrie impuritățile magnetice înglobate în metale.

Apoi vom crea Circuit-urile Trotter cu 16 Qubit-uri folosite pentru a genera subspațiul Krylov cuantic.

# 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

Apoi, vom genera subspațiul Krylov cuantic cu un set de Circuit-uri cuantice Trotterizate. Aici creăm funcții ajutătoare pentru generarea stării inițiale (de referință), precum și a evoluției temporale pentru componentele cu unul și două corpuri ale Hamiltonianului. Pentru o descriere mai detaliată a acestui model și a modului în care sunt concepute Circuit-urile, consultă lucrarea.

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

Generează d stări evoluate temporal care specifică subspațiul Krylov cuantic. Aici am ales d=8. Eroarea din eșantionarea stărilor din baza Krylov converge odată cu creșterea lui d. Rețineți că particularitățile acestei instanțe de problemă permit o compilare eficientă a evoluției cu un corp folosind OrbitalRotationJW; totuși, în general, se poate folosi PauliEvolutionGate din Qiskit pentru a implementa evoluția temporală Trotterizată a Hamiltonianului 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)

Pasul 2: Optimizarea problemei

După ce am creat Circuit-urile Trotterizate, le vom optimiza pentru un hardware țintă. Trebuie să alegem dispozitivul hardware de utilizat înainte de optimizare. Vom folosi un Backend fals de 127 de Qubit-uri din qiskit_ibm_runtime pentru a emula un dispozitiv real. Pentru a rula pe un QPU real, înlocuiește pur și simplu Backend-ul fals cu unul real. Consultă documentația Qiskit IBM Runtime pentru mai multe informații.

from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke

backend = FakeSherbrooke()

Apoi, vom transpila Circuit-urile pe Backend-ul țintă folosind 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)

Pasul 3: Execuția experimentelor

După optimizarea Circuit-urilor pentru execuția pe hardware, suntem gata să le rulăm pe hardware-ul țintă și să colectăm eșantioane pentru estimarea energiei stării fundamentale. Aici folosim SamplerV2 din qiskit-ibm-runtime pentru a simula eșantioane zgomotoase preluate din Backend-ul ibm_sherbrooke. Combinăm apoi contorizările din fiecare stare a bazei Krylov într-un singur dicționar de contorizări și reprezentăm grafic cele mai frecvente 20 de biți eșantionați.

Notă: Qiskit Aer este necesar pentru simularea eșantioanelor din Circuit-uri transpilate.

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)

Pasul 4: Post-procesarea rezultatelor

Acum rulăm algoritmul SQD folosind funcția diagonalize_fermionic_hamiltonian. Consultă documentația API pentru explicații despre argumentele acestei funcții.

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

Acum reprezentăm grafic rezultatele. Primul grafic arată că după câteva iterații obținem energia exactă a stării fundamentale.

Acest exemplu este suficient de mic încât putem explora întreg spațiul Hilbert, după cum se vede în instrucțiunile de afișare de mai sus. Reamintim că spațiul Hilbert complet are dimensiunea (num_orbitals choose nelec_a) * (num_orbitals choose nelec_b). Astfel pentru această problemă: (8 choose 4)**2 = 4900. Dimensiunile subspațiului cresc datorită recuperării îmbunătățite a configurației, dar și faptului că algoritmul SQD transferă configurațiile importante de la o iterație la alta. La ultima iterație, diagonalizăm în întreg spațiul Hilbert.

Al doilea grafic arată ocupanța medie a fiecărui orbital spațial pentru soluțiile tuturor loturilor. Observăm că cu probabilitate mare fiecare orbital conține un electron.

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