# NPAHierarchy

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
 Other toolboxes required NPAHierarchy Determines whether or not a set of probabilities satisfy the conditions of the NPA hierarchy CVX BellInequalityMaxXORGameValue Nonlocality and Bell inequalities yes (concave)

NPAHierarchy is a function that determines whether or not a given set of probabilities satisfy the constraints specified by the NPA hierarchy, which are a set of necessary conditions that probabilities arising from quantum mechanics must satisfy. That is, given a set of probabilities $\{p(a,b|x,y)\}$, this function tries to determine whether or not there is a quantum state $\rho$ shared between Alice and Bob, and sets of measurement operators $\{M_{a|x}\}$ for Alice and $\{M_{b|y}\}$ for Bob such that $p(a,b|x,y) = \mathrm{Tr}(\rho(M_{a|x} \otimes M_{b|y}))$ for all $a,b,x,y$ (here $x$ indexes the different measurement settings for Alice, while $a$ indexes the different measurement outcomes within that setting, and similarly for $y$ and $b$ for Bob).

## Syntax

• IS_NPA = NPAHierarchy(P)
• IS_NPA = NPAHierarchy(P,K)

## Argument descriptions

• P: A 4D array such that P(a,b,x,y) is the probability that Alice and Bob obtain measurement result (a,b) when they use measurement setting (x,y) (i.e., this is equal to $p(a,b|x,y)$ from above).
• K (optional, default 1): A non-negative integer that determines the level of the NPA hierarchy to use. If K = 0 then the NPA hierarchy isn't use at all, and this function just checks whether or not P is a valid probability array that satisfies non-signalling constraints (i.e., all entries are non-negative, each matrix face sums to 1, and the marginal probabilities for Alice and Bob are consistent with each other). Alternatively, K can be a string of a form like '1+ab+aab', which indicates that an intermediate level of the hierarchy should be used, where this example uses all products of 1 measurement, all products of one Alice and one Bob measurement, and all products of two Alice and one Bob measurement. Use plus signs to separate the different categories of products, as above. The first character of this string should always be a number, indicating the base level to use.

## Examples

### Can be used within CVX

Within CVX, this function is treated as any other concave function (so you can use it as an equality or >= constraint, or maximize it subject to other constraints). Maximizing it has the interpretation of trying to find a set of probabilities that satisfy the constraints of the NPA hierarchy, as well as whatever other constraints you specify. For example, the following code finds a $2$-measurement, $2$-outcome probability array satisfying the constraints of the second level of the hierarchy and also has $p(1,1|1,1) = p(1,2|1,1) + 0.1$ and $p(1,2|1,2) = p(2,2|1,2) - 0.2$ (these last two constraints have no special significance; they are just there for illustrative purposes):

>> cvx_begin sdp quiet
variable p(2,2,2,2); % 2 measurements, 2 outcomes each
maximize NPAHierarchy(p,2) % second level of the hierarchy
subject to
% these constraints do not have any particular significance
p(1,1,1,1) == p(1,2,1,1) + 0.1;
p(1,2,1,2) == p(2,2,1,2) - 0.2;
cvx_end
cvx_optval

cvx_optval =

1 % a value of 1 indicates that a probability array satisfying the constraints exists

>> p % let's have a look at it

p(:,:,1,1) =

0.2637    0.1637
0.2698    0.3029

p(:,:,2,1) =

0.2667    0.2333
0.2667    0.2333

p(:,:,1,2) =

0.2692    0.1582
0.2145    0.3582

p(:,:,2,2) =

0.2418    0.2582
0.2418    0.2582

## Source code

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

