Meningkatkan estimasi energi model kisi Fermionik dengan SQD
Dalam tutorial ini kita mengimplementasikan sebuah pola Qiskit yang menunjukkan cara memproses sampel kuantum yang bising untuk menemukan pendekatan ground state dari Hamiltonian kisi Fermionik yang dikenal sebagai model Anderson impuritas tunggal. Kita akan mengikuti pendekatan diagonalisasi kuantum berbasis sampel untuk memproses sampel yang diambil dari sekumpulan 16-Qubit keadaan basis Krylov pada interval waktu yang semakin meningkat. Keadaan-keadaan ini direalisasikan pada perangkat kuantum menggunakan Trotterisasi evolusi waktu. Untuk memperhitungkan efek noise kuantum, teknik pemulihan konfigurasi digunakan. Dengan asumsi keadaan awal yang baik dan sparsitas dari ground state, pendekatan ini terbukti konvergen secara efisien.
Polanya dapat digambarkan dalam empat langkah:
- Langkah 1: Petakan ke masalah kuantum
- Hasilkan sekumpulan keadaan basis Krylov (yaitu, Circuit evolusi waktu Trotterisasi) pada interval waktu yang semakin meningkat untuk mengestimasi ground state
- Langkah 2: Optimalkan masalah
- Transpile Circuit untuk Backend
- Langkah 3: Jalankan eksperimen
- Ambil sampel dari Circuit menggunakan primitif
Sampler
- Ambil sampel dari Circuit menggunakan primitif
- Langkah 4: Pasca-proses hasil
- Loop pemulihan konfigurasi mandiri
- Pasca-proses seluruh kumpulan sampel bitstring, menggunakan pengetahuan awal tentang jumlah partikel dan rata-rata hunian orbital yang dihitung pada iterasi terbaru
- Buat batch subsampel secara probabilistik dari bitstring yang dipulihkan
- Proyeksikan dan diagonalisasikan Hamiltonian kisi Fermionik pada setiap subruang yang diambil sampelnya
- Simpan energi ground state minimum yang ditemukan di semua batch dan perbarui rata-rata hunian orbital
- Loop pemulihan konfigurasi mandiri
Langkah 1: Petakan masalah ke Circuit kuantum​
Pertama, kita akan membuat Hamiltonian satu-badan dan dua-badan yang menggambarkan model Anderson impuritas tunggal (SIAM) satu dimensi dengan 7 situs bath (8 elektron dalam 8 orbital). Model ini digunakan untuk menggambarkan impuritas magnetik yang tertanam dalam logam.
Kemudian kita akan membuat Circuit Trotter 16-Qubit yang digunakan untuk menghasilkan subruang Krylov kuantum.
# 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
Selanjutnya, kita akan menghasilkan subruang Krylov kuantum dengan sekumpulan Circuit kuantum Trotterisasi. Di sini kita membuat fungsi pembantu untuk menghasilkan keadaan awal (referensi) serta evolusi waktu untuk bagian satu-badan dan dua-badan dari Hamiltonian. Untuk deskripsi yang lebih rinci tentang model ini dan bagaimana Circuit dirancang, silakan merujuk ke paper.
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],
)
Hasilkan d keadaan yang berevolusi terhadap waktu yang menentukan subruang Krylov kuantum. Di sini, kita telah memilih d=8. Error dari pengambilan sampel keadaan basis Krylov konvergen seiring meningkatnya d. Perlu dicatat bahwa ciri khas instance masalah ini memungkinkan kompilasi efisien dari evolusi satu-badan menggunakan OrbitalRotationJW; namun, secara umum, seseorang dapat menggunakan PauliEvolutionGate dari Qiskit untuk mengimplementasikan evolusi waktu Trotterisasi dari Hamiltonian lengkap.
# 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)
Langkah 2: Optimalkan masalah​
Setelah kita membuat Circuit Trotterisasi, kita akan mengoptimalkannya untuk perangkat keras target. Kita perlu memilih perangkat keras yang akan digunakan sebelum optimasi. Kita akan menggunakan Backend 127-Qubit palsu dari qiskit_ibm_runtime untuk mengemulasi perangkat nyata. Untuk menjalankan pada QPU nyata, cukup ganti Backend palsu dengan Backend nyata. Lihat dokumentasi Qiskit IBM Runtime untuk info lebih lanjut.
from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke
backend = FakeSherbrooke()
Selanjutnya, kita akan men-transpile Circuit ke Backend target menggunakan 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)
Langkah 3: Jalankan eksperimen​
Setelah mengoptimalkan Circuit untuk eksekusi perangkat keras, kita siap menjalankannya pada perangkat keras target dan mengumpulkan sampel untuk estimasi energi ground state. Di sini kita menggunakan SamplerV2 dari qiskit-ibm-runtime untuk mensimulasikan sampel bising yang diambil dari Backend ibm_sherbrooke. Kemudian kita menggabungkan counts dari setiap keadaan basis Krylov menjadi satu kamus counts dan memplot 20 bitstring yang paling sering diambil sampelnya.
Catatan: Qiskit Aer diperlukan untuk mensimulasikan sampel dari Circuit yang telah di-transpile.
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)
Langkah 4: Pasca-proses hasil​
Sekarang, kita menjalankan algoritma SQD menggunakan fungsi diagonalize_fermionic_hamiltonian. Lihat dokumentasi API untuk penjelasan argumen fungsi ini.
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
Sekarang, kita memplot hasilnya. Plot pertama menunjukkan bahwa setelah beberapa iterasi, kita mendapatkan energi ground state yang tepat.
Contoh ini cukup kecil sehingga kita dapat mengeksplorasi seluruh ruang Hilbert, seperti terlihat pada pernyataan print di atas. Ingat, ruang Hilbert penuh berdimensi (num_orbitals pilih nelec_a) * (num_orbitals pilih nelec_b). Jadi untuk masalah ini: (8 choose 4)**2 = 4900. Dimensi subruang meningkat karena pemulihan konfigurasi yang ditingkatkan, dan juga karena algoritma SQD membawa konfigurasi penting dari satu iterasi ke iterasi berikutnya. Pada iterasi terakhir, kita melakukan diagonalisasi pada seluruh ruang Hilbert.
Plot kedua menunjukkan rata-rata hunian setiap orbital spasial di semua solusi batch. Kita melihat bahwa dengan probabilitas tinggi setiap orbital mengandung satu elektron.
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
