Convergence rates of efficient global optimization algorithms

Introduction

Many standard global optimization algorithms exist, including genetic algorithms, multistart, and simulated annealing (Pardalos and Romeijn, 2002), but these algorithms are designed for functions that are cheap to evaluate. When $f$ is expensive, we need an efficient algorithm, one which will choose its observations to maximize the information gained.

We can consider this a continuum-armed-bandit problem (Srinivas et al., 2010, and references therein), with noiseless data, and loss measured by the simple regret (Bubeck et al., 2009). At time $n$, we choose a design point $x_{n}\in X$, make an observation $z_{n}=f(x_{n})$, and then report a point $x_{n}^{*}$ where we believe $f(x_{n}^{*})$ will be low. Our goal is to find a strategy for choosing the $x_{n}$ and $x_{n}^{*}$, in terms of previous observations, so as to minimize $f(x_{n}^{*})$.

We would like to find a strategy which can guarantee convergence: for functions $f$ in some smoothness class, $f(x_{n}^{*})$ should tend to $\min f$, preferably at some fast rate. The simplest method would be to fix a sequence of $x_{n}$ in advance, and set $x^{*}_{n}=\arg\min\hat{f}_{n}$, for some approximation $\hat{f}_{n}$ to $f$. We will show that if $\hat{f}_{n}$ converges in supremum norm at the optimal rate, then $f(x_{n}^{*})$ also converges at its optimal rate. However, while this strategy gives a good worst-case bound, on average it is clearly a poor method of optimization: the design points $x_{n}$ are completely independent of the observations $z_{n}$.

We may therefore ask if there are more efficient methods, with better average-case performance, that nevertheless provide good guarantees of convergence. The difficulty in designing such a method lies in the trade-off between exploration and exploitation. If we exploit the data, observing in regions where $f$ is known to be low, we will be more likely to find the optimum quickly; however, unless we explore every region of $X$, we may not find it at all (Macready and Wolpert, 1998).

Initial attempts at this problem include work on Lipschitz optimization (summarized in Hansen et al., 1992) and the DIRECT algorithm (Jones et al., 1993), but perhaps the best-known strategy is expected improvement. It is sometimes called Bayesian optimization, and first appeared in Močkus (1974) as a Bayesian decision-theoretic solution to the problem. Contemporary computers were not powerful enough to implement the technique in full, and it was later popularized by Jones et al. (1998), who provided a computationally efficient implementation. More recently, it has also been called a knowledge-gradient policy by Frazier et al. (2009). Many extensions and alterations have been suggested by further authors; a good summary can be found in Brochu et al. (2010).

Vazquez and Bect (2010) show that when $\pi$ is a fixed Gaussian process prior of finite smoothness, expected improvement converges on the minimum of any $f\in\mathcal{H}$, and almost surely for $f$ drawn from $\pi$. Grunewalder et al. (2010) bound the convergence rate of a computationally infeasible version of expected improvement: for priors $\pi$ of smoothness $\nu$, they show convergence at a rate $O^{*}(n^{-(\nu\wedge 0.5)/d})$ on $f$ drawn from $\pi$. We begin by bounding the convergence rate of the feasible algorithm, and show convergence at a rate $O^{*}(n^{-(\nu\wedge 1)/d})$ on all $f\in\mathcal{H}$. We go on to show that a modification of expected improvement converges at the near-optimal rate $O^{*}(n^{-\nu/d})$.

For practitioners, however, these results are somewhat misleading. In typical applications, the prior is not held fixed, but depends on parameters estimated sequentially from the data. This process ensures the choice of observations is invariant under translation and scaling of $f$, and is believed to be more efficient (Jones et al., 1998, §2). It has a profound effect on convergence, however: Locatelli (1997, §3.2) shows that, for a Brownian motion prior with estimated parameters, expected improvement may not converge at all.

We extend this result to more general settings, showing that for standard priors with estimated parameters, there exist smooth functions $f$ on which expected improvement does not converge. We then propose alternative estimates of the prior parameters, chosen to minimize the constants in the convergence rate. We show that these estimators give an automatic choice of parameters, while retaining the convergence rates of a fixed prior.

In Section 2, we briefly describe the expected-improvement algorithm, and detail our assumptions on the priors used. We state our main results in Section 3, and discuss implications for further work in Section 4. Finally, we give proofs in Appendix A.

Expected Improvement

and our goal, given $\pi$, is to choose the strategy $u$ to minimize this quantity.

For $N>1$ this problem is very computationally intensive (Osborne, 2010, §6.3), but we can solve a simplified version of it. First, we restrict the choice of $x_{n}^{*}$ to the previous design points $x_{1},\dots,x_{n}$. (In practice this is reasonable, as choosing an $x_{n}^{*}$ we have not observed can be unreliable.) Secondly, rather than finding an optimal strategy for the problem, we derive the myopic strategy: the strategy which is optimal if we always assume we will stop after the next observation. This strategy is suboptimal (Ginsbourger et al., 2008, §3.1), but performs well, and greatly simplifies the calculations involved.

In this setting, given $\mathcal{F}_{n}$, if we are to stop at time $n$ we should choose $x^{*}_{n}\coloneqq x_{i_{n}^{*}}$, where $i_{n}^{*}\coloneqq\arg\min_{1,\dots,n}z_{i}$. (In the case of ties, we may pick any minimizing $i_{n}^{*}$.) We then suffer a loss $z_{n}^{*}-\min f$, where $z_{n}^{*}\coloneqq z_{i_{n}^{*}}$. Were we to observe at $x_{n+1}$ before stopping, the expected loss would be

