Convergence Rates of Epsilon-Greedy Global Optimization Under Radial Basis Function Interpolation

1. Introduction

Consider the optimization problem ${\min }_{\mathbf{x}\in \mathcal{X}}f(x)$, where no closed-form expression for f is available. We can observe the function values $f({\mathbf{x}}_n)$ at individual design points ${\{{\mathbf{x}}_n\}}_{n\ge 1}$ of our choosing (we assume that the observations are noiseless), but we have no information about the derivative of f at these points. This problem class is also known as “global optimization” and “derivative-free optimization” (Conn et al. 2009b) and is often applied to tune parameters in engineering simulators (Giuliani and Camponogara 2015) or machine learning models (Eitrich and Lang 2006).

When the domain $\mathcal{X}\subseteq {\mathbb{R}}^d$ is compact and connected, we construct a function ${\widehat{f}}_N$ (often called a “metamodel” or “surrogate model”) that interpolates the observed function values $f({\mathbf{x}}_1),\dots ,f({\mathbf{x}}_N)$ in some way. This allows us to predict values at points we have not yet observed and to approximate the optimal solution ${\mathbf{x}}^{*}=\arg {\min }_{\mathbf{x}\in \mathcal{X}}f(x)$ by calculating ${\widehat{\mathbf{x}}}_N^{*}=\arg {\min }_{\mathbf{x}\in \mathcal{X}}{\widehat{f}}_N(\mathbf{x})$. We can also use the interpolation to guide the selection of new design points; that is, x_n can be allowed to depend on ${\widehat{f}}_{n-1}$. Thus, the quality of our estimate of the optimal solution (i.e., the difference $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$) is determined by two factors: (1) the particular interpolation method used to construct ${\widehat{f}}_N$ and 2) the policy used to determine $\{{\mathbf{x}}_n\}$ based on previous observations.

In this paper, we derive new results on the convergence rate of $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$ under specific choices for the policy and interpolation. We assume that the metamodel ${\widehat{f}}_N$ is constructed using the method of radial basis functions (RBFs), which is widely used in global optimization, and is closely related to Gaussian process regression (itself a very popular technique). As for the choice of design points, we focus on a relatively simple sequential policy known as ε-greedy (Sutton and Barto 2018): at each time stage n, we either sample uniformly from a small neighborhood of ${\widehat{\mathbf{x}}}_{n-1}^{*}$ (with probability ε), or we sample uniformly from the entire domain $\mathcal{X}$ (with probability $1-\epsilon$); note that this policy is randomized. We will explain the reasons for this choice of policy further down, but first we will state the two main results of the paper: the pathwise convergence rate

where b is determined from ε and k is a parameter of the interpolation model, and the concentration inequality where $c,c^{\prime }$ are problem-specific constants. The second result complements the first by bounding the probability that ${\widehat{\mathbf{x}}}_N^{*}$ has not yet entered the asymptotic regime at a fixed time N.

Although the global optimization literature has a long history and many sophisticated sampling procedures have been developed, these types of results remain fairly rare: most studies focus on empirical performance and/or on weaker theoretical guarantees such as convergence to a first-order critical point. Among papers that do study convergence rates, many require additional structure on f, such as convexity (Bauschke et al. 2015, Duchi et al. 2015) or strong convexity (Berahas et al. 2019). Among the very few papers that do not require such assumptions, we highlight several that study convergence rates in Gaussian process regression, namely Srinivas et al. (2010), Bull (2011), and Vakili et al. (2020). There are some differences between the metrics and sampling procedures whose convergence rates are studied in these papers, but roughly speaking, our results are stronger than those of Bull (2011) because of the way we shrink the local search region (to be discussed further down); without this improvement, we would have obtained similar rates. Regarding the other papers, the results in Vakili et al. (2020) are comparable to those of Bull (2011), and stronger than those of Srinivas et al. (2010). The setting of our paper is more general than any of the three, as RBF is a more general class of interpolation models than Gaussian process regression, allowing fewer assumptions on the kernel function used to build the interpolant (in fact, Vakili et al. 2020 focuses on a single class of kernels, namely the Matérn family).

Several other results are worth mentioning. Recent work by Calvin et al. (2018) obtains a strong rate of $O(e^{-c\sqrt{N}})$ on the optimality gap ${\min }_{1\le n\le N}f({\mathbf{x}}_n)-f({\mathbf{x}}^{*})$ but requires a computationally expensive multilinear interpolation model and the numerical evaluation of complicated integrals; the constant c also vanishes very quickly in the dimension d. The uncertainty quantification literature has studied the convergence rate of the estimation error across the entire domain without an explicit optimization problem; an example is Teckentrup (2020), which studies Gaussian process regression with Matérn kernels. Last, Tikhomirov (2006) derived a bound on the time required to reach a certain accuracy using randomized direct search (without any interpolation model).

In light of this discussion, there is value in studying the ε-greedy policy, which has had a long history in reinforcement learning and is still actively used in applications such as recommender systems (Kamishima and Akaho 2011) and crowdsourcing (Raykar and Agrawal 2014). This simple policy captures the key tradeoff between local and global search, governed by the parameter ε. It enables a tractable analysis of convergence rates under RBF interpolation, and produces bounds that are competitive against the theoretical state of the art. We also present numerical evidence that ε-greedy can also be competitive in some practical situations where the function f is high-dimensional, multimodal, and badly scaled, creating challenges for more sophisticated methods. Furthermore, the rates that we derive also hold for generalizations of ε-greedy where global search is conducted by sampling from an arbitrary density on $\mathcal{X}$ (this changes the multiplicative constant, but not the order of the rate), so in that sense our choice of policy is not restrictive. However, our analysis cannot improve the rates by using nonuniform sampling, because our proof technique relies on a connection between the estimation error $\Vert {\widehat{f}}_N-f\Vert$ under RBFs and the so-called “mesh norm,” which measures how evenly the design points are spread out over $\mathcal{X}$.

We do, however, obtain an insight into the optimal size of the local search region. Namely, we find that the size of the neighborhood around ${\widehat{x}}_N^{*}$ should shrink over time, at a rate proportional to ${\left(\frac{\log (bN)}{bN}\right)}^{\frac{k}{2d}}$. Essentially, if the local search region shrinks too slowly, we will be wasting design points that should have been used to explore the domain; however, if the local search region shrinks too quickly, there is a risk that it will no longer cover $x^{*}$ even when N is very large. Shrinking the local search region at the rate indicated above improves the convergence rate by a factor of ${\left(\frac{\log (bN)}{bN}\right)}^{\frac{k^2}{4d}}$ that otherwise would not be there.

2. Literature Review

There is a large class of global optimization methods that either do not require a metamodel at all or can be applied very generally (with virtually any metamodel). These include heuristics such as evolutionary algorithms (Back 1996), simulated annealing (Corana et al. 1987), and particle swarm (Hu et al. 2004) algorithms. Such approaches have shown promise in global optimization (Schutte and Groenwold 2005, Yang 2010). To give some examples of the available theory, Van den Bergh and Engelbrecht (2006) proved convergence of particle swarm to stationary points, while Vaz and Vicente (2007) proved the existence of a subsequence of design points that converges to a first-order critical point. Orosz and Jacobson (2002) studied the expected number of samples required by simulated annealing to identify a suboptimal solution within some fixed tolerance level; the resulting bounds, however, are difficult to compute and have to be evaluated numerically.

Direct search methods (Torczon 1997) also do not require an interpolation model but rather move toward $x^{*}$ by a sequence of local directional searches. This methodology can handle extensions such as constrained problems (Lewis and Torczon 1999, 2000); see Kolda et al. (2003) for a review of various extensions and improvements. A major advance in this literature was the development of mesh adaptive direct search (Audet and Dennis 2006), which allows substantially more flexibility in the choice of direction. Again, many extensions are possible, for example, to nonsmooth optimization (Audet et al. 2008) or multiobjective optimization (Audet et al. 2010). The theory generally focuses on convergence to first-order critical points, with Abramson and Audet (2006) proving convergence to second-order stationary points.

Trust-region methods conduct local search on a suitably defined region using a metamodel, for example, linear (Powell 1994, Conn et al. 1997), quadratic (Powell 2002), or polynomial (Shashaani et al. 2018) interpolation. The practical potential of RBF interpolation within the trust-region framework was investigated by Wild et al. (2008). With regard to theory, fast convergence rates can be derived when the derivative of f is observable (Shi and Guo 2008), but in the derivative-free setting, the main focus has been on global convergence to first-order (Wild and Shoemaker 2011) and sometimes second-order (Conn et al. 2009a) critical points.

By contrast, methods based on Gaussian process regression typically do not explicitly distinguish between local and global search, but rather accomplish this tradeoff through a stochastic metamodel with built-in uncertainty quantification. The most popular algorithmic approach in this stream is expected improvement (Jones et al. 1998) and its many variants (Sasena et al. 2002, Huang et al. 2006). The theory primarily focused on the pointwise consistency of the metamodel (Vazquez and Bect 2010) until the convergence rate analysis of Bull (2011) and Vakili et al. (2020), which were discussed earlier. Closely related is the probability of improvement criterion (Zhigljavsky and Zilinskas 2008), which motivated the rate analysis of Calvin et al. (2018), also discussed previously.

Last, RBF interpolation has had a long history in numerical analysis (Buhmann 2003) outside the setting of global optimization. Our analysis draws on this literature, specifically theory by Wu and Schaback (1993) characterizing the convergence rate of the estimation error (under RBF interpolation) given an arbitrary collection of design points. The first RBF-based global optimization procedure was proposed by Gutmann (2001), with later improvements by Regis and Shoemaker (2007) and Holmström (2008). Sampling in these papers is based on a measure of the smoothness of the interpolation, with additional logic for balancing global and local search. Other sampling criteria have also been considered: for example, Regis and Shoemaker (2005) aims to spread out the design points to avoid excessive clustering. Extensions include parallelized methods (Regis and Shoemaker 2009) and hybrid methods combining RBF with coordinate search (Regis and Shoemaker 2013). Much of this work is computationally oriented and focuses on complex engineering applications, with the theory mostly limited to global convergence. In the computer science community, Srinivas et al. (2010) derived rate results that apply to RBFs, but the setting there is online learning, where one optimizes cumulative error over time, rather than the offline setting more typical of global optimization (global search plays a much greater role in offline algorithms).

Overall, the approach and results presented here are intended, not to supplant the existing work on global optimization with RBFs, but to complement it from a theoretical viewpoint. Our paper adds to a very small number of prior studies of convergence rates for derivative-free optimization.

3. Problem Statement and RBF Interpolation

Let f be a function defined on a compact and connected domain $\mathcal{X}\subseteq {\mathbb{R}}^d$, and suppose that ${\mathbf{x}}^{*}=\arg {\min }_{\mathbf{x}\in \mathcal{X}}f(\mathbf{x})$ is the unique global minimizer of f. Let ${\{{\mathbf{x}}_n\}}_{n=1}^N$ be a finite sequence of design points in ${\mathbb{R}}^d$, where for each x_n we observe $f({\mathbf{x}}_n)$ without noise. The design points can be predetermined by the decision maker or chosen adaptively; for the moment, however, suppose that they are simply given and that all the observations have been made. Using these observations, we construct a RBF interpolation ${\widehat{f}}_N$ of f and use ${\widehat{\mathbf{x}}}_N^{*}=\arg {\min }_{\mathbf{x}\in \mathcal{X}}{\widehat{f}}_N(\mathbf{x})$ as our estimate of ${\mathbf{x}}^{*}$.

