Extract effective hamiltonian from LocalMPO

Question

Dear Miles,

I am trying to implement a TDVP algorithm with ITensor. I am following the structure of the DMRG algorithm in the tutorial since the logic of TDVP is the same.
One needs to exponentiate the effective Hamiltonian for the TDVP time evolution. In the DMRG code the effective Hamiltonian is computed as

auto Heff = LocalMPO(H);

psi.position(b);
Heff.position(b,psi);

The question is: how do I extract the actual tensor giving the effective Hamiltonian from Heff? Ideally then one should exponentiate this with expHermitian or some other method (do you have suggestions on this?)

Thanks a lot for your time!

Jacopo.

Answer 1

Hi Jacopo,
Glad to hear that. If you get the TDVP code working well, please consider sharing it, like on Github, and we'd be happy to advertise it on the ITensor website (the "Codes" page). Just a thought but take your time. We're planning to make a TDVP of our own too, but it might take some months to complete with the other things we're doing simulatenously.

To answer your question, the LocalMPO class lets you retrieve the projections of the MPO H into the left- and right-hand MPS basis by calling the .L() and .R() methods. So if you call Heff.position(j), then Heff.L() should be all MPO tensor from 1 up to j-1 projected into the MPS tensors (bra and ket). And likewise Heff.R() should be the MPO tensors from j+2 to N projected.

The reason that .R() goes from j+2 is that it's set up for doing 2-site DMRG. You may want 1-site for TDVP so you may need to modify the logic of the LocalMPO class for your own purposes, or if you'd like it's something we can probably easily do (and should, since we want to add a 1-site DMRG to ITensor eventually).

To get the middle MPO tensors on site j and j+1 you can just get them from H directly, if you have access to H, or through the LocalMPO class's .H() method which gives you access to H. So like Heff.H().A(j) to get the jth tensor of the MPO inside the LocalMPO.

Also I think you might benefit from reading through the LocalMPO code to make sure you understand how it's working and to confirm everything I said is correct. And you can try the things above which I said and print out the .L() and .R() tensors to make sure they are what you expect in terms of their index structure.

Finally, you might want to consider making your own class that's based on LocalMPO but which is simplified and adapted for your own needs, just using LocalMPO as an starting point.

Best regards,
Miles

Comments

Dear Miles,
thanks a lot for your prompt and detailed answer. Still regarding the TDVP implementation, one needs to apply the exponent of the effective Hamiltonian (which is a 6 legs tensor) to the MPS tensor (3 -legs tensor). I guess the most direct way is to reshape the 6 legs of the Heff into 2 legs and the 3 legs of the MPS into 1 legs and use expHermitian to exponentiate the Heff. What is the best or more direct way to reshape tensors in ITensor ?  Or do you think there are other ways to apply the exponential of Heff to the local MPS matrix ?
Thanks!
Best,
Jacopo.

---

Hi Jacopo,