so the myopic strategy should choose $x_{n+1}$ to minimize this quantity. Equivalently, it should maximize the expected improvement over the current loss,

So far, we have merely replaced one optimization problem with another. However, for suitable priors, $EI_{n}$ can be evaluated cheaply, and thus maximized by standard techniques. The expected-improvement algorithm is then given by choosing $x_{n+1}$ to maximize (1).

2 Gaussian Process Models

We still need to choose a prior $\pi$ for $f$. Typically, we model $f$ as a stationary Gaussian process: we consider the values $f(x)$ to be jointly Gaussian, with mean and covariance

for an underlying kernel $K$ with $K(0)=1$. (Note that we can always satisfy this condition by suitably scaling $K$ and $\sigma$.) The $\theta_{i}>0$ are the length-scales of the process: two values $f(x)$ and $f(y)$ will be highly correlated if each $x_{i}-y_{i}$ is small compared with $\theta_{i}$. For now, we will assume the parameters $\sigma$ and $\theta$ are fixed in advance.

For (2) and (3) to define a consistent Gaussian process, $K$ must be a symmetric positive-definite function. We will also make the following assumptions.

and by Bochner’s theorem, $\widehat{K}$ is non-negative and integrable.

$\widehat{K}$ is isotropic and radially non-increasing.

In other words, $\widehat{K}(x)=\widehat{k}(\lVert x\rVert)$ for a non-increasing function $\widehat{k}:[0,\infty)\to[0,\infty)$; as a consequence, $K$ is isotropic.

$\widehat{K}(x)=\Theta(\lVert x\rVert^{-2\nu-d})$ for some $\nu>0$; or

$\widehat{K}(x)=O(\lVert x\rVert^{-2\nu-d})$ for all $\nu>0$ (we will then say that $\nu=\infty$).

Note the condition $\nu>0$ is required for $\widehat{K}$ to be integrable.

$K$ is $C^{k}$, for $k$ the largest integer less than $2\nu$, and at the origin, $K$ has $k$-th order Taylor approximation $P_{k}$ satisfying

When $\alpha=0$, this is just the condition that $K$ be $2\nu$-Hölder at the origin; when $\alpha>0$, we instead require this condition up to a log factor.

The rate $\nu$ controls the smoothness of functions from the prior: almost surely, $f$ has continuous derivatives of any order $k<\nu$ (Adler and Taylor, 2007, §1.4.2). Popular kernels include the Matérn class,

where $k_{\nu}$ is a modified Bessel function of the second kind, and the Gaussian kernel,

Having chosen our prior distribution, we may now derive its posterior. We find

for $z=(z_{i})_{i=1}^{n}$, $V=(K_{\theta}(x_{i}-x_{j}))_{i,j=1}^{n}$, and $v=(K_{\theta}(x-x_{i}))_{i=1}^{n}$ (Santner et al., 2003, §4.1.3). Equivalently, these expressions are the best linear unbiased predictor of $f(x)$ and its variance, as given in Jones et al. (1998, §2). We will also need the reduced sum of squares,

3 Expected Improvement Strategies

Under our assumptions on $\pi$, we may now derive an analytic form for (1), as in Jones et al. (1998, §4.1). We obtain

and $\Phi$ and $\varphi$ are the standard normal distribution and density functions respectively.

For a prior $\pi$ as above, expected improvement chooses $x_{n+1}$ to maximize (8), but this does not fully define the strategy. Firstly, we must describe how the strategy breaks ties, when more than one $x\in X$ maximizes $EI_{n}$. In general, this will not affect the behaviour of the algorithm, so we allow any choice of $x_{n+1}$ maximizing (8).

Secondly, we must say how to choose $x_{1}$, as the above expressions are undefined when $n=0$. In fact, Jones et al. (1998, §4.2) find that expected improvement can be unreliable given few data points, and recommend that several initial design points be chosen in a random quasi-uniform arrangement. We will therefore assume that until some fixed time $k$, points $x_{1},\dots,x_{k}$ are instead chosen by some (potentially random) method independent of $f$. We thus obtain the following strategy.

initial design points $x_{1},\dots,x_{k}$ independently of $f$; and

further design points $x_{n+1}\ (n\geq k)$ from the maximizers of (8).

So far, we have not considered the choice of parameters $\sigma$ and $\theta$. While these can be fixed in advance, doing so requires us to specify characteristic scales of the unknown function $f$, and causes expected improvement to behave differently on a rescaling of the same function. We would prefer an algorithm which could adapt automatically to the scale of $f$.

A natural approach is to take maximum likelihood estimates of the parameters, as recommended by Jones et al. (1998, §2). Given $\theta$, the MLE $\hat{\sigma}^{2}_{n}=\hat{R}_{n}^{2}(\theta)/n$; for full generality, we will allow any choice $\hat{\sigma}^{2}_{n}=c_{n}\hat{R}_{n}^{2}(\theta)$, where $c_{n}=o(1/\log n)$. Estimates of $\theta$, however, must be obtained by numerical optimization. As $\theta$ can vary widely in scale, this optimization is best performed over $\log\theta$; as the likelihood surface is typically multimodal, this requires the use of a global optimizer. We must therefore place (implicit or explicit) bounds on the allowed values of $\log\theta$. We have thus described the following strategy.

Let $\hat{\pi}_{n}$ be a sequence of priors, with parameters $\hat{\sigma}_{n}$, $\hat{\theta}_{n}$ satisfying:

$\hat{\sigma}^{2}_{n}=c_{n}\hat{R}_{n}^{2}(\hat{\theta}_{n})$ for constants $c_{n}>0$, $c_{n}=o(1/\log n)$; and

An $EI(\hat{\pi})$ strategy satisfies 1, replacing $\pi$ with $\hat{\pi}_{n}$ in (8).

Convergence Rates

which holds for all $f\in\mathcal{E}(S)$. See Aronszajn (1950), Berlinet and Thomas-Agnan (2004), Wendland (2005) and van der Vaart and van Zanten (2008).

and there is a unique $g$ minimizing this expression.

is finite. Thus, for the kernel $K$ with Fourier transform $\widehat{K}(\xi)=(1+\lVert\xi\rVert^{2})^{s/2}$, this is just the RKHS $\mathcal{H}(D)$. More generally, if $K$ satisfies our assumptions with $\nu<\infty$, these spaces are equivalent in the sense of normed spaces: they contain the same functions, and have norms $\lVert\,\cdot\,\rVert_{1},\lVert\,\cdot\,\rVert_{2}$ satisfying

If $\nu<\infty$, $\mathcal{H}_{\theta}(\bar{D})$ is equivalent to the Sobolev Hilbert space $H^{\nu+d/2}(D)$.

If $\nu=\infty$, $\mathcal{H}_{\theta}(\bar{D})$ is continuously embedded in $H^{s}(D)$ for all $s$.

Thus if $\nu<\infty$, and $X$ is, say, a product of intervals $\prod_{i=1}^{d}[a_{i},b_{i}]$, the RKHS $\mathcal{H}_{\theta}(X)$ is equivalent to the Sobolev Hilbert space $H^{\nu+d/2}(\prod_{i=1}^{d}(a_{i},b_{i}))$, identifying each function in that space with its unique continuous extension to $X$.

2 Fixed Parameters

We will say that $u$ converges on the optimum at rate $r_{n}$, if

for all $R>0$. Note that we do not allow $u$ to vary with $R$; the strategy must achieve this rate without prior knowledge of $\lVert f\rVert_{\mathcal{H}_{\theta}(X)}$.

We begin by showing that the minimax rate of convergence is $n^{-\nu/d}$.

and this rate can be achieved by a strategy $u$ not depending on $R$.

The upper bound is provided by a naive strategy as in the introduction: we fix a quasi-uniform sequence $x_{n}$ in advance, and take $x_{n}^{*}$ to minimize a radial basis function interpolant of the data. As remarked previously, however, this naive strategy is not very satisfying; in practice it will be outperformed by any good strategy varying with the data. We may thus ask whether more sophisticated strategies, with better practical performance, can still provide good worst-case bounds.

One such strategy is the $EI(\pi)$ strategy of 1. We can show this strategy converges at least at rate $n^{-(\nu\wedge 1)/d}$, up to log factors.

For $\nu\leq 1$, these rates are near-optimal. For $\nu>1$, we are faced with a more difficult problem; we discuss this in more detail in Section 3.4.

3 Estimated Parameters

First, we consider the effect of the prior parameters on $EI(\pi)$. While the previous result gives a convergence rate for any fixed choice of parameters, the constant in that rate will depend on the parameters chosen; to choose well, we must somehow estimate these parameters from the data. The $EI(\hat{\pi})$ strategy, given by 2, uses maximum likelihood estimates for this purpose. We can show, however, that this may cause the strategy to never converge.

The counterexamples constructed in the proof of the theorem may be difficult to minimize, but they are not badly-behaved (Figure 1). A good optimization strategy should be able to minimize such functions, and we must ask why expected improvement fails.

We can understand the issue by considering the constant in Theorem 2. Define

From the proof of Theorem 2, the dominant term in the convergence rate has constant

for $C>0$ not depending on $R$ or $\sigma$. In Appendix A, we will prove the following result.

$\hat{R}_{n}(\theta)$ is non-decreasing in $n$, and bounded above by $\lVert f\rVert_{\mathcal{H}_{\theta}(X)}$.

Hence for fixed $\theta$, the estimate $\hat{\sigma}_{n}^{2}=\hat{R}_{n}^{2}(\theta)/n\leq R^{2}/n$, and thus $R/\hat{\sigma}_{n}\geq n^{1/2}$. Inserting this choice into (10) gives a constant growing exponentially in $n$, destroying our convergence rate.

To resolve the issue, we will instead try to pick $\sigma$ to minimize (10). The term $R+\sigma$ is increasing in $\sigma$, and the term $\tau(R/\sigma)/\tau(-R/\sigma)$ is decreasing in $\sigma$; we may balance the terms by taking $\sigma=R$. The constant is then proportional to $R$, which we may minimize by taking $R=\lVert f\rVert_{\mathcal{H}_{\theta}(X)}$. In practice, we will not know $\lVert f\rVert_{\mathcal{H}_{\theta}(X)}$ in advance, so we must estimate it from the data; from 1, a convenient estimate is $\hat{R}_{n}(\theta)$.

