Yalmip - Baron error

[Aras]

Hi Dr. Lofberg,

I have a convex problem which looks like:

w = sdpvar(n,1);
w0 = sdpvar(1,1);
Constraints = [w>=0, sum(w)==1];
Objective = - (A(:,i).'*w - D(:,i).'*(u + V*w + r*(entropy(w))));
optimize(Constraints,Objective,sdpsettings('debug',1,'solver','baron'));

The data I use is hard to share here so I will try to explain my problem. Knitro, or IPOPT works for this optimization. However, Baron gives the following error:

Index exceeds the number of array elements (0).


Error in baron


Error in baron


Error in baron


Error in callbaron (line 85)
[x,fval,exitflag,info,allsol] = baron(obj,A,rl,ru,lb,ub,con,cl,cu,xtype,x0,opts);


Error in solvesdp (line 361)
    eval(['output = ' solver.call '(interfacedata);']);


Error in optimize (line 31)
[varargout{1:nargout}] = solvesdp(varargin{:});

Can this be an issue about YALMIP? Thank you for your time.

[Johan Löfberg]

why would you use baron for a convex problem?

but no, this appears to be a baron issue, but since you do not supply a reproducible example, we cannot know (just use random but similiar data)

[Erling D. Andersen]

Warning: Shameless advetisement ahead.

This can be formulated using the exponential cone and hence solved by Mosek v9.

See

   https://docs.mosek.com/cheatsheets/conic.pdf 

I would be surprised if that does not beat the other optimizers.

[Johan Löfberg]

Absolutely, YALMIP automatically selects mosek as first choice solver here (assuming r has the correct sign)

[Aras]

Dear Both,

Thank you for your nice comments.

I was just benchmarking all the solvers and since the objective is convex (linear constraints) I was hoping BARON to work ok. Since it is the only solver which does not work I specifically wanted to ask about your opinion. I am going to share the data soon (I am trying to reduce the size for the easy debugging). I hope this will not take your time, otherwise, I can skip BARON.

Also, I never reformulated entropy functions as exponential cones. That would be a very nice method for me, and I am now trying! Thank you!

[Johan Löfberg]

YALMIP does the exponential cone formulation for you automatically if you have mosek 9 installed

I suspect there is a bug in baron. This fails on my installation

K>> obj = @(x)(x(1)-(-x(1)*log(x(1))-x(2)*log(x(2))))

K>> baron(obj,[],[],[],[],[],[],[],[],['CC'])

while this works

K>> obj = @(x)(-(-x(1)*log(x(1))-x(2)*log(x(2)))+x(1))

K>> baron(obj,[],[],[],[],[],[],[],[],['CC'])

[Johan Löfberg]

..although that error doesn't appear in the latest version of baron after updating the matlab interface, so still interested in a reproducible example for your problem

[Aras]

I attached the data. I tried with Mosek 9, works very well :) So I am still interested in the reason of BARON error. Can you please tell me if you have an idea? Many thanks.

[Johan Löfberg]

runs without problems here, so you are experiencing a bug in the baron binary running on your machine

having said that, solving this problem using baron would be crazy inefficient

[Aras]

Ah, that's interesting. I had downloaded and installed BARON's newest version and except this problem it works in all problems. I will try to fix this issue, good to hear it doesn't crash at your PC.

And yes, I am sticking to convex solvers :) Observed already that MOSEK is the fastest.

However, why do you think BARON is inefficient? Doesn't it reduce to a local solver when the problem is convex?

[Johan Löfberg]

no,  it runs the whole machinery for all problems to close the gap from below, it simply assumes it i a hard nonconvex problem as you never would use it otherwise

[Aras]

That's so good to know, thanks.

[Mark L. Stone]

Section 4,4 of the BARON manual (2019.3.22) https://minlp.com/downloads/docs/baron%20manual.pdf

CONVEX_EQUATIONS <list equation names>;
This equation declaration can be used to specify constraints that are convex. This is optional and must follow after the EQUATIONS declaration and before the equation definitions.

Using epigraph formulation, convex objective function could de declared (as a constraint) as well (I'm pretty sure BARON will recognize a linear objective as convex)

As far as I know, this functionality is not available under the MATLAB interface.

-----------------------

Aras, are you using matbar 1.88, which is the latest version of the MATLAB BARON interface?

[Aras]

Hi Dr. Stone,

In matbar folder version.txt says v1.88, so I guess I am using the latest version. 

So, I am concluding that unless I can make CONVEX_EQUATIONS work there is no sense to use BARON to benchmark a convex problem. For general convex programming, I am currently benchmarking IPOPT, Knitro, MOSEK right now.

By the way, IPOPT works a bit worse than the others with the entropy functions. I think this is about the algorithm of IP when the variable gets close to 0. Knitro was the fastest, but now MOSEK 9 works very well too :)

[Johan Löfberg]

Even if you did that, it would still make no sense, as you would time Barons internal nonlinear local solver (which ever that is) + a lot of overhead. It would be similar to specifying the solver as YALMIP global solver bmibnb as a very expensive way to time, e.g., fmincon.

Note that when you time fmincon/ipopt etc via YALMIP, you are timing the time it takes for these solvers to solve the problem as coded by YALMIP. This means for instance, that the solvers have no Hessian information as yalmip doesn't support this. ipopt suffers enormously from this, and your model has a trivial hessian, so a manually coded call to ipopt would most likely lead to much faster solves.

[Aras]

I got more information than I could ever imagine under one forum question. Thanks :)

[Erling D. Andersen]

How big are your problems i.e. how long time does take it take to solve a problem for Mosek and Knitro.

