DMRG energy oscillating rather than converging

DMRG and Numerical Methods
1

Hello, I am trying to use ITensor’s DMRG to find the ground state of a 2D fermionic system.
I am finding unexpected behaviour in that as optimisation progresses, the energy of the DMRG sweep is oscillating / jumping around, even in the limit of giving it a much larger than necessary bond dimension. Due to this stochastic behaviour the final energy returned from the final sweep is not the lowest energy the optimiser found over the course of the sweeps which I find odd…

I am also suspicious of the behaviour due to the bond dimension going from 2 up to 256, without a significant decrease in the energy, whilst the truncation error goes from 0.1 to 1E-12,

Please can you suggest an explanation for this behaviour. I have already gone through the DMRG FAQ’s to try the tips for improving convergence:

  • I have allowed a bond dimension of up to 1024, however it seems to only need a bond dim of 256 at it has a truncation error of around 1e-12 at this point.

  • I have given it a large number of sweeps and it still hasn’t converged.

  • I have added a noise term which decrease in magnitude as the sweeps progress.

Hamiltonian taken from Eq (1) of [1911.06368] Quantum Computing for Neutrino-nucleus Scattering.

Code:

using ITensors
using ITensors.HDF5
using NPZ

function main()
    
    N = 16
    
    D=2
    A = 3
    
    t = 1.0
    C0 = 2 * -7
    D0 = 2 * 28
    
    sites = siteinds("Fermion",N)

    os = OpSum()
    
    # Const
    os += 2*D*t*A, "Id", 1
    
    # Kinetic
    Kinetic_terms = [[(i, i+4), (i, i+8), (i+4, i+12), (i+8, i+12)] for i in 1:4]
        
    for terms in Kinetic_terms
        for term in terms
            os += -t, "Cdag", term[1], "C", term[2]
            os += -t, "C", term[1], "Cdag", term[2]
        end
    end
    
    
    # 2-body
    twobody_terms = [[(i,i+1), (i,i+2), (i,i+3), (i+1,i+2), (i+1,i+3), (i+2,i+3)] for i in [1,5,9,13]]
    
    for terms in twobody_terms
        for term in terms
            os += C0/2, "N", term[1], "N", term[2]
        end
    end
    
    # 3-body
    threebody_terms = [[(i,i+1,i+2), (i,i+1,i+3), (i,i+2,i+3), (i+1,i+2,i+3)] for i in [1,5,9,13]]
    
    for terms in threebody_terms
        for term in terms
            os += D0/2, "N", term[1], "N", term[2], "N", term[3]
        end
    end
    
    H = MPO(os,sites)
    
    psi0 = randomMPS(sites, 16)
    
    nsweeps = 50
    maxdim = [2^i for i in 1:10]
    
    cutoff = [1E-12]

    noise = [1E-6, 1E-6, 1E-7, 1E-7, 1E-8, 1E-8, 1E-10, 1E-10, 1E-12, 1E-12, 0.0]

    energy, psi = dmrg(H,psi0; nsweeps, maxdim, cutoff, eigsolve_krylovdim=5, eigsolve_maxiter=3)
        
    println("Energy: $(energy)")
    
    return energy, psi
end

energy, psi2 = main();

Output:

