Questions about DMRG

Question

Hi Miles,

I have some questions for DMRG in ITensor. Hope you could help.

1. What kind of optimization algorithm does it use, single-site or twosite?

2. PBC is difficult. Will the "Efficient matrix-product state method for periodic boundary conditions" be built in this Library?

3. Can it run in parallel?

4. I tried some simulation with DMRG. I get different numbers with different sweeps.maxm(). The code is:

   #include "itensor/all.h"
   using namespace itensor;
   int main()
   {
   int N = 32;
   
   auto sites = SpinOne(N);
   
   auto ampo = AutoMPO(sites);
   for(int j = 1; j < N; ++j)
       {
       ampo += 10.0,"S+",j,"S-",j+1;
       ampo += 10.0,"S-",j,"S+",j+1;
       ampo +=  "Sz*Sz",j;
       }
   ampo += 10.0,"S+",N,"S-",1;
   ampo += 10.0,"S-",N,"S+",1;
   ampo +=  "Sz*Sz",N;
   
   auto H = IQMPO(ampo);
   
   auto state = InitState(sites);
   for(int i = 1; i <= N; ++i)
   {
   if(i%2 == 1)
       state.set(i,"Up");
   else
       state.set(i,"Dn");
   }
   
   auto psi = IQMPS(state);
   
   auto sweeps = Sweeps(200);
   sweeps.maxm() = 20,40,100,100,200;
   sweeps.cutoff() = 1E-10;
   
   auto energy = dmrg(psi,H,sweeps,{"Quiet=",true});
   
   return 0;
   }
   

   For this case, the vN Entropy at center bond converges at 1.834080476998. If I use "sweeps.maxm() = 20,40,60,80,100,120,140,160,180,200;", the vN Entropy at center bond converges at 1.834080201341. The difference is not big. But when I set "N=64", the two vN entropy become 2.059443246913 and 2.059446044498. It seems this difference will be bigger with larger N. Is this normal? Which one I should trust? Is there a rule to set "sweeps.maxm()"?

Many thanks.

Best regards,

Jin

Answer 1

Hi Jin,
Thanks for the questions. Here are some answers for you:

1. Right now the code uses two-site optimization. Single-site can in principle be faster, though one has to account for the extra overhead needed to ensure convergence (things like the noise trick or the Hubig et al. expansion) when making the comparison. We'd like to offer both in the future. The differences are pretty minor for most purposes though.

2. Right now you can do PBC the "poor man's" way by using AutoMPO, say, to connect a long bond from site 1 to site N. This has the disadvantage of increasing the entanglement of the MPS, though. We can't really support the Pippan et al. PBC proposal right now because it's still not a standard method; more research is needed to determine the cases where that method works well and when it runs into trouble. For example, it may not work well enough for critical systems, which are one of the few cases where PBC can really be helpful. For gapped systems often OBC are all you need (just take the system size longer than the correlation length). Also do consider infinite DMRG. If it's necessary to study periodic systems, one can take infinite MPS and put them onto a finite cylinder; for example see arxiv:1208.2623. Also see my article about boundary conditions here: http://itensor.org/docs.cgi?page=articles/periodic

3. ITensor sends tensor contractions to BLAS dgemm, so if your BLAS has multithreading enabled (e.g. if you have set OMP_NUM_THREADS or MKL_NUM_THREADS or VECLIB_NUM_THREADS environment variables greater than 1) then you will see greater than 100% cpu usage for large enough tensor contractions. Other than that right now there's no other multithreading or paralellism happening in ITensor, though we'd like to improve this soon. If you are asking about the real-space parallel DMRG algorithm I wrote about in 2013, that's not publicly available because the code needs to be polished, but if you email me I can send you an old copy of the code if you want to try fixing it up yourself.

4. The way DMRG works, you will always see small differences in the numbers you get depending on the accuracy parameters (maxm, cutoff, etc.). These have a non-trivial and complicated effect on the entanglement of the MPS at each bond. So the behavior you're seeing is normal. For cases where you need to get reliablt converged observables for a given number of states m, say, the standard practice is to take a small cutoff, then do two or more sweeps at the same m. So you would do something like maxm() = 50,50,100,100,200,200; That way the second sweep at each m is closer to being the optimal MPS for that particular number of states.