SkOperatorNorm

 Other toolboxes required SkOperatorNorm Computes the S(k)-norm of an operator cvx IsBlockPositiveSkVectorNorm Norms 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. N. Johnston. Norms and Cones in the Theory of Quantum Entanglement. PhD thesis, University of Guelph, 2012. E-print: arXiv:1207.1479 [quant-ph]