{ "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, also cite the paper. You can find the citation in the\n", "`CITATION.bib` file, with citation key `merz2026crossing`.\n", "\n", "For more details on the SQD code used in this example, see the [chemistry Hamiltonian tutorial](https://quantum.cloud.ibm.com/docs/tutorials/sample-based-quantum-diagonalization)." ] }, { "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", "\n", "One trim SQD iteration does the following:\n", "\n", "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.\n", "2. **Merged round.** Diagonalize the merged survivors once. This is the result the iteration reports, and its own trimmed strings seed the next iteration.\n", "\n", "The loop's `sci_solver` performs both rounds, and defaults to this package's PySCF solver.\n", "\n", "Enable trim SQD 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": [ "## Check 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. `TrimPolicy` manages this by modifying the behavior of `diagonalize_fermionic_hamiltonian` to draw one pool of `samples_per_batch * num_batches` bitstrings and partitioning it among the batches, so no bitstring appears in more than one." ] }, { "cell_type": "markdown", "id": "b5c6d7eb", "metadata": {}, "source": [ "## Bound the subspace\n", "\n", "A subspace here is spanned by the *Cartesian product* of a spin-alpha and a spin-beta CI string array.\n", "\n", "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.\n", "\n", "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.\n", "\n", "Four arguments bound it, and it is useful to set 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. 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, including the merged one. It acts as a backstop and applies to whatever the policy returns, after the carryover is merged with fresh samples.\n", "\n", "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." ] }, { "cell_type": "markdown", "id": "c04f21a7", "metadata": {}, "source": [ "## Use an external eigensolver\n", "\n", "You can use trim SQD with an external eigensolver by passing the `sci_solver` argument:\n", "\n", "```python\n", "from sbd.sbd_solver import create_sbd_solver\n", "\n", "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", " 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", "\n", "The prior snippet requires the `sbd-eigensolver` package to be installed from PyPI." ] }, { "cell_type": "markdown", "id": "c6d7e8fc", "metadata": {}, "source": [ "## Distribute 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", "Two levels of parallelism are potentially available, and they are independent:\n", "\n", "- **Within a diagonalization.** A solver such as SBD can 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, the solver determines how to delegate work 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 }