The RBF interpolation follows the method in Wu and Schaback (1993). Let ${\mathcal{P}}_q$ be a space of polynomial functions on ${\mathbb{R}}^d$ with total order not exceeding q; when q > 0, suppose that, if for some $p\in {\mathcal{P}}_q$ we have $p({\mathbf{x}}_n)=0$ for all $n=1,\dots ,N$, then $p\equiv 0$. Now let $\varphi :{\mathbb{R}}_{+}\to \mathbb{R}$ be a function chosen to make the mapping $r\mapsto \varphi (\sqrt{r})$ conditionally positive definite of order q. This means that the N × N kernel matrix, whose (i, j)th entry is $\varphi (\Vert x_i-{\mathbf{x}}_j\Vert )$, is positive definite on the set of $u\in {\mathbb{R}}^d$ satisfying

where $Q=\left(\begin{array}{c}q+d-1\\ d\end{array}\right)$ is the dimension of ${\mathcal{P}}_q$, and $(p_1,\dots ,p_Q)$ is any basis of ${\mathcal{P}}_q$. The interpolation ${\widehat{f}}_N$ then has the form where the coefficients $(a_1,\dots ,a_N)$ and $(a_1^{\prime },\dots ,a_Q^{\prime })$ constitute the solution to the linear system

The conditional positive definiteness assumptions on $\varphi$ guarantee that this system is nonsingular.

Given the RBF $\varphi$, let $\psi (\mathbf{x})=\varphi (\Vert \mathbf{x}\Vert )$ and take $\tilde{\psi }$ to be the Fourier transform of ψ (analogously, $\tilde{f}$ is the Fourier transform of f ). We define

and require f to satisfy $c_{f,\varphi }^2<\infty$, a condition that is also imposed in Wu and Schaback (1993). Moreover, although f itself is not required to be convex, we assume local strong convexity of f around its global minimizer: that is, there exists a closed ball centered at ${\mathbf{x}}^{*}$ with strictly positive radius where f is strongly convex.

The condition $c_{f,\varphi }^2<\infty$ is the only assumption in our paper that relates f to the RBF $\varphi$. There may be many choices of $\varphi$ that satisfy this condition for a given f. Similarly, for a given $\varphi$, there may be many functions f for which the condition holds. This is an important difference between our model and Gaussian process regression, as studied by Vakili et al. (2020), Teckentrup (2020), and others. These papers model f as a realization from a distribution that depends on some RBF, which gives rise to a misspecification problem when this true RBF is unknown. In our setting, there is no true RBF, and our theory will hold for any $\varphi$ satisfying $c_{f,\varphi }^2<\infty$.

For our analysis, we will also require the constant $k=\frac{1}{2}{s}_{\infty }$, where $s_{\infty }$ satisfies

for $\Vert t\Vert \to \infty$. This condition is satisfied by many commonly used kernels. For example, if $\varphi (r)=e^{-\alpha r^2}$ is the Gaussian kernel (for some $\alpha >0$), we may have arbitrarily large values of $s_{\infty }$, leading to arbitrarily large values of k, which in effect causes ${\widehat{\mathbf{x}}}_n^{*}$ to converge to ${\mathbf{x}}^{*}$ even more quickly than the rate we derive in Section 5. For other types of kernels, k can often be computed explicitly. For example, Teckentrup (2020) can be used to compute it for the Matérn kernel; several other examples are given in Wu and Schaback (1993). In general, smaller values of $s_{\infty }$ allow a wider range of functions f that satisfy $c_{f,\varphi }^2<\infty$, but slow down the convergence.

To measure the local density of the design points, denote by

the mesh norm of the compact subset $\mathcal{D}$. Letting $\mathbf{\mu }\in {\mathcal{N}}^d$ with $|\mathbf{\mu }|={\sum }_j{\mathbf{\mu }}_j$, we use the standard multi-index notation $f^{(\mathbf{\mu })}$ for the function obtained after sequentially applying to f the μ_jth-order partial derivative with respect to x_j. We also use the notation $B(\mathbf{x},r)=\{\mathbf{y}:\Vert \mathbf{x}-\mathbf{y}\Vert \le r\}$ to represent the ball of radius r centered at $\mathbf{x}\in {\mathbb{R}}^d$. With these preliminaries, we can now state a result from Wu and Schaback (1993) that will be referenced and applied throughout this paper.

Lemma 1.

(Wu and Schaback 1993). Suppose that f is C² on $\mathcal{X}$, and the RBF $\varphi$ is C² on $(0,\infty )$ and C⁴ in a neighborhood of zero, with $k\ge 2$. Then, for $\rho >0$, there exists $C\in {\mathbb{R}}_{+}$ such that, for any set of design points and any $\mathbf{x}\in \mathcal{X}$ satisfying $h_{B(\mathbf{x},\rho )\cap \mathcal{X}}<h_0$, with h₀ being a constant whose value depends on k, we have the inequality

$$\left|{\widehat{f}}_N^{(\mathbf{\mu })}(\mathbf{x})-f^{(\mathbf{\mu })}(\mathbf{x})\right|\le c_{f,\varphi }C{h}_{B(\mathbf{x},\rho )\cap \mathcal{X}}^{k-|\mathbf{\mu }|}$$

for all μ satisfying $|\mathbf{\mu }|\le 2$. Above, $h_{B(\mathbf{x},\rho )\cap \mathcal{X}}^{k-|\mathbf{\mu }|}={(h_{B(\mathbf{x},\rho )\cap \mathcal{X}})}^{k-|\mathbf{\mu }|}$, that is, the mesh norm is raised to the power $k-|\mathbf{\mu }|$.

The assumptions of Lemma 1 on $\varphi$ will also be carried over to our analysis. All the existing work on convergence rates of RBF models makes some assumptions on the smoothness of $\varphi$. Many of these are arguably stronger than ours: for example, Bull (2011) requires the mapping $\mathbf{x}\mapsto \varphi (\Vert \mathbf{x}\Vert )$ to be $C^{2k-1}$, so the needed order of smoothness grows with k.

4. Properties of Local Mesh Norms

In the following, we prove several useful properties of the local mesh norm defined in (1) on a particular class of domains. These results will motivate certain algorithmic design choices in Section 5 and will also be used in the subsequent analysis. To prove the results, we build on the properties of several other measures of local data density, which we show to be equivalent to (1).

Definition 1.

A compact set $\mathcal{D}\subseteq {\mathbb{R}}^d$ is shape-regular if there exists a continuously differentiable bijection $L_{\mathcal{D}}$, mapping either ${[0,1]}^d$ or $B(0,1)$ (either can be chosen as the domain) onto $\mathcal{D}$, whose Jacobian has nonzero determinant everywhere on the domain.

Definition 2.

Let ${\mathcal{X}}_N={\{{\mathbf{x}}_n\}}_{n=1}^N$ and define

$$\begin{array}{l}{\breve{h}}_{\mathcal{D}}({\mathcal{X}}_N)=\underset{\mathbf{y}\in \mathcal{D}}{\sup }\underset{\mathbf{x}\in {\mathcal{X}}_N\cap \mathcal{D}}{\inf }\Vert \mathbf{x}-\mathbf{y}{\Vert }_2,\\ {\overline{h}}_{\mathcal{D}}({\mathcal{X}}_N)=\underset{\mathbf{y}\in \mathcal{D}}{\sup }\underset{\mathbf{x}\in {\mathcal{X}}_N\cup \partial \mathcal{D}}{\inf }\Vert \mathbf{x}-\mathbf{y}{\Vert }_2,\\ {\tilde{h}}_{\mathcal{D}}({\mathcal{X}}_N)=\underset{\mathbf{y}\in \mathcal{D}}{\sup }\underset{\mathbf{x}\in {\mathcal{X}}_N\cup \partial \mathcal{D}}{\inf }\Vert \mathbf{x}-\mathbf{y}{\Vert }_{\infty }.\end{array}$$

For simplicity, we may omit the explicit dependence of these quantities on ${\mathcal{X}}_N$ from the notation when there is no ambiguity.

To compare these and other quantities, we introduce the following notation. For two positive sequences $\{F_n^1\}$ and $\{F_n^2\}$, we write $F_n^1\mathrm{\lesssim }{F}_n^2$ if there exists a constant c, independent of n, such that $F_n^1\le c{F}_n^2$ for all n. We write $F_n^1\cong F_n^2$ when ${\lim }_{n\to \infty }\frac{F_n^1}{F_n^2}=1$ (this is stronger than having both $F_n^1\mathrm{\lesssim }{F}_n^2$ and $F_n^2\mathrm{\lesssim }{F}_n^1$).

Lemma 2.

Let $\mathcal{D}\subseteq \mathcal{X}$. For any N and any ${\mathcal{X}}_N$, we have ${\overline{h}}_{\mathcal{D}}\mathrm{\lesssim }{\tilde{h}}_{\mathcal{D}}$ and ${\tilde{h}}_{\mathcal{D}}\mathrm{\lesssim }{\overline{h}}_{\mathcal{D}}$. We also have ${\overline{h}}_{\mathcal{D}}\le h_{\mathcal{D}}\le {\breve{h}}_{\mathcal{D}}$.

Proof.

The first statement follows from the equivalence of the l²-norm and $l^{\infty }$-norm. The second statement follows from the relationship ${\mathcal{X}}_N\cap \mathcal{D}\subseteq {\mathcal{X}}_N\subseteq {\mathcal{X}}_N\cup \partial \mathcal{D}$. □

From Lemma 1, we know that the interpolation error is bounded by some power of the local mesh norm. Therefore, to study the convergence rate of the interpolation error under any sampling scheme, we essentially require an appropriate decreasing rate for the local mesh norm. Janson (1987) derived such a rate for the specific case where the design points are sampled from a uniform distribution on $\mathcal{X}$. For convenience, we state this result here.

Lemma 3

(Janson 1987). Suppose that $\mathcal{X}={[0,1]}^d$ or $\mathcal{X}=B(0,1)$. Suppose also that the design points ${\{{\mathbf{x}}_n\}}_{n=1}^N$ are sampled independently from a uniform distribution on $\mathcal{X}$. Then,

$${\tilde{h}}_{\mathcal{X}}=O\left({\left(\frac{\log N}{N}\right)}^{\frac{1}{d}}\right)$$(2)

almost surely. Furthermore, the multiplicative constant in (2) is nonrandom, that is, the limit superior ${\lim \sup }_{N\to \infty }{\tilde{h}}_{\mathcal{X}}{\left(\frac{\log N}{N}\right)}^{-\frac{1}{d}}$ is almost surely (a.s.) equal to a deterministic quantity.

Building on this result, we prove an analogous rate for a more general case where the design points are sampled from an arbitrary probability distribution. This result helps motivate our subsequent focus on the epsilon-greedy policy, which uses uniform sampling for global search, because the convergence rate of the mesh norm is faster when the essential infimum of the sampling density is higher.

Lemma 4.

Suppose that $\mathcal{X}={[0,1]}^d$ or $\mathcal{X}=B(\mathbf{0},1)$. Suppose also that the design points ${\{{\mathbf{x}}_n\}}_{n=1}^N$ are sampled independently from a probability distribution with density $g:\mathcal{X}\to \mathbb{R}$. Let $g_{\min }= \text{ess}\hspace{0.17em}{\text{inf}}_{\mathit{\text{x}}\in \mathcal{X}}\mathit{\text{g}}(\mathit{\text{x}})$ and suppose that $g_{\min }>0$. Then,

$$h_{\mathcal{X}}=O\left({\left(\frac{\log \left(g_{\min }N\right)}{g_{\min }N}\right)}^{\frac{1}{d}}\right)$$

almost surely.

Proof.

