Function Repository Resource:

MatrixPartialTrace

Calculate a partial trace of a matrix

Contributed by: Jaroslav Kysela

ResourceFunction["MatrixPartialTrace"][mat,n,{d₁,d₂,…}]

calculates the partial trace of matrix mat over the n^th subspace, where mat is assumed to lie in a space constructed as a tensor product of subspaces with dimensions {d₁,d₂,…}.

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

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

ResourceFunction["MatrixPartialTrace"][mat,{n₁,n₂,…},dims]

calculates the partial trace of mat over subspaces in positions {n₁,n₂,…}.

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[{n₁,n₂,…}],dims]

calculates the partial trace of mat over all subspaces except those in positions {n₁,n₂,…}.

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:

Method

Automatic

method to use for the calculation of the partial trace

"Verbose"

False

whether 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

Automatic

use "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]:=

Out[1]=

Calculate the partial trace over the second subspace :

In[2]:=

Out[2]=

Calculate a partial trace over multiple subspaces at once:

In[3]:=

Out[3]=

In[4]:=

Out[4]=

Calculate a partial trace of a numeric matrix:

In[5]:=

Out[5]=

Scope (10)

Data types (3)

Calculate a partial trace of a numeric matrix:

In[6]:=

Out[6]=

Calculate a partial trace of a symbolic matrix:

In[7]:=

Out[7]=

Calculate a partial trace of a SparseArray:

In[8]:=

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]:=

Out[9]=

In[10]:=

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]:=

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]:=

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]:=

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]:=

Out[14]=

More indices can be specified:

In[15]:=

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]:=

Out[16]=

In[17]:=

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]:=

Out[19]=

In[20]:=

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]:=

In[22]:=

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]:=

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]:=

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]:=

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]:=

Out[27]=

In[28]:=

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]:=

Out[30]=

In[31]:=

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]:=

In[33]:=

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]:=

Out[34]=

In[35]:=

Out[35]=

Verbose (3)

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

In[36]:=

Out[36]=

In[37]:=

Out[37]=

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

In[38]:=

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

In[39]:=

Out[39]=

Applications (2)

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

Let ρ₁₂ 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 ρ₁ 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 ρ₂ 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 ρ₁₂₃ 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 ρ₁₂₃ 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]:=

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]:=

In[52]:=

rmat1 = ResourceFunction["MatrixPartialTrace"][mat, {3}, {2, 3, 2}];
rmat2 = ResourceFunction["MatrixPartialTrace"][rmat1, {2}, {2, 3}];
MatrixForm[rmat2]

Out[52]=

In[53]:=

Out[53]=

In[54]:=

Out[54]=

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

In[55]:=

In[56]:=

Out[56]=

In[57]:=

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]:=

Out[59]=

In[60]:=

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]:=

Out[62]=

In[63]:=

Out[63]=

In[64]:=

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]:=

Out[65]=

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

In[66]:=

In[67]:=

Out[67]=

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

In[68]:=

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

In[69]:=

In[70]:=

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]:=