{ "cells": [ { "cell_type": "markdown", "id": "0a1f5c30", "metadata": {}, "source": [ "# Leverage trim SQD for improved configuration selection\n", "\n", "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.\n", "\n", "The trim SQD implementation is based on joint research from\n", "the Cleveland Clinic Foundation, RIKEN, and IBM. If you use this feature in\n", "your research, please also use the citation for the paper, which can be found\n", "in the `CITATION.bib` file with citation key `merz2026crossing`.\n", "\n", "For more details on the SQD code used in this example, check out [tutorial 1](https://qiskit.github.io/qiskit-addon-sqd/tutorials/01_chemistry_hamiltonian.html)." ] }, { "cell_type": "markdown", "id": "2c3d7e52", "metadata": {}, "source": [ "## Set up the molecule" ] }, { "cell_type": "code", "execution_count": 1, "id": "3d4e8f63", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "WARN: Unable to to identify input symmetry using original axes.\n", "Different symmetry axes will be used.\n", "\n", "converged SCF energy = -108.835236570774\n", "CASCI E = -109.046671778080 E(CI) = -32.8155692383188 S^2 = 0.0000000\n" ] } ], "source": [ "import numpy as np\n", "import pyscf\n", "import pyscf.mcscf\n", "from qiskit_addon_sqd.counts import generate_bit_array_uniform\n", "from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian\n", "from qiskit_addon_sqd.trim import TrimPolicy\n", "\n", "# Build N2 molecule\n", "mol = pyscf.gto.Mole()\n", "mol.build(\n", " atom=[[\"N\", (0, 0, 0)], [\"N\", (1.0, 0, 0)]],\n", " basis=\"6-31g\",\n", " symmetry=\"Dooh\",\n", ")\n", "\n", "# Define active space\n", "n_frozen = 2\n", "active_space = range(n_frozen, mol.nao_nr())\n", "\n", "# Get molecular integrals\n", "scf = pyscf.scf.RHF(mol).run()\n", "num_orbitals = len(active_space)\n", "n_electrons = int(sum(scf.mo_occ[active_space]))\n", "num_elec_a = (n_electrons + mol.spin) // 2\n", "num_elec_b = (n_electrons - mol.spin) // 2\n", "cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))\n", "mo = cas.sort_mo(active_space, base=0)\n", "hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n", "eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)\n", "\n", "# Compute exact energy\n", "exact_energy = cas.run().e_tot\n", "\n", "# Create a seed to control randomness throughout this workflow\n", "rng = np.random.default_rng(24)\n", "\n", "# Generate random samples\n", "bit_array = generate_bit_array_uniform(10_000, num_orbitals * 2, rand_seed=rng)" ] }, { "cell_type": "markdown", "id": "5f60a185", "metadata": {}, "source": "## Run the loop\n\nOne trim SQD iteration does the following:\n\n1. **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.\n2. **Merged round.** Diagonalize the merged survivors once. This is the result the iteration reports, and its own trimmed strings seed the next iteration.\n\nBoth rounds are performed by the loop's `sci_solver`, which defaults to this package's PySCF solver.\n\nTrimSQD is enabled by passing the `TrimPolicy` as follows:" }, { "cell_type": "code", "execution_count": 2, "id": "60718296", "metadata": {}, "outputs": [], "source": "result = diagonalize_fermionic_hamiltonian(\n hcore,\n eri,\n bit_array,\n samples_per_batch=300,\n norb=num_orbitals,\n nelec=(num_elec_a, num_elec_b),\n num_batches=5,\n max_iterations=5,\n policy=TrimPolicy(trim_ratio=0.5, max_strings_per_trim=250),\n symmetrize_spin=True,\n max_dim=250,\n seed=rng,\n)" }, { "cell_type": "code", "execution_count": 3, "id": "8293a4b8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Merged subspace: 148 x 148 CI strings\n", "Exact energy: -109.04667177808024\n", "Estimated energy: -109.034369147063\n" ] } ], "source": [ "sci_state = result.sci_state\n", "print(f\"Merged subspace: {len(sci_state.ci_strs_a)} x {len(sci_state.ci_strs_b)} CI strings\")\n", "print(f\"Exact energy: {exact_energy}\")\n", "print(f\"Estimated energy: {result.energy + nuclear_repulsion_energy}\")" ] }, { "cell_type": "markdown", "id": "a4b5c6da", "metadata": {}, "source": [ "## Checking whether screening is helping\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "markdown", "id": "b5c6d7eb", "metadata": {}, "source": "## Bounding the subspace\n\nA subspace here is spanned by the *Cartesian product* of a spin-alpha and a spin-beta CI string array.\n\nTrimming 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.\n\nTrimming 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.\n\nFour arguments bound it, and it is worth setting at least one:\n\n- **`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.\n- **`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.\n- **`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.\n- **`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.\n\nThe 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." }, { "cell_type": "markdown", "id": "c04f21a7", "metadata": {}, "source": "## Using an external eigensolver\n\nTrimSQD can be used with an external eigensolver by passing the `sci_solver` argument:\n\n```python\nfrom sbd.sbd_solver import create_sbd_solver\n\nresult = diagonalize_fermionic_hamiltonian(\n hcore,\n eri,\n bit_array,\n samples_per_batch=300,\n norb=num_orbitals,\n nelec=(num_elec_a, num_elec_b),\n num_batches=5,\n policy=TrimPolicy(trim_ratio=0.1),\n sci_solver=create_sbd_solver(sbd_config={\"method\": 0, \"eps\": 1e-10}),\n seed=rng,\n)\n```\n\nThe prior snippet requires the `sbd-eigensolver` package to be installed from PyPI." }, { "cell_type": "markdown", "id": "c6d7e8fc", "metadata": {}, "source": [ "## Distributing the work\n", "\n", "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](hpc_acceleration.ipynb) invokes collectively on every MPI rank.\n", "\n", "There are potentially two levels of parallelism available, and they are independent:\n", "\n", "- **Within a diagonalization.** A solver such as SBD is able to spread a single diagonalization across ranks.\n", "- **Across the batches.** A solver could, in principle, diagonalize the batches concurrently.\n", "\n", "Given a set of batches, it is up to the solver to determine the way in which work is delegated to its worker processes." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.7" } }, "nbformat": 4, "nbformat_minor": 5 }