Optimizer and solvesdp give different results when using optimizer with sedumi, but not with gurobi?

[Rasmus Brandt]

Hello,

I have stumbled upon some weird behaviour with optimizer in conjunction with sedumi. The following code will yield different results from optimizer and solvesdp, when using sedumi. When using gurobi, the results are the same! Am I doing something wrong here, or is this a bug?

I guess that my code uses the beta version of general variable substitution in optimizer, but both 'gurobi' and '+gurobi' seem to work equally fine.

Great work with Yalmip otherwise!

best,

Rasmus

==

% Parameters, and generate input data

Mt = 3;

Nd = 2;

Pt = 100;

alpha = 5/100;

Wsqrtm = crandn(Nd,Nd); W = Wsqrtm*Wsqrtm';

G = crandn(Mt,Nd);

t = crandn(Mt,Mt); T = t*t' + G*G'; Tsqrtm = sqrtm(T);

Ssqrtm = crandn(Mt,Mt); S = Ssqrtm*Ssqrtm';

% Solve in standard way

Vs = sdpvar(Mt,Nd,'full','complex');

ctsq = sdpvar(Mt,1);

constraints = [norm(Vs,'fro')^2 + sum(ctsq) <= Pt];

for m = 1:Mt

constraints = constraints + [ctsq(m) >= (alpha^2)*norm(Vs(m,:),'fro')^2];

end

options = sdpsettings('solver','sedumi','verbose',1);

