Bayesian Optimization via Continual Variational Last Layer Training
Authors: Paul Brunzema
Mikkel Jordahn
John Willes
Sebastian Trimpe
Jasper Snoek
James Harrison
Venue: ICLR 2025 Spotlight (2025)
Keywords: Bayesian deep learning
bayesian optimization
uncertainty
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We develop an efficient and expressive Bayesian neural network surrogate for Bayesian optimization

Abstract: Gaussian Processes (GPs) are widely seen as the state-of-the-art surrogate models for Bayesian optimization (BO) due to their ability to model uncertainty and their performance on tasks where correlations are easily captured (such as those defined by Euclidean metrics) and their ability to be efficiently updated online. However, the performance of GPs depends on the choice of kernel, and kernel selection for complex correlation structures is often difficult or must be made bespoke. While Bayesian neural networks (BNNs) are a promising direction for higher capacity surrogate models, they have so far seen limited use due to poor performance on some problem types. In this paper, we propose an approach which shows competitive performance on many problem types, including some that BNNs typically struggle with. We build on variational Bayesian last layers (VBLLs), and connect training of these models to exact conditioning in GPs. We exploit this connection to develop an efficient online training algorithm that interleaves conditioning and optimization. Our findings suggest that VBLL networks significantly outperform GPs and other BNN architectures on tasks with complex input correlations, and match the performance of well-tuned GPs on established benchmark tasks.

Published as a conference paper at ICLR 2025
BAYESIAN OPTIMIZATION VIA CONTINUAL VARIATIONAL LAST LAYER TRAINING
Paul Brunzema', Mikkel Jordahn', John Willes?, Sebastian Trimpe',
Jasper Snoek*, James Harrison*
'RWTH Aachen University, Technical University of Denmark, Vector Institute,
*Google DeepMind
Contact: brunzema @dsme.
de, jamesharrison@google..com
ABSTRACT
Gaussian Processes (GPs) are widely seen as the state-of-the-art surrogate models for Bayesian optimization (BO) due their ability model uncertainty and their performance on tasks where correlations are easily captured (such as those defined by Euclidean metrics) and their ability to be efficiently updated online. However, the performance of GPs depends on the choice of kernel, and kernel selection for complex correlation structures often difficult or must be made bespoke. While Bayesian neural networks (BNNs) are promising direction for higher capacity surrogate models, they have so far seen limited use due to poor performance on some problem types. In this paper, propose approach which shows compet- itive performance on many problem types, including some that BNNs typically struggle with. We build on variational Bayesian last layers (VBLLs), and connect training of these models to exact conditioning in GPs. We exploit this connection develop an efficient online training algorithm that interleaves conditioning and optimization. Our findings suggest that VBLL networks significantly outperform GPs and other BNN architectures on tasks with complex input correlations, and match the performance of well-tuned GPs on established benchmark tasks.
INTRODUCTION
Bayesian optimization (BO) has become an immensely popular method for optimizing black-box functions that are expensive to evaluate and has seen large success in a variety of applications (Garnett
et al., 2010; Snock et al., 2012: Calandra et al., 2016; Marco et al., 2016; Frazier & Wang, 2016; Berkenkamp et al., 2017: Chen et al., 2018; Neumann-Brosig et al., 2019; Griffiths & Hernández- Lobato, 2020; Colliandre& & Muller, 2023). In BO, the goal is to optimize some black-box objective f: X' > R K (where X R²) in few samples possible whilst only having access to sequentially sampled, potentially noisy, data points from possibly multiple objectives.
Gaussian processes (GPs) have long been the de facto surrogate models in BO due to their well- calibrated uncertainty quantification and strong performance in small-data regimes. However, their application becomes challenging high-dimensional, non-stationary, and structured data environ- ments such as drug-discovery (Colliandre & Muller, 2023; Griffiths & Hernández-Lobato, 2020) and materials science (Frazier & Wang, 2016). Here, often prohibitively expensive bespoke kernels are necessary to capture meaningful correlations between data points. Furthermore, the scaling of GPs large datasets typically associated with high- dimenssinlii spaces can be limiting- especially limiting-expeciilly combined with online hyperparameter estimation. To address these challenges, integrating Bayesian Neural Networks (BNNs) into BO as alternative surrogate models has gained increasing attention (Snock et al., 2015; Springenberg et al., 2016; Foldager et al., 2023; Li et al., 2024). While BNNs inherently scale with data, challenges like efficient conditioning on new data and consistency across tasks persist. Moreover, the performance of BNNs many tasks has so far not matched the perfor- mance of GPs due to issues like under-fitting (Li et al., 2024; Wenzel et al., 2020; Ovadia et al., 2019; Izmailov et al., 2021).
In this work, we develop an approach for surrogate modeling in Bayesian optimization that combines the strengths of GPs (efficient conditioning, strong predictive performance, simple and effective

---

Published as a conference paper at ICLR 2025
True objective * Observations
VBLL posterior predictive
VBLL posterior samples
BO Step
ВО Step 2
ВО Step
-2
0.0
0.5
1.00.0
0.5 Parameter x
1.0
0.5
1.0
Parameter X'
Parameter X
Figure I: Variational Bayesian last layer model as surrogate model for BO on a toy example. The VBLL model can capture in-between uncertainty and analytic posterior samples are easily obtained through its parametric form making it suitable surrogate for BO.
uncertainty quantification) with the strengths of BNNs (scalability and ability handle highly non- Euclidean correlation structures). Our approach builds on Variational Bayesian Last Layer (VBLL) neural networks (Harrison et al., 2024), which are similar to parameteric GPs with learned kernel.
Contributions: Our main contributions and findings are:
We show VBLL surrogate model outperform wide variety of baselines on diverse set of problems, including those with discrete inputs and multi-bbjective problems.
We present connection between model conditioning (in Bayesian sense) and model optimiza- tion, and develop an efficient online optimization scheme that exploits this connection.
We experimentally connect the hyperparameters of the VBLL model to problem features, such as stochasticity, and discuss implications for other BNN models.
RELATED WORK AND BACKGROUND
Various forms of Bayesian or partially-Bayesiie neural networks have been explored for BO, including mean field BNNs (Foldager et al., 2023), networks trained via stochastic gradient Hamiltonian Monte al.. 2016; Kim et al., 2021), and last layer Laplace approximation BNNs (Kristiadi et al., 2021; Daxherger et al., 2021). Li et al. (2024) find that infinite-width BNNs (I-BNNs) (Lee et al., 2017; Adlam et al., 2021: 2023) perform particularly well especially on high-dimensional. non-stationary and -on-Euclidean problems, setting for which standard GP priors are inappropriate or not well specified, and for which designing suitable kernels is difficult.
Carlo (Springenberg
While BNNs are promising, they have often proven to be challenging to train and complex to use in practice. Bayesian last layer networks- provide simple (and often much casier to train) Purtially-Baessiad neural network model (Snock et al., 2015; Azizzadenesheli et al., 2018; Harrison et al., 2018; Weber et al., 2018; Riquelme et al., 2018; Fiedler & Lucia, 2023). Concretely, the standard model for regression with Bayesian last layer networks and a one-dimensional! output is
consider uncertainty only over the output layer-
20 + +,,
where w 1.th the last layer for which Bayesian inference is performed, and de are features learned by neural network backbone with (point estimated) parameters 0. The noise ~ o?) assumed be independent and identically distributed.
With this observation model, fixed features Φg. and Gaussian prior on the weights as p(w) = A (W. Sp). posterior inference for the weights analytically tractable via Bayesian linear regression. In particular, Bayesian linear regression is possible via a set of recursive updates in the natural
For multi-variate modeling. we assume each element of is independent, and so the N-dimensional problem separable into N independent inference problems with shared features Φb. Relaxing this assumption is possible and leads matrix normal-distribuedd posterior.
2

---

