DMRG3S: strictly single-site DMRG subspace expansion (Hubig et al. 2015)

[VinceNeede]

Summary

Implements DMRG3S — strictly single-site DMRG with subspace expansion — from Hubig, McCulloch, Schollwöck & Wolf, Phys. Rev. B 91, 155115 (2015). Unlike CBE (alg_expand), which enriches the bond before the local eigensolve, DMRG3S enriches it after, using a perturbation built from the just-optimized tensor, the local Hamiltonian, and the environment. This lets single-site DMRG escape local minima caused by missing quantum-number sectors that CBE-style pre-expansion can't reach.

This PR also targets TDVP, pending design agreement

DMRG only for now. TDVP support belongs in this same PR, but I'd like design agreement on the interface below first, since it directly shapes how it plugs into TDVP too. Docs are also on hold for the same reason — no point documenting an interface that might still change.

What's here

• Refactored find_groundstate!'s left-to-right and right-to-left sweeps from two separate, near-duplicate loops into a single loop over direction in (:right, :left), dispatching on Val. No behavioral change; this just made room for the new expansion step without doubling its code.
• alg_post_expand field on DMRG, dispatching a new post_expand! step after the eigensolve, alongside the existing pre-eigensolve alg_expand. Default (NoExpand()) reproduces current behavior exactly.
• DMRG3S, with a composable NoiseSchedule for the mixing factor (ExponentialDecay, Warmup, ∘ for combining them).
• alg_expand and alg_post_expand are independent and can be used together or separately.
• A warning fires if either expansion mechanism is active without a truncating alg_gauge, since the bond would otherwise grow unboundedly each sweep.
• Tests: a bond-growth smoke test alongside the existing CBE tests, and a reproduction of the paper's Sec. VII A example (plain single-site DMRG gets stuck at E ≈ -6.35, DMRG3S recovers E ≈ -8.68).

Open design questions for maintainers

1. alg_post_expand sentinel: I used a dedicated NoExpand() type rather than nothing (which alg_expand uses). This lets NoExpand carry its own post_expand! method for dispatch, but it's an inconsistency with alg_expand's convention worth flagging.
2. Expand+gauge+neighbor-pad are fused into one post_expand! call rather than kept separate, because splitting them left the MPS briefly inconsistent between steps (caught via a real bug during development — happy to share details). This does mean post_expand! needs alg_gauge threaded in as a parameter, coupling it to the gauge step more tightly than alg_expand is.

TODO before merge

• Docs
• TDVP support (Decided to add it in another PR)
• Decide on Q1–Q2 above

[github-actions]

Your PR no longer requires formatting changes. Thank you for your contribution!

[codecov]

Codecov Report

❌ Patch coverage is 87.50000% with 10 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/algorithms/post_expand/dmrg3s.jl	88.67%	6 Missing ⚠️
src/algorithms/groundstate/dmrg.jl	78.94%	4 Missing ⚠️

Files with missing lines	Coverage Δ
src/MPSKit.jl	100.00% <ø> (ø)
src/algorithms/approximate/fvomps.jl	91.30% <100.00%> (ø)
src/algorithms/post_expand/post_expand.jl	100.00% <100.00%> (ø)
src/utility/utility.jl	59.76% <100.00%> (ø)
src/algorithms/groundstate/dmrg.jl	92.78% <78.94%> (-2.46%)	⬇️
src/algorithms/post_expand/dmrg3s.jl	88.67% <88.67%> (ø)

🚀 New features to boost your workflow:

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

[lkdvos]

Thanks for this, I had a brief look and I think it already looks quite well.
I'll try and make some time for a proper review either tomorrow or the day after, but mostly I'm wondering if it might make sense to wait for #455, as that has a similar kind of refactor of the DMRG machinery to avoid duplication.

