Skip to content

Preconditioned LOBPCG solver for large non-periodic problems (v0.4.0) - #8

Merged
galjos merged 4 commits into
mainfrom
feature/lobpcg-solver
Aug 11, 2026
Merged

Preconditioned LOBPCG solver for large non-periodic problems (v0.4.0)#8
galjos merged 4 commits into
mainfrom
feature/lobpcg-solver

Conversation

@galjos

@galjos galjos commented Jul 18, 2026

Copy link
Copy Markdown
Member

Adds solver = lobpcg (input files and API, opt-in — arpack stays default): LOBPCG with a tensor-product kinetic preconditioner built from per-dimension 1D stencil operators, applied exactly via small per-dimension eigenbasis transforms — no full-operator factorization, so no 3D fill-in.

Solver time (3D harmonic oscillator)

grid arpack lobpcg speedup agreement
35³ 16.2 s 4.2 s 3.8× 1e-12
51³ 133.8 s 13.8 s 9.7× 4e-12

Accuracy safeguards

  • Every result's residual is checked against the real Hamiltonian; near-zero eigenvectors count as infinitely loose (catches silent LOBPCG collapse).
  • Escalation ladder: breakdown → retry → Arpack fallback; verification failure → Arpack re-solve. Worst case matches current accuracy/speed.
  • Periodic (complex Hermitian) problems are rejected with a clear error.
  • The preconditioner is exact for 1D/2D, an approximation for 3D — doesn't affect correctness (residual check is solver/stencil-agnostic), only convergence speed. Docstring and tests reflect this.

Adds IterativeSolvers to deps. Version 0.4.0.

New opt-in solver = lobpcg (input files and API): LOBPCG with a
tensor-product kinetic preconditioner - the kinetic operator is
approximated by the Kronecker sum of the per-dimension 1D stencil
operators, whose shifted inverse is applied exactly through small
per-dimension eigenbasis transforms, so no factorization of the full
operator and no 3D fill-in.

Measured on the 3D harmonic oscillator (solver time): 9.7x faster
than Arpack shift-invert at 51^3 (13.8 s vs 133.8 s), 3.8x at 35^3,
with eigenvalue agreement at 1e-12; end-to-end 3.7x at 51^3.

Accuracy safeguards, since correctness outranks speed:
- every iterative solve is verified after the fact: the kept
  eigenpairs' relative residuals are computed against the actual
  Hamiltonian, near-zero eigenvectors count as infinitely loose (a
  silently collapsed block can never pass), and lobpcg results that
  fail verification are re-solved with Arpack shift-invert
- lobpcg breakdowns (its internal factorizations can fail) retry
  with a fresh random block, then fall back to Arpack - never worse
  than the previous solver in either accuracy or robustness
- complex Hermitian (periodic) problems are rejected with a clear
  error; arpack remains the default everywhere

Adds IterativeSolvers as a dependency; version 0.4.0.
@galjos
galjos requested a review from 97gamjak July 18, 2026 13:37
@codecov

codecov Bot commented Jul 18, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 99.12281% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 97.20%. Comparing base (a3867e1) to head (f3974a9).
⚠️ Report is 5 commits behind head on main.

Files with missing lines Patch % Lines
src/solve.jl 97.77% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main       #8      +/-   ##
==========================================
+ Coverage   97.02%   97.20%   +0.18%     
==========================================
  Files          32       33       +1     
  Lines        1645     1752     +107     
==========================================
+ Hits         1596     1703     +107     
  Misses         49       49              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@github-actions

github-actions Bot commented Jul 18, 2026

Copy link
Copy Markdown

Benchmark Results (Julia v1)

Time benchmarks
main f3974a9... main / f3974a9...
assembly/laplacian_3D_15 24.7 ± 2.5 ms 24.7 ± 3.2 ms 1 ± 0.16
bandstructure/1D_kronigpenney_10k 0.0751 ± 0.013 s 0.0776 ± 0.014 s 0.968 ± 0.24
solve/1D_harmonic_201 4.21 ± 0.11 ms 4.21 ± 0.14 ms 1 ± 0.043
solve/2D_water 0.0939 ± 0.0027 s 0.0988 ± 0.0033 s 0.95 ± 0.042
solve/3D_harmonic_15 0.444 ± 0.15 s 0.438 ± 0.17 s 1.01 ± 0.53
time_to_load 1.47 ± 0.0071 s 1.78 ± 0.0051 s 0.826 ± 0.0046
Memory benchmarks
main f3974a9... main / f3974a9...
assembly/laplacian_3D_15 28.5 k allocs: 0.0561 GB 28.5 k allocs: 0.0561 GB 1
bandstructure/1D_kronigpenney_10k 0.776 M allocs: 0.0818 GB 0.776 M allocs: 0.0825 GB 0.993
solve/1D_harmonic_201 0.0376 M allocs: 4.49 MB 0.0376 M allocs: 4.52 MB 0.993
solve/2D_water 0.425 M allocs: 0.0926 GB 0.426 M allocs: 0.0929 GB 0.997
solve/3D_harmonic_15 0.544 M allocs: 0.237 GB 0.544 M allocs: 0.238 GB 0.999
time_to_load 0.149 k allocs: 11.1 kB 0.145 k allocs: 11 kB 1.02

