Suggestion for development of optimization solvers under MPT

I am afraid, that this suggestion is a bit out of scope for MCT. All above mentioned commercial solvers are realized typically at C/C++ and available Matlab toolboxes (all supports only double-precision class) plays only a role of interface between these low-level solvers (implemented as C/C++ MEX files) and Matlab environment.

The only way, in this case, is incorporation of multi-precision computing libraries (like MPIR  C Library for Multiple Precision Integers and Rationals, MPFR C Library for Multiple-Precision Floating-point computations with correct Rounding, MPC C Library for the Arithmetic of Complex numbers with arbitrary high precision, MPFR C++ multi-precision floating point number class, etc.), to these low-level solvers source code, which is definitely not simple task at all.

By my best knowledge of MCT, the only possible way, how to perform multi-precision optimization via MCT is modification of suitable optimization solvers, which are completely implemented at MATLAB to use mp class. But, I am sure (!!!), that for above mentioned commercial solvers are not available any strictly MATLAB coded solvers.

Dear Mark,

Thank you very much for your suggestion and for summary on optimization solvers. I will use it as a reference.

In order to support multiple precision, core C/C++ code of each solver must be modified. 

The minimal modification would be just changing the data types, constants, tolerances, etc. However in most of the cases these codes also rely on third-party libraries like Intel MKL/LAPACK, SuiteSparse or even some outdated FORTAN libraries (like NAG). This part must be modified to call analogous (if present) functions from MCT (also on C/C++ level). 

Overall this work requires tight collaboration with those vendors with access to C/C++ codes on both sides. 

We are open for such collaboration, but from my experience, other companies are not ready to invest because market is really tiny and it won't pay off.

For example, I talked to Maple top managers few years ago. They said that fast multiple precision is low priority task because of low demand.

P.S.

To my knowledge, developers of some of the libraries you mentioned (Mosek, CVX) are using MCT in their development work.

Well, I didn't think it would be easy, but perhaps much easier than writing a high quality optimization solver from scratch. There are a tremendous amount of non-obvious subtleties which make a huge difference in actual performance and reliability, and that is not easy to do.

SDPA is available in the variant SDPA-GMP (SDPA with arbitrary precision arithmetic);) http://sdpa.sourceforge.net/index.html . However, SDPA is not the most reliable LMI solver (not in Mosek's league).

quadMINOS is a quadrupke (not arbitrary) precision version of MINOS which is available commercially https://www.gams.com/latest/docs/S_MINOS.html .  See also https://web.stanford.edu/group/SOL/talks/14beijingQuadMINOS.pdf , including FAQ on the last page, which as of Dec 2016, states there is no MATLAB version of quadMINOS, and also that there is quadSNOPT (in f90 snopt9 we can change 1 line), i..e, quad precision version of the SNOPT continuous variable smooth nonlinear optimization solver.

