Function Repository Resource:

MatrixPartialTrace

Source Notebook

Calculate a partial trace of a matrix

Contributed by: Jaroslav Kysela

ResourceFunction["MatrixPartialTrace"][mat,n,{d1,d2,}]

calculates the partial trace of matrix mat over the nth subspace, where mat is assumed to lie in a space constructed as a tensor product of subspaces with dimensions {d1,d2,}.

ResourceFunction["MatrixPartialTrace"][mat,n,d]

calculates the partial trace of mat over the nth subspace, where all subspaces have dimension d.

ResourceFunction["MatrixPartialTrace"][mat,{n1,n2,},dims]

calculates the partial trace of mat over subspaces in positions {n1,n2,}.

ResourceFunction["MatrixPartialTrace"][mat,m;;n,dims]

calculates the partial trace of mat over subspaces in consecutive positions {m,m+1,,n}.

ResourceFunction["MatrixPartialTrace"][mat,Except[{n1,n2,}],dims]

calculates the partial trace of mat over all subspaces except those in positions {n1,n2,}.

ResourceFunction["MatrixPartialTrace"][mat,All,dims]

calculates the partial trace of mat over all subspaces.

Details and Options

The partial trace is an operation that is widely used in quantum theory. The state of a subsystem can be retrieved from the state of the composite system by taking an appropriate partial trace.
The partial trace acts on square matrices and returns square matrices.
ResourceFunction["MatrixPartialTrace"] always returns a two-dimensional array, even in the case when the trace is taken over all subspaces.
ResourceFunction["MatrixPartialTrace"] works for numeric as well as symbolic matrices.
ResourceFunction["MatrixPartialTrace"] works also for SparseArray objects.
The following options can be given:
MethodAutomaticmethod to use for the calculation of the partial trace
"Verbose"Falsewhether to print an additional summary of preprocessed input parameters
The Method option can be set to one of the following values:
"TensorContract"use built-in function TensorContract to calculate the partial trace; usually faster for numeric matrices
"Sum"use built-in function Sum to calculate the partial trace; usually faster for symbolic matrices
Automaticuse "TensorContract" when the input matrix is numeric, and "Sum" when the matrix is symbolic
In the special case of an empty list in the second argument, the trace is not evaluated and the input matrix is returned unaltered.
In the special case of a trace over all subspaces, the function Tr is used internally for both numeric and symbolic matrices.

Examples

Basic Examples (3) 

Consider a 4×4 matrix, which can be interpreted as an element of space , i.e. of the tensor product of two subspaces :

