Bayesian Optimization with Exponential Convergence

Introduction

The general global optimization problem is known to be intractable if we make no further assumptions . The simplest additional assumption to restore tractability is to assume the existence of a bound on the slope of $f$. A well-known variant of this assumption is Lipschitz continuity with a known Lipschitz constant, and many algorithms have been proposed in this setting . These algorithms successfully guaranteed certain bounds on the regret. However appealing from a theoretical point of view, a practical concern was soon raised regarding the assumption that a tight Lipschitz constant is known. Some researchers relaxed this somewhat strong assumption by proposing procedures to estimate a Lipschitz constant during the optimization process .

Bayesian optimization is an efficient way to relax this assumption of complete knowledge of the Lipschitz constant, and has become a well-recognized method for solving global optimization problems with non-convex black-box functions. In the machine learning community, Bayesian optimization—especially by means of a Gaussian process (GP)—is an active research area . With the requirement of the access to the $\delta$-cover sampling procedure (it samples the function uniformly such that the density of samples doubles in the feasible regions at each iteration), de Freitas et al. recently proposed a theoretical procedure that maintains an exponential convergence rate (exponential regret). However, as pointed out by Wang et al. , one remaining problem is to derive a GP-based optimization method with an exponential convergence rate without the $\delta$-cover sampling procedure, which is computationally too demanding in many cases.

In this paper, we propose a novel GP-based global optimization algorithm, which maintains an exponential convergence rate and converges rapidly without the $\delta$-cover sampling procedure.

Gaussian Process Optimization

where $\mu(x_{N+1}|{\cal D}_{N})=m(x_{N+1})+\mathbf{k}^{T}\mathbf{K}^{-1}(\mathbf{f(x}_{1:N})-\mathbf{m(x}_{1:N}))$ and $\sigma^{2}(x_{N+1}|{\cal D}_{N})=\kappa(\mathbf{x}_{N+1},\mathbf{x}_{N+1})-\mathbf{k}^{T}\mathbf{K}^{-1}\mathbf{k}$. One advantage of GP is that this closed-form solution simplifies both its analysis and implementation.

To use a GP, we must specify the mean function and the covariance function. The mean function is usually set to be zero. With this zero mean function, the conditional mean $\mu(x_{N+1}|{\cal D}_{N})$ can still be flexibly specified by the covariance function, as shown in the above equation for $\mu$. For the covariance function, there are several common choices, including the Matern kernel and the Gaussian kernel. For example, the Gaussian kernel is defined as $\kappa({x,x^{\prime}})=\text{exp}\left(-\frac{1}{2}{(x-x^{\prime})}^{T}{\Sigma^{-1}(x-x^{\prime})}\right)$ where $\Sigma^{-1}$ is the kernel parameter matrix. The kernel parameters or hyperparameters can be estimated by empirical Bayesian methods ; see for more information about GP.

For deterministic function, de Freitas et al. recently presented a theoretical procedure that maintains exponential convergence rate. However, their own paper and the follow-up research point out that this result relies on an impractical sampling procedure, the $\delta$-cover sampling. To overcome this issue, Wang et al. combined GP-UCB with a hierarchical partitioning optimization method, the SOO algorithm , providing a regret bound with polynomial dependence on the number of function evaluations. They concluded that creating a GP-based algorithm with an exponential convergence rate without the impractical sampling procedure remained an open problem.

Infinite-Metric GP Optimization

The GP-UCB algorithm can be seen as a member of the class of bound-based search methods, which includes Lipschitz optimization, A* search, and PAC-MDP algorithms with optimism in the face of uncertainty. Bound-based search methods have a common property: the tightness of the bound determines its effectiveness. The tighter the bound is, the better the performance becomes. However, it is often difficult to obtain a tight bound while maintaining correctness. For example, in A* search, admissible heuristics maintain the correctness of the bound, but the estimated bound with admissibility is often too loose in practice, resulting in a long period of global search.

