12Representing and compressing the Hamiltonian

Keep the model fixed while changing its numerical representation

The nitrogen story now has a precisely identified molecular geometry, orbital basis, active-electron partition, and correlated reference strategy. We can ask a new engineering question: how should we store and manipulate the Hamiltonian of that chosen model?

Five operations must remain distinct. An exact change of storage preserves the numerical operator. Truncating a factorization approximates that operator. Selecting active orbitals changes which electron occupations the model permits. Downfolding would incorporate excluded-space effects through a specifically defined effective operator; as Lesson 11 established, its QDK interface currently has no bundled concrete implementation. Qubit encoding, taught in Lesson 14, represents a fermionic operator with qubit operators.

A low-rank tensor is not automatically a better chemical model. It can be a cheaper representation of the same model, or a controlled numerical approximation that needs its own error assessment. It does not restore correlation omitted by a six-orbital active space or repair an inappropriate spin reference.

All API examples in this lesson are unexecuted. Algebraic examples are illustrative, not measured nitrogen factor ranks or energies. We keep the 1.85-Å tutorial separate from the 1.270025-Å notebook and other regression molecules.

The Hamiltonian contains numbers with several different index meanings

For real spatial orbitals, the one-electron integral $h_{pq}$ describes kinetic energy and the external potential between orbitals $p$ and $q$. In a frozen-core active Hamiltonian it also includes the inactive-electron field. The two-electron integral in chemist notation is

The first pair shares electron coordinate one; the second pair shares coordinate two. Physicist notation instead groups a bra pair and a ket pair: $v_{ijkl}=\langle ij|kl\rangle=(ik|jl)=g_{ikjl}$. Switching between these notations permutes indices; it is not a different Coulomb interaction.

For example, with two orbitals, chemist entry $g_{0011}$ is the same number as physicist entry $v_{0101}$. If each four-index array is flattened with its final index fastest, the first entry occupies offset three while the second occupies offset five. Reading the same flat buffer under the wrong notation therefore changes the interpreted Hamiltonian.

The concrete canonical container uses chemist-ordered (i,j,k,l) entries with offset

where $K$ is its active spatial-orbital count. The element accessor implements that mapping directly; it does not automatically translate an arbitrary caller’s physicist-order request. Some older RST text says physicist notation is the default, so this is a case where implementation matters more than the overview.

Indices in the active integral arrays are positions within the active list. If the full MO indices are [4, 5, 6, 7, 8, 9], active position zero refers to original orbital four. get_two_body_element(0,0,5,5) then refers to the first and last orbitals of that active list. Passing original index nine into a six-orbital tensor is not a request for automatic remapping.

Spin blocks, symmetries, and constants are part of the representation

The common one-body getter returns alpha and beta matrices. The common two-body getter returns three arrays labeled aaaa, aabb, and bbbb. For a restricted common spatial-orbital representation these can share data. For unrestricted orbitals, different alpha/beta coefficients give different same-spin and mixed-spin integrals, even though the underlying Coulomb interaction is spin independent.

The native four-center construction describes its mixed array as $(\alpha\alpha|\beta\beta)$; the Cholesky reconstruction and its dense-equivalence tests use the corresponding left-alpha/right-beta factors. The low-level transformation call can pass beta coefficient blocks first because of its internal layout, so call order is not a reliable statement of returned spin order. The constructor, actual flattened access, and matching producer/consumer tests must agree. Do not infer a universal beta-first storage convention from a legacy contraction comment.

For real Coulomb integrals, swapping members inside either orbital pair or exchanging the two pairs gives the familiar permutation symmetries. For complex orbitals some identities require conjugation and are not interchangeable real-array equalities. Antisymmetrized interaction tensors also have different signs and structure. The molecular containers and factorization algorithms traced here store real matrices/vectors; their mathematical assumptions must not be extended automatically to every complex or effective interaction.

get_core_energy() is the zero-body constant: it shifts every eigenvalue in the model by the same amount. For an all-active all-electron model it includes nuclear repulsion. With frozen inactive orbitals it includes their energy; with an ECP it uses the appropriate effective-charge convention and ECP one-electron contribution. It is not necessarily “the energy of the nuclei alone.” An inactive Fock matrix, when supplied, records the field of the frozen partition and can have full MO dimensions, even when active integrals have smaller dimensions.

