# Parallel optimization algorithms for a problem with very expensive objective function

I am optimizing a function of 10-20 variables. The bad news is that each function evaluation is expensive, approx 30 min of serial computation. The good news is that I have a cluster with a few dozen computational nodes at my disposal.

Thus the question: are there optimization algorithms available that would allow me to use all that computational power efficiently?

On one side of the spectrum would be an exhaustive search: subdivide the whole search space into a fine grid and compute the function at each gridpoint independently. This is certainly a very parallel computation, but the algorithm is horribly inefficient.

On the other side of the spectrum would be quasi-Newton algorithms: intelligently update the next estimate of the parameters based on the prior history. This is an efficient algorithm, but I don't know how to make it parallel: the concept of the "estimate of the parameters based on the prior history" sounds like a serial computation.

Quadratic algorithms seem to be somewhere in the middle: one can build the initial "surrogate model" by computing a bunch of values in parallel, but I don't know whether the remaining iterations are parallelizable.

Any suggestions on what kind of gradient-free optimization methods would perform well on a cluster? Also, are there any parallel implementations of optimization algorithms currently available?

• You can always calculate the gradient in parallel (for quasi-Newton methods using finite differences) and get a speedup proportional to the number of parameters i.e., 10x-20x. – stali Oct 17 '14 at 1:27
• @stali: You need the Hessian for quasi-Newton methods in optimization. Computing the Hessian via finite differences of function evaluations is really not a good idea. Computing finite difference approximations of the gradient for optimization is also generally not a good idea. – Geoff Oxberry Oct 17 '14 at 2:35
• Many quasi-Newton methods such as BFGS don't require the Hessian explicitly. I think by using gradients, in combination with L-BFGS the OP can quickly achieve what he wants. – stali Oct 17 '14 at 3:25
• @stali: I pointed out why using a finite difference approximation to the gradient would be a bad idea in my answer. It will degrade convergence by introducing error into the right-hand side of the quasi-Newton iteration. Also, it wastes function evaluations because there's no opportunity for reusing old evaluations (unlike surrogate methods). Using BFGS only addresses half the issues with your proposed approach. – Geoff Oxberry Oct 17 '14 at 9:37
• This is more appropriately a comment not an answer. But I have no choice, since I don't have enough rep to post a comment. Michael, I have very similar type of problem: expensive function evaluations that involve complex simulations running on a cluster. Did you ever found a code appropriate for running optimization when the function evaluation involves a simulation on a cluster? – MoonMan May 24 '16 at 19:37

With 10-20 variables and expensive function evaluations, the tendency is to recommend derivative-free optimization algorithms. I am going to disagree strongly with Paul's advice: you generally need a machine-precision gradient unless you're using some sort of special method (for instance, stochastic gradient descent in machine learning will exploit the form of the objective to come up with reasonable gradient estimates).

Each quasi-Newton step is going to be of the form:

\begin{align} \tilde{H}(x_{k})d_{k} = -\nabla{f}(x_{k}), \end{align}

where $\tilde{H}$ is some approximation of the Hessian matrix, $d_{k}$ is the search direction, $x_{k}$ is the value of the decision variables at the current iterate, $f$ is your objective function, and $\nabla{f}$ is the gradient of your objective, and the decision variables are updated like $x_{k+1} = x_{k} + \alpha_{k}d_{k}$, where $\alpha_{k}$ is a step size determined in some fashion (like a line search). You can get away with approximating the Hessian in certain ways and your iterations will converge, although if you use something like finite difference approximations of the Hessian via exact gradients, you might suffer from issues due to ill-conditioning. Typically, the Hessian is approximated using the gradient (for instance, BFGS type methods with rank-1 updates to the Hessian).

Approximating the Hessian and the gradient both via finite differences is a bad idea for a number of reasons:

• you're going to have error in the gradient, so the quasi-Newton method you're applying is akin to finding the root of a noisy function
• if function evaluations are expensive and you're trying to evaluate a gradient with respect to $N$ variables, it's going to cost you $N$ function evaluations per iteration
• if you have error in the gradient, you're going to have more error in your Hessian, which is a big problem in terms of the conditioning of the linear system
• ...and it's going to cost you $N^{2}$ function evaluations per iteration

So to get one bad iteration of quasi-Newton, you're doing something like up to 420 function evaluations at 30 minutes per evaluation, which means you're either going to be waiting a while for each iteration, or you're going to need a big cluster just for the function evaluations. The actual linear solves are going to be of 20 by 20 matrices (at most!), so there's no reason to parallelize those. If you can get gradient information by, for instance, solving an adjoint problem, then it might be more worthwhile, in which case, it might be worth looking at a book like Nocedal & Wright.

If you're going to do a lot of function evaluations in parallel, you should look instead at surrogate modeling approaches or at generating set search methods before considering quasi-Newton approaches. The classic review articles are the one by Rios and Sahinidis on derivative-free methods, which was published in 2012 and provides a really good, broad comparison; the benchmarking article by More and Wild from 2009; the 2009 textbook "Introduction to Derivative-Free Optimization" by Conn, Scheinberg, and Vicente; and the review article on generating set search methods by Kolda, Lewis, and Torczon from 2003.