We prove this lemma for $\mathcal{X}={[0,1]}^d$, as the proof for the closed ball is similar. Define random variables $z_n\sim \mathit{Bernoulli}(g_{\min })$. Then, the distribution of the design points can be rewritten as follows: if z_n = 1, draw ${\mathbf{x}}_n\sim U({[0,1]}^d)$, and if z_n = 0, draw x_n from a distribution with density $\mathbf{x}\mapsto \frac{g(\mathbf{x})-g_{\min }}{1-g_{\min }}$. We can also denote by ${\mathcal{X}}_N^U=\{{\mathbf{x}}_n:z_n=1\}$ the subset of the design points coming from the uniform density, with $N^U={\sum }_{n=1}^N{z}_n$ being the size of this subset. By the strong law of large numbers, $N^U\cong g_{\min }N$ almost surely.

Now, observe that for general $\mathcal{D}\subseteq \mathcal{X}$ and ${\mathcal{Y}}_1\subseteq {\mathcal{Y}}_2\subseteq \mathcal{X}$, we have ${\tilde{h}}_{\mathcal{D}}({\mathcal{Y}}_1)\ge {\tilde{h}}_{\mathcal{D}}({\mathcal{Y}}_2)$. Therefore, on $\mathcal{D}={[0,1]}^d$, we have

$$\begin{array}{l}{\tilde{h}}_{{[0,1]}^d}\left({\mathcal{X}}_N\right)\le {\tilde{h}}_{{[0,1]}^d}\left({\mathcal{X}}_N^U\right)\\ =O\left({\left(\frac{\log N^U}{N^U}\right)}^{\frac{1}{d}}\right)\hspace{1em}\hspace{1em} \text{a}.\text{s}.\\ =O\left({\left(\frac{\log \left(g_{\min }N\right)}{g_{\min }N}\right)}^{\frac{1}{d}}\right)\hspace{1em}\hspace{1em} \text{a}.\text{s}.\end{array}$$

It can be easily shown that

$${\breve{h}}_{{[0,1]}^d}\left({\mathcal{X}}_N\right)\le \left(1+\sqrt{d}\right){\overline{h}}_{{[0,1]}^d}\left({\mathcal{X}}_N\right)$$

because of the geometry of the unit cube and the fact that ${\overline{h}}_{{[0,1]}^d}$ is the radius of the largest ball inside ${[0,1]}^d$ with no design points in its interior. Combining this with Lemma 2, we obtain $h_{{[0,1]}^d}\lesssim {\overline{h}}_{{[0,1]}^d}\mathrm{\lesssim }{\tilde{h}}_{{[0,1]}^d}$, which completes the proof. □

Using Lemma 4, we can extend the result of Lemma 3 to a more general class of domains, while keeping the assumption of uniformly distributed design points.

Lemma 5.

Let $\mathcal{X}$ be compact and shape-regular, and suppose that the design points ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_N$ are sampled independently from a uniform distribution on $\mathcal{X}$. Then,

$$h_{\mathcal{X}}\left({\mathcal{X}}_N\right)=O\left({\left(\frac{\log N}{N}\right)}^{\frac{1}{d}}\right)\hspace{1em}\hspace{1em} \mathit{\text{a}}.\mathit{\text{s}}.$$

Proof.

By Definition 1, there exists a continuously differentiable bijection $L_{\mathcal{X}}:{[0,1]}^d\to \mathcal{X}$ whose Jacobian has nonzero determinant everywhere on the domain. Here we use ${[0,1]}^d$ as the domain of $L_{\mathcal{X}}$, but the proof is similar if $B(0,1)$ is used instead.

Because $L_{\mathcal{X}}$ is continuously differentiable on a compact set, it is Lipschitz. Then,

$$\begin{array}{l}{h}_{\mathcal{X}}\le {\breve{h}}_{\mathcal{X}}\\ =\underset{\mathbf{y}\in \mathcal{X}}{\sup }\underset{\mathbf{x}\in {\mathcal{X}}_N}{\inf }\Vert \mathbf{x}-\mathbf{y}{\Vert }_2\\ =\underset{\mathbf{y}\in {[0,1]}^d}{\sup }\underset{\mathbf{x}\in L_{\mathcal{X}}^{-1}\left({\mathcal{X}}_N\right)}{\inf }\Vert L_{\mathcal{X}}\left(\mathbf{x}\right)-L_{\mathcal{X}}\left(\mathbf{y}\right){\Vert }_2\\ \le \underset{\mathbf{y}\in {[0,1]}^d}{\sup }\underset{\mathbf{x}\in L_{\mathcal{X}}^{-1}\left({\mathcal{X}}_N\right)}{\inf }{c}_{L_{\mathcal{X}}}\Vert \mathbf{x}-\mathbf{y}{\Vert }_2,\end{array}$$

where $c_{L_{\mathcal{X}}}$ is the Lipschitz constant of $L_{\mathcal{X}}$. Using similar arguments as in the proof of Lemma 4, we obtain

$$h_{\mathcal{X}}\mathrm{\lesssim }{h}_{{[0,1]}^d}\left(L_{\mathcal{X}}^{-1}\left({\mathcal{X}}_N\right)\right).$$

Let $|\mathcal{X}|$ be the volume of the domain $\mathcal{X}$ under Lebesgue measure. The design points $L_{\mathcal{X}}^{-1}({\mathcal{X}}_N)$ are drawn from a distribution with density

$$g_0\left(\mathbf{y}\right)=\frac{1}{|\mathcal{X}|}\left|\det \left(\frac{d{L}_{\mathcal{X}}\left(\mathbf{y}\right)}{d\mathbf{y}}\right)\right|$$

whose essential infimum satisfies

$$g_{0,\min }= \underset{\text{y}\in {[0,1]}^{\text{d}}}{\text{ess}\hspace{0.17em}\text{inf}}{\mathit{\text{g}}}_0\left(\mathit{\text{y}}\right)>0$$

because of the assumptions on $L_{\mathcal{X}}$. Then, by Lemma 4, for the mesh norm on ${[0,1]}^d$ with design points $L_{\mathcal{X}}^{-1}({\mathcal{X}}_N)$, we have

$$h_{{[0,1]}^d}\left(L_{\mathcal{X}}^{-1}\left({\mathcal{X}}_N\right)\right)=O\left({\left(\frac{\log \left(g_{0,\min }N\right)}{g_{0,\min }N}\right)}^{\frac{1}{d}}\right)=O\left({\left(\frac{\log N}{N}\right)}^{\frac{1}{d}}\right)$$

almost surely. This completes the proof. □

5. The Epsilon-Greedy Algorithm

Section 5.1 gives a formal statement of the sequential algorithm used to select design points and states our main theoretical results on its convergence rate. The proofs of these results are given later in Sections 6 and 7. Before proceeding to these proofs, we present a numerical example in Section 5.2 of a difficult problem where ε-greedy is demonstrably useful.

5.1. Algorithm and Main Results

The theoretical foundation for our analysis comes from interpolation theory: Lemma 1 connects the error $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$ to the mesh norm, which can be made tractable under uniform sampling as seen in Lemma 5. This motivates the structure of our algorithm. Epsilon-greedy can be seen as a mixture of two uniform distributions, one on the domain $\mathcal{X}$ (global search) and one on a neighborhood of the estimated optimum (local search). The specific rate of the mesh norm obtained in Lemma 5 further motivates how we set the size of the local search region, which is crucial for improving the performance of the algorithm.

Let $S(\mathbf{x},r)=\{\mathbf{y}\hspace{0.17em}:\hspace{0.17em}\Vert \mathbf{x}-\mathbf{y}{\Vert }_{\infty }\le \frac{r}{2}\}$ be the hypercube centered at x with sides of length r parallel to the coordinate axes. The ε-greedy algorithm, as defined in this paper, will randomly choose between uniform sampling inside this hypercube centered at the current estimated minimizer ${\widehat{\mathbf{x}}}_n^{*}$ (local search) and uniform sampling on $\mathcal{X}$ (global search). Formally, let z_n be a Bernoulli random variable with success probability $b·|\mathcal{X}|$, where $b\in \left(0,\frac{1}{|\mathcal{X}|}\right)$ is a constant and $|\mathcal{X}|$ is, again, the volume of $\mathcal{X}$ under Lebesgue measure. The success probability $b·|\mathcal{X}|$ corresponds to $1-\epsilon$ in the ε-greedy policy and governs the frequency of global search.

The distribution of the nth design point x_n is determined adaptively, after ${\mathbf{x}}_0,\dots ,{\mathbf{x}}_{n-1}$ and $f({\mathbf{x}}_0),\dots ,f({\mathbf{x}}_{n-1})$ have been observed, in the following way: when z_n = 1, x_n is sampled from a uniform distribution on $\mathcal{X}$. When z_n = 0, x_n is sampled from a uniform distribution on $S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)$, where the side length r_n will be discussed later, but must be small enough to make $S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)\subseteq \mathcal{X}$. In the situation where ${\widehat{\mathbf{x}}}_{n-1}^{*}$ is exactly on the boundary of $\mathcal{X}$, local exploration will sample uniformly on $S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)\cap \mathcal{X}$ instead. Treating z_n as a latent variable, x_n follows a distribution whose density g_n is a weighted average of two uniform distributions defined on distinct regions:

where t_n is a constant satisfying $t_n>\frac{1}{|\mathcal{X}|}>b$, whose value can be determined by the normalization condition from the values of b and r_n. Note that $\text{ess}\hspace{0.17em}{\text{inf}}_{\mathit{\text{x}}\in \mathcal{X}}{\mathit{\text{g}}}_{\mathit{\text{n}}}(\mathit{\text{x}})=b$ for all n.

In our analysis, the sequence $\{r_n\}$ of side lengths for the local search region is set according to

where c_r is a large enough constant, and k is the constant in Lemma 1, which is computed from the choice of kernel, as discussed in Section 3. The specific rate chosen in (5) is related to the decreasing rate of the mesh norm from Lemmas 1 and 5 and plays a major role in proving the main results of this paper. The first major result bounds the convergence rate of $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$ on almost every sample path.

Theorem 1.

Assume that $\mathcal{X}$ is shape-regular, f is strongly convex on $B({\mathbf{x}}^{*},{\rho }_0)$ for some constant radius ${\rho }_0>0$, both f and $\varphi$ satisfy the assumptions of Lemma 1, and $\{r_n\}$ is chosen according to (5). Then,

$$\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert =O\left({\left(\frac{\log N}{N}\right)}^{\frac{k}{2d}}{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k^2}{4d}}\right)\hspace{1em}\hspace{1em} \mathit{\text{a}}.\mathit{\text{s}}.$$(6)

In (6), the factor ${\left(\frac{\log N}{N}\right)}^{\frac{k}{2d}}$ is a kind of “baseline” obtained from the theoretical framework of Lemma 1. As will be shown in the proof of Theorem 1, this factor represents the convergence rate achieved when the size of the local search region is constant, i.e., $r_n\equiv r_0$ for some $r_0>0$. The second factor represents the additional improvement achieved by shrinking r_n. The specific rate in (5) leads to the best possible improvement that can be achieved using our techniques.

The second major result is a concentration inequality for the rate in Theorem 1. Essentially, for any given N there may be a nonnegligible set of sample paths on which the rate in (6) does not hold. The probability measure of this set is bounded as follows.

Theorem 2.

Suppose that we are in the situation of Theorem 1. Then, there exists a deterministic $N^{\prime }$, which depends only on b and the domain $\mathcal{X}$, such that for all $N\ge N^{\prime }$, we have

$$P\left(\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert >c{\left(\frac{\log N}{N}\right)}^{\frac{k}{2d}}{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k^2}{4d}}\right)\le \frac{c^{\prime }}{N},$$

where c, $c^{\prime }$ are constants.

We also considered a variant of the algorithm where b was also allowed to vary over n (recall from (3) that b is the value of the sampling density outside the local search region), as well as over $\mathcal{X}$. However, using nonuniform global search does not help the bound on the convergence rate, because our analysis relies on a connection between the estimation error and the mesh norm, and the convergence rate of the mesh norm for an arbitrary density g depends on $\text{ess}\hspace{0.17em}{\text{inf}}_{\mathit{\text{x}}}\mathit{\text{g}}(\mathbf{x})$ as was seen in Lemma 4. If we then use uniform global search, but allow b to vary over time, we find that the bound becomes worse when ${\lim \inf }_n{b}_n=0$ (in fact, if b_n vanishes too quickly, we may not even have ${\widehat{\mathbf{x}}}_N^{*}\to {\mathbf{x}}^{*}$). On the other hand, when b_n varies between constant, strictly positive lower and upper bounds, the order of the rate does not change, and these bounds only contribute to the multiplicative constant. For this reason, we decided not to overcomplicate the presentation with these details and have simply used a constant b in the following.

5.2. Numerical Example

Before we proceed to the proofs of Theorems 1 and 2, we first present a numerical example of a problem in which ε-greedy is an effective solution strategy. Some context is required for this statement, because (as discussed in Section 2) global optimization is a mature field with many different algorithmic approaches. Nonetheless, even from a purely practical standpoint, ε-greedy can be useful if the problem is sufficiently large and difficult. In such a setting, even the most sophisticated methods will inevitably require a certain amount of trial and error, making ε-greedy a reasonable alternative.

We consider one such difficult problem, namely the Rastrigin function

in d = 30 dimensions. The function has been scaled so that its domain is the unit cube centered at the origin; the global minimum is also at the origin, and the optimal value is zero. Figure 1 provides a visualization in d = 2 dimensions, showing that the function is highly multimodal and badly scaled. These issues are exacerbated as d grows.

We compared ε-greedy against the methods of Regis and Shoemaker (2005) and Regis and Shoemaker (2007) on this function. Both methods (referred to, respectively, as RS05 and RS07) are based on RBF interpolation, and thus are natural benchmarks for the framework of this paper. For context, we note that the extensive computational study in Regis and Shoemaker (2007) does not consider any problem with more than 6 dimensions. Thus, the 30-dimensional Rastrigin function is a very large and difficult problem by the standards of this literature.

All three methods were implemented with the Gaussian kernel $\varphi (r)=e^{-0.073r^2}$. The RBF interpolant was constructed as laid out in Section 3, but for simplicity we did not use any polynomial components. All three methods were initialized using a single point uniformly sampled on the domain; the initialization strategies used by RS05 and RS07 as originally proposed are not practical in 30 dimensions because they would require too many points. All three methods were also allowed to use the well-known practical strategy, dating back to Björkman and Holmström (2000), of replacing large function values by the median of all sampled values. Finally, all three methods used the same subroutine, based on trust-region optimization (chapter 14 of Conn et al. 2000), for computing the minimizer ${\widehat{\mathbf{x}}}_n^{*}$ of ${\widehat{f}}_n$.

We made some adjustments to the benchmarks to make their computational cost roughly comparable to that of ε-greedy (to the extent possible). In every iteration, RS07 chooses from a certain number of points generated uniformly at random. This number is set to 1,000d in the method as originally proposed, but that is obviously not practical in 30 dimensions, so we used 2,250 samples per iteration instead. RS05 solves two auxiliary optimization problems on a discretization of the domain; we used 1,000 points in the discretization. The settings of ε-greedy were set to c_r = 50, k = 50, and b = 0.3.

The focus of our theoretical study is the metric $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$. This metric, averaged over 10 runs, is reported in Figure 2(a) with 95% confidence intervals for up to n = 600 function evaluations. We see that, by the end of this time frame, ε-greedy is much more effective at estimating the location of the global minimizer. RS05 is marginally better up to about n = 320 but then slows down. RS07 never really improves on its initial estimate.

To better understand what is happening, let ${\tilde{\mathbf{x}}}_N=\arg {\min }_{n\le N}f({\mathbf{x}}_n)$ be the best among the first N points observed. Figure 2(c) reports $|f({\tilde{\mathbf{x}}}_N)-f({\mathbf{x}}^{*})|$, a measure of the quality of the observed points. Because ε-greedy chooses many points uniformly at random, one would not expect it to do well under this metric. Yet, by n = 600, we find some overlap between the CIs for all three methods, indicating that ε-greedy has largely caught up with the benchmarks. One can also improve the result for ε-greedy by increasing the probability of global search, though this will make the primary metric $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$ worse.

Figure 2(d) reports $\Vert {\tilde{\mathbf{x}}}_N-{\mathbf{x}}^{*}\Vert$. Here we see very clearly that RS07 is actually moving away from ${\mathbf{x}}^{*}$ over time, toward a local optimum that is producing seemingly low values of $f({\tilde{\mathbf{x}}}_N)$. RS05 gets stuck at a local optimum very early on (also seen in Figure 2(c)) but is actually building a more accurate RBF model with the points it samples, explaining the improvement in Figure 2(a). The wide CIs seen in Figure 2(c) for RS05 are due to the instability of the auxiliary problem solved by this procedure; because RS05 is stuck after approximately 50 iterations, the confidence bounds do not shrink.

Figure 2(b) shows that all three methods incur higher computational cost over time, which is to be expected because the complexity of the RBF model increases with the number of design points. There is some fluctuation in the cost (most notably for RS05) because of the numerical instability of the auxiliary optimization problems solved by the methods. Overall, ε-greedy and RS07 have roughly similar cost, whereas RS05 is significantly slower.

In conclusion, ε-greedy shows itself to be more effective in estimating the global minimizer than the two benchmarks. At the same time, it takes longer to find a good solution, as measured by $f({\tilde{\mathbf{x}}}_n)$. In our experience, the tradeoff between these goals can be tuned using the parameter b. Of course, we do not claim that such a simple algorithm as ε-greedy will “always” outperform two of the leading computational benchmarks in the field. Rather, the conclusion that should be drawn is that there are certain settings where even the most sophisticated methods will struggle, and very competitive performance can then be obtained from very simple algorithms.

6. Proof of Theorem 1

We now move to the proofs of the main results of our paper. The pathwise rate in Theorem 1 can be obtained once ${\widehat{\mathbf{x}}}_N^{*}$ is in a neighborhood $B({\mathbf{x}}^{*},{\rho }_0)$ of the global optimum, in which we can make use of the local strong convexity of f. Thus, the proof of the theorem consists of two parts. Section 6.1 proves that ${\widehat{\mathbf{x}}}_N^{*}$ will remain in this neighborhood for all sufficiently large N. Then, in Section 6.2 we can assume that ${\widehat{\mathbf{x}}}_N^{*}\in B({\mathbf{x}}^{*},{\rho }_0)$ and study the convergence rate.

6.1. Convergence of ${\widehat{\mathbf{x}}}_N^{*}$ to a Neighborhood of ${\mathbf{x}}^{*}$

The main goal in this section is to prove

which means that, asymptotically, ${\widehat{\mathbf{x}}}_N^{*}\in B({\mathbf{x}}^{*},{\rho }_0)$ with probability (w.p.) 1. The proof consists of several steps, but is self-contained; that is, in Section 6.2 we will assume that N is large enough so that ${\widehat{\mathbf{x}}}_N^{*}\in B({\mathbf{x}}^{*},{\rho }_0)$, but otherwise we will not refer to any of the intermediate results in this section. The flow of the proof is fairly straightforward, with each new step building on the previous one.

Our analysis proceeds by deriving a lower bound on $P({\widehat{\mathbf{x}}}_N^{*}\in B({\mathbf{x}}^{*},{\rho }_0))$. This bound will eventually be shown to converge to 1. Similarly to the proof of Lemma 4, we rewrite the distribution of x_n as follows. For any n, let $z_n\sim \mathit{Bernoulli}(b·|\mathcal{X}|)$ be a latent variable; then, if z_n = 1, draw ${\mathbf{x}}_n\sim U(\mathcal{X})$, and if z_n = 0, draw x_n from a distribution with density $\mathbf{x}\mapsto \frac{g_n-b}{1-b|\mathcal{X}|}$, where g_n is as in (3). We then denote by ${\mathcal{X}}_N^U=\{{\mathbf{x}}_n:z_n=1\}$ the subset of the design points coming from the uniform density, with $N^U={\sum }_{n=1}^N{z}_n$ being the number of such samples.

The first step in the proof, stated in the following technical lemma, is to bound the probability $P({\widehat{\mathbf{x}}}_N^{*}\in B({\mathbf{x}}^{*},{\rho }_0))$ of interest by the sum of two more tractable probabilities, which will also be bounded in subsequent steps. The proof is moved to the Appendix.

Lemma 6.

Define $\tilde{\mathbf{x}}=\arg {\min }_{\mathbf{x}\in \text{cl}\left(\mathcal{X}\setminus B({\mathbf{x}}^{*},{\rho }_0)\right)}f(\mathbf{x})$ and

$$\mathcal{D}=\left\{\mathbf{x}\in \mathcal{X}\hspace{0.17em}:\hspace{0.17em}f\left(\mathbf{x}\right)<\frac{1}{2}\left(f\left(\tilde{\mathbf{x}}\right)+f\left({\mathbf{x}}^{*}\right)\right)\right\}.$$

There exists a positive constant c_w, independent of the sampling policy and the design points, such that

$$P\left({\widehat{\mathbf{x}}}_N^{*}\in B\left({\mathbf{x}}^{*},{\rho }_0\right)\right)\ge P\left(h_{\mathcal{X}}\left({\mathcal{X}}_N^U\right)<c_w\right)+P\left({\mathcal{X}}_N^U\cap \mathcal{D}\ne \varnothing \right)-1.$$(8)

The next lemma, whose proof is also moved to the Appendix, bounds the first term on the right-hand side of (8). After that, in Lemma 8, we will bound the second term and combine this with the result of Lemma 7 to complete the proof of the theorem.

Lemma 7.

Let c_w be the constant obtained from Lemma 6. There exists another positive constant ${\overline{c}}_w$ such that

$$P\left(h_{\mathcal{X}}\left({\mathcal{X}}_N^U\right)<c_w\right)\ge 1-{\overline{c}}_w{e}^{-\frac{1}{2}{b}^2{|\mathcal{X}|}^{2}N}.$$

Lemma 8.

There exists a positive constant ${\widehat{c}}_w$ such that

$$P\left({\widehat{\mathbf{x}}}_N^{*}\in \mathcal{X}\setminus B\left({\mathbf{x}}^{*},{\rho }_0\right)\right)\mathrm{\lesssim }{e}^{-{\widehat{c}}_{w}N}.$$

Proof.

Let $\tilde{\mathbf{x}}$ and $\mathcal{D}$ be as in the statement of Lemma 6. We calculate

$$\begin{array}{c}P\left({\mathcal{X}}_N^U\cap \mathcal{D}\ne \varnothing \hspace{0.17em}|\hspace{0.17em}{N}^U\right)=1-P\left({\mathcal{X}}_N^U\subseteq \mathcal{X}\setminus \mathcal{D}\hspace{0.17em}|\hspace{0.17em}{N}^U\right)\\ =1-{\left(1-\frac{|\mathcal{D}|}{|\mathcal{X}|}\right)}^{N^U}.\end{array}$$

Taking the expectation over the distribution of N^U, we obtain

$$1-P\left({\mathcal{X}}_N^U\cap \mathcal{D}\ne \varnothing \right)=\mathbb{E}\left({\left(1-\frac{|\mathcal{D}|}{|\mathcal{X}|}\right)}^{N^U}\right).$$

By the independence of $\{z_n\}$, we obtain

$$\begin{array}{l}\mathbb{E}\left({\left(1-\frac{|\mathcal{D}|}{|\mathcal{X}|}\right)}^{N^U}\right)=\mathbb{E}{\displaystyle \prod _{n=1}^N{\left(1-\frac{|\mathcal{D}|}{|\mathcal{X}|}\right)}^{z_n}}\\ ={\left(\left(1-\frac{|\mathcal{D}|}{|\mathcal{X}|}\right)b|\mathcal{X}|+1-b|\mathcal{X}|\right)}^N.\end{array}$$

Then,

$$1-P\left({\mathcal{X}}_N^U\cap \mathcal{D}\ne \varnothing \right)=e^{N\log \left(1-b|\mathcal{D}|\right)}.$$(9)

Combining (9) with Lemmas 6 and 7, we obtain

$$\begin{array}{l}P\left({\widehat{\mathbf{x}}}_N^{*}\in B\left({\mathbf{x}}^{*},{\rho }_0\right)\right)\ge P\left(h_{\mathcal{X}}\left({\mathcal{X}}_N^U\right)<c_w\right)+P\left({\mathcal{X}}_N^U\cap \mathcal{D}\ne \varnothing \right)-1\\ \ge 1-{\overline{c}}_w{e}^{-\frac{1}{2}{b}^2{|\mathcal{X}|}^{2}N}-e^{N\log \left(1-b|\mathcal{D}|\right)},\end{array}$$

whence the desired result follows. □

From Lemma 8, it follows that

We then obtain (7) by a direct application of the Borel-Cantelli lemma.

6.2. Convergence Rate Around ${\mathbf{x}}^{*}$

From the preceding, we know that there is an a.s. finite random number $\overline{N}$ such that ${\widehat{\mathbf{x}}}_N^{*}\in B({\mathbf{x}}^{*},{\rho }_0)$ whenever $N\ge \overline{N}$. Other than this, we will not refer to any result from Section 6.1. The flow of the proof is somewhat complicated; the relationships between various intermediate steps are given in Figure 3. However, one can think of the proof as consisting of two major tasks. The first task is to prove a rate on $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$ inside the ball $B({\mathbf{x}}^{*},{\rho }_0)$. In Figure 3, the end result of this task is represented by (17), and the components that lead into it in the diagram represent intermediate results mainly focusing on the properties of the Hessian of f obtained from local strong convexity.

The second task, represented in Figure 3 by Lemma 10 and (18), improves on the rate obtained from the first task by showing more explicitly how it depends on the local search radius r_n. Theorem 1 is then concluded by setting r_n as indicated in (5).

We now begin the first task. Keeping the notation introduced previously, we begin by applying the strong law of large numbers to N^U. Because

we have $\frac{N^U-\mathbb{E}(N^U)}{N}\to 0$ as $N\to \infty$, that is, $N^U\cong |\mathcal{X}|bN$ a.s. Applying Lemma 5, we obtain

Then,

almost surely as $N\to \infty$.

Our immediate goal now is to relate $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$ to $h_{B({\mathbf{x}}^{*},{\rho }_0)}({\mathcal{X}}_N)$, thus yielding a rate through (10). This connection can be made through the function f, because applying Lemma 1 on $B({\mathbf{x}}^{*},{\rho }_0)$ gives

for $\mathbf{\mu }\in {\mathcal{N}}^d$ with $|\mathbf{\mu }|\le 2\le k$. We can quantify the approximation error of ${\widehat{f}}_N$ and its Hessian as

and

We will now work to connect these results to $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$ using a second-order Taylor expansion. Letting $H_f(\mathbf{x})$ and $H_{{\widehat{f}}_N}(\mathbf{x})$ be the Hessian matrices of (respectively) f and ${\widehat{f}}_N$ at x, we can write

for some $\breve{\mathbf{x}}$ on the segment joining ${\widehat{\mathbf{x}}}_N^{*}$ and ${\mathbf{x}}^{*}$, that is,

Similarly, we have

for some ${\breve{\mathbf{x}}}_N$ on the segment joining ${\widehat{\mathbf{x}}}_N^{*}$ and ${\mathbf{x}}^{*}$. Adding (11) and (12), we obtain the upper bound

A lower bound involving $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$ can be obtained via

The constant

satisfies ${\lambda }_0>0$ by the assumption that f is locally strongly convex inside $B({\mathbf{x}}^{*},{\rho }_0)$. A further lower bound can be obtained from the following technical lemma (whose proof is deferred to the Appendix).

Lemma 9.

For a positive definite d × d matrix A satisfying $|A_{ij}|\le t_A$ for all i, j, then ${\lambda }_{\max }(A)\le d·t_A$.

Applying Lemma 9 to (14) leads to

Combining (13) with (15) yields

almost surely when $N\ge \overline{N}$ and $h_{B({\mathbf{x}}^{*},{\rho }_0)}$ is small enough. Considering the decreasing rate of the mesh norm obtained from (10), this in turn implies

The rate in (17) can be improved by narrowing the local search region over time, as long as both ${\mathbf{x}}^{*}$ and ${\widehat{\mathbf{x}}}_N^{*}$ are elements of each region in the sequence, and we collect infinitely many samples from these regions. The decay rate of $\{r_n\}$ begins to play an important role in ensuring that these conditions hold. If r_n decays too slowly, we will undersample in the local regions; if r_n decays too quickly, our local regions may fail to cover ${\mathbf{x}}^{*}$.

It turns out that the best possible rate for r_n is $c_r{\left(\frac{\log (bn)}{bn}\right)}^{\frac{k}{2d}}$ for some large enough constant $c_r>0$. With this specific choice, the following result is obtained.

Lemma 10.

Suppose that x_n is sampled from the density g_n defined in (3). Let $r_n=c_r{\left(\frac{\log (bn)}{bn}\right)}^{\frac{k}{2d}}$ for $c_r>0$ and define

$$S_n=S\left({\mathbf{x}}^{*},2c_b{\left(\frac{\log \left(bn\right)}{bn}\right)}^{\frac{k}{2d}}\right).$$

Then, for c_r large enough, there exists a constant $c_b>0$ and a random integer $n_r>\max \{\overline{N},1\}$ such that, when $N>n_r$,

$$S_N\subseteq {\cap }_{n=n_r}^{N}S\left({\widehat{\mathbf{x}}}_{n-1}^{*},r_n\right)$$

almost surely.

Proof.

We have already obtained (17). From this result, there exist a constant c_b and an integer n_r such that, for all $n\ge n_r$,

$$\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert \le c_b{\left(\frac{\log \left(bn\right)}{bn}\right)}^{\frac{k}{2d}}$$

and (by the equivalence of vector norms),

$$\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}{\Vert }_{\infty }\le c_b{\left(\frac{\log \left(bn\right)}{bn}\right)}^{\frac{k}{2d}}.$$

Let c_r be large enough to satisfy, uniformly over the possible values of n_r, the inequality

$$r_n>c_b{\left(\frac{\log \left(b\left(n-1\right)\right)}{b\left(n-1\right)}\right)}^{\frac{k}{2d}}+c_b{\left(\frac{\log \left(bn\right)}{n}\right)}^{\frac{k}{2d}},\hspace{1em}\hspace{1em}\forall n\ge n_r.$$

Then, we have $r_n>\Vert {\widehat{\mathbf{x}}}_{n-1}^{*}-{\mathbf{x}}^{*}{\Vert }_{\infty }+\Vert {\widehat{\mathbf{x}}}_n^{*}-{\mathbf{x}}^{*}{\Vert }_{\infty }$ a.s., whence ${\mathbf{x}}^{*}\in S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)$ and $S_n\subseteq S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)$. Similarly, for all n, N satisfying $N\ge n\ge n_r>\max \{\overline{N},1\}$, we have $r_n>\Vert {\widehat{\mathbf{x}}}_{n-1}^{*}-{\mathbf{x}}^{*}{\Vert }_{\infty }+\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}{\Vert }_{\infty }$ a.s., whence we obtain $S_N\subseteq S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)$, as required. □

From Lemma 10, it follows that the ball

inscribed in S_N satisfies $B_N\subseteq {\cap }_{n=n_r}^{N}S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)$. Additionally, from the results of Section 6.1, we have ${\widehat{\mathbf{x}}}_N^{*}\in B_N$ almost surely whenever $N\ge \overline{n}$, where $\overline{n}\ge n_r$ is an a.s. finite random number.

From (3), we have $g_N(\mathbf{x})=t_N$ a.s. for all $\mathbf{x}\in B_N$ and $N\ge \overline{n}\ge n_r$. Let

Any $\mathbf{x}\in {\mathcal{X}}^{B_N}$ was drawn from a uniform distribution on B_N. The volume of B_N satisfies

for some constant c_B.

Conditional on $\overline{n}$, take $N\ge \overline{n}$ and let $z_{\overline{n}}^{\prime },\dots ,z_N^{\prime }$ be independent Bernoulli random variables, with the success probability of $z_n^{\prime }$ being $t_n|B_N|$ for $n=\overline{n},\dots ,N$. Sampling x_n from the density g_n is equivalent to sampling uniformly on B_N if $z_n^{\prime }=1$, and sampling from a distribution with density function $\mathbf{x}\mapsto g_n(\mathbf{x})-t_n{1}_{B_N}(\mathbf{x})$ if $z_n^{\prime }=0$. Then, because $N^{B_N}={\sum }_{n=\overline{n}}^N{z}_n^{\prime }$, it follows by the strong law of large numbers and the a.s. finiteness of $\overline{n}$ that

Note that B_N is shape-regular because we can use a shifting and scaling mapping on $B(0,1)$ into B_N. The scaling factor is precisely $c_b{\left(\frac{\log (bN)}{bN}\right)}^{\frac{k}{2d}}$. The mesh norm of a set $\mathcal{D}$ will be invariant if $\mathcal{D}$ is translated and the sampling density is changed accordingly. If $\mathcal{D}$ is scaled isotropically and the sampling density is changed accordingly, the mesh norm will be scaled by the same scaling factor.

Conditional on $N^{B_N},\hspace{0.17em}{\mathcal{X}}^{B_N}$ contains $N^{B_N}$ points sampled uniformly on B^N. By Lemma 4 and the properties of the bijection $L^{B_N}$ obtained from shape regularity, we have

with each relation holding a.s., as long as $t_N·|B_N|·N$ does not converge to zero. We will now prove this last assertion.

Because ${\widehat{\mathbf{x}}}_N^{*},{\mathbf{x}}^{*}\in B_N$ when $N\ge \overline{n}$, with $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert =O\left({\left(\frac{\log (bN)}{bN}\right)}^{\frac{k}{2d}}\right)$ almost surely, we can obtain the a.s. inequality

similarly to (16), for $h_{B_N}$ small enough and $N\ge \overline{n}$, by repeating the steps of the analysis done for $B({\mathbf{x}}^{*},{\rho }_0)$ using the strong convexity of f on B_N. The quantity λ_N in (18) is defined as ${\lambda }_N={\inf }_{\mathbf{x}\in B_N}{\lambda }_{\min }(H_f(\mathbf{x}))$.

Recall that $h_{B_N}({\mathcal{X}}_N)\le h_{B_N}({\mathcal{X}}^{B_N})$, and observe that ${\lambda }_N\to {\lambda }_{\min }(H_f({\mathbf{x}}^{*}))$ as $N\to \infty$. It follows, analogously to (17), that

With $r_N=c_r{\left(\frac{\log (bN)}{bN}\right)}^{\frac{k}{2d}}$, the rate of t_N follows from (4), and is given by

From this it follows that

which means that $t_N·|B_N|·N=O(N)$, as required. Finally, substituting (20) into (19) yields

as claimed by Theorem 1.

Equation (19) illustrates the benefit of shrinking r_N. If the size of the local search region is constant (i.e., $r_N\equiv r$), we can see that t_N will also be a constant, and the convergence rate will reduce to (17). The bound is tightened when t_N becomes larger, with the tightest bound obtained when r_N follows (5). At the same time, we cannot make r_N vanish more quickly than the rate obtained in (17) using the properties of the mesh norm.

