As promised in the comments on my first answer, here is an implementation of an all-compiled-code Nelder-Mead minimizer, which hopefully represents a more useful response to the question. The algorithm used here corresponds to that given by Lagarias et al. in SIAM J. Optim. 9 (1), 112 (1998) (abridged .pdf). It is compatible with Mathematica versions 6, 7, 8, and 9, but not 5.2 or any previous version. This is not only due to the use of the new-in-6 functions OptionsPattern, FilterRules, and OptionValue, but also to apparent limitiations of the compiler--in particular, the robustness of the type inference mechanism was not entirely satisfactory prior to version 6.
The code is many times faster than NMinimize in all versions, although I would recommend using Mathematica 8 or 9 if possible. The performance of the compiled code is much better here than in versions 6 and 7, and many more functions are supported for compilation. Compilation to native code via C can also result in substantially improved performance. In fact, LibraryLink and/or the Mathematica runtime seem to have gained additional performance improvements in version 9, so this seems to be the optimal version to use as of this posting, being about 25% faster even than version 8.
A very important consideration is that, if the minimand is not compilable, performance will suffer due to calls out of compiled code to perform top-level evaluations. Indeed, these calls are so expensive that the resulting compiled function may easily be slower than the equivalent top-level code. It's also worth noting that FindMinimum possesses a very efficient implementation, so if only local optimization is needed, that function is likely to remain the best choice. For global optimization, an advantageous strategy might consist of using this package to quickly explore large parts of the parameter space (perhaps trying many different initial simplices) followed by the use of the optimized values as a starting point for FindMinimum, which will provide tight convergence to the final result.
Unlike for NMinimize, constrained optimization is not supported, because the Nelder-Mead algorithm is fundamentally an unconstrained method. For constrained problems, NMinimize performs a sort of regularization of the minimand such that the minimum of the resulting function fulfils the Karush-Kuhn-Tucker conditions, thus allowing unconstrained methods to continue be used. I may include this in a future update, but currently it is not implemented. Another difference relative to NMinimize is the convergence criterion: the one used here can more easily distinguish slow convergence from the minimizer having stalled without finding a minimum, which is useful for poorly behaved minimands. Instead of PrecisionGoal/AccuracyGoal, one specifies a tolerance, "ConvergenceTolerance", for the minimum allowed change in the average function value (sampled at the vertices of the simplex) within a given number of iterations. The default settings typically result in tighter convergence than NMinimize achieves, while still terminating the optimization if it genuinely does not converge.
This latest update contains several fixes and improvements:
- Option handling for
NelderMeadMinimize has been fixed--options given for this function were incorrectly being overridden by those of NelderMeadMinimize`Dump`CompiledNelderMead in previous versions, which would have been confusing to the user.
- It is now possible to refer to
NelderMeadMinimize`Private`dimension in options. This value represents the dimension of the problem and allows one to specify this parameter in an abstract way. An application of this will be demonstrated below.
- The interpretation of values given for the
"InitialPoints" option has been improved.
- Diagnostics in
NelderMeadMinimize`Dump`CompiledNelderMead can now be enabled for any return type. When disabled (as by default), the operation counts will no longer be maintained and no reference to them will appear in the compiled code. When enabled, these values will be given along with their descriptions in NelderMeadMinimize`Dump`OperationCounts on return.
- The code has undergone some general clean-up and should be easier to read as a result.
- An additional test function, the rotated (nonseparable) hyperellipsoid, has been provided. The rotation should not present much of a hindrance for the Nelder-Mead algorithm, but this is not necessarily the case for other approaches, particularly when scaled to hundreds of dimensions, whereupon e.g. differential evolution begins to encounter difficulties with it. This function is therefore useful for comparative purposes.
The package can no longer be presented in a code block in this post because it is too long, so please download it from its GitHub repository here.
Because the question involves performing a large number of similar minimizations, and in order to avoid expensive calls out of compiled code, I decided to inline the minimand into the Nelder-Mead algorithm itself. For each function minimized with a given set of options, compiled code is generated on the first call and memoized in order to amortize the compilation overhead over subsequent calls. The minimization can be run again with a different starting simplex (or some other set of initial points, specified via the "InitialPoints" option), or with different settings for "RandomSeed", "ConvergenceTolerance", or MaxIterations without re-compilation. Changing any other parameters or options will result in a new minimizer being generated.
To further reduce overheads, very little error checking is done. In fact, only the forms of some of the arguments are verified. As a result, if incorrect parameters or options are specified, the resulting errors will be returned to the top level.
For testing, I've included a few simple problems: the n-dimensional shifted hyperellipsoid function and its rotated counterpart (Schwefel's problem 1.2) and the n-dimensional generalized Rosenbrock's function. Contrary to common belief, the latter is not a unimodal function for all n: as shown by Yun-Wei Shang and Yu-Huang Qiu in Evol. Comp. 14 (1), 119-126 (2006) (link), there are actually two minima for n $\ge$ 4, and the Nelder-Mead algorithm (which is not strictly a global optimization algorithm) might converge to either of them. While these problems are not very difficult, I think they serve well enough for expository purposes. So, let's test the code. First, a simple usage example:
NelderMeadMinimize[x^2 + y^2, {x, y}]
(* or, equivalently, *)
NelderMeadMinimize[Function[{x, y}, x^2 + y^2], {x, y}]
(* or even: *)
With[{cf = Compile[{{x, _Real, 0}, {y, _Real, 0}}, x^2 + y^2]},
NelderMeadMinimize[cf, {x, y}]
]
(* -> {4.53016*10^-20, {x -> 1.90885*10^-10, y -> 9.41508*10^-11}} *)
(Note that when the minimand is passed as a pure or compiled function, the names of the variables are not actually important; they can be anything, as we demonstrate below. Note also that NelderMeadMinimize has HoldAll--although this can safely be removed if you prefer consistency with NMinimize to the convenience of not having to Block your variables.)
Now, a performance comparison:
(* Generate some variables *)
vars = Block[{x}, Unique[ConstantArray[x, 10], Temporary]];
(* First let's try NMinimize: *)
NMinimize[
NelderMeadMinimize`Dump`Hyperellipsoid @@ vars, vars,
Method -> {"NelderMead", "PostProcess" -> False}, MaxIterations -> 10000
] // Timing
(* -> {0.515625, {8.34607*10^-9, {
x$405 -> 0.999988, x$406 -> 1.000010, x$407 -> 1.,
x$408 -> 1.000030, x$409 -> 0.999995, x$410 -> 1.00001,
x$411 -> 0.999999, x$412 -> 1.000020, x$413 -> 1.00001,
x$414 -> 1.00001}}} *)
(* Now NelderMeadMinimize: *)
NelderMeadMinimize[
NelderMeadMinimize`Dump`Hyperellipsoid, Evaluate[vars],
CompilationTarget -> "C"
] // Timing
(* -> {0.391375, {1.73652*10^-16, {
x$405 -> 1., x$406 -> 1., x$407 -> 1., x$408 -> 1.,
x$409 -> 1., x$410 -> 1., x$411 -> 1., x$412 -> 1.,
x$413 -> 1., x$414 -> 1.}}} *)
We've achieved much better convergence, somewhat faster than NMinimize. But this includes the time taken for compilation to C! Trying again now that the minimizer has already been generated reveals that almost all of the above timing is in fact due to the compilation step:
Do[
NelderMeadMinimize[
NelderMeadMinimize`Dump`Hyperellipsoid, Evaluate[vars],
CompilationTarget -> "C"
], {100}
] // Timing
(* -> {1.296875, Null} *)
On the second and subsequent minimizations, we beat NMinimize by a factor of around 40, despite a tighter convergence tolerance. Let's now even the odds, as it's well known that the Nelder-Mead algorithm is quite slow to converge to very tight tolerances:
Do[
NelderMeadMinimize[
NelderMeadMinimize`Dump`Hyperellipsoid, Evaluate[vars],
"ConvergenceTolerance" -> 10^-9,
CompilationTarget -> "C"
], {100}
] // Timing
(* -> {0.953125, Null} *)
That's a better than 50-fold improvement over NMinimize in a more or less fair test, with each minimization of this 10-dimensional function taking under 10 ms. It may be of interest in some cases to record the number of function evaluations and the types of steps taken by the Nelder-Mead algorithm. If so, we may set the option "Diagnostics" -> True: after re-running the optimization we then find the relevant information recorded in the value of NelderMeadMinimize`Dump`OperationCounts:
NelderMeadMinimize`Dump`OperationCounts
(* {"Function evaluations" -> 2441,
"Reflections" -> 1289, "Expansions" -> 80,
"Contractions" -> 347, "Shrinkages" -> 0} *)
If absolutely minimum overhead is required, the compiled minimizer can be called directly, with the requirements that the minimand is given as a Function or CompiledFunction and the starting simplex is fully specified (taking the form of an array of real numbers having dimensions {d + 1, d}, where d is the dimension of the problem). Also, to specify MaxIterations -> Infinity, the third parameter should be a negative integer. This works as follows:
Do[
NelderMeadMinimize`Dump`CompiledNelderMead[
NelderMeadMinimize`Dump`Hyperellipsoid, vars,
CompilationTarget -> "C"
][RandomReal[{0, 1}, {Length[vars] + 1, Length[vars]}], 10^-9, -1],
{100}
] // Timing
(* -> {0.734375, Null} *)
This was a bit more work, but we have now achieved a 70-fold improvement over NMinimize. However, it should be noted that timings are generally much more sensitive to the initial simplex (and thus the number of iterations performed before convergence) than to the method in which the code is called. Working with the compiled minimizer directly is therefore perhaps better thought of as a means to incorporate it as a building block into other code (as shown below, where many minimizations are performed in parallel) than a means of achieving higher performance in its own right.
Now, we try to minimize the 50-dimensional Rosenbrock's function, even though the performance of the Nelder-Mead algorithm is usually worse (both slower and less reliable) than other methods (e.g. Storn-Price differential evolution) for high-dimensional minimization:
vars = Block[{x}, Unique[ConstantArray[x, 50], Temporary]];
NelderMeadMinimize[
NelderMeadMinimize`Dump`Rosenbrock, Evaluate[vars],
"RandomSeed" -> 10, CompilationTarget -> "C"
] // Timing
(* -> {24.109375, {2.44425*10^-15, {
x$567 -> 1., x$568 -> 1., x$569 -> 1., x$570 -> 1., x$571 -> 1.,
x$572 -> 1., x$573 -> 1., x$574 -> 1., x$575 -> 1., x$576 -> 1.,
x$577 -> 1., x$578 -> 1., x$579 -> 1., x$580 -> 1., x$581 -> 1.,
x$582 -> 1., x$583 -> 1., x$584 -> 1., x$585 -> 1., x$586 -> 1.,
x$587 -> 1., x$588 -> 1., x$589 -> 1., x$590 -> 1., x$591 -> 1.,
x$592 -> 1., x$593 -> 1., x$594 -> 1., x$595 -> 1., x$596 -> 1.,
x$597 -> 1., x$598 -> 1., x$599 -> 1., x$600 -> 1., x$601 -> 1.,
x$602 -> 1., x$603 -> 1., x$604 -> 1., x$605 -> 1., x$606 -> 1.,
x$607 -> 1., x$608 -> 1., x$609 -> 1., x$610 -> 1., x$611 -> 1.,
x$612 -> 1., x$613 -> 1., x$614 -> 1., x$615 -> 1., x$616 -> 1.}}} *)
We found the minimum in a reasonable time, although it required a non-default random seed to do so. As a performance comparison, a differential evolution minimizer I wrote in Python can minimize this function in about 21 seconds, which is certainly better considering that Python is interpreted while the result shown here is after compilation to C. However, this is still a huge improvement over NMinimize, which cannot detect convergence properly in this case, taking over 7 times as long trying (and ultimately failing) to find the minimum. In fact, we can do better if we employ the modified ("adaptive") scale parameters proposed by Fuchang Gao and Lixing Han in Comput. Optim. Appl. 51 (1), 259-277 (2012) (.pdf available from Gao's website):
With[{dim := NelderMeadMinimize`Private`dimension},
NelderMeadMinimize[
NelderMeadMinimize`Dump`Rosenbrock, Evaluate[vars],
"ReflectRatio" -> 1, "ExpandRatio" :> 1 + 2/dim,
"ContractRatio" :> 3/4 - 1/(2 dim), "ShrinkRatio" :> 1 - 1/dim,
CompilationTarget -> "C"
] // Timing
]
(* -> {16.781250, {1.5829*10^-15, { identical result omitted }}} *)
As described in the paper, these parameter values improve the efficacy of the expansion and contraction steps for high-dimensional problems and help to prevent the simplex from degenerating into a hyperplane, which otherwise would lead to failure of the Nelder-Mead algorithm. Convergence is thus achieved more reliably, to tighter tolerances, and without having to adjust the random seed. What is not so clear from this example is that, in favorable cases, the performance improvements can also be very dramatic: for the 35-dimensional hyperellipsoid function, for instance, the modified parameters yield a tenfold reduction in execution timing versus the classical values. I would thus strongly recommend at least trying the modified settings for larger problems, hence the incorporation of the symbol NelderMeadMinimize`Private`dimension to represent the problem dimension when doing so.
Edit: response to Ajasja's edits
Re-compilation of the minimizer on every call for a CompiledFunction objective function was a result of my inadequate testing of this case, so thanks go to Ajasja for noticing and reporting this issue, as well as the bug in handling specified initial points. Until it was pointed out to me by Leonid, I had somehow managed to overlook the fact that CompiledFunctions contain (in their second argument) patterns specifying the types of the arguments they accept. Pattern matching a CompiledFunction against itself will therefore not produce a match unless Verbatim is used to wrap the CompiledFunction being treated as the pattern, and memoization of function values similarly will not work where a CompiledFunction appears in an argument unless Verbatim is used in defining the applicable DownValue. This issue has now been fixed and compiled objective functions will be properly memoized by the updated version (both posted code and downloadable files) above.
The second point regards the question of how to incorporate parameters in the objective function other than the values actually being minimized. In fact this was possible without any additional modifications right from the first posted version of the code, although I didn't make this explicit or specify how it can be done. I hope to rectify this omission now, alongside describing how to obtain the best performance from this approach.
Let's take an example related to the scenario described in the question: namely, fitting a model to data. Here we will perform least-squares fitting of a cubic polynomial, i.e. the minimand is,
Norm[data - (a + b x + c x^2 + d x^3)]
with data (ordinate values only) and x (abscissae) given, and a, b, c, and d being the values under optimization. (In principle, using the monomials as basis functions is less than ideal because of its numerical instability, which combines poorly with the Nelder-Mead algorithm's tendency to get trapped in local minima. Practically speaking, it works well enough as an example.) This can be given to NelderMeadMinimize essentially directly:
fitter = Block[{data = #},
NelderMeadMinimize[
Block[{x = Range@Length[data]}, Norm[data - (a + b x + c x^2 + d x^3)]],
{a, b, c, d}
]
] &
The point to note here is that data appears lexically in the minimand, but not as a variable as far as NelderMeadMinimize is concerned. The first time this is called with actual data, compiled code will be generated that is a closure over the non-localized symbol data; where data is referenced, code is generated for a call into the main evaluator to retrieve its value. (The Block inside the minimand isn't relevant to this; it simply generates abscissae suitable for the given data and will be compiled completely since x is localized.) As it's the symbol data that appears inside the minimand and not actual data, compilation occurs only once rather than for every dataset fitted.
We try it:
datasets = Accumulate /@ RandomReal[{-1, 1}, {5, 100}];
fits = fitter /@ datasets;
fittedmodels = a + b x + c x^2 + d x^3 /. fits[[All, 2]];
Show[
{ Plot[fittedmodels, {x, 1, Last@Dimensions[datasets]}],
ListLinePlot[datasets] }, PlotRange -> All
]
Giving:
which seems like it was reasonably successful. So, this is a fine proof of principle, but quite slow ($\approx$ 100 ms/fit) due to the expensive call out of compiled code on every objective function evaluation. Obviously, we can do much better.
Since the aim is to completely eliminate calls into the main evaluator while fitting an arbitrary number of datasets, the new option "ReturnValues" of NelderMeadMinimize`Dump`CompiledNelderMead will come in useful. This enables compiled minimizers to be generated that produce any of several different return values, facilitating their use as building blocks in other compiled code. The possible option values are:
"OptimizedParameters": return a list of the values of the variables that minimize the objective function."AugmentedOptimizedParameters": as for "OptimizedParameters", but with the corresponding (minimal) value of the objective function prepended."Simplex": like "OptimizedParameters", but now returning a list of all d + 1 points of the final simplex obtained by the Nelder-Mead algorithm."AugmentedSimplex": as for "Simplex", but with each point having the corresponding value of the objective function prepended to it.
It seems to me that "ReturnValues" -> "OptimizedParameters" is the most suitable for the present application, so let's proceed as such. We now turn to the question of the parameter value accesses.
As Leonid has noted here, if compiled closures are inlined (using CompilationOptions -> {"InlineCompiledFunctions" -> True}) into other compiled code containing the values they close over, calls to the main evaluator can be eliminated entirely:
With[{
minimizer = NelderMeadMinimize`Dump`CompiledNelderMead[
Function[{a, b, c, d},
Block[{x = Range@Length[data]}, Norm[data - (a + b x + c x^2 + d x^3)]]
], {a, b, c, d}, "ReturnValues" -> "OptimizedParameters"
],
epsilon = $MachineEpsilon
},
serialFitter = Compile[{{datasets, _Real, 2}},
Table[minimizer[RandomReal[{0, 1}, {4 + 1, 4}], epsilon, -1], {data, datasets}],
CompilationOptions -> {"InlineCompiledFunctions" -> True},
RuntimeOptions -> {"Speed", "EvaluateSymbolically" -> False},
CompilationTarget -> "C"
];
parallelFitter = Compile[{{data, _Real, 1}},
minimizer[RandomReal[{0, 1}, {4 + 1, 4}], epsilon, -1],
CompilationOptions -> {"InlineCompiledFunctions" -> True},
RuntimeOptions -> {"Speed", "EvaluateSymbolically" -> False},
CompilationTarget -> "C",
Parallelization -> True, RuntimeAttributes -> {Listable}
];
];
Here the same minimand as used above is enclosed in a Function to make it suitable for NelderMeadMinimize`Dump`CompiledNelderMead, which is then called with this objective function and the option "ReturnValues" -> "OptimizedParameters" to generate a compiled minimizer that can be used from within other compiled code. Two callers are defined: serialFitter simply loops over each dataset given in its argument, while parallelFitter is Listable and automatically parallelizes over multiple datasets. Let's check them:
Block[{x = Range[50]},
Round@serialFitter[{3 + x - 4 x^2 + 2 x^3, -8 - 2 x + 7 x^2 - x^3}] ==
Round@parallelFitter[{3 + x - 4 x^2 + 2 x^3, -8 - 2 x + 7 x^2 - x^3}] ==
{{3, 1, -4, 2}, {-8, -2, 7, -1}}
]
As expected, we get True, so these can both correctly fit cubic polynomials. What about performance?
datasets = Accumulate /@ RandomReal[{-1, 1}, {1000, 100}];
serialFitter[datasets]; // Timing
(* 7.813 seconds *)
parallelFitter[datasets]; // AbsoluteTiming
(* 2.141 seconds *)
So, we have 7.8ms and 2.1ms per dataset, respectively. While these datasets, the model being fitted, and the convergence tolerances are admittedly all different to those in Ajasja's problem, that's still not too bad at all in my opinion. Furthermore, if you have a computer with support for simultaneous multithreading (SMT, e.g. Intel HT), the performance of parallelFitter can be further improved by evaluating SetSystemOptions["ParallelOptions" -> "ParallelThreadNumber" -> n] (where the value n depends on the number of logical processors available). Evidently, Mathematica's compiled code is not quite optimal, even after being translated to C and compiled to native code, since I found that this setting provided about 20% better performance on an Intel i7-2600 CPU.
AbsoluteTimingincludes everything going on on your computer until the calculation finished. What results do you get withTiminginstead ofAbsoluteTiming? $\endgroup$Do[Cos[0.3], {LOOPCOUNT}]? Is it just filler to represent something in your real problem? $\endgroup$Pause[]. So that one evaluation ofHi2ptakes about 40 microseconds. $\endgroup$