partial trace of MPO

[cgh]

Hi miles,
I want to ask more details from the question :https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79

For instance, I have a MPO rho:

n = 6
s = siteinds("S=1/2",n)
rho = randomMPO(s)

I want to get another MPO M with partial trace:

I have tried code in above link but I do not know what I can get through them.
How can I get M?

[miles]

Hi cgh, here is the code for the tr of MPO function inside of ITensor:

function tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair=ts"" => ts"")
  N = length(M)
  L = tr(M[1]; plev=plev, tags=tags)
  for j in 2:N
    L *= M[j]
    L = tr(L; plev=plev, tags=tags)
  end
  return L
end

What you could do to adapt this to your case is the following:

1. change the range in the for loop to only go from 1:2 or 1:NL where NL is the number of tensors you want to trace over on the left
2. multiply L into the next tensor of the MPO, so the one at location NL+1
3. define a new MPO of length (N-NL) and assign the MPO tensors from NL+1:N to it

Do you think that will work for your case?

[cgh]

Yeah. I change the function and it works well. Thank you so much.

[vyason]

@miles How do I change this code if I divide my sites (1 – N) into these parts A B C D E.

Region A: 1–n1
Region B: n1+1 – n2
Region C: n2+1 – n3
Region D: n3+1 – n4
Region E: n4+1 – N

And I want to trace out regions A C and E, and find the reduced MPO for sites B and D only.

Thanks a lot for your code snippets! They are really helpful

[cgh]

I think you can finish this via ITensor, but using partial trace three times.
First tracing out A as

    p1=outer(psi',psi)

    k1=n1
    L1 = tr(p1[1])
    for i = 2:k1
        L1 *= p1[i]
        L1 = tr(L1)
    end

    p1[k1+1]=L1*p1[k1+1]

Do this for C and E similarly.