[ITensors] [ENHANCEMENT] Bump default value for eigsolve_krylovdim in DMRG to 4

[jtschneider]

Description

In light of the spectacularly improved convergence of the ground state energy of a critical spin chain, I would suggest increasing the default dimension of the Krylov subspace in the eigensolver from KrylovKit to 4

Minimal runnable code

using ITensors
# from  Eq. (13) in T W Burkhardt and I Guim (1985), J. Phys. A: Math. Gen. 18 L33, https://doi.org/10.1088/0305-4470/18/1/006
exact_crit_ising_energy_OBC(N::Int) =  1.0 - csc(π/(2.0*(2.0*N+1)))

@time begin 
  J = 1.0
  h = 1.0

  N = 80
  χ = 150
  SVDcutoff = 1e-16

  sites = siteinds("S=1/2",N)

  ampo = OpSum()
  for j=1:N-1
    ampo += -J,"X",j,"X",j+1
    ampo += -h,"Z",j
  end
  ampo += -h,"Z",N

  H = MPO(ampo,sites)

  psi0 = randomMPS(sites,10);

  sweeps = Sweeps(12)
  mindim!(sweeps, 5, 10, 20)
  maxdim!(sweeps, floor.(Int, LinRange(10, χ, 10))...)
  cutoff!(sweeps, SVDcutoff)


  energy_3, psi_3 = dmrg(H,psi0, sweeps, eigsolve_krylovdim=3)

  energy_4, psi_4 = dmrg(H,psi0, sweeps, eigsolve_krylovdim=4)

  @show exact_crit_ising_energy_OBC(N)

  @show energy_3
  @show energy_4

  @show abs(energy_3 - exact_crit_ising_energy_OBC(N))
  @show abs(energy_4 - exact_crit_ising_energy_OBC(N))

  @show maxlinkdim(psi_3)
  @show maxlinkdim(psi_4)
  
end

Output

Output of minimal runnable code

After sweep 1 energy=-101.27197171871956  maxlinkdim=10 maxerr=1.41E-03 time=0.085
After sweep 2 energy=-101.42825239181545  maxlinkdim=20 maxerr=1.10E-09 time=0.122
After sweep 3 energy=-101.46263187509906  maxlinkdim=30 maxerr=4.92E-11 time=0.182
After sweep 4 energy=-101.4749475856145  maxlinkdim=40 maxerr=1.18E-12 time=0.283
After sweep 5 energy=-101.48106043338078  maxlinkdim=50 maxerr=5.16E-13 time=0.409
After sweep 6 energy=-101.48516199496554  maxlinkdim=60 maxerr=4.17E-13 time=0.523
After sweep 7 energy=-101.48852231537329  maxlinkdim=69 maxerr=1.94E-13 time=0.693
After sweep 8 energy=-101.49114036933987  maxlinkdim=80 maxerr=8.77E-14 time=0.855
After sweep 9 energy=-101.49322776755494  maxlinkdim=90 maxerr=4.16E-14 time=1.056
After sweep 10 energy=-101.49479749446336  maxlinkdim=100 maxerr=2.71E-14 time=1.274
After sweep 11 energy=-101.49584825672818  maxlinkdim=100 maxerr=2.90E-14 time=1.308
After sweep 12 energy=-101.49651752362102  maxlinkdim=100 maxerr=3.24E-14 time=1.303
After sweep 1 energy=-101.38857321946286  maxlinkdim=10 maxerr=5.60E-04 time=0.104
After sweep 2 energy=-101.47636860931595  maxlinkdim=20 maxerr=4.03E-09 time=0.148
After sweep 3 energy=-101.48875353720699  maxlinkdim=30 maxerr=7.41E-11 time=0.230
After sweep 4 energy=-101.4936643693508  maxlinkdim=40 maxerr=5.15E-11 time=0.333
After sweep 5 energy=-101.49628517973699  maxlinkdim=50 maxerr=9.02E-12 time=0.476
After sweep 6 energy=-101.49717247505477  maxlinkdim=60 maxerr=1.52E-12 time=0.649
After sweep 7 energy=-101.49736173965282  maxlinkdim=69 maxerr=3.05E-13 time=0.855
After sweep 8 energy=-101.49739979242437  maxlinkdim=80 maxerr=3.24E-14 time=1.055
After sweep 9 energy=-101.4974080482929  maxlinkdim=90 maxerr=4.05E-14 time=1.167
After sweep 10 energy=-101.49740929033591  maxlinkdim=100 maxerr=1.74E-15 time=1.122
After sweep 11 energy=-101.49740944279279  maxlinkdim=93 maxerr=1.02E-15 time=0.962
After sweep 12 energy=-101.4974094520415  maxlinkdim=72 maxerr=9.98E-16 time=0.682
exact_crit_ising_energy_OBC(N) = -101.49740945239314
energy_3 = -101.49651752362102
energy_4 = -101.4974094520415
abs(energy_3 - exact_crit_ising_energy_OBC(N)) = 0.0008919287721198543
abs(energy_4 - exact_crit_ising_energy_OBC(N)) = 3.516475999276736e-10
maxlinkdim(psi_3) = 100
maxlinkdim(psi_4) = 72
 15.897904 seconds (11.27 M allocations: 42.475 GiB, 27.51% gc time)
