Mickaël Binois
email: mickael.binois@inria.fr

Nathan Wycoff

A survey on high-dimensional Gaussian process modeling with application to Bayesian optimization

Keywords: Black-box optimization, low-intrinsic dimensionality, additivity, variable selection, active subspace

Bayesian Optimization, the application of Bayesian function approximation to finding optima of expensive functions, has exploded in popularity in recent years. In particular, much attention has been paid to improving its efficiency on problems with many parameters to optimize. This attention has trickled down to the workhorse of high dimensional BO, high dimensional Gaussian process regression, which is also of independent interest. The great flexibility that the Gaussian process prior implies is a boon when modeling complicated, low dimensional surfaces but simply says too little when dimension grows too large. A variety of structural model assumptions have been tested to tame high dimensions, from variable selection and additive decomposition to low dimensional embeddings and beyond. Most of these approaches in turn require modifications of the acquisition function optimization strategy as well. Here we review the defining structural model assumptions and discuss the benefits and drawbacks of these approaches in practice.

Introduction

A large number of parameters (high dimensionality) is regularly mentioned as a critical challenge for black-box optimization. The root of the issue is the so-called curse of dimensionality, as coined by [START_REF] Bellman | Dynamic programming[END_REF]: the exponential dependence of complexity on input to BO. Optimization of the acquisition function to select the next design point benefits from (comparatively) quick evaluation and gradient evaluation, but still suffers from these same high-dimensional effects. Depending on the difficulty of the optimization problem, these issues can occur when reaching ten variables, or manifest for dozens of them. Still, dedicated methods have been shown to work for billions of variables [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF]though under rather limiting assumptions.

Indeed, stronger structural model assumptions are needed as the dimension increases, with three main categories. One idea is to reduce the dimension by removing variables with little or no impact on the output. Another is to define a few new variables based on linear or non-linear combinations of the original ones. The last direction is to assume additivity of effects of variables, or groups of them. With more structure to learn, model inference becomes harder, adding estimation risk to the difficulties above. Indisputably, estimation risk can be reduced by avoiding estimation, and this has been the approach taken by researchers relying on randomly defined structures (though obviously this is at the cost of variance and bias possibly gained from the randomization procedure). If we want to infer the structure, however, we also must face the fact that in the sequential design framework, only limited data is at first available to estimate characteristics of the function.

Either the structure can be learnt and updated "online" as the data are observed, or we can break the sequential process into two consecutive stages, the first of which is focused on estimating structure and the second on exploiting it. In this latter case, the optimal balance of budget to dedicate to each stage can be quite problem dependent.

The chosen structure affects the acquisition phase in several ways. It can be exploited for the modeling effort, as in reduced dimension or additivity assumptions, or it can constrain the optimization domain to aid the acquisition function search. Independently, strategies such as deploying trust regions as in [START_REF] Eriksson | Scalable Global Optimization via Local Bayesian Optimization[END_REF] are also available to limit the effect of the curse of dimensionality by limiting the volume of the search region. The resulting optimization approach is more local, but is complemented by restarts for globalization.

There exists many additional challenges and refinements for GPs and BO that are out of scope here, for which the corresponding techniques may need to be adapted for high dimension. One such issue is the scaling in terms of number of design points. Logically, more designs are needed to learn in larger dimensions, but one can not hope to match the exponential dependence on the dimension. Techniques for coping with large data have also attracted a lot of attention, see for instance [START_REF] Hensman | Gaussian processes for Big data[END_REF]; [START_REF] Heaton | A case study competition among methods for analyzing large spatial data[END_REF]; [START_REF] Wang | Exact Gaussian processes on a million data points[END_REF]; [START_REF] Kleijnen | Prediction for big data through Kriging: small sequential and one-shot designs[END_REF], where some are independent of the input dimension as is the case for, say, local models. In general we will assume here that running the black-box remains limiting compared to running the BO framework. Other refinements for GP and BO to cope with complex noise modeling (e.g., non Gaussian noise, input dependent variance) or non-stationarity could be adapted, but high-dimension exacerbates the difficulty of the learning task. While we focus on unconstrained optimization, batch (or parallel) optimization, constrained optimization and others could be entertained as well. We refer the interested reader to [START_REF] Garnett | Active learning of linear embeddings for Gaussian processes[END_REF] for a broader and more introductory overview of BO.

The remainder of this paper is as follows. First, key notions and notations are introduced in Section 2. Next, structural assumptions for high dimensional GP modeling are detailed in Section 3, before consequences and adaptations for the acquisition function optimization are presented in Section 4. Section 5 includes a list of possible test functions. Finally, some practical guidelines and a summary on promising research directions is given in Section 6.

Background

Let us consider an expensive-to-evaluate black-box simulator f : X ⊂ R d → R that we want to globally optimize:

find x * ∈ argmin x∈X f (x).
(1)

By black-box, we mean that nothing is assumed known about the functional form of f : f (x) can only be queried at any given input point x (sometimes gradients are also assumed available).

The Gaussian stochastic process

Given an index set X , in our case typically a closed and bounded subset of R d , a stochastic process is simply a rule for assigning to B members of that set x (1) , x (2) , . . . , x (B) a joint distribution of random variables; maybe it has a density δ y(x (1) ), y(x (2) ), . . . , y(x (B) ) with respect to some dominating measure. If the joint distribution is assumed to be Gaussian, we simply need a way of deciding what the mean vector and covariance matrix are going to be for any possible set of B points. These are naturally referred to as the prior mean function µ and covariance function k (also kernel function). Such a mathematical construct taken all together is called a Gaussian process (GP), see, e.g., [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF].

If we observe y(x * ) (possibly corrupted by Gaussian noise), we can imagine that there is some stochastic process that maps x (1) , x (2) , . . . , x (B) to δ y(x (1) ), y(x (2) ), . . . , y(x (B) )|y(x * ) , that is, that maps any set of points to its conditional distribution given y(x * ). Happily, this new stochastic process also happens to be a GP (conjugacy with respect to the normal likelihood), and its new mean and covariance functions are available in closed form. In particular, if we observed y(x (1) ), . . . , y(x (n) ), then:

m n (x) = µ + k(x) K -1 (y -µ1), s 2 n (x) = k(x, x) -k(x) K -1 k(x)
where y := (y 1 , . . . , y n ), k(x) := (k(x,

x (i) )) 1≤i≤n , K := (k(x (i) , x (j) ) + τ 2 1 i=j ) 1≤i,j≤n . τ 2 is
the noise hyperparameter, when assuming

y i = f (x (i) ) + ε i , with ε i = N (0, τ 2 ).
Often, the prior mean function µ is chosen to be zero (after centering observed y): the GP's focus is on defining and exploiting spatial covariance. However, it can definitely be useful to have different µ. If the objective function is known to be nonstationary, adding in polynomial trend terms can be a way to tame this. However, if the proper trend is not a priori known, it pays to be cautious when using one in high dimension as it's easy to accidentally extrapolate rather severely and perniciously [START_REF] Journel | Geostatistics for conditional simulation of ore bodies[END_REF]. Sparse trend selection might be preferable, as performed for instance by [START_REF] Kersaudy | A new surrogate modeling technique combining Kriging and polynomial chaos expansions-Application to uncertainty analysis in computational dosimetry[END_REF] in a combination with polynomial chaos expansion, as may be using basis elements informed via cross-validation [START_REF] Liang | Using Cross-Validation to Design Trend Function in Kriging Surrogate Modeling[END_REF] or a Bayesian framework [START_REF] Joseph | Blind Kriging: a new method for developing metamodels[END_REF].

The choice of covariance function family determines the qualitative properties of the GP. Of course, we want to make sure that when we're done applying the kernel function to all our data, the covariance matrix that comes out is positive semi-definite. Functions which guarantee this are themselves called positive semi-definite. Attention is often further restricted to stationary kernels, where the kernel function k is a function of x-x , k(x, x ) = k(x-x ). Within this class, the most popular kernel is probably the squared exponential one

in product form, k(x, x ) = σ 2 d i=1
exp -

(x i -x i ) 2 2θ i
. This kernel function is defined only up to its free parameters θ 1 , . . . , θ d , the lengthscales which determines at what distance covariance begins to drop off between points, and σ 2 which converts correlation to covariance based on the scale of y. This kernel induces an infinitely differentiable metamodel, which may or may not be desirable, depending on how well that matches our understanding of the black-box. For precise control over differentiability1 , consider the Matérn class indexed by smoothness parameter ν, which is defined in terms of the modified Bessel function of the second kind (though can be written in terms of simple functions when ν + 0.5 is a whole number). See Rasmussen and Williams (2006, Chapter 4) for more on the Matérn and other covariance function classes. Such GPs can be differentiated at most ν -0.5 times.