7. Proof of Theorem 2

The flow of the proof is quite complicated and requires many interlocking parts, whose relationships are illustrated in Figure 4. Three main tasks, corresponding to Sections 7.1 to 7.3, can be identified:

• In Section 7.1, the main goal is to derive a concentration inequality for the mesh norm defined on ${[0,1]}^d$ under uniform sampling. This task is represented by the top part of Figure 4, above the first dashed line. The analysis is self-contained; only the final result, stated in Lemma 15, is used for subsequent tasks.

• In Section 7.2, we then derive a similar concentration inequality for the mesh norm, but this time on a general domain $\mathcal{X}$ and under ε-greedy sampling. The final result, shown in Figure 4 above the second dashed line, is Lemma 18.

• Finally, in Section 7.3, the concentration inequalities obtained for the mesh norm are converted into analogous results for the estimation error $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$, arriving at Theorem 2.

7.1. Concentration of Mesh Norm on ${[0,1]}^d$ Under Uniform Sampling

In this section, we derive a concentration inequality on the mesh norm $h_{{[0,1]}^d}({\mathcal{X}}_0^U)$ with ${\mathcal{X}}_0^U$ being a set of N uniformly sampled design points. Specifically, we bound the probability that, at a fixed time N, the mesh norm does not yet obey the asymptotic pathwise rate obtained from Lemma 3.

To study the behavior of the mesh norm in finite time, we partition the domain into a finite number of subsets, then grow the size of the partition as more design points are sampled at a suitably chosen rate. The growth rate of the partition size can be related to the declining behavior of the mesh norm. This analysis is based on combinatorics and largely self-contained, with only the final result being used in further tasks.

Specifically, we consider $m=K^d$ sets ${\zeta }_{m,1},\dots ,{\zeta }_{m,m}$ of the form $\left[\frac{i_1}{K},\frac{i_1+1}{K}\right]\times \left[\frac{i_2}{K},\frac{i_2+1}{K}\right]\times \dots \times \left[\frac{i_d}{K},\frac{i_d+1}{K}\right]$ with $i_j=0,1,\dots ,K-1$ for $j=1,\dots d$. We then have ${\displaystyle {\cup }_i{\zeta }_{m,i}}={[0,1]}^d$ and ${\zeta }_{m,i}\cap {\zeta }_{m,j}$ has zero Lebesgue measure for any $i\ne j$ (also under the measure induced by the uniform distribution). Each subset ${\zeta }_{m,i}$ also has the same volume under either measure.

Given N design points, let M be the maximum number of subsets into which the domain can be partitioned (according to the method described previously) such that each subset contains at least one design point in its interior. The probability that any design points will fall exactly on the joint boundary of two adjacent subsets is zero. Formally,

Thus, M is a random variable that takes positive integer values and is dependent on N. We are interested in the increasing rate of M as N becomes large.

First, for fixed positive integers m and N, let $\mathbf{W}(N,m)=(W_1,\dots ,W_m)$ be a random vector following a multinomial distribution with parameters (N, m) and probability vector $\frac{1}{m}·\mathbf{\delta }$, where ${\delta }_j=1$ for each component j. The following result (proved in the Appendix) calculates the probability that no subsets will be empty.

Lemma 11.

For any N and m,

$$P\left(W\left(N,m\right)\succcurlyeq \delta \right)=\frac{1}{m^N}{\displaystyle \sum _{j=0}^m\left(\begin{array}{c}m\\ j\end{array}\right)}{\left(-1\right)}^{m-j}{j}^N.$$(21)

Now observe that

because the maximum partition size is at least m if and only if, in the partition whose size is m, each subset contains at least one design point. In another way of writing, $P(M\ge m)=\frac{m!}{m^N}\left\{\begin{array}{c}N\\ m\end{array}\right\}$, where $\left\{\begin{array}{c}N\\ m\end{array}\right\}$ represents a Stirling number of the second kind, defined for general m, n as

By following the asymptotic analysis of Stirling numbers (Temme 1993), we can obtain the following increasing rate of M. The next result (whose proof is deferred to the Appendix) then derives a more useful equivalent form for this approximation.

Lemma 12.

(Temme 1993). The Stirling number of the second kind $\left\{\begin{array}{c}n\\ m\end{array}\right\}$ can be approximated (uniformly in m, and asymptotically as $n\to \infty$) as

$$\left\{\begin{array}{c}n\\ m\end{array}\right\}\cong e^T{m}^{n-m}\tau \left(\frac{n}{m}-1\right)\left(\begin{array}{c}n\\ m\end{array}\right),$$(22)

where

$$\begin{array}{l}T=\gamma \left({\mathbf{y}}_0\right)-n+m+\left(n-m\right)\log \left(\frac{n}{m}-1\right),\\ \tau \left(s\right)=\sqrt{\frac{ms}{n\left({\mathbf{y}}_0-s\right)}},\\ \gamma \left(\mathbf{y}\right)=-n\log \mathbf{y}+m\log \left(e^{\mathbf{y}}-1\right),\end{array}$$

and y₀ is the solution of $\frac{m}{n}\mathbf{y}=1-e^{-\mathbf{y}}$.

Lemma 13.

The Stirling number of the second kind $\left\{\begin{array}{c}n\\ m\end{array}\right\}$ can be approximated (uniformly in m, and asymptotically as $n\to \infty$) as

$$\left\{\begin{array}{c}n\\ m\end{array}\right\}\cong e^{m-n}{\left(\frac{\frac{n}{m}-1}{\frac{n}{m}-w}\right)}^{n-m}\frac{m^{n-m}}{w^m}\sqrt{\frac{\frac{n}{m}-1}{\frac{n}{m}\left(1-w\right)}}\left(\begin{array}{c}n\\ m\end{array}\right),$$

where $w=-W_0\left(-\frac{n}{m}{e}^{-\frac{n}{m}}\right)$ and W₀ is the upper branch of the Lambert W function.

Now, we can apply Lemma 13 to obtain an asymptotic rate for $P(M\ge m)$ that is expressible in closed form as a function of N and m. Consequently, as N grows large, we can choose a suitable growth rate for m that would allow us to achieve a partition of size m w.p.1. The following bound (proved in the Appendix) provides a sufficient condition for this growth rate to be suitable.

Lemma 14.

If $m\le N$ and $m,N\to \infty$ with $\frac{N}{m}\to \infty$, then $P(M<m)\le c_t^{\prime }\frac{1}{N-m}$ for some constant $c_t^{\prime }>0$.

With these technical results, we return to the decreasing rate of the mesh norm on ${[0,1]}^d$ under uniform sampling. The final result of this section follows fairly straightforwardly from the preceding.

Lemma 15.

Let ${\mathcal{X}}_0^U$ be a set of N design points sampled independently from the uniform distribution on ${[0,1]}^d$. There exist constants $c_h,c_t>0$ such that, for any n > 1,

$$P\left(h_{{[0,1]}^d}\left({\mathcal{X}}_0^U\right)>c_h{\left(\frac{\log N}{N}\right)}^{\frac{1}{d}}\right)\le c_t\frac{1}{N}.$$

Proof.

There exists a constant $c_h^{\prime }$ such that the statement $M\ge m$ (for any m) implies that $h_{{[0,1]}^d}({\mathcal{X}}_0^U)\le c_h^{\prime }{m}^{-\frac{1}{d}}$. To match the almost sure convergence rate of the mesh norm, choose $m^{*}\cong \frac{N}{\log N}$. Then, there exist constants $c_h,c_t>0$ such that, for any N > 1,

$$\begin{array}{l}P\left(h_{{[0,1]}^d}\left({\mathcal{X}}_0^U\right)>c_h{\left(\frac{\log N}{N}\right)}^{\frac{1}{d}}\right)\le P\left(M<m^{*}\right)\\ \le c_t^{\prime }\frac{1}{N-m^{*}}\\ \le c_t\frac{1}{N},\end{array}$$(23)

with (23) following from Lemma 14. □

7.2. Concentration of Mesh Norm on $\mathcal{X}$ Under ε-Greedy Sampling

The final result of this section is an analog of Lemma 15 from the previous section with the general domain $\mathcal{X}$ and the design points ${\mathcal{X}}_N$ obtained from ε-greedy sampling. We move to this more general case in two steps. First, we consider nonuniform sampling while keeping the domain as ${[0,1]}^d$, and then extend the result to a general domain.

Let ${\mathcal{X}}_N^{\prime }=\{{\mathbf{x}}_1^{\prime },\dots ,{\mathbf{x}}_N^{\prime }\}$ where each ${\mathbf{x}}_n^{\prime }$ is sampled independently from a density $g_n^{\prime }$ with support ${[0,1]}^d$. Let $b_n^{\prime }={\inf }_{\mathbf{x}\in {[0,1]}^d}{g}_n^{\prime }(\mathbf{x})$ and suppose that $b_n^{\prime }>0$ and ${\sum }_{n=1}^N{b}_n^{\prime }\to \infty$ as $N\to \infty$. Let $z_n^{\prime }$ be independent Bernoulli random variables with success probabilities $b_n^{\prime }$ so that, if $z_n^{\prime }=1$, then ${\mathbf{x}}_n^{\prime }$ is sampled from a uniform distribution on ${[0,1]}^d$, and if $z_n^{\prime }=0$, then ${\mathbf{x}}_n^{\prime }$ is sampled from the density $\mathbf{x}\mapsto \frac{g_n^{\prime }(\mathbf{x})-b_n^{\prime }}{1-b_n^{\prime }}$.

Let $N^{U^{\prime }}={\sum }_{n=1}^N{z}_n^{\prime }$ and ${\overline{b}}_N^{\prime }=\frac{1}{N}{\sum }_{n=1}^N{b}_n^{\prime }$. A direct application of Hoeffding’s inequality yields

In later proofs, this will be combined with the following technical result (proved in the Appendix).

Lemma 16.

If $|N^{U^{\prime }}-{\overline{b}}_N^{\prime }N|\le \frac{1}{2}{\overline{b}}_N^{\prime }N-1$, then

$${\left(\frac{\log N^{U^{\prime }}}{N^{U^{\prime }}}\right)}^{\frac{1}{d}}\le 2^{\frac{1}{d}}{\left(\frac{\log \left({\overline{b}}_N^{\prime }N\right)}{{\overline{b}}_N^{\prime }N}\right)}^{\frac{1}{d}}.$$

With the previous preliminaries, we may establish the following concentration inequality for the mesh norm on ${[0,1]}^d$ under nonuniform sampling. The main idea is to apply Lemma 15 to the part of the data that comes from a uniform distribution.

Lemma 17.

There exist constants $c_{h,1},c_{t,1}$ such that, for large enough N,

$$P\left(h_{{[0,1]}^d}\left({\mathcal{X}}_N^{\prime }\right)>c_{h,1}{\left(\frac{\log \left({\overline{b}}_N^{\prime }N\right)}{{\overline{b}}_N^{\prime }N}\right)}^{\frac{1}{d}}\right)\le c_{t,1}\frac{1}{{\overline{b}}_N^{\prime }N}.$$

Proof.

For notational compactness, denote the event

$$E=\left\{\left|N^{U^{\prime }}-{\overline{b}}_N^{\prime }N\right|\le \frac{1}{2}{\overline{b}}_N^{\prime }N-1\right\}.$$

For any arbitrary positive value of $c_{h,1}$, we derive

