SymmetricProjection

From QETLAB
Revision as of 13:27, 6 April 2015 by Nathaniel (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
SymmetricProjection
Produces the projection onto the symmetric subspace

Other toolboxes required none
Related functions AntisymmetricProjection
PermutationOperator
SwapOperator
Function category Permutations and symmetry of subsystems

SymmetricProjection is a function that computes the orthogonal projection onto the symmetric subspace of two or more subsystems. The output of this function is always a sparse matrix.

Syntax

  • PS = SymmetricProjection(DIM)
  • PS = SymmetricProjection(DIM,P)
  • PS = SymmetricProjection(DIM,P,PARTIAL)
  • PS = SymmetricProjection(DIM,P,PARTIAL,MODE)

Argument descriptions

  • DIM: The dimension of each of the subsystems.
  • P (optional, default 2): The number of subsystems.
  • PARTIAL (optional, default 0): If PARTIAL = 1 then PS isn't the orthogonal projection itself, but rather a matrix whose columns form an orthonormal basis for the symmetric subspace (and hence PS*PS' is the orthogonal projection onto the symmetric subspace).
  • MODE (optional, default -1): A flag that determines which of two algorithms is used to compute the symmetric projection. If MODE = -1 then this script chooses which algorithm it thinks will be faster based on the values of DIM and P. If you wish to force the script to use a specific one of the algorithms (not recommended!), they are as follows:
    • MODE = 0: Computes the symmetric projection by explicitly constructing an orthonormal basis of the symmetric subspace. The details of how to construct such a basis can be found in [1]. This method is typically fast when DIM is small compared to P.
    • MODE = 1: Computes the symmetric projection by averaging all P! permutation operators (in the sense of the PermutationOperator function). Because P! grows very quickly, this method is only practical when P is small.

Examples

Two subsystems

To compute the symmetric projection on two-qubit space, the following code suffices:

>> SymmetricProjection(2)
 
ans =
 
   (1,1)       1.0000
   (2,2)       0.5000
   (3,2)       0.5000
   (2,3)       0.5000
   (3,3)       0.5000
   (4,4)       1.0000

Note that the output of this function is always sparse. If you want a full matrix (not recommended for even moderately large DIM or P), you must explicitly convert it (as in the following example).

Three subsystems

To compute a matrix whose columns form an orthonormal basis for the symmetric subspace of three-qubit space, set PARTIAL = 1:

>> PS = full(SymmetricProjection(2,3,1))
 
PS =
 
    1.0000         0         0         0
         0         0   -0.5774         0
         0         0   -0.5774         0
         0         0         0   -0.5774
         0         0   -0.5774         0
         0         0         0   -0.5774
         0         0         0   -0.5774
         0    1.0000         0         0

Note that PS is an isometry from the symmetric subspace to the full three-qubit space. In other words, PS'*PS is the identity matrix and PS*PS' is the orthogonal projection onto the symmetric subspace, which we can verify as follows:

>> PS'*PS
 
ans =
 
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1
 
>> PS*PS'
 
ans =
 
    1.0000         0         0         0         0         0         0         0
         0    0.3333    0.3333         0    0.3333         0         0         0
         0    0.3333    0.3333         0    0.3333         0         0         0
         0         0         0    0.3333         0    0.3333    0.3333         0
         0    0.3333    0.3333         0    0.3333         0         0         0
         0         0         0    0.3333         0    0.3333    0.3333         0
         0         0         0    0.3333         0    0.3333    0.3333         0
         0         0         0         0         0         0         0    1.0000