This dependence of GPs on distances between points will be problematic in high dimension. And, consequently, many of the approaches we review in this article will define a kernel function not on the original input space, but on some transformation thereof. Of course, it is possible to view the transformation and the kernel function thereon together as a novel kernel function. But we will instead take the perspective in this article of separating dimension-reducing transformations and the subsequent purely distance-based kernel, as we view this modularity as helpful in model design and comparison.

Insofar as a Gaussian process is the infinite dimensional generalization of the multivariate Gaussian distribution, it retains many of the analytical properties that have propelled Gaussianity to its prominence. Of particular use to us will be the closed form evaluation of certain exceedance probabilities and related metrics which form the basis of deploying these Gaussian processes to optimization problems.

We refer to [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]; [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences[END_REF] for additional details on GP regression, and, e.g., to [START_REF] Kanagawa | Gaussian processes and kernel methods: A review on connections and equivalences[END_REF] for a discussion of connections to other kernel methods.

Bayesian optimization

Popularized with the Efficient Global Optimization method of [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF], Bayesian optimization, see, e.g., [START_REF] Maus | Local Latent Space Bayesian Optimization over Structured Inputs[END_REF], relies on the GP probabilistic prior of f based on a few initial observations to define an acquisition function α : X → R. This is later optimized to select new points to evaluate sequentially. The initial stage, known as the Design of Experiments (DoE), typically starts by evaluating designs based on optimized Latin hypercube samples, for their space filling properties, see, e.g., [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences[END_REF]. A pseudo-code is given in Algorithm 1.

Algorithm 1 Pseudo-code for BO Require: N max (total budget), GP model trained on initial DoE of size n = n 0 : (x (i) , y i ) 1≤i≤n 1: while n ≤ N max do 2:

Choose x (n+1) ∈ arg max x∈X α(x)

3:

Update the GP model by conditioning on (x (n+1) , y n+1 )

4: n ← n + 1 5: end while
The expected improvement (Mockus et al., 1978, EI) and upper confidence bound (Srinivas et al., 2010, UCB) criteria are typically used due to their analytical tractability in both value and gradient. Other choices are possible, see e.g., Shahriari et al. (2016b). The generic acquisition function optimization problem is:

find x * ∈ argmax x∈X α(x).
(2) Some acquisition functions are natively suitable for when the black-box is noisy, while some need adaptation like EI, see e.g., [START_REF] Letham | Constrained Bayesian optimization with noisy experiments[END_REF] for a discussion. A pinch of noise might be beneficial for its regularization effect, see, e.g., [START_REF] Gramacy | Cases for the nugget in modeling computer experiments[END_REF], but the effort of finding a solution increases drastically as the signal to noise ratio sinks.

Not to mention input-varying and non-Gaussian noises that are hard to deal with even in low dimension. We refer to [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF]; [START_REF] Roustant | DiceKriging, DiceOptim: Two R Packages for the Analysis of Computer Experiments by Kriging-Based Metamodeling and Optimization[END_REF]; [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences[END_REF]; [START_REF] Garnett | Active learning of linear embeddings for Gaussian processes[END_REF] for more details on how to apply BO in practice.

The curse of dimensionality and its effects

An increase in input dimension is first felt in distance computation, which is at the core of most common covariance kernels. The difficulty is that points are relatively further away in high dimension, such that it becomes harder to learn with a covariance based on distances. As is illustrated in Fig. 1, for large d, the squared distance between uniformly sampled points concentrates increasingly further away. In particular, the expected squared distance approaches to d 6 with standard deviation relative to the width of the interval of possible distances [0, √ d] decreasing as 7 180d

2 . See [START_REF] Köppen | The Curse of Dimensionality[END_REF] for a discussion of other implications and properties of high dimensional cubes. Except for isotropic kernels, the second effect is the increase of the number of hyperparameters of the models, typically with d many lengthscale parameters to tune for product covariance kernels. Maximizing the (non-convex) likelihood thus becomes increasingly hard, since the curse of dimensionality also affects the gradient based optimization routines typically used. Grouping variables to share the same lengthscales alleviate this, which has been done either manually based on problem-specific considerations as done by [START_REF] Binois | Uncertainty quantification on Pareto fronts and high-dimensional strategies in Bayesian optimization, with applications in multi-objective automotive design[END_REF] or using the Bayesian information criterion (BIC) as done by [START_REF] Blanchet-Scalliet | A specific kriging kernel for dimensionality reduction: Isotropic by group kernel[END_REF].

Similarly, accurately optimizing the acquisition function is complex in high dimension, especially since it is a multi-modal optimization task, with flat regions. Furthermore, given that most of the volume is on the boundary of the domain and that GP-based predictive variance increases with the distance to the designs -a desirable property at least in low dimension -the downside is that the optimum is generally found on one of the exponentially many vertices or sides of the d-hypercube, tilting the balance towards blind exploration.

The curse of dimensionality also manifests itself when random sampling, required for Monte Carlo estimation of certain integral quantities, among them entropy, the optimizing point of a posterior realization (as is required in Thompson sampling), and certain global sensitivity measures. The issue is that a uniformly sampled point will be increasingly distant, on average, to any given point as dimension increases. Hence the challenge is not only to scale in terms of modeling accuracy but also to keep inference manageable and avoid pitfalls in acquisition function optimization.

3 High-dimensional Gaussian process modeling

Structural model assumptions are the only way to avoid the exponential dependence on the dimension in modeling, for which various options are summarized for instance by [START_REF] Bach | Breaking the curse of dimensionality with convex neural networks[END_REF], who also gives the corresponding generalization bounds. These different structural model assumptions have been adapted in the GP context, plus others, as we detail next.

One exception is to use isotropic kernels and thus always a single lengthscale, for which increasing the dimension influences only the behavior of distances.

Variable selection or screening

Whenever possible, limiting the number of variables based on expert knowledge on the optimization problem is recommended. If expert knowledge is unavailable, a simple, datadriven idea is to perform variable selection, or screening, before optimization, e.g., using the Morris technique [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF] or via hierarchical diagonal sampling as in [START_REF] Chen | Joint optimization and variable selection of high-dimensional Gaussian processes[END_REF]. Other global sensitivity analysis techniques can be used to select variables individually and we refer to [START_REF] Iooss | A review on global sensitivity analysis methods[END_REF] for an entry point to this topic.

Hence, one of the early attempts to tackle high-dimensions is to assume that most of the variables have no effect:

model: f (x) ≈ g(x I ) with I ⊂ {1, . . . , d}, |I| d (3)
and then identify those influential variables in the set I. For the Gaussian and other stationary kernels in tensor product form, this can be performed by looking at the lengthscale values: the covariance varies less for very large values of the lengthscale (when parameterized as above), whose corresponding input variables can thus be removed. This is also known under the term automatic relevance determination (ARD), as discussed, e.g., by Rasmussen and Williams (2006, Chapter 5). [START_REF] Salem | Sequential dimension reduction for learning features of expensive black-box functions[END_REF] shows that this indeed holds asymptotically. Hence the idea is to rank the variables based on their lengthscales to determine the more influential from the less. These influential variables are then used to build a GP model and optimize expected improvement.

Nevertheless, θ 1 > θ 2 does not necessarily imply that x 1 is less important than x 2 , see for example [START_REF] Wycoff | Sequential Learning of Active Subspaces[END_REF]; [START_REF] Lin | Transformation and additivity in Gaussian processes[END_REF]. Hence rather than looking only only at the lengthscales values, [START_REF] Linkletter | Variable selection for Gaussian process models in computer experiments[END_REF] propose to compare the posterior distributions of lengthscales corresponding to both real and artificially added inert variables in a fully Bayesian GP framework. This is complemented by a local variable selection in [START_REF] Winkel | Sequential Optimization in Locally Important Dimensions[END_REF], relying of predictions made when ignoring some dimensions. [START_REF] Eriksson | High-Dimensional Bayesian Optimization with Sparse Axis-Aligned Subspaces[END_REF] choose instead to put a Horseshoe prior [START_REF] Carvalho | Handling sparsity via the horseshoe[END_REF] on the inverse lengthscales and perform gradient-based numerical posterior simulation in order to sample from this posterior of functions defined on subsets of variables. With a more incremental approach, [START_REF] Marrel | An efficient methodology for modeling complex computer codes with Gaussian processes[END_REF] propose to sequentially add variables to the regression (resp. covariance) element based on the Akaike information criterion (resp.

predictivity coefficient Q 2 ).
Rather than definitively deciding on which variables to keep, which can be timeconsuming, [START_REF] Ulmasov | Bayesian optimization with dimension scheduling: Application to biological systems[END_REF] propose to sample a few variables at each iteration,

where the weight vector is determined with principal component analysis (PCA) on the (x) 1≤i≤n . This may not be efficient for small budget as indicated in [START_REF] Li | High Dimensional Bayesian Optimization using Dropout[END_REF], who prefer to select variables uniformly and to fill in values for the remaining variables.

In practice, most variables typically have a possibly limited but non-null influence on the output. Though this can mean that choosing the number of variables to keep is somewhat arbitrary, the simplicity, both in terms of implementation and interpretation, of variable selection makes it a compelling approach to dimension reduction, particularly for a first pass. Of course, the variable selection itself becomes only harder as the dimension increases, but this is to a lesser extent than some of the more sophisticated approaches discussed below. Nevertheless, the approach's strength is also its limitation, and it is common for the complicated functions which make up black-box problems to truly and strongly depend on all input parameters. In the case where all variables have the same influence, variable selection would even fail to reduce the dimension by one no matter how simple the relationship. Other structural assumptions can overcome these limitations, as well as exploit interactions between variables.

Additive and ANOVA models

One set of structural assumptions amenable to keeping all variables but limiting their interaction is that of additivity:

model: f (x) ≈ µ + d i=1 g i (x i ) (4)
with univariate functions g i . This has been transposed to the GP framework by [START_REF] Neal | Monte Carlo implementation of Gaussian process models for Bayesian regression and classification[END_REF]; [START_REF] Plate | Accuracy versus interpretability in flexible modeling: Implementing a tradeoff using Gaussian process models[END_REF]; [START_REF] Durrande | Additive Kernels for Gaussian Process Modeling[END_REF]; [START_REF] Duvenaud | Additive Gaussian processes[END_REF], initially via the summation of univariate kernels:

k(x, x ) = d i=1 σ 2 i k i (x i , x i ),
which results in a valid covariance function just as the product form. A useful property for interpretability and visualization is that the GP predictive mean can be decomposed into a sum of univariate components:

m n (x) = k(x) K -1 y = d i=1 k i (x i )K -1 y = d i=1 m n,i (x i ) with k i (x i ) := (k i (x i , x (j) i )) 1≤j≤n .
A more surprising property is that the covariance can become non invertible due to linear relationships appearing between distinct observations -see e.g., [START_REF] Durrande | Additive Kernels for Gaussian Process Modeling[END_REF] for examples-which can be alleviated by adding a noise term. Consequently, the predictive variance can be zero at an unobserved design point, a somewhat detrimental side effect for exploration and thus optimization. Another difficulty is the estimation of the hyperparameters, with the need of an additional variance parameter σ i per coordinate (thus approximately doubling the number of hyperparameters). Still, compared to the tensor product form whose values quickly go to zero in high dimension, the sum form scales much better. Making the black-box more additive by applying a transformation of the outputs is explored by [START_REF] Lin | Transformation and additivity in Gaussian processes[END_REF].

Higher-order models can be defined in the same way [START_REF] Duvenaud | Additive Gaussian processes[END_REF], usually restricted to order two or selecting higher order components only if all lower orders are already selected. [START_REF] Plate | Accuracy versus interpretability in flexible modeling: Implementing a tradeoff using Gaussian process models[END_REF] instead proposes adding all interactions at once on top of first order ones. Directly identifying groups of variables is also possible, see e.g., [START_REF] Kandasamy | High dimensional Bayesian optimisation and bandits via additive models[END_REF]; [START_REF] Gardner | Discovering and Exploiting Additive Structure for Bayesian Optimization[END_REF]; [START_REF] Wang | Batched High-dimensional Bayesian Optimization via Structural Kernel Learning[END_REF][START_REF] Wang | Batched Large-scale Bayesian Optimization in High-dimensional Spaces[END_REF]:

model: f (x) ≈ µ + M i=1 g i (x A i ) (5)
with multivariate g i functions acting on A i disjoint subsets of {1, . . . , d} such that M i A i = {1, . . . , d}. The restriction of disjoint subsets of variables have been further lifted in subsequent works, see e.g., [START_REF] Rolland | High-Dimensional Bayesian Optimization via Additive Models with Overlapping Groups[END_REF]; [START_REF] Hoang | Decentralized highdimensional Bayesian optimization with factor graphs[END_REF].

Similar in their form but rooted in global sensitivity analysis are additive models based on the functional analysis of variance (fANOVA, a.k.a. Sobol-Hoeffding) decomposition [START_REF] Efron | The jackknife estimate of variance[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF]:

model: f (x) ≈ c + d i=1 g i (x i ) + j<k g jk (x j , x k ) + • • • + g 12...d (x 1 , x 2 , . . . , x d ) (6)
with elementary functions g ... that are required to be centered and orthogonal for the uniqueness of the decomposition. The benefit is that sensitivity analysis can be conducted and interpreted as in a regular ANOVA. [START_REF] Muehlenstaedt | Data-driven Kriging models based on FANOVA-decomposition[END_REF] rely on this decomposition up to second order interaction to build their model, choosing their components on the basis of Sobol indices that appear in this formulation (and estimated by a first-pass GP with an anisotropic tensor-product kernel). These are the so-called main effects (one variable only) and total interaction effects (effect of two variables at any order). Selecting which interaction to remove to form cliques requires a thresholding scheme. [START_REF] Ulaganathan | High dimensional Kriging metamodelling utilising gradient information[END_REF] proposes a similar approach with the addition of cut-points when gradient observations are available. [START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis[END_REF] go further with a dedicated kernel whose The main downside of these techniques is that inference is challenging in this context, with combinatorially many terms. Hence various techniques to estimate the hyperparameters have been applied: coordinate-ascent-like [START_REF] Durrande | Additive Kernels for Gaussian Process Modeling[END_REF] or quasi-Newton [START_REF] Duvenaud | Additive Gaussian processes[END_REF] methods for likelihood maximization. Relying on randomness is sometimes preferred to bypass the cost of a full optimization, like by [START_REF] Kandasamy | High dimensional Bayesian optimisation and bandits via additive models[END_REF] where random decompositions are sampled and the best one for the likelihood is selected (fixing the order and number of terms), or by [START_REF] Wang | Batched Large-scale Bayesian Optimization in High-dimensional Spaces[END_REF]. [START_REF] Journel | Geostatistics for conditional simulation of ore bodies[END_REF].

form is k ANOVA (x, x ) = D i=1 (1 + k i (x i , x i )),
In terms of advantages, on the other hand, this method maintains the interpretability of variable selection, particularly if the selected model includes mostly first or second order terms and is somewhat sparse. Learning which variables interact by observing which pairs are selected by an inference-driven procedure can be scientifically interesting in and of itself.

Some works are dedicated to scaling to many observations as well, such as Mutny and Krause ( 2018) with a basis expansion of GP kernels. [START_REF] Sung | Multiresolution functional anova for large-scale, many-input computer experiments[END_REF] also propose a basis expansion, further complementing it with a multi-resolution scheme with a group lasso estimation procedure. [START_REF] Wang | Batched Large-scale Bayesian Optimization in High-dimensional Spaces[END_REF] scale to many observations with a random partition of the input space, random additive approximation and random features decomposition of the kernel.

Here, one underlying hypothesis is that high order interaction components are negligible (because they are hard to estimate in this formalism), though this is a long-standing tradition in experimental design. Another concern is that ANOVA cannot detect non-linearity and multi-modality [START_REF] Palar | Exploiting active subspaces in global optimization: how complex is your problem?[END_REF]. These downsides are mitigated by the following framework.

Linear embeddings

One approach to designing high dimensional kernels is to avoid doing so directly, and instead run the data through a dimension reducing map first. The image of these data under this map is called its embedding, and an obvious class of functions to use for dimension reduction is the class of linear functions. Indeed, as noted by [START_REF] Marcy | Bayesian Gaussian Process Models for Dimension Reduction Uncertainties[END_REF], this idea of using a linear mapping dates back to at least [START_REF] Matérn | Spatial variation. Stochastic models and their application to some problems in forest surveys and other sampling investigations[END_REF]. The question is then simply which linear function to use in particular.

Denote this generic mapping by z = Ax with A ∈ R r×d . In the case of [START_REF] Vivarelli | Discovering hidden features with Gaussian processes regression[END_REF], r = d, so the mapping serves to rotate the space rather than actually reduce the dimension, but we are more interested in the r d case. When r = 1, this is a popular dimension reduction technique called the single index model:

model: f (x) ≈ g(a x) with a ∈ R d . (7)
We refer to [START_REF] Gramacy | Gaussian process single-index models as emulators for computer experiments[END_REF] for the GP treatment, with a treated as an additional kernel hyperparameter. a could instead simply be randomly chosen at each iteration, as suggested by [START_REF] Kirschner | Adaptive and safe Bayesian optimization in high dimensions via one-dimensional subspaces[END_REF]. Unfortunately, extending to r > 1 is by no means trivial.

The assumption is:

model: f (x) ≈ g(A x) (8)
where such functions are called ridge functions in the literature. It corresponds to the observation, sometimes backed by theoretical evidence [START_REF] Constantine | Many physical laws are ridge functions[END_REF], that the variation of high dimensional functions can be concentrated around a few but unknown directions. There are several approaches to choosing A, which we'll review in turn.

For GP regression, arguably the most direct way is simply to treat A as just another kernel hyperparameter to learn, e.g., by marginal likelihood optimization. Fixing r, [START_REF] Garnett | Active learning of linear embeddings for Gaussian processes[END_REF] provide an approximately Bayesian scheme to do so, with a Laplace approximation to the likelihood followed by an approximate marginalization over hyperparameters. [START_REF] Tripathy | Gaussian processes with built-in dimensionality reduction: Applications to high-dimensional uncertainty propagation[END_REF] rely on a two stage approach: first using the likelihood to learn an orthogonal A, then finding the rest of the hyperparameters (and repeat). BIC is used to identify r. The orthogonality constraint requires some special consideration, since these matrices lie on the Stiefel manifold.

:= X ∇(f (x)) ∇(f (x))λ(dx)
where λ is any well-behaved measure on the design space (typically the Lebesgue one on hypercubic domains). Then a Monte Carlo estimator can be used. Loosening the assumption that some eigenvalues will be exactly zero leads to the active subspace (AS) methodology [START_REF] Constantine | Active subspaces: Emerging ideas for dimension reduction in parameter studies[END_REF], analyzed for dimension reduction and visualization. In particular, the presence of an AS is hinted at by gaps in the eigenvalues of C [START_REF] Constantine | Active subspaces: Emerging ideas for dimension reduction in parameter studies[END_REF]. If the gradient can be evaluated, this is a good pre-processing step before GP modeling [START_REF] Eriksson | Scaling Gaussian Process Regression with Derivatives[END_REF]. Without the gradient information, finite differencing is generally too costly, such that using a GP to estimate C might be more appropriate, e.g., as in [START_REF] Fukumizu | Gradient-based kernel dimension reduction for regression[END_REF]; [START_REF] Palar | Exploiting active subspaces in global optimization: how complex is your problem?[END_REF]. [START_REF] Djolonga | High-Dimensional Gaussian Process Bandits[END_REF] uses directional derivatives with finite difference to recover A with low rank matrix recovery before BO. Other compressed sensing techniques can be used, as done initially by Carpentier and Munos (2012) and later by [START_REF] Groves | Efficient and scalable batch Bayesian optimization using K-means[END_REF]. The MC estimators of A are applied in a two stages approach, due to the assumption of iid sampling of the design points. [START_REF] Wycoff | Sequential Learning of Active Subspaces[END_REF] showed that an estimate of the C matrix of a GP is tractable for standard stationary kernels, relieving the need for this sampling assumption. Other works in the vein of sensitivity analysis on slice inverse regression or partial least squares to recover A and possibly update it during optimization [START_REF] Bouhlel | Improving kriging surrogates of high-dimensional design models by Partial Least Squares dimension reduction[END_REF][START_REF] Zhang | High dimensional Bayesian optimization via supervised dimension reduction[END_REF]Chen et al., 2020a). [START_REF] Lee | Modified Active Subspaces Using the Average of Gradients[END_REF] advocate for a modified AS matrix which exaggerates the influence of the average gradient (over the input space).

And finally, rather than learning A as part of inference on GP hyperparameters or as the result of a sensitivity analysis, a third option is to simply select A randomly, either by generating a single A matrix before the modeling step or using different As, say, at each iteration. For example, in random embedding BO (REMBO), [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF][START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF]) use a fixed and randomly sampled A. One justification is the stability of random projection of the L 2 norm due to the Johnson-Lindenstrauss lemma [START_REF] Johnson | Extensions of Lipschitz mappings into a Hilbert space 26[END_REF], as noticed, e.g., by [START_REF] Letham | Re-Examining Linear Embeddings for High-Dimensional Bayesian Optimization[END_REF]. In the context of BO, the rationale is further that, at least for an unbounded domain, there exists a solution to the problem on the low dimensional embedding, explaining the success of random optimization [START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF] on some hyperparameter tuning problems.

Even with a candidate A in hand, whether through optimization, sensitivity analysis, or random sampling, we still have modeling decisions to make (not to mention acquisition decisions; Section 4.3). The seemingly innocuous choice of a hypercubic domain can give headaches when combined with a linear dimension reduction, as the possible solution space is no longer simply a cube but a polytope defined by A. If our GP is fit to the low dimensional space, it may have an easier time finding a next candidate point to optimize, but we have to figure out which point in the original space corresponds to this low-dimensional optimum. This means that not only are preimages non-unique, but the backprojection of a point is not even guaranteed to live in the original unit hypercube, and a convex projection is possibly required to regain feasibility. This means that our backprojection is no longer linear. The remediation options proposed in [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF] are to ignore the corresponding non-injectivity issue, or to use a kernel defined on X , losing the benefits of low dimensional GP modeling. [START_REF] Binois | A Warped Kernel Improving Robustness in Bayesian Optimization Via Random Embeddings[END_REF] proposed to include high-dimensional information in Y with a warping to address the non-injectivity issue, while Binois et al. Though in general any GP can be fit in the reduced space, links to additive assumptions are natural, and [START_REF] Gilboa | Scaling multidimensional Gaussian processes using projected additive approximations[END_REF] combine these in the projection-pursuit style, leading to projected additive approximations. In this context, tuning the r parameter is performed by adding dimensions iteratively and halting when the accuracy stops increasing. [START_REF] Li | High Dimensional Bayesian Optimization via Restricted Projection Pursuit Models[END_REF] further discuss how to extend this approach to optimization, in a restricted-projection pursuit setup. This can be viewed as somewhere in between fitting a full low dimensional model and using first order additive model. These ideas can be combined by estimating several low-dimensional subspaces, see e.g., Yenicelik (2020); [START_REF] Wong | Embedded ridge approximations[END_REF].

The estimation risk is strong when using a linear embedding (which has p×r parameters in general), with unresolved questions: i) How well must the matrix A be learned before it is better than directly fitting the high dimensional problem? ii) How much of the budget should be dedicated to this task? iii) When is it better to learn A dynamically, leading to a noisy GP (whenever design points are not all on the same embedding as for REMBO)? Plus finding an appropriate r remains difficult unless expert information is given. Fortunately, taking larger values than necessary is not detrimental -except that it quickly reintroduces high-dimensional challenges [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF][START_REF] Cartis | A dimensionality reduction technique for unconstrained global optimization of functions with low effective dimensionality[END_REF].

