# PartialTrace

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
 Other toolboxes required PartialTrace Computes the partial trace of a matrix none PartialMapPartialTranspose Superoperators

PartialTrace is a function that computes the partial trace of 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 matrix to have its partial trace 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. 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.

## Source code

Click on "expand" to the right to view the MATLAB source code for this function.

1. %%  PARTIALTRACE    Computes the partial trace of a matrix
2. %   This function has one required argument:
3. %     X: a square matrix
4. %
5. %   XPT = PartialTrace(X) is the partial trace of the matrix X,
6. %   where it is assumed that length(X) is a perfect squares and both
7. %   subsystems have equal dimension. The trace is taken over the second
8. %   subsystem.
9. %
10. %   This function has three optional arguments:
11. %     SYS (default 2)
12. %     DIM (default has all subsystems of equal dimension)
13. %     MODE (default -1)
14. %
15. %   XPT = PartialTrace(X,SYS,DIM,MODE) gives the partial trace of the matrix X,
16. %   where the dimensions of the (possibly more than 2) subsystems are given
17. %   by the vector DIM and the subsystems to take the trace on are given by
18. %   the scalar or vector SYS. MODE is a flag that determines which of two
19. %   algorithms is used to compute the partial trace. If MODE = -1 then this
20. %   script chooses whichever algorithm it thinks will be faster based on
21. %   the dimensions of the subsystems being traced out and the sparsity of
22. %   X. If you wish to force one specific algorithm, set either MODE = 0
23. %   (which generally works best for full or non-numeric matrices, or sparse
24. %   matrices when most of the subsystems are being traced out) or MODE = 1
25. %   (which generally works best when X is large and sparse, and the partial
26. %   trace of X will also be large).
27. %
28. %   URL: http://www.qetlab.com/PartialTrace
29. 
30. %   requires: opt_args.m, PermuteSystems.m
31. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
32. %   package: QETLAB
33. %   last updated: November 22, 2012
34. 
35. function Xpt = PartialTrace(X,varargin)
36. 
37. lX = length(X);
38. 
39. % set optional argument defaults: sys=2, dim=round(sqrt(length(X))), mode=-1
40. [sys,dim,mode] = opt_args({ 2, round(sqrt(lX)), -1 },varargin{:});
41. 
42. num_sys = length(dim);
43. 
44. % allow the user to enter a single number for dim
45. if(num_sys == 1)
46.     dim = [dim,lX/dim];
47.     if abs(dim(2) - round(dim(2))) >= 2*lX*eps
48.         error('PartialTrace:InvalidDim','If DIM is a scalar, DIM must evenly divide length(X).');
49.     end
50.     dim(2) = round(dim(2));
51.     num_sys = 2;
52. end
53. 
54. sp = issparse(X);
55. isnum = isnumeric(X);
56. prod_dim = prod(dim);
57. prod_dim_sys = prod(dim(sys));
58. 
59. % Determine which of two computation methods to use (i.e., guess which
60. % method will be faster).
61. if(mode == -1)
62.     mode = (isnum && sp && prod_dim_sys^2 <= prod_dim);
63. end
64. 
65. % If the matrix is sparse and the amount we are tracing over is smaller
66. % than the amount remaining, just do the naive thing and manually add up
67. % the blocks.
68. if(mode)
69.     sub_sys_vec = prod_dim*ones(1,prod_dim_sys)/prod_dim_sys;
70.     perm = [sys,setdiff(1:num_sys,sys)];
71. 
72.     X = mat2cell(PermuteSystems(X,perm,dim), sub_sys_vec, sub_sys_vec); % permute the subsystems so that we just have to do the partial trace on the first (potentially larger) subsystem
73.     Xpt = sparse(sub_sys_vec(1),sub_sys_vec(1));
74.     for j = 1:prod_dim_sys
75.         Xpt = Xpt + X{j,j};
76.     end
77. 
78. % Otherwise, do a clever trick with mat2cell or reshaping, which is almost always faster.
79. else
80.     sub_prod = prod_dim/prod_dim_sys;
81.     sub_sys_vec = prod_dim*ones(1,sub_prod)/sub_prod;
82.     perm = [setdiff(1:num_sys,sys),sys];
83. 
84.     Xpt = PermuteSystems(X,perm,dim); % permute the subsystems so that we just have to do the partial trace on the second (potentially larger) subsystem
85. 
86.     if(isnum) % if the input is a numeric matrix, perform the partial trace operation the fastest way we know how
87.         Xpt = cellfun(@(x) full(trace(x)), mat2cell(Xpt, sub_sys_vec, sub_sys_vec)); % partial trace on second subsystem
88.         if(sp) % if input was sparse, output should be too
89.             Xpt = sparse(Xpt);
90.         end
91.     else % if the input is not numeric (such as a variable in a semidefinite program), do a slower method that avoids mat2cell (mat2cell doesn't like non-numeric arrays)
92.         Xpt = reshape(permute(reshape(Xpt,[sub_sys_vec(1),sub_prod,sub_sys_vec(1),sub_prod]),[2,4,1,3]),[sub_prod,sub_prod,sub_sys_vec(1)^2]);
93.         Xpt = sum(Xpt(:,:,1:sub_sys_vec(1)+1:sub_sys_vec(1)^2),3);
94.     end
95. end