# How to apply a partial transpose to a list of matrices

$\rho= \begin{pmatrix} 0&0&0&0\\ 0&1&0&0\\ 0&0&0&0\\ 0&0&0&0 \end{pmatrix}$

and evolve this according to the unitary evolution operator

$U=e^{-iHt}$

where $H$ is the Hamiltonian

$H=0.03(\sigma_x\otimes\sigma_x+\sigma_y\otimes\sigma_y+\sigma_z\otimes\sigma_z)$

such that the output matrix is given by

$\rho_{out}(t)=U\rho U^\dagger$

I can plot a table of values of $\rho$ using just

rho={{0,0,0,0},{0,1,0,0},{0,0,0,0},{0,0,0,0}}

H=0.03(KroneckerProduct[PauliMatrix[1],PauliMatrix[1]]
+KroneckerProduct[PauliMatrix[2],PauliMatrix[2]]
+KroneckerProduct[PauliMatrix[3],PauliMatrix[3]])

U[t_]=MatrixExp[-I*H*t]

rhoout[t_]=U[t].rho.U[t]

rhooutlist=Table[rhoout[t],{t,0,500}].


However, this output is not a list of matrices, but one matrix with a list of values for each element. The next thing I want to do is to take the partial transpose of each output matrix, probably using the method in How does one transpose in place the off-diagonal 16 x 16 blocks of a 32 x 32 matrix?.

I'm not sure whether the best way would be to somehow convert the output to a list of matrices or whether to use a different method for the partial trace.

• "However, this output is not a list of matrices, but one matrix with a list of values for each element. " I don't get it. rhooutlist is a list containing 501 $4 \times 4$ matrices. Apr 27, 2018 at 10:00
• Well I mean if I try to do some operation, for example rhoout.ConjugateTranspose[rhoout] it sees rhoout as one matrix that is the wrong shape rather than applying this to each matrix in the list
– JJH
Apr 27, 2018 at 10:21
• The simplest way to deal with that would be to learn about Map. Define a function that processes a $4 \times 4$ matrix and Map it over the list rhooutlist. Apr 27, 2018 at 10:24
• I'm still not sure how this would apply to the operation rhoout.ConjugateTranspose[rhoout]. The description of how Map works seems to be for operations that are applied to the list, not where the list is a list of operations.
– JJH
Apr 27, 2018 at 10:35
• Map[m \[Function] m.ConjugateTranspose[m], rhoout] Apr 27, 2018 at 10:38

I am not sure, but maybe this does what you want:

r1 = Map[
m \[Function]
Flatten[Transpose[Partition[m, {2, 2}]], {{1, 3}, {2, 4}}],
rhooutlist
]; // AbsoluteTiming // First


0.056507

For each matrix in the list rhooutlist, this swaps the upper right $2 \times 2$ block with the lower left $2 \times 2$ block.

The following does the same but is almost twice as fast:

r2 = Flatten[
Transpose[
Flatten[
Partition[rhooutlist, {1, 2, 2}],
{{1}, {2}, {3}, {4, 5}, {6}}
],
{1, 3, 2, 4, 5}],
{{1}, {2, 4}, {3, 5}}
];// AbsoluteTiming // First

r1 == r2


0.028758

True

And even faster is using Compile which enables parallelization:

r3 = Compile[{{m, _Complex, 2}},
Partition[Flatten[Transpose[Partition[m, {2, 2}], {3, 1, 2, 4}]], 4],
CompilationTarget -> "WVM",
RuntimeAttributes -> {Listable},
Parallelization -> True
][rhooutlist]; // AbsoluteTiming // First

r1 == r2 == r3


0.007143

True

Henrik's answer is hardcoded for transposing specific parts only, but for anyone working with quantum mechanics it would help a lot to have a general partial transpose available. Here is mine:

First I need two helper functions to map "matrix indexes" to the "indexes in the subsystems" and back:

indexCombinedToIndexSystems::usage = "
indexCombinedToIndexSystems[indexCombined,systemsDimensions]

indexCombined: index of the combined system
systemsDimensions: list of subsystem dimensions

automatically threads over indexCombined, if list

Example for 3 2-state systems:
indexSystemsToIndexCombined[6,{2,2,2}] == {1,1,0}
";
indexCombinedToIndexSystems[indexCombined : {__Integer},
systemsDimensions : {__Integer}] :=
indexCombinedToIndexSystems[#, systemsDimensions] & /@ indexCombined;
indexCombinedToIndexSystems[indexCombined_Integer,
systemsDimensions : {__Integer}] := Module[
{reminder = indexCombined - 1, indexSystem, indexSystems},
indexSystems = Table[
{indexSystem, reminder} = QuotientRemainder[reminder, base];
indexSystem,
{base, Reverse@FoldList[Times][systemsDimensions[[;; -2]]]}
];
Append[indexSystems, reminder] + 1
];
indexSystemsToIndexCombined::usage = "
indexSystemsToIndexCombined[indexSystems, systemsDimensions]

indexSystems: list of indexes in each system
systemsDimensions: list of subsystem dimensions

automatically threads over indexSystems, if list of lists

Example for 3 2-state systems:
indexSystemsToIndexCombined[{1,1,0},{2,2,2}] == 6
";
indexSystemsToIndexCombined[indexSystems : {{__Integer} ...},
systemsDimensions : {__Integer}] /;
Dimensions[indexSystems][[2]] === Length[systemsDimensions] :=
indexSystemsToIndexCombined[#, systemsDimensions] & /@ indexSystems;
indexSystemsToIndexCombined[indexSystems : {__Integer},
systemsDimensions : {__Integer}] /;
Length[indexSystems] === Length[systemsDimensions] :=
Total[Reverse@
FoldList[Times][
systemsDimensions[[;; -2]]] (indexSystems[[;; -2]] - 1)] +
indexSystems[[-1]]


These functions can be used a follows:

In[284]:=
Table[indexCombinedToIndexSystems[i, {2, 2, 2}], {i, 8}]
indexSystemsToIndexCombined[#, {2, 2, 2}] & /@ %

Out[284]=
{{1, 1, 1}, {1, 1, 2}, {1, 2, 1}, {1, 2, 2}, {2, 1, 1}, {2, 1, 2},
{2, 2, 1}, {2, 2, 2}}

Out[285]=
{1, 2, 3, 4, 5, 6, 7, 8}

In[294]:=
indexCombinedToIndexSystems[Range[1, 8], {2, 2, 2}]
indexSystemsToIndexCombined[%, {2, 2, 2}]

Out[294]=
{{1, 1, 1}, {1, 1, 2}, {1, 2, 1}, {1, 2, 2}, {2, 1, 1}, {2, 1, 2},
{2, 2, 1}, {2, 2, 2}}

Out[295]=
{1, 2, 3, 4, 5, 6, 7, 8}


Now it is easy to describe the partial transpose as it swaps the "indexes" in a subspace only:

partialTranspose[a_?MatrixQ, systemsDimensions : {__Integer},
systemsTransposed : {___?BooleanQ}] /;
Times @@ systemsDimensions === Length[a] :=
Table[Module[{k, l},
(* Swap the indexes for transposed subsystems *)
If[#2, Reverse@#1, #1]&,
{  indexCombinedToIndexSystems[{i, j},
systemsDimensions]\[Transpose], systemsTransposed}
]\[Transpose]~indexSystemsToIndexCombined~systemsDimensions;
a[[k, l]]],
{i, Length[a]},
{j, Length[a]}
]

comb = Array[Subscript[a, #1, #2] &, {2, 2}]~KroneckerProduct~
Array[Subscript[b, #1, #2] &, {2, 2}];
MatrixForm[comb]


MatrixForm@partialTranspose[comb, {2, 2}, {True, False}]