TDVP (2-site) stuck in the initial state

[adityab999]

Hi, for the given Hamiltonian below, 2-site TDVP (v0.4 as well as v0.3) is always stuck in the initial product state in its real-time evolution. It’s happening not just with the specific initial state written below but also with a few other product states that are relevant to me for this model. It’s the first time I am seeing 2-TDVP being unable to get out of the initial state. Please let me know if it’s something I am doing trivially wrong, and if not so then if there’s some way to resolve this issue.

N = 40  
τ = 0.01
T = 1

s = siteinds("Qubit",2*N; conserve_qns=false)  

mu1 = 1 
mu2 = 1 
x = 1 

os1 = OpSum()      

for j=1:2*N-2

    os1 += -x/4,"X + iY",j,"Z",j+1,"X - iY",j+2
    os1 += -x/4,"X - iY",j,"Z",j+1,"X + iY",j+2

end

for j=1:N

    os1 += mu1/2*(-1)^(j+1),"Id",2*j-1
    os1 += mu1/2*(-1)^(j+1),"Z",2*j-1
    os1 += mu2/2*(-1)^(j+1),"Id",2*j
    os1 += mu2/2*(-1)^(j+1),"Z",2*j
    
end


for j=1:N-1

    osj = OpSum()

    for k=1:j

        osj += (-1)^k,"Id",k
        osj += 1/2,"Z",2*k-1
        osj += 1/2,"Z",2*k

    end    

    os1 += Ops.expand(osj*osj)

end

H = MPO(os1,s)


function step(; sweep, bond, half_sweep)
  if bond == 1 && half_sweep == 2
    return sweep
  end
  return nothing
end


function current_time(; current_time, bond, half_sweep)
  if bond == 1 && half_sweep == 2
    return current_time
  end
  return nothing
end


function return_state(; state, bond, half_sweep)
  if bond == 1 && half_sweep == 2
    return state
  end
  return nothing
end


obs = observer("steps" => step, "times" => current_time, "psis" => return_state)

ψ0 = MPS(s, n -> isodd(n) ? "Dn" : "Up")