Suppose, then, that we make some bounded estimate $\hat{\theta}_{n}$ of $\theta$, and set $\hat{\sigma}_{n}^{2}=\hat{R}_{n}^{2}(\hat{\theta}_{n})$. As Theorem 3 holds for any $\hat{\sigma}^{2}_{n}$ of faster than logarithmic decay, such a choice is necessary to ensure convergence. (We may also choose $\theta$ to minimize (10); we might then pick $\hat{\theta}_{n}$ minimizing $\hat{R}_{n}(\theta)\prod_{i=1}^{d}\theta_{i}^{-\nu/d},$ but our assumptions on $\hat{\theta}_{n}$ are weak enough that we need not consider this further.)

If we believe our Gaussian-process model, this estimate $\hat{\sigma}_{n}$ is certainly unusual. We should, however, take care before placing too much faith in the model. The function in Figure 1 is a reasonable function to optimize, but as a Gaussian process it is highly atypical: there are intervals on which the function is constant, an event which in our model occurs with probability zero. If we want our algorithm to succeed on more general classes of functions, we will need to choose our parameter estimates appropriately.

To obtain good rates, we must add a further condition to our strategy. If $z_{1}=\dots=z_{n}$, $EI_{n}(\,\cdot\,;\hat{\pi}_{n})$ is identically zero, and all choices of $x_{n+1}$ are equally valid. To ensure we fully explore $f$, we will therefore require that when our strategy is applied to a constant function $f(x)=c$, it produces a sequence $x_{n}$ dense in $X$. (This can be achieved, for example, by choosing $x_{n+1}$ uniformly at random from $X$ when $z_{1}=\dots=z_{n}$.) We have thus described the following strategy.

we instead set $\hat{\sigma}_{n}^{2}=\hat{R}_{n}^{2}(\hat{\theta}_{n})$; and

we require the choice of $x_{n+1}$ maximizing (8) to be such that, if $f$ is constant, the design points are almost surely dense in $X$.

We cannot now prove a convergence result uniform over balls in $\mathcal{H}_{\theta}(X)$, as the rate of convergence depends on the ratio $R/\hat{R}_{n}$, which is unbounded. (Indeed, any estimator of $\lVert f\rVert_{\mathcal{H}_{\theta}(X)}$ must sometimes perform poorly: $f$ can appear from the data to have arbitrarily small norm, while in fact having a spike somewhere we have not yet observed.) We can, however, provide the same convergence rates as in Theorem 2, in a slightly weaker sense.

4 Near-Optimal Rates

So far, our rates have been near-optimal only for $\nu\leq 1$. To obtain good rates for $\nu>1$, standard results on the performance of Gaussian-process interpolation (Narcowich et al., 2003, §6) then require the design points $x_{i}$ to be quasi-uniform in a region of interest. It is unclear whether this occurs naturally under expected improvement, but there are many ways we can modify the algorithm to ensure it.

Perhaps the simplest, and most well-known, is an $\varepsilon$-greedy strategy (Sutton and Barto, 1998, §2.2). In such a strategy, at each step with probability $1-\varepsilon$ we make a decision to maximize some greedy criterion; with probability $\varepsilon$ we make a decision completely at random. This random choice ensures that the short-term nature of the greedy criterion does not overshadow our long-term goal.

The parameter $\varepsilon$ controls the trade-off between global and local search: a good choice of $\varepsilon$ will be small enough to not interfere with the expected-improvement algorithm, but large enough to prevent it from getting stuck in a local minimum. Sutton and Barto (1998, §2.2) consider the values $\varepsilon=0.1$ and $\varepsilon=0.01$, but in practical work $\varepsilon$ should of course be calibrated to a typical problem set.

We therefore define the following strategies.

chooses initial design points $x_{1},\cdots,x_{k}$ independently of $f$;

with probability $1-\varepsilon$, chooses design point $x_{n+1}\ (n\geq k)$ as in $EI(\,\cdot\,)$; or

with probability $\varepsilon$, chooses $x_{n+1}\ (n\geq k)$ uniformly at random from $X$.

We can show that these strategies achieve near-optimal rates of convergence for all $\nu<\infty$.

Let $EI(\,\cdot\,,\varepsilon)$ be one of the strategies in 4. If $\nu<\infty$, then for any $R>0$,

while if $\nu=\infty$, the statement holds for all $\nu<\infty$.

Conclusions

We have shown that expected improvement can converge near-optimally, but a naive implementation may not converge at all. We thus echo Diaconis and Freedman (1986) in stating that, for infinite-dimensional problems, Bayesian methods are not always guaranteed to find the right answer; such guarantees can only be provided by considering the problem at hand.

We might ask, however, if our framework can also be improved. Our upper bounds on convergence were established using naive algorithms, which in practice would prove inefficient. If a sophisticated algorithm fails where a naive one succeeds, then the sophisticated algorithm is certainly at fault; we might, however, prefer methods of evaluation which do not consider naive algorithms so successful.

Vazquez and Bect (2010) and Grunewalder et al. (2010) consider a more Bayesian formulation of the problem, where the unknown function $f$ is distributed according to the prior $\pi$, but this approach can prove restrictive: as we saw in Section 3.3, placing too much faith in the prior may exclude functions of interest. Further, Grunewalder et al. find the same issues are present also within the Bayesian framework.

A more interesting approach is given by the continuum-armed-bandit problem (Srinivas et al., 2010, and references therein). Here the goal is to minimize the cumulative regret,

