Simulated Annealing proof of convergence

Question

I implemented downhill simplex simulated annealing algorithm. Algorithm is very hard to tune, w.r.t. parameters including cooling schedule, starting temperature...

My first question is about convergence. In general, independently from choice of initial parameters, is it sure that algorithm converges to global optimal solution?

More precisely: I use SA with downhill simplex based on Num.recipes, however, I am adding some specific actions for observing constraints. Do you think my additional conditions for constraints are cause for bad convergence?

Then, I would like your advice on how to go on with solver conception: for instance, my algo is tested on data in 10-D space (10 parameters). It goes as follows: a preprocessing (monte carlo) picks up a good starting point. Then I run simplex downhill --> does not converge. I run simulated annealing (with simplex downhill to go to next step), and it remains far to solution (RMSE == 0.001). Excel solver on these data does the job very good. Excel method is generalized reduced gradient. Should i go on with fine tuning for SA, or go to non-derivative-free methods s.a. GRG ?

Finally, to speak the same language: with my example, convergence is not reach AT ALL. i don't have time to reach a local minimum either, since I have put some exit conditions when the simplex gets stuck. what does one mean in general with 'solver does not converge'? does it mean that termination condition (s.a. go below a tolerance value) is met and solution is sub optimal, or can it mean that termination condition is never met?

Thanks & regards.

EDIT

My objective function is given explicitly. It can, for instance, look like

$F(p) = \sum_{j=1}^{n}\Bigg(y_{j} - (p_{1} + p_{2}\exp{\Big(-\frac{(x_{i}-p_{3})^{2}}{p_{4}}\Big)} + p_{5}\exp{\Big(-\frac{(x_{j}-p_{6})^{2}}{p_{7}}\Big)} + p_{8}\exp{\Big(-\frac{(x_{j}-p_{9})^{2}}{p_{10}}\Big))}\Bigg)^{2}$

so, the squared residuals with a 'three Gaussians + step' model.

This is an example, solver should solve for any combination of functions among gaussians, exponential, power, polynomial, linear, rational with 'reasonable' number of parameters.

Box constraints on parameters are specified, with form $l_{i} \leq p_{i} \leq u_{i}$.

Answer 1

There is a theorem that syas that a black box algorithm is guaranteed to find the global minimum of an arbitrary smooth (i.e., twice continuously differentiable) function if and only if it samples points densely in the search space.

Here dense is meant in the topological sense, i.e., it must sample points in arbirarily small neighborhoods of every point.

In this case, the worst case complexity is $O(\Big(\frac{d}{\delta}\Big)^{-n})$, [the Latex parser chokes with nested big brackets] where $d$ is the diameter of the search space and $\delta$ the guaranteed error in $x$.

Edit: While this looks like being the worst case complexity for locating an $\delta$-accurate $x$ rather than for locating a point for a $\epsilon$-accurate $f$, the complexity for the latter is as bad (even for the rather big value of $\epsilon=(f_{first}-f_{global})/2$), as one can easily construct a smooth function which interpolates all data points so far but takes much smaller values in the center of the largest open ball not containing one of the points evaluated.

Thus guaranteed convergence to a global minimum is worthless in practice. (For example, uniformly random search has a convergence guarantee to the global minimizer, whereas most practically fast algorithms don't have one.)

Note that simulated annealing usually performs much worse than modern methods. Rather use code recommended in:

Comparison of derivative-free optimization algorithms (2012, by Nick Sahinidis)

Black-Box Optimization Benchmarking (BBOB) 2012 (by Auger, Hansen, et al.)

Answer 2

No, it is not. Very few methods are provably convergent for nonconvex optimization problems. If you're looking for such a method, you might look into branch-and-bound methods.

Answer 3

Do you have box constraints, or do your constraints at least define a convex region? If not, they certainly can cause bad convergence.

I looked up SA with downhill simplex in my copy of Num.recipes. This is not a normal SA implementation, and probably behaves more like an improved downhill simplex than like an "improved" SA implementation.

In my own experience, the convergence of a simple SA implementation can be surprisingly robust, even so you will probably need at least 100000 steps for 10 parameters. A simple SA implementation would randomly select an axis, and then do a "small" random move along that axis. My own attempts to do something more clever than this always slowed down convergence significantly or even completely ruined it. Perhaps something similar happens for SA with downhill simplex, even so it should still converge better than a pure downhill simplex.

To check for convergence, I found it helpful to always do three independent SA optimization runs, and see whether they produce similar results. I have to admit that most of the time when the three runs produced significantly different results, I actually had a bug somewhere in my implementation of the target function.

The (implicit) suggestion of Arnold Neumaier to try more than one optimization algorithm on your problem (for example from the DAKOTA project) is probably also a good idea. But these algorithms also have tuning parameters, so being familiar with the real world behavior of some optimization algorithms and their tuning parameters still seems useful.

Answer 4

You're probably better off switching to a deterministic method instead of simulated annealing. Provided you can import the data to a text file, you can learn the basics of the algebraic modeling language GAMS without too much effort, and use the BARON solver to attempt to solve your problem. BARON is a very good (though not foolproof) nonconvex nonlinear programming solver developed by Nick Sahinidis and Mohit Tawarmalani. Ten decision variables and box constraints should still be within the capabilities of a free GAMS license, but if that is not the case, you can submit a GAMS job (using the BARON solver) on the NEOS optimization server and obtain results that way.

Stochastic optimization solvers, while useful in some circumstances, can be very finicky when it comes to tuning parameters. You may even be better off attempting to solve your least-squares problem with a traditional convex or least-squares solver; even though it may get stuck in a local minimum, using judicious initial guesses or multistart could give you a sufficiently good solution for your purposes.