# How to compute the partial trace of a 4x4 matrix?

I want to compute a partial trace using formula

$\rho_A=\sum\langle B|\rho_{AB}|B\rangle$ .

Example,

$\rho_{AB}= $$\begin{pmatrix} a & b & c & d\\ e & f & g & h\\ i & j & k & l\\ m & n & o & p \end{pmatrix} I try to evaluate MatrixForm[ Sum[ TensorProduct[{{a, b, c, d}, {e, f, g, h}, {i, j, k, l}, {m, n, o, p}}, (q*r)], {1, 2}]]  but Mathematica returns the expression unevaluated. The output is wrong. I am not sure whether the operator should be TensorProduct or Cross. I expected to get \rho_{A}=$$ \begin{pmatrix} a +f & c+h\\ i+n & k+p\\ \end{pmatrix}$

• Not familiar with partial trace, but does this package help? Jun 23, 2016 at 8:36
• Thank you again. Finally. I know. I explore it first Jun 23, 2016 at 10:17
• Your expression for Sum is ill-formed. Jun 23, 2016 at 13:46
• I will check it .tq] Jun 23, 2016 at 21:00

Your example can be achieved using Map with a level specification, Partition to generate the sub-matrices and Tr to calculate the traces.

ClearAll[a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p]
r = {{a, b, c, d}, {e, f, g, h}, {i, j, k, l}, {m, n, o, p}};

Map[Tr, Partition[r, {2, 2}], {2}]


{{a + f, c + h}, {i + n, k + p}}

• I dont get the same value Jun 23, 2016 at 10:17
• Do you mean that evaluating the expressions above don't give you the same result? Jun 23, 2016 at 10:28
• yes. but i change the r in term in matriks, may you explain the command Map, Tr, Partition work ? tq Jun 23, 2016 at 12:00
• I'm guessing that this Q&A may help you. To learn how this works try experimenting, e.g. evaluate Partition[r, {2, 2}] on its own and perhaps something like Map[dummy, Partition[r, {2, 2}], {2}] and read the documentation for these functions. Jun 23, 2016 at 12:11
• thank you very much MikeLimaOscar Jun 23, 2016 at 20:59

You can achieve this by partitioning your matrix and performing a tensor contraction. To be clear, having: $$\rho_{AB} = \sum_{ijkl} \rho_{ij}^{kl} | ij \rangle \langle kl |$$ the partial trace over $B$ may be defined as: $$\mathrm{tr}_B \rho_{AB} = \sum_{jl} |j\rangle \langle l| \sum_m \rho_{mj}^{ml},$$ which is a tensor contraction. You have written $\rho$ as a matrix (or a two-indexed tensor), so first you need to partition it:

rhop = ArrayReshape[rho, {2, 2, 2, 2}];


Now rhop[[i,j,k,l]] will give you $\rho_{ij}^{kl}$. Then the tensor contraction is easily done by Mathematica's libraries:

rhoc = TensorContract[rhop, {2, 4}]; rhoc//TeXForm


This yields $$\left(\begin{array}{cc} a+f & c+h \\ i+n & k+p \\\end{array}\right).$$ Similarly, if you wanted to take the partial trace over $A$, you just need to contract the first indices:

rhod = TensorContract[rhop, {1, 3}]; rhod // TeXForm


Which returns:

$$\left( \begin{array}{cc} a+k & b+l \\ e+o & f+p \\ \end{array} \right).$$

Note that using Partition yields a slightly different tensor than ArrayReshape, but you can convert one in another by performing a TensorTranspose over the indices 2 and 3.