Questions about random initial MPS under QN conservation ("randomMPS(sites, state, linkdim)" in Julia)

Question

I am exploring the "randomMPS(sites, state, linkdim)" function in the Julia version of the code, used to initialize a random MPS under QN conservation. I toyed around with the following basic code:

N = 100
linkdim = 1
sites = Index[]
sites = siteinds("S=1/2", N, conserve_qns=true)
ampo_H = AutoMPO()
for j = 1:N-1
  add!(ampo_H, 1, "S+", j, "S-", j+1)
end
H = MPO(ampo_H, sites)
sweeps = Sweeps(10)
maxdim!(sweeps, 10,20,100,100,200)
cutoff!(sweeps, 1E-10)
init_state = [isodd(n) ? 1 : 2 for n = 1:N]
psi0 = randomMPS(sites, init_state, linkdim)
energy, psi = dmrg(H, psi0, sweeps)
println("Ground state energy = $energy")

When I used linkdim = 0 or linkdim = 1, I got the following (the times were different, but the energies after each sweep were identical):

After sweep 1 energy=-20.093813857715 maxlinkdim=3 time=46.208
After sweep 2 energy=-19.645716775053 maxlinkdim=6 time=0.633
After sweep 3 energy=-23.704160971099 maxlinkdim=21 time=0.837
After sweep 4 energy=-21.976581439869 maxlinkdim=50 time=1.242
After sweep 5 energy=-17.979053940609 maxlinkdim=113 time=2.240
After sweep 6 energy=-12.781825304274 maxlinkdim=200 time=4.182
After sweep 7 energy=-18.735327685101 maxlinkdim=200 time=6.456
After sweep 8 energy=-17.554817334205 maxlinkdim=200 time=6.428
After sweep 9 energy=-11.592260452089 maxlinkdim=200 time=6.489
After sweep 10 energy=-18.150972193173 maxlinkdim=200 time=6.039
Ground state energy = -18.150972193173168

When I used linkdim = 2, I got the following:

After sweep 1 energy=-17.138433694197 maxlinkdim=7 time=16.842
After sweep 2 energy=-19.329082230479 maxlinkdim=20 time=0.785
After sweep 3 energy=-14.965351818888 maxlinkdim=55 time=1.271
After sweep 4 energy=-14.037715046585 maxlinkdim=100 time=2.089
After sweep 5 energy=-18.253177305167 maxlinkdim=200 time=4.063
After sweep 6 energy=-17.967870954570 maxlinkdim=200 time=6.248
After sweep 7 energy=-19.513691584655 maxlinkdim=200 time=6.625
After sweep 8 energy=-18.245438077285 maxlinkdim=200 time=6.188
After sweep 9 energy=-17.493978866991 maxlinkdim=200 time=6.266
After sweep 10 energy=-17.804622748715 maxlinkdim=200 time=5.827
Ground state energy = -17.804622748714813

Finally, when I used linkdim = 10, I got the following:

After sweep 1 energy=-18.432676648899 maxlinkdim=10 time=17.006
After sweep 2 energy=-19.905158282446 maxlinkdim=20 time=0.990
After sweep 3 energy=-23.774621818681 maxlinkdim=65 time=1.544
After sweep 4 energy=-16.035532730501 maxlinkdim=100 time=2.895
After sweep 5 energy=-10.002545372991 maxlinkdim=200 time=5.878
After sweep 6 energy=-25.311587050851 maxlinkdim=200 time=7.119
After sweep 7 energy=-17.279971718254 maxlinkdim=200 time=7.238
After sweep 8 energy=-13.021013397538 maxlinkdim=200 time=7.013
After sweep 9 energy=-18.874581289730 maxlinkdim=200 time=7.202
After sweep 10 energy=-28.130306372441 maxlinkdim=200 time=6.768
Ground state energy = -28.13030637244088