Already, the linear embeddings approach jettison much of the interpretability that selection and ANOVA approaches offer (it will get worse when we consider nonlinear embeddings in the next section). In the case when r ∈ {1, 2}, it is possible to create visualizations of the function, which can be great sources of intuition. But in larger dimensional cases, we are reduced to squinting at the loadings of each variable in the retained directions. On the other hand, there are indeed functions which can be reduced to much smaller dimension when using linear embeddings compared to the axis-aligned variable selection3 . Which of these two approaches will be more appropriate is quite problem-dependent.

Non-linear embeddings and structured spaces

Dropping the linearity assumption adds even more flexibility to the model at the cost of requiring even more data to fit it. Recovering such a manifold suitable for regression has been proposed e.g., by [START_REF] Guhaniyogi | Compressed Gaussian process for manifold regression[END_REF], when the data is on a low dimensional manifold. A simple way to extend linear dimension reduction to the nonlinear case is to use a locally linear approach, as does [START_REF] Wycoff | Gradient-Based Sensitivity Analysis with Kernels[END_REF] GPs. In this line of thought, it is also possible to include generative topographic mapping [START_REF] Viswanath | Dimension reduction for aerodynamic design optimization[END_REF], GP latent variable models, e.g., [START_REF] Lawrence | Probabilistic non-linear principal component analysis with Gaussian process latent variable models[END_REF][START_REF] Titsias | Bayesian Gaussian process latent variable model[END_REF], or deep GPs [START_REF] Damianou | Deep Gaussian processes[END_REF][START_REF] Hebbal | Multiobjective optimization using Deep Gaussian Processes: Application to aerospace vehicle design[END_REF][START_REF] Sauer | Active Learning for Deep Gaussian Process Surrogates[END_REF]. Deep Gaussian processes refer to the modeling strategy of assuming that the input is linked to the output via a chain of Gaussian processes, the output of one serving as the input for the next. This is not to be confused with the concept of deep kernels [START_REF] Wilson | Deep kernel learning[END_REF][START_REF] Huang | Scalable Gaussian process regression using deep neural networks[END_REF], which involves only a single Gaussian process, but whose kernel function is parameterized by a neural network. Typically, inference is conducted on the GP and neural net weights via optimization of (an estimate of) the log marginal likelihood, but in the event of structured inputs or partial labeling the neural network can be initialized as an autoencoder (i.e., initialized to reconstruct the input). An orthogonal direction to avoid the larger optimization budgets is via multi-fidelity, when cheaper but less accurate version(s) of the black-box are available, as exploited e.g., by [START_REF] Ginsbourger | Distance-based kriging relying on proxy simulations for inverse conditioning[END_REF]; [START_REF] Falkner | BOHB: Robust and efficient hyperparameter optimization at scale[END_REF].

