Convergence rates of efficient global optimization algorithms

1 Introduction

Suppose we wish to minimize a continuous function $f : X \to R$ , where $X$ is a compact subset of $R^{d}$ . Observing $f (x)$ is costly (it may require a lengthy computer simulation or physical experiment), so we wish to use as few observations as possible. We know little about the shape of $f$ ; in particular we will be unable to make assumptions of convexity or unimodality. We therefore need a global optimization algorithm, one which attempts to find a global minimum.

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 {^f}_{n}$ , for some approximation ${^f}_{n}$ to $f$ . We will show that if ${^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).

Expected improvement performs well in experiments (Osborne, 2010, §9.5), but little is known about its theoretical properties. The behaviour of the algorithm depends crucially on the Gaussian process prior $π$ chosen for $f$ . Each prior has an associated space of functions $H$ , its reproducing-kernel Hilbert space. $H$ contains all functions $X \to R$ as smooth as a posterior mean of $f$ , and is the natural space in which to study questions of convergence.

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

Table 1 summarizes the notation used in this paper. We say $f : R^{d} \to R$ is a bump function if $f$ is infinitely differentiable and of compact support, and $f : R^{d} \to C$ is Hermitian if $¯ ¯¯¯¯¯¯¯¯ ¯ f (x) = f (- x)$ . We use the Landau notation $f = O (g)$ to denote $lim sup | f / g | < \infty$ , and $f = o (g)$ to denote $f / g \to 0$ . If $g = O (f)$ , we say $f = Ω (g)$ , and if both $f = O (g)$ and $f = Ω (g)$ , we say $f = Θ (g)$ . If further $f / g \to 1$ , we say $f \sim g$ . Finally, if $f$ and $g$ are random, and $P (sup | f / g | \leq M) \to 1$ as $M \to \infty$ , we say $f = O_{p} (g)$ .

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.

Section 1
$f$	unknown function $X \to R$ to be minimized
$X$	compact subset of $R^{d}$ to minimize over
$d$	number of dimensions to minimize over
$x_{n}$	points in $X$ at which $f$ is observed
$z_{n}$	observations $z_{n} = f (x_{n})$ of $f$
$x_{n}^{*}$	estimated minimum of $f$ , given $z_{1}, \dots, z_{n}$
Section 2.1
$π$	prior distribution for $f$
$u$	strategy for choosing $x_{n}$ , $x_{n}^{*}$
$F_{n}$	filtration $F_{n} = σ (x_{i}, z_{i} : i \leq n)$
$z_{n}^{*}$	best observation $z_{n}^{*} = {min}_{i = 1, \dots, n} z_{i}$
$E I_{n}$	expected improvement given $F_{n}$
Section 2.2
$μ$ , $σ^{2}$	global mean and variance of Gaussian-process prior $π$
$K$	underlying correlation kernel for $π$
$K_{θ}$	correlation kernel for $π$ with length-scales $θ$
$ν$ , $α$	smoothness parameters of $K$
${^μ}_{n}$ , ${^f}_{n}$ , $s_{n}^{2}$ , ${^R}_{n}^{2}$	quantities describing posterior distribution of $f$ given $F_{n}$
Section 2.3
$E I (π)$	expected improvement strategy with fixed prior
${^σ}_{n}^{2}$ , ${^θ}_{n}$	estimates of prior parameters $σ^{2}$ , $θ$
$c_{n}$	rate of decay of ${^σ}_{n}^{2}$
$θ^{L}$ , $θ^{U}$	bounds on ${^θ}_{n}$
$E I (^π)$	expected improvement strategy with estimated prior
Section 3.1
$H_{θ} (S)$	reproducing-kernel Hilbert space of $K_{θ}$ on $S$
$H^{s} (D)$	Sobolev Hilbert space of order $s$ on $D$
Section 3.2
$L_{n}$	loss suffered over an RKHS ball after $n$ steps
Section 3.3
$E I (~ π)$	expected improvement strategy with robust estimated prior
Section 3.4
$E I (\cdot, ε)$	$ε$ -greedy expected improvement strategies

Table 1: Notation

2 Expected Improvement