objective = trace(W) - 2*real(trace(Wsqrtm*G'*Vs)) + norm(Tsqrtm*Vs,'fro')^2  + diag(T)'*ctsq + (alpha^2)*norm(Ssqrtm*Vs,'fro')^2;

solvesdp(constraints,objective,options);

double(Vs)

% Solve using optimizer

Vs_opt = sdpvar(Mt,Nd,'full','complex');

ctsq_opt = sdpvar(Mt,1);

constraints_opt = [norm(Vs_opt,'fro')^2 + sum(ctsq_opt) <= Pt];

for m = 1:Mt

constraints_opt = constraints_opt + [ctsq_opt(m) >= (alpha^2)*norm(Vs_opt(m,:),'fro')^2];

end

W_in = sdpvar(Nd,Nd,'full','complex');

Wsqrtm_in = sdpvar(Nd,Nd,'full','complex');

G_in = sdpvar(Mt,Nd,'full','complex');

Tsqrtm_in = sdpvar(Mt,Mt,'full','complex');

T_in = sdpvar(Mt,Mt,'full','complex');

Ssqrtm_in = sdpvar(Mt,Mt,'full','complex');

options_opt = sdpsettings('solver','sedumi','verbose',2);

objective_opt = trace(W_in) - 2*real(trace(Wsqrtm_in*G_in'*Vs_opt)) + norm(Tsqrtm_in*Vs_opt,'fro')^2  + diag(T_in)'*ctsq_opt + (alpha^2)*norm(Ssqrtm_in*Vs_opt,'fro')^2;

solver = optimizer(constraints_opt,objective_opt,options_opt,{W_in,Wsqrtm_in,G_in,Tsqrtm_in,T_in,Ssqrtm_in},Vs_opt);

optimizer_solution = solver{{W,Wsqrtm,G,Tsqrtm,T,Ssqrtm}};

double(optimizer_solution)

[Johan Löfberg]

I get the same objective using both approaches with SeDuMi. The problem doesn't have a unique minimizer, and any small change can lead to a different solution. Gurobi simply happens to be more stable with respect to changes in the model setup.

[Rasmus Brandt]

The problem should have a unique minimizer, this can be seen by solving it analytically. I have created a smaller example, where the same problem arises for sedumi and sdpt3, but not gurobi. See results below and attached code. I tried the code both with Release 20131002, as well as the latest checkin from github. Both exhibit the same behaviour.

Two other things I noticed:

1) When removing the constraints in the attached code, optimizer fails with error "Nonlinear replacement in optimizer object only supported in MATLAB R2012A or later", even though I am running Matlab 2012a.

2) When having several complex optimization variables being output from optimizer, they are returned as real matrices with doubled first dimension. This does not happen when asking for only one complex optimization variable to be returned from optimizer.

best,

Rasmus

=====

Sedumi

--

SOLVESDP

Objective value

   -0.5000

Objective value, relative error

   2.5074e-11

Solution

    0.1944

    0.0278

   -0.1389

Solution, relative error

   5.0740e-06

--

OPTIMIZER

Objective value

   4.4537e+03

Objective value, relative error

   8.9085e+03

Solution

    8.0178

    5.3452

    2.6726

Solution, relative error

   41.0222

--

SDPT3

--

SOLVESDP

Objective value

   -0.5000

Objective value, relative error

   1.4154e-10

Solution

    0.1944

    0.0278

   -0.1389

Solution, relative error

   1.1979e-05

--

OPTIMIZER

Objective value

   4.4537e+03

Objective value, relative error

   8.9085e+03

Solution

    8.0178

    5.3452

    2.6726

Solution, relative error

   41.0222

--

Gurobi

--

SOLVESDP

Objective value

   -0.5000

Objective value, relative error

   2.1010e-09

Solution

    0.1945

    0.0278

   -0.1389

Solution, relative error

   5.6557e-05

--

OPTIMIZER

Objective value

   -0.5000

Objective value, relative error

   2.1010e-09

Solution

    0.1945

    0.0278

   -0.1389

Solution, relative error

   5.6557e-05

--

[Johan Löfberg]

Now I see the problem better.

The problem here is that you have high-level constructs which aren't reformulated when data is fixed. The objective is quartic and for fixed data it simplifies to a quadratic. However, the optimizer layer does not expect unsupported stuff to show up in the instantiated model (sedumi does not support quadratic objectives), hence, it disappears. You should reformulate your model using the low-level cone operator instead (i.e., instead of using (T*x-G)'*(T*x-G) as the objective, introduce a new variable t, and define e = T*x-G;e = e(:); and add the cone constraint cone(e,t). This involves a quadratic expression inside the cone operator, but for fixed T, it is a simple SOCP constraint which is handled by SeDuMi

[Rasmus Brandt]

Great, now at least the simple example seems to work fine with sedumi!

Is it generally a good idea to use the 'cone' command, viz. norms, if I will end up using gurobi as the solver, as well?

Thanks,

// Rasmus

[Johan Löfberg]

Generally, you never use the cone command since it forces you to do the modelling manually, which is against the whole idea of a modelling language.

In the optimizer scenario though, when the expression isn't directly in a form suitable for the targeted solver, you have to use the cone operator. That is the price you pay for speed.

Hence, (decisionvariable*parameter-constant)^2 has to be modeled using cone in optimizer, since the expression not is SOCP representable before parameters are fixed, hence YALMIP does not setup a SOCP structure, and simply assumes that the whole expression will boil down to something linear when parameter is fixed.

(decisionvariable*constant - parameter)^2 can be modelled using standard norm (or simply using the square as here) since it is SOCP representable with unfixed parameters, and YALMIP thus precompiles an SOCP structure in which data can be fixed later on.

[Rasmus Brandt]

I made a coding error, so the simple example used numerical matrices, not sdpvars... So I'm still not able to get it to work properly, even when using the cone construct for the "(decisionvariable*parameter-constant)^2" case. (In the end, I need to do  "(decisionvariable*parameter-parameter)^2" though). 

For the simple example, I'm getting the following error, which I don't really understand:

==

Error using error

Not enough input arguments.

Error in lmi/categorizeproblem (line 252)

                                error

Error in compileinterfacedata (line 251)

[ProblemClass,integer_variables,binary_variables,parametric_variables,uncertain_variables,semicont_variables,quad_info]

= categorizeproblem(F,logdetStruct,h,options.relax,parametric,evaluation_based,

Error in export (line 136)

    [interfacedata,recoverdata,solver,diagnostic,F] =

    compileinterfacedata(F,[],logdetStruct,h,options,findallsolvers);

Error in optimizer (line 185)

    [aux1,aux2,aux3,model] = export(set(x ==

    repmat(pi,nIn*mIn,1))+Constraints,Objective,options,[],[],0);

Error in yalmip_export_error (line 18)

solver = optimizer(constraints_opt,objective_opt,options_opt,Tsqrtm_in,Vs_opt);

==

when trying to do the following

==

% Data

Mt = 3;

G = [3;2;1];

t = [3 2 1; 1 3 2; 2 1 3]; T = t*t' + G*G'; Tsqrtm = sqrtm(T);

 

% Optimization variables

Vs_opt = sdpvar(Mt,1);

t_opt = sdpvar(1);

 

% Input

Tsqrtm_in = sdpvar(Mt,Mt);

 

% Set up

objective_opt = t_opt;

e = Tsqrtm_in*Vs_opt - G; constraints_opt = cone(e,t_opt);

options_opt = sdpsettings('solver','sedumi','verbose',2);

 

solver = optimizer(constraints_opt,objective_opt,options_opt,Tsqrtm_in,Vs_opt);

Vs_opt_sol = solver{Tsqrtm};

===

[Johan Löfberg]

OK, see where it fails. Easiest is to move all nonlinearities to an equality instead

e = sdpvar(3,1);
constraints_opt = [e == Tsqrtm_in*Vs_opt - G,cone(e,t_opt)];

[Johan Löfberg]

FYI, the discovered bugs have all been fixed, and the code now crashes if quadratic objectives are discovered after fixing parameters, when using a solver which doesn't support quadratic objectives.

[Rasmus Brandt]

Ok, that sounds great.

I also managed to get my non-simplified optimization problem working fine with the equality constraint trick. Using 'optimizer', I am seeing speedups on the order of 20x compared to my original CVX implementation. This speedup really saves me, since the full simulation will probably only take around 10-20 hours to run now!

Thanks for the excellent support!

hälsningar

Rasmus, KTH signalbehandling