# Optimization when Cost Function Slow to Evaluate

Gradient descent and many other methods are useful for finding local minima in cost functions. They can be efficient when the cost function can be evaluated quickly at each point, whether numerically or analytically.

I have what appears to me to be an unusual situation. Each evaluation of my cost function is expensive. I am attempting to find a set of parameters that minimize a 3D surface against ground truth surfaces. Whenever I change a parameter, I need to run the algorithm against the entire sample cohort to measure its effect. In order to calculate a gradient, I need to change all 15 parameters independently, meaning I have to regenerate all the surfaces and compare against the sample cohort way too many times per gradient, and definitely way too many times over the course of optimization.

I have developed a method to circumvent this problem and am currently evaluating it, but I am surprised that I have not found much in the literature regarding expensive cost function evaluations. This makes me wonder if I am making the problem harder than it is and that there might be a better way already available.

So my questions are basically this: Does anyone know of methods for optimizing cost functions, convex or not, when evaluation is slow? Or, am I doing something silly in the first place by rerunning the algorithm and comparing against the sample cohort so many times?

• Have you heard of stochastic gradient descent? For deep neural networks applied to large Training sets you have a similar issue (but can Eval gradient analytically) and the standard solution is to make gradient descent based only on a single sample(stochastic) vs the whole cohort (batch) Jan 31, 2016 at 10:06
• I'm only vaguely familiar with the area, so hence this being a comment rather than an answer. But what you are discussing sounds a lot like the topic of Uncertainty Quantification, often faced by engineers, where a single evaluation of the target function made take weeks to evaluate (at least in the problems faced by my engineering coworkers). My very limited understanding of how this is handled is by making a surrogate approximation that is much easier to evaluate based on past evaluations and simpler engineering models and then use these surrogate models to chose the next evaluation... Jan 31, 2016 at 22:45
• ... of the more expensive target function. I hate to say it, but I don't know anymore about the topic at the moment; I've only been told about it briefly while discussing research topics with said engineers. Interestingly, it seems like a very challenging research area: I believe good models require both a good understanding of the physics and statistics. Jan 31, 2016 at 22:47
• @seanv507 Yes, thanks, but I avoided it for a similar reason. It takes about 30 seconds to a minute to run one sample. If I have 15 parameters, I am looking at nearly 8 minutes per gradient calculation even if I only use one sample. If the space is large, it may take a very long time. Please correct me if you had other ideas in mind. Feb 1, 2016 at 1:52
• "what appears to me to be an unusual situation. Each evaluation of my cost function is expensive.". This is not at all an unusual situation, in general. It shows up all over the place, for example when ever your cost function comes from running a simulation (Eg in this paper: white.ucc.asn.au/publications/White2015PsoTransistorSizing.pdf we simulate a circuit in SPICE taking 10s). More boardly, in experimental science, evaluations can take ages. One of my friends Masters project is basically optimizing 5 parameters to find the best way of inserting DNA. Each evaluation takes 24 hours. Feb 1, 2016 at 8:55

# TL;DR

I recommend using LIPO. It is provably correct and provably better than pure random search (PRS). It is also extremely simple to implement, and has no hyperparameters. I have not conducted an analysis that compares LIPO to BO, but my expectation is that the simplicity and efficiency of LIPO imply that it will out-perform BO.

# LIPO and its Variants

This is an exciting arrival which, if it is not new, is certainly new to me. It proceeds by alternating between placing informed bounds on the function, and sampling from the best bound, and using quadratic approximations. I'm still working through all the details, but I think this is very promising. This is a nice blog write-up, and the paper is Cédric Malherbe and Nicolas Vayatis "Global optimization of Lipschitz functions."

LIPO is most useful when the number of hyper-parameters that you are searching over is small.

# Bayesian Optimization

Bayesian Optimization-type methods build Gaussian process surrogate models to explore the parameter space. The main idea is that parameter tuples that are closer together will have similar function values, so the assumption of a co-variance structure among points allows the algorithm to make educated guesses about what best parameter tuple is most worthwhile to try next. This strategy helps to reduce the number of function evaluations; in fact, the motivation of BO methods is to keep the number of function evaluations as low as possible while "using the whole buffalo" to make good guesses about what point to test next. There are different figures of merit (expected improvement, expected quantile improvement, probability of improvement...) which are used to compare points to visit next.

Contrast this to something like a grid search, which will never use any information from its previous function evaluations to inform where to go next.

Incidentally, this is also a powerful global optimization technique, and as such makes no assumptions about the convexity of the surface. Additionally, if the function is stochastic (say, evaluations have some inherent random noise), this can be directly accounted for in the GP model.

On the other hand, you'll have to fit at least one GP at every iteration (or several, picking the "best", or averaging over alternatives, or fully Bayesian methods). Then, the model is used to make (probably thousands) of predictions, usually in the form of multistart local optimization, with the observation that it's much cheaper to evaluate the GP prediction function than the function under optimization. But even with this computational overhead, it tends to be the case that even nonconvex functions can be optimized with a relatively small number of function calls.

A downside to GP is that the number of iterations to get a good result tends to grow with the number of hyper-parameters to search over.

A widely-cited paper on the topic is Jones et al (1998), "Efficient Global Optimization of Expensive Black-Box Functions." But there are many variations on this idea.

# Random Search

Even when the cost function is expensive to evaluate, random search can still be useful. Random search is dirt-simple to implement. The only choice for a researcher to make is setting the the probability $$p$$ that you want your results to lie in some quantile $$q$$; the rest proceeds automatically using results from basic probability.

Suppose your quantile is $$q = 0.95$$ and you want a $$p=0.95$$ probability that the model results are in top $$100\times (1-q)=5$$ percent of all hyperparameter tuples. The probability that all $$n$$ attempted tuples are not in that window is $$q^n = 0.95^n$$ (because they are chosen independently at random from the same distribution), so the probability that at least one tuple is in that region is $$1 - 0.95^n$$. Putting it all together, we have

$$1 - q^n \ge p \implies n \ge \frac{\log(1 - p)}{\log(q)}$$

which in our specific case yields $$n \ge 59$$.

This result is why most people recommend $$n=60$$ attempted tuples for random search. It's worth noting that $$n=60$$ is comparable to the number of experiments required to get good results with Gaussian Process-based methods when there are a moderate number of parameters.

Unlike Gaussian Processes, for random search, the number of queried tuples does not grow with the number of hyper-parameters to search over. Indee,d the dimension of the problem does not appear in the expression that recommends attempting $$n=60$$ random values. However, this does not mean that random search is "immune" to curse of dimensionality. Increasing the dimension of the hyperparameter search space can mean that the average result drawn from among the "best 5% of values" is still very poor. More information: The "Amazing Hidden Power" of Random Search? The intuition is that if we increase the volume of the search space, then we are naturally also increasing the volume of 5% of the search space.

Since you have a probabilistic characterization of how good the results are, this result can be a persuasive tool to convince your boss that running additional experiments will yield diminishing marginal returns.

• This seems like a modern variant of response surface methods! May 16, 2017 at 17:06
• Actually, random search can work remarkably well: argmin.net/2016/06/20/hypertuning
– Tim
Jun 21, 2018 at 15:09
• @Tim Yes, your point is well-taken. I didn't want to "decide" the issue of which is better in this post, since there are essentially endless permutations on BO, each claiming to be "the best" black-box optimizer -- making it impossible to be definitive. I agree that random search can work quite well, but I would actually recommend LIPO over PRS. LIPO is provably correct and strongly out-performs PRS (on average) in all of my experiments. LIPO also has minimal estimation cost: if you can minimize a QP, then you can use LIPO, and LIPO has zero hyperparameters (in contrast to BO).
– Sycorax
Jun 21, 2018 at 15:14
• Does the recommendation n=60 depend on the number of dimensions of space to be searched? I'm dealing with a hyperparameter space with 5 dimensions and 60 sounds small to me. Dec 10, 2022 at 0:34
• @DavidSilverberg It doesn't; dimension is not a number that appears in the equations that recommend $n=60$, but see caveats here: stats.stackexchange.com/questions/561164/… Essentially, you can make the search space so large that the "best 5% of hyperparameters" is still too poor quality on average to be useful.
– Sycorax
Dec 10, 2022 at 0:35