As for optimization solvers written in MATLAB, there is the free solver PENLAB, which can handle nonlinear SDPs,(none of the solvers listed in my post above address nonlinear SDPs, to include Bilinear Matrix Inequalities (BMIs), and nonlinear constraints such as FMINCON addresses. It is a free, lower performing version of the commercial solver PENNON, now distributed via NAG in the NAG C library.. Unfortunately, PENLAB is riddled with bugs, and frequently hangs, unlike PENNON, or its other sibling, PENBMI (addresses BMIs rather than fully general nonlinear SDPs), bath written in C/C++, and having many fewer bugs than PENLAB. Supposedly a new version of PENLAB is under development, or is at least being considered for such, but I don't know its statues,, I wouldn't invest much effort in trying to make an arbitrary precision version of PENLAB until there's a relatively bug-free, or at least "usable" version available as a starting point. A version of PENNON packaged to be called from MATLAB will supposedly be available in NAG Toolbox Mark 27 in summer 2019 - however as of 2 years ago, it was going to come out 1.3.4 years ago, so who knows. Although PENLAB and PENNON address standard nonlinear continuous nonlinear programming (in addition to SDP constraints, they can't hold a candle to KNITRO when solving problems addressable by KNITRO (or FMINCON)., so they are not a good choice except for problems having SDP constraints, possibly in addition to other constraints,

Also I forgot to mention in my post above that KNITRO also has an SR1 Quasi-Newton option, in addition to BFGS and BFGS-L (LBFGS), whereas FMINCON does not have offer an SR1 option. The SR1 update, when used in a rust region, as in KNITRO, has better capability to handle sustained non-convexiities in the Lagrangian on the path to the optimum,on account of not enforcing (maintaining) positive semoidefiniteness of the approximation of the Hessian of the Lagrnagian as does BFGS, SR1 can quickly adjust its Hessian eigenvalues to match the curvature it encounters. But if your problem is convex, you're probably better off with BFGS.

Thank you for sharing your deep knowledge - very helpful.

I quickly checked the source code of SDPA-GMP.

Indeed, there are many weak points. It uses MPACK which is very buggy and unstable conversion of some LAPACK functions to arbitrary precision. I tested it myself few years ago - that was nightmare. E.g. half of the indices in the SVD code was hand-translated from FORTRAN to 0-based indices, another half still uses FORTRAN 1-based indices..... Not sure how people use it at all. I ended up writing my own MP-LAPACK with parallelism - it is part of toolbox now.

Another issue is that SDPA-GMP relies on Spooles for sparse matrices factorization. This is way too suboptimal, something like UMFPACK/CHOLMOD would be much better alternative. Also I would use something like iterative refinement to make sure our matrix solution has relative accuracy comparable to eps.

As for algorithm itself, I am not sure why code relies on Cholesky decomposition which is unstable in a presence of ill-conditioning. It already breaks when cond(A) ~ 1/sqrt(eps). Rank revealing (or at least with decent column pivoting) QR would be much more stable, it  can handle ill-conditioning up to cond(A) ~ 1/eps. If our goal is accuracy, then algorithm must seek the most reliable way to handle ill-conditioning. But I don't see this in SDPA-GMP.

QUADMINOS looks interesting but I am not sure if I can access the source code.

Thanks for sharing all the links - I will study them (wasn't aware of SR1 option in KNITRO).

What solvers would you suggest to have in arbitrary precision (by priority)?

In some sense of decreasing priority (but other people have different priorities I'm sure):

I'd pick KNITRO.  I suppose arbitrary precision need not be implemented for all algorithm options of which there are several https://www.artelys.com/tools/knitro_doc/2_userGuide/algorithms.html ,. combined with derivative options https://www.artelys.com/tools/knitro_doc/2_userGuide/derivatives.html .

Also, Mosek which can address LMIs and other conic problems (in forthcoming Mosek version 9, any product of Second Order Cones, semidefinite cones, exponential cones, and power cones, with exponential cones and power cones being new for Mosek 9). As you know with your eigenvalue work, sometimes you have an inherently ill-conditioned problem, and need to solve semidefinjite problems to high precision. I really don't like it when solvers cheat on the semidefinite constraints allowing slightly negative eigenvalue for psd constraint, and I don't like having to build in buffer against that by constraining to be >= small_number_as buffer_against_cheating). There actually is VSDP,http://www.ti3.tu-harburg.de/jansson/vsdp/VSDP2012Guide.pdf which is built on top of INTLAB, but I haven';t heard of anyone using VSDPt to practical effect.

CPLEX ,might be a good choice, but IBM is bureaucratic as hell, so good luck with that. Gurobi (smaller problem class than CPLEX, since it will not do non-convex QP or non-convex MIQP), also has top reputation for what it does, and was started by the founder of CPLEX after CPLEX was sold.

Arbitrary precision seems like it might not be useful for tough mixed integer problems (except when continuous relaxations are extremely ill-conditioned) because you are usually limited by long run-time to specify a non-trivial gap, such as 1% or higher,(as termination criterion)( between upper bound from best integer feasible solution and lower bound from continuous relaxation. But with arbitrary precision, you can reduce the integrality tolerance, allowing Big M logic constraints to have larger M (bound on variable) and have logic still function correctly. Many naive users specify a very large Big M value, which given default integrality tolerance means that their Big M constraints might not function as intended, thereby producing the wrong answer. But for certain very tough problems, the M needs to be very big, which means that the integrality tolerance needs to be very small, which can be achieved with arbitrary precision,.

Of course, there will be a huge performance penalty from doing floating point in SW.  I think some folks have looked at hybrid implementations which implement higher precision only in vital places, and use HW floating point for everything else.

Any progress on this topic?

No, there is no progress on this. High-quality solvers are commercial without access to source code (CPLEX, KNITRO, etc.). So this thread was just a discussion of wishful possibilities, not a realistic commitment/promise for development.

Besides, good solver enabled with extended precision capability would take years of efforts to design and implement. In fact, it would be comparable to toolbox itself by amount of work required. So we need some serious investor to come by and fund this work.

I fully understand your final respond. I expected something like this from the beginning of this discussion.

But, there is still open question: What should MATLAB users do in a case of need to perform any kind of optimization with multi-precision support?

The general recommendation is: use suitable optimization solver in pure MATLAB code and then modify it for MCT. But good optimization solvers are mainly not written in pure MATLAB.

So, situation is really bad for anybody who need to solve optimization problem with multi-precision accuracy.

In my case I am solving linear constrained optimization problem, where constraints are defined by exponential functions of optimized parameters. In MATLAB double class I am reaching always overflow problem.

I came across YALMIP project (MATLAB optimization environment), and there is REFINER: Built-in layer for arbitrary precision linear programming (https://yalmip.github.io/solver/refiner/), which is based on GEM library (https://github.com/jdbancal/gem), which is of course not so sophisticated as MCT.

Could be possible to create modified version of REFINER based on MCT???

Most likely it is possible. If YALMIP code is precision-independent then porting might be easy. Please try to port the functionality you need to MCT. 

Or maybe you can suggest authors of YALMIP to check/add compatibility with MCT. All they need is to write precision-independent code, e.g.: https://www.advanpix.com/2016/07/21/how-to-write-precision-independent-code-in-matlab/

This will make sure YALMIP can work with any toolbox.

The probably most simple and fully functional way how to solve LP and NLP optimization problem with full multi-precision support is based on Matlab toolbox for Maple (MTM) + Maple 2020 + Matlab 2020a. Everything works smoothly and surprisingly effective.

Communication with YILMAP development community was true nightmare, there is absolutely no motivation to change anything... :(   

Michal,

Can you pease provide more information on how that combination is used for multi-precision optimization, for example, constrained nonlinear optimization? Which optimization solvers can be used, and are all the under the hood parameter settings in the optimizer really changed for multi precision, or are some still based on settings the developer thought appropriate for a double precision solver?

Paval,

YALMIP is not a solver. It is an optimization modeling tool. It transforms problems, and calls any of a large number of supported solvers (all having MATLAB interfaces) to solve the problem. Similar to CVX, it can do conic reformulations and call conic solvers, but it can also handle non-convex problems. The developer is a Professor in Sweden, Johan Löfberg https://liu.se/en/employee/johlo46 , and he runs the YALMIP forum at https://groups.google.com/forum/?fromgroups#!forum/yalmip . YALMIP comes bundled with a few solvers, such as the BMIBNB branch and bound global optimizer, but other than that (and even for BMIBNB) requires other solvers to be supplied and called.In fact, YALMIP can call some higher precision solvers, although I guess thee is the matter of ensuring it does not muck anything up with double precision pre-processing.  Actually, presolving is one of the most "dangerous" operations performed bty double precision solvers - roundoff error can cause problems which are trivially feasible to be infeasible after presolve (large coefficients in the objective function can couple into the constrains during presolve).