Leverage trim SQD for improved configuration selection

The default behavior of `diagonalize_fermionic_hamiltonian <../apidocs/qiskit_addon_sqd.fermion.rst>`__ is to diagonalize each subsampled batch and keep the best result. Trim SQD instead treats the batches as candidates to be screened: each is diagonalized, trimmed down to the configurations carrying the largest weight, and the survivors of all the batches are merged and diagonalized together. That merged diagonalization is the one whose energy the iteration reports, and its own trimmed output seeds the next iteration.

The trim SQD implementation is based on joint research from the Cleveland Clinic Foundation, RIKEN, and IBM. If you use this feature in your research, please also use the citation for the paper, which can be found in the CITATION.bib file with citation key merz2026crossing.

For more details on the SQD code used in this example, check out tutorial 1.

Set up the molecule

[1]:

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

Run the loop

One trim SQD iteration does the following:

  1. Screening round. Partition the bitstrings into batches. Diagonalize each batch and keep the trim_ratio fraction of its CI strings, per spin sector, ranked by weight.

  2. Merged round. Diagonalize the merged survivors once. This is the result the iteration reports, and its own trimmed strings seed the next iteration.

Both rounds are performed by the loop’s sci_solver, which defaults to this package’s PySCF solver.

TrimSQD is enabled by passing the TrimPolicy as follows:

[2]:

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

[3]:

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

Checking whether screening is helping

Screening has the most to offer when the sampled configurations greatly outnumber what can be diagonalized at once, so that ranking a larger pool genuinely finds configurations a single batch would have missed. It has the least to offer when a subspace of the affordable size already captures most of the wavefunction.

Screening depends on the batches differing from one another, and TrimPolicy takes care of that: it modifies the behavior of diagonalize_fermionic_hamiltonian to draw one pool of samples_per_batch * num_batches bitstrings and partition it among the batches, so no bitstring appears in more than one.

Bounding the subspace

A subspace here is spanned by the Cartesian product of a spin-alpha and a spin-beta CI string array.

Trimming ranks configurations, but the carryover is two per-spin string lists. Keeping the 150 highest-weight configurations of a subspace typically involves close to 150 distinct alpha strings and 150 distinct beta strings — whose product spans roughly 22,500 configurations, not 150. So the subspace handed to the next iteration is generally much larger than the set that was ranked, and it grows quadratically in what you retain.

Trimming and carrying over are bounded separately, because they happen at different points: trimming shrinks each batch of the screening round, while the carryover is what crosses into the next iteration.

Four arguments bound it, and it is worth setting at least one:

  • ``trim_ratio`` is the fraction of each batch’s strings, per spin sector, that survive screening. Lower values shrink the merged subspace but discard configurations the next iteration must rediscover by sampling.

  • ``max_strings_per_trim`` caps the strings each batch contributes, per sector, regardless of the ratio. Useful with symmetrize_spin, which merges both spin sectors into one list and so roughly doubles what each batch offers. Note that this bound is per batch: the merged subspace can hold up to num_batches times as many strings, less whatever the batches had in common, which in practice is not much.

  • ``carryover_threshold``, ``carryover_ratio`` and ``max_carryover`` bound the strings that seed the next iteration, applied to the merged subspace as a whole rather than per batch. carryover_threshold is the amplitude cutoff that StandardPolicy uses, with the same default; the other two bound the count, which the cutoff alone does not.

  • ``max_dim``, on diagonalize_fermionic_hamiltonian itself, bounds each spin sector of every subspace the loop diagonalizes, the merged one included. This is the backstop: it applies to whatever the policy returns, after the carryover is merged with fresh samples.

The run above sets max_strings_per_trim=250 and max_dim=250 together, so neither what each batch contributes nor the subspace built from it can run away.

Using an external eigensolver

TrimSQD can be used with an external eigensolver by passing the sci_solver argument:

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

The prior snippet requires the sbd-eigensolver package to be installed from PyPI.

Distributing the work

With TrimPolicy, the loop calls sci_solver twice per iteration: once with all the batches, and once with the merged subspace. Both calls go through the ordinary sci_solver interface, which this package’s multi-process support invokes collectively on every MPI rank.

There are potentially two levels of parallelism available, and they are independent:

  • Within a diagonalization. A solver such as SBD is able to spread a single diagonalization across ranks.

  • Across the batches. A solver could, in principle, diagonalize the batches concurrently.

Given a set of batches, it is up to the solver to determine the way in which work is delegated to its worker processes.