After sweep 1 energy=-15.994966526037368  maxlinkdim=2 maxerr=1.99E-01 time=4.714
After sweep 2 energy=-15.946111844575244  maxlinkdim=4 maxerr=1.20E-03 time=0.062
After sweep 3 energy=-15.907545825668999  maxlinkdim=8 maxerr=4.13E-04 time=0.033
After sweep 4 energy=-15.976404553400846  maxlinkdim=16 maxerr=9.42E-06 time=0.201
After sweep 5 energy=-16.00435608400482  maxlinkdim=32 maxerr=3.07E-08 time=0.241
After sweep 6 energy=-16.00163306719081  maxlinkdim=64 maxerr=4.54E-09 time=0.493
After sweep 7 energy=-15.999956574113618  maxlinkdim=128 maxerr=1.24E-11 time=0.973
After sweep 8 energy=-15.930308909897244  maxlinkdim=240 maxerr=9.94E-13 time=1.378
After sweep 9 energy=-15.986788492008595  maxlinkdim=228 maxerr=9.80E-13 time=1.600
After sweep 10 energy=-15.999067539235353  maxlinkdim=228 maxerr=9.60E-13 time=1.385
After sweep 11 energy=-15.999491299654544  maxlinkdim=217 maxerr=9.55E-13 time=1.207
After sweep 12 energy=-16.00156450927571  maxlinkdim=204 maxerr=9.58E-13 time=1.764
After sweep 13 energy=-15.966392923185586  maxlinkdim=234 maxerr=9.72E-13 time=1.869
After sweep 14 energy=-15.992740923416052  maxlinkdim=213 maxerr=9.93E-13 time=1.747
After sweep 15 energy=-15.9968859735475  maxlinkdim=235 maxerr=9.18E-13 time=1.791
After sweep 16 energy=-15.999757147621862  maxlinkdim=202 maxerr=9.46E-13 time=1.602
After sweep 17 energy=-15.999927240170724  maxlinkdim=174 maxerr=9.83E-13 time=1.475
After sweep 18 energy=-16.041650727369614  maxlinkdim=126 maxerr=9.70E-13 time=0.996
After sweep 19 energy=-15.949825936439298  maxlinkdim=238 maxerr=9.55E-13 time=1.394
After sweep 20 energy=-16.001678839550845  maxlinkdim=213 maxerr=9.53E-13 time=1.775
After sweep 21 energy=-15.960162484984036  maxlinkdim=206 maxerr=9.85E-13 time=1.494
After sweep 22 energy=-16.169776183102805  maxlinkdim=237 maxerr=8.36E-13 time=1.762
After sweep 23 energy=-15.997452012014097  maxlinkdim=214 maxerr=9.58E-13 time=1.726
After sweep 24 energy=-15.998318845231353  maxlinkdim=190 maxerr=9.88E-13 time=1.414
After sweep 25 energy=-15.998871001788038  maxlinkdim=206 maxerr=9.99E-13 time=1.134
After sweep 26 energy=-15.99908876476493  maxlinkdim=189 maxerr=9.69E-13 time=1.393
After sweep 27 energy=-16.015694891103937  maxlinkdim=212 maxerr=9.71E-13 time=1.283
After sweep 28 energy=-16.02500335412863  maxlinkdim=220 maxerr=9.81E-13 time=2.031
After sweep 29 energy=-16.000950066608734  maxlinkdim=202 maxerr=9.75E-13 time=1.524
After sweep 30 energy=-15.9272280497569  maxlinkdim=243 maxerr=9.62E-13 time=1.546
After sweep 31 energy=-15.963063748119833  maxlinkdim=246 maxerr=9.30E-13 time=1.483
After sweep 32 energy=-15.97998893816897  maxlinkdim=248 maxerr=8.36E-13 time=1.946
After sweep 33 energy=-15.998227656006312  maxlinkdim=250 maxerr=9.95E-13 time=1.800
After sweep 34 energy=-15.988286229947642  maxlinkdim=253 maxerr=9.52E-13 time=1.829
After sweep 35 energy=-15.993651358743104  maxlinkdim=249 maxerr=6.87E-13 time=1.900
After sweep 36 energy=-15.999961041496704  maxlinkdim=199 maxerr=9.67E-13 time=1.623
After sweep 37 energy=-15.999997034276024  maxlinkdim=154 maxerr=9.88E-13 time=1.787
After sweep 38 energy=-15.999941889935851  maxlinkdim=187 maxerr=9.74E-13 time=1.359
After sweep 39 energy=-16.006805349032923  maxlinkdim=232 maxerr=9.26E-13 time=1.810
After sweep 40 energy=-15.932823364379395  maxlinkdim=250 maxerr=8.83E-13 time=1.565
After sweep 41 energy=-15.993328621575355  maxlinkdim=247 maxerr=9.31E-13 time=1.742
After sweep 42 energy=-15.998973781688106  maxlinkdim=231 maxerr=8.85E-13 time=1.905
After sweep 43 energy=-15.994502940992493  maxlinkdim=245 maxerr=8.47E-13 time=1.672
After sweep 44 energy=-15.999774505213233  maxlinkdim=233 maxerr=9.40E-13 time=1.949
After sweep 45 energy=-15.991341761050778  maxlinkdim=246 maxerr=9.33E-13 time=1.493
After sweep 46 energy=-15.99654891122973  maxlinkdim=243 maxerr=9.85E-13 time=1.763
After sweep 47 energy=-16.003838675957997  maxlinkdim=231 maxerr=8.92E-13 time=1.909
After sweep 48 energy=-15.994796958951014  maxlinkdim=237 maxerr=9.94E-13 time=1.734
After sweep 49 energy=-15.919778594313923  maxlinkdim=214 maxerr=9.94E-13 time=2.083
After sweep 50 energy=-16.015028045282456  maxlinkdim=238 maxerr=9.59E-13 time=1.913
Energy: -16.015028045282456

Thanks,
Joe

2

Hi Joe,
I believe your hopping term is defined incorrectly, i.e. has the wrong sign in the second term. To be Hermitian, it should be either:

            os += -t, "Cdag", term[1], "C", term[2]
            os += +t, "C", term[1], "Cdag", term[2]

or else

            os += -t, "Cdag", term[1], "C", term[2]
            os += -t, "Cdag", term[2], "C", term[1]
3

Thanks Miles, that was well spotted, seems to have solved the issue!

4

Yes when DMRG acts up that much, it’s usually the Hamiltonian not being Hermitian. Can be subtle!

6

Hi miles
How do we deal with this oscillation if the Hamiltonian is not Hermitian, or is there an algorithm like dmrg？

7

The good news is you can use a non-Hermitian Hamiltonian in the ITensor dmrg code. You just need to pass the keyword argument ishermitian = false. There are also options to specify which eigenvalue DMRG should target such as the one with largest absolute value, or largest real part, and so on. The documentation on our DMRG function provides more details about these options.