1
$\begingroup$

If I start with some 4x4 input matrix

$\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.

$\endgroup$
  • $\begingroup$ "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. $\endgroup$ – Henrik Schumacher Apr 27 '18 at 10:00
  • $\begingroup$ 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 $\endgroup$ – JJH Apr 27 '18 at 10:21
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Apr 27 '18 at 10:24
  • $\begingroup$ 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. $\endgroup$ – JJH Apr 27 '18 at 10:35
  • $\begingroup$ Map[m \[Function] m.ConjugateTranspose[m], rhoout] $\endgroup$ – Henrik Schumacher Apr 27 '18 at 10:38
3
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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},
   {k, l} = MapThread[
       (* 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]

product

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

transposed

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.