PartialTrace

PartialTrace
Computes the partial trace of a vector or matrix
Other toolboxes required	none
Related functions	PartialMap PartialTranspose
Function category	Superoperators

PartialTrace is a function that computes the partial trace (i.e., reduced density matrix) of a vector or a matrix. The trace may be taken on any subset of the subsystems on which the matrix acts.

Syntax

• XPT = PartialTrace(X)
• XPT = PartialTrace(X,SYS)
• XPT = PartialTrace(X,SYS,DIM)
• XPT = PartialTrace(X,SYS,DIM,MODE)

Argument descriptions

• X: A vector or a square matrix to have its partial trace (i.e., reduced density matrix) returned.
• SYS (optional, default 2): A scalar or vector containing the indices of the subsystems on which the trace is to be applied.
• DIM (optional, by default has all subsystems of equal dimension): A specification of the dimensions of the subsystems that X lives on. DIM can be provided in one of two ways:
  • If DIM is a scalar, it is assumed that X lives on the tensor product of two spaces, the first of which has dimension DIM and the second of which has dimension length(X)/DIM.
  • If $X \in M_{n_1} \otimes \cdots \otimes M_{n_p}$ then DIM should be a row vector containing the dimensions (i.e., DIM = [n_1, ..., n_p]).
• MODE (optional, default -1): A flag that determines which of two algorithms is used to compute the partial trace. Is only used if X is a matrix (not a vector). If MODE = -1 then this script chooses whichever algorithm it thinks will be faster based on the dimensions of the subsystems being traced out and the sparsity of X. If you wish to force the script to use a specific one of the algorithms (not recommended!), they are generally best used in the following situations:
  • MODE = 0: Best used when X is full or non-numeric. Sometimes also appropriate when X is sparse, particularly when XPT will be quite small compared to X (i.e., when most of the subsystems are traced out).
  • MODE = 1: Best used when X is sparse, particularly when XPT will not be significantly smaller than X itself (i.e., when only a few of the subsystems are traced out).

Examples

A bipartite square matrix

By default, the PartialTrace function takes the trace over the second subsystem:

>> X = reshape(1:16,4,4)'

X =

     1     2     3     4
     5     6     7     8
     9    10    11    12
    13    14    15    16

>> PartialTrace(X)

ans =

     7    11
    23    27

By specifying the SYS argument, you can take the trace over the first subsystem instead:

>> PartialTrace(X,1)

ans =

    12    14
    20    22

Taking the trace over both the first and second subsystems results in the standard trace of X:

>> PartialTrace(X,[1,2])

ans =

    34

A large sparse matrix

This function has no trouble dealing with large sparse matrices. The following line of code generates a random sparse matrix in $M_2 \otimes M_2 \otimes M_5 \otimes M_5 \otimes M_{10} \otimes M_{100} \otimes M_2 \otimes M_5$ and traces out the first, third, fourth, fifth, sixth, and eighth subsystems (resulting in an operator living in $M_2 \otimes M_2$) in under 1/2 of a second on a standard desktop computer:

>> PartialTrace(sprand(1000000,1000000,0.000001),[1,3,4,5,6,8],[2,2,5,5,10,100,2,5])

ans =

   (2,1)       0.9638
   (1,2)       0.5716
   (4,3)       0.7358
   (4,4)       0.7300

Notice that, because input to PartialTrace was sparse, so it its output.

Partial trace of a vector

If a pure state vector is provided as the input X, its reduced density matrix (i.e., the partial trace of X*X') will be returned:

>> X = MaxEntangled(2)

X =

    0.7071
         0
         0
    0.7071

>> PartialTrace(X)

ans =

    0.5000         0
         0    0.5000

Applying the PartialTrace function to a vector X is significantly faster and less memory-intensive than applying it to the pure state matrix X*X', especially in high-dimensional systems:

>> X = RandomStateVector(10*20*30);
>> tic;PartialTrace(X*X',3,[10,20,30]);toc
Elapsed time is 2.029789 seconds.

>> tic;PartialTrace(X,3,[10,20,30]);toc
Elapsed time is 0.001600 seconds.

Source code

Click here to view this function's source code on github.