trace and partial trace of MPO

ITensor Julia Questions
,
1

Hi,

Suppose I have a MPO, how to get the trace/ partial trace of it?

For the trace part, I tried

n = 10
s = siteinds("Qubit",n)
A = randomMPO(s)
A_trace = tr(A)

Or

n = 10
s = siteinds("Qubit",n)
A = randomMPO(s)
i = MPO(s,"Id")
A_trace=inner(A,i)

These seem to work and give me very close result. But when I tried

n = 10
s = siteinds("Qubit",n)
A = randomMPO(s)
for i = 1:n
       A *= delta(s[i],s[i]')
end

like the MPS case, I failed.

Also, for the partial trace part, I tried

n = 10
j = 6
s = siteinds("Qubit",n)
A = randomMPO(s)
for i = j:n
       A *= delta(s[i],s[i]')
end

I also failed. Then what is the correct way to get the partial trace?

Thanks

2

It is a good question. We should add an option to the tr(::MPO) function to just do a partial trace; that would be a good expansion of the uses of that function.

For now, you are very close to a solution with your codes that use the delta tensor approach. But you cannot just contract a whole MPO with a delta tensor like that. To do the partial trace, you need to contract the individual MPO tensors such as A[j] with delta(s[j],s[j]') and then contract the results of those traces with the neighboring MPO tensors so that you will get a smaller MPO. There also needs to be some code logic to define a smaller MPO with fewer tensors to hold the results, and assign the results of the partial trace to that smaller MPO. I think if you diagram out the operation of a partial trace you will see what I mean.

It would be easier to do in the case where the sites you want to trace over are either at the beginning (1:n for some n) or a the end (m:length(A)). Is your application one of these cases?

1 Like
3

Hi Miles, thanks a lot for the answer. For the partial trace, I guess what you mean is something like

n = 10
s = siteinds("Qubit",n)
rho = randomMPS(s)
for i = 1:k  #k is the site I stopped. or i = m:length(rho)
rho[i] = rho[i]*delta(s[i],s[i]')
end

right? But I am not sure how to transfer this into a smaller MPO. This MPO still has n sites. Can you give more hints?

Yes, the sites I want to trace out are either 1:n or m:length(rho). The question I am considering these days is calculate the Renyi entropy or entanglement negativity or a quantum circuit with both haar random evolution and quantum channels. Due to the existence of quantum channels, I need a mixed state in general. But basically I just want to trace out either left or right half of the sites.

Thanks

4

Yes you’re very close to getting the right code written for it. I was about to take some time making that improvement to the tr function earlier tonight, but then I realized it might take some time to really do properly so thought I’d just answer with a sketch of the idea.

Here is a diagram showing the remaining steps you should do after your code above:

IMG_40AD7BE68506-1
IMG_40AD7BE68506-12271×1000 86.6 KB

Basically you can make a loop like this:

L = ITensor(1.0)
for i = 1:k
  L *= rho[i]
end

then do rho[k+1] *= L to ‘absorb’ L into the (k+1)'st MPO tensor. After doing that, rho will no longer be a proper MPO at least in terms of what ITensor expects MPOs to be like. So it would be best to finally define a new MPO like M = MPO(n-k) and assign M[i] = rho[k+i] for i=1:(n-k) (though please check that I have my loop logic right there!).

4 Likes
5

Hi miles, If the above rho is a MPO which is from
rho=outer(psi',psi)
and the psi is a MPS.
I used

L = ITensor(1.0)
for i = 1:k
 L *= rho[i]
end
rho[k+1] *= L 
M = MPO(n-k)
 for i=1:(n-k)
   M[i] = rho[k+i]
 end

I can get M and I can’t make sure whether the M is reduced density matrix from rho.
Thanks

6

Hi, thanks for your comment, but I am not sure if you are asking a question and, if so, what is the question?

7

Can I get after partial trace MPO ‘M’ from MPO ‘rho’ using above code?

8

I think your question is whether the MPO obtained by your code is the correct partial density operator. Personally I think the answer is yes.

9

Ok. Thanks for your help.

10

Hi cgh: I don’t see where in your code you are actually tracing over the rho[i] tensors for i=1:k. If you just contract them with L then the result will be a tensor with many indices.

For creating codes like this, a very useful thing to do is to print out the indices of the tensors at each step, so inside your for loop put:

@show inds(L)

If you want to slow down the output so you can see it at each step, you can do:

import ITensors: pause
...
for i = 1:k
 L *= rho[i]
 @show inds(L)
  pause()
end
1 Like