…r accuracy claim

Codecov flagged src/preconditioner.jl and src/solve.jl at ~87% patch
coverage on PR #8: the uncovered lines were the lobpcg retry/fallback
and residual-escalation branches, which had only been exercised by
chance (RNG-dependent) rather than deliberately, so coverage varied
by platform/session.

- Expose lobpcg_maxiter and krylov_maxiter as tunable keywords on
  solveWrapper and use them to deterministically force the three
  failure modes these safety nets exist for, using genuine solver
  behavior rather than mocks:
    * IterativeSolvers.lobpcg refuses to run (throws) when the
      matrix is smaller than 3x the block size - reliably fails
      both retry attempts, forcing the arpack fallback
    * lobpcg_maxiter = 1 produces an under-converged (not thrown)
      result - caught by the residual verifier, re-solved with arpack
    * krylov_maxiter = 1 similarly under-converges KrylovKit's
      solver, exercising the generic (non-escalating) residual
      warning for solvers with no further fallback
- Add test_KineticPreconditioner, verifying the preconditioner
  against an independently-derived dense reference (not calling any
  of its own internals) across 1D/2D/3D and a periodic dimension,
  through all four ldiv! dispatches

While building that reference, discovered that the preconditioner's
Kronecker-sum kinetic operator is only an EXACT match to the true
production operator for 1D and 2D (buildLaplace_2d's Laplacian is
separable - dividing by 2^(dimension-1) recovers the Kronecker sum
bit-for-bit); buildLaplace_3d uses a more elaborate, non-separable
stencil, so for 3D the preconditioner is a heuristic approximation.
This does not affect correctness - every result is independently
residual-verified regardless of preconditioner quality - but the
docstring previously implied a close approximation to the true
operator in all cases, which was wrong for 3D. Corrected the wording
and added a regression assertion for the 1D/2D exactness.

Local coverage of both flagged files is now 0 uncovered lines.
…input validation gaps

- solveWrapper: re-check the residual after the lobpcg->arpack escalation
  solve, instead of returning it unconditionally - the rescue itself can
  land above tolerance for ill-conditioned Hamiltonians and was previously
  returned without any warning.
- setupSystem: reject solver=lobpcg for reciprocal (periodic k-point) runs
  as soon as reciprocal is known, before any output file is touched -
  previously this only threw deep in solve(), after main.jl had already
  deleted a pre-existing eigenvalues.dat.
- api.jl: extend the n_eigenvalues+5 < grid size guard to :lobpcg and
  :krylov, not just :arpack - lobpcg transparently falls back to arpack and
  was raising an opaque BoundsError instead of a clear ArgumentError.
- preconditioner.jl: correct the docstring's claim that the preconditioner
  is exact for all 2D problems - it is only exact for the 5-point stencil,
  not the package's default (9); broaden the regression test accordingly.
- test_3Dsmoke: seed the RNG and assert no fallback/escalation warning
  fires in the lobpcg-vs-arpack degenerate-cluster check, so the test
  actually proves lobpcg converged rather than merely that the pipeline's
  answer is correct (measured a 1/20-7/20 unseeded fallback rate that the
  old assertions couldn't distinguish from genuine convergence).
@galjos

galjos commented Jul 21, 2026

Copy link
Copy Markdown
Member Author

Pushed a follow-up commit after a deeper self-review pass:

  • The lobpcg→arpack escalation (triggered by a loose residual, not an exception) wasn't re-verified before being returned - fixed, with a deterministic regression test reproducing a case where the rescue itself lands above tolerance.
  • solver=lobpcg + periodic k-points now gets rejected in setupSystem, before any output file is touched - previously the ArgumentError only fired deep in solve(), after main.jl had already deleted a pre-existing eigenvalues.dat.
  • n_eigenvalues+5 < grid size is now validated for :lobpcg and :krylov too, not just :arpack - lobpcg was raising an opaque BoundsError for the same misconfiguration arpack already rejects cleanly.
  • Corrected the preconditioner docstring: it's only exact for 2D's 5-point stencil, not 2D in general (the default stencil is 9) - broadened the regression test to cover that case.
  • test_3Dsmoke's lobpcg-vs-arpack degenerate-cluster check now seeds the RNG and asserts no fallback/escalation warning fires, so it actually proves lobpcg converged (unseeded, it silently accepted an arpack rescue 1-in-20 to 7-in-20 runs depending on n_eigenvalues, since only the final energies were checked).

Full suite (1000 tests) passes.

Heads-up for whoever merges this and #9: both branches independently add Random = "1.10" to [compat] in Project.toml at different insertion points, so a plain merge of both is conflict-free but leaves a duplicate TOML key that breaks Pkg.instantiate(). One-line fix (delete the duplicate) after merging both.

@galjos
galjos merged commit b8a89d8 into main Aug 11, 2026
7 checks passed
@galjos
galjos deleted the feature/lobpcg-solver branch August 11, 2026 08:20
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant