forked from nathanieljohnston/QETLAB
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBellInequalityMax.m
193 lines (160 loc) · 6.76 KB
/
BellInequalityMax.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
%% BELLINEQUALITYMAX Computes the maximum value of a Bell inequality
% This function has three required input arguments:
% COEFFICIENTS: a matrix or tensor that specifies the coefficients of the
% Bell inequality in either full probability, full correlator,
% or Collins-Gisin notation
% DESC: a vector [oa ob ma mb] that specifies the number of outputs and inputs
% of Alice and Bob
% NOTATION: a character vector that specifies the notation the coefficients are
% written in. the possible values are 'fp', 'fc', and 'cg', for
% full probability, full correlator, and Collins-Gisin
%
% BMAX = BellInequalityMax(COEFFICIENTS,DESC,NOTATION) is the
% maximum value that the specified Bell inequality can take on in
% classical mechanics. For the maximum quantum or no-signalling value,
% see the optional arguments described below.
%
% This function has two optional input arguments:
% MTYPE (default 'classical'): one of 'classical', 'quantum', or
% 'nosignal', indicating what type of Bell inequality maximum should
% be computed. IMPORTANT NOTE: if MTYPE='quantum' then only an upper
% bound on the Bell inequality is computed, not its exact value (see
% the argument K below).
% K (default 1): if MYTPE='quantum', then K is a non-negative integer
% or string indicating what level of the NPA hierarchy to use to
% bound the Bell inequality (higher values give better bounds, but
% require more computation time). See the NPAHierarchy function for
% details.
%
% BMAX = BellInequalityMax(COEFFICIENTS,DESC,NOTATION,MTYPE,K) is
% the maximum value that the specified Bell inequality can take on in the
% setting (classical, quantum, or no-signalling) specified by MTYPE.
%
% URL: http://www.qetlab.com/BellInequalityMax
% requires: CVX (http://cvxr.com/cvx/)
% author: Nathaniel Johnston ([email protected]) and Mateus Araújo
% package: QETLAB
% last updated: September 16, 2022
function bmax = BellInequalityMax(coefficients,desc,notation,varargin)
% set optional argument defaults: MTYPE='classical', K=1
[mtype,k] = opt_args({ 'classical', 1 },varargin{:});
% Get some basic values
oa = desc(1);
ob = desc(2);
ma = desc(3);
mb = desc(4);
% The no-signalling maximum is just a LP. We use the Collins-Gisin notation,
% so the only constraint left to impose is positivity of the probabilities
if(strcmpi(mtype,'nosignal'))
% Set up the Bell inequality.
if (strcmpi(notation,'cg'))
M = coefficients;
elseif (strcmpi(notation,'fp'))
M = fp2cg(coefficients);
elseif (strcmpi(notation,'fc'))
M = fc2cg(coefficients);
end
cvx_begin quiet
variable p_cg((oa-1)*ma+1,(ob-1)*mb+1);
p_fp = cg2fp(p_cg,desc,1);
bmax = M(:)'*p_cg(:);
maximize bmax;
subject to
p_cg(1,1) == 1;
p_fp >= 0;
cvx_end
% Compute the quantum maximum value of the Bell inequality
elseif(strcmpi(mtype,'quantum'))
% Set up the Bell inequality.
if (strcmpi(notation,'cg'))
M = coefficients;
elseif (strcmpi(notation,'fp'))
M = fp2cg(coefficients);
elseif (strcmpi(notation,'fc'))
M = fc2cg(coefficients);
end
cvx_begin quiet
variable p((oa-1)*ma+1,(ob-1)*mb+1);
bmax = M(:)'*p(:);
maximize bmax;
subject to
p(1,1) == 1;
NPAHierarchy(p,desc,k) == 1;
cvx_end
% Deal with error messages.
if(strcmpi(cvx_status,'Inaccurate/Solved'))
warning('BellInequalityMax:NumericalProblems','Minor numerical problems encountered by CVX. Consider adjusting the tolerance level TOL and re-running the script.');
elseif(strcmpi(cvx_status,'Inaccurate/Infeasible'))
warning('BellInequalityMax:NumericalProblems','Minor numerical problems encountered by CVX. Consider adjusting the tolerance level TOL and re-running the script.');
elseif(strcmpi(cvx_status,'Unbounded') || strcmpi(cvx_status,'Inaccurate/Unbounded') || strcmpi(cvx_status,'Failed'))
error('BellInequalityMax:NumericalProblems',strcat('Numerical problems encountered (CVX status: ',cvx_status,'). Please try adjusting the tolerance level TOL.'));
end
elseif(strcmpi(mtype,'classical'))
% Compute the classical maximum value by going through all of Bob's strategies and picking
% Alice's corresponding optimal strategy
if oa == 2 && ob == 2 %if we have only 2 outcomes convert to full correlation notation to use the faster algorithm
if (strcmpi(notation,'fc'))
M = coefficients;
elseif (strcmpi(notation,'fp'))
M = fp2fc(coefficients);
elseif (strcmpi(notation,'cg'))
M = cg2fc(coefficients);
end
if (ma < mb) %if Alice has fewer inputs than Bob we swap them
M = M';
ma = desc(4);
mb = desc(3);
end
if (mb <= 21) %for few inputs it's faster not to parallelize
parallel_threads = 0;
else
parallel_threads = Inf;
end
Bmarginal = M(1,2:end);
Amarginal = M(2:end,1);
correlations = M(2:end,2:end);
bmax = -Inf;
parfor (b=0:2^mb-1,parallel_threads)
b_vec = 1-2*integer_digits(b,2,mb);
temp_bmax = Bmarginal*b_vec + sum(abs(Amarginal + correlations*b_vec));
bmax = max(bmax,temp_bmax);
end
bmax = bmax + M(1,1);
else
if (strcmpi(notation,'fp')) %if we have more than 2 output we convert to full probability notation
M = coefficients;
elseif (strcmpi(notation,'cg'))
M = CG2FP(coefficients,desc);
end
if (oa^ma < ob^mb) % we choose Bob as the party with the fewest strategies
M = permute(M,[2,1,4,3]);
[oa ob ma mb] = size(M);
end
M = permute(M,[1 3 2 4]);
M = reshape(M,oa*ma,ob*mb);
offset = 1+ob*[0:mb-1]';
if (ob^mb <= 10^6) %for few strategies it's faster not to parallelize
parallel_threads = 0;
else
parallel_threads = Inf;
end
bmax = -Inf;
parfor (b=0:ob^mb-1,parallel_threads)
ind = integer_digits(b,ob,mb) + offset;
Ma = sum(M(:,ind),2);
temp_bmax = sum(max(reshape(Ma,oa,ma)));
bmax = max(bmax,temp_bmax);
end
end
else
error('BellInequalityMax:InvalidMTYPE','MTYPE must be one of ''classical'', ''quantum'', or ''nosignal''.');
end
end
% converts number "number" to base "base" with digits "digits"
function dits = integer_digits(number,base,digits)
dits = zeros(digits,1);
for i=1:digits
dits(digits+1-i) = mod(number,base);
number = floor(number/base);
end
end