sdp and socp formulation

[Peter]

I have an issue with the SOCP and SDP formulation


SOCP formulation
x=sdpvar(n,1);
obj=(v'*x);
f= [a'*x - b  >= c*norm(sqrt(E)*x,2)];
f=f+[ones(n,1)'*x==1.2];
f=f+[-.1<=x<=1.3];

SDP fomulation
x=sdpvar(n,1);
obj=(v'*x);
X=sdpvar(n,n,'di');
f1= [trace(X*a*a') - 2*a'*x*b + b^2 >= (c^2)*trace(X*E) ];
f1=f1+[ones(n,1)'*x==1.2];
f1=f1+[-.1<=x<=1.3];
f1=f1+[[1 x'; x X]>=0];


b and c are known constants
E is a diagonal matrix
a and v are known vectors

The objective value is different in both the cases.
Tried without the 'di' in X but there is no change.
Is there any fault in my formulation ?


Regards,
Peter

[Johan Löfberg]

If I understand you correctly, you are applying an SDP relaxation to the following model, which you think is equivalent to the original constraint 

(a'*x - b)^2  >= c^2*x'*E*x

To begin with, I hope you understand the SDP formulation is vastly inefficient, since you're trying to lift a convex SOCP to a semidefinite relaxation, thus introducing a large number of variables which you really don't need. The SOCP can much easier be reformulated to the folllowing SDP by simply applying a Schur complement (if you really have to reformulate as an SDP, typically you solve these problems using SOCP capable solvers such as SeDuMi, SDPT3, Mosek, CPLEX etc)

[(a'*x - b)*eye(n) c*sqrt(E)*x; x'*sqrt(E)*c (a'*x - b)] >=0

Without looking any deeper at your relaxation, I would suspect your issues are related to the fact that a'*x-b has to be non-negative. I.e., the correct equivalence is

[a'*x - b  >= c*norm(sqrt(E)*x,2)] <=> [(a'*x-b)^2 >= c^2*x'*E*x, a'*x-b >=0]

[Johan Löfberg]

Additionally, even if you would fix this, the first SDP relaxation is not tight, i.e., for this model you need a higher order relaxation. You can see that by using the built-in module for semidefinite relaxation

x=sdpvar(n,1);

obj=(v'*x);

f= [a'*x - b  >= c*norm(sqrt(E)*x,2)];

f=f+[ones(n,1)'*x==1.2];

f=f+[-.1<=x<=1.3];

solvesdp(f,obj)

% True solution

double(obj)

% Solve SDP relaxation instead

Model = [(a'*x-b)^2 >= c^2*x'*E*x, a'*x-b>=0,sum(x)==1.2,-.1 <= x <= 1.3];

solvemoment(Model,obj,sdpsettings('moment.order',1));

% Standard SDP relaxation is not tight, lower bound much lower than true

relaxdouble(obj)

% First higher-order (include quartics) is almost tight (albeit with numerical problems)

solvemoment(Model,obj,sdpsettings('moment.order',2))

relaxdouble(obj)
% Include 6th order monomials in relaxation yields tight bound

solvemoment(Model,obj,sdpsettings('moment.order',3))

relaxdouble(obj)