Convex optimization with arbitrary precision

[O S]

Dear all,

I'm new to YALMIP and will be happy for your help.

The optimization I have involves a convex objective function that involves relative entropy and linear terms, the constraints are linear.

Several other solvers gave me results that are pretty good in terms of tolerance (~e^-12), but I need further precision for my research problem.

Is it possible to set a tolerance smaller than e^-12 in YALMIP?

Thanks a lot,

Oron

[Johan Löfberg]

YALMIP is not a solver but a modelling language.

However, solvers typically cannot work that that high precision. In almost all cases where someone requires this kind of precision, it is a sign of a badly scaled model/ or generally flawed model/idea/methodology.

You will probably have to develop your own solver using higher precision floating-point libraries etc to get this.

Maybe Mark has some insight on this, I know he's interested in these issues.

[Johan Löfberg]

BTW, I read that as 10^-12

e^-12 is 6*10^-6 which is pretty standard precision, and solvers typically go perhaps a couple of digits further.

[Johan Löfberg]

As a side comment, you've posted this on the gurobi forum, and there you referenced cvx. Perhaps you've used cvx successive approximation schemes then to solve the problem. That is really not the way to solve the problem if you want high precision. Simply solve it using a standard nonlinear solver (such as fmincon/ipopt/knitro/ etc that you can call from YALMIP) and see what kind of precision you get with a real nonlinear solver. There are various entropy-operators you can use to improve the computations of gradient callbacks etc, instead of simply writing the exp/log functions https://yalmip.github.io/tags/#exponential-and-logarithmic-functions, and if your problem is an exponential cone progra, these operators will be applicable also for exponential cone solvers

[Mark L. Stone]

quadMINOS (feasibility and optimality tolerance of 1e-17) or quadSNOPT , but I don't believe under MATLAB.  I've never used them. Quad precision will be in software, so not fast.
https://web.stanford.edu/group/SOL/talks/14informsQuadMINOS.pdf
https://web.stanford.edu/group/SOL/talks/16meridaQuadMINOS.pdf

[Joachim Dahl]

How do you assert that your current solution is not accurate enough?  Is it because you encounter numerical problems in the solution process, or does the solver converge and you then check that the constraints are not satisfied with high accuracy?

As an aside, the up-coming  version of MOSEK 9.0 handles exponential cones natively,  but only with standard double-precision. 

[O S]

Dear Joachim, 

There are instances of problems that I could solve analytically and therefore I know what is the exact optimal solution.

If the the solvers within the YALMIP supporting solvers will not suffice, I will try the new MOSEK when it is released, thanks!

[O S]

Thanks Johan,

Formalizing the problem using the YALMIP language was quite simple and the problem has been 'successfully solved.' after setting an initial guess using the 'assign' function. 

Is it possible to change the stopping criteria or the tolerance of the fmincon solver through the YALMIP, so I can get better accuracy for the result? I added a figure with the results I got and I will be happy to learn more about the accuracy of the solution I got.

One more question, how can I extract the optimization results (optimal value, optimal arguments etx.) as variables? When I write R = optimize(...) it returns a structure including yalmiptime, solvertime, info and problem without the variables I mentioned.

Thanks again!

[Johan Löfberg]

All explained in the very first example here

https://yalmip.github.io/tutorial/basics/

[O S]

Thank you Johan,

following your suggestion, I implemented the YALMIP code in order to check the results accuracy:

1. The fmincon gave me a result which is not good compared to the sedumi (using CVX).

2. The IPOPT solver returns "EXIT: Problem has too few degrees of freedom." Do you have any idea why is that? (my problem is formulated well since the fmincon gave me the correct result +- eps)

3. Additionally, are there more solvers that I can try using the YALMIP for my problem that involves the relative entropy function?

 

[Johan Löfberg]

if you've solved an exponential problem using cvx approximations and sedumi, and nonlinear solvers such as fmincon performs worse, it is really odd and you should post your code.

[O S]

You are correct, I measured the wrong quantity (the tolerance that the CVX and the YALMIP returned), while the absolute error (with respect to the analytic optimum which I could derive for this instance) of the fmincon is 2e-8 and the CVX error is 8e-8 so the fmincon performs better.

