DMRG does nothing in certain cases if initial state chosen deterministically? (Julia)

Question

I am finding a problem with something seemingly very simple. Basically, the following Julia code implements a Hamiltonian that is just Sz on each of N spin-1/2 sites. Clearly the ground state energy is -N/2.

N = 100
sites = siteinds("S=1/2", N)
ampo_H = AutoMPO()
for j = 1:N
  add!(ampo_H, 1, "Sz", j)
end
H = MPO(ampo_H, sites)
sweeps = Sweeps(5)
maxdim!(sweeps, 10,20,100,100,200)
cutoff!(sweeps, 1E-10)
init_state = [n % 3 == 0 ? "Up" : "Dn" for n=1:N]
psi0 = productMPS(sites, init_state)
# psi0 = randomMPS(sites)
energy, psi = dmrg(H, psi0, sweeps)
println("Ground state energy = $energy")

However, if you run the code, you will not get a ground state energy of -50. DMRG will simply stay put at an energy of -17, which you may notice is just the energy of the initial state. Feel free to try other deterministic choices of psi0; you should see the same thing happen.

In contrast, if you remove the deterministic choice of psi0 and uncomment the line that says "psi0 = randomMPS(sites)", you should now see that DMRG immediately converges to the right answer, as one would expect.

I know that I have run DMRG on a lattice formulation of the Schwinger model (Kogut-Susskind staggered lattice with Jordan-Wigner transformation and open boundary conditions) with successful results. I used a deterministic initial state for that, and that lattice was also on spin-1/2 sites, so I'm not sure why this one isn't working.

What exactly is going on here? If DMRG is looking for the ground state energy, why should it matter whether the initial state was chosen randomly or deterministically? Any clarification on this would be greatly appreciated.

Comments

On a completely unrelated note, there also seems to be a funny output error in this code? I've never seen it before except here. Here is what the code mentioned in the original question prints out when executed:

    After sweep 1 energy=-17.000000000000 maxlinkdim=1 time=0.078
    2 energy=-17.000000000000 maxlinkdim=1 time=0.086
    After sweep 3 energy=-17.000000000000 maxlinkdim=1 time=0.100
    After sweep 4 energy=-17.000000000000 maxlinkdim=1 time=0.084
    After sweep 5 energy=-17.000000000000 maxlinkdim=1 time=0.085
    Ground state energy = -17.0

For some reason, the words "After sweep " just don't show up for the second line of output! This happens with remarkable consistency. Also, this happens even if you use a different deterministic choice of psi0 or use the "psi0 = randomMPS(sites)" line instead. There was one exception in the case of choosing psi0 randomly--in that one case, this happened:

    0.000000000000 maxlinkdim=1 time=0.124
    After sweep 2 energy=-50.000000000000 maxlinkdim=1 time=0.088
    After sweep 3 energy=-50.000000000000 maxlinkdim=1 time=0.093
    After sweep 4 energy=-50.000000000000 maxlinkdim=1 time=0.107
    After sweep 5 energy=-50.000000000000 maxlinkdim=1 time=0.082
    Ground state energy = -50.0

In this case, the text "After sweep 1 energy=-5" got cut from the first line of output. But this happened exactly once and then never again.

I also tried changing the value of N. It seems that if you increase N, then it has a change of behavior at around N=135 and starts messing up the first line of output instead, as shown below:

1 energy=68.000000000000 maxlinkdim=1 time=0.133
After sweep 2 energy=68.000000000000 maxlinkdim=1 time=0.107
After sweep 3 energy=68.000000000000 maxlinkdim=1 time=0.108
After sweep 4 energy=68.000000000000 maxlinkdim=1 time=0.138
After sweep 5 energy=68.000000000000 maxlinkdim=1 time=0.112
Ground state energy = 68.0

In addition, when N starts going lower, then the one line that gets messed up moves later and later, until eventually the output comes out normally. My best guess is that the output gets messed up for the sweep where the total time first exceeds some amount, which I would estimate at about 0.13 seconds. But then again, I haven't observed this on any code other than this one.

Given that I've never seen this before today and that it has no ties to the actual computations, as far as I can tell, I'm pretty sure this is just some utterly insignificant bug. But I thought I should still mention it. I also wasn't sure whether this merited its own question; at first it seemed like a one-time thing, so I figured it did not, and then when I tested it further it ballooned into a more consistent trend on this particular code.

Answer 1

Hi Sujay,
Thanks for the question. This is a good example of the DMRG "sticking problem" which is a long-standing issue with how the DMRG algorithm works, barring additional convergence tricks. Many people new to DMRG can find it surprising, so it's a question we get often and I mention it in the Guidelines for Asking Questions: http://itensor.org/support/2003/guidelines-for-asking-questions

First of all, as a formal statement DMRG is not rigorously guaranteed to converge and return the true ground state of the system. So care must be taken when working with challenging or unusual Hamiltonians, or with systems with a high degree of symmetry, to check that DMRG is working correctly. Examples of good checks include comparing to exact limits or comparing to exact diagonalization of small systems.

At a more intuitive level, it's not a bad cartoon of DMRG to think of it as a fancy version of the "power method" for finding dominant eigenvectors of matrices. So just like in the power method, if the initial guess vector is orthogonal to the eigenvector one seeks, then convergence will not happen. You can think of DMRG as taking the initial vector and enhacing the part that overlaps with the ground state while suppressing all of the other components. But if there's no overlap with the ground state at all, this can't be done.

So that's why random MPS are actually in some sense the best possible starting state: they have a non-zero overlap with the ground state with probability approaching 1.0 in the limit of the random MPS having a large bond dimension. In fact, one option we provide in the randomMPS function is the ability to make randomMPS of different bond dimensions by doing randomMPS(sites,chi) where chi is the requested bond dimension. If you choose it larger, you will generally see faster convergence of DMRG in terms of number of sweeps.

Best regards,
Miles

Comments

OK, that's very helpful, in particular the comment about orthogonality to the desired state. I didn't realize that this was a limitation of DMRG, so that's good to know. Thank you so much! I do have a few follow-ups though:

(1) Is there a way to choose the initial state randomly under QN conservation in the Julia code? And in particular, choose the initial state randomly so that it falls within a given sector? Because I believe this was not possible two months ago, but I don't see whether a method to this effect was ever added.

(2) Also, do you have any idea about the weird printout mistakes? Have you ever seen anything like this before? Again, I hadn't until last night.

---

Hi Sujay,
Yes, there is a version of randomMPS that additionally takes a "state vector" which is either an array of integers or array of strings specifying the initial product state to use as the starting point for randomization. The total QN of the random MPS will be the same as whatever the state vector's is.

For example, to make a random MPS of bond dimension 4 for a spin 1/2 model that has total Sz of zero:

    state = [isodd(n) ? 1 : 2 for n=1:N]
    psi = randomMPS(sites,state,4)

It's documented here:
https://itensor.github.io/ITensors.jl/dev/MPSandMPO.html#MPS-Constructors-1

Regarding (2), I ran the code you posted above, except with one change which is that I changed the input of siteinds to N instead of 2*N. (Otherwise I was getting an error.) But when I ran the code I didn't get any strange printout, so I'm not sure how to reproduce that issue. If you could post an updated code below that reproduces that printing issue I'd be happy to know about it and fix it -  thanks!

Miles