HamiltonianType.Hermitian or NonHermitian is part of the data description. A metadata label is not a proof that every numerical array satisfies every symmetry or that a chosen solver supports non-Hermitian input. The energy and RDM conventions from Lesson 9 remain in force.

The facade preserves one API while containers choose the storage

Hamiltonian delegates to a concrete container. get_container_type() identifies it, and Python get_container() gives access to container-specific information. The shared get_one_body_integrals, get_two_body_integrals, get_one_body_element, get_two_body_element, and availability methods are convenient, but their computational costs can differ enormously.

These are existing concrete containers, not four interchangeable algorithms that all consume the same constructor arguments. In particular, a sparse model container constructs ModelOrbitals; simply moving molecular arrays into it does not preserve all molecular-basis metadata automatically.

The defining data are supplied at construction. The current interfaces do not provide a general “set this one integral in place” workflow like some older setter-oriented prose suggests. Construct a new representation deliberately. Internal lazy reconstruction caches are different from user mutation of the defining operator.

Canonical four-center storage is simple and expensive

With $K$ active spatial orbitals, one full two-electron block has $K^4$ real entries. For the six-orbital model, that is 1296 entries; for the original 28-orbital space it is 614656. Those counts are arithmetic, not measured allocations for a complete calculation. Multiple spin blocks, intermediate transforms, solver matrices and cached properties add other costs.

Real permutation symmetry can reduce the number of independent values, but the canonical container exposes full flattened blocks rather than a universal eightfold-minimal packed tensor. Its symmetry-blocked storage can alias equivalent spin blocks for a restricted representation. That saves duplicated spin data without making every four-index entry disappear.

The advantage is direct lookup and simple contractions. Many established CI and correlation codes naturally consume these arrays. The disadvantage is quartic storage and expensive AO-to-MO transformation or repeated many-index contractions. The fact that a Hamiltonian is small in active orbital count does not imply that its original AO integral construction was equally small.

Cholesky uses the Coulomb tensor’s pair structure

Combine each pair of orbital indices into one composite index, $I=(p,q)$, and form the supermatrix $G_{I,J}=g_{pqrs}$, with $J=(r,s)$. A supermatrix here is just an ordinary matrix whose labels happen to be index pairs.

For the real Coulomb kernel, this pair matrix is positive semidefinite: the Coulomb quadratic form of any real pair-density combination is nonnegative. Positive semidefinite (PSD) means $x^{\mathsf T}Gx\ge0$ for every real vector $x$. It can have zero eigenvalues, especially because reversed orbital pairs are redundant.

A pivoted Cholesky decomposition approximates

$R$ is the retained Cholesky rank, not electron count or number of atoms. The factors can be stored as a $K^2$-by-$R$ matrix. The term “three-center” in this container describes two orbital indices plus an auxiliary index; a Cholesky auxiliary need not be a physical third atom or a separately selected chemical basis function.

The algorithm chooses the largest remaining residual diagonal as its next pivot. If previous columns are known and pivot $k$ has residual diagonal $D_k$, the next column is

The column removes the strongest remaining pair direction; repeated updates shrink the residual. A threshold controls when the remaining diagonal is judged negligible. This is very different from selecting determinants or deleting orbitals: every active orbital can remain present in the factors.

A purely illustrative two-orbital example takes

In pair order (00,01,10,11), its vector is (1,0,0,2). The outer product has nonzero entries $g_{0000}=1$, $g_{0011}=g_{1100}=2$, and $g_{1111}=4$, and has rank one. One Cholesky column reconstructs the tensor exactly. These invented values demonstrate structure, not a molecular integral set obtained in this session.

Construct Cholesky factors directly from molecular orbitals

The registered implementation hamiltonian_constructor/qdk_cholesky takes Orbitals, rather than an already constructed Hamiltonian. Its implementation evaluates needed AO Coulomb shell blocks through Libint while performing pivoted decomposition. It avoids first storing the complete AO four-index tensor.

The AO factors are transformed to active MOs by

Unrestricted input uses separate alpha and beta coefficient transforms with a shared auxiliary space. The constructor requires nonempty active spaces, matching alpha/beta active dimensions, and valid basis/index information; it handles noncontiguous selected columns explicitly. Native orbital/basis conversion restrictions from earlier lessons still apply.

