Variational Search Distributions

## Abstract

Abstract We develop variational search distributions (VSD) , a method for finding discrete, combinatorial designs of a rare desired class in a batch sequential manner with a fixed experimental budget. We formalize the requirements and desiderata for this problem and formulate a solution via variational inference. In particular, VSD uses off-the-shelf gradient based optimization routines, can learn powerful generative models for designs, and can take advantage of scalable predictive models. We derive asymptotic convergence rates for learning the true conditional generative distribution of designs with certain configurations of our method. After illustrating the generative model on images, we empirically demonstrate that VSD can outperform existing baseline methods on a set of real sequence-design problems in various biological systems.

## 1 Introduction

We consider a variant of the *active search* problem (Garnett et al., ; Jiang et al., ; Vanchinathan et al., ), where we wish to find as many members (designs) of a rare desired class in a batch sequential manner with a fixed experimental budget. Examples of this are compounds that could be useful pharmaceutical drugs, or highly active enzymes for catalyzing chemical reactions. We assume the design space is discrete or partially discrete, high-dimensional, and practically *innumerable*.
For example, the number possible configurations of a single protein is $20^{{\mathcal{{O}}}(100)}$ (see, e.g., Sarkisyan et al., ).

We are interested in this objective for a variety of reasons. We may wish to study the properties of the “fitness landscape” (Papkou et al., ) to gain a better scientific understanding of a phenomenon such as natural evolution. Or, we may not be able to completely specify the constraints and objectives of a task, but we would like to characterize the space of feasible designs. For example, we want enzymes that can degrade plastics in an industrial setting, but we may not yet know the exact conditions (e.g. temperature, pH), some of which may be anti-correlated with enzyme catalytic activity.

Assuming we can take advantage of a prior distribution over designs, we formulate the search problem as inferring the posterior distribution over rare, desirable designs. Importantly, this posterior can be used for generating new designs. Specifically, we use (black-box) variational inference (VI) (Ranganath et al., ), and so refer to our method as variational search distributions (VSD). Our major contributions are: (1) we formulate the batch active search objective over a (practically) innumerable discrete design space as an instance of variational inference, (2) we present an algorithm, VSD, which solves this objective, and (3) we show that
VSD satisfies well-defined requirements and desiderata specific to our problem. For example, it uses off-the-shelf gradient based optimization routines, is able to learn powerful generative models, and can take advantage of scalable predictive models. In our experiments we show that VSD can outperform existing baseline methods on a set of real applications.
Finally, we evaluate our approach on the related sequential black-box optimization (BBO) problem, where we want to find the globally optimal design for a specific objective and show competitive performance when compared with state-of-the-art methods, e.g., based on Bayesian optimization (BO) (Garnett, ).

## 2 Method

In this section we formalize our problem and describe its requirements and desiderata. We also develop our proposed solution, based on variational inference, which we will refer to as variational search distributions (VSD).

### 2.1 Problem Formulation

We are given a design space $\mathcal{X}$, which can be discrete or mixed discrete-continuous and high dimensional, and where for each instance that we choose $\mathbf{x}\in\mathcal{X}$, we measure some corresponding property of interest (so-called fitness) $y\in$. For example, in our motivating application of DNA/RNA or protein sequences (henceforth referred to as just sequences), $\mathcal{X}=\mathcal{V}^{M}$ where $\mathcal{V}$ is the sequence vocabulary (e.g., amino acid labels, $|\mathcal{V}|=20$) and $M$ is the length of the sequence. However, we do not limit the application of our method to sequences. Using this framing, a real world experiment (for example, measuring the activity of an enzyme) can be modeled as an unknown relationship,

$$ $\displaystyle y={f\!\centerdot}(\mathbf{x})+\epsilon,$ (1) $$

for some black-box function (e.g. the experiment), ${f\!\centerdot}$, and measurement error $\epsilon\in$, distributed according to $p({\epsilon})$ with $\mathbb{E}_{p({\epsilon})}\!\left[{\epsilon}\right]=0$.
Instead of wanting to model the whole space, we are only interested in a set of events which we choose based on fitness $y$. For instance (refer to Figure 1), we may be interested in the fittest measurable design; all designs above a minimum level of feasibility (e.g. a wild-type sequence), $\tau\in$; the distribution of these feasible designs; or the shape of the black-box function for these feasible designs,

$$ $\displaystyle\mathbf{x}^{*}=\operatorname*{arg\,\!max}_{\mathbf{x}}{f\! \centerdot}(\mathbf{x})+\epsilon,\quad\mathcal{S}:=\{\mathbf{x}:y>\tau\},\quad p ({\mathbf{x}}|{y>\tau}),\quad\textrm{or}\quad\mathcal{F}:=\{{f\!\centerdot}( \mathbf{x}):\mathbf{x}\in\mathcal{S}\}.$ (2) $$

Our primary focus in this work is to estimate the super level-set distribution $p({\mathbf{x}}|{y>\tau})$ in a sequential manner.
We assume that $\mathcal{S}$ are rare events in a high dimensional space, and that we have access to a prior belief, $p({\mathbf{x}})$, which helps narrow in on this subset of $\mathcal{X}$.
We are given a dataset, $\mathcal{D}_{N}:=\{(y_{n},\mathbf{x}_{n})\}_{n=1}^{N}$, which may contain only a few instances of $y_{n}>\tau$. Given $p({\mathbf{x}})$ and $\mathcal{D}_{N}$ we aim to recommend batches of unique candidates, $\{\mathbf{x}_{bt}\}^{B}_{b=1}$, for experimental evaluation (Equation 1) in a series of rounds, $t\in\{1,\ldots,T\}$, where $B={\mathcal{{O}}}(1000)$ and we desire $\mathbf{x}_{bt}\in\mathcal{S}$. Each round, $\mathcal{D}_{N}$ is augmented with the experimental results of the previous batch, so $N\leftarrow N+B$.
Estimating this super level-set distribution of $\mathbf{x}$ is computationally and statistically challenging and, therefore, we cast this as a *variational inference* problem. As we shall see later, our solution allows us to satisfy the following requirements and additional desiderata for our problem.

###### Requirements & Desiderata .

Problem requirements (R) and other desiderata (D).

feasible designs, $y>\tau$, are rare events in $\mathcal{X}$ that need to be identified

non-myopic candidate generation, $\mathbf{x}\in\mathcal{S}\subset\mathcal{X}$, for sequential experiments

search over (combinatorially) large design spaces, e.g. $\mathbf{x}\in\mathcal{X}=\mathcal{V}^{m}$

generation of up to $\mathcal{O}(1000)$ *diverse* candidate designs per round

convergence for certain choices of priors, variational distributions and predictive models

based optimization strategies for candidate searching

models, $\mathbf{x}^{(s)}\sim q({\mathbf{x}})$, that are task-specific for fit designs

predictive models that enable high-throughput experiments.

Figure: (a) $\operatorname*{arg\,\!max}_{\mathbf{x}}y(\mathbf{x})$
Refer to caption: extracted/5897346/figs/argmaxx.png

###### Requirements & Desiderata .

### 2.2 Variational Search Distributions

We cast the estimation of $p({\mathbf{x}}|{y>\tau})$ as a sequential optimization problem. A suitable objective for a round, $t$, is to minimize a divergence,

$$ $\displaystyle\phi_{t}^{*}=\operatorname*{arg\,\!min}_{\phi}\mathbb{D}\!\left[{ p({\mathbf{x}}|{y>\tau})}\|{q({\mathbf{x}}|{\phi})}\right]$ (3) $$

where $q({\mathbf{x}}|{\phi})$ is a parameterized distribution from which we sample experimental candidate designs $\mathbf{x}_{bt}$, (D3) Generative, and which we aim to match to $p({\mathbf{x}}|{y>\tau})$.
The difficulty is that we cannot directly evaluate or empirically sample from $p({\mathbf{x}}|{y>\tau})$. However, if we consider the reverse Kullback-Leibler (KL) divergence,

$$ $\displaystyle\operatorname*{arg\,\!min}_{\phi}\,\mathbb{D}_{\mathrm{KL}}\! \left[{q({\mathbf{x}}|{\phi})}\|{p({\mathbf{x}}|{y>\tau})}\right]= \operatorname*{arg\,\!min}_{\phi}\mathbb{E}_{q({\mathbf{x}}|{\phi})}\!\left[{ \log\frac{q({\mathbf{x}}|{\phi})}{p({\mathbf{x}})}-\log p({y>\tau}|{\mathbf{x} })}\right],$ (4) $$

where we have expanded $p({\mathbf{x}}|{y>\tau})$ using Bayes rule and dropped the constant term $p({y>\tau})$, we note that we no longer require evaluation of $p({\mathbf{x}}|{y>\tau})$ directly. We recognize the right hand side of Equation 4 as the well known (negative) variational evidence lower bound (ELBO),

$$ $\displaystyle\mathcal{L_{\textrm{ELBO}}}\!\left({\phi}\right):=\mathbb{E}_{q({ \mathbf{x}}|{\phi})}\!\left[{\log p({y>\tau}|{\mathbf{x}})}\right]-\mathbb{D}_ {\mathrm{KL}}\!\left[{q({\mathbf{x}}|{\phi})}\|{p({\mathbf{x}})}\right].$ (5) $$

For this we assume access to a prior distribution over the space of designs, $p({\mathbf{x}})$, that may be informed from the data at hand. Henceforth, as we will develop a sequential algorithm, we will denote this prior with $p({\mathbf{x}}|{\mathcal{D}_{0}})$.
We note the relationship between $\log p({y>\tau}|{\mathbf{x}})$ and the probability of improvement (PI) acquisition function from BO (Kushner, ),

$$ $\displaystyle\log p({y>\tau}|{\mathbf{x}}):=\log\mathbb{E}_{p({y}|{\mathbf{x}, \mathcal{D}_{N}})}\!\left[{\mathbbm{1}[{y>\tau}]}\right]=\log\alpha_{PI}\! \left({\mathbf{x},\mathcal{D}_{N},\tau}\right).$ (6) $$

Here $\mathbbm{1}:\{\textrm{false},\textrm{true}\}\to\{0,1\}$ is the indicator function and $p({y}|{\mathbf{x},\mathcal{D}_{N}})$ is typically estimated using the posterior predictive distribution of a Gaussian process (GP) given data, $\mathcal{D}_{N}$. So $p({y>\tau}|{\mathbf{x},\mathcal{D}_{N}})=\Psi((\mu_{N}(\mathbf{x})-\tau)/
\sigma_{N}(\mathbf{x}))$, where $\Psi(\cdot)$ is a cumulative standard normal distribution function, and $\mu_{N}(\mathbf{x})$, $\sigma_{N}^{2}(\mathbf{x})$ are the posterior predictive mean and variance, respectively, of the GP. We refer to this estimation strategy as GP-PI, and rewrite the ELBO accordingly,

$$ $\displaystyle\mathcal{L_{\textrm{ELBO}}}\!\left({\phi,\tau,\mathcal{D}_{N}} \right)=\mathbb{E}_{q({\mathbf{x}}|{\phi})}\!\left[{\log\alpha_{PI}\!\left({ \mathbf{x},\mathcal{D}_{N},\tau}\right)}\right]-\mathbb{D}_{\mathrm{KL}}\! \left[{q({\mathbf{x}}|{\phi})}\|{p({\mathbf{x}}|{\mathcal{D}_{0}})}\right].$ (7) $$

The method that maximizes the objective in Equation 7 we call variational search distributions (VSD), since we are using the variational posterior distribution as a means of searching the space of fit designs, satisfying (R1) Rare, (R2) Sequential and (R4) Batch. It is well known that when the true posterior is a member of the variational family indexed by $\phi$, the above variational inference procedure has the potential to recover the exact posterior distribution. To recommend experimental candidates we sample a set of designs from our search distribution each round,

$$ $\displaystyle\{\mathbf{x}_{bt}\}^{B}_{b=1}\sim\prod^{B}_{b=1}q({\mathbf{x}}|{ \phi_{t}^{*}}),\quad\textrm{where}\quad\phi_{t}^{*}=\operatorname*{arg\,\!max} _{\phi}\mathcal{L_{\textrm{ELBO}}}\!\left({\phi,\tau,\mathcal{D}_{N}}\right).$ (8) $$

We discuss the relationship between VSD and BO in Appendix F. In general, because of the discrete combinatorial nature of our problem, we cannot use the re-parameterization trick (Kingma & Welling, ) to estimate the gradients of the ELBO. Instead, we use the score function gradient estimator (Williams, ; Mohamed et al., ) with standard gradient descent methods (D2) Gradient,

$$ $\displaystyle\nabla_{\phi}\mathcal{L_{\textrm{ELBO}}}\!\left({\phi,\tau, \mathcal{D}_{N}}\right)$ $\displaystyle=\mathbb{E}_{q({\mathbf{x}}|{\phi})}\!\left[{\left(\log\alpha_{PI }\!\left({\mathbf{x},\mathcal{D}_{N},\tau}\right)-\log\frac{q({\mathbf{x}}|{ \phi})}{p({\mathbf{x}}|{\mathcal{D}_{0}})}\right)\nabla_{\phi}\log q({\mathbf{ x}}|{\phi})}\right],$ (9) $$

where we use Monte-Carlo sampling to approximate this expectation with a suitable variance reduction scheme, such as using a control variate or baseline (Mohamed et al., ). We find that the exponentially smoothed average of the ELBO works well in practice, and is the same strategy employed in Daulton et al. (). Effectively, VSD implements black-box variational inference (Ranganath et al., ) for parameter estimation, and despite the high-dimensional nature of $\mathcal{X}$, we find we only need ${\mathcal{{O}}}(1000)$ samples to estimate the required expectations for ELBO optimization on problems with $M=\mathcal{O}(100)$, satisfying (R3) Discrete. Note that Equation 7 – 9 do not involve any data ($\mathcal{D}_{N}$) directly, only indirectly through the acquisition function. Hence the scalability of VSD is dependent on the complexity of training the underlying estimator of $p({y}|{\mathbf{x},\mathcal{D}_{N}})$.

### 2.3 Class Probability Estimation

So far our method indirectly computes the PI by transforming the predictions of a GP surrogate model, $p({y}|{\mathbf{x},\mathcal{D}_{N}})$, as in Equation 6. Instead we may choose to follow the reasoning used by Bayesian optimization by density-ratio estimation (BORE) in Tiao et al. (); Oliveira et al. (); Song et al. (), and directly estimate the quantity we care about, $p({y>\tau}|{\mathbf{x},\mathcal{D}_{N}})$. We do this with class probability estimation (CPE) on the labels $z:=\mathbbm{1}[{y>\tau}]\in\{0,1\}$ so
$p({y>\tau}|{\mathbf{x},\mathcal{D}_{N}})=p({z=1}|{\mathbf{x},\mathcal{D}_{N}})
\approx\pi_{\theta}({\mathbf{x}})$, where $\pi_{\theta}:\mathcal{X}\to[0,1]$. We can recover the class probability
estimates using a proper scoring rule (Gneiting & Raftery, ) such as
Brier score or log-loss on training data, $\mathcal{D}_{N}^{z}=\{(z_{n},\mathbf{x}_{n})\}^{N}_{n=1}$, e.g.,

$$ $\displaystyle\mathcal{L_{\textrm{CPE}}}\!\left({\theta,\mathcal{D}_{N}^{z}} \right):=-\frac{1}{N}\sum_{n=1}^{N}z_{n}\log\pi_{\theta}({\mathbf{x}_{n}})+(1- z_{n})\log(1-\pi_{\theta}({\mathbf{x}_{n}})).$ (10) $$

The VSD objective using CPE becomes,

$$ $\displaystyle\mathcal{L_{\textrm{ELBO}}}\!\left({\phi,\theta}\right)=\mathbb{E }_{q({\mathbf{x}}|{\phi})}\!\left[{\log\pi_{\theta}({\mathbf{x}})}\right]- \mathbb{D}_{\mathrm{KL}}\!\left[{q({\mathbf{x}}|{\phi})}\|{p({\mathbf{x}}|{ \mathcal{D}_{0}})}\right],$ (11) $$

into which we plug $\theta_{t}^{*}=\operatorname*{arg\,\!min}_{\theta}\mathcal{L_{\textrm{CPE}}}\!
\left({\theta,\mathcal{D}_{N}^{z}}\right)$. Using a CPE also opens up the choice of estimators that are more scalable than GP-PI, satisfying our desiderata (D4) Scalable. This may be crucial if we choose to run more than a few rounds of experiments with $B={\mathcal{{O}}}(1000)$. Additionally, since VSD is a black-box method, we can choose to use CPEs that are non-differentiable, such as decision tree ensembles.
The complete VSD algorithm is given in Algorithm 1, in which we have allowed for a threshold function, $\tau_{t}=f_{\tau}\!\left({\{y:y\in\mathcal{D}_{N}\},\gamma_{t}}\right)$. This function can be used to modify the threshold each round, e.g. following Tiao et al. (), an empirical quantile function $\tau_{t}=\hat{Q}_{y}(\gamma_{t})$ where $\gamma_{t}\in(0,1)$, or a constant $\tau$ in the case of estimating the distribution of the super level-set.

Figure: Algorithm 1 VSD optimization loop with CPE.

### 2.4 Theoretical Analysis

In this section we summarize the main theoretical results concerning VSD and its estimates. We show that VSD sampling distributions converge to a target distribution that characterizes the level set given by $\tau$, satisfying (D1) Guaranteed in two general settings.
We first derive results assuming ${f\!\centerdot}$ is drawn from a Gaussian process, i.e., ${f\!\centerdot}\sim\mathcal{GP}(0,k)$, with a positive-semidefinite covariance (or kernel) function $k:\mathcal{X}\times\mathcal{X}\to$ (Appendix D) and then using GP-PI as the CPE for VSD. These results are then extended to probabilistic classifiers based on wide neural networks (Appendix E) by means of the neural tangent kernel (NTK) for the given architecture (Jacot et al., ).
For simplicity, we set $B=1$ and $N=t$, though sampling with $B>1$ should improve rates by a multiplicative factor.
{restatable}theoremmainthm
Let assumptions D.1 to D.5 hold. Then VSD equipped with GP-PI approaches the level-set distribution at the following rate:

$$ $\begin{split}\mathbb{D}\!\left[{p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})} \|{p({\mathbf{x}}|{y>\tau_{t},{f\!\centerdot}})}\right]\in{\mathcal{{O}}}_{ \mathbb{P}}(t^{-1/2})\end{split}\,.$ $$

This result is based on showing that the GP posterior variance vanishes at an optimal rate of ${\mathcal{{O}}}(t^{-1})$ in our setting (D.5).
We also analyze the rate at which VSD finds feasible designs, or “hits”, compared to an oracle with full knowledge of ${f\!\centerdot}$.
After $T$ rounds, the number of hits found by VSD is
$H_{T}=\sum_{t=1}^{T}\mathbbm{1}[{y_{t}>\tau_{t-1}}]$, where $y_{t}$ is generated from Equation 1 and $\mathbf{x}_{t}\sim p({\mathbf{x}}|{y>\tau_{t-1},\mathcal{D}_{t-1}})$. The number of hits, $H_{T}^{*}$, from an agent that fully knows ${f\!\centerdot}$ is the same but for generating conditioned on ${f\!\centerdot}$ with $\mathbf{x}_{t}\sim p({\mathbf{x}}|{y>\tau_{t-1},{f\!\centerdot}})$.
Using this definition and Sec. 2.4 we show the following.
{restatable}corollaryhitscor
Under the settings in Sec. 2.4, we also have that:

$$ $\mathbb{E}[|H_{T}-H_{T}^{*}|]\in{\mathcal{{O}}}(\sqrt{T})\,.$ $$

$\mathbb{E}[H_{T}]$ is related to the empirical recall measure in Equation 16 up to the normalization constant, but it does not account for repeated hits, which are treated as false discoveries (false positives) under recall. Lastly, for NN-based CPEs, we obtain convergence rates dependent on the spectrum of the NTK (E.2), which we instantiate for ReLU networks below. For the full results and proofs, please see Appendix D for the GP-based analysis and Appendix E for the NTK results.
{restatable}corollaryrelukl
Let $\pi_{\theta}$ be modeled via a fully connected ReLU network. Then, under the assumptions in E.2, VSD achieves:

$$ $\mathbb{D}\!\left[{p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})}\|{p({\mathbf{ x}}|{y>\tau,{f\!\centerdot}})}\right]\in\widetilde{{\mathcal{{O}}}}_{\mathbb{P }}\left(t^{-\frac{1}{2(M+1)}}\right),$ (12) $$

which asymptotically vanishes for all finite sequence lengths $M$.

## 3 Related Work

VSD can be viewed as one of many methods that makes use of the bound (Staines & Barber, ),

$$ $\displaystyle\max_{\mathbf{x}}{f\!\centerdot}(\mathbf{x})\geq\max_{\phi} \mathbb{E}_{q({\mathbf{x}}|{\phi})}\!\left[{{f\!\centerdot}(\mathbf{x})}\right].$ (13) $$

The maximum is always greater than or equal to the expected value of a random variable. This bound is useful for black-box optimization (BBO)of ${f\!\centerdot}$, and becomes tight if $q({\mathbf{x}}|{\phi})\to\delta(\mathbf{x}^{*})$, see Appendix F for more detail and VSD’s relation. Other well known methods that make use of this bound are evolution strategies (ES) and natural evolution strategies (NES) (Wierstra et al., ), variational optimization (VO) (Staines & Barber, ; Bird et al., ), estimation of distribution algorithms (EDA) (Larrañaga & Lozano, ), and Bayesian optimization with probabilistic reparameterisation (BOPR) (Daulton et al., ). For learning the parameters of the variational distribution, $\phi$, they variously make use of maximum likelihood estimation or the score function gradient estimator (REINFORCE) (Williams, ).
Algorithms that modify Equation 13 to stop the collapse of $q({\mathbf{x}}|{\phi})$ to a point mass for batch design include design by adaptive sampling (DbAS) (Brookes & Listgarten, ) and conditioning by adaptive sampling (CbAS) (Brookes et al., ). They use fixed samples $\mathbf{x}^{(s)}$ from $q({\mathbf{x}}|{\phi^{*}_{t-1}})$ for approximating the expectation, and then optimize $\phi$ using a weighted maximum-likelihood or variational style procedure. We can take a unifying view of many of these algorithms by recognizing the general gradient estimator,

$$ $\displaystyle\mathbb{E}_{q({\mathbf{x}}|{\phi^{\prime}})}\!\left[{w(\mathbf{x} )\nabla_{\phi}\log q({\mathbf{x}}|{\phi})}\right],$ (14) $$

where we give each component in Table 1. BORE ^∗ has been adapted to discrete $\mathcal{X}$ by using the score function gradient estimator and CbAS and DbAS have been adapted to use a CPE– their original derivations use the equivalent of a PI acquisition function.

**Table 1: How related methods can be adapted from Equation 14. VSD, CbAS and DbAS may also use a cumulative distribution representation of $\alpha_{\textrm{PI}}\!\left({\mathbf{x},\mathcal{D}_{N},\tau}\right)$ in place of $\pi_{\theta^{*}}({\mathbf{x}})$.**
| Method | $w(\mathbf{x})$ | $\phi^{\prime}$ | Fixed $\mathbf{x}^{(s)}\sim q({\mathbf{x}}|{\phi^{\prime}})$? |
| --- | --- | --- | --- |
| VSD | $\log\pi_{\theta^{*}}({\mathbf{x}})+\log p({\mathbf{x}}|{\mathcal{D}_{0}})-\log<br>q<br>({\mathbf{x}}|{\phi})$ | $\phi$ | No |
| CbAS | $\pi_{\theta^{*}}({\mathbf{x}})p({\mathbf{x}}|{\mathcal{D}_{0}})/q({\mathbf{x}}<br>|{\phi^{*}_{t-1}})$ | $\phi^{*}_{t-1}$ | Yes |
| DbAS | $\pi_{\theta^{*}}({\mathbf{x}})$ | $\phi^{*}_{t-1}$ | Yes |
| BORE ^∗ | $\pi_{\theta^{*}}({\mathbf{x}})$ | $\phi$ | No |
| BOPR | $\alpha\!\left({\mathbf{x},\mathcal{D}_{N}}\right)$ | $\phi$ | No |

A number of finite horizon methods have been applied to biological sequence BBO tasks, such as Amortized BO (Swersky et al., ), GFlowNets (Jain et al., ), and the reinforcement learning based DynaPPO (Angermueller et al., ). Heuristic stochastic search methods such as AdaLead (Sinai et al., ) and proximal exploration (PEX) (Ren et al., ) have also demonstrated strong empirical performance on these tasks. We compare the properties of the most relevant methods to our problem in Table 2.

Other methods that explicitly attempt to estimate the super level-set distribution, $p({\mathbf{x}}|{y>\tau})$, include CbAS, which optimizes the forward KL divergence, $\mathbb{D}_{\mathrm{KL}}\!\left[{p({\mathbf{x}}|{y>\tau})}\|{q({\mathbf{x}}|{
\phi})}\right]$, in a sequential fashion using importance weighted cross entropy estimation (Rubinstein, ). Batch-BORE (Oliveira et al., ) also optimizes the reverse KL divergence and uses CPE, but with Stein variational inference (Liu & Wang, ) for continuous and diverse batch candidates $\{\mathbf{x}_{tb}\in^{M}\}^{B}_{b=1}$ (the authors optimize on the simplex, and then choose the maximum coordinate for their discrete experiments). There is a rich literature on the related task of active learning and BO for level set estimation (LSE) (Bryan et al., ; Gotovos et al., ; Bogunovic et al., ; Zhang et al., ). However we focus on the task of learning a generative model over a practically innumerable discrete space.

**Table 2: Feature table of competing methods: ✓ has feature, ✗ does not have feature, – partially has feature, or requires simple modification. We follow Swersky et al. () in their definition of amortization referring to the ability to use $q({\mathbf{x}}|{\phi^{*}_{t-1}})$ for warm-starting the optimization of $\phi_{t}$.**
| Method | Rare $\mathbf{x}\in\mathcal{S}$ (R1) Rare | Sequential (R2) Sequential | Discrete $\mathcal{X}$ (R3) Discrete | Batch $\{\mathbf{x}_{bt}\}^{B}_{b=1}$ (R4) Batch | Guaranteed (D1) Guaranteed | Gradient descent (D2) Gradient | Generative $q({\mathbf{x}}|{\phi})$ (D3) Generative | Scalable (D4) Scalable | General acq./reward fn. | Amortization |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| BOPR (Daulton et al., ) | ✗ | ✓ | ✓ | ✗ | ✓ | ✓ | – | ✗ | ✓ | – |
| BORE (Tiao et al., ) | ✗ | ✓ | – | ✗ | ✓ | ✓ | – | ✓ | ✗ | – |
| Batch BORE (Oliveira et al., ) | ✓ | ✓ | ✗ | ✓ | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ |
| DbAS (Brookes & Listgarten, ) | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | ✓ | ✗ | ✓ |
| CbAS (Brookes et al., ) | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | ✓ | ✗ | ✓ |
| Amortized BO (Swersky et al., ) | ✗ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | ✓ | ✓ | ✓ |
| GFlowNets (Jain et al., ) | ✗ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | ✓ | ✓ | ✓ |
| DynaPPO (Angermueller et al., ) | ✗ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | – | ✓ | ✓ |
| AdaLead (Sinai et al., ) | ✗ | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | – | ✗ | ✗ |
| PEX (Ren et al., ) | ✗ | ✓ | ✓ | ✓ | ✗ | ✗ | ✗ | – | ✗ | ✗ |
| GGS (Kirjner et al., ) | ✗ | ✗ | ✓ | ✓ | ✗ | ✗ | ✓ | ✗ | ✗ | ✗ |
| VSD (ours) | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ |

## 4 Experiments

Firstly we test our method, VSD, on its ability to generate complex, structured candidates, $\mathbf{x}$, in a single round by training it to generate a subset of handwritten digits from flattened MNIST images (LeCun et al., ) in Sec. 4.1.
We then compare VSD on two sequence design tasks against existing baseline methods.
The first of these tasks (Sec. 4.2) is to discover as many unique, fit sequences as possible using the datasets DHFR (Papkou et al., ), TrpB (Johnston et al., ) and TFBIND8 (Barrera et al., ). These datasets contain near complete evaluations of $\mathcal{X}$, and to our knowledge DHFR and TrpB are novel in the machine learning literature. The second (Sec. 4.3) is a more traditional black-box optimization task of finding the maximum of an unknown function; using datasets AAV (Bryant et al., ) and GFP (Sarkisyan et al., ). A complete combinatorial assessment is infeasible for these latter datasets, and so we use the convolution neural network oracle presented in Kirjner et al. () as *in-silico* ground truth.
The corresponding datasets involve $|\mathcal{V}|\in\{4,20\}$, $4\leq M\leq 237$ and $65,000<|\mathcal{X}|<20^{237}$. We discuss the settings and properties of these datasets in greater detail in Appendix B.

For the biological sequence experiments we run a predetermined number of experimental rounds, $T=10$, we set the batch size to $B=128$, and use five different seeds for random initialization. For the first set of real-sequence experiments in Sec. 4.2, we use a fixed threshold, $\tau$, with the aim of finding $\mathbf{x}\in\mathcal{S}$ using both probabilistic and non probabilistic models. For the next set of experiments in Sec. 4.3, we set $\tau_{t}$ adaptively each round for testing VSD’s ability to find the maximizer, $\mathbf{x}^{*}$, of the black-box function, ${f\!\centerdot}$. We compare against DbAS (Brookes & Listgarten, ), CbAS (Brookes et al., ), AdaLead (Sinai et al., ), PEX (Ren et al., ), BORE (Tiao et al., ) adapted to use the score function gradient estimator, and a naïve baseline that uses random samples from the prior, $p({\mathbf{x}}|{\mathcal{D}_{0}})$. To reduce confounding, all methods share the same surrogate model and other components where possible (acquisition functions, priors and variational distributions). We optimize VSD, CbAS, DbAS and BORE for a minimum of 5000 iterations each round. When we use a CPE, AdaLead’s $\kappa$ parameter is set to 0.5 since the CPE already incorporates the appropriate threshold.

### 4.1 Conditional Generation of Handwritten Digits

Our motivating application for VSD is to model the space of fit DNA and protein sequences, which are string-representations of complex 3-dimensional structures. In this experiment we aim to demonstrate, by analogy, that VSD can operate with the types of generative models required for representing these complex structures.
For this task, we have chosen to ‘unroll’ (reverse the order of every odd row, and flatten) down-scaled ($14\times 14$ pixel, 8-bit) MNIST (LeCun et al., ) images into sequences, $\mathbf{x}$, where $M$ = 196 and $|\mathcal{V}|=8$. We then train long short-term memory (LSTM) recurrent neural network (RNN) and decoder-only causal transformer generative models on the entire MNIST training set by maximum likelihood (ML). These generative distributions are used as our prior models, $p({\mathbf{x}}|{\mathcal{D}_{0}})$, for VSD and we detail their form in Sec. B.1.

Figure: (a) LSTM Prior
Refer to caption: extracted/5897346/figs/cde/prior_digit_samples_LSTM.png

The task is then to use VSD in one round to estimate the posterior $p({\mathbf{x}}|{y\in\{3,5\}})$ using a CPE trained on labels $z_{n}=\mathbbm{1}[{y_{n}\in\{3,5\}}]$. We use a convolutional architecture for the CPE given in Sec. B.2, and it achieves a test balanced accuracy score of $\sim 99\%$. We parameterize the variational distributions, $q({\mathbf{x}}|{\phi})$, in the same way as the priors, and initialize these distribution parameters from the prior distribution parameters. During training with ELBO the prior distribution parameters are locked, and we run training for 5000 iterations. This is exactly lines 8 and 9 in Algorithm 1. Samples are visualized from the resulting variational distributions with the corresponding priors in Figure 2.
We can see that the prior LSTM and transformer are able to generate convincing digits once the sampled sequences are ‘re-rolled’, and that VSD is able to effectively refine these distributions, even though it does not have access to any data directly – only scores from the CPE. Both the LSTM and transformer yield qualitatively similar results, and have similar mean scores from the CPE.

### 4.2 Fitness Landscapes

In this setting we wish to find $\mathbf{x}\in\mathcal{S}$, so we fix $\tau$ over all rounds for all competing methods. We only consider the combinatorially (near) complete datasets to avoid any pathological behavior from relying on machine learning oracles (Surana et al., ). Results are presented in Figure 3. The primary measures by which we compare methods are precision, recall and performance (the last adapted from Jain et al. ()),

$$ $\displaystyle\textrm{Precision}_{t}$ $\displaystyle=\frac{1}{\min\{tB,|\mathcal{S}|\}}\sum^{t}_{r=1}\sum^{B}_{b=1} \mathbbm{1}[{y_{br}>\tau}]\cdot\mathbbm{1}[{\mathbf{x}_{br}\notin\mathcal{X}^{ q}_{b-1,r}}],$ (15) $\displaystyle\textrm{Recall}_{t}$ $\displaystyle=\frac{1}{\min\{TB,|\mathcal{S}|\}}\sum^{t}_{r=1}\sum^{B}_{b=1} \mathbbm{1}[{y_{br}>\tau}]\cdot\mathbbm{1}[{\mathbf{x}_{br}\notin\mathcal{X}^{ q}_{b-1,r}}],$ (16) $\displaystyle\textrm{Performance}_{t}$ $\displaystyle=\sum^{t}_{r=1}\sum^{B}_{b=1}y_{br}\cdot\mathbbm{1}[{\mathbf{x}_{ br}\notin\mathcal{X}^{q}_{b-1,r}}].$ (17) $$

Here $\mathcal{X}^{q}_{br}\subset\mathcal{X}$ is the set of experimentally queried sequences by the $b$th batch member of the $r$th round, including the initial training set.
These measures are comparable among probabilistic and non probabilistic methods. Precision and recall measure the ability of a method to efficiently explore $\mathcal{S}$, where $\min\{tB,|\mathcal{S}|\}$ is the size of the selected set at round $t$ (bounded by the number of good solutions), and $\min\{TB,|\mathcal{S}|\}$ is the number of positive elements possible in the experimental budget. Strictly, recall should be normalized by $|\mathcal{S}|$, but we use $TB$ here since it may not be realistic to have the experimental budget to fully explore $\mathcal{S}$. Performance measures the cumulative fitness of the unique batch members, and unlike Jain et al. () we do not normalize this measure.

For the DHFR and TrpB experiments we set maximum fitness in the training dataset to be that of the wild type, and $\tau$ to be slightly below the wild type fitness value (so we have $\sim\!10$ positive examples to train the CPE with). We use a randomly selected $N_{\textrm{train}}=2000$ below the wild-type fitness to initially train the CPE, we also explicitly include the wild-type. The thresholds and wild-type fitness values are; DHRF: $\tau=-0.1$, $y_{\textrm{wt}}=0$, TrpB: $\tau=0.35$, $y_{\textrm{wt}}=0.409$. We follow the same procedure for the TFBIND8 experiment, however, there is no notion of a wild-type sequence in this data, and so we set $\tau=0.75$, and $y_{\textrm{train max}}=0.85$.
We use a uniform prior over sequences, $p({\mathbf{x}})=\prod^{M}_{m=1}\mathrm{Categ}({x_{m}}|{\mathbf{1}\cdot|
\mathcal{V}|^{-1}})$, since these are relatively small search spaces, and the sub-sequences of nucleic/amino acids have been specifically selected for their task. Similarly, we find that relatively simple independent (mean-field) variational distributions of the form in Equation 20 and MLP based CPE s work best for these experiments (details in Sec. B.2).

Figure: (a) DHFR
Refer to caption: extracted/5897346/figs/fl/dhfr/Precision.png