The literature on evaluation of expensive black-box function is quite vast and it is usually based on surrogate-model methods, as other people pointed out. Black-box here means that little is known about the underlying function, the only thing you can do is evaluate $f(x)$ at a chosen point $x$ (gradients are usually not available).

I would say that the current gold standard for evaluation of (very) costly black-box function is (global) Bayesian optimization (BO). Sycorax already described some features of BO, so I am just adding some information that might be useful.

As a starting point, you might want to read this overview paper 1. There is also a more recent one .

Bayesian optimization has been growing steadily as a field in the recent years, with a series of dedicated workshops (e.g., BayesOpt, and check out these videos from the Sheffield workshop on BO), since it has very practical applications in machine learning, such as for optimizing hyper-parameters of ML algorithms -- see e.g. this paper  and related toolbox, SpearMint. There are many other packages in various languages that implement various kinds of Bayesian optimization algorithms.

As I mentioned, the underlying requirement is that each function evaluation is very costly, so that the BO-related computations add a negligible overhead. To give a ballpark, BO can be definitely helpful if your function evaluates in a time of the order of minutes or more. You can also apply it for quicker computations (e.g. tens of seconds), but depending on which algorithm you use you may have to adopt various approximations. If your function evaluates in the time scale of seconds, I think you're hitting the boundaries of current research and perhaps other methods might become more useful. Also, I have to say, BO is rarely truly black-box and you often have to tweak the algorithms, sometimes a lot, to make it work at full potential with a specific real-world problem.

BO aside, for a review of general derivative-free optimization methods you can have a look at this review  and check for algorithms that have good properties of quick convergence. For example, Multi-level Coordinate Search (MCS) usually converges very quickly to a neighbourhood of a minimum (not always the global minimum, of course). MCS is thought for global optimization, but you can make it local by setting appropriate bound constraints.

Finally, you are interested in BO for target functions that are both costly and noisy, see my answer to this question.

References:

1 Brochu et al., "A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning" (2010).

 Shahriari et al., "Taking the Human Out of the Loop: A Review of Bayesian Optimization" (2015).

 Snoek et al., "Practical Bayesian Optimization of Machine Learning Algorithms", NIPS (2012).

 Rios and Sahinidis, "Derivative-free optimization: a review of algorithms and comparison of software implementations", Journal of Global Optimization (2013).

• @ lacerbi : I really like your answer! Can you please explain what you mean by costly and noisy? My guess is that "costly" means when you have a loss function with many model parameters (i.e. the weight connections in the neural network) and repeatedly evaluating the derivatives of such a loss function is very "costly". On the other hand, does "noisy" refer to loss functions that involve non-deterministic (i.e. random) variables? Why do non-gradient based optimization algorithms present advantages for "noisy" functions? can you please explain this further? thank you so much! Jan 17, 2022 at 6:16
• @stats555: Your guess is inexact in that you are describing a very specific case: that's not what "costly" generally means. First, here we are talking about black-box optimization, so there are no derivatives involved. Costly means costly: there is a cost involved in the evaluation of the objective function, be it time, energy, money, other resources, etc. and so you cannot afford to perform an arbitrary large number of function evaluations. Functions can be costly without the need to have many parameters (e.g., complex simulators). Jan 18, 2022 at 8:46
• @stats555: Noisy means stochastic, i.e. evaluating the target at the same input values twice may yield different output values for the target functions. That may happen when computation of the target involves randomness of some kind, broadly speaking. Jan 18, 2022 at 8:49