I am bit surprised Knitro should be faster. Would you be willing to give us a problem.

[Johan Löfberg]

I don't see knitro being the fastest. It is a very small problem, but mosek appears 20x faster (mosek 9, knitro 11, once again with the caveat that knitro is called from yalmip and thus not necessarily done in the best way)

>> optimize(Constraints,Objective,sdpsettings('solver','knitro','verbose',0))

ans = 

  struct with fields:

    yalmipversion: '20190425'

       yalmiptime: 1.559811467071738e-01

       solvertime: 9.301885329282590e-02

             info: 'Successfully solved (KNITRO)'

          problem: 0

>> optimize(Constraints,Objective,sdpsettings('solver','mosek','verbose',0))

ans = 

  struct with fields:

    yalmipversion: '20190425'

       yalmiptime: 1.066808792425237e-01

       solvertime: 4.319120757476922e-03

             info: 'Successfully solved (MOSEK)'

          problem: 0

[Erling D. Andersen]

In any case it is hard to do reliable timing for such small problems. 

[Johan Löfberg]

btw, just so you don't think anything else, entropy(w) is -sum(w*log(w)) (so you don't think it is an elementwise operator)

[Johan Löfberg]

Yes. My guess is that Aras hasn't turned off display, so basically compared the iteration display functions of knitro and mosek (and still, the time 4.319120757476922e-03 in mosek is so low that it basically is down to the time it takes matlab to make a call to mosek. Solving an LP with 1 variable and 1 constraint takes 1.3e-3, which thus essentially is the call overhead

>> sdpvar x

>> optimize(x>=0,x,sdpsettings('solver','mosek','verbose',0))

ans = 

  struct with fields:

    yalmipversion: '20190425'

       yalmiptime: 1.037132350146018e-01

       solvertime: 1.286764985395122e-03

[Erling D. Andersen]

Logging is very expensive in MATLAB and should be turned off in serious benchmarks.

[Aras]

Hi both,

I don't think Knitro is the fastest. I was not using MOSEK before for this particular case. However please see my full code below:

load('upper_soltn.mat')

w = sdpvar(n,1);
solver_time= [];

Constraints = [w>=0, sum(w)==1];

tic
%Objective = -(d'*u + d'*V*w - (1+d'*r)*sum(w.*log(w)));
Objective = -(d'*u + d'*V*w + max(0,(1+d'*r))*entropy(w) - tau);
P = optimize(Constraints,Objective,sdpsettings('solver','mosek','verbose',0));
solver_time = [solver_time P.solvertime];
wsave = value(w);
%yalmip clear;

for i=1:n

w = sdpvar(n,1);

Constraints = [w>=0, sum(w)==1];

%Objective = - (A(:,i).'*w - D(:,i).'*(u + V*w + r*(-sum(w.*log(w)))));

Objective = - (A(:,i).'*w - D(:,i).'*(u + V*w + r*(entropy(w))));

P = optimize(Constraints,Objective,sdpsettings('solver','mosek','verbose',0));
solver_time = [solver_time P.solvertime];
wsave = [wsave value(w)];
yalmip clear;
end

toc

save('save_w', 'wsave')

for i=1:q
w = sdpvar(n,1);

Constraints = [w>=0, sum(w)==1];

%Objective = u(i) + V(i,:)*w + r(i)*(-sum(w.*log(w)));
Objective = u(i) + V(i,:)*w + r(i)*(entropy(w));
P = optimize(Constraints,Objective,sdpsettings('solver','mosek','verbose',0));
solver_time = [solver_time P.solvertime];
wsave = [wsave value(w)];
yalmip clear;
end
toc
save('save_w', 'wsave')

So I am solving similar problems n+q+1 (20 + 20 + 1) times. The overall tic-toc time takes 6 secs for Knitro and 5.3 secs for MOSEK. I am aware that it may take some extra time for Yalmip to convert problems to the exponential conic formulation for Mosek, so further I check the solver_time. The mean P.solvertime is 0.0024 for MOSEK,  and 0.0195 for Knitro. Also, the summation of solver times is 0.0981 to 0.7980. This shows that MOSEK in my problem is faster.

I turned off the logging, I was screening the outcomes. The reason is, my problem is a scalable adjustable robust optimization problem. This was the n=20 case, however if I can obtain my 'fastest' form, then I will solve n=100,200,500 ... cases.

Finally, I didn't get this sentence: "btw, just so you don't think anything else, entropy(w) is -sum(w*log(w)) (so you don't think it is an elementwise operator)". I believe -sum(w.*log(w))) and entropy(w) give the same result. So I am using a higher level operator as you suggested in another post.

Ps: I am aware that I don't need to re-define 

[Aras]

Ps: I am aware that I don't need to re-define Constraints and variable each time. I will change these :)

 

And, I used max(0,(1+d'*r)) instead of (1+d'*r) because the value of the latter is  -5.6428e-07, and I obtain this from another model which I use MOSEK with a constraint of 1+d'*r >= 0. Clearly, this is close to 0 and does not make an issue in the way I obtain it. However, in the script I shared above, MOSEK doesn't accept the convexity (so no rounding to zero since it is a coefficient). Thus, I made it max(0,(1+d'*r)) but I am really open for suggestions for these kind of problems... 

[Johan Löfberg]

yes, regarding that comment, I misread your code, thinking it indicated you though entropy(x) was a vector

[Aras]

Also Dr. Andersen, when I said "Knitro was the fastest, but now MOSEK 9 works very well too" I meant before trying MOSEK. Now, as I showed above, MOSEK works the fastest.

[Erling D. Andersen]

:-).