Published as a conference paper at ICLR 2025
parameters
q: 0 0 Dilt + 9t-1
S. " +
where w, Siq. is the vector of precision-mean and Φi == Φ (x)). For dataset Dr
{(,,, 31)))=111111 this recursive update yields posterior p(w Dr) = -V ST) and posterior predic- tive
ply x. Dr, 0) (o
(x (r) 0²)
Since the predictive distribution is Gaussian, it pairs ftraight-forwardly with conventional acquisition functions Bayesian optimization and bandit tasks (Snoek et al., 2015; Riquelme et al., 2018).
Usually, such BLL models are trained using gradient descent on the exact (log) marginal likelihood over all data points. This can be computed either via the marginal likelihood or by backpropagating through the recursive last layer update (Harrison et al., 2018). The cost of this is prohibitive, as it requires iterating over the full dataset for each gradient computation. Moreover, it often numerically unstable and can result in pathological behavior of the learned features (Thakur et al., 2020; Ober et al., 2021). 1). An alternative strategy to train on mini-batches (Snoek et al., 2015) to learn features and condition the last layer the full dataset after training. However, this yields biased gradients and often results in a severely over-comcentrated final model (Ober & Rasmussen, 2019).
To increase efficiency, recent work (Harrison al., 2024; Watson et al., 2021) developed determin- istic variational lower bound to the exact marginal likelihhod and proposed to optimize this instead, resulting in the variational Bayesian last layer (VBLL.) model. In our setup, fix variational posterior ((w) = NN S). The lack of subscript denotes that the last layer parameters belong to the variational posterior, and we will write n (w, S) for convenience. Following Harrison al. (2024, Theorem 1), the variational lower bound for regression with BLLs (under the prior defined previously) is
logp(Y Xx,0 2
Σ
1 Log.N(sn wT... o. 0²) - -0-2--0 S00002) (L(g,(w) || ((w)) =: C(n, 0).
The variational posterior over the last layer is trained with the network weights 0 in standard neural network training loop. yielding lightweight Bayesian formulation. We refer to Appendix A for how VBLLs relate to other BNN methods we benchmark against in this work.
3 TRAINING AND MODEL SPECIFICATION
In this section, we discuss the training procedure we use for VBLLs in the BO loop. We first identify a relationship between variational training Bayesian last layer models (as used in the VBLL class of models) and recursive last layer computation (as is used in standard BLL models). We further present different methods for choosing whether to perform a recursive last layer update or full model retraining for newly-observed datapoint. We then describe the feature training procedure, and diseuss early stopping and previous model re-use. Finally, we present the overall training loop that combines these two training approaches, and relevant hyperparameters. We expand on all design decisions in the Appendix.
3.1 VARIATIONAL INFERENCE AND RECURSIVE LAST LAYER UPDATING
In this subsection, we describe connection between the VBLL variational objective (5)-which is optimized through standard mini-bated optimization-ann the recursive updating associated with Bayesian linear regression and BLL models.
Theorem 1. Fix O. Then, the variational posterior parameterized by
"in" : 17" % arg max L(n, 0)
equivalent to the posterior computed by the recursive least squares inferential procedure described by (2) and (3), iterated over the full dataset.
3

---

Published as a conference paper at ICLR 2025
This result follows from the fact that the true posterior for Bayesian linear regression is contained in our chosen variational family, and that the variational posterior is tight. The full proof is provided in Appendix B.1.
This equivalence provides a bridge between two different interpretations of model training: optimization- based training (as is used in NNs) and conditioning (as is used in GPs). This equivalence allows us to consider an online optimization procedure that interleaves two steps: full model training (including feature weights 0) on the variational objective, and last layer-only training via recursive conditioning. While the former optimizes all parameters and the latter only optimizes the last layer parameters, they are optimizing the same objective which stabilizes the interleaving of these steps.
We parameterize a dense precision (in contrast to standard VBLL models which parameterized the covariance) for the variational posterior by Cholesky decomposition, S-1 LL (with L lower triangular). The recursive precision update is then
ب = Φ ---1L
The updated Cholesky factor Lt can be computed efficiently via rank-1 update that preserves trian- gularity (Krause & Igel, 2015). The mean update can be computed via the natural parameterization
qi 0iye + qt-1
(8)
from which the last layer mean be computed. We do not store the computation graph associated with these updates and do not backpropagate through model updates. This update has quadratic complexity in the dimensionality of Φ, and enables quadratic complexity for (efficient, randomized) computation of all terms in training. We discuss the numerical complexity of this procedure and the overall training loop in Appendix B.5.
3.2 FULL MODEL TRAINING
In full model training, directly train the neural network weights 0, the variational posterior parameters y7, and the MAP estimate of the noise covariance o² via gradient descent on (5).
To improve training efficiency, we perform early stopping based on the training loss. Whereas standard neural network training will typically do early stopping based on validation loss, we directly terminate training when training loss does not improve for more than set number of epochs. Whereas standard neural networks would typically substantially overfit with this procedure, we find that VBLL networks are not prone severe overfitting. We additionally experiment with a more advanced form of continual learning in which we initialize network weights with previous training parameters. However, find that the benefits from this are relatively minor and there often exploration advantages to re-initializing the network weights at the start of training. Details of early stopping and model re-use are provided in Appendix B.2.
Training efficiency:
Optimization: For each iteration of model training, we will re-initialize the optimizer state. Because are interleaving recursive updates, continued use of previous optimizer states such momentum and scale normalization (as in e.g. Adam (Kingma & Ba, 2015)) actively harm performance. However, we do find benefits in using adaptive optimizers for training in general and use AdamW (Loshchilov & Hutter, 2017) throughout our experiments (with no regularization on the variational last layer). All training details are described Appendix B.4.
3.3 CONTINUAL LEARNING TRAINING LOOP AND HYPERPARAMETERS
To alleviate VBLL surrogate fit times, propose continual learning training loop (Algorithm 1) which interleaves two training steps: full model training via minibatch gradient descent, and last layer conditioning.
Choosing when to re-train model features: There are several possible methods for choosing (for each new data point) whether to perform full model re-training or a last layer recursive update. We consider three methods in this paper. First, we consider a simple scheme re-training the network every M 2 steps, and otherwise doing recursive updates. Practically this performs poorly, and allocating compute to model training earlier in the episode performs substantially better. Therefore, also consider scheduling scheme in which re-training randomized, with the probability of
4

---

Published as a conference paper at ICLR 2025
Algorithm 1 VBLL Bayesian Optimization Loop with Continual Variational Last Layer Training
Require: Model train frequency model and total number of evaluations T: Wishart prior scale V
1: D <- Din obtained from evaluations at pscudo-random input locations from a Sobol sequence
2: for 0 to T do
3:
Preinit <- RE-INIT CRITERION(D, > Periodic, scheduled, or event-triggent initialization
4: if Treinit then
5:
Init. model parameters 7, 0, V and train via minibatch optimization on (5) using D
6:
else
7:
Update 77 via (rank-1 Cholesky) recursive updates (2) and (3) using (II-1, 3t-1) and Σ
8: end
9: Select at, via optimizing acquisition function; query objective function with œy and receive y/
10: D - DU ((x,)))))
end for
re-training declining via sigmoid curve over the course of the episode. Lastly, we consider loss- triggered method for re-training, in which we re-train the full model if the log predictive likelihood of the data under the model is below some threshold, and otherwise perform recursive updates. In the main body of the paper, we will only present the latter approach when performing continual learning. Details re-initialization schemes and results for other scheduling methods are provided in Appendix B.3.
Hyperparameters: The VBLL models have several hyperparameters, some of which similar to GPs and some of which are substantially different. We include hyperparameter studies in Appendix We specifically investigate the hyperparameter sensitivity of the noise covariance prior, network reinitialization rate and neural network architecture. Overall we find that the VBLL models are highly robust these hyperparameters, but that proper tuning of them for the problem at hand can further improve performance.
4
ACQUISITION FUNCTIONS
In this section, we describe acquisition functions for our VBLL models in both the single objective and multi-objective case.
4.1 SINGLE OBJECTIVE ACQUISITION FUNCTIONS
VBLLs yield a Gaussian predictive distribution and thus most acquisition functions that are straight- forward to compute for GPs are also straightforward for VBLLs. However, parametric VBLLs are also especially well suited for Thompson sampling compared to non- con-paraettric models like GPs?. For Thompson sampling, we simply sample from the variational posteríor of w iteration and then construct a sample from the predictive (cf. Fig. 1) as a generalized linear model as
i
on w))
2
f(x) = ib Φo(x))
3
E++1 arg max f(x) #EX
(9)
This sample of the predictive can then be optimized analytically, which differs from classic Thompson sampling methods used for non non-paramen paraamemtriic GPs. This can similarly be done for LLLA BNNs (cf. Sec. 5), although Laplace approximated posterior samples are known to suffer from underfitting on training data (Lawrence, 2001; Roy et al., 2024). For the analytic optimization of the parametric sample, we use L-BFGS-B (7hu et al., 1997) leveraging the fact that we can also easily obtain the gradient of f. We initialize the optimization at multiple random initial conditions and choose the best argmax as the next query location.
TThompson sampling for GPs often involves drawing samples from gigh-dimensionll posterior distributions generated at pseudo-randon input locations (e.g.. using Sobol sequences) and then selecting the argmax of the discrete samples the next query locations (Eriksson et al., 2019). It is worth noting that while it is possible construct analytic approximate posterior samples for GPs (Wilson al., 2020; 2021), this approach not yet. all kerness.
commonly adopted current practice, and is not possible to
In the following, the number of random initial conditions will be set the same number of random initial conditions for the standard optimization of the acquisition functions for fair comparison.
5

---

Published as a conference paper at ICLR 2025
True Pareto front
Pareto front of T'hompson sample P
Training points
Reference point
Index for next point
Iteraton
100
Objective (x)
Tteration 7
Tteration 15
Tteration 23
Tteration 100
75
-10
25
Current HV: 0.00
Current HV: 36.99
Current HV: 54.50
Final HV; 58.36
-15
-100
100
0
100
0
- 1100
Objective fi (wx
Objective fi (x)
Objective (x)
Objective fx
Figure 2: Multi-objective Thompson sampling on BraninCurrin. At each iteration, we optimize the multi-objective Thompson sample and choose as the next point the index that increases the predicted hypervolume the most. At the end of the optimization (right with colormap), we can observe that the true Pareto front is nicely approximated.
4.2 MULTI-OBIECTIVE ACQUISITION FUNCTIONS
In multi-objective BO it is common model each objective with separate GP (Zitzler al., 2003; Daulton et al., 2020; Ament et al., 2024). While it is often preferable to model correlations between outputs (Swersky et al., 2013), this requires correlation kernels for the objectives, which may be difficult to specify. With neural network based surrogate models, one can simply set the number of outputs to the number of objectives, thus sharing feature learning between objectives. For VBLL networks, we can use the multivariate regression formulation to obtain multivariate normal prediction for each point. Since this again is Gaussian, it is straightforward use popular acquisition functions for multi-objective optimization such as expected hypervolume improvement (Zitzler et al., 2003: Daulton et al., 2020; Ament al., 2024).
Also for multi-objective optimization problems, we can leverage the parametric form of the VBLLS to do efficient Thompson sampling. It also possible to do Thompson sampling with GPs for multi-objective optimization, but this usually involves approximations with high number of Fourier bases functions (Bradford et al., 2018; Belakaria et al., 2019; Paria et al., 2020). For Thompson sampling with VBLLs, we first sample NN realization of the VBLL network as in (9) and optimize the sample with NSGA-II (Deb et al., 2002) using Pymoo (Blank & Deb, 2020). This yields a predicted Pareto front P == (91,.... ýr) with size equal to the population size P of the NSGA-II solver. In summary, the Thompson sampling procedure for multi-objective optimization is
W 2 (W)
2
f(x)
W (x()))
P max f(x). SE
(10)
To choose the next candidate, greedily select the index of the point in the predicted Pareto front that maximizes the improvement of the hypervolume. The hypervolume (HV), denoted by HVPP, r). a commonly used metric multi-objective optimization, which measures the volume of the region dominated by the Pareto front P with respect to predefined reference point r in the objective space. Our goal in multi-objective BO is to expand this dominated region by adding new candidate to the Pareto front. At each iteration k the current Pareto front is denoted by Pt. To select the next candidate, we maximize the hypervolume improvement, which is the increase in hypervolume when adding a predicted new point ij P to the existing Pareto front Px as
#k++ E arg max HV(PP r)). gEP
(11)
no points in the predicted Pareto front improve the hypervolume or the solution non-unique, we randomly sample from the set of maximizers. Figure 2 shows Thompson sampling with VBLLS on the BraninCurrin benchmark (cf. Sec. 5). After the initialization, none of the points in the Pareto front P improve the HV. Therefore random index from the list of solutions in chosen for the next query location. After some iteration, the P partially includes the true Pareto front. The algorithm continuously samples points that improve the HV until after 100 iterations P is well approximated.

---

Published as a conference paper at ICLR 2025
In all experiments in this paper, we assume the reference point to be known but we should note that it is also possible to estimate or adaptively choose the reference point т, which can improve the effectiveness of the optimization process (Bradford al., 2018). Additionally, this strategy of maximizing hypervolume can be combined with other acquisition functions; of particular interest are information-theoree acquisition functions in the single objective case. For instance, max-value entropy search for multi-objective optimization (Belakaria et al., 2019) relies on random Fourier features approximate samples from the GP posterior; this can easily be switched using Thompson sampling with VBLL networks.
RESULTS
We evaluate the performance of the VBLL surrogate model on various standard benchmarks and three more complex optimization problems, where the optimization landscape is non-stationary. For experimental details and ablations of hyperparameters, we refer the reader to Appendix and E.
5.1 SURROGATE MODELS
We use the following baselines throughout the experiments in addition to VBLLs: Gaussian Processes (GPs), Infinite- -iith Bayesian Neural Networks (I-BNNs) (Lee al., 2018: Li al., 2024), Deep Kernel Learning (DKL) (Wilson et al., 2016), Last Layer Laplace Approximations (LLLA) (Daxberger et al., 2021; Kristiadi et al.), and Deep Ensembles (DE). For the configuration of neural network backbone of the BNN baselines as well as the VBLLLs, we closely follow the setup in Li et al. (2024). For details all baselines and their configurations, refer Appendix D.1.
In all subsequent experiments, we select the number of initial points for the single objective bench- marks equal to the imput dimensionality D and for the multi-objective benchmarks we use 2 (D + 1) initial points (Daulton et al., 2020; Balandat et al., 2020). We further use batch size of one as in the classic BO setting.* We compare the performance of all surrogates for the following acquisition functions: (i) log expected improvement (IogEI) (Ament et al., 2024), numerically more stable version of standard expected improvement, (ii) Thompson sampling (TS) (Thompson, 1933; Russo et al., 2018), (ili) and log expected hypervolume improvement 1oTEHVI (Ament et al., 2024) for the multi-objective benchmarks.
5.2 BENCHMARK PROBLEMS
We begin by examining a set of standard benchmark problems commonly used to assess the perfor- mance of BO algorithms (Eriksson al., 2019; Ament et al., 2024). Figure 3 (top) illustrates the performance of all surrogates on these benchmark problems. It can be observed that, as expected, GPs perform well. The BNN baselines also demonstrate strong performance on lower-dimensional problems, although they do not match the performance of GPs on the Hartmann function. VBLLs, however, generally perform on par or better than most other BNN surrogates on the classic bench- marks. Interestingly, for TS, we notice that on the Ackley20 2D and AckleySD benchmark, the VBLLs with analytic optimization of the Thompson samples even surpass the performance of GPs as well as all other baselines. The continual learning baselines show slightly worse performance than the standard VBLLs but at the benefit of reduced computing.
5.3
AND
PROBLEMS
GPs without tailored kernels often struggle in high-dimensional and non-stationary environments (Li al., 2024); areas where deep learning approaches are expected to excel. Our results on the 200.D NNdraw benchmark (Li et al., 2024), the real-world 25D cont rol benchmark (Oh et al., 2019), and the 12D Lunarl ander benchmark (Eriksson al., 2019) are shown in Figure 3 (bot- tom). On these benchmarks, all BNN surrogates show strong performance; especially LLLA, DEs and the VBLL networks. While GPs perform well with LogEI on NNdras aw and the I-BNNs show good performance on restcontrol.. the VBLLs and LLLA models are consistently the best-performing surrogates for TS. Similar the classic benchmarks, the continual learning version of the VBLLS
Note that this differs from the experimental setup by Li (2024) which use larger hatch sizes.
7