Its actual settings are cholesky_tolerance=1e-8, eri_threshold=1e-12, store_ao_cholesky_vectors=False, and cholesky_gemm_batch_cols=20. The first controls the Cholesky residual stopping scale; the second independently screens small integral shell contributions. The implementation raises an excessively small Cholesky tolerance to a floor proportional to the square root of the maximum AO pair rank times machine epsilon. A nominal zero tolerance is therefore not unlimited exact arithmetic.

cholesky_gemm_batch_cols controls batching of matrix multiplication work and must be positive; it is not a desired rank or chemistry-error tolerance. store_ao_cholesky_vectors=True retains the original AO factors for possible reuse but increases stored data. The actual auxiliary rank after transforming to a smaller active space can be larger than the minimal rank of that active tensor; projection does not automatically recompress the auxiliary list.

The constructor evaluates the one-electron kinetic/attraction terms and adds ECP integrals where present. It also uses its Cholesky factors to build inactive-electron Coulomb/exchange fields and the inactive constant. Consequently, finite Cholesky error can affect the active one-body term and core energy as well as the active two-body tensor. This differs from factoring an already constructed Hamiltonian, which can preserve its existing one-body/constant terms unchanged.

The Hamiltonian tests cover restricted/unrestricted Cholesky construction and verify the correct effective nuclear constant for ECP-containing AgH. These are source test cases, not calculations run here.

A factor getter and an integral getter can have very different costs

A CholeskyHamiltonianContainer exposes get_three_center_integrals() for the alpha/beta factor matrices and get_ao_cholesky_vectors() for the optional AO data. These are the useful interfaces when a consumer can work directly with factors.

Its shared get_two_body_integrals() interface reconstructs and caches full four-index arrays on first access. More subtly, the current Cholesky get_two_body_element also builds that full cache before returning one entry. A request that looks like one scalar lookup can therefore trigger quartic memory. The container implementation makes this explicit.

Low-rank storage uses roughly $K^2R$ entries, advantageous when $R$ is substantially smaller than $K^2$. It can be less attractive for a tiny active tensor that inherits many AO factors. Once a legacy consumer requests the reconstructed tensor, peak memory may include both the factors and $K^4$ cache. Native CI methods discussed earlier call dense integral getters, so changing the container alone does not guarantee their full calculation stays low rank.

Other consumers have specialized paths. The current native qubit mapper forwards Cholesky factors to a C++ provider that contracts one pair row at a time, avoiding a dense four-index intermediate. Its sparse path likewise iterates stored entries. This is a concrete algorithmic advantage of a matched consumer, not a universal benefit of the class name. The equivalence tests compare those paths term by term for the same underlying integrals.

Hamiltonian factorization is a separate registered algorithm family

HamiltonianFactorization takes an existing Hamiltonian and returns a new one. The current factory registers double_factorization as its default concrete implementation. This is a working registered algorithm, unlike the interface-only effective-Hamiltonian constructor in Lesson 11.

The first stage produces Cholesky columns of the active pair tensor. The second stage treats each column as a symmetric $K$-by-$K$ matrix and diagonalizes it:

The rows of $U^t$ in this storage convention are the fragment eigenvectors. Each $L^t$ can have positive and negative eigenvalues even though the outer pair matrix $G$ is PSD. The Cholesky column already contains its magnitude, so the second-stage eigenvalues are used unscaled; multiplying by an additional square root of an imagined first-stage eigenvalue would double count scaling.

For the toy diagonal fragment above, the second stage simply gives eigenvalues one and two with the identity rotation. More generally it expresses each two-electron fragment using rotated one-electron directions, useful for later quantum constructions. It does not optimize one common new molecular-orbital basis for the whole Hamiltonian; different fragments have different rotations.

The exact current truncation semantics

The double-factorization code requires a non-null restricted Hamiltonian with a two-body term. For a dense input it works in the symmetric pair basis of dimension $K(K+1)/2$, averages the eight real chemist permutations, and performs pivoted Cholesky. It still receives an existing dense $K^4$ tensor and allocates a reduced-pair matrix whose size is quartic asymptotically. This route is not an on-the-fly AO integral algorithm.

truncation_threshold defaults to 1e-12 and applies to the largest remaining residual diagonal of that first-stage decomposition. The actual stopping level is at least its floating-point noise floor. It is not a cutoff applied independently to every original integral, every second-stage eigenvalue, or the final energy error. The current second stage diagonalizes all retained Cholesky matrices without an independent inner-rank truncation setting.

