Khac Vu
email: vu@lix.polytechnique.fr

Claudia D ' Ky

Youssef Ambrosio
email: dambrosio@lix.polytechnique.fr

Leo Hamadi

Liberti
email: liberti@lix.polytechnique.fr

Surrogate-based methods for black-box optimization

In this paper, we survey methods that are currently used in black-box optimization, i.e. the kind of problems whose objective functions are very expensive to evaluate and no analytical or derivative information are available. We concentrate on a particular family of methods, in which surrogate (or meta) models are iteratively constructed and used to search for global solutions.

Introduction

In this paper, we survey about a special method for solving the following optimization problem:

z = min {f (x) | x ∈ D}, (1)
in which f : R d → R is a continuous black-box and D is a compact subset of R d . Here, black-boxes refer to functions that have no analytical form and are often expensive to evaluate. Black-box functions appear in a variety of different settings. The most notable example is computer simulation programs, which may consist of thousands of code lines and take hours, or even days to complete a single run. Since the computer programs implementing black-box functions are often very complicated, it is difficult, impractical, or downright impossible, to obtain an explicit description of these functions. Thus, simulation programs can be treated as black-box functions.

The methods we will present are based on the construction of surrogate models to approximate black-box functions and using them to search for (global) optimal solutions. They are often called surrogate-based method. In this paper, we will mainly focus on continuous variables. The continuity and compactness of f are technical assumptions which ensure the existence of a global minimum of this problem. Later on we will briefly discuss about the case when there are discrete variables.

Motivational examples

One classic example of real world application of black-box optimization is reservoir engineering, e.g., the so-called wells placement problem. In this problem, we aim at deciding the location and geometry (i.e. trajectory) of wells. The independent variables are: the number of wells, their trajectory, and their type (injector or producer). The main constraints are represented by minimum distance between wells to be respected while the objective is typically to either maximize the oil production or the Net Present Value (NPV) function. For a given solution, the oil production or the NPV is given as result of a simulation that is heavy in practice (minutes for simplified cases to several hours for the more realistic ones), see PumaFlow [START_REF] Franlab | Pumaflow reservoir simulator reference manual[END_REF]. In the literature, evolutionary/metaheuristic algorithms can be found, see, for example, Bouzarkouna et al. [START_REF] Bouzarkouna | Using evolution strategy with meta-models for well placement optimization[END_REF] and Onwunalu et al. [START_REF] Onwunalu | Application of a p article swarm optimization algorithm for determining optimum well location and type[END_REF]. More recently, hybrid approaches have been proposed, see, for example, Lizon et al. [START_REF] Lizon | A mixed-integer nonlinear optimization approach for well placement and geometry[END_REF].

INTRODUCTION

2

Building design is another field to which black-box optimization is often applied. One of the examples is the SmartBuilding project of Microsoft Research [START_REF] Fialho | A multi-objective approach to balance buildings construction cost and energy efficiency[END_REF], with the aim to devise new methods for optimizing smart buildings design by minimizing energy consumption, construction cost, and so on. The energy consumption of a building (in one year period) is estimated by EnergyPlus, a simulation program that is computationally expensive, which takes up to 1 hour for a single simulation. Decision variables are selected to be both continuous and discrete: the orientation angle, the scaling factors of window sizes, and the choices of material for each inner insulation layer of the external walls. In [START_REF] Fialho | A multi-objective approach to balance buildings construction cost and energy efficiency[END_REF], the problem is solved using an evolutionary algorithm called HYPE. A similar problem in architectural designs [START_REF] Costa | RBFOPT: an open-source library for black-box optimization with costly function evaluations[END_REF] is solved using a Radial Basis Functions based method.

Heuristical methods

There are many difficulties associated with black-box optimization. The most obvious is that the problems has no structure, so we can not exploit their analytical and derivative information. Due to this lack of structure, the most natural way to attack black-box optimization is to used heuristics such as evolutionary algorithms, simulated annealing, tabu search... This approach has the advantage that it only relies on function evaluations and is relatively easy to implement. However, due to the computational cost of function evaluation, these methods are only applicable for problem with small dimension since they often require a large number of function evaluations. Moreover, sometimes they are not desirable since they lack of mathematical convergence guarantees: generally, we do not know how close our solutions are from the optimum.

Derivative-free methods

The lack of derivatives motivates us to use derivative-free optimization (DFO) methods to solve blackbox optimization problems. One of these methods is directional direct-search, which is often studied in a search-and-poll framework [START_REF] Audet | A survey on direct search methods for blackbox optimization and their applications[END_REF][START_REF] Audet | Mesh adaptive direct search algorithms for constrained optimization[END_REF][START_REF] Torczon | On the convergence of pattern search algorithms[END_REF]. This method makes use of the fact that if a solution is not locally optimal, then we can improve it by moving along one of the directions in some positive bases (or spanning sets). Based on this observation, from a given point x k , we will look at all vectors u k in a positive basis and select an appropriate step size α k such that f (x k + α k u k ) < f (x k ). Since the number of vectors in a positive basis is large and feasible step sizes are unknown, in general, it might be expensive to find an updated better point. Moreover, most of the directional direct-search algorithms can only find local solutions. Mesh adaptive direct search (MADS), proposed by Audet and Dennis [START_REF] Audet | Mesh adaptive direct search algorithms for constrained optimization[END_REF], is one of the few algorithms in this class that can converge to local solutions in smooth and non-smooth cases. This algorithm is implemented in NOMAD [START_REF] Abramson | The NOMAD project[END_REF][START_REF] Digabel | Algorithm 909: NOMAD: Nonlinear optimization with the MADS algorithm[END_REF], an open-source software for solving BB optimization, which can also deal with black-box constraints, multiple objectives and parallelism, in both continuous and discrete cases.

Another DFO method is simplicial direct-search. In this method, new points are searched in the direction that is as far as possible from the vertex with the worst function value. One of the most popular methods in this class is proposed by Nelder-Mead [START_REF] Nelder | A simplex method for function minimization[END_REF], in which a simplex consists of d + 1 points are maintained at each iteration. Then, new points are found by performing either reflection, expansion, or contraction. The worst points in the simplex will be replaced by a point constructed from one of the 3 points. This method converses, but sometimes not to a stationary point. Several additional modifications are taken into consideration to overcome this drawback.

The two above methods may require substantial computing resources since the number of function evaluations is sometimes very large. This occurs because the inherent smoothness of the blackbox function is not exploited appropriately. A method which overcomes this disadvantage to a certain extent is the trust-region algorithms. In trust-region methods, at each step, the objective function is approximated by a linear (or quadratic) model using current data points. The model is solved within a local region (also called trust-region), which is basically a ball around the best current point (for certain norms). The

INTRODUCTION

3 trust-region is then adaptively changed, either by increasing or decreasing the radius or by moving its center, based on how much improvement of the newly found solution compared to the current best point. With the assumption that the model is fully linear (or fully quadratic), trust-region method converges to stationary solutions [START_REF] Conn | Introduction to derivative-free optimization[END_REF]. However, similar to the two previous methods, it fails to converge to a global solution.

Surrogate-based methods

In this paper, we will concentrate on another approach which appears to be very successful at finding global solutions of black-box optimization problems. The main idea behind these methods is to iteratively construct surrogate models to approximate the black-box functions (globally) and use them to search for optimal solutions [START_REF] Sóbester | Engineering design applications of surrogate-assisted optimization techniques[END_REF]. A common approach (in its simplest form) for surrogate-based methods is as follows

• Phase 1 (design): Let k := 0. Select and evaluate a set S 0 of starting points.

• While some given stopping criteria are not met:

-Phase 2 (model): From the data {(x, f (x)) | x ∈ S k }),
:= k + 1.
Phase 1 is commonly referred to as sampling or design of experiment (DOE) [START_REF] Forrester | Engineering Design via Surrogate Modelling, a practical guide[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF]. Its purpose is to find a set of points uniformly spread over the domain, so that, if we evaluate the function at these points, we might obtain a global picture of its range. To better achieve this a uniformity measure of each point set is employed. Our expected sampling is then found by maximizing uniformity (this is known as the design problem). In Section 2, various such measures as well as methods for solving associated design problems are introduced. We also mention the definition of stratified sample plans. Among those Latin hypercube [START_REF] Morris | Exploratory designs for computational experiment[END_REF] is the most popular.

In Phase 2, various models can be used to approximate black-box functions. In Section 3, we will introduce some of the most popular models, including polynomials, radial basis functions (RBF) and kriging. While polynomials are well-studied, their use as global models for high-dimensional black-box functions is only limited to linear and quadratic cases. Other models like RBF and kriging can fit more complicated functions while still being relatively simple to build. For this reason, many recent surrogatebased methods for black-box optimization make use either one of the two as interpolation models. It is also useful to use multiple surrogate models at the same time, because it can prevent us from choosing poorly fitted models [START_REF] Viana | Multiple surrogates: how cross-validation errors can help us to obtain the best predictor[END_REF].

