Infeasibility of SDP in MATLAB CVX (Works in Python)

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

You haven’t provided input data and you haven’t shown the CVX and solver output (nor the CVXPY output). That greatly limits the ability of forum readers to assess what’s going on.

If CVX provided the dual to Mosek, then DUAL_INFEASIBLE means the original problem is primal infeasible. But you haven’t shown the CVX output, so I don’t know that that is the case. If it is, all but section 1 of Debugging infeasible models - YALMIP also applies to CVX.

Thank you for your response!
I’ve now included the CVX output, and the input values are d=2,oa=2,ma=2 which I’ve stated above the optimization problem.
The goal is to obtain a ** negative** value of \mu for this specific assemblage:
sigma_a_x(:,:,1,1) = [0.2500, 0.2500; 0.2500, 0.2500];
sigma_a_x(:,:,2,1) = [0.2500, -0.2500; -0.2500, 0.2500];
sigma_a_x(:,:,1,2) = [0.5000, 0.0000; 0.0000, 0.0000];
sigma_a_x(:,:,2,2) = [0.0000, 0.0000; 0.0000, 0.5000];
This would indicate the impossibility of a local hidden state decomposition, aligning with known results. However, while CVXPY (Python) correctly finds a negative \mu, CVX (MATLAB) reports infeasibility.

Have you tried the advice in the link in my previous post?

Yes, I have checked, and the issue seems to be with the first constraint (the unsteerability condition). I tried relaxing it by allowing some tolerance, but the problem only becomes feasible when the tolerance is above 0.1, which is far too large to be acceptable, and even then the value of \mu is not negative.

In addition to constraining sig_lam(:, :, lam) - mu * eye(d) to be semidefinite, you are also constraining sig_lam(:,:,lam) to be semidefinite (by virtue of the variable declaration). If mu is negative, the latter is more constraining than the former. I don’t really understand what you’re doing, so I ll let you follow up on that. In any event, that semidefintie variable declaration is imposing Ndet semidefinte constraints, and they should be considered as such for purposes of infeasibility debugging.

Thanks a lot! You’re right, defining sig_lam as semidefinite was adding unnecessary constraints, which made the problem infeasible when mu was negative. It works now.

Really appreciate your help!