The implementation rejects NaN/non-finite tensor data and indefiniteness exposed through a sufficiently negative residual diagonal. A supplied tensor must genuinely have the PSD Coulomb structure needed by Cholesky. A generic Hermitian many-electron Hamiltonian, an antisymmetrized tensor, or a prospective downfolded interaction need not satisfy that pair-matrix condition. Small hidden negative directions can be truncated before a diagonal check exposes them. Averaging permutation symmetry also means an arbitrary nonsymmetric input is modified, not certified correct.

If the input is already a CholeskyHamiltonianContainer with the expected shape, the first stage is skipped. The stored factors are consumed directly, truncation_threshold is ignored, and each stored auxiliary column produces a second-stage fragment. Even redundant stored columns are retained. The tests demonstrate this by supplying five columns spanning only rank three and verifying that five fragments are returned, including with an enormous truncation threshold.

An asymmetric stored Cholesky matrix is symmetrized for diagonalization, with a warning for appreciable asymmetry. Discarding all dense-input fragments is an error rather than a valid zero-interaction output. The one-electron arrays, orbital object, core constant, inactive Fock matrix, and Hamiltonian type are otherwise carried over unchanged. Increasing a threshold cannot undo truncation already present in an input Cholesky representation.

The factorized container is more general than the default factorizer

FactorizedHamiltonianContainer stores U flattened as [R,B,K], W as [R,B,C], and WB as [R,C]. Here $R$ counts groups/fragments, $B$ counts basis rows per group, and $C$ counts copies per group; this $C$ is a dimension label, not the earlier MO coefficient matrix.

Define a leaf matrix

Its reconstructed tensor is

This is a plain sum of outer-product squares. WB does not enter the reconstructed two-electron tensor. The general container can hold a double-factorized tensor-hypercontraction-style representation with multiple copies and normalized basis rows; it does not require all those rows to be mutually orthogonal. Thus arbitrary W entries should not all be called eigenvalues of a fragment. They are eigenvalues in the particular orthogonal second-stage decomposition produced by DoubleFactorization.

That default producer sets $B=K$, one copy per fragment, W to the fragment eigenvalues, and WB=0. It does not implement a separate nonlinear fitting optimizer that searches all possible general U/W/WB tensor-hypercontraction representations. Having a storage format does not imply every fitting method for it is registered.

The constructor validation checks dimensions, finite factors, active-space consistency, and unit norm of each U basis row within 1e-8. Scaling a row requires compensating its W values, not silently changing the operator. The container is restricted-only and supports real factors.

get_u_matrices, get_w_matrices, get_wb_matrix, and the rank/basis/copy count getters expose the representation directly. get_two_body_integrals() reconstructs a quartic cache. Unlike the Cholesky scalar getter, its get_two_body_element() can contract just the requested pair entries from the factors when that cache is absent. reconstruct_two_body_integrals() explicitly requests a new full tensor. Those choices should be deliberate.

Adjusted one-body terms and normalization are not new chemistry energies

The factorized container also supplies get_h1_prime() and get_lambda() for downstream factor-based operator constructions. They are not replacements for the original stored one-electron integrals or total molecular energy.

Writing $h$ for the stored one-body matrix and $w_B^{tc}$ for WB, the implemented adjusted matrix is

The first correction includes the one-body term produced when products of one-electron operators are converted to properly normal-ordered two-electron interactions. The trace and identity-weight terms belong to the later decomposition convention. The fermionic algebra behind that distinction is developed in Lesson 14; simply treating a numerical square as an already normal-ordered two-electron operator would miss the correction.

The reported normalization is

This is a factorization-dependent algorithmic normalization, not the ground-state energy, not an entropy, and not necessarily the Pauli coefficient norm of a different mapped representation. Its implementation requires a symmetric adjusted one-body matrix and throws if that condition is not met. The separate core energy is not simply another term inserted into this formula.

Changing WB leaves the ordinary reconstructed two-electron tensor unchanged while changing the adjusted matrix and normalization. The factorized-container tests verify that distinction and the exact one-body correction formula. Do not tune representation weights and report a lower $\Lambda$ as a lower physical energy. These data will matter for sum-of-squares and quantum-walk constructions later, but their circuit and workload costs still need analysis.

