Strange .orthogonalize behaviour when operators are multiplied

Question

Hi,
while digging deeper into my previous question, I found the following behaviour which totally confuses me. It might be a bug I can't find a workaround for.

Try to run following code:

  Hubbard sites(4);

  IQMPO testA(sites);
  testA.Anc(2) = sites.op("Cdagdn", 2) *
                 setElt(leftLinkInd(testA, 2)(1)) *
                 setElt(rightLinkInd(testA, 2)(1));

  testA.Anc(3) = sites.op("Cdn", 3) *
                 setElt(leftLinkInd(testA, 3)(1)) *
                 setElt(rightLinkInd(testA, 3)(1));
  testA.orthogonalize();
  std::cout << "Checking testA";
  checkQNs(testA);
  std::cout << " successfull" << std::endl;

  IQMPO testB(sites);
  testB.Anc(2) = multSiteOps(sites.op("Cdagdn", 2), sites.op("Id", 2)) *
                 setElt(leftLinkInd(testB, 2)(1)) *
                 setElt(rightLinkInd(testB, 2)(1));
  testB.Anc(3) = sites.op("Cdn", 3) *
                 setElt(leftLinkInd(testB, 3)(1)) *
                 setElt(rightLinkInd(testB, 3)(1));
  std::cout << "Checking testB";
  testB.orthogonalize();
  checkQNs(testB);
  std::cout << " successfull" << std::endl;

You will find something like

Checking testA succesful
Checking testBcheckQNs: At site 4 to the right of the OC, Left side Link not pointing In
From line 250, file mps/mpo.cc

Incorrect Arrow in IQMPO

Incorrect Arrow in IQMPO
Aborted (core dumped)

This means if an IQMPO is constructed from a siteset and you put an operator at some place it works (testA). But as soon as you multiply this operator by one (or any other site operator), the orthogonalization does not work anymore and the IQMPO breaks.

I already checked the multSiteOps function and tried to use putMPOLinks instead of the setElt(left....) thing. Both approaches do not solve the problem. The single workaround I figured out so far is to multiply the site operators by hand and put them together in a new created operator, but I hope there is something better.

I hope you have an idea.

Best,
Lars

Comments

Hi Lars,
I finally got the new .orthogonalize() code written and put into the MPS/IQMPS class. I checked and it fixes the bug in the sample code you provided above. It's an entirely new algorithm from the previous .orthogonalize(). The new one is simpler in terms of taking fewer steps and is guaranteed to be globally optimal in terms of the truncations it does.

Please let me know if you encounter more related issues and I'll try to be quicker in getting to them.

Best regards,
Miles

Answer 1

Thanks for reporting this. I do find the same bug.

I've been meaning to put in a new algorithm for .orthogonalize() so I am working on that currently. It will fix this bug plus it ought to be faster too.

Comments

Hi Miles,
Sorry for my delayed reply and thanks for the update.

Are you sure you pushed the update? I am seeing the latest commit to https://github.com/ITensor/ITensor/tree/master to be from December 4th.

Anyway, I compared https://github.com/ITensor/ITensor/tree/master to my fork https://github.com/lfrahm/ITensor/tree/cmplmntryoprtrs, which has the old .orthogonalize() and my implementation of the plussers function.

When executing: (I'm seeing there is little point of using code like this, but just as an example since it reproduces the problem I have)

int main(int argc, char* argv[])
    {
    int bS = 10;
    Hubbard sites;
    sites = Hubbard(bS);

    IQMPO H;

    auto T = MPOt<IQTensor>(sites);

    for (int i = 1; i <= bS; i++)
        {
        for (int j = 1; j <= bS; j++)
            {
            for (int k = 1; k <= bS; k++)
                {
                if ((k == i) && (k == j))
                    T.Anc(k) = sites.op("Nup", k);
                else if (k == i)
                    T.Anc(k) = sites.op("Cdagup", k);
                else if (k == j)
                    T.Anc(k) = sites.op("Cup", k);
                else
                    T.Anc(k) = sites.op("Id", k);
                }
            putMPOLinks(T);

            if (!H)
                H = T;
            else
                H.plusEq(T);
            }
        }

     return 0;
}

The latest code needs 18.390 seconds, where the fork needs just 3.821 seconds on my machine (using the debug library and standard Lapack, same options.mk file). I am curious whether you can reproduce this.

I tried to profile the code and I'm seeing the .orthogonalize() needs most of the time in both cases. I hope this helps.

Thanks again and just let me know if you need anything else.

Best,
Lars

---

Hi Lars,
Yes I may have mistakenly not pushed the patch I made after Dec 4. But I just pushed an even newer patch which includes the previous fix for the truncation issue I mentioned, as well as an additional check that uses a theoretical upper bound on the possible MPS/MPO dimension when truncating within the new orthogonalize code.

Can you try this new code out (pushed to master a few minutes ago) and see if it restores your code's speed back to near what it was before?

Best,
Miles

---

Wow, that was fast. Thanks for your prompt reply.

I just tried the patch and it reduced the calculation time of the example above to 14s. So it's still not really close to the old code.

Do you have any other ideas what could cause this?
Lars

---

Hi Lars,
I just tried the code you provided with the newest version and compared to the version from Dec 4 (before the latest two patches). I found that the code was really fast for both library version (took about 1.4 seconds with both versions).

I see above that you said you were timing with the debug build of the library? If so, does the speed difference persist if you time with the optimized build? It's possible that the debug build speed could be much more sensitive to small changes in the code than the optimized build.

---

Yeah, I am also seeing the slowdown when using the optimized build. So the version you compared to is the one with the old .orthogonalize() implementation but with the changes you had to do to make plusEq working for IQMPOs?

Would you be so kind and compare to my fork at https://github.com/lfrahm/ITensor/tree/cmplmntryoprtrs ? Of course you might have to increase the bS value in the code to make the time representative.

If you don't see a difference there either it must be something at my end.
Thanks,
Lars