Phase 3 is the crucial step in the procedure. Given the information from the current surrogate model, we need to decide which point(s) should be evaluated in the subsequent step. There are many strategies to do that, and indeed they represent the main feature distinguishing different surrogate-based optimization methods. The most common idea of these strategies is to define a merit (or cost) function of candidate point that predicts the objective and/or model accuracy improvement if we evaluate the black-box there. We select the next point for evaluation to be the one that maximizes (or minimizes) the merit function (resp. cost function). Different merit functions and methods for solving associated subproblems will be introduced in Section 4.

With the emphasis on expensive functions, the whole process for solving black-box optimization problems is often dominated by the number of function evaluations. Therefore, the performance of an algorithm is often measured by the number of evaluations needed until an acceptable (global) solution is found. In fact, it may also be used as a stopping criterion in black-box optimization: we stop when we reach a certain number of evaluations.

The above procedure is only a simplification of various surrogate-based methods for black-box optimization. In the literature, there are more complex frameworks. For example, surrogate models can be used in the poll step of a directional direct-search: Before the poll trial points are evaluated on the true functions, they are evaluated on the surrogates. This list of points is sorted according to the surrogate values so that the most promising points are evaluated first during the true evaluations [START_REF] Booker | A rigorous framework for optimization of expensive functions by surrogates[END_REF][START_REF] Conn | Use of quadratic models with mesh-adaptive direct search for constrained black box optimization[END_REF][START_REF] Digabel | Algorithm 909: NOMAD: Nonlinear optimization with the MADS algorithm[END_REF]. This is called surrogate management framework [START_REF] Booker | A rigorous framework for optimization of expensive functions by surrogates[END_REF], various formulations for it is proposed in [START_REF] Talgorn | Statistical surrogate formulations for simulation-based design optimization[END_REF]. In [START_REF] Regis | A quasi-multistart framework for global optimization of expensive functions using response surface models[END_REF], a quasi-multistart framework called AQUARS is proposed, which runs a multistart method to identify the local minima of the current surrogate model. At every iteration, AQUARS checks if a local minimum has been "fully explored", if this is the case, we do not need to carry out any further exploration in its vicinity. AQUARS also decides on which type of local search (breadth first or depth first) should be performed. Another framework that combines local and global methods is presented in [START_REF] Ji | A new framework for combining global and local methods in black box optimization[END_REF]. It consists of two procedures: the first constructs a global model to approximate the black-box function and explores an unvisited area to search for a global solution; the other identifies a promising local region and conducts a local search to ensure local optimality. At each iteration, we divide all existing data points into 2 sets: a full information set, which consists of all evaluated points so far, and a partial information set, which is the set of only evaluated points since the most recent Restart. Surrogate models are used to fit both these two types of sets to serve different purposes: we use full information model to switch back and forth between local and global searches and use partial information model to generate candidate points for global search.

Dealing with constraints and multiple objectives

Black-box optimization problems become much more difficult with the present of multiple black-box functions, such as in the case of multiple black-box objectives or constraints (we remark that constraints having inexpensive evaluation, such as when the analytical form is given, are implicitly included in the feasible set D). The most natural idea to deal with multiple expensive constraints is to combine them all into a weighted penalty term added to the objective function. Another idea is to aggregate all these nonlinear constraints into a single constraint violation function, and possibly even treat the optimization problem as an unconstrained bi-objective problem [START_REF] Audet | A pattern search filter method for nonlinear programming without derivatives[END_REF] or treat it as a single constraint [START_REF] Audet | A progressive barrier for derivative-free nonlinear programming[END_REF]. However, a more direct approach for solving these problems is to construct one surrogate model for each black-box function [START_REF] Akhtar | Multi objective optimization of computationally expensive multimodal functions with rbf surrogates and multi-rule selection[END_REF][START_REF] Digabel | Efficient global optimization algorithm assisted by multiple surrogate techniques[END_REF]. Using surrogate models for constraint functions helps us to predict feasible points, so we can avoid spending time to evaluate points that might be infeasible.

The rest of this paper is constructed as follows: In Section 2, we give an overview of the experimental design problem. In Section 3, we introduce some popular surrogate models used in practice. In Section 4, we review common choices of merit functions, associated search domains and their uses for solving blackbox optimization. Section 5 concludes the paper by summarizing the idea of surrogate-based methods and suggesting a few possible research directions.

Overview of experimental designs

The first step in any surrogate-based method is to find a set of starting points for evaluation. This step is commonly referred as design of experiments (DOE) [START_REF] Forrester | Engineering Design via Surrogate Modelling, a practical guide[END_REF][START_REF] Johnson | Minimax and maximin distance designs[END_REF][START_REF] Morris | Exploratory designs for computational experiment[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF]. Function values given at these points provide us the very first information about the black-box function, which is then used to construct suitable surrogate models. Good choices of experimental designs can help us to identify inadequacies in proposed models, as well as to avoid model bias [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF].

