SU(2) DMRG performance unexpectedly slower than U(1) ITensor version

[drichg]

Hi, thank you for the great package!
I'm benchmarking DMRG calculations using MPSKit.jl with SU(2) symmetry, and I encountered a significant performance discrepancy compared to a U(1)-symmetric calculation done with ITensor.
Here are my benchmarks:
ITensors (Julia):
Symmetry: U(1)
Bond dimension: D = 4800
Time per sweep: 20.92 hours.
output:
After sweep 1 energy=-337.23217296396973 maxlinkdim=4800 maxerr=3.46E-05 time=70248.585
After sweep 2 energy=-337.2389290839439 maxlinkdim=4800 maxerr=3.78E-05 time=75324.747

MPSKit (0.12.6)
Symmetry: SU(2)
Max bond dimension: D = 4000
Time per sweep: 31.66 hours.
output:
[ Info: DMRG2 1: obj = -3.164484013797e+02 err = 9.1771680577e-01 time = 4.42 min
[ Info: DMRG2 2: obj = -3.292824077866e+02 err = 1.7325899602e-02 time = 17.07 min
[ Info: DMRG2 3: obj = -3.348343234046e+02 err = 4.5124942262e-03 time = 5.39 hr
[ Info: DMRG2 4: obj = -3.360135495584e+02 err = 8.0470138898e-04 time = 31.66 hr

Both calculations were performed on the same hardware, using a 64-core, 128-thread CPU machine. This seems counterintuitive, as I expected that the SU(2) symmetry would help reduce the computational cost due to block structure. Is this behavior expected, or might I have misconfigured something in my setup?
Thanks again for your hard work on this project!

[lkdvos]

Thanks for sharing these results, this is definitely very valuable! Would you mind sharing a bit more details about what model you are trying out? I agree that we are expecting that the symmetry should lead to an increase in performance, so I'm happy to check in more detail what's going on.

If you prefer to share in a more confident setting, also do feel free to email me about the details, I'm definitely also happy to try and debug what's going on that way.

[drichg]

Thanks for sharing these results, this is definitely very valuable! Would you mind sharing a bit more details about what model you are trying out? I agree that we are expecting that the symmetry should lead to an increase in performance, so I'm happy to check in more detail what's going on.

If you prefer to share in a more confident setting, also do feel free to email me about the details, I'm definitely also happy to try and debug what's going on that way.

Thank you very much for your quick and positive response!

I’m studying a J1-J2-J3 Heisenberg model on a 12×36 square lattice cylinder geometry.

Below is the main part of my MPSKit code for setting up the calculation:
`
function main(J, J2, J3, L, W, Dcut, filename; symmetry=[], state=[])

initial_bond_dim = 10

if symmetry == SU2Irrep && state == []

    physical_space = SU2Space(1 // 2 => 1)

    virtual_space = SU2Space(0 => initial_bond_dim, 1 // 2 => initial_bond_dim, 1 => initial_bond_dim)

    state = FiniteMPS(L * W, physical_space, virtual_space)

end
 lattice = FiniteCylinder(L, L * W)
H = hamiltonian(J, J2, J3, lattice; symmetry=symmetry)

groundstate, envs, delta = find_groundstate(state, H,
    DMRG2(; trscheme=truncdim(Dcut), maxiter=5,
        tol=1e-6))

save(filename, "groundstate", groundstate)

return groundstate, envs

end`
If you think it would help, I can also share more details about the lattice or Hamiltonian construction.

Thank you very much for your kind offer to help debug.

Looking forward to your insights.

[lkdvos]

I'm happy to help, but it would be easier for me to have an idea of what's going on if I have some code that could reproduce the example. In this case, I'm looking for the definition of the Hamiltonian function, just so I can run a profiler to see what's going on.

EDIT: it would also be nice just to see what the ITensor code is you are comparing this too, I'm happy to investigate the differences!

[drichg]

I'm happy to help, but it would be easier for me to have an idea of what's going on if I have some code that could reproduce the example. In this case, I'm looking for the definition of the Hamiltonian function, just so I can run a profiler to see what's going on.

EDIT: it would also be nice just to see what the ITensor code is you are comparing this too, I'm happy to investigate the differences!

Thank you for your response!

I've just sent the code to your email for your reference. Please let me know if you didn't receive it. Looking forward to your feedback, and many thanks for your support!

[ZongYongyue]

Based on my usage experience of the two packages, the difference is likely caused by the different default parameter settings, such as the Lanczos default Krylov space dimension of ITensor is 3, while in MPSKit it is 30 and so on.

[drichg]

Based on my usage experience of the two packages, the difference is likely caused by the different default parameter settings, such as the Lanczos default Krylov space dimension of ITensor is 3, while in MPSKit it is 30 and so on.

I've just set the MPSKit.Defaults.krylovdim=3, and test the Heisenberg model on the 8*24 cylinder. But it seems not work.
MPSKit output:
[ Info: DMRG2 1: obj = -1.213500932450e+02 err = 9.0701102162e-01 time = 2.43 min
[ Info: DMRG2 2: obj = -1.264925681821e+02 err = 6.0410417120e-03 time = 1.42 min
[ Info: DMRG2 3: obj = -1.269778764295e+02 err = 4.2692908963e-04 time = 4.22 min
[ Info: DMRG2 4: obj = -1.269852659285e+02 err = 1.1090712938e-05 time = 4.82 min
ITensor output:
After sweep 1 energy=-126.92640345731412 maxlinkdim=400 maxerr=1.01E-04 time=114.704
After sweep 2 energy=-127.04513516791893 maxlinkdim=400 maxerr=3.19E-05 time=85.454
After sweep 3 energy=-127.04801351477134 maxlinkdim=400 maxerr=3.29E-05 time=83.365
After sweep 4 energy=-127.04833641221003 maxlinkdim=400 maxerr=3.33E-05 time=83.451
After sweep 5 energy=-127.04841561129761 maxlinkdim=400 maxerr=3.34E-05 time=83.187
After sweep 6 energy=-127.04844590430321 maxlinkdim=400 maxerr=3.36E-05 time=83.553.

[lkdvos]

Let me recap some of my findings and answers here:

As @ZongYongyue mentioned, there are some important differences in the default settings of the different packages that make direct comparisons a bit harder.

The first important caveat is that the way the eigensolvers are handled is very different. In principle, in order to obtain a similar treatment of the local eigenvalue solvers you would have to use:

alg = DMRG2(; ..., alg_eigsolve=(; krylovdim=3, maxiter=1, eager=false))
groundstate, envs, delta = find_groundstate(state, H, alg)

I.e. the default ITensors settings just apply the Hamiltonian 3 times and then simply return.
The main issue with simply doing this is that MPSKit is set up to print a warning whenever the eigensolver claims that it is not converged, which in this case will be every time, by design.
We could try and discuss whether or not this is necessary, and have the option to switch these warnings off, but I'm not entirely sure if that is fair.

The second major difference is that the ITensors version of OpSum automatically "compresses" the resulting MPO, by truncating near zero singular values in the spectrum. In contrast, MPSKit works with sparse MPO objects which, depending on the model and geometry can be more or less efficient.
Unfortunately, here we really have one such case where our MPO constructions leads to an MPO which has a bit too many terms, further hindering the performance.

I'm working on at least supporting some version of this compression scheme right now, I'll link the PR below, and it is definitely interesting to compare the results. Some preliminary testing shows indeed that the timings may vary greatly, on my machine for the 8 * 24 cylinder, spin 1//2 Heisenberg I end up with the following timings, where I locally disabled the fixedpoint warnings for MPSKit:

julia> alg = DMRG2(; trscheme=truncdim(Dcut), maxiter=5, tol=1e-6, alg_eigsolve=(; krylovdim=3, maxiter=1, verbosity=0, eager=false, dynamic_tols=false))
DMRG2{KrylovKit.Lanczos{KrylovKit.ModifiedGramSchmidt2, Float64}, SDD, typeof(MPSKit.Defaults._finalize)}:
 * tol: 1.0e-6
 * maxiter: 5
 * verbosity: 3
 * alg_eigsolve: Lanczos{ModifiedGramSchmidt2, Float64}(ModifiedGramSchmidt2(), 3, 1, 1.0e-10, false, 0)
 * alg_svd: SDD()
 * trscheme: TruncationDimension(400)
 * finalize: _finalize


julia> groundstate, envs, delta = find_groundstate(state, H, alg);
[ Info: DMRG2   1:      obj = -1.202223551735e+02       err = 5.1562064104e-01  time = 4.55 sec
[ Info: DMRG2   2:      obj = -1.264430753877e+02       err = 5.7545598623e-03  time = 23.12 sec
[ Info: DMRG2   3:      obj = -1.269688219386e+02       err = 3.7078659017e-04  time = 35.38 sec
[ Info: DMRG2   4:      obj = -1.269749195564e+02       err = 9.6834750497e-06  time = 31.81 sec
┌ Warning: DMRG2 cancel 5:      obj = -1.269751605557e+02       err = 4.9052342919e-06  time = 2.10 min

Bonus, with the compressed hamiltonian:

julia> groundstate2, envs2, delta = find_groundstate(state, H2, alg);
[ Info: DMRG2   1:      obj = -1.202223551735e+02       err = 5.1562064104e-01  time = 3.42 sec
[ Info: DMRG2   2:      obj = -1.264430753877e+02       err = 5.7545598623e-03  time = 17.85 sec
[ Info: DMRG2   3:      obj = -1.269688219386e+02       err = 3.7078659017e-04  time = 30.08 sec
[ Info: DMRG2   4:      obj = -1.269749195564e+02       err = 9.6834750496e-06  time = 27.00 sec
┌ Warning: DMRG2 cancel 5:      obj = -1.269751605557e+02       err = 4.9052342923e-06  time = 1.74 min

[drichg]

Let me recap some of my findings and answers here:

As @ZongYongyue mentioned, there are some important differences in the default settings of the different packages that make direct comparisons a bit harder.

The first important caveat is that the way the eigensolvers are handled is very different. In principle, in order to obtain a similar treatment of the local eigenvalue solvers you would have to use:

alg = DMRG2(; ..., alg_eigsolve=(; krylovdim=3, maxiter=1, eager=false))
groundstate, envs, delta = find_groundstate(state, H, alg)
I.e. the default ITensors settings just apply the Hamiltonian 3 times and then simply return. The main issue with simply doing this is that MPSKit is set up to print a warning whenever the eigensolver claims that it is not converged, which in this case will be every time, by design. We could try and discuss whether or not this is necessary, and have the option to switch these warnings off, but I'm not entirely sure if that is fair.

The second major difference is that the ITensors version of OpSum automatically "compresses" the resulting MPO, by truncating near zero singular values in the spectrum. In contrast, MPSKit works with sparse MPO objects which, depending on the model and geometry can be more or less efficient. Unfortunately, here we really have one such case where our MPO constructions leads to an MPO which has a bit too many terms, further hindering the performance.

I'm working on at least supporting some version of this compression scheme right now, I'll link the PR below, and it is definitely interesting to compare the results. Some preliminary testing shows indeed that the timings may vary greatly, on my machine for the 8 * 24 cylinder, spin 1//2 Heisenberg I end up with the following timings, where I locally disabled the fixedpoint warnings for MPSKit:

julia> alg = DMRG2(; trscheme=truncdim(Dcut), maxiter=5, tol=1e-6, alg_eigsolve=(; krylovdim=3, maxiter=1, verbosity=0, eager=false, dynamic_tols=false))
DMRG2{KrylovKit.Lanczos{KrylovKit.ModifiedGramSchmidt2, Float64}, SDD, typeof(MPSKit.Defaults._finalize)}:

• tol: 1.0e-6
• maxiter: 5
• verbosity: 3
• alg_eigsolve: Lanczos{ModifiedGramSchmidt2, Float64}(ModifiedGramSchmidt2(), 3, 1, 1.0e-10, false, 0)
• alg_svd: SDD()
• trscheme: TruncationDimension(400)
• finalize: _finalize

julia> groundstate, envs, delta = find_groundstate(state, H, alg);
[ Info: DMRG2 1: obj = -1.202223551735e+02 err = 5.1562064104e-01 time = 4.55 sec
[ Info: DMRG2 2: obj = -1.264430753877e+02 err = 5.7545598623e-03 time = 23.12 sec
[ Info: DMRG2 3: obj = -1.269688219386e+02 err = 3.7078659017e-04 time = 35.38 sec
[ Info: DMRG2 4: obj = -1.269749195564e+02 err = 9.6834750497e-06 time = 31.81 sec
┌ Warning: DMRG2 cancel 5: obj = -1.269751605557e+02 err = 4.9052342919e-06 time = 2.10 min
Bonus, with the compressed hamiltonian:

julia> groundstate2, envs2, delta = find_groundstate(state, H2, alg);
[ Info: DMRG2 1: obj = -1.202223551735e+02 err = 5.1562064104e-01 time = 3.42 sec
[ Info: DMRG2 2: obj = -1.264430753877e+02 err = 5.7545598623e-03 time = 17.85 sec
[ Info: DMRG2 3: obj = -1.269688219386e+02 err = 3.7078659017e-04 time = 30.08 sec
[ Info: DMRG2 4: obj = -1.269749195564e+02 err = 9.6834750496e-06 time = 27.00 sec
┌ Warning: DMRG2 cancel 5: obj = -1.269751605557e+02 err = 4.9052342923e-06 time = 1.74 min

With this setup, I can confirm that the performance has improved significantly compared to my earlier runs with the default parameters. In particular, for models with long-range interactions, MPSKit now runs noticeably faster than ITensors in my tests, which is really impressive! Thanks again for your time and support.

[ZongYongyue]

@lkdvos Hi Lukas , the second point is interesting. What do you think are the respective advantage ranges of these two ways? For example, in large-scale calculations, due to the large number of terms and the linear operations between MPO and MPS take longer, does using automatically compressed MPO have greater advantages in both performance and memory overhead?

[lkdvos]

Well, I would have to experiment a bit with it myself, but it also somewhat depends on how things are implemented, and technically I'm still missing a couple algorithms to have the maximal performance. In what follows, I always mean lossless compression.

To elaborate: when creating the MPOHamiltonian object, currently we are doing the most naive thing and simply creating a big sparse representation of the finite state machine, without trying to do anything smart. This is in general suboptimal.

The first alternative, which is what I'm now implementing and which more closely resembles ITensors implementation is to create the large object, and then find some reasonable way of compressing it (losslessly), typically through a combination of QR and SVD decompositions. This is in some sense optimal, giving the MPO with the smallest bond dimension for a given hamiltonian. The main drawbacks here are that this can be quite expensive for "large" hamiltonians (think quantum chemistry models etc), and that this will spoil the sparsity that is naturally present in these operators.
See also https://journals.aps.org/prb/abstract/10.1103/PhysRevB.102.035147

The second alternative is based on more graph-theoretic approaches, and explicitly performs resummations without making the hamiltonian dense. From what I gather the state-of-the-art algorithms here is https://arxiv.org/pdf/2502.18630.

I'm not entirely sure when either of these options will be the most beneficial.
It will probably depend a lot on the symmetry, the exact values of the bond dimension, and how sparse the resulting operator can remain.

I definitely want to eventually support both, unfortunately I just haven't gotten around to it yet. All help and contributions obviously always welcome.

[ZongYongyue]

Thank you very much for your detailed explanation and sharing of cutting-edge insights in this area!