VUMPS - 2D models and memory consumption

I’ve been trying to reproduce the results of the VUMPS paper for the 2D Heisenberg model, so far without success since my code consumes surprisingly large amounts of RAM: more than 100 GB for the relatively modest N=W=6, D=256 (the run ended prematurely since I only allocated 100GB).

Here’s my definition of the model (rotated to restore TI in the GS):

function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenberg2D", W = 6)
  opsum = OpSum()

  for i in 1:W

    # up

    opsum += -0.5, "S+", i, "S-", mod(i+1, W)
    opsum += -0.5, "S-", i, "S+", mod(i+1, W)
    opsum += "Sz", i, "Sz", mod(i+1, W)

    # right

    opsum += -0.5, "S+", i, "S-", i+W
    opsum += -0.5, "S-", i, "S+", i+W
    opsum += "Sz", i, "Sz", i+W
  end

  return opsum
end

Other than the model definition and setting N=6 & maxdim=256 I’ve used the Heisenberg example as-is. Did I make a mistake in defining the model, thus causing the abnormal memory consumption? Or is it normal after all?

Hello,do you solve this problem now?

Sadly, no. Did you also run into it?

Thanks for the report. It is very possible there is a performance bug in the VUMPS implementation. At some point I found some performance issues which I fixed in these pull requests:

but I mostly tested for models with next nearest neighbor interactions and not longer range interactions like ones that show up in 2D systems.

It may also depend on the Hamiltonian representation you are using, are you using InfiniteSum{MPO}?

Thanks for the response!

I’ve started working with ITensor before these fixes, I’ll update and try again (though they seem unlikely to cause significantly increased memory consumption, no?)

Yes. I’m using the VUMPS Heisenberg example as a template.

They improved both the memory consumption and computation time, since they improved some contraction sequences that weren’t being performed in the optimal ordering. Often with tensor network calculations memory and time go hand-in-hand (tensor contractions which involve larger tensors/aren’t done in an optimal sequence use more memory and take longer to perform).

You can see that before optimizing, the example their used 365.922 GiB of memory, while by the end the example was using 8.782 GiB. Note that is the output from @time, which tracks all memory allocations throughout the calculation, so would be much larger than the peak memory usage.

Hi,mtfishman

l have run the code following the Model"heisenberg2D" as Vladimir has post,with the optimizations eager=true and InfiniteSum{MPO} and N=W=6

But I still found that the error “Out of memory” will be reported, and the program will be interrupted. When using ITensor.jl to calculate two-dimensional problems, it seems that you can use write to disk to avoid this situation. I would like to ask whether there is such a function in ITensorInfiniteMPS, or is there any other way to avoid such error?

I guess a relevant question is how does the VUMPS (peak) memory usage scale. My original post was based on the assumption that VUMPS’s memory usage is similar to iDMRG, which in turn should be as good as or better than finite DMRG in typical use cases.

Hi, mtfishman

As discussed above, I expect to repeat the two dimensional models on infinite cylinders, such as the 2d heisenberg model, from this VUMPS paper. But there are some difficulties. As Vladimir said, these problems have been plaguing us for a long time. Next I will show my error code, it may have some naive mistakes.But I truly hope to get your help.

using MKL
using ITensorInfiniteMPS
using ITensors

include(
  joinpath(
    pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src", "vumps_subspace_expansion.jl"
  ),
)

##############################################################################
# VUMPS parameters
#

maxdim = 256 # Maximum bond dimension
cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
max_vumps_iters = 100 # Maximum number of iterations of the VUMPS algorithm at each bond dimension
vumps_tol = 1e-6
conserve_qns = false
outer_iters = 10 # Number of times to increase the bond dimension
eager = true # Introduced in this PR

##############################################################################
# CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
#

N = 6 # Number of sites in the unit cell

initstate(n) = isodd(n) ? "↑" : "↓"
s = infsiteinds("S=1/2", N; conserve_qns, initstate)
ψ = InfMPS(s, initstate)

function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenberg2D", W = 6)
    opsum = OpSum()
  
    for i in 1:W
  
      # up
  
      opsum += -0.5, "S+", i, "S-", mod(i+1, W)
      opsum += -0.5, "S-", i, "S+", mod(i+1, W)
      opsum += "Sz", i, "Sz", mod(i+1, W)
  
      # right
  
      opsum += -0.5, "S+", i, "S-", i+W
      opsum += -0.5, "S-", i, "S+", i+W
      opsum += "Sz", i, "Sz", i+W
    end
  
    return opsum
  end
model = Model("heisenberg2D")

# Form the Hamiltonian
H = InfiniteSum{MPO}(model, s)

# Check translational invariance
# println("\nCheck translation invariance of the initial VUMPS state")
# @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))

vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, eager)
subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

ψ = vumps_subspace_expansion(H, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs)

# Check translational invariance
# println()
# println("==============================================================")
# println()
# println("\nCheck translation invariance of the final VUMPS state")
# @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))

# Sz = [expect(ψ, "Sz", n) for n in 1:N]

energy_infinite = expect(ψ, H)
@show energy_infinite

The output before the error is

And the error is

ERROR: LoadError: OutOfMemoryError()
Stacktrace:
  [1] Array
    @ ./boot.jl:459 [inlined]
  [2] similar
    @ /public2/kevinH/packages/packages/NDTensors/y9JPS/src/similar.jl:45 [inlined]
  [3] similar
    @ /public2/kevinH/packages/packages/NDTensors/y9JPS/src/similar.jl:48 [inlined]
  [4] similar
    @ /public2/kevinH/packages/packages/NDTensors/y9JPS/src/dense/dense.jl:91 [inlined]
  [5] _similar_from_dims
    @ /public2/kevinH/packages/packages/NDTensors/y9JPS/src/tensor.jl:213 [inlined]
  [6] _similar_from_dims
    @ /public2/kevinH/packages/packages/NDTensors/y9JPS/src/tensor.jl:207 [inlined]
  [7] similar
    @ /public2/kevinH/packages/packages/NDTensors/y9JPS/src/tensor.jl:180 [inlined]
  [8] contraction_output
    @ /public2/kevinH/packages/packages/NDTensors/y9JPS/src/dense/dense.jl:564 [inlined]
  [9] contraction_output
    @ /public2/kevinH/packages/packages/NDTensors/y9JPS/src/generic_tensor_operations.jl:28 [inlined]
 [10] contract(T1::NDTensors.DenseTensor{Float64, 7, NTuple{7, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, labelsT1::NTuple{7, Int64}, T2::NDTensors.DenseTensor{Floa                t64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, labelsT2::Tuple{Int64, Int64, Int64}, labelsR::NTuple{8, Int64})
    @ NDTensors /public2/kevinH/packages/packages/NDTensors/y9JPS/src/generic_tensor_operations.jl:66
 [11] contract(::Type{NDTensors.CanContract{NDTensors.DenseTensor{Float64, 7, NTuple{7, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, NDTensors.DenseTensor{Float64, 3,                 Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}}}, t1::NDTensors.DenseTensor{Float64, 7, NTuple{7, Index{Int64}}, NDTensors.Dense{Float                64, Vector{Float64}}}, labels_t1::NTuple{7, Int64}, t2::NDTensors.DenseTensor{Float64, 3, Tuple{Index{Int64}, Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}                }, labels_t2::Tuple{Int64, Int64, Int64})
    @ NDTensors /public2/kevinH/packages/packages/NDTensors/y9JPS/src/generic_tensor_operations.jl:51
 [12] contract
    @ /public2/kevinH/packages/packages/SimpleTraits/l1ZsK/src/SimpleTraits.jl:331 [inlined]
 [13] _contract(A::NDTensors.DenseTensor{Float64, 7, NTuple{7, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, B::NDTensors.DenseTensor{Float64, 3, Tuple{Index{Int64}, I                ndex{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}})
    @ ITensors /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:1870
 [14] _contract(A::ITensor, B::ITensor)
    @ ITensors /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:1876
 [15] contract(A::ITensor, B::ITensor)
    @ ITensors /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:1980
 [16] #233
    @ /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:2066 [inlined]
 [17] BottomRF
    @ ./reduce.jl:81 [inlined]
 [18] afoldl
    @ ./operators.jl:549 [inlined]
 [19] _foldl_impl
    @ ./tuple.jl:277 [inlined]
 [20] foldl_impl
    @ ./reduce.jl:48 [inlined]
 [21] mapfoldl_impl
    @ ./reduce.jl:44 [inlined]
 [22] #mapfoldl#258
    @ ./reduce.jl:162 [inlined]
 [23] mapfoldl
    @ ./reduce.jl:162 [inlined]
 [24] #foldl#259
    @ ./reduce.jl:185 [inlined]
 [25] foldl
    @ ./reduce.jl:185 [inlined]
 [26] contract(As::NTuple{4, ITensor}; sequence::String, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:2066
 [27] contract
    @ /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:2062 [inlined]
 [28] #contract#237
    @ /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:2076 [inlined]
 [29] contract
    @ /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:2076 [inlined]
 [30] #*#239
    @ /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:2086 [inlined]
 [31] *(::ITensor, ::ITensor, ::ITensor, ::ITensor)
    @ ITensors /public2/kevinH/packages/packages/ITensors/PisGD/src/itensor.jl:2086
 [32] generate_twobody_nullspace(ψ::InfiniteCanonicalMPS, H::InfiniteSum{MPO}, b::Tuple{Int64, Int64}; atol::Float64)
    @ ITensorInfiniteMPS /public2/kevinH/packages/packages/ITensorInfiniteMPS/Ogcvb/src/subspace_expansion.jl:84
 [33] subspace_expansion(ψ::InfiniteCanonicalMPS, H::InfiniteSum{MPO}, b::Tuple{Int64, Int64}; maxdim::Int64, cutoff::Float64, atol::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tu                ple{}, NamedTuple{(), Tuple{}}})
    @ ITensorInfiniteMPS /public2/kevinH/packages/packages/ITensorInfiniteMPS/Ogcvb/src/subspace_expansion.jl:195
 [34] subspace_expansion(ψ::InfiniteCanonicalMPS, H::InfiniteSum{MPO}; kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol}, NamedTuple{(:cutoff, :maxdim), Tuple{Float64, Int64}}                })
    @ ITensorInfiniteMPS /public2/kevinH/packages/packages/ITensorInfiniteMPS/Ogcvb/src/subspace_expansion.jl:321
 [35] macro expansion
    @ /public2/kevinH/packages/packages/ITensorInfiniteMPS/Ogcvb/examples/vumps/src/vumps_subspace_expansion.jl:14 [inlined]
 [36] macro expansion
    @ ./timing.jl:263 [inlined]
 [37] tdvp_subspace_expansion(H::InfiniteSum{MPO}, ψ::InfiniteCanonicalMPS; time_step::Float64, outer_iters::Int64, subspace_expansion_kwargs::NamedTuple{(:cutoff, :maxdim), Tuple{                Float64, Int64}}, vumps_kwargs::NamedTuple{(:tol, :maxiter, :eager), Tuple{Float64, Int64, Bool}})
    @ Main /public2/kevinH/packages/packages/ITensorInfiniteMPS/Ogcvb/examples/vumps/src/vumps_subspace_expansion.jl:9
 [38] #vumps_subspace_expansion#9
    @ /public2/kevinH/packages/packages/ITensorInfiniteMPS/Ogcvb/examples/vumps/src/vumps_subspace_expansion.jl:25 [inlined]
 [39] top-level scope
    @ /public2/kevinH/entanglement/2dheisenberg.jl:94
 [40] include(fname::String)
    @ Base.MainInclude ./client.jl:476
 [41] top-level scope
    @ REPL[3]:1
in expression starting at /public2/kevinH/entanglement/2dheisenberg.jl:94

It seems like I can’t get the right result because

  1. OutOfMemoryError(), program running stopped
  2. The energy in the unit cell does not converge, which is different from the previous result of one-dimensional heisenberg dmrg

So , in order to get the correct 2d heisenberg model on infinite cylinders. What should I do? Any answers and suggestions from you will be of great help to me!!!

The peak memory of VUMPS should scale the same with the MPS bond dimension as DMRG/iDMRG, the issues being reported here are likely due to a performance bug in the VUMPS implementation in ITensorInfiniteMPS.jl (not inherent to the VUMPS algorithm).

@kevinh thanks for sharing a code, that is helpful.

I can confirm that I see the same issue even with the latest version of the ITensorInfiniteMPS.jl package, so there must be some lingering performance bug in the VUMPS implementation that wasn’t addressed by my optimizations in Optimizations for VUMPS by mtfishman · Pull Request #64 · ITensor/ITensorInfiniteMPS.jl · GitHub.

It looks like there is a performance issue in the subspace expansion step, where the infinite MPS bond dimension is increased. That is where I’m seeing the code crash with an out-of-memory error, and also seeing the code using more memory and taking more time to run than I would expect.

From profiling just the subspace expansion step (using ProfileView.jl), I see it spending most of the time in this contraction: ITensorInfiniteMPS.jl/subspace_expansion.jl at db941c6c97ba006a2af441f74949811c05f03bb5 · ITensor/ITensorInfiniteMPS.jl · GitHub

Likely there is a relatively simple fix by improving the order of the tensor contractions in the function generate_twobody_nullspace along the lines of the optimizations I did in Optimizations for VUMPS by mtfishman · Pull Request #64 · ITensor/ITensorInfiniteMPS.jl · GitHub.

@kevinh @Vladimir I tracked down the performance bug in the subspace expansion, it should be fixed on the latest main branch.

The fix was implemented here: Optimize subspace expansion for long range interactions by mtfishman · Pull Request #67 · ITensor/ITensorInfiniteMPS.jl · GitHub.

2 Likes

Hi,Vladimir

@Vladimir
l think you may know about this algorithm, do you know how to solve the last question I asked mtfishman?

Hi,mtfishman

@mtfishman

Thank you very much for fixing this bug. Now the code is running and working properly.

My next question is supposed to have nothing to do with code. Take the result in the paper as an example. The parameter I refer to is Table IV,width=6,D=256. The result of the code is

while in paper:

Table

So I’m confused about the following question

  1. Different from the one-dimensional code result (as shown in the figure energy1d), the two dimension result is not uniform. Is this the correct result?

  1. How can I get the results from the code that correspond to this paper in Table IV

Thank you again for your help, which is very useful for learners!!!

I don’t yet have any experience with VUMPS, no. I do think there might be a bit of a physics issue here, though.

I defined the antiferromagnetic Heisenberg Hamiltonian in a slightly unusual basis wherein all the spins on one of the two sublattices are flipped; this is useful since in this basis the ground state is ferromagnetic (translationally invariant) rather than antiferromagnetic (doubled unit cell), allowing us to use N=W rather than N=2W.

The (potential) issue is that the initial state you used is a good starting point in the original basis, but not in the rotated one. If you were to run the simulation with quantum number conservation this initial state would definitely be a problem since it doesn’t belong to the ground state’s total magnetization sector. It could still be a problem even with quantum number conservation turned off, by making the algorithm find a local minimum of the energy rather than the global one.

In fact, running your code with quantum number conservation turned on seems to give the exact same results as you’ve posted. Not sure how to interpret this.

if you dont mind me to ask in your question…

i’ve come across the memory comsumption problem.
I have use the qn to reduce the dimension but it allocate much more memory than it actually uses.
Is this allocation necessary?

Precision error 9.707990079574844e-6 reached tolerance 1.0e-5, stopping VUMPS after 53 iterations (of a maximum 100).
ERROR: LoadError: OutOfMemoryError()
Stacktrace:
  [1] Array
    @ ./boot.jl:448 [inlined]
  [2] similar
    @ ~/.julia/packages/NDTensors/y9JPS/src/similar.jl:45 [inlined]
  [3] similar
    @ ~/.julia/packages/NDTensors/y9JPS/src/dense/dense.jl:80 [inlined]
  [4] zeros
    @ ~/.julia/packages/NDTensors/y9JPS/src/dense/dense.jl:100 [inlined]
  [5] _zeros(TensorT::Type{NDTensors.DenseTensor{ComplexF64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.Dense{ComplexF64, Vector{ComplexF64}}}}, inds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}})
    @ NDTensors ~/.julia/packages/NDTensors/y9JPS/src/dense/dense.jl:211
  [6] zeros
    @ ~/.julia/packages/NDTensors/y9JPS/src/dense/dense.jl:215 [inlined]
  [7] dense(T::NDTensors.DiagBlockSparseTensor{ComplexF64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.DiagBlockSparse{ComplexF64, Vector{ComplexF64}, 2}})
    @ NDTensors ~/.julia/packages/NDTensors/y9JPS/src/blocksparse/diagblocksparse.jl:387
  [8] array(T::NDTensors.DiagBlockSparseTensor{ComplexF64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.DiagBlockSparse{ComplexF64, Vector{ComplexF64}, 2}})
    @ NDTensors ~/.julia/packages/NDTensors/y9JPS/src/blocksparse/diagblocksparse.jl:299
  [9] array(T::ITensor)
    @ ITensors ~/.julia/packages/ITensors/PisGD/src/itensor.jl:2772
 [10] top-level scope
    @ ~/ITensorInfiniteMPS.jl/imps/zigzagsu3/vumps_su3.jl:111
in expression starting at /data/home/jhzhang/ITensorInfiniteMPS.jl/imps/zigzagsu3/vumps_su3.jl:111
260.626271 seconds (1.51 G allocations: 378.015 GiB, 14.20% gc time, 0.04% compilation time)
456.286747 seconds (2.76 G allocations: 542.244 GiB, 13.48% gc time, 15.94% compilation time)


I’m sorry for not reading the relevant discussions carefully.
Thanks for your effort in solving the issue.
Waiting for good news, thank you.

If you use a 2-site unit cell, there is no guarantee the results will be exactly uniform up to a single site, what you are seeing is normal. In general errors in translational invariance seem to correspond to errors in observables of the state. Running for more iterations and using larger bond dimensions can make the state and observables more uniform within the unit cell.

I’m not sure why the results aren’t matching that table. In case you don’t know, probably you should take an average of the energies that you printed, which correspond to the bond energies of each bond in a single unit cell of the uniform MPS, but it looks like that average would be a bit off from the number reported in that table. Perhaps you need to converge your state better with more iterations.

Have you tried updating to the latest version of the package and seeing if that fixes our of memory errors you are seeing?

I have re-installed the packages and Julia today.
But the problem seems to be the same.
My system has three sites and has a next-neighbor hamiltonian.

Precision error 9.32681722119151e-6 reached tolerance 1.0e-5, stopping VUMPS after 31 iterations (of a maximum 100).
 60.058589 seconds (457.97 M allocations: 72.982 GiB, 23.60% gc time, 0.59% compilation time)

Increase bond dimension 6 out of 6, starting from dimension 99
cutoff = 1.0e-8, maxdim = 100
  0.067100 seconds (327.43 k allocations: 63.284 MiB, 12.00% gc time)

Run VUMPS with new bond dimension 99
Using VUMPS solver with time step -Inf
Running VUMPS with multisite_update_alg = parallel
VUMPS iteration 1 (out of maximum 100). Bond dimension = 99, energy = ([-1.703077624471701, -1.7030776267160406, -1.70307761886441], [-1.7030775761440335, -1.703077577785663, -1.70307757020592]), ϵᵖʳᵉˢ = 7.941985298259432e-6, tol = 1.0e-5, iteration time = 2.016377081 seconds
Precision error 7.941985298259432e-6 reached tolerance 1.0e-5, stopping VUMPS after 1 iterations (of a maximum 100).
  2.016616 seconds (15.33 M allocations: 2.442 GiB, 24.27% gc time)
137.447451 seconds (825.87 M allocations: 118.390 GiB, 18.22% gc time, 30.51% compilation time)

Real memory occupation grows as time goes on.
It grows 0.05G/s to the limit when doing a update with a constant bond dimension 928.

and I run some extra updating on same node, they both survives.
until now I don’t receive the out_of_memory error