Optimize v3 native KPM energy truncation, benchmark against untruncated KPM - #98
Conversation
…d benchmark Chain::kpm_local_krylov_projection() (mpscpp3/chain_session.h) built its dense per-site Krylov Hamiltonian Hk in two passes: a Gram-Schmidt extension loop that computes H*V_j via LocalMPO::product() while orthogonalizing against existing vectors, then a second, separate loop that called LocalMPO::product() again for every vector just to rebuild Hk from scratch. That second pass is redundant -- the orthogonalization coefficients the first loop already computes and discards (summed across its two reorthogonalization passes) are exactly Hk's off/diagonal entries to working precision, since a Hermitian operator's projected matrix satisfies Hk(j,i) = conj(Hk(i,j)). Only the last accepted vector's own row/column is genuinely never produced by the extension loop (which stops one short by design), so exactly one extra matvec covers it, not k. This cuts LocalMPO::product() calls -- the dominant cost of this method, a full local-effective-Hamiltonian contraction -- from ~2k-1 down to ~k per site, with no change to the numerics (same dense Hermitian Hk, same Hermitization/diagHermitian/thresholding). Verified against tests/test_kpm_energy_truncation_v3.py, tests/test_kpm_energy_truncation_v3_accuracy.py, and tests/test_kpm_divergence_guard_catchable.py (all still pass), plus a full pytest tests/ run (293 passed, 0 failed). Added examples/dynamical_correlator/kpm_energy_truncation_v3_benchmark, turning what was previously an ad hoc, unsaved timing observation into a saved, rerunnable regression. Measured on this system (4- and 6-site Heisenberg chains, dK=30/nsweeps=10, back-to-back before/after the optimization): n=4: truncated KPM 0.93s -> 0.56s (~1.7x faster) n=6: truncated KPM 13.9s -> 7.2s (~1.9x faster) Energy truncation remains a net slowdown vs. plain KPM on these small systems (untruncated: 0.06s / 0.14s) -- this optimization narrows that gap, it doesn't eliminate it; the resolution win truncation buys still doesn't pay for itself until a system is large enough for a correlator's spectral weight to genuinely separate from the full bandwidth (see the example's own module docstring). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01DoZKWzhMo7qKoiP12gbcvv
…_wide Found while running the full suite to verify an unrelated v3 KPM optimization: this test builds an excited state via _excited_state(sc), which consumes draws from Python's global np.random state through pyitensor's randomized excited-state search -- but unlike its two sibling tests in the same file, it never called np.random.seed(0) first. Whichever tests ran earlier in the same pytest process therefore determined its random draws, making it pass or fail depending on test collection order rather than on its own logic (reproduced deterministically twice in a row against the same commit: identical failure value 0.10780681457731987 both times, matching pytest's fixed, non-randomized default collection order). Added the same np.random.seed(0) call its sibling tests already use, in the same position (before any chain/DMRG construction). Verified with a full pytest tests/ run (293 passed, 0 failed) and a separate run of this file after a broad slice of other random-draw-consuming test files (190 passed, 0 failed). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01DoZKWzhMo7qKoiP12gbcvv
Follow-up: 20-site benchmark (larger, more extensive system)Tested whether the truncation cost/benefit crossover starts to favor truncation at a larger, more extensive system size, per the earlier note that this "may flip favorable only for genuinely large systems (not benchmarked)." Setup: 20-site S=1/2 Heisenberg chain, same
Timing (delta=0.08): untruncated 38.4s avg. Truncated: dK=30/ns=10 (paper's Table I) 422.0s — 11.0x slower, down from 48-91x slower at n=6 with the same dK/nsweeps. The slowdown is shrinking with system size as the paper's own regime assumption predicts, but there's still no crossover to a real speedup at n=20. Accuracy: truncated (dK≥20) lineshape correlates well with untruncated (Pearson r≈0.96). Peak locations look different at first glance (untruncated argmax=0.6 vs truncated argmax=0.4), but this isn't real disagreement — untruncated y(0.4)=0.324 vs y(0.6)=0.327 is a <1% tie; at this size the correlator is a broad finite-size continuum (level spacing ≈ delta's own resolution), not the isolated single peak the 4/6-site tests rely on, so a strict peak-match check is the wrong accuracy metric here. Separate, notable finding: the library's own hardcoded defaults ( 🤖 Generated with Claude Code |
Summary
LocalMPO::product()(matvec) calls inChain::kpm_local_krylov_projection()(mpscpp3/chain_session.h) from ~2k-1 down to ~k per site: the Gram-Schmidt extension loop already computes the dense projected HamiltonianHk's entries as a byproduct of orthogonalization (summed across its two reorthogonalization passes), but the code discarded them and recomputedHkfrom scratch via a second full matvec pass. No change to the numerics — same dense HermitianHk, same Hermitization/diagHermitian/thresholding.examples/dynamical_correlator/kpm_energy_truncation_v3_benchmark, turning a previously ad hoc, unsaved timing observation (from PR Add KPM energy truncation (pyitensor): narrower-than-bandwidth dynamical-correlator windows #96) into a saved, rerunnable regression comparing untruncated vs. truncated v3 KPM.test_energy_truncate_noop_when_window_is_wideintests/test_kpm_energy_truncation.py) found while verifying this change against the full suite: it was missing thenp.random.seed(0)call its sibling tests already have, so its outcome depended on pytest's test collection order rather than its own logic.Benchmark (this system, 4- and 6-site Heisenberg chains, dK=30/nsweeps=10, back-to-back before/after)
Energy truncation remains a net slowdown vs. plain KPM on these small systems — this optimization narrows that gap, it doesn't eliminate it. The resolution win truncation buys doesn't pay for itself until a system is large enough for a correlator's spectral weight to genuinely separate from the full bandwidth (see PR #96 and the new example's module docstring).
Test plan
tests/test_kpm_energy_truncation_v3.py,tests/test_kpm_energy_truncation_v3_accuracy.py,tests/test_kpm_divergence_guard_catchable.pyall passpytest tests/(v2+v3 compiled): 293 passed, 0 failed (run twice for determinism)🤖 Generated with Claude Code
https://claude.ai/code/session_01DoZKWzhMo7qKoiP12gbcvv