# Finding a global minimum of a smooth, bounded, non-convex 2D function that is costly to evaluate

I have a bounded non-convex 2-D function which I'd like to find the minimum of. The function is quite smooth. Evaluating it is costly. An acceptable error is about 3% of the function's domain in each axis.

I tried running the implementation of the DIRECT algorithm in the NLOPT library, but it didn't give a considerable improvement over brute force search in terms of the amount of function evaluations needed for the required accuracy and there were some outliers.

Which other global optimization solvers should I consider?

-
Can you compute gradients, or would you need to approximate them by difference quotients? – Arnold Neumaier May 17 '12 at 12:25
I need to approximate them by difference quotients. – Victor May May 17 '12 at 12:41
In this case, Newton's method cannot be recommended, as numerical second derivatives are numerically very unstable, and difficult to tune to work safely. – Arnold Neumaier May 17 '12 at 12:45
@Victor May, what did you end up with ? (If you could post a function similar to yours, that would really help people to compare and tune different algorithms.) – denis Oct 16 '12 at 10:15
@Denis, I was trying to get more speed out of an algorithm for tracking an object in video. The output of the algorithm was a likelihood estimate for each image location to contain the tracked object. The image containing these likelihood estimates is the function I was trying to optimize. I ended up with brute-forcing at several resolution steps. For more information on the tracking algorithm in question read the paper "Robust Fragments-based Tracking using the Integral Histogram". – Victor May Dec 15 '13 at 19:10

I would like to suggest a somewhat different approach compared to the other answers, though @barron has indirectly discussed the same thing.

Instead of optimizing your function directly, i.e. by evaluating it at a series of points $\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_k$ points that (hopefully) converge to a (local) optimum, you could use the concept of $\textit{surrogate modelling}$, which is very well suited for problems of the type you describe (high cost, smooth, bounded, low dimensional i.e. less than 20 unknowns).

Specifically, surrogate modelling works by setting up a model function $c \in \mathbb{R}^d \rightarrow \mathbb{R}$ of your true function $f \in \mathbb{R}^d \rightarrow \mathbb{R}$. The key is that while $c$ of course does not perfectly represent $f$, it is far cheaper to evaluate.

So, a typical optimization process would be as follows:

1. Evaluate $f$ at a set of $j$ initial points $\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_j$. Note that derivatives are not needed. Also note that these points should be distributed evenly throughout the search space, e.g. by Latin Hypercube Sampling or a similar space-filling design.
2. Based on this original dataset, create a model function $c$. You could use cross validation to validate your model (i.e. use only a subset of the original $j$ points to create $c$, and then use the remainder of the dataset to check how well $c$ predicts those values)
3. Use a criterion such as the Expected Improvement (EI) criterion to find out where to ''fill in'' more samples to make $c$ more accurate by sampling $f$. This is actually far better studied theoretically than it might seem, and the EI criterion is very well researched. The EI criterion is also not a greedy criterion, so you both get good overall improvement of the model accuracy, whilst prioritizing accuracy near potential optima.
4. If your model is not accurate enough, repeat step 3, else use your favourite optimization routine to find the optimum of $c$, which will be very cheap to evaluate (so you could use any routine you want, even ones that requires derivatives, or simply evaluate the function in a fine mesh).

In general, this is what is meant by EGO, Efficient Global Optimization, as @barron suggested. I would like to stress that for your application, this seems perfectly suitable — you get a surprisingly accurate model based on relatively few evaluations of $f$, and can then use any optimization algorithm you want. What's often also interesting is that you can now evaluate $c$ on a mesh and plot it, thereby gaining insight into the general appearance of $f$. Another interesting point is that most surrogate modelling techniques also provide statistical error estimates, thereby allowing uncertainty estimation.

How to construct $c$ is of course an open question, but often Kriging or so-called space-mapping models are used.

Of course, this is all quite a bit of coding work, but a lot of other people have done very good implementations. In Matlab, I only know of the DACE software toolbox DACE is free. TOMLAB might also offer a Matlab package, but costs money — however, I believe it also works in C++ and has far more capabilities than DACE will ever have. (Note: I am one of the developers of the new version of DACE, soon to be released, which will offer added support for EGO.)

Hope that this rough overview has helped you, please ask questions if there are points that can be made more clear or stuff I've missed, or if you would like further material on the subject.

-
Fwiw, google surrogate-model brings up a Surrogate Modeling Lab at Ghent University, and a book Engineering Design via Surrogate Modelling , 2008 228p 0470770791. A problem with any very general approach is that you soon have a kitchen sink full of method variants, more than real test functions. – denis Oct 17 '12 at 17:05

See

L.M. Rios and N.V. Sahinidis, Derivative-free optimization: A review of algorithms and comparison of software implementations

for a very useful recent comparison of solvers.

-

For a smooth function, the Efficient Global Optimization method should perform quite well and be dramatically more efficient than DIRECT. Implementations are available in TOMLAB (haven't used it myself) and DAKOTA (which I've had some success with).

-

Since the function is smooth, Newton's method will be the overwhelmingly most efficient method to find a minimum. Since the function is not convex you will have to apply the usual tricks to make Newton's method converge (Levenberg-Marquardt modification, line search or trust region to globalize). If you can't get derivatives of your function try either computing it via finite differences or using a BFGS update. If you suspect that the problem has more than one local minimum, one would simply start Newton's method from a bunch of randomly or not quite so randomly chosen points and see where they converge.

-
My problem does indeed have local minima. What methods are there for choosing the starting points? – Victor May May 17 '12 at 10:26
Unless you know anything about the problem, statistical sampling is essentially your only choice. – Wolfgang Bangerth May 17 '12 at 13:28
@Wolfgang: Any ideas how to approach "statistical sampling"? Just try 10, 100, ... random initial guesses? Are there "more rigorous" approaches? I ask, because I have more or less a similar problem (see scicomp.stackexchange.com/q/4708/1789) – Andre Dec 30 '13 at 8:01
It all depends on what you know about the function. If you know something like a "typical length scale" for your function that would give an indication of how far local extrema would be separated. This will also give you an indication how many points you may have to start with, and how far apart they should be chosen of each other. – Wolfgang Bangerth Jan 8 '14 at 2:12