I am solving an SDP problem in MATLAB using CVX to check for the existence of a local hidden state decomposition of a given assemblage sigma_a|x(d,d,oa,ma) (a 4-D array where first two dimensions contain the d x d quantum states, while the last two dimensions are outcomes and measurements)
The same problem in Python using CVXPY has feasible solutions, but I need it to be working in MATLAB. And somehow in MATLAB, the SDP is feasible if I set:
sigma_a_x = repmat(eye(d)/4, [1, 1, oa, ma]) but otherwise, it becomes infeasible, which is unexpected since the problem is strictly feasible for any valid sigma_a|x.
Problem status : DUAL_INFEASIBLE
Why does the MATLAB version become infeasible?
d=2; % dimension for Alice and Bo
ma = d; % number of measurement settings for Alice
oa = d; % number of outcomes of Alice for each setting
Ndet = oa^ma; % number of deterministic strategies for Alice
SingleParty = genSinglePartyArray(oa,ma); % generate array containing single party strategies
Optimization Problem:
cvx_begin sdp
variable sig_lam(d, d, Ndet) hermitian semidefinite
variable mu
maximize mu
subject to
% unsteerability condition
for x = 1:ma
for a = 1:oa
sigma_sum = cvx(0);
for lam = 1:Ndet
sigma_sum = sigma_sum + SingleParty(a, x, lam) * sig_lam(:, :, lam);
end
sigma_sum - sigma_a_x(:, :, a, x) == semidefinite(d,d);
end
end
% normalization constraint on sig_lam
trace_sum = cvx(0);
for lam = 1:Ndet
trace_sum = trace_sum + trace(sig_lam(:,:,lam));
end
trace_sum == 1;
% positivity constraint on sig_lam
for lam = 1:Ndet
sig_lam(:, :, lam) - mu * eye(d) == semidefinite(d,d);
end
cvx_end
cvx output
For improved efficiency, Mosek is solving the dual problem.
MOSEK Version 9.1.9 (Build date: 2019-11-21 11:27:13)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Linux/64-X86
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 17
Cones : 12
Scalar variables : 49
Matrix variables : 0
Integer variables : 0
Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 1 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.00
Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 17
Cones : 12
Scalar variables : 49
Matrix variables : 0
Integer variables : 0
Optimizer - threads : 10
Optimizer - solved problem : the primal
Optimizer - Constraints : 17
Optimizer - Cones : 13
Optimizer - Scalar variables : 50 conic : 50
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 0.00 dense det. time : 0.00
Factor - ML order time : 0.00 GP order time : 0.00
Factor - nonzeros before factor : 137 after factor : 153
Factor - dense dim. : 0 flops : 2.83e+03
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 3.5e+00 1.4e+00 1.0e+00 0.00e+00 -2.000000000e+00 0.000000000e+00 1.0e+00 0.00
1 1.0e+00 4.1e-01 9.5e-02 1.23e+00 -4.238967825e-01 -8.684944020e-02 2.9e-01 0.00
2 2.6e-01 1.1e-01 1.7e-02 9.09e-01 -2.549009470e-01 -1.305726897e-01 7.6e-02 0.00
3 7.0e-02 2.8e-02 1.1e-02 -4.09e-01 -1.388270035e+00 -9.519483161e-01 2.0e-02 0.00
4 4.3e-04 1.7e-04 1.0e-03 -8.64e-01 -3.252340044e+02 -2.588986352e+02 1.2e-04 0.00
5 1.7e-06 6.8e-07 6.3e-05 -9.99e-01 -8.326647970e+04 -6.633433373e+04 4.8e-07 0.00
6 9.6e-12 1.0e-11 1.4e-10 -1.00e+00 -2.476368758e-01 -1.972807881e-01 1.1e-12 0.00
Optimizer terminated. Time: 0.00
Interior-point solution summary
Problem status : DUAL_INFEASIBLE
Solution status : DUAL_INFEASIBLE_CER
Primal. obj: -2.4763687582e-01 nrm: 3e+00 Viol. con: 2e-11 var: 0e+00 cones: 2e-23
Optimizer summary
Optimizer - time: 0.00
Interior-point - iterations : 6 time: 0.00
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00
Status: Infeasible
Optimal value (cvx_optval): -Inf
SingleParty function
function SingleParty = genSinglePartyArray(oa,ma)
Ndet = oa^ma; %number of deterministic behaviours
SingleParty = zeros(oa,ma,Ndet); % initialise array
for lam = 0:Ndet-1
lamdec = dec2base(lam,oa,ma)-‘0’; % generates the string of outcomes a
%(for each x), for the given variable lam
for x = 1:ma
for a = 0:oa-1
SingleParty(a+1,x,lam+1) = (lamdec(x) == a);
% probability = 1 if a = lamdec(x), 0 otherwise
end
end
end
end