Fastest way to create and apply diagonal exponential operator

Question

Hello Miles,

I am trying to time evolve a state using a Hamiltonian, which changes every time step. The Hamiltonian consists of a constant two-site term (for which I evolve using Trotter gates), and a varying single-site term, which is diagonal.
Since the sequence will be run many times, the total time of the process of exponentiating the Hamiltonian and applying it needs to be as fast as possible.
Currently I create the exponent through toExpH, and apply the tensors directly to the MPS right after the forward sweep of Trotter gates, like this:

if (ni1 >= i2)
    {
    psi.svdBond(i1,AA,Fromleft,args);
    auto Ap = psi.A(i1)*expHt.A(i1);
    Ap.mapprime(1,0,Site);
    psi.setA(i1,Ap);
    psi.position(ni1);
    }

Although this method is fairly fast, I need it to be even faster. Is there any way I can effectively truncate the MPS after applying the expHt, or can I somehow utilize the diagonal form of the time dependent Hamiltonian Ht?

Thank you!

Answer 1

Hi, interesting set of questions but it's hard to give a totally general answer, as it could depend on many things.

The first and most important question is: have you done careful profiling/timing of the code to be sure about which part is taking the most time? Sometimes it can be surprising as to which part is taking the longest if you measure.

Assuming you've now done that, then for the two likely slowest parts here are some of the things you can try:

(0) you can just try lowering the accuracy you are asking for (larger cutoff, limiting maxm) or using a bigger time step, and see if you still get acceptable results

(1) if a slow part is the decomposition of the local piece of the wavefunction after changing it (the call to svdBond), then you can make extra sure that the method used to factor the wavefunction tensor is the "denmatDecomp" method and not the "svd" method. The .svdBond method (confusingly) switches between these depending on the particular accuracy parameters you set; if you are setting a rather small cutoff, the code may be using the more accurate svd routine versus the faster denmatDecomp routine.

(2) if a slow part is exponentiating your operators, and you can justify separately exponentiating your diagonal operators without causing too much of an additional Trotter / time-step error (there may not be any extra error; just something to make sure to think about) then yes there would be an analytical expression for the exponential of a diagonal operator of course, so you can just directly set the elements of the associated ITensor/IQTensor you want to construct by using the .set method.

Regarding the .set method, you may be interested to know we just added a more convenient way to set multiple elements:

auto T = ITensor(i,j,k);
T.order(k,j,i); //ensures that the order of the indices is k,j,i
T.set(2,1,3,0.5235); //set the k=2,j=1,i=3 element to 0.5235

Best regards,
Miles

Comments

Hi Miles,

Thank you very much for your response. Just to make clear:

(1) the "denmatDecomp" method requires me to create two new tensors "A" and "B" from the two-site tensor "AA" and then manually set their correct location in the MPS?
I have tried using this approach by setting

    IQTensor A,B;
    denmatDecomp(AA,A,B,Fromleft,args);
    A.mapprime(1,0,Site);
    B.mapprime(1,0,Site);
    psi.setA(i1,A);
    psi.setA(i2,B);
    
but this results in a error regarding arrows.
Should I instead modify the .svdBond() method, such that it is forced to use the denmatDecomp routine?

(2) how does one use the set method for operators. Say I want exponentiate a number operator "N" for time evolution, which of course i diagonal. Should you do
T .set(1,1,1 , exp(0) )
T .set(2,2,2 , exp(1*dt) )
T .set(3,3,3 , exp(2*dt) )
.... and so on ?

---

Sure, good questions.

(1) Yes, after doing a denmatDecomp call, you do have to place the factors A and B back into the MPS similar to what you did. The missing piece that may be causing the bug you're getting is that you have to make sure that either A or B have some of the indices of "AA", namely the ones that you want to end up on A or B. So for example you usually want to put the leftward facing link/virtual index of AA as well as the left-hand site index of AA into the set of indices of A, which tells the denmatDecomp routine to make sure those indices end up on A and the rest of the indices end up on B (and there's a new index connecting A to B). It's confusing to explain in words but pretty intuitive once you get the hang of it. I'd suggest doing Print(A); and Print(B); then PAUSE; right after the call to denmatDecomp to see if A and B are coming out with the indices that you expect, then you can tweak things until they do.

