How do we properly add two IQMPO's?

Question

(I'll use the language of ITensorV2 for this question, but I don't believe my question is version specific.)

As discussed by Schollwock (pg 58 and 42), if I wish to add two MPO's, say
Ha = A1 * A2 * ... *AN
and
Hb = B1 * B2 * ... *BN
for some N-site system; then I'll need to form the direct product out of all the local tensors Ai and Bi. In other words, form a larger tensor with Ai and Bi on the diagonal: Mi = diagonal(Ai,Bi), where Ha+Hb = M1 * M2 * ... * MN. From which it follows the new tensor would have dimensions dim(Mi) = dim(Ai) + dim(Bi).

It's not so clear to me how one would do this if Ai and Bi are block sparse IQTensors. If LinkA and LinkB represent the sets of Link indices for A and B, then would I need to create a new set LinkM, where the dimensions of like-QN sectors are summed? What is the best way of doing this, or are there functions which already handle this?

I'm aware there is a function add(Ha,Hb), however, doing this fails with the complaint:

From line 271, file itensor_operators.cc
ITensorT::operator+=: different index structure

In this particular run, Ha and Hb are nearly identical, and have identical quantum number. (I'm happy to share my code if that will help.)

Comments

To add two MPOs, they need to have the same pairs of site indices (i.e. same ids, prime levels and tags). Were they created from the same SiteSets? Posting your code may help, if it is not too long.

---

Yes, I had a similar question for you Eugenio, about when you say Ha and Hb are "nearly" identical. Are they identical in terms of their uncontracted indices and number of tensors?

---

Y'all are really amazing. Thank you for the responses.

Let me be more precise. When I Print() MPOs Ha and Hb, I see the following output:

Ha = 
/--------------IQTensor--------------
r=4 div=(Sz=0,Nf=0) log(scale)=0
IQIndex("site=1",4,Site|515) <In>
("Emp 1",1,Site|493) (Sz=0,Nf=0)
("Up 1",1,Site|707) (Sz=1,Nf=1)
("Dn 1",1,Site|313) (Sz=-1,Nf=1)
("UpDn 1",1,Site|681) (Sz=0,Nf=2)

IQIndex("site=1",4,Site|515)' <Out>
("Emp 1",1,Site|493)' (Sz=0,Nf=0)
("Up 1",1,Site|707)' (Sz=1,Nf=1)
("Dn 1",1,Site|313)' (Sz=-1,Nf=1)
("UpDn 1",1,Site|681)' (Sz=0,Nf=2)

IQIndex("l1",7,Link|557) <In>
("q0_0",7,Link|565) (Sz=0,Nf=0)

IQIndex("l1",7,Link|244) <Out>
("q0_1",7,Link|600) (Sz=0,Nf=0)

for the first site tensor of Ha, and

Hb = 
/--------------IQTensor--------------
r=4 div=(Sz=0,Nf=0) log(scale)=0
IQIndex("site=1",4,Site|515) <In>
("Emp 1",1,Site|493) (Sz=0,Nf=0)
("Up 1",1,Site|707) (Sz=1,Nf=1)
("Dn 1",1,Site|313) (Sz=-1,Nf=1)
("UpDn 1",1,Site|681) (Sz=0,Nf=2)

IQIndex("site=1",4,Site|515)' <Out>
("Emp 1",1,Site|493)' (Sz=0,Nf=0)
("Up 1",1,Site|707)' (Sz=1,Nf=1)
("Dn 1",1,Site|313)' (Sz=-1,Nf=1)
("UpDn 1",1,Site|681)' (Sz=0,Nf=2)

IQIndex("l1",7,Link|557) <In>
("q0_0",7,Link|565) (Sz=0,Nf=0)

IQIndex("l1",7,Link|244) <Out>
("q0_1",7,Link|600) (Sz=0,Nf=0)

for the first site tensor of Hb. The other site tensors follow this pattern -- they appear to have identical index structure.

I've created this Github page for my code, which is written for ITensorV2. I hope this helps y'all reproduce the issue. In my code, Ha=Hupup and Hb=Hdndn. The background: I'm building a long-range Hubbard model as an exercise in constructing MPOs. The Hamiltonian Hupup is the inter-site density-density interaction between the up-spin particles, which has this Hamiltonian and these local tensors (where the interaction is truncated to range=2). Hdndn is the same as Hupup, save the local density operators have the opposite spin; which is what I meant by "nearly identical".

It would really be amazing if sum(Hupup,Hdndn) works for two MPOs. If it doesn't, my backup plan is to use directSum (which I believe Matt wrote). I've had some success summing local IQTensors this way.

I hope I made clear the issue without diverging too much into background!

---

It is a strange error, since it does look like the ITensors of the MPOs have the same indices. If you take a look at the test code here:

https://github.com/ITensor/ITensor/blob/v2.1.1/unittest/mpo_test.cc#L96

you can see a working example of sum using IQMPOs. Is there a reason why you aren't making your MPOs with AutoMPO? Could you try using AutoMPO to make your MPOs and add those to see if it works? The only thing that stands out is that your MPOs seem to have extra link indices on the edges, perhaps sum was not made with that case in mind.

---

You may be right about sum() not permitting that case. I removed the dandling bonds by contracting with a pair of edge vectors, and sum(Ha,Hb) works.

However, I actually want to do idmrg, which I don't believe works with autoMPO. So I'm building finite MPO's from scratch then checking their results compared to what autoMPO produces, then flipping a switch to put the edge vectors in the correct place for idmrg. This being said, even if I place the edge vectors to the left and right of the chain (as described in Heisenberg.h), the sum() function still produces a "different index structure" error. Does sum() not support MPO's constructed for idmrg?

Using your directsum.cc, I coded a custom sum function which does work for both cases -- finite and infinite -- though without doing any compression (unlike ITensor's sum). I'll include it here as part of an example, in case anyone is interested.

(Also, for anyone searching these forums with a related problem: maybe this is a relevant discussion.)

---

That sounds right that it doesn't support the MPOs constructed for iDMRG, though I would have to look into the details of the sum() function (I haven't looked at it in a while). In general the infinite MPS support in C++ ITensor is a bit ad hoc, and I think it really requires devoted iMPS and iMPO classes with their own methods (particularly in light of improved algorithms like VUMPS, which require their own functionality for gauging, etc.).

One thing you could try is explicitly make an MPO of size `N+2`, where the first and last MPO tensors are the edge tensors (possibly with some "dummy" one-dimensional site indices). That should work in ITensor's sum() function.

---

Yes, I would say that officially we don't have a function for adding infinite MPOs. While it's true that the existing functions might do that because we wrote the algorithm in a general enough way, it's not really a tested or guaranteed behavior. It's more like a future area of ITensor where we want to add support for iMPS and iMPO types as a distinct feature, which we plan to do.