The `expHermitian` function does not require any reshaping. It accepts an ITensor with pairs of primed and unprimed indices (i.e. {i,j,k,...,i',j',k',...}), and treats the ITensor as a matrix where the rows are the unprimed indices and the columns are the primed indices. So you should just be able to input `Heff` directly into `expHermitian`. See this documentation page for more details: http://www.itensor.org/docs.cgi?vers=cppv2&page=classes/decomp

Cheers,
Matt

---

Hello,

still regarding the construction of the effective Hamiltonian, I noticed that while it can be easily built when sweeping from left to right, when going from right to left there are issues. For example the simple code

auto psi = MPS(sites); //random starting state

for(int b = N; b>1; b--)  {
                          
                    psi.position(b);
                    Heff.position(b,psi);
}
give error when it reach site b=N-1 :

libc++abi.dylib: terminating with uncaught exception of type std::out_of_range: vector
Abort trap: 6

Do you maybe know what could be the source of error here?
Thanks a lot!
All the best,
Jacopo.

---

Hi, so I'm not sure precisely what the error is just from seeing that you got an "out_of_range vector" error, since that could be caused by many things.

But if I had to guess, I would say that it's because you are mixing the psi.position and Heff.position methods. The LocalMPO position method *assumes* that the MPS used to project the MPO has only changed on sites b and b+1. All of the other tensors must stay unchanged especially their index structure. But if you call psi.position it could make big changes to the MPS and make things incompatible. So you have to really understand the exact changes being made to the MPS. Inside of the DMRG function we carefully only update the MPS on two sites and we don't call position on it at any point. So that's probably the same approach you'd want to take here.

Miles

---

Dear Miles,

you are right, indeed that was the issue. Now that I only modify the site b and b + 1 there is no error produced. Thanks a lot! When you say "LocalMPO position method *assumes* that the MPS used to project the MPO has only changed on sites b and b+1" is this because it keeps the memory of the previous effective Hamiltonian? Could you point me out to the piece of code that does that? Also, do you think there could be possible issues in calling Heff.position(b,psi) at each site b and then modifying site b and b+1 or the command Heff.position(b,psi)  is tailored to be called only at sites b ,b+2, b+4 and so on (as in your 2-site dmrg)?
Thanks so much!
Best,
Jacopo.

---

Hi Jacopo, when I say it assumes the MPS hasn't changed, this just has to do with the details of how the LocalMPO::position method works. As you can see, it takes the MPS "psi" as one of its arguments. Then it uses the MPS to advance the left or right projections of the MPO into the MPS basis. If the indices from the previous left or right projections don't match anymore with the next MPS tensor, then things will not work correctly. The piece of code for this is essentially the functions makeL and makeR inside of LocalMPO.

Hope that helps!

Miles

---

Dear Miles,

I am computing the 1-site effective Hamiltonian at site b as following

 Heff = LocalMPO<ITensor>(H);

 for(int b = 1; b<=N-1; b++)  {
                Heff.position(b,psi);
                ITensor G1=  Heff.L() * (  Heff.H().A(b)*Heff.H().A(b+1) ) *  dag(prime(psi.A(b+1))) * psi.A(b+1) * Heff.R() ;

....

However the tensor contraction  necessary to compute G1 seems to be extremely slow  as long as bond dimension is higher than 5 or 6. Are there suggestions to improve the efficiency of the ITensor contractions?
I noticed that when ITensor builds Heff (for example in each DMRG sweeps) is extremely fast, so I must be doing something wrong in the contraction above ...

Thanks a lot!

Best,
Jacopo

---

Hi Jacopo,
Usually the best thing to do is to break down the contraction that you have on a single line there (which is fine once your code is tested) into separate steps, and print the intermediate tensor that results after each step. Take a careful look at each one and make sure it has the correct number of indices and specific indices that you expect. If it's extremely slow, probably there is a mistake with the contraction logic and one of the intermediate tensors has the wrong indices.

Best regards,
Miles

---

Hi,

thanks a lot for your quick reply.
The indices were correct (the code was giving the right result) but indeed it appears that by splitting the contraction this way the code runs much faster

            ITensor G1;
            ITensor AA;
            
            Heff = LocalMPO<ITensor>(H);
            
            
            for(int b = 2; b<=N-2; b++)  {
                
               
                Heff.position(b,psi);
                
                AA =  dag(prime(psi.A(b+1)))*Heff.H().A(b+1)* psi.A(b+1);
                AA = AA * Heff.R();
                G1 =    Heff.L() * Heff.H().A(b);
                G1 = G1 *  AA;
                
So I guess the general rule of thumb to have an efficient contraction is to try to not contract tensors with large ranks?
Thanks!
All the Best,
Jacopo.

---

Hi Jacopo, glad you got it to run much faster. Not only is avoiding contracting tensors with large ranks a good rule of thumb, it is the most crucial consideration to keep in mind when devising efficient tensor network algorithms. In more detail, there is a straightforward theory of how to calculate the computational cost of each tensor contraction: it is just the product of the dimensions of all of the indices participating in the contraction (you count the contracted indices only once in this calculation). For more information see this lecture slide:
https://itensor.org/miles/BrazilLectures/TNAndApplications02.pdf
(see slide number 70 and the slides after)

Answer 2

The effective Ham support in modified localmpo.h and localop.h can be found here.
With LocalMPO object PH

PH.position(b,psi); Tenosr heff; PH.localh(heff);
gives the two-site (b,b+1) effective Hamiltonian held in heff

The most challenging problem is how to calculate e^(-i Heff t)|psi>, a possible way is using the Keylov subspace algorithm. In my practice, it is difficult to implemented this algorithm under itensor, because it needs expnonHermitian, or eigensystem of non-Hermitian matrix.

Comments

Hi, thank you for your answer. Regarding exponentiating the Heff, obviusly the easiest approach is to use expHermitian but there is also a related discussion about this here http://itensor.org/support/1337/krylov-subspace-method-to-the-matrix-exponential-for-tdvp?show=1337#q1337.