(2) Yes the code you wrote for the .set operator is correct, as long as the ordering of the indices is set to the ordering you want. (Of course for the diagonal case it won't matter.)

---

Hi Miles, sorry for all the questions ...
I did the trick with directly building the exponentiated operator, and it works great! However, I   am still having issues with denmatDecomp. When running the code above I get

Tensor AA:

/--------------IQTensor--------------
r=3 div=QN(5) log(scale)=2.22045e-16
IQIndex("d",11,Link|774) <Out>
  ("d_0",1,Link|168) QN(0)
  ("d_1",2,Link|428) QN(1)
  ("d_2",2,Link|891) QN(2)
  ("d_3",3,Link|536) QN(3)
  ("d_4",2,Link|571) QN(4)
  ("d_5",1,Link|534) QN(5)

IQIndex("Boson 2",5,Site|972) <Out>
  ("Emp 2",1,Site|519) QN(0)
  ("Occ1 2",1,Site|162) QN(1)
  ("Occ2 2",1,Site|213) QN(2)
  ("Occ3 2",1,Site|60) QN(3)
  ("Occ4 2",1,Site|696) QN(4)

IQIndex("Boson 1",5,Site|113) <Out>
  ("Emp 1",1,Site|845) QN(0)
  ("Occ1 1",1,Site|715) QN(1)
  ("Occ2 1",1,Site|695) QN(2)
  ("Occ3 1",1,Site|867) QN(3)
  ("Occ4 1",1,Site|439) QN(4)
\------------------------------------

Tensor A:

/--------------IQTensor--------------
r=4 div=QN(0) log(scale)=0
IQIndex("d",11,Link|774) <Out>
  ("d_0",1,Link|168) QN(0)
  ("d_1",2,Link|428) QN(1)
  ("d_2",2,Link|891) QN(2)
  ("d_3",3,Link|536) QN(3)
  ("d_4",2,Link|571) QN(4)
  ("d_5",1,Link|534) QN(5)

IQIndex("Boson 2",5,Site|972) <Out>
  ("Emp 2",1,Site|519) QN(0)
  ("Occ1 2",1,Site|162) QN(1)
  ("Occ2 2",1,Site|213) QN(2)
  ("Occ3 2",1,Site|60) QN(3)
  ("Occ4 2",1,Site|696) QN(4)

IQIndex("Boson 1",5,Site|113) <Out>
  ("Emp 1",1,Site|845) QN(0)
  ("Occ1 1",1,Site|715) QN(1)
  ("Occ2 1",1,Site|695) QN(2)
  ("Occ3 1",1,Site|867) QN(3)
  ("Occ4 1",1,Site|439) QN(4)

IQIndex("d",1,Link|694) <In>
  ("d_5",1,Link|224) QN(5)
\------------------------------------

Tensor B:

/--------------IQTensor--------------
r=1 div=QN(5) log(scale)=4.44089e-16
IQIndex("d",1,Link|694) <Out>
  ("d_5",1,Link|224) QN(5)
\------------------------------------

It appears like the decomp creates a sort of dummy index, d_5, which is turned into the tensor B, while A gets to keep all the indices of the combined tensor. I am unsure how this happens, and how I can prevent this.

---

Hi, so if you can tell me which indices of "AA" you want to end up on A versus on B then I could write you a quick code example

---

Excellent!
I want the denmatDecomp method to achieve the same as the svdBond method.
Hence, tensor A should have the "Boson 1" Site index, while B should have the "Boson 2" Site index and the "d" link index from AA - keep in mind that this is from the left end of the tensor network. In the general case, A should have an additional link index to the left.
Also, how come the "d" index reaches QN(5), while all the site indices only go up to QN(4)?

Thank you very much for your help!

---

Hi, sorry about the slow reply.

So to make sure the "Bosoy 1" index ends up on tensor A, you need to make sure that tensor A has this index prior to calling denmatDecomp. The denmatDecomp method inspects all of the indices of tensor A, and whichever ones match an index that the input AA also has, stay on tensor A.

So I think if you call it as follows it should do what you want:

auto i1 = sites(1); //get index "Boson 1"; there may be a different way you need to get it

ITensor A(i1); //make sure A has index "i1"
ITensor B;
denmatDecomp(AA,A,B,Fromleft,{"Cutoff=",cutoff,"Maxm=",maxm});

Print(A);
Print(B);