SkOperatorNorm

From QETLAB
Jump to: navigation, search
SkOperatorNorm
Computes the S(k)-norm of an operator

Other toolboxes required cvx
Related functions IsBlockPositive
SkVectorNorm
Function category Norms
Usable within CVX? no

SkOperatorNorm is a function that computes the S(k)-norm of an operator[1][2]: $$ \|X\|_{S(k)} := \sup_{|v\rangle , |w\rangle } \Big\{ \big| \langle w| X |v \rangle \big| : SR(|v \rangle), SR(|v \rangle) \leq k, \big\||v \rangle\big\| = \big\||w \rangle\big\| = 1 \Big\}, $$ where $SR(\cdot)$ refers to the Schmidt rank of a pure state.

Syntax

  • LB = SkOperatorNorm(X)
  • LB = SkOperatorNorm(X,K)
  • LB = SkOperatorNorm(X,K,DIM)
  • LB = SkOperatorNorm(X,K,DIM,STR)
  • LB = SkOperatorNorm(X,K,DIM,STR,TARGET)
  • LB = SkOperatorNorm(X,K,DIM,STR,TARGET,TOL)
  • [LB,~,UB] = SkOperatorNorm(X,K,DIM,STR,TARGET,TOL)
  • [LB,LWIT,UB,UWIT] = SkOperatorNorm(X,K,DIM,STR,TARGET,TOL)

Argument descriptions

