perf: DMRG(2) with smarter adaptive tolerances#455
Conversation
8e22df4 to
8ab7f54
Compare
|
Here I've added some results on the comparisons for the sweeptimes, but also on the energy reached after a given walltime, in an attempt to not just measure what the fastest way is to get 10 sweeps done. I'd also say that it seems really hard to make generalizations about the actual walltimes it takes to get to a given energy, since the faster sweeps get compensated with the necessity for more sweeps. In general, this is not faster everywhere, but I think it at least feels faster because the sweeps take a little less time, and it makes a naive benchmark no longer favor the ITensor version that much. 
|
borisdevos
left a comment
borisdevos
left a comment
There was a problem hiding this comment.
This is amazing work, Lukas! The benchmarks also look very promising.
I was just wondering whether we should be more verbose about how AdaptiveKrylov is now the default everywhere. Also, this is not being compared explicitly in our tests by name, though I'm not sure what the criteria would be.
The refactoring indeed makes the code much easier to follow! The title of the PR could also added "performance" arguably ;)
You should definitely let others review this PR as well, since I'm no expert really, but hopefully I brought up some things with substance.
| function instantiate_algorithm( | ||
| alg::DynamicTol, decay_rate::Real, g_local::Real, g_global::Real, eps_trunc::Real | ||
| ) | ||
| tol = clamp(g_global * alg.tol_factor, alg.tol_min, alg.tol_max) | ||
| return _updatetol(alg.alg, tol) | ||
| end |
There was a problem hiding this comment.
This behavior is very different to what the user would previously expect from DynamicTol, no? I think someone can write DMRG(; alg_eigsolve = DynamicTol(Lanczos())), and then expect line 65's formula to still hold (using 1/sqrt(iter)).
There was a problem hiding this comment.
Good point, I somehow forgot about that since my intention was to completely supersede the DynamicTol implementation. I'll make sure to restore this though, thanks!
| HAC2 = normalize!(Heff * ac2) | ||
| AC2′ = copy(HAC2) | ||
| project_complement!(AC2′, ψ.AL[pos]) | ||
| project_complement_right!(AC2′, _transpose_tail(ψ.AR[pos + 1])) | ||
| ϵ_local = norm(AC2′) | ||
|
|
||
| # 1. local two-site update | ||
| alg_eigsolve = instantiate_algorithm(alg.alg_eigsolve, decay_rate, ϵ_local, ϵ_global, ϵ_trunc) | ||
| newA2center, info = @timeit timeroutput "AC2_eigsolve" begin | ||
| _, newA2center, info = fixedpoint(Heff, HAC2, :SR, alg_eigsolve) | ||
| (newA2center, info) | ||
| end |
There was a problem hiding this comment.
Since H * AC2 is used as an initial vector to the fixed point equation, this one application should also count to the number of matrix-vector multiplications, but this isn't being logged. Maybe this 1 matvec doesn't matter in the grand scheme of things though 🤷
There was a problem hiding this comment.
It definitely does, we are now doing on the order of 5-10 of them so it is 10-20%. However, the debug information is mostly for inspecting the eigsolve behavior, since that is also what is configured with the adaptive tolerances (krylovdim=3 would end up doing 4 matvecs if we count the initial one), so I think debug-wise this is probably fine?
There was a problem hiding this comment.
Is it useful to have H * AC2 instead of just AC2 as a starting vector? For transfer matrix fixed points, this should be a strictly better initial guess, but I am a bit worried for Hamiltonians which are such that the ground state energy is (close to) zero, and so H * AC2 might actually project out the most relevant contribution to the eigenvector when you are close to convergence.
There was a problem hiding this comment.
Oh you are right, I was definitely thinking of the transfer matrix case.
|
Thanks a lot @borisdevos for taking the time to go over this, the feedback is definitely appreciated! Some more things I would love inputs on:
|
| krylovdim = iszero(ρ) ? alg.krylovdim_min : alg.krylovdim_max | ||
| maxiter = iszero(ρ) ? alg.iter_min : alg.iter_max | ||
| return (T ? Lanczos : Arnoldi)(; alg.orth, krylovdim, maxiter, tol, alg.eager, alg.verbosity) | ||
| end |
There was a problem hiding this comment.
With decay_rate = clamp((first(info.normres) / ϵ_local)^(1 / max(1, info.numops)), 1.0e-3, 0.999) above, can ρ ever fall outside (or even reach) the limits 0.0 or 1.0?
There was a problem hiding this comment.
yes, the very first iteration I don't have any history so I set it to 0 to detect this. I'll add a comment here to clarify!
lkdvos
changed the title
Switch the DMRG/DMRG2 stop test from the overlap-based convergence error `ϵ_conv = 1 - |⟨old|new⟩|` to the Galerkin error (already computed each local update as `ϵ_local`/`ϵ_global`), matching VUMPS/iDMRG. Drop the now-unused `ϵ_conv` bookkeeping from both `local_update!` methods and `find_groundstate!`. Perturb the default `FiniteExcited` initial guess with a small per-site random component. Under the Galerkin criterion an unperturbed ground-state guess can already be (near) an eigenvector of the shifted operator, so convergence would trigger on the first sweep and return the ground state instead of climbing to the excited state. The perturbation preserves spaces (symmetry sector) and bond dimensions. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
DMRG2 reused `HAC2 = Heff * ac2` (needed for the Galerkin `ϵ_local`) as the eigensolver seed. Seeding with H*AC is power-iteration intuition (it amplifies the largest-magnitude eigenvector), which is wrong for `:SR` (smallest algebraic): it de-weights the ground state, and when the effective ground eigenvalue is 0 it annihilates the ground-state component exactly. Since H preserves the ground state's orthogonal complement, the Krylov solve then converges to the smallest nonzero eigenvalue instead — the wrong state. This bites frustration-free/parent Hamiltonians and any H - E0 shift. Seed with the plain two-site center `ac2` instead, consistent with single-site DMRG's `ac_old`. `HAC2` is still computed for the Galerkin gradient only. Verified on a finite TFIM (L=16, g=1) shifted so E0=0: old seeding gives a non-convergent 4.73 (ϵ_conv≈0.88), fixed seeding gives 5e-9≈0 (ϵ_conv≈2e-15); unshifted DMRG2 still matches single-site DMRG to 1e-14. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…port different metrics and ignore ones that aren't present
`gauge2!` unconditionally promoted the SVD center matrix to complex (`complex(c)`), which for a real-valued state forced a complexify then coerce-back-to-real roundtrip when installing the tensors. Only complexify when the state's scalartype is itself complex. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
3 participants
This PR features a somewhat large refactor of the DMRG implementations, and finally fixes #311 hopefully in a more systematic way.
The main idea is that I want to use loose bounds on the krylov solvers, since currently everyone has the wrong assumption that ITensors is a lot faster than MPSKit, which is just due to the default settings.
Nevertheless, while for ITensorMPS you just have to know and read the docs to figure out that if you want to actually converge something for a non-gapped state, I wanted to keep some guarantees on stability and convergence rates.
Therefore, this PR introduces a bunch more knobs to determine adaptive tolerance and krylov settings based on the information that is available:
A secondary part of this PR is just me refactoring the DMRG implementations into a more readable and consolidated form, which I think overall just improves the code quailty.
Will still attempt to post some actual benchmarks here comparing main/this branch, as well as some comparisons with ITensorMPS.