The GP-UCB algorithm has the same problem. The bound in GP-UCB is represented by UCB, which has the following property: $f(x)\leq{\cal U}(x|{\cal D})$ with some probability. We formalize this property in the analysis of our algorithm. The problem is essentially due to the difficulty of obtaining a tight bound ${\cal U}(x|{\cal D})$ such that $f(x)\leq{\cal U}(x|{\cal D})$ and $f(x)\approx{\cal U}(x|{\cal D})$ (with some probability). Our solution strategy is to first admit that the bound encoded in GP prior may not be tight enough to be useful by itself. Instead of relying on a single bound given by the GP, we leverage the existence of an unknown bound encoded in the continuity at a global optimizer.

2 Description of Algorithm

Figure 1 illustrates how the algorithm works with a simple 1-dimensional objective function. We employ hierarchical partitioning to maintain hyperintervals, as illustrated by the line segments in the figure. We consider a hyperrectangle as our hyperinterval, with its center being the evaluation point of $f$ (blue points in each line segment in Figure 1). For each iteration $t$, the algorithm performs the following procedure for each interval size:

Select the interval with the maximum center value among the intervals of the same size.

Keep the interval selected by (i) if it has a center value greater than that of any larger interval.

Keep the interval accepted by (ii) if it contains a UCB greater than the center value of any smaller interval.

If an interval is accepted by (iii), divide it along with the longest coordinate into three new intervals.

For each new interval, if the UCB of the evaluation point is less than the best function value found so far, skip the evaluation and use the UCB value as the center value until the interval is accepted in step (ii) on some future iteration; otherwise, evaluate the center value.

Repeat steps (i)–(v) until every size of intervals are considered

3 Technical Detail of Algorithm

We define $h$ to be the depth of the hierarchical partitioning tree, and $c_{h,i}$ to be the center point of the $i^{th}$ hyperrectangle at depth $h$. $N_{gp}$ is the number of the GP evaluations. Define $depth(\mathcal{T})$ to be the largest integer $h$ such that the set $\mathcal{T}_{h}$ is not empty. To compute UCB $\mathcal{U}$, we use $\varsigma_{M}=\sqrt{2\log(\pi^{2}M^{2}/12\eta)}$ where $M$ is the number of the calls made so far for ${\cal U}$ (i.e., each time we use ${\cal U}$, we increment $M$ by one). This particular form of $\varsigma_{M}$ is to maintain the property of $f(x)\leq{\cal U}(x|{\cal D})$ during an execution of our algorithm with probability at least $1-\eta$. Here, $\eta$ is the parameter of IMGPO. $\Xi_{max}$ is another parameter, but it is only used to limit the possibly long computation of step (iii) (in the worst case, step (iii) computes UCBs $3^{\Xi_{max}}$ times although it would rarely happen).

The pseudocode is shown in Algorithm 1. Lines 8 to 23 correspond to steps (i)-(iii). These lines compute the index $i_{h}^{*}$ of the candidate of the rectangle that may contain a global optimizer for each depth $h$. For each depth $h$, non-null index $i_{h}^{*}$ at Line 24 indicates the remaining candidate of a rectangle that we want to divide. Lines 24 to 33 correspond to steps (iv)-(v) where the remaining candidates of the rectangles for all $h$ are divided. To provide a simple executable division scheme (line 29), we assume $\Omega$ to be a hyperrectangle (see the last paragraph of section 4 for a general case).

Lines 8 to 17 correspond to steps (i)-(ii). Specifically, line 10 implements step (i) where a single candidate is selected for each depth, and lines 11 to 12 conduct step (ii) where some candidates are screened out. Lines 13 to 17 resolve the the temporary dummy values computed by GP. Lines 18 to 23 correspond to step (iii) where the candidates are further screened out. At line 21, $\mathcal{T}^{\prime}_{h+\xi}(c_{h,i^{*}_{h}})$ indicates the set of all center points of a fully expanded tree until depth $h+\xi$ within the region covered by the hyperrectangle centered at $c_{h,i^{*}_{h}}$. In other words, $\mathcal{T}^{\prime}_{h+\xi}(c_{h,i^{*}_{h}})$ contains the nodes of the fully expanded tree rooted at $c_{h,i^{*}_{h}}$ with depth $\xi$ and can be computed by dividing the current rectangle at $c_{h,i^{*}_{h}}$ and recursively divide all the resulting new rectangles until depth $\xi$ (i.e., depth $\xi$ from $c_{h,i^{*}_{h}}$, which is depth $h+\xi$ in the whole tree).

4 Relationship to Previous Algorithms

The idea of considering a set of infinitely many bounds was first proposed by Jones et al. . Their DIRECT algorithm has been successfully applied to real-world problems , but it only maintains the consistency property (i.e., convergence in the limit) from a theoretical viewpoint. DIRECT takes an input parameter $\epsilon$ to balance the global and local search efforts. This idea was generalized to the case of an unknown semi-metric and strengthened with a theoretical support (finite regret bound) by Munos in the SOO algorithm. By limiting the depth of the search tree with a parameter $h_{max}$, the SOO algorithm achieves a finite regret bound that depends on the near-optimality dimension.

Analysis

In this section, we prove an exponential convergence rate of IMGPO and theoretically discuss the reason why the novel idea underling IMGPO is beneficial. The proofs are provided in the supplementary material. To examine the effect of considering infinitely many possible candidates of the bounds, we introduce the following term.

(Infinite-metric exploration loss). The infinite-metric exploration loss $\rho_{t}$ is the number of intervals to be divided during iteration $t$.

In Theorem 1, we show that the exponential convergence rate $O\left(\lambda^{N+N_{gp}}\right)$ with $\lambda<1$ is achieved. We define $\Xi_{n}\leq\Xi_{max}$ to be the largest $\xi$ used so far with $n$ total node expansions. For simplicity, we assume that $\Omega$ is a square, which we satisfied in our experiments by scaling original $\Omega$.

Assume Assumptions 1 and 2. Let $\beta=\sup_{x,x^{\prime}\in\Omega}\frac{1}{2}\|x-x^{\prime}\|_{\infty}$. Let $\lambda=3^{-\frac{\alpha}{2CD\bar{\rho}_{t}}}<1$. Then, with probability at least $1-\eta$, the regret of IMGPO is bounded as

We note that $\lambda$ can get close to one as input dimension $D$ increases, which suggests that there is a remaining challenge in scalability for higher dimensionality. One strategy for addressing this problem would be to leverage additional assumptions such as those in .

(The effect of the tightness of UCB by GP) If UCB computed by GP is “useful” such that $N/\bar{\rho}_{t}=\Omega(N)$, then our regret bound becomes $O\left(\exp\left(-\frac{N+N_{gp}}{2CD}\alpha\ln 3\right)\right)$. If the bound due to UCB by GP is too loose (and thus useless), $\bar{\rho}_{t}$ can increase up to $O(N/t)$ (due to $\bar{\rho}_{t}\leq\sum^{t}_{i=1}{i}/{t}\leq O(N/t)$), resulting in the regret bound of $O\left(\exp\left(-\frac{t(1+N_{gp}/N)}{2CD}\alpha\ln 3\right)\right)$, which can be bounded by $O\left(\exp\left(-\frac{N+N_{gp}}{2CD}\max(\frac{1}{\sqrt{N}},\frac{t}{N})\alpha\ln 3\right)\right)$This can be done by limiting the depth of search tree as $depth(T)=O(\sqrt{N})$. Our proof works with this additional mechanism, but results in the regret bound with $N$ being replaced by $\sqrt{N}$. Thus, if we assume to have at least “not useless” UCBs such that $N/\bar{\rho}_{t}=\Omega(\sqrt{N})$, this additional mechanism can be disadvantageous. Accordingly, we do not adopt it in our experiments.. This is still better than the known results.

One may improve the algorithm with different division procedures than one presented in Algorithm 1 as discussed in the supplementary material.

Experiments

In this section, we compare the IMGPO algorithm with the SOO, BaMSOO, GP-PI and GP-EI algorithms . In previous work, BaMSOO and GP-UCB were tested with a pair of a handpicked good kernel and hyperparameters for each function . In our experiments, we assume that the knowledge of good kernel and hyperparameters is unavailable, which is usually the case in practice. Thus, for IMGPO, BaMSOO, GP-PI and GP-EI, we simply used one of the most popular kernels, the isotropic Matern kernel with $\nu=5/2$. This is given by $\kappa({x,x^{\prime}})=g(\sqrt{5||{x-x^{\prime}}||^{2}/l})$, where $g(z)={\sigma^{2}}(1+z+z^{2}/3)\exp(-z)$. Then, we blindly initialized the hyperparameters to $\sigma=1$ and $l=0.25$ for all the experiments; these values were updated with an empirical Bayesian method after each iteration. To compute the UCB by GP, we used $\eta=0.05$ for IMGPO and BaMSOO. For IMGPO, $\Xi_{max}$ was fixed to be $2^{2}$ (the effect of selecting different values is discussed later). For BaMSOO and SOO, the parameter $h_{max}$ was set to $\sqrt{n}$, according to Corollary 4.3 in . For GP-PI and GP-EI, we used the SOO algorithm and a local optimization method using gradients to solve the auxiliary optimization. For SOO, BaMSOO and IMGPO, we used the corresponding deterministic division procedure (given $\Omega$, the initial point is fixed and no randomness exists). For GP-PI and GP-EI, we randomly initialized the first evaluation point and report the mean and one standard deviation for 50 runs.

The experimental results for eight different objective functions are shown in Figure 2. The vertical axis is log${}_{10}(f(x^{*})-f(x^{+}))$, where $f(x^{*})$ is the global optima and $f(x^{+})$ is the best value found by the algorithm. Hence, the lower the plotted value on the vertical axis, the better the algorithm’s performance. The last five functions are standard benchmarks for global optimization . The first two were used in to test SOO, and can be written as $f_{sin1}(x)=(\sin(13x)\sin+1)/2$ for Sin1 and $f_{sin2}(x)=f_{sin1}(x_{1})f_{sin1}(x_{2})$ for Sin2. The form of the third function is given in Equation (16) and Figure 2 in . The last function is Sin2 embedded in 1000 dimension in the same manner described in Section 4.1 in , which is used here to illustrate a possibility of using IMGPO as a main subroutine to scale up to higher dimensions with additional assumptions. For this function, we used REMBO with IMGPO and BaMSOO as its Bayesian optimization subroutine. All of these functions are multimodal, except for Rosenbrock2, with dimensionality from 1 to 1000.

As we can see from Figure 2, IMGPO outperformed the other algorithms in general. SOO produced the competitive results for Rosenbrock2 because our GP prior was misleading (i.e., it did not model the objective function well and thus the property $f(x)\leq{\cal U}(x|{\cal D})$ did not hold many times). As can be seen in Table 1, IMGPO is much faster than traditional GP optimization methods although it is slower than SOO. For Sin 1, Sin2, Branin and Hartmann3, increasing $\Xi_{max}$ does not affect IMGPO because $\Xi_{n}$ did not reach $\Xi_{max}=2^{2}$ (Figure 2). For the rest of the test functions, we would be able to improve the performance of IMGPO by increasing $\Xi_{max}$ at the cost of extra CPU time.

Conclusion

We have presented the first GP-based optimization method with an exponential convergence rate $O\left(\lambda^{N+N_{gp}}\right)$ ($\lambda<1$) without the need of auxiliary optimization and the $\delta$-cover sampling. Perhaps more importantly in the viewpoint of a broader global optimization community, we have provided a practically oriented analysis framework, enabling us to see why not relying on a particular bound is advantageous, and how a non-tight bound can still be useful (in Remarks 1, 2 and 3). Following the advent of the DIRECT algorithm, the literature diverged along two paths, one with a particular bound and one without. GP-UCB can be categorized into the former. Our approach illustrates the benefits of combining these two paths.

As stated in Section 3.1, our solution idea was to use a bound-based method but rely less on the estimated bound by considering all the possible bounds. It would be interesting to see if a similar principle can be applicable to other types of bound-based methods such as planning algorithms (e.g., A* search and the UCT or FSSS algorithm ) and learning algorithms (e.g., PAC-MDP algorithms ).

The authors would like to thank Dr. Remi Munos for his thoughtful comments and suggestions. We gratefully acknowledge support from NSF grant 1420927, from ONR grant N00014-14-1-0486, and from ARO grant W911NF1410433. Kenji Kawaguchi was supported in part by the Funai Overseas Scholarship. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of our sponsors.

Bayesian Optimization with Exponential Convergence: Supplementary Material

In this supplementary material, we provide the proofs of the theoretical results. Along the way, we also prove regret bounds for a general class of algorithms, the result of which may be used to design a new algorithm.

We first provide a known property of the upper confidence bound of GP.

(Bound Estimated by GP) According to the belief encoded in the GP prior/posteriorThus, the probability in this analysis should be seen as that of the subjective view. If we assume that $f$ is indeed a sample from the GP, we have the same result with the objective view of probability. , for any $x$, $f(x)\leq{\cal U}(x|{\cal D})$ holds during the execution of Algorithm 1 with probability at least $1-\eta$.

It follows the proof of lemma 5.1 of . From the property of the standard gaussian distribution, $\Pr(f(x)>{\cal U}(x|{\cal D}))<\frac{1}{2}e^{-\varsigma^{2}_{M}/2}$. Taking union bound on the entire execution of Algorithm 1, $\Pr(f(x)>{\cal U}(x|{\cal D})\;\forall M\geq 1)<\frac{1}{2}\sum_{M=1}^{\infty}e^{-\varsigma^{2}_{M}/2}$. Substituting $\varsigma_{M}=\sqrt{2\log(\pi^{2}M^{2}/12\eta)}$, we obtain the statement. ∎

Our algorithm has a concrete division procedure in line 27 of Algorithm 1. However, one may improve the algorithm with different division procedures. Accordingly, we first derive abstract version of regret bound for the IMGPO (Algorithm 1) under a family of division procedures that satisfy Assumptions 3 and 4. After that, we provide a proof for the main results in the paper.

A With Family of Division Procedure

In this section, we modify the result obtained by . Let $x_{h,i}$ to be any point in the region covered by the $i$th hyperinterval at depth $h$, and $x_{h,i}^{*}$ be the global optimizer that may exist in the $i$th hyperinterval at depth $h$. The previous work provided the regret bound of the SOO algorithm with a family of division procedure that satisfies the following two assumptions.

Thus, in this section, hyperinterval is not restricted to hyperrectangle. We now revisit the definitions of several terms and variables used in . Let the $\epsilon$-optimal space $X_{\epsilon}$ be defined as $X_{\epsilon}:=\{x\in\Omega:f(x)+\epsilon\geq f({x^{*}})\}$. That is, the $\epsilon$-optimal space is the set of input vectors whose function value is at least $\epsilon$-close to the global optima. To bound the number of hyperintervals relevant to this $\epsilon$-optimal space, we define a near-optimality dimension as follows.

Finally, we define the set of $\delta$-optimal hyperintervals $I_{\delta(h)}$ as ${I_{\delta(h)}}:=\{c\in\Omega:f(c)+\delta(h)\geq f({x^{*}}),c\text{ is the center point of the interval},\omega_{h,i},\text{for some }(h,i)\}$. The $\delta$-optimal hyperinterval ${I_{\delta(h)}}$ is used to relate the hyperintervals to the $\epsilon$-optimal space. Indeed, the $\delta$-optimal hyperinterval ${I_{\delta(h)}}$ is almost identical to the $\delta(h)$-optimal space $X_{\delta(h)}$, except that ${I_{\delta(h)}}$ is focused on the center points whereas ${X_{\delta(h)}}$ considers the whole input vector space. In the following, we use $|{I_{\delta(h)}}|$ to denote the number of ${I_{\delta(h)}}$ and derive its upper bound.

(Lemma 3.1 in ) Let $d$ be the near-optimality dimension and $C$ denote the corresponding constant in Definition 1. Then, the number of $\delta$-optimal hyperintervals is bounded by $|{I_{\delta(h)}}|\leq C\delta{(h)^{-d}}$.

We are now ready to present the main result in this section. In the following, we use the term optimal hyperinterval to indicate a hyperinterval that contains a global optimizer $x^{*}$. We say a hyperinterval is dominated by other intervals when it is rejected or not selected in step (i)-(iii). In Lemma 3, we bound the maximum size of the optimal hyperinterval. From Assumption 1, this can be translated to the regret bound, as we shall see in Theorem 2.

Let $\Xi_{n}\leq\min(\Xi,\Xi_{max})$ be the largest $\xi$ used so far with $n$ total node expansions. Let $h_{n}^{*}$ be the depth of the deepest expanded node that contains a global optimizer ${x^{*}}$ after $n$ total node expansions (i.e., $h_{n}^{*}\leq n$ determines the size of the optimal hyperinterval). Then, with probability at least $1-\eta$, $h_{n}^{*}$ is bounded below by some $h^{\prime}$ that satisfies

Let $T_{h}$ denote the time at which the optimal hyperinterval is further divided. We prove the statement by showing that the time difference $T_{h+1}-T_{h}$ is bounded by the number of $\delta$-optimal hyperintervals. To do so, we first note that there are three types of hyperinterval that can dominate an optimal hyperinterval $c_{h+1,*}$ during the time $[T_{h},T_{h+1}-1]$, all of which belong to $\delta$-optimal hyperintervals $I_{\delta}$. The first type has the same size (i.e., same depth $h$), $c_{h+1,i}$. In this case,

where the first inequality is due to line 10 (step (i)) and the second follows Assumptions 1 and 2. Thus, it must be $c_{h+1,i}\in I_{h+1}$. The second case is where the optimal hyperinterval may be dominated by a hyperinterval of larger size (depth $l<h+1$), $c_{l,i}$. In this case, similarly,

where the first inequality is due to lines 11 to 12 (step (ii)) and thus $c_{l,i}\in I_{l}$. In the final scenario, the optimal hyperinterval is dominated by a hyperinterval of smaller size (depth $h+1+\xi$), $c_{h+1+\xi,i}$. In this case,

with probability at least $1-\eta$ where $z(\cdot,\cdot)$ is defined in line 21 of Algorithm 1. The first inequality is due to lines 19 to 23 (step (iii)) and the second inequality follows Lemma 1 and Assumptions 1 and 3. Hence, we can see that $c_{h+1+\xi,i}\in I_{h+1+\xi}$.

For all of the above arguments, the temporarily assigned ${\cal U}$ under GP has no effect. This is because the algorithm still covers the above three types of $\delta$-optimal hyperintervals $I_{\delta}$, as ${\cal U}\geq f$ with probability at least $1-\eta$ (Lemma 1). However, these are only expanded based on $f$ because of the temporary nature of ${\cal U}$. Putting these results together,

Since if one of the $I_{\delta}$ is divided during $[T_{h},T_{h+1}-1]$, it cannot be divided again during another time period,

where on the right-hand side, we could combine the summation $\sum^{h_{n}^{*}}_{h=0}$ and $\sum_{\tau=1}^{\sum_{l=1}^{h+1+\Xi_{n}}|I_{\delta(l)}|}$ into the one, because each $h$ in the summation refers to the same $\delta$-optimal interval $I_{\delta(l)}$ with $l\leq h_{n}^{*}+1+\Xi_{n}$, and should not be double-counted. As $\sum_{h=0}^{h_{n}^{*}}T_{h+1}-T_{h}=T_{h_{n}^{*}+1}-T_{0}$, $T_{0}=1$ and $|I_{\delta(0)}|=1$,

As $T_{{h_{n}^{*}}+1}>n$ by definition, for any $h^{\prime}$ such that $\sum_{\tau=1}^{\sum_{l=0}^{h^{\prime}+\Xi_{n}}|I_{l}|}\rho_{\tau}\leq n<\sum_{\tau=1}^{\sum_{l=0}^{h_{n}^{*}+1+\Xi_{n}}|I_{l}|}\rho_{\tau}$, we have $h_{n}^{*}>h^{\prime}$. ∎

With Lemmas 2 and 3, we are ready to present a finite regret bound with the family of division procedures.

Assume Assumptions 1, 3, and 4. Let $h(n)$ be the smallest integer $h$ such that

Then, with probability at least $1-\eta$, the regret of the IMGPO with any general division procedure is bounded as

Let $c(n)$ and $c_{h_{n}^{*},*}$ be the center point expanded at the $n$th expansion and the optimal hyperinterval containing a global optimizer $x^{*}$, respectively. Then, from Assumptions 1, 3, and 4, $f(c(n))\geq f(c_{h_{n}^{*},*})\geq f^{*}-\delta(h_{n}^{*})$, where $f^{*}$ is the global optima. Hence, the regret bound is $r_{h}\leq\delta(h_{n}^{*})$. To find a lower bound for the quantity $h_{n}^{*}$, we first relate $h(n)$ to Lemma 3 by

where the first inequality comes from the definition of $h(n)$, and the second follows from Lemma 2. Then, from Lemma 3, we have $h_{n}^{*}\geq h(n)-1$. Therefore, $r_{n}\leq\delta(h_{n}^{*})\leq\delta(h(n)-1)$. ∎

(Decreasing diameter revisit) The decreasing diameter defined in Assumption 3 can be written as $\delta(h)=c_{1}\gamma^{h/D}$ for some $c_{1}>0$ and $\gamma<1$ with a division procedure that requires $c_{2}$ function evaluations per node expansion.

Assume Assumptions 1, 3, 4, and 5. Then, if $d=0$, with probability at least $1-\eta$,

If $d>0$, with probability at least $1-\eta$,

For the case $d=0$, we have $n\leq\sum_{\tau=1}^{C\sum_{l=0}^{h(n)+\Xi_{n}}\delta(l)^{-d}}\rho_{\tau}\leq\sum_{\tau=1}^{C(h(n)+\Xi_{n}+1)}\bar{\rho}_{t}$, where the first inequality follows from the definition of $h(n)$, and the second comes from the definition of $\bar{\rho}_{t}$ and the assumption $d=0$. The second inequality holds for $\bar{\rho}_{t}$ that only considers $\rho_{\tau}$ with $\tau\leq t$. This is computable, because $\tau\leq t$ by construction. Indeed, the condition of Lemma 3 implies $t\geq\sum_{l=0}^{h^{\prime}+\Xi_{n}}|I_{l}|$. Therefore, the two inequalities hold, and we can deduce that $h(n)\geq\frac{n}{C\bar{\rho}_{t}}-\Xi_{n}-1$ by algebraic manipulation. By Assumption 5, $n=(N+N_{gp})/c_{2}$. With this, substituting the lower bound of $h(n)$ into the statement of Theorem 2 with Assumption 5,

and hence $c\gamma^{\frac{h(n)+\Xi_{n}}{D}}\leq\left(\frac{n(1-\gamma^{d/D})}{C\bar{\rho}_{t}}\right)^{-1/d}$ by algebraic manipulation. Substituting this into the result of Theorem 2, we arrive at the desired result. ∎

B With a Concrete Division Procedure

In this section, we prove the main result in the paper. In Theorem 1, we show that the exponential convergence rate bound $O\left(\lambda^{N+N_{gp}}\right)$ with $\lambda<1$ is achieved without Assumptions 3, 4 and 5 and without the assumption that $d=0$.

Assume Assumptions 1 and 2. Let $\beta=\sup_{x,x^{\prime}\in\Omega}\frac{1}{2}\|x-x^{\prime}\|_{\infty}$. Let $\lambda=3^{-\frac{\alpha}{2C\bar{\rho}_{t}D}}<1$. Then, without Assumptions 3, 4 and 5 and without the assumption on $d$, with probability at least $1-\eta$, the regret of IMGPO with the division procedure in Algorithm 1 is bounded as

To prove the statement, we show that Assumptions 3, 4, and 5 can all be satisfied while maintaining $d=0$.

From Assumption 2 (i), and based on the division procedure that the algorithm uses,

This upper bound corresponds to the diagonal length of each hyperrectangle with respect to $p$-norm, where $3^{-\lfloor h/D\rfloor}\beta$ corresponds to the length of the longest side. We fix the form of $\delta$ as $\delta(h)=L3^{\alpha}D^{\alpha/p}3^{-h\alpha/D}\beta^{\alpha}\geq L(3^{-\lfloor h/D\rfloor}\beta D^{1/p})^{\alpha}$, which satisfies Assumption 3.

This form of $\delta(h)$ also satisfies Assumption 5 with $\gamma=3^{-\alpha}$ and $c_{1}=L3^{\alpha}D^{\alpha/p}\beta^{\alpha}$.

Finally, we show that $d=0$. The set of $\delta$-optimal hyperintervals $I_{\delta(h)}$ is contained by the $\delta(h)$-optimal space $X_{\delta(h)}$ as

Now that we have satisfied Assumptions 3, 4, and 5 with $d=0$, $\gamma=3^{-\alpha}$, and $c_{1}=L3^{\alpha}D^{\alpha/p}\beta^{\alpha}$, we follow the proof of Corollary 1 and deduce the desired statement. ∎

References