1. %%  NPAHIERARCHY    Determines whether or not a set of probabilities satisfy the conditions of the NPA hierarchy
2. %   This function has one required input argument:
3. %     P: a 4D array of probabilities, where P(a,b,x,y) is the probability
4. %        that Alice and Bob get measurement outcome (a,b) when they use
5. %        measurement setting (x,y)
6. %
7. %   IS_NPA = NPAHierarchy(P) is either 1 or 0, indicating that the
8. %   probabilities in the array P do or do not satisfy the conditions of the
9. %   first level of the NPA heierarchy.
10. %
11. %   This function has one optional input argument:
12. %     K (default 1): a non-negative integer, or a string indicating which
13. %       level of the hierarchy to check (see below for details)
14. %
15. %   IS_NPA = NPAHierarchy(P,K) is as above, but the K-th level of the NPA
16. %   hierarchy is checked. If K = 0 then it is just verified that P is a
17. %   valid probability array (i.e., each entry is non-negative, each matrix
18. %   "face" adds up to 1, and Alice's and Bob's marginal probabilities are
19. %   consistent with each other).
20. %
21. %   If K is a string, it must be a string of a form like '1+ab+aab', which
22. %   indicates that the intermediate level of the hierarchy should be used,
23. %   which uses all products of 1 measurement, all products of one Alice and
24. %   one Bob measurement, and all products of two Alice and one Bob
25. %   measurement. Use plus signs to separate the different categories of
26. %   products, as above. The first character of this string should always be
27. %   a number, indicating the base level to use.
28. %
29. %   URL: http://www.qetlab.com/NPAHierarchy
30. 
31. %   requires: CVX (http://cvxr.com/cvx/), opt_args.m
32. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
33. %   package: QETLAB
34. %   last updated: January 22, 2015
35. 
36. function is_npa = NPAHierarchy(p,varargin)
37. 
38.     % set optional argument defaults: K=1
39.     [k] = opt_args({ 1 },varargin{:});
40. 
41.     % Parse the input argument K to determine which measurement operators
42.     % to use. BASE_K is the integer part of the input (i.e., we use all
43.     % operators that are the product of BASE_K or fewer measurements) and K
44.     % is the maximum number of products that we ever use (e.g., '1+ab+aab'
45.     % would result in BASE_K = 1, K = 3).
46.     if(isnumeric(k))
47.         base_k = k;
48.         num_k_compon = 1;
49.     elseif(isa(k,'char'))
50.         k_types = textscan(lower(k),'%s','Delimiter','+'); % works with old versions of MATLAB, unlike strsplit
51.         k_types = k_types{1};
52.         base_k = str2double(k_types{1});
53. 
54.         num_k_compon = length(k_types);
55.         if(num_k_compon > 1)
56.             k = base_k;
57.             for j = 2:num_k_compon
58.                 k_types{j} = strtrim(k_types{j});
59.                 k = max(k,length(k_types{j}));
60.             end
61.         else
62.             k = base_k;
63.         end
64.     else
65.         error('NPAHierarchy:InvalidK','K must be a positive integer or a string.');
66.     end
67. 
68.     % Start by computing the number measurement settings for Alice and Bob (MA
69.     % and MB) and the number of outcomes for each measurement setting (OA and
70.     % OB).
71.     [oa,ob,ma,mb] = size(p);
72.     o_vec = [oa;ob];
73.     m_vec = [ma;mb];
74.     tot_dim = 1 + ((o_vec-1)'*m_vec)^k; % upper bound on the dimension of the (compact version of the) matrix GAMMA used in the NPA SDP
75.     tol = eps^(3/4); % numerical tolerance used
76. 
77.     % Make sure that P really is a probability array and that its marginal
78.     % probabilities are consistent with each other (but only if P is not a
79.     % CVX variable... if it is a CVX variable we will enforce these
80.     % constraints within the SDP below).
81.     if(~isa(p,'cvx'))
82.         % Require that P is a probability matrix (all entries are
83.         % non-negative and its faces sum to 1).
84.         if(min(min(min(min(p)))) < -tol^(3/4))
85.             is_npa = 0;
86.             return;
87.         end
88. 
89.         for i = 1:ma
90.             for j = 1:mb
91.                 if(abs(sum(sum(p(:,:,i,j))) - 1) > tol^(3/4))
92.                     is_npa = 0;
93.                     return;
94.                 end
95.             end
96.         end
97. 
98.         % Let's check Bob's marginal probabilities.
99.         for i = 1:ob
100.             for j = 1:mb
101.                 marg = sum(squeeze(p(:,i,:,j)),1);
102.                 if(max(abs(marg - marg(1))) > tol^(3/4))
103.                     is_npa = 0;
104.                     return;
105.                 end
106.             end
107.         end
108. 
109.         % Now check Alice's marginal probabilities.
110.         for i = 1:oa
111.             for j = 1:ma
112.                 marg = sum(squeeze(p(i,:,j,:)),1);
113.                 if(max(abs(marg - marg(1))) > tol^(3/4))
114.                     is_npa = 0;
115.                     return;
116.                 end
117.             end
118.         end
119.     end
120. 
121.     % Check the NPA SDP (if K is large enough).
122.     if(k >= 1 || isa(p,'cvx'))
123.         i_ind = [zeros(1,k);-ones(1,k)];
124.         j_ind = [zeros(1,k);-ones(1,k)];
125. 
126.         if(k >= 1)
127.             % Start by generating all of the product of measurements that
128.             % you need.
129.             ind_catalog = cell(0);
130.             for j = 1:tot_dim
131.                 [res,res_type] = product_of_orthogonal(j_ind,m_vec);
132.                 res_fnd = find_in_cell(res,ind_catalog);
133. 
134.                 % Make sure that this measurement is (1) new, and (2) valid
135.                 % given the user input.
136.                 if(res_fnd == 0 && res_type ~= 0)
137.                     is_valid_res = (size(res,2) <= base_k);
138.                     if(~is_valid_res && num_k_compon >= 2)
139.                         num_a_res = sum(res(1,:) < m_vec(1));
140.                         num_b_res = size(res,2) - num_a_res;
141.                         for i = 2:num_k_compon
142.                             num_a_fnd = length(find(k_types{i}=='a'));
143.                             num_b_fnd = length(find(k_types{i}=='b'));
144.                             if(num_a_res <= num_a_fnd && num_b_res <= num_b_fnd)
145.                                 is_valid_res = true;
146.                                 break;
147.                             end
148.                         end
149.                     end
150. 
151.                     if(is_valid_res)
152.                         ind_catalog{end+1} = res;
153.                     end
154.                 end
155. 
156.                 j_ind = update_ind(j_ind,k,m_vec,o_vec-1);
157.             end
158.             real_dim = length(ind_catalog);
159.         end
160. 
161.         cvx_begin quiet
162.             cvx_precision(tol);
163.             % We only enforce the actual NPA constraints if K >= 1... if
164.             % K = 0 we are just verifying marginals are consistent (i.e.,
165.             % no-signalling).
166.             if(k >= 1)
167.                 variable G(real_dim,real_dim) symmetric
168.             end
169.             subject to
170.                 if(k >= 1)
171.                     res_catalog = cell(0);
172.                     res_loc = cell(0);
173. 
174.                     % The following double loop loops over all entries of G and
175.                     % enforces entry-by-entry the (somewhat complicated) set of NPA
176.                     % constraints.
177.                     for i = 1:real_dim
178.                         for j = i:real_dim
179.                             % First determine what "type" of product of
180.                             % measurement operators the given matrix entry
181.                             % corresponds to (see product_of_orthogonal function
182.                             % below for details).
183.                             [res,res_type] = product_of_orthogonal([fliplr(ind_catalog{i}),ind_catalog{j}],m_vec);
184. 
185.                             % Entry is 0 if S_i^dagger*S_j = 0.
186.                             if(res_type == 0)
187.                                 G(i,j) == 0;
188. 
189.                             % Entry is a single probability from the P array if
190.                             % S_i^dagger*S_j measures on both Alice and Bob's
191.                             % sytems.
192.                             elseif(res_type == 2)
193.                                 G(i,j) == p(res(2,1)+1,res(2,2)+1,res(1,1)+1,res(1,2)-m_vec(1)+1);
194. 
195.                             % Entry is a sum of probabilities from the P array if
196.                             % S_i^dagger*S_j measures on just one system.
197.                             elseif(res_type == 1)
198.                                 if(isequal(res,[0;-1]))
199.                                     G(i,j) == 1; % identity measurement
200.                                 elseif(res(1) >= m_vec(1)) % measure on Bob's system
201.                                     G(i,j) == sum(p(:,res(2)+1,1,res(1)-m_vec(1)+1));
202.                                 else % measure on Alice's system
203.                                     G(i,j) == sum(p(res(2)+1,:,res(1)+1,1));
204.                                 end
205. 
206.                             % Entry is a product of non-commuting
207.                             % measurement operators. We can't specify its
208.                             % value, but we can specify that it is equal to
209.                             % other entries that are the *same* product of
210.                             % measurement operators.
211.                             else % res_type == -1
212.                                 % Check to see if we have run into this
213.                                 % particular RES before.
214.                                 res_fnd = find_in_cell(res,res_catalog);
215.                                 if(res_fnd == 0)
216.                                     res_fnd = find_in_cell(product_of_orthogonal(fliplr(res),m_vec),res_catalog);
217.                                 end
218. 
219.                                 % No, this RES is new to us.
220.                                 if(res_fnd == 0)
221.                                     res_catalog{end+1} = res;
222.                                     res_loc{end+1} = [i,j];
223. 
224.                                 % Yes, we have seen this RES before.
225.                                 else
226.                                     G(i,j) == G(res_loc{res_fnd}(1),res_loc{res_fnd}(2));
227.                                 end
228.                             end
229.                         end
230.                     end
231.                     G == semidefinite(real_dim);
232.                 end
233. 
234.                 % Now enforce that P is a probability array and that its
235.                 % marginals are consistent with each other.
236.                 if(isa(p,'cvx'))
237.                     % First, P is an array of probabilities, so it must be
238.                     % non-negative and its faces must sum to 1.
239.                     p >= 0;
240.                     for i = 1:ma
241.                         for j = 1:mb
242.                             sum(sum(p(:,:,i,j))) == 1;
243.                         end
244.                     end
245. 
246.                     % Bob's marginal probabilities must be consistent.
247.                     for i = ob:-1:1
248.                         for j = mb:-1:1
249.                             marg_a{i,j} = sum(squeeze(p(:,i,:,j)),1);
250.                             marg_a{i,j} == marg_a{i,j}(1)*ones(1,ma);
251.                         end
252.                     end
253. 
254.                     % Alice's marginal probabilities must be consistent.
255.                     for i = oa:-1:1
256.                         for j = ma:-1:1
257.                             marg_b{i,j} = sum(squeeze(p(i,:,j,:)),1);
258.                             marg_b{i,j} == marg_b{i,j}(1)*ones(1,mb);
259.                         end
260.                     end
261.                 end
262.         cvx_end
263. 
264.         is_npa = 1-min(cvx_optval,1); % CVX-safe way to map (0,Inf) to (1,0)
265.         if(~isa(p,'cvx')) % make the output prettier if it's not a CVX input
266.             is_npa = round(is_npa);
267. 
268.             % Deal with error messages.
269.             if(strcmpi(cvx_status,'Inaccurate/Solved'))
270.                 warning('NPAHierarchy:NumericalProblems','Minor numerical problems encountered by CVX. Consider adjusting the tolerance level TOL and re-running the script.');
271.             elseif(strcmpi(cvx_status,'Inaccurate/Infeasible'))
272.                 warning('NPAHierarchy:NumericalProblems','Minor numerical problems encountered by CVX. Consider adjusting the tolerance level TOL and re-running the script.');
273.             elseif(strcmpi(cvx_status,'Unbounded') || strcmpi(cvx_status,'Inaccurate/Unbounded') || strcmpi(cvx_status,'Failed'))
274.                 error('NPAHierarchy:NumericalProblems',strcat('Numerical problems encountered (CVX status: ',cvx_status,'). Please try adjusting the tolerance level TOL.'));
275.             end
276.         end
277.     else
278.         is_npa = 1;
279.     end
280. end
281. 
282. 
283. function ind = find_in_cell(val,A)
284.     ind = 0;
285.     for i = 1:numel(A)
286.         if(isequal(A{i},val))
287.             ind = i;
288.             return;
289.         end
290.     end
291. end
292. 
293. % This is a function that computes the next index matrix, given an old
294. % one. Index matrices have 2 rows, and keep track of the measurement
295. % operators that are being multiplied together, from left to right. The
296. % first row contains the index of the measurement, and the second row
297. % contains the index of the result of that measurement. For example,
298. % the index matrix [2 3 3;0 1 4] refers to the product of 3 measurement
299. % operators. The first measurement operator is result 0 of measurement
300. % 2, the second is result 1 of measurement 3, and the third is result 4
301. % of measurement 3. Note that the entries of this matrix are indexed
302. % starting at 0 (i.e., there is a measurement 0 and a measurement
303. % result 0).
304. %
305. % Note that we need this function, rather than the simpler update_odometer
306. % function, because the upper limit in the second row of an index matrix
307. % depends on the value in the first row.
308. function new_ind = update_ind(old_ind,k,m_vec,o_vec)
309.     % Do we have the identity measurement right now? Go to the first
310.     % non-identity one.
311.     if(min(min(old_ind)) == -1)
312.         new_ind = zeros(2,k);
313.         return;
314.     end
315. 
316.     % Start by increasing the last index by 1.
317.     new_ind = old_ind;
318.     new_ind(2,k) = new_ind(2,k)+1;
319. 
320.     % Now we work the "odometer": repeatedly set each digit to 0 if it
321.     % is too high and carry the addition to the left until we hit a
322.     % digit that *isn't* too high.
323.     for l = k:-1:1
324.         % If we've hit the end of the outcomes for this particular
325.         % measurement, move to the next measurement setting.
326.         if(new_ind(2,l) >= o_vec(min(floor(new_ind(1,l)/m_vec(1)),1)+1))
327.             new_ind(2,l) = 0;
328.             new_ind(1,l) = new_ind(1,l) + 1;
329.         else
330.             return; % always return if the odometer doesn't turn over
331.         end
332. 
333.         % If we've hit the end of the measurement settings within
334.         % this particular measurement operator, move to the previous
335.         % measurement operator.
336.         if(new_ind(1,l) >= sum(m_vec))
337.             new_ind(1,l) = 0;
338.             if(l >= 2)
339.                 new_ind(2,l-1) = new_ind(2,l-1) + 1;
340.             else % if L = 1 too then we are completely done! It doesn't matter what we do from here; just exit.
341.                 return;
342.             end
343.         else
344.             return;
345.         end
346.     end
347. end
348. 
349. % This function determines the nature of the operator specified by the
350. % index matrix IND. If IND corresponds to something that is not (generally)
351. % a measurement, then -1 is returned. If it corresponds to the zero
352. % operator, 0 is returned. If it corresponds to the identity operator,
353. % [0;0] is returned. If it corresponds to a measurement then the index
354. % matrix of that measurement is returned.
355. function [res,res_type] = product_of_orthogonal(ind,m_vec)
356.     res_type = -1;
357.     res = ind;
358.     len = size(ind,2);
359. 
360.     % IND is the product of just one measurement operator.
361.     if(len == 1)
362.         res_type = 1;
363.         return;
364. 
365.     % IND is the product of two commuting non-identity measurement
366.     % operators.
367.     elseif(len == 2 && ind(2,1) >= 0 && ind(2,2) >= 0 && min(floor(ind(1,1)/m_vec(1)),1) ~= min(floor(ind(1,2)/m_vec(1)),1))
368.         res = sortrows(res.').'; % sort so that Alice's measurement comes first
369.         res_type = 2;
370.         return;
371.     end
372. 
373.     % IND is more complicated. Recursively figure out how much it can be
374.     % simplified.
375.     for i = 1:len-1
376.         for j = i+1:len
377.             % These two measurements are next to each other and are
378.             % orthogonal!
379.             if(ind(2,i) >= 0 && ind(2,j) >= 0 && ind(1,i) == ind(1,j) && ind(2,i) ~= ind(2,j))
380.                 res_type = 0;
381.                 return; % one is enough; break out of the loop when one is found
382. 
383.             % These two measurements are next to each other and are the
384.             % same! Merge them and then start over.
385.             elseif(ind(1,i) == ind(1,j) && ind(2,i) == ind(2,j))
386.                 [res,res_type] = product_of_orthogonal(ind(:,[1:j-1,j+1:len]),m_vec);
387.                 return;
388. 
389.             % The first of these two measurement operators is the identity.
390.             % Merge them and then start over.
391.             elseif(ind(2,i) == -1)
392.                 [res,res_type] = product_of_orthogonal(ind(:,[1:i-1,i+1:len]),m_vec);
393.                 return;
394. 
395.             % The second of these two measurement operators is the
396.             % identity. Merge them and then start over.
397.             elseif(ind(2,j) == -1)
398.                 [res,res_type] = product_of_orthogonal(ind(:,[1:j-1,j+1:len]),m_vec);
399.                 return;
400. 
401.             % These two measurements act on the same party, but are not
402.             % orthogonal. Stop increasing J and increase I again.
403.             elseif(min(floor(ind(1,i)/m_vec(1)),1) == min(floor(ind(1,j)/m_vec(1)),1))
404.                 break;
405. 
406.             % These two measurements act on different parties and are both
407.             % non-identity. They commute; move Alice's before Bob's
408.             elseif(ind(1,i) > ind(1,j))
409.                 [res,res_type] = product_of_orthogonal(ind(:,[1:i-1,j,i+1:j-1,i,j+1:len]),m_vec);
410.                 return;
411.             end
412.         end
413.     end
414. end