# trace and partial trace of MPO

**URL:** <https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79>\
**Category:** ITensor Julia Questions\
**Tags:** julia, mpo\
**Created:** [April 24, 2022, 11:28pm UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79 "2022-04-24T23:28:32Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![jisutich](https://avatars.discourse-cdn.com/v4/letter/j/e79b87/32.png) [@jisutich](https://itensor.discourse.group/u/jisutich)\
**Post date:** [April 24, 2022, 11:28pm UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/1 "2022-04-24T23:28:32Z")

</div>

Hi,

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

For the trace part, I tried

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

```

Or

```auto
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

```auto
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

```auto
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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [April 25, 2022, 2:18am UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/2 "2022-04-25T02:18:16Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![jisutich](https://avatars.discourse-cdn.com/v4/letter/j/e79b87/32.png) [@jisutich](https://itensor.discourse.group/u/jisutich)\
**Post date:** [April 25, 2022, 2:26am UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/3 "2022-04-25T02:26:16Z")

</div>

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

```auto
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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [April 25, 2022, 2:54am UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/4 "2022-04-25T02:54:16Z")

</div>

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](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/63f8e2b8f9dc706ae914c224d284700add5efa27.jpeg)

Basically you can make a loop like this:

```auto
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!).

---

<div class="post-metadata">

**Author:** ![cgh](https://avatars.discourse-cdn.com/v4/letter/c/e95f7d/32.png) [@cgh](https://itensor.discourse.group/u/cgh)\
**Post date:** [April 29, 2022, 12:36am UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/5 "2022-04-29T00:36:12Z")

</div>

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

```auto
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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [April 29, 2022, 12:57am UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/6 "2022-04-29T00:57:29Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![cgh](https://avatars.discourse-cdn.com/v4/letter/c/e95f7d/32.png) [@cgh](https://itensor.discourse.group/u/cgh)\
**Post date:** [April 29, 2022, 3:47am UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/7 "2022-04-29T03:47:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jisutich](https://avatars.discourse-cdn.com/v4/letter/j/e79b87/32.png) [@jisutich](https://itensor.discourse.group/u/jisutich)\
**Post date:** [April 29, 2022, 4:03am UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/8 "2022-04-29T04:03:49Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![cgh](https://avatars.discourse-cdn.com/v4/letter/c/e95f7d/32.png) [@cgh](https://itensor.discourse.group/u/cgh)\
**Post date:** [April 29, 2022, 6:15am UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/9 "2022-04-29T06:15:47Z")

</div>

Ok. Thanks for your help.

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [April 29, 2022, 7:35pm UTC](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/10 "2022-04-29T19:35:14Z")

</div>

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:

```auto
@show inds(L)

```

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

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

```