in general observing the function $f$ under noise. Algorithms controlling the cumulative regret at rate $r_{n}$ also solve the optimization problem, at rate $r_{n}/n$ (Bubeck et al., 2009, §3). The naive algorithms above, however, have poor cumulative regret. We might, then, consider the cumulative regret to be a better measure of performance, but this approach too has limitations. Firstly, the cumulative regret is necessarily increasing, so cannot establish rates of optimization faster than $n^{-1}$. (This is not an issue under noise, where typically $r_{n}=\Omega(n^{1/2})$, see Kleinberg and Slivkins, 2010.) Secondly, if our goal is optimization, then minimizing the regret, a cost we do not incur, may obscure the problem at hand.

Bubeck et al. (2010) study this problem with the additional assumption that $f$ has finitely many minima, and is, say, quadratic in a neighbourhood of each. This assumption may suffice in practice, and allows the authors to obtain impressive rates of convergence. For optimization, however, a further weakness is that these rates hold only once the algorithm has found a basin of attraction; they thus measure local, rather than global, performance. It may be that convergence rates alone are not sufficient to capture the performance of a global optimization algorithm, and the time taken to find a basin of attraction is more relevant. In any case, the choice of an appropriate framework to measure performance in global optimization merits further study.

Finally, we should also ask how to choose the smoothness parameter $\nu$ (or the equivalent parameter in similar algorithms). van der Vaart and van Zanten (2009) show that Bayesian Gaussian-process models can, in some contexts, automatically adapt to the smoothness of an unknown function $f$. Their technique requires, however, that the estimated length-scales $\hat{\theta}_{n}$ to tend to 0, posing both practical and theoretical challenges. The question of how best to optimize functions of unknown smoothness remains open.

We would like to thank the referees, as well as Richard Nickl and Steffen Grunewalder, for their valuable comments and suggestions.

Appendix A Proofs

and as $\lVert\widehat{K}\rVert_{\infty}\leq\lVert K\rVert_{1}$,

so $\widehat{f}\in L^{1}\cap L^{2}$. $\widehat{f}$ is thus the Fourier transform of a real continuous $f\in L^{2}$, satisfying the Fourier inversion formula everywhere.

If $\nu<\infty$, by assumption $\widehat{K}(\xi)=\widehat{k}(\lVert\xi\rVert)$, for a finite non-increasing function $\widehat{k}$ satisfying $\widehat{k}(\lVert\xi\rVert)=\Theta(\lVert\xi\rVert^{-2\nu-d})$ as $\xi\to\infty$. Hence

If $\nu=\infty$, by a similar argument $\mathcal{H}_{\theta}(\bar{D})$ is continuously embedded in all $H^{s}(D)$. ∎

From 1, we can derive results on the behaviour of $\lVert f\rVert_{\mathcal{H}_{\theta}(S)}$ as $\theta$ varies. For small $\theta$, we obtain the following result.

If $f\in\mathcal{H}_{\theta}(S)$, then $f\in\mathcal{H}_{\theta^{\prime}}(S)$ for all $0<\theta^{\prime}\leq\theta$, and

Let $C=\prod_{i=1}^{d}(\theta^{\prime}_{i}/\theta_{i})$. As $\widehat{K}$ is isotropic and radially non-increasing,

Likewise, for large $\theta$, we obtain the following.

If $\nu<\infty$, $f\in\mathcal{H}_{\theta}(S)$, then $f\in\mathcal{H}_{t\theta}(S)$ for $t\geq 1$, and

for a $C^{\prime\prime}>0$ depending only on $K$ and $\theta$.

As in the proof of 3, we have constants $C,C^{\prime}>0$ such that

and we may argue as in the previous lemma. ∎

We can also describe the posterior distribution of $f$ in terms of $\mathcal{H}_{\theta}(S)$; as a consequence, we may deduce 1.

Suppose $f(x)=\mu+g(x)$, $g\in\mathcal{H}_{\theta}(S)$.

$\hat{f}_{n}(x;\theta)=\hat{\mu}_{n}+\hat{g}_{n}(x)$ solves the optimization problem

with minimum value $\hat{R}_{n}^{2}(\theta)$.

with equality for some $g\in\mathcal{H}_{\theta}(S)$.

Let $W=\operatorname{span}(k_{x_{1}},\dots,k_{x_{n}})$, and write $\hat{g}=\hat{g}^{\parallel}+\hat{g}^{\perp}$ for $\hat{g}^{\parallel}\in W$, $\hat{g}^{\perp}\in W^{\perp}$. $\hat{g}^{\perp}(x_{i})=\langle\hat{g}^{\perp},k_{x_{i}}\rangle=0$, so $\hat{g}^{\perp}$ affects the optimization only through $\lVert\hat{g}\rVert$. The minimal $\hat{g}$ thus has $\hat{g}^{\perp}=0$, so $\hat{g}=\sum_{i=1}^{n}\lambda_{i}k_{x_{i}}$. The problem then becomes

The solution is given by (4) and (5), with value (7).

By symmetry, the prediction error does not depend on $\mu$, so we may take $\mu=0$. Then

for $e_{n,x}=k_{x}-\sum_{i=1}^{n}\lambda_{i}k_{x_{i}}$, and

Now, $\lVert e_{n,x}\rVert^{2}_{\mathcal{H}_{\theta}(S)}=s_{n}^{2}(x;\theta)$, as given by (6); this is a consequence of Loève’s isometry, but is easily verified algebraically. The result then follows by Cauchy-Schwarz. ∎