In[1]:=
MatrixForm[mat = Array[Subscript[a, ##] &, {4, 4}]]
Out[1]=

Calculate the partial trace over the second subspace :

In[2]:=
ResourceFunction["MatrixPartialTrace"][mat, {2}, 2] // MatrixForm
Out[2]=

Calculate a partial trace over multiple subspaces at once:

In[3]:=
MatrixForm[mat = Array[Subscript[a, ##] &, {8, 8}]]
Out[3]=
In[4]:=
ResourceFunction["MatrixPartialTrace"][mat, {2, 3}, 2] // MatrixForm
Out[4]=

Calculate a partial trace of a numeric matrix:

In[5]:=
mat = RandomComplex[1 + I, {12, 12}];
ResourceFunction["MatrixPartialTrace"][
  mat, {3}, {2, 2, 3}] // MatrixForm
Out[5]=

Scope (10) 

Data types (3) 

Calculate a partial trace of a numeric matrix:

In[6]:=
mat = RandomComplex[1 + I, {8, 8}];
ResourceFunction["MatrixPartialTrace"][
  mat, {3}, {2, 2, 2}] // MatrixForm
Out[6]=

Calculate a partial trace of a symbolic matrix:

In[7]:=
mat = Array[a, {8, 8}];
ResourceFunction["MatrixPartialTrace"][
  mat, {3}, {2, 2, 2}] // MatrixForm
Out[7]=

Calculate a partial trace of a SparseArray:

In[8]:=
mat = SparseArray@Array[a, {8, 8}];
ResourceFunction["MatrixPartialTrace"][mat, {3}, {2, 2, 2}]
Out[8]=

If the input is a SparseArray, so is the output. This is true even if the resulting array is just a 1×1 array:

In[9]:=
mat = SparseArray@Array[Subscript[a, ##] &, {12, 12}];
rmat = ResourceFunction["MatrixPartialTrace"][mat, All, {2, 3, 2}]
Out[9]=
In[10]:=
Normal[rmat]
Out[10]=

Subspace specification (5) 

The second argument can contain negative numbers. In that case, the counting goes from the end of the list, in complete analogy to Take specification:

In[11]:=
mat = Array[Subscript[a, ##] &, {8, 8}];
ResourceFunction["MatrixPartialTrace"][mat, {1, -1}, 2] // MatrixForm
Out[11]=

If the partial trace is taken only over one subspace, the second argument can be specified as the index of this subspace:

In[12]:=
mat = Array[Subscript[a, ##] &, {8, 8}];
ResourceFunction["MatrixPartialTrace"][mat, 3, 2] // MatrixForm
Out[12]=

If the partial trace is to be taken over a range {a,…,b} of successive subspaces, the second argument can be specified as a;;b:

In[13]:=
mat = Array[Subscript[a, ##] &, {16, 16}];
ResourceFunction["MatrixPartialTrace"][mat, 1 ;; 3, 2] // MatrixForm
Out[13]=

If the partial trace is to be taken over all but one subspace with index k, the second argument can also be specified as Except[k]:

In[14]:=
mat = Array[Subscript[a, ##] &, {16, 16}];
ResourceFunction["MatrixPartialTrace"][mat, Except[4], 2] // MatrixForm
Out[14]=

More indices can be specified:

In[15]:=
mat = Array[Subscript[a, ##] &, {16, 16}];
ResourceFunction["MatrixPartialTrace"][mat, Except[{1, 4}], 2] // MatrixForm
Out[15]=

When the second argument is set to All, the partial trace is taken over all subspaces, which corresponds to the standard trace:

In[16]:=
mat = Array[Subscript[a, ##] &, {16, 16}];
ResourceFunction["MatrixPartialTrace"][mat, All, 2] // MatrixForm
Out[16]=
In[17]:=
Tr[mat]
Out[17]=

Dimension specification (2) 

If all subspaces have the same dimension, only this one dimension can be entered instead of the whole list:

In[18]:=
mat = Array[Subscript[a, ##] &, {8, 8}];
rmat1 = ResourceFunction["MatrixPartialTrace"][mat, {1, 3}, {2, 2, 2}];
rmat2 = ResourceFunction["MatrixPartialTrace"][mat, {1, 3}, 2];
In[19]:=
rmat1 === rmat2
Out[19]=
In[20]:=
MatrixForm[rmat1]
Out[20]=

Dimensions can be different for different subspaces. However, their product has to be equal to the dimension of the input matrix:

In[21]:=
mat = Array[Subscript[a, ##] &, {2 3 4, 2 3 4}];
rmat = ResourceFunction["MatrixPartialTrace"][mat, {1, 3}, {2, 3, 4}];
In[22]:=
MatrixForm[rmat]
Out[22]=

Options (6) 

Method (3) 

Calculate the partial trace using Mathematica’s built-in function TensorContract. This routine is usually faster for numeric matrices:

In[23]:=
mat = Array[Subscript[a, ##] &, {8, 8}];
ResourceFunction["MatrixPartialTrace"][mat, 2, 2, Method -> "TensorContract"] // MatrixForm
Out[23]=

Calculate the partial trace using block-wise summation implemented by Mathematica’s built-in function Sum. This routine is usually faster for symbolic matrices:

In[24]:=
mat = Array[Subscript[a, ##] &, {8, 8}];
ResourceFunction["MatrixPartialTrace"][mat, 2, 2, Method -> "Sum"] // MatrixForm
Out[24]=

By default, the "TensorContract" method is used, when the input matrix is numeric, and the "Sum" method is used when it is symbolic:

In[25]:=
mat = Array[Subscript[a, ##] &, {8, 8}];
ResourceFunction["MatrixPartialTrace"][mat, 2, 2, Method -> Automatic] // MatrixForm
Out[25]=

Compare the speed of different methods for numeric and symbolic matrices:

In[26]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/4413703f-d8a6-4492-b85b-5fd5be2b152f"]
In[27]:=
mat1 == mat2 == mat3
Out[27]=
In[28]:=
{time1, time2, time3}
Out[28]=
In[29]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/05b121a7-fd50-4424-a57a-8ae8addf8fdf"]
In[30]:=
mat1 === mat2 === mat3
Out[30]=
In[31]:=
{time1, time2, time3}
Out[31]=

There are two special cases that do not follow the previously stated rules. The first is when the list of subspaces to be traced over is empty. In that case, the input matrix is returned:

In[32]:=
mat = Array[a, {16, 16}];
rmat = ResourceFunction["MatrixPartialTrace"][mat, {}, 2];
In[33]:=
mat === rmat
Out[33]=

The second case is when the partial trace is to be taken over all subspaces. In that case, Mathematica’s built-in function Tr is used internally instead:

In[34]:=
mat = Array[Subscript[a, ##] &, {16, 16}];
ResourceFunction["MatrixPartialTrace"][mat, All, 2]
Out[34]=
In[35]:=
Tr[mat]
Out[35]=

Verbose (3) 

When set to True, there is an additional summary of preprocessed values of parameters printed:

In[36]:=
mat = RandomComplex[1 + I, {16, 16}];
ResourceFunction["MatrixPartialTrace"][mat, {-3, 4}, 2, "Verbose" -> True]
Out[36]=
In[37]:=
mat = Array[Subscript[a, ##] &, {16, 16}];
ResourceFunction["MatrixPartialTrace"][mat, Except[2], 2, "Verbose" -> True]
Out[37]=

When the second argument is an empty list, the "method" item returns None:

In[38]:=
mat = RandomComplex[1 + I, {16, 16}];
ResourceFunction["MatrixPartialTrace"][mat, {}, {4, 4}, "Verbose" -> True];

When the second argument is the list of all subspaces, the routine internally uses Tr and the "method" item returns "Tr":

In[39]:=
mat = Array[Subscript[a, ##] &, {16, 16}];
ResourceFunction["MatrixPartialTrace"][mat, 1 ;; 4, 2, "Verbose" -> True]
Out[39]=

Applications (2) 

In quantum mechanics, a partial trace of a density matrix gives the quantum state of a subsystem.

Let ρ12 be a density matrix of a system of two maximally entangled photons in their polarization:

In[40]:=
Subscript[\[Rho], 12] = {{1/2, 0, 0, -(1/2)}, {0, 0, 0, 0}, {0, 0, 0, 0}, {-(1/2), 0, 0, 1/2}};

The quantum state of one of the photons is obtained as a partial trace over the state of the other photon. In this case, the resulting state ρ1 is maximally mixed:

In[41]:=
Subscript[\[Rho], 1] = ResourceFunction["MatrixPartialTrace"][Subscript[\[Rho], 12], 2, 2]
Out[41]=

The same result is obtained when tracing over the first photon to get the state ρ2 of the second photon:

In[42]:=
Subscript[\[Rho], 2] = ResourceFunction["MatrixPartialTrace"][Subscript[\[Rho], 12], 1, 2]
Out[42]=

Generate a random quantum state of three spin-1/2 particles:

In[43]:=
mat = RandomComplex[1 + I, {8, 8}];
mat = mat\[ConjugateTranspose] . mat;
Subscript[\[Rho], 123] = mat/Tr[mat];

Matrix ρ123 is a positive semidefinite matrix with unit trace and thus represents a valid density operator:

In[44]:=
PositiveSemidefiniteMatrixQ[Subscript[\[Rho], 123]]
Out[44]=
In[45]:=
Tr[Subscript[\[Rho], 123]]
Out[45]=

The state of the second particle is obtained as the partial trace of ρ123 over the first and third particles:

In[46]:=
Subscript[\[Rho], 2] = ResourceFunction["MatrixPartialTrace"][Subscript[\[Rho], 123], {1, 3}, 2];
MatrixForm[Subscript[\[Rho], 2]]
Out[46]=

This matrix again represents a valid density operator:

In[47]:=
PositiveSemidefiniteMatrixQ[Subscript[\[Rho], 2]]
Out[47]=
In[48]:=
Tr[Subscript[\[Rho], 2]]
Out[48]=

Properties and Relations (3) 

The order in which the subspaces are specified is irrelevant:

In[49]:=
mat = Array[a, {32, 32}];
rmat1 = ResourceFunction["MatrixPartialTrace"][mat, {1, 3, 5}, 2];
rmat2 = ResourceFunction["MatrixPartialTrace"][mat, {5, 3, 1}, 2];
In[50]:=
rmat1 === rmat2
Out[50]=

Successive applications of a partial trace for different subspaces gives the same results as a single application of a multi-subspace partial trace:

In[51]:=
mat = Array[Subscript[a, ##] &, {12, 12}];
In[52]:=
rmat1 = ResourceFunction["MatrixPartialTrace"][mat, {3}, {2, 3, 2}];
rmat2 = ResourceFunction["MatrixPartialTrace"][rmat1, {2}, {2, 3}];
MatrixForm[rmat2]
Out[52]=
In[53]:=
rmat12 = ResourceFunction["MatrixPartialTrace"][
   mat, {2, 3}, {2, 3, 2}];
MatrixForm[rmat12]
Out[53]=
In[54]:=
rmat12 === rmat2
Out[54]=

The partial trace over all subspaces corresponds to the standard trace:

In[55]:=
mat = Array[Subscript[a, ##] &, {12, 12}];
In[56]:=
ResourceFunction["MatrixPartialTrace"][mat, {1, 2, 3}, {2, 3, 2}]
Out[56]=
In[57]:=
Tr[mat]
Out[57]=

Possible Issues (5) 

MatrixPartialTrace always returns a matrix, even when the trace is taken over all subspaces:

In[58]:=
mat = Array[Subscript[a, ##] &, {2^3, 2^3}];
tr1 = ResourceFunction["MatrixPartialTrace"][
  mat, {1, 2, 3}, {2, 2, 2}]
Out[58]=

The trace over all subspaces can be calculated also using Tr, which returns the sum of diagonal elements:

In[59]:=
tr2 = Tr[mat]
Out[59]=
In[60]:=
tr1[[1, 1]] == tr2
Out[60]=

The same is true also when the input matrix is a SparseArray:

In[61]:=
mat = SparseArray@Array[Subscript[a, ##] &, {2^3, 2^3}];
In[62]:=
tr1 = ResourceFunction["MatrixPartialTrace"][
  mat, {1, 2, 3}, {2, 2, 2}]
Out[62]=
In[63]:=
tr2 = Tr[mat]
Out[63]=
In[64]:=
tr1[[1, 1]] == tr2
Out[64]=

When both positive and negative indices are used in the second argument, it might happen that multiple indices refer to the same subspace. This case is nevertheless not allowed:

In[65]:=
mat = Array[Subscript[a, ##] &, {12, 12}];
ResourceFunction["MatrixPartialTrace"][mat, {1, -3}, {2, 3, 2}]
Out[65]=

For some settings of the initial index, final index, and the step in the Span specification, undesired behavior may occur:

In[66]:=
mat = Array[a, {32, 32}];
rmat = ResourceFunction["MatrixPartialTrace"][mat, 3 ;; 1, 2];
In[67]:=
rmat === mat
Out[67]=

Use "Verbose"True to see how the input parameters are interpreted:

In[68]:=
mat = Array[a, {32, 32}];
rmat = ResourceFunction["MatrixPartialTrace"][mat, 3 ;; 1, 2, "Verbose" -> True];

When the indices to be excluded from the tracing are out of the valid range, the trace is done over all subspaces:

In[69]:=
mat = Array[a, {32, 32}];
rmat = ResourceFunction["MatrixPartialTrace"][mat, Except[30], 2];
In[70]:=
rmat === {{Tr[mat]}}
Out[70]=

Note that when tracing over multiple subspaces one by one, the position and dimensions of subspaces in each new application of the partial traces may change:

In[71]:=
mat = Array[Subscript[a, ##] &, {12, 12}];
In[72]:=
rmat12 = ResourceFunction["MatrixPartialTrace"][
   mat, {1, 2}, {2, 3, 2}];
MatrixForm[rmat12]
Out[72]=
In[73]:=
rmat1 = ResourceFunction["MatrixPartialTrace"][mat, {1}, {2, 3, 2}];
rmat2 = ResourceFunction["MatrixPartialTrace"][rmat1, {1}, {3, 2}];
MatrixForm[rmat2]
Out[73]=
In[74]:=
rmat12 === rmat2
Out[74]=

When the tracing proceeds from the end of the list of positions, here {1,2}, no recalculation of subspace positions is necessary:

In[75]:=
rmat1 = ResourceFunction["MatrixPartialTrace"][mat, {2}, {2, 3, 2}];
rmat2 = ResourceFunction["MatrixPartialTrace"][rmat1, {1}, {2, 2}];
MatrixForm[rmat2]
Out[75]=
In[76]:=
rmat12 === rmat2
Out[76]=

Neat Examples (1) 

Plot the partial trace over an increasing number of subspaces:

In[77]:=
n = 4;
mat = Table[
   Boole@PossibleZeroQ@
     Total[IntegerDigits[i, 2, n] - IntegerDigits[j, 2, n]], {i, 2^
    n}, {j, 2^n}];
MatrixPlot[mat]
Out[77]=
In[78]:=
MatrixPlot[ResourceFunction["MatrixPartialTrace"][mat, #, 2], PlotLabel -> ("Trace over subspace: " <> ToString[#])] & /@ Range[n]
Out[78]=
In[79]:=
MatrixPlot[ResourceFunction["MatrixPartialTrace"][mat, #, 2], PlotLabel -> ("Trace over subspaces: " <> ToString[#])] & /@ Subsets[Range[n], {2}]
Out[79]=
In[80]:=
MatrixPlot[ResourceFunction["MatrixPartialTrace"][mat, #, 2], PlotLabel -> ("Trace over subspaces: " <> ToString[#])] & /@ Subsets[Range[n], {3}]
Out[80]=

Publisher

Jaroslav Kysela

Version History

  • 1.0.0 – 03 May 2021

License Information