refactor: DMRG(2) with smarter adaptive tolerances by lkdvos · Pull Request #455 · QuantumKitHub/MPSKit.jl

lkdvos · 2026-07-10T01:04:49Z

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:

The tolerance of the krylov solver depends on the global convergence, and oversolving locally (especially for finite systems) has very little benefit.
The tolerance of the krylov solver should not go far below the truncation error incurred in adaptive DMRG schemes, since that error will dominate the eigenvector precision anyways.
The gap of the local effective system can be predicted by investigating the decay rate of the previous sweep, which is theoretically predicted to be geometric for Lanczos solvers. Using this, we can make an educated guess for how gapped the system is and adapt accordingly.

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.

codecov · 2026-07-10T01:28:53Z

Codecov Report

❌ Patch coverage is 88.07947% with 18 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/utility/dynamictols.jl	75.86%	7 Missing ⚠️
src/states/orthoview.jl	82.85%	6 Missing ⚠️
src/algorithms/fixedpoint.jl	40.00%	3 Missing ⚠️
src/algorithms/groundstate/dmrg.jl	98.48%	1 Missing ⚠️
src/states/finitemps.jl	83.33%	1 Missing ⚠️

Files with missing lines	Coverage Δ
src/MPSKit.jl	`100.00% <ø> (ø)`
src/algorithms/approximate/fvomps.jl	`91.30% <100.00%> (ø)`
src/algorithms/changebonds/optimalexpand.jl	`98.73% <100.00%> (ø)`
src/algorithms/changebonds/randexpand.jl	`64.38% <100.00%> (ø)`
src/algorithms/changebonds/sketchedexpand.jl	`98.43% <100.00%> (ø)`
src/operators/jordanmpotensor.jl	`84.93% <100.00%> (+0.19%)`	⬆️
src/utility/defaults.jl	`100.00% <100.00%> (ø)`
src/algorithms/groundstate/dmrg.jl	`95.55% <98.48%> (+3.59%)`	⬆️
src/states/finitemps.jl	`93.75% <83.33%> (-0.26%)`	⬇️
src/algorithms/fixedpoint.jl	`70.00% <40.00%> (+0.76%)`	⬆️
... and 2 more

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

lkdvos · 2026-07-12T03:19:41Z

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.
What's interesting is that we now consistently outperform ITensorMPS, even without symmetries, but that probably we will have to investigate the single-site version a bit more, since that seems to not perform the way I'd expect.

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

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.

borisdevos · 2026-07-13T09:05:57Z

+                        ϵ_global, ϵ_truncs[pos], decay_rates[pos],
+                        timeroutput
+                    )
+                    ϵ_global = maximum(ϵ_locals)


What's the role of the local and global errors now? It seems that ϵ_conv is the only one that matters now

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

borisdevos · 2026-07-13T09:22:32Z


-# 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)


Same comment as below (oops) on the dispatch design choice.

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

function AC2(ψ::InfiniteMPS, i::Integer; kind = :ACAR)

if kind == :ACAR

return ψ.AC[i] * _transpose_tail(ψ.AR[i + 1])

elseif kind == :ALAC

return ψ.AL[i] * _transpose_tail(ψ.AC[i + 1])

else

throw(ArgumentError("Invalid kind: $kind"))

end

end

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.

borisdevos · 2026-07-13T09:29:13Z

+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


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)).

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!

borisdevos · 2026-07-13T09:34:38Z

+    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


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 🤷

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?

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.

Oh you are right, I was definitely thinking of the transfer matrix case.

lkdvos · 2026-07-13T13:28:20Z

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:

Looking at the benchmarks right now, it is actually not entirely clear that these looser defaults are really helping. Of course, if the metric is to get fast sweeps this makes them a lot faster, but if the metric is to get to the groundstate faster, this doesn't seem to really be the case, especially in the gapped case or for higher bond dimensions. I might redo an actual hyperparameter sweep in the current implementation and see if I can actually get some more quantitative data on this.
I really don't know what to do with the convergence criterion. The galerkin error seems more principled to me, but especially in the twosite context I'm not even sure if I'm defining this correctly, since the correct thing might be to project onto the twosite left nullspace? Maybe @Jutho can say something more about this. Ultimately I'm happy to in the future work towards making this configurable, but I don't want to tackle that in this PR.
I have a bunch more tricks up my sleeve to make this faster, which mostly have to do with the way the environments are constructed and stored, which arguably will make a much larger difference across the board than this hyperparameter tweaking, so I don't want to spend too much time dwelling on this. I'm probably going to try and add some of these benchmark results to the documentation, where for the comparison with ITensor I will just use the "fixed krylov" settings, which is also what TenPy settled on for their benchmarks.
The conv vs galerkin errors currently have a different scale - the former is 1 - dot(x, y) while the latter is norm = sqrt(dot(...)). However, if I naively square root the former the behavior doesn't really fix my excitations immediately converging issue...

Jutho · 2026-07-14T09:47:53Z

-                    @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′)))


Are the norms of ac_old and AC' necessary? Aren't they 1 automatically?

Jutho · 2026-07-14T10:57:20Z

+_num_updates(::DMRG2, ψ) = length(ψ) - 1
+
+_sweep_ranges(::DMRG, ψ) = (1:(length(ψ) - 1), length(ψ):-1:2)
+_sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1)


It doesn't really matter, but this is a bit inconsistent between DMRG and DMRG2, I would have expected:

Suggested change

_sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1)

_sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 2), (length(ψ) - 1):-1:2)

Also, for both DMRG and DMRG2, it is indeed true that the first (1) and last site (length(ψ)forDMRGandlength(ψ)-1forDMRG2`) only needs to be visited in either forward or backward, e.g. backwards only goes to site 2, and then the next forward starts again at site 1. However, in the very last iteration, you do want to also visit site 1 after the backward sweep. (In most cases this doesn't matter, since anyway if Dl * d = Dr, there is actually no point in optimizing that site, but there might be problems where d, the physical dimension is very large, and where this would have an effect).

While it is definitely possible to include site 1 in the last sweep, I don't really see how this makes any real change, except for where the tolerance is checked. If we do that 1 extra update, doesn't that then mean that now site 2 might need updating?
I definitely might be missing something, but somehow I always kind of assumed that this is completely arbitrary, and we could also just check convergence after every local update and stop at whatever site reaches that?

lkdvos force-pushed the ld-dynamictols branch from 5ebd191 to 8e22df4 Compare July 12, 2026 00:05

lkdvos added 5 commits July 11, 2026 20:06

refactor fixedpoint to avoid warnings

2e061fe

refactor DMRG implementation

7bc6056

refactor gauging

8a4cc2b

add AdaptiveKrylov solver

2a8ade7

refactor DMRG2 to be in line with DMRG

8ab7f54

lkdvos force-pushed the ld-dynamictols branch from 8e22df4 to 8ab7f54 Compare July 12, 2026 00:06

lkdvos marked this pull request as ready for review July 12, 2026 00:06

lkdvos added 3 commits July 11, 2026 23:21

small type stability improvements

a15f465

small fixes

36bd652

restore old DMRG convergence behavior

f426dd4

borisdevos reviewed Jul 13, 2026

View reviewed changes

Jutho reviewed Jul 14, 2026

View reviewed changes

Comment thread src/algorithms/groundstate/dmrg.jl

Jutho reviewed Jul 14, 2026

View reviewed changes

	function AC2(ψ::InfiniteMPS, i::Integer; kind = :ACAR)
	if kind == :ACAR
	return ψ.AC[i] * _transpose_tail(ψ.AR[i + 1])
	elseif kind == :ALAC
	return ψ.AL[i] * _transpose_tail(ψ.AC[i + 1])
	else
	throw(ArgumentError("Invalid kind: $kind"))
	end
	end

	_sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 1), (length(ψ) - 2):-1:1)
	_sweep_ranges(::DMRG2, ψ) = (1:(length(ψ) - 2), (length(ψ) - 1):-1:2)

Uh oh!

Conversation

lkdvos commented Jul 10, 2026

Uh oh!

codecov Bot commented Jul 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

lkdvos commented Jul 12, 2026

Uh oh!

borisdevos left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lkdvos commented Jul 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

codecov Bot commented Jul 10, 2026 •

edited

Loading

lkdvos commented Jul 13, 2026 •

edited

Loading