A.2 Fixed Parameters

We first establish the lower bound. Suppose we have $2n$ functions $\psi_{m}$ with disjoint supports. We will argue that, given $n$ observations, we cannot distinguish between all the $\psi_{m}$, and thus cannot accurately pick a minimum $x_{n}^{*}$.

but on that event, $u$ cannot distinguish between and $\psi_{m}$ before time $n$, so

As the minimax loss is non-increasing in $n$, for $(2(k-1))^{d}/2\leq n<(2k)^{d}/2$ we conclude

For general $X$ having non-empty interior, we can find a hypercube $S=x_{0}+[0,\varepsilon]^{d}\subseteq X$, with $\varepsilon>0$. We may then proceed as above, picking functions $\psi_{m}$ supported inside $S$.

For the upper bound, consider a strategy $u$ choosing a fixed sequence $x_{n}$, independent of the $z_{n}$. Fit a radial basis function interpolant $\hat{f}_{n}$ to the data, and pick $x_{n}^{*}$ to minimize $\hat{f}_{n}$. Then if $x^{*}$ minimizes $f$,

so the loss is bounded by the error in $\hat{f}_{n}$.

From results in Narcowich et al. (2003, §6) and Wendland (2005, §11.5), for suitable radial basis functions the error is uniformly bounded by

To prove Theorem 2, we first show that some observations $z_{n}$ will be well-predicted by past data.

If $\nu\leq\frac{1}{2}$, then by assumption

as $x\to 0$. If $\nu>\frac{1}{2}$, then $K$ is differentiable, so as $K$ is symmetric, $\nabla K(0)=0$. If further $\nu\leq 1$, then

Similarly, if $\nu>1$, then $K$ is $C^{2}$, so

for a constant $C>0$ depending only on $X$, $K$ and $\theta$.

We next show that most design points $x_{n+1}$ are close to a previous $x_{i}$. $X$ is bounded, so can be covered by $k$ balls of radius $O(k^{-1/d})$. If $x_{n+1}$ lies in a ball containing some earlier point $x_{i}$, $i\leq n$, then we may conclude

for a constant $C^{\prime}>0$ depending only on $X$, $K$ and $\theta$. Hence as there are $k$ balls, at most $k$ points $x_{n+1}$ can satisfy

Next, we provide bounds on the expected improvement when $f$ lies in the RKHS.

If $s=0$, then by 6, $\hat{f}_{n}(x;\theta)=f(x)$, so $EI_{n}(x;\pi)=I$, and the result is trivial. Suppose $s>0$, and set $t=(f(x_{n}^{*})-f(x))/s$, $u=(f(x_{n}^{*})-\hat{f}_{n}(x;\theta))/s$. From (8) and (9),

and by 6, $\lvert u-t\rvert\leq R$. As $\tau^{\prime}(z)=\Phi(z)\in$, $\tau$ is non-decreasing, and $\tau(z)\leq 1+z$ for $z\geq 0$. Hence

If $I=0$, then as $EI$ is the expectation of a non-negative quantity, $EI\geq 0$, and the lower bounds are trivial. Suppose $I>0$. Then as $EI\geq 0$, $\tau(z)\geq 0$ for all $z$, and $\tau(z)=z+\tau(-z)\geq z$. Thus

Combining these bounds, and eliminating $s$, we obtain

We may now prove the theorem. We will use the above bounds to show that there must be times $n_{k}$ when the expected improvement is low, and thus $f(x_{n_{k}}^{*})$ is close to $\min f$.

holds at most $k$ times. Furthermore, $z_{n}^{*}-z_{n+1}^{*}\geq 0$, and for $\lVert f\rVert_{\mathcal{H}_{\theta}(X)}\leq R$,

so $z_{n}^{*}-z_{n+1}^{*}>2Rk^{-1}$ at most $k$ times. Since $z_{n}^{*}-f(x_{n+1})\leq z_{n}^{*}-z_{n+1}^{*}$, we have also $z_{n}^{*}-f(x_{n+1})>2Rk^{-1}$ at most $k$ times. Thus there is a time $n_{k}$, $k\leq n_{k}\leq 3k$, for which $s_{n_{k}}(x_{n_{k}+1};\theta)\leq Ck^{-(\nu\wedge 1)/d}(\log k)^{\beta}$ and $z_{n_{k}}^{*}-f(x_{n_{k}+1})\leq 2Rk^{-1}$.

Let $f$ have minimum $z^{*}$ at $x^{*}$. For $k$ large, $x_{n_{k}+1}$ will have been chosen by expected improvement (rather than being an initial design point, chosen at random). Then as $z_{n}^{*}$ is non-increasing in $n$, for $3k\leq n<3(k+1)$ we have by 8,

This bound is uniform in $f$ with $\lVert f\rVert_{\mathcal{H}_{\theta}(X)}\leq R$, so we obtain

A.3 Estimated Parameters

To prove Theorem 3, we first establish lower bounds on the posterior variance.

uniformly in the sequences $x_{n}$, $\theta_{n}$.

Given $n$ design points $x_{1},\dots,x_{n}$, there must be some $\psi_{m}$ such that $\psi_{m}(x_{i})=0$, $1\leq i\leq n$. By 6, the posterior mean of $\psi_{m}$ given these observations is the zero function. Thus for $x\in T$ minimizing $\psi_{m}$,

