Neural-BO: A Black-box Optimization Algorithm using Deep Neural Networks
Bayesian Optimization (BO) is an effective approach for global optimization
of black-box functions when function evaluations are expensive. Most prior
works use Gaussian processes to model the black-box function, however, the use
of kernels in Gaussian processes leads to two problems: first, the kernel-based
methods scale poorly with the number of data points and second, kernel methods
are usually not effective on complex structured high dimensional data due to
curse of dimensionality. Therefore, we propose a novel black-box optimization
algorithm where the black-box function is modeled using a neural network. Our
algorithm does not need a Bayesian neural network to estimate predictive
uncertainty and is therefore computationally favorable. We analyze the
theoretical behavior of our algorithm in terms of regret bound using advances
in NTK theory showing its efficient convergence. We perform experiments with
both synthetic and real-world optimization tasks and show that our algorithm is
more sample efficient compared to existing methods.
Dat Phan-Trong, Hung Tran-The, Sunil Gupta
# Neural-BO: A Black-box Optimization Algorithm using Deep Neural Networks

###### Abstract

Bayesian Optimization (BO) is an effective approach for global optimization of black-box functions when function evaluations are expensive. Most prior works use Gaussian processes to model the black-box function, however, the use of kernels in Gaussian processes leads to two problems: first, the kernel-based methods scale poorly with the number of data points and second, kernel methods are usually not effective on complex structured high dimensional data due to curse of dimensionality. Therefore, we propose a novel black-box optimization algorithm where the black-box function is modeled using a neural network. Our algorithm does not need a _Bayesian_ neural network to estimate predictive uncertainty and is therefore computationally favorable. We analyze the theoretical behavior of our algorithm in terms of regret bound using advances in NTK theory showing its efficient convergence. We perform experiments with both synthetic and real-world optimization tasks and show that our algorithm is more sample efficient compared to existing methods.

keywords: Black-box Optimization, Neural Tangent Kernel +
Footnote †: journal: Neurocomputing

## 1 Introduction

Optimizing a black-box function is crucial in many real-world application domains. Some well-known examples include hyper-parameter optimization
of machine learning algorithms [1; 2], synthesis of short polymer fiber materials, alloy design, 3D bio-printing, molecule design, etc [3; 4]. Bayesian Optimization (BO) has emerged as an effective approach to optimize such expensive-to-evaluate black-box functions. Most BO algorithms work as follows. First, a probabilistic regression model e.g. a Gaussian Process (GP) is trained on available function observations. This model is then used to construct an acquisition function that trades off between two potentially conflicting objectives: exploration and exploitation. The acquisition function is then optimized to suggest a point where the black-box functions should be next evaluated. There are several choices for acquisition functions, e.g., improvement-based schemes (Probability of Improvement [5], Expected improvement [6]), Upper Confidence Bound scheme [7], entropy search (ES [8], PES [9]), Thompson Sampling [10] and Knowledge Gradient [11].

While GPs are a natural and usually the most prominent probabilistic model choice for BO, unfortunately they are not suitable for optimizing functions in complex structured high dimensional spaces such as functions involving optimization of images and text documents. This is because GPs use kernels to model the functions and the popular kernels e.g. Square Exponential kernel, Matern kernel cannot accurately capture the similarity between complex structured data. This limitation prevents BO with GP from being applied to computer vision and natural language processing applications. Further, GPs do not scale well computationally due to kernel inversion requiring cubic complexity.

Recently, Deep Neural Networks (DNNs) have emerged as state-of-the-art models for a wide range of tasks such as speech recognition, computer vision, natural language processing, and so on. The emergence of DNNs proves their power to scale and generalize to most kinds of data, including high dimensional and structured data as they can extract powerful features and scale linearly with the number of data points unlike GPs. Leveraging these advantages, some efforts have been made to utilize DNNs for optimization e.g., black-box function optimization in continuous search spaces [12; 13; 14] and contextual bandit optimizations in discrete search spaces [15; 16; 17; 18]. However, these works have several drawbacks: (a) some of these approaches are computationally expensive because of relying on Bayesian Neural Networks (BNNs), making them as unscalable as GPs. (b) other works in continuous search space setting do not provide any theoretical analysis of the respective optimization algorithms. Therefore, the problem of developing a DNN-based black-box optimization algorithm that is computationally effi
cient and comes with a mathematical guarantee of its convergence and sample efficiency remains open.

Addressing the above challenge, we propose a novel black-box optimization algorithm where the unknown function is modeled using an over-parameterized deep neural network. Our algorithm does not need a _Bayesian Neural Network_ to model predictive uncertainty and is therefore computationally favorable. We utilize the recent advances in neural network theory using Neural Tangent Kernels (NTK) to estimate a confidence interval around the neural network model and then use it for Thompson Sampling to draw a random sample of the function, which is maximized to recommend the next function evaluation point. We analyze the theoretical behavior of our algorithm in terms of its regret bound and show that our algorithm is convergent and more sample efficient compared to existing methods. In summary, our contributions are:

* A deep neural network based black-box optimization (Neural-BO) that can perform efficient optimization by accurately modeling the complex structured data (e.g. image, text, etc) and also being computationally favorable, i.e. scaling linearly with the number of data points.
* A theoretical analysis of our proposed Neural-BO algorithm showing that our algorithm is convergent with a sub-linear regret bound. To the best of our knowledge, our work is the first to achieve sub-linear bound in any neural network based bandit or black-box optimization algorithms.
* Experiments with both synthetic benchmark optimization functions and real-world optimization tasks showing that our algorithm outperforms current state-of-the-art black-box optimization methods.

## 2 Related Works

### Bayesian optimization

Bayesian optimization is a well-established approach for performing the global optimization of noisy, expensive black-box functions [6]. Traditionally, Bayesian optimization relies on the construction of a probabilistic surrogate model that places a prior on objective functions and updates the model using the existing function evaluations. To handle the promise and uncertainty of a new point, an acquisition function is used to capture the trade-off between
exploration and exploitation through the posterior distribution of the functions. By performing a proxy optimization over this acquisition function, the next evaluation point is determined. The Expected Improvement (EI) [6] and Upper Confidence Bound (UCB) [7] criteria are typically used due to their analytical tractability in both value and gradient. Beside, the information-theoretic acquisition functions e.g., (ES [8], PES [19], MES [9]) guide our evaluations to locations where we can maximise our learning about the unknown minimum rather than to locations where we expect to obtain higher function values. For a review on Bayesian optimization, see [3; 4; 20; 21].

### Bayesian Optimization using GPs

