Single site SVD - QR decomposition

Question

Hello,

for many different reasons, i.e. 1-site TDVP algorithm, it is useful to QR decompose the MPS tensor A=QR of at position b. Since QR is not implemented in ITensor one may still use SVD, namely decompose psi.A(b) into USV and identify Q=U and SV= R. However most of the examples in ITensor are for two-site svd decomposition. I tried to implement it for one site as follow

  auto ir1 = commonIndex(psi.A(b-1),psi.A(b),Link);
  auto ir2 = findtype(psi.A(b) , Site);
  ITensor U(ir1,ir2),V,S;
  svd(phi,U,S,V);

However when giving

    psi.setA(b,U);

a segmentation fault is produced. The reason is that the r=3 tensor U has not its indices ordered as in psi.A(b). U has index (Link,Site,Link) while psi.A(b) has (Site, Link,Link).
Do you know how to quickly solve the problem? Or maybe there is a better way to do the decomposition. Thanks a lot!
Best,
Jacopo.

Answer 1

Hi Jacopo,

That sounds like a good approach for doing what you want to do. In the code above, is phi a copy of psi.A(b)?

I am not sure why it would be seg faulting. The reason would not be that U and psi.A(b) have different ordering of indices, since ITensor doesn't care what the ordering of the indices are. Could you run your code in debug mode and see if it outputs any more specific errors?

In order to help debug it would be helpful to have a minimal working example of the code that is failing.

Here is a minimal example of how this approach would work (this runs correctly for me):

auto N = 10;
auto sites = SpinHalf(N);
auto psi = MPS(sites);
auto b = 4;
// Set the center site gauge with ITensor
psi.position(b);
// Move the center gauge manually
auto lb = linkInd(psi,b);
ITensor U,S,V(lb);
svd(psi.A(b),U,S,V);
psi.setA(b,U);
psi.Aref(b+1) *= S*V;

linkInd(psi,b) is a shorthand for commonIndex(psi.A(b),psi.A(b+1),Link). Also note that if you put indices only on the V tensor, the svd function will use those to determine how to reshape the input tensor.

Note that a more efficient alternative to using svd, if you don't need extremely high accuracy of your singular values, is to use denmatDecomp ( http://itensor.org/docs.cgi?page=classes/decomp ), which directly computes T = U*B where B = S*V.

Cheers,
Matt

Comments

Dear Matt,
thanks a lot, the code you wrote works indeed very well. I guess the error I was getting was because I was imposing the index lb on the matrix U instead of the matrix V (and yes, phi was a copy of psi.A(b) ). Thanks a lot for your time and apologies for the numerous questions!
Jacopo.

---

Glad to hear it works!

It would also be fine specifying the indices of `U`, I was just showing a slightly simpler way since there is only one index that goes onto `V`. Perhaps there was another small bug in your code that rewriting managed to fix.

Cheers,
Matt