$$\begin{array}{l}\hspace{0.17em}P\left(h_{{\left[0,1\right]}^d}\left({\mathcal{X}}_N^{\prime }\right)>c_{h,1}{\left(\frac{\log \left({\overline{b}}_N^{\prime }N\right)}{{\overline{b}}_N^{\prime }N}\right)}^{\frac{1}{d}}\right)\\ \le \mathbb{E}\left[P\left(h_{{\left[0,1\right]}^d}\left({\mathcal{X}}_N^{\prime }\right)>c_{h,1}{\left(\frac{\log \left({\overline{b}}_N^{\prime }N\right)}{{\overline{b}}_N^{\prime }N}\right)}^{\frac{1}{d}}\hspace{0.17em}|\hspace{0.17em}{N}^{U^{\prime }},E\right)\hspace{0.17em}|\hspace{0.17em}E\right]+P\left(E^c\right)\\ \le \mathbb{E}\left[P\left(h_{{\left[0,1\right]}^d}\left({\mathcal{X}}_N^{\prime }\right)>c_{h,1}{2}^{-\frac{1}{d}}{\left(\frac{\log \left(N^{U^{\prime }}\right)}{N^{U^{\prime }}}\right)}^{\frac{1}{d}}\hspace{0.17em}|\hspace{0.17em}{N}^{U^{\prime }},E\right)\hspace{0.17em}|\hspace{0.17em}E\right]+o\left(\frac{1}{{\overline{b}}_N^{\prime }N}\right),\end{array}$$(25)

where (25) is due to Lemma 16 and (24).

When N satisfies

$$N\ge \min \{n\hspace{0.17em}:\hspace{0.17em}{\overline{b}}_{n^{\prime }}^{\prime }{n}^{\prime }>e\hspace{0.17em}\forall n^{\prime }\ge n\},$$

we have $N^{U^{\prime }}\ge \frac{1}{2}{\overline{b}}_N^{\prime }N+1>1$ on the event E. Consequently, letting $c_{h,1}=2^{\frac{1}{d}}{c}_h$, where c_h is the constant obtained from Lemma 15, we obtain

$$\begin{array}{l}P\left(h_{{[0,1]}^d}\left({\mathcal{X}}_N^{\prime }\right)>c_{h,1}{\left(\frac{\log \left({\overline{b}}_N^{\prime }N\right)}{{\overline{b}}_N^{\prime }N}\right)}^{\frac{1}{d}}\right)\le \mathbb{E}\left(c_t\frac{1}{N^{U^{\prime }}}\hspace{0.17em}|\hspace{0.17em}E\right)+o\left(\frac{1}{{\overline{b}}_N^{\prime }N}\right)\\ \le c_t\frac{1}{\frac{1}{2}{\overline{b}}_N^{\prime }N+1}+o\left(\frac{1}{{\overline{b}}_N^{\prime }N}\right)\\ \le c_{t,1}\frac{1}{{\overline{b}}_N^{\prime }N},\end{array}$$(26)

where (26) applies Lemma 15 to (25), and $c_{t,1}$ is suitably chosen to dominate the second term in (26). □

The concentration inequality can now be generalized to the domain $\mathcal{X}$. Because we assume that $\mathcal{X}$ is shape regular, there exists a continuously differentiable bijection $L_{\mathcal{X}}:{[0,1]}^d\to \mathcal{X}$ whose Jacobian has nonzero determinant everywhere on its domain. The proof works with the inverse image of ${\mathcal{X}}_N$ under $L_{\mathcal{X}}$, similar to the proof of Lemma 5.

Lemma 18.

There exist constants $c_{h,2},c_{t,2}>0$ such that, for all large enough N,

$$P\left(h_{\mathcal{X}}\left({\mathcal{X}}_N\right)>c_{h,2}{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{1}{d}}\right)\le c_{t,2}\frac{1}{bN}.$$(27)

Proof.

As in the proof of Lemma 5, we observe that the mapping $L_{\mathcal{X}}$ is Lipschitz. The Lipschitz condition implies

$$h_{\mathcal{X}}\left({\mathcal{X}}_N\right)\le c_{L_{\mathcal{X}}}{h}_{{[0,1]}^d}\left(L_{\mathcal{X}}^{-1}\left({\mathcal{X}}_N\right)\right).$$

For each n, define a density

$$g_n^{\prime }\left(\mathbf{y}\right)=g_n\left(L_{\mathcal{X}}\left(\mathbf{y}\right)\right)\left|\det \left(\frac{d{L}_{\mathcal{X}}\left(\mathbf{y}\right)}{d\mathbf{y}}\right)\right|.$$

Then, letting $b_n^{\prime }={\inf }_{\mathbf{y}\in {[0,1]}^d}{g}_n^{\prime }(\mathbf{y})$, we have $b_n^{\prime }\ge b{c}_J$, where $c_J={\inf }_{\mathbf{y}\in {[0,1]}^d}\left|\det \left(\frac{d{L}_{\mathcal{X}}(\mathbf{y})}{d\mathbf{y}}\right)\right|$. We also have ${\overline{b}}_N^{\prime }\ge b{c}_J$, where ${\overline{b}}_N^{\prime }=\frac{1}{N}{\sum }_{n=1}^N{b}_n^{\prime }$.

For any arbitrary positive value of $c_{h,2}$, the inequalities

$$\begin{array}{l}\hspace{0.17em}P\left(h_{\mathcal{X}}\left({\mathcal{X}}_N\right)>c_{h,2}{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{1}{d}}\right)\\ \le P\left(h_{{[0,1]}^d}\left(L_{\mathcal{X}}^{-1}\left({\mathcal{X}}_N\right)\right)>\frac{c_{h,2}}{c_{L_{\mathcal{X}}}}{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{1}{d}}\right)\\ \le P\left(h_{{[0,1]}^d}\left(L_{\mathcal{X}}^{-1}\left({\mathcal{X}}_N\right)\right)>\frac{c_{h,2}}{c_{L_{\mathcal{X}}}}{c}_J^{\prime }{\left(\frac{\log \left({\overline{b}}_N^{\prime }N\right)}{{\overline{b}}_N^{\prime }N}\right)}^{\frac{1}{d}}\right),\end{array}$$

where $c_J^{\prime }$ is some constant, hold for all N satisfying

$$N\ge \min \left\{n\hspace{0.17em}:\hspace{0.17em}b{n}^{\prime }>\max \left\{e,\frac{e}{c_J}\right\}\hspace{0.17em}\forall n^{\prime }\ge n\right\}.$$(28)

Now choose $c_{h,2}=c_{h,1}\frac{c_{L_{\mathcal{X}}}}{c_J^{\prime }}$, where $c_{h,1}$ is the constant obtained from Lemma 17. It follows that, for all N satisfying (28), we have

$$P\left(h_{\mathcal{X}}\left({\mathcal{X}}_N\right)>c_{h,2}{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{1}{d}}\right)\le c_{t,1}\frac{1}{{\overline{b}}_N^{\prime }N}\le \frac{c_{t,1}}{c_J}\frac{1}{bN},$$

where $c_{t,1}$ is the constant obtained from Lemma 17. Letting $c_{t,2}=\frac{c_{t,1}}{c_J}$ yields the desired result. □

If the Inequality (27) holds for a certain value of $c_{h,2}$, it also holds for any larger value. The same is true of the other concentration inequalities on the mesh norm that we derived throughout this section.

7.3. Concentration Inequality for Estimation Error

Finally, we connect the previously obtained results for the mesh norm to the estimation error $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$. We first take care of the situation where ${\widehat{\mathbf{x}}}_N^{*}$ is converging to ${\mathbf{x}}^{*}$ too slowly, by deriving a bound on the probability of this event. This is done by building on Lemma 8.

Then, conditional on the event $\{{\widehat{\mathbf{x}}}_N^{*}\in B({\mathbf{x}}^{*},{\rho }_0)\}$, we apply (16), which implies the existence of a constant $c_b^{\prime }$ such that, for all N,

At this point, Lemma 18 provides a probabilistic bound for the right-hand side of (29). We thus obtain the following concentration inequality on $\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert$, although it is still not the strongest one possible because the right-hand side of (29) only considers a ball $B({\mathbf{x}}^{*},{\rho }_0)$ with fixed radius. The situation is analogous to Section 6.2, where we also first derived a rate within that ball, and later improved it by incorporating the shrinking radius.

Lemma 19.

There exist constants $c_b,c_{t,3}>0$ such that, for all large enough N,

$$P\left({\widehat{\mathbf{x}}}_N^{*}\in \mathcal{X}\setminus B\left({\mathbf{x}}^{*},\min \left\{{\rho }_0,c_b{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k}{2d}}\right\}\right)\right)\le c_{t,3}\frac{1}{bN}.$$

Proof.

We derive

$$\begin{array}{l}\hspace{0.17em}P\left({\widehat{\mathbf{x}}}_N^{*}\in \mathcal{X}\setminus B\left({\mathbf{x}}^{*},\min \left\{{\rho }_0,c_b{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k}{2d}}\right\}\right)\right)\\ =P\left(\left\{\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert >c_b{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k}{2d}}\right\}\hspace{0.17em}\cup \hspace{0.17em}\left\{{\widehat{\mathbf{x}}}_N^{*}\in \mathcal{X}\setminus B\left({\mathbf{x}}^{*},{\rho }_0\right)\right\}\right)\\ \le P\left(\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert >c_b{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k}{2d}}\hspace{0.17em}|\hspace{0.17em}{\widehat{\mathbf{x}}}_N^{*}\in B\left({\mathbf{x}}^{*},{\rho }_0\right)\right)+P\left({\widehat{\mathbf{x}}}_N^{*}\in \mathcal{X}\setminus B\left({\mathbf{x}}^{*},{\rho }_0\right)\right)\\ \le P\left(\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert >c_b{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k}{2d}}\hspace{0.17em}|\hspace{0.17em}{\widehat{\mathbf{x}}}_N^{*}\in B\left({\mathbf{x}}^{*},{\rho }_0\right)\right)+o\left(\frac{1}{bN}\right),\end{array}$$

where the last line follows from Lemma 8. We then derive

$$\begin{array}{l}\hspace{0.17em}P\left(\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert >c_b{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k}{2d}}\hspace{0.17em}|\hspace{0.17em}{\widehat{\mathbf{x}}}_N^{*}\in B\left({\mathbf{x}}^{*},{\rho }_0\right)\right)\\ \le P\left(c_b^{\prime }{h}_{\mathcal{X}}^{\frac{k}{2}}\left({\mathcal{X}}_N\right)>c_b{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k}{2d}}\hspace{0.17em}|\hspace{0.17em}{\widehat{\mathbf{x}}}_N^{*}\in B\left({\mathbf{x}}^{*},{\rho }_0\right)\right)\end{array},$$(30)

where $c_b^{\prime }$ in (30) is the same value as in (29). Now, if we choose $c_b=c_b^{\prime }{c}_{h,2}^{\frac{k}{2}}$, where $c_{h,2}$ is the same as in Lemma 18, we obtain

$$P\left(\Vert {\widehat{\mathbf{x}}}_N^{*}-{\mathbf{x}}^{*}\Vert >c_b{\left(\frac{\log \left(bN\right)}{bN}\right)}^{\frac{k}{2d}}\hspace{0.17em}|\hspace{0.17em}{\widehat{\mathbf{x}}}_N^{*}\in B\left({\mathbf{x}}^{*},{\rho }_0\right)\right)\le c_{t,2}\frac{1}{bN},$$

where $c_{t,2}$ is the same as in Lemma 18. The desired result follows. □

Now, analogously to Section 6.2, we will strengthen this result by considering the local search radius r_n set according to (5) with some sufficiently large c_r. We first show that, using such a radius, the local sampling regions $S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)$ for $n\le N$ will cover

where c_b is the value obtained from Lemma 19, with sufficiently high frequency as $N\to \infty$. In other words, there will be sufficiently many iterations n in which a rectangle centered at ${\widehat{\mathbf{x}}}_n^{*}$ will cover a rectangle centered at ${\mathbf{x}}^{*}$, with the sizes of both rectangles shrinking as n grows large.