Because of their flexibility, well-calibrated uncertainty, and analytical properties, GPs have been a popular choice to model the distribution over functions for Bayesian optimization. Formally, a function \(f\) is assumed to be drawn from GP, i.e., \(f\sim\text{GP}(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime}))\), where \(\mu(\mathbf{x})\) is a mean function and \(k(\mathbf{x},\mathbf{x}^{\prime})\) is a covariance function. Common covariance functions include linear kernel, squared exponential (SE) kernel and Matern kernel. Thus, given a set of \(t\) observed points \(\mathcal{D}_{1:t}=\{\mathbf{x}_{i},y_{i}\}_{i=1}^{t}\), the posterior mean and variance at step \(t+1\) are computed as: \(\mu_{t+1}(\mathbf{x})=\mathbf{k}^{\top}\mathbf{K}^{-1}\mathbf{y}\) and variance \(\sigma_{t+1}^{2}(\mathbf{x})=k(\mathbf{x},\mathbf{x})-\mathbf{k}^{\top} \mathbf{K}^{-1}\mathbf{k}\), where \(\mathbf{K}=\left[k(\mathbf{x}_{i},\mathbf{x}_{j})\right]_{1\leq i,j\leq t}\), \(\mathbf{k}=\left[k(\mathbf{x},\mathbf{x}_{1})\ \ k(\mathbf{x},\mathbf{x}_{2})\ \ldots\ \ k( \mathbf{x},\mathbf{x}_{t})\right]\) and \(\mathbf{y}=[y_{1},y_{2},\cdots,y_{t}]\). Nonetheless, a major drawback of GP-based Bayesian optimization is that the computational complexity grows cubically in the number of observations, as it necessitates the inversion of the covariance matrix. Further, the covariance functions are usually too simplistic to model the function effectively when the inputs are high dimensional and complex structured data such as images and text.

### Black-box Optimization using models other than GPs

There are relatively few works in black-box optimization that have used models other than GPs. For example, random forest (RF) and neural networks (NNs) have been used to estimate black-box functions. The first work that utilized RF as surrogate model was proposed in SMAC [22] and later refined in [23]. Although random forest performs well by making good predictions in the neighbourhood of training data, it faces trouble with extrapolation [4]. Over the past decade, neural networks have emerged as a potential candidate to play the role of surrogate model in black-box optimization. The first work on black-box optimization using neural networks [12] utilizes
a deep neural network as a mapping to transform high-dimensional inputs to lower-dimensional feature vectors before fitting a Bayesian linear regressor surrogate model. Another attempt has been made in [13] where the objective function has been modeled using a Bayesian neural network, which has high computational complexity due to the need to maintain an ensemble of models to get the posterior distribution over the parameters. Although these early works show the promise of neural networks for black-box optimization tasks through empirical results, they are still computationally unscalable (due to using Bayesian models) and no theoretical guarantees have been provided to ensure the convergence of their algorithms.

Very recently, [14] also explores the use of deep neural networks as an alternative to GPs to model black-box functions. However, this work only provides theoretical guarantees for their algorithm under _infinite width_ of DNN setting, which is unrealistic in practice, while we provide theoretical guarantees for over-parameterized DNN, which is more general than their setting. Furthermore, [14] simply used a noise on the target does not correctly account for epistemic uncertainty and thus this method performs exploration randomly like \(\epsilon\)-greedy, which is known to be sub-optimal. Moreover, their experimental results are limited to only synthetic functions, which is not enough to demonstrate the effectiveness of an algorithm in real world. In contrast, we demonstrate the effectiveness of their algorithm on various real applications with high dimensional structured and complex inputs.

### Contextual Bandits with Neural Networks

There is a research line of contextual bandits which also uses deep neural networks to model reward functions. This problem is defined as follows: at each round \(t\), the agent is presented with a set of \(K\) actions associated with \(K\) feature vectors called "contexts". After choosing an arm, the agent receives a scalar reward generated from some unknown function with respect to the corresponding chosen feature vector. The goal of the agent is to maximize the expected cumulative reward over \(T\) rounds. For theoretical guarantees, most of works in this line (e.g., [15; 17; 16; 18]) utilize a _over-parametrized_ DNN to predict rewards. While NeuralUCB [15] and NeuralTS [17] use UCB and Thompson Sampling acquisition functions to control the trade-off between exploration and exploitation by utilizing gradient of the network in algorithms, [16] performs LinearUCB on top of the representation of the last layer of the neural network, motivated by [12].
Similar to the above contextual bandit algorithms, we attempt to develop an effective algorithm with theoretical analysis for black-box optimization using over-parametrized neural networks with _finite_ width. The key challenge we need to address when extending the contextual bandit algorithms to black-box optimization is that of dealing with a continuous action space instead of a discrete space. Having continuous action space requires to bound the function modeling error at arbitrary locations in the action space. More specifically, from neural networks' viewpoint, this requirement is equivalent to characterizing predictive output on unseen data, which requires using neural network generalization theory different from the problem setting of contextual bandits.

## 3 Preliminaries

### Problem Setting

We consider a global optimization problem to maximize \(f(\mathbf{x})\) subject to \(\mathbf{x}\in\mathcal{D}\subset\mathbb{R}^{d}\), where \(\mathbf{x}\) is the input with \(d\) dimensions. The function \(f\) is a noisy black-box function that can only be evaluated via noisy evaluations of \(f\) in the form \(y=f(x)+\epsilon\), where \(\epsilon\) is a sub-Gaussian noise which we will discuss more clearly in Assumption 2 in our regret analysis section. In our work, we consider input space \(\mathcal{D}\) of the form: \(a\leq\left\|\mathbf{x}\right\|_{2}\leq b\), where \(a,b\) are absolute constants and \(a\leq b\). We note that this assumption is mild compared to current works in neural contextual bandits [15; 17; 16; 18] that required that \(\left\|\mathbf{x}\right\|_{2}=1\).

To imply the smoothness of function \(f\), we assume the objective function \(f\) is from an RKHS corresponding to a positive definite Neural Tangent Kernel \(k_{\mathrm{NTK}}\). In particular, \(\left\|f\right\|_{\mathcal{H}_{k_{\mathrm{NTK}}}}\leq B\), for a finite \(B>0\). These assumptions are regular and have been typically used in many previous works [10; 7; 24].

### RKHS with Neural Tangent Kernel

