Pagpapabuti ng energy estimation ng isang Fermionic lattice model gamit ang SQD

Sa tutorial na ito ay nag-iimplementa tayo ng isang Qiskit pattern na nagpapakita kung paano i-post-process ang noisy quantum samples upang makahanap ng aproksimasyon sa ground state ng isang Fermionic lattice Hamiltonian na kilala bilang single-impurity Anderson model. Susundan natin ang isang sample-based quantum diagonalization approach upang i-process ang samples na kinuha mula sa isang set ng 16-qubit Krylov basis states sa lumalaking time intervals. Ang mga state na ito ay isinasakatuparan sa quantum device gamit ang Trotterization ng time evolution. Upang isaalang-alang ang epekto ng quantum noise, ginagamit ang configuration recovery technique. Sa pag-aakalang isang magandang initial state at sparsity ng ground state, ang approach na ito ay napatunayang ma-converge nang efficient.

Ang pattern ay maaaring ilarawan sa apat na hakbang:

1. Hakbang 1: I-map sa quantum problem
   • Bumuo ng set ng Krylov basis states (i.e., Trotterized time-evolution circuits) sa lumalaking time intervals para sa pag-estimate ng ground state
2. Hakbang 2: I-optimize ang problem
   • I-transpile ang circuits para sa backend
3. Hakbang 3: Isagawa ang mga eksperimento
   • Kumuha ng samples mula sa circuits gamit ang Sampler primitive
4. Hakbang 4: I-post-process ang resulta
   • Self-consistent configuration recovery loop
     • I-post-process ang buong set ng bitstring samples, gamit ang prior knowledge ng particle number at ang average orbital occupancy na kinakalkula sa pinakahuling iteration
     • Probabilistically gumawa ng batches ng subsamples mula sa recovered bitstrings
     • I-project at i-diagonalize ang Fermionic lattice Hamiltonian sa bawat sampled subspace
     • I-save ang minimum ground state energy na natagpuan sa lahat ng batches at i-update ang avg orbital occupancy

Hakbang 1: I-map ang problem sa isang quantum circuit

Una, gagawa tayo ng one- at two-body Hamiltonians na naglalarawan sa one-dimensional single-impurity Anderson model (SIAM) na may 7 bath sites (8 electrons sa 8 orbitals). Ang model na ito ay ginagamit upang ilarawan ang magnetic impurities na naka-embed sa mga metal.

Pagkatapos, gagawa tayo ng 16-qubit Trotter circuits na ginagamit upang bumuo ng quantum Krylov subspace.

# 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

Susunod, bubuo tayo ng quantum Krylov subspace gamit ang isang set ng Trotterized quantum circuits. Dito ay gumagawa tayo ng helpers para sa pagbuo ng initial (reference) state pati na rin ng time evolution para sa one- at two-body parts ng Hamiltonian. Para sa mas detalyadong paglalarawan ng model na ito at kung paano nilikha ang circuits, mangyaring sumangguni sa the 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],
        )

Bumuo ng d time-evolved states na tumutukoy sa quantum Krylov subspace. Dito, pinili natin ang d=8. Ang error mula sa pag-sample ng Krylov basis states ay nagco-converge habang tumataas ang d. Tandaan na ang mga particularidad ng problema instance na ito ay nagpapahintulot ng efficient compilation ng one-body evolution gamit ang OrbitalRotationJW; gayunpaman, sa pangkalahatan, maaaring gamitin ang PauliEvolutionGate ng Qiskit upang i-implementa ang Trotterized time evolution ng buong Hamiltonian.

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

Hakbang 2: I-optimize ang problem

Pagkatapos nating gumawa ng Trotterized circuits, i-ooptimize natin ang mga ito para sa target hardware. Kailangan nating pumili ng hardware device na gagamitin bago i-optimize. Gagamit tayo ng fake 127-qubit backend mula sa qiskit_ibm_runtime upang gayahin ang isang totoong device. Para tumakbo sa isang totoong QPU, palitan lamang ang fake backend ng totoong backend. Tingnan ang Qiskit IBM Runtime docs para sa karagdagang impormasyon.

from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke

backend = FakeSherbrooke()

Susunod, ita-transpile natin ang circuits sa target backend gamit ang 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)

Hakbang 3: Isagawa ang mga eksperimento

Pagkatapos i-optimize ang circuits para sa hardware execution, handa na tayong patakbuhin ang mga ito sa target hardware at mangolekta ng samples para sa ground state energy estimation. Dito ay gumagamit tayo ng SamplerV2 mula sa qiskit-ibm-runtime upang gayahin ang noisy samples na kinuha mula sa ibm_sherbrooke backend. Pagkatapos ay pinagsasama natin ang counts mula sa bawat isa sa Krylov basis states sa isang single counts dictionary at ipla-plot ang nangungunang 20 na pinakakaraniwang sampled bitstrings.

Tandaan: Ang Qiskit Aer ay kinakailangan upang gayahin ang samples mula sa transpiled circuits.

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)

Hakbang 4: I-post-process ang resulta

Ngayon, pinapatakbo natin ang SQD algorithm gamit ang diagonalize_fermionic_hamiltonian function. Tingnan ang API documentation para sa mga paliwanag ng arguments sa function na ito.

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

Ngayon, ipla-plot natin ang mga resulta. Ang unang plot ay nagpapakita na pagkatapos ng ilang iterations, nakukuha natin ang exact ground state energy.

Ang halimbawang ito ay sapat na maliit para makapag-explore tayo ng buong Hilbert space, tulad ng nakikita sa print statements sa itaas. Tandaan, ang buong Hilbert space ay may dimension na (num_orbitals choose nelec_a) * (num_orbitals choose nelec_b). Kaya para sa problema na ito: (8 choose 4)**2 = 4900. Ang subspace dimensions ay tumataas dahil sa enhanced configuration recovery, at gayundin sa katotohanan na ang SQD algorithm ay nagdadala ng mahahalagang configurations mula sa isang iteration patungo sa kasunod. Sa huling iteration, dia-diagonalize na natin sa buong Hilbert space.

Ang pangalawang plot ay nagpapakita ng average occupancy ng bawat spatial orbital sa lahat ng solusyon ng batches. Nakikita natin na may mataas na probability ang bawat orbital ay naglalaman ng isang 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