Sk iterate

From QETLAB
Jump to: navigation, search
sk_iterate
Computes a lower bound of the S(k)-norm of an operator

Other toolboxes required none
Related functions SchmidtDecomposition
SchmidtRank
SkOperatorNorm
Function category Helper functions
This is a helper function that only exists to aid other functions in QETLAB. If you are an end-user of QETLAB, you likely will never have a reason to use this function.

sk_iterate is a function that iteratively computes a lower bound on 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. The method used to compute this lower bound is described here.

Syntax

  • SK = sk_iterate(X)
  • SK = sk_iterate(X,K)
  • SK = sk_iterate(X,K,DIM)
  • SK = sk_iterate(X,K,DIM,TOL)
  • SK = sk_iterate(X,K,DIM,TOL,V0)
  • [SK,V] = sk_iterate(X,K,DIM,TOL,V0)

Argument descriptions

Input arguments

  • X: A square positive semidefinite matrix to have its S(k)-norm bounded.
  • K (optional, default 1): A positive integer, the Schmidt rank to optimize over.
  • 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.
  • TOL (optional, default \(10^{-5}\)): The numerical tolerance used when determining whether or not the iterative procedure has converged.
  • V0 (optional, default is randomly-generated): The vector to begin the iterative procedure from.

Output arguments

  • V (optional): A vector with Schmidt rank at most K such that V'*X*V == SK.

Examples

A two-qubit example

In [3], it was shown that the density matrix $$ \rho = \frac{1}{8}\begin{bmatrix}5 & 1 & 1 & 1\\1 & 1 & 1 & 1\\1 & 1 & 1 & 1\\1 & 1 & 1 & 1\end{bmatrix} $$ has S(1)-norm equal to $(3+2\sqrt{2})/8 \approx 0.7286$. The following code shows that this quantity is indeed a lower bound of the S(1)-norm:

>> rho = [5 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1]/8;
>> sk_iterate(rho)
 
ans =
 
    0.7286