A setup where these highly flexible models may be more amenable is in the case where some additional information is available about the structure of the problem. One such case is the optimization of geometric shapes, typically airfoils. AS is performed in this context by [START_REF] Lukaczyk | Active subspaces for shape optimization[END_REF]. As many options to parameterize these shapes are available, some are more adapted to optimization. Independently of this choice, [START_REF] Gaudrie | Modeling and optimization with Gaussian processes in reduced eigenbases[END_REF] propose to work in the shape space, defined with shape eigenvectors and their values, which is not costly since the computation of the geometry is in general negligible next to the subsequent simulation. A non-linear computation of the shape basis is used in (2016), that do not scale well). The underlying assumption is that the solution is close to the center of the domain, if suitably chosen. Transforming the coordinates amounts to separating radial and angular components:

T (x) =    (||x|| 2 , x/||x|| 2 ) for ||x|| 2 = 0
(0, a arbitrary ) for ||x|| 2 = 0 , while the inverse transformation is T -1 (r, a) = ra. The corresponding covariance kernel is et al., 2014) with a user defined P .

k cyl (x, x ) = k r (r, r ) × k a (a, a ). The 1D radius component k r (r, r ) = k(|(1 -r α ) β -(1 - r α ) β |), α, β > 0 is
The nonlinear approaches have seen great success in high dimension, particularly when the space is structured/non-Euclidean. However, this comes at the cost of added complexity, both in computation (e.g., when using deep neural networks) and interpretation.