Considering the expansion algorithm itself - without having looked in a lot of detail so definitely feel free to tell me I'm wrong - I think it could be nice to avoid having both an alg_gauge and an alg_postexpand, and simply merge the two together. As always, naming things is really the hardest part, but it feels a bit clumsy to have the if-guard in the generic local step, and to have to pass 2 algorithms to the perturb! body.
Do you think it might be possible to simply put the alg_gauge inside of the DMRG3S, and then keep only one of these around that we dispatch on? I'm not sure how blasphemous it would be to call these both alg_gauge...
We might also try and discuss more carefully about the invalid state though, the FiniteMPS struct and the interaction with the environments is quite subtle, and it might still be that it is still fixable to separate the two cleanly as well.

For the TDVP implementation, I definitely agree to flush this one out first, and honestly it's probably easier to tackle and review that in a separate PR.

Considering your other question, I'm definitely happy to add the NoExpand, and since we are in the middle of a breaking cycle, I would be happy to change the alg_expand to use that too.

[VinceNeede]

Hi @lkdvos,

Pushed a sketch of merging alg_gauge/alg_post_expand along the lines you suggested. The dispatch side is in place: NoExpand/DMRG3S both wrap an inner gauge algorithm, gauge!/gauge2! dispatch on whichever one alg_gauge is, and local_update! ends up with a single unconditional gauge!(...) call — no if-guard, no second algorithm argument. The noise schedule advances once per outer iteration via _update_alg_gauge, threaded through find_groundstate!.

I haven't wired up DMRG's constructor yet, so this isn't runnable end-to-end through the public API — that's next, and here's the direction I'm planning to take it.

Since DMRG's trscheme/alg_svd/alg_orth also determine the inner gauge, a user could in principle construct DMRG3S with its own gauge already set (DMRG3S(noise, schedule, some_gauge)) and pass a different trscheme to DMRG, giving two conflicting sources for the same thing. My plan is to make trscheme/alg_svd/alg_orth keyword arguments of NoExpand/DMRG3S themselves, building the inner gauge there — rather than have them accessible from DMRG's constructor directly. That way each gauge object is self-contained, and the conflict can't arise in the first place. Let me know if you'd rather keep trscheme etc. as DMRG-level kwargs instead (with some other rule for resolving the conflict), otherwise I'll go ahead and wire it up this way.

[lkdvos]

Thanks for rebasing onto that other work!

Considering the DMRG constructor, I think I agree that we should exclude being able to set conflicting schemes, but unless I'm missing something that would already not be possible - the gauging algorithm is either a (truncated) orth algorithm or DMRG3S, and in the latter case there is only a single (truncated) orth algorithm in the DMRG3S.
If at all possible, I would like to design the DMRG constructor to be able to take the trscheme still, and have that routed through automatically, or error if both a trscheme and DMRG3S are provided. My experience is that while these "algorithm structs" are very convenient as a developer, it is somewhat harder as a newcomer to the package to keep track of everything, and a simple keyword-argument approach, which possibly does not expose all settings but handles the common cases, helps alleviate some of these issues.

[lkdvos]

I think this contraction might be the same as our AC_projection contractions, and it could be beneficial to reuse that - the main reason is that we also support density-matrix MPS, which have 2 physical legs and would break this contraction, and reusing that function automatically fixes that.

[VinceNeede]

I already took a look at AC_projection, and from my understanding, it requires to supply all the enviroments, both left and right, here instead only one is required at a time. I'm not sure it can be adapted to this use case

[lkdvos]

yes, but the combiner can be thought of as an environment in this case, so while the code routing is a bit more complex, the contraction is the same as this one:

MPSKit.jl/src/algorithms/derivatives/mpo_derivatives.jl

Lines 69 to 70 in 5bbab08

	@plansor y[-1 -2; -3] ≔ h.leftenv[-1 5; 4] * x[4 2; 1] * h.operators[1][5 -2; 2 3] *
	h.rightenv[1 3; -3]

[VinceNeede]

I have tried using MPO_AC_Hamiltonian, but the assumption of this line breaks:

MPSKit.jl/src/algorithms/derivatives/mpo_derivatives.jl

Line 71 in 5bbab08

return y isa AbstractBlockTensorMap ? only(y) : y

