Paglimita sa dimensyon ng subspace

Sa walkthrough na ito, ipapakita namin ang epekto ng dimensyon ng subspace sa self-consistent configuration recovery technique.

A priori, hindi natin alam kung ano ang tamang dimensyon ng subspace para makuha ang target na antas ng katumpakan. Gayunpaman, alam nating pinapataas ng pagtaas ng dimensyon ng subspace ang katumpakan ng pamamaraan. Kaya, maaari nating pag-aralan ang katumpakan ng mga prediksyon bilang function ng dimensyon ng subspace. Tukuyin ang molecule at ang mga katangian nito.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "pyscf": "pyscf", "qiskit_addon_sqd": "qiskit-addon-sqd"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import warnings

import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)

# Compute exact energy
exact_energy = cas.run().e_tot

converged SCF energy = -108.835236570775
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

Gumawa ng ilang random na bitstring bilang proxy ng mga sample ng QPU.

import numpy as np
from qiskit_addon_sqd.counts import generate_bit_array_uniform

# Create a seed to control randomness throughout this workflow
rng = np.random.default_rng(24)

# Generate random samples
bit_array = generate_bit_array_uniform(
    10_000, num_orbitals * 2, rand_seed=rng
)

Tawagin ang SQD na may pataas na laki ng batch.

from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian

list_samples_per_batch = [50, 200, 400, 600]

# SQD options
max_iterations = 5

# Eigenstate solver options
num_batches = 10
max_davidson_cycles = 200

energies = []
subspace_dimensions = []

for samples_per_batch in list_samples_per_batch:
    result = diagonalize_fermionic_hamiltonian(
        hcore,
        eri,
        bit_array,
        samples_per_batch=samples_per_batch,
        norb=num_orbitals,
        nelec=(num_elec_a, num_elec_b),
        num_batches=num_batches,
        max_iterations=max_iterations,
        symmetrize_spin=True,
        seed=rng,
    )
    energies.append(result.energy)
    subspace_dimensions.append(np.prod(result.sci_state.amplitudes.shape))

Ipinapakita ng plot na ito na ang pagtaas ng dimensyon ng subspace ay humahantong sa mas tumpak na mga resulta.

import matplotlib.pyplot as plt

# Data for energies plot
x1 = subspace_dimensions
y1 = np.array(energies) + nuclear_repulsion_energy

fig, axs = plt.subplots(1, 1, figsize=(12, 6))

# Plot energies
axs.plot(x1, y1, marker=".", markersize=20, label="Estimated")
axs.set_xticks(x1)
axs.set_xticklabels(x1)
axs.axhline(y=exact_energy, color="red", linestyle="--", label="Exact")
axs.set_title("Approximated Ground State Energy vs subspace dimension")
axs.set_xlabel("Subspace dimension")
axs.set_ylabel("Energy (Ha)")
axs.legend()

plt.tight_layout()
plt.show()

Pagtatakda ng tahasang limitasyon gamit ang max_dim

Sa sweep sa itaas, hindi direktang kinontrol ang dimensyon ng subspace. Itinatakda ng samples_per_batch kung ilang bitstring ang kukunin sa bawat batch, ngunit ang dimensyon ng nagreresultang subspace ay natutukoy ng kung ilang magkakaibang spin-alpha at spin-beta string ang nilalaman ng mga bitstring na iyon, at ang pool ng mga configuration na magagamit para sa subsampling ay lumalaki mula sa isang iteration ng configuration recovery patungo sa susunod. Iyan ang dahilan kung bakit hindi regular na mga numero ang mga dimensyon ng subspace sa horizontal axis sa itaas na hindi natin pinili.

Direktang nililimitahan ng argumentong max_dim ng diagonalize_fermionic_hamiltonian ang dimensyon. Nililimitahan nito ang bilang ng mga CI string na pinapanatili sa bawat spin sector, at dahil ang subspace ay binubuo ng Cartesian product ng dalawang sector, ang dimensyon nito ay ang product ng dalawang dimensyon ng sector:

• Nililimitahan ng max_dim=d ang parehong sector sa d na mga string, kaya ang dimensyon ng subspace ay hindi hihigit sa d ** 2.

