# Sk iterate

 Other toolboxes required sk_iterate Computes a lower bound of the S(k)-norm of an operator none SchmidtDecompositionSchmidtRankSkOperatorNorm 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]