For a proper environment, the resulting MPS does have a single block, but in the case of a combiner this breaks. For an hamiltonian of bond dimension 3, the right environment has blocks $1 \times 3 \times 1$ while the combiner has $1\times 3 \times 3$, which give y as a $1\times1\times3$ blocked.
I also tried JordanMPO_AC_Hamiltonian, but it requires indexing the right environment

MPSKit.jl/src/algorithms/derivatives/hamiltonian_derivatives.jl

Line 83 in 5bbab08

A = MPO_AC_Hamiltonian(GL[2:(end - 1)], WA, GR[2:(end - 1)])

And this does not work for combiner, because getindex only works witth zero or 1 non trivial dimensions https://github.com/QuantumKitHub/BlockTensorKit.jl/blob/215f9283301cc4d964d5819754320d6c4a25a064/src/vectorspaces/sumspaceindices.jl#L63-L76

I don't think I'm in a good position to safely change either of these, so I'd rather get your take on how you'd want this handled before I try anything myself. Is there a way you'd recommend getting the raw, uncollapsed contraction out of this framework, or is this a case where the current abstraction just isn't meant to support it and I should go another route entirely?

[lkdvos]

Thanks for checking and explaining why it fails, I definitely also missed that.
Looking into this in a little more detail, I guess we probably don't want to get a
$1 \times 1 \times 3$ object out anyways, since you want to concatenate it with the MPSTensor right afterwards.
In particular, it should be more efficient to immediately make the fuser that collapses these levels together. I would have to look at the exact implementation requirements, but I think you can call TensorKit.oplus on the fused spaces to combine them. (Unfortunately, that indeed now means we can no longer reuse the fuser function, but that's fine).

Let me know if this makes sense, I definitely might be missing something, and I'm happy to explain more of what I'm thinking of.
The idea would definitely be to try and reuse the MPO_AC_Hamiltonian, not the JordanMPO one.

[VinceNeede]

Hi @lkdvos,

I was trying to use oplus to flatten a leg, but ran into a loss of precision in doing so (my second DMRG3S test, which checks proper subspace exploration, started failing). Tracked it down to unitary/isomorphism producing a non-unitary map when the domain is a SumSpace and the codomain is its flattened (oplus(...)) counterpart — even when every sector's dimension matches exactly between the two. Minimal example:

using TensorKit, BlockTensorKit, LinearAlgebra

V1 = Rep[U₁](0 => 2, 1 => 1)
V2 = Rep[U₁](0 => 2, 1 => 1)
Vsum = V1 ⊞ V2
Vflat = oplus(Vsum)

u = unitary(ComplexF64, Vflat, Vsum)
@show norm(u' * u - id(Vsum))          # 1.73..., should be ~0
@show norm(only(u * u') - id(Vflat))   # 1.73..., should be ~0

Is this a known limitation, or should I open a proper issue (and if so, here or on TensorKit)? Happy to give more detail / try a fix there.

Thanks!

[lkdvos]

This is definitely not a known limitation and almost certainly is a bug. Feel free to open an issue on BlockTensorKit.jl, where I assume this is originating from, and I'll try and have a look asap.

[VinceNeede]

Update

• DMRG constructor has been updated to accept alg_gauge. The inner algorithm used in the svd/ortho step is constructed as before, ensuring back compatibility. trscheme is initially defaulted to nothing to prevent supplying the svd/ortho algorithm twice, it is then repristinated to the previous default as trscheme = something(trscheme, notrunc()). The docstring has been updated to explain the behavior and DMRG3S role.
• fuser type restriction has been removed, and is now used to construct the combiner. The @plansor contraction has been reordered to mimic MPO_AC_Hamiltonian, even though it still can't be used directly for the above mentioned problems.

TODO

• Move update_alg_gauge inside local_update!, verifying no significative change in performance/allocations is present
• Understand how MPO_AC_Hamiltonian can be used for the contraction step

[lkdvos]

Thanks for the update! I commented on the reuse of the AC_projection, hopefully that makes sense, and I think it is safe to assume that the alg_gauge computation per site is not hurting performance, indeed after these two this looks ready to go for me.
Already thanks for putting in this effort, this is a very nice PR.