Therefore, our recommendation is to reserve these techniques for the toughest problems, and consider simpler solutions to simpler ones, particularly if insight into function, and not only an optimal solution, is desired.

There are many hoops and hurdles to construct a satisfying GP model in high dimension.

Yet, this only gets one half the way, since it remains to optimize the acquisition function.

High-dimensional acquisition function optimization

A proposed high-dimensional GP model is usually combined with a strategy for acquisition function optimization, adapted to the specific model structure and the curse of dimensionality. Many of these approaches can be extricated from the modeling framework in which they were proposed, leaving many possible combinations unexplored. We now discuss them on their own merits.

Faster evaluation time and the availability of gradients provide limited help towards globally optimizing the multi-modal acquisition function. Global optimality guarantees are out of reach, since branch-and-bound methods or DIRECT do not scale [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]. Aside from local minima, another pitfall is the presence of large plateaus where the acquisition function value is flat, while local optima may be peaked. To decrease these two effects, [START_REF] Rana | High dimensional Bayesian optimization with elastic Gaussian process[END_REF] Regardless of the quality of the fitted GP model, the expanding boundary volume as a function of the number of variables has a large impact on the acquisition's optimization (Equation 2). If based on uniform sampling, e.g., with multi-start gradient optimization or even evolutionary algorithms, the search will focus in these boundary areas. This is reinforced by the larger variance on the boundary that would locate the acquisition function's optimum there, with no hope of evaluating all the faces or vertices before focusing on the interior. Commenting on and discussing these effects, [START_REF] Oh | BOCK: Bayesian Optimization with Cylindrical Kernels[END_REF] proposed the use of cylindrical coordinates with BOCK to upweight the volume of the interior, imparting the prior knowledge that the optimum is close to the center. Adding virtual derivative observations on the boundary like in [START_REF] Siivola | Correcting boundary over-exploration deficiencies in Bayesian optimization with virtual derivative sign observations[END_REF] is hardly feasible in high dimension, as derivative information is only analytically enforceable at finitely many collocation points in GPs (and indeed computation scales cubically in the same). An infinite version could possibly be entertained via spectral methods, e.g., based on [START_REF] Gauthier | Spectral approach for kernel-based interpolation[END_REF]. In such a case, the use of a trust-region (TR) can limit the size of the search space drastically. Trust region methods focus the optimization locally within the neighborhood of a current best solution (TR center), whose size increases if the new candidate point improves sufficiently over the TR center, or decreases instead, see e.g., [START_REF] Larson | Derivative-free optimization methods[END_REF] for a review. Combined with BO, this has been shown to be quite beneficial in high-dimension [START_REF] Regis | Trust regions in Kriging-based optimization with expected improvement[END_REF][START_REF] Eriksson | Scalable Global Optimization via Local Bayesian Optimization[END_REF][START_REF] Diouane | TREGO: a Trust-Region Framework for Efficient Global Optimization[END_REF]Zhou et al., 2021;[START_REF] Daulton | Multi-Objective Bayesian Optimization over High-Dimensional Search Spaces[END_REF], perhaps at the cost of a less global search (possibly compensated for with restarts or parallel TR). To further avoid the attraction of the boundary of the TR, ( 2) is only optimized over a discrete set in [START_REF] Eriksson | Scaling Gaussian Process Regression with Derivatives[END_REF], with some coordinates randomly kept at the TR center.

While the strategies above are mostly independent of the GP model (e.g., used with isotropic or anisotropic product kernels), we next detail strategies for GPs with structural assumptions. Note that these assumptions, never perfectly fulfilled or estimated in practice, require the introduction of a noise component to account for the introduced approximations, with a few exceptions. Hence this has to be taken into account in the optimization of (2).

Additive case

Additivity provides the opportunity to rely on the same decomposition as the GP for solving (2), hence reducing the search to several lower-dimensional searches, possibly in parallel. That is, rather than focusing on the posterior of f , the idea is to look at those of the additive components instead, the g i from models (4, 5, 6):

N (m n,i (x I ), s 2 n,i (x I )) where s 2 n,i (x I ) = k i (x I , x I ) -k i (x i ) K -1 k i (x i
) for a general index I. Then partial acquisition functions are defined on each g i model, as a sum, e.g., as in [START_REF] Kandasamy | High dimensional Bayesian optimisation and bandits via additive models[END_REF]. More care is needed for overlapping subsets, say relying on message passing [START_REF] Rolland | High-Dimensional Bayesian Optimization via Additive Models with Overlapping Groups[END_REF][START_REF] Hoang | Decentralized highdimensional Bayesian optimization with factor graphs[END_REF] but it remains more efficient than optimizing in the original space.

The major advantage of this approach is that it allows for solving d many one dimensional optimization problems instead of one d dimensional problem, which ceteris paribus is many times easier in the nonconvex case. Nevertheless, the extent to which the zero variance at unobserved locations with additive models affects the globality of the search remains unknown. Adding a kernel a component that helps mitigate incorrect assumptions as for variable selection or AS as below may be interesting. Additive GPs should be considered for problems where the high dimensionality challenges the traditional acquisition function optimization and the budget is small enough that we cannot expect to learn a high-fidelity representation of the function.

Variable selection

If variables are completely removed beforehand in a preliminary stage, then a fixed value may be used, returning to a low dimensional problem. Observations with other values are generally discarded to keep the problem deterministic, the alternative being to add some noise if the inert variables can vary. Otherwise, after optimizing (2) on few variables, values for the screened variables must be determined to evaluate f . Alternatives include fixing these to a constant value, taking the values at the best design sampled so far for these coordinates, a random sampling, or a combination of these, see e.g., [START_REF] Li | High Dimensional Bayesian Optimization using Dropout[END_REF]; [START_REF] Spagnol | Global sensitivity analysis for optimization with variable selection[END_REF]. In [START_REF] Salem | Sequential dimension reduction for learning features of expensive black-box functions[END_REF], alternative lengthscales are estimated for the remaining variables by finding the most different values still passing a likelihood ratio test.

Then, values for these variables are selected where the difference between the predictive means are the most different between the two sets of hyperparameters, to challenge the initial split.

When considering the acquisition function, the main advantage of variable selection approaches is to reduce the dimension of the search space, which lessens the burden of the acquisition optimization. On the other hand, the problem doesn't decompose, as in the additive case, and may be using suboptimal values for the "inactive" variables.

Embedding case

With embeddings, additional difficulties arise with bounded domains. Let W = [A W 2 ] be a basis of R d . Splitting between active and inactive (or less active) variables: ∀x ∈ R d ,