I don't know the algorithms myself, but I believe the kind of optimization algorithm that you are looking for is derivative-free optimization, which is used when the objective is costly or noisy.

For example, take a look at this paper (Björkman, M. & Holmström, K. "Global Optimization of Costly Nonconvex Functions Using Radial Basis Functions." Optimization and Engineering (2000) 1: 373. doi:10.1023/A:1011584207202) whose abstract seems to indicate this is exactly what you want:

The paper considers global optimization of costly objective functions, i.e. the problem of finding the global minimum when there are several local minima and each function value takes considerable CPU time to compute. Such problems often arise in industrial and financial applications, where a function value could be a result of a time-consuming computer simulation or optimization. Derivatives are most often hard to obtain, and the algorithms presented make no use of such information.

• Can you please explain what you mean by costly and noisy? My guess is that "costly" means when you have a loss function with many model parameters (i.e. the weight connections in the neural network) and repeatedly evaluating the derivatives of such a loss function is very "costly". On the other hand, does "noisy" refer to loss functions that involve non-deterministic (i.e. random) variables? Why do non-gradient based optimization algorithms present advantages for "noisy" functions? can you please explain this further? thank you so much! Jan 17, 2022 at 6:18
• @stats555: "Costly" just means "too expensive to compute". Like imagine optimizing the shape of an airplane wing, or the folding of a protein, or some other process that might take hours to simulate (...or possibly require you to run a physical experiment!) and/or potentially a ton of money. Simply taking one gradient is going to be incredibly expensive at that point, never mind N of them. Jan 17, 2022 at 9:57
• Thank you so much! Can you please describe "noisy" as well? Jan 17, 2022 at 10:02
• @stats555: For "noisy", I guess you can think of it as random variables, but I think it's more about the idea of not being able to compute the exact gradient (forcing you to approximate the gradient, which would add "noise"). Like think of stochastic gradient descent. You might not have enough resources to compute the full gradient, so you'd have to sample the dataset and compute a noisy approximation to it. Which might be fine—maybe that approach will work for your problem—but you can also imagine problems where a different approach that doesn't use gradients would be superior to SGD. Jan 17, 2022 at 10:04
• @stats555: As a trivial example, if your problem happens to be monotonic across the variables of interest, you might be able to just binary search (or interpolation search) instead of trying to approximate the gradient. Jan 17, 2022 at 10:06

You are not alone.

Expensive-to-evaluate systems are very common in engineering, such as finite element method (FEM) models and computational fluid dynamics (CFD) models. Optimization of these computatationaly expensive models is very needed and challenge because evoluationary algorithms often needs tens of thouands of evaluations of the problem which is not an option for expensive-to-evaluate problems. Fortunately, there are lots of methods (algorithms) avaiable to solve this problem. As far as I know, most of them are based on surrogate models (metamodels). Some are listed below.

• Efficient Global Optimization (EGO, also known as Bayesian optimization)  . The EGO algorithm has been mentioned above and may be the most famous surrogate-based optimization algorithm. It is based on the Kriging model and an infill criterion called expected improvement function (EI). R packages including the EGO algorithm are DiceOptim and DiceKriging.
• Mode-pursuing sampling (MPS) method . The MPS algorithm is built on the RBF model and an adptive sampling strategy is used to pick up candidate points. The MATLAB code is publised by the authors at http://www.sfu.ca/~gwa5/software.html. The MPS algorithm may need more evaluations to get the optimum, but can handle more complicated probelms than the EGO algorithm from my personal experience.
• Ensemble surrogates models by Juliane Müller . She used multiple surrogates to enhance the searching ability. The MATLAB toolbox MATSuMoTo is avaiable at https://github.com/Piiloblondie/MATSuMoTo.