This takes me back to my initial question on which solvers might provide better accuracy for an exponential problem.

1. I will be happy to know if the fmincon accuracy can be improved.

2. When I try to use the IPOPT for the same problem it returns "EXIT: Problem has too few degrees of freedom." Do you have any idea why it occurs?

3. Are there more solvers that can handle exponential problems?

Thanks a lot, Oron

[Joachim Dahl]

Perhaps the difference is error so small in this example, that the comparison is inconclusive.

[O S]

Dear Joachim, 

Indeed, it may be inconclusive. Still, my current goal is to improve the accuracy beyond the 1e-8. 

The reason for this is that I use the the convex optimization to derive an upper bound on a desired quantity (called channel capacity, in information theory), and then I check if the arguments of the solution satisfy a sufficient condition. If the condition is satisfied, the upper bound becomes an equality.

In several instances, I could show analytically that the optimal solution of the convex optimization satisfies the sufficient condition but in the numerical optimization, the condition is very sensitive to the arguments so I need better accuracy.

[Johan Löfberg]

without actual examples, this is a discussion in complete darkness

[O S]

I didn't want to overload with the code details, but many thanks for your support!

A simplifies example is below and attached as main.m; it needs the attached functions and 'variables. mat'. (All files are also zipped)

Please let me know if more explanations or other examples are needed.

clear all

load 'variables.mat'

%%The first part of the code is fmincon using YALMIP 

x = sdpvar(1,24);

%This function puts the variables x in a new matrix that corresponds to the non-zero entries of SQ. 

J = Joint(x,M);

%The following two functions return linear combinations of the variables x

YQ = MargYQ(squeeze(sum(J,3)),4,2);

Q_DUP = MargQ(YQ,2);

Constraints = [A*reshape(J,[1 256])' == b, x<ones(1,24), -x < zeros(1,24)];

assign(x,ones(1,24)); 

Obj = -kullbackleibler(reshape(YQ,[1,8]),Q_DUP) + sum(P.*reshape(J,[1 256]));

ops = sdpsettings('solver','fmincon','debug',1,'usex0',1);

ops.fmincon.TolCon = 1e-30;

ops.fmincon.TolFun = 1e-30;

ops.fmincon.TolFunValue = 1e-30;

ops.fmincon.TolFunValue = 1e-30;

optimize(Constraints,-Obj,ops)

x =  value(x)

% In this part, the optimal objective is computed

J = Joint(x,M);

YQ = MargYQ(squeeze(sum(J,3)),4,2);

Q_DUP = MargQ(YQ,2);

Obj = -kullbackleibler(reshape(YQ,[1,8]),Q_DUP) + sum(P.*reshape(J,[1 256]));

R = Obj/log(2);

%This part of the code is the cvx implementation of the same problem

cvx_begin quiet

    cvx_solver sedumi

    variable x(24) nonnegative

    J = Joint(x,M);

    YQ = MargYQ(squeeze(sum(J,3)),4,2);

    Q_DUP = MargQ(YQ,2);

    maximize -sum(rel_entr(reshape(YQ,[1,8]),Q_DUP)) + sum(P.*reshape(J,[1 256]))

    A*reshape(J,[1 256])' == b

cvx_end

RR = cvx_optval/log(2);

%The analytic solution for this problem (after division by log2) is log2(1+sqrt(5))-1

ERR_fmincon = abs(R-log2(1+sqrt(5))+1);

ERR_CVX = abs(RR-log2(1+sqrt(5))+1);

[O S]

I have now noted that I omitted from my code the cvx precision line: 

cvx_precision best 

If you add this line, then the sedumi outperforms with error 1.67e-11, while the fmincon is ~1e-8

[Johan Löfberg]

I would not draw any kind of conclusions from the numbers you report (as you only report the objective R, but not how the how the solver has satisfied the constraints on x). Also, you are comparing numbers beyond the intended precisions of the solvers so you basically look at primarily random effects

For instance, switch 'fmincon.alg' to 'sqp' and all of a sudden you have that you magic number is 10^-13 instead

BTW, you cannot use strict inequalities, and you could just as well write 0<=x<=1