x = WW x = AA x + W 2 W 2 x = Ay + W 2 z, y ∈ R r , z ∈ R d-r . If f has a true active subspace, the problem becomes find y * ∈ argmin y∈Y⊆R r f (p X (Ay)) (9)
with p X the convex projection onto X ; otherwise, the problem is:

find (y * , z * ) ∈ argmin y∈Y⊆R r ,z∈Z⊆R d-r f (Ay + W 2 z) s.t. Ay + W 2 z ∈ X .
As discussed by [START_REF] Constantine | Active subspaces: Emerging ideas for dimension reduction in parameter studies[END_REF], perhaps the optimization over z is secondary, which can be rewritten in the form

find y * ∈ argmin y∈Y⊆R r min z∈Z⊆R d-r f (Ay + W 2 z) s.t. Ay ∈ X
where less effort can be put on solving the high-dimensional subproblem involving z. This can be performed by random sampling [START_REF] Constantine | Active subspaces: Emerging ideas for dimension reduction in parameter studies[END_REF]. Alternatively, we can take a page from optimization in the context of variable selection and fix the inactive directions, see e.g., [START_REF] Cartis | Constrained global optimization of functions with low effective dimensionality using multiple random embeddings[END_REF]. Innocuous for unbounded domains4 , the constraint to belong to the domain for evaluation is complex for a compact domain X . The optimization problem is noisy whenever A changes over iterations as design points need to be projected on the embedding (that is, when a z component is present but ignored in the modeling of f ).

Discussing solely the (centered) [-1, 1] d hypercubic case, the intersection with a linear embedding is a convex polytope, defined by Y = {y ∈ R r s.t. -1 ≤ Ay ≤ 1}. This is the domain advocated by [START_REF] Letham | Re-Examining Linear Embeddings for High-Dimensional Bayesian Optimization[END_REF] to solve (9), at the extra cost of handling the linear constraints in optimizing α. With random A, a simpler option is to optimize y within a hypercubic domain Y = [-l, l] r , with l chosen based on the probability of finding a solution, as derived in [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF][START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF]; [START_REF] Qian | Derivative-Free Optimization of High-Dimensional Non-Convex Functions by Sequential Random Embeddings[END_REF]; [START_REF] Cartis | Constrained global optimization of functions with low effective dimensionality using multiple random embeddings[END_REF].

These results are simpler if the true active subspace is axis aligned, and depend on the difference between r and the true low dimension (if it exists). Nonetheless, these choices for Y may not even contain y * . The smallest compact set for this is a star-shaped polygon described in Binois et al. (2020). The reason to focus on smaller Y is two-fold. So when it comes to embeddings, they more often complicate than simplify the acquisition function search. This can still be worth the extra effort when a sufficient evaluation budget exists to learn the embedding.

These problems of domain selection are amplified for non-linear dimension reduction as illustrated by [START_REF] Siivola | Good practices for Bayesian optimization of high dimensional structured spaces[END_REF]. [START_REF] Li | High Dimensional Bayesian Optimization via Restricted Projection Pursuit Models[END_REF] faces both additivity and embedding issues, restricting the extent of the projection term to remain mostly in the domain. Next we provide pointers to benchmark problems to empirically assess the various methods detailed above.

Synthetic test optimization problems

Due to the multiplicity of possible methods and combinations, a thorough empirical comparison is out of scope in this work. Besides, implementations are not always available or even compatible. Inference techniques may differ as well, complicating comparisons. A less obvious difficulty is the absence of standard benchmark functions for high-dimensional surrogate-based optimization. We thus list here options we have come across.

First, some of the standard benchmark functions from global optimization, say as in [START_REF] Hansen | COCO: A platform for comparing continuous optimizers in a black-box setting[END_REF], can be scaled to high dimension. These synthetic test functions may be grouped into those which are separable, unimodal with moderate (resp. high) conditioning and finally multi-modal with global structure. Only the last one is in the realm of BO, as stated by [START_REF] Diouane | TREGO: a Trust-Region Framework for Efficient Global Optimization[END_REF]. But even these are often pathological, for instance with exponentially many local optima as d increases, and consequently demand a greater budget than could be reasonably afforded to find the global solution, regardless of how quickly the simulator can be evaluated. Defining reasonable, sub-globally-optimal targets for these scalable problems under limited budgets could be an option. Another option is to add artificial, inactive variables to classical low dimensional test functions, making the problem nominally high dimensional yet intrinsically low dimensional. While this scales easily to billions of variables as in [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF], it may not be realistic. Repeated versions of these low dimensional toy functions are also proposed, e.g., by [START_REF] Oh | BOCK: Bayesian Optimization with Cylindrical Kernels[END_REF]. This approach simply sums many low-dimensional versions of the same function together, allowing for more active variables. Another option is to rotate the low-dimensional function via a linear transformation, with two snags: the initial hypercubic domain may not fill the entire high-dimensional one, or parts may be mapped outside, including known optima. Analytical functions can be extended easily to fill the gaps, but the difficulty of the problem becomes rotation dependent. A set of fixed rotations could be used for ease of comparison, or embedding matrices like the class proposed in Nayebi et al.

(2019) that do not have the same domain issues.

GP realizations are an easy option to test the estimation difficulty of a given structure, and may also be used for comparisons between structures as in [START_REF] Ginsbourger | On ANOVA decompositions of kernels and Gaussian random field paths[END_REF].

The extent to which these results apply to real applications remains a topic of research. To this end, a list of problems which may be reproduced using publicly available software is provided in Table 1. Some of these problems arise from hyperparameter tuning applications that are popular in the machine learning community, and consist of tuning various neural networks properties with BO. The choice of the neural networks structure seems to often be hand tuned, the results are noisy and platform dependent, but the code is sometimes available to reproduce results. The most realistic examples comes from engineering and simulation, but the pipeline to reproduce high-fidelity simulation is rarely available in these cases. Besides the ones in Table 1 for which this is the case, examples include the design of airfoils, wings or fans [START_REF] Chen | Airfoil Design Parameterization and Optimization Using Bézier Generative Adversarial Networks[END_REF][START_REF] Gaudrie | Modeling and optimization with Gaussian processes in reduced eigenbases[END_REF][START_REF] Palar | Exploiting active subspaces in global optimization: how complex is your problem?[END_REF][START_REF] Viswanath | Dimension reduction for aerodynamic design optimization[END_REF][START_REF] Lukaczyk | Active subspaces for shape optimization[END_REF][START_REF] Seshadri | Dimension reduction via Gaussian ridge functions[END_REF], automotive industry test cases [START_REF] Binois | Uncertainty quantification on Pareto fronts and high-dimensional strategies in Bayesian optimization, with applications in multi-objective automotive design[END_REF], alloy design [START_REF] Li | High Dimensional Bayesian Optimization using Dropout[END_REF][START_REF] Rana | High dimensional Bayesian optimization with elastic Gaussian process[END_REF], biology [START_REF] Ulmasov | Bayesian optimization with dimension scheduling: Application to biological systems[END_REF], physics [START_REF] Kirschner | Adaptive and safe Bayesian optimization in high dimensions via one-dimensional subspaces[END_REF][START_REF] Mutny | Efficient high dimensional Bayesian optimization with additivity and quadrature Fourier features[END_REF], and electronics [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]. High dimensional data sets are another proper option to test high-dimensional GP regression, and can be extended to BO by allowing the BO algorithm to select not arbitrary points in the input space, but only those points corresponding to observed data. However, this discounts the effect of continuous optimization of the acquisition function, a fully nonneglible component of BO on real black-boxes. Still, some common data sets for regression tasks that have been used are listed in Table 2. One option, following the example of [START_REF] Jones | Large-scale multi-disciplinary mass optimization in the auto industry[END_REF] presenting at MOPTA, is to fit an interpolant to a given dataset and then treat that interpolant as the black-box to be optimized. We should remark that if the interpolant used to generate the "black-box" is the same as that used for BO, the results may be unrealistically rosy.

From the set of examples listed above, we see that the definition of "high dimensional"

in BO varies extensively, from the low tens to several thousands. The optimization budgets share the same range of variation, but each analysis is conducted only for a fixed budget.

There is thus room for a thorough study on which high-dimensional structure is the most adapted given the dimension and budget, or even better, enabling more complexity as more evaluations become available.

Conclusion and perspectives

Often rightly presented as one of the top challenges for GP-based BO, high-dimensionality has generated many different ideas to address each aspect of the manifestation of the curse of dimensionality. While the structural assumptions seem to cluster into one or more of variable selection, additive decomposition or linear embeddings, no consensus is present on the best way to tackle the problem. Even if this is largely a problem-dependent issue, the inference methods used also restrict which structures are possible to assume, whether relying on a random structure or actually trying to infer it. Even besides these specifics, the particular strategy used for acquisition function optimization can overshadow the rest.

