Combined re-orthogonalization of an MPS and repositioning of a ProjMPO

ITensor Julia Questions
, ,
1

Hi all!

I have an MPS called state and an MPO called H; I’m trying to calculate the product of the projection, using ProjMPO, of H with some tensors obtained from state.
This is partly like a typical TDVP1 sweep, and in fact I based my code on that algorithm.
I need to compute the following.

  1. The product of H’s one-site projection on the first site with the first tensor of the MPS (with its orthogonality center on the first site):
PH = ProjMPO(H)
PH.nsite = 1
orthogonalize!(state, 1)
ITensors.position!(PH, state, 1)
H1 = PH(state[1])
  1. then, I need to decompose the first site of the MPS in a left-orthogonal tensor A_1 and a “residue” C, and multiply this by the zero-site projection of H:
_, S, V = svd(state[1], uniqueinds(state[1], state[2]))
C = S * V
PH.nsite = 0
ITensors.position!(PH, state, 2)
K = PH(C)
  1. The same as point 1, but on the second bond, i.e. the product of H’s one-site projection on the second site with the second tensor of state with the orthogonality center moved on the second bond as well:
PH.nsite = 1
orthogonalize!(state, 2)
ITensors.position!(PH, state, 2)
H2 = PH(state[2]) # <---- error

When I run this code (in sequence) I get the following error, referring to the line I pointed out:

ERROR: LoadError: The order of the ProjMPO-ITensor product P*v is not equal to the order of the ITensor v, this is probably due to an index mismatch.
Common reasons for this error: 
(1) You are trying to multiply the ProjMPO with the 1-site wave-function at the wrong position.
(2) `orthogonalize!` was called, changing the MPS without updating the ProjMPO.

P*v inds: ((dim=1|id=245|"Link,l=1"), (dim=1|id=875|"Link,l=1"), (dim=1|id=875|"Link,l=1")', (dim=4|id=788|"Site,n=2,vS=1/2")', (dim=1|id=898|"Link,l=2")') 

v inds: ((dim=4|id=788|"Site,n=2,vS=1/2"), (dim=1|id=898|"Link,l=2"), (dim=1|id=245|"Link,l=1"))

I can see that there are two spurious indices on P*v, id=875 and its primed version, that weren’t present in v, but I don’t get where they come from.
I did call orthogonalize! on the state, but also I updated the ProjMPO as well. Or so I think. Does my code not do what I think it does?

2

By the way, here’s a minimal (not) working example:

using ITensors

let
    s = siteinds("S=1/2", 6)
    z = randomMPS(s; linkdims=4)
    h = randomMPO(s)

    ph = ProjMPO(h)
    
    ph.nsite = 1
    orthogonalize!(z, 1)
    position!(ph, z, 1)
    ph(z[1])

    ph.nsite = 0
    _, S, V = svd(z[1], uniqueinds(z[1], z[2]))
    position!(ph, z, 2)
    ph(S * V)

    ph.nsite = 1
    orthogonalize!(z, 2)
    position!(ph, z, 2)
    ph(z[2])
end

I find that if I remove the position!(ph, z, 2) then it doesn’t give any error, but I don’t understand why. I thought that I needed to re-project the ProjMPO object on the newly-orthogonalized state.

3

After some trial and error I managed to figure out the culprit. I’ll post my solution.

Clearly the SVD factorization doesn’t affect the MPS z so it can be left out. Actually, if I remove the whole block containing the factorization I encounter no issues, i.e. the following code works:

using ITensors

let
    s = siteinds("S=1/2", 6)
    z = randomMPS(s; linkdims=4)
    h = randomMPO(s)

    ph = ProjMPO(h)
    
    ph.nsite = 1
    orthogonalize!(z, 1)
    position!(ph, z, 1)
    ph(z[1])

    orthogonalize!(z, 2)
    position!(ph, z, 2)
    ph(z[2])
end

The following, however, does not:

using ITensors

let
    s = siteinds("S=1/2", 6)
    z = randomMPS(s; linkdims=4)
    h = randomMPO(s)

    ph = ProjMPO(h)
    
    ph.nsite = 1
    orthogonalize!(z, 1)
    position!(ph, z, 1)
    ph(z[1])

    position!(ph, z, 2) # <-------- !!

    orthogonalize!(z, 2)
    position!(ph, z, 2)
    ph(z[2])
end

By looking at the code for position! I see that there are some checks on the site parameter, such that if it hasn’t changed since the last time the function was called (and ph.nsite is not changed either) then the function does nothing. This means that the second position!(ph, z, 2) call doesn’t change ph at all, even if the MPS z has been changed in the meantime.
As a result, when ph(z[2]) is called, z has a new Index object (i.e. with a different id) on the (1,2) link due to the re-orthogonalization process, while ph is still using the old set of indices from z, causing an index mismatch.

My previous example works if I switch the last two blocks—then position! uses the new information on ph.nsite and the projections are correctly updated after the re-orthogonalization of the MPS.