---

Published as a conference paper at ICLR 2025
LBNN
DKL
LLLA
DE
VBLL
VBLL (ET CL)
Best value
logEl
-10
20
10
Branin (2D)
Ackley (2D)
Ackley (5D)
Hartmann (6D)
30
Best value
10
TS
20
Branin (2D)
10
Ackley (2D)
10
Ackley (6D)
Hartmann (6D)
-30
20
40
20
40
50
100
50
100
Iteration
Iteration
Iteration
Iteration
Best value
7500
12.5
logel
200
5000
15.0
2500
NN draw (200D).
Pest; control (25D)
Lunar lander (12D).
Best value
7500
12.5
200
5000
15.0
2500
(2008)
Pest control (25D)
Lunar lander (12D)
200
400
600
50
100
200
Iteration
Iteration
Iteration
Figure 3: Classic benchmarks (top) and high-dimensionel and non- statiiati benchmarks (bottom). Performance of all surrogates for logEI (top) and TS (bottom).
shows similar performance to the VBLLs. Additionally, we note that by adjusting the Wishart scale or changing the re-initialization strategy we can further enhance the performance of continual learning for LogEI (cf. Appendix E.2 and Appendix B.3.3, respectively).
5.4 IMPACT OF RE-INITIALIZATION METHODS
VBLL
VBLL (ET CL)
VBLL (SG CL)
We compare our two strategies for determin- ing when to re-initialize the model: event- triggered and scheduled re-initialization (ef. Ap- pendix B.3). For this, we define the indica- tor function I(t), which is if the model is re-initialized at iteration t and 0 if updated re- cursively. We show the expected value of this
TS Best value
-10
-20
Branin (2D)
-10
Ackley (5D)
-30
ET: -45 RT
ET: -6% RT
Indicator function to releara model
SG: -60.73% RT
SG: -67.35% RT
indicator function as well as the performance (top) ten seeds (bottom) in Fig. 4. We can observe that both strategies perform similarly for TS, achieving Rear-VBLL performance with large reductions in runtime as highlighted by the percentage decreases. Further, the event trigger
1.0
1.0
0.5
0.5
0.0
0.0
20
40
50
100
adapts more flexibly than the pre-defined sig- moid schedule, particularly for Branin, where it reduces runtime early on. In Appendix B.3.3, we included longer discussion on the compar- ison of the re-in i--itialization strategies. We find that the event- trint-triggeeed approach provides greater adaptivity compared to the sigmoid schedule but can lead to stagnation in high-dimensional settings, as observed for NNdraw with logET in Fig. 3.
Iteration
Iteration
Figure 4: Impact of event-triggered (ET) and sig- moid scheduled (SG) re-initialization the per- formance and runtime (RT) saving.
8

---

Published as a conference paper at ICLR 2025
GP
L-BNN
DKL
LLLA
VBLL
VBLL+TS
Hypervolume diff. log (HV HV)
BraninCurrin (2D, 2K)
DTL21 (6D, 2K)
DTLZ2 (5D, 2K)
10
Dil Sorbent
10
(7D, 3K)
10
10-1
10a
10
100
10a
103
1.6 x10
X
Numerical issues of logEHVI
x10
Hypervolume HV
60
0.4
40
0.3
1.0
1.5
(2D,
2K)
(6D,
2K)
(6D, 2K)
0.5
(7D, 3K)
BraninCurrin
DTLZ1
0.1
DTL22-
0i1
Sorbent
50
100
50
100
100
100
200
300
Iteration
Tteration
Tteration
Iteration
Figure 5: Multi-objectivve benchmarks. Performance of all surrogate models using logENVI and VBLLs with TS. A cross indicates the crash surrogate's s furthest seed due to a numerically unstable acquisition function. VBLL+TS successfully navigates the numerically unstable HV plateau in OilSorbent, enabling a more accurate approximation its three-dimensional Pareto front.
GP
T-BNN
DKL
LLLA
VBLL
DE
VBLL (ET+CL)
12
TS
7500
logEf
Best value
5000
-20
-16
2500
Bramin (2D)
-10
ey (5D)
BERLOOUT, (25D)
2DD au (200D)
-30
10
10a
10a
10
10a
10a
10
Total fit time
[s]
Total fit time
[s]
Total fit time [s]
Total fit time [s]
Figure 6: Performance vs. accumulated surrogate fit time. VBLLs are the most expensive train
surrogate model. Using the proposed continual learning scheme can significantly reduce runtime while maintaining good performance.
Interestingly, this stagnation does not early when using TS, suggesting that the additional randomness helps mitigate this issue.
5.5 MULTT-OBJECTIVE PROBLEMS
We evaluate our approach on multi objectib optimization problems. Here, we consider the standard
benchmarks BraninCurrin (D 2, K 2), DTLZI (D = 5, K = 2). DTLZ2 (D = 5. K 2), and the real world benchmark Oil Sorbent (D = 7, K 3) (Wang et al., 2020; Li et al., 2024). Figure shows the performance of all surrogate in terms of the obtained HV with respect to fixed reference points, as well as the logarithmic HV difference between the maximum HVS and the obtained HV as common metries in multi-objective BO (Belakaria al., 2019; Daulton al., 2020). The performance of all BNN models is similar. It is however notable that OilSorbent all baselines using 1 LogEHVI crash due to numerical issues while optimizing the acquisition function. Such numerical issues are known for expected improvement type acquisition functions due to vanishing gradients (Ament et al., 2024). Thompson sampling with VBLLs does not result numerical issues, we observe that with this combination of model and acquisition function, the Pareto front can be further refined.
For BraninCurrin, DTLZ71, and DTT.72 we the values provided by BoTorch, and for OilSorbent estimate the maximum HV based the asymptotic performance of the best performing surrogate.

---