More systematic comparisons are necessary, as is the definition of suitable benchmarks.

Some general guidelines

As mentioned at the beginning of Section 3, most structural model assumptions can be seen as instances of the more general model:

model: f (x) ≈ κ i=1 g i (A i x) (10)
with A i ∈ R r×d , κ ∈ N * , which allows for sums of low dimensional effects, This includes at one end the additive Gaussian process of [START_REF] Durrande | Additive Kernels for Gaussian Process Modeling[END_REF] which separates each variable individually into its own bin, at another end linear embeddings of the form f (x) = h(Ax) which ignore certain input directions, and beyond lie the functional ANOVA kernels of [START_REF] Muehlenstaedt | Data-driven Kriging models based on FANOVA-decomposition[END_REF] or in combining several random low dimensional subspaces such as in [START_REF] Delbridge | Randomly projected additive Gaussian processes for regression[END_REF]. There are still promising research directions at the intersections of these methods, and we argue in particular that "default" kernel components, for instance a simple isotropic kernel, could be useful in complementing methods which search a particular part of the input space. The various structural modeling options are presented in Figure 2, highlighting the links between options.

Rather than starting directly from the most general model, whose direct inference would make a very steep initial step, we recommend to try the simpler models first. That is, begin Alternatively, if the budget of optimization is limited, a simple filling strategy for the "inactive" variables combined with optimization on some restricted subspace could be preferable. Finally the use of the trust region framework has shown promising results

for optimization in high dimensional problems, and is amenable to parallelization.

Future work

There is potential for hybrids between models, inference methods and mathematical opti- The ability to doubt the current structural assumption may also be quite helpful in order to avoid negative feedback loop effects appearing in adaptive design, see e.g., [START_REF] Gramacy | Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences[END_REF].

