Is expressed as \lVert \mathbf y - \mathbf A'\mathbf p(\mathbf x)\rVert_2^2, where \mathbf y is a simple vector, \mathbf A is a simple matrix (lower-triangular), and p_1,...,p_m are nasty polynomials.
Polynomials p_1,...,p_m are computed through a recursive implementation, say a function p(m,x) that computes all of them at once and which should be considered a black box (the number of coefficients of the polynomials are in the millions if we try to express them explicitelyā¦).
Moreover:
The loss is written in pure julia and automatic differentiation can pass through
But the very polynomial nature of the loss gives it a huge amount of local minimas, and any gradient descent (e.g. Optim.LBFGS()) gets stuck on the nearest local minimum.
For all these reasons, I am currently optimizing through Optim.ParticleSwarm(), which is a global algorithm for non-convex non-differentiable functions. It works well but itās quite slow.
Is there somewhere a global optimization algorithm that is more adapted to problems where:
The loss function is clearly non-convex with a lot of local minima.
But it is āsmoothā and automatic differentiation can pass through, and therefore local gradient descent are possibles.
It is required that the optimisation routines are implemented in pure julia as i use BigFloats. Furthermore if linear equality and bound contraints are posible itās a bonus
I tend to use MLSL for this sort of problem: it is a multi-start algorithm where you repeatedly use a local optimizer (which can be e.g. BFGS) from different starting points, and which includes some techniques to prevent it from repeatedly searching the same local optima. e.g. NLopt includes an MLSL implementation. It requires some experimentation to choose good termination tolerances for the local search.
You didnāt say how many parameters you have, however. In high dimensions global optimization becomes pretty hard by any method!
It is required that the optimisation routines are implemented in pure julia as i use BigFloats.
NLopt is calling an external C library and only handles Float64; youād have to re-implement MLSL or some similar algorithm.
You should really try to see whether you can avoid BigFloats. For example, there are lots of ways of working with polynomials that are numerically unstable (e.g. using a monomial basis, Lagrange interpolation, etcetera) and might at first seem to require BigFloat, but which can be reformulated in a stable fashion (Chebyshev-polynomial bases, barycentric interpolation, etcetera) that works fine in double precision.
Edit: I plan to add these to Nonconvex.jl, automatically extracting the functionsā expressions but when I get sometime to integrate Nonconvex with ModelingToolkit.
Can a Jump model use a standard julia function as an objective function ? In which case, it might be easy to write it as a Jump model and iāll try those two optimizers, for sure
No we already did a lot of investigations, and found out that multiple precision is indeed needed for this problem. The polynomials themselves are dense, in the sense that they have millions/billions of coefficients in the standard monomial basis (I cannot even compute them). The algo i have for the loss is currently stable with BigFloats, and thatās the reason iām on Julia. MultiFloats.jl also works well.
If there are no MLSL that are Julia-native then this is not a solution for me
The work will be in https://github.com/mohamed82008/Nonconvex.jl. There is no description other than that I want to run ModelingToolkit on the function, extract its expression and pass it to MathOptInterface for JuMP-based solvers to have access to it.
Looks perfect for what I want to do: allows me to specify my objective function as a genuine Julia function. Definitely tell me when you have a first working version
So many solvers and packages ! Everyday i found a new one at least.
So iām trying EAGO through JuMP, but it seems like BigFloats are not possibleā¦ I do not find how to specify that i want my variables to be BigFloats.
Iām sorry but i did not found an example without JuMP. As JuMP standards and syntax is quite new to me, I assumed that EAGo could not be used by itself. Do you have some kind of example ? I found nothing in the documentation (by a quick glance).
As I understood this relies on Nlopt.jl for local searches, so no Bigfloats ? Is there the possibility to use native Julia local searches (like Optim.LBFGS for example), to allow the whole thing to be type-agnostic ?