SymmetricExtension

From QETLAB
Jump to: navigation, search
SymmetricExtension
Determines whether or not an operator has a symmetric extension

Other toolboxes required cvx
Related functions IsPPT
IsSeparable
SymmetricInnerExtension
SymmetricProjection
Function category Entanglement and separability

SymmetricExtension is a function that determines whether or not a given positive semidefinite operator has a symmetric extension. This function is extremely useful for showing that quantum states are entangled (see the Examples section). Various types of symmetric extensions (such as Bosonic and/or PPT extensions) can be looked for by specifying optional arguments in the function.

Syntax

  • EX = SymmetricExtension(X)
  • EX = SymmetricExtension(X,K)
  • EX = SymmetricExtension(X,K,DIM)
  • EX = SymmetricExtension(X,K,DIM,PPT)
  • EX = SymmetricExtension(X,K,DIM,PPT,BOS)
  • EX = SymmetricExtension(X,K,DIM,PPT,BOS,TOL)
  • [EX,WIT] = SymmetricExtension(X,K,DIM,PPT,BOS,TOL)

Argument descriptions

Input arguments

  • X: A positive semidefinite operator.
  • K (optional, default 2): The number of copies of the second subsystem in the desired symmetric extension.
  • DIM (optional, by default has both subsystems of equal dimension): A 1-by-2 vector containing the dimensions of the two subsystems that X acts on.
  • PPT (optional, default 0): A flag (either 1 or 0) that indicates whether or not the desired symmetric extension must have positive partial transpose.
  • BOS (optional, default 0): A flag (either 1 or 0) that indicates whether or not the desired symmetric extension must be Bosonic (i.e., have its range contained within the symmetric subspace).
  • TOL (optional, default eps^(1/4)): The numerical tolerance used throughout this script. It is recommended that this is left at the default value unless numerical problems arise and the script has difficulty determining whether or not X has a symmetric extension.

Output arguments

  • EX: A flag (either 1 or 0) indicating that X does or does not have a symmetric extension of the desired type.
  • WIT (optional): A witness that verifies that the answer provided by EX is correct. If EX = 1 (i.e., X has a symmetric extension) then WIT is such a symmetric extension. If EX = 0 (i.e., no symmetric extension exists) then WIT is an entanglement witness with trace(WIT*X) = -1 but trace(WIT*Y) >= 0 for all symmetrically extendable Y.

Examples

2-qubit symmetric extension

It is known[1] that a 2-qubit state $\rho_{AB}$ has a (not necessarily PPT) symmetric extension if and only if ${\rm Tr}(\rho_B^2) \geq {\rm Tr}(\rho_{AB}^2) - 4\sqrt{\det(\rho_{AB})}$. The following code verifies that one such state does indeed have a symmetric extension.

>> rho = [1 0 0 -1;0 1 1/2 0;0 1/2 1 0;-1 0 0 1];
>> [trace(PartialTrace(rho)^2), trace(rho^2) - 4*sqrt(det(rho))] % if the first number is >= the second number, rho has a symmetric extension
 
ans =
 
    8.0000    6.5000
 
>> SymmetricExtension(rho) % verify that rho has a symmetric extension
 
ans =
 
     1

Notes

If your goal is to detect entanglement in an operator, then it is always in your best interest to set the optional arguments PPT and BOS to be 1. Setting BOS = 1 increases the effectiveness of the entanglement test without any computational overhead at all. Setting PPT = 1 slows down the computation quite a bit, but increases the effectiveness as an entanglement test considerably.