As $s_{n}(x;\theta)$ is non-increasing in $n$, for $\frac{1}{2}(2(k-1))^{d}<n\leq\frac{1}{2}(2k)^{d}$ we obtain

Next, we bound the expected improvement when prior parameters are estimated by maximum likelihood.

Let $\lVert f\rVert_{\mathcal{H}_{\theta^{U}}(X)}\leq R$, $x_{n},y_{n}\in X$. Set $I_{n}(x)=z_{n}^{*}-f(x)$, $s_{n}(x)=s_{n}(x;\hat{\theta}_{n})$, and $t_{n}(x)=I_{n}(x)/s_{n}(x)$. Suppose:

for some $T_{n}\to-\infty$, $t_{n}(x_{n+1})\leq T_{n}$ whenever $s_{n}(x_{n+1})>0$;

for some $C>0$, $s_{n}(y_{n+1})\geq e^{-C/c_{n}}$.

Then for $\hat{\pi}_{n}$ as in 2, eventually $EI_{n}(x_{n+1};\hat{\pi}_{n})<EI_{n}(y_{n+1};\hat{\pi}_{n})$. If the conditions hold on a subsequence, so does the conclusion.

Let $\hat{R}_{n}^{2}(\theta)$ be given by (7), and set $\hat{R}_{n}^{2}=\hat{R}_{n}^{2}(\hat{\theta}_{n})$. For $n\geq j$, $\hat{R}_{n}^{2}>0$, and by 4 and 1,

Thus $0<\hat{\sigma}^{2}_{n}\leq S^{2}c_{n}$. Then if $s_{n}(x)>0$, for some $\lvert u_{n}(x)-t_{n}(x)\rvert\leq S$,

If $s_{n}(x_{n+1})=0$, then $x_{n+1}\in\{x_{1},\dots,x_{n}\}$, so

When $s_{n}(x_{n+1})>0$, as $\tau$ is increasing we may upper bound $EI_{n}(x_{n+1};\hat{\pi}_{n})$ using $u_{n}(x_{n+1})\leq T_{n}+S$, and lower bound $EI_{n}(y_{n+1};\hat{\pi}_{n})$ using $u_{n}(y_{n+1})\geq-S$. Since $s_{n}(x_{n+1})\leq 1$, and $\tau(x)=\Theta(x^{-2}e^{-x^{2}/2})$ as $x\to-\infty$ (Abramowitz and Stegun, 1965, §7.1),

If the conditions hold on a subsequence, we may similarly argue along that subsequence. ∎

Finally, we will require the following technical lemma.

Given $\varepsilon>0$, fix $m\geq n/\varepsilon$, and pick disjoint open sets $U_{1},\dots,U_{m}\subset S$. Then

We may now prove the theorem. We will construct a function $f$ on which the $EI(\hat{\pi})$ strategy never observes within a region $W$. We may then construct a function $g$, agreeing with $f$ except on $W$, but having different minimum. As the strategy cannot distinguish between $f$ and $g$, it cannot successfully find the minimum of both.

We work conditional on the event $A$, having probability at least $1-\varepsilon$, that $z^{*}_{k}=0$, and thus $z^{*}_{n}=0$ for all $n\geq k$. Suppose $x_{n}\in V_{1}$ infinitely often, so the $z_{n}$ are not all equal. By 7, $s_{n}(x_{n+1};\hat{\theta}_{n})\to 0$, so on a subsequence with $x_{n+1}\in V_{1}$, we have

whenever $s_{n}(x_{n+1};\hat{\theta}_{n})>0$. However, by 9, there are points $y_{n}\in V_{0}$ with $z_{n}^{*}-f(y_{n+1})=0$, and $s_{n}(y_{n+1};\hat{\theta}_{n})=\Omega(n^{-\nu/d})$. Hence by 10, $EI_{n}(x_{n+1};\hat{\pi}_{n})<EI_{n}(y_{n+1};\hat{\pi}_{n})$ for some $n$, contradicting the definition of $x_{n+1}$.

Construct a smooth function $g$ by adding to $f$ a $C^{\infty}$ function which is 0 outside $W$, and has minimum $-2$. Then $\min g=-1$, but on the event $C$, $EI(\hat{\pi})$ cannot distinguish between $f$ and $g$, and $g(x_{n}^{*})\geq 0$. Thus for $\delta=1$,

As the behaviour of $EI(\hat{\pi})$ is invariant under rescaling, we may scale $g$ to have norm $\lVert g\rVert_{\mathcal{H}_{\theta}(X)}\leq R$, and the above remains true for some $\delta>0$. ∎

As in the proof of Theorem 2, we will show there are times $n_{k}$ when the expected improvement is small, so $f(x_{n_{k}})$ must be close to the minimum. First, however, we must control the estimated parameters $\hat{\sigma}^{2}_{n}$, $\hat{\theta}_{n}$.

If the $z_{n}$ are all equal, then by assumption the $x_{n}$ are dense in $X$, so $f$ is constant, and the result is trivial. Suppose the $z_{n}$ are not all equal, and let $T$ be a random variable satisfying $z_{T}\neq z_{i}$ for some $i<T$. Set $U=\inf_{\theta^{L}\leq\theta\leq\theta^{U}}\hat{R}_{T}(\theta)$. $\hat{R}_{T}(\theta)$ is a continuous positive function, so $U>0$. Let $S^{2}=R^{2}\prod_{i=1}^{d}(\theta^{U}_{i}/\theta^{L}_{i})$. By 4, $\lVert f\rVert_{\mathcal{H}_{\hat{\theta}_{n}}(X)}\leq S$, so by 1, for $n\geq T$,