72

As you can see, the DMRG converges up to 4 orders of magnitude when using 3 Krylov vectors, while it converges to a spectacular 6 orders of magnitude more accurately (1e-10 absolutely) to the true ground state energy when using 4 Krylov vectors. Note also that this is not a matter of bond dimension or entanglement, as both states have the same SVD cutoff and maximum link dimension. Bizarrely, the ground state of the algorithm with 4 krylov vectors has a lower bond dimension 🤔

In this vain, one might decide on the layout as you discussed in ITensor/ITensorMPS.jl#49

Note that there is a performance cost of ≈ 16% when running the code away from the critical point. I tested h/J=0.5 and got the following timer values: (note that H is now different), while there is arguably no accuracy improvement...

julia> @time energy_3, psi_3 = dmrg(H,psi0, sweeps, eigsolve_krylovdim=3, outputlevel=0)
  1.508486 seconds (4.85 M allocations: 2.802 GiB, 25.16% gc time, 0.72% compilation time)
(-84.21288175017561, MPS

julia> @time energy_4, psi_4 = dmrg(H,psi0, sweeps, eigsolve_krylovdim=4; outputlevel=0)
  1.750208 seconds (6.20 M allocations: 3.402 GiB, 26.02% gc time)
(-84.21288175017594, MPS

julia> energy_4 - energy_3
-3.268496584496461e-13

Now, clearly this is a design decision on your side. I suppose "only" using 3 times as many Krylov basis vectors as eigenvalues one is asking for (in this case the first eigenvalue, i.e. of the ground state) is typically on the lower end of the recommended factor of Krylov vectors to use in an Arnoldi eigensolver [citation needed].

Let me know what you think of this design choice as well as the one in ITensor/ITensorMPS.jl#49

Cheers,

Version information

• Output from versioninfo():

julia> versioninfo()
Julia Version 1.9.3
Commit bed2cd540a1 (2023-08-24 14:43 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × AMD Ryzen 9 7950X 16-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, znver3)
  Threads: 1 on 32 virtual cores
Environment:
  JULIA_DEPOT_PATH = /home/jan/.julia

• Output from using Pkg; Pkg.status("ITensors"):

using Pkg; Pkg.status("ITensors")
Status `~/.julia/environments/v1.9/Project.toml`
  [9136182c] ITensors v0.3.49

[mtfishman]

Hmm maybe we could just add a note in the documentation about it? Choosing defaults is tough and always requires some tradeoffs.

[jtschneider]

I would be in favour of the increased Krylov dimension, as this is not a standard parameter to be adjusted in the narrative of better convergence of DMRG algorithms such as the increased bond dimension or the reduced SVD cutoff (which is what newcomers adhere to), and this parameter feels a bit more buried behind a few layers of abstraction and is hard to find (particularily for newcomers).

Surely, a note in the documentation is also helpful!

[jtschneider]

In fact, the whooping 6 orders of magnitude improvement are not possible (in the critical Ising chain) with improving the bond dimension or the SVD cut off (or both)!

So I think this has to be documented at least.

[mtfishman]

@emstoudenmire any opinion on this?

I think we would just need a more systematic analysis of a wider range of models to analyze the tradeoffs of accuracy vs. speed before making a decision. But adding a documentation note that increasing eigsolve_krylovdim for models near criticality is a more benign decision and achieves many of the same goals.

I do see that it is already mentioned here: https://itensor.github.io/ITensors.jl/dev/faq/DMRG.html#Ensuring-a-DMRG-calculation-is-converged

[jtschneider]

Ah very good, I have missed that part of the documentation!

The tradeoff vs accuracy certainly needs a much wider analysis than just this one example.
I was thinking the same that one may emphasise the better convergence close to criticality.

[emstoudenmire]

Hi Jan, thanks for starting this discussion and providing a sample code for us to use.

On my machine (M1 Macbook running Julia 1.10-rc1) I get rather different results, however, that are more in line with what we would hope, namely that results are quite good with Krylov dim of 3 and do get better with 4, but not dramatically so.

Here's what I get:

exact_crit_ising_energy_OBC(N) = -101.49740945239314
energy_3 = -101.49740580485742
energy_4 = -101.49740944736686
abs(energy_3 - exact_crit_ising_energy_OBC(N)) = 3.647535720574524e-6
abs(energy_4 - exact_crit_ising_energy_OBC(N)) = 5.026279836783942e-9
maxlinkdim(psi_3) = 20
maxlinkdim(psi_4) = 20

So I wonder what's accounting for the difference? Would you be willing to try Julia 1.10-rc1 to see if perhaps that's it? Otherwise it could be something like the available BLAS or LAPACK on your hardware versus mine (are you using the default BLAS?). It would be good to understand this especially if there's a set of defaults we could shift to that would give more consistent results across hardware and Julia/ITensor versions.

I agree with Matt's point that we should gather more info first. And thanks for your PR improving the docs.

[jtschneider]

Hi Miles, thanks for your reply! I have just ran the code on my PC with Julia 1.10-rc1 and the output is this:

exact_crit_ising_energy_OBC(N) = -101.49740945239314
energy_3 = -101.4965602937405
energy_4 = -101.49740941575797
abs(energy_3 - exact_crit_ising_energy_OBC(N)) = 0.0008491586526417905
abs(energy_4 - exact_crit_ising_energy_OBC(N)) = 3.6635171341004025e-8
maxlinkdim(psi_3) = 150
maxlinkdim(psi_4) = 150

What is absolutely bewildering to me is that your two ground states have maxlinkdim(psi_3) = 20 = maxlinkdim(psi_4) while mine both returned maxlinkdim(psi_3) = 150 = maxlinkdim(psi_4). I can see you used the same system size N=80 as you print exact_crit_ising_energy_OBC(N) = -101.49740945239314 but did you change the SVD cut off? I don't believe so, as your energy values are close to the theoretical value...

I wonder what the hell is going on here 🤔

I read about some tricks about exchanging the BLAS system on an AMD CPU to MKL, and tested the following code where I used MKL on the AMD CPU I have by previously exporting MKL_DEBUG_CPU_TYPE=5, which sadly did not improve the convergence of the DMRG algorithm with eigsolve_krylovdim=3:

Code with MKL

ENV["MKL_DEBUG_CPU_TYPE"] = "5"
using LinearAlgebra, MKL
using ITensors
# from  Eq. (13) in T W Burkhardt and I Guim (1985), J. Phys. A: Math. Gen. 18 L33, https://doi.org/10.1088/0305-4470/18/1/006
exact_crit_ising_energy_OBC(N::Int) =  1.0 - csc(π/(2.0*(2.0*N+1)))

@time begin 
  J = 1.0
  h = 1.0

  N = 80
  χ = 150
  SVDcutoff = 1e-16

  sites = siteinds("S=1/2",N)

  ampo = OpSum()
  for j=1:N-1
    ampo += -J,"X",j,"X",j+1
    ampo += -h,"Z",j
  end
  ampo += -h,"Z",N

  H = MPO(ampo,sites)

  psi0 = randomMPS(sites,10);

  sweeps = Sweeps(12)
  mindim!(sweeps, 5, 10, 20)
  maxdim!(sweeps, floor.(Int, LinRange(10, χ, 10))...)
  cutoff!(sweeps, SVDcutoff)


  energy_3, psi_3 = dmrg(H,psi0, sweeps, eigsolve_krylovdim=3)

  energy_4, psi_4 = dmrg(H,psi0, sweeps, eigsolve_krylovdim=4)

  @show exact_crit_ising_energy_OBC(N)

  @show energy_3
  @show energy_4

  @show abs(energy_3 - exact_crit_ising_energy_OBC(N))
  @show abs(energy_4 - exact_crit_ising_energy_OBC(N))

  @show maxlinkdim(psi_3)
  @show maxlinkdim(psi_4)
  
end

julia> LinearAlgebra.BLAS.get_config()
LinearAlgebra.BLAS.LBTConfig
Libraries: 
├ [ILP64] libmkl_rt.so
└ [ LP64] libmkl_rt.so

julia> versioninfo()
Julia Version 1.10.0-rc1
Commit 5aaa9485436 (2023-11-03 07:44 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × AMD Ryzen 9 7950X 16-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
  Threads: 1 on 32 virtual cores
Environment:
  JULIA_DEPOT_PATH = /home/jan/.julia

Which outputted the same behaviour as before when using OpenBLAS:

exact_crit_ising_energy_OBC(N) = -101.49740945239314
energy_3 = -101.49671243345877
energy_4 = -101.49740945210006
abs(energy_3 - exact_crit_ising_energy_OBC(N)) = 0.0006970189343746824
abs(energy_4 - exact_crit_ising_energy_OBC(N)) = 2.930846676463261e-10
maxlinkdim(psi_3) = 150
maxlinkdim(psi_4) = 95

While testing, I found that I can occasionally get a good result even with 3 krylov vectors because of an advantageous random initial state, I suppose. I ran the @time begin ... end bracket for a couple of times and got one time the following output:

exact_crit_ising_energy_OBC(N) = -101.49740945239314
energy_3 = -101.49739657222555
energy_4 = -101.49740944490975
abs(energy_3 - exact_crit_ising_energy_OBC(N)) = 1.2880167588491531e-5
abs(energy_4 - exact_crit_ising_energy_OBC(N)) = 7.48339346046123e-9
maxlinkdim(psi_3) = 150
maxlinkdim(psi_4) = 121

which makes sense in the context that this result is also achievable by increasing the number of sweeps for example up to 20, so I have to completely retract my above statement that a better result is not so easily achieved with egisolve_krylovdim=3 (note the persistent 8 orders of magnitude difference)

The below output was generated only changing sweeps = Sweeps(12) to sweeps = Sweeps(20) :

After sweep 1 energy=-101.40339650748034  maxlinkdim=10 maxerr=1.55E-03 time=0.050
After sweep 2 energy=-101.46031134289957  maxlinkdim=25 maxerr=2.36E-11 time=0.079
After sweep 3 energy=-101.4704113569282  maxlinkdim=41 maxerr=2.31E-12 time=0.343
After sweep 4 energy=-101.4764805892804  maxlinkdim=56 maxerr=7.78E-13 time=0.532
After sweep 5 energy=-101.4807265475735  maxlinkdim=72 maxerr=5.01E-13 time=0.858
After sweep 6 energy=-101.48403543732168  maxlinkdim=87 maxerr=1.97E-13 time=1.351
After sweep 7 energy=-101.48694783009688  maxlinkdim=103 maxerr=1.89E-13 time=1.614
After sweep 8 energy=-101.48951174716731  maxlinkdim=118 maxerr=1.06E-13 time=2.193
After sweep 9 energy=-101.4917851523758  maxlinkdim=134 maxerr=5.12E-14 time=2.472
After sweep 10 energy=-101.49368195192015  maxlinkdim=150 maxerr=2.78E-14 time=3.150
After sweep 11 energy=-101.49508127970724  maxlinkdim=150 maxerr=3.90E-14 time=3.381
After sweep 12 energy=-101.49601481073685  maxlinkdim=150 maxerr=3.46E-14 time=3.324
After sweep 13 energy=-101.49660647651778  maxlinkdim=150 maxerr=2.33E-14 time=3.042
After sweep 14 energy=-101.49696216413135  maxlinkdim=150 maxerr=1.61E-14 time=2.808
After sweep 15 energy=-101.4971729760927  maxlinkdim=150 maxerr=1.05E-14 time=2.698
After sweep 16 energy=-101.49728451701536  maxlinkdim=150 maxerr=6.33E-15 time=2.784
After sweep 17 energy=-101.49734165321406  maxlinkdim=150 maxerr=3.45E-15 time=2.518
After sweep 18 energy=-101.49737273784434  maxlinkdim=150 maxerr=1.86E-15 time=2.760
After sweep 19 energy=-101.49738987827969  maxlinkdim=150 maxerr=1.17E-15 time=2.542
After sweep 20 energy=-101.4973992694482  maxlinkdim=150 maxerr=7.37E-16 time=2.313
After sweep 1 energy=-101.37291163356562  maxlinkdim=10 maxerr=8.95E-04 time=0.063
After sweep 2 energy=-101.46779740896534  maxlinkdim=25 maxerr=6.49E-10 time=0.268
After sweep 3 energy=-101.4924910367722  maxlinkdim=41 maxerr=6.91E-11 time=0.448
After sweep 4 energy=-101.49653748673155  maxlinkdim=56 maxerr=4.15E-12 time=0.521
After sweep 5 energy=-101.49728144111938  maxlinkdim=72 maxerr=3.06E-13 time=0.810
After sweep 6 energy=-101.49738814570195  maxlinkdim=87 maxerr=4.45E-14 time=1.394
After sweep 7 energy=-101.49740498773751  maxlinkdim=103 maxerr=5.22E-15 time=1.305
After sweep 8 energy=-101.49740853693427  maxlinkdim=118 maxerr=1.01E-15 time=1.675
After sweep 9 energy=-101.49740929514599  maxlinkdim=133 maxerr=9.99E-17 time=1.571
After sweep 10 energy=-101.49740943000546  maxlinkdim=124 maxerr=9.99E-17 time=1.294
After sweep 11 energy=-101.49740945060206  maxlinkdim=106 maxerr=9.96E-17 time=1.153
After sweep 12 energy=-101.49740945231196  maxlinkdim=72 maxerr=9.95E-17 time=0.647
After sweep 13 energy=-101.49740945239259  maxlinkdim=64 maxerr=1.00E-16 time=0.434
After sweep 14 energy=-101.49740945239311  maxlinkdim=50 maxerr=9.98E-17 time=0.347
After sweep 15 energy=-101.49740945239294  maxlinkdim=50 maxerr=9.93E-17 time=0.309
After sweep 16 energy=-101.49740945239277  maxlinkdim=50 maxerr=9.77E-17 time=0.305
After sweep 17 energy=-101.49740945239287  maxlinkdim=50 maxerr=9.76E-17 time=0.317
After sweep 18 energy=-101.49740945239294  maxlinkdim=50 maxerr=9.76E-17 time=0.304
After sweep 19 energy=-101.49740945239249  maxlinkdim=50 maxerr=9.76E-17 time=0.555
After sweep 20 energy=-101.49740945239267  maxlinkdim=50 maxerr=9.76E-17 time=0.387
exact_crit_ising_energy_OBC(N) = -101.49740945239314
energy_3 = -101.4973992694482
energy_4 = -101.49740945239267
abs(energy_3 - exact_crit_ising_energy_OBC(N)) = 1.018294494770089e-5
abs(energy_4 - exact_crit_ising_energy_OBC(N)) = 4.689582056016661e-13
maxlinkdim(psi_3) = 150
maxlinkdim(psi_4) = 50

It becomes evident with more sweeps that the SVD compression reduces the critical Ising state further down to smaller bond dimension with a more efficient representation, but it is only able to do so in my case for eigsolve_krylovdim=4 as you can see in the above output.
In fact, you can see the energy oscillate quite strongly in the case of 3 krylov vectors once it approached the energy up to 5 orders of magnitude, which seems to be its
This could explain your result of only maxlinkdim(psi_gs)=20... I am not so much bewildered by the low bond dimension in this "small" system size of N=80 as by the fact that OpenBLAS seems to only yield this lower and more efficient representation when using 4 krylov vectors, while your BLAS could find it easier...

[jtschneider]

Given that the energy converges for more sweeps and that we put a note into the documentation, I would close this issue.

Perhaps interestingly, I have observed essentially the same behaviour of the convergence of the energy on the nvidia A100, meaning it only converges significantly closer to the true value with more than 15 sweeps and 4 Krylov vectors:

Code

using ITensors
using CUDA

# from  Eq. (13) in T W Burkhardt and I Guim (1985), J. Phys. A: Math. Gen. 18 L33, https://doi.org/10.1088/0305-4470/18/1/006
exact_crit_ising_energy_OBC(N::Int) =  1.0 - csc(π/(2.0*(2.0*N+1)))

platform = NDTensors.cu
# platform = identity

begin 
  J = 1.0
  h = 1.0

  N = 80
  χ = 150
  SVDcutoff = 1e-16

  sites = siteinds("S=1/2",N)

  ampo = OpSum()
  for j=1:N-1
    ampo += -J,"X",j,"X",j+1
    ampo += -h,"Z",j
  end
  ampo += -h,"Z",N

  H = platform(MPO(ampo,sites))

  psi0 = platform(randomMPS(sites,10))

  sweeps = Sweeps(20)
  mindim!(sweeps, 5, 10, 20)
  maxdim!(sweeps, floor.(Int, LinRange(10, χ, 10))...)
  cutoff!(sweeps, SVDcutoff)


  energy_3, psi_3 = dmrg(H,psi0, sweeps, eigsolve_krylovdim=3)

  energy_4, psi_4 = dmrg(H,psi0, sweeps, eigsolve_krylovdim=4)

  @show exact_crit_ising_energy_OBC(N)

  @show energy_3
  @show energy_4

  @show abs(energy_3 - exact_crit_ising_energy_OBC(N))
  @show abs(energy_4 - exact_crit_ising_energy_OBC(N))

  @show maxlinkdim(psi_3)
  @show maxlinkdim(psi_4)
  
end

Output:

exact_crit_ising_energy_OBC(N) = -101.49740945239314
energy_3 = -101.4973974694192
energy_4 = -101.49740945239256
abs(energy_3 - exact_crit_ising_energy_OBC(N)) = 1.1982973944668629e-5
abs(energy_4 - exact_crit_ising_energy_OBC(N)) = 5.826450433232822e-13
maxlinkdim(psi_3) = 150
maxlinkdim(psi_4) = 50

This is obviously nothing dramatic, so I would just close this issue, given that the documentation also contains the note to be careful around critical points.
It is however interesting to note that the BLAS vendor on the Apple chip is apparently better at finding a more compressed ground state than both the OpenBLAS and the CUDA vendors.
Perhaps that is also something you want to note somewhere in the documentation?

Technical output

julia> LinearAlgebra.BLAS.get_config()
LinearAlgebra.BLAS.LBTConfig
Libraries: 
└ [ILP64] libopenblas64_.so

julia> versioninfo()
Julia Version 1.9.3
Commit bed2cd540a1 (2023-08-24 14:43 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 56 × Intel(R) Xeon(R) Gold 6330 CPU @ 2.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, icelake-server)
  Threads: 5 on 56 virtual cores
Environment:
  LD_LIBRARY_PATH = /dragofs/sw/foss/0.2/software/ncurses/6.2-GCCcore-11.2.0/lib:/dragofs/sw/foss/0.2/software/binutils/2.37-GCCcore-11.2.0/lib:/dragofs/sw/foss/0.2/software/zlib/1.2.11-GCCcore-11.2.0/lib:/dragofs/sw/foss/0.2/software/GCCcore/11.2.0/lib64
  JULIA_NUM_THREADS = 4

julia> CUDA.current_device()
CuDevice(0): NVIDIA A100-PCIE-40GB