Trim SQD für eine verbesserte Konfigurationsauswahl nutzen

Das Standardverhalten von diagonalize_fermionic_hamiltonian besteht darin, jeden subgesampelten Batch zu diagonalisieren und das beste Ergebnis zu behalten. Trim SQD behandelt die Batches stattdessen als Kandidaten, die geprüft werden: Jeder wird diagonalisiert und auf die Konfigurationen mit dem größten Gewicht getrimmt, und die Überlebenden aller Batches werden zusammengeführt und gemeinsam diagonalisiert. Diese zusammengeführte Diagonalisierung liefert die Energie, die die Iteration meldet, und ihre eigene getrimmte Ausgabe dient als Ausgangspunkt der nächsten Iteration.

Die Trim-SQD-Implementierung beruht auf gemeinsamer Forschung der Cleveland Clinic Foundation, von RIKEN und IBM®. Wenn du diese Funktion in deiner Forschung verwendest, zitiere bitte auch die Arbeit. Die Referenz findest du in der Datei CITATION.bib unter dem Zitationsschlüssel merz2026crossing.

Weitere Details zum in diesem Beispiel verwendeten SQD-Code findest du im Tutorial zum chemischen Hamiltonoperator.

Das Molekül einrichten

# 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

Die Schleife ausführen

Eine Trim-SQD-Iteration führt Folgendes aus:

1. Screening-Runde. Die Bitstrings in Batches aufteilen. Jeden Batch diagonalisieren und pro Spinsektor den Anteil trim_ratio seiner CI-Strings behalten, geordnet nach Gewicht.

2. Zusammengeführte Runde. Die zusammengeführten Überlebenden einmal diagonalisieren. Das ist das Ergebnis, das die Iteration meldet, und ihre eigenen getrimmten Strings dienen als Ausgangspunkt der nächsten Iteration.

Der sci_solver der Schleife führt beide Runden aus und ist standardmäßig der PySCF-Solver dieses Pakets.

Aktiviere Trim SQD, indem du die TrimPolicy wie folgt übergibst:

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

Prüfen, ob das Screening hilft

Das Screening bringt am meisten, wenn die gesampelten Konfigurationen die Zahl derer, die auf einmal diagonalisiert werden können, weit übersteigen, sodass das Ranking eines größeren Pools tatsächlich Konfigurationen findet, die ein einzelner Batch übersehen hätte. Es bringt am wenigsten, wenn ein Unterraum der bezahlbaren Größe bereits den Großteil der Wellenfunktion erfasst.

Das Screening hängt davon ab, dass sich die Batches voneinander unterscheiden. TrimPolicy sorgt dafür, indem sie das Verhalten von diagonalize_fermionic_hamiltonian so ändert, dass ein Pool von samples_per_batch * num_batches Bitstrings gezogen und auf die Batches aufgeteilt wird, sodass kein Bitstring in mehr als einem Batch vorkommt.

Den Unterraum begrenzen

Ein Unterraum wird hier vom kartesischen Produkt eines Spin-Alpha- und eines Spin-Beta-CI-String-Arrays aufgespannt.

Das Trimmen ordnet Konfigurationen nach Rang, aber der Übertrag besteht aus zwei String-Listen pro Spin. Die 150 Konfigurationen mit dem höchsten Gewicht eines Unterraums zu behalten, umfasst typischerweise nahezu 150 verschiedene Alpha-Strings und 150 verschiedene Beta-Strings, deren Produkt etwa 22.500 Konfigurationen aufspannt, nicht 150. Der Unterraum, der an die nächste Iteration weitergegeben wird, ist daher im Allgemeinen deutlich größer als die Menge, die geordnet wurde, und er wächst quadratisch mit dem, was du behältst.

Trimmen und Übertragen werden getrennt begrenzt, weil sie an verschiedenen Stellen stattfinden: Das Trimmen verkleinert jeden Batch der Screening-Runde, während der Übertrag das ist, was in die nächste Iteration übergeht.

Vier Argumente begrenzen ihn, und es ist sinnvoll, mindestens eines zu setzen:

• trim_ratio ist der Anteil der Strings jedes Batches pro Spinsektor, der das Screening überlebt. Niedrigere Werte verkleinern den zusammengeführten Unterraum, verwerfen aber Konfigurationen, die die nächste Iteration durch Sampling erneut entdecken muss.

• max_strings_per_trim begrenzt die Strings, die jeder Batch pro Sektor beiträgt, unabhängig vom Verhältnis. Nützlich mit symmetrize_spin, das beide Spinsektoren in einer Liste zusammenführt und damit ungefähr verdoppelt, was jeder Batch anbietet. Diese Grenze gilt pro Batch: Der zusammengeführte Unterraum kann bis zu num_batches-mal so viele Strings enthalten, abzüglich dessen, was die Batches gemeinsam hatten, was in der Praxis nicht viel ist.

• carryover_threshold, carryover_ratio und max_carryover begrenzen die Strings, die die nächste Iteration anstoßen, und werden auf den zusammengeführten Unterraum als Ganzes angewendet, nicht pro Batch. carryover_threshold ist der Amplitudenschwellenwert, den StandardPolicy verwendet, mit demselben Standardwert; die anderen beiden begrenzen die Anzahl, was der Schwellenwert allein nicht leistet.

• max_dim an diagonalize_fermionic_hamiltonian selbst begrenzt jeden Spinsektor jedes Unterraums, den die Schleife diagonalisiert, einschließlich des zusammengeführten. Es wirkt als Absicherung und gilt für alles, was die Policy zurückgibt, nachdem der Übertrag mit frischen Samples zusammengeführt wurde.

Der obige Lauf setzt max_strings_per_trim=250 und max_dim=250 gemeinsam, sodass weder das, was jeder Batch beiträgt, noch der daraus aufgebaute Unterraum außer Kontrolle geraten kann.

Einen externen Eigensolver verwenden

Du kannst Trim SQD mit einem externen Eigensolver verwenden, indem du das Argument sci_solver übergibst:

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

Das vorherige Snippet setzt voraus, dass das Paket sbd-eigensolver von PyPI installiert ist.

Die Arbeit verteilen

Mit TrimPolicy ruft die Schleife sci_solver pro Iteration zweimal auf: einmal mit allen Batches und einmal mit dem zusammengeführten Unterraum. Beide Aufrufe laufen über die gewöhnliche sci_solver-Schnittstelle, die die Multiprozess-Unterstützung dieses Pakets kollektiv auf jedem MPI-Rank aufruft.

Potenziell stehen zwei Ebenen der Parallelität zur Verfügung, und sie sind voneinander unabhängig:

• Innerhalb einer Diagonalisierung. Ein Solver wie SBD kann eine einzelne Diagonalisierung auf mehrere Ranks verteilen.

• Über die Batches hinweg. Ein Solver könnte im Prinzip die Batches gleichzeitig diagonalisieren.

Bei gegebener Menge von Batches bestimmt der Solver, wie er die Arbeit an seine Worker-Prozesse delegiert.