As in the proof of Theorem 2, we have a constant $C>0$, and some $n_{k}$, $k\leq n_{k}\leq 3k$, for which $z_{n_{k}}^{*}-f(x_{n_{k}+1})\leq 2Rk^{-1}$ and $s_{n_{k}}(x_{n_{k}+1};\hat{\theta}_{n_{k}})\leq Ck^{-\alpha}(\log k)^{\beta}$. Then for $k\geq T$, $3k\leq n<3(k+1)$, arguing as in Theorem 2 we obtain

We thus have a random variable $C^{\prime}$ satisfying $z_{n}^{*}-z^{*}\leq C^{\prime}n^{-(\nu\wedge 1)/d}(\log n)^{\beta}$ for all $n$, and the result follows. ∎

A.4 Near-Optimal Rates

To prove Theorem 5, we first show that the points chosen at random will be quasi-uniform in $X$.

Let $x_{n}$ be i.i.d. random variables, distributed uniformly over $X$, and define their mesh norm,

For any $\gamma>0$, there exists $C>0$ such that

We will partition $X$ into $n$ regions of size $O(n^{-1/d})$, and show that with high probability we will place an $x_{i}$ in each one. Then every point $x$ will be close to an $x_{i}$, and the mesh norm will be small.

For $n$ large, $\mu_{n}\leq 1$, so by the generalized Chernoff bound of Panconesi and Srinivasan (1997, §3.1),

On the event $\sum_{m}I_{m}<1$, $I_{m}=0$ for all $m$. For any $x\in X$, we then have $x\in X_{m}$ for some $m$, and $x_{j}\in X_{m}$ for some $1\leq j\leq\lfloor\gamma n\log n\rfloor$. Thus

As this bound is uniform in $x$, we obtain $h_{\lfloor\gamma n\log n\rfloor}\leq\sqrt{d}k^{-1}$. Thus for $n=k^{d}$,

and as $h_{n}$ is non-increasing in $n$, this bound holds also for $k^{d}\leq n<(k+1)^{d}$. By a change of variables, we then obtain

and the result follows by choosing $\gamma$ large. For general $X$, as $X$ is bounded it can be partitioned into $n$ regions of measure $\Theta(n^{-1/d})$, so we may argue similarly. ∎

We may now prove the theorem. We will show that the points $x_{n}$ must be quasi-uniform in $X$, so posterior variances must be small. Then, as in the proofs of Theorems 2 and 4, we have times when the expected improvement is small, so $f(x_{n}^{*})$ is close to $\min f$.

First suppose $\nu<\infty$. Let the $EI(\,\cdot\,,\varepsilon)$ choose $k$ initial design points independent of $f$, and suppose $n\geq 2k$. Let $A_{n}$ be the event that $\lfloor\frac{\varepsilon}{4}n\rfloor$ of the points $x_{k+1},\dots,x_{n}$ are chosen uniformly at random, so by a Chernoff bound,

Let $B_{n}$ be the event that one of the points $x_{n+1},\dots,x_{2n}$ is chosen by expected improvement, so

Finally, let $C_{n}$ be the event that $A_{n}$ and $B_{n}$ occur, and further the mesh norm $h_{n}\leq C(n/\log n)^{-1/d}$, for the constant $C$ from 12. Set $r_{n}=(n/\log n)^{-\nu/d}(\log n)^{\alpha}$. Then by 12, since $C_{n}\subset A_{n}$,

for a constant $C^{\prime}>0$ not depending on $f$.

Let $EI(\,\cdot\,,\varepsilon)$ have prior $\pi_{n}$ at time $n$, with (fixed or estimated) parameters $\sigma_{n}$, $\theta_{n}$. Suppose $\lVert f\rVert_{\mathcal{H}_{\theta^{U}}(X)}\leq R$, and set $S^{2}=R^{2}\prod_{i=1}^{d}(\theta^{U}_{i}/\theta^{L}_{i})$, so by 4, $\lVert f\rVert_{\mathcal{H}_{\theta_{n}}(X)}\leq S$. If $\alpha=0$, then by Narcowich et al. (2003, §6),

uniformly in $\theta$, for $M(\theta)$ a continuous function of $\theta$. Hence on the event $C_{n}$,

for a constant $C^{\prime\prime}>0$ depending only on $X$, $K$, $C$, $\theta^{L}$ and $\theta^{U}$. If $\alpha>0$, the same result holds by a similar argument.

On the event $C_{n}$, we have some $x_{m}$ chosen by expected improvement, $n<m\leq 2n$. Let $f$ have minimum $z^{*}$ at $x^{*}$. Then by 8,

for a constant $T>0$. (Under $EI(\pi,\varepsilon)$, we have $T=2S+\sigma$; otherwise $\sigma_{m-1}\leq S$ by 1, so $T=3S$.) Thus, rearranging,

On the event $C_{n}^{c}$, we have $z_{2n}^{*}-z^{*}\leq 2\lVert f\rVert_{\infty}\leq 2R$, so

As this bound is uniform in $f$ with $\lVert f\rVert_{\mathcal{H}_{\theta^{U}}(X)}\leq R$, the result follows. If instead $\nu=\infty$, the above argument holds for any $\nu<\infty$. ∎

References