Whether the ability to keep a kernel defined on the whole initial domain is superior to defining one directly on a low dimensional embedding (which while appealing is associated with domain selection issues) remains to be seen. If a kernel defined on an embedding is indeed preferred, keeping kernels sufficiently expressive (e.g., not restricted to be first order additive like in [START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis[END_REF]) by adding an isotropic component may be improved by considering ortho-additivity [START_REF] Ginsbourger | On ANOVA decompositions of kernels and Gaussian random field paths[END_REF] instead. These orthogonality constraints are likely to help inference by reducing identifiability issues.

Additionally, the success of trust region BO approaches which restrict the search space to some subset of the input space motivates dimension reduction approaches which are themselves local. For instance Wycoff (2021a) propose to define an active subspace with respect to probability measures with non-uniform densities in order to emphasize certain regions of interest. Namely, the focus is restricted to a trust region, but other schemes could be considered. Of course, a dimension reduction which is locally linear is in fact globally nonlinear, and such local linear models may prove to be a tractable approach to developing nonlinear dimension reduction for BO.

Furthermore, there have already been useful results born of the influence of mathematical programming on approaches to high dimensional Bayesian optimization. In particular, use of a rectangular trust region seems to be of benefit in [START_REF] Eriksson | Scalable Global Optimization via Local Bayesian Optimization[END_REF] and a

Gaussian process analogue of the augmented Lagrangian method is developed in [START_REF] Gramacy | Modeling an Augmented Lagrangian for Blackbox Constrained Optimization[END_REF]; [START_REF] Picheny | Bayesian optimization under mixed constraints with a slack-variable augmented Lagrangian[END_REF]. Helping to develop new algorithms which combine benefits from low dimensional modeling with trust regions or other approaches successful in high dimensional mathematical programming broadly could help BO finally make the leap to practical, widespread usage on complex, high dimensional problems.

A Synthetic benchmark functions

Some global optimization benchmark functions are used for testing high-dimensional BO, of which a list is provided in Table 3.

Figure 1 :

 1 Figure 1: Distances Concentrate in High Dimension: Randomly sampling 10,000 points according to the uniform measure in [0, 1] d and then calculating squared interpoint distances reveals that these concentrate within the bounds of possible values in high dimension: [0, d].

  as introduced by[START_REF] Stitson | Support vector regression with ANOVA decomposition kernels[END_REF]. There the sensitivity indices of the ANOVA representation are analytically tractable. In[START_REF] Ginsbourger | On ANOVA decompositions of kernels and Gaussian random field paths[END_REF], the ANOVA decomposition is performed directly on the kernel and propagated to the corresponding random field under appropriate orthogonality conditions. Looking for sparsity to avoid estimating the entirety of the 2 d components, they define projectors that allow separation of the additive components (with cross-covariance between main effects) from their complement.

(

  2020) defines an alternative mapping from Y to X to avoid it.[START_REF] Nayebi | A framework for Bayesian optimization in embedded subspaces[END_REF] also bypasses this problem by choosing a sparse random matrix with {-1, 0, 1} elements only, in essence selecting from embeddings on diagonals of the hypercube. These ideas can be extended to affine embeddings, as proposed for global optimization by[START_REF] Cartis | Constrained global optimization of functions with low effective dimensionality using multiple random embeddings[END_REF];[START_REF] Cartis | A dimensionality reduction technique for unconstrained global optimization of functions with low effective dimensionality[END_REF].Then comes the question of how to use high dimensional data that are not on the linear embedding (unlike for REMBO). Unless A is perfectly recovered, introducing a noise term is necessary to account for the discrepancy. To this end,[START_REF] Moriconi | High-dimensional Bayesian optimization with projections using quantile Gaussian processes[END_REF] uses GP regression on quantiles for axis-aligned projections. Another question is how to select the lengthscales in this context since the product kernels are not preserved in the embedded space, for which[START_REF] Letham | Re-Examining Linear Embeddings for High-Dimensional Bayesian Optimization[END_REF] show that a specific parameterization (Mahalanobis distance based kernel) is preferable and handles distortions. Still, inference is complicated whereas the REMBO program needs only to fit a low dimensional GP. What distribution to sample A from has also been briefly investigated in[START_REF] Binois | Uncertainty quantification on Pareto fronts and high-dimensional strategies in Bayesian optimization, with applications in multi-objective automotive design[END_REF];Binois et al. (2020);[START_REF] Letham | Re-Examining Linear Embeddings for High-Dimensional Bayesian Optimization[END_REF].

  chosen to further focus on the center[START_REF] Oh | BOCK: Bayesian Optimization with Cylindrical Kernels[END_REF] while the angular component is a continuous radial kernel K a (a, a ) = P p=0 c p (a a ) p , c p > 0, ∀p (Jayasumana

  suggest to artificially inflate the gradient by taking a larger lengthscale within an isotropic GP. Then the lengthscale is successively decreased while tracking subsequent optima (warm started by previous values). Applying the random linear embedding of REMBO (Section 3.3) at the level of acquisition function optimization, Tran-The et al. (2019) optimize over a finite set of subspaces.

  First, most of the difference corresponds to points on the boundary of X , with increasingly large volume compared to the intersection, bringing back the curse of dimensionality. Second is that the corresponding observations are distorted with the convex projection and are thus harder to model. These issues are by-passed by[START_REF] Nayebi | A framework for Bayesian optimization in embedded subspaces[END_REF];[START_REF] Moriconi | High-dimensional Bayesian optimization with projections using quantile Gaussian processes[END_REF], at the cost of the stricter diagonal or axis-aligned embedding assumptions. Changing the domain, as is implicitly done by choosing l relates to another strategy not explored much in BO: to reduce the optimization space, as discussed by[START_REF] Stork | Open issues in surrogate-assisted optimization[END_REF] or employed with TRs.A weighted-PCA estimation of A is proposed in[START_REF] Raponi | High dimensional Bayesian optimization assisted by principal component analysis[END_REF], and they handle the domain issue with a penalty for infeasibility.Chen et al. (2020a) use a semi-supervised version of sliced inverse regression (SIR) to find important input directions, using both labeled (evaluated) and unlabeled (unevaluated) designs. Namely, they collect points which have high acquisition function values but were not chosen for evaluation into an "unlabeled dataset", which they incorporate into the estimation of the embedding. The "inverse" in SIR indicates the goal of swapping the role of x and y, finding which values of the designs lead to given values of the outputs. Selecting a domain suffers the same difficulties as above, plus those related to updating the SIR model.Most existing work fixes the embedding, or does not propagate the uncertainty when estimated iteratively. Exploring the interplay between acquisition functions defined over AS functions, e.g., in[START_REF] Garnett | Active learning of linear embeddings for Gaussian processes[END_REF];[START_REF] Wycoff | Gradient-Based Sensitivity Analysis with Kernels[END_REF] and those defined for optimization is a promising direction to enhance both of them.

  benchmarking too: Borehole (8d), OTL circuit (6d), piston (7d), robot arm (8d), and wingweight (10d). The power circuit function (13d)[START_REF] Lee | Modified Active Subspaces Using the Average of Gradients[END_REF]) makes one additional such example. Several functions designed for sensitivity analysis have been experimented on as well: the Sobol g-function (any d), the Ishigami function (3d), and the 8d flood model in[START_REF] Iooss | A review on global sensitivity analysis methods[END_REF]. Assessing the adequacy of these analytical examples with the structural model assumptions is interesting future work.

  with a standard GP to build an alternative one based on the ARD principle, or to compute sensitivity indices like Sobol indices. Then build a fully additive GP model or a single index one as the basic models of the additive and linear embedding families. From comparing the prediction given by these simple models, building more advanced models from the best performing one is possible: functional ANOVA or block additive from the additive one; linear embedding on top of single index. If none is performing well, projection pursuit or non-linear embeddings could be entertained, if n is sufficiently large to allow a reliable inference of the corresponding model.More precisely for the inference, depending on d, the budget n and the application dependent complexity of f , inference by maximum likelihood (or, even better, a Bayesian alternative) is preferred. If too difficult or slow when d reaches the hundreds, relying on a random structure to bypass the full-dimensional model inference remains possible. The simplest models scale well in computation with increasing d, while the more complex ones may benefit from recent advances brought by the increasing use of GPUs and automatic differentiation frameworks. This is especially true as more complex models require larger values of n to show improvements.

  mization techniques. For global optimization, a research avenue pointed out by Spagnol et al. (2019) is to focus on variables that are important to reach low function values, rather than evaluating their influence on the whole range of f . Chen et al. (2020a) go some way to focusing on low values, as do Guhaniyogi and Dunson (2016) on manifold-bound inputs.

  [START_REF] Hoang | Decentralized highdimensional Bayesian optimization with factor graphs[END_REF] also rely on random groups while[START_REF] Rolland | High-Dimensional Bayesian Optimization via Additive Models with Overlapping Groups[END_REF] use a dependence graph and Gibbs sampling to perform inference. In this vein of methods, there is[START_REF] Delbridge | Randomly projected additive Gaussian processes for regression[END_REF] that reconstruct the high dimensional kernel by a sum of univariate ones in random directions, in the spirit of the turning band method from geostatistics

[START_REF] Gardner | Discovering and Exploiting Additive Structure for Bayesian Optimization[END_REF]

attempts to elicit the structure via a dedicated Metropolis-Hastings algorithm. For overlapping subsets,

  Running a sensitivity analysis independent of the GP model to select A as a kind of preprocessing step is also an option. When gradients of the black-box are available, recovering the matrix A is relatively straightforward. It corresponds, up to a rotation, to the eigenvectors with non-zero eigenvalues of the matrix C

	Seshadri et al. (2019) reinterpreted a similar problem,
	drawing connections between regression and approximation. Also building on Tripathy
	et al. (2016), Yenicelik (2020) observed that the likelihood may not always select the best
	matrix among alternatives. A full Bayesian treatment was proposed by Marcy (2018), with
	priors on matrix manifolds and the dimension, requiring advanced Monte Carlo methods.

  in the context of AS. Also defined by analogy to AS is Bridges et al. (2019)'s 1D active manifold, not yet applied to

  Table 3 of the Appendix enumerates these global optimization functions. Additional ones, created for BO, are tested by Siivola et al. (2018). The five physically-motivated analytical functions from https://www.sfu.ca/ ~ssurjano/optimization.html have been the subject of BO

Table 1 :

 1 Examples of reproducible, realistic test problems.

	Name(s)	Field/Domain	d	Example reference
	Bayesian neural net Hyperparameter tuning	5	Falkner et al. (2018)
	Cartpole swing-up task Hyperparameter tuning	8	Falkner et al. (2018)
	Covid-19 model	Epidemiology	8	Wycoff et al. (2021)
	Elliptic PDE	Simulation	100	Tripathy et al. (2016)
	Groundwater remediation	Environment	6	Gramacy (2020)
	Neural net Hyperparameter tuning	6	Falkner et al. (2018)
	Robot pushing	Robotics	14	Wang et (2018)
	Rover trajectory	Robotics	60	Wang et al. (2018)
	Satellite drag	Astrodynamics	7	Gramacy (2020)
	SDSS	Cosmology	9	Kandasamy et al. (2015)
	Shallow Net Weights Hyperparameter tuning 100-500 Oh et al. (2018)
	Three-link walk robot	Robotics	25	Zhang et al. (2019)
	VJ Cascade classifier Hyperparameter tuning	22	Kandasamy et al. (2015)

Table 2 :

 2 Example of data sets for regression used for GPs and BO.

	Name(s)	d	n	Example reference
	Airfoil	5	1,503 Francom and Sansó (2019)
	Amazon commerce reviews 10,000 1,500 Fukumizu and Leng (2014)
	Autos	24	159	Delbridge et al. (2020)
	Communities and crime	96	1,993 Garnett et al. (2014)
	Concrete	8	1,030 Garnett et al. (2014)
	Crash	9	20	Gramacy (2020)
	CT slices	318	3,071 Garnett et al. (2014)
	Elevators	17	8,752 Gilboa et al. (2013)
	Engine block and head joint sealing	8	27	Joseph et al. (2008)
	Heat exchanger	4	64	Lin and Joseph (2020)
	Housing	13	506	Delbridge et al. (2020)
	Isomap faces 4,096	698	Guhaniyogi and Dunson (2016)
	Kink40k	8	10,000 Gilboa et al. (2013)
	LGBB	3	3,167 Gramacy (2020)
	MARTHE	20	300	Marrel et al. (2008)
	Olivetti faces data set (modified) 1,935	400	Delbridge et al. (2020)
	PIC	9	45,730 Hoang et al. (2018)
	Pumadyn	8	7,168 Gilboa et al. (2013)
	Pumadyn	32	7,168 Gilboa et al. (2013)
	Sarcos robot	21	44,484 Winkel et al. (2021)
	Servo	3	167	Duvenaud et al. (2011)
	Tecator	100	215	Plate (1999)
	Temperature	106	10,675 Garnett et al. (2014)
	Yacht	6	308	Garnett et al. (2014)

  Figure2: A taxonomy of structural model assumptions, going from simpler (top) to more general ones (bottom). Arrows marks models that can be generalized, ( †) indicates that the estimation of a noise term is needed and ( ‡) cases for which a filling strategy is needed unless optimization is conducted only on a low dimensional manifold.

	Variable		1st order		Single
	selection † ‡ (3)	additive † (4)	index † ‡ (7)
			fANOVA † (6)	High	Linear
				order	embedding † ‡ (8)
			additive † (5)
	Anisotropic	Isotropic	Sum of linear	Projection	Non-linear
	with all	with all	embeddings † ‡	pursuit † ‡	embedding † ‡
	variables	variables	(10)

Table 3 :

 3 Synthetic problems for global optimization benchmarks tested on for highdimensional BO. See e.g., https://www.sfu.ca/ ~ssurjano/optimization.html for further details, expressions and codes, or[START_REF] Cartis | A dimensionality reduction technique for unconstrained global optimization of functions with low effective dimensionality[END_REF]. Cola is used byBinois et al. (2020) while the Beach, Dream and Simba functions can be found inWinkel et al.

	(2021).
	Name(s)	d	Name(s)	d
	Ackley Any	Hartmann 3	3
	Beach	6	Hartmann 6	6
	Beale	2	Levy Any
	Branin	2	Michalewicz Any
	Brent	2	Perm 4, 0.5	4
	Bukin N.6	2	Quadratic Form Any
	Camel	2	Rastrigin any
	Cola 17	Rosenbrock Any
	Colville	4	Simba	6
	Dream	6	Shekel 5,7,10	4
	Easom	2	Shubert	2
	Franke	2	Styblinski Tang Any
	Friedman	5	Trid any
	Giunta	2	Welsh 20
	Goldstein-Price	2	Zakharov any
	Griewank Any	Zettl	2

For more details on what exactly is meant by the derivative of a stochastic process see(Papoulis, 1965, Chapter 10).

The results follow from the fact that the density of the squared distance d between two uniform variables on [0, 1] is given by 1 √ d -1.

A toy example is sin(

4π d 1 x).

At least in the formulation, unbounded domains are another topic in BO, see e.g.,Shahriari et al. (2016a).

Zhou, J., Yang, Z., Si, Y., Kang, L., Li, H., Wang, M., and Zhang, Z. (2021). "A Trust-Region Parallel Bayesian Optimization Method for Simulation-Driven Antenna Design." Antennas and Propagation,69,7,[3966][3967][3968][3969][3970][3971][3972][3973][3974][3975][3976][3977][3978][3979][3980][3981]

IEEE Transactions on