Input arguments

  • X: A square matrix acting on bipartite space. Generally, X should be positive semidefinite – the bounds produced otherwise are quite poor.
  • K (optional, default 1): A positive integer.
  • DIM (optional, by default has both subsystems of equal dimension): A 1-by-2 vector containing the dimensions of the subsystems that X acts on.
  • STR (optional, default 2): An integer that determines how hard the script should work to compute the lower and upper bounds (STR = -1 means that the script won't stop working until the bounds match each other). Other valid values are 0, 1, 2, 3, .... In practice, if STR >= 4 then most computers will run out of memory and/or the sun will explode before computation completes.
  • TARGET (optional, default -1): A value that you wish to prove that the norm is above or below. If, at any point while this script is running, it proves that LB >= TARGET or that UB <= TARGET, then the script will immediately abort and return the best lower bound and upper bound computed up to that point. This is a time-saving feature that can be avoided by setting TARGET = -1.
  • TOL (optional, default eps^(3/8)): The numerical tolerance used throughout the script.

Output arguments

  • LB: A lower bound of the S(k)-operator norm of X.
  • LWIT: A witness that verifies that LB is indeed a lower bound of the S(k)-operator norm of X. More specifically, LWIT is a unit vector such that SchmidtRank(LWIT,DIM) <= K and LWIT'*X*LWIT = LB.
  • UB: An upper bound of the S(k)-operator norm of X.
  • UWIT: A witness that verifies that UB is indeed an upper bound of the S(k)-operator norm of X. More specifically, UWIT is a feasible point of the dual semidefinite program presented in Section 5.2.3 of [3].

Examples

Exact computation in small dimensions

When X lives in $M_2 \otimes M_2$, $M_2 \otimes M_3$, or $M_3 \otimes M_2$ (i.e., when prod(DIM) <= 6), the script is guaranteed to compute the exact value of $\|X\|_{S(1)}$:

>> X = [5 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1]/8;
>> SkOperatorNorm(X)
 
ans =
 
    0.7286

The fact that this computation is correct is illustrated in Example 5.2.11 of [3], where it was shown that the S(1)-norm is exactly $(3 + 2\sqrt{2})/8 \approx 0.7286$. However, if we were still unconvinced, we could request witnesses that verify that 0.7286 is both a lower bound and an upper bound of the S(1)-norm:

>> [lb,lwit,ub,uwit] = SkOperatorNorm(X)
 
lb =
 
    0.7286
 
 
lwit =
 
   0.8536 + 0.0000i
   0.3536 - 0.0000i
   0.3536 + 0.0000i
   0.1464          
 
 
ub =
 
    0.7286
 
 
uwit =
 
   0.0518 + 0.0000i  -0.0625 + 0.0000i  -0.0625 - 0.0000i  -0.1250 - 0.0000i
  -0.0625 - 0.0000i   0.3018 + 0.0000i   0.0000 + 0.0000i  -0.0625 - 0.0000i
  -0.0625 + 0.0000i   0.0000 - 0.0000i   0.3018 + 0.0000i  -0.0625 + 0.0000i
  -0.1250 + 0.0000i  -0.0625 + 0.0000i  -0.0625 - 0.0000i   0.3018 + 0.0000i
 
>> lwit'*X*lwit % verify that the lower bound is correct
 
ans =
 
   0.7286 + 0.0000i
 
>> norm(X + uwit) % verify that the upper bound is correct
 
ans =
 
    0.7286

Only interested in the lower and upper bounds; not the witnesses

If all you want are the lower and upper bounds, but don't require the witnesses LWIT and UWIT, you can use code like the following. Note that in this case, $\|X\|_{S(1)}$ is computed exactly, as the lower and upper bound are equal (though this will not happen for all X!). However, all we know about $\|X\|_{S(2)}$ is that it lies in the interval [0.3522, 0.3546]. It is unsurprising that $\|X\|_{S(3)}$ is the usual operator norm of X, since this is always the case when K >= min(DIM).

>> X = RandomDensityMatrix(9);
>> [lb,~,ub] = SkOperatorNorm(X,1)
 
lb =
 
    0.2955
 
 
ub =
 
    0.2955
 
>> [lb,~,ub] = SkOperatorNorm(X,2)
 
lb =
 
    0.3522
 
 
ub =
 
    0.3546
 
>> [lb,~,ub] = SkOperatorNorm(X,3)
 
lb =
 
    0.3770
 
 
ub =
 
    0.3770
 
>> norm(X)
 
ans =
 
    0.3770

Source code

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

  1. %%  SKOPERATORNORM    Bounds the S(k)-norm of an operator
  2. %   This function has one required input argument:
  3. %     X: a square matrix
  4. %
  5. %   LB = SkOperatorNorm(X) is a lower bound of the S(1)-norm of the
  6. %   operator X, which is assumed to act on two systems of the same
  7. %   dimension. Note that X should be positive semidefinite in almost all
  8. %   cases in which this function is used -- the bounds provided are
  9. %   generally not very good otherwise.
  10. %
  11. %   This function has five optional input arguments:
  12. %     K (default 1)
  13. %     DIM (default has both subsystems of equal dimension)
  14. %     STR (default 2)
  15. %     TARGET (default -1)
  16. %     TOL (default eps^(3/8))
  17. %
  18. %   [LB,LWIT,UB,UWIT] = SkOperatorNorm(X,K,DIM,STR,TARGET,TOL) provides a
  19. %   lower bound (LB) and an upper bound (UB) of the S(K)-norm of the
  20. %   operator X, which acts on two subsystems, whose dimensions are given in
  21. %   the 1-by-2 vector DIM.
  22. %
  23. %   K is the "index" of the norm -- that is, it is the Schmidt rank of the
  24. %   vectors that are multiplying X on the left and right in the definition
  25. %   of the norm.
  26. %
  27. %   STR is the amount of computation that you want to devote to computing
  28. %   the bounds. Higher values of STR correspond to more computation time
  29. %   (and thus better bounds).
  30. %
  31. %   TARGET is a target value that you wish to prove that the norm is above
  32. %   or below. Once the script has proved that LB >= TARGET or UB <= TARGET,
  33. %   the script immediately aborts. This is a time-saving feature that can
  34. %   be avoided by setting TARGET = -1.
  35. %
  36. %   TOL is the numerical tolerance used throughout the script.
  37. %
  38. %   LWIT and UWIT are witnesses that verify that the bounds LB and UB are
  39. %   correct. More specifically, LWIT is a vector with Schmidt rank <= K
  40. %   such that LWIT'*X*LWIT = LB, and UWIT is the optimal positive
  41. %   semidefinite operator Y in the dual semidefinite program (5.2.3) in [3]
  42. %   (see online documentation for more details).
  43. %
  44. %   Most of the lower bounds and basic facts about these norms were derived
  45. %   in [1]. The semidefinite program method for computing upper bounds was
  46. %   derived in [2]. See [3] for a summary of results and more information.
  47. %
  48. %   URL: http://www.qetlab.com/SkOperatorNorm
  49. %
  50. %   References:
  51. %   [1] N. Johnston and D. W. Kribs. A Family of Norms With Applications in
  52. %       Quantum Information Theory. Journal of Mathematical Physics,
  53. %       51:082202, 2010.
  54. %   [2] N. Johnston and D. W. Kribs. A Family of Norms With Applications in
  55. %       Quantum Information Theory II. Quantum Information & Computation,
  56. %       11(1 & 2):104-123, 2011.
  57. %   [3] N. Johnston. Norms and Cones in the Theory of Quantum Entanglement.
  58. %       PhD thesis, University of Guelph, 2012.
  59.  
  60. %   requires: cvx (http://cvxr.com/cvx/), iden.m, IsPSD.m, kpNorm.m,
  61. %             MaxEntangled.m, normalize_cols.m, opt_args.m, PartialMap.m,
  62. %             PartialTrace.m, PartialTranspose.m, PermuteSystems.m,
  63. %             Realignment.m, SchmidtDecomposition.m, SchmidtRank.m,
  64. %             sk_iterate.m, SkVectorNorm.m, sporth.m, Swap.m
  65. %
  66. %   author: Nathaniel Johnston (nathaniel@njohnston.ca), based on joint
  67. %           work with David W. Kribs
  68. %   package: QETLAB
  69. %   last updated: September 22, 2014
  70.  
  71. function [lb,lwit,ub,uwit] = SkOperatorNorm(X,varargin)
  72.  
  73.     X = full(X);
  74.     dX = length(X);
  75.     sdX = round(sqrt(dX));
  76.     lwit = 0;
  77.     uwit = 0;
  78.  
  79.     % Set optional argument defaults: k=1, dim=sqrt(length(X)), str=2,
  80.     % target=-1, tol=eps^(3/8)
  81.     [k,dim,str,target,tol] = opt_args({ 1, [sdX sdX], 2, -1, eps^(3/8) },varargin{:});
  82.     if(str == -1)
  83.         str = 1/eps; % keep going forever!
  84.     end
  85.  
  86.     % allow the user to enter a single number for dim
  87.     if(length(dim) == 1)
  88.         dim = [dim,dX/dim];
  89.         if abs(dim(2) - round(dim(2))) >= 2*dX*eps
  90.             error('SkOperatorNorm:InvalidDim','If DIM is a scalar, X must be square and DIM must evenly divide length(X); please provide the DIM array containing the dimensions of the subsystems.');
  91.         end
  92.         dim(2) = round(dim(2));
  93.     end
  94.  
  95.     % some useful, repeatedly-used, values
  96.     prod_dim = prod(dim);
  97.     op_norm = norm(X);
  98.     x_rank = rank(X);
  99.     X = X/op_norm; % rescale X to have unit norm
  100.  
  101.     % The S(k)-norm is just the operator norm if k is large enough.
  102.     if k >= min(dim)
  103.         lb = op_norm;
  104.         ub = op_norm;
  105.         return
  106.     end
  107.  
  108.     % If X is rank 1 then the S(k)-norm is easy to compute via Proposition
  109.     % 10 of [1].
  110.     if x_rank == 1
  111.         [U,~,V] = svds(X,1);
  112.         lb = op_norm*SkVectorNorm(U,k,dim)*SkVectorNorm(V,k,dim);
  113.         ub = lb;
  114.         return
  115.     end
  116.  
  117.     isprojection = 0;
  118.     ishermitian = 0;
  119.     ispositive = 0;
  120.     isnormal = (max(max(abs(X'*X-X*X'))) <= 2*prod_dim*eps);
  121.     if isnormal
  122.         ishermitian = (max(max(abs(X-X'))) <= 100*eps);
  123.         if ishermitian
  124.             X = (X+X')/2; % avoids some numerical problems later
  125.             ispositive = IsPSD(X);
  126.             if ispositive
  127.                 isprojection = (max(max(abs(X-X*X))) <= eps*prod_dim^2);
  128.             end
  129.         end
  130.     end
  131.     is_trans_exact = (min(dim) == 2 && max(dim) <= 3);
  132.  
  133.     % Compute some more simple bounds. We will multiply these by op_norm before
  134.     % we leave this function.
  135.     lb = (k/min(dim));    % comes from Theorem 17 in [1]
  136.     ub = 1; % our most basic upper bound
  137.  
  138.     % break out of the function if the target value has already been met
  139.     if(MetTarget(lb,ub,op_norm,tol,target))
  140.         lb = op_norm*lb;
  141.         ub = op_norm*ub;
  142.         return
  143.     end
  144.  
  145.     if ~(ispositive && is_trans_exact && k == 1 && str >= 1) % if the exact answer won't be found by SDP, compute bounds via other methods first
  146.         if(ishermitian)
  147.             [eig_vec,eig_val] = eig(X);
  148.             [eig_val,ord] = sort(real(diag(eig_val)),'ascend');
  149.             eig_vec = eig_vec(:,ord);
  150.  
  151.             % use the lower bound of Proposition 4.14 of [1]
  152.             for r = k:min(dim)
  153.                 t_ind = prod(dim) - prod(dim - r);
  154.                 lb = max([(k/r)*eig_val(t_ind),lb]);
  155.             end
  156.  
  157.             % use the lower bound of Theorem 4.2.15 of [3]
  158.             if(k==1)
  159.                 lb = max([(trace(X) + sqrt((prod_dim*trace(X^2)-trace(X)^2)/(prod_dim-1)))/prod_dim,lb]);
  160.             end
  161.         end
  162.  
  163.         if(ispositive && nargout > 2)
  164.             % Use the upper bound of Proposition 15 of [1].
  165.             t_ub = 0;
  166.             for i = 1:prod_dim
  167.                 t_ub = t_ub + abs(eig_val(i))*SkVectorNorm(eig_vec(:,i),k,dim)^2;
  168.             end
  169.             ub = min(t_ub,ub);
  170.  
  171.             % Use the upper bound of Proposition 4.2.11 of [3].
  172.             ub = min(kpNorm(Realignment(X,dim),k^2,2),ub);
  173.         end
  174.  
  175.         % Use the lower bound of Theorem 4.2.17 of [3].
  176.         if(isprojection)
  177.             lb = max(min([1,k/ceil((dim(1)+dim(2) - sqrt((dim(1)-dim(2))^2 + 4*x_rank - 4))/2)]),lb);
  178.             lb = max((min(dim)-k)*(x_rank+sqrt((prod_dim*x_rank-x_rank^2)/(prod_dim-1)))/(prod_dim*(min(dim)-1)) + (k-1)/(min(dim)-1),lb);
  179.         end
  180.  
  181.         % break out of the function if the target value has already been met
  182.         if(MetTarget(lb,ub,op_norm,tol,target))
  183.             lb = op_norm*lb;
  184.             ub = op_norm*ub;
  185.             return
  186.         end
  187.  
  188.         % Use a randomized iterative method to try to improve the lower
  189.         % bound.
  190.         if(ispositive)
  191.             for j=5^str:-1:1
  192.                 [tlb,twit] = sk_iterate(X,k,dim,max(tol^2,eps^(3/4)));
  193.                 if(tlb >= lb)
  194.                     lb = tlb;
  195.                     lwit = twit;
  196.  
  197.                     % break out of the function if the target value has already been met
  198.                     if(MetTarget(lb,ub,op_norm,tol,target))
  199.                         lb = op_norm*lb;
  200.                         ub = op_norm*ub;
  201.                         return
  202.                     end
  203.                 end
  204.             end
  205.         end
  206.     end
  207.  
  208.     % Start the semidefinite programming approach for getting upper bounds.
  209.     if(str >= 1 && ((lb + tol < ub && ispositive && nargout > 2) || (ispositive && is_trans_exact && k == 1)))
  210.         cvx_begin sdp quiet
  211.             variable rho(prod_dim,prod_dim) hermitian
  212.             dual variable uwit
  213.             maximize trace(X*rho)
  214.             subject to
  215.                 rho >= 0;
  216.                 trace(rho) <= 1;
  217.                 if(k == 1)
  218.                     uwit : PartialTranspose(rho,2,dim) >= 0;
  219.                 else
  220.                     uwit : k*kron(PartialTrace(rho,2,dim),eye(dim(2))) >= rho;
  221.                 end
  222.         cvx_end
  223.  
  224.         % format the output of cvx into a slightly more user-friendly form
  225.         if(strcmpi(cvx_status,'Solved'))
  226.             ub = min(ub,real(cvx_optval));
  227.             if(k == 1)
  228.                 uwit = op_norm*PartialTranspose(uwit,2,dim);
  229.             else
  230.                 uwit = op_norm*(k*kron(PartialTrace(uwit,2,dim),eye(dim(2))) - uwit);
  231.             end
  232.         else
  233.             error('SkOperatorNorm:NumericalProblems',strcat('Numerical problems encountered (cvx: ',cvx_status,').'));
  234.         end
  235.  
  236.         % In small dimensions, the transpose map gets the result exactly.
  237.         if(is_trans_exact && k == 1 && str >= 1)
  238.             lb = ub;
  239.             [lwit,twit] = eig(rho);
  240.             [~,twit] = max(diag(twit));
  241.             lwit = lwit(:,twit(1));
  242.         elseif(k == 1)
  243.              % we can also get decent lower bounds from the SDP results when k=1
  244.             gs = min(1-roots(jacobi_poly(dim(2)-2,1,1)));
  245.             xmineig = min(real(eig(X)));
  246.             tlb = real(cvx_optval)*(1 - dim(2)*gs/(2*dim(2)-1)) + xmineig*gs/(2*dim(2)-2);
  247.  
  248.             if(tlb > lb)
  249.                 lb = tlb;
  250.                 lwit = 0; % unfortunately, we don't have a lower bound witness anymore
  251.             end
  252.  
  253.             % done the str = 1 SDP, now get better upper bounds via symmetric
  254.             % extensions if str >= 2
  255.             for j = 2:str
  256.                  % break out of the function if the target value has already been met
  257.                 if(MetTarget(lb,ub,op_norm,tol,target))
  258.                     lb = op_norm*lb;
  259.                     ub = op_norm*ub;
  260.                     return
  261.                 end
  262.  
  263.                 sym_dim = [dim(1),dim(2)*ones(1,j)];
  264.                 prod_sym_dim = dim(1)*dim(2)^j;
  265.                 P = kron(eye(dim(1)),SymmetricProjection(dim(2),j));
  266.  
  267.                 cvx_begin sdp quiet
  268.                     variable rho(prod_sym_dim,prod_sym_dim) hermitian
  269.                     dual variable suwit
  270.                     maximize trace(X*PartialTrace(rho,3:j+1,sym_dim))
  271.                     subject to
  272.                         rho >= 0;
  273.                         trace(rho) <= 1;
  274.                         P*rho*P' == rho;
  275.                         suwit : PartialTranspose(rho,1:ceil(j/2)+1,sym_dim) >= 0;
  276.                 cvx_end
  277.  
  278.                 % format the output of cvx into a slightly more user-friendly form
  279.                 if(strcmpi(cvx_status,'Solved'))
  280.                     if(real(cvx_optval) < ub)
  281.                         ub = real(cvx_optval);
  282.                         uwit = op_norm*PartialTranspose(suwit,1:ceil(j/2)+1,sym_dim);
  283.                     end
  284.  
  285.                     gs = min(1-roots(jacobi_poly(dim(2)-2,mod(j,2),floor(j/2)+1)));
  286.                     tlb = real(cvx_optval)*(1 - dim(2)*gs/(2*dim(2)-1)) + xmineig*gs/(2*dim(2)-2);
  287.  
  288.                     if(tlb > lb)
  289.                         lb = tlb;
  290.                         lwit = 0; % unfortunately, we don't have a lower bound witness anymore
  291.                     end
  292.                 else
  293.                     error('SkOperatorNorm:NumericalProblems',strcat('Numerical problems encountered (cvx: ',cvx_status,').'));
  294.                 end
  295.             end
  296.         end
  297.     end
  298.  
  299.     lb = op_norm*lb;
  300.     ub = op_norm*ub;
  301. end
  302.  
  303.  
  304. % This function checks whether or not the lower bound or upper bound
  305. % already computed meets the desired target value (within numerical error)
  306. % and thus we can abort the script early.
  307. function mt = MetTarget(lb,ub,op_norm,tol,target)
  308.     mt = (op_norm*(lb + tol) >= op_norm*ub || (target >= 0 && (op_norm*(lb - tol) >= target || op_norm*(ub + tol) <= target)));    
  309. end

References

  1. N. Johnston and D. W. Kribs. A Family of Norms With Applications in Quantum Information Theory. J. Math. Phys., 51:082202, 2010. E-print: arXiv:0909.3907 [quant-ph]
  2. N. Johnston and D. W. Kribs. A Family of Norms With Applications in Quantum Information Theory II. Quantum Information & Computation, 11(1 & 2):104–123, 2011. E-print: arXiv:1006.0898 [quant-ph]
  3. 3.0 3.1 N. Johnston. Norms and Cones in the Theory of Quantum Entanglement. PhD thesis, University of Guelph, 2012. E-print: arXiv:1207.1479 [quant-ph]