Suppose we wish to minimize an unknown function $f$ , choosing design points $x_{n}$ and estimated minima $x_{n}^{*}$ as in the introduction. If we pick a prior distribution $π$ for $f$ , representing our beliefs about the unknown function, we can describe this problem in terms of decision theory. Let $(Ω, F, P)$

be a probability space, equipped with a random process

f

having law

π

. A strategy

u

is a collection of random variables

(x_{n})

(x_{n}^{*})

taking values in

X

. Set

zn\coloneqqf(xn)

, and define the filtration

Fn\coloneqqσ(xi,zi:i≤n)

. The strategy

u

is valid if

x_{n}

is conditionally independent of

f

given

F_{n - 1}

, and likewise

x_{n}^{*}

given

F_{n}

. (Note that we allow random strategies, provided they do not depend on unknown information about

f

When taking probabilities and expectations we will write $P_{π}^{u}$ and $E_{π}^{u}$ , denoting the dependence on both the prior $π$ and strategy $u$ . The average-case performance at some future time $N$ is then given by the expected loss,

E_{π}^{u} [f (x_{N}^{*}) - min f],

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

2.1 Bayesian Optimization

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 $F_{n}$ , if we are to stop at time $n$ we should choose $x∗n\coloneqqxi∗n$ , where $i∗n\coloneqqargmin1,…,nzi$ . (In the case of ties, we may pick any minimizing $i_{n}^{*}$ .) We then suffer a loss $z_{n}^{*} - min f$ , where $z∗n\coloneqqzi∗n$ . Were we to observe at $x_{n + 1}$ before stopping, the expected loss would be

E_{π}^{u} [z_{n + 1}^{*} - min f ∣ F_{n}],

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

EIn(xn+1;π)\coloneqqEuπ[z∗n−z∗n+1∣Fn]=Euπ[(z∗n−zn+1)+∣Fn],

(1)

where $x^{+} = max (x, 0)$ .

So far, we have merely replaced one optimization problem with another. However, for suitable priors, $E I_{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.2 Gaussian Process Models

We still need to choose a prior $π$ 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

E_{π} [f (x)] = μ, C {o v}_{π} [f (x), f (y)] = σ^{2} K_{θ} (x - y) .

(2)

$μ \in R$ is the global mean of $f$ ; we place a flat prior on $μ$ , reflecting our uncertainty over the location of $f$ .

$σ > 0$ is the global scale of variation of $f$ , and $K_{θ} : R^{d} \to R$ its correlation kernel, governing the local properties of $f$ . In the following, we will consider kernels

Kθ(t1,…,td)\coloneqqK(t1/θ1,…,td/θd),

(3)

for an underlying kernel $K$ with $K (0) = 1$ . (Note that we can always satisfy this condition by suitably scaling $K$ and $σ$ .) The $θ_{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 $θ_{i}$ . For now, we will assume the parameters $σ$ and $θ$ 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.

Assumption 1.

$K$ is continuous and integrable.

$K$

thus has Fourier transform

ˆK(ξ)\coloneqq∫Rde−2πi⟨x,ξ⟩K(x)dx,

and by Bochner’s theorem, $ˆ K$ is non-negative and integrable.

Assumption 2.

$ˆ K$ is isotropic and radially non-increasing.

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

Assumption 3.

As $x \to \infty$ , either:

$ˆ K (x) = Θ (∥ x ∥^{- 2 ν - d})$ for some $ν > 0$ ; or
$ˆ K (x) = O (∥ x ∥^{- 2 ν - d})$ for all $ν > 0$ (we will then say that $ν = \infty$ ).

Note the condition $ν > 0$ is required for $ˆ K$ to be integrable.

Assumption 4.

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

| K (x) - P_{k} (x) | = O (∥ x ∥^{2 ν} (- log ∥ x ∥)^{2 α})

as $x \to 0$ , for some $α \geq 0$ .

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

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

Kν(x)\coloneqq21−νΓ(ν)(√2ν∥x∥)νkν(√2ν∥x∥),ν∈(0,∞),

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

K∞(x)\coloneqqe−12∥x∥2,

obtained in the limit $ν \to \infty$ (Rasmussen and Williams, 2006, §4.2). Between them, these kernels cover the full range of smoothness $0 < ν \leq \infty$ . Both kernels satisfy Assumptions 1–4 for the $ν$ given; $α = 0$ except for the Matérn kernel with $ν \in N$ , where $α = \frac{1}{2}$ (Abramowitz and Stegun, 1965, §9.6).

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

f (x) ∣ z_{1}, \dots, z_{n} \sim N ({^f}_{n} (x; θ), σ^{2} s_{n}^{2} (x; θ)),

where

	${^μ}_{n} (θ)$	$\coloneqq1TV−1z1TV−11,$	(4)
	${^f}_{n} (x; θ)$	$\coloneqq^μn+vTV−1(z−^μn1),$	(5)
and
	$s_{n}^{2} (x; θ)$	$\coloneqq1−vTV−1v+(1−1TV−1v)21TV−11,$	(7)

for $z = (z_{i})_{i = 1}^{n}$ , $V = (K_{θ} (x_{i} - x_{j}))_{i, j = 1}^{n}$ , and $v = (K_{θ} (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,

^R2n(θ)\coloneqq(z−^μn1)TV−1(z−^μn1).

(8)

2.3 Expected Improvement Strategies

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

E I_{n} (x_{n + 1}; π) = ρ (z_{n}^{*} - {^f}_{n} (x_{n + 1}; θ), σ s_{n} (x_{n + 1}; θ)),

(9)

where

ρ(y,s)\coloneqq{yΦ(y/s)+sφ(y/s),s>0,max(y,0),s=0,

(10)

and $Φ$ and $φ$

are the standard normal distribution and density functions respectively.

For a prior $π$ as above, expected improvement chooses $x_{n + 1}$ to maximize (9), 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 $E I_{n}$ . In general, this will not affect the behaviour of the algorithm, so we allow any choice of $x_{n + 1}$ maximizing (9).

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.

Definition 1.

An $E I (π)$ strategy chooses:

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 (9).

So far, we have not considered the choice of parameters $σ$ and $θ$ . 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 $θ$ , the MLE ${^σ}_{n}^{2} = {^R}_{n}^{2} (θ) / n$ ; for full generality, we will allow any choice ${^σ}_{n}^{2} = c_{n} {^R}_{n}^{2} (θ)$ , where $c_{n} = o (1 / log n)$ . Estimates of $θ$ , however, must be obtained by numerical optimization. As $θ$ can vary widely in scale, this optimization is best performed over $log θ$ ; 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 θ$ . We have thus described the following strategy.

Definition 2.

Let ${^π}_{n}$ be a sequence of priors, with parameters ${^σ}_{n}$ , ${^θ}_{n}$ satisfying:

${^σ}_{n}^{2} = c_{n} {^R}_{n}^{2} ({^θ}_{n})$ for constants $c_{n} > 0$ , $c_{n} = o (1 / log n)$ ; and
$θ^{L} \leq {^θ}_{n} \leq θ^{U}$ for constants $θ^{L}$ , $θ^{U} \in R_{+}^{d}$ .

An $E I (^π)$ strategy satisfies 1, replacing $π$ with ${^π}_{n}$ in (9).

3 Convergence Rates

To discuss convergence, we must first choose a smoothness class for the unknown function $f$ . Each kernel $K_{θ}$ is associated with a space of functions $H_{θ} (X)$ , its reproducing-kernel Hilbert space (RKHS) or native space. $H_{θ} (X)$ contains all functions $X \to R$ as smooth as a posterior mean of $f$ , and is the natural space to study convergence of expected-improvement algorithms, allowing a tractable analysis of their asymptotic behaviour.

3.1 Reproducing-Kernel Hilbert Spaces

Given a symmetric positive-definite kernel $K$ on $R^{d}$ , set $k_{x} (t) = K (t - x)$ . For $S \subseteq R^{d}$ , let $E (S)$ be the space of functions $S \to R$ spanned by the $k_{x}$ , for $x \in S$ . Furnish $E (S)$ with the inner product defined by

⟨kx,ky⟩\coloneqqK(x−y).

The completion of $E (S)$ under this inner product is the reproducing-kernel Hilbert space $H (S)$ of $K$ on $S$ . The members $f \in H (S)$ are abstract objects, but we can identify them with functions $f : S \to R$ through the reproducing property,

f (x) = ⟨ f, k_{x} ⟩,

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

We will find it convenient also to use an alternative characterization of $H (S)$ . We begin by describing $H (R^{d})$ in terms of Fourier transforms. Let $ˆ f$ denote the Fourier transform of a function $f \in L^{2}$ . The following result is stated in Parzen (1963, §2), and proved in Wendland (2005, §10.2); we give a short proof in Appendix A.

Lemma 1.

$H (R^{d})$ is the space of real continuous $f \in L^{2} (R^{d})$ whose norm

∥f∥2H(Rd)\coloneqq∫|ˆf(ξ)|2ˆK(ξ)dξ

is finite, taking $0 / 0 = 0$ .

We may now describe $H (S)$ in terms of $H (R^{d})$ .

Lemma 2 (Aronszajn, 1950, §1.5).

$H (S)$ is the space of functions $f = g |_{S}$ for some $g \in H (R^{d})$ , with norm

∥f∥H(S)\coloneqqinfg|S=f∥g∥H(Rd),

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

These spaces are in fact closely related to the Sobolev Hilbert spaces of functional analysis. Say a domain $D \subseteq R^{d}$ is Lipschitz if its boundary is locally the graph of a Lipschitz function (see Tartar, 2007, §12, for a precise definition). For such a domain $D$ , the Sobolev Hilbert space $H^{s} (D)$ is the space of functions $f : D \to R$ , given by the restriction of some $g : R^{d} \to R$ , whose norm

∥f∥2Hs(D)\coloneqqinfg|D=f∫|ˆg(ξ)|2(1+∥ξ∥2)s/2dξ

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

C ∥ f ∥_{1} \leq ∥ f ∥_{2} \leq C^{'} ∥ f ∥_{1},

for constants $0 < C \leq C^{'}$ .

Lemma 3.

Let $H_{θ} (S)$ denote the RKHS of $K_{θ}$ on $S$ , and $D \subseteq R^{d}$ be a Lipschitz domain.

If $ν < \infty$ , $H_{θ} (¯ D)$ is equivalent to the Sobolev Hilbert space $H^{ν + d / 2} (D)$ .
If $ν = \infty$ , $H_{θ} (¯ D)$ is continuously embedded in $H^{s} (D)$ for all $s$ .

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

3.2 Fixed Parameters

We are now ready to state our main results. Let $X \subset R^{d}$ be compact with non-empty interior. For a function $f : X \to R$ , let $P_{f}^{u}$ and $E_{f}^{u}$ denote probability and expectation when minimizing the fixed function $f$ with strategy $u$ . (Note that while $f$ is fixed, $u$ may be random, so its performance is still probabilistic in nature.) We define the loss suffered over the ball $B_{R}$ in $H_{θ} (X)$ after $n$ steps by a strategy $u$ ,

Ln(u,Hθ(X),R)\coloneqqsup∥f∥Hθ(X)≤REuf[f(x∗n)−minf].

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

L_{n} (u, H_{θ} (X), R) = O (r_{n})

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 $∥ f ∥_{H_{θ} (X)}$ .

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

Theorem 1.

If $ν < \infty$ , then for any $θ \in R_{+}^{d}$ , $R > 0$ ,

inf u L_{n} (u, H_{θ} (X), R) = Θ (n^{- ν / 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 $E I (π)$ strategy of 1. We can show this strategy converges at least at rate $n^{- (ν \land 1) / d}$ , up to log factors.

Theorem 2.

Let $π$ be a prior with length-scales $θ \in R_{+}^{d}$ . For any $R > 0$ ,

L_{n} (E I (π), H_{θ} (X), R) = {\begin{matrix} O (n^{- ν / d} (log n)^{α}), & ν \leq 1, O (n^{- 1 / d}), & ν > 1. \end{matrix}

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

3.3 Estimated Parameters

First, we consider the effect of the prior parameters on $E I (π)$ . 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 $E I (^π)$ strategy, given by 2, uses maximum likelihood estimates for this purpose. We can show, however, that this may cause the strategy to never converge.

Theorem 3.

Suppose $ν < \infty$ . Given $θ \in R_{+}^{d}$ , $R > 0$ , $ε > 0$ , there exists $f \in H_{θ} (X)$ satisfying $∥ f ∥_{H_{θ} (X)} \leq R$ , and for some fixed $δ > 0$ ,

P_{f}^{E I (^π)} (inf n f (x_{n}^{*}) - min f \geq δ) > 1 - ε .

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.

Figure 1: A counterexample from Theorem 3

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

τ(x)\coloneqqxΦ(x)+φ(x).

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

C (R + σ) \frac{τ (R / σ)}{τ (- R / σ)},

(11)

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

Corollary 1.

${^R}_{n} (θ)$ is non-decreasing in $n$ , and bounded above by $∥ f ∥_{H_{θ} (X)}$ .

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

To resolve the issue, we will instead try to pick $σ$ to minimize (11). The term $R + σ$ is increasing in $σ$ , and the term $τ (R / σ) / τ (- R / σ)$ is decreasing in $σ$ ; we may balance the terms by taking $σ = R$ . The constant is then proportional to $R$ , which we may minimize by taking $R = ∥ f ∥_{H_{θ} (X)}$ . In practice, we will not know $∥ f ∥_{H_{θ} (X)}$ in advance, so we must estimate it from the data; from 1, a convenient estimate is ${^R}_{n} (θ)$ .

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

If we believe our Gaussian-process model, this estimate ${^σ}_{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}$ , $E I_{n} (\cdot; {^π}_{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.

Definition 3.

An $E I (~ π)$ strategy satisfies 2, except:

we instead set ${^σ}_{n}^{2} = {^R}_{n}^{2} ({^θ}_{n})$ ; and
we require the choice of $x_{n + 1}$ maximizing (9) 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 $H_{θ} (X)$ , as the rate of convergence depends on the ratio $R / {^R}_{n}$ , which is unbounded. (Indeed, any estimator of $∥ f ∥_{H_{θ} (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.

Theorem 4.

For any $f \in H_{θ^{U}} (X)$ , under $P_{f}^{E I (~ π)}$ ,

f (x_{n}^{*}) - min f = {\begin{matrix} O_{p} (n^{- ν / d} (log n)^{α}), & ν \leq 1, O_{p} (n^{- 1 / d}), & ν > 1. \end{matrix}

3.4 Near-Optimal Rates

So far, our rates have been near-optimal only for $ν \leq 1$ . To obtain good rates for $ν > 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 $ε$ -greedy strategy (Sutton and Barto, 1998, §2.2). In such a strategy, at each step with probability $1 - ε$ we make a decision to maximize some greedy criterion; with probability $ε$ 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 $ε$ controls the trade-off between global and local search: a good choice of $ε$ 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 $ε = 0.1$ and $ε = 0.01$ , but in practical work $ε$ should of course be calibrated to a typical problem set.

We therefore define the following strategies.

Definition 4.

Let $\cdot$ denote $π$ , $^π$ or $~ π$ . For $0 < ε < 1$ , an $E I (\cdot, ε)$ strategy:

chooses initial design points $x_{1}, \dots, x_{k}$ independently of $f$ ;
with probability $1 - ε$ , chooses design point $x_{n + 1} (n \geq k)$ as in $E I (\cdot)$ ; or
with probability $ε$ , 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 $ν < \infty$ .

Theorem 5.

Let $E I (\cdot, ε)$ be one of the strategies in 4. If $ν < \infty$ , then for any $R > 0$ ,

L_{n} (E I (\cdot, ε), H_{θ^{U}} (X), R) = O ((n / log n)^{- ν / d} (log n)^{α}),

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

Note that unlike a typical $ε$ -greedy algorithm, we do not rely on random choice to obtain global convergence: as above, the $E I (π)$ and $E I (~ π)$ strategies are already globally convergent. Instead, we use random choice simply to improve upon the worst-case rate. Note also that the result does not in general hold when $ε = 1$ ; to obtain good rates, we must combine global search with inference about $f$ .

4 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 $π$ , 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,

Rn\coloneqqn∑i=1(f(xi)−minf),

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} = Ω (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 $ν$ (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 ${^θ}_{n}$ to tend to 0, posing both practical and theoretical challenges. The question of how best to optimize functions of unknown smoothness remains open.

Acknowledgements

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

Appendix A Proofs

We now prove the results in Section 3.

a.1 Reproducing-Kernel Hilbert Spaces

Proof of 1.

Let $V$ be the space of functions described, and $W$ be the closed real subspace of Hermitian functions in $L^{2} (R^{d}, {ˆ K}^{- 1})$ . We will show $f \mapsto ˆ f$ is an isomorphism $V \to W$ , so we may equivalently work with $W$ . Given $ˆ f \in W$ , by Cauchy-Schwarz and Bochner’s theorem,

\int | ˆ f | \leq {(\int ˆ K)}^{1 / 2} {(\int | ˆ f |^{2} / ˆ K)}^{1 / 2} < \infty,

and as $∥ ˆ K ∥_{\infty} \leq ∥ K ∥_{1}$ ,

\int | ˆ f |^{2} \leq ∥ ˆ K ∥_{\infty} \int | ˆ f |^{2} / ˆ K < \infty,

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

$f \mapsto ˆ f$ is hence an isomorphism $V \to W$ . It remains to show that $V = H (R^{d})$ . $W$ is complete, so $V$ is. Further, $E (R^{d}) \subset V$ , and by Fourier inversion each $f \in V$ satisfies the reproducing property,

f (x) = \int e^{2 π i ⟨ x, ξ ⟩} ˆ f (ξ) d ξ = \int \frac{ˆ f (ξ) ¯ ¯¯¯¯¯¯¯¯¯¯ ¯ {ˆ k}_{x} (ξ)}{ˆ K (ξ)} d ξ = ⟨ f, k_{x} ⟩,

so $H (R^{d})$ is a closed subspace of $V$ . Given $f \in H (R^{d})^{⊥}$ , $f (x) = ⟨ f, k_{x} ⟩ = 0$ for all $x$ , so $f = 0$ . Thus $V = H (R^{d})$ . ∎

Proof of 3.

By 1, the norm on $H_{θ} (R^{d})$ is

∥ f ∥_{H_{θ} (R^{d})}^{2} = \int \frac{| ˆ f (ξ) |^{2}}{{ˆ K}_{θ} (ξ)} d ξ,

and $K_{θ}$ has Fourier transform

{ˆ K}_{θ} (ξ) = \frac{ˆ K (ξ_{1} / θ_{1}, \dots, ξ_{d} / θ_{d})}{\prod_{i = 1}^{d} θ_{i}} .

If $ν < \infty$ , by assumption $ˆ K (ξ) = ˆ k (∥ ξ ∥)$ , for a finite non-increasing function $ˆ k$ satisfying $ˆ k (∥ ξ ∥) = Θ (∥ ξ ∥^{- 2 ν - d})$ as $ξ \to \infty$ . Hence

C (1 + ∥ ξ ∥^{2})^{- (ν + d / 2)} \leq {ˆ K}_{θ} (ξ) \leq C^{'} (1 + ∥ ξ ∥^{2})^{- (ν + d / 2)},

for constants $C, C^{'} > 0$ , and we obtain that $H_{θ} (R^{d})$ is equivalent to the Sobolev space $H^{ν + d / 2} (R^{d})$ .

From 2, $H_{θ} (D)$ is given by the restriction of functions in $H_{θ} (R^{d})$ ; as $D$ is Lipschitz, the same is true of $H^{ν + d / 2}$ . $H_{θ} (D)$ is thus equivalent to $H^{ν + d / 2} (D)$ . Finally, functions in $H_{θ} (¯ D)$ are continuous, so uniquely identified by their restriction to $D$ , and

H_{θ} (¯ D) ≃ H_{θ} (D) ≃ H^{ν + d / 2} (D) .

If $ν = \infty$ , by a similar argument $H_{θ} (¯ D)$ is continuously embedded in all $H^{s} (D)$ . ∎

From 1, we can derive results on the behaviour of $∥ f ∥_{H_{θ} (S)}$ as $θ$ varies. For small $θ$ , we obtain the following result.

Lemma 4.

If $f \in H_{θ} (S)$ , then $f \in H_{θ^{'}} (S)$ for all $0 < θ^{'} \leq θ$ , and

∥ f ∥_{H_{θ^{'}} (S)}^{2} \leq (d \prod i = 1 θ_{i} / θ_{i}^{'}) ∥ f ∥_{H_{θ} (S)}^{2} .

Proof.

Let $C = \prod_{i = 1}^{d} (θ_{i}^{'} / θ_{i})$ . As $ˆ K$ is isotropic and radially non-increasing,

{ˆ K}_{θ^{'}} (ξ) = C {ˆ K}_{θ} ((θ_{1}^{'} / θ_{1}) ξ_{1}, \dots, (θ_{d}^{'} / θ_{d}) ξ_{d}) \geq C {ˆ K}_{θ} (ξ) .

Given $f \in H_{θ} (S)$ , let $g \in H_{θ} (R^{d})$ be its minimum norm extension, as in 2. By 1,

∥f∥2Hθ′(S)≤∥g∥2Hθ′(Rd)=∫|ˆg|2ˆKθ′≤∫|ˆg|2CˆKθ=C−1∥f∥2Hθ(S).\qed

Likewise, for large $θ$ , we obtain the following.

Lemma 5.

If $ν < \infty$ , $f \in H_{θ} (S)$ , then $f \in H_{t θ} (S)$ for $t \geq 1$ , and

∥ f ∥_{H_{t θ} (S)}^{2} \leq C^{''} t^{2 ν} ∥ f ∥_{H_{θ} (S)}^{2},

for a $C^{''} > 0$ depending only on $K$ and $θ$ .

Proof.

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

C (1 + ∥ ξ ∥^{2})^{- (ν + d / 2)} \leq {ˆ K}_{θ} (ξ) \leq C^{'} (1 + ∥ ξ ∥^{2})^{- (ν + d / 2)} .

Thus for $t \geq 1$ ,

	${ˆ K}_{t θ} (ξ) = t^{d} {ˆ K}_{θ} (t ξ)$	$\geq C t^{d} (1 + t^{2} ∥ ξ ∥^{2})^{- (ν + d / 2)}$
		$\geq C t^{- 2 ν} (1 + ∥ ξ ∥^{2})^{- (ν + d / 2)}$
		$\geq C C^{' - 1} t^{- 2 ν} {ˆ K}_{θ} (ξ),$

and we may argue as in the previous lemma. ∎

We can also describe the posterior distribution of $f$ in terms of $H_{θ} (S)$ ; as a consequence, we may deduce 1.

Lemma 6.

Suppose $f (x) = μ + g (x)$ , $g \in H_{θ} (S)$ .

${^f}_{n} (x; θ$

Convergence rates of efficient global optimization algorithms

Authors

Convergence of Gaussian process regression: Optimality, robustness, and relationship with kernel ridge regression

A hierarchical expected improvement method for Bayesian optimization

Analysis of Thompson Sampling for Gaussian Process Optimization in the Bandit Setting

Regularized Contextual Bandits

Finding Global Minima via Kernel Approximations

Histogram Transform Ensembles for Large-scale Regression

A simple parameter-free and adaptive approach to optimization under a minimal local smoothness assumption

This week in AI

1 Introduction

2 Expected Improvement

2.1 Bayesian Optimization

2.2 Gaussian Process Models

Assumption 1.

Assumption 2.

Assumption 3.

Assumption 4.

2.3 Expected Improvement Strategies

Definition 1.

Definition 2.

3 Convergence Rates

3.1 Reproducing-Kernel Hilbert Spaces

Lemma 1.

Lemma 2 (Aronszajn, 1950, §1.5).

Lemma 3.

3.2 Fixed Parameters

Theorem 1.

Theorem 2.

3.3 Estimated Parameters

Theorem 3.

Corollary 1.

Definition 3.

Theorem 4.

3.4 Near-Optimal Rates

Definition 4.

Theorem 5.

4 Conclusions

Acknowledgements

Appendix A Proofs

a.1 Reproducing-Kernel Hilbert Spaces

Proof of 1.

Proof of 3.

Lemma 4.

Proof.

Lemma 5.

Proof.

Lemma 6.