NPAHierarchy

From QETLAB
Revision as of 17:40, 22 January 2015 by Nathaniel (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
NPAHierarchy
Determines whether or not a set of probabilities satisfy the conditions of the NPA hierarchy

Other toolboxes required CVX
Related functions BellInequalityMax
XORGameValue
Function category Nonlocality and Bell inequalities
Usable within CVX? yes (concave)

NPAHierarchy is a function that determines whether or not a given set of probabilities satisfy the constraints specified by the NPA hierarchy[1][2], 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

References

  1. M. Navascués, S. Pironio, and A. Acín. Bounding the set of quantum correlations. Phys. Rev. Lett., 98:010401, 2007. E-print: arXiv:quant-ph/0607119
  2. M. Navascués, S. Pironio, and A. Acín. A convergent hierarchy of semidefinite programs characterizing the set of quantum correlations. New J. Phys., 10:073013, 2008. E-print: arXiv:0803.4290 [quant-ph]