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

Not sure whether you already aware of what causes this problem, but I figured out that if I replace "multSiteOps(A, B)" by "prime(A,Site)*B.mapprime(2,1,Site)" it works.

---

Very helpful to know. The new algorithm will take a while. This workaround makes it sound like the bug may be something really minor, so maybe I can get it fixed quickly before I get the newer code ready.

---

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

---

Hi Miles,

thanks for your feedback and the new .orthogonalize() code. I tried it and you are right it fixes the problem I described above. However, I'm seeing a massive slowdown in problem I use it for (adding IQMPOs). Do you have any idea where this slowdown comes from? Are you sure the new code is more efficient in terms of CPU time?

Best,
Lars

---

Hi Lars,
Glad it fixes the bug. I did briefly also have a problem myself where it was much slower (mainly due to a minor code issue where the Maxm parameter wasn't being passed through correctly). However when I test it now it doesn't seem to have any speed or scaling issues (I just summed 50 random product states together and it was very fast). Did you get the very latest version from Github (Nov 28th version)? If the bug persists, could you send me some sample code that I can use to reproduce the issue?

Thanks,
Miles

---

Thanks for the reply. I think I got the newest version but I will see what happened there on Monday (It's already late here in Germany). Don't invest time in this one. I will get back to you early next week.

Thanks for the great support!
Lars

---

Hi Lars,
I just encountered a massive slowdown when adding MPS with the new code (same code is used underneath for adding MPS and MPO). The culprit in my case was some faulty logic in our svd routine that was not truncating correctly in the case that all the singular values were exactly zero. I just pushed a fix for this.

It's unlikely that it's the same problem you are having, but perhaps it was (it could occur if one of your MPO's was actually zero somehow).

Otherwise please let me know if the bug persists and let me know of your case. I have one other improvement I'm going to make which might fix this kind of bug more generally that I plan to push soon by the way.

Thanks, Miles

---

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