@show inner(ψ0', H, ψ0) / inner(ψ0, ψ0)

@time ψf = tdvp(H,-im*T,ψ0; time_step=-im*τ, reverse_step=true, normalize=true, maxdim=20, cutoff=1e-10, (observer!)=obs, outputlevel=1, updater_kwargs=(; krylovdim=6, tol=1e-10, maxiter=200), nsite=2)

@show norm(ψ0)
@show norm(ψf)

steps = obs.steps
psis = obs.psis
times = obs.times

[mtfishman]

Have you tried the new expand function (GitHub - ITensor/ITensorTDVP.jl: Time dependent variational principle (TDVP) of MPS based on ITensors.jl.)?

[adityab999]

Thanks for pointing it out, should have tested it out earlier.

Just tried it with both ψ0 = expand(psi0, H; alg="global_krylov", krylovdim=2, cutoff=1e-9) (also with krylovdim=4,8,10), as well as ψ0 = expand(psi0, [psi01,psi02]; alg="orthogonalize", cutoff=1e-9). Unfortunately, still stuck at the initial (or at some) product state with maxlinkdim=1 being shown all the way.
Here, psi0 = MPS(s, n -> isodd(n) ? "Dn" : "Up"),
psi01 = MPS(s, n -> isodd(n) ? "Up" : "Up") ,
psi02 = MPS(s, n -> isodd(n) ? "Up" : "Dn")

Will input other reference states later in case non-product reference states are actually needed, such as by doing a quick DMRG on this Hamiltonian and supplying the resultant state as a reference state, if that makes any sense.

[adityab999]

Update - the issue is resolved for a certain other type of initial product states that are physically-plausible configurations in some extremities of parameter values… The previously mentioned Neel product state and some related ones that I had tried earlier (not mentioned above for brevity) are all rather highly excited states in the parameter regimes that I am concerned with, so that’s probably a reason why TDVP had a hard time with evolving them (perhaps a very high value of krylovdim to be supplied to expand()) ?

It is worth mentioning that without using the expand function, the issue still persisted. It is resolved only with this. So thanks indeed for insisting on updating ITensorTDVP and using expand().

[mtfishman]

Thanks for the detailed report, it is very helpful to get feedback about successes and failures of TDVP so we can know if the tools we are providing are sufficient for most use cases or if the code or algorithms need to be improved. EDIT: The algorithm is still under active development and to this day there are still new methods being developed to efficiently improve the stability and accuracy of the algorithm, see for example [2403.00562] Comment on "Controlled Bond Expansion for Density Matrix Renormalization Group Ground State Search at Single-Site Costs" (Extended Version).

Would you mind sharing a (formatted) math formula for the Hamiltonian you are studying? That will be useful as a reference for us and others.

Something you could try is performing the initial time evolution (up to some short time) with TEBD and then switch over to TDVP, which was used in [2009.10736] Stripes, Antiferromagnetism, and the Pseudogap in the Doped Hubbard Model at Finite Temperature (see Fig. 9).

If the starting state is very highly excited, indeed maybe that causes issues for global Krylov, since the intuition behind the method is based on projecting out highly excited states. I don’t have much personal experience and intuition about the method, mostly just anecdotal reports from colleagues. However what I have heard is that global Krylov seems to work in cases when “all else fails”, for example local subspace expansion methods, so is generally considered a robust option so if it is struggling in your case you must have a particularly challenging Hamiltonian evolution to perform.

Another thing is that there are a few hidden/experimental options in expand you can try, for example:

expand(state, operator; alg="global_krylov", krylovdim=2, cutoff=1e-12, apply_kwargs=(; maxdim=20, cutoff=1e-10))

cutoff controls the cutoff of the truncation of the states being expanded (setting it lower will potentially include more states), and apply_kwargs controls how the MPO operator is applied to the MPS state to generate the Krylov vectors in the global Krylov method.

[adityab999]

Hi,

Sure, the following is the Hamiltonian (slightly different in notation at one place compared to what I posted in the code snippet above) -

\small H = -x \sum_{j=1}^{2N-2} [ \sigma_j^+ \sigma_{j+1}^z \sigma_{j+2}^- + hc] + \sum_{j=1}^{N}[ \frac{\mu_1{(-1)^j}}{2}(1+\sigma_{2j-1}^z)] + \small \sum_{j=1}^N [\frac{\mu_2{(-1)^j}}{2}(1+\sigma_{2j}^z)] + \sum_{j=1}^N [\sum_{k=1}^j ((-1)^k + \sigma_{2k-1}^z + \sigma_{2k}^z ) ]^2 .

This results from Jordan-Wigner of 2-flavour Schwinger model (with staggered fermions), see the appendix of https://arxiv.org/pdf/1611.00705.

As one can see, the last term (with the double summations) results in arbitrary long-distance 2-qubit terms (with constant coefficients, not separation-dependent) which result from having integrated out the U(1) gauge field which can always be done in 1d with OBC making use of the Gauss’ law. Because of this term, TEBD is hard to naively implement and also will be inefficient, if not impossible, to do. Hence my almost entire reliance on TDVP for time-evolving this and other (qubitized-) lattice gauge theories in 1d.

I thought more and experimented with a few more initial product states and have come to the following conclusion - the issue is not that the algorithm is getting stuck in initial states that are highly excited ones for a given set of parameters, but that it’s getting stuck in initial states that are “unphysical” due to violating an underlying Gauss law (that exists explicitly in the original form of the model before integrating out the gauge degree of freedom) though that may not be directly apparent from the spin Hamiltonian written above. For some other initial product states that are also fairly highly excited and satisfy Gauss’ law, it works fine. It’s cool to see that the TDVP engine seemingly wouldn’t accept (and can recognize?) Gauss law-violating initial states, thus raising an alert to make the user look harder into it

So all in all, I think the newer TDVP with expand() is working as expected I think. Thanks for pointing out the additional options within expand().

[mtfishman]

Thanks for the update, glad to hear tdvp and expand appear to be working well and you’ve found a physical justification for why some initial states are failing to evolve.