Before we proceed, we slightly relax the constant $c_{h,2}$ in Lemma 18. From (10), we have

almost surely. From Lemma 3, we know that there exists a nonrandom constant $c_{h,3}$ such that

for all large enough N (the exact threshold value of N at which this happens may be random, however). Thus, we can let $c_{h,4}=\max \{c_{h,2},c_{h,3}\}$ and replace $c_{h,2}$ in Lemma 18 by $c_{h,4}$ without changing the result. Similarly, if N is large enough for (32) to hold, then we will still have Inequality (32) if we replace $c_{h,3}$ by $c_{h,4}$.

Now, consider the set ${\mathcal{X}}_n$ of the first n data points. We observe that $h_{\mathcal{X}}({\mathcal{X}}_n)$ is decreasing in n. Thus, if ${\widehat{\mathbf{x}}}_N^{*}\in B({\mathbf{x}}^{*},{\rho }_0)$ and $h_{\mathcal{X}}({\mathcal{X}}_{n^{\prime }})\le c_{h,4}{\left(\frac{\log (b{n}^{\prime })}{b{n}^{\prime }}\right)}^{\frac{1}{d}}$ for some $n^{\prime }$, then

for all $n^{\prime }\le n\le N$. The value of $c_b^{\prime }$ in (33) is the same as in (29), and the value of c_b in (34) is obtained from Lemma 19 with $c_{h,4}$ replacing $c_{h,3}$ as discussed previously. Using similar arguments as in the proof of Lemma 10, for suitably chosen c_r in (5), we have

where $S_N,S_{n^{\prime }+1}$ are as in (31). It follows that the number of time stages $n\le N$ in which S_N is covered by $S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)$ is bounded below by the number of times that

is achieved.

The remainder of our analysis is based on this concept of coverage. Whenever S_N is covered by $S({\widehat{\mathbf{x}}}_{n-1}^{*},r_n)$, there is a certain probability (varying with n) that a sample drawn from the local search region falls in S_N. Recalling that the ε-greedy algorithm conducts local search by sampling uniformly inside the local region, one can then interpret the distribution of that sample as a mixture of uniforms. In other words, whenever S_N is contained in the local search region, there is a certain probability that the sample will be drawn from a uniform distribution on S_N. One can then make use of Lemma 15 to study the convergence rate of the mesh norm in S_N, a neighborhood of the global optimum.

To make these arguments, we first need to ensure that coverage occurs, that is, the probability that S_N is never covered in the above-described manner is sufficiently small. Let n be some fixed integer large enough to satisfy Inequality (28) and define

with $N^S={\sum }_{n^{\prime }=n}^N{z}_{n^{\prime }}^S$. Thus, N^S counts the number of time stages $n\le n^{\prime }\le N$ in which $S({\widehat{\mathbf{x}}}_{n^{\prime }-1}^{*},r_{n^{\prime }})$ covers S_N. Recall that (32) holds for all large enough $n^{\prime }$, although the exact threshold after which this occurs may be random. Therefore, $z_{n^{\prime }}^S=1$ for all large enough $n^{\prime }$, whereupon

holds almost surely. We will discard a suitable set of measure zero from the outcome space so that (35) can be assumed to always hold, and we do not have to keep conditioning on this event in our analysis.

With this, let $N\ge n$ be a fixed value satisfying (28). Then,

where the last line follows by Lemma 18.

Now, we study the portion of the design points that can be viewed as originating from a uniform distribution on S_N. Let $z_{n^{\prime }}^{\prime }$, for $n\le n^{\prime }\le N$, be independent Bernoulli random variables with success probabilities $t_{n^{\prime }}|S_N|$. Then, $N^{S,U}={\sum }_{n^{\prime }=n}^N{z}_{n^{\prime }}^S{z}_{n^{\prime }}^{\prime }$ is the cardinality of the subset ${\mathcal{X}}_N^{S,U}$ of the data that was sampled uniformly on S_N. Note that $z_{n^{\prime }}^S$ and $z_{n^{\prime }}^{\prime }$ are independent. We introduce the notation n_k, for $k=1,2,\dots ,N^S$, to represent those time stages $n^{\prime }$ for which $z_{n^{\prime }}^S=1$, that is,

and use the notation $\mathcal{N}={\{n_k\}}_{k=1}^{N^S}$ to denote the entire sequence of such time stages (all of which are random variables).

Recall from Lemma 15 that, if ${\mathcal{X}}_0^U$ is a set of independent and identically distributed (i.i.d.) samples from the uniform density on ${[0,1]}^d$, then

for some c_h, c_t. As we have observed previously in Section 6.2, when any set $\mathcal{D}$ is scaled isotropically and the sampling density is also appropriately scaled, the mesh norm will be changed by the same scaling factor. For this reason, we can apply Inequality (37) to the mesh norm on S_N by treating it as a scaled mesh norm on ${[0,1]}^d$. That is, if we suppose that the value of $N^{S,U}$ is given, and let

for notational convenience, we then have

Recall from (20) that $t_n=b+\frac{1-b|\mathcal{X}|}{r_n^d}=O\left({\left(\frac{\log (bn)}{bn}\right)}^{-\frac{k}{2}}\right)$. Now, we define

and repeat the proof of Lemma 16 to obtain the inequality

under the condition that $N^{S,U}>\frac{1}{2}{\overline{t}}^S·|S_N|·N^S$. This allows us to improve the lower bound on $h_{S_N}({\mathcal{X}}_N^{S,U})$ in the event whose probability is computed in (38) to match the asymptotic a.s. rate of $h_{S_N}$. The resulting inequality can then be directly connected to the estimation error.

First, we give a technical lemma characterizing the probability that the condition required for (39) is not satisfied.

Lemma 20.

The following statement holds:

$$P\left(N^{S,U}<\frac{1}{2}{\overline{t}}^S|S_N|N^S+1\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right)=O\left(\frac{1}{{\left(t_N|S_N|\right)}^{2}N}\right).$$

Proof.

Conditional on $\mathcal{N},\hspace{0.17em}{N}^{S,U}$ is a sum of independent Bernoulli random variables. Applying Hoeffding’s concentration inequality, we derive

$$\begin{array}{l}\hspace{0.17em}P\left(N^{S,U}<\frac{1}{2}{\overline{t}}^S|S_N|N^S+1\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right)\\ \le P\left(|N^{S,U}-{\overline{t}}^S·|S_N|·N^S|>\frac{1}{2}{\overline{t}}^S|S_N|N^S-1\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right)\\ =\mathbb{E}\left[P\left(|N^{S,U}-{\overline{t}}^S·|S_N|·N^S|>\frac{1}{2}{\overline{t}}^S|S_N|N^S-1\hspace{0.17em}|\hspace{0.17em}\mathcal{N},E_N^{\mathbf{x}},N^S\ge 1\right)\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right]\\ \le 2\mathbb{E}\left[e^{-\frac{2}{N^S}{\left(\frac{1}{2}{\overline{t}}^S|S_N|N^S-1\right)}^2}\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right]\\ \mathrm{\lesssim }\mathbb{E}\left[\frac{1}{{\left({\overline{t}}^S|S_N|\right)}^2{N}^S}\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right]\\ =\mathbb{E}\left[\frac{1}{{\left({\overline{t}}^S\right)}^2{N}^S}\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right]\frac{1}{{|S_N|}^2}.\end{array}$$(40)

Provided that $N^S\ge 1$ and ${\lim \inf }_{N\to \infty }\frac{N^S}{N}=1$, the inequality $N^S\le n_{N^S}$ implies that $n_{N^S}\cong N$, whence, applying (20), we obtain

$$t_N\cong t_{n_{N^S}}\cong {\overline{t}}^S,$$

whence ${\lim \sup }_{N\to \infty }\frac{t_N}{{\overline{t}}^S}=1$ and ${\lim \sup }_{N\to \infty }\frac{N}{N^S}=1$. Therefore,

$$\underset{N\to \infty }{\lim \sup }{\left(\frac{t_N}{{\overline{t}}^S}\right)}^2\frac{N}{N^S}\le \underset{N\to \infty }{\lim \sup }{\left(\frac{t_N}{{\overline{t}}^S}\right)}^2·\underset{N\to \infty }{\lim \sup }\frac{N}{N^S}=1.$$

Because ${\overline{t}}^S$ and N^S are all strictly positive random variables when $N^S\ge 1$, we can apply Fatou’s lemma to obtain

$$\begin{array}{l}\underset{N\to \infty }{\lim \sup }\mathbb{E}\left[\frac{t_N^{2}N}{{\left({\overline{t}}^S\right)}^2{N}^S}\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right]\le \mathbb{E}\left[\underset{N\to \infty }{\lim \sup }\frac{t_N^{2}N}{{\left({\overline{t}}^S\right)}^2{N}^S}\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right]\le 1,\end{array}$$

which yields a further bound on (40) because of the relation

$$\mathbb{E}\left[\frac{1}{{\left({\overline{t}}^S\right)}^2{N}^S}\hspace{0.17em}|\hspace{0.17em}{E}_N^{\mathbf{x}},N^S\ge 1\right]\mathrm{\lesssim }\frac{1}{t_N^{2}N},$$

for any N satisfying (28). This completes the proof. □

Now, letting $c_{h,5}=2^{\frac{1}{d}}{c}_h$, where c_h is the value obtained from (38), and also letting $E_N^S=\{N^S\ge 1\}$ for notational convenience, we can derive

where (41) is obtained by applying (39), and (42) follows from the definition of $c_{h,5}$.

We then derive the bound

In this derivation, (44) is obtained by applying (42) to the first term of (43), and Lemma 20 to the second term. Then, (45) is due to (38).

By repeating the proof of Lemma 20, we can derive

for any N satisfying (28). Because $t_N|S_N|=O(1)$, this yields

Similarly to (19), when ${\widehat{\mathbf{x}}}_N^{*}\in B({\mathbf{x}}^{*},{\rho }_0)\cap S_N$, we can bound the error of ${\widehat{\mathbf{x}}}_N^{*}$ by the local mesh norm on S_N as

It follows from the preceding that there exists a constant $c_{b,1}$ such that

Using the fact that $t_N\cong {\overline{t}}^S$ and $N^S\cong N$, there exists a constant $c_{b,2}$ such that

shown by applying Lemma 19 together with (36) and (46). By combining this bound with the preceding analysis on the rate of t_N, we obtain constants $c,c^{\prime }$ for which

for all N satisfying (28). This condition on N depends only on b and the constant c_J, which itself depends only on the domain $\mathcal{X}$ and not on f.

8. Conclusion

We presented new convergence rate results for global optimization using radial basis function interpolation and an ε-greedy algorithm that randomly chooses between uniform sampling on the entire domain versus uniform sampling on a local neighborhood of the current estimated minimizer. This sampling method is simple to implement but captures the key distinction between local and global search that is present in many other algorithms that are not amenable to theoretical analysis. We find that convergence rates are improved when the size of the local search region is made to shrink over time at a suitable rate, that is, local search concentrates around the estimated minimizer over time.

The theory of RBF interpolation relies on a bound connecting the interpolation error to the distribution of the design points on the domain. The overall structure of our proof could potentially be repeated for other interpolation models as long as those models also had bounds of this kind (not necessarily the same order) available. The bounds improve when the design points are more evenly spread out; thus, although it is possible to obtain very similar rates for the case where global search is conducted using nonuniform sampling, it is not possible to improve the rates using this analytical technique. To obtain such improvements, it would be necessary to develop new theory that makes a closer connection between the error of RBF interpolation and the shape of the underlying function; we leave this problem for future work.