In the following, such a set of starting points will be called a design and the entire domain D will be called the design space. Moreover, we will mostly discuss computer experiments. The main difference between a physical and a computer experiment is that the latter is deterministic, i.e. replicate observations at the same input yield the same results [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. Computer experiments therefore lack random errors. Because of this, each input corresponds to a unique response. A word of caution about this statement is in order, however: in practice, computers can (and usually do) behave nondeterministically, as hardware glitches provide a constant source of errors -only the error correction codes save us from total chaos.

If no detail about the functional behavior of the response is available, we need to explore the entire design space. Therefore, we will need to find a design set so that a variety of models can be fitted and it must spreads over the experimental region [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. Moreover, if one of the input parameters has little influence on the function value, two points that differ only in this parameter should be considered the same. To avoid this situation, we require that no two design points share common value in any dimension.

From the above discussion, we summarize the two most important requirements for a good experimental design [START_REF] Husslage | Space-filling latin hypercube designs for computer experiments[END_REF][START_REF] Stinstra | Constrained maximin designs for computer experiments[END_REF]:

• space-fill: The design points should be uniformly spread over the entire design space.

• Non-collapse: Two design points should not share any coordinate value if we do not know a priori which dimensions are important.

Note that when the number of variables is high, preprocessing the problem to reduce its dimension is crucial. If all the parameters in the problem are selected to be important, the non-collapse requirement becomes less significant. In this case, designs with onlyspace-fil are easy to construct. However, many practical applications are often subject to several constraints, therefore, we list an additional requirement (which is also hard to deal with) for a design set

• Constraints-fil: The design points have to satisfy a set of linear and/or nonlinear constraints.

The concept of space-fill is based on spreading points uniformly over the design space. There are many definitions which result in different quality measures for a design [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. Some of the most popular measures in the literature, often based on either geometrical or statistical properties, are minimax, maximin, IMSE, maximum entropy, Audze-Eglais, and so on. The associated design problems consist in maximizing (or minimizing) these measures.

There is a trade-off between space-fill and non-collapse properties: a design with good space-fill properties is often collapsing (a striking example is full factorial designs). To obtain a good design in terms of both space-fill and non-collapse, we often solve the design problem over a reduced class of noncollapsing designs. Among such classes, Latin hypercube designs (LHD) [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] are the best known. On the other hand, some point generation techniques, such as Sobol' sequence, avoid this trade-off by uniformly filling the space with non-collapsing point sets.

In fact, finding an optimal design is a very complicated problem, especially in high dimension. For example, with a relatively simple measure among them, the maximin LHD problem is believed to be NP-hard [START_REF] Fang | Centered ℓ 2 -discrepancy of random sampling and latin hypercube design, and construction of uniform designs[END_REF] (to the best of our knowledge, the question is still open).

In the next subsection, we will introduce various criteria to measure the space-fill and/or non-collapse of a designs. After that, we overview common methods to find optimal designs with respect to such criteria.

Space-filling and non-collapsing designs

In this subsection, for simplicity, we assume that the design space D is a compact subset of R d and is defined by box constraints. More specifically, let

D = {x ∈ R d | L i ≤ x i ≤ U i }.
In most cases, we assume that all design parameters are equally important, therefore we can scale each interval [L i , U i ] linearly to [0, 1] for all dimensions i. This subsection reviews the most common space-filling and non-collapsing designs over D. Various criteria based on the same idea for selecting optimal designs are introduced: first, define a measure Φ(.) of space-fill and then optimize Φ(S) over all S ⊆ D with |S| = n. When we 6 restrict the search space to a certain class of non-collapsing designs, we can find designs that both satisfy the two requirements. There are many choices for Φ(.) and each choice corresponds to a different type of designs. We categorize them into two main groups: geometrical designs and statistical designs.

Geometrical designs

Factorial designs: Full factorial designs might be the simplest and most straightforward way to sample points in the uniform manner [START_REF] Forrester | Engineering Design via Surrogate Modelling, a practical guide[END_REF]. They are constructed by dividing each dimension i of the design space D ⊆ R d into n i equal intervals and then selecting centers of all n = n 1 n 2 . . . n d resulting cells. Obviously, full factorial designs are highly uniform but lack of the non-collapse property. Another drawback is their inflexibility regarding the number of sampling points. In particular, a full factorial design can only be used if the number of sampling points is of the form n = n 1 n 2 . . . n d , which is often very large (n ≥ 2 d ). This almost prohibits the use of full factorial designs for sampling sets of size that is small or given in advance. Sometimes, we use a variation called fractional factorial design, in which only a fraction of the design points specified by the full factorial designs are considered.

Latin hypercube designs: Latin hypercube designs (LHD) [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] are the most popular non-collapsing designs in the literature. An LHD with n points is constructed as follows: we first divide each axis into n equal intervals: now the design space D will consist of n d identical cells. An LHD is then obtained by assigning n points to centers of these n d cells such that there are no two points with the same coordinate in any dimension. One point can be conveniently represented as an element in {1, 2, . . . , n} d (not to be confused with its coordinate). Therefore, we can interpret an LHD as a set of n points x 1 , . . . , x n ∈ {1, 2, . . . , n} d with the property that, for each j ∈ {1, . . . , d}, the set {x 1j , . . . , x nj } is a permutation of {1, 2, . . . , n}. Formally, let S = {x 1 , . . . , x n } be a set of n points in {1, 2, . . . , n} d , we define the design matrix of S as follows: Then S is an LHD if each column of M (S) is a permutation of {1, 2, . . . , n}. Note that an LHD is not necessary space-filling. In Figure 1, two different LHDs are presented but in the first design, all points are distributed along the main diagonal (so they do not spread uniformly over the design space). The definition of minimax distance designs is very intuitive: we want each point in D to be presented (i.e. be close to) by at least one of the selected points. Therefore, the information about that point can be exploited through its representative. In particular, in the minimax design problem, we want to choose a minimal radius r > 0 and a set of centers x 1 , . . . , x n such that the design space D is covered by the union of all the spheres B(x i , r), i = 1, . . . , n. Therefore, minimax distance designs closely relate to the covering problems. They are, however, are computationally difficult to find, because for a given design, an optimization step is required to calculate the value of the measure max x∈D d(x, S) [START_REF] Stinstra | Constrained maximin designs for computer experiments[END_REF].

Maximin (Mm) distance design [START_REF] Johnson | Minimax and maximin distance designs[END_REF]: Similarly, a maximin distance design S * of cardinality n is the one that solves max

|S|=n min x,y∈S d(x, y).
The definition of maximin distance designs is also very intuitive: since the number of design points is often relatively small and the design space is often very large, we do not want to evaluate two points that are too close to each other. The maximin distance design problem is similar to the packing problem, in which we need to choose a maximal radius r > 0 and a set of centers x 1 , . . . , x n such that all the spheres B(x i , r), i = 1, . . . , n are pairwise disjoint. Moreover, the space-fill measure in this case, i.e. min x,y∈S d(x, y), is quite simple. For that reason, M m distance designs are the most well-studied among all experimental designs [START_REF] Forrester | Engineering Design via Surrogate Modelling, a practical guide[END_REF][START_REF] Husslage | Space-filling latin hypercube designs for computer experiments[END_REF][START_REF] Stinstra | Constrained maximin designs for computer experiments[END_REF].

Φ p designs [START_REF] Morris | Exploratory designs for computational experiment[END_REF]: One problem with maximin distance designs is it only considers the smallest distance of a point set and neglects other distances. Another design [START_REF] Morris | Exploratory designs for computational experiment[END_REF], makes all distances of the point set counted. In this design, we have to minimize:

Φ p (S) = x,y∈S 1 d(x, y) p 1/p
, where p is a parameter specified by users. Note that the choice of p is very important. When p is large enough, the contributions of largest distances become the dominating factors and Φ p designs are equivalent to minimax distance designs. However, when p gets larger, the optimization problem becomes more difficult. For lower values of p, the measure Φ p might not exactly be the same as the minimax design, but it is more amenable to optimization methods [START_REF] Forrester | Engineering Design via Surrogate Modelling, a practical guide[END_REF]. Note that when p = 2, we have the Audze-Eglais distance design (also called potential energy) [START_REF] Audze | New approach for planning out of experiments[END_REF]:

S * ∈ arg min |S|=n x,y∈S 1 d(x, y) 2 .
This design is favorable because when dealing with euclidean metrics, the quantities d(x, y)2 will cancel out all square roots. Therefore, the associated objective function becomes analytically simple. Moreover, when minimizing this measure over the reduced class of LHDs, we can obtain designs with very good properties and useful for many applications [START_REF] Bates | Formulation of the optimal latin hypercube design of experiments using a permutation genetic algorithm[END_REF][START_REF] Husslage | Space-filling latin hypercube designs for computer experiments[END_REF][START_REF] Stocki | A method to improve design reliability using optimal latin hypercube sampling[END_REF].

Statistical designs

In statistical designs, we treat a deterministic response of a black-box function as if it were a realization of a stochastic process X. We assume X has mean zero and correlation function R with the property that for all x i , x j ∈ D, R(x i , x j ) increases when d(x i , x j ) decreases (i.e. when x i and x j are close, their correlation is high and vice versa). Many other functions of this kind can be found in [START_REF] Johnson | Minimax and maximin distance designs[END_REF]. One example of such a correlation function is R(x i , x j ) = exp -θ x i -x j for some θ > 0. In [START_REF] Sacks | Design and analysis of computer experiments[END_REF], a different choice of R is given by

R(x i , x j ) = exp d k=1 -θ k |x ik -x jk | t , 1 ≤ i, j ≤ n,
for some θ k > 0 and 1 ≤ t ≤ 2. These functions are more convenient than the former, since we can express them as products of one-dimensional correlations, while still provide enough flexibility for adequate prediction in most cases. We will explain more about this function later. Now we want to choose the design set S such that the response X can be best predicted. This is the subject of classical optimal design. Many different criteria for choosing optimal designs which exploit the information about the assumed correlation R(., .) have been proposed. We will not go into details about all of them, instead we discuss some relatively intuitive criteria being proposed by [START_REF] Johnson | Minimax and maximin distance designs[END_REF]. These criteria rely on the predictors X(x) at each x ∈ D, which are found based on the assumed data at {x 1 , . . . , x n }. Details on how to find these predictors will be presented in Section 3.

I-optimality:

We are given a weight function φ(s) and let Var X(s) be the variance of the prediction X(s) (note that both X(s) and X(s) are random). A design S * is called I-optimal if it solves min |S|=n D Var X(s)φ(s) ds.

In [START_REF] Sacks | Design and analysis of computer experiments[END_REF], the function being minimized is called integrated mean squared error (IMSE). When D is finite, we often select φ(s) ≡ 1, so the objective becomes si∈D Var X(s i ), which is the total error associated with all the prediction points.

G-optimality:

Similar to I-optimality, a design S * is called G-optimal if it solves min |S|=n max s∈D Var X(s) ,
i.e. minimizes the maximum error of all predictions. In [START_REF] Sacks | Design and analysis of computer experiments[END_REF], the function being minimized is called maximum mean squared error (MMSE).

D-optimality:

For any design set S = {x 1 , x 2 , . . . , x n }, let

R(S) =   Corr (X(x 1 ), X(x 1 )) . . . Corr (X(x 1 ), X(x n )) . . . . . . . . . Corr (X(x n ), X(x 1 )) . . . Corr (X(x n ), X(x n ))   be the correlation matrix of S. A design S * is called D-optimal if it solves max S det R(S),
i.e. maximizes the determinant of the correlation matrix R(S). Note that large values of det R(S) imply the low correlations of all random variables X(x i ). Since R(x, y) decreases when the distance d(x, y) increases, in general, D-optimal designs are uniformly distributed. Moreover, D-optimality is much easier to deal with than G-optimality and I-optimality since it does not involve predictions at any x ∈ D.

Another way to explain D-optimality is via the concept of entropy: we want to find design sets that minimize the expected posterior entropy, i.e. the sets of points at which we have the least information. Shewry and Wynn [START_REF] Shewry | Maximum entropy sampling[END_REF] showed that this is equivalent to minimizing -log (det R(S)), which is obviously equivalent to finding D-optimal designs.

Methods for finding optimal designs

From the previous subsection, we know that a general design problem can be written as min S⊆D,|S|=n where Φ(.) is a measure that quantifies the non-uniformity of a set. To get the balance between spacefilling and non-collapsing properties, we also consider Latin hypercube design problems, in which the search space for S is restricted to the class LD of all Latin hypercube designs: min S⊆LH, |S|=n Φ(S).

(3)

In the previous subsection, we have introduced many design measures Φ(.) but have not mentioned how to solve the problems (2), [START_REF] Audet | A survey on direct search methods for blackbox optimization and their applications[END_REF]. Before moving on to discuss methods for solving them, we should note that it is not obvious which measures should be used in design problems. A design found by optimizing one criterion is not necessarily good at the others. Moreover, design problems are often very hard to solve. Moreover, optimality guarantees are not required. Therefore, heuristic methods, such as evolutionary algorithms, may be more practical. For example, in [START_REF] Bates | Formulation of the optimal latin hypercube design of experiments using a permutation genetic algorithm[END_REF], the authors apply a permutation genetic algorithm to find good Audze-Eglais LHDs. Results for eight different combinations of n and k are also reported [START_REF] Bates | Formulation of the optimal latin hypercube design of experiments using a permutation genetic algorithm[END_REF][START_REF] Husslage | Space-filling latin hypercube designs for computer experiments[END_REF].

In principle, we can apply standard nonlinear programming (NLP) methods to find optimal designs. Most of the previously mentioned design problems are already in the form of NLPs, so we can simply plug them into a solver. For example, in the maximin distance design problem, we exploit the idea of maximizing the radius of n identical balls centered in D which form a packing to obtain the following NLP:

max {r | x i ∈ D, d(x i , x j ) ≥ 2r for all i, j}. (4)
If n is small, this problem can be handled by powerful NLP solvers. However, when n gets larger, these NLPs become very difficult due to their nonconvexity and the number of constraints [START_REF] Stinstra | Constrained maximin designs for computer experiments[END_REF].

Local improvements

In experimental design problems, we are required to optimize the locations of n points. Therefore, it is natural to use the idea of local improvement: we start with a particular design and iteratively improve it by optimizing the location of one point while fixing n -1 others. The problem of nd variables will now be converted into subproblems of only n variables.

This idea can be implemented for many design problems. For example, for Φ p family, we can consider the following subproblems

min xi∈D j =i 1 d(x i , x * j ) p
in which x * 1 , . . . , x * n is the current design set. In [START_REF] Stinstra | Constrained maximin designs for computer experiments[END_REF], the NLP (4) is replaced by the following NLP subproblems, which are still non-convex, but generally much easier to solve:

max xi∈D {r | d(x i , x * j ) ≥ 2r for all j = i}.
Note that at each step we have to select a suitable index i to process. We can sequentially select i from any permutation of {1, 2, . . . , n}, but a smarter choice will probably better perform. For example, in [START_REF] Stinstra | Constrained maximin designs for computer experiments[END_REF] an index i will be skipped if the previous iteration shows no significant improvement with respect to x i .

Column exchange

Column exchange is a special local technique applied for Latin hypercube design problems. Instead of locally modifying one point from the current design, we look for improvements by perturbing one or several columns of the design matrix. Recall that a LHD matrix is a matrix in which each column is a permutation of {1, 2, . . . , n}. Therefore if we replace one column by another permutation of {1, 2, . . . , n}, we still obtain a LHD matrix.

In [START_REF] Morris | Exploratory designs for computational experiment[END_REF], the authors use a simulated annealing method to optimize LHDs. The method starts with a certain Latin hypercube design S. At each iteration, a column of the design matrix M (S) is selected and a new design S ′ is generated from S by interchanging two random entries in that column. If the objective function Φ is improved, S will be replaced by S ′ . Otherwise, we decide whether to keep S or update S to S ′ by computing a probability given by exp Φ(D ′ )-Φ(D) t , where t is an algorithm parameter often known as temperature. In [START_REF] Grosso | Finding maximin latin hypercube designs by iterated local search heuristics[END_REF], an Iterated Local Search (ILS) heuristics is proposed, in which column exchanges are used in both local search steps and perturbation steps. In the local search step, we define a set of critical points and the swaps are only made between components of these points. In the perturbation step, cyclic order exchanges are performed to get rid of local optimum.

Since the class of all LHDs is quite large (there are (n!) d-1 possibilities), a restricted class of symmetric Latin hypercube designs (SLHD) with desirable properties is introduced in [START_REF] Ye | Algorithmic construction of optimal symmetric latin hypercube designs[END_REF]. A design S is called symmetric if the sums of the i th and (n + 1 -i) th row vectors are always equal to (n + 1, . . . , n + 1) for all i. SLHDs are claimed to work better than regular LHDs in the cases of entropy or maximin measures [START_REF] Ye | Algorithmic construction of optimal symmetric latin hypercube designs[END_REF]. In [START_REF] Ye | Algorithmic construction of optimal symmetric latin hypercube designs[END_REF], an columnwise-pairwise algorithm is proposed for finding optimal designs within this class. At each iteration, two pairs of entries (instead of one) are exchanged to ensure that the resulting design is still symmetric. These points are not randomly chosen as in [START_REF] Morris | Exploratory designs for computational experiment[END_REF], but they are selected as the best two simultaneous exchanges.

Surrogate models

This section introduces some of the most popular surrogate models being used in engineering designs. We begin with polynomials, a class of functions that are well-studied. They are, for example, used in trust-region methods to provide approximation of the true function in local areas. However, they are unsuitable as global models for highly nonlinear, multidimensional functions. Other models that are better performed, including radial basis functions and kriging will be introduced in the subsequences.

Polynomials

Let P m be the space of all polynomials of degree at most m on x ∈ R d . Let {φ 1 , φ 2 , . . . , φ M } be a basis of P m , then each polynomial on P m can be written as

p(x) = M i=1 a i φ i (x) (5)
where a i ∈ R, i = 1, 2, . . . , M. An example of such a basis is the natural basis φ α (x) = x α := x α1 1 . . . x αd d , where α j ∈ N for all 1 ≤ j ≤ d and α 1 + . . . + α d ≤ m. Assume that we want to find a polynomial of the form (5) that interpolates the data (x 1 , y 1 ), . . . , (x n , y n ), then we need to solve the system of equations:

p(x i ) = y i
for all i = 1, 2, . . . , n.

Denote by

M φ,x =   φ 1 (x 1 ) . . . φ M (x 1 ) . . . . . . . . . φ 1 (x n ) . . . φ M (x n )   , a =   a 1 . . . a M   and y =   y 1 . . . y n   ,
then the above system can be written as M φ,x a = y. This system has a solution if the matrix M is full-row rank. It is well-known that when n = M , the existance and uniqueness of the interpolation does not depend on particular choices of the basis {φ i } i (see the proof in [START_REF] Conn | Introduction to derivative-free optimization[END_REF]).

Lemma

Let x 1 , . . . , x n be such that M φ * ,x is a non-singular square matrix for some basis φ * = {φ * 1 , . . . , φ * M } in P m . Then for any choice of basis φ, the system M φ,x a = y has an unique solution.

The above lemma, which requires n = M , offers no flexibility on choosing the number of evaluation points n. In general, when n < M , we need to restrict the search to a subspace V ⊆ P d in order to ensure the unique existance of the interpolated polynomial. The theory of polynomial interpolation goes beyond the scope of this paper, so we recommend the survey [START_REF] Gasca | Polynomial interpolation in several variables[END_REF] by M. Gasca and T. Sauer for interested readers. The primary weakness of polynomials are their inability to fit smooth functions of any shapes [START_REF] Hussain | Metamodeling: Radial basis functions versus polynomials[END_REF]. Moreover, they require a large number of parameters to be estimated (due to the high dimension of the space P d ). Polynomial model, however, is used for interpolation in practice due to the wide availability of techniques and softwares for computing polynomials.

Radial basis function

A radial function φ is the one whose value at each point depends only on the distance between that point and the origin, i.e. φ(x) = ξ( x ) where ξ is a real function. Examples of such functions are

• Linear: φ(x) = x • Cubic: φ(x) = x 3 • Thin plate spline: φ(x) = x 2 log x • Multiquadric: φ(x) = x 2 + γ 2
• Gaussian: φ(r) = e -γ x 2 where γ is a prescribed positive constant.

Given a set of points S = {x 1 , x 2 , . . . , x n }, x i ∈ R d and their corresponding responses y 1 , y 2 , . . . , y n ∈ R, we want to interpolate the data (x 1 , y 1 ), (x 2 , y 2 ), . . . , (x n , y n ) by functions of the form

s n (x|S) = n i=1 λ i φ( x -x i ) + p(x) (6)
where p is a polynomial and φ is one of the radial basis functions listed above. Let denote by p = (p(x 1 ), . . . , p(x n )), y = (y 1 , . . . , y n ), λ = (λ 1 , . . . , λ n ) and Φ = φ( x i -x j ) 1≤i,j≤n , then the interpolation condition system can be written as Φλ + p = y. The polynomial p appears here to ensure the existence of at least one solution for this system (otherwise, when Φ is nonsingular, the system Φλ = y has no solution).

In [START_REF] Gutmann | A radial basis function method for global optimization[END_REF], all RBFs listed above are discussed (using the concept of conditionally positive/negative definiteness), but later the author claimed that the multiquadric and Gaussian cases are disappointing (according to [START_REF] Björkman | Global optimization of costly nonconvex functions using radial basis functions[END_REF]). If φ(x) is either linear, cubic or thin plate spline, the polynomial p can be taken in the form of p(x) = b T x + a (or just a constant p(x) = a in the case φ(x) is linear). The unknown parameters λ 1 , . . . , λ n , a, b 1 , . . . , b d are obtained by solving the following system of linear equations

Φ P P T 0 λ c = y 0 , where P =   x T 1 1 . . . . . . x T n 1   and c =     b 1 . . . b d a     .
When rank (P ) = d+1, the matrix Φ P P T 0 is nonsingular, so the above equation has a unique solution.

Let Π 1 be the set of linear polynomials and V be the space of all λ ∈ R n that satisfy

n i=1 λ i p(x i ) = 0 for all p ∈ Π 1 .
For all RBFs being listed above, there is a nonnegative integer m 0 such that (-1) m0 λ T Φλ > 0 for all λ ∈ V. Let A φ be the set of functions that are defined as in [START_REF] Audet | Mesh adaptive direct search algorithms for constrained optimization[END_REF] with λ ∈ V. For any s and u in A φ , say

s(x) = n1 i=1 λ i φ( x -x i ) + p(x) and u(x) = n2 j=1 µ j φ( x -z j ) + q(x),
we define a semi-inner product of s and u by

s, u := (-1) m0 n1 i=1 n2 j=1 λ i µ j φ( x i -z j ).
As usual, s = s, s also defines a semi-norm on A φ and we have

s 2 = s, s = (-1) m0 λ T Φ(x) λ.
When a lower bound on the black-box function is known, the RBF model can be modified to justify this information. If we know that y is a lower bound of function f (x), then any surrogate model whose minimum value drops below y should be avoided. A procedure for generating lower-bounded interpolants is proposed in [START_REF] Cassioli | Global optimization of expensive black box problems with a known lower bound[END_REF]. The procedure will augment the sample set S n by a set of additional points S f that are discarded once a new point x n+1 is selected and f (x n+1 ) evaluated.

Assume that x is a global solution for the current model s n (x|S n ). If s n (x|S n ) ≥ y, then our constructed model satisfies the required lower bound y and we can use it without any further modification. Otherwise, x is temporarily grouped into our design set S n . Now, s n (x|S n ) is replaced by a new model, which is obtained as follows: First, we iteratively solve min λ,c,ǫ

(-1) m0 λ T Φ(x) λ subject to: Φ(x)λ + P c - 0 ǫ = y y P λ = 0 ǫ ≥ 0.
This quadratic convex problem is well-studied and can be solved efficiently by many current methods.

When a new surrogate model is constructed, we can find its minimal value and re-check whether that value is at least y or not. If it is not the case, then the process is repeated by inserting an additional point and constraint to increase the lower bound of the interpolation model. We continue this procedure until we reach a maximum number of iterations, and then replace s n (x|S n ) by the last interpolation model.

Kriging

Kriging is a special case of Gaussian processes. The use of kriging in black-box optimization begins with the famous paper by Sacks et al. [START_REF] Sacks | Design and analysis of computer experiments[END_REF]. In that paper, each deterministic function value y(x) is treated as a realization of a stochastic regression process Y (x) = h j=1 β j f j (x) + Z(x), in which the random process Z(.) is assumed to have mean zero and covariance

Cov (Z(x), Z(s)) = Φ 2 R(x, s),
where Φ 2 is the process variance and R(x, s) is the correlation. For each untried point s, the value of Y (s) is predicted as the best linear unbiased predictor (BLUP). This is equivalent to minimizing the mean squared error (MSE) subject to the unbiasedness constraint, i.e.

minimize c(x) E[c(x) T Y S -Y (x)] 2 subject to: E[c(x) T Y S ] = E[Y (x)].
By introducing Lagrange multipliers λ(x) for the unbiasedness constraint, we can easily solve this problem. Denote by

f (x) =    f 1 (x) . . . f h (x)    , F =    f 1 (x 1 ) . . . f h (x 1 ) . . . . . . . . . f 1 (x n ) . . . f h (x n )    , r(x) =    R(x 1 , x) . . . R(x n , x)    , R = [R(x i , x j )] 1≤i,j≤n ,
then a convenient representation for the mean squared error of a predictor ŷ(s) is given by [START_REF] Sacks | Design and analysis of computer experiments[END_REF]:

MSE (ŷ(s)) = Φ 2 1 -f (x) T , r(x) T 0 F T F R -1 f (x) r(x) (7)
The correlation R(x, s) has to be specified to compute this MSE. In [START_REF] Sacks | Design and analysis of computer experiments[END_REF], the authors restrict their attention to a special correlation class of the form

R(x, s) = d j=1 exp(-θ j |x j -s j | pj ) = exp - d j=1 θ j |x j -s j | (8)
where θ j > 0 and 0 < p j ≤ 2. According to [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF], these correlations have the property that if x = s then R(x, s) = 1. Similarly, when x i -x j → ∞, they tend to zero. The parameter θ j determines how fast the correlation decreases in the j th coordinate direction. The parameter p j determines how smooth the function is in the j th coordinate direction.

In [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF], a more natural way that leads to kriging is introduced. In that paper, the authors treats y(x) as a realization of a normally distributed random variable Y (x) with mean µ and variance Φ 2 . The correlation R(x, s) is specified the same as in equation ( 8). However, the parameters µ, Φ 2 , θ j , p j are chosen to maximize the likelihood of the observed data, which is

1 (2π) n 2 (Φ 2 ) n 2 (det R) 1 2 exp -(y -1µ) T R -1 (y -1µ) 2Φ 2
where y = (y 1 , y 2 , . . . , y n ) T . Taking logarithm and using first order conditions, we can estimate μ and Φ2 as functions of det R and y. To predict the response at an untried point x * ∈ D, we choose a predicted value y * to maximize the augmented likelihood function. It is shown to be equivalent to minimizing

y -1μ y * -μ T R r r T 1 -1 y -1μ y * -μ .
By applying the partitioned inverse formula, this function can be written as a quadratic function of y * . Therefore we obtain a standard formula for the kriging predictor

y * = μ + r T R -1 (y -1μ). (9)

Support Vector Machine

Another surrogate model that is also used for engineering design and black-box optimization is Support Vector Machine (SVM). SVM is the model used very often in supervised learning to help analyzing classification and linear regression problems. Intuitively, it constructs a hyperplane with maximal margin that separates the training data. The use of SVM for black-box optimization is suggested in [START_REF] Nakayama | Using support vector machines in optimization for black-box objective functions[END_REF], in which the authors report promising computational results.

In the simplest case, a SVM model is constructed as the following linear regression:

f (x) = a T x + b,
in which a, b are the parameters that minimize the quadratic lost function

1 2 a 2 + C n i=1 a T x i + b -y i 2 ε
where z ε := max{0, |z| -ε}.

In order to fit nonlinear functions, the data (x i ) 1≤i≤n are embedded into a higher dimensional space by a mapping φ (e.g. polynomials, Gaussian RBF, Hyperbolic tangent, Kriging, . . . ). Then we can rewrite the SVM problem in the following form:

min 1 2 a 2 + C n i=1 (µ 2 i + ξ 2 i ) subject to a T φ(x i ) + b -y i ≤ ε + µ i y i -a T φ(x i ) -b ≤ ε + ξ i µ i , ξ i ≥ 0.
By rewriting this problem in the dual form and using what is called the kernel trick K(x i , x j ) = φ(x i ) T φ(x j ), we obtain the following max 1 2 n i,j=1

(α + i -α - i )K(x i , x j )(α + j -α - j ) - ε 2 n i=1 (α + i + α - i ) + n i=1 y i (α + i -α - i ) subject to n i=1 (α + i -α - i ) = 0 0 ≤ α + i , α - i ≤ C n .
This is a quadratic problem and can be solved efficient by current methods such as coordinate descent.

Mixed surrogate models

The models presented above have their own advantages and disadvantages, depending on the specific applications that we consider. A certain model that is a perfect choice for some problems might perform poorly on others. It is natural to ask, given a black-box problem, what is the best surrogate model to use. In [START_REF] Müller | Mixture surrogate models based on dempster-shafer theory for global optimization problems[END_REF], the authors study the influence of combining multiple surrogate models on the performance of methods for solving black-box optimization. They define a mixture surrogate model as a convex combination of individual surrogate models, i.e.

s mix (x) = r∈M w r s r (x), where r∈M w r = 1, and w r ≥ 0, in which M is the set of all available surrogate models. For each r ∈ M , the large w r means that s r (x) is a good choice as a surrogate model for the specific problem. They propose to apply Dempster-Shafer theory in order to find the most suitable model as well as the best combination of weights. In [START_REF] Müller | Influence of ensemble surrogate models and sampling strategy on the solution quality of algorithms for computationally expensive black-box global optimization problems[END_REF], they examine various surrogate models and their combinations within the so-called SO-M-c and SO-M-s frameworks.

A Matlab toolbox called MATSuMoTo [START_REF] Mueller | Matsumoto: The matlab surrogate model toolbox for computationally expensive blackbox global optimization problems[END_REF] is provided so that users can freely choose various models and sampling strategies from its library to test the algorithm. Moreover, this toolbox also enable users to choose experimental design strategies and make use of MATLAB Parallel Computing Toolbox to perform multiple function evaluations in parallel.

Surrogate-based subproblems

We have mentioned in the introduction part that, at each iteration of surrogate-based algorithms, we have to find the next candidate point(s) for evaluation. Their selection is crucial. Even if these points are not better than the previous ones, they provide us more information about the black-box function. This information is of great importance since in black-box optimization, all we know is the values at a very limited number of data points. Thus, the more data points we have, the more we know about the function.

Since we want the function value at each candidate to be improved, we have to focus the search around the current best point. But we also want to search in areas that are far from all previous points to ensure that we do not overlook a global solution. The two types of searches above are respectively called local and global ones. In black-box optimization, they are also mentioned as exploitation and exploration. Both these two processes are indispensable: if we spend too much time on exploitation, we might get trapped at a local solution; if we spend too much time on exploration, we might not reach any optimal solution at all, even a local one. The management of the two processes (i.e. how to balance them), is at the center of surrogate-based algorithms.

There are several ideas for the selection of candidate points to achieve this balance. The first idea is to keep exploiting until there is no sign of improvement and then move to a new unexplored area. The second idea is to separate in advance how many points we will use to search locally and globally at each iteration. We can also look for compromises, i.e. the searches that are neither local nor global, but somewhere "in between".

In black-box optimization, new candidate points are often chosen by solving auxiliary subproblems. We define these problems by introducing a merit function and looking for the maximum or minimum of this function. But an auxiliary subproblem defined by a single merit function over the whole domain D is usually not enough for searching both locally and globally. In the literature, the most common approach is to use not only one but several merit functions and/or adaptively change the domain of consideration (which we will call search regions). However, instead of using different types of merit functions, we can use a unique type and parameterize it as σ(x|λ), in which a small λ correspond to a local search and when λ increases, the searches become more global. In this section, we will introduce various choices of such (parameterized) merit functions and their associated search regions.

Choices of merit functions

In the following, we will review some popular merit functions. We first consider the class of merit functions that make use of standard error in statistical models. We then consider RBF-based merit functions like bumpiness and its variations. Lastly, general merit functions which do not depend on any particular surrogate model are introduced.

• Surrogate models as merit functions. The most intuitive way for choosing a new candidate point is to choose it as a minimizer of the current surrogate model. One advantage of this choice is that the model itself is normally relatively simple to optimize, so we can quickly solve auxiliary subproblems. This choice works well if the model is reliable, but it can be very misleading if the model is biased. In fact, this method does not converge to a global solution. Furthermore, if one of the data points is a minimizer of the current surrogate model, the method might converge to a point that is not even a local optimum. In this case, a minor modification has to be made: in order to require the gradient of surrogate model to agree with that of the black-box function, we need to sample in a small neighborhood of the current solution [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF][START_REF] Koehler | Computer experiments[END_REF]. Note that with proper choices of search regions, we can still use surrogate models as merit functions to find global solutions. For example, the CORS algorithm in [START_REF] Regis | Constrained global optimization of expensive black-box functions using radial basis functions[END_REF] looks for candidate points by minimizing surrogate models over search regions of the form:

D n := D \ n i=1 B(x i , r i ).
In that algorithm, B(x, r) = [x -r, x + r] and all radius r i have the same value r i = β n . Here β n are determined by cycles of N + 1 iterations where each cycle represent a range of values which starts at a large value close to 1 (global search) and gradually decreases to β n = 0 (local search). This method is able to find global minimizers since the iterates generated by the algorithm are proved to be dense [START_REF] Regis | Constrained global optimization of expensive black-box functions using radial basis functions[END_REF].

• Statistical lower confidence bound. Using kriging, we can measure the standard error between the true black-box and the surrogate functions. In particular, we can associated s n (x) with a statistical lower confidence bound of the form:

Lcb (x) := s n (x) -bσ M SE(x)
where σ is the process variance being estimated as in [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF] and b is a factor determined by users. This merit function is used in the so-called Sequential Design for Optimization (SDO) algorithm [START_REF] Cox | SDO: A statistical method for global optimization[END_REF][START_REF] Cox | A statistical method for global optimization[END_REF].

The Lcb has the property that it coincides the black-box function at all data points (since the mean squared errors at those points are all equal to zeros); the values of Lcb at other points become smaller and smaller comparing to the values of the surrogate function as the errors increase (since we hope the black-box function values to be small at uncertain points). Therefore, the Lcb function seems to be more promising than the surrogate model itself.

Since Lcb is quite difficult to minimize, in the SDO algorithm, we only compare the values of Lcb among a finite number of prediction points. The algorithm proceeds until a (user-specified) maximum number of function evaluations or when we reach the convergence given by y min < min i Lcb (x i ). Here y min is the current best function value and the minimum is taken over all possible remaining candidate points [START_REF] Cox | SDO: A statistical method for global optimization[END_REF]. The SDO algorithm, however, concentrates too much on local searches and can potentially delete some regions of search space, so it might not be able to find global minimizers [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF]. Therefore, if we want to use it we have to appropriately manage the search regions (as mentioned previously).

• Probability of improvement. Assume that T is a real number that is strictly smaller than the current best function value. We associate each x ∈ D with the probability that the prediction Y (x) of x (which is a normally distributed random variable) is smaller than T . This is called the probability of improvement (PI) at x and is formally defined as follows [START_REF] Kushner | A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise[END_REF][START_REF] Mockus | Application of bayesian approach to numerical methods of global and stochastic optimization[END_REF][START_REF] Žilinskas | A review of statistical models for global optimization[END_REF]:

P I(x) := P [Y (x) ≤ T ] = Φ T -μ(x) σ(x) ,
here Φ is the normal cumulative distribution function. The next candidate point is chosen to maximize this PI function. More specifically, [START_REF] Mockus | Application of bayesian approach to numerical methods of global and stochastic optimization[END_REF] discusses the general Bayesian approach in which the distribution is taken arbitrarily and [START_REF] Žilinskas | A review of statistical models for global optimization[END_REF] introduces an axiomatic treatment of the method, which he calls the P -algorithm.

The key advantage of PI over Lcb is that, under certain mild assumptions, the iterates are dense [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF].

The reason is that, when more and more points are sampled around the current best point, the standard error σ(x) gets smaller, therefore Φ T -μ(x) σ(x)

will become so small that we have to switch to unexplored regions. However, this PI function is extremely sensitive to the selection of targets T [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF]. If T is too close to the current best value y min (i.e. T is large), the search will be highly local and only focus in a small neighborhood of the current minimizer. Otherwise, if T is too far from y min (i.e. T is too small), the search will be highly global. Therefore, proper selections of T will help us to balance between local and global searches. A simple but powerful idea is to use different values of T , ranging from large to -∞, and find all candidate points from resulting subproblems.

• Expected Improvement. Assume that we use kriging to interpolate the black-box function. For any x ∈ D, let Y (x) be the prediction of the black-box function value at x. The improvement at x is then defined by I(x) = max y min -Y (x), 0 , where y min is the current best value. Therefore, the improvement at each point is a nonnegative random variable, thus we can compute its expectation as follows:

EI(x) := E(I(x)) = E max (y min -Y (x), 0) = [y min -μ(x)]Φ y min -μ(x) σ(x) + σ(x)φ y min -μ(x) σ(x) ,
where Φ and φ are the normal cumulative distribution and density functions, respectively. EI(x) is called the expected improvement at x [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]. The next candidate point is chosen to maximize the EI function.

Expected improvement is used in the EGO algorithm by Jones et al. [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]. In EGO, auxiliary subproblems are solved using a branch-and-bound algorithm. They use the fact that the EI function is monotonic in μ and in σ to simplify the computation of its lower and upper bounds in each rectangular. Under mild assumptions, the iterates from EGO are dense [START_REF] Vazquez | Convergence properties of the expected improvement algorithm with fixed mean and covariance functions[END_REF], so it can be used to find a global optimum. However, the method requires fairly exhaustive search around the initial best point before the algorithm begins to search more globally (because we treat the estimated standard error as if it is correct). Therefore EGO runs relatively slowly compared to modern methods.

• Bumpiness of functions. The RBF-based algorithm proposed by [START_REF] Gutmann | A radial basis function method for global optimization[END_REF] is one of the most popular methods in black-box optimization. It is based on the observation that, among different surrogate functions to approximate a black-box, in practice, the smoothest one is often the most accurate. Assume that we use certain radial basis functions for interpolation. Given data points (x 1 , y 1 ), . . . , (x n , y n ), we have for each additional (x * , y * ) ∈ D × R a different radial basis function.

The idea is, if we assign y * to a fixed value T , we will look for x * such that the resulting function is the "least bumpy". Here T can be regarded as an estimate of the optimal value and we expect that it is attained at x (i.e. f (x * ) = T ). Figure 2: An example of two surrogate models that interpolate 4 data points. The solid-line model is more likely than the dashed-line one since it is less bumpy [START_REF] Costa | RBFOPT: an open-source library for black-box optimization with costly function evaluations[END_REF].

For each w / ∈ {x 1 , . . . , x n }, let s w be the radial basis function defined by the following interpolation conditions:

s w (x i ) = y i for all i = 1, 2, . . . , n.

s w (w) = T. (10)
Then the next candidate point is selected by

x n+1 ∈ arg min σ(s w ) | w / ∈ {x 1 , . . . , x n } ,
where σ is a measure of the bumpiness of functions, which is defined by the semi-norm σ(s w ) = s w = s w , s w (see Section 3.2). In [START_REF] Gutmann | A radial basis function method for global optimization[END_REF], the bumpiness function σ(s w ) is simplified to

g n (w) = (-1) m0 [T -s n (w)] 2 µ n (w),
where µ n (w) is the coefficient corresponding to variable w of the Lagrangian function L, which satisfies L(x i ) = 0, for i = 1, . . . , n and L(w) = 1. We can see that µ n and g n are not defined at x 1 , . . . , x n and µ n (w) = ∞ if w ∈ {x 1 , . . . , x n }, so there might be computational difficulties if µ n is evaluated at any points closed to x 1 , . . . , x n . Therefore, instead of minimizing g n (w), we maximize a differentiable everywhere function h n (w) which is defined by

h n (w) = 1 gn(w) if w = x 1 , . . . , x n 0 else.
It is not too difficult to solve this auxiliary subproblem, but we need to be careful when selecting the target T since the problem is very sensitive to the choice of T . If T is large and close to y min , the search is highly local; when T gets smaller, the search will be more global. To handle it, we select T by performing cycles of iterations, in which each cycle starts with a low value (global search) and ends with a high value close to min s n (x) (local search). Then we go back to a global search, starting the cycle again [START_REF] Gutmann | A radial basis function method for global optimization[END_REF]. With certain choices of the targets T , the iterates from this algorithm are dense [START_REF] Gutmann | A radial basis function method for global optimization[END_REF]. rbfSolve [START_REF] Björkman | Global optimization of costly nonconvex functions using radial basis functions[END_REF] is a Matlab implementation of the RBF-Gutmann method that uses several strategies for selecting the targets T . The first strategy is to perform a cycle of length N + 1, starting at some n = ñ and for ñ ≤ n ≤ ñ + N -1, we choose

T = min w∈D s n (w) -W n max i∈Sn f (x i ) -min w∈D s n (w) , in which W n = N -n+ñ N 2
and S n is form by {x 1 , . . . , x n } after removing n -n max points with largest function values. Note that in a cycle, the values of W n are decreasing, therefore the corresponding target gradually get closer to the minimum value s * n of s n (x). The second strategy is to choose T as the optimal value of the following auxiliary problem:

min x∈D T (x) subject to: µ n (x) s n (x) -T (x) 2 ≤ α 2 n
and then perform a cycle of length N +1 on the choice of α n (in the paper, the authors use Cholesky factorization to get update for the next interpolations. It helps us to reduce computation time for solving linear system of equations). Another adaptive choice of the targets T is introduced in the ARBF method [START_REF] Holmström | An adaptive radial basis algorithm (arbf) for expensive black-box global optimization[END_REF], in which a large set of targets T j are used:

T j = s min n -β .W j . f △ ,
where β is an adaptive factor and the range f △ is selected as in [START_REF] Holmström | An adaptive radial basis algorithm (arbf) for expensive black-box global optimization[END_REF]. If for different target values, we obtain optimal solutions that are very close to each other, we need to increase β. On the other extreme, if the optimal solution regarding the second target is too far from the one obtained with the first target, we need to decrease β because it seems like the target values are too spread out.

V D n (x) = △ max n -△n(x) △ max n -△ min n if △ max n = △ min n 1
otherwise.

The MRS weighted-score merit function is then defined by:

W λ n (x) := λ V R n (x) + (1 -λ)V D n (x)
, where λ n ∈ [0, 1] is the weighted score specified by the user. The new candidate point will be found by minimizing W λ n . Note that when λ n = 1, the subproblem is equivalent to minimizing s n (x) (a local search). On the other hand, when λ n = 0, the subproblem is equivalent to maximizing △ n (x), i.e. to force x to stay away from all existing data points (a global search). Different choices of λ n will help us to manage the balance of local and global searches. Under mild assumptions, the MSRS method is shown to convergent almost surely. In [START_REF] Regis | A stochastic radial basis function method for the global optimization of expensive functions[END_REF], two specific implementations are discussed: Global Metric Stochastic Response Surface (GMSRS) and Local Metric Stochastic Response Surface (LMSRS). In GMSRS, the set Ω n of test points are uniformly generated throughout D. To achieve a balance between global and local search, the weight λ n is allowed to cycle through the iterations, starting from a high value (global search) and ending with a low value (local search). In LMSRS, the set of test points is generated by adding random perturbations to the current best solution x * n . These perturbations are chosen to be normally distributed with zero mean and with covariance matrix σ 2 n I d , in which σ n is the step size defined by the user. Global search can be done by restarting the algorithm whenever it appears to have converged to a local solution, i.e. the number of consecutively failed attempts is large enough.

An extension of the LMSRS method when some of the constraints are black-boxes is given in [START_REF] Regis | Stochastic radial basis function algorithms for large-scale optimization involving expensive black-box objective and constraint functions[END_REF]. In this case, the authors use multiple surrogate models to approximate the objective and expensive constraint functions. Starting with a design having at least one feasible point, we then perform the following loop: First, we construct surrogate models for all these black-box functions. Next, test points are generated randomly by perturbing the current best feasible point. These test points are then filtered so that they consist of only the elements that do not violate (or violate the minimum number of) surrogate-constraints. The next candidate point is selected by minimizing the MRS weighted-score function. During the process, we continuously update counters for consecutive failures and successes. The purpose is to check whether we are trapped at a local minimum or not, so we can appropriately adjust the step-size σ n (see LMSRS method) to switch back and forth between local and global search.

• Weighted Model Uncertainty. In the qualSolve algorithm, proposed by Jakobsson et al. [START_REF] Jakobsson | A method for simulation based optimization using radial basis functions[END_REF], the strategy for selecting new points for evaluation is to minimize the total model uncertainty weighted against the surrogate function value. Similar to [START_REF] Regis | Constrained global optimization of expensive black-box functions using radial basis functions[END_REF][START_REF] Regis | A stochastic radial basis function method for the global optimization of expensive functions[END_REF], the measure of uncertainty at point x is defined as

U Sn (x) := min xi∈Sn d(x, x i ).
Let ω(S(x)) be a weight function which gives a point x a low weight value when S(x) is high and gives x a high weight value when S(x) is low. Then we consider the following quality function

Q(y) := D (U Sn (x) -U Sn∪{y} (x))ω(S(x))dV (x).
The point x n+1 to evaluate next is obtained by solving

x n+1 ∈ arg max y∈D Q(y).
The weight function ω should satisfy the following demands -no area are completely disregarded, i.e. ω(z) > 0 for all z ∈ R.

areas with low surrogate function values are weighted higher than areas with high surrogate function values, i.e ω is strictly decreasing.

In [START_REF] Jakobsson | A method for simulation based optimization using radial basis functions[END_REF], the following weight function ω is proposed

ω(x) = γ exp - z -s min n s max n -s min n ,
where γ ≥ 1 is the parameter which controls the balance between local and global search.

Choices of search regions

Search regions are the domains where we solve auxiliary subproblems, i.e. to minimize or maximize the merit functions. In many cases, they are simply chosen to be the whole design space D. However, proper choices of search regions can also help us to manage the balance between local and global searches.

There are several ways to select a search region. For example, we can choose it as a neighborhood of some promising point (e.g. a minimizer of the surrogate function, the current best point). Normally, we use closed balls (with l 1 norm, for simplicity) to define a neighborhood. In [START_REF] Regis | Improved strategies for radial basis function methods for global optimization[END_REF], the authors realize that the original implementation of RBF-Gutmann algorithm can sometimes converge slowly to the global minimum on some test problems due to its failure to perform local search. This is because the candidate point for evaluation might be far from previous evaluated points and the global minimizers of s n (x). Another situation might cause the failure is when a global minimizer of s n (x) coincides or is too close to one of previously evaluated points. When the above cases occur, the strategy of choosing small W n in [START_REF] Gutmann | A radial basis function method for global optimization[END_REF] does not guarantee local search. Therefore, they suggest that, whenever W n is small, we restrict the global minimization of the bumpiness function to a small hyper-rectangle centered at a global minimizer of s n (x), i.e. we solve min

[x * n -rn,x * n +rn]∩D B n (x),
where x * n is a global minimizer of s n (x), B n (x) is the bumpiness function. Radius r n is an increasing function of the parameter W n , which tends to 0 when W n is small and vice versa. In particular, they use the following choice for r n

r n = v √ W n (b -a) if 0 ≤ W n ≤ U b -a otherwise,
where 0 < v < 1 and U > 0 are parameters to be specified and a, b are given from the definition of design space

D = {x ∈ R n | a ≤ x ≤ b}.
We can also define a search region of certain distance from existing data points. The idea is to push the next candidate point away from the points that have already been evaluated, thus enable us to move to unexplored regions. In this case, the search regions can be selected as

D n := D \ n i=1 B(x i , r i )
One of the examples, the CORS algorithm, proposed by Regis and Schoemaker [START_REF] Regis | Constrained global optimization of expensive black-box functions using radial basis functions[END_REF], uses these search regions. Let S n be the set of all existing data points, then the next candidate point x n+1 is chosen to minimizes the surrogate model over the set

x ∈ D, x -S n ≥ β n △ n ,
where β n is decided by the user and △ n is defined by

△ n = max x∈D x -S n .
The parameters β n are set by performing cycles of N + 1 iterations where each cycle employs a range of values for β n , starting with a high value close to 1 (global search) and ending with β n = 0 (local search). The iterates generated by the algorithm are dense therefore the method are able to find a global solution.

Solving black-box optimization problems

We know that by minimizing (or maximizing) a merit function over search regions, we will obtain next candidate points for evaluation. In the previous subsections, we have introduced various choices for such merit functions and search regions. Now we will discuss how to solve these auxiliary subproblems.

In principle, these auxiliary subproblems can be written explicitly as NLPs, therefore we can solve them by a generic solver. Moreover, we can exploit derivative information, which is available in most of cases, to make use of classical methods. However, due to the high nonlinearity and multi-modality of some merit functions (e.g. Expected Improvement), these subproblems are sometimes still difficult to solve. Moreover, it is questionable whether we should solve these problems to optimality or only approximate solutions are sufficient. Merit functions are constructed based on approximate models of the true function, which might be inaccurate. For this reason, it is risky if we spend too much time to solve the auxiliary subproblems, since in many cases, it is the quality of the model that matters. Therefore, approximate solutions (even solutions that are fairly good ) are sufficient.

Taking advantages of the fact that merit functions are given explicitly and very cheap to evaluate, we can use heuristics like genetic algorithms, simulated annealing and so on, to solve them. The simplest way is to sample a large number of points in the search region, and select the point with the best merit function values. This idea is implemented in [START_REF] Regis | A stochastic radial basis function method for the global optimization of expensive functions[END_REF], in which test points are randomly generated using either uniform and normal distributions and the MRS weighted score is used as the merit function. Under mild assumptions, the corresponding method converges to a global minimum almost surely, provided that the algorithm is allowed to run indefinitely.

In many real-life problems, the black-box function can be affected by noises. This makes the surrogate models far from accurate and it may lead to many difficulties for finding optimal solutions. Therefore, we have to adjust surrogate models to deal with the present of noises before going further. In [START_REF] Jakobsson | A method for simulation based optimization using radial basis functions[END_REF], a surrogate model that deviates from the data points is chosen by solving the following subproblem:

min µ S 2 + (1 -µ) e 2 2 s.t S(x i ) = f i + e i e i ∈ R n .
This problem is formed by adding a penalty term to the objective. In [START_REF] Jakobsson | A method for simulation based optimization using radial basis functions[END_REF], the authors show that the optimal error term e can be obtained by solving an explicit linear system:

µ -1 µ I n×n -B T AB e = B T AB f .
A similar approach to deal with noise is introduced in [START_REF] Costa | RBFOPT: an open-source library for black-box optimization with costly function evaluations[END_REF]. In this paper, the range within which function values are allowed to vary is required to specify. Assume that in addition to the black-box f (x) we have access to f (x) such that f (x) = f (x)(1+ ǫ r )+ ǫ a , where ǫ r , ǫ a are random variables with bounded support and unknown distribution. To determine the surrogate model, the authors introduce a vector of slack variables ξ ∈ R k and solve the problem: There are also cases when there is no noise in the black-box, but some points might have much larger function values comparing to others. In these cases, the surrogate model tends to oscillate strongly. It might also happen that min s n (y) is much lower than the best known function value, which leads to a choice of control parameters that overemphasizes global search. To handle these problems, it is suggested that we replace large function values at each iteration by the median of all computed function values [START_REF] Björkman | Global optimization of costly nonconvex functions using radial basis functions[END_REF] (see more discussion about it in [START_REF] Gutmann | A radial basis function method for global optimization[END_REF][START_REF] Ji | A new framework for combining global and local methods in black box optimization[END_REF]).

Since black-box functions are often expensive to evaluate, we can use parallel computation to speed up the algorithms. The idea of parallel computing is to break the main computation task into independent tasks so that each can be executed in a separate worker computer (or processing unit) simultaneously with the others. In black-box optimization, it is relatively easy to design a parallel algorithm. At each iteration, we find different candidate points (by solving multiple auxiliary subproblems) and assign each point to one worker computer for evaluation. In [START_REF] Regis | Improved strategies for radial basis function methods for global optimization[END_REF], parallel versions of RBF-Gutmann and CORS-RBF methods are implemented. While in RBF-Gutmann, the set of candidate points are obtained by solving subproblems for different values of the targets T n , in CORS-RBF, they are obtained by choosing different factors β n . All these subproblems are solved by the master computer, however we can modify the algorithm to empower worker computers to do that. It is easy to notice that parallel computation can be implemented for many other surrogate-based methods in the same way. For example, different auxiliary subproblems can be found by choosing different values for the radius r n (in RBF-CG method), the score λ n (in MSRS method), the control parameter γ (in qualSolve) and so on.

Dealing with discrete variables

In many applications, we need to deal with integer and binary variables, in addition to continuous ones. It is suggested that we can modify the surrogate-based methods to address these cases. In [START_REF] Müller | So-mi: A surrogate model algorithm for computationally expensive nonlinear mixed-integer black-box global optimization problems[END_REF], a method called SO-MI for solving black-box optimization problems with mixed integer variables is proposed. It starts with some feasible point and try to iteratively improve it. In order to do that, four potential new points are chosen for expensive black-box evaluations. They are selected from 4 different groups of points which are generated by perturbing the current solutions in different ways (either continuous or discrete variables or both) or by sampling at unexplored regions. By constructing RBF surrogate models, we determine the best point from each group with respect to the MRS weighted scores. A similar method called SO-I is developed in [START_REF] Müller | So-i: a surrogate model algorithm for expensive nonlinear integer programming problems including global optimization applications[END_REF] to deal with pure integer problem. Note that both of them start with a symmetric Latin hypercube design and round to the nearest integers. Moreover, they can handle blackbox constraints by incorporating them into the objective functions under penalty terms. The starting feasible point, in those cases, are found by iteratively solving some auxiliary problems.

Conclusions

In this paper, we present surrogate-based methods for solving black-box optimization problems. In these methods, the key to success is the ability to construct good approximations for the black-box functions. As long as we have a good approximation of the objective, in principle, it is not so difficult to find global minimal solutions. It is because approximate models are often quite simple and easy to work with. Moreover, even if they are analytically complicated, they are much faster to evaluate and we can apply heuristics such as genetic algorithms to solve. However, it is often impractical to expect that we can accurately approximate the black-box function everywhere in the domain: such an approximation requires us to evaluate of a huge number of data points, which is extremely expensive. In fact, for our purpose, we only need to find the models that approximate the black-box fairly well in the neighborhoods of (local or global) optima. Therefore, we should concentrate more on areas where optima are likely to locate. This observation is fundamental to all the black-box algorithms we introduced in the previous sections. A good model in this case needs not to be a good approximation everywhere but has to fit the black-box function well in the areas that are intensively evaluated. This is very intuitive since the more we explore in one area, the better we understand the function there. Kriging and radial basis functions are two such models that are very popular in the literature.

Since it is unknown where optima are located, we will look at the two types of potential areas: (i) neighborhoods of the best data points so far and (ii) areas that are least explored. Surrogate-based methods work by iteratively updating data points from one of these such areas. One idea to do that is to construct merit functions so that their optima are either close to the best current points or not explored (or both). Many such merit functions have been introduced in the previous sections.

Another idea is to properly manage the search regions, i.e. areas where we look for candidate points. Some authors have combined the both two ideas in their algorithms, but a systematic test of all combinations (models + merit functions + search regions) has not yet implemented. It might be a good idea to do such a test to compare all current black-box algorithms.

Another possible research direction is to devise parallel implementations for black-box optimization algorithms. This is already done in [START_REF] Regis | Parallel radial basis function methods for the global optimization of expensive functions[END_REF] for RBF with relatively simple ideas. Modifications of this method to apply for other merit functions are possible. Moreover, hopefully more complicated designs, in which the relations between candidate points found by different worker computers are taken into consideration, will speed up the algorithms.

x

  11 x 12 . . . x 1d x 21 x 22 . . . x 2d . . . . . . . . . . . . x n1 x n2 . . . x nd

Figure 1 :

 1 Figure 1: Latin square examples (reproduced from [Viana13]).

min (- 1 )

 1 dmin+1 λ T Φλ s.t Φλ + P h + ξ = F P T λ = 0 -ǫ r | f (x i )| -ǫ a ≤ ǫ i ≤ ǫ r | f (x i )| + ǫ a for all i ∈ L ǫ i = 0 for all i / ∈ L.

  construct a surrogate model s k (.) that approximates the black-box function.

-Phase 3 (search): Use s k (.) to search for a new point to evaluate.

Evaluate the new chosen point, update the data set S k . Assign k

•

  Metric Response Surface (MRS) weighted score. The Metric Response Surface (MRS) weighted score is a merit function proposed by Regis et al.[START_REF] Regis | A stochastic radial basis function method for the global optimization of expensive functions[END_REF]. In order to find the next candidate for evaluation, the authors do not minimize (or maximize) that merit function over a continuous set. Instead, they select the best among a finite set of N randomly generated points. Merit functions are made up from two criteria: the value of the surrogate model at each point and its minimum distance to existing data points. Explicitly, let Ω n be the set of N randomly generated points at the iteration n. Let s max n , △ max

	For each x ∈ Ω n , we identify two scores
	V R n (x) =	sn(x)-s min n s max n -s min n 1	if s max n otherwise = s min n
	and
	n s n (x) and △ n (x) over Ω n , where	(corr. s min n , △ min n ) be the maximum (corr. minimum) values of
		△ n (x) := min 1≤i≤n	d(x, x i ).

OVERVIEW OF EXPERIMENTAL DESIGNS

SURROGATE-BASED SUBPROBLEMS

Acknowledgements

The first author of this work is supported by a Microsoft Research PhD Fellowship. The second and fourth authors are partly supported by the FP7 Marie-Curie Initial Training Network "MINO". We thank an anonymous referee for many useful suggestion.