PartialTrace
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.
Contents
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.
%% PARTIALTRACE Computes the partial trace of a matrix
% This function has one required argument:
% X: a square matrix
%
% XPT = PartialTrace(X) is the partial trace of the matrix X,
% where it is assumed that length(X) is a perfect squares and both
% subsystems have equal dimension. The trace is taken over the second
% subsystem.
%
% This function has three optional arguments:
% SYS (default 2)
% DIM (default has all subsystems of equal dimension)
% MODE (default -1)
%
% XPT = PartialTrace(X,SYS,DIM,MODE) gives the partial trace of the matrix X,
% where the dimensions of the (possibly more than 2) subsystems are given
% by the vector DIM and the subsystems to take the trace on are given by
% the scalar or vector SYS. MODE is 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 one specific algorithm, set either MODE = 0
% (which generally works best for full or non-numeric matrices, or sparse
% matrices when most of the subsystems are being traced out) or MODE = 1
% (which generally works best when X is large and sparse, and the partial
% trace of X will also be large).
%
% URL: http://www.qetlab.com/PartialTrace
% requires: opt_args.m, PermuteSystems.m
% author: Nathaniel Johnston (nathaniel@njohnston.ca)
% package: QETLAB
% last updated: November 22, 2012
function Xpt = PartialTrace(X,varargin)
lX = length(X);
% set optional argument defaults: sys=2, dim=round(sqrt(length(X))), mode=-1
[sys,dim,mode] = opt_args({ 2, round(sqrt(lX)), -1 },varargin{:});
num_sys = length(dim);
% allow the user to enter a single number for dim
if(num_sys == 1)
dim = [dim,lX/dim];
if abs(dim(2) - round(dim(2))) >= 2*lX*eps
error('PartialTrace:InvalidDim','If DIM is a scalar, DIM must evenly divide length(X).');
end
dim(2) = round(dim(2));
num_sys = 2;
end
sp = issparse(X);
isnum = isnumeric(X);
prod_dim = prod(dim);
prod_dim_sys = prod(dim(sys));
% Determine which of two computation methods to use (i.e., guess which
% method will be faster).
if(mode == -1)
mode = (isnum && sp && prod_dim_sys^2 <= prod_dim);
end
% If the matrix is sparse and the amount we are tracing over is smaller
% than the amount remaining, just do the naive thing and manually add up
% the blocks.
if(mode)
sub_sys_vec = prod_dim*ones(1,prod_dim_sys)/prod_dim_sys;
perm = [sys,setdiff(1:num_sys,sys)];
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
Xpt = sparse(sub_sys_vec(1),sub_sys_vec(1));
for j = 1:prod_dim_sys
Xpt = Xpt + X{j,j};
end
% Otherwise, do a clever trick with mat2cell or reshaping, which is almost always faster.
else
sub_prod = prod_dim/prod_dim_sys;
sub_sys_vec = prod_dim*ones(1,sub_prod)/sub_prod;
perm = [setdiff(1:num_sys,sys),sys];
Xpt = PermuteSystems(X,perm,dim); % permute the subsystems so that we just have to do the partial trace on the second (potentially larger) subsystem
if(isnum) % if the input is a numeric matrix, perform the partial trace operation the fastest way we know how
Xpt = cellfun(@(x) full(trace(x)), mat2cell(Xpt, sub_sys_vec, sub_sys_vec)); % partial trace on second subsystem
if(sp) % if input was sparse, output should be too
Xpt = sparse(Xpt);
end
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)
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]);
Xpt = sum(Xpt(:,:,1:sub_sys_vec(1)+1:sub_sys_vec(1)^2),3);
end
end