Published as a conference paper at ICLR 2025
5.6 TRAINING TIME COMPARISON
Lastly, we consider surrogate fit time comparison in which we keep constant BO budget and track accumulated fit time to highlight that recursive updates can significantly speed up run time. It should be noted that in the usual BO setting, we assume that the black-box function is expensive to evaluate, and hence, fit times play only a minor role in many applications, such as drug- or materials discovery. Still, there may be applications where fit times are important. In Figure 6, observe that VBLLS are the most expensive surrogates to train on smaller problems, while on NNDraw, DKL and DE the most expensive. The event-triggered re-initializing combined with recursive updates of the variational posterior significantly reduces the runtime while maintaining good performance.
DISCUSSION AND CONCLUSIONS
In this paper, we have developed bridge between classical conjugate Bayesian methods for BO and BNN-based methods. We built upon VBLLs (Harrison et al., 2024) to obtain the expressivity of neural networks, together with the simple and robust uncertainty characterization associated with GPs.
Continual learning: We introduced an online training scheme that enables online complexity comparable to GPs. While we have developed an efficient training scheme, further developments in combining continual learning with our architecture are possible. For our experiments, we use a simple threshold the log- legliikli from model selection can likely be used to make even more informed decisions. While our approach was effective, many methods in continual learning (e.g. Nguyen et al. (2018)) have been developed to improve training efficiency. While the combination of seeviously-developee schemes with our last layer updating may be non-trivial, this is likely promising direction for improving efficiency. Balancing the functional randomness afforded by network re- re-initializaaion with efficient continual learning will likely be necessary focus of future work this direction.
decide re-initializing the model. Here, further concepts
continual learning
Additionally, as VBLLs are new model class, we believe that there exist various options beyond the approaches developed herein that can be used to speed up runtime further. our experiments, we used the default initialization provided in Harrison et al. (2024). Further investigating the initialization of the VBLLS through form of stochastic feature re-use and its effect on fit time can likely improve wall clock performance. Also, standard deep learning best practices, such learning rate scheduling, could be beneficial, and exploring these directions is promising if wall clock time is important.
VBLLS perform on par with GPs on standard low-dimensional bench- marks, yet significantly outperform GPs in high-dimensional and non-stationary problems. Fur- thermore, VBLLs outperform wide array of Bayesian neural network models, including I-BNNs (Adlam et al., 2021), Laplace approximation-bon (Hansen & Salamon, 1990; Lakshminarayanan et al., 2017) and deep kernel learning (Wilson et al., 2016). Especially the combination of VBLLs with TS shows strong performance throughout. While laplace-approximation methods perform well in our experiments, find that they appear more sensitive to observation noise (see Appendix C.3 for experimental results highlighting this).
Performance and baselines:
models (Daxberger et al., 2021), 1), deep ensembles
Because of the parametric architecture of VBLLs, they can efficiently be combined with Thompson sampling as an acquisition function. This yields reasonable benefits the univariate case (as shown by our full set of results) but yields substantial benefits in the multi-objective case due to better numerical conditioning. Moreover, the ease of use with Thompson sampling also enable easy use of TS-based ieformation-thoorett acquisition functions such as max-value entropy search (Wang& & Jegelka, 2017; ; Belakaria et al., 2019).
Acquistion functions:
ACKNOWLEDGMENTS
The authors thank Friedrich Solowjow, Alexander von Rohr, Kevin Murphy, and Zi Wang for their helpful comments and discussions. Paul Brunzema is partially funded by the Deutsche Forschungs- gemeinschaft (DFG, German Research Foundation)-RTT 2236/2 (UnRAVeL). Simulations were performed in part with computing resources granted by RWTH Aachen University under projects wwth1579 and p0022034.
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Ben Adlam, Jaehoon Lee, Lechao Xiao, Jeffrey Pennington, and Jasper Snoek. Exploring the uncertainty properties of neural networks' implicit priors in the infinite-width limit. International Conference on Learning Representations (ICLR), 2021.
Ben Adlam, Jaehoon Lee, Shreyas Padhy, Zachary Nado, and Jasper Snoek. Kernel regression with infinite-width neural networks on millions of examples. arXiv:2303.05420 2023.
Sebastian Ament, Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bakshy. Unex- pected improvements to expected improvement for Bayesian optimization. Neural Information Processing Systems (NeurIPS), 36, 2024.
Kamyar Azizzadenesheli, Emma Brunskill, and Animashree Anandkumar. Efficient exploration through Bayesian deep q q-nttworks. arXiv:1802.00 1802.04412, 2018.
Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, An- drew Gordon Wilson, and Eytan Bakshy. BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization. In Neural Information Processing Systems (NeurIPS), 2020.
Syrine Belakaria, Aryan Deshwal, and Janardhan Rao Doppa. Max-value entropy search for multi- objective Bayesian optimization. Neural Information Processing Systems (NeurIPS), 32, 2019.
Felix Berkenkamp, Matteo Turchetta, Angela Schoellig, and Andreas Krause. Safe model-based reinforcement learning with stability guarantees. Neural Information Processing Systems (NeurIPS), 2017.
Julian Blank and Kalyanmoy Deb. Pymoo: Multi-objective optimization in python. IEEE Access, 89497-89509, 2020.
Eric Bradford, Artur M Schweidtmann, and Alexei Lapkin. Efficient multiobjective optimization employing Gaussian processes, spectral sampling and a genetic algorithm. Journal of Global Optimization, (1(2):407-438, 2018.
Roberto Calandra, André Seyfarth, Jan Peters, and Mare Peter Deisenroth. Bayesian optimization for learning gaits under uncertainty: An experimental comparison on dynamic bipedal walker. Annals of Mathematics and Artificial Intelligence. 76:5-23, 2016.
Yutian Chen, Aja Huang, Ziyu Wang, Ioannis Antonoglou, Julian Schrittwieser, David Silver, and Nando Freitas. Bayesian optimization alphago. arXiv: 06855,, 2018.
Lionel Colliandre and Christophe Muller. Bayesian optimization in drug discovery. High Performance Computing for Drug Discovery and Biomedicine, 2023.
Samuel Daulton, Maximilian Balandat, and Eytan Bakshy. Differentiable expected hypervolume
improvement for parallel multi-oojectiive Bayesian optimization. Neural Information Processing Systems (NeurIPS), 33:9851-9864, 2020.
Erik Daxberger, Agustinus Kristiadi, Alexander Immer, Runa Eschenhagen, Matthias Bauer, and Philipp Hennig. Laplace refux-effortless Bayesian deep learning. Neural Information Processing Systems (NeurIPS), 2021.
Kalyanmoy Deb, Amrit Pratap, Sameer Agarwal, and TAMT Meyarivan. A fast and elitist multi- objective genetic algorithm: NSGA-IL. IEEE Transactions on Evolutionary Computation, 6(2): 182-197, 2002.
David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization. Neural Information Processing Systems (NeurIPS), 2019.
Felix Fiedler and Sergio Lucia. Improved uncertainty quantification for neural networks with Bayesian last layer. JEEE Access, 2023.
Jonathan Foldager, Mikkel Jordahn, Lars Kai Hansen, and Michael Riis Andersen. On the role of model uncertainties in Bayesian optimization. In Uncertainty in Artificial Intelligence (UAI), 2023.
11

---

Published as a conference paper at ICLR 2025
Peter Frazier and Jialei Wang. Bayesian optimization for materials design. Information science for materials discovery and design, 2016.
Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wilson. Gpytorch: Blackbox matrix-matrinn gaussian process inference with gpu acceleration. In Neural Information Processing Systems (NeurIPS), 2018.
Roman Garnett, Michael A Osborne, and Stephen J Roberts. Bayesian optimization for sensor set selection. International conference on information processing in sensor networks, 2010.
Ryan-Rhys Griffiths and José Miguel Hernández-Lobato, Constrained Bayesian optimization for automatic chemical design using variational autoencoders. Chemical science, 11(2)5777-586
2020.
Raia Hadsell, Dushyant Rao, Andrei A Rusu, and Razvan Pascanu. Embracing change: Continual learning in deep neural networks. Trends cognitive sciences, 2020.
L.K. Hansen and P. Salamon. Neural network ensembles. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(10):993-100, 1990. doi: 10.1109/34.5881.
James Harrison, Apoorva Sharma, and Marco Pavone. Meta-learning priors efficient online Bayesian regression. Workshop the Algorithmic Foundations of Robotics (WAFR), 2018.
James Harrison, John Willes, and Jasper Snoek. Variational Bayesian last layers. International Conference on Learning Representations (ICLR), 2024.
Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in craistiis-Simmlatioon and Computation, 1989.
Carl Hvarfner, Erik Orm Hellsten, and Luigi Nardi. Vanilla bayesian optimization performs great in high dimensions. In International Conference on Machine Learning (ICML), 2024.
Pavel Izmailov, Sharad Vikram, Matthew D Hoffman, and Andrew Gordon Gordon Wilson. What are bayesian neural network posteriors really like? In International Conference Machine Learning (ICML), 2021.
Samuel Kim, Peter Y Lu, Charlotte Loh, Jamie Smith, Jasper Snoek, and Marin Soljačić. Deep learning for Bayesian optimization of scientific problems with high-dimensional structure. orXiv:2104.11665 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. International Conference on Learning Representations (ICLR), 2015.
Jeremias Knoblauch, Jack Jewson, and Theodoros Damoulas. Generalized variational inference: Three arguments for deriving new posteriors. arXiv preprint arXiv: 0rXiv:/904.02063. 2019.
Oswin Krause and Christian Igel. A more efficient rank-one covariance matrix update for evolution strategies. In ACM Conference on Foundations of GGGetic Algorithms, 2015.
Agustinus Kristiadi, Alexander Immer, Runa Eschenhagen, and Vincent Fortuin. Promises and pitfalls of the linearized laplace in Bayesian optimization. In Fifth Symposium on Advances in Approximate Bayesian Inference.
Agustinus Kristiadi, Matthias Hein, and Philipp Hennig. Learnable uncertainty under laplace approximations. In Uncertainty in Artificial Intelligence (UAT), 2021.
Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. Neural Information Processing Systems (NeurIPS), 2017.
Neil David Lawrence. Variational inference in probabilistic models. PhD thesis, Citeseer, 2001.
Gilwoo Lee, Siddhartha S Srinivasa, and Matthew T Mason. GP-iLQG: Data-driven robust optimal control for uncertain nonlinear dynamical systems. arXiv:1705.05344 2017.
12

---

Published as a conference paper at ICLR 2025
Jachoon Lee, Yasaman Bahri, Roman Novak, Samuel S Schoenholz, Jeffrey Pennington, and Jascha Sohl-Dickstein. Deep neural networks as Gaussian processes. International Conference on Learning Representations (ICLR), 2018.
Yucen Lily Li, Tim G. J. Rudner, and Andrew Gordon Wilson. A study of Bayesian neural network surrogates for Bayesian optimization. In International Conference on Learning Representations (ICLR), 2024.
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. arXiv:1711.0510. 2017.
David JC MacKay. A practical bayesian framework for backpropagation networks. Neural Computa-
fion, 1992.
Alonso Marco, Philipp Hennig, Jeannette Bohg. Stefan Schaal, and Sebastian Trimpe. Automatic lqr tuning based on gaussian process global optimization. IEEE International Conference on Robotics and Automation (ICRA), 2016.
Raphael A Meyer, Cameron Musco, Christopher Musco, and David P Woodruff. Hutch++: Optimal stochastic trace estimation. In SIAM Symposium Simplicity in Algorithms (SOSA), 2021.
Matthias Neumann-Browe Brosig, Alonso Marco, Dieter Schwarzmann, and Sebastian Trimpe. Data- efficient autotuning with Bayesian optimization: An industrial control study. IEEE Transactions on Control Systems Technology. 28(3):730-740, 2019.
Cuong V Nguyen, Yingzhen Li, Thang D Bui, and Richard E Turner. Variational continual learning. In International Conference on Learning Representations (ICLR), 2018.
Sebastian W Ober and Carl Edward Rasmussen. Benchmarking the neural linear model for regression. arXiv preprint arXiv:1912.088. 1912.0441,, 2019.
Sebastian W Ober, Carl E Rasmussen, and Mark van der Wilk. The promises and pitfalls of deep kernel learning. In Uncertainty in Artificial Intelligence (UAI), 2021.
Changyong Oh, Jakub Tomezak, Efstratios Gavves, and Max Welling. Combinatorial bayesian optimization using the graph cartesian product. Neural Information Processing Systems (NeurIPS), 32, 2019.
Yaniv Ovadia, Emily Fertig, Jie Ren, Zachary Nado, Sebastian Nowozin, Joshua Dillon, Balaji Lakshminarayanan, and Jasper Snoek. Can you trust your model's uncertainty? evaluating predictive uncertainty under dataset shift. In Neural Information Processing Systems (NeurIPS), 2019.
Biswajit Paria, Kirthevasan Kandasamy, and Barnabás Póezos. flexible framework for multi- objective Bayesian optimization using random scalarizations. In Uncertainty in Artificial Intelli- gence, pp. 766-776. PMLR, 2020.
Carlos Riquelme, George Tucker, and Jasper Snoek. Deep Bayesian bandits showdown. In Interna- lional Conference on Learning Representations (ICLR), 2018.
Hrittik Roy, Marco Miani, Carl Henrik Ek, Philipp Hennig, Marvin Pförtner, Lukas Tatzel, and Spren Hauberg, Reparameterization invariance in approximate Bayesian inference, 2024.
Daniel Russo, Benjamin Van Roy, Abbas Kazerouni, Ian Osband, Zheng Wen, et al. A tutorial on Thompson sampling. Foundations and Trends in Machine Learning, 2018.
Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical Bayesian optimization of machine learning algorithms. Neural Information Processing Systems (NeurIPS), 2012.
Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary, Prabhat, and Ryan Adams. Scalable Bayesian optimization using deep neural networks. International Conference on Machine Learning (ICML), 2015.
Jost Tobias Springenberg, Aaron Klein, Stefan Falkner, and Frank Hutter. Bayesian optimization with robust bayesian neural networks. Neural Information Processing Systems (NeurIPS), 2016.
13

---

Published as a conference paper at ICLR 2025
Kevin Swersky, Jasper Snock, and Ryan P. Adams. Multi-task bayesian optimization. 2013.
Sujay Thakur, Cooper Lorsung, Yaniv Yacoby, Finale Doshi-Velez, and Weiwei Pan. Uncertainty- aware (una) bases for Bayesian regression using multi-headed auxiliary networks. arXiv preprint arXiv:2006. 1695, 2020.
William R Thompson. On the likelihood that one unknown probability exceeds another view of the evidence of two samples. Biometrika, 1933.
Boqian Wang, Jiacheng Cai, Chuangui Liu, Jian Yang, and Xianting Ding. Harnessing novel -eatiei--ttnnn---sssseeddd assisted evolutionary algorithm to co-optimize three characteristics an electro- spun oil sorbent. ACS Applied Materials & Interfaces, 12(38):42842 42849, 2020.
machine-learming -assisted
Zi Wang and Stefanie Jegelka. Max-value entropy search for efficient Bayesian optimization. In International Conference on Machine Learning (ICML), pp. 3627-3635. PMLR, 2017.
Joe Watson, Jihao Andreas Lin, Pascal Klink, Joni Pajarinen, and Jan Peters. Latent derivative Bayesian last layer networks. In Artificial Intelligence and Statistics (AISTATS), 2021.
Noah Weber, Janez Stare, Arpit Mittal, Roi Blanco, and Lluís Mårquez. Optimizing over a Bayesian last layer. In NeurlPS workshop Bayesian Deep Learning, 2018.
Florian Wenzel, Kevin Roth, Bastiaan S Veeling, Jakub Swiatkowski, Linh Tran, Stephan Mandt, Jasper Snoek, Tim Salimans, Rodolphe Jenatton, and Sebastian Nowozin. How good is the bayes posterior in deep neural networks really? In International Conference on Machine Learning (ICML), 2020.
Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning. In Artificial Intelligence and Statistics (AISTATS), 2016.
James Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Mare Deisenroth. Efficiently sampling functions from Gaussian process posteriors. In International Conference Machine Learning (ICML). PMLR, 2020.
James T Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter Deisenroth. Pathwise conditioning of Gaussian processes. Journal of Machine Learning Research, 22(105):1-47, 2021.
Friedemann Zenke, Ben Poole, and Surya Ganguli. Continual learning through synaptic intelligence. International Conference on Machine Learning (ICML), 2017.
Ciyou Zhu, Richard H Byrd, Peihuang Lu, and Jorge Nocedal. Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale software (TOMS), 23(4):550-560, 1997.
optimization. ACM Transactions mathematical
Eckart Zitzler, Lothar Thiele, Marco Laumanns, Carlos M Fonseca, and Viviane Grunert Da Fonseca. Performance assessment of multiobjective optimizers: An analysis and review. IEEE Transactions on Evolutionary Computation, 7(2):117-132, 2003.
14

---

Published as a conference paper at ICLR 2025
A Related Work Comparisons
16
A.1 VBLL Comparison to Last Layer Laplace Approximations
16
A.2 VBLL Comparison Deep Kernel Learning
16
Algorithmic Details
16
B.1
Proof of Theorem
16
B.2
Model Re-Use and Early Stopping
17
B.3
Recursive Update Scheduling.
18
B.3.1
Model re-initialization through event triggering
18
B.3.2
Model re-initialization through scheduling
18
B.3.3
Comparison of re-initialization strategies
19
B.4
Training Details
20
B.5
Computational Complexity
20
B.6
Further Details the Setup and Acquisition Functions
20
Discussion and Baseline Ablations
C.1 On Neural Network Thompson Sampling
21
C.2 Comparison to D-scaled Gaussian process priors
22
C.3 Comparison of Sensitivity to Noise for LLLA and VBLLs
23
Experimental Details
23
D.1 Baselines
................................
23
D.2 Experiments
24
Hyperparameter Sensitivity
25
E.1 Wishart Scale and Continual Learning Sensitivity
25
E.2 Influence of the Wishart Scale on the Performance of Continual Learning
25
E.3 Model Width Ablation
25
E.4 Model Performance in the Presence of Noise
26
15

---

Published as a conference paper at ICLR 2025
A RELATED WORK COMPARISONS
A.1 VBLL COMPARISON TO LAST LAYER LAPLACE APPROXIMATIONS
Laplace methods approximate the predictive covariance with the (inverse) Hessian. In the line of work on Laplace approximations for neural networks (e.g. Daxberger et al. (2021)), the Fisher information matrix-which consists of the outer- outer-product of network gradients, summed over all data is used as an approximation to the Hessian. For last layer Laplace in particular, this reduces to the outer product of features, equivalent covariance computation for Bayesian last layer methods. The primary difference is Daxberger et al. (2021) also typically optimize the hyperparameters (including the last layer prior and hyperparameters) via an empirical Bayes objective. Thus, the last layer Laplace method in this setting can be interpreted as a form of the (non- Variational) Bayesian last layer model. While we optimize the noise covariance we do not optimize the last layer prior, which is the primary difference between our VBLL-based approach and last layer Laplace approximations. Another key difference is that in the presented VBLL models, the variational posterior learned jointly with the features (or basis functions), which is not the case for last layer Laplace where the approximate posterior distribution is learned post-hoc fitting of the features.
A.2 VBLL COMPARISON TO DEEP KERNEL LEARNING
Deep kernel learning (DKL) (Wilson al., 2016) aims learn an input transform kernel with deep neural network to effectively capture non-stationary and complex correlations in the input space. Specifically, uses neural network ge as input wrapping to the kernel as ÁDKL = (z)) ge (z')). The kernel's 's hyperparameters and the neural network are trained simultaneously by optimizing the exact (log) marginal likelihood of the GP posterior on the full data set. With the learned kernel, inference is identical to standard non-parametric GPs through explicit conditioning on past data and therefore scales as ( This explicit conditioning can become especially problematic for larger datasets. Here, parametric models such VBLLs and LLLA, which are trained mini-batches are preferable. Intuitively, VBLLs can be considered a parametric version of DKL, also learning correlations by optimizing the variational posterior of the last layer in tandem with the feature- extracting backbone network. Optimizing DKL is challenging primarily because of the computational complexity of gradient computation with marginal likelihood objective and overparameterizeed models (Ober et al., 2021). However, it can still yield good performance on many benchmarks which also demonstrate in Figure 3. Still, similar GPs, it struggles in high-damensional Thompson sampling due to its non-parametric nature.
B ALGORITHMIC DETAILS
B.1 PROOF OF THEOREM I
We restate the theorem for completeness.
Theorem 1. Fix O. Then, the variational posterior parameterized by
(w"
arg maxx(mp
(12)
is equivalent the posterior computed by the recursive least squares inferential procedure described by (2) and (3), iterated over the full dataset.
Proof. There are two methods to prove this result. First, we can note that for the standard ELBO constructed from single application of Jensen's inequality, exchanging the log and the expectation the variational lower bound is tight if the chosen variational family contains the true posterior (see e.g. Knoblauch et al. (2019) for discussion). Our chosen variational posterior formulation consists of independent posteriors for each output dimension, with dense covariances. This contains the exact posterior for multivariate Bayesian linear regression with diagonal Σ and thus the variational lower bound is tight under our modeling assumptions. Moreover, because the variational lower bound (for fixed 0) is strictly concave for appropriate choice of priors, the maximizer is unique.
16

---

Published as a conference paper at ICLR 2025
Early stop training
103
70
-0.5
60
Best value
Loss
50
-1.0
10000 Epochs
2500
5000
7500
10000 Epochs
1.55
w. early stopping
10a
1000 Epochs
0
2500
5000
7500
10000
2500
5000
7500
Epoch
Runtime [8]
(a) Training loss with early stopping for three differ- ent model initializations.
(b) Comparison of runtime and performance for dif- ferent fixed epochs and using patience on Ackley5D.
Figure Illustration of early stopping in VBLL training (left) and its impact on runtime and performance (right). By introducing the early stopping, significantly reduce runtime while maintaining good performance.
second proof approach constructive: note that the variational lower bound is strictly concave and compute the maximizer by computing the gradient and using first order necessary conditions for optimality. By inspection one can see that this is equivalent to the recursive least squares posterior. O
B.2 MODEL RE-USE AND EARLY STOPPING
Early stopping: The training process for BNNs typically follows standard neural network procedures and relies on held-out validation set to determine when to stop training (Watson et al., 2021; Harrison et al., 2024). However, this approach is impractical in the context of BO with limited data, particularly the early stages of optimization when the training data size very small. To address this, we implement early stopping for VBI.I networks based on the training loss. Specifically, we track the average loss for each training epoch and stop the average loss does not improve for M consecutive epochs, using the model parameters that produced the lowest training loss.
Figure 7a shows the training of the same VBLL network with different initialization but using identical data and optimizer settings. The dashed vertical lines indicate that the training loss reaches a plateau at significantly different epochs across these initializations. Such a large variance in the optimal number of training epochs was also observed in the original VBLL paper for regression tasks (Harrison et al., 2024). By applying an early stopping criterion, we are able to specify a high maximum number of epochs while still reducing computational costs, each initialization is only trained as long as necessary to achieve good empirical performance. In Figure 7b, we can see that by using early stopping, the runtime can be reduced while maintaining performance. Only training for small fixed number of epochs has similar savings in runtime but at the cost of final performance.
is important emphasize that this approach is not exclusive to VBLLs; we believe that applying a similar early stopping criterion could prove beneficial for training neural network- based surrogate models BO more broadly. our experiments, we therefore applied early stopping all applicable BNN baselines.
Feature re-use: The naive approach to training requires training new network for each step of full model training. While we mitigate the cost of training via the last layer recursive update, we further explore using continual feature learning for faster convergence. For this, we initialize the VBLLS at each (full model training) iteration with the feature weights from the last full model training iteration, and the variational posterior from the previous recursive update iteration. This warm start, combined with early stopping, yields substantial speed-ups in training. This feature re-use approach follows conventional continual learning, in which networks are trained for streaming data (Zenke et al., 2017; Hadsell et al., 2020). However, discussed in Appendix C.I, re- re-initializing networks plays an important role that has strong connections to Thompson sampling over the features, and thus we will interleave continual (feature) training with occasional model re-initialization.
17

---

Published as a conference paper at ICLR 2025
VBLL
VBLL (Periodic CT)
VBLL (ET CT)
Best value
-10
TS
2
-0
10
Branin (2D)
Ackley (2D)
Ackley (6D)
Hartmann (6D)
-30
20
40
20
40
50
100
50
100
Iteration
Iteration
Iteration
Iteration
Figure 8: Comparison of periodic retraining (Periodic CL) and event-triggered retraining (ET CL).
B.3 RECURSIVE UPDATE SCHEDULING
We next discuss the different approaches for deciding when to re-initialize the model. There are different heuristics approach would be to fix a re-initialization rate M 2 1 a-priori.
apply to decide when re-initialize the network. A straightforward
B.3.1 MODEL RE- INITIALIZATION THROUGH EVENT TRIGGERING
simple rule for deciding online whether to re-initialize the network is to check how consistent a new observation is with the current variational posterior. For this, simply compute the log-likelhhood of the incoming point and compare it against a threshold A as
log log.Nog, (ye v..) < A O re-initialize model
(13)
Here, A is a user-defined threshold and Nog, is the predictive under the current variation posterior. In all our experiments, we set A = 0 and did not further tune this. Intuitively, this gives practitioners a tuning nob trading off between computational efficiency and predictive accuracy. In Figure 9, show the influence of A on both final performance and runtime.
In Figure 8, we compare the event-triggered strategy to periodic retraining (M = 5) demonstrating an improved performance on all benchmarks. However,
it should be noted that with periodic retraining, the total runtime can be better estimated as fraction of retraining at every time step, whereas this is not possible with the event-triggered strategy.
Bost value
A -10
A 0
=
A = 0.1
0.5
6
0..1
B.3.2 MODEL
RE-INITIALIZATION THROUGH SCHEDULING
2000
4000
6000
Runtime
[+]
Another idea to decide whether to re-initialize the model perform recursive update is through scheduling. Various scheduling approaches can be leveraged here, but scheduling re-initialization with a sigmoid arises naturally: in the beginning, when data still spare, we want a high retrain rate, whereas the later stages, the basis functions of the model expressive enough model the function. We can then use efficient recursive updates. Specifically, we parameterize the probability of re-training with sigmoid taking as input the current time step
Figure 9: Performance and runtime for vari- ous thresholds A on Ackley5D with TS. The threshold defines trade-off between final per- formance and total runtime.
p(t)
Z Bernulli(p(t)) re-initialize model.
(14)
+ (e-(c)-
Here, is the center of the sigmoid and is the stretch parameter. We set the stretch parameter as 2 2ln(9) where w is the transition window ratio.
T.
Remark. Setting the stretch parameter
2 1n(9)
T ensures that the sigmoid transition from p(t) = 0.9 We
to p(t) = 0.1 occurs
time steps. With the sigmoid as p(t) I+e
18

---

Published as a conference paper at ICLR 2025
p(t1) 0.9 and p(t2) = 0.1, and solve for t1 and t2 t In(9) and 82 length of the transition window is l2 t 2ln(0) Setting this equal to the desired window w T, we have 2ln(9) w. which implies s 21n(9) T-w
c+ In(9) The
As defaults, we set the center to 2/2 and w = 0.5. However, similar to the threshold for the event trigger, practitioners can use these tuning knobs prioritize computational efficiency or predictive accuracy.
B.3.3 COMPARISON RE-INITIALIZATION STRATEGIES
We next compare the two main strategies for determining when to re-initialize the model. For this, we define the following indicator function:
{
if model is re-initialized at iteration t,
T(() 0 if model updated recursively iteration f.
(15)
We compare the performance and re-initialization behavior of the approaches on selected subset of experiments in Figure 10. We can observe that for TS, both approaches perform similarly and almost
VBLL
VBLL (ET CL)
VBLL (Sigmoid CL)
Best value
-10
TS
-5
2
-0
Branin (2D)
10
Ackley (2D)
10
Ackley (6D)
Hartmann (6D)
30
ET: -85. 14% RT
ET: -64.77% RT
ET: -74.16% RT
ET: -78.04% RT
SG: -60. 73% RT
SG: -54.97% RT
SG: -67 47.35% RT
SG: -68..2%% RI
Indicstor function to relearn model
1.0
1.0
1.0
1.0
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
20
20
40
100
100
7500
12
log Best value
200
5000
14
- 16
2500
NN draw (200D)
Pest control (25D)
Lunar lander (12D)
ET: -61 61.42% RT
ET: 77777.9%% RT
ET: -27 27% RT
SG: -2.62% RT
SG: -50.12% RT
SG: -30.82% RT
Indicator function
to relearn model
1.0
1.0
1.0
0.5
0.5
0.0
0.0
0.0
200
400
600
50
100
200
Iteration
Iteration
Iteration
Figure 10: Comparison of wriggeerggeeed and scheduled re-initialization on a selected subset of problems. The top always shows the performance, and the bottom row shows [E[[(t) across seeds. In between, we report the improvement in runtime compared to the standard VBLL baseline.
on par with the VBLLs but at a significant reduction total runtime. Through all experiments, we observe the more adaptive behavior of the event trigger compared to the pre-defined sigmoid schedule. Specifically for Branin, we can observe that the event trigger can significantly reduce the runtime already early
19

---

Published as a conference paper at ICLR 2025
For -igh-dimensional problems using logEI, observe a notable difference in performance for NNdraw. When employing the event trigger. the algorithm stagnates resulting in poor performance. In contrast, sigmoid scheduling allows the algorithm to achieve performance levels that approximate those of the VBL.s. Interestingly, we do not observe the same behavior when using TS (cf. Figure 3). Here, the additional exploration introduced by Thompson sampling appears mitigate these issues. Based on these observations, we hypothesize that especially the interplay of TS with an online strategy leveraging event-triggered re-learning offers significant potential for practical applications.
Finally, it is important to acknowledge that the comparisons in this section are based on a single set of parameters, limiting the scope for direct comparisons, particularly due differences in runtime savings. Moreover, while both strategies arise naturally, an intriguing direction for future work would be investigate more advanced approaches, such as those from model selection.
B.4 TRAINING DETAILS
For training the VBLL models, we closely follow Harrison et al. (2024). For all experiments, we use AdamW (Loshchilov & Hutter, 2017) as our optimizer with a learning rate of 10 33 set the weight decay for the backbone (not including the parameters of the VBLL. to 10-1, and use norm-based gradient clipping with value of 1. For the VBI.I, we set the prior scale to and the the Wishart scale to 0.01. A sensitivity analysis of the Wishart scale on the performance and run time is in Appendix E. As mentioned in Sec. 3, we employ early stopping for all VBLL models based on the training loss. We track the average loss of training epoch and if this average loss does not improve for 100 epochs in a row, we stop training and use the model parameters that yielded the lowest training loss.
B.5 COMPUTATIONAL COMPLEXITY
The update formula for the last layer (defined by the rank-1 Cholesky update for the covariance and (8)) has quadratic computational complexity in the feature dimension. Computing wù can similarly be computed with quadratic complexity as the product L L, qu- This can be computed via efficient linear solves due to the triangularity of L, We also note that computing the log determinant - loog det S, is equivalent to computing log det S, = 2 log det L,, which is equal to the sum of the log diagonal elements, and thus has linear complexity.
The KL. penalty on the last layer variational posterior requires computing the trace of the covariance, tr(L, L, )) This can be expressed as ||? and complexity for this is cubic. We can turn instead to randomized algorithms- in particular, Hutchinson estimators (Hutchinson, 1989) to yield an approximate solution with quadratic complexity. This estimator is
solve L,, 1,,,)= tr(L, 1 = S
N E solve((] ,G)²
1
(16)
for appropriately sampled (for example, ~ NN T)). Here, solve(::.) denotes linear solve. Thus this estimator is quadratic in the feature dimension and linear in N. The MSE convergence rate for this estimator is linear, but modifications exist that achieve quadratic convergence (Meyer et al., 2021) by exploiting quasi-Monte Carlo methods.
With the exception of the trace computation, the computational complexity is quadratic in the feature dimension. Therefore, for an exact solution the complexity is cubic (due to the cubic complexity trace computation), and rapidly convergent randomized solution is possible with quadratic complexity.
B.6 FURTHER DETAILS ON THE SETUP AND ACQUISITION FUNCTIONS
We implement the VBLLS within BoTorch (Balandat et al., 2020). We further build on the imple- mentation of Li et al. (2024) for the different baselines which are also based on BoTorch as well as GPyTorch (Gardner et al., 2018). As best-practice in BO, we standardize the data to mean zero and a variance of one at each iteration. We further transform the input space specified by the problem intro the hypercube X' E [0, 1). For completeness, we again list all the baselines below and then discuss the optimization of the acquisition functions.
Acquisition functions: We use 10 restarts and 512 raw samples for optimizing the acquisition functions UCB and 1 logEI for all models. For TS optimize the analytic sample of the VBLLS
20

---

Published as a conference paper at ICLR 2025
with 10 random restart using L-BFGS-B (Zhu et al., 1997) the optimization method. For TS from the non-parametric models, we use the same heuristic as Eriksson et al. (2019) and generate min (5000, (200, 200 }}} pseudo-random input points from Sobol sequence and then sample from the high-dimensional multi-variate normal. The next location is then the argmax of the sample path. For the multi-objective problems, use 1 logEHVI for all models with consistent reference points, For IS with VBLLs, we follow the procedure described in Sec. 4. For NSGA-II, we use population size of 100 and terminate after 200 generations.
C DISCUSSION AND BASELINE ABLATIONS
C.I ON NEURAL NETWORK THOMPSON SAMPLING
We further investigated the use of using neural network based Thompson samples. With VBLLs, effectively maintain distribution over plausible deterministic NNs that are in accordance with the current data set and noise level. Instead of maintaining such distribution and then sampling from the variational posterior of the weights w, a straight-forward idea would be directly train a deterministic neural network using the same optimizer (including the same early stopping etc.), and L2 loss, as well as L2 regularization for the backbone directly generating MAP Thompson sample. Note that with such an approach, one would no longer be able to leverage acquisition functions that rely on good uncertainty quantification such as logEI or more sophisticated such as ieformation-theoree
acquisition functions based on entropy search such as MES Wang & Jegelka (2017), which may require large ensembles (up to 100 NNs), making it computationally expensive. Still, we tested this baseline and the results are summarized in Fig. 11.
GP
LBNN
DKL
LLLA
VBLL
MAP
Best value
TS
-20
-5
-40
Branin (2D)
-10
Ackley (2D)
-0
Ackley (BD)
Hartmann (6D)
20
40
20
40
50
100
50
100
Tteration
Iteration
Tteration
Iteration
Figure 11: Comparison of neural network Thompson sampling synthetic benchmarks. The deterministic MAP Thompson sample shows surprisingly good performance however also yields large variance on simple benchmarks which is undesirable.
Here MAP refers to the MAP Thompson sample baseline. We can observe that this simple baseline performs surprisingly well but cannot match the performance of the other baselines. We also observe that in relatively simple problems, such as Branin n and Ack lel2De, the variance is significantly larger compared the other baselines, which is undesirable in BO applications. We also compared this baseline on the high-dimensional and real-world benchmarks and the results are shown in Fig. 12.
The MAP approach demonstrates superior performance on the NNdraw task. This is likely because the Wishart scale the VBLL baselines not well-tuned for the NNdraw problem (cf. Sec. E). The data is normalized to zero mean and a standard deviation of one; however, with Wishart scale greater than zero, the assumed noise in this normalized space affects the accuracy the correlations between data points-especially for the large value range in NNdraw. In contrast, the MAP baseline does not, by design, account for noise, which may contribute to its better performance in this context. Additionally, the MAP baseline exhibits slightly faster convergence Pestcontza For Lunarl Ietaaraanderr the MAP approach shows considerable variance and fails to match the final performance of the VBLLs.
21

---

Published as a conference paper at ICLR 2025
GP
I-BNN
DKL
LLLA
VBLL
MAP
8000
12
300
Best valuee
6000
200
14
1S
4000
100
-16
2000
(2009)
Fest control (25D)
Lunar lander (12D)
200
400
600
50
100
200
Iteration
Tteration
Iteration
Figure 12: Comparison of neural network Thompson sampling high-dimensional and real-world benchmarks. On these benchmarks, the MAP Thompson sample baseline shows mixed performance. performs better on NNdraw but exibits large variance for Lunarlander
GP
GP (with D-scaled hyperprior)
VBLL
-2
7500
Best valuue
TS
-5
5000
-10
Ackley (5D)
2500
NW draw (206D)
-16
Pest control (25D)
0
-12
log Best value
7500
-5
114
5000
-16
-10
Ackley (5D)
2500
NN draw (200D)
Pest control (25D)
50
100
200
400
600
50
100
Iteration
Iteration
Iteration
Figure 13: Comparison of standard GP with and without D-scaled hyperprior to VBLLS.
C.2 COMPARISON TO D-SCALED GAUSSIAN PROCESS PRIORS
We further compare the VBLLs to Gaussian process with so-called D-scaled hyperpriors on all the lengthscales which has demonstrated promising results in high high-dimensional settings (Hvarfner et al., 2024). With this, we aim to test the hypothesis that the VBLL's improvement in performance over the GPs tasks such Pestcontrol does not come from bad GP hyperparameters but rather from its ability to capture complex non- -on-Euclidean input correlations. For the hyperprior, we closely follow Hvarfner et al. (2024) and define log-normal distribution the hyperprior for all lengthscales as
log D
p(es) LN
N4
2 000
(17)
and set the constants to μo = /2 and ou = /3. As in Hvarfner et al. (2024), we no longer use box constraints on the lengthscales for the D-scaled GP. The lengthscales are initialized with the mean of the hyperprior at each iteration and then optimized by minimizing the marginal log-likelihood. In the following, we focus only on Ackley (5D), Pestcontrol, and NNDraw. The results are in Figure 13
We can observe that for Ackley, the additional hyperprior does not have an influence on the per- formance. For the NNDraw example, the inductive bias introduced by the hyperprior benefits performance when using logEI, but it is negligible for Thompson sampling. In high-dimensionall settings, the critical factor for Thompson sampling is not the kernel choice but whether the surrogate model is parametric performance but also negatively impacts it when using logEI. Even with the additional hyperprior, the GP model is not able to capture the input correlations, which prevents it from yielding good
non-pann parameeriic. On Pestcontrol, the hyperprior not only fails to improve
22

---

Published as a conference paper at ICLR 2025
performance. Here change in kernel is necessary to achieve competitive performance with GP surrogate model. The VBLLs and the other parametric BNN models can capture this complex input correlation and converge fairly quickly to the optimum, given that the problem has 25 dimensions. Lastly, it is important to note that also for the VBLLS and the other baselines, different hyperparameters such network architecture also improve performance (see Appendix E.3).
C.3 COMPARISON OF SENSITIVITY TO NOISE FOR LLLA AND VBLLS
We conduct additional experiments to compare the robustness of the LLLA and VBLL surrogates to objective noise. Specifically, consider Ackley2D and Ackley5D and set the output noise standard deviation to a E [0.00 0.001, 0.01, 0.1, 1] to test a broad spectrum. For each σ, we run 10 seeds. The results are summarized in Figure 14. The left plot displays the mean and the 10th and 90th percent quantiles across all σ.. We can observe that the performance is more consistent for VBLLs on the lower dimension problem. The center and right plots show mean and quantiles for each for LLLA and VBLLs, respectively. As reference, we overlay the performance of the noise-free case blue. Also here we can observe that while the mean of the best observed values remains similar for Ackley2D, VBLLs display higher consistency compared to LLLA considering the quantiles. A similar but less amplified trend can be observed on Ackley5D.
0= 0.1
a 0.01
0.001
0.0
TS (Ackley2D) Best value
-5
-5
LLLA
-10
VBLL
-10
LLLA
-10
VBLL
20
40
20
40
20
40
0
0
TS (Ackley Best value
-5
-5
VBLI.
LLLA
-10
-10
LLLA
VBLL
50
100
50
100
50
100
Iteration
Iteration
Iteration
Figure 14: Sensitivity to noise on Ackley2D and Ackley5D for LLLA and VBLLs
D EXPERIMENTAL DETAILS
D. BASELINES
As the de-facto standard in BO, we compare against GPs. As kernel, we choose a Matérn kernel with v = 2.5 and use individual lengthscales f; for all input dimensions that are optimized within box constraints following recommended best practices (Eriksson et al., 2019; Balandat et al., 2020). We use box constraints f E [0.005, (Eriksson et al., 2019)-we refer to Appendix C.2 for further discussion on this choice as well as a comparison to D-scaled GPs (Hvarfner et al., 2024). We expect the performance of GPs to be particularly good stationary benchmarks.
GPs:
We compare against infinite-width Bayesian neural networks (I-BNNs) (Lee et al., 2018), which have shown promising results in recent work (Li et al., 2024). As in Li et al. (2024), we set the depth to and initialize the weight variance to 10 and the bias variance to 1.6. Note that the I-BNN expressed as a kernel and therefore the model is still non non-parametric and does not learn features.
I-BNNs:
'Note that we track the best observed value. Therefore, values greater zero possible (see 1 for the VBLLS is D)) though max [[)()))) for Ackley.
23

---

Published as a conference paper at ICLR 2025
DKL: Deep kernel learning (DKL) (Wilson al., 2016) combines feature extracting with neural networks with GPs. It uses a neural network ge as imput wrapping to the kernel as İDKL = ({)) (o))))"(()))) (o to allow for non-stationary modeling and exact inference. For the neural network, we use the same architecture as in (Li et al., 2024), i.e., 3 layers with 128 neurons. We further use ELU activations for all layers.
Last layer Laplace approximations (LLLA) are a computationally efficient way obtain uncertainty estimates after training a neural network (MacKay, 1992; Daxherger et al., 2021). With this, they are a well suited BNN surrogate model for BO (Kristiadi et al.; Li et al., 2024). As NN, we also use layers with 128 neurons and ELU activations. Note that for TS, can also optimize the parametric neural network function numerically with LLLA.
LLLA:
cheap and efficient way to obtain predictions with uncertainty. Each member of the ensemble is initialized with a different random seed, and all are trained on the same data, minimizing an MSE loss using weight decay. We use 5 ensemble members. We parameterize the predictive as a normal distribution with mean and variance from the ensemble. We note that an ensemble of VBLL networks is also possible, although we do not include this approach.
DE: Deep ensembles
VBLL:
For the VBLLs, we use the same architecture as for DKL and LLLA. In the body of the paper, we include two VBLL models: training the features from scratch every iteration (performant but expensive) and re- re-training the features based on an event trigger and otherwise doing recursive last layer updates (VBLL (ET CL)).
D.2 EXPERIMENTS
Branin: A standard two dimensional optimization benchmark with three global optima.
Ackley: A standard optimization benchmark with various local optima (depending on the dimen- sionality) and one global optimum. In out experiments. we compare the surrogates on 2D and 5D version and set the feasible set to the hypercube X' [-,, 10j4 as in Eriksson et al. (2019).
Hartmann: A standard six dimensional benchmark with six local optima and one global optimum.
NN draw: In this optimization problem, our goal is to find the global optimum of a function defined by sample from neural network within the hypereube X' [0, 1]d. This benchmark was also employed in Li et al. (2024). We use a fully connected neural network with two hidden layers, each containing 50 nodes, and ReLU activation functions. The input size corresponds to the dimensionality of the optimization problem (in our case, 200), and the output size is one. To generate a function, we sample all weights from the standard normal distribution (0, 1). For a fair comparison, we use the same fixed seed across all baselines ensuring that the same objective function is used.
This optimization problem was also in Li et al. (2024) and aims to minimizing the spread of pests while minimizing the prevention costs of treatment and was introduced in Oh et al. (2019). In this experiment, we define the setting as categorical optimization problem with 25 categorical variables corresponding to stages of intervention, with different values at each stage. As mentioned in Oh et al. (2019), dynamics behind this problem are highly complex resulting in involved correlations between the inputs.
Pest control:
Lunar lander is an environment from OpenAl gym. The objective is to maximize the average final reward over 50 randomly generated environments. For this, 12 continuous parameters of a controller have to be tuned as in Eriksson et al. (2019).
Lunar lander:
Branin Currin: A standard two dimensional benchmark with two objectives.
DTLZ1: A standard benchmark of which we choose five dimensional version with two objectives.
DTL.22: A standard benchmark of which we choose a five dimensional version with two objectives.
Oil Sorbent:
Building the implemention by Li et al. (2024), optimize the properties of material to maximize its performance as a sorbent for marine oil spills as introduced by Wang et al. (2020). The problem consists of five ordinal parameters and two continuous parameters which control
24

---

Published as a conference paper at ICLR 2025
the manufacturing process of electrospun materials, and the three objectives water contact angle, oil absorption capacity, and mechanical strength.
E HYPERPARAMETER SENSITIVITY
The parametric VBLI surrogate has hyperparameters that have to be specified a-priori. In the following, we present results on the hyperparameter sensitivity of the VBLL surrogate model and demonstrate that tuning hyperparameters can improve empirical performance but also that the VBLL surrogate model is rather robust for a wide range of specifications. In Sec. E.I, we will first consider the sensitivity with respect to the Wishart scale and the reinitialization rate for continual learning. Following this, Sec. E.3 then studies the sensitivity regarding the width of the neural network backbone. Lastly, Sec. E,4 considers the robustness to different noise levels.
E.1 WISHART SCALE AND CONTINUAL LEARNING SENSITIVITY
We sweep number of hyperparameters in the VBLLs in order to experiment with the hyperparameter sensitivity of the VBLL models. In particular we sweep the Wishart scale and the re-initialization rate of the model. The re-initialization rate determines how often the VBLL model is re-initialized rather than using CL on the backbone and the variational posterior.
The results of the hyperparameter sweep of the Wishart scale and reinitialization rate on Ackley (SD) can be seen in Fig. 15 and on Pestcontrol in Fig. 16. Please note that we have computed running averages on these figures to make qualitative assessment easier. We find that both of these hyperparameters have impact on BO performance, but that the VBLL models are not exceedingly brittle to the values of these hyperparameters. These results also indicate that continual learning for VBLLs is an area of interest, not only to reduce fitting time, but also because there are indications that the VBLL surrogate benefits from not always being re-initialized (see e-g., Fig. 15 (b) for reinitialization rates 3 and 5). Based on these results we also hypothesise that tuning the Wishart scale appropriately for the problem at hand may lead to increased model performance.
E.2 INFLUENCE OF THE WISHART SCALE ON THE PERFORMANCE OF CONTINUAL LEARNING
We can leverage the insights from the previous section also to improve the performance of the continual learning baseline for logEI. We saw that decreasing the Wishart scale can improve performance. The default value the Wishart scale used in all experiments in Section 5.2 and Section 5.3 was V = 0.01. This induces a hyperprior on the noise with mass centered at o² 0.01. However, expected improvement is known to get stuck for noisy objectives. We can observe that combined with recursive updates this effect is increased compared the standard VBLL baseline which still converges well due to the additional randomness of a new initialization (cf. Figure 3). In Figure 17, we set the Wishart scale to the standard VBLLS as a reference. We can observe that with this tuned Wishart scale, we can significantly improve the performance on benchmarks such as Pestcontrol. This further shows that similar to tuning hyperprior for GP, tuning hyperparameters can have benefit.
10 for the continual learning baseline. We include
E.3 MODEL WIDTH ABLATION
To evaluate the impact of model width on the performance of both the MAP and VBLL Thompson sampling methods, we conducted a series of experiments on the Ackley5D and rest benchmarks varying the model width.
As illustrated in Fig. (a) and 19 (a), increasing the model capacity (width) of the MAP baseline results in a significant increase in variance, especially pronounced in the Ackley5D benchmark. This high variance suggests that the MAP method is highly sensitive to changes in model width, making it challenging tune effectively for consistent performance across different tasks.
In comparison, the VBLL method exhibits more robustness to model capacity, as shown in Fig. 18 (b) and 19 (b). Despite increasing the model width, VBLL does not suffer from the high variance observed in the MAP baseline. Additionally, the next x distance metric plots in both Fig. 18 and 19 indicates that later BO stages the VBLL models return exploring uncertain regions of the input
25

---

Published as a conference paper at ICLR 2025
1c-05
000000
0.001
0
Next distance
1.0
Accum. Fit Time
Best valuu
40
-5
0.5
20
-10
0.0
0
Next distance
1.0
Accum. Fit Time
50
Best value
-5
0.5
110
0.0
Iteration
Iteration
Itratiion
(a) Sensitivity to the Wishart scale (noise free (top), noisy (bottom))
0
Next distance
Accum. Fit Time
40
Best value
-5
0.5
20
10
0.0
0
0
Next distance
Accum. Fit Time
Best value
40
5
0.5
20
-10
0.0
0
50
50
50
Iteration
Iteration
Iteration
(b) Sensitivity to network reinitialization rate (noise free (top), noisy (bottom))
Figure 15: Hyperparameter sensitivity on the Ackley5D benchmark.
space whereas the MAP models get stuck in local optima accounting for the continous improvement best values for the VBL.I models. This robustness is advantageous practical application where extensive model tuning is unfeasible and hints that VBL. may also perform better when scaling to larger model sizes.
E.4 MODEL PERFORMANCE IN THE PRESENCE OF NOISE
Lastly, also benchmark the different surrogates different noise levels. We again only consider Ack Ee (Fig. 20) and Pest.ont.rs (Fig. 21). For these experiments, we use the same Wishart scale of 0.01 for the VBLL baseline. We can observe that all models, besides the MAP baseline in Ackley5D, are rather robust the change in noise level.
26

---

Published as a conference paper at ICLR 2025
1c:05
12
Next a distance
Accum. Fit Time
Best valtue
16
12
Next æ distance
2
Accum. Fit Time
Best valute
14
-16
Iteration
Iteration
Iteration
(a) Sensitivity to the Wishart scale (noise free (top), noisy (bottom))
-22
Next or distance
2
Accum. Fit Time
Best value
50
-14
25
16
0
-12
Next ir distance
2
Accum. Fit Time
50
Best valute
-14
25
-16
0
50
50
50
Iteration
Iteration
Iteration
(b) Sensitivity to network reinitialization rate (noise free (top), noisy (bottom))
Figure 16: Hyperparameter sensitivity on Pestcontrol benchmark.
27

---

Published as a conference paper at ICLR 2025
VBLL
VBLL (ET CL), 0.01
VBLL (ET CL),
Best value
logel
5
Branin (2D)
-10
Ackley (2D)
10
Ackley (5D)
Hartmann (6D)
-30
20
20
40
100
100
12
Best value
6000
200
log
-14
4000
100
-16
2000
NN draw (200D)
Peat control (25D)
Lunar lander (12D)
200
400
600
100
200
Iteration
Iteration
Iteration
Figure 17: Classic benchmarks (top) and high-dimensional and hon-stationary benchmarks (bottom). Performance of all surrogates for logEI (top) and TS (bottom).
64
128
512
0
Next x distance
1.0
Best. value
0.5
-10
Ackley (5D)
0.0
Ackley (5D)
25
50
75
25
50
75
Iteration
Iteration
(a) MAP Thompson sampling with varying model width
64
128
512
Next distance
1.0
Best value
5
0.5
-10
Ackley (5D)
0.0
Ackley ey (5D)
0
25
50
75
0
25
50
75
Iteration
Iteration
(b) VBLL Thompson sampling with varying model width
Figure 18: Comparison of neural network Thompson sampling methods on the Ackley5D bench- mark with varying model width. The models were trained with width of 64, 128 and 512 neurons.
28

---

Published as a conference paper at ICLR 2025
64
128
512
12.5
Next: x distance
Best value
-15.0
Pest control (25D)
Past control (25D)
17.5
25
75
50
75
Iteration
Iteration
(a) MAP Thompson sampling with varying model width
64
128
512
Next a distance
2
12.5
Best value
15.00
Pest control (25D)
Pest control (25D)
17.5
25
50
75
25
50
75
Iteration
Iteration
(b) VBLL Thompson sampling with varying model width
Figure 19: Comparison of neural network Thompson sampling methods on the restcontrol benchmark with varying model width. The models were trained with width of 64, 128 and 512 neurons.
GP
I-BN
VBLL
MAP
Ackley (5D) Best valtee
5
-5
-5
-10
-10
-10
UCB
logEI
TS
50
100
50
100
50
100
Iteration
Iteration
Iteration
(a) Noisy objective with noise standard deviation of 0.01.
0
Ackley (5D) Best tte
-5
-10
-0
UCB
logEI
TS
50
50
100
50 Iteration
100
Iteration
Iteration
(b) Noisy objective with noise standard deviation 0.1.
Figure 20: Performance comparison of baseline methods Ackley5D benchmark with noise.
29

---

Published as a conference paper at ICLR 2025
GP
I-NN
VBLL
MAP
Pest control (25D)
-22
12
12
Best value
-14
14
14
-16
-16
UCB
logEI
TS
50
100
50
100
50
100
Iteration
Iteration
Iteration
Noisy objective with noise standard deviation of 0.01.
Pest control (25D
Best value
-14
14
14
-16
-16
16
UCB
logEI
TS
50
100
50
100
50
100
Iteration
Iteration
Iteration
(b) Noisy objective with noise standard deviation of 0.1.
Figure 21: Performance comparison of baseline methods on Pestcontrol benchmark with noise.
30


Reviews:
[{"review_id": "9673IFZhVr", "reviewer": "Reviewer_wtu3", "summary": "The authors propose VBLL networks as a surrogate for Bayesian optimization and identify a relationship between optimizing the variational posterior and recursive Bayesian linear regression. They demonstrate competitive BO results with VBLL on diverse single and multi-objective benchmarks.", "questions": "- Could you further elaborate on the differences between LLLA and VBLL? LLLA outperforms VBLL on many of the non-stationary benchmarks, and from line 530, it appears that VBLL-based approaches may be less flexible than last-layer Laplace approximations.", "limitations": "", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- The method is well-explained and is theoretically justified, and there are additional modifications which can be made to increase the efficiency such as feature re-use and sparse full model retraining. This flexibility enables practitioners to balance the tradeoff between model performance and computational cost.\n- The authors use a diverse setting of test objectives, specifically demonstrating performance on instances with high-dimensionality and non-stationarity. \n- VBLL appears to be robust to hyperparameter choices and can be used as a drop-in surrogate model, unlike typical GPs which require careful kernel selection.", "weaknesses": "- Although it appears that one of the primary motivations behind this work is the increased efficiency compared to other BNN surrogates, there is no measure of runtime or computational cost within the paper. It would be helpful to understand how these methods perform as a function of computational budget. This could also help clarify the difference in performances between VBLL and VBLL CL. \n- There is also currently no demonstration of why this approximation would be preferred over the using the exact marginal likelihood with approaches like DNGO [1]. Without these baselines, there is minimal evidence that the proposed method has practical merit over existing work. Furthermore, it would also be useful to compare last-layer methods like VBLL to more expensive BNN surrogates like deep ensembles so we can assess the tradeoff between computation and performance. \n\n[1] Snoek et al, Scalable Bayesian Optimization Using Deep Neural Networks, 2015", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ieWNxVjDDI", "reviewer": "Reviewer_1qWo", "summary": "This paper presents a combination of two ideas: (!) Bayesian last-layer training of neural networks with (2) using a parametric Bayesian linear model for black-box function optimization. Using natural parameterization of the last layer Gaussian and assuming independent Gaussian noise, it is possible to use a continual update which is only $O(N^2)$ in the last-layer number of neurons $N$. As with every parametric Bayesian function model, the acquisition function can be directly and analytically optimized using Thompson sampling. The empirical results on a wide variety of Bayesian Optimization tasks are promising.", "questions": "* In line 110 und 111 it helps to point out that $\\mathbf{w} = S^{-1} \\mathbf{q}$ is the vector of precision-means (for those unfamiliar with natural parameters of the Normal distribution)\n* The proof of Theorem 1 in the appendix relies on the simple observation that the approximating family contains the true distribution. I would have preferred to see this in the main body of the text; it's less \"mechanical\" than I expected and key to the reasoning of the paper (line 178-189 can be significantly shortened b/c it uses the well known Cholesky decomposition for approximating the inverse and log-determinant of the precision)\n* Line 217: Where is the parameter $V$ (Wishart prior scale) necessary?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The paper combines two well studied ideas in an elegant way and the presentation is (relatively) easy to follow (though I would have wished a bit more emphasis on the natural parameterization of the Normal distributions as this is key to the computational efficiency). The empirical studies are extensive and well discussed.", "weaknesses": "One aspect that is disregarded by the paper is how to chose the network architecture for all but the last-layer; I have no idea how sensitive the quality of the proposed approach is to this. In essence, the complexity of choosing a kernel function for GPs has been shifted to the network architecture of the underlying neural network. This is not discussed in sufficient detail. Also, only at the end the difference to Laplace approximation of the last layer is discussed; I would have expected this in the Related Work section.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bjJ39n87cw", "reviewer": "Reviewer_StZZ", "summary": "Although Gaussian processes (GPs) are widely used for Bayesian Optimisation (BO), they are not always well suited to modelling functions with complex correlations across inputs, and are often limited by the choice of kernel function. On the other hand, Bayesian neural networks (BNNs) can better handle complex non-Euclidean data, but are computationally expensive and challenging to condition on new data. In this work, the authors extend recent work on Variational Bayesian Last Layer (VBLL) models specifically for BO. They show how VBLLs can be adapted as efficient surrogate models for BO tasks through modified training procedures, which enable continual, online training with recursive conditioning, improving scalability. Additionally, the authors demonstrate how VBLL’s parametric structure enables effective Thompson sampling for both single- and multi-objective acquisition functions, offering more stability and numerical efficiency compared to GPs. Experiments compare VBLL’s performance against other techniques such as baseline GPs, and BNNs, and show that VBLL performs especially well on complex tasks, such as the Oil Sorbent benchmark, where other approaches struggle due to numerical instability.", "questions": "Please refer to comments in \"Weaknesses\".", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-written and a pleasure to read. The problem statement is clear from the outset, and the connections to related work are extensive. I also appreciated how the paper’s focus on practical aspects such as improving training efficiency via continual learning. Having the method implemented in BoTorch is also appealing to practitioners wanting to experiment using this method in real-world settings. \n\nThe experiments demonstrate that VBLL performs well in the targeted settings having complex input correlations and non-Euclidean structures. Showing that VBLL outperforms competing techniques on real-world datasets such as the Oil Sorbent and Pest Control datasets adds further credence to how VBLL is suited to multi-objective settings prone to numerical instability.\n\nWhile the contributions may initially appear incremental, adapting VBLL to BO introduces challenges that require non-trivial solutions. The need for efficient, online training in BO necessitated the development of recursive conditioning and continual learning updates, which are distinct from standard regression tasks. Addressing the requirements of multi-objective and high-dimensional settings also required effective workarounds to address numerical stability issues.", "weaknesses": "1. Deep kernel learning was the first method to come to mind when reading the motivation for this work. While I appreciated its inclusion in the experimental section, I would have liked more discussion in the earlier sections on why DKL might be less ideal than VBLL. To my understanding, DKL’s computational complexity, especially in high-dimensional settings, might be a key differentiator, but additional detail on this would help clarify VBLL’s practical advantages right from the outset.\n2. While the experiments are quite extensive, I would appreciate more insight on the cases where the method is expected to underperform compared to other approaches. Although dedicated experiments are provided in the supplementary material, high-level insights on possible sensitivity to hyper-parameters could also be included in the main text.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "lGopSQgMyt", "reviewer": "Reviewer_asqH", "summary": "Gaussian Process (GP) models are widely regarded as state-of-the-art surrogates for Bayesian Optimization, thanks to their strong predictive performance, efficient updates for small datasets, and straightforward uncertainty quantification. However, they can be challenging to apply to search problems with non-stationary or complex correlation structures. In contrast, Bayesian Neural Networks (BNNs) are more flexible and can handle such problems with minimal adaptation, but they have traditionally been associated with high computational costs and unreliable uncertainty estimates.\n\nThis paper introduces a new method that leverages variational last-layer BNNs, combining advantages of both GPs and BNNs. The proposed approach demonstrates superior performance over GPs and several other BNN architectures in tasks with complex input correlations while achieving comparable results to GPs on standard benchmark problems.", "questions": "1. I think the runtime advantage of the suggested algorithm must be more clearly presented, since it is the motivation of much of the methodology. Specifically, be accompanied by results showcasing a benefit (in e.g. regret/performance per wallclock time), especially compared to LLLA which in terms of regret per iteration performs similarly. \n\n2. There are a couple of sections in the methodology which I think unjustifiably and unnecessarily make claims without backing (see Weaknesses). I recommend the authors to look over the claims and make sure they have backing; either by adding relevant proofs, experiments or references for claims important to the paper, or lessening/removing claims which may be unnecessary. It is okay that not every design decision in a larger algorithm (or system or model) is fully backed, but then those design decisions should arguably not be presented as central parts of the methodology accompanied by unbacked claims. \n\nOverall, the paper addresses a meaningful gap in the literature. If the concerns outlined above are addressed, I would be inclined to raise my evaluation score.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 1, "presentation": 3, "contribution": 3, "strengths": "The paper makes a valuable contribution to the field by highlighting the capabilities of a class of BNN models for Bayesian Optimization (BO) and introducing a practical technique for efficient updates in online settings, including BO scenarios. The writing is clear and well-structured, and the findings are substantiated by experiments conducted in both single-objective and multi-objective settings.", "weaknesses": "The presented method appears to perform comparably to Last Layer Laplace Approximations (LLLA) without clearly demonstrating new advantages. The paper emphasises computational concerns — and presumably could provide a favourable runtime-performance trade-off—but that is not empirically validated and the method is seemingly not compared against the baselines in this respect.\n\nThe claim that Laplace approximations are more sensitive to noise lacks sufficient support, as there is no accompanying experiment or reference to substantiate it (see Section 6, \"Performance and Baselines\"). Providing evidence here would strengthen the argument.\n\nThe discussion around early-stopping based on training loss (see Section 3.2, \"Early Stopping\") makes significant claims, such as training \"only as long as necessary to achieve optimal performance\" and suggesting that applying a similar criterion could benefit training neural network-based surrogate models in BO more broadly. While it is reasonable to argue that stopping training before full convergence improves runtime efficiency and can serve as a regularisation heuristic, the lack of experimental results of the trade-off in the setting when presented as a methodological contribution is a major omission. The effect of early stopping on the quality of the fitted model should be demonstrated through empirical evaluation, such as predictive error on relevant functions or/and assessing its impact on BO performance.\n\nThe choice of length scales [0.005, 4] for the GP model (see Section 5.1, \"Surrogate Models\") appears to be unsuitable for the high-dimensional benchmarks considered. As demonstrated in (1) \"Vanilla Bayesian Optimization Performs Great in High Dimensions\" (ICML 2024), length scales around \\sqrt{D} are generally more effective for Bayesian Optimization in high-dimensional settings. Using more appropriate lengthscales (specifically adopting a suitable lengthscale prior with mass concentrated near \\sqrt{D}) could potentially dramatically enhance the model's performance, making it a more informative comparison. \n(1) https://arxiv.org/pdf/2402.02229", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Spotlight)