PartialTrace

From QETLAB
Jump to: navigation, search
PartialTrace
Computes the partial trace of a matrix

Other toolboxes required none
Related functions PartialMap
PartialTranspose
Function category 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