Source code

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

  1. %%  SK_ITERATE    Computes a lower bound of the S(k)-norm of an operator
  2. %   This function has one required argument:
  3. %     X: a square positive semidefinite matrix
  4. %
  5. %   SK = sk_iterate(X) is a lower bound of the S(1)-norm of X (i.e., the
  6. %   maximum overlap of X with a separable pure state -- see references
  7. %   [1,2,3]). The two subsystems on which X acts are assumed to have equal
  8. %   dimension in this case (specify the optional DIM parameter if they are
  9. %   of unequal dimension). The lower bound is computed via a randomized
  10. %   method that searches for local maxima. X can be full or sparse.
  11. %
  12. %   This function has four optional arguments:
  13. %     K (the Schmidt rank to optimize, default 1)
  14. %     DIM (default has both subsystems of equal dimension)
  15. %     TOL (default 10^-5)
  16. %     V0 (default is a randomly-chosen vector)
  17. %   
  18. %   [SK,V] = sk_iterate(X,K,DIM,TOL,V0) computes a lower bound of the
  19. %   S(K)-norm of X, as above. The search for a local maximum starts with
  20. %   the vector V0, and the dimensions of the subsystems that X acts on are
  21. %   provided by the 1-by-2 vector DIM. The algorithm terminates when two
  22. %   iterations result in values that are within TOL of each other. V is the
  23. %   local maximizing vector of Schmidt rank <= K.
  24. %
  25. %   URL: http://www.qetlab.com/sk_iterate
  26. %
  27. %   References:
  28. %   [1] N. Johnston and D. W. Kribs. A Family of Norms With Applications in
  29. %       Quantum Information Theory. Journal of Mathematical Physics,
  30. %       51:082202, 2010.
  31. %   [2] N. Johnston and D. W. Kribs. A Family of Norms With Applications in
  32. %       Quantum Information Theory II. Quantum Information & Computation,
  33. %       11(1 & 2):104-123, 2011.
  34. %   [3] N. Johnston. Norms and Cones in the Theory of Quantum Entanglement.
  35. %       PhD thesis, University of Guelph, 2012.
  36.  
  37. %   requires: iden.m, MaxEntangled.m, normalize_cols.m, opt_args.m,
  38. %             PermuteSystems.m, SchmidtDecomposition.m, SchmidtRank.m,
  39. %             sporth.m, Swap.m
  40. %             
  41. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  42. %   package: QETLAB
  43. %   last updated: November 14, 2014
  44.  
  45. %% WHEN UPDATING FOR NON-HERMITIAN, ALSO UPDATE HOW IT'S USED IN SKOPERATORNORM
  46.  
  47. function [Sk,v] = sk_iterate(X,varargin)
  48.  
  49. dX = length(X);
  50. sdX = round(sqrt(dX));
  51.  
  52. % Set optional argument defaults: k=1, dim=sqrt(length(X)), tol=10^-5, v0
  53. % is a random initial vector (set to -1 for now).
  54. [k,dim,tol,v0] = opt_args({ 1, sdX, 10^(-5), -1 },varargin{:});
  55.  
  56. % allow the user to enter a single number for dim
  57. if(length(dim) == 1)
  58.     dim = [dim,dX/dim];
  59.     if abs(dim(2) - round(dim(2))) >= 2*dX*eps
  60.         error('sk_iterate: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.');
  61.     end
  62.     dim(2) = round(dim(2));
  63. end
  64. da = dim(1);
  65. db = dim(2);
  66.  
  67. % Some of this preparation is unnecessary when k = 1, but it's cheap so we
  68. % don't really care.
  69. psi = MaxEntangled(k,1,0);
  70. psiI = Swap(kron(psi,speye(da*db)),[2,3],[k,k,da,db],1);
  71. PSI = psiI*psiI';
  72. PSIX = Swap(kron(psi*psi',X),[2,3],[k,k,da,db],0);
  73.  
  74. % If the user specified a starting guess v0, parse it; otherwise randomly
  75. % generate one.
  76. randv0 = 1;
  77. if(max(size(v0)) > 1)
  78.     [s0,a0,b0] = SchmidtDecomposition(v0,dim);
  79.     sr = length(s0);
  80.     if(sr > k)
  81.         warning('sk_iterate:SchmidtRankMismatch','The Schmidt rank of the initial vector v0 is %d, which is larger than k=%d. Using a randomly-generated intial vector instead.',sr,k);
  82.     else
  83.         randv0 = 0;
  84.         vp(:,1) = padarray(reshape(a0*diag(s0),da*sr,1),da*(k-sr),'post');
  85.         vp(:,2) = padarray(reshape(b0,db*sr,1),db*(k-sr),'post');
  86.     end
  87. end
  88. if randv0 % generate a random starting vector v0, if appropriate
  89.     vp = randn(max(dim)*k,2) + 1i*randn(max(dim)*k,2);
  90.     v = psiI'*kron(vp(1:k*da,1),vp(1:k*db,2));
  91.     v = v/norm(v);
  92. else
  93.     v = v0;
  94. end
  95. vp = normalize_cols(vp);
  96.  
  97. % Preparation is done; now do the actual iteration.
  98. it_err = 2*tol+1;
  99. Sk = v'*X*v;
  100. while it_err > tol
  101.     it_err = 0;
  102.  
  103.     % Loop through the 2 parties.
  104.     for p = 0:1
  105.         % If Schmidt rank is not full, we will have numerical problems; go to
  106.         % lower Schmidt rank iteration.
  107.         sr = SchmidtRank(v,dim,1000*eps);
  108.         if(sr < k)
  109.             [Sk,v] = sk_iterate(X,sr,dim,tol,v);
  110.             break;
  111.         end
  112.  
  113.         % Fix one of the parties and optimize over the other party.
  114.         if(p==1)
  115.             V0 = kron(vp(1:da*k,1),speye(db*k));
  116.         else
  117.             V0 = kron(speye(da*k),vp(1:db*k,2));            
  118.         end
  119.         V1 = V0'*PSI*V0;
  120.  
  121.         try
  122.             [vp(1:size(V0,2),p+1),NSk] = eigs(V0'*PSIX*V0,V1,1,'LR');
  123.         catch err
  124.             % In case of ARPACK errors, try to compute the maximal
  125.             % eigenvalue by converting to full and using pinv, if it is
  126.             % reasonable to do so. Otherwise, just ignore it altogether and
  127.             % return the best lower bound found so far.
  128.             if(strcmpi(err.identifier,'MATLAB:eigs:ARPACKroutineError'))
  129.                 %if(sdX <= 50)
  130.                 %    [tmp,NSk] = eig(pinv(full(V1))*full(V0'*PSIX*V0));
  131.                 %    [NSk,ind] = max(diag(NSk));
  132.                 %    vp(:,p+1) = tmp(:,ind);
  133.                 %else
  134.                     return
  135.                 %end
  136.             else
  137.                 rethrow(err);
  138.             end
  139.         end
  140.         vp(:,p+1) = vp(:,p+1)/norm(vp(:,p+1));
  141.  
  142.         it_err = it_err + abs(Sk-NSk);
  143.         Sk = real(NSk);
  144.         v = psiI'*kron(vp(1:k*da,1),vp(1:k*db,2));
  145.         v = v/norm(v);
  146.     end
  147. 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. N. Johnston. Norms and Cones in the Theory of Quantum Entanglement. PhD thesis, University of Guelph, 2012. E-print: arXiv:1207.1479 [quant-ph]