• Hiwalay na nililimitahan ng max_dim=(d_a, d_b) ang mga spin-alpha at spin-beta sector, kaya ang dimensyon ng subspace ay hindi hihigit sa d_a * d_b. Kapag symmetrize_spin=True, dapat pantay ang dalawang limitasyon.

• Ang max_dim=None, ang default, ay hindi nagtatakda ng anumang limitasyon.

Kapag ang isang batch ay may mas maraming string kaysa sa pinapayagan ng limitasyon, itinatapon ang mga pinakakaunting mahalaga: unang pinapanatili ang mga string na ipinasa sa pamamagitan ng include_configurations, kasunod ang mga string na dinala mula sa nakaraang iteration (tingnan ang carryover_threshold), at pinupunan ang natitirang slot ng mga na-sample na string ayon sa pababang pagkakasunod-sunod ng dalas ng pag-sample sa mga ito.

Dahil nililimitahan ng max_dim ang laki ng pinakamalaking diagonalization na hihilingin sa eigenstate solver na gawin, isa itong maginhawang paraan para panatilihin ang solver sa loob ng badyet ng runtime at memory: piliin nang sagana ang samples_per_batch, at hayaang ang max_dim ang magpasya kung gaano kalaki ang maaaring lumaki ang subspace.

Sa ibaba, uulitin natin ang pag-aaral, ngunit sa pagkakataong ito ay pinananatiling nakapirmi ang samples_per_batch habang pinapataas ang max_dim.

list_max_dim = [100, 200, 400, 800]

# Use a generous number of samples per batch, so that max_dim is what limits
# the dimension of the subspace
samples_per_batch = 600

energies_max_dim = []
subspace_dimensions_max_dim = []

for max_dim in list_max_dim:
    result = diagonalize_fermionic_hamiltonian(
        hcore,
        eri,
        bit_array,
        samples_per_batch=samples_per_batch,
        norb=num_orbitals,
        nelec=(num_elec_a, num_elec_b),
        num_batches=num_batches,
        max_iterations=max_iterations,
        symmetrize_spin=True,
        max_dim=max_dim,
        seed=rng,
    )
    energies_max_dim.append(result.energy)
    subspace_dimensions_max_dim.append(
        int(np.prod(result.sci_state.amplitudes.shape))
    )

print(f"max_dim:            {list_max_dim}")
print(f"subspace dimension: {subspace_dimensions_max_dim}")

max_dim:            [100, 200, 400, 800]
subspace dimension: [10000, 40000, 160000, 640000]

Ang dimensyon ng subspace ay eksaktong max_dim ** 2 na ngayon sa bawat takbo: pinili ito sa halip na matuklasan. Ipinapakita ng pagpa-plot sa dalawang sweep nang magkasama na ang katumpakan ay nakadepende sa mismong dimensyon ng subspace, anuman ang argumentong ginamit para marating ito. Gayunpaman, hindi kailangang eksaktong magkapareho ang dalawang curve. Ang mga takbo ng max_dim ay kumukuha ng 600 bitstring bawat batch at pagkatapos ay pinapanatili ang mga pinakamadalas na na-sample na string, kaya sa isang partikular na dimensyon ay maaari silang magpanatili ng mas magandang hanay ng mga configuration kaysa sa takbong umabot sa dimensyong iyon sa pamamagitan ng pag-sample ng mas kaunting bitstring sa simula pa lang.

fig, axs = plt.subplots(1, 1, figsize=(12, 6))

axs.plot(
    subspace_dimensions,
    np.array(energies) + nuclear_repulsion_energy,
    marker=".",
    markersize=20,
    label="Varying samples_per_batch",
)
axs.plot(
    subspace_dimensions_max_dim,
    np.array(energies_max_dim) + nuclear_repulsion_energy,
    marker="^",
    markersize=12,
    label="Varying max_dim",
)
axs.axhline(y=exact_energy, color="red", linestyle="--", label="Exact")
axs.set_xscale("log")
axs.ticklabel_format(axis="y", useOffset=False)
axs.set_title("Approximated Ground State Energy vs subspace dimension")
axs.set_xlabel("Subspace dimension")
axs.set_ylabel("Energy (Ha)")
axs.legend()

plt.tight_layout()
plt.show()