refactor: DMRG(2) with smarter adaptive tolerances#455
Conversation
|
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.
| ϵ_global, ϵ_truncs[pos], decay_rates[pos], | ||
| timeroutput | ||
| ) | ||
| ϵ_global = maximum(ϵ_locals) |
There was a problem hiding this comment.
What's the role of the local and global errors now? It seems that ϵ_conv is the only one that matters now
There was a problem hiding this comment.
The ϵ_global is used in the adaptive scheme, and signifies more or less the norm of the gradient projected onto the tangent space. On the other hand, the ϵ_conv gives how much the tensors have drifted.
I really wanted to avoid to use this latter quantity, but that was what DMRG2 was using before, and the problem is that the former can become zero even if the state is not converged, as it becomes zero for any eigenstate.
This broke the excitation tests, because FiniteExcited basically adds a penalty term to the ground state, which does not change the eigenvectors, and therefore the algorithm immediately terminates.
A different solution would be to just add a miniter as well, which I'm definitely also open to
There was a problem hiding this comment.
Does it make sense for ϵ_global = norm(ϵ_locals). I think the total projected norm is in fact the 2-norm of the local ones.
As for using ϵ_conv and the issue you describe, this sounds like a deeper problem. If ϵ_global = 0, then all local eigenvalue problems should also be satisfied by their initial guess, and should return immediately. And if they don't, as in the case you describe where you add a penalty term but still use the ground state as initial guess, that is really relying on some numerical noise to mix in new eigenvector components in the Krylov subspace, and it does not sound very robust. For finding excited states in this way, I think you really want to explicitly perturb the state away from the ground state as initial guess.
There was a problem hiding this comment.
I think what is subtle about this is that the local errors are computed in a non-simultaneous way, so really only the one that got computed last is fully correct with the current state, since we are updating as we go along the sweep. That's also why I'm using maximum here rather than taking the norm. A different reason that is more pragmatic is that I want something that is size-independent, as otherwise the various tolerances and settings would have to scale inversely with system size. This is less of a problem in the infinite settings, but here I think it makes more sense to have a 50-site and a 500-site have the same local tolerances, not a global one?
There was a problem hiding this comment.
Yes, in particular the size dependence is also something that occurred to me after writing this comment, so I think maximum is a suitable choice.
|
|
||
| # TODO: check where gauge center is to determine efficient kind | ||
| AC2(psi::FiniteMPS, site::Int) = psi.AC[site] * _transpose_tail(psi.AR[site + 1]) | ||
| function AC2(psi::FiniteMPS, site::Int; kind = :ACAR) |
There was a problem hiding this comment.
Same comment as below (oops) on the dispatch design choice.
There was a problem hiding this comment.
This is in line with the function that already existed for the infinite case:
MPSKit.jl/src/states/infinitemps.jl
Lines 237 to 245 in 4e0e891
The main difference is that the return type does not change for different values of this keyword argument, so the Val dispatch is not really necessary and we can just use runtime to distinguish between the two. The left_gauge and right_gauge functions on the other hand return their arguments in different orders, so to avoid type-instability I added the Val dispatch. In hindsight, I guess at the gauge! level this is not really relevant since it just returns the state, so I agree that it might be reasonable to change that to a keyword argument too, to take some strain from the compiler.
| 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:
|
| @timeit timeroutput "gauge" left_gauge!(ψ, pos, AC′, alg.alg_gauge; normalize = true) | ||
| end | ||
| # convergence error: how much the eigensolve moved the center, `1 - |⟨old|new⟩|`. | ||
| ϵ_conv = abs(1 - abs(dot(ac_old, AC′)) / (norm(ac_old) * norm(AC′))) |
There was a problem hiding this comment.
Are the norms of ac_old and AC' necessary? Aren't they 1 automatically?
| 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!
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.