These results were quite perplexing. Furthermore, in a more complicated setup, I ended up getting an error "MPS center bond dim less than requested". (I was trying to set up the Schwinger model, in line with http://itensor.org/support/2140/implementing-local-gauge-symmetries. I can post the code if you want, but I'm not sure how helpful it would be here.) I got this error with several values of linkdim, including 4, 10, and 100, although I did not get it with the default linkdim = 1.

As a result of all of these things, I had a number of questions:

(1) As you can see in all of the above examples, the energy is not decreasing with each sweep as one would naturally expect. Instead it is jumping all over the place. Am I doing something fundamentally wrong here?

(2) As you can see above, one can set linkdim to 0, and that it gives the same answer as linkdim=1. Why is a linkdim of 0 even allowed? Maybe I'm totally misunderstanding what this is representing here, but I thought it had to be at least 1. I assume an answer to this will also explain why the code gives identical results for linkdim=0 and linkdim=1.

(3) What is the source of the error "MPS center bond dim less than requested", and how does one make sure it doesn't happen? I looked for it in the code, but I really didn't understand it.

(4) More broadly, are there some general pointers for how one should choose the linkdim value for a given situation? It seems based that the default of linkdim=1 will just return the state you specified without introducing any randomness, so I am assuming you should not pick this for sure.

Answer 1

Hi Sujay,
Thanks for the detailed feedback: this is really helpful to make the code better. Some of the issues you are running into are typical of newer code, where we have not ironed out every bug or hardened it against every possible input. Let me address each of these as far as I can without running tests of my own:

(1) the reason for the energy jumping around is probably that your Hamiltonian MPO is not Hermitian. You only put S+ S- terms and no S- S+ terms. I’ve seen this kind of behavior before for non-Hermitian inputs to DMRG. Ideally we’d have a check that throws an error if the Hamiltonian is non-Hermitian, as we do not support this case.

(2) An input of linkdim=0 was not intended to be allowed, and is not a meaningful input. I’ll update the code to throw an error if linkdim < 1. I’m not sure what linkdim of zero would mean.

(3) The error about “MPS center bond dim” is something I need to work on. What the randomMPS algorithm currently does is act on the input product state with a random quantum circuit. It keeps applying layers of the circuit until the bond dimension grows to linkdim. So apparently there’s a case where this isn’t happening properly and I’ll look into it.

(4) For this code indeed linkdim = 1 just returns the product state you input. For the QN conserving case it’s hard to think of anything else it could do very easily. If returned state was a superposition of other product states of the same total QN then the bond dimension would have to be greater than 1. On the other hand if it was a totally random product state the total QN guarantee wouldn’t be satisfied, which would be a failure of the purpose of this function. So this behavior is the only consistent one. Really this function becomes more useful for linkdim > 1 and is not particularly useful for linkdim = 1, though I’m not sure we should have that case throw an error since the behavior is meaningful and well defined. I’m open to any suggestions you may have though.

Thanks,
Miles

Comments

Thank you so much Miles!

Regarding (1), you were right. In my desire to make the test case as simple as possible, I forgot that having S+S- but not S-S+ would make the Hamiltonian not Hermitian. I added the conjugate term (S- on site j and S+ on site j+1, for j = 1 to N-1), and I can confirm that the DMRG now decreases the energy monotonically.

Regarding (4), I realize that I worded my question to put too much emphasis on the trivial linkdim=1 case. More broadly, do you have any pointers on what order of magnitude one should put for linkdim, e.g. choosing whether to put 3 or 10 or 30 or 100 (other than just trial and error, I guess)?

---

Hi Sujay,
Glad it's working better. I'm a bit embarrassed about the "MPS center bond dim" errors though I filed an issue and will work on fixing that. It's a very new code after all.

I don't have a general guideline about a best value of linkdim, because it depends on the application and how much need there is for help converging a certain calculation, but I would think about it like this: just the cost alone of constructing the random MPS of some linkdim=chi and then inputting it into DMRG will have a similar cost to doing one DMRG sweep at the same chi. The difference is that it will (hopefully) quickly give an energy which is close to the best one can reach for that chi, and not get stuck in some local minimum because of the random properties (here I think one could actually develop a rather deep theoretical statement about concentration of measure and random SVD algorithms etc. but that's only been done in the literature somewhat for MPS and DMRG).

So if you are aiming to spend a certain "budget" on your DMRG calculation in terms of the max dimension you want to reach, I'd choose linkdim to be something like one-tenth of that max dimension (last sweep maxdim) and also start the maxdim of the first sweep already at the same dimension as the randomMPS you use.

I'm planning myself to do some more experiments with randomMPS, to see how quickly properties converge in the first one or two sweeps of DMRG compared with just random product states.

Best,
Miles