In summery, these surrogate-based optimization algorithms try to find the global optimum of the problem using as few evaluations as possible. This is achieved by making the full use of the informations that the surrogate (surrogates) provides. Reviews on optimization of compuationally expensive problems are in [4-6].

Reference:

1. D. R. Jones, M. Schonlau, and W. J. Welch, "Efficient global optimization of expensive black-box functions," Journal of Global Optimization, vol. 13, pp. 455-492, 1998.
2. L. Wang, S. Shan, and G. G. Wang, "Mode-pursuing sampling method for global optimization on expensive black-box functions," Engineering Optimization, vol. 36, pp. 419-438, 2004.
3. J. Müller, "Surrogate Model Algorithms for Computationally Expensive Black-Box Global Optimization Problems," Tampere University of Technology, 2012.
4. G. G. Wang and S. Shan, "Review of metamodeling techniques in support of engineering design optimization," Journal of Mechanical Design, vol. 129, pp. 370-380, 2007.
5. A. I. Forrester and A. J. Keane, "Recent advances in surrogate-based optimization," Progress in Aerospace Sciences, vol. 45, pp. 50-79, 2009.
6. F. A. C. Viana, T. W. Simpson, V. Balabanov, and V. Toropov, "Metamodeling in Multidisciplinary Design Optimization: How Far Have We Really Come?," AIAA Journal, vol. 52, pp. 670-690, 2014/04/01 2014.
• Could you compare the methods you mentioned to the Bayesian optimization method mentioned in other answers?
– Hans
Nov 22, 2020 at 18:34
• @Hans Yes. The mentioned EGO and the Bayesian optimization are in fact the same algorithm. It is called EGO in the engineering design optimization field and is called Bayesian optimization in the statistics and machine learning fields. Other mentioned algorithms have very similar framework as the Bayesian optimization with different surrogate models and infill criteria. Dec 2, 2020 at 1:54
• Thank you. Would you care to also compare them to LIPO?
– Hans
Dec 2, 2020 at 14:45
• @Hans I am afraid not because I have no knowledge about LIPO. Dec 12, 2020 at 3:52

The two simple strategies that I have successfully used in the past are:

1. If possible, try to find a simpler surrogate function approximating your full cost function evaluation -- typical an analytical model replacing a simulation. Optimize this simpler function. Then validate and fine tune the resulting solution with your exact cost function.
2. If possible, try to find a way to evaluate an exact "delta-cost" function -- exact as opposed to be an approximation from using the gradient. That is, from an initial 15-dimensional point for which you have the full cost evaluated, find a way to derive how the cost would change by making a small change to one (or several) of the 15 components of your current point. You would need to exploit localization properties of a small perturbation if any in your particular case and you would likely need to define, cache, and update an internal state variable along the way.

Those strategies are very case specific, I don't know whether they can be applicable in your case or not, sorry if they are not. Both could be applicable (as it was in my use cases): apply the "delta-cost" strategy to a simpler analytical model -- performance may improve by several orders of magnitudes.

Another strategy would be to use a second order method that typically tends to reduce the number of iterations (but each iteration is more complex) -- e.g., Levenberg–Marquardt algorithm. But considering you don't seem to have a way to directly and efficiently evaluate the gradient, this is probably not a viable option in this case.

There are many tricks used in stochastic gradient descent that can be also applied to objective function evaluation. The overall idea is trying to approximate the objective function using a subset of data.

My answers in these two posts discuss why stochastic gradient descent works: the intuition behind it is to approximate the gradient using a subset of data.

How could stochastic gradient descent save time comparing to standard gradient descent?

How to run linear regression in a parallel/distributed way for big data setting?

The same trick applies to the objective function.

Let's still use linear regression as an example: suppose the objective function is $$\|Ax-b\|^2$$. If $$A$$ is huge, say a trillion rows, evaluating it once will take a very long time. We can always use a subset of $$A$$ and $$b$$ to approximate the objective function, which is the squared loss on the subset of the data.