Source code

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

  1. %%  SYMMETRICEXTENSION    Determines whether or not an operator has a symmetric extension
  2. %   This function has one required argument:
  3. %     X: a positive semidefinite matrix
  4. %
  5. %   EX = SymmetricExtension(X) is either 1 or 0, indicating that X does or
  6. %   does not have a 2-copy positive partial transpose Bosonic symmetric
  7. %   extension. The extension is always taken on the second subsystem of X.
  8. %
  9. %   This function has five optional arguments:
  10. %     K (default 2)
  11. %     DIM (default has both subsystems of equal dimension)
  12. %     PPT (default 0)
  13. %     BOS (default 0)
  14. %     TOL (default eps^(1/4))
  15. %
  16. %   [EX,WIT] = SymmetricExtension(X,K,DIM,PPT,BOS,TOL) determines whether
  17. %   or not X has a symmetric extension and provides a witness WIT that
  18. %   verifies the answer. If a symmetric extension of X exists
  19. %   (i.e., EX = 1) then WIT is such a symmetric extension. If no symmetric
  20. %   extension exists (i.e., EX = 0) then WIT is an entanglement witness
  21. %   with trace(WIT*X) = -1 but trace(WIT*Y) >= 0 for all symmetrically
  22. %   extendable Y.
  23. %
  24. %   K is the desired number of copies of the second subsystem. DIM is a
  25. %   1-by-2 vector containing the dimensions of the subsystems on which X
  26. %   acts. PPT is a flag (either 0 or 1) indicating whether the desired
  27. %   symmetric extension must have positive partial transpose. BOS is a flag
  28. %   (either 0 or 1) indicating whether the desired symmetric extension must
  29. %   be Bosonic (i.e., be supported on the symmetric subspace). TOL is the
  30. %   numerical tolerance used when determining whether or not a symmetric
  31. %   extension exists.
  32. %
  33. %   URL: http://www.qetlab.com/SymmetricExtension
  34.  
  35. %   requires: cvx (http://cvxr.com/cvx/), IsPPT.m, IsPSD.m, opt_args.m,
  36. %             PartialTrace.m, PartialTranspose.m, PermutationOperator.m,
  37. %             PermuteSystems.m, sporth.m, SymmetricProjection.m
  38. %
  39. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  40. %   package: QETLAB
  41. %   last updated: December 12, 2014
  42.  
  43. function [ex,wit] = SymmetricExtension(X,varargin)
  44.  
  45. lX = length(X);
  46.  
  47. % set optional argument defaults: k=2, dim=sqrt(length(X)), ppt=1, bos=1, tol=eps^(1/4)
  48. [k,dim,ppt,bos,tol] = opt_args({ 2, round(sqrt(lX)), 0, 0, eps^(1/4) },varargin{:});
  49.  
  50. % allow the user to enter a single number for dim
  51. if(length(dim) == 1)
  52.     dim = [dim,lX/dim];
  53.     if abs(dim(2) - round(dim(2))) >= 2*lX*eps
  54.         error('SymmetricExtension:InvalidDim','If DIM is a scalar, it must evenly divide length(X); please provide the DIM array containing the dimensions of the subsystems.');
  55.     end
  56.     dim(2) = round(dim(2));
  57. end
  58.  
  59. % In certain situations, we don't need semidefinite programming.
  60. if(k == 1 || (lX <= 6 && ppt && nargout <= 1))
  61.     if(k == 1 && ~ppt) % in some cases, the problem is *really* trivial
  62.         if(nargout > 1)
  63.             [ex,wit] = IsPSD(X,tol);
  64.             if(ex == 1)
  65.                 wit = X;
  66.             end
  67.         else
  68.             ex = IsPSD(X,tol);
  69.         end
  70.  
  71.         return
  72.     end
  73.  
  74.     % In this case, all they asked for is a 1-copy PPT symmetric extension
  75.     % (i.e., they're asking if the state is PPT).
  76.     if(nargout > 1)
  77.         [ex,wit] = IsPPT(X,2,dim,tol);
  78.         if(ex)
  79.             wit = X;
  80.         else
  81.             wit = PartialTranspose(wit*wit');
  82.             wit = -wit/trace(wit*X);
  83.         end
  84.     else
  85.         ex = (IsPPT(X,2,dim,tol) && IsPSD(X,tol));
  86.     end
  87.  
  88. % In the 2-qubit case, an analytic formula is known for whether or not a
  89. % state has a (2-copy, non-PPT) symmetric extension that is much
  90. % faster to use than semidefinite programming.
  91. elseif(~isa(X,'cvx') && k == 2 && ~ppt && dim(1) == 2 && dim(2) == 2 && nargout <= 1) % we don't need "&& ~bos" thanks to a lemma of Myhr and Lutkenhaus
  92.     ex = (trace(PartialTrace(X,1)^2) >= trace(X^2) - 4*sqrt(det(X)));
  93.  
  94. % otherwise, use semidefinite programming to find a symmetric extension
  95. elseif(k > 1)
  96.     sdp_dim = [dim(1),dim(2)*ones(1,k)];
  97.     % For Bosonic extensions, it suffices to optimize over the symmetric
  98.     % subspace, which is smaller.
  99.     if(bos)
  100.         sdp_prod_dim = dim(1)*nchoosek(k+dim(2)-1, dim(2)-1);
  101.     else
  102.         sdp_prod_dim = dim(1)*dim(2)^k;
  103.     end
  104.  
  105.     cvx_begin sdp quiet
  106.         cvx_precision(tol);
  107.         variable rho(sdp_prod_dim,sdp_prod_dim) hermitian
  108.         if(nargout > 1) % don't want to compute the dual solution in general (especially not if this is called within CVX)
  109.             dual variable W
  110.         end
  111.         subject to
  112.             rho >= 0;
  113.             if(bos) % check for a Bosonic extension (faster *and* more effective)
  114.                 P = kron(speye(dim(1)),SymmetricProjection(dim(2),k,1));
  115.                 if(nargout > 1)
  116.                     W : PartialTrace(P*rho*P',3:k+1,sdp_dim) == X;
  117.                 else
  118.                     PartialTrace(P*rho*P',3:k+1,sdp_dim) == X;
  119.                 end
  120.                 if(ppt)
  121.                     PartialTranspose(P*rho*P',1:ceil(k/2)+1,sdp_dim) >= 0;
  122.                 end
  123.             else % check for a non-Bosonic extension (slower but perhaps sometimes useful)
  124.                 % It's a bit messy that the code from above is largely
  125.                 % repeated here, but the goal is to not slow down the
  126.                 % Bosonic optimization at all, so we don't want to
  127.                 % introduce a new variable sigma = P*rho*P' in that case.
  128.                 if(nargout > 1)
  129.                     W : PartialTrace(rho,3:k+1,sdp_dim) == X;
  130.                 else
  131.                     PartialTrace(rho,3:k+1,sdp_dim) == X;
  132.                 end
  133.                 for j = 3:k+1
  134.                     PartialTrace(rho,setdiff(2:k+1,j),sdp_dim) == X;
  135.                 end
  136.                 if(ppt)
  137.                     PartialTranspose(rho,1:ceil(k/2)+1,sdp_dim) >= 0;
  138.                 end
  139.             end
  140.     cvx_end
  141.  
  142.     ex = 1-min(cvx_optval,1); % CVX-safe way to map (0,Inf) to (1,0)
  143.     if(~isa(X,'cvx')) % make the output prettier if it's not a CVX input
  144.         ex = round(ex);
  145.  
  146.         % Deal with error messages and witnesses.
  147.         if(nargout > 1 && strcmpi(cvx_status,'Solved'))
  148.             if(bos)
  149.                 wit = P*rho*P';
  150.             else
  151.                 wit = rho;
  152.             end
  153.         elseif(strcmpi(cvx_status,'Inaccurate/Solved'))
  154.             if(bos)
  155.                 wit = P*rho*P';
  156.             else
  157.                 wit = rho;
  158.             end
  159.             warning('SymmetricExtension:NumericalProblems','Minor numerical problems encountered by CVX. Consider adjusting the tolerance level TOL and re-running the script.');
  160.         elseif(nargout > 1 && strcmpi(cvx_status,'Infeasible'))
  161.             wit = -W;
  162.         elseif(strcmpi(cvx_status,'Inaccurate/Infeasible'))
  163.             if(nargout > 1)
  164.                 wit = -W;
  165.             end
  166.             warning('SymmetricExtension:NumericalProblems','Minor numerical problems encountered by CVX. Consider adjusting the tolerance level TOL and re-running the script.');
  167.         elseif(strcmpi(cvx_status,'Unbounded') || strcmpi(cvx_status,'Inaccurate/Unbounded') || strcmpi(cvx_status,'Failed'))
  168.             error('SymmetricExtension:NumericalProblems',strcat('Numerical problems encountered (CVX status: ',cvx_status,'). Please try adjusting the tolerance level TOL.'));
  169.         end
  170.     end
  171. else
  172.     error('SymmetricExtension:InvalidK','K must be a positive integer.');
  173. end

References

  1. J. Chen, Z. Ji, D. Kribs, N. Lütkenhaus, and B. Zeng. Symmetric extension of two-qubit states. E-print: arXiv:1310.3530 [quant-ph]