Sparse storage means explicit missing entries, not inferred interactions

The SparseHamiltonianContainer primarily supports lattice/model Hamiltonians. It stores a sparse one-body matrix and a map from four-index keys to two-body values. A missing two-body entry is interpreted as zero; the map is not an instruction to fill all orbital permutations automatically. The common restricted constructors alias the same interaction across appropriate spin channels.

Converting a dense two-body array to this map keeps entries whose value is not exactly zero. It is not a user-tunable molecular integral screening algorithm. The dense-to-sparse one-body path uses Eigen’s sparse conversion. If one deliberately drops small nonzero interactions to obtain sparsity, that is an additional approximation whose symmetry and error must be checked. Duplicate operator contributions need to be combined consistently; a map of unique keys is not a substitute for that algebra.

Sparse-specific accessors allow direct iteration over stored data. The base facade still retains dense one-body information, so this implementation is not sparse-only in every array. A full two-body getter lazily allocates a zero-filled $K^4$ tensor and inserts the stored entries; even a one-body-only model can incur that allocation if a consumer asks for the dense two-body tensor.

has_two_body_integrals() checks whether a sparse-map object is present, not whether a materialized tensor has nonzero norm. An absent map and an explicitly supplied empty map can therefore have different availability metadata while representing zero interaction. The general rule is to inspect the concrete contract rather than infer cost or nonzero physics from one boolean.

Sparse storage is attractive when locality or a model definition makes most entries genuinely zero. Low rank is attractive when many dense entries can be expressed through common factors. A tensor can be dense but low rank, sparse but not very low rank, both, or neither. This is why no single container is always the most compact.

A source-backed representation workflow

Using the consistently selected orbital object from Lesson 10, the following describes three actual native requests without executing them:

from qdk_chemistry.algorithms import create

selected_orbitals = selected.get_orbitals()
canonical = create("hamiltonian_constructor", "qdk").run(selected_orbitals)
cholesky = create(
    "hamiltonian_constructor", "qdk_cholesky",
    cholesky_tolerance=1e-8, eri_threshold=1e-12,
    store_ao_cholesky_vectors=False,
).run(selected_orbitals)
factorized = create(
    "hamiltonian_factorization", "double_factorization",
    truncation_threshold=1e-12,
).run(cholesky)
factors = factorized.get_container()
rank = factors.get_num_ranks()

The last threshold does not re-truncate the supplied Cholesky factors. To study the dense-input cutoff, give the factorizer canonical instead. Neither example chooses a different active space. The selected carrier is used only to supply orbitals, not falsely presented as the correlated CI state whose amplitudes were preserved.

For a small validation problem, compare reconstructed tensors and core/one-body data, or evaluate the same normalized reference against both operators with consistent RDMs. Re-solving CI at matching active counts tests an eigenvalue consequence but may materialize dense tensors and incur substantial cost. The Python double-factorization tests use explicitly PSD synthetic tensors to check rank and reconstruction; they are not measured nitrogen compression ratios.

Serialization preserves the defining representation through container_type: the loader dispatch knows canonical, Cholesky, sparse, and factorized containers and rejects unknown types. Factors, orbital metadata, constants, and versions belong in that record; a lazy dense cache is not the scientific definition. No serialized files were written during preparation.

Budget factorization error separately from chemistry error

For a fixed normalized state, the energy change caused by integral approximation is a contraction of its density information with changes in the one-body, two-body, and constant terms. Schematically,

The number and magnitude of these terms matter; a pivot threshold is not automatically an absolute molecular energy error. For Hermitian operators in the same state space, an operator-norm bound on $\delta H$ bounds eigenvalue shifts, but this implementation does not turn every requested tolerance into that certificate.

A useful assessment varies factorization/screening thresholds while keeping geometry, basis, active partition, spin restriction, and solver precision fixed. Separate reconstruction error from error introduced when one constructs the inactive field with approximate factors. Preserve signs, pair ordering and constants before blaming a discrepancy on rank truncation.

Finally, factorization can lower memory or provide more useful input to a quantum algorithm without proving quantum advantage. A smaller rank may alter a later normalization and circuit decomposition, but the complete cost also includes fragment rotations, state preparation, precision, shots/repetitions, and error correction. Before analyzing those costs, we need a concrete vocabulary for the quantum registers and operations that will consume these representations.