As linked above, the DAKOTA software package will implement some of those methods, and so will NLOPT, which implements DIRECT, and a few of Powell's surrogate modeling methods. You might also take a look at MCS; it's written in MATLAB, but maybe you can port the MATLAB implementation to the language of your choice. DAKOTA's essentially a collection of scripts you can use to run your expensive simulation and collect data for optimization algorithms, and NLOPT has interfaces in a large number of languages, so choice of programming language shouldn't be a huge issue in using either software package; DAKOTA does take a while to learn, though, and has a massive amount of documentation to sift through.

• Its a pleasure for me to be completely wrong and learn something new and useful in the process:). – Paul Oct 17 '14 at 4:45
• Thanks! Just one more clarification: which of those algorithms are capable of performing function evaluations in parallel? For example, on k-way grid performing iterations n+1,...,n+k based only on information gained from iterations 1,...,n? – Michael Oct 17 '14 at 15:53
• @Michael: Many of the algorithms will have a step like "evaluate function at each of $k$ specified sample points", and that step can be embarrassingly parallel, usually with a scatter before that step to distribute sample points to processors, and a gather after that step to collect the results. Sometimes, evaluations can be cached (for instance, with surrogate modeling) and used for later iterations. – Geoff Oxberry Oct 18 '14 at 8:44

Perhaps surrogate-based optimization algorithms are what you are looking for. These algorithms use surrogate models to replace the real computationally expensive models during the optimization process and try to get a suitable solution using as few evaluations of the computationally expensive models as possiable.

I think the Mode Pursuing Sampling method may be used to solve your problem. This algorithm uses the RBF surrogate model to approximate the expensive objective function and can handle nonlinear constraints. More importantly, it selects multiple candidates to do the expensive function evaluations so that you can distribute these candidates for parallel computing to further speed up the search process. The code is open-source and written in MATLAB.

Reference

Wang, L., Shan, S., & Wang, G. G. (2004). Mode-pursuing sampling method for global optimization on expensive black-box functions. Engineering Optimization, 36(4), 419-438.

I am not sure a parallel algorithm is truly what you're looking for. It's your function evaluations which are very costly. What you want to do is parallelize the function itself, not necessarily the optimization algorithm.

If you can't do that, then there's a middle ground between exhaustive search and Newton algorithm, it's Monte Carlo methods. You can, on a bunch of different cores/nodes, start the same algorithm which is prone to fall to local optima (say quasi-Newton algorithms), but all with random initial conditions. Your best guess then for the true optima is the minimum of the minimums. This is trivial to parallelize and can be used to extend any method. While not perfectly efficient, if you have enough computing power at your disposal it can definite win the programming productivity vs algorithm performance battle (if you have lots of computing power, this can finish before you'd ever be finished making a fancier algorithm).

The choice of optimization algorithm (and thus its parallelization) highly depends on the properties of the objective function and constraints. Without knowing more about the problem, it's hard to give any kind of meaningful advice.

But from your considerations of newton methods, I infer that your objective function is differentiable. If possible, your problem would benefit greatly from parallelizing the function evaluation. If not possible, then you may also consider an inexact newton method which replaces the exact gradients/hessians by finite difference approximations. Then, you can use all those processors at your disposal to compute each non-zero element of the jacobian, as @stali suggests.

For more information, read Nocedal & Wright's Numerical Optimization, Chapter 7. There are many optimization software packages that implement this in parallel. Among the most widely used freeware is the DAKOTA software package (Sandia National Labs).

• Unless you have machine-precision gradients available (analytically, through adjoint computations, through some sort of forward sensitivity analysis), this approach is really not a good idea; it would require a huge number of simulations per Hessian evaluation, and you'd be better off taking those function evaluations and using them to build a surrogate model (for instance, like BOBYQA; ORBIT could build a surrogate model in $N$ function evaluations using radial basis functions). – Geoff Oxberry Oct 17 '14 at 2:38

Here is a solution to your problem.

Description of a mathematical method is provided in this paper.

• Welcome to SciComp.SE. Can you provide details about the approach described in the paper and implemented in the software? What is the method used? Why is it good? What is provided in this approach that the other answers don't cover? – nicoguaro May 18 '16 at 22:48
• Also, it seems that this is your own work. If that is true, please state that explicitly in your answer. – nicoguaro May 18 '16 at 22:49
• @nicoguaro: thank you, but I know how to click links. – Michael May 19 '16 at 2:33
• @Michael, it is not for you. The philosophy of this site is to be a collection of answers. You are obtaining your answer today, but in the future somebody else can need the same help. That is why there are standards de facto of what a good answer is. – nicoguaro May 19 '16 at 2:47
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 1, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 1, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.4548717439174652, "perplexity": 551.9851263660481}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2019-43/segments/1570987836368.96/warc/CC-MAIN-20191023225038-20191024012538-00344.warc.gz"}