Given a training set \(\mathcal{X}_{t}=\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{t}\subset\mathcal{D}\times \mathbb{R}\), that is a set of \(t\) arbitrary points picked from \(\mathcal{D}\) with its corresponding observations, and consider training a neural network \(h(\mathbf{x};\mathbf{\theta})\) with gradient descent using an infinitesimally small learning rate. The Neural Tangent Kernel (NTK) \(k_{\mathrm{NTK}}(\mathbf{x},\mathbf{x}^{\prime})\) is defined as: \(k_{\mathrm{NTK}}(\mathbf{x},\mathbf{x}^{\prime})=\langle\mathbf{g}(\mathbf{x };\mathbf{\theta}),\mathbf{g}(\mathbf{x}^{\prime};\mathbf{\theta})\rangle/m\), where \(\mathbf{g}(\mathbf{x};\mathbf{\theta})\) is the gradient of \(h(\mathbf{x};\mathbf{\theta})\) as defined. When the network width \(m\) increases, \(\langle\mathbf{g}(\mathbf{x};\mathbf{\theta}),\mathbf{g}(\mathbf{x}^{\prime},\bm {\theta})\rangle\) gradually converges to \(\langle\mathbf{g}(\mathbf{x},\mathbf{\theta}_{0}),\mathbf{g}(\mathbf{x}^{\prime}, \mathbf{\theta}_{0})\rangle\) and the corresponding NTK matrix is defined in Definition 1.
**Definition 1** ([25]): _Consider the same set \(\mathcal{X}_{t}\) as in Section 3.2. Define_

\[\widetilde{\mathbf{H}}_{i,j}^{(1)}=\mathbf{\Sigma}_{i,j}^{(1)}=\left\langle \mathbf{z}_{i},\mathbf{z}_{j}\right\rangle\,\mathbf{A}_{i,j}^{(l)}=\begin{pmatrix} \mathbf{\Sigma}_{i,i}^{(l)}&\mathbf{\Sigma}_{i,j}^{(l)}\\ \mathbf{\Sigma}_{i,j}^{(l)}&\mathbf{\Sigma}_{j,j}^{(l)}\end{pmatrix}\]

\[\mathbf{\Sigma}_{i,j}^{(l+1)}=2\mathbb{E}_{(u,v)\sim\mathbf{N}(\mathbf{0}, \mathbf{A}_{i,j}^{(l)})}[\phi(u)\phi(v)]\]

\[\widetilde{\mathbf{H}}_{i,j}^{(l+1)}=2\widetilde{\mathbf{H}}_{i,j}^{(l)} \mathbb{E}_{(u,v)\sim\mathbf{N}(\mathbf{0},\mathbf{A}_{i,j}^{(l)}}[\phi^{\prime }(u)\phi^{\prime}(v)]+\mathbf{\Sigma}_{i,j}^{(l+1)},\]

where \(\phi(\cdot)\) is the activation function of the neural network \(h(\mathbf{x};\boldsymbol{\theta})\). Then \(\mathbf{H}_{t}=\frac{\widehat{\mathbf{H}}^{(L)}+\mathbf{\Sigma}^{(L)}}{2}\) is called the NTK matrix on the set \(\mathbf{X}\).

We assume \(f\) to be a member of the Reproducing Kernel Hilbert Space (RKHS) of real-valued functions on \(\mathcal{D}\) with Neural Tangent Kernel (NTK, see Section 3.2), and a bounded norm, \(\|f\|_{\mathcal{H}_{k_{\text{NTK}}}}\leq B\). This RKHS, denoted by \(\mathcal{H}_{k_{\text{NTK}}}(\mathcal{D})\), is completely specified by its kernel function \(k_{\text{NTK}}(\cdot,\cdot)\) and vice-versa, with an inner product \(\left\langle\cdot,\cdot\right\rangle\) obeying the reproducing property: \(f(\mathbf{x})=\left\langle f,k_{\text{NTK}}(\cdot,\mathbf{x})\right\rangle\) for all \(f\in\mathcal{H}_{k_{\text{NTK}}}(\mathcal{D})\). Equivalently, the RKHS associated with \(k_{\text{NTK}}\) is then given by:

\[\mathcal{H}_{k_{\text{NTK}}}=\{f:f=\sum_{i}\alpha_{i}k_{\text{NTK}}(\mathbf{x} _{i},\cdot),\alpha_{i}\in\mathbb{R}\}\]

The induced RKHS norm \(\|f\|_{\mathcal{H}_{k_{\text{NTK}}}}=\sqrt{\left\langle f,f\right\rangle_{ \mathcal{H}_{k_{\text{NTK}}}}}\) measures the smoothness of \(f\), with respect to the kernel function \(k_{\text{NTK}}\), and satisfies: \(f\in\mathcal{H}_{k_{\text{NTK}}}(\mathcal{D})\) if and only if \(\|f\|_{k_{\mathcal{H}_{\text{NTK}}}}<\infty\).

### Maximum Information Gain

Assume after \(t\) steps, the model receives an input sequence \(\mathcal{X}_{t}=(\mathbf{x}_{1},\mathbf{x}_{2},\ldots\mathbf{x}_{t})\) and observes noisy rewards \(\mathbf{y}_{t}=(y_{1},y_{2},\ldots,y_{t})\). The _information gain_ at step \(t\), quantifies the reduction in uncertainty about \(f\) after observing \(\mathbf{y}_{t}\), defined as the mutual information between \(\mathbf{y}_{t}\) and \(f\):

\[I(\mathbf{y}_{t};f):=H(\mathbf{y}_{t})-H(\mathbf{y}_{t}|f),\]

where \(H\) denotes the entropy. To obtain closed from the expression of information gain, one needs to introduce a GP model where \(f\) is assumed to be a zero mean GP indexed on \(\mathcal{X}\) with kernel \(k\). Let \(\mathbf{f}_{t}=(f(\mathbf{x}_{1}),f(\mathbf{x}_{1}),\ldots,f(\mathbf{x}_{t})])\)
be the corresponding true function values. From [26], the mutual information between two multivariate Gaussian distributions is:

\[I(\mathbf{y}_{t};f)=I(\mathbf{y}_{t};\mathbf{f}_{t})=\frac{1}{2}\log\det\bigl{(} \mathbf{I}+\lambda^{-1}\mathbf{H}_{t}\bigr{)},\]

where \(\lambda>0\) is a regularization parameter and \(\mathbf{H}_{t}\) is the kernel matrix. In our case, \(\mathbf{H}_{t}\) can be referred to as the NTK matrix associated with Neural Tangent Kernel.

Assuming GP model with NTK kernel to calculate closed-form expression of maximum information gain is popular in NTK-related works [18; 27].

To further achieve input-independent and kernel-specific bounds, we define the _maximum information gain_ as follows: The _maximum information gain_\(\gamma_{t}\), after \(t\) steps, defined as:

\[\gamma_{t}:=\max_{\mathcal{A}\subset\mathcal{X},|A|=t}I(\mathbf{y}_{A}; \mathbf{f}_{A})\]

As stated in Lemma 3 of [10], the maximum information gain \(\gamma_{t}\) can be lower bounded as:

\[\gamma_{t}\geq\frac{1}{2}\log\det\bigl{(}\mathbf{I}+\lambda^{-1}\mathbf{H}_{t }\bigr{)}\]

## 4 Proposed Neural-BO Algorithm

```
1:Set \(\mathbf{U}_{0}=\lambda\mathbf{I}\)
2:for\(t=1\) to \(T\)do
3:\(\sigma_{t}^{2}(\mathbf{x})=\lambda\mathbf{g}(\mathbf{x};\boldsymbol{\theta}_ {0})^{\top}\mathbf{U}_{t-1}^{-1}\mathbf{g}(\mathbf{x};\boldsymbol{\theta}_{0 })/m\)
4:\(\widetilde{f}_{t}(\mathbf{x})\sim\mathcal{N}(h(\mathbf{x};\boldsymbol{\theta}_ {t-1}),\nu_{t}^{2}\sigma_{t}^{2}(\mathbf{x}))\)
5: Choose \(\mathbf{x}_{t}=\operatorname*{argmax}_{\mathbf{x}\in\mathcal{D}}\widetilde{f}_ {t}(\mathbf{x})\) and receive observation \(y_{t}=f(\mathbf{x}_{t})+\epsilon_{t}\)
6: Update \(\boldsymbol{\theta}_{t}\leftarrow\operatorname*{TrainNN}(\{\mathbf{x}_{i}\}_{ i=1}^{t}\{y_{i}\}_{i=1}^{t},\boldsymbol{\theta}_{0})\)
7:\(\mathbf{U}_{t}\leftarrow\mathbf{U}_{t-1}+\frac{\mathbf{g}(\mathbf{x}_{t}; \boldsymbol{\theta}_{0})\mathbf{g}(\mathbf{x}_{t};\boldsymbol{\theta}_{0})^{ \top}}{m}\)
8:endfor
```

**Algorithm 1** Neural Black-box Optimization (Neural-BO)
In this section, we present our proposed neural-network based black-box algorithm: Neural-BO. Following the principle of BO, our algorithm consists of two main steps: (1) Build a model of the black-box objective function, and (2) Use the black-box model to choose the next function evaluation point at each iteration. For the first step, unlike traditional BO algorithms which use Gaussian process to model the objective function, our idea is to use a fully connected neural network \(h(\mathbf{x};\mathbf{\theta})\) to learn the function \(f\) as follows:

\[h(\mathbf{x};\mathbf{\theta})=\sqrt{m}\mathbf{W}_{L}\phi(\mathbf{W}_{L-1}\phi( \cdots\phi(\mathbf{W}_{1}\mathbf{x})),\]

where \(\phi(x)=\max(x,0)\) is the rectified linear unit (ReLU) activation function, \(\mathbf{W}_{1}\in\mathbb{R}^{m\times d},\mathbf{W}_{i}\in\mathbb{R}^{m\times m },2\leq i\leq L-1,\mathbf{W}_{L}\in\mathbb{R}^{1\times m}\), and \(\mathbf{\theta}=(vec(\mathbf{W}_{1}),\cdots,vec(\mathbf{W}_{L}))\in\mathbb{R}^{p}\) is the collection of parameters of the neural network, \(p=md+m^{2}(L-2)+m\). To keep the analysis convenient, we assume that the width \(m\) is the same for all hidden layers. We also denote the gradient of the neural network by \(\mathbf{g}(\mathbf{x};\mathbf{\theta})=\nabla_{\mathbf{\theta}}h(\mathbf{x};\mathbf{ \theta})\in\mathbb{R}^{p}\).

For the second step, we choose the next sample point \(\mathbf{x}_{t}\) by using a Thompson Sampling strategy. In particular, given each \(\mathbf{x}\), our algorithm maintains a Gaussian distribution for \(f(\mathbf{x})\). To choose the next point, it samples a random function \(\widetilde{f_{t}}\) from posterior distribution \(\mathcal{N}(h(\mathbf{x},\mathbf{\theta}_{t-1}),\nu_{t}^{2}\sigma_{t}^{2}(\mathbf{ x}))\), where \(\sigma_{t}^{2}(\mathbf{x})=\lambda\mathbf{g}(\mathbf{x};\mathbf{\theta}_{0}) \mathbf{U}_{t-1}^{-1}\mathbf{g}(\mathbf{x};\mathbf{\theta}_{0})/m\) and \(\nu_{t}=\sqrt{2}B+\frac{R}{\sqrt{\lambda}}(\sqrt{2\log(1/\alpha)}\). From here, \(\sigma_{t}(\mathbf{x})\) and \(\nu_{t}\) construct a confidence interval, where for all \(\mathbf{x}\in\mathcal{D}\), we have \(|f(\mathbf{x})-f(\mathbf{x},\mathbf{\theta})|\leq\nu_{t}\sigma_{t}(\mathbf{x})\) that holds with probability \(1-\alpha\). We also note the difference from [17] that uses \(\mathbf{g}(\mathbf{x};\mathbf{\theta}_{t})/\sqrt{m}\) as dynamic feature map, we here use a fixed feature map \(\mathbf{g}(\mathbf{x};\mathbf{\theta}_{0})/\sqrt{m}\) which can be viewed as a finite approximation of the feature map \(\phi(\mathbf{x})\) of \(k_{\text{NTK}}\). This allows us to eliminate \(\sqrt{m}\) away the regret bound of our algorithm. Our algorithm chooses \(\mathbf{x}_{t}\) as \(\mathbf{x}_{t}=\text{argmax}_{\mathbf{x}\in\mathcal{D}}\widetilde{f_{t}}( \mathbf{x})\).
Besides, to initialize the network \(h(\mathbf{x};\mathbf{\theta}_{0})\), we randomly generate each element of \(\mathbf{\theta}_{0}\) from an appropriate Gaussian distribution: for each \(1\leq l\leq L-1\), each entry of \(\mathbf{W}\) is generated independently from \(\mathcal{N}(0,2/m)\); while each entry of the last layer \(\mathbf{W}_{L}\) is set to zero to make \(h(\mathbf{x},\mathbf{\theta}_{0})=0\). Here we inherit so-called _He initialization_[28], which ensures the stability of the expected length of the output vector at each hidden layer.

The proposed Neural-BO is summarized using Algorithm 1 and Algorithm 2. In section 5, we will demonstrate that our algorithm can achieve a sub-linear regret bound, and works well in practice (See our section 7).

DiscussionA simple way to have a neural network based black-box algorithm is to extend the existing works of [15; 17] or [18] to our setting where the search space is continuous. However, this is difficult because these algorithms are designed for a finite set of actions. A simple adaptation can yield non-vanishing bounds on cumulative regret. For example, the regret bound in [18] depends on the size of finite action spaces \(|\mathcal{A}|\). If \(|\mathcal{A}|\to\infty\) then the regret bound goes to infinity. In [15; 17], the regret bounds are built on gram matrix \(\mathbf{H}\) which is defined through the set of actions. However, such a matrix cannot even be defined in our setting where the search space (action space) is infinite. We solve this challenge by using a discretization of the search space at each iteration as mentioned in Section 5. In addition, for our confidence bound calculation, using \(\mathbf{g}(\mathbf{x};\mathbf{\theta}_{0})/\sqrt{m}\) instead of \(\mathbf{g}(\mathbf{x};\mathbf{\theta}_{t})/\sqrt{m}\) as in [15; 17] allows us to reduce a factor \(\sqrt{m}\) from our regret bound.

Our algorithm is also different from [10] in traditional Bayesian optimization which uses a Gaussian process to model the function with posterior computed in closed form. In contrast, our algorithm computes posterior by gradient descent.

## 5 Regret Analysis

In this section, we provide a regret bound for the proposed NeuralBO algorithm. Our regret analysis is built upon the recent advances in NTK theory [29; 30] and proof techniques of GP-TS [10]. Unlike most previous neural-based contextual bandits works [17; 15] that restrict necessarily the search space to a set \(\mathbb{S}^{d-1}=\{\mathbf{x}\in\mathbb{R}^{d}\colon\left\|\mathbf{x}\right\| _{2}=1\}\), we here derive our results for a flexible condition, where the norm of inputs are bounded: \(a\leq\left\|\mathbf{x}\right\|_{2}\leq b,\) where \(0<a\leq b\).

For theoretical guarantees, we need to use the following assumptions on the observations \(\{\mathbf{x}_{i}\}_{i=1}^{T}\).
**Assumption 1**.: There exists \(\lambda_{0}>0\), \(\mathbf{H}_{T}\geq\lambda_{0}\mathbf{I}_{T}\).

This is a common assumption in the literature [15, 16, 17], and is satisfied as long as we sufficiently explore the input space such that no two inputs \(\mathbf{x}\) and \(\mathbf{x}^{\prime}\) are identical.

**Assumption 2**.: We assume the noises \(\{\epsilon_{i}\}_{i=1}^{T}\) where \(\epsilon_{i}=y_{i}-f(\mathbf{x}_{i})\), are i.i.d and sub-Gaussian with parameter \(R>0\).

This assumption is mild and similar to the work of [31] where the assumption is used to prove an optimal order of the regret bound for Gaussian process bandits. We use it here to provide a regret bound with the sublinear order for our proposed algorithm.

Now we are ready to bound the regret of our proposed NeuralBO algorithm. To measure the regret of the algorithm, we use the cumulative regret which is defined as follows: \(R_{T}=\sum_{t=1}^{T}r_{t}\) after \(T\) iterations, where \(\mathbf{x}^{*}=\operatorname*{argmax}_{\mathbf{x}\in\mathcal{D}}f(\mathbf{x})\) is the maximum point of unknown function \(f\) and \(r_{t}=f(\mathbf{x}^{*})-f(\mathbf{x}_{t})\) is instantaneous regret incurred at time \(t\). We present our main theoretical result.

**Theorem 1**.: _Let \(\alpha\in(0,1)\). Assume that the true underlying \(f\) lies in the RKHS \(\mathcal{H}_{k_{\mathrm{NTK}}}(\mathcal{D})\) corresponding to the NTK kernel \(k_{\mathrm{NTK}}(\mathbf{x},\mathbf{x}^{\prime})\) with RKHS norm bounded by \(B\). Set the parameters in Algorithm 1 as \(\lambda=1+\frac{1}{T}\), \(\nu_{t}=\sqrt{2}B+\frac{R}{\sqrt{\lambda}}(\sqrt{2\log(1/\alpha)}\), where \(R\) is the noise sub-Gaussianity parameter. Set \(\eta=(m\lambda+mLT)^{-1}\) and \(J=(1+LT/\lambda)\left(1+\log(T^{3}L\lambda^{-1}\log\left(1/\alpha\right))\right)\). If the network width m satisfies:_

\[m\geq\operatorname{poly}\left(\lambda,T,\log(1/\alpha),\lambda_{0}^{-1}\right),\]

_then with probability at least \(1-\alpha\), the regret of Algorithm 1 is bounded as_

\[R_{T} \leq C_{1}(1+c_{T})\nu_{T}\sqrt{L}\sqrt{\frac{\lambda BT}{\log(B+ 1)}(2\gamma_{T}+1)}\] \[+(4+C_{2}(1+c_{T})\nu_{T}L)\sqrt{2\log(1/\alpha)T}+\frac{(2B+1) \pi^{2}}{6}+\frac{b(a^{2}+b^{2})}{a^{3}}\]

_where \(C_{1},C_{2}\) are absolute constants, \(a\) and \(b\) are lower and upper norm bounds of \(\mathbf{x}\in\mathcal{D}\) and \(c_{T}=\sqrt{4\log T+2\log\ |\mathcal{D}_{t}|}\)._
**Remark 1**.: There exists a component \(\frac{b(a^{2}+b^{2})}{a^{3}}\) in the regret bound \(R_{T}\). This component comes from the condition of the search space. If we remove the influence of constants \(a,b,B\) and \(R\), Theorem 1 implies that our regret bound follows an order \(\widetilde{O}(\sqrt{T\gamma_{T}})\). This regret bound has the same order as the one of the MVR algorithm in [31], but is significantly better than traditional BO algorithms (e.g., [7; 10]) and existing neural contextual bandits (e.g., [15; 17]) that have the order \(\widetilde{O}(\gamma_{T}\sqrt{T})\). We note that the regret bounds in [15; 17] are expressed through the effective dimension \(\widetilde{d}=\frac{\log\det(\mathbf{I}+\mathbf{H}/\lambda)}{\log(1+TK/\lambda)}\). However, \(\gamma_{T}\) and the effective dimension are of the same order up to a ratio of \(1/(2\log(1+TK/\lambda))\). The improvement of factor \(\sqrt{\gamma_{T}}\) is significant because for NTK kernel as shown by [18], \(\gamma_{T}\) follows order \(\mathcal{O}(T^{1-1/d})\), where \(d\) is the input dimension. For \(d\geq 2\), existing bounds become vacuous and are no longer sub-linear. In contrast, our regret bound remains sub-linear for any \(d>0\).

## 6 Proof of the Main Theorem

To perform analysis for continuous actions as in our setting, we use a discretization technique. At each time \(t\), we use a discretization \(\mathcal{D}_{t}\subset\mathcal{D}\) with the property that \(|f(\mathbf{x})-f([\mathbf{x}]_{t})|\leq\frac{1}{t^{2}}\) where \([\mathbf{x}]_{t}\in\mathcal{D}_{t}\) is the closest point to \(\mathbf{x}\in\mathcal{D}\). Then we choose \(\mathcal{D}_{t}\) with size \(|\mathcal{D}_{t}|=\left((b-a)BC_{\text{lip}}t^{2}\right)^{d}\) that satisfies \(\|\mathbf{x}-[\mathbf{x}]_{t}\|_{1}\leq\frac{b-a}{(b-a)BC_{\text{lip}}t^{2}}= \frac{1}{BC_{\text{lip}}t^{2}}\) for all \(\mathbf{x}\in\mathcal{D}\), where \(C_{\text{lip}}=\underset{\mathbf{x}\in\mathcal{D}j\in[d]}{\text{supsup}}\left( \frac{\partial^{2}k_{\text{NTK}}(\mathbf{p},\mathbf{q})}{\partial\mathbf{p}_{j }\mathbf{q}_{j}}|\mathbf{p}=\mathbf{q}=\mathbf{x}\right)\). This implies, for all \(\mathbf{x}\in\mathcal{D}\):

\[|f(\mathbf{x})-f([\mathbf{x}]_{t})|\leq\|f\|_{\mathcal{H}_{k_{\text{NTK}}}}C_{ \text{lip}}\|\mathbf{x}-[\mathbf{x}]_{t}\|_{1}\leq BC_{\text{lip}}\frac{1}{BC _{\text{lip}}t^{2}}=1/t^{2}, \tag{1}\]

is Lipschitz continuous of any \(f\in\mathcal{H}_{k_{\text{NTK}}}(\mathcal{D})\) with Lipschitz constant \(BC_{\text{lip}}\), where we have used the inequality \(\|f\|_{\mathcal{H}_{k_{\text{NTK}}}}\leq B\) which is our assumption about function \(f\). We bound the regret of the proposed algorithm by starting from the instantaneous regret \(r_{t}=f(\mathbf{x}^{*})-f(\mathbf{x}_{t})\) can be decomposed to \(r_{t}=[f(\mathbf{x}^{*})-f([\mathbf{x}^{*}]_{t})]+[f([\mathbf{x}^{*}]_{t}-f( \mathbf{x}_{t})]\). While the first term is bounded by Eq(1), we bound the second term in Lemma 7 provided in our Appendix. The proof of Lemma 7 requires a few steps - we introduce a saturated set \(\mathcal{S}_{t}\) (see Definition A.2), we then combine it with the results of Lemma 2 (deriving a bound on \(|\widetilde{f}_{t}(\mathbf{x})-h(\mathbf{x};\boldsymbol{\theta}_{t-1})|\)) necessitated due to the Thompson sampling and Lemma 3 (deriving a bound on \(|f(\mathbf{x})-h(\mathbf{x};\boldsymbol{\theta}_{t-1})|\)) utilizing
the generalization result of the over-parameterized neural networks. It is noted that in our proof, we consider the effect of our general search space setting (where input \(\mathbf{x}\in\mathcal{D}\) is norm-bounded, i.e., \(0<a\leq\left\|\mathbf{x}\right\|_{2}\leq b\)) on the output of the over-parameterized deep neural network at each layer in Lemma 1 and on the gradient of this network utilized in Lemma 3.2. These results contribute to the appearance of the lower and the upper norm bounds \(a\) and \(b\), the first in the generalization property of over-parameterized neural networks in Lemma 3 and the later affect the cumulative regret bound in Theorem 1. Our proof follows the proof style of Lemma 4.1 of [30], Lemma B.3 of [30] and Theorem 5 of [29].

On the other hand, Lemma 3.3, based on Assumption 2, provides a tighter bound on the confidence interval that eliminated the factor \(\sqrt{\gamma_{T}}\), in comparison with previous relevant works [10, 15], leads to the sub-linear regret bound in Theorem 1. By these arguments, we achieve a cumulative regret bound for our proposed algorithm in terms of the maximum information gain \(\gamma_{T}\) (Lemma 6, 7, 8). **Our detailed proof is provided in Appendix A**.

Proof Sketch for Theorem 1.: With probability at least \(1-\alpha\), we have

\[R_{T} =\sum_{t=1}^{T}f(\mathbf{x}^{*})-f(\mathbf{x}_{t})\] \[=\sum_{t=1}^{T}\left[f(\mathbf{x}^{*})-f([\mathbf{x}^{*}]_{t}) \right]+\left[f([\mathbf{x}^{*}]_{t})-f(\mathbf{x}_{t})\right]\] \[\leq 4T\epsilon(m)+\frac{(2B+1)\pi^{2}}{6}+\bar{C_{1}}(1+c_{T}) \nu_{T}\sqrt{L}\sum_{i=1}^{T}\min(\sigma_{t}(\mathbf{x}_{t}),B)\] \[\quad+(4+\bar{C_{2}}(1+c_{T})\nu_{T}L+4\epsilon(m))\sqrt{2\log(1 /\alpha)T}\] \[\leq\bar{C_{1}}(1+c_{T})\nu_{T}\sqrt{L}\sqrt{\frac{\lambda BT}{ \log(B+1)}(2\gamma_{T}+1)}+\frac{(2B+1)\pi^{2}}{6}+4T\epsilon(m)\] \[\quad+4\epsilon(m)\sqrt{2\log(1/\alpha)T}+\left(4+\bar{C_{2}}(1+ c_{T})\nu_{T}L\right)\sqrt{2\log(1/\alpha)T}\] \[=\bar{C_{1}}(1+c_{T})\nu_{t}\sqrt{L}\sqrt{\frac{\lambda BT}{\log( B+1)}(2\gamma_{T}+1)}+\frac{(2B+1)\pi^{2}}{6}\] \[\quad+\epsilon(m)(4T+\sqrt{2\log(1/\alpha)T})+(4+\bar{C_{2}}(1+c_ {T})\nu_{t}L)\sqrt{2\log(1/\alpha)T}\]
The first inequality is due to Lemma 7, which provides the bound for cumulative regret \(R_{T}\) in terms of \(\sum_{t=1}^{T}\min(\sigma_{t}(\mathbf{x}_{t}),B)\). The second inequality further provides the bound of term \(\sum_{t=1}^{T}\min(\sigma_{t}(\mathbf{x}_{t}),B)\) due to Lemma 8, while the last equality rearranges addition. Picking \(\eta=(m\lambda+mLT)^{-1}\) and \(J=(1+LT/\lambda)\left(\log(C_{\epsilon,2})+\log(T^{3}L\lambda^{-1}\log(1/ \alpha))\right)\), we have

\[\frac{b}{a}C_{\epsilon,2}(1-\eta m\lambda)^{J}\sqrt{TL/\lambda} \left(4T+\sqrt{2\log(1/\alpha)T}\right)\] \[= \frac{b}{a}C_{\epsilon,2}\left(1-\frac{1}{1+LT/\lambda}\right)^{ J}\left(4T+\sqrt{2\log(1/\alpha)T}\right)\] \[= \frac{b}{a}C_{\epsilon,2}e^{-\left(\log(C_{\epsilon,2})+\log \left(T^{3}L\lambda^{-1}\log(1/\alpha)\right)\right)}\left(4T+\sqrt{2\log(1/ \alpha)T}\right)\] \[= \frac{b}{a}\frac{1}{C_{\epsilon,2}}.T^{-3}L^{-1}\lambda\log^{-1} (1/\alpha)\left(4T+\sqrt{2\log(1/\alpha)T}\right)\leq\frac{b}{a}\]

Then choosing \(m\) that satisfies:

\[\left(\frac{b}{a}C_{\epsilon,1}m^{-1/6}\lambda^{-2/3}L^{3}\sqrt{ \log m}+\left(\frac{b}{a}\right)^{3}C_{\epsilon,3}m^{-1/6}\sqrt{\log m}L^{4}T ^{5/3}\lambda^{-5/3}(1+\sqrt{T/\lambda})\right)\] \[\left(4T+\sqrt{2\log(1/\alpha)T}\right)\leq\left(\frac{b}{a} \right)^{3}\]

We finally achieve the bound of \(R_{T}\) as:

\[R_{T} \leq\bar{C}_{1}(1+c_{T})\nu_{T}\sqrt{L}\sqrt{\frac{\lambda BT}{ \log(B+1)}(2\gamma_{T}+1)}\] \[\quad+(4+\bar{C}_{2}(1+c_{T})\nu_{T}L)\sqrt{2\log(1/\alpha)T}+ \frac{(2B+1)\pi^{2}}{6}+\frac{b(a^{2}+b^{2})}{a^{3}}\]

## 7 Experiments

In this section, we demonstrate the effectiveness of our proposed NeuralBO algorithm compared to traditional BO algorithms and several other BO algorithms on various synthetic benchmark optimisation functions and various real world optimisation problems which have complex and structured inputs. The real world optimization problems consists of (1) generating sensitive samples to detect tampered model trained on MNIST [32] dataset; (2) optimizing control parameters for robot pushing considered in [9]; and
(3) optimizing the number of evaluations needed to find a document that resembles an intended (unknown) target document.

For all experiments, we compared our algorithm with common classes of surrogate models used in black-box optimization, including Gaussian Processes (GPs), Random Forests and Deep Neural Networks (DNNs). For GPs, we have three baselines: GP-UCB [7], GP-TS [10] and GP-EI [33] with the common Squared Exponential (SE) Kernel. Our implementations for GP-based Bayesian Optimization baselines utilize public library GPyTorch 1 and BOTorch 2. Next, we consider SMAC [22] as a baseline using RF as surrogate model. Lastly, we also include two recent DNNs-based works for black-box optimization: DNGO [12] and Neural Greedy [14]. We further describe the implementations for RF and DNN-based baselines as follows:

Footnote 1: [https://gpytorch.ai/](https://gpytorch.ai/)

Footnote 2: [https://botorch.org/](https://botorch.org/)

* DNGO [12] models black-box function by a Bayesian linear regressor, using low-dimensional features (extracted by DNN) as inputs. We run DNGO algorithm with the implementation of AutoML 3 with default settings. Footnote 3: [https://github.com/automl/pybnn](https://github.com/automl/pybnn)
* NeuralGreedy [14] fits a neural network to the current set of observations, where the function values are randomly perturbed before learning the neural network. The learnt neural network is then used as the acquisition function to determine the next query point. Our re-implementation of NeuralGreedy follows the setting described in Appendix F.2[14].
* Random forest (RF) based surrogate model is implemented using public library GPyOpt 4 and optimization is performed with EI acquisition function. Footnote 4: [https://github.com/SheffieldML/GPyOpt](https://github.com/SheffieldML/GPyOpt)

For our proposed NeuralBO algorithm, we model the functions using two-layer neural networks \(f(\mathbf{x};\boldsymbol{\theta})=\sqrt{m}\mathbf{W}_{2}\sigma(\mathbf{W}_{1} \mathbf{x})\) with network width \(m=500\). The weights of the network are initialized with independent samples from a normal distribution \(\mathcal{N}(0,1/m)\). To train the surrogate neural
network models, we use SGD optimizer with batch size \(50\), epochs \(50\), learning rate \(\eta=0.001\) and regularizer \(\lambda=0.01\). We set exploration parameter in Algorithm 1: \(\nu_{t}=\nu\) and do a grid search over \(\{0.1,1,10\}\).

### Synthetic Benchmark Functions

We optimized several commonly used synthetic benchmark functions: Ackley, Levy, and Michalewics. The exact expression of these functions can be found at: [http://www.sfu.ca/~ssurjano/optimization.html](http://www.sfu.ca/~ssurjano/optimization.html). For each function, we perform optimization in two different dimensions. The latter dimension is twice as high as the former. This is done to evaluate the methods across a range of dimensions. More precisely, we choose dimensions \(15\) and \(30\) for Michalewics, \(20\) and \(40\) for Levy, \(50\) and \(100\) for Ackley. The function evaluation noise follows a normal distribution with zero mean and the variance depends on objective function.
All reported experiments are averaged over ten runs, where each run uses a random initialization. All the methods start with the same set of initial points. As seen from Figure 1, Neural-BO outperforms all other baselines, including GP-based BO algorithms, both NN-based BO algorithms (DNGO

Figure 1: The plots show the minimum true value observed after optimizing several synthetic functions over 2000 iterations of our proposed algorithm and 6 baselines. The dimension of each function is shown in the parenthesis.
and NeuralGreedy) and one algorithm using random forest. Moreover, our algorithm is also promising for high dimensions.

### Real Applications

#### 7.2.1 Designing Sensitive Samples for Detection of Model Tampering

We consider an attack scenario where a company offering Machine Learning as a service (MLaaS) hosts its model on a cloud. However, an adversary with backdoor access may tamper with the model to change its weight. This requires the detection of model tampering. To deal with this problem, [34] suggests generating a set of test vectors named _Sensitive-Samples_\(\{v_{i}\}_{i=1}^{n}\), whose outputs predicted by the compromised model will be different from the outputs predicted by the original model. As formalized in [34], suppose we suspect a pre-trained model \(f_{\theta}(x)\) of having been modified by an attacker after it was sent to the cloud service provider. Finding sensitive

Figure 2: The plot shows the **detection rates** corresponding to the number of samples on MNIST dataset. The larger the number of sensitive samples, the higher the detection rate. As shown in the figure, NeuralBO can generate sensitive samples that achieves nearly 90% of the detection rate with at least 8 samples.
samples for verifying the model's integrity is equivalent to optimizing task: \(v=\operatorname*{argmax}_{x}\left\|\frac{\partial f_{\theta}(x)}{\partial\theta} \right\|_{F}\), where \(\left\|\cdot\right\|_{F}\) is the Frobenius norm of a matrix. A _successful detection_ is defined as "given \(N_{S}\) sensitive samples, there is at least one sample, whose top-1 label predicted by the compromised model is different from the top-1 label predicted by the correct model". Clearly, optimizing this expensive function requires a BO algorithm to be able to work with high-dimensional structured images, unlike usual inputs that take values in hyper-rectangles.

We used a hand-written digit classification model (pre-trained on MNIST data) as the original model and tampered it by randomly adding noise to each weight of this model. We repeat this procedure 1000 times to generate 1000 different tampered models. The top-1 accuracy of the MNIST original model is 93% and is reduced to \(87.73\%\pm 0.08\%\) after modifications.

To reduce the computations, we reduce the dimension of images from \(28\times 28\) to \(7\times 7\) and do the optimization process in 49-dimensional space. After finding the optimal points, we recover these points to the original dimension by applying an upscale operator to get sensitive samples. We compare our method with all other baselines by the average detection rate with respect to the number of sensitive samples. From Figure 2, it can be seen that our method can generate samples with better detection ability than other baselines. This demonstrates the ability of our method to deal with complex structured data such as images.
#### 7.2.2 Unknown target document retrieval

Next, we evaluate our proposed method on a retrieval problem where our goal is to retrieve a document from a corpus that matches user's preference. The optimization algorithm works as follows: it retrieves a document, receives user's feedback score and then updates its belief and attempts again until it reaches a high score, or its budget is depleted. The objective target function is defined as the user's evaluation of each document, which usually has a complex structure. It is clear that evaluations are expensive since the user must read each suggestion. Searching for the most relevant document is considered as finding document \(d\) in the dataset \(S_{text}\) that maximizes the match to the target document \(d_{t}\): \(d=\operatorname*{argmax}_{d\in S_{text}}\operatorname{Match}(d,d_{t})\), where \(\operatorname{Match}(d,d_{t})\) is a matching score between documents \(d\) and \(d_{t}\). We represent each document by a word frequency vector \(x_{n}=(x_{n1},\cdots,x_{nJ})\), where \(x_{nj}\) is the number of occurrences of the \(j\)-th vocabulary term in the \(n\)-th document,

Figure 3: We search for the most related document for a specified target document in **Amazon product reviews** dataset and report the maximum **hierarchical F1 score** found by all baselines. All methods show similar behaviour and NeuralBO performs comparably and much better than GP-based baselines.
and \(J\) is the vocabulary size.

Our experiment uses Amazon product reviews dataset 5, which are taken from users' product reviews from Amazon's online selling platform. This dataset has hierarchical categories, where the category classes are sorted from general to detail. The dataset are structured as: 6 classes in "level 1", 64 classes in "level 2" and 464 classes in "level 3". The number of users' reviews was originally 40K, which was reduced to 37738 after ignoring reviews with "unknown" category. We choose the size of vocabulary to be 500 and use hierarchical F1-score introduced in [35] as a scoring metric for the target and retrieved documents. We report the mean and standard deviation of hierarchical F1-score between target and retrieved documents over ten runs for all methods in Figure 3. Figures 3 indicate that our method shows a better performance for the Amazon product review dataset in comparison with other approaches.

Footnote 5: [https://www.kaggle.com/datasets/kashnitsky/hierarchical-text-classification](https://www.kaggle.com/datasets/kashnitsky/hierarchical-text-classification)
#### 7.2.3 Optimizing control parameters for robot pushing

Finally, we evaluate our method on tuning control parameters for robot pushing problem considered in [9]. We run each method for a total of 6000 evaluations and repeat ten times to take average optimization results. Neural-BO and all other methods are initialized with 15 points. Figure 4 summarizes the median of the best rewards achieved by all methods. It can be seen that Neural-BO achieves the highest reward after 5K optimization steps.

We also use this task to demonstrate the GPU memory used by three practical baselines: GP-UCB, Neural-BO, and NeuralGreedy in Figure 5. The NeuralGreedy and NeuralBO GPU memory requirements do not appear to have changed as the memory used by these two methods remains almost constant, unlike GP. For GP, GPU memory requirement grows fast with the iterations due to the need to compute the inverse of larger and larger covariance matrices. We report here GP-UCB, but similar behavior is also noted for GP-EI and GP-TS. Neural Greedy is slightly lesser memory-consuming

Figure 4: Optimization results for control parameters of 14D robot pushing problem. The X-axis shows iterations, and the y-axis shows the median of the best reward obtained.
because it only needs memory to store the neural network model and for the training process. Neural-BO used a little bit more memory due to the calculation of matrix \(\mathbf{U}\) described in Algorithm 1. However, this matrix has a fixed size; therefore, the memory consumption does not increase over time. We do not report memory used by RF and DNGO because their implementations are CPU-based, and their running speeds are too slow compared to all the above baselines.

## 8 Conclusion

We proposed a new algorithm for Bayesian optimization using deep neural networks. A key advantage of our algorithm is that it is computationally

Figure 5: The GPU memory used by three practical baselines for 14D robot pushing control parameters optimization task.
efficient and performs better than traditional Gaussian process based methods, especially for complex structured design problems. We provided rigorous theoretical analysis for our proposed algorithm and showed that its cumulative regret converges with a sub-linear regret bound. Using both synthetic benchmark optimization functions and a few real-world optimization tasks, we showed the effectiveness of our proposed algorithm.

Our proposed method currently employs over-parametric neural networks, which can be analyzed using neural tangent kernels. As a future work, it would be interesting to investigate if it is possible to extend this work beyond over-parametric neural networks to allow a broader class of neural networks that may have smaller widths.

## 9 Acknowledgement

This research was partially funded by the Australian Government through the Australian Research Council (ARC). Prof Venkatesh is the recipient of an ARC Australian Laureate Fellowship (FL170100006).
