Využití trim SQD pro lepší výběr konfigurací

Výchozí chování diagonalize_fermionic_hamiltonian je diagonalizovat každou podvzorkovanou dávku a ponechat nejlepší výsledek. Trim SQD místo toho zachází s dávkami jako s kandidáty k prověření: každá se diagonalizuje, ořízne se na konfigurace s největší vahou a přeživší ze všech dávek se sloučí a diagonalizují společně. Tato sloučená diagonalizace je ta, jejíž energii iterace hlásí, a její vlastní oříznutý výstup slouží jako základ další iterace.

Implementace trim SQD vychází ze společného výzkumu Cleveland Clinic Foundation, RIKEN a IBM®. Pokud tuto funkci použiješ ve svém výzkumu, cituj také příslušný článek. Citaci najdeš v souboru CITATION.bib pod citačním klíčem merz2026crossing.

Další podrobnosti o kódu SQD použitém v tomto příkladu najdeš v tutoriálu o chemickém Hamiltoniánu.

Nastavení molekuly

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

_needed = {"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 numpy as np
import pyscf
import pyscf.mcscf
from qiskit_addon_sqd.counts import generate_bit_array_uniform
from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
from qiskit_addon_sqd.trim import TrimPolicy

# 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

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

WARN: Unable to to identify input symmetry using original axes.
Different symmetry axes will be used.

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

Spuštění smyčky

Jedna iterace trim SQD provede následující:

Prověřovací kolo. Rozděl bitové řetězce do dávek. Diagonalizuj každou dávku a ponech zlomek trim_ratio jejích řetězců CI, pro každý spinový sektor zvlášť, seřazených podle váhy.
Sloučené kolo. Diagonalizuj sloučené přeživší jednou. Toto je výsledek, který iterace hlásí, a její vlastní oříznuté řetězce slouží jako základ další iterace.

sci_solver smyčky provádí obě kola a ve výchozím nastavení používá řešič PySCF z tohoto balíčku.

Trim SQD povolíš předáním TrimPolicy takto:

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=300,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=5,
    max_iterations=5,
    policy=TrimPolicy(trim_ratio=0.5, max_strings_per_trim=250),
    symmetrize_spin=True,
    max_dim=250,
    seed=rng,
)

sci_state = result.sci_state
print(
    f"Merged subspace: {len(sci_state.ci_strs_a)} x {len(sci_state.ci_strs_b)} CI strings"
)
print(f"Exact energy:     {exact_energy}")
print(f"Estimated energy: {result.energy + nuclear_repulsion_energy}")

Merged subspace: 148 x 148 CI strings
Exact energy:     -109.04667177808024
Estimated energy: -109.034369147063

Ověř, zda prověřování pomáhá

Prověřování nabízí nejvíce tehdy, když vzorkované konfigurace výrazně převyšují to, co lze diagonalizovat najednou, takže seřazení většího souboru skutečně najde konfigurace, které by jedna dávka minula. Nabízí nejméně tehdy, když podprostor dosažitelné velikosti již zachycuje většinu vlnové funkce.

Prověřování závisí na tom, že se dávky navzájem liší. TrimPolicy to zajišťuje úpravou chování diagonalize_fermionic_hamiltonian tak, aby vylosoval jeden soubor samples_per_batch * num_batches bitových řetězců a rozdělil ho mezi dávky, takže se žádný bitový řetězec neobjeví ve více než jedné dávce.

Omezení podprostoru

Podprostor je zde generován kartézským součinem pole řetězců CI se spinem alfa a pole se spinem beta.

Trimování řadí konfigurace, ale přenos tvoří dva seznamy řetězců pro jednotlivé spiny. Ponechání 150 konfigurací s nejvyšší vahou z podprostoru typicky zahrnuje téměř 150 různých řetězců alfa a 150 různých řetězců beta, jejichž součin generuje zhruba 22 500 konfigurací, nikoli 150. Podprostor předaný další iteraci je tedy obecně mnohem větší než množina, která byla seřazena, a roste kvadraticky s tím, kolik ponecháš.

Trimování a přenos jsou omezeny samostatně, protože se dějí v různých bodech: trimování zmenšuje každou dávku prověřovacího kola, zatímco přenos je to, co přechází do další iterace.

Omezují ho čtyři argumenty a je užitečné nastavit alespoň jeden:

trim_ratio je zlomek řetězců každé dávky, pro každý spinový sektor, které prověřováním přežijí. Nižší hodnoty zmenšují sloučený podprostor, ale zahazují konfigurace, které musí další iterace znovu objevit vzorkováním.
max_strings_per_trim omezuje počet řetězců, které každá dávka přispěje, pro každý sektor, bez ohledu na poměr. Hodí se se symmetrize_spin, které sloučí oba spinové sektory do jednoho seznamu, a tak zhruba zdvojnásobí to, co každá dávka nabízí. Toto omezení platí pro každou dávku: sloučený podprostor může obsahovat až num_batches krát tolik řetězců, minus to, co dávky měly společné, což v praxi není mnoho.
carryover_threshold, carryover_ratio a max_carryover omezují řetězce, které slouží jako základ další iterace, a aplikují se na sloučený podprostor jako celek, nikoli na každou dávku. carryover_threshold je mez amplitudy, kterou používá StandardPolicy, se stejnou výchozí hodnotou; další dva omezují počet, což samotná mez nedokáže.
max_dim v samotné diagonalize_fermionic_hamiltonian omezuje každý spinový sektor každého podprostoru, který smyčka diagonalizuje, včetně sloučeného. Působí jako záchranná pojistka a vztahuje se na vše, co politika vrátí, poté co je přenos sloučen s novými vzorky.

Výše uvedený běh nastavuje současně max_strings_per_trim=250 a max_dim=250, takže se nemůže vymknout kontrole ani to, co každá dávka přispívá, ani podprostor z toho sestavený.

Použití externího eigensolveru

Trim SQD můžeš použít s externím eigensolverem předáním argumentu sci_solver:

from sbd.sbd_solver import create_sbd_solver

result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=300,
    norb=num_orbitals,
    nelec=(num_elec_a, num_elec_b),
    num_batches=5,
    policy=TrimPolicy(trim_ratio=0.1),
    sci_solver=create_sbd_solver(sbd_config={"method": 0, "eps": 1e-10}),
    seed=rng,
)

Předchozí úryvek vyžaduje instalaci balíčku sbd-eigensolver z PyPI.

Rozdělení práce

S TrimPolicy smyčka volá sci_solver dvakrát za iteraci: jednou se všemi dávkami a jednou se sloučeným podprostorem. Obě volání procházejí běžným rozhraním sci_solver, které vícenásobné procesní podpora tohoto balíčku vyvolává kolektivně na každém ranku MPI.

K dispozici jsou potenciálně dvě úrovně paralelismu a jsou nezávislé:

V rámci jedné diagonalizace. Řešič jako SBD může rozložit jedinou diagonalizaci mezi ranky.
Napříč dávkami. Řešič by mohl v principu diagonalizovat dávky souběžně.

Pro danou sadu dávek řešič určuje, jak práci přidělí svým pracovním procesům.

Nastavení molekuly​

Spuštění smyčky​

Ověř, zda prověřování pomáhá​

Omezení podprostoru​

Použití externího eigensolveru​

Rozdělení práce​