A Global Optimization Algorithm Worth Using
blog.dlib.net
blog.dlib.net
These solvers (BARON [1], ANTIGONE [2], Couenne [3] etc.) are capable of solving nonconvex global optimization problems that are structurally more complicated (i.e. non-ML applications) than the ones presented in the article linked above.
They employ a bunch of mathematically-rigorous methods such as interval analysis, branch-and-cut and McCormick relaxations. There's a treasure trove of little ideas from decades of global optimization research that I feel the Deep Learning community could benefit from taking a peek at.
[0] https://neos-server.org/neos/solvers/index.html
[1] http://archimedes.cheme.cmu.edu/?q=baron
I recently did an "Intro to OR" class and most of the tools we used were proprietary. I've been looking for more open source tools (and more info on OR in general, I think it's fascinating!).
Edit: NLOpt (https://nlopt.readthedocs.io/en/latest/) is also a nice suite of optimization algorithms for non-linear optimization.
The only issue is that much of the code on that site is very cutting-edge and academic. You may need to read a few benchmarks/publications to figure which solvers are production quality and which are not.
That said, a handful of these free solvers are so well-written that they exceed the performance of many commercial solvers for certain types of optimization problems.
(A notable exception is MIPs. Commercial MIP solvers like Gurobi or CPLEX are orders of magnitude better than the open-source CBC or GLPK. CBC is still pretty ok though, and I bundle it when I need to include a free solver in my code to solve smaller problems.)
For general optimisation, i've used NLopt. It was very easy to use, and has a pretty sensible set of algorithms.
One thing i've learnt is that comparing optimisation algorithms is really hard. One might be better at one problem, one at another. One pretty general axis of variation is problem size: some algorithms scale to orders of magnitude more variables (or constraints etc) than others.
In particular, i get the impression that a lot of the cutting-edge research is about being able to solve huge problems at all, or in hours instead of days. Perhaps those algorithms will be more useful for ML hyperparameter optimisation, but unfortunately, my needs are about solving much simpler problems in a fraction of a second. I get a lot more excited about Michael Powell's algorithms, which are cunning but simple, and work well on my problems.
It still boggles my mind how simplex, something that theoretically runs worse that interior point methods is competitive (at least based on what I've read online). I guess with these problems, the instances of the problem has a huge effect on how long approaches take.
>In particular, i get the impression that a lot of the cutting-edge research is about being able to solve huge problems at all
This is also the impression I get. When I look at benchmarks (btw, does anyone know an updated publication for these? A lot of what I find seems to be from years ago. I'm not sure how fast development of these are, but it would be nice to be updated once in a while) it always has some measure of time along with number of instances solved.
>Perhaps those algorithms will be more useful for ML hyperparameter optimisation
I thought about this, and maybe the reason why they largely don't use these have to do with getting results that are good enough (generalize well). A global optimum might not be worth the effort, so they stick to some variant of gradient descent. The measure they look at, after all, is performance on the test set. Aside from that, there may be something specific about a well defined problem that they can use to speed computations up, and a more general approach probably can't assume these for other instances of NLP for example.
Hans Mittelmann (also a wonderful individual) has you covered: http://plato.la.asu.edu/bench.html
Well, theoretical worst cases are just that. They are worst case bounds. People think NP-hard problem are intractable, but this is a misunderstanding. Many NP-hard problems can actually be solved with fairly good average case performance.
Case in point: the Simplex method is worst-case exponential time, but its average case performance is actually pretty good.
When Khachiyan first came up with his elipsoid method, it was polynomial time but in practice it was too slow. Karmarkar's interior point algorithm was polynomial time too, but performs much efficiently.
These days though, solve times are predominantly affected by solver heuristics and randomness. The choice of Simplex vs Interior Point does not make a significant difference in many cases.
Agreed. CBC/Clp are the best open-source LP/MIP solvers out there. SCIP is better, but it's not fully free.
> For general optimisation, i've used NLopt. It was very easy to use, and has a pretty sensible set of algorithms.
I noticed the omission of IPOPT support in NLopt. IPOPT is one of the best large-scale nonconvex NLP solvers out there (if you have exactly 2nd order derivative information, typically from AD).
> One thing i've learnt is that comparing optimisation algorithms is really hard. One might be better at one problem, one at another. One pretty general axis of variation is problem size: some algorithms scale to orders of magnitude more variables (or constraints etc) than others.
Very true. The academic world has standardized on Dolan-More charts for comparing optimizers on a standard corpus of problems, and those are the closest we have to quantifying general performance. But the no-free lunch theorem applies. There's also a lot of randomness and chance involved, e.g. initialization points can drastically affect nonconvex solver performance on the same problem. So can little things like ordering of constraints or the addition of redundant non-binding constraints (!) for MIPs (this was demonstrated by Fischetti). That's why modern MIP solvers include a random perturbation parameter.
> ML hyperparameter optimisation
Most types of hyperparameter optimization are black-box optimization problem with no analytic derivative information, which much of research in deterministic optimization does not address. Black-box problems typically rely on DFO (Derivative-Free Optimization) methods. DFO methods are even harder to quantify in terms of performance. The Nelder-Mead Simplex method is as good a method as any.
I've had good results from simplex; i'm a bit sad to hear there isn't some magical algorithm that will give me a massive speedup over it.
I haven't tried automatic differentiation. That seems extremely daunting given the volume of code in my problem. Someone has done some groundwork on AD of the main library i use in my model, but it would still be a serious bit of work.
Optimal knot placement, though seemingly simple, can be a gnarly problem due to degeneracy in many cases. Anecdotally I've found quadratic splines to still be somewhat tractable with exact methods, but polynomial splines with cubic or higher orders can be tricky, especially when they are large-scale.
For derivative-free optimization, you can sometimes get speedups from just adopting a embarassingly parallel algorithm like Particle Swarm. These types of stochastic algorithms aren't necessarily algorithmically better than Nelder-Mead Simplex, but they can explore a larger terrain in parallel, which may help. If you can find some parallel stochastic code, you should be able to spin up a VM on the cloud with lots of cores to help you search the space.
I've been looking at open source solvers for a while now to solve an MIP, and this one seems to be the best.
There is this annoying thing about the whole ecosystem though, aside from the proprietary bits.
I started by looking at GLPK and wrote my MIP in GMPL. For some reason a lot of tools I looked at don't "just work". I have to export it to some other format. It works, but feels like a work-around.
I wish there were a lingua franca in this domain and sort of have support from all major solver implementations, but each one seems to have their own way, or just provide an API (ehem google/or-tools)
I am happy I can make it work, but it may not be as easy for others.
JuMP is kind of that in terms of modeling software - http://www.juliaopt.org/
Pyomo is also pretty useful if you'd rather use Python than Julia, though it's not as well-documented and doesn't make installation of the open-source solvers quite as easy for you.
If you want to get a feel for it, here's an example in GMPL, which is a subset of AMPL used in GLPK:
What I was thinking of was a common language that can be used by all solvers. Similar to how C can be compiled by gcc or clang, this language can be used directly with the solver binary. Something like `export CC=my_solver` and `CC -d data.dat -m model.mod`, if you will.
What a lot of these projects do feels like calling a solver by translating it into some format that's different for everyone, whereas I wanted for the solvers to all agree on some language.
Or at the very least, some standard similar to how languages have libraries for something like XML. It would be nice if JuMP and PyOmO could read these standardized formats.
Then again, it's probably a pain to implement and even C compilers implement different extensions. But one can dream.
It's the common language that matter more though. I've used GAMS in my undergrad, so I'm under the impression that it's had at least some adoption, but it's proprietary AFAIK. I've know some people who use AMPL too, which also seems popular, but I'm not sure about it's licenses.
Umm, they (mostly) do. Pyomo supports the AMPL .nl format, which means you can call any AMPL solver from Pyomo (most COIN-OR solvers have AMPL interfaces). The AMPL .nl format is the de facto standard for open-source optimization solvers because it's the only open standard, I believe [0].
So you would define your mathematical program in Pyomo (which isn't a terribly challenging syntax compared to AMPL). When you invoke solve, Pyomo will spit out an .nl file and call a solver (which reads the .nl file, solves the problem, and outputs a .sol file which can be read back into Pyomo).
Incidentally this is exactly how it works in AMPL as well -- all the interop is done through .nl and .sol files. I recommend reading [0] to understand how this works.
As for Pyomo, while at first glance I dislike the syntax (it's not as easy to share with others, at least compared to something like GNU MathProg, where the receiver doesn't need to concern themselves with python or julia), I like python enough to maybe give it a try. I'll have to learn more about it though. Thanks to you and the others!
The proprietary part of AMPL is the program that turns human-readable .mod and .dat files into low-level .nl encodings of the problem representation. The .nl format is documented, and there's an MIT licensed library on netlib (and now github too) that reads it and evaluates functions, constraints, gradients, Jacobians, Hessians, etc. (AD is super useful and has been around way longer than the deep learning popularization of it.)
Julia has a low-level solver-independent API called MathProgBase that sits between solver wrappers and modeling languages like JuMP / Convex.jl. It's in the process of being rewritten as MathOptInterface, but it would allow people to write solvers in pure Julia and be usable right away with JuMP.
The MiniZinc Challenge (http://www.minizinc.org/challenge2017/results2017.html) is a competition between various solvers where the organizers enter some MIP solver back-ends also.
Personally I enjoy using MiniZinc as a language for prototyping models in, and enabling me to experiment and compare various different solver technologies easily.
The (sole, afaik) library author and maintainer Davis King is also extremely responsive and helpful. dlib and Davis deserve a lot more recognition.
[0] http://dlib.net
However I just wish it's Python binding were published as binary wheels. Dlib is very large and the compilation times suck when installing from pip.
In fact, I might as well contact him and help setup this!
Another method for global optimization is interval arithmetic, but it is not a black-box method.
By the way, the parent mentioned interval arithmetic. If anyone is looking for optimization algorithm look for "interval Newton". It's a cool algorithm.
Instead, the algorithm promises to actually match that lower bound, while always being at least as good as random search, and offering exponentially large improvements in certain cases. Not exactly the kind of claim that has implications for P =? NP.
Why? Continuous global optimization on convex problems can be solved in polynomial time. Did you mean on nonconvex problems?
Also, x*(x-1)=0 is a complementarity constraint and it is mathematically degenerate (it violates constraint qualifications) when solved as a continuous optimization problem. There are many ways to work-around it [0, 1, 2] and even solvers built for this of problems (MPECs = Mathematical Programs with Equilibrium Constraints) which I suspect only work well for narrow classes of problems like Linear Complementarity Problem (LCPs).
I have found that ultimately for nonlinear non-convex problems, the solution technique for continuous optimization of complementarity constraint problems has been empirically unreliable or the solution is often low equality (even if it can be shown that the complementarity penalty method does not violate constraint qualifications).
Bilinear constraints generally add nonconvexities and add a lot of difficult terrain for a continuous optimizer.
For general nonconvex problems, it's much more reliable to cast the problem as an NP-hard MINLP.
And if you're solving an LCP, you might as well recast it as an MIP rather than mucking around with an LCP. MIP solvers like Gurobi and CPLEX are so amazingly performant these days that it's just a no-brainer
[0] http://www.tandfonline.com/doi/abs/10.1080/10556780410001654...
[1] http://epubs.siam.org/doi/abs/10.1137/S1052623403429081
[2] http://epubs.siam.org/doi/abs/10.1137/040621065
[4] https://www.artelys.com/tools/knitro_doc/2_userGuide/complem...
https://mathoverflow.net/questions/92939/is-that-true-all-th...
About Lipschitz interval methods, as stated here (first googled example, this from 2001): https://link.springer.com/referenceworkentry/10.1007%2F0-306...
"Interval algorithms for constrained and unconstrained optimization are based on adaptive, exhaustive search of the domain. Their overall structure is virtually identical to Lipschitz optimization as in [4], since interval evaluations of an objective function ϕ over an interval vector x correspond to estimation of the range of ϕ over x with Lipschitz constants. However, there are additional opportunities for acceleration of the process with interval algorithms, and use of outwardly rounded interval arithmetic gives the computations the rigor of a mathematical proof."
edit: Any grad students looking for a paper idea might consider extending this result to allow stochastic objective functions satisfying a stochastic Lipschitz condition.
> However, if you want to use LIPO in practice there are some issues that need to be addressed. The rest of this blog post discusses these issues and how the dlib implementation addresses them. First, if f(x) is noisy or discontinuous even a little it's not going to work reliably since k will be infinity.
> ... Now each sample from f(x) has its own noise term, σi, which should be 0 most of the time unless xi is really close to a discontinuity or there is some stochasticity.
Is that not what you were pointing out?
[1] https://link.springer.com/article/10.1007/s10898-016-0482-9 / http://ssamot.me/papers/global-optimisation.pdf
I'd like to share another global optimization algorithm 'worth using' that I've recently discovered for myself: Hyperband [0]. It's extremely simple to implement (i.e. a dozen lines of code), is based on successive halving principle and is particularly suitable for iterative machine learning algorithms where early stopping is applicable (i.e. NN, GBM).
[0] https://people.eecs.berkeley.edu/~kjamieson/hyperband.html
this formula defines a function U(x). This means that it gives a recipe that takes a number x and produces a number U(x). The recipe depends on an already prepared set of ingredients: a set of numbers x1, x2, ..., xt, k, and another function f. Now, given a number x, to compute U(x), you first compute all the numbers f(xi)-k*|x-xi|, and then you pick the smallest one. This smallest value is then U(x).
The graph below the formula displays its interpretation. The graph of the function f is the red curve, the points (xi, f(xi)) are the black dots, and U(x) is in green.
f(x) - f(y) <= k|x-y|
hence
f(x) <= f(y) + k|x-y|.
Therefore the function:
u(x) = f(y) + k|x-y|
is an upper bound for the function f(x). Repeating this for multiple points x1, x2, x3,... you can construct a stricter upper bound by taking the minimum:
U(x) = min_i f(xi) + k|x - xi|
The equation says that U(X) is equal to the minimum value of the expression f(x_i) + k * (Euclidean distance from x to x_i), where i ranges over integers 1 to t.
In Go pseudo-code:
func U(X []float64, x [][]float64, f func([]float64) float64, k float64) float64 {
min := math.Inf(1)
for _, x_i := range x {
e := f(x_i) + k * EuclideanDistance(X, x_i)
if e < min {
min := e
}
}
return min
}It seems like even in his demo video, k starts off underestimated which made the upperbound not a upperbound (f is sometimes above the green line U) and is only slowly corrected towards the end. Why does it still work?