Results are presented in Figure 3. VSD is the best performing method by most of the measures, with the related method CbAS also often performing well. We have found the AdaLead and PEX evolutionary-search based methods to be effective on lower-dimensional problems (TFBIND8 being the lowest here), however we consistently observe their performance degrading as the dimension of the problem increases. We suspect this is a direct consequence of their random mutation strategies being suited to exploration in low dimensions, but less efficient in higher dimensions compared to the learned generative models employed by VSD, CbAS, and DbAS. Our modified version of BORE (which is just the expected log-likelihood component of Equation 11) performs badly in all cases, and this is a direct consequence of its variational distribution collapsing to a point mass. In a non-batch setting this behavior is not problematic, but shows the importance of the KL divergence of VSD in this batch setting.
We replicate these experiments in Sec. C.1 using GP-PI, also backed by our guarantees. In all cases VSD’s results remain similar or improve slightly, whereas the other methods results remain similar or degrade. We also report on batch diversity scores in Sec. C.2.

### 4.3 Black-Box Optimization

In this experiment we aim to find the global maximizers, $(y^{*}\!,\mathbf{x}^{*})$, of the black-box function, ${f\!\centerdot}$. For this, we set $\tau_{t}$ adaptively by specifying it as an empirical quantile of the observed target values,

$$ $\displaystyle\tau_{t}=\tilde{Q}^{t}_{y}(\gamma_{t}{=}p_{t-1}^{\eta})$ (18) $$

where $\tilde{Q}^{t}_{y}$ is the empirical quantile function of targets at round $t$, $p_{t-1}$ is a percentile from the previous round, and $\eta\in[0,1]$ is a parameter that controls an annealing-like schedule for $\tau_{t}$ that prioritizes exploration of the fitness landscape in earlier rounds and exploitation of known fit regions in later rounds. This is a strategy loosely-similar to Srinivas et al. (). The measure of interest for these experiments is simple/instantaneous regret $\mathrm{r}_{t}$ which quantifies how close the methods get to finding the globally fittest sequence,

$$ $\displaystyle\mathrm{r}_{t}=y^{*}-\max_{y}\{y_{bi}\}^{B,t}_{b=1,i=1},$ (19) $$

where $y^{*}$ is the fitness value of the fittest sequence $\mathbf{x}^{*}$.
We use the higher dimensional AAV ($y^{*}{=}19.54$) and GFP ($y^{*}{=}4.12$) datasets to show that VSD can scale to higher dimensional problems. However, the $\mathcal{X}$ of these experiments is completely intractable to fully explore experimentally, and so we use a predictive oracle trained on all of the original experimental data as the ground-truth black-box function. This is the same strategy used in Brookes et al. (); Jain et al. (); Trabucco et al. (); Kirjner et al. () among others, and we use the exact CNN-based oracles from Kirjner et al. () for these experiments. However, we note here that some of the oracles used in these experiments do not predict well out-of-distribution (Surana et al., ), which limits their real-world applicability.

We follow Kirjner et al. () in the experimental settings for the AAV and GFP datasets, but we modify the maximum fitness training point and training dataset sizes to make them more amenable to a sequential optimization setting. The initial percentiles, schedule, and max training fitness values are; AAV: $p_{0}=0.8$, $\eta=0.7$, $y_{\textrm{max}}=5$, GFP: $p_{0}=0.8$, $\eta=0.7$ $y_{\textrm{max}}=1.9$. The edit distance between $\mathbf{x}^{*}$ and the fittest sequence in the CPE training data is 8 for GFP, and 13 for AAV. We again use a random $N_{\textrm{train}}=2000$ for training the CPE s, which in this case are CNNs – architecture specifics are in Sec. B.2.

Figure: (a) Independent
Refer to caption: extracted/5897346/figs/opt/gfp/Simple_regret_mc.png

In these higher dimensional settings, we find that performance of the methods heavily relies on using an informed prior (in the case of VSD and CbAS), or initial variational distribution (in the case of DbAS and BORE). To this end, we follow Brookes et al. () and fit the initial variational distribution to the CPE training sequences (regardless of fitness), but we use maximum likelihood. Then for VSD and CbAS we copy this distribution and fix its parameters for the remainder of the experiment for use as a prior. We also use this prior for the Random method, and AdaLead and PEX use alternative generative heuristics. For these experiments we use the simple independent variational distribution from before, and also the same LSTM and causal decoder-only transformer models from Sec. 4.1. We present additional results in Sec. C.3 with uninformed priors and transition variational distributions of the form $q({\mathbf{x}_{t}}|{\mathbf{x}_{t-1},\phi})$.

The results are summarized in Figure 4. VSD is among the leading methods for all experiments. VSD appears better take advantage of the more complex variational distributions than the competing methods CbAS and DbAS since it can sample from the adapted variational distribution while learning it. We can see that AdaLead, PEX and often BORE all perform worse than random for reasons previously mentioned. Simple regret can drop below zero for these experiments since an oracle is used as the black box function, but the global maximizer is taken from the experimental data. This potentially highlights some of the overconfidence issues when using oracles as outlined in Surana et al. (). Batch diversity scores for these experiments are presented in Sec. C.2.

## 5 Conclusion

We have presented variational search distributions (VSD), a method for efficiently finding designs of a rare class sequentially under some experimental constraints. VSD is underpinned by variational inference, which allows it to satisfy several critical requirements and important desiderata specific to this problem, including learning generative models for fit sequences and batch candidate generation. We show that VSD converges asymptotically to the true level-set distribution at the same rate as a Monte-Carlo estimator with full knowledge of the true distribution. We have also showcased the benefits of our method empirically on a set of combinatorially complete and high dimensional sequential-design biological problems and show that it can effectively learn powerful generative models of fit designs. Finally, our framework can be generalized to more complex application scenarios, potentially involving other challenging combinatorial optimization problems (Bengio et al., ), such as graph structures (Annadani et al., ), and mixed discrete-continuous variables, which are worth investigating as future work directions.

## Appendix A Acronyms

## Acronyms

## Appendix B Experimental Details

We use three well established datasets; a green fluorescent protein (GFP) from Aequorea Victoria (Sarkisyan et al., ), an adeno-associated virus (AAV) Bryant et al. (); and DNA binding activity to a human transcription factor (TFBIND8) (Trabucco et al., ; Barrera et al., ). These datasets have been used variously by Brookes & Listgarten (); Brookes et al. (); Angermueller et al. (); Kirjner et al. (); Jain et al. () among others. The GFP task is to maximize fluorescence, this protein consists of 238 amino acids, of which 237 can mutate. The AAV task us to maximize the genetic payload that can be delivered, and the associated protein has 28 amino acids, all of which can mutate. A complete combinatorial assessment is infeasible for these tasks, and so we use the convolution neural network oracle presented in Kirjner et al. () as *in-silico* ground truth. TFBIND8 contains a complete combinatorial assessment of the effect of changing 8 nucleotides on binding to human transcription factor SIX6 REF R1 (Barrera et al., ). The dataset we use contains all 65536 sequences prepared by Trabucco et al. (). We also use two novel datasets from recent works that experimentally assess the (near) complete combinatorial space of short sequences. The first dataset measures the antibiotic resistance of Escherichia coli metabolic gene folA, which encodes dihydrofolate reductase (DHFR) (Papkou et al., ). Only a sub-sequence of this gene is varied (9 nucleic acids which encode 3 amino acids), and so a near-complete (99.7%) combinatorial scan is available. For variants that have no fitness (resistance) data available, we give a score of $-1$. The next dataset is near-complete combinatorial scan of four interacting amino acid residues near the active site of the enzyme tryptophan synthase (TrpB) (Johnston et al., ), with 159,129 unique sequences and fitness values, we use $-0.2$ for the missing fitness values (we do not use the authors’ imputed values). These residues are explicitly shown to exhibit epistasis – or non-additive effects on catalytic function – which makes navigating this landscape a more interesting challenge from an optimization perspective. The properties of these datasets are presented in Table 3.

**Table 3: Alphabet size, sequence length, and number of availbale sequences for each of the datasets we use in this work.**
| Dataset | $|\mathcal{V}|$ | $M$ | $|\mathcal{X}_{\textrm{available}}|$ | $|\mathcal{X}|$ |
| --- | --- | --- | --- | --- |
| TFBIND8 | 4 | 8 | 65,536 | 65,536 |
| TrpB | 20 | 4 | 159,129 | 160,000 |
| DHFR | 4 | 9 | 261,333 | 262,144 |
| AAV | 20 | 28 | 42,340 | $20^{28}$ |
| GFP | 20 | 237 | 51,715 | $20^{237}$ |

### B.1 Variational distributions

In this section we summarize the main variational distribution architectures considered for VSD, BORE, CbAS and DbAS, and the sampling distributions for the Random baseline method. Somewhat surprisingly, we find that we obtain consistently good results for the biological sequence experiments using a simple independent (or mean-field) variational distribution,

$$ $\displaystyle q({\mathbf{x}}|{\phi})=\prod^{M}_{m=1}\mathrm{Categ}({x_{m}}|{ \mathrm{softmax}\!\left({\phi_{m}}\right)}),$ (20) $$

where $x_{m}\in\mathcal{V}$ and $\phi_{m}\in^{|\mathcal{V}|}$. However, this simple mean-field distribution was not capable of generating convincing handwritten digits. We have also tested a variety of transition variational distributions,

$$ $\displaystyle q({\mathbf{x}_{t}}|{\mathbf{x}_{t-1},\phi})=\prod^{M}_{m=1} \mathrm{Categ}({x_{tm}}|{\mathrm{softmax}\!\left({\textrm{NN}_{m}(\mathbf{x}_{ t-1},\phi)}\right)}),$ (21) $$

where $\textrm{NN}_{m}(\mathbf{x}_{t-1},\phi)$ is the $m^{\textrm{th}}$ vector output of a neural network that takes a sequence from the previous round, $\mathbf{x}_{t-1}$, as input. We have implemented multiple neural net encoder/decoder architectures for $\textrm{NN}_{m}(\mathbf{x}_{t-1},\phi)$, but we did not consider architectures of the form $\textrm{NN}_{m}(\phi)$ since the variational distribution in Equation 20 can always learn a $\phi_{m}=\textrm{NN}_{m}(\phi^{\prime})$. We found that none of these transition architectures significantly outperformed the mean-field distribution (Equation 20) when it was initialized well (e.g. fit to the CPE training sequences), see Sec. C.3 for results. We also implemented auto-regressive variational distributions of the form,

$$ $\displaystyle q({\mathbf{x}}|{\phi})$ $\displaystyle=\mathrm{Categ}({x_{1}}|{\mathrm{softmax}\!\left({\phi_{1}}\right )})\prod^{M}_{m=2}q({x_{m}}|{x_{1:m-1},\phi_{1:m}})\quad\textrm{where,}$ (22) $\displaystyle q({x_{m}}|{x_{1:m-1},\phi_{1:m}})$ $\displaystyle=\begin{cases}\mathrm{Categ}({x_{m}}|{\mathrm{softmax}\!\left({ \mathrm{LSTM}(x_{m-1},\phi_{m-1:m})}\right)}),\\ \mathrm{Categ}({x_{m}}|{\mathrm{softmax}\!\left({\mathrm{DTransformer}(x_{1:m- 1},\phi_{1:m})}\right)}).\end{cases}$ $$

For a LSTM RNN and a decoder-only transformer with a causal mask, for the latter see Phuong & Hutter (, Algorithm 10 & Algorithm 14) for maximum likelihood training and sampling implementation details respectively. For the digits experiment, the LSTM uses 5 layers of dimension 128, and the transformer uses 8-attention heads, with 4 layers and a feed-forward network size of 256. For the AAV and GFP experiments, the LSTM uses 4 layers of dimension 32, and the transformer uses 4 attention heads for GFP, 2 attention heads for AAV, 1 layer, and a feed-forward network size of 64. Any larger than this, and we found these models would over-fit. We use additive positional encoding for both of these models.

### B.2 Class Probability Estimator Architectures

For the fitness landscape experiments on the smaller combinatorially complete datasets we use a two-hidden layer MLP, with an input embedding layer. The architecture is given in Figure 5 (a). For the larger dimensional AAV and GFP datasets, we use the convolutional architecture given in Figure 5 (b). Five fold cross validation is used to select the hyper parameters before the CPE s are trained on the whole training set for use in the subsequent experimental rounds. Model updates are performed by retraining on the whole query set.

Figure: Figure 5: CPE architectures used for the experiments in PyTorch syntax. $\texttt{A}=|\mathcal{V}|$, $\texttt{M}=M$, GFP uses a max pooling kernel size and stride of 4, all other datasets use 2.

## Appendix C Additional Experimental Results

### C.1 Fitness Landscapes – Gaussian Process Probability of Improvement

Here we present additional fitness landscape experimental results, where we have used a GP as a surrogate model for $p({y}|{\mathbf{x},\mathcal{D}_{N}})$ in conjunction with a complementary Normal CDF as the PI acquisition function. This is the one of the main frameworks supported by our theoretical analysis.
VSD, DbAS, CbAS and BORE make use of the GP-PI acquisition function, and so BORE is BOPR in this instance since we are not using a CPE. PEX and AdaLead only use the GP surrogate, as per their original formulation.
The GP uses a simple categorical kernel with automatic relevance determination from Balandat et al. (),

$$ $\displaystyle k(\mathbf{x},\mathbf{x}^{\prime})=\sigma\exp\!\left(-\frac{1}{M} \sum^{M}_{m=1}\frac{\mathbbm{1}[{x_{m}=x^{\prime}_{m}}]}{l_{m}}\right),$ (23) $$

where $\sigma$ and $l_{m}$ are hyper-parameters controlling scale and length-scale respectively. See Figure 6 for the results.

Figure: (a) DHFR
Refer to caption: extracted/5897346/figs/fl/dhfr_gp/Precision.png

### C.2 Diversity Scores

The diversity of batches of candidates is a common thing to report in the literature, and to that end we present the diversity of our results here. We have taken the definition of diversity from (Jain et al., ) as,

$$ $\displaystyle\mathrm{Diversity}_{t}=\frac{1}{B(B-1)}\sum_{\mathbf{x}_{i}\in \mathcal{D}_{Bt}}\sum_{\mathbf{x}_{j}\in\mathcal{D}_{Bt}\setminus\{\mathbf{x}_ {i}\}}\mathrm{Lev}(\mathbf{x}_{i},\mathbf{x}_{j}),$ (24) $$

where $\mathrm{Lev}:\mathcal{X}\times\mathcal{X}\to\mathbb{N}_{0}$ is the Levenshtein distance. We caution the reader as to the interpretation of these results however, as more diverse batches often do not lead to better performance, precision, recall or simple regret (as can be seen from the Random method results). Though insufficient diversity can also explain poor performance, as in the case of BORE. Results for the fitness landscape experiment are presented in Figure 7, and black-box optimization in Figure 8.

Figure: (a) DHFR
Refer to caption: extracted/5897346/figs/fl/dhfr/Diversity.png

Figure: (a) Independent
Refer to caption: extracted/5897346/figs/opt/gfp/Diversity_mc.png

### C.3 Ablations – Variational and Prior Distributions

In Figure 9 we present ablation results for VSD using different priors and variational distributions. We use the BBO experimental datasets for this task as they are higher-dimensional and so more sensitive to these design choices. We test the following prior and variational posterior distributions:

Independent categorical variational posterior distribution of the form in Equation 20, and a uniform prior distribution, $p({\mathbf{x}})=\prod^{M}_{m=1}\mathrm{Categ}({x_{m}}|{\mathbf{1}\cdot|
\mathcal{V}|^{-1}})$.

Independent categorical prior and variational posterior of the form in Equation 20. The prior is fit using ML on the initial CPE training data.

LSTM prior and variational posterior of the form Equation 22. The prior is fit using ML on the initial CPE training data.

Decoder-only causal transformer prior and variational posterior of the form Equation 22. The prior is fit using ML on the initial CPE training data.

Independent categorical prior and a transition-style auto-encoder variational posterior of the form Equation 21, where we use two-hidden layer MLPs for the encoder and decoder. The prior is fit using ML on the initial CPE training data.

Independent categorical prior and a transition-style convolutional auto-encoder variational posterior of the form Equation 21, where we use a convolutional encoder, and transpose convolutional decoder. The prior is fit using ML on the initial CPE training data.

We use the informed-independent priors with the transition variational distributions since they are somewhat counter-intuitive to use as priors themselves.

Figure: (a) AAV
Refer to caption: extracted/5897346/figs/ablations/Simple_regret_aav.png

From Figure 9 we can see that while using an uninformative prior works in the lower-dimensional fitness landscape experiments, using an informative prior is crucial for these higher dimensional problems. We found a similar result when using this uninformative prior with CbAS, or using a uniform initialization with DbAS and BORE. The methods are not able to make any significant progress within the experimental budget given. The independent and transition variational distributions achieve similar performance, whereas the auto-regressive models generally outperform all others. This is because of the LSTM and transformer’s superior generalization performance when generating sequences – measured both when training the priors (on held-out sequences) and during VSD adaptation.

## Appendix D Theoretical analysis for GP-based CPEs

In this section, we present theoretical results concerning VSD and its estimates when equipped with Gaussian process regression models (Rasmussen & Williams, ). We show that VSD sampling distributions converge to a target distribution that characterizes the level set given by $\tau$. The approximation error mainly depends on the predictive uncertainty of the probabilistic model with respect to the true underlying function ${f\!\centerdot}$. For the analysis, we will assume that ${f\!\centerdot}$ is drawn from a Gaussian process, i.e., ${f\!\centerdot}\sim\mathcal{GP}(0,k)$, with a positive-semidefinite covariance (or kernel) function $k:\mathcal{X}\times\mathcal{X}\to$. In this case, we can show that the predictive uncertainty of the model converges (in probability) to zero as the number of observations grows. From this result, we prove asymptotic convergence guarantees for VSD equipped with GP-PI-based CPEs. These results form the basis for our analysis of CPEs based on neural networks (Appendix E).

### D.1 Gaussian process posterior

Let ${f\!\centerdot}\sim\mathcal{GP}(0,k)$ be a zero-mean Gaussian process with a positive-semidefinite covariance function $k:\mathcal{X}\times\mathcal{X}\to$. Assume that we are given a set $\mathcal{D}_{N}:=\{(\mathbf{x}_{i},y_{i})\}_{i=1}^{N}$ of $N\geq 1$ observations $y_{i}={f\!\centerdot}(\mathbf{x}_{i})+\epsilon_{i}$, where $\epsilon\sim\mathcal{N}\!\left({0,\sigma_{\epsilon}^{2}}\right)$ and $\mathbf{x}_{i}\in\mathcal{X}$. The GP posterior predictive distribution at any $\mathbf{x}\in\mathcal{X}$ is then given by (Rasmussen & Williams, ):

$$ $\displaystyle{f\!\centerdot}(\mathbf{x})|\mathcal{D}_{N}$ $\displaystyle\sim\mathcal{N}\!\left({\mu_{N}(\mathbf{x}),\sigma_{N}^{2}( \mathbf{x})}\right)$ (25) $\displaystyle\mu_{N}(\mathbf{x})$ $\displaystyle={\boldsymbol{\mathbf{k}}}_{N}(\mathbf{x})^{\top}({\boldsymbol{ \mathbf{{K}}}}_{N}+\sigma_{\epsilon}^{2}\mathbf{I})^{-1}{\boldsymbol{\mathbf{y }}}_{N}$ (26) $\displaystyle k_{N}(\mathbf{x},\mathbf{x}^{\prime})$ $\displaystyle=k(\mathbf{x},\mathbf{x}^{\prime})-{\boldsymbol{\mathbf{k}}}_{N}( \mathbf{x})^{\top}({\boldsymbol{\mathbf{{K}}}}_{N}+\sigma_{\epsilon}^{2} \mathbf{I})^{-1}{\boldsymbol{\mathbf{k}}}_{N}(\mathbf{x}^{\prime})$ (27) $\displaystyle\sigma_{N}^{2}(\mathbf{x})$ $\displaystyle=k_{N}(\mathbf{x},\mathbf{x})\,,$ (28) $$

where ${\boldsymbol{\mathbf{k}}}_{N}(\mathbf{x}):=[k(\mathbf{x},\mathbf{x}_{i})]_{i=1
}^{N}\in^{N}$, ${\boldsymbol{\mathbf{{K}}}}_{N}:=[k(\mathbf{x}_{i},\mathbf{x}_{j})]_{i,j=1}^{N
,N}\in^{N\times N}$, and ${\boldsymbol{\mathbf{y}}}_{N}:=[y_{i}]_{i=1}^{N}\in^{N}$.

#### Batch size.

In the following, we will assume a batch of size $B=1$ to keep the proofs simple. With this assumption, at every iteration $t\geq 1$, we have $N=t$ observations available in the dataset. We would, however, like to emphasize that sampling a batch of multiple observations, instead of a single observation, per iteration should only improve the convergence rates by a constant (batch-size-dependent) multiplicative factor. Therefore, our results remain valid as an upper bound for the convergence rates of VSD in the batch setting.

### D.2 Background

We will consider an underlying probability space $(\Omega,{\mathfrak{{A}}},\mathbb{P})$, where $\Omega$ is the sample space, ${\mathfrak{{A}}}$ denotes the $\sigma$-algebra of events, and $\mathbb{P}$ is a probability measure. For any event ${\mathcal{{A}}}\in{\mathfrak{{A}}}$, we have that $\mathbb{P}\!\left[{\mathcal{{A}}}\right]\in[0,1]$ quantifies the probability of that event. For events involving a random variable, e.g., $\chi:(\Omega,{\mathfrak{{A}}})\to(,{\mathfrak{{B}}}_{)}$, where ${\mathfrak{{B}}}$ denotes the Borel $\sigma$-algebra of the real line with its usual topology, we will let:

$$ $\mathbb{P}\!\left[\chi>0\right]=\mathbb{P}\!\left[\{\omega\in\Omega:\chi( \omega)>0\}\right]\,.$ (29) $$

We will also use conditional expectations, i.e., given a $\sigma$-sub-algebra ${\mathfrak{{S}}}$ of ${\mathfrak{{A}}}$, the conditional expectation $\mathbb{E}[\chi|{\mathfrak{{S}}}]$ is a ${\mathfrak{{S}}}$-measurable random variable such that:

$$ $\forall{\mathcal{{A}}}\in{\mathfrak{{S}}}\,,\quad\int_{\mathcal{{A}}}\mathbb{E }[\chi|{\mathfrak{{S}}}]{\mathop{}\operatorname{d}}\mathbb{P}=\int_{\mathcal{{ A}}}\chi{\mathop{}\operatorname{d}}\mathbb{P}=\mathbb{E}[\chi|{\mathcal{{A}}}]\,.$ (30) $$

We will denote by $\{{\mathfrak{{F}}}_{t}\}_{t=0}^{\infty}$ an increasing filtration on ${\mathfrak{{A}}}$. For instance, we could set ${\mathfrak{{F}}}_{t}$ as the $\sigma$-algebra generated by the random variables in the algorithm (i.e., the candidates, target observations, etc.) at time $t$. For more details on the measure-theoretic definition of probability, we refer the reader to classic textbooks in the area (e.g. Bauer, ; Durrett, )

We will use the following well known notation for asymptotic convergence results. For a given strictly positive function $g:\mathbb{N}\to$, we define ${\mathcal{{O}}}(g(t))$ as the set of functions asymptotically bounded by $g$ (up to a constant factor) as:

$$ ${\mathcal{{O}}}(g(t)):=\left\{h:\mathbb{N}\to\>\middle|\>\limsup_{t\to\infty} \frac{|h(t)|}{g(t)}<\infty\right\},$ (31) $$

and for convergence in probability we use its stochastic counterpart:

$$ ${\mathcal{{O}}}_{\mathbb{P}}(g(t)):=\left\{\rho:\mathbb{N}\times(\Omega,{ \mathfrak{{A}}})\to(,{\mathfrak{{B}}}_{)}\>\middle|\>\lim_{C\to\infty}\limsup_ {t\to\infty}\mathbb{P}\!\left[\frac{|\rho(t)|}{g(t)}>C\right]=0\right\}\,,$ (32) $$

which is equivalent to:

$$ $\forall\varepsilon>0,\quad\exists C_{\varepsilon}\in(0,\infty):\quad\mathbb{P} \!\left[|\rho_{t}|>C_{\varepsilon}\right]\leq\varepsilon,\quad\forall t\geq T_ {\varepsilon}\,,$ (33) $$

for some $T_{\varepsilon}\in\mathbb{N}$. For almost sure convergence, we may also say that a sequence of random variables $\rho_{t}$, $t\in\mathbb{N}$, is almost surely ${\mathcal{{O}}}(g(t))$ if $\mathbb{P}\!\left[\rho_{t}\in{\mathcal{{O}}}(g(t))\right]=1$. A deeper overview on these notations and their properties can be found in García-Portugués ().

### D.3 Auxiliary results

We start with a few technical results which will form the basis for our derivations. The following recursive relations allow us to derive convergence rates for the variance of a GP posterior by analyzing how much it reduces per iteration.

###### Lemma D.1 ( Chowdhury & Gopalan ( 2017 , Appendix F) ) .

The posterior mean and covariance functions of a Gaussian process given $t\geq 1$ observations obey the following recursive identities:

$$ $\displaystyle\mu_{t}(\mathbf{x})$ $\displaystyle=\mu_{t-1}(\mathbf{x})+\frac{k(\mathbf{x},\mathbf{x}_{t})}{\sigma _{\epsilon}^{2}+\sigma_{t-1}^{2}(\mathbf{x}_{t})}(y_{t}-\mu_{t-1}(\mathbf{x}))$ (34) $\displaystyle k_{t}(\mathbf{x},\mathbf{x}^{\prime})$ $\displaystyle=k_{t-1}(\mathbf{x},\mathbf{x}^{\prime})-\frac{k_{t-1}(\mathbf{x} ,\mathbf{x}_{t})k_{t-1}(\mathbf{x}_{t},\mathbf{x}^{\prime})}{\sigma_{\epsilon} ^{2}+\sigma_{t-1}^{2}(\mathbf{x}_{t})}$ (35) $\displaystyle\sigma_{t}^{2}(\mathbf{x})$ $\displaystyle=\sigma_{t-1}^{2}(\mathbf{x})-\frac{k_{t-1}^{2}(\mathbf{x}, \mathbf{x}_{t})}{\sigma_{\epsilon}^{2}+\sigma_{t-1}^{2}(\mathbf{x}_{t})}\,,$ (36) $$

for $\mathbf{x},\mathbf{x}^{\prime}\in\mathcal{X}$.

We will also make use of the following version of the second Borel-Cantelli lemma adapted from Durrett (, Thr. 4.5.5) and its original statement in Dubins & Freedman ().

###### Lemma D.2 (Second Borel-Cantelli lemma) .

Let $\{{\mathcal{{A}}}_{t}\}_{t=1}^{\infty}$ be a sequence of events where ${\mathcal{{A}}}_{t}\in{\mathfrak{{F}}}_{t}$, for all $t\in\mathbb{N}$, and let $\chi_{t}:\omega\mapsto\mathbbm{1}[{\omega\in{\mathcal{{A}}}_{t}}]$, for $\omega\in\Omega$. Then the following holds with probability 1:

$$ $\lim_{T\to\infty}\frac{\sum_{t=1}^{T}\chi_{t}}{\sum_{t=1}^{T}\mathbb{P}\!\left [{\mathcal{{A}}}_{t}|{\mathfrak{{F}}}_{t-1}\right]}=L<\infty\,,$ (37) $$

assuming $\mathbb{P}\!\left[{\mathcal{{A}}}_{1}|{\mathfrak{{F}}}_{0}\right]>0$. In addition, if $\lim_{T\to\infty}\sum_{t=1}^{T}\mathbb{P}\!\left[{\mathcal{{A}}}_{t}|{
\mathfrak{{F}}}_{t-1}\right]=\infty$, then $L=1$.

The next result provides us with an upper bound on the posterior variance of a Gaussian process which is valid for any covariance function.

###### Lemma D.3 .

Let $k:\mathcal{X}\times\mathcal{X}\to$ be any positive-semidefinite kernel on $\mathcal{X}$, and let $\tilde{k}:\mathcal{X}\times\mathcal{X}\to$ be a kernel defined as:

$$ $\tilde{k}(\mathbf{x},\mathbf{x}^{\prime})=\begin{cases}k(\mathbf{x},\mathbf{x} ),\quad&\mathbf{x}=\mathbf{x}^{\prime}\\ 0,\quad&\mathbf{x}\neq\mathbf{x}^{\prime},\end{cases}$ (38) $$

for $\mathbf{x},\mathbf{x}^{\prime}\in\mathcal{X}$. Given any set of observations $\{\mathbf{x}_{i},y_{i}\}_{i=1}^{t}$, for $t\geq 1$, denote by $\sigma_{t}^{2}$ the predictive variance of a GP model with prior covariance given by $k$, and let $\tilde{\sigma}_{t}^{2}$ denote the predictive variance of a GP model configured with $\tilde{k}$ as prior covariance function, where both models are given the same set of observations. Then the following holds for all $t\geq 0$:

$$ $\sigma_{t}^{2}(\mathbf{x})\leq\tilde{\sigma}_{t}^{2}(\mathbf{x})=\frac{\sigma_ {\epsilon}^{2}\tilde{\sigma}_{0}^{2}(\mathbf{x})}{\sigma_{\epsilon}^{2}+N_{t}( \mathbf{x})\tilde{\sigma}_{0}^{2}(\mathbf{x})},\quad\forall\mathbf{x}\in \mathcal{X}\,,$ (39) $$

where $N_{t}(\mathbf{x})$ denotes the number of observations at $\mathbf{x}$, and $\tilde{\sigma}_{0}^{2}(\mathbf{x})=\sigma_{0}^{2}(\mathbf{x}):=k(\mathbf{x},
\mathbf{x})$, for $\mathbf{x}\in\mathcal{X}$.

###### Proof.

It is not hard to show that $\tilde{k}$ defines a valid positive-semidefinite covariance function whenever $k$ is positive semidefinite. We will then focus on proving the main statement by an induction argument. The proof that the statement holds for the base case at $t=0$ is trivial given the definition:

$$ $\sigma_{0}^{2}(\mathbf{x})=k(\mathbf{x},\mathbf{x})=\tilde{k}(\mathbf{x}, \mathbf{x})=\tilde{\sigma}_{0}^{2}(\mathbf{x}),\quad\forall\mathbf{x}\in \mathcal{X}\,.$ (40) $$

Now assume that, for a given $t>0$, it holds that $\sigma_{t}^{2}(\mathbf{x})\leq\tilde{\sigma}_{t}^{2}(\mathbf{x})$, for all $\mathbf{x}\in\mathcal{X}$. We will then check if the inequality remains valid at $t+1$. By D.1, we have that:

$$ $\begin{split}\sigma_{t+1}^{2}(\mathbf{x})&=\sigma_{t}^{2}(\mathbf{x})-\frac{k_ {t}^{2}(\mathbf{x},\mathbf{x}_{t+1})}{\sigma_{t}^{2}(\mathbf{x}_{t+1})+\sigma_ {\epsilon}^{2}}\\ \end{split}$ (41) $$

For any $\mathbf{x}\in\mathcal{X}$ such that $\mathbf{x}\neq\mathbf{x}_{t+1}$, we know that $\tilde{k}_{t}(\mathbf{x},\mathbf{x}_{t+1})\geq 0$, so that (again by D.1):

$$ $\tilde{k}_{t}^{2}(\mathbf{x},\mathbf{x}_{t+1})\leq\tilde{k}^{2}(\mathbf{x}, \mathbf{x}_{t+1})=0\,,$ (42) $$

which shows that:

$$ $\forall\mathbf{x}\neq\mathbf{x}_{t+1},\quad\sigma_{t+1}^{2}(\mathbf{x})\leq \sigma_{t}^{2}(\mathbf{x})\leq\tilde{\sigma}_{t}^{2}(\mathbf{x})=\tilde{\sigma }_{t+1}^{2}(\mathbf{x})\,.$ (43) $$

At $\mathbf{x}=\mathbf{x}_{t+1}$, we can rewrite $\sigma_{t+1}^{2}(\mathbf{x})=\sigma_{t+1}^{2}(\mathbf{x}_{t+1})$ as:

$$ $\begin{split}\sigma_{t+1}^{2}(\mathbf{x}_{t+1})=\frac{\sigma_{\epsilon}^{2} \sigma_{t}^{2}(\mathbf{x}_{t+1})}{\sigma_{t}^{2}(\mathbf{x}_{t+1})+\sigma_{ \epsilon}^{2}}\,.\end{split}$ (44) $$

We then check the difference:

$$ $\begin{split}\sigma_{t+1}^{2}(\mathbf{x}_{t+1})-\tilde{\sigma}_{t+1}^{2}( \mathbf{x}_{t+1})&=\frac{\sigma_{\epsilon}^{2}\sigma_{t}^{2}(\mathbf{x}_{t+1}) }{\sigma_{t}^{2}(\mathbf{x}_{t+1})+\sigma_{\epsilon}^{2}}-\frac{\sigma_{ \epsilon}^{2}\tilde{\sigma}_{t}^{2}(\mathbf{x}_{t+1})}{\tilde{\sigma}_{t}^{2}( \mathbf{x}_{t+1})+\sigma_{\epsilon}^{2}}\\ &=\frac{\sigma_{\epsilon}^{2}\sigma_{t}^{2}(\mathbf{x}_{t+1})(\tilde{\sigma}_{ t}^{2}(\mathbf{x}_{t+1})+\sigma_{\epsilon}^{2})-\sigma_{\epsilon}^{2}\tilde{ \sigma}_{t}^{2}(\mathbf{x}_{t+1})(\sigma_{t}^{2}(\mathbf{x}_{t+1})+\sigma_{ \epsilon}^{2})}{(\sigma_{t}^{2}(\mathbf{x}_{t+1})+\sigma_{\epsilon}^{2})( \tilde{\sigma}_{t}^{2}(\mathbf{x}_{t+1})+\sigma_{\epsilon}^{2})}\\ &=\frac{\sigma_{\epsilon}^{4}(\sigma_{t}^{2}(\mathbf{x}_{t+1})-\tilde{\sigma}_ {t}^{2}(\mathbf{x}_{t+1}))}{(\sigma_{t}^{2}(\mathbf{x}_{t+1})+\sigma_{\epsilon }^{2})(\tilde{\sigma}_{t}^{2}(\mathbf{x}_{t+1})+\sigma_{\epsilon}^{2})}\\ &\leq 0\,,\end{split}$ (45) $$

since $\sigma_{t}^{2}(\mathbf{x}_{t+1})\leq\tilde{\sigma}_{t}^{2}(\mathbf{x}_{t+1})$ by our assumption for time $t$. Therefore, we have shown that:

$$ $\sigma_{t}^{2}(\mathbf{x})\leq\tilde{\sigma}_{t}^{2}(\mathbf{x})\implies\sigma _{t+1}^{2}(\mathbf{x})\leq\tilde{\sigma}_{t+1}^{2}(\mathbf{x})\,,\quad\forall \mathbf{x}\in\mathcal{X}\,.$ (46) $$

From the conclusion above and the base case, the inequality in the main result follows by induction.

Now we derive an explicit form for $\tilde{\sigma}_{t}^{2}$. Note that this case corresponds to an independent Gaussian model, i.e., ${f\!\centerdot}(\mathbf{x})\perp\!\!\!\perp{f\!\centerdot}(\mathbf{x}^{\prime})$ whenever $\mathbf{x}\neq\mathbf{x}^{\prime}$, for ${f\!\centerdot}\sim\mathcal{GP}(0,\tilde{k})$. For any $t\geq 1$, this model's predictive variance at any $\mathbf{x}\in\mathcal{X}$ is given by:

$$ $\tilde{\sigma}_{t}^{2}(\mathbf{x})=\begin{dcases}\tilde{\sigma}_{t-1}^{2}( \mathbf{x}),&\mathbf{x}\neq\mathbf{x}_{t}\\ \frac{\sigma_{\epsilon}^{2}\tilde{\sigma}_{t-1}^{2}(\mathbf{x}_{t})}{\sigma_{ \epsilon}^{2}+\tilde{\sigma}_{t-1}^{2}(\mathbf{x}_{t})}=\left(\frac{1}{\tilde{ \sigma}_{t-1}^{2}(\mathbf{x}_{t})}+\frac{1}{\sigma_{\epsilon}^{2}}\right)^{-1} ,&\mathbf{x}=\mathbf{x}_{t}\end{dcases}$ (47) $$

Looking at the reciprocal, we have that:

$$ $\forall t\geq 1,\quad\frac{1}{\tilde{\sigma}_{t}^{2}(\mathbf{x})}=\frac{1}{ \tilde{\sigma}_{t-1}^{2}(\mathbf{x}_{t})}+\frac{\mathbbm{1}[{\mathbf{x}_{t}= \mathbf{x}}]}{\sigma_{\epsilon}^{2}},\quad\forall\mathbf{x}\in\mathcal{X}.$ (48) $$

Therefore, every observation at $\mathbf{x}$ is simply adding a factor of $\sigma_{\epsilon}^{-2}$ to $\tilde{\sigma}_{t}^{-2}(\mathbf{x})$. Unwrapping this recursion leads us to:

$$ $\forall t\geq 1,\quad\frac{1}{\tilde{\sigma}_{t}^{2}(\mathbf{x})}=\frac{1}{ \tilde{\sigma}_{0}^{2}(\mathbf{x})}+\frac{1}{\sigma_{\epsilon}^{2}}\sum_{i=1}^ {t}\mathbbm{1}[{\mathbf{x}_{i}=\mathbf{x}}]\,,\quad\forall\mathbf{x}\in \mathcal{X}\,.$ (49) $$

The result in D.3 then follows as the reciprocal of the above, which concludes the proof.
∎

###### Lemma D.4 .

Let ${f\!\centerdot}\sim\mathcal{GP}(0,k)$ for a given $k:\mathcal{X}\times\mathcal{X}\to$, where $\sigma_{\mathcal{X}}^{2}:=\sup_{\mathbf{x}\in\mathcal{X}}k(\mathbf{x},\mathbf{
x})<\infty$, and $|\mathcal{X}|<\infty$. Then ${f\!\centerdot}$ is almost surely bounded, and:

$$ $\mathbb{E}\!\left[\sup_{\mathbf{x}\in\mathcal{X}}\lvert{f\!\centerdot}(\mathbf {x})\rvert\right]\leq\sigma_{\mathcal{X}}\sqrt{2\log|\mathcal{X}|}\,.$ (50) $$

###### Proof.

The result follows by an application of a concentration inequality for the maximum of a finite collection of sub-Gaussian random variables (Boucheron et al., , Sec. 2.5). Note that $\{{f\!\centerdot}(\mathbf{x})\}_{\mathbf{x}\in\mathcal{X}}$ is a collection of $|\mathcal{X}|$ Gaussian, and therefore sub-Gaussian, random variables with sub-Gaussian parameter given by $\sigma_{\mathcal{X}}^{2}\geq\sigma_{t}^{2}(\mathbf{x})$, for all $\mathcal{X}$. Applying the maximal inequality for a finite collection sub-Gaussian random variables (Boucheron et al., , Thr. 2.5), we have that:

$$ $\mathbb{E}\!\left[\max_{\mathbf{x}\in\mathcal{X}}{f\!\centerdot}(\mathbf{x}) \right]\leq\sigma_{\mathcal{X}}\sqrt{2\log|\mathcal{X}|}<\infty\,.$ (51) $$

By symmetry, we know that $-{f\!\centerdot}(\mathbf{x})$ is also sub-Gaussian with the same parameter, so that the bound remains valid for $\max_{\mathbf{x}\in\mathcal{X}}-{f\!\centerdot}(\mathbf{x})$. As a consequence, the expected value of the maximum of $\rvert{f\!\centerdot}(\mathbf{x})\rvert$ is upper bounded by the same constant. On a finite set, the maximum and the supremum coincide. As the expected value of the supremum is finite, the supremum must be almost surely finite by Markov's inequality, and therefore ${f\!\centerdot}$ is almost surely bounded.
∎

###### Lemma D.1 ( Chowdhury & Gopalan ( 2017 , Appendix F) ) .

###### Lemma D.2 (Second Borel-Cantelli lemma) .

###### Lemma D.3 .

###### Proof.

###### Lemma D.4 .

###### Proof.

### D.4 Asymptotic convergence

The main assumption we will be working with in this section is the following.

###### Assumption D.1 .

The objective function is a sample from a Gaussian process ${f\!\centerdot}\sim\mathcal{GP}(0,k)$, where $k:\mathcal{X}\times\mathcal{X}\to$ is a bounded positive-semidefinite kernel on $\mathcal{X}$.

The next result allows us to derive a convergence rate for the posterior variance of a GP as a function of the sampling probabilities. This result might also be useful by itself for other sampling problems involving GP-based approximations.

###### Lemma D.5 .

Let $\{\mathbf{x}_{t}\}_{t\geq 1}$ be a sequence of $\mathcal{X}$-valued random variables adapted to the filtration $\{{\mathfrak{{F}}}_{t}\}_{t\geq 1}$. For a given $\mathbf{x}\in\mathcal{X}$, assume that the following holds:

$$ $\exists T_{*}\in\mathbb{N}:\quad\forall T\geq T_{*},\quad\sum_{t=1}^{T}\mathbb {P}\!\left[\mathbf{x}_{t}=\mathbf{x}\mid{\mathfrak{{F}}}_{t-1}\right]\geq B_{T }>0\,,$ (52) $$

for a some sequence of lower bounds $\{B_{t}\}_{t\in\mathbb{N}}$. Then, under D.1, given observations at $\{\mathbf{x}_{i}\}_{i=1}^{t}$, the following holds with probability 1:

$$ $\sigma_{t}^{2}(\mathbf{x})\in{\mathcal{{O}}}(B_{t}^{-1}).$ (53) $$

In addition, if $B_{t}\to\infty$, then $\lim_{t\to\infty}B_{t}\sigma_{t}^{2}(\mathbf{x})\leq\sigma_{\epsilon}^{2}$.

###### Proof.

At any iteration $t$, the posterior variance $\sigma_{t}^{2}$ of a GP model is upper bounded by a worst case assumption of no correlation between observations (see D.3). In this case, we have that:

$$ $\sigma_{t}^{2}(\mathbf{x})\leq\tilde{\sigma}_{t}^{2}(\mathbf{x})=\frac{\sigma_ {\epsilon}^{2}\tilde{\sigma}_{0}^{2}(\mathbf{x})}{\sigma_{\epsilon}^{2}+N_{t} \tilde{\sigma}_{0}^{2}(\mathbf{x})}\,,$ (54) $$

where $\tilde{\sigma}_{0}^{2}(\mathbf{x}):=\tilde{k}(\mathbf{x},\mathbf{x})=k(\mathbf
{x},\mathbf{x})$, and $N_{t}:=N_{t}(\mathbf{x})\leq t$ denotes the total number of observations taken at $\mathbf{x}$ as of iteration $t$. Without loss of generality, assume that $\tilde{\sigma}_{0}^{2}(\mathbf{x})=1$.

The only random variable to be bounded in Equation 54 is $N_{t}$. Let $\chi_{t}:=\mathbbm{1}[{\mathbf{x}_{t}=\mathbf{x}}]$, so that:

$$ $N_{t}=\sum_{i=1}^{t}\chi_{i}=\sum_{i=1}^{t}\mathbbm{1}[{\mathbf{x}_{t}=\mathbf {x}}]\,,\quad t\geq 1.$ (55) $$

We now apply the second Borel-Cantelli lemma (D.2) to $N_{t}$. Namely, let $\widehat{N}_{t}$ denote the sum of conditional expectations of $\{\chi_{i}\}_{i=1}^{t}$ given available data, i.e.:

$$ $\widehat{N}_{t}:=\sum_{i=1}^{t}\mathbb{E}[\chi_{i}\mid{\mathfrak{{F}}}_{i-1}]= \sum_{i=1}^{t}\mathbb{E}[\mathbbm{1}[{\mathbf{x}_{t}=\mathbf{x}}]\mid{ \mathfrak{{F}}}_{i-1}]=\sum_{i=1}^{t}\mathbb{P}\!\left[\mathbf{x}_{i}=\mathbf{ x}\mid{\mathfrak{{F}}}_{i-1}\right]\,.$ (56) $$

By D.2, we know that the following holds for some $L\in$:

$$ $\lim_{t\to\infty}\frac{N_{t}}{\widehat{N}_{t}}=L<\infty\,.$ (57) $$

Hence, $N_{t}$ is asymptotically equivalent to $\widehat{N}_{t}$. Applying this fact to $\tilde{\sigma}_{t}^{2}$, we have that:

$$ $\begin{split}\lim_{t\to\infty}B_{t}\tilde{\sigma}_{t}^{2}(\mathbf{x})&=\lim_{t \to\infty}\frac{B_{t}\sigma_{\epsilon}^{2}}{\sigma_{\epsilon}^{2}+N_{t}}\\ &=\lim_{t\to\infty}\frac{B_{t}\sigma_{\epsilon}^{2}}{\sigma_{\epsilon}^{2}+L \widehat{N}_{t}}\\ &\leq\lim_{t\to\infty}\frac{B_{t}\sigma_{\epsilon}^{2}}{\sigma_{\epsilon}^{2}+ LB_{t}}\\ &\leq\frac{1}{L}\lim_{t\to\infty}\min\{LB_{t},\sigma_{\epsilon}^{2}\}\\ &<\infty\,,\end{split}$ (58) $$

which holds with probability 1. Lastly, note that, if $B_{t}\to\infty$, then $L=1$ by D.2, and the last limit above becomes $\sigma_{\epsilon}^{2}$. The main result then follows by an application of D.3 and the definition of the big-${\mathcal{{O}}}$ notation (see Equation 31).(^1^11Recall that for convergent sequences $\lim$ and $\limsup$ coincide.)
∎

We assume a finite search space, which is the case for spaces of discrete sequences of bounded length. However, we conjecture that our results can be extended to continuous or mixed discrete-continuous search spaces via a discretization argument under further assumptions on the kernel $k$ (e.g., ensuring that ${f\!\centerdot}$ is Lipschitz continuous, as in Srinivas et al. ()).

###### Assumption D.2 .

The search space $\mathcal{X}$ is finite, $|\mathcal{X}|<\infty$.

We assume that our family of variational distributions is rich enough to be able to represent the PI-based distribution $p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})$, which is the optimum of our variational objective when the optimal classifier is given by GP-PI. Although this assumption could be seen as strong, note that, due to Gaussian noise, the classification probability $p({y>\tau_{t}}|{\mathbf{x},\mathcal{D}_{t}})$ should be a reasonably smooth function of $\mathbf{x}$, which facilitates the approximation of the resulting posterior by a generative model.

###### Assumption D.3 .

For every $t\geq 0$, $p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})$ is a member of the variational family, i.e.:

$$ $\exists\phi^{*}_{t}:\quad\mathbb{D}\!\left[{q({\mathbf{x}}|{\phi_{t}^{*}})}\|{ p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})}\right]=0.$ (59) $$

The next assumption is a technical one to ensure that the thresholds will not diverge to infinity.

###### Assumption D.4 .

The sequence of thresholds is almost surely bounded:(^2^22We do not require $\tau_{*}$ to be known, only finite.)

$$ $\sup_{t\in\mathbb{N}}|\tau_{t}|\leq\tau_{*}<\infty\,.$ (60) $$

We can now state our main result regarding the GP-based approximations learned by VSD.

###### Theorem D.1 .

Let assumptions D.1 to D.4 hold. Then the following holds with probability 1 for VSD equipped with GP-PI:

$$ $\sigma_{t}^{2}(\mathbf{x})\in{\mathcal{{O}}}(t^{-1})\,,$ (61) $$

at every $\mathbf{x}\in\mathcal{X}$ such that $p(\mathbf{x})>0$.

###### Proof.

Let $\ell_{t}(\mathbf{x}):=p({y>\tau_{t}}|{\mathbf{x},\mathcal{D}_{t}})$. For any given $\mathbf{x}\in\mathcal{X}$ where $p({\mathbf{x}})>0$, by D.2, we have that the next candidate will be sampled according to:

$$ $\begin{split}\forall t\geq 0,\quad\mathbb{P}\!\left[\mathbf{x}_{t+1}=\mathbf{x }\mid{\mathfrak{{F}}}_{t}\right]&=p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}}) \\ &=\frac{\ell_{t}(\mathbf{x})p({\mathbf{x}})}{\mathbb{E}_{p({\mathbf{x}})}\! \left[{\ell_{t}(\mathbf{x})}\right]}\\ &\geq\ell_{t}(\mathbf{x})p({\mathbf{x}}),\end{split}$ (62) $$

where we used the fact that $\mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}(\mathbf{x})}\right]\leq 1$, since $\ell_{t}(\mathbf{x})\leq 1$, for all $\mathbf{x}\in\mathcal{X}$. As $p({\mathbf{x}})>0$, we only have to derive a lower bound on $\ell_{t}(\mathbf{x})$ to apply D.5 and derive a convergence rate.

A lower bound on $\ell_{t}(\mathbf{x})$ is given by:

$$ $\begin{split}\forall t\geq 0,\quad\ell_{t}(\mathbf{x})=\Psi\!\left(\frac{\mu_{ t}(\mathbf{x})-\tau_{t}}{\sqrt{\sigma_{t}^{2}(\mathbf{x})+\sigma_{\epsilon}^{2 }}}\right)\geq\Psi\!\left(-\frac{\lVert\mu_{t}\rVert_{\infty}+\tau_{*}}{\sigma _{\epsilon}}\right),\end{split}$ (63) $$

where $\Psi(\cdot)$ denotes the cumulative distribution function of a standard normal random variable, and $\lVert\cdot\rVert_{\infty}$ denotes the essential supremum of a function under $\mathbb{P}$ (the probability measure of the underlying abstract probability space). Therefore, if $\lim_{t\to\infty}\lVert\mu_{t}\rVert_{\infty}<\infty$, we will have that $\lim_{t\to\infty}\ell_{t}(\mathbf{x})>0$, and the sum in D.5 will diverge.

By Jensen's inequality for conditional expectations, we have that:

$$ $\forall t\geq 0,\quad\lVert\mu_{t}\rVert_{\infty}=\lVert\mathbb{E}[{f\! \centerdot}\mid{\mathfrak{{F}}}_{t}]\rVert_{\infty}\leq\mathbb{E}[\lVert{f\! \centerdot}\rVert_{\infty}\mid{\mathfrak{{F}}}_{t}].$ (64) $$

As $\mathbb{E}[\mathbb{E}[\lVert{f\!\centerdot}\rVert_{\infty}\mid{\mathfrak{{F}}}
_{t}]]=\mathbb{E}[\lVert{f\!\centerdot}\rVert_{\infty}]<\infty$ (cf. D.4), an application of Markov's inequality implies that:

$$ $\lim_{a\to\infty}\mathbb{P}\!\left[\mathbb{E}[\lVert{f\!\centerdot}\rVert_{ \infty}\mid{\mathfrak{{F}}}_{t}]\geq a\right]\leq\lim_{a\to\infty}\frac{1}{a} \mathbb{E}[\lVert{f\!\centerdot}\rVert_{\infty}]=0.$ (65) $$

Furthermore, $m_{t}:=\mathbb{E}[\lVert{f\!\centerdot}\rVert_{\infty}\mid{\mathfrak{{F}}}_{t}]$ also defines a non-negative martingale, and by the martingale convergence theorem (Durrett, , Thr. 4.2.11), $\lim_{t\to\infty}m_{t}=m_{\infty}:=\mathbb{E}[\lVert{f\!\centerdot}\rVert_{
\infty}\mid{\mathfrak{{F}}}_{\infty}]$ is well defined and $\mathbb{E}[\mathbb{E}[\lVert{f\!\centerdot}\rVert_{\infty}\mid{\mathfrak{{F}}}
_{\infty}]]=\mathbb{E}[\lVert{f\!\centerdot}\rVert_{\infty}]<\infty$. Again, by Markov's inequality, for any $a>0$, we have that:

$$ $\mathbb{P}\!\left[\lim_{t\to\infty}\lVert\mu_{t}\rVert_{\infty}\geq a\mathbb{E }[\lVert{f\!\centerdot}\rVert_{\infty}]\right]\leq\frac{\mathbb{E}\left[\lim_{ t\to\infty}\lVert\mu_{t}\rVert_{\infty}\right]}{a\mathbb{E}[\lVert{f\! \centerdot}\rVert_{\infty}]}\leq\frac{\mathbb{E}\left[\lim_{t\to\infty}\mathbb {E}[\lVert{f\!\centerdot}\rVert_{\infty}\mid{\mathfrak{{F}}}_{t}]\right]}{a \mathbb{E}[\lVert{f\!\centerdot}\rVert_{\infty}]}=\frac{1}{a}\,.$ (66) $$

Therefore, for any $a>0$ and any given $\mathbf{x}\in\mathcal{X}$, with probability at least $1-\frac{1}{a}$, the following holds:

$$ $\begin{split}\lim_{t\to\infty}\mathbb{P}\!\left[\mathbf{x}_{t}=\mathbf{x}\mid{ \mathfrak{{F}}}_{t-1}\right]&\geq p({\mathbf{x}})\lim_{t\to\infty}\ell_{t-1}( \mathbf{x})\\ &\geq p({\mathbf{x}})\lim_{t\to\infty}\Psi\!\left(-\frac{\lVert\mu_{t-1}\rVert _{\infty}+\tau_{*}}{\sigma_{\epsilon}}\right)\\ &\geq p({\mathbf{x}})\Psi\!\left(-\frac{a\mathbb{E}[\lVert{f\!\centerdot} \rVert_{\infty}]+\tau_{*}}{\sigma_{\epsilon}}\right)\\ &=:b_{\infty}(a)>0\,.\end{split}$ (67) $$

Hence, for any $\varepsilon_{a}\in(0,b_{\infty}(a))$, there is $N_{a}\in\mathbb{N}$, such that $\mathbb{P}\!\left[\mathbf{x}_{t}=\mathbf{x}\mid{\mathfrak{{F}}}_{t-1}\right]
\geq b_{\infty}(a)-\varepsilon_{a}>0$, for all $t\geq N_{a}$. As a result, $\sum_{t^{\prime}=1}^{t}\mathbb{P}\!\left[\mathbf{x}_{t^{\prime}}=\mathbf{x}
\mid{\mathfrak{{F}}}_{t^{\prime}-1}\right]\geq(b_{\infty}(a)-\varepsilon_{a})(
t-N_{a})$, for all $t\geq N_{a}$, which asymptotically diverges at a rate proportional to $t$. By D.5 and the definition of the big-${\mathcal{{O}}}$ notation, for any $\mathbf{x}\in\mathcal{X}$, we then have that:

$$ $\forall a>0,\quad\mathbb{P}\!\left[\limsup_{t\to\infty}\left\lvert t\sigma_{t} ^{2}(\mathbf{x})\right\rvert\leq\sigma_{\epsilon}^{2}<\infty\right]\geq 1- \frac{1}{a}\,.$ (68) $$

Taking the limit as $a\to\infty$, we can finally conclude that:

$$ $\mathbb{P}\!\left[\limsup_{t\to\infty}\left\lvert t\sigma_{t}^{2}(\mathbf{x}) \right\rvert<\infty\right]=1,$ (69) $$

i.e., $\sigma_{t}^{2}$ is almost surely ${\mathcal{{O}}}(t^{-1})$, which concludes the proof.
∎

###### Remark D.1 .

The convergence rate in Theorem D.1 is optimal and cannot be further improved. As shown by previous works in the online learning literature (Mutný & Krause, ; Takeno et al., ), a lower bound on the GP variance at each iteration $t\geq 1$ is given by $\sigma_{t}^{2}(\mathbf{x})\geq\sigma_{\epsilon}^{2}(\sigma_{\epsilon}^{2}+t)^{
-1}$ (assuming $k(\mathbf{x},\mathbf{x})=1$), which is the case when every observation in the dataset was collected at the same point $\mathbf{x}\in\mathcal{X}$ (see Takeno et al., , Lem. 4.2). Therefore, the lower and upper bounds on the asymptotic convergence rates for the GP variance differ by only up to a multiplicative constant.

The result in Theorem D.1 now allows us to derive a convergence rate for VSD's approximations to the level-set distributions. To do so, however, we will require the following mild assumption, which is satisfied by any prior distribution which has support on the entire domain $\mathcal{X}$.

###### Assumption D.5 .

The prior distribution is such that $p({\mathbf{x}})>0$, for all $\mathbf{x}\in\mathcal{X}$.

*

###### Proof.

We first prove an upper bound for the KL divergence in terms of the PI approximation error. We then derive a bound for this term and apply Theorem D.1 to obtain a convergence rate.

*KL bound formulation.* Let $\ell_{t}(\mathbf{x}):=p({y>\tau_{t}}|{\mathbf{x},\mathcal{D}_{t}})$ and $\ell_{t}^{*}(\mathbf{x}):=p({y>\tau_{t}}|{\mathbf{x},{f\!\centerdot}})$, for $\mathbf{x}\in\mathcal{X}$. From the definition of the KL divergence, we have that:

$$ $\begin{split}\mathbb{D}\!\left[{p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})} \|{p({\mathbf{x}}|{y>\tau_{t},{f\!\centerdot}})}\right]&=\mathbb{E}_{p({ \mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})}\!\left[{\log p({\mathbf{x}}|{y>\tau _{t},\mathcal{D}_{t}})-\log p({\mathbf{x}}|{y>\tau_{t},{f\!\centerdot}})} \right]\\ &=\mathbb{E}_{p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})}\!\left[{\log\ell_{ t}(\mathbf{x})-\log\ell_{t}^{*}(\mathbf{x})}\right]\\ &\quad+\log\mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}^{*}(\mathbf{x})} \right]-\log\mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}(\mathbf{x})}\right] \\ &=\mathbb{E}_{p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})}\!\left[{\log\left( \frac{\ell_{t}(\mathbf{x})}{\ell_{t}^{*}(\mathbf{x})}\right)}\right]+\log\left (\frac{\mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}^{*}(\mathbf{x})}\right]}{ \mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}(\mathbf{x})}\right]}\right).\end {split}$ (70) $$

For logarithms, we know that $\log(1+a)\leq a$, for all $a>-1$, which shows that:

$$ $\displaystyle\log\left(\frac{\ell_{t}(\mathbf{x})}{\ell_{t}^{*}(\mathbf{x})}\right)$ $\displaystyle=\log\left(1+\frac{\ell_{t}(\mathbf{x})-\ell_{t}^{*}(\mathbf{x})} {\ell_{t}^{*}(\mathbf{x})}\right)\leq\frac{\ell_{t}(\mathbf{x})-\ell_{t}^{*}( \mathbf{x})}{\ell_{t}^{*}(\mathbf{x})}$ (71) $\displaystyle\log\left(\frac{\mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}^{*} (\mathbf{x})}\right]}{\mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}(\mathbf{x} )}\right]}\right)$ $\displaystyle=\log\left(1+\frac{\mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}^ {*}(\mathbf{x})-\ell_{t}(\mathbf{x})}\right]}{\mathbb{E}_{p({\mathbf{x}})}\! \left[{\ell_{t}(\mathbf{x})}\right]}\right)\leq\frac{\mathbb{E}_{p({\mathbf{x} })}\!\left[{\ell_{t}^{*}(\mathbf{x})-\ell_{t}(\mathbf{x})}\right]}{\mathbb{E}_ {p({\mathbf{x}})}\!\left[{\ell_{t}(\mathbf{x})}\right]}\,.$ (72) $$

Combining the above into Equation 70 yields:

$$ $\begin{split}\mathbb{D}\!\left[{p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})} \|{p({\mathbf{x}}|{y>\tau_{t},{f\!\centerdot}})}\right]&\leq\mathbb{E}_{p({ \mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})}\!\left[{\frac{\ell_{t}(\mathbf{x})- \ell_{t}^{*}(\mathbf{x})}{\ell_{t}^{*}(\mathbf{x})}}\right]+\frac{\mathbb{E}_{ p({\mathbf{x}})}\!\left[{\ell_{t}^{*}(\mathbf{x})-\ell_{t}(\mathbf{x})}\right] }{\mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}(\mathbf{x})}\right]}\end{split }\,.$ (73) $$

The denominator in the expression above is such that:

$$ $\forall t\geq 0,\quad\ell_{t}^{*}(\mathbf{x})=p({y>\tau_{t}}|{\mathbf{x},{f\! \centerdot}})=\Psi\!\left(\frac{{f\!\centerdot}(\mathbf{x})-\tau_{t}}{\sigma_{ \epsilon}}\right)\geq\Psi\!\left(-\frac{\lVert{f\!\centerdot}\rVert_{\infty}+ \tau_{*}}{\sigma_{\epsilon}}\right),\quad\forall\mathbf{x}\in\mathcal{X}\,.$ (74) $$

By D.4, we know that $\mathbb{E}[\lVert{f\!\centerdot}\rVert_{\infty}]<\infty$, which implies that $\mathbb{P}\!\left[\lVert{f\!\centerdot}\rVert_{\infty}<\infty\right]=1$ by Markov's inequality. Next, we derive a bound for the approximation error term.

*Error bound.* We now derive an upper bound for the difference ${\Delta\ell}_{t}(\mathbf{x}):=\ell_{t}(\mathbf{x})-\ell_{t}^{*}(\mathbf{x})$ and then show that it asymptotically vanishes.
Applying Taylor's theorem to $\Psi$, we can bound ${\Delta\ell}_{t}$ as a function of the approximation error between the mean $\mu_{t}$ and the true function ${f\!\centerdot}$ as:

$$ $\begin{split}\forall t\geq 0,\quad|{\Delta\ell}_{t}(\mathbf{x})|&=\left\lvert \Psi\!\left(\frac{\mu_{t}(\mathbf{x})-\tau_{t}}{\sqrt{\sigma_{t}^{2}(\mathbf{x })+\sigma_{\epsilon}^{2}}}\right)-\Psi\!\left(\frac{{f\!\centerdot}(\mathbf{x} )-\tau_{t}}{\sigma_{\epsilon}}\right)\right\rvert\\ &\leq\frac{1}{\sqrt{2\pi}}\left\lvert\frac{\mu_{t}(\mathbf{x})-\tau_{t}}{\sqrt {\sigma_{t}^{2}(\mathbf{x})+\sigma_{\epsilon}^{2}}}-\frac{{f\!\centerdot}( \mathbf{x})-\tau_{t}}{\sigma_{\epsilon}}\right\rvert\\ &=\frac{1}{\sqrt{2\pi}}\left\lvert\frac{\sigma_{\epsilon}\mu_{t}(\mathbf{x})-{ f\!\centerdot}(\mathbf{x})\sqrt{\sigma_{t}^{2}(\mathbf{x})+\sigma_{\epsilon}^{ 2}}+\tau_{t}(\sqrt{\sigma_{t}^{2}(\mathbf{x})+\sigma_{\epsilon}^{2}}-\sigma_{ \epsilon})}{\sigma_{\epsilon}\sqrt{\sigma_{t}^{2}(\mathbf{x})+\sigma_{\epsilon }^{2}}}\right\rvert\\ &\leq\frac{\lvert\sigma_{\epsilon}\mu_{t}(\mathbf{x})-{f\!\centerdot}(\mathbf{ x})\sqrt{\sigma_{t}^{2}(\mathbf{x})+\sigma_{\epsilon}^{2}}\rvert+|\tau_{t}| \sigma_{t}(\mathbf{x})}{\sigma_{\epsilon}^{2}\sqrt{2\pi}}\\ &\leq\frac{\sigma_{\epsilon}\lvert\mu_{t}(\mathbf{x})-{f\!\centerdot}(\mathbf{ x})\rvert+\sigma_{t}(\mathbf{x})(|{f\!\centerdot}(\mathbf{x})|+|\tau_{t}|)}{ \sigma_{\epsilon}^{2}\sqrt{2\pi}},\quad\forall\mathbf{x}\in\mathcal{X},\\ \end{split}$ (75) $$

since $\sup_{\epsilon\in}\left\lvert\frac{{\mathop{}\operatorname{d}}\Psi(\epsilon)}{
{\mathop{}\operatorname{d}}\epsilon}\right\rvert=\frac{1}{\sqrt{2\pi}}<1$, and we used the fact that $\sigma_{\epsilon}\leq\sqrt{\sigma_{t}^{2}(\mathbf{x})+\sigma_{\epsilon}^{2}}
\leq\sigma_{t}(\mathbf{x})+\sigma_{\epsilon}$ to obtain the last two inequalities.

*Convergence rate.* To derive a convergence rate, given any $\mathbf{x}\in\mathcal{X}$ and $t\geq 0$, we have that:

$$ $\mathbb{E}[|{\Delta\ell}_{t}(\mathbf{x})|\mid{\mathfrak{{F}}}_{t}]\leq\frac{ \sigma_{\epsilon}\mathbb{E}[\lvert\mu_{t}(\mathbf{x})-{f\!\centerdot}(\mathbf{ x})\rvert\mid{\mathfrak{{F}}}_{t}]+\sigma_{t}(\mathbf{x})(\mathbb{E}[|{f\! \centerdot}(\mathbf{x})|\mid{\mathfrak{{F}}}_{t}]+|\tau_{t}|)}{\sigma_{ \epsilon}^{2}\sqrt{2\pi}}.$ (76) $$

We know that $\mathbb{E}[|{f\!\centerdot}(\mathbf{x})|\mid{\mathfrak{{F}}}_{t}]$ is almost surely bounded, and by Jensen's inequality, it also holds that:

$$ $\mathbb{E}[\lvert\mu_{t}(\mathbf{x})-{f\!\centerdot}(\mathbf{x})\rvert\mid{ \mathfrak{{F}}}_{t}]\leq\sigma_{t}(\mathbf{x}).$ (77) $$

Applying Theorem D.1, we then have that:

$$ $|{\Delta\ell}_{t}(\mathbf{x})|\in{\mathcal{{O}}}_{\mathbb{P}}(t^{-1/2}).$ (78) $$

Since $\lVert\mu_{t}\rVert_{\infty}\leq\mathbb{E}[\lVert{f\!\centerdot}\rVert_{\infty
}|{\mathfrak{{F}}}_{t}]\in{\mathcal{{O}}}_{\mathbb{P}}(1)$, we also have that:

$$ $\frac{1}{\mathbb{E}_{p({\mathbf{x}})}\!\left[{\ell_{t}(\mathbf{x})}\right]}\in {\mathcal{{O}}}_{\mathbb{P}}(1)\,.$ (79) $$

Lastly, we know that $\frac{1}{\ell_{t}^{*}(\mathbf{x})}\in{\mathcal{{O}}}_{\mathbb{P}}(1)$ by Equation 74 and the observation that $\lVert{f\!\centerdot}\rVert_{\infty}\in{\mathcal{{O}}}_{\mathbb{P}}(1)$. The main result then follows by combining the rates above into Equation 73.
∎

###### Assumption D.1 .

###### Lemma D.5 .

###### Proof.

###### Assumption D.2 .

###### Assumption D.3 .

###### Assumption D.4 .

###### Theorem D.1 .

###### Proof.

###### Remark D.1 .

###### Assumption D.5 .

###### Proof.

### D.5 Performance analysis

At every iteration $t\geq 1$, VSD samples $\mathbf{x}_{t}$ from (an approximation to) the target $p({\mathbf{x}}|{y>\tau_{t-1},\mathcal{D}_{t-1}})$ and obtains an observation $y_{t}\sim p({y}|{\mathbf{x}_{t}})$. A positive hit consists of an event $y_{t}>\tau_{t-1}$, where $\tau_{t-1}$ is computed based on the data available in $\mathcal{D}_{t-1}$ or a constant. Therefore, we can compute the probability of a positive hit for a given realization of ${f\!\centerdot}$ as:

$$ $\mathbb{P}\!\left[y_{t}>\tau_{t-1}\mid\mathcal{D}_{t-1},{f\!\centerdot}\right] =\mathbb{E}_{p({\mathbf{x}}|{y>\tau_{t-1},\mathcal{D}_{t-1}})}\!\left[{p({y> \tau_{t-1}}|{\mathbf{x},{f\!\centerdot}})}\right].$ (80) $$

Then the expected number of hits $H_{T}$ after $T\geq 1$ iterations is given by:

$$ $\mathbb{E}[H_{T}\mid{f\!\centerdot}]=\sum_{t=1}^{T}\mathbb{E}_{p({\mathbf{x}}| {y>\tau_{t-1},\mathcal{D}_{t-1}})}\!\left[{p({y>\tau_{t-1}}|{\mathbf{x},{f\! \centerdot}})}\right].$ (81) $$

We will compare this quantity with the expected number of hits $H_{T}^{*}$ obtained by a sampling distribution with full knowledge of the objective function ${f\!\centerdot}$:

$$ $\mathbb{E}[H_{T}^{*}\mid{f\!\centerdot}]=\sum_{t=1}^{T}\mathbb{E}_{p({\mathbf{ x}}|{y>\tau_{t-1},{f\!\centerdot}})}\!\left[{p({y_{t}>\tau_{t-1}}|{\mathbf{x}, {f\!\centerdot}})}\right].$ (82) $$

The next result allows us to bound the difference between these two quantities.

*

###### Proof.

For all $T\geq 1$, we have that:

$$ $\begin{split}\mathbb{E}[H_{T}-H_{T}^{*}]&=\mathbb{E}\!\left[\sum_{t=1}^{T} \mathbb{E}_{p({\mathbf{x}}|{y>\tau_{t-1},\mathcal{D}_{t-1}})}\!\left[{p({y> \tau_{t-1}}|{\mathbf{x},{f\!\centerdot}})}\right]-\mathbb{E}_{p({\mathbf{x}}|{ y>\tau_{t-1},{f\!\centerdot}})}\!\left[{p({y_{t}>\tau_{t-1}}|{\mathbf{x},{f\! \centerdot}})}\right]\right]\\ &=\mathbb{E}\!\left[\sum_{t=1}^{T}\sum_{\mathbf{x}\in\mathcal{X}}p({y>\tau_{t- 1}}|{\mathbf{x},{f\!\centerdot}})\left(p({\mathbf{x}}|{y>\tau_{t-1},\mathcal{D }_{t-1}})-p({\mathbf{x}}|{y>\tau_{t-1},{f\!\centerdot}})\right)\right]\\ &=\mathbb{E}\!\left[\sum_{t=0}^{T-1}\sum_{\mathbf{x}\in\mathcal{X}}p({y>\tau_{ t}}|{\mathbf{x},{f\!\centerdot}})p({\mathbf{x}})\left(\frac{\ell_{t}(\mathbf{x })}{\mathbb{E}_{p({\mathbf{x}^{\prime}})}\!\left[{\ell_{t}(\mathbf{x}^{\prime} )}\right]}-\frac{\ell_{t}^{*}(\mathbf{x})}{\mathbb{E}_{p({\mathbf{x}^{\prime}} )}\!\left[{\ell_{t}^{*}(\mathbf{x}^{\prime})}\right]}\right)\right]\\ &\leq\mathbb{E}\!\left[\sum_{t=0}^{T-1}\sum_{\mathbf{x}\in\mathcal{X}}p({ \mathbf{x}})\left(\frac{|{\Delta\ell}_{t}(\mathbf{x})|}{\min\{\mathbb{E}_{p({ \mathbf{x}^{\prime}})}\!\left[{\ell_{t}(\mathbf{x}^{\prime})}\right],\mathbb{E }_{p({\mathbf{x}^{\prime}})}\!\left[{\ell_{t}^{*}(\mathbf{x}^{\prime})}\right] \}}\right)\right],\end{split}$ (83) $$

since $p({y>\tau_{t-1}}|{\mathbf{x},{f\!\centerdot}})\leq 1$, for all $t\geq 1$. As both $\lVert\mu_{t}\rVert_{\infty}$ and $\lVert{f\!\centerdot}\rVert_{\infty}$ are in ${\mathcal{{O}}}_{\mathbb{P}}(1)$, $\min\{\mathbb{E}_{p({\mathbf{x}^{\prime}})}\!\left[{\ell_{t}(\mathbf{x}^{
\prime})}\right],\mathbb{E}_{p({\mathbf{x}^{\prime}})}\!\left[{\ell_{t}^{*}(
\mathbf{x}^{\prime})}\right]\}$ is lower bounded by some constant. As ${\Delta\ell}_{t}(\mathbf{x})\in{\mathcal{{O}}}_{\mathbb{P}}(t^{-1/2})$, for $T$ large enough and some $C>0$, we then have that:

$$ $\begin{split}\mathbb{E}[|H_{T}-H_{T}^{*}|]\leq C\sum_{t=1}^{T}\frac{1}{\sqrt{t }}\leq 2C\sqrt{T}\in{\mathcal{{O}}}(\sqrt{T}),\end{split}$ (84) $$

which follows by an application of the Euler-Maclaurin formula, since $\int_{1}^{T}\frac{1}{\sqrt{t}}{\mathop{}\operatorname{d}}t=2\sqrt{T}-2$ and the remainder term asymptotically vanishes.
∎

###### Remark D.2 .

If the oracle achieves $\mathbb{E}[H_{T}^{*}]=T$, the error bound in Sec. 2.4 suggests an increasing rate of positive hits by VSD as $\frac{1}{T}\mathbb{E}[H_{T}]\geq 1-CT^{-1/2}$, for some constant $C>0$ and large enough $T$. Therefore, VSD should asymptotically achieve a full rate of 1 positive hit per iteration in the single-point batch setting we consider. Note, however, that the results above do not discount for repeated samples, though should still indicate that VSD achieves a high discovery rate over the course of its execution.

###### Proof.

###### Remark D.2 .

## Appendix E VSD with neural network CPEs

In this section, we consider VSD with class probability estimators that are not based on GP regression, which was the case for the previous section, while specifically focusing on neural network models. We will, however, show that with a kernel-based formulation we are able to capture the classification models based on neural networks which we use. This is possible by analyzing the behavior of infinite-width neural networks (Jacot et al., ; Lee et al., ), whose approximation error with respect to the finite-width model can be bounded (Liu et al., ; Eldan et al., ).

Although our classifiers are learned by minimizing the cross-entropy (CE) loss, we can connect their approximations with theoretical results from the infinite-width neural network (NN) literature, which are mostly based on the mean squared error (MSE) loss. Recall that, given a dataset $\mathcal{D}_{N}^{z}:=\{(\mathbf{x}_{n},z_{n})\}_{n=1}^{N}$ with binary labels $z_{n}\in\{0,1\}$, the cross-entropy loss for a probabilistic classifier $\pi_{\theta}:\mathcal{X}\to[0,1]$ parameterized by $\theta$ is given by(^3^33We implicitly assume that $0<\pi_{\theta}({\mathbf{x}_{n}})<1$, for $n\in\{1,\dots,N\}$, so that the CE loss is well defined. This assumption can, however, be relaxed when dealing with the MSE loss, which remains well defined otherwise.):

$$ $\displaystyle\mathcal{L_{\textrm{CPE}}}\!\left({\theta,\mathcal{D}_{N}^{z}} \right):=-\frac{1}{N}\sum_{n=1}^{N}z_{n}\log\pi_{\theta}({\mathbf{x}_{n}})+(1- z_{n})\log(1-\pi_{\theta}({\mathbf{x}_{n}}))\,.$ (85) $$

The MSE loss for the same model corresponds to:

$$ $\mathcal{L}_{\mathrm{MSE}}(\theta,\mathcal{D}_{N}^{z}):=\frac{1}{N}\sum_{n=1}^ {N}(z_{n}-\pi_{\theta}({\mathbf{x}_{n}}))^{2}\,.$ (86) $$

The following result establishes a connection between the two loss functions.

###### Proposition E.1 .

Given a binary classification dataset $\mathcal{D}_{N}^{z}$ of size $N\geq 1$, the following holds for the cross-entropy and the mean-square error losses:

$$ $\mathcal{L_{\textrm{CPE}}}\!\left({\theta,\mathcal{D}_{N}^{z}}\right)\geq \mathcal{L}_{\mathrm{MSE}}(\theta,\mathcal{D}_{N}^{z}),\quad\forall N\in \mathbb{N}\,.$ (87) $$

###### Proof.

Applying the basic logarithmic inequality $\log(1+a)\leq a$, for all $a>-1$, to the cross-entropy loss definition yields:

$$ $\begin{split}\mathcal{L_{\textrm{CPE}}}\!\left({\theta,\mathcal{D}_{N}^{z}} \right)&:=-\frac{1}{N}\sum_{n=1}^{N}z_{n}\log\pi_{\theta}({\mathbf{x}_{n}})+(1 -z_{n})\log(1-\pi_{\theta}({\mathbf{x}_{n}}))\\ &\geq-\frac{1}{N}\sum_{n=1}^{N}z_{n}(\pi_{\theta}({\mathbf{x}_{n}})-1)-(1-z_{n })\pi_{\theta}({\mathbf{x}_{n}})\\ &=-\frac{1}{N}\sum_{n=1}^{N}2z_{n}\pi_{\theta}({\mathbf{x}_{n}})-z_{n}-\pi_{ \theta}({\mathbf{x}_{n}})\\ &=\frac{1}{N}\sum_{n=1}^{N}z_{n}-2z_{n}\pi_{\theta}({\mathbf{x}_{n}})+\pi_{ \theta}({\mathbf{x}_{n}})\,.\end{split}$ (88) $$

Now note that $z_{n}=z_{n}^{2}$, for $z_{n}\in\{0,1\}$, and $\pi_{\theta}({\mathbf{x}_{n}})\geq\pi_{\theta}({\mathbf{x}_{n}})^{2}$, as $\pi_{\theta}({\mathbf{x}_{n}})\in[0,1]$, for all $n\in\{1,\dots,N\}$. Making these substitutions in Equation 88, we obtain:

$$ $\mathcal{L_{\textrm{CPE}}}\!\left({\theta,\mathcal{D}_{N}^{z}}\right)\geq\frac {1}{N}\sum_{n=1}^{N}z_{n}^{2}-2z_{n}\pi_{\theta}({\mathbf{x}_{n}})+\pi_{\theta }({\mathbf{x}_{n}})^{2}=\mathcal{L}_{\mathrm{MSE}}(\theta,\mathcal{D}_{N}^{z})\,,$ (89) $$

which concludes the proof.
∎

The result in E.1 suggests that minimizing the cross-entropy loss will lead us to minimize the MSE loss as well, since the latter is upper bounded by the former. This result provides us with theoretical justification to derive convergence results based on the MSE loss, which has been better analyzed in the NN literature (Jacot et al., ; Lee et al., ), as a proxy to establish convergence guarantees for the CE-based VSD setting.

###### Proposition E.1 .

###### Proof.

### E.1 Linear approximations via the neural tangent kernel

Let $\pi^{*}$ denote the unknown true classifier, i.e., $\pi(\mathbf{x}):=p({y>\tau}|{\mathbf{x},{f\!\centerdot}})$, for $\mathbf{x}\in\mathcal{X}$. In the following, we will assume that $\pi^{*}$ is an unknown, fixed element of a reproducing kernel Hilbert space (RKHS) associated with a given kernel (Schölkopf & Smola, ). In the case of infinite-width neural networks, we know that under certain assumptions the NN trained via gradient descent under the MSE loss will asymptotically converge to a kernel ridge regression solution whose kernel is given by the neural tangent kernel (NTK, Jacot et al. ()). This asymptotic solution is equivalent to the posterior mean of a Gaussian process that assumes no observation noise. For a finite amount of training steps with a non-infinitesimal learning rate, however, the literature has shown that gradient-based training provides a form of implicit regularization. In that case, we recover a regularized kernel ridge regression solution which can be robust to label noise (Hu et al., ) for us to use. Lastly, although our analysis will be based on NTK results, the approximation error between the infinite-width and the finite-width NN vanishes with the square root of the network width for most popular NN architectures (Liu et al., ). Therefore, we can assume that these approximation guarantees will remain useful for wide-enough, finite-width NN models.

#### Implicit regularization.

Several results in the literature have shown that training overparameterized neural networks via gradient descent provides a form of implicit regularization on the learned model (Fleming, ; Yao et al., ; Soudry et al., ; Barrett & Dherin, ), with some of the same behavior extending to the stochastic gradient setting (Smith et al., ). In particular, Fleming () showed a direct equivalence between an early stopped gradient-descent linear model and the solution of a regularized least-squares problem with a penalty on the parameters vector Euclidean norm. Therefore, as wide deep neural networks behave as linear models in the infinite-width limit when trained via gradient descent (Jacot et al., ; Lee et al., ), it is reasonable to model our least-squares problem with the MSE loss via its regularized version:

$$ $\hat{\pi}_{N}\in\operatorname*{arg\,\!min}_{\pi\in{\mathcal{{F}}}_{k}}\sum_{i= 1}^{N}(\pi(\mathbf{x}_{i})-z_{i})^{2}+\rho\lVert\pi\rVert_{k}^{2},$ (90) $$

where $\rho>0$ is a possibly unknown regularization factor, ${\mathcal{{F}}}_{k}$ denotes the RKHS associated with the neural tangent kernel $k$ for the given CPE NN architecture, and $\lVert\cdot\rVert_{k}$ denotes the RKHS norm. The problem above is equivalent to regularized kernel ridge regression (Shawe-Taylor & Cristianini, ), whose solution is given by:

$$ $\hat{\pi}_{N}(\mathbf{x})=\pi_{0}(\mathbf{x})+{\boldsymbol{\mathbf{k}}}_{N}( \mathbf{x})^{\top}({\boldsymbol{\mathbf{{K}}}}_{N}+\rho\mathbf{I})^{-1}({ \boldsymbol{\mathbf{z}}}_{N}-\pi_{0}(\mathcal{X}_{N})),$ (91) $$

where ${\boldsymbol{\mathbf{k}}}_{N}$ and ${\boldsymbol{\mathbf{{K}}}}_{N}$ are defined in the same way as for the GP case, $\pi_{0}$ denotes the untrained NN classifier at initialization, and $\pi_{0}(\mathcal{X}_{N}):=[\pi_{0}(\mathbf{x}_{i})]_{i=1}^{N}$. For our analysis, we may assume that $\pi_{0}=0$ at times, noting that the least-squares problem can always be solved for the residuals $z-\pi_{0}(\mathbf{x})$ and then have $\pi_{0}$ added back to the solution. We refer the reader to Jacot et al. () for further discussion on the effect of the network initialization.

#### Approximation for finite-width networks.

For fully connected, convolutional or residual networks equipped with smooth activation functions (e.g., sigmoid or $\tanh$), Liu et al. () showed that the approximation error between the linear model and the finite-width NN is $\widetilde{{\mathcal{{O}}}}(m^{-1/2})$, where $m$ denotes the minimum layer width, and the $\widetilde{{\mathcal{{O}}}}$ notation corresponds to the ${\mathcal{{O}}}$-notation with logarithmic factors suppressed. NTK results for other activation functions, e.g., ReLU (Chen & Xu, ), and different neural network architectures, such as multi-head attention (Hron et al., ), are also available in the literature.

### E.2 Assumptions

In the following, we present a series of mild technical assumptions needed for our theoretical analysis of NN-based CPEs. Firstly, we assume a bounded NTK.

###### Assumption E.1 .

The NTK $k$ corresponding to the network architecture in $\pi_{\theta}$ is bounded:

$$ $\sup_{\mathbf{x}\in\mathcal{X}}k(\mathbf{x},\mathbf{x})<\infty\,,$ (92) $$

for some constant $b_{k}>0$.

We will also assume that the threshold is fixed to simplify the analysis. However, our results should still be applicable to the time-varying threshold setting after minor adjustments.

###### Assumption E.2 .

The threshold is fixed, i.e., $\tau_{t}=\tau\in$, for all $t\geq 1$.

The following assumption on label noise should always hold for Bernoulli random variables (Boucheron et al., ). Any upper bound on the sub-Gaussian parameter should suffice for the analysis (e.g., $\sigma_{z}\leq 1$ for Bernoulli variables).

###### Assumption E.3 .

For all $t\in\mathbb{N}$ and all $\mathbf{x}\in\mathcal{X}$, label noise $\zeta=\mathbbm{1}[{y>\tau}]-\pi^{*}(\mathbf{x})$, with $y\sim p({y}|{\mathbf{x},{f\!\centerdot}})$, is $\sigma_{z}$-sub-Gaussian:

$$ $\forall a\in,\quad\mathbb{E}\left[\exp\left(a\zeta\right)\right]\leq\exp\left( \frac{a^{2}\sigma_{z}^{2}}{2}\right),$ (93) $$

for some $\sigma_{z}\geq 0$.

For this analysis, we mainly assume that the true classifier $\pi^{*}(\mathbf{x})=p({y>\tau}|{\mathbf{x},{f\!\centerdot}})$ is a fixed, though unknown, element of the RKHS given by the NTK, which is formalized by the following assumption.

###### Assumption E.4 .

There is $\pi^{*}\in{\mathcal{{F}}}_{k}$ such that:

$$ $\pi^{*}(\mathbf{x})=p({y>\tau}|{\mathbf{x},{f\!\centerdot}}),\quad\forall \mathbf{x}\in\mathcal{X}.$ (94) $$

For a rich enough RKHS, such assumption is mild, especially given that most popular NN architectures possess universal approximation guarantees (Hornik et al., ). Finally, the next assumption ensures enough sampling asymptotically over the domain $\mathcal{X}$, which we still assume is finite.

###### Assumption E.5 .

For any $t\geq 1$, the variational family is such that sampling probabilities are bounded away from 0, i.e.:

$$ $\exists b>0:\quad\forall t\in\mathbb{N},\quad q({\mathbf{x}}|{\phi_{t}})\geq b\,.$ (95) $$

The assumption above only imposes mild constraints on the generative models $q({\mathbf{x}}|{\phi})$, so that probabilities for all candidates $\mathbf{x}\in\mathcal{X}$ are never exactly 0, though still allowed to be arbitrarily small.

###### Assumption E.1 .

###### Assumption E.2 .

###### Assumption E.3 .

###### Assumption E.4 .

###### Assumption E.5 .

### E.3 Approximation error for NN-based CPEs

Similar to the GP-PI setting, we will assume a batch size of 1, so that we can simply use the iteration index $t\geq 0$ for our estimators. We recall that convergence rates for the batch setting should only be affected by a batch-size-dependent multiplicative factor, preserving big-${\mathcal{{O}}}$ convergence rates. We start by defining the following *proxy* variance:

$$ $t\geq 1,\quad\hat{\sigma}_{t}^{2}(\mathbf{x})=k(\mathbf{x},\mathbf{x})-{ \boldsymbol{\mathbf{k}}}_{t}(\mathbf{x})^{\top}({\boldsymbol{\mathbf{{K}}}}_{t }+\rho\mathbf{I})^{-1}{\boldsymbol{\mathbf{k}}}_{t}(\mathbf{x})\,,\quad\mathbf {x}\in\mathcal{X},$ (96) $$

which is equivalent to a GP posterior variance when observation noise is assumed to be zero-mean Gaussian with variance given by $\rho$. Given its similarities, we have that if enough sampling is asymptotically guaranteed, we can apply the same convergence results available for the GP-PI-based CPE, i.e., $\hat{\sigma}_{t}^{2}\in{\mathcal{{O}}}(t^{-1/2})$ almost surely. We then invoke E.5 to derive the following result.

###### Lemma E.1 .

Let E.1 and E.5 hold. Then the following almost surely holds for the proxy variance:

$$ $\hat{\sigma}_{t}^{2}\in{\mathcal{{O}}}(t^{-1})\,.$ (97) $$

###### Proof.

The proof follows by verifying that the sum of sampling probabilities at any point $\mathbf{x}\in\mathcal{X}$ diverges as $t\to\infty$ by E.5, and then by D.5 the result follows.
∎

###### Lemma E.2 .

Let assumptions E.1 to E.4 hold. Then, given any $\delta\in(0,1]$, the following holds with probability at least $1-\delta$ for the approximation error between $\hat{\pi}_{t}$ and $\pi^{*}$:

$$ $\forall t\geq 1,\quad|\hat{\pi}_{t}(\mathbf{x})-\pi^{*}(\mathbf{x})|\leq\beta_ {t}(\delta)\hat{\sigma}_{t}(\mathbf{x}),\quad\mathbf{x}\in\mathcal{X},$ (98) $$

where $\beta_{t}(\delta):=\lVert\pi^{*}\rVert_{k}+\sigma_{\zeta}\sqrt{2\rho^{-1}\log(
\det(\mathbf{I}+\rho^{-1}{\boldsymbol{\mathbf{{K}}}}_{t})^{1/2}/\delta)}$.

###### Proof.

The result above is a direct application of Theorem 3.5 in Maillard () which provides an upper confidence bound on the kernelized least-squares regressor approximation error (another version of the same result is also available in Durand et al. (, Thr. 1)).
∎

For the next result, we need to define the following quantity:

$$ $\xi_{T}:=\max_{\mathcal{X}_{T}\subset\mathcal{X}:|\mathcal{X}_{T}|\leq T}\frac {1}{2}\log\det(\mathbf{I}+\rho^{-1}{\boldsymbol{\mathbf{{K}}}}(\mathcal{X}_{T} ))\,,$ (99) $$

where ${\boldsymbol{\mathbf{{K}}}}(\mathcal{X}_{T}):=[k(\mathbf{x},\mathbf{x}^{\prime
})]_{\mathbf{x},\mathbf{x}^{\prime}\in\mathcal{X}_{T}}\in^{|\mathcal{X}_{T}|
\times|\mathcal{X}_{T}|}$. Note that $\xi_{T}$ corresponds to the maximum information gain of a GP model (Srinivas et al., ) with covariance function given by the NTK, assuming Gaussian observation noise with variance given by $\rho$. Then $\xi_{T}$ is mainly dependent on the eigenvalue decay of the kernel under its spectral decomposition (Vakili et al., ). For the spectrum of the NTK, a few results are available in the literature (Murray et al., ).

###### Proposition E.2 .

Let assumptions E.1 to E.5 hold. Then, given $\delta\in(0,1]$, the following holds with probability at least $1-\delta$ for VSD equipped with a wide enough NN-based CPE model $\hat{\pi}_{t}$:

$$ $\mathbb{D}\!\left[{p({\mathbf{x}}|{y>\tau_{t},\mathcal{D}_{t}})}\|{p({\mathbf{ x}}|{y>\tau_{t},{f\!\centerdot}})}\right]\in{\mathcal{{O}}}\left(\sqrt{\frac{ \xi_{t}}{t}}\right)\,.$ (100) $$

###### Proof.

The result follows by applying the same steps as in the proof of Sec. 2.4. We note that $\ell_{t}^{*}(\mathbf{x})=\pi^{*}(\mathbf{x})>0$, due to observation noise, so that $\ell_{t}^{*}(\mathbf{x})^{-1}\in{\mathcal{{O}}}_{\mathbb{P}}(1)$. Similarly, E.2 implies that $|\hat{\pi}_{t}(\mathbf{x})-\pi^{*}(\mathbf{x})|\leq\beta_{t}(\delta)\sigma_{t}
(\mathbf{x})$ with probability at least $1-\delta$ simultaneously over all $\mathbf{x}\in\mathcal{X}$, so that ratio-dependent terms in Sec. 2.4 should remain bounded in probability. The upper bound in the result then follows by noticing that in our case $|{\Delta\ell}_{t}(\mathbf{x})|\leq\beta_{t}(\delta)\hat{\sigma}_{t}(\mathbf{x})$ with high probability, where $\hat{\sigma}_{t}\in{\mathcal{{O}}}(t^{-1/2})$ by E.1, and $\beta_{t}(\delta)\in{\mathcal{{O}}}(\sqrt{\xi_{t}})$ by E.2 and the definition of $\xi_{t}$ in Equation 99.
∎

The result above tells us that VSD is able to recover a similar asymptotic convergence guarantee to the one we derived for the GP-PI case, depending on the choice of NN architecture and more specifically on the spectrum of its associated NTK. In the case of a fully connected multi-layer ReLU network, for example, Chen & Xu () showed an equivalence between the RKHS of the ReLU NTK and that of the Laplace kernel $k(\mathbf{x},\mathbf{x}^{\prime})=\exp(-C\lVert\mathbf{x}-\mathbf{x}^{\prime}\rVert)$. As the latter is equivalent to a Matérn kernel with smoothness parameter set to $0.5$ (Rasmussen & Williams, ), the corresponding information gain bound is $\xi_{t}\in\widetilde{\mathcal{{O}}}(t^{\frac{d}{1+d}})$, where $d$ here denotes the dimensionality of the domain $\mathcal{X}$ (Vakili & Olkhovskaya, ). In the case of discrete sequences of length $M$, the dimensionality of $\mathcal{X}$ is determined by $M$. Therefore, in this case, we have proven Sec. 2.4.(^4^44Here $\widetilde{{\mathcal{{O}}}}_{\mathbb{P}}$ suppresses logarithmic factors, as in $\widetilde{{\mathcal{{O}}}}$.)

*

Similar steps can be applied to derive convergence guarantees for VSD with other neural network architectures based on the eigenspectrum of their NTK (Murray et al., ) and following the recipe in, e.g., Vakili et al. () or Srinivas et al. ().

###### Lemma E.1 .

###### Proof.

###### Lemma E.2 .

###### Proof.

###### Proposition E.2 .

###### Proof.

## Appendix F VSD as a Black-Box Optimization Lower Bound

A natural question to ask is how VSD relates to the BO objective for probability of improvement (Garnett, , Ch.7),

$$ $\displaystyle\mathbf{x}^{*}_{t}=\operatorname*{arg\,\!max}_{\mathbf{x}}\log \alpha_{PI}\!\left({\mathbf{x},\mathcal{D}_{N},\tau}\right).$ (101) $$

Firstly, we can see that the expected log-likelihood of term of Equation 7 lower-bounds this quantity.

###### Proposition F.1 .

For a parametric model, $q({\mathbf{x}}|{\phi})$, given $\phi\in\Phi\subseteq^{m}$ and $q\in\mathcal{P}:\mathcal{X}\times\Phi\to[0,1]$,

$$ $\displaystyle\max_{\mathbf{x}}\log\alpha_{PI}\!\left({\mathbf{x},\mathcal{D}_{ N},\tau}\right)\geq\max_{\phi}\mathbb{E}_{q({\mathbf{x}}|{\phi})}\!\left[{\log \alpha_{PI}\!\left({\mathbf{x},\mathcal{D}_{N},\tau}\right)}\right],$ (102) $$

and the bound becomes tight as $q({\mathbf{x}}|{\phi^{*}_{t}})\to\delta(\mathbf{x}^{*}_{t})$, a Dirac delta function at the maximizer $\mathbf{x}^{*}_{t}$.

Taking the argmax of the RHS will result in the variational distribution collapsing to a delta distribution at $\mathbf{x}^{*}_{t}$ for an appropriate choice of $q({\mathbf{x}}|{\phi})$. The intuition for Equation 102 is that the expected value of a random variable is always less than or equal to its maximum. The proof of this is in Daulton et al. (); Staines & Barber (). Extending this lower bound, we can show the following.

###### Proposition F.2 .

For a divergence $\mathbb{D}:\mathcal{P}(\mathcal{X})\times\mathcal{P}(\mathcal{X})\to[0,\infty)$, and a prior $p_{0}\in\mathcal{P}(\mathcal{X})$,

$$ $\displaystyle\max_{\mathbf{x}}\log\alpha_{PI}\!\left({\mathbf{x},\mathcal{D}_{ N},\tau}\right)\geq\max_{\phi}\mathbb{E}_{q({\mathbf{x}}|{\phi})}\!\left[{\log \alpha_{PI}\!\left({\mathbf{x},\mathcal{D}_{N},\tau}\right)}\right]-\mathbb{D} \!\left[{q({\mathbf{x}}|{\phi})}\|{p_{0}(\mathbf{x})}\right].$ (103) $$

We can see that this bound is trivially true given the range of divergences, and this covers VSD as a special case. However, this bound is tight if and only if $p_{0}$ concentrates as a Dirac delta at $\mathbf{x}_{t}^{*}$ with an appropriate choice of $q({\mathbf{x}}|{\phi})$. In any case, the lower bound remains valid for any choice of informative prior $p_{0}$ or even a uninformed prior, which allows us to maintain the framework flexible to incorporate existing prior information whenever that is available.

###### Proposition F.1 .

###### Proposition F.2 .

## Peer Reviews

confidence: 5
contribution: 3
initial_rating: 3
presentation: 3
questions: The following questions are sincere:

- Who is the audience for this paper? 

- What questions is this paper answering?

- What does the variational inference framing get us in the end? Access to a set of tools for theoretical analysis?
soundness: 4
strengths: This paper demonstrates a clarity of thought and composition that is commendable, I particularly enjoyed the related work section.

Likewise I do not have any major concerns regarding the technical soundness of the results presented.

As a good conceptual introduction to the topic, I think this draft could be useful to researchers new to the topic with some revisions.
summary: This paper casts sequential black-box optimization as a variational inference (i.e. amortized optimization) problem, and uses this perspective to unify a collection of different black-box optimization algorithms under a common theoretical framework and presents some proof of concept results on easy sequence optimization tasks.

#### UPDATE - 12/02/2024 ####
After extended discussion and an exceptionally thorough response from the authors, my concerns have been addressed and I am recommending acceptance. I believe this paper presents a very nice conceptual view of the topic of active generation and the empirical results raise interesting questions. I encourage my fellow reviewers to review the updated manuscript and re-evaluate their scores. I have left my original review unaltered for any interested readers.
weaknesses: I have two general impressions of this paper. 

First, it seems like the authors have not really chosen a direction for the paper. There are at least three different directions here, A) a unifying view of sequential black box optimization algorithms, B) a practical algorithm for sequential BBO, and C) theoretical analysis of convergence rates of a particular sequential BBO algorithm under strong assumptions. I would suggest you pick no more than two directions, preferably one. I actually think this particular subfield could really benefit from a more holistic perspective of the work that has been done, as I constantly see minor variations of these algorithms in my social media feed and review stack with no apparent awareness of the relationships between them. From what I can tell from this draft, it seems that A and C likely play more to your strengths.


Second, the authors seem blissfully unaware of a substantial body of work on this topic. To be quite candid, the paper reads like it was written circa September 2021. This is not mere rhetoric. The most recent baseline the authors consider was published at ICML 2021. It is also odd that two of the baselines you did include, DbAS and CbAS, are not even designed for the sequential setting. As a very active researcher in this exact area, I struggle to understand who this paper is for and how the authors pictured their place in the broader dialogue on this topic. I am sure you worked very hard on this paper and I commend your effort, but I honestly believe the best advice I can give you is to talk to more people working on this topic, preferably from outside your immediate academic circle. While it is difficult to hear this feedback, one of the functions of peer review is to reveal "unknown unknowns". I want to be sure this review is constructive, so I will provide some key references if you are serious about diving into this topic. You should also consider making use of tools like [Connected Papers](https://www.connectedpapers.com/) to improve your literature review process and avoid this situation in the future. 

You can start with [A survey and benchmark of high-dimensional Bayesian optimization of discrete sequences](https://arxiv.org/abs/2406.04739). This work is the most up-to-date complete survey on the topic I have seen, and the benchmarking rigor is notably good. This paper is associated with two repositories, [poli](https://github.com/MachineLearningLifeScience/poli) and [poli-baselines](https://github.com/MachineLearningLifeScience/poli-baselines). The former contains a suite of test functions that are much more up to date than the combinatorially complete landscapes considered in this paper, and the latter contains a suite of baseline solvers. You may even want to consider contributing your method as a solver to poli-baselines at some point.

Some key axes of variation to consider: 

How is the optimization problem solved? Most fall into one of three categories, directed evolution (which you seem to be familiar with based on your inclusion of AdaLead and PEX), generative search with explicit guidance, e.g. [2, 3, 4, 5, 6], and generative search with implicit guidance [7, 8], which can also be seen as a kind of amortized search. I could cite more papers but I believe I have made my point. Algorithms also differ in their handling of constraints, and their approach to managing the feedback covariate shift induced by online active data collection by an agent. 

In particular I will draw your attention to [a tutorial for LaMBO-2](https://github.com/prescient-design/cortex/blob/main/tutorials/4_guided_diffusion.ipynb) if you want to start considering more up to date baselines, however I would recommend using the solver interface provided in poli-baselines for actual experiments. You may also be interested in Ehrlich functions if you would like a convenient test function that is much more difficult to solve than small combinatorially complete landscapes but still easy to work with [9]. Ehrlich functions are available in [a small standalone package](https://github.com/prescient-design/holo-bench) or [as part of the poli package](https://machinelearninglifescience.github.io/poli-docs/using_poli/objective_repository/ehrlich_functions.html).

While I'm sure this is not the outcome you hoped for, science is a dialogue, and good science requires awareness of what is happening outside your academic niche. Hopefully my feedback is clear and actionable enough to benefit this work and your progression as a scientist.

References

- [1] González-Duque, M., Michael, R., Bartels, S., Zainchkovskyy, Y., Hauberg, S., & Boomsma, W. (2024). A survey and benchmark of high-dimensional Bayesian optimization of discrete sequences. arXiv preprint arXiv:2406.04739.
- [2] Tripp, A., Daxberger, E., & Hernández-Lobato, J. M. (2020). Sample-efficient optimization in the latent space of deep generative models via weighted retraining. Advances in Neural Information Processing Systems, 33, 11259-11272.
- [3] Stanton, S., Maddox, W., Gruver, N., Maffettone, P., Delaney, E., Greenside, P., & Wilson, A. G. (2022, June). Accelerating bayesian optimization for biological sequence design with denoising autoencoders. In International Conference on Machine Learning (pp. 20459-20478). PMLR.
- [4] Gruver, N., Stanton, S., Frey, N., Rudner, T. G., Hotzel, I., Lafrance-Vanasse, J., ... & Wilson, A. G. (2023). Protein design with guided discrete diffusion. Advances in neural information processing systems, 36.
- [5] Maus, N., Jones, H., Moore, J., Kusner, M. J., Bradshaw, J., & Gardner, J. (2022). Local latent space bayesian optimization over structured inputs. Advances in neural information processing systems, 35, 34505-34518.
- [6] Maus, N., Wu, K., Eriksson, D., & Gardner, J. (2022). Discovering many diverse solutions with bayesian optimization. arXiv preprint arXiv:2210.10953.
- [7] Tagasovska, N., Gligorijević, V., Cho, K., & Loukas, A. (2024). Implicitly Guided Design with PropEn: Match your Data to Follow the Gradient. arXiv preprint arXiv:2405.18075.
- [8] Chen, A., Stanton, S. D., Alberstein, R. G., Watkins, A. M., Bonneau, R., Gligorijevi, V., ... & Frey, N. C. (2024). LLMs are Highly-Constrained Biophysical Sequence Optimizers. arXiv preprint arXiv:2410.22296.
- [9] Stanton, S., Alberstein, R., Frey, N., Watkins, A., & Cho, K. (2024). Closed-Form Test Functions for Biophysical Sequence Optimization Algorithms. arXiv preprint arXiv:2407.00236.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: * What's 'x' in the title of Figure 1?
* What are the limitations of this approach?
* How is diversity within a batch enforced? 
* The reverse KLD is known to result in mode collapse. Why wasn't this an issue?
* Which variation reduction method did you use for the gradient estimator?
soundness: 3
strengths: * The problem is important as it has applications in pharmaceutical drugs/enzyme design.
* The paper paper is well written and the method is sound
* Experimental results on high dimensional datasets demonstrate superiority of the approach
summary: The authors propose a black box varoiational inference approachfor discrete designs generation. The authors derive asymptotic convergence rates for learning the true conditional generative distribution of designs. Compelling results on high dimensional sequence-design problems are demonstrated.
weaknesses: * The method lacks novelty, it's based on putting together blocks that have already been proposed in the litterature
* The paper clarity can be improved with an overview plot of the method

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 2
questions: - How robust are the results to the selection of the threshold $\tau$ and the batch size $B$?
- While the reviewer is not familiar with the field, could the authors give some intuitions about the difference between VSD and active learning approaches like Bayesian optimization, and why VSD is better?
soundness: 3
strengths: - The paper formulates the batch active search problem in the variational inference framework and provides theoretical guarantees to the learned distribution based on the sequentially attained data.
- Experimental results on real-world biological datasets demonstrate the practical use of the algorithm and its effectiveness to solve the problem.
summary: The paper develops the variational search distribution method to solve the active search problem in biological design. VSD estimates the super level-set distribution in a sequential manner by generating batches of data points for sequential experiments. Empirical results on optimizing protein fitness in several datasets showcase the effectiveness of VSD.
weaknesses: - The precision of VSD and most other methods is decreasing with more rounds in TrpB and TFBIND8 datasets while the recall values are in general low. However, an ideal method should achieve a better estimation of the ground truth super level-set distribution as more samples are collected. This may be due to the initial training set size being too large or the fitness landscape being easy to model. How do the models perform with a smaller initial training set size?
- How is VSD compared with the simple and commonly used directed evolution method?

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1.Since your algorithm heavily relies on another model ($p(x | D)$), I would be highly interested in better understanding the influence of a good prior on your variational distribution.
2. Regarding the GFP experiments, do you sample already existing sequences ? What is the influence of the relative poor performance of the oracle on ood data on the interpretation of the results  ?
3. How can you explain that only a very simple prior such as a mean field performs on average better ?  It seems quite logical for GFP for instance where a wild type exists, however it is less intuitive for datasets without wild type.

Typo: the recall and precision have the same expression.
soundness: 4
strengths: The paper is clear, well-written, and aligns with well-established benchmarks in the field, such as CBAS (Brookes et al.). 
The model is supported by convergence analysis and an extensive set of well-handled experiments.
summary: The authors present a novel variational method for learning to sample from rarely observed events, aiming to minimize a distance between the distribution of interest, namely $p(x∣y>t)$, and its parametric variational counterpart $q(x|\phi)$. The problem is reformulated to leverage the “whole dataset,” not just rarely observed events, and is expressed as Equation (5), which comprises two terms: $log p (y>t∣x)$ and the negative KL divergence between$q(x|\phi)$ and $p(x)$. The authors' final proposal is to estimate $p(y>t∣x)$ using a parametric function instead of a simple PI estimate. The variational distribution is optimized by a REINFORCE gradient estimator.
weaknesses: While the model description is clear, the model comprise a parametric distribution $p(x|D_0)$ which might be the biggest model shortcoming originating from the model own formulation. 

Its major impact is that it reweights the gradient estimates of $q(x|\phi)$. Intuitively, how would that compare simply to the iterative strategy of Cbas ?

Meta-review: This paper considers an active search problem and presents a method to generate new designs of rare desired class under budget constraints. The target problem is of discrete and combinatorial nature. The proposed algorithm is based on variational inference. Theoretical analysis is conducted to understand its performance. All reviewers agree the paper contains interesting and potentially impactful results. One major concern is the lacking of comparison with strong baselines. One reviewer also points out the paper misses some important related work. Overall, this is a nice contribution to active search algorithms.

Decision: Accept (Poster)