# Simulated Annealing proof of convergence

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}$.

• Could you write out the formulation of the problem you're trying to solve? Is the objective function a black box (I presume it is), or is its mathematical expression available to you? What does the feasible set look like? This information may or may not be useful in answering your questions. – Geoff Oxberry Sep 26 '12 at 12:19
• @GeoffOxberry thanks Geoff! I have explicit objective function (for instance, the model in my example is the sum of three gaussians + one step), but function may or not be derivable. Explicit form of derivatives if any is unknown. Constraints are box constraints on parameters (l <= PARAMi <= u). does it help you to figure out what is appropriate, or should I give more details? – kiriloff Sep 26 '12 at 12:22
• If you could write all of that in the post, that'd be great. Stack Exchange supports Markdown with MathJax, which means you can use LaTeX-style notation to write out the math so that it looks like, well, math. – Geoff Oxberry Sep 26 '12 at 12:31
• I'm guessing that $p_4$, $p_7$, and $p_{10}$ are all bounded away from zero, because they all signify standard deviations? – Geoff Oxberry Sep 26 '12 at 12:53
• @GeoffOxberry bounds are well chosen, whatever the mathematical model is. here, indeed, parameters are more than 0. since data are available, 'visual' estimation of parameters' value may be possible, so that bounds can frame precisely the true solution. if you prefer, I want to reimplement an Excel like solver. – kiriloff Sep 26 '12 at 13:02

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.)

• Your worst case bound is misleading, because your $\epsilon$ is with respect to distance in the search space ($\delta$ might be more appropriate for this), not with respect to difference of function values between real optimum and approximate optimum. What I also find slightly misleading is that you refer to more modern optimization methods, without first addressing the question whether the OP uses SA correctly. SA certainly has its shortcomings, but it's strange that it doesn't converge at all. There might be some easy to fix bug somewhere. – Thomas Klimpel Sep 25 '12 at 21:24
• @ThomasKlimpelL I know of no implementation of SA that performed well in a comparison of more than a few methods. (Any algorithm can look nice if compared against other not very efficient methods.) I challenge you to test your favorite SA algorithm against one of the public codes that performed significantly better than the average algorithm on one of the two benchmarks cited. – Arnold Neumaier Sep 26 '12 at 7:58
• @ThomasKlimpel: The complexity for distance in $f$-space is not better. See the edited answer. – Arnold Neumaier Sep 26 '12 at 8:16
• Please, I only requested to use $\delta$ for the guaranteed error in $x$ and $\epsilon$ for the guaranteed error in $f$. The edited answer might contain the sketch of a proof, but I also find it misleading for the following reasons: 1) The global convergence guarantee of uniform random search (and simulated annealing) wouldn't hold, if such hand tailored evaluation point dependent target function were allowed. 2) The meaning of smooth (or continuous) is completely lost when each evaluation point is considered in isolation with the remark that we can interpolate between all given eval. points.. – Thomas Klimpel Sep 26 '12 at 11:13

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.

• @fonjibe: I agree with david. Simulated annealing is an approximation method, and is not guaranteed to converge to the optimal solution in general. It can avoid stagnation at some of the higher valued local minima, but in later iterations it can still get stuck at some lower valued local minimum that is still not optimal. – Paul Sep 25 '12 at 13:58
• @Paul thanks! 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 ? what is your advice? – kiriloff Sep 25 '12 at 14:18
• @Paul 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? thanks for any advice. – kiriloff Sep 25 '12 at 14:21
• @Paul and a last question: I mean that 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!!! – kiriloff Sep 25 '12 at 15:29
• @fonjibe: it looks like your question is much more complex than you wrote initially. Please incorporate all this information into your question. – Paul Sep 25 '12 at 15:38

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.

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.

• thanks Thomas, I do restrict vrariables to a convex region. since I am beginning with a MC sampling of the feasible variables space (which we assume is not too bad chosen), and SA parameters are adaptative (among other) with the starting point, I obtain very similar results with severals runs. However, my criterion 'quantity below tolerance threshold' is never met. results are closer to global optimum than with simplex downhill alone. However, solution found by Excel solver (GRG2 algo) is never met, and I am far from it. – kiriloff Sep 26 '12 at 6:11
• also, would you advise to go on with tuning (here each paper seems to have one own magic recipe, and experimental tuning is extremely difficult since tuning for one dataset may fail on next ...), or go to another method (GRG, CMA-ES)? – kiriloff Sep 26 '12 at 6:21
• It depends on what you want to achieve. If you want to learn more about SA and get some feeling for its real world behavior, I suggest to implement or download an implementation of a simple SA (not the downhill simplex variant), and convince yourself that it actually converges when you cool down slow enough. If you just need a fast optimization algorithm for your problem at hand, the better option is to try out some existing implementations of optimization algorithms. – Thomas Klimpel Sep 26 '12 at 11:25