How do you separate product of two ITensor objects if the two tensors do not have any common indices? (Julia)

Question

As an example, suppose the first and last sites of an MPS M are combined into an ITensor as phi = M[1] * M[N]. Then you perform some optimization that changes phi, for instance eigsolve using some other tensor as the "matrix" and phi as the "vector." Now I want to put the updated value of phi back into M[1] and M[N], but this requires me to somehow split the ITensor so that the two resulting ITensor objects have the same index structure as M[1] and M[N], respectively. How do I do this?

In regular DMRG, this is done using the factorize method, but I'm pretty sure this only works when the two things you want to separate have an index in common. For instance, if you do phi = M[1] * M[2], then since M[1] and M[2] share a common index, then the appropriate factorization will "sever" the appropriate bond between the two tensors, thus giving you one tensor with the same index structure as M[1] and one tensor with the same index structure as M[2].

However, when I tried this, it doesn't work. In particular, the code seems to introduce a new index to bind the two tensors together. To show an example using code, I did phi = M[1] * M[2] * M[3] * M[6] * M[7] * M[8] on an 8-site MPS M. Then, after doing some changes to phi using eigsolve, I wanted to separate it into M[1] * M[2] * M[3] and M[6] * M[7] * M[8] (each of which can then be separated using further factorize commands). I tried to do so using the following command:

index_list = Index[]
for j=1:3
  append!(index_list, inds(M[j]))
end
index_set = IndexSet(index_list)
L,R,spec = factorize(phi, index_set; which_decomp = which_decomp,
                                       tags = tags(linkind(M,3)),
                                       kwargs...)

But the result I got was the following:

left side of phi:
IndexSet{5} (dim=2|id=16|"S=1/2,Site,n=1") (dim=2|id=978|"S=1/2,Site,n=2") 
(dim=2|id=311|"S=1/2,Site,n=3") (dim=1|id=305|"Link,l=3") (dim=8|id=372|"Link,l=3") 
right side of phi:
IndexSet{5} (dim=8|id=372|"Link,l=3") (dim=1|id=608|"Link,l=5") 
(dim=2|id=514|"S=1/2,Site,n=6") (dim=2|id=558|"S=1/2,Site,n=7") 
(dim=2|id=515|"S=1/2,Site,n=8")

Notice the new index (dim=8|id=372|"Link,l=3") used to connect the two tensors. I don't want to do this--I just want to separate the two tensors so that they have the same index structure as before. How do I do this? It's possible that factorization isn't the right term to describe this.

Comments

Hi Sujay,
I think I see what you are asking but have a few questions. First of all, if you modify the tensor `phi` going into factorize then it really just may be the case that it cannot be factorized into an outer product again without some loss of accuracy, and possibly a very significant loss. This is just a mathematical issue as I'm sure you understand, that in general most matrices do not admit a rank-1 decomposition (outer product decomposition) unless they have, well, a rank of one (approximately).

So does that property hold of your phi, namely that it does have approximate matrix rank of 1? It certainly does if you form phi by just multiplying M[1] times M[N] and they don't share any indices before you contract them, but you are modifying phi before re-factorizing it correct?

Or are you saying you want to force the decomposition to be rank 1 no matter what error it might incur?

If you are ok with a rank-1 factorization, then just set maxdim=1 in the keyword argument to factorize. There will still be an Index connecting the resulting factor tensors but you can remove it by doing the following trick:

    L,R,spec = factorize(...) #I'm just omitting the arguments for brevity here
    link = commonind(L,R)
    L *= setelt(link=>1)
    R *= setelt(link=>1)

(I didn't test this though so you may have to tweak it a bit.) The `setelt` function makes an ITensor with the specified tensor element set to 1.0 and all other elements set to 0.0 (though here there is only one element if `dim(link)==1`).

Does that answer your question? If so I'll paste it as the answer below but if not we can discuss more.

Best,
Miles

---

Hi Miles, yes, this does answer the question. I was not aware that deleting this bond is essentially forcing the decomposition to be rank-1, so thank you for pointing that out to me. I am currently playing with including these types of optimizations (along with the standard optimization steps) within my own custom version of DMRG. Much remains to be seen about whether these will work as intended, but it is starting to seem as though there is no unproblematic way to implement what I want to implement. Thank you so much!

---

Ok glad it was a helpful reply! Thanks for letting me know.

Answer 1

(Please see answer to question in discussion above.)