Optimizations for VUMPS (#64)

* Add `eager=true` to `eigsolve` for reducing the number of matrix-vector multiplications needed for `AC` and `C`.

* Start improving tensor contraction sequences when computing `AC` and `C` in `InfiniteSum{MPO}`.

* Make sure the phases of `AC` and `C` match when computing the precision error `|AC - AL * C|`.

* Use `splitblocks` when making `MPO`s as an additional optimization.

Author: mtfishman

examples/vumps/vumps_hubbard_extended.jl

@@ -5,11 +5,12 @@ using ITensorInfiniteMPS
 # VUMPS parameters
 #
 
-maxdim = 30 # Maximum bond dimension
+maxdim = 50 # Maximum bond dimension
 cutoff = 1e-6 # Singular value cutoff when increasing the bond dimension
 max_vumps_iters = 200 # Maximum number of iterations of the VUMPS algorithm at each bond dimension
 outer_iters = 5 # Number of times to increase the bond dimension
-localham_type = ITensor # or MPO
+localham_type = MPO # or ITensor
+eager = true
 
 model_params = (t=1.0, U=10.0, V=0.0)
 
@@ -47,7 +48,7 @@ println("\nCheck translational invariance of initial infinite MPS")
 @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...))
 
 outputlevel = 1
-vumps_kwargs = (tol=1e-8, maxiter=max_vumps_iters, outputlevel=outputlevel)
+vumps_kwargs = (tol=1e-8, maxiter=max_vumps_iters, outputlevel=outputlevel, eager)
 subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
 
 # For now, to increase the bond dimension you must alternate
@@ -60,7 +61,7 @@ println("\nRun VUMPS on initial product state, unit cell size $N")
 
 @time for _ in 1:outer_iters
   println("\nIncrease bond dimension")
-  global ψ = subspace_expansion(ψ, H; subspace_expansion_kwargs...)
+  global ψ = @time subspace_expansion(ψ, H; subspace_expansion_kwargs...)
   println("Run VUMPS with new bond dimension")
   global ψ = @time vumps(H, ψ; vumps_kwargs...)
 end
@@ -107,7 +108,7 @@ Hfinite = MPO(model, sfinite; model_params...)
 ψfinite = randomMPS(sfinite, initstate; linkdims=10)
 println("\nQN sector of starting finite MPS")
 @show flux(ψfinite)
-sweeps = Sweeps(30)
+sweeps = Sweeps(10)
 maxdims =
   min.(maxdim, [2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 16, 16, 16, 16, 32, 32, 32, 32, 50])
 @show maxdims
@@ -120,7 +121,7 @@ println("\nEnergy density")
 
 nfinite = Nfinite ÷ 2 - 1
 bsfinite = [(nfinite, nfinite + 1), (nfinite + 1, nfinite + 2)]
-hfinite(b) = ITensor(model, sfinite[b[1]], sfinite[b[2]]; model_params...)
+hfinite(b) = ITensor(model, [sfinite[b[1]], sfinite[b[2]]]; model_params...)
 energy_finite = map(b -> expect_two_site(ψfinite, hfinite(b), b), bsfinite)
 
 Nup_finite = ITensors.expect(ψfinite, "Nup")[nfinite:(nfinite + 1)]

examples/vumps/vumps_ising_extended.jl

@@ -10,6 +10,7 @@ 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
 outer_iters = 10 # Number of times to increase the bond dimension
+eager = true
 
 ##############################################################################
 # CODE BELOW HERE DOES NOT NEED TO BE MODIFIED
@@ -33,7 +34,7 @@ H = InfiniteSum{MPO}(model, s; J=J, J₂=J₂, h=h);
 
 @show norm(contract(ψ.AL[1:N]..., ψ.C[N]) - contract(ψ.C[0], ψ.AR[1:N]...));
 
-vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters)
+vumps_kwargs = (tol=vumps_tol, maxiter=max_vumps_iters, eager)
 subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
 ψ_0 = @time vumps(H, ψ; vumps_kwargs...)
 
@@ -52,18 +53,17 @@ sfinite = siteinds("S=1/2", Nfinite; conserve_szparity=true)
 Hfinite = MPO(model, sfinite; J=J, J₂=J₂, h=h)
 ψfinite = randomMPS(sfinite, initstate; linkdims=10)
 @show flux(ψfinite)
-sweeps = Sweeps(15)
+sweeps = Sweeps(10)
 setmaxdim!(sweeps, maxdim)
 setcutoff!(sweeps, cutoff)
 energy_finite_total, ψfinite = dmrg(Hfinite, ψfinite, sweeps)
 
 nfinite = Nfinite ÷ 2
 Sz_finite = expect(ψfinite, "Sz")[nfinite:(nfinite + N - 1)]
 
-@show (
-  energy_infinite,
-  energy_finite_total / Nfinite,
-  reference(model, Observable("energy"); J=J, h=h, J₂=J₂),
-)
+@show energy_infinite
+@show energy_finite_total / Nfinite
+@show reference(model, Observable("energy"); J=J, h=h, J₂=J₂)
 
-@show (Sz, Sz_finite)
+@show Sz
+@show Sz_finite

src/ITensorInfiniteMPS.jl

@@ -43,7 +43,6 @@ include("models/hubbard.jl")
 include("models/xx.jl")
 include("orthogonalize.jl")
 include("infinitemps_approx.jl")
-include("nullspace.jl")
 include("subspace_expansion.jl")
 include("vumps_generic.jl")
 include("vumps_localham.jl")

src/models/fqhe13.jl

@@ -26,7 +26,7 @@ function ITensors.MPO(::Model"fqhe_2b_pot", s::CelledVector, n::Int64; Ly, Vs, p
   end
   opsum = generate_Hamiltonian(opt)
 
-  return MPO(opsum, [s[x] for x in n:(n + range_model)])
+  return splitblocks(linkinds, MPO(opsum, [s[x] for x in n:(n + range_model)]))
 end
 
 #Please contact Loic Herviou before using this part of the code for production

src/models/heisenberg.jl

@@ -15,7 +15,7 @@ function ITensors.MPO(::Model"heisenberg", s)
     os .+= 0.5, "S-", j, "S+", j + 1
     os .+= "Sz", j, "Sz", j + 1
   end
-  return MPO(os, s)
+  return splitblocks(linkinds, MPO(os, s))
 end
 
 """

src/models/hubbard.jl

@@ -31,7 +31,7 @@ function ITensors.MPO(::Model"hubbard", s; t, U, V)
       opsum .+= U, "Nupdn", n
     end
   end
-  return MPO(opsum, s)
+  return splitblocks(linkinds, MPO(opsum, s))
 end
 
 function ITensors.ITensor(model::Model"hubbard", s; kwargs...)

src/models/ising.jl

@@ -8,7 +8,7 @@ function ITensors.MPO(::Model"ising", s; J, h)
   for n in 1:N
     a .+= -h, "Z", n
   end
-  return MPO(a, s)
+  return splitblocks(linkinds, MPO(a, s))
 end
 
 # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ

src/models/ising_extended.jl

@@ -18,7 +18,7 @@ function ITensors.MPO(::Model"ising_extended", s; J=1.0, h=1.0, J₂=0.0)
       a += -h, "Z", n
     end
   end
-  return MPO(a, s)
+  return splitblocks(linkinds, MPO(a, s))
 end
 
 # H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ - J₂ Σⱼ XⱼZⱼ₊₁Xⱼ₊₂

src/models/models.jl

@@ -24,7 +24,7 @@ end
 
 function InfiniteSum{MPO}(model::Model, s::CelledVector; kwargs...)
   N = length(s)
-  mpos = [MPO(model, s, n; kwargs...) for n in 1:N] #slightly improved version. Note: the current implementation does not really allow for staggered potentials for example
+  mpos = [splitblocks(linkinds, MPO(model, s, n; kwargs...)) for n in 1:N] #slightly improved version. Note: the current implementation does not really allow for staggered potentials for example
   return InfiniteSum{MPO}(mpos, translator(s))
 end
 
@@ -38,7 +38,7 @@ end
 function ITensors.MPO(model::Model, s::CelledVector, n::Int64; kwargs...)
   n1, n2 = 1, 2
   opsum = OpSum(model, n1, n2; kwargs...)
-  return MPO(opsum, [s[x] for x in n:(n + nrange(model) - 1)]) #modification to allow for more than two sites per term in the Hamiltonians
+  return splitblocks(linkinds, MPO(opsum, [s[x] for x in n:(n + nrange(model) - 1)])) #modification to allow for more than two sites per term in the Hamiltonians
 end
 
 # Version accepting IndexSet

src/models/xx.jl

@@ -6,7 +6,7 @@ function ITensors.MPO(::Model"xx", s)
     os .+= 1, "X", n, "X", n + 1
     os .+= 1, "Y", n, "Y", n + 1
   end
-  return MPO(os, s)
+  return splitblocks(linkinds, MPO(os, s))
 end
 
 # H = X₁X₂ + Y₁Y₂