שיפור אמידת אנרגיה של מודל סריג פרמיוני עם SQD
במדריך זה נממש תבנית Qiskit שמדגימה כיצד לעבד בדיעבד דגימות קוונטיות רועשות כדי למצוא קירוב למצב היסוד של המילטוניאן סריג פרמיוני הידוע בשם מודל אנדרסון עם פגם בודד. נעקוב אחר גישת אלכסון קוונטי מבוסס-דגימות כדי לעבד דגימות שנלקחו מתוך קבוצת 16 מצבי בסיס Krylov של Qubit, על פני מרווחי זמן גדלים. מצבים אלו ממומשים על המכשיר הקוונטי באמצעות Trotterization של אבולוציית הזמן. כדי להתחשב בהשפעת הרעש הקוונטי, נעשה שימוש בטכניקת שחזור הקונפיגורציה. בהנחה של מצב התחלתי טוב ודלילות של מצב היסוד, הוכח שגישה זו מתכנסת ביעילות.
התבנית מתוארת בארבעה שלבים:
- שלב 1: מיפוי לבעיה קוונטית
- יצירת קבוצת מצבי בסיס Krylov (כלומר, Circuit של אבולוציית זמן מסוג Trotterized) על פני מרווחי זמן גדלים לצורך אמידת מצב היסוד
- שלב 2: אופטימיזציה של הבעיה
- Transpile של ה-Circuit לפי ה-Backend
- שלב 3: ביצוע ניסויים
- שליפת דגימות מה-Circuit באמצעות הפרימיטיב
Sampler
- שליפת דגימות מה-Circuit באמצעות הפרימיטיב
- שלב 4: עיבוד בדיעבד של התוצאות
- לולאת שחזור קונפיגורציה עצמית-עקבית
- עיבוד בדיעבד של מלוא קבוצת דגימות ה-bitstring, תוך שימוש ב ידע מוקדם על מספר החלקיקים ועל תפוסת המסלולים הממוצעת שחושבה בחזרה האחרונה
- יצירה הסתברותית של אצוות תת-דגימות מ-bitstring ששוחזרו
- הקרנה ואלכסון של המילטוניאן הסריג הפרמיוני על כל תת-מרחב שנדגם
- שמירת אנרגיית מצב היסוד המינימלית שנמצאה בכל האצוות ועדכון תפוסת המסלולים הממוצעת
- לולאת שחזור קונפיגורציה עצמית-עקבית
שלב 1: מיפוי הבעיה ל-Circuit קוונטי
ראשית, ניצור את המילטוניאנים חד-גופיים ודו-גופיים המתארים את מודל אנדרסון חד-ממדי עם פגם בודד (SIAM) עם 7 אתרי אמבטיה (8 אלקטרונים ב-8 מסלולים). מודל זה משמש לתיאור פגמים מגנטיים משובצים במתכות.
לאחר מכן ניצור את ה-Circuit של Trotter עם 16 Qubit, שישמשו ליצירת תת-מרחב Krylov הקוונטי.
# 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
בשלב הבא, ניצור את תת-המרחב הקוונטי של Krylov עם קבוצת Circuit קוונטיים מסוג Trotterized. כאן ניצור פונקציות עזר ליצירת המצב ההתחלתי (הייחוס), וכן לאבולוציית הזמן עבור החלקים חד-גופיים ודו-גופיים של המילטוניאן. לתיאור מפורט יותר של מודל זה ושל אופן עיצוב ה-Circuit, ניתן לפנות למאמר.
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],
)
ניצור d מצבים מאופיינים בזמן, המגדירים את תת-המרחב הקוונטי של Krylov. כאן בחרנו d=8. השגיאה הנובעת מדגימת מצבי בסיס Krylov מתכנסת עם עליית d. שים לב שהפרטים של מופע בעיה זה מאפשרים קומפילציה יעילה של האבולוציה החד-גופית באמצעות OrbitalRotationJW; עם זאת, באופן כללי, ניתן להשתמש ב-PauliEvolutionGate של Qiskit כדי לממש את אבולוציית הזמן מסוג Trotterized של המילטוניאן המלא.
# 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)
שלב 2: אופטימיזציה של הבעיה
לאחר שיצרנו את ה-Circuit מסוג Trotterized, נבצע אופטימיזציה שלהם עבור חומרה מסוימת. צריך לבחור את מכשיר החומרה לפני האופטימיזציה. נשתמש ב-Backend מדומה בן 127 Qubit מ-qiskit_ibm_runtime כדי לדמות מכשיר אמיתי. להרצה על QPU אמיתי, פשוט החלף את ה-Backend המדומה ב-Backend אמיתי. ניתן למצוא מידע נוסף בתיעוד Qiskit IBM Runtime.
from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke
backend = FakeSherbrooke()
בשלב הבא, נבצע Transpile של ה-Circuit אל ה-Backend היעד באמצעות 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)
שלב 3: ביצוע ניסויים
לאחר אופטימיזציה של ה-Circuit לביצוע חומרתי, אנחנו מוכנים להריצם על החומרה היעד ולאסוף דגימות לאמידת אנרגיית מצב היסוד. כאן אנחנו משתמשים ב-SamplerV2 מ-qiskit-ibm-runtime כדי לדמות דגימות רועשות שנלקחות מה-Backend של ibm_sherbrooke. לאחר מכן אנחנו משלבים את הספירות מכל אחד ממצבי בסיס Krylov למילון ספירות אחד ומשרטטים את 20 ה-bitstring הנפוצים ביותר.
הערה: Qiskit Aer נדרש לדמות דגימות מ-Circuit מתורגמים.
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)
שלב 4: עיבוד בדיעבד של התוצאות
כעת נריץ את אלגוריתם SQD באמצעות הפונקציה diagonalize_fermionic_hamiltonian. ראה את תיעוד ה-API להסבר על הארגומנטים של פונקציה זו.
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
כעת נשרטט את התוצאות. הגרף הראשון מראה שלאחר מספר חזרות, אנחנו מקבלים את אנרגיית מצב היסוד המדויקת.
דוגמה זו קטנה מספיק כדי שנוכל לחקור את מרחב הילברט המלא, כפי שנראה בהצהרות ההדפסה לעיל. זכור, מרחב הילברט המלא הוא בממד (num_orbitals choose nelec_a) * (num_orbitals choose nelec_b). כך עבור בעיה זו: (8 choose 4)**2 = 4900. ממדי תת-המ רחב גדלים בשל שחזור קונפיגורציה משופר, וגם בשל העובדה שאלגוריתם SQD מעביר קונפיגורציות חשובות מחזרה לחזרה. בחזרה האחרונה, אנחנו מבצעים אלכסון בכל מרחב הילברט.
הגרף השני מראה את תפוסת המסלולים המרחביים הממוצעת על פני הפתרונות של כל האצוות. אנחנו רואים שבהסתברות גבוהה כל מסלול מכיל אלקטרון אחד.
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