Source code

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

  1. %%  SYMMETRICPROJECTION    Produces the projection onto the symmetric subspace
  2. %   This function has one required argument:
  3. %     DIM: the dimension of the local systems
  4. %
  5. %   PS = SymmetricProjection(DIM) is the orthogonal projection onto the
  6. %   symmetric subspace of two copies of DIM-dimensional space. PS is always
  7. %   a sparse matrix.
  8. %
  9. %   This function has five optional arguments:
  10. %     P (default 2)
  11. %     PARTIAL (default 0)
  12. %     MODE (default -1)
  13. %
  14. %   PS = SymmetricProjection(DIM,P,PARTIAL,MODE) is the orthogonal
  15. %   projection onto the symmetric subspace of P copies of DIM-dimensional
  16. %   space. If PARTIAL = 1 then PS isn't the orthogonal projection itself,
  17. %   but rather a matrix whose columns form an orthonormal basis for the
  18. %   symmetric subspace (and hence PS*PS' is the orthogonal projection onto
  19. %   the symmetric subspace). MODE is a flag that determines which of two
  20. %   algorithms is used to compute the symmetric projection. If MODE = -1
  21. %   then this script chooses whichever algorithm it thinks will be faster
  22. %   based on the values of DIM and P. If you wish to force one specific
  23. %   algorithm, set either MODE = 0 (which generally works best when DIM is
  24. %   small) or MODE = 1 (which generally works best when P is small).
  25. %
  26. %   URL: http://www.qetlab.com/SymmetricProjection
  27.  
  28. %   requires: opt_args.m, PermutationOperator.m, PermuteSystems.m, sporth.m
  29. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  30. %   package: QETLAB
  31. %   last updated: December 15, 2014
  32.  
  33. function PS = SymmetricProjection(dim,varargin)
  34.  
  35.     % set optional argument defaults: p=2, partial=0, mode=-1
  36.     [p,partial,mode] = opt_args({ 2, 0, -1 },varargin{:});
  37.     dimp = dim^p;
  38.  
  39.     if(p == 1)
  40.         PS = speye(dim);
  41.         return
  42.     elseif(mode == -1)
  43.         mode = (dim >= p-1);
  44.     end
  45.  
  46.     % The symmetric projection is the normalization of the sum of all
  47.     % permutation operators. This method of computing the symmetric
  48.     % projection is reasonably quick when the local dimension is large
  49.     % compared to the number of spaces.
  50.     if(mode == 1)
  51.         plist = perms(1:p);
  52.         pfac = factorial(p);
  53.         PS = sparse(dimp,dimp);
  54.  
  55.         for j = 1:pfac
  56.             PS = PS + PermutationOperator(dim*ones(1,p),plist(j,:),0,1);
  57.         end
  58.         PS = PS/pfac;
  59.  
  60.         if(partial)
  61.             PS = sporth(PS);
  62.         end
  63.  
  64.     % When the local dimension is small compared to the number of spaces,
  65.     % it is quicker to just explicitly construct an orthonormal basis for
  66.     % the symmetric subspace.
  67.     else
  68.         lim = nchoosek(dim+p-1,dim-1);
  69.         PS = spalloc(dimp,lim,dimp); % allocate dim^p non-zero elements to PS for speed/memory reasons
  70.         slist = sum_vector(p,dim);
  71.         Id = speye(dim);
  72.  
  73.         for j = 1:lim
  74.             ind = cell2mat(arrayfun(@(x,y) repmat(y,1,x), slist(j,:),1:dim,'un',0));
  75.             plist = unique_perms(ind);
  76.             v = spalloc(dimp,1,nchoosek(p,floor(p/2)));
  77.             sp = size(plist,1);
  78.             for k = 1:sp
  79.                 vt = Id(:,plist(k,1));
  80.                 for m = 2:p
  81.                     vt = kron(vt,Id(:,plist(k,m)));
  82.                 end
  83.                 v = v + vt/sqrt(sp);
  84.             end
  85.             PS(:,j) = v;
  86.         end
  87.  
  88.         if(~partial)
  89.             PS = PS*PS';
  90.         end
  91.     end
  92. end
  93.  
  94. % We need some helper functions to help us through the MODE = 0 algorithm.
  95. function slist = sum_vector(dim,p)
  96.     if p <= 1
  97.         slist = dim;
  98.     else
  99.         k = 0;
  100.         slist = zeros(nchoosek(dim+p-1,p-1),p);
  101.         for j = 0:dim
  102.             cs = nchoosek(j+p-2,p-2);
  103.             t = [sum_vector(j,p-1),(dim-j)*ones(cs,1)];
  104.             slist(k+1:k+cs,:) = t;
  105.             k = k + cs;
  106.         end
  107.     end
  108. end

References

  1. John Watrous. Lecture 21: The quantum de Finetti theorem, Theory of Quantum Information Lecture Notes, 2008.