Deep Learning for Bayesian Optimization of Scientific Problems with High-Dimensional Structure
Authors: Samuel Kim
Peter Y Lu
Charlotte Loh
Jamie Smith
Jasper Snoek
Marin Soljacic
Venue: Accepted by TMLR ()

Abstract: Bayesian optimization (BO) is a popular paradigm for global optimization of expensive black-box functions, but there are many domains where the function is not completely a black-box. The data may have some known structure (e.g.\ symmetries) and/or the data generation process may be a composite process that yields useful intermediate or auxiliary information in addition to the value of the optimization objective. However, surrogate models traditionally employed in BO, such as Gaussian Processes (GPs), scale poorly with dataset size and do not easily accommodate known structure. Instead, we use Bayesian neural networks, a class of scalable and flexible surrogate models with inductive biases, to extend BO to complex, structured problems with high dimensionality. We demonstrate BO on a number of realistic problems in physics and chemistry, including topology optimization of photonic crystal materials using convolutional neural networks, and chemical property optimization of molecules using graph neural networks. On these complex tasks, we show that neural networks often outperform GPs as surrogate models for BO in terms of both sampling efficiency and computational cost.

Published in Transactions on Machine Learning Research (09/2022)
Deep Learning for Bayesian Optimization of Scientific Prob- lems with High-Dimensione Structure
Samuel Kim
samkim@@mit.c edu
Deportment of Electrical Engineering and Computer Science Massachusetts Institute of Technology
Peter Y. Lu
Department Physics
Massachusetts Institute of Technology
Charlotte Loh
Department of Electrical Engineering and Computer Science Massachusetts Institute of Technology
Jamie Smith
Google Research
Jasper Snoek Google Research
Marin Soljačié
soljacic@mit..cm edu
Department of Physics
Massachusells Institute of Technology
Reviewed on OpenReview: https: // /ttp://pppeee optpr/veew.eeee..o net/ forun? id= tPMQ6Je2r8
Abstract
Bayesian optimization (BO) a popular paradigm for global optimization of expensive black-box functions, but there are many domains where the function is not completely black-box. The data may have some known structure (e.g. symmetries) and/or the data generation process may be composite process that yields useful intermediate auxiliary information in addition to the value of the optimization objective. However, surrogate models traditionally employed in BO, such as Gaussian Processes (GPs), scale poorly with dataset size and do not easily accommodate known structure. Instead, we use Bayesian neural networks, class of scalable and flexible surrogate models with inductive biases, to extend BO to complex, structured problems with high dimensionalily. We demonstrate BO number of realistic problems in physics and chemistry, including topology optimization of photonic crystal materials using convolutional neural networks, and chemical property optimization of molecules using graph neural networks. On these complex tasks, we show that neural networks often outperform GPs as surrogate models for BO in terms of both sampling efficiency and computational cost.
1 Introduction
Bayesian optimization (BO) is a methodology well-suited for global (as opposed to local) optimization of expensive, black-box (e.g. derivative-free) functions and has been successfully applied to wide range of problems in science and engineering (Ueno al., 2016 Griffiths & Hernándzz-Loath. 2020 Korovina al. 2020) as well as hyperparaneter tuning of machine learning models (Snock al. 2012 Swersky et al., 2014
Klein et al. 2017 Turner et. al.] 2020 Ru et. al.] 2021 BO works by iteratively deciding the next data

---

Published in Transactions on Machine Learning Research (09/2022)
point to label in order to maximize sampling efficiency and minimize the number of data points required to optimize a function, which is critical in many contexts where experiments or simulations can be costly or time-consuming.
However, in many domains, the system not a complete black box. For example, certain types of high- dimensional input spaces such as images or molecules have some known structure, symmetrics and invari- ances. In addition, the function may be decomposed into other functions; rather than directly outputting the value of the objective, the data collection process may provide intermediate auxiliary information from which the objective function can be cheaply computed. For example, scientific experiment or simulation may produce a high-dimensional observation or multiple mcasurements simultancously, such as the optical scattering spectrum of a nanoparticle over a range of wavelengths, or multiple quantum chemistry properties of molecule from a single density functional theory (DFT) calculation. All of these physically-infore insights into the system are potentially useful and important factors for designing surrogate models through inductive blases, but they are often not fully exploited in existing methods and applications.
BO relies specifying a surogate model which captures distribution over potential functions to incorpo- rate uncertainty in its predictions. These surrogate models are typically Gaussian Processes (GPs), as the posterior distribution of GPs can be expressed analytically. However, (1) inference in GPs scales cubically in time with the number of observations ions and output dimensionality, limiting their use to smaller datasets or to problems with low output dimensionality without the use of kernel approximations, and (2) GPs oper- ate most naturally over continuous low-dimensional input spaces, so kernels for high-dimensional data with complex structure must be carefully formulated by hand for each new domain. Thus, encoding inductive biases can be challenging.
Neural networks (NNs) and Bayesian neural networks (BNNs) have been proposed as alternative to GPs due to their scalability and flexibility (Snoek et al. 2015 Springenberg et al.) 2016). Alternatively, neural networks have also been used to create continuous latent spaces so that BO with vanilla GPs can be more easily applied (Kusner et al. 2017 Tripp et al., 2020). The ability to incorporate a variety of constraints, symmetries, and inductive biases into BNN architectures offers the potential for BO to be applied to more complex tasks with structured data.
This work demonstrates the use of deep learning to enable BO for complex, real-world scientific datasets, without the need for pre-trained models. In particular:
We take advantage of auxiliary or intermediate information to improve BO for tasks with high- dimensional observations.
We demonstrate BO on complex input spaces including images and molecules using convolutional and graph neural networks, respectively.
We apply BO to several realistic scientific datasets, including the optical scattering of nanopar- licle, lopology optimization of a photonic crystal material, and chemical properly optimization of molecules from the QM9 dataset.
We show that neural networks are often able to significantly outperform GPs as surrogate models on these problems, and we believe that these strong results will also generalize to other contexts and enable BO to be applied to wider range of problems. We note that while our methods are based on existing methods, we use novel combination of these methods to tailor existing BO frameworks to real-world, complex applications.
1.1 Related Work
Various methods have been formulated to scale GPs to larger problems. For example, Bruinsma et al. (2020) proposes a framework for multi-output GPs that scale linearly with m, where m is the dimensionality of a low-dimensional sub-space of the data. Maddox et al.) (2021a) uses multi-task GPs to perform BO over problems with large output dimensionalities. Additionally, GPs have been demonstrated on extremely large datasets through the use of GPUs and intelligent. preconditioners Gardner et al. (2018): Wang et al. (019
2

---

Published in Transactions on Machine Learning Research (09/2022)
through the use of various approximations Rahimi & Recht (007 Wang et al. (2018) Liu et al. (020
Maddox et al. (2021b).
Another approach to scaling BO to larger problems is by combining it with other methods such that the surrogate model does not need to train on the entire dataset. For example, TuRBO uses a collection of independent probabilistic models in different trust regions, iteratively deciding in which trust. region to perform BO and thus reducing the problem to a set of local optimizations (Briksson et al. 2019). Methods such LA-MCTS build upon TuRBO and dynamically learn the partition function separating different regions Wang al. 2020
GPs have been extended to complex problem settings to enable BO on a wider variety of problems. Astudillo & Frazier (2019) decompose synthetic problems as a composition of other functions, and lake advanlage of the additional structure to improve BO. However, the multi-output GP they use scales poorly with output dimensionality, and so this approach
limited to simpler problems. This work has also been extended (2021a). GP kernels have also been formulated for complex input [Van der Wilk et al., 2017 Novak et al. 2020 Wilson et al. 2016)
Balandat al. (2020); Maddox et al.
spaces including convolutional kernels
and graph kernels (Shervashidze et al.
2011 Walker & Glocker 2019). The graph kernels have been used to
apply BO to neural architecture scarch (NAS) where the architecture and connectivity of a neural network
itself can be optimized (Ru et. al. 2021
Deep learning has been used as scalable and flexible surrogate model for BO. In particular, Snoek et al. (2015) uses neural networks as an adaptive basis function for Bayesian linear regression, which allows BO scale to large datasets. This approach also enables BO in more complex settings including transfer learning of the adaptive basis across multiple tasks, and modcling of auxiliary signals to improve performance (Perrone et al. 2018 Additionally, Bayesian neural networks (BNNs) that use Hamiltonian Monte Carlo to sample the posterior have been used for single-task -ask and multi-task BO for hyperparameter optimization (Springenberg
ct al. 2016
Another popular approach for BO on high-dimensional spaces latent-space approaches. Here, an autoen- coder such a VAE is trained on a dataset to create a continuous latent space that represents the data, From here, more conventional optimization algorithm, such BO using GPs, can be used to optimize over the continuous latent space. This approach has been applied to complex tasks such arithmetic ex- pression optimization and chemical design Kusner et al. 2017 Gómez-Bombarelli et al. 2018 Griffiths & Hernández--Lobato 22200 Tripp et al. 2020 Deshwal & Doppa 2021j. Note that these approaches focus on both data generation and optimization simultancously, whereas our work focuses on just the optimization
process.
Random forests have also been used for iterative optimization such as sequential model-basee algorithm configuration (SMAC) as they do not face scaling challenges (Hutter et al. 2011). Treo-structured Parzen Estimators (TPE) are also popular choice for syper-parameter tuning (Bergstra et al. 2013). However, these approaches still face the same issues with encoding complex, structured inputs such as images and graphs.
Deep learning has also been applied to improve tasks other than BO. For example, active learning is a similar scheme to BO that, instcad of optimizing an objective function, aims to optimize the predictive ability of model with few data points as possible. The inductive biases of neural networks has enabled active learning on a variety of high-dimensional data including images (Gal et al.] 2017), language (Siddhant & Lipton 2018), and partial differential equations (Zhang et al. 2019ah. BNNs have also been applied to the contextual bandits problem, where the model chooses between discrete actions to maximize expected reward
[Blundell et al. 2015 Riquelme et al. 2018).
2 Bayesian Optimization
2.1
Prerequisites
Now, we briefly introduce the BO methodology: more details can be found in the literature (Brochu et al., 2010 Shahriari et al. 2015 Garnett 2022). We formulate our optimization task maximization problem
3

---

Published in Transactions on Machine Learning Research (09/2022)
in which we wish to find the imput x" X that maximizes some function f such that x" arg maxx fx). The input z is most simply a real-valued continuous vector, but can be generalized to categorical variables. images, or even discrete objects such as molecules. The function f returns the value of the objective f(z) (which we also refer to as the "label" z), and can represent some performance metric that we wish to maximize. In general f can be a noisy function.
A key ingredient in BO the surrogate model that produces a distribution of predictions as opposed to a single point estimate for the prediction. Such surrogate models are ideally Bayesian models, but in practice, variety of approximate Bayesian models even frequentist (i.e. empirical) distributions have been used. In iteration N, a Bayesian surrogate model M is trained on labeled dataset Perain ((Xn:Bn))ee An acquisition function then uses M to suggest the next data point XN+1 X to label, where
XN+1 = arg max a M, Derain). xEX'
(1)
The new data is evaluated to get WN+1 f(xN+1), and (XN+1, UN+1) is added to Dursin
2.2 Acquisition Function
An important consideration within BO is how to choose the next data point X E X given the model M and labelled dataset Dirain- This is parameterized through the "acquisition function" " a, which we maximize to get the next data point to label as shown Equation
We choose the expected improvement (EI) acquisition function ogr (Jones et al. 1998 When the postcrior
predictive distribution of the surrogate model M is normal distribution N (N((x), of (x)). analytically
El can be expressed
(x) = (x)
((7(x))),
(2)
where 7(x) = (u(x) - //res)//(()) Ybest = max( (onx((nn))) ) is the best value of the objective function so far, and á and &b are the PDF and CDF of the standard normal (O0 1), respectively, For surrogate models that do not give analytical form for the posterior predietive distribution, sample from the posterior NMc times and use Monte Carlo (MC) approximation of EI:
NMc
OBLMM(() NMC
Σ max (x) - -Ybestt 0)
(3)
where ( is a prediction sampled from the posterior of M Wilson et al. 2018 While works such as Lakshminarayanaan et al. (2017) fit the output of the surrogate model a Gaussian use Eq. 2 for acquisition, this is not valid when the model prediction for y not Gaussian, which is generally the case for composite functions (see Section 2.4)
EI has the advantage over other acquisition functions in that the MC approximation (1) remains differentiable to facilitate opiimizaiion of the acquisition function in the inner loop (i.e. the MC approximation ion of upper confidence bound (UCB) is not differentiable and can result in ties) and (2) is inexpensive (i.e. naive Thompson sampling for ensembles would require re-training a model from scratch in each iteration).
2.3 Continued Training with Learning Rate Annealing
One challenge is that training surrogate model Dirain from scratch in every optimization loop adds large computational cost that limits the applicability of BO, especially since neural networks ideally trained for long time until convergence. To minimize the training time of BNNs in each optimization loop, we use the model that has been trained in the Nth oplimization loop iteration as the initialization ion (also known as a "warm start") for the (N + 1)th iteration, rather than training from a random initialization. In particular, we use the cosine anncaling learning rate proposed in Loshchilov Hutter (2016) which starts with a large learning rate and drops the learning rate to O. For more details, refer to Section A.3 in the A Appendix.
4

---

Published in Transactions on Machine Learning Research (09/2022)
2.4 Auxiliary Information
Typically we assume I' is black box function, we train M: X + y to model f. Here we consider the case where the experiment or observation may provide some intermediate auxiliary information z E 2, such that f can be decomposed as
f(x) = h(g(x)),
(4)
where g: X + E is the expensive labeling process, and h: 2 + y is a known objective function that can be cheaply computed. Note that this is also known as "composite functions" (Astudillo & Frazier 2019 Balandat et al. 2020 In this case, we train M: X + Z to model 9, and the approximate El acquisition function bocomes
NMc
Σ max (n (x())x))) - V/best, 0).
O'I-MC-aux (x) NMC
(5)
which can be seen as Monte Carlo version of the acquisition function presented in AAtudilloo & Frazier (201)). We denote models trained using auxiliary information with the suffix "-aux?" Because is not necessarily linear, (x)) is not general Gaussian even if (() itself may be, which makes the MC approximation convenient or even necessary.
3 Surrogate Models
Bayesian models are able to capture uncertainty associated with both the data and the model parameters in the form of probability distributions. To do this, there is a prior probability distribution P(0) placed upon the model parameters, and the posterior belief of the model parameters can be calculated using Bayes' theorem upon observing new data. Fully Bayesian neural networks have been studied in small architectures, but are impractical for realistically-sizer neural networks as the nonlinearities between layers render the posterior intractable, thus requiring the use of MCMC methods to sample the posterior. In the last decade, however, there have been numerous proposals for approximate Bayesian neural networks that are able to capture some of the Bayesian properties and produce a predictive probability distribution. In this work, we compare several different options for the BNN surrogate model. In addition, we compare against other non-BNN basclines. We list some of the more notable models here, and model details and results can be found in Section A.4.1 of the Appendix.
Ensembles combine multiple models into model to improve predictive performance by averaging the results of the single models Ensembles of neural networks have been reported to be more robust than other BNNs (Ovadia et al. 2211),, and we use "ENSEMBLE" to denote an ensemble of neural networks with identical architectures but different random initializations, which provide enough variation for the individual models to give different predictions. Using the individual models can be interpreted as sampling from posterior distribution, and so we use Eq. for acquisition. Our ensemble size NMC = 10.
Other BNNS. We also compare to variational BNNs including Bayes by Backprop (BBB) Blundell et al. 2015 and Multiplicative Normalizing Flows (MNF) Louuioos & Welling 2017]; BOHAMIANN (Springenberg et al. 2016): ; and NEURAILLINAR (Sooek et al. 2015 For BBB, we also experiment with KT. annealing, denoted by "-Anneal."
GP Baselines. GPs are largely defined by their kernel (also called "covariance functions") which determines the prior and posterior distribution, how different data points relate to each other, and the type of data the GP can operate on. In this work, we will use "GP" to refer a specific, standard specification that uses a Matérn 5/2 kernel, a popular kernel that operates over real-valued continuous spaces. To operate on images, use convolutional kernel, labeled as "ConvGP" which is implemented using the infinite-width limit of convolutional neural network (Novak et al. 2020). Finally, to operate directly on graphs, we use the Weisfeiler-Lehman (WL) kernel as implemented by (Ru et. al., (021), which we label as "GRAPHGP". The WL kernel is able to operate on undirected graphs containing node and edge features making it appropriate for chemical molocule graphs, and was used by Ru et al. (2021) to optimize neural network architectures a method they call NAS-BOWI, Additionally, we compare against "GP-AUX" which use multi-output GPs for problems with auxiliary information (also known composite functions) (Astudillo & Frazier 2019). Tn
5

---

Published in Transactions on Machine Learning Research (09/2022)
(a)
(b)
(c)
(d)
Layer thicknesses
Layer thicknesses X E I2
cross-seciion aftr"
x R
Scattering
3.5
GP
BNN-aux
30
Scattering spectrum 22
X1
X3
Wavelength [nm]]
Objective y R
Objective y E R
Figure 1: (a) A cross-section of three-layer nanoparticle parameterized by the layer thicknesses. (b) An example of the scattering cross-soction spectrum of six-layer nanoparticle. (c) Whereas GPs are trained to directly predict the objective function, (d) multi-oulput BNNs can be trained with auxiliary information, which here is the scattering spectrum.
the Appendix, we also look at GPs that use infinite-width and infinite-ensemble neural network limits as the kernel (Novak et al. 2020 as well as TuRBO which combines GP-based BO with trust regions (Eriksson et al. 2019
unlabelled dataset representative of X. This allows to encode complex input spaces, such as chemical molecules, into a continuous latent space over which con- ventional GP-based BO methods can be applied, even enabling generation and discovery of novel molecules that were not contained in the original datasct. Here, we modifiod the implementation provided by (Tripp et al. 2020) which they use a junction tree VAE (JTVAE) to encode chemical molecules (Jin et al. 2018 More details can be found in the Appendix.
VAE-GP uses a VAE trained ahead of time
Other Basclines. We compare against two variations of Bayesian optimization, TuRBO Eriksson et al, 2019) and TPE (Bergstra et al. 2013). We also compare against several global optimization algorithms that do not use surrogate models and are cheap to run, including LIPO (Malherbe & Vayatis 2017 DIRECT-L (Gablonsky & Kelley 2001), and CMA-ES.
We emphasize that ensembles and variational methods can easily scale up to high-dimensional outputs with minimal increase in computational cost by simply changing the output layer size. Neural Linear and GPs scale cubically with output dimensionality (without the use of covariance approximations), making them difficult to train on high-dimensional auxiliary intermediate information.
4 Results
We now look at three real-world scientific optimization tasks all of which provide intermediate or auxiliary information that can be leveraged. In the latter two lasks, the structure of the data also becomes important and hence BNNs with various inductive biases significantly outperform GPs and other baselines. For sim- plicity, we only highlight results from select architectures (sce Appendix for full results along with datasct and hyperparameter details). All BO results are averaged over multiple trials, and the shaded area in the plots represents + one standard error over the trials.
4.1 Multilayer Nanoparticle
We first consider the simple problem of light scattering from a multilayer nanoparticle, which has wide variety of applications that demand a tailored optical response (Ghosh Chaudhuri & Paria, 2012) including biological imaging (Saltsberger et al.) 2012). improved solar cell efficiency (Ho et al., 2012 Shi et al. 2013) and catalytic materials (Tang & Henkelman 2200). In particular, the nanoparticle we consider consists of lossless silica core and spherical shells of alternating TiO, and silica. The nanoparticle is parameterized by the core radius and layer thicknesses shown in Figure [a], which restrict. to the range 30 nm to

---

Published in Transactions on Machine Learning Research (09/2022)
Narrowband
Highpass
0.1600
Random
1.05
Random
0.1575
Veest
0.1550
GP-aux
BBB-Anneal-aux
1.02
GP-aux
Ensemble
0.1525
BBB-Anneal
Ensemble-aux
0.98
0.1475
0.36
300
400
800
800
1000
200
400
500
800
1000
N
Figure 2: BO results for two different objective functions for the nanoparticle scattering problem. Training with auxiliary information (where M is trained to predict z) is denoted with "-aux". Adding auxiliary information to BNNs significantly improves performance.
70 nm. Because the size of the nanoparticle is on the order of the wavelength of light, its optical properties can be tuned by the number and thicknesses of the layers. The scattering spectrum can be calculated semi-analytically, as detailed in Section A.1.T of the Appendix.
We wish to optimize the scattering cross-section spectrum over a range of visible wavelengths, an example of which is shown in Figure [bb. In particular, we compare two different objective functions: the narrowband objective that aims to maximize scattering in the small wavelength range 600 to 640 nm and minimize it elsewhere, and the highpass objective that aims to maximize scattering above 600 nm and minimize it elsewhere. While conventional GPs train using the objective function as the label directly, BNNs with auxiliary information can be trained to predict the full scattering spectrum, i.e. the auxiliary information 2 E 2201 which is then used to calculate the objective function, as shown in Figure (e.d).
BO results are shown Figure 2 Adding auxiliary information significantly improves BO performance for BNNs. Additionally, they are competitive with GPs, making BNNs a viable approach for scaling BO large datasets. In Appendix A.5] we see similar trends for other types of BNNs. Due to poor scaling of multi-output GPs with respect to output dimensionality,
are only able to run GP-AUX for small number of iterations in a reasonable time. Within these few iterations, GP-AUX performs poorly, only slightly better than random sampling. We also see in the Appendix that BO with either GPs or BNNs are comparable with, or outperform other global optimization algorithms, including DIRECT-L and CMA-ES.
4.2 Photonic Crystal Topology
Next we look at a more complex, high-dimensional domain that contains symmetries not easily exploitable by GPs. Photonic crystals (PCs) are nanostructured materials that are engineered to exhibit exotic optical properties not found in bulk materials, including photonic band gaps, negative index of refraction, and angular selective transparency (John 1987: Yablonovitch 1987 Joannopoulos et al. 2008 Shen et al. 2014). As advanced ion techniques are enabling smaller and smaller Teature sizes, there has been growing interest in inverse design and topology optimization to design even more sophisticated PCs (Jensen & Sigmund] 2011 1: Men et al.] 2014) for applications in photonic integrated circuits, flat lenses, and sensors
Piggott et al. 2015: Lin et al. 2019
Here we consider 2D PCs consisting of periodic unit cells represented by a 32 pixel image, as shown in Figure 3(a), with white and black regions representing vacuum (or air) and silicon, respectively. Because optimizing over raw pixel values may lead to pixel-sized features or intermediate pixel values that cannot be fabricated, we have parameterized the PCs with level-set function Φ: X V that converts 51-dimensional feature vector x = [c1, C2,..., Coo, A] E IR$1 representing the level-set parameters into an image v E 3232332 that represents the PC. More details can be found in Section A.1.2 in the Appendix.
We test BO on two different data distributions, which are shown Figure B (bcc). In the PC-A distribution, x spans e, E [-1,1], , E [-3,3], 3). In the PC-B distribution, we arbitrarily restrict the domain c; E [0, 1]. The PC-A data distribution is translation invariant, meaning that any PC with translational shift will also
7

---

Published in Transactions on Machine Learning Research (09/2022)
(a)
(b)
(e)
(f)
Level set parameters R
Level set parameters 12
Unit cell image v E Q
(c)
(d)
GP
BCNN-aux
20 8
Density of States z E 5000
10-
frequency
Objective y IR
Objective y E
Figure 3: (a) A 2D photonic crystal (PC). The black and white regions represent different materials, and the periodic unit cells are outlined in red. Examples of PC unit cells drawn from the (b) PC-A distribution and (c) the PC-B distributions. The PC-A data distribution translation invariant, whereas unit cells drawn from the PC-B distribution all have white regions in the middle of the unit cell, so the distribution is not translation invariant.
(d) Example of a PC density of states (DOS). (e, f) Comparison of the process flow for training the surrogate model in the case of (e) GPs and (f) Bayesian Convolutional NNs (BCNN), The BCNN can train directly on the images to take advantage of the structure and symmetries in the data, and predict the multi-dimennsinal DOS.
be in the data distribution. However, the PC-B data distribution is not translation invariant, as shown by the white regions in the center of all the examples in Figure Te).
The optical properties of PCs can be characterized by their photonic density of states (DOS), c.g. see Figure Bd)), We choose an objective function that aims to minimize the DOS in a certain frequency range while maximizing it everywhere else, which corresponds to opening up a photonic band gap in said frequency range. As shown in Figure 3(e,f,, we train GPs directly on the level-set parameters X', whereas we train the Bayesian convolutional NNs (BCNNs) on the more natural unit cell image space V. BCNNs can also be trained to predict the full DOS as auxiliary information E R500
The BO results, seen in Figure aa), show that BCNNs outperform GPs by a significant margin both datasets, which is due to both the auxiliary information and the inductive bias of the convolutional layers, as shown in Figure b). Because the behavior of PCs is determined by their topology rather than individual pixel values or level-sct parameters, BCNNs are much better suited to analyze this datasct compared to GPs. Additionally, BCNNs can be made much more data--fficient since they directly encode translation invariance and thus learn the behavior of a whole class of translated images from a single image. Because GP-AUX is extremely expensive compared to GP (500x longer on this datasct), we only able to run GP-AUX for small number of iterations, where it performs comparably to random sampling. We also compare to GPs using a convolutional kernel (s ("ConvGP-NNGP") in Figure [a). ConyGP-NNGP only performs slightly better than random sampling, which is likely due to a lack of auxiliary information and inflexibility to learn the most suitable representation for this dataset.
For our main experiments with BCNNs, we use an architecture that respects translation invariance. To demonstrate the effect of another commonly used deep learning training technique, we also experiment with incorporating translation invariance into translation dependent (i.e. not translation invariant) architecture using a data augmentation scheme in which each image is randomly translated, flipped, and rotated during training. We expect data augmentation to improve performance when the data distribution exhibits the corresponding symmetries: in this case, we focus on translation invariance. As shown in Figure c), we indeed find that data augmentation improves the BO performance of the translation dependent architecture when trained on the translation invariant PC-A datasot, even matching the performance of a translation invariant architecture on PC-A. However, on the translation dependent PC-B dataset, data augmentation initially hurts the BO performance of the translation dependent architecture because the model is unable to
8

---

Published in Transactions on Machine Learning Research (09/2022)
(a)
(b)
(c)
(d)
4000
3000
PC-A
2000
1000
Yoest
PC-B
3000
2000
1000
250
500
750
1000
250
500
750
1000
250
500
750
1000
0
250
500
750
1000
N
Rand
Ensemble
ConvBBB
T1
Ensemble-aux
ConvBBB-aux
TD TD-augment
ConvGP-NNGP
CorvEnsemble
NeuralLinear
GP-aux ConvEnsemble-aux
ConvEnsemble-add
ConvNeuralLinear
DIRECT-L
Figure 4: Three sets of comparisons for BO results on the (top row) PC-A and (bottom row) PC-B datasets. (a) BNNs with inductive biases outperform all other GP baselines and the random baseline. Note that GP-AUX is comparable to random sampling.
(b) The inductive bias of convolutional layers and the addition of suxiliary information siguificantly improve performance BCNNs. (c) Additional comparisons. (d) Data augmentation boosts performance if the augmentations reflect symmetry present in the dataset but not enforced by the model architecture. "TT" refers to translation invariant BCNN architecture, whereas "TD" refers to translation dependent architecture. "-augment" signifies that data augmentation of the photonic crystal image is applied, which includes periodic translations, flips, and rotations.
quickly specialize the more compact distribution of PC-B, putting BO performance more on par with models trained on PC-A. These results show that techniques used improve generalization performance (such as data augmentation or invariant architectures) for training deep learning architectures can also be applied to BO surrogate models and, when used appropriately, directly translate into improved BO performance. Note that data augmentation would not be feasible for GPs without a hand-crafted kernel as the increased size of the dataset would cause inference to become computationally intractable.
4.3 Organic Molecule Quantum Chemistry
Finally, we optimize the chemical properties of molecules. Chemical optimization is of significant interest in both academia and industry with applications in drug design and materials optimization (Hughes al. 2011). This is difficult problem where computational approaches such as density functional theory (DFT) can take days for simple molecules and are intractable for larger molecules; synthesis is expensive and time- consuming, and the space of synthesizable molecules is large and complex. There have been many approaches for molecular optimization that largely revolve around finding continuous latent space of molecules (Gómez- Bombarelli et al. 2018) or hand-crafting kernels to operate on molecules (Korovina et al. 2020).
Here we focus on the QM9 dataset (Ruddigkeit et al. 2012 Ramakrishnan et al. 2014). which consists of 133,885 small organic molecules along with their geometric, electronic, and thermodynamics quantities that have been calculated with DFT. Instead of optimizing of available molecules and iteratively select the next molecule to add to D This is problem setting especially common to materials design where databases are incomplete and the space of experimentally- feasible materials is small.
continuous space, draw from the fixed pool
We use a Bayesian graph neural network (BGNN) for our surrogate model, as GNNs have bocome popular for chemistry applications due the natural encoding of molecule as graph with atoms and bonds nodes and edges, respectively. For baselines that operate over continuous spaces (i.e. GPs and simple neural

---

Published in Transactions on Machine Learning Research (09/2022)
(a)
(b)
(c)
Egep
Molecule
Molecule
160
0.8
150
140
Random
130
GP
GraphEP
SOAP descriptor R2
Molecule graph
120
VAE-GP
X E Ruoka, A E mumm E E gooomsss
110
93-
nremhinsembic
160
GraphEnsembiieaax
02
9.1
BGNN-aux
100
220
100
600
540
100
200
400
500
properties E 12%
Time per iteration [+]
GraphGP GraphiPr Popl-1et GraphEnsemble GraphEnserble Pasfmled
Objective yER
Objective y e IR
10*
100 200
Figure Quantum chemistry task and results. (a) The GP is trained on the SOAP descriptor, which is precomputed for each molccule. (b) The BGNN operates directly graph representation of the molccule, where atoms and bonds are represented by nodes and edges, respectively. The BGNN can be trained on multiple properties given in the QM9 dataset. (c) BO results for various properties. Note that GraphEnsem- ble is a type of BGNN. (d) Time per BO iteration on the GPU. (Note the logarithmic scale on the y-axis.) GRAPHGP takes orders of magnitudes longer than BGNNs for moderate N.
networks), we use the Smooth Overlap of Atomic Positions (SOAP) descriptor to produce a fixed-length feature vector for each molecule, as shown in Figure 5(a) De et al. 2016 Himanen el. al. 2020
We compare two different optimization objectives derived from the QM9 dataset: the isotropic polarizability a and (a (gap) where Egap is the HOMO-LUMO energy gap. Other objectives are included in Appendix A.5.4 Because many of the chemical properties in the QM9 dataset can be collectively computed by single DFT molecular dynamics calculation, we treat group of labels from QM9 auxiliary information and train BGNN to predict this entire group simultaneously. The objective function h then simply picks out the property of interest.
As shown in Figure 5 c), GRAPHGP and the BGNN variants significantly outperform GPs, showing that the inductive bias the graph structure leads to a much more natural representation the molecule and properties. In the case of maximizing the polarizability a; including the auxiliary information improves BÖ performance, showing signs of positive transfer. However, it does not have a significant impact on the other objectives, which may be due to the small size of the available auxiliary information (only handful of chemical properties from the QM dataset) as compared with the nanoparticle and photonic crystal tasks. In more realistic online setting, we would have significantly more physically-ine innormaalii information available from a DFT calculation, e.g. we could easily compute the electronic density of states (the electronic analogue of the auxiliary information used the photonics task).
As scen Figure d), we also note that the GRAPHGP relatively computationally expensive (15x longer than GPs for small N and 800x longer than BGNNs for N 100) and so we are only able run it for limited N in reasonable time frame. We see that BGNNs perform comparably or better than GRAPHGPs despite incurring fraction of the computational cost.
VAE-GP uses a modified version of the latent-space optimization method implementation provided by Tripp et al. (2020). Rather than optimizing over a contimuous latent space of the VAE, we feed the data pool through the VAE encoder to find their latent space representation, and then apply the acquisition function to the latent points to pick out the best unlabeled point to sample. We keep as many hyper- parameters the same the original implementation possible, with the exception of the weighted re-
10

---

Published in Transactions on Machine Learning Research (09/2022)
training which we forgo since we have fixed data pool that was used to train the VAE. This setup is similar to GRAPHNRUUA.LLNNAN in that deep learning architecture is used to encode the molecule as a continuous vector, although GRAPHNEURALLINAA is only trained on the labelled data. The results for this experiment show that VAE-GP performs worse than BNNs on two of the three objective functions we tested and slightly better on one objective. We also note that the performance of VAE-GP depends very heavily on the pre-training of the VAE, as choosing different hyper-parameters or even different random seed can significantly deteriorate performance (see Figure 15 in the Appendix).
5 Discussion
Introducing physics-informed priors (in the form of inductive biases) into the model are critical for their performance. Well-known inductive biases in deep learning include convolutional and graph neural networks for images and graph structures, respectively, which significantly improve BO performance. Another indue- tive bias that we introduce is the addition of auxiliary information present in composite functions, which significantly improves the performance of BQ for the nanoparticle and photonic crystal tasks. We conjecture that the additional information forces the BNN to learn a more consistent physical model of the system since it must learn features that are shared across the multi-dimennional auxiliary information, thus enabling the BNN to generalize better. For example, the scattering spectrum of the multilayer particle consists of mul- tiple resonances (sharp peaks), the width and location of which are determined by the material properties and layer thicknesses. The BNN could potentially learn these more abstract features, and thus, the deeper physics, to help it interpolate more efficiently, akin to data augmentation (Peuriloy et al.] 2018). Auxiliary information can also be interpreted as a form of data augmentation. Indeed, tracking the prediction error on a validation set shows that models with auxiliary information tend to have a lower loss than those without (see Appendix A.5). It is also possible that the loss landscape for the auxiliary information is smoother than that of the objective function and that the auxiliary information acts as an implicit regularization that improves generalization performance.
Interestingly, GP-AUX performs extremely poorly on the nanoparticle and photonic crystal tasks. One possible reason is that we are only able to run GP-AUX for a few iterations, and it is not uncommon for GP-based BÖ to require some critical number of iterations to reach convergence especially in the case of high-dimensional systems where the size of the covariance matrix scales with the square of the dimensionality. It may also be possible that GP-AUX only works on certain types of decompositions of functions and cannot be applied broadly to all composite functions, as the inductive biases in GPs are often hard-coded.
There is an interesting connection between how well BNNs are able to capture and explore a multi-modal posterior distribution and their performance in BO. For example, we have noticed that larger batch sizes tend to significantly hurt BO performance. On the one hand, larger batch sizes may be rosulting in poorer generalization as the model finds sharper local minima in the loss landscape. Another explanation is that the stochasticity inherent in smaller batch sizes allows the BNN to more easily explore the posterior distribution, which is known to be highly multi-modal (Fort al. 2 Indeed, BÖ often underperforms for very small dataset sizes N but quickly catches up as N increases, indicating that batch size is an important hyperparameter which must be balanced with computational cost.
All our results use continued training (or warm restart) to minimize training costs. We note that re-initializing M and training from scratch in every iteration performs better than continued training on some tasks (results the Appendix), which points to how BNNs may not sufficiently represent a multi-modal posterior distribution or that continued training may skew the training distribution that the BNN sces. Future work will consider using stochastic training approaches such as SG-MCMC methods for exploring posterior distributions (Welling& & Teh 2011 Zhang et al., 2019b) well as other continual learning techniques to further minimize training costs, especially for larger datasots (Parisi et al. 2019).
When comparing BNN architectures, we find that ensembles lend to consistently perform among the best, which is supported by previous literature showing that ensembles capture uncertainty much better than variational methods (Ovadia et al. 2019; Gusstffson et al.] 2020) espocially in multi-modal loss landscapes Fortt et al.] (2019). Ensembles are also attractive because they require no additional hyperparameters and they are simple to implement. Although training costs increase linearly with the size of the ensemble, this
11

---

Published in Transactions on Machine Learning Research (09/2022)
can be easily parallelized on modern computing infrastructures. Furthermore, recent work that aims to model efficient, ensembles that minimize computational cost could be an interesting future direction (Havasi et al. 2020 Wen et al.) 2020j. NEURALLINEAR variants are also quite powerful and cheap, making them very promising for tasks without high-dimensional auxiliary information. Integrating Neural Linear with multi-oulpul GPs is an interesting direction for future work. The other BNNs either require extensive hyper-parameter tuning or perform poorly, making them difficult to use in practice. Additional discussion can be found in Appendix A.5.5
As seen in Appendix A.5.4 VAE-GP performs worse than our method on two of the chemistry objectives and better on one objective. While latent-space optimization methods are often applied to domains where one wants to simultaneoushy generate data and optimize over the data distribution, these methods can also be applied to the cases in this work, where a data pool (e.g. QM9 dataset for the chemistry task) or separate data generation process (e.g. level-set process for the photonic crystal task) is already available. In these cases, the VAE is not used as a generative model, but rather as a way to learn appropriate representations. While latent- space approaches are able to take advantage of well-developed and widely available optimization algorithms, they also require unsupervised pre-training the necessary inductive biases. Such models are available in chemistry where there have been significant development, but are more limited in other domains such as photonics. On the other hand, our method is able to incorporate the data structure or domain knowledge in an end-to-end manner during training, although future work is needed to more carefully evaluate how much of an advantage this is and whether it depends on specific dataset or domain characteristics. For settings where we do not need a generative model, it would also be interesting to replace the autoencoder with self-supervised model (Hendrycks et al., 2019: semi-supervised model (Kingma et al. 2014) to create suitable latent space.
sizeable dataset and a suitable autoencoder model with
Loh al. 2021
6 Conclusion
We have demonstrated global optimization on multiple tasks using a combination of deep learning and BO. In particular, we have shown how BNNs can be used as surrogate models in BO, which enables the scaling of BO to large datasets and provides the flexibility to incorporate a wide variety of constraints, data augmentation techniques, and inductive biases. We have demonstrated that integrating domwin-kkoowledge on the structure and symmetries of the data into the surrogate model as well as exploiting intermediate or auxiliary information significantly improves BO performance, all of which can be interpreted as physics- informed priors. Intuitively, providing the BNN surrogate model with all available information allows the BNN to learn a more faithful physical model of the system of interest, thus enhancing the performance of BO. Finally, we have applied BO to real-world, high-dimensional scientific datasets, and our results show that BNNs can outperform our best-effort GPs, even with strong domain-dependent structure encoded in the covariance functions. We note that our method is not necessarily tied to any particular application domain, and can lower the barrier of entry for design and optimization.
Future work will investigate more complex BNN architectures with stronger inductive biases. For example, output constraints can be placed through unsupervised learning Karpatne et al. 2017) or by variationally fitting a BNN prior (Yang et al. 2020). Custom architectures have also been proposed for partial differential
equations (Raissi et al.) 2017 Lu et al. 2020), body systems (Cranmer et al. 2020), and generalized symmetries (Hutchinson et al. 20200) which will enable effective BO on a wider range of tasks. The methods and experiments presented here enable BO to be effectively applied in wider variely of settings. There are also variants of BO including TuRBO which perform extremely well on our tasks, and so future work will also include incorporating INNs into these variants.
We make our datasets and code publicly available at https: httpa//w...oooommmmeeee/ame///ee
Acknowledgements
The authors would like acknowlodge Rodolphe Jenatton, Thomas Christensen, Andrew Ma, Rumen Dan- govski, Joy Zeng, Tin D. Nguyen, Charles Roques-Carmes, Carmes, and Mohammed Benzaonia for fruitful conversa- tions. The authors acknowledge the MIT SuperCloud and Lincoln Laboratory Supercomputing Center for
12

---

Published in Transactions on Machine Learning Research (09/2022)
providing HPC resources that have contributed to the research results reported within this paper. This work is supported in part by the the National Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, http::/iaaiii research was also sponsored in part by the Department of Defense through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program. This material is based upon work partly sup- portod by the Air Force Office of Scientific Research under the award number FA9550-21-1-0317, as well partly supported by the US Office of Naval Research (ONR) Multidisciplinary University Research Initiative Robust Photonic Materials with High-Order Topological Protection It
This
(MURI) grant NO -01-220--22325 is also based upon work supported in part by the U.S. Army Research Office through the Institute for Sol- dier Nanotechnologies at MIT, under Collaborative Agreement Number W911NN-18-2-00445 Research was sponsored by the United States Air Force Research Laboralory and the United States Air Force Artificial Intelligence Accelerator and was accomplished under Cooperative Agreement Number FA8750-19-2- 2-000 The views and conclusions contained in this document are those of the authors and should not be inter- preted as representing the official policies, either expressed implied, of the United States Air Force or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.
References
Raul Astudillo and Peter Frazier. Bayesian optimization of composite functions. In International Conference on Machine Learning, pp. 354 363. PMLR, 2019.
Maximilian Balandat, Brian Karrer, Daniel Jiang, Samuel Daulton, Ben Letham, Andrew G Wilson, and Eytan Bakshy. Botorch: a framework for efficient monte-carlo bayesian optimization. Advances in newral information processing systems, 33:21524- 2--5522--1538, 2020.
James Bergstra, Daniel Yamins, and David Cox. Making a science of model search: Hyperparameter op- timization in hundreds dimensions for vision architectures. In International conference on machine learning, pp. 123. PMLR, 2013.
Charles Blundell, Julien Cornebise, Koray Kavukeuoglu, and Daan Wierstra. Weight uncertainty in neural networks. Xiv preprint 154...5440, 2015.
Eric Brochu, Vlad M Cora, and Nando De Freitas. A futorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning. arXiv preprint arXiv:1012.2500 2010.
Wessel Bruinsma, Eric Perim, William Tebbutt, Scott Hosking, Arno Solin, and Richard Turner. Scalable exact inference in multi-output gaussian processes. In International Conference on Machine Learning, pp. 1190 1201. PMLR, 2020.
Miles Cranmer, Alvaro Banchez-Gonzalee Peter Battaglia, Rui Xu, Kyle Cranmer, David Spergel, and Shirley Ho. Discovering symbolic models from deep learning with inductive biases. ar preprinl arXiv:206..1 11287, 2020.
Sandip De, Albert P Bartók, Gábor Csányi, and Michele Ceriouti. Comparing molecules and solids across structural and alchemical space. Physical Chemistry Chemical Physics, 18(20):13754-137000 2016.
Aryan Deshwal and Jana Doppa. Combining latent space and structured kernels for bayesian optimization over combinatorial spaces. Advances in Neural Information Processing Systems, 34, 2021.
David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local bayesian optimization. Advances in neural information processing systems, 32, 2019.
Stanislav Fort, Huiyi Hu, and Balaji Lakshminarayanan. Deep ensembles: A loss landscape perspective. arXiv preprint ar
2019.
Joerg M Gablonsky and Carl T Kelley. A locally-biased form of the direct algorithm. Journal of Global Optimization, 21(1):77335 37, 2001.
13

---

Published in Transactions on Machine Learning Research (09/2022)
Yarin Gal, Riashat Islam, and Zoubin Ghahramani. Deep bayesian active learning with image data. arXiv preprinl Xiv: 1700. 02910, 2017.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. Advances in newral information processing systems, 31, 2018.
Roman Garnett. Bayesian Optimization. Cambridge University Press, 2022. in preparation.
Rajib Ghosh Chaudhuri and Santamu Paria. Core/shell nanoparticles: classes, properties, synthesis mecha- nisms, characterization, and applications. Chemical revicus, 112(4))2373-2438 2012.
Rafael Homezz-oombarell. Jennifer N Wei, David Duvenaud, José Miguel Hernández-Lobato, Benjamín Jänchez-Lengeling Dennis Sheberla, Jorge Aguilera-Iparraguinng Timothy D Hirzel, Ryan P Adams, and Alán Aspuru-Guzik. Automatic chemical design using data-driven continuous representation of molecules. ACS central science, 4(2):268-27) 2018.
Daniele Grattarola and Cesare Alippi. Graph neural networks in tensorflow and keras with spektral. arXiv preprint Xix::2006. 12138, 2020.
Ryan-Rhys Griffiths and José Miguel HeráándozLLobato. Constrained bayesian optimization for automatic chemical design using variational autoencoders. Chemical science, 11(2):577 586, 2020.
Chuan Guo, Gcoff Plciss, Yu Sun, and Kilian Q Weinberger. On calibration of modern neural networks. In International conference on machine learning, pp. 1321-1330. PMLR, 2017.
Fredrik K Gustalsson, Martin Danelljan, and Thomas B Schon. Evaluating scalable bayesian deep learning methods for robust computer vision. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pp. 318-319, 2020.
Nikolaus Hansen, Youhei Akimoto, and Petr Baudis. CMA-FS/pycma on Github. Zenodo, D01::00.2811/eeooo02255333 25559644,, February 2019. URL https://doi.org// http:///oo...ogg111.52222zzoooo 2559634
Marton Havasi, Rodolphe Jenatton, Stanislav Fort, Jeremiah Zhe Liu, Jasper Snook, Balaji Lakshmi- narayanan, Andrew Dai, and Dustin Tran. Training independent subnetworks for robust prediction. arXiv preprinl ar xrXi:220000000 0666100, 2020.
Dan Hendryeks, Mantas Mazeika, Saurav Kadavath, and Dawn Song. Using self-supervised learning can improve model robustness and uncertainty. Advances in neural information processing systems, 32, 2019.
Lauri Himanen, Marc O. J. Jäger, Fiaki V. Morooka, Filippo Federici Canova, Yashasvi S. Ranawat, David Z. Gao, Patrick Rinke, and Adam S. Foster. DScribe; Library of descriptors for machine learning in materials science. Computer Physics Communications, 247:106949, 2020. ISSN 0010-4655. doi: 10.1016/j..pc.2019. 106949. URL https:///oi. /ttp://oi...rr/11.1111111/1/1/11 106994 106949
Chung-I Ho, Dan-Ju Yeh, Vin-Cent Su, Chieh-Hung Yang, Po-Chuan Yang, Ming-Yi Yi Pu, Chieh-Hsiung Kuan, I-Chun Cheng, and Si-Chen Lee. Plasmonic multilayer nanoparticles enhanced photocurrent in thin film hydrogenated amorphous silicon solar cells. Journal of Applied Physics, 112(2):023113, 2012.
Jilin Hu, Jianbing Shen, Bin Yang, and Ling Shao. Infinitely wide graph convolutional networks: semi-
supervised learning via gaussian processes. arXiv preprint 121668,0 2020.
Gao Huang, Yixuan Li, Geoff Pleiss, Zhuang Liu, John E Hoperoft, and Kilian Q Weinberger. Snapshot ensembles: Train 1, get m for free. arXiv preprint arXiv:1 704000109,, 2017.
James P Hughes, Stephen Rees, S Barrett Kalindjian, and Karen L Philpott. Principles of early drug discovery. British journal of pharmacology, (219))112991245 2011.
Michael Hutchinson, Charline Le Lan, Sheheryar Zaidi, Emilien Dupont, Yee Whye Tch, and Hyunjik Kim. Lietransformer: Equivariant self-attention for lie groups. arXiv preprint ar 20v: 10885, 2020.
14

---

Published in Transactions on Machine Learning Research (09/2022)
Frank Hutter, Holger H Hoos, and Kevin Leyton- Leyton-Brown. algorithm configuration. In International conference learning and intelligent oplimization, pp. 507-523 Springer, 2011.
Brown. Sequential model-based optimization for general
Jakob Sondergaard Jensen and Ole Sigmund. Topology optimization for nano-photonics. Laser & Photonics Reviews, 5(2):308-322 2011.
Wengong Jin, Regina Barzilay, and Tommi Jaakkola. Junction tree variational autoencoder for molecular graph generation. In International conference on machine learning, pp. 2323-2332. PMLR, 2018.
John D. Joannopoulos, Steven G. Johnson, Joshua N. Winn, and Robert D. Meade. Photonic Crystals: Molding the Flow of Light (Second Edition). Princeton University Press, edition, 2008. ISBN 0691124566.
Sajeev John. Strong localization of photons in certain disordered dielectric superlattices. Physical review letters, 58(23)22486, 1987.
Steven G. Johnson. nlopt
The nlopt noulhnea--optimization package, 2010. URL http:://gattbb..cm//ttttt/////
Steven G Johnson and John D Joannopoulos. Block-iterative frequency-domain methods for maxwell's equations in a planewave basis. Optics express, 8(3):173-199) 190, 2001.
Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13( ():4)) 492, Dec 1998. ISSN 1573-2916. doi:
10.1028/A:100030641114 URL https://doll org/ 10. 100012222/000030441114
Amuj Karpatne, William Watkins, Jordan Read, and Vipin Kumar. Physics-guided neural networks (pgnn): An application in lake temperature modeling. arXiv preprint arXi::7110.1143 2017.
Davis E. King. Dlib-ml: A machine learning toolkit. Journal of Machine Learning Research, 10:1755 17588 2009.
Durk P Kingma, Shakir Mohamed, Danilo Jimenez Rezende, and Max Welling. Semi-supervised learning with deep generative models. Advances in newral information processing systems, 27, 2014.
Aaron Klein, Stefan Falkner, Simon Bartels, Philipp Hennig, and Frank Hutter. Fast bayesian optimization of machine learning hyperparameters on large datasets. In Artificial Intelligence and Statistics, pp. 538-536. PMLR, 2017.
Ksenia Korovina, Sailun Xu, Kirthevasan Kandasamy, Willie Neiswanger, Barnabas Poczos, Jeff Schneider, and Eric Xing. Chembo: Bayesian optimization of small organic molecules with synthesizable recom- mendations. In International Conference on Artificial Intelligence and Statistics, pp. 2020.
PMLR,
Volodymyt Kuleshov, Nathan Fenner, and Stefano Ermon. Accurate uncertainties for deep learning using calibrated regression. In International conference machine learning, pp. 2796-2804. PMLR, 2018.
Matt J Kusner, Brooks Paige, and José Miguel Hernández-Lobato. Grammar variational autoencoder. In International conference on machine learning, pp. 1945 1954. PMLR, 2017.
Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive un- certainty estimation using deep ensembles. In Advances in neurol information processing systems, pp. 6402-6413, 2017.
Zin Lin, Victor Liu, Raphaël Pestourie, and Steven G Johnson. Topology optimization of freeform large-area metasurfaces. Optics express, 2711))15663-11775 2019.
Boyuan Liu, Steven G Johnson, John Joannopoulos, and Ling Lu. Generalized gilat raubenheimerr method for density-of-states calculation in photonic crystals. Journal of Optics, 20(4):044005, 2018.
15

---

Published in Transactions on Machine Learning Research (09/2022)
Haitao Liu, Yew-Soon Ong, Xiaobo Shen, and Jianfei Cai. When gaussian process meets big data: A review of scalable gps. IEEE transactions neural networks and learning systems, (1(1)):4405 4423, 2020.
Charlotte Loh, Thomas Christensen, Rumen Dangovski, Samuel Kim, and Marin Soljacic. Surrogate- and invariance-boosted contrastive learning for data-scarce applications in science. arXiv preprint arXiv:21 0et:21100.88. 10,00000, 2021.
Iya Loshchilov and Frank Hutter. Sgdr: Stochastic gradient descent with warm restarts. arXiv preprint arXiv:1688.0300 2016.
Christos Louizos and Max Welling. Multiplicative normalizing flows for variational bayesian neural networks. arXiv preprint ar X 707.770.100 01961, 2017.
Peter Y Lu, Samuel Kim, and Marin Soljačić. Extracting interpretable physical parameters from spatiotem- poral systems using unsupervised learning. Physical Review X, 10330003 2020.
Wesley Maddox, Maximilian Balandat, Andrew G Wilson, and Eytan Bakshy. Bayesian optimization with high-dimensional outputs. Advances in Neural Information Processing Systems, 34, 2021a.
Wesley J Maddox, Samuel Stanton, and Andrew G Wilson. Conditioning sparse variational gaussian processes for online decision-making Advances in Neurol Information Processing Systems, 34, 2021b.
Cédric Malherbe and Nicolas Vayatis. Global optimization of lipschitz functions. arXiv:170s D.XX::2033.0628, 02628, 2017.
arXiv preprint
Han Men, Karen YK Lee, Robert M Freund, Jaime Peraire, and Steven G Johnson. Robust topology optimization of chree-dimensional photonic-crystal band-gap structures. Optics express, 22(19):22632 22648, 2014.
Roman Novak, Lechao Xiao, Jiri Hron, Jaehoon Lee, Alexander A. Alemi, Jascha Sohl-Dickstein, and Samuel S. Schoenholz. Neural tangents: Fast and easy infinite neural networks in python. In International Conference on Learning Representations, 2020. URL https: t/gttbb www.o/oogoooooooonoooooooooo
Yaniv Ovadia, Emily Fertig, Jie Ren, Zachary Nado, David Sculley, Sebastian Nowozin, Joshua Dillon, Balaji Lakshminarayanan, and Jasper Snoek. Can you trust your model's uncertainty? evaluating predictive uncertainty under dataset shift. In Advances in Neural Information Processing Systems, pp. 13991 14002, 2019.
German I Parisi, Ronald Kemker, Jose L Part, Christopher Kanan, and Stefan Wermter. Continual lifelong learning with neural networks: A review. Neural Networks, 113:54-71, 2019.
Valerio Perrone, Rodolphe Jenatton, Matthias Seeger, and Cédric Archambeau. Scalable hyperparameter transfer learning. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pp. 6846-6856, 2018
John Peurifoy, Yichen Shen, Jing, Yi Yang, Fidel Cano-Renteria, Brendan G DeLacy, John Joannopou- los, Max Tegmark, and Marin Soljačíé. Nanophotonic particle simulation and inverse design using artificial neural networks. Science advances, 4(6):caar4206, 2018.
Alexander Y Piggott, Jesse Lu, Konstantinos G Lagoudakis, Jan Petykiewicz, Thomas M Babinec, and Jelena Vučković. Inverse design and demonstration of compact and broadband on-chip wavelength demultiplexer. Nature Photonics, 976):374 37777,77- 2015.
Wenjun Qiu, Brendan G DeLacy, Steven G Johnson, John D Joannopoulos, and Marin Soljačíć. Optimization of broadband optical response of multilayer nanospheres. Optics express, 20(16):18494 18504, 2012.
Ali Rahimi and Benjamin Recht. Random features for largo-scale kernel machines. Advances in neural information processing systems, 20, 2007.
16

---

Published in Transactions on Machine Learning Research (09/2022)
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations. or Nix preprint or 1ix: 1...11111106,, 2017.
Raghunathan Ramakrishnan, Pavlo O Dral, Matthias Rupp, and O Anatole Von Lilienfeld. Quantum chemistry structures and properties of 134 kilo molecules. Scientific data, 1(1): 7, 2014.
Carlos Riquelme, George Tucker, and Jasper Snoek. Deep bayesian bandits showdown: An empirical com- parison of bayosian deep networks for thompson sampling. arXiv preprint arXiv: 18X2.0888209127 2018.
Binxin Ru, Xingehen Wan, Xiaowen Dong, and Michael Osborne. Interpretable neural architecture search via bayesian optimisation with weisfeiler-lehman kernels. In International Conference on Learning Repre- senlalions, 2021. URL
Lars Ruddigkeit, Ruud Van Deursen, Lorenz C Blum, and Jean-Louis Reymond. Enumeration of 166 billion organic small molecules in the chemical universe database gdb-17. Journal of chemical information and modeling, 52 (1211:228442 2012.
S Saltsborger, I Steinberg, and Israel Gannot. Multilayer mic scattering model for investigation of intracellular structural changes in the nucleolus and cytoplasm. International Journal of Optics, 2012, 2012.
Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando De Freitas. Taking the human out of the loop: A review of bayesian optimization. Proceedings of the IEEE, 104(1):148-I1) 175, 2015.
Yichen Shen, Dexin Ye, Ivan Celanovic, Steven G Johnson, John Joannopoulos, and Marin Soljačić. Optical broadband angular selectivity. Science, 343(6178):1449 -501, 2014.
Nino Shervashidze, Pascal Schweitzer, Erik Jan Van Leeuwen, Kurt Mehlhorn, and Karsten M Borgwardt. Weisfeiler-lcmman graph kernels. Journal of Machine Learning Research, 12(9), 2011.
Yanpeng Shi, Xiaodong Wang, Wen Liu, Tianshu Yang, Rui Xu, and Fuhua Yang. Multilayer silver nanopar- ticles for light trapping in thin film solar cells, 2013.
Aditya Siddhant and Zachary C Lipton. Deep bayesian active learning for natural language processing:
Results of a large-scale empirical study. arXiv preprinl urXi::1808. 05690, 2018.
Martin Simonovsky and Nikos Komodakis. Dynamic edge-conditioned filters in convolutional neural networks on graphs In Proceedings of the IEEE conference computer vision and pattern recognition, pp. 3693-
3702, 2017.
Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical bayesian optimization of machine learning algorithms. Advances in neural information processing systems, 25:2951 2959, 2012.
Jasper Snoek, Oren Rippel, Kevin Swersky, Ryan Kiros, Nadathur Satish, Narayanan Sundaram, Mostofa Patwary, Mr Prabhat, and Ryan Adams. Scalable bayesian optimization using deep neural networks. In International conference on machine learning, pp. 2171-21881 2015.
Jost Tobias Springenberg, Aaron Klein, Stefan Falkner, and Frank Hutter. Bayesian optimization with robust bayesian neural networks. Advances in neural in information processing systems, 20:4134-4142, 2016.
Shengyang Sun, Guodong Zhang, Jiaxin Shi, and Roger Grosse. Functional variational bayesian neural networks. arXiv preprint XrX3::1903.057 2019.
Kevin Swersky, Jasper Snoek, and Ryan Prescott Adams, Freeze-thaw bayesian optimization. arXiv preprint 1rXiv:1406.3800 2014.
Wenjic Tang and Graeme Henkelman. Charge redistribution in corc-shell nanoparticles to promote oxygen reduction. The Journal of chemical physics, 130(19):194504, 2009.
The GPyOpt authors. GPyOpt: A bayesian optimization framework in python. http://github. com/
Shef 2016.
17

---

Published in Transactions on Machine Learning Research (09/2022)
Sebastian Thrun. Lifelong learning algorithms. In Learning to learn, pp. 181-209. Springer, 1998.
Austin Tripp, Erik Daxberger, and José Miguel Hernández-Lobato Sample-effcient optimization in the latent space of deep generative models via weighted retraining. Advances in Newral Information Processing Systems, 33, 2020.
Ryan Turner, David Eriksson, Michael J. McCourt, Juha Kiili, Eero Laaksonen, Zhen Xu, and Isabelle Guyon. Bayosian optimization is superior to random scarch for machine learning hyperparamoter tuning: Analysis of the black-bol box optimization challenge 2020. NewrIPS, 2020.
Tsuyoshi Ueno, Trevor David Rhone, Zhufeng Hou, Teruyasu Mizoguchi, and Koji Tsuda. Combo: efficient bayesian optimization library for materials science. Materials discovery, 4:18-20 21,-21, 2016.
Mark Van der Wilk, Carl Edward Rasmussen, and James Hensman. Convolutional gaussian processes. arXiv preprint or Xiv:1 709.0000000 01894, 2017.
Tan Walker and Ben Glocker. Graph convolutional gaussian processes. In International Conference on Machine Learning, pp. 6495 6504. PMLR, 2019.
Ke Wang, Geoff Pleiss, Jacob Gardner, Stephen Tyree, Kilian Q Weinberger, and Andrew Gordon Wilson. Exact gaussian processes on million data points. Advances in Newral Information Processing Systems, 32, 2019.
Linnan Wang, Rodrigo Fonseca, and Yuandong Tian, Learming search space partition for black- bbx optimiza- tion using monte carlo tree search. Advances in Neural Information Processing Systems, 33:19511- 2020.
Zi Wang, Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka. Batched large-scale bayesian optimization in high-dimensional spaces. In International Conference on Artificial Intelligence and Statistics, pp. 745- 754. PMLR, 2018.
Max Welling and Yee W Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th international conference on machine learning (ICML-11), pp. 681-688, 2011.
Yeming Wen, Dustin Tran, and Jimmy Ba. Batchensemble: an alternative approach to efficient ensemble and lifelong learning. arXiv preprint ar arXie:2002.06715 2020.
Florian Wenzel, Kevin Roth, Bastiaan S Veeling, Jakub Swiatkowski, Linh Tran, Stephan Mandt, Jasper Snock, Tim Salimans, Rodolphe Jenatton, and Scbastian Nowozin. How good the bayes postcrior in deep neural networks really? arXiv preprint arXiv:2002.02405, 2020.
Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning. In Artificial intelligence and statistics, pp. 370-378. PMLR, 2016.
James Wilson, Frank Hutter, and Mare Deisenroth. Maximizing acquisition functions for bayesian optimiza- tion. Advances in Neural Information Processing Systems, 31:9884-988 9895, 2018.
Eli Yablonovitch. Inhibited spontaneous emission solid-state physics and electronics. Physical review letters, 20220) 1987.
Wanqian Yang, Lars Lorch, Moritz A Graule, Himabindu Lakkaraju, and Finale Doshi-Velez. Incorporating interpretable output constraints in bayesian noural networks. arXiv preprint (XXx:2010.10960 2020.
Dongkun Zhang, Lu Lu, Ling Guo, and George Em Karniadakis. Quantifying total uncertainty in physics- informed neural networks for solving forward and inverse stochastic problems. Journal of Compulational Physics, 397:108850, 2019a.
Ruqi Zhang, Chunyuan Li, Jianyi Zhang, Changyou Chen, and Andrew Gordon Wilson. Cyclical stochastic gradient meme for bayesian deep learning, orXiv preprint Xiv:: 039322......2.......0022..222,22 2019b.
18

---

Published in Transactions on Machine Learning Research (09/2022)
Table 1: Summary of dataset dimensionalities. Note that alternate inputs for photonic crystal and organic molecule datasets are binary images and molocule graphs, respectively.
CONTINUOUS INPUT ALTERNATE INPUT
AUXILIARY DIMENSION
DIMENSION
DIMENSION
NANOPARTICLE SCATTERING
6
N/A
201
PHOTONIC CRYSTAL DOS
x 32 1024
500
MOLECULE QUANTUM CHEMISTRY
480
+9999+9x99 171
A Appendix
A.1 Datasets
The dimensionalities of the datasets are summarized in table [ The continuous input dimension for chemical molecules refers to the SOAP descriptor. While the space of chemical molccule graphs in general do not have a well-defined dimensionality size of molecules by only sampling from the QM9 dataset, and can define the dimensionality the sum of the adjacency, node, and edge matrix dimensionalities.
chemical molecules can be arbitrarily large and complex, we limit the
The high dimensionalities of all of these problems make Bayesian neural networks well-suited as surrogate models to enable scaling. Note that the nanoparticle scattering problem can be adjusted to be less or more difficult by either changing the input dimensionality (i.e. the number of nanoparticle layers) or the auxiliary dimension (i.e. the resolution or range of wavelengths that are sampled).
A.1.1 Nanoparticle Scattering
The multilayer nanoparticle consists of a lossless silica core surrounded by alternating spherical layers of lossless TiO₂ and lossless silica. The relative permittivity of silica is ésilica = 2.04. The relative permiltivity of TiO2 dispersive and depends on the wavelength of light:
0.2441
ETIO, 5.913 + 10 x2 0.0803
(6)
where is the wavelength given in units of nm. The entire particle is surrounded by water, which has a relative permittivity of Ewater 1.77.
For given set of thicknesses, we analytically solve for the scattering spectrum, i.e. the scattering cross- section (A) as function of wavelength A, using Mie scattering as described in Qiu et al. (0012). The code for computing was adapted from Peurifoy et al. (2018).
The objective functions for the narrowband and highpass objectives are:
SxEUD a (A) dd
lab(z)) Selsewhere o(A) dv
21=126
(7)
+ Σ
201
iil
146 24
SXCUp o (A) dd
201
hhp Teleewhere (A)) dX
(8)
120
1== 24
where z 2001 is the discretized scattering secsinn o(A) from 350 mm to 750
A.1.2 Photonic Crystal
The photonic crystal (PC) consists of periodic unit cells with periodicity a = 1 au, where each unit cell is depicted as a "two-Lone" image, with the white regions representing silicon with permittivity e = 11.4 and black regions representing vacuum (or air) with permittivity Eo =
The photonic crystal (PC) structure defined by spatially varying permittivity e(x, y) (eo. over a 2D periodic unit cell with spatial coordinates x, E [0, a). To parameterize e, we choose level set of
19

---

Published in Transactions on Machine Learning Research (09/2022)
Fourier sum function o, defined as linear combination of plane waves with frequencies evenly spaced in the reciprocal lattice space up to a maximum cutofT. Intuitively, the upper limit on the frequencies roughly corresponds to a lower limit on the feature size such that the photonic crystal remains within reasonable fabrication constraints. Here we set the cutoff such that there are 25 complex frequencies corresponding to real coefficients c= (e1, C2. .2,, (50).
Explicitly, we have
< v) = R
(ck + ick+25) e enniaaa
(9)
where each exponential term is composed from the 25 different pairs (Na, with nay, By E 11,,,1,,,1,2,, We then choose a level-set offset A to determine the PC structure, where regions with Ø > A are assigned to be silicon and regions where ø I A are vacuum. Thus, the photonic crystal unit cell topology is parameterized by a 51-dimensional vector, [c1, C2, .... C50, A E R More specifically,
dfc v) A E1
<(x,y) = =/c, A|(
(10)
Eo
dlc] y) S
which is discretized to result in a 32 x 32 pixel image v f ] This formulation also has the advantage of enforcing periodic boundary conditions.
For each unit cell, we use the MIT Photonics Bands (MPB) software (Johnson Joannopoulos 2001) to compute the band structure of the photonic crystal, (k), up to the lowest 10 bands, using a 32 x 32 spatial resolution (or equivalently, 32 x 32 k-points over the Brillouin zone - < k < 5). We also extract the group velocities at each k-point and compute the density-of-states (DOS) via an extrapolative technique, adapted from Liu et al. (2018). The DOS is computed at resolution of 20,000 points, and a Gaussian filter of kernel size 100 is used to smooth the DOS spectrum. To normalize the frequency scale across the different unit cells, the froquency is rescaled via w /V/arg + Whorm: where Eang = a² Jo Jo E(x, y) dx dy = 1 is the average permittivity over all pixels. Finally, the DOS spectrum is truncated at Whorm = 1.2 and interpolated using 500 points to give z E R2000
(32)2 Σ.g V.
The objective fonction aims to minimize the DOS in a small frequency range and maximize it elsewhere. We use the following:
>300
300
550 500
(z) = 1 Ei=301 2/
2; +
Zi-351 351
2
(11)
where the is added in the denominator to avoid singular values.
A.1.3 Organic Molecule Quantum Chemistry
The Smooth Overlap of Atomic Positions (SOAP) descriptor (De et al. 2016) uses smoothed atomic densities to describe local environments for each atom in the molecule through fixed-length feature vector, which can then be averaged over all the atoms in the molecule to produce a fixed-length feature vector for the molecule. This descriptor is invariant to translations, rotations, and permutations. We use the SOAP descriptor implemented by DScribe (Himanen et al.) 2020) using the parameters: local cutoff rcut = 5, number of radial basis functions nmax = 3, and maximum degree of spherical harmonics 1max = 3. We use outer averaging, which averages over the power spectrum of different sites.
The graph representation of each molecule is processed by the Spektral package (Grattarola & Alippi 2020). Each graph is represented by node feature matrix X E Rexd,, an adjacency matrix A R*Xs and an edge matrix E E Rexd, where s is the number of atoms in the molecule, e is the number of bonds, and d,, de are the number of features for nodes and edges, respectively.
The properties that we use from the QM9 dataset are listed in Table 2 two categories: (1) the ground state quantities which are calculated from a single DFT calculation of the molecule and include geometric, energetic, and electronic quantities, and (2) the thermodynamic quantities which are typically calculated from molecular dynamics simulation.
We separate these properties into
20

---

Published in Transactions on Machine Learning Research (09/2022)
Table List of properties from the QM9 dataset used as labels PROPERTY
UNIT
DESCRIPTION
GROUND STATE QUANTITIES
A
GHz
ROTATIONAL CONSTANT
B
GHz
ROTATIONAL CONSTANT
C
GIl
ROTATIONAL CONSTANT
DIPOLE MOMENT
ISOTROPIC POLARIZABILITY
Ha
ENERGY OF HOMO
CLUMO
Ha
ENERGY LUMO
COAP
Ha
GAP (CLUMO - CHоМО)
(R*)
ab
ELECTRONIC SPATIAL EXTENT
THERMODYNAMIC
QUANTITIES AT 298.15 K
U
Ha
INTERNAL ENERGY
H
Ha
ENTHALPY
G
Ha
FREE ENERGY
Cy
cal
HEAT CAPACITY
molk
Algorithm Bayesian optimization with auxiliary information
1: Input: Labelled dataset Dirain {(Xn, Ano vn) 55 )
2: for = 5 to 1000 do
3: Train M: + Z on Drrain
4: Form an unlabelled dataset, Kpool
Find XN+1 arg max Nponl a (x; M, Dtrain) Label the data 2N+1 ((x++1), UN+1 h(2N41)
Derain Duain U (N++1, #N41)
end for
The auxiliary information for this task consist of the properties listed in Table 2 that are in the same category as the objective property, as these properties would be caleulated together. The objective function then simply picks out the corresponding feature from the auxiliary information. More precisely, for the ground state objectives, the auxiliary information is
[A, B, C, Ot, éHOMO, eLUMO: Egaps (R²) E R
and the objective functions are
ha (z) = 25
ha- Coop (z)
25 28
0.02
191
0.6
where the quantities for the latter objective normalized that they have the same magnitude.
A.2 Bayesian Optimization and Acquisition Function
Our algorithm for Bayesian optimization using auxiliary information is shown in Algorithm 1 This algorithm reduces to the basic BO algorithm in the case where h is the identity function and Z such that we can ignore mention of z in Algorithm
As mentioned in the main text, the inner optimization loop in line 5 of Algorithm is performed by finding the maximum value of a over a pool of [Xpooll randomly sampled points. We can see in Figure 0 that increasing |pool in the acquisition step tends to improve BÒ performance. Thus, there is likely further room for improvement of the inner optimization loop using more sophisticated algorithms, possibly using the gradient information provided by BNNs Unless otherwise stated, optimize the inner loop of Bayesian
21

---

Published in Transactions on Machine Learning Research (09/2022)
0.362
Bx256
0.160
164512-aux
0.159
0.157
0.156
10
10
165
[xpod
Figure Eiffet m |Xpooll used in the inner optimization loop to maximize the acquisition function on overall BO performance. Ybest is taken from the narrowband objective function using the ensemble architecture. The "aux" in the legend denotes using auxiliary information and the numbers represent the architecture (i.e. 8 layers of 256 units or 16 layers of 512 units).
Narrowband
Highpass
0.162
1.00
B8B-aux
BBB-qux-norestart
1.06
Ensemb le aux-restart
0.156
Ensemble Ensembleeauux t-aux-restart--ydee
1.04
0.154
le-aux-norestart
Ensemb ey-au--norestat-cccee
3.03
0.152
0.150
N
N
Figure 7: Effect of restarting the BNN training from scratch in each BO iteration.
optimization to choose the next data point to label by maximizing El on a pool of [Xpooil = 105 randomly sampled points.
A.3 Continued Training
As mentioned Section 2.3 of the main text, the BNN ideally trained from scratch until convergence each iteration loop, although this comes at great computational cost. An alternative is the warm restart method of continuing the training from the previous ileration which enables the model's training loss to converge in only few epochs. However, shown in Figure 7 we have found that naive continued training can result in poor BO performance. This is likely because (a) training does not converge for the new data point Dnew = (NN+1,9N+11 relative to the rest of the data under a limited computational budget, resulting in the acquisition function possibly labeling similar points in consecutive iterations, and (b) the BNN gets trapped in a local minima the loss landscape that is not ideal for learning future data points. To mitigate this, we use the cosine anncaling learning rate proposed in Loshchilov & Hutter (2016). The large lcarning rale at the start of training allows the model to more easily escape local minima and explore a multimodal posterior (Huang et al. 2017,, while the small learning rate towards the end of the annealing cycle allows the model to converge more easily. Note that the idea of warm restart is similar to "continual learning,"
which is an open and active sub-problem in machine learning research (Trun 1998 Parisi et al. 2019 In particular, we re-train the BNN using 10 epochs.
A.4 Models and Hyperparameters
A.4.1 Additional Surrogate Models
Variational BNNs model a prior and posterior distribution over the neural network weights, but use some approximation the distributions to make the BNN tractable. In particular, Bayes by Backprop
22

---

Published in Transactions on Machine Learning Research (09/2022)
(BBB) (also referred to as the "mean field" approximation), which approximates the posterior over the neural network weights with independent normal distributions (Blundell et al. 2015). We also compare Multiplicative Normalizing Flows (MNF),, which uses normalizing flows on top of each layer output for more expressive posterior distributions Louizos & Welling 2017).
BOHAMIANN proposed to use BNNs in BO by using stochastic gradient Hamiltonian Monte Carlo (SGHMC) to approximately sample the BNN posterior, combined with scale adaptation adapt it for an iterative setting (Springenberg et al. 2016].
NeuralLinear trains conventional neural network on the data, but then replaces the last layer with Bayesian linear regression such that the neural network serves as an adaptive basis for the linear regression Snoek et al. 2015).
TuRBO (trust region Bayesian Optimization) a method that maintains M trust regions and performs Bayesian optimization within each trust region, maintaining M local surrogate models, to scale BO to high- dimensional problems that require thousands of observations (Eriksson et al. 2019). We use M and
M 5, labelled "TuRBO-1" and "TuRBO-5", respectively.
TPE (Tree Parzen Estimator) is a method that instead of modeling p(y|z), models p(zly) and p(v) for the surrogate model and fits into the BO framework (Bergstra et al. 2013). The tree-structure of the surrogate model allows it to define leaf variables only when node variables take particular values, which makes well- suited for hyper-parameter search (e.g. the learning rate momentum is only defined for momentum-basee gradient descent methods).
LIPO is a paramcter-free algorithm that assumes the underlying function is Lipschitz function and osti- mates the bounds of the function (Malherbe & Vayatis 2017). We use the implementation provided by the
dlib library (King) 2009).
DIRECT-L (DIviding RECTanglcs-Local) systematically divides the search domain into smaller and smaller hyperreclangles to efficiently search the space (Gablonsky & Kelley 2001). We use the implementation provided by the NLopt library (Johnson 2010).
CMA-ES (covariance matrix adaptation evolution strategy) is an evolutionary algorithm that samples new data based on mullivariate normal distribution and refines the parameters of this distribution until reaching convergence. We us the implementation provided by the pyema library (Hansen et al. 2019
A.4.2 Implementation Details
Unless otherwise stated, we set NMC = 30. All BNNs other than the infinitely-wide networks are implemented in TensorFlow v1. Models are trained using the Adam optimizer using the cosine anncaling learning rate with a base learning rate of 10-3 (Loshchilov & Hutter 2016). All hidden layers use ReLU as the activation function, and no activation function is applied to the output layer.
Infinite-wddth neural networks are implemented using the Neural Tangents library (Novak et al. 2020 We
use two different types of infinite networks: (1) "GP-" refers to a closed form expression for Gaussian process inference using the infinite- width neural network as kernel, and (2) "INF-" refers to an infinite ensemble of infinite-width networks that have been "trained" with continuous gradient descent for an infinite time. We compare NNGP and NTK kernels well as the parameterization of the layers. By default, we use the NTK parameterization, but also use the standard parameterization, denoted by "-STD".
We implement BO using GPs with a Matérn kernel using the GPyOpt library The GPyOpt authors 2016 The library optimizes over the acquisition function in the inner loop using the T-BFGS algorithm.
LIPO (Malherbe & Vayatis 2017) is implemented in the dlib library King 2009 D. DIRECT-L (Gablonsky &
Kelley 2001) implemented in the NLopt library (Johnson 2010). CMA- ES implemented in the pycma
library (Hansen et al. 2019).
23

---

Published in Transactions on Machine Learning Research (09/2022)
Branin
Hartmann-6
-0.22
50
-1.22
GP
Ensemble
-1.75
BBB
-2.00
05
100
220
302
400
300
200
200
202
400
300
Figure 8: BO results for the Branin and Hartmann-6 functions.
Narrowband
Highpass
0.1600
Ensemble-aux
Ensemble
Ybust
BBB-aux
0.1550
1.06
Neural Linear
0.1500
1.05
BOHAMIANN
0.1475
1.04
800
1000
200
400
0.36
1.08
Ensemble-aax
Yoest
DIRECT-L
1.04
CMA
1.00
200
400
800
800
1000
0
200
400
600
800
1.000
Figure 9: Additional optimization result curves for the nanoparticle scattering task. (Top) Various BNNs. Note that results using auxiliary information are denoted by solid line, while those that do not are denoted by dashed line. Also note that the y-axis zoomed in differentiate the curves. (Bottom) Various non-BO algorithms. ENSEMBLE-AUX is replicated here for case of comparison.
A.5 Additional Results
A.5.1 Test Functions
We test BO several common synthetic functions used for optimization, namely the Branin and 6- dimensional Hartmann functions. We use BNNs with hidden layers and 256 units in each hidden layer, where each hidden layer is followed by ReLU activation function. Plots of the best value Shest at each BO iteration are shown in Figure S As expected, GPs perform the best. Ensembles and BBB also perform com- petitively and much better than random sampling, showing that deop BO is viable even for low-dimensional black-box functions.
A.5.2 Nanoparticle Scattering
Detailed BO results for the nanoparticle scattering problem are shown in Table B
All the BNNs used for the nanoparticle scattering problem use an architecture consisting hidden layers with 256 units each, with the exception of BOHAMIANN where we used the original architccture consisting of hidden layers with 50 units each. The infinite-width neural networks for the nanoparticle task consist of 8 hidden layers of infinite width, each of which followed by ReLU activation functions.
24

---

Published in Transactions on Machine Learning Research (09/2022)
Table 3: BO results for the nanoparticle scattering problem. denotes that Ybest is measured at N 100 due to computational constraints
MODEL
NARROWBAND
HIGHPASS
Yaeer AT
250
yEeS7 AT N 1000
YBEST AT N 250
UBEST AT N
100
MEAN
SE
MEAN
SE
MEAN
SE
MEAN
SE
GP
0.1606
0.0005
0.1621
0.0001
1.0839
0.0017
1.0851
0.0008
GP-AUX
*0.1541
0.0019
*1.0110
0.0234
ENSEMBLE
0.1558
0.0011
0.1607
0.0003
1.0729
0.0025
1.077
0.0021
ENSEMBLE-AUX
0.1578
0.0014
0.1593
0.0013
1.0783
0.0003
1.0822
0.001
BBB
0.1596
0.0006
0.1596
0.0006
1.0753
0.0005
1.0753
0.0005
BBB-AUX
0.1601
0.001
0.1601
0.001
1.076
0.0028
1.076
0.0028
BBB-ANNEAL
0.1598
0.001
0.1611
0.0001
1.0813
0.0003
1.0821
0.0005
BBB-AUX- LUX-ANNEAL
0.1613
0.0003
0.1619
0
1.0826
0.0008
1.0834
0.0005
MNF
0.15
0.0005
0.1547
0.0004
1.027
0.005
1.0312
0.0036
MNF-AUX
0.1549
0.0014
0.1569
0.0006
0.9957
0.0168
1.028
0.0157
NEURAL LINEAR
0.1543
0.002
0.1579
0.0015
1.0798
0.0007
1.0836
0.0007
BOHAMIANN
0.1616
0.0001
1.0717
0.0031
INF-NNGP
0.1541
0.0011
0.157
0.0009
1.055
0.0036
1.0653
0.0022
INF-NTK
0.1536
0.0008
0.1571
0.001
1.041
0.004
1.0612
0.0011
INF-NNGP-STD
0.1551
0.0006
0.1598
0.0006
1.0615
0.0043
1.069
0,0018
IN--NTK-STD
0.1564
0.0006
0.1607
0.0001
1.0607
0.0039
1.0761
0.0014
GP-NNGP
0.1582
0.0007
0.1609
0.0001
1.0621
0.0027
1.0694
0.0019
GP-NTK
0.1573
0.001
0.1611
0.0001
1.0667
0.0032
1.0732
0.0012
GP-NNGP-STD
0.1562
0.0008
0.1595
0.001
1.0615
0.0058
1.0718
0.0024
GP-NTK-STD
0.1592
0.0011
0.1608
0.0002
1.0641
0.0033
1.0704
0.0017
TuRBO-1
0.1572
0.0017
0.1619
0.0011
1.0831
0.022
1.0871
0.0005
TuRBO-1
0.1605
0.0011
0.1619
0.0001
1.0867
0.0016
1.0890
0.0003
TPE
0.1561
0.0007
0.1615
0.0001
1.0517
0.0035
1.0794
0.0010
RANDOM
0.1527
0.0008
0.1555
0.0006
1.0053
0.0063
1.0362
0.0047
LIPO
0.1604
0.0016
0.1619
0.0006
1.0792
0.0066
1.087
0.0034
DIRECT-L
0.1544
0
0.156
0
1.0777
1.0801
0
CMA
0.1424
0.0046
0.143
0.0048
1.059
0.0117
1.076
0.0127
25

---

Published in Transactions on Machine Learning Research (09/2022)
1.050
1.025
layers4
layers4-aux
Yiett
layersa
layersB-aux
1.0750
layers16
1.0725
layers16-eux
1,0700
64
Layer width
Figure 10: Comparison of Ybest at N 1000 for the nanoparticle narrowband objective function for variety of neural network sizes. All results are ensembles, and "aux" denotes using auxiliary information.
= 30.09,, 39.76, 39.76,30.1169.60
[67:44.36.04669902222 69.06]
3.0-
4.5
4.0
3.5
Scattering
Www
25
500
700
300 Weveeength [nm]
Wovelength [nm]
Figure 11: Examples of optimized nanoparticles and their scattering spectrum using the "ENSEMBLE-AUX" architecture for the (a) narrowband and (c) highpass objectives. Orange shaded regions mark the range over which we wish to maximize the scattering.
We also experiment with KL annealing BBB, a proposed method to improve the performance of variational methods for BNNs in which the weight of the KL term in the loss function is slowly increased throughout exponentially anneal the KT term with weight function of epoch i when training from scratch; during the continued training, the weight is held constant at ORt. = 10-3
training Wenzel et al. (020)). For these experiments, (t) 10
KL annealing in the BBB architecture significantly improves performance for the narrowband objective, although results are mixed for the highpass objective. Additionally, KI. annealing has the downside of introducing more parameters that must be carefully tuned for optimal performance. MNF performs poorly, especially on the highpass objective where is comparable random sampling, and we have found that MNF is quite sensitive to the choice of hyperparameters for uncertainty estimates even simple regression problems.
The different variants infinite-width neural networks do not perform as well as the BNNs on both objective functions, despite the hyper-parameter search.
LIPO seems perform as well as GPs on both objective functions, which is impressive given the computa- tional speed of the LIPO algorithm. Interestingly DIRECT-L does not perform as well as LIPO or GPs on the narrowband objective, and actually performs comparably random sampling the highpass objective. Additionally, CMA performs poorly both objectives, likely due to the highly multimodal nature of the objective function landscape.
We also look al the effect of model size in terms of number of layers and units in Figure 10 for ensembles. While including auxiliary information clearly improves performance across all architectures, there is not clear trend of performance with respect to the model size. Thus, the performance of BO seems to be somewhat robust to the exact architecture as long as the model is large enough accurately and efficiently train the data.
26

---

Published in Transactions on Machine Learning Research (09/2022)
Table 4: Various architectures for BNNs and BCNNS used in the PC problem. Numbers represent the number of channels and units for the convolutional and fully-connected layers, respoctively. All convolutional layers use 3 x 3-sized filters with stride 11, and periodic boundaries. "MP" denotes max-pooling layers of size 2 x 2 with stride (2, 2), and "AP" denotes wverage-pooling layers of size 2 x 2 with stride (1,1). "CONV" denotes BCNNs whereas "FC" denotes BNNs (containing only fully-connected layers) that act on the level- set parameterization x rather than on the image v. "TI" denotes translation invariant architectures, whereas "TD" denotes translation dependent architectures (i.e. not translation invariant).
CONVOLUTIONAL
FULLY-CONNECTED
ARCHITECTURE
LAYERS
LAYERS
CoNv-TI
16-MP-30 12--32-M-- MP-6PP-4--M- 28-MMP226- 256
256 256 256-256
CONv-TD
A-PP8-MP-16-AM AP-332 APP----P--------- MP-32-AP
256 256 256-226622222 256
FC
N/A
256- 25--2556-2222225
(a)
(b) 25
(c)
(d)
X
30-
10
10
0.0
0.2 0.0 1.0 1.2
frequency
Figure 12: Examples of optimized photonic crystal unit cells over multiple trials for (a) PC-A distribution and (c) PC-B distribution. (b,d) Examples of the optimized DOS. Note that the DOS has been minimized to nearly zero thin frequency range. Orange shaded regions mark the frequency range which we wish minimize the DOS, All results were optimized by the "ENSEMBLE- AUX" architecture.
Examples of the optimized structures by the "ENSEMBLE-AUX" architecture are shown in Figure 11 We can see that the scattering spectra peak in the shaded region of interest, as desired by the respective objective functions.
A.5.3 Photonic Crystal
The BNN and BCNN architectures that we use for the PC task are listed in Table The size of the "FC" architectures are chosen to have similar number of parameters as their convolutional counterparts. Unless otherwise stated, all results in the main text and here use the "CoNv-TI" and "FC" architectures for BCNNs and BNNs, respectively.
The infinito-width convolutional neural networks (which act convolutional kernels for GPs) in the PC task consist of 5 convolutional layers followed by 4 fully-connected layers of infinite width. Because the pooling layers in the Neural Tangents library are currently too slow for use in application, we increased the size of the filters to 5 x 5 to increase the receptive field of each filter.
Detailed BO results for the PC problem are shown in Table 5 For algorithms that optimize over the level set parameterization R51, we see that GPs perform consistently well, although BNNs using auxiliary information (e.g. Ensemble-Aux) can outperform GPs. DIRECT-L and CMA perform extremely well on the PC-A distribution but performs worse than GP on the PC-B distribution.
Adding convolutional layers and auxiliary information improves performance such that BCNNs significantly outperform GPs. Interestingly, the infinite-width networks perform extremely poorly, although this may be due to a lack pooling layers in their architecturs ure which limits the receptive field of the convolutions.
Examples of the optimized structures by the "ENSEMBLE-AUX" architecture are shown Figure 12 The photonic crystal unit cells generally converged the same shape: square lattice of silicon posts with periodicity V2a.
27

---

Published in Transactions on Machine Learning Research (09/2022)
Table 5: Select BO results for the PC problem. denotes that Wbest is measured at N 130 due to computational constraints. denotes that Ybest is measured at N = 750 due to computational constraints. MODEL
PC-A
PC-B
Yues AT
250
Your AT
1000
Yueer AT
250
Doesr AT N 100
MEAN
SE
MEAN
SE
MEAN
SE
MEAN
GP
548
450
2109
448
781
394
3502
49
GP-AUX
16
99
ENSEMBLE
30
841
448
216
145
1318
465
ENSEMBLE-AUX
305
217
1310
509
2909
408
3633
130
CONVENSEMBLE
1140
2375
390
263
2070
505
CONVENSEMBLE-AUX
2623
558
3468
120
3752
106
4002
92
BBB
75
350
207
704
502
780
485
BBB-AUX
39
7
413
313
554
371
1605
541
CONVBBB
712
416
1486
490
928
600
930
599
CONYBBB-AN AUX
2109
583
3124
43
1761
724
1928
711
NEURALLINEAR
1009
549
1235
481
685
488
2853
291
CONVNEURALLLNEAR
1160
540
2524
479
1643
596
2722
647
CONv-IN--NNGP
29
322
181
21
157
42
CONV-INP-NTI
49
425
322
28
907
711
CONv-GP-NNGP
15
2
221
118
37
830
533
CONY-GP-NIK
20
12
45
CONV-INF-NNGP-STT
17
3
732
432
66
15
889
442
CONVIINE-NTK-STD
8
12
CONV-GP-NNGPA
20
7
'101
59
100
55
1124
49
CON--GP-NTK-STTD
13
1332
77
7
1686
575
RANDOM
141
61
402
184
471
398
485
395
TuRBO-1
1150
638
4451
20
3865
289
4476
16
TERBO-5
3738
4456
37
4128
49
4466
26
TPE
1001
648
3901
140
3045
571
4119
156
LIPO
940
1073
1280
1073
1837
1792
2266
1626
DIRECT-L
20
0
4351
8
2525
38
CMA
9
4078
126
1777
969
28

---

Published in Transactions on Machine Learning Research (09/2022)
(8)
1090
100
MSE
MAE
NLL
107
Calibration
107
107
10-3-
400
620
800
600
230
N
N
N
N
[6]
Nm50
N=200
N=500
.... Ideal
GP-aux
0.6
Ensambla-aux
ConvErsemble
0.2
ConvEnsemble-aux
NTK
1.0
Figure 13: (a) Various metrics tracked during BO of the PC-A dataset distribution validation dataset of 1000 datapoints. (b) Uncertainty calibration curves measured at various points during BO. Note that the calibration curve for GP-AUX is only shown for N = 50, as it becomes computationally intractable for larger N.
Validation Metrics
To explore more deeply why certain surrogate models perform well while others do not, track various metrics of the model during BO on validation dataset with 1000 randomly sampled data points. In particular, we look at the mean squared error (MSE), the mean absolute error (MAE), the negative log- likelihood (NLL), and the calibration error on the PC-A data distribution. Results are shown in Figure 13).).
The calibration error is quantitative measure of the uncertainty of the model, which is important for the performance of BO as the acquisition function the uncertainty to balance exploration and exploitation. Intuitively, expect that 50% confidence interval contains the correct answer 50% of the time. In particular, we use the forecast calibration as proposed by Kuleshov et al. (208)), which is an extension of the calibration error proposed by Guo et al. (0177) to regression tasks:
cal(F1, V1, Fr,
Σ (p) -,)²
where F, is the CDF of the predictive distribution, Py is the confidence level, and py is the empirical frequency. We choose to measure the error along the confidence levels P5 = (j 1)/10 for 1,,, ....,,,,. 11. The CDF
F2 (15) an be analytically calculated for models that have an analytical predictive distribution. For models that do not have an analytical predictive distribution, we use the empirical CDE:
1 Σ 4
F(y)
==
where is the indicator function. We also plot the calibration, (D);,;))M in Figure 13 (). Perfectly calibrated predictions correspond to a straight line.
Figure 13 shows that the infinite neural network kernel (NTK) has the highest prediction error, which is likely a contributing factor to its poor BO performance. Interestingly, vanilla GPs have the lowest MSE, so the prediction error is the only indicator for BO performance. Looking at the calibration, the infinite neural network kernel has the highest calibration error, and from Figure (b) that tends to be
29

---

Published in Transactions on Machine Learning Research (09/2022)
Egap
Egop
0.6
140
Ybert
GrashEnsemble-auu BBB FC
NeuralLinear
100 0.2
0.4
VAE-GP
400
200
N
N
N
Figure 14: Additional BO results for several different objective functions on the chemistry dataset. GP and GRAPHENSEMBLR-ACX curves are replicated from the main text for convenience.
Table 6: BO results for the four different quantum chemistry objective functions. denotes that Phest is measured at N 100 due to computational constraints.
YBEST AT
500
MODEL
CoAP
Cuar
C
Ct GGAr
MEAN
SD
MEAN
SD
MEAN
SD
MEAN
SD
GP
0.41
0.04
-0.10
0.02
101.08
1.05
0.29
0.07
GRAPHGP
*0662
0.00
0.10
0.02
1131.9
14.59
*.224
0.03
ENSEMBLE
0.62
0.00
-0.08
0.00
86.56
0.28
0.00
ENSEMBLE-AUX
0.62
0.00
-0.10
0.02
83.86
4.45
0.13
0.05
GRAPHENSEMBLE
0.62
0.00
-0.10
0.00
143.53
0.00
0.49
0.00
GRAPHENSEMBLE-UX
0.62
0.00
-0.10
0.00
143.53
0.00
0.49
0.00
GRAPHBBB
0.38
0.01
-0.11
0.01
94.46
1.16
0.25
0.01
GRAPHBBB FC
0.62
0.00
-0.10
0.00
135.64
13.67
0.39
0.14
GRAPHINEURALLINEA
0.62
0.00
0..0
0.01
143.53
0.00
0.46
0.09
VAE-GP
0.62
0.06
123.33
13.02
0.61
0.34
VAE-GP-2
110.84
16.68
0.56
0.35
AAE-GP-LATENTIU
154.66
35.96
0.40
0.10
133.66
13.25
0.42
0.13
PAE-GP-LATENT22
114.83
14.64
0.38
RANDOM
0.38
0.02
-0.11
0.03
105.19
7.87
0.29
0.07
overconfident in its predictions. GPs have a higher calibration error than the ensemble neural network methods, and tends be significantly underconfident in prodictions. GP-AUX has higher validation loss, calibration error, and NLL than most, if not all, of the other methods, which explain its poor performance.
The ensemble NN methods tend to be reasonably well-calibrated. Within the ensemble NNs, the "-aux" methods have lower MSE and calibration error than their respective counterparts, and ConvEnscmble-aux has the lowest NLL calibration error out of all the methods, although interestingly Ensemble-aux seems to have the lowest MSE and MAE out the ensemble NNs.
These results together show that calibration of Bayesian models extremely important for use surrogate models in BO.
A.5.4 Organic Molecule Quantum Chemistry
The Bayesian graph neural networks (BGNNs) used for the chemical property optimization task consist of 4 edge-conditioned graph convolutional layers with 32 channels each, followed by a global average pooling oper- ation, followed by 4 fully-connected hidden layers of 64 units each. The edge-conditioned graph convolutional layers Simonovsky & Komodakis (2017) are implemented by Spektral Grattarola & Alippi (020).
More detailed results for the quantum chemistry dataset are shown in Table and Figure 14 The archi-
tecture with the Bayes by Backprop variational approximation applied to every layer including the graph
30

---

Published in Transactions on Machine Learning Research (09/2022)
Egup
0.6
160
0.6
VAE-GP
Youst
140
0.5
VAE-GP-atten1128
0.4
120
VAE-GP-atent22--b6600000
0.3
la-tePlatent322
100
0.2
0.1
80
100
400
500
100
500
N
N
Figure 15: Additional BO results for VAE-GP using different pre-trained VAEs.
convolutional layers (BBBB"), performs extremely poorly, even worse than random sampling in some cases. However, only making the fully-connected layers Bayesian (*BBB-FC") performs surprisingly well.
Ensembles trained with auxiliary information ("ENSEMBLE-AUX") and neural linear ("NEURALLINEAR") perform the best on all objective functions. Adding auxiliary information to ensembles helps for the a objective function, and neither helps nor hurts for the other objective functions. Additionally, BNNs perform at least as well or significantly better than GPs in all cases. GPs perform comparably or worse than random sampling in several cases.
As noted in the main text, the performance of VAE-GP depends on the quality of the pro-trained VAE, shown in Figure 15 The VAE-GP benchmark uses the same pre-trained VAE, and "VAE-GP-2" refers to the same method using a different random seed for the VAE. Even with the exact same method, VAE-GP-2 performs significantly worse on both objective functions. We also increase the latent space dimensionality benchmark, which performs even worse on the a Egap
from 52 to 128 in the "VAE-GP-LATENT122 benchmark although it performs significantly better on the a benchmark. We also adjust the learning rate momentum to 8 0.001 in "VAE-GP-LATES" "VAE-GP-LATENT322 There is no clear trend with the different ayper-parameters, which may point to the random seed of the VAE pre-training being greater factor in BO performance than the
and the latent space dimensionality to 32 in
Validation Metrics
As in Appendix A.5.3 we track the MSE, NLL, and calibration error during optimization on the chemistry task. Results are shown in Figure 16 The various metrics correlate with the respective methods' peformances during BO. For example, VAE-GP has an extremely high MSE and calibration error on the a objective, where it performs poorly, but has an MSE and calibration error more comparable with that of other methods well as extremely low NLL on the metrics for GRAPHGP are very high on the a Cgap objective, where it performs poorly. GraphEnsemble tends to be among the better methods in terms of these metrics, which translates into good BO performance.
Égap objective, where it performs extremely well. Likewise, the
A.5.5 Additional Discussion
BBB performs reasonably well and is competitive with or even better than ensembles on some tasks, but it requires significant. hyperparameter tuning. The tendency of variational methods such as BBB to underesti- mate uncertainty is likely detrimental to their performance in BO. Additionally, [Sun al. (2019) shows that BBB has trouble scaling to larger network sizes, which may make them unsuitable for more complex lasks such those in our work. BOHAMIANN performs very well on the nanoparticle narrowband objective and comparable to other BNNs without auxiliary information on the nanoparticle highpass objective. This is likely due to its effectiveness in exploring a multi-modal posterior. However, the need for SGHMC to sample the posterior makes this method computationally expensive, and so we were only able to run it for a limited number of iterations using a small neural network architecture.
Infinitely wide neural networks are another interesting research direction, as the ability to derive infinitely wide versions of various neural network architectures such as convolutions, and more recently graph con-
31

---

Published in Transactions on Machine Learning Research (09/2022)
MSE
Calibration error
NLL
Calibration Curve
100
10"
10a
10
10
10a
100
0.2
10a
0.0
10
1.0
3x 10
1011
20-1
10a
10-4
107
0.2
0.0
100
200
300
400
100
200
300
400
100
200
300
400
0.00
0.25
0.50
0.75
1.00
p
Ideal
GraphGP
Ensemble
GraphEnsemble
VAE-GP
Figure 16: (a) Various metrics tracked during BO of the PC-A dataset distribution on a validation dataset of 1000 datapoints. (b) Uncertainty calibration curves measured at various points during BO.
volutional layers (Hu et al. 2020) could potentially bring the power of GPs and BO to complex problems in low-data regimes. However, find they perform relatively poorly in BO, are quite sensitive to hyper- parameters (e.g. kernel and parameterization)) and current implementations of certain operations such pooling are too slow for practical use in an iterative setting, In particular, BO using an infinite ensemble of infinitewwddth networks performs poorly compared to normal ensembles, suggesting that the infinitewiddth formulations do not fully capture the dynamics of their finite-width counterparts.
Non- -on-Bayesian global optimization methods such as LIPO and DIRECT-L are quite powerful in spite of their small computational overhead and can outperform BO some simpler tasks. However, they are not consistent BO, performing more comparably to random sampling other tasks. CMA-ES performs poorly all the lasks here. Also, like GPs, these non-Bayesian algorithiss assume continuous input space and cannot be effectively applied to structured, high-dimensional problems.
A.6 Compute
All experiments were carried out on systems with NVIDIA Volta V100 GPUs and Intel Xeon Gold 6248 CPUs. All training and inference using neural botwork-based models, graph kernels, and infinito-width neural network approximations are carried out on the GPUs. All other models are carried out on the CPUs.
32


Reviews:
[{"review_id": "xvSs58OrOl", "reviewer": "Reviewer_bRGJ", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\n### Main claims are supported by experiments\n\nThe first claim of the paper is that auxiliary information helps to speed-up the optimization process. On all three benchmarks the proposed BO strategy with neural networks ensembles achieves the same or faster convergence when it is augmented with auxiliary information. The only exception is for the Organic Molecule Quantum Chemistry with a HOMO-LUMO energy gap, where the GraphEnsemble seems to perform slightly worse when augmented with auxiliary performance.\n\nSecond, the paper claims that BO with convolutional or graph neural networks achieves a better performance on complex input spaces such as images or molecules. This is shown for the Photonic Crystal Topology benchmark. However, as mentioned below, I miss another baseline to support this claim for the Organic Molecule Quantum Chemistry that uses standard feed forward neural networks instead of graph neural networks.\n\nLastly, the paper claims that neural networks outperform Gaussian processes on these problems. This is also supported by empirical evidence at least for two benchmarks. However, on the Multilayer Nanoparticle GPs without auxiliary data seem to work slightly better than the proposed BO strategy with neural networks.\n\n\n### Relevance to the community\n\nI think the paper could be interesting for the BO community because it a) provides some interesting and challenging benchmarks that encourage further research into this direction and b) establishes some sensible baselines for neural network based Bayesian optimization.\n\n## Weaknesses\n\n\n### Empirical evaluation can be improved\n\nApart from adding the missing baseline for the Organic Molecule Quantum Chemistry, I think the empirical evaluation of the paper can be further improved by a comparison to other Bayesian optimization method that use other probabilistic models, such as random forests or kernel density estimators, which have shown competitive performance on high dimensional structured spaces before. Also, I don't understand why the VAE based GP method has only be considered for the Organic Molecule Quantum Chemistry - where it achieved competitive performance -  and not for the other two tasks.\n\n\n### Further clarifications\n\nSome parts of the paper require further clarifications:\n\n- Equation 3 seems wrong. Normally you first compute the empirical mean \\mu = \\sum_i \\mu^(i) and variance \\sigma = \\sum_i (\\mu^(i) - \\mu)^2) of your ensemble and use \\mu and \\sigma to compute Equation 2.\n\n- Section 2.2:  How exactly do you use cosine annealing here? For how many steps? I also don't understand how cosine annealing here helps to learn the new data point? Do you also use consine annealing if you train from scratch?\n\n- How do you optimize the acquisition function for neural network based Bayesian optimization methods?\n\n- Figure 5 (d): Do you train both GraphGP and GraphEnsemble on GPU or CPU?\n\n", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper discusses three different scientific optimization problems with high-dimensional search spaces form physics and chemistry applications, that exhibit known structures or auxiliary information and how to solve them with Bayesian optimization. Due to the high dimensionality of the data, the paper proposes deep neural networks as probabilistic model for Bayesian optimization, which scale better with the number of data points and dimensions than more commonly used Gaussian processes.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Section 4.3: It is unclear how much the graph neural network actually helps for this task compared to using a ensemble based on MLPs that operates on the SOAP descriptor directly. It would be great if the authors could add this as additional baseline\n\nPlease also add a comparison to Bayesian optimization with other probabilistic models such as random forests or kernel density estimators. Furthermore, it would be nice if the VAE-GP baseline could be also applied for the other two tasks.\n\n\n## Minor comments:\n\n- Typo Section 1.1 drop the word  'Ref' in the second and third sentence.\n\n- Can you increase the figure size to make the plots more visible? You waste a lot of white space for the margins. Alternatively, you can move some plots from the appendix to the main paper, to fill these white spaces.\n\n- Figure 5: Does 'GraphEnsemble' in the subplot (c) corresponds to BGNN? If so, could you use the same label?\n\n", "broader_impact_concerns": "I do not see any ethical concerns with this work", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "pc0xKjvVn0", "reviewer": "Reviewer_Cgt7", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "As a full disclosure, I’ve reviewed this work previously, and do appreciate many of the changes made. The figures are of much higher quality now, and the study on VAE – GPs is much clearer now. \n\n**Strengths:**\n\nThe description of the tasks and their relevance to a scientific audience is a selling point. I also appreciate the charts of their descriptions.\n\nExperimental design is generally quite broad and incorporates good baselines:\n-\tExperiments using CMA-ES (and other non BO approaches) (Fig. 9), even if LIPO tends to perform quite well. \n\n-\tExperiments on tabular data, even if standard BO performs quite well (Fig. 8). \n\n-\tKnockout / ablation experiments on several of the test problems. (Figs. 6, 10).\n\n\n**Weaknesses:**\n\nThe scientific and experimental claims are still a bit confusing to me (see Model fits, TuRBO comparison, and VAE-GP proposals for fixing).\n\n-\tThe abstract + intro could be re-worked a bit more to make the claims / experiments in the main text stand out. For example, VAE-GPs are mentioned, but they don’t appear anywhere in the main text.\n\n-\tThe takeaways from the paper proposals could be editorialized a bit more – what exactly did you demonstrate (that inductive biases and NN ensembles help in BayesOpt, e.g. well designed surrogate models)?\n\nCode isn’t provided (but I do hope that this would be remedied on acceptance).\n", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies Bayesian optimization using Bayesian neural networks and ensembles of neural networks as the surrogate models for optimizing three complicated scientific design tasks (multilayer nanoparticles, crystal topology, and quantum chemistry). The authors also demonstrate the success of using auxiliary information, e.g. information outside the existing objective, as well as the strong performance of simple ensembles of neural networks compared to simple Gaussian process models.\n\nThe authors use both cyclical learning rates for continued retraining of their networks (an advance over BOHIMIANN, which instead uses scale adaptation, a weaker type of retraining) as well as incorporate auxiliary information (like image outputs from the simulators) into their NNs (much like Astudillo & Frazier, ’19 and Balandat et al, ’20 did with GPs).\nOverall, the authors demonstrate the importance of high performing surrogate models, presumably in terms of accuracy, for having highly performant Bayesian optimization loops for scientific design problems.\n\n**Would some individuals in TMLR's audience be interested in the findings of this paper?**\n\nYes, Bayesian optimization for scientific design and discovery is an extremely interesting set of applications for machine learning. The findings in this paper should certainly interest people who read TMLR.\n\n**Are the claims made in the submission supported by accurate, convincing and clear evidence?**\n\nPossibly, I still think that the interplay between Bayesian neural nets, ensembles, and Gaussian processes could be better studied and understood in this work. Please see the following sections.\n\n", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "_Model fits (medium / high):_\n\nI think that showing model fits, describing the performance on a heldout set on each task (as much as possible), for a couple of surrogate models would vastly clean up the story of this paper. Judging from the experiments, it’s sort of easy to tell when various models are likely performing well in prediction tasks (e.g. the Graph nets on Fig. 5 or the GP on highpass in Fig. 2), but it would be really nice to tell this type of performance for sure. \n\nAfter all, the experimental results seem to demonstrate to me that one should use the best performing surrogate model for a given problem (as long as it has reasonable uncertainty – show the ensembles have reasonable uncertainty please) and be somewhat agnostic about using GPs or any other model class. However, this type of message would be a lot more clear with these types of model fits. \nThis would also be helpful in diagnosing the underperformance of the GP-aux model as well (Fig. 4). A similar argument could be made to study the weak performance of the infinite NNs throughout (less flexibility produces worse performance).\n\nIt also would help parse out when / how auxiliary information helps – does it help in the objective itself, or just in uncertainty. How exactly?\n\nCode (medium / high): Providing the code would be helpful. I’d personally like to see the implementation of the graphGP, as it seems very strange that it both uses the GPU and is somehow 15x slower than a GP on the CPU. \n\n_Comparison with TurBO (medium):_\n\nI still think that a strong baseline, especially in these sorts of large n problems would be helpful. \n\n-\tUseful reference is Eriksson et al, ’19 and the implementations in BoTorch (https://botorch.org/tutorials/turbo_1) as well as LaMCTS (Wang et al, ’20) (https://github.com/facebookresearch/LaMCTS/) which certainly scales even farther. \n\n-\tGP-EI is a pretty weak baseline for a GP enabled BayesOpt loop these days (UCB, KG, MVES, etc..), given the vast amount of BayesOpt tooling that exists, particularly from BoTorch.\n\n_More detailed study of VAE-GP (medium):_\n\n Fig. 13 is really interesting and probably deserves space in the main text given the popularity of deep models + GPs for BayesOpt. Perhaps the authors could cut down on sum of the ending discussion in order to move it into the main text and to flesh out the experiment a bit more, with a bit more detail there.\n\n_Claims (medium):_\n\n“Neural Linear and GPs scale cubically with output dimensionality (without the use of covariance approximations),”\n-\tThis claim is still a bit misleading. Various multi-task Gaussian processes aren’t exactly “approximate GPs” or even “covariance approximations”, they’re different model classes than a GP fit over all nt data points, that make different independence assumptions than standard GPs. All covariances are exact in a MTGP implementation, see Bonilla et al, ’08 and also Maddox et al, ’21 for a BayesOpt application. \n\nReference Wang et al, ’18 is all about scaling BayesOpt to tens of thousands of data points, not just pure prediction, and so probably deserves a bit more detail / comparison to your approaches. Even if, in my understanding, it is quite computationally intensive.\n\nDoes performance depend on ensemble size?\n\n_Small fixes:_\n\nSection 1.1: don’t use “Ref”. Strike those words.\n\nFig 13: the comparison with the VAE-GP still probably needs a bit more comparison and explanation. Why is the VAE precisely so helpful on only one task, which seems to be the most difficult to model?\n-\tAlso, why no error bars on the GraphEnsemble (also seems to apply to Fig. 5)?\n\nI feel like I’m missing this somewhere, but how many trials were each experiment run over?\n\n\nSpace after graphs at second to last paragraph of Section 1.\n\n_References:_\n\nEriksson, David, et al. \"Scalable global optimization via local bayesian optimization.\" Advances in Neural Information Processing Systems 32 (2019).\n\nWang, Linnan, Rodrigo Fonseca, and Yuandong Tian. \"Learning search space partition for black-box optimization using monte carlo tree search.\" Advances in Neural Information Processing Systems 33 (2020): 19511-19522.\n\nBonilla, Edwin V., Kian Chai, and Christopher Williams. \"Multi-task Gaussian process prediction.\" Advances in neural information processing systems 20 (2007).\n\nMaddox, Wesley J., et al. \"Bayesian optimization with high-dimensional outputs.\" Advances in Neural Information Processing Systems 34 (2021).\n\n", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "2U27jyFBbc", "reviewer": "Reviewer_Trh3", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "- The problems in the experimental section are interesting (catering to scientific applications) and will be useful for the community if released (currently the attached repo is empty https://github.com/placeholder/deepBO)\n- In terms of methodology, unfortunately there is limited new knowledge in the paper. BNNs have been used earlier in multiple papers for BO and there doesn't seem to be any new technique to train them introduced here. One point about using cosine annealing learning rate (Loshchilov and Hutter) is mentioned when retraining the model on acquiring a new point instead of re-training from scratch. However, Figure 7 shows that the gains from this technique are not very effective.\n- One other big concern is about the baselines. There is no comparison/discussion with any high-dimensional BO approach using GPs like TurBO, SAASBO, BOCK etc. TurBO can easily be scaled to 10000 points much larger than the sizes in the experiment of the paper (1000). \n\nTurBO: \"Scalable global optimization via local bayesian optimization.\" Advances in Neural Information Processing Systems 32 (2019)\n\nSAASBO: David Eriksson and Martin Jankowiak. High-dimensional bayesian optimization with sparse axis- aligned subspaces. In Uncertainty in Artificial Intelligence, pages 493–503. PMLR, 202\n\nBOCK: ChangYong Oh, Efstratios Gavves, and Max Welling. Bock: Bayesian optimization with cylindrical kernels. In International Conference on Machine Learning, pages 3868–3877. PMLR, 2018\n\nThere is no comparison with any of the latent space VAE style approaches as well.\n\nLocal Latent Space Bayesian Optimization over Structured Inputs. CoRR abs/2201.11872 (2022)\n\nCombining Latent Space and Structured Kernels for Bayesian Optimization over Combinatorial Spaces. NeurIPS 2021: 8185-8200\n\nSample-Efficient Optimization in the Latent Space of Deep Generative Models via Weighted Retraining. NeurIPS 2020", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper studies the problem of optimizing expensive black-box functions in high dimensions with BO framework when we have additional auxiliary information for e.g. a-priori known structure of the input domain (e.g. symmetries or invariances). The key idea is to employ Bayesian neural networks (BNNs) surrogate models with model architectures capturing the auxiliary information for e.g. like convolutional NNs to exploit translation invariance in images, or graph NNs to exploit symmetries in molecule.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Overall, I think the paper is at a workshop level right now. However, adding new insights and new qualitative/quantitative discussion with respect to a large body of existing work (as mentioned in the review) will be useful and strengthen the paper.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper demonstrates through numerical examples that Bayesian Optimization (BO) can be applied in relatively high dimensional settings when it is augmented with auxiliary information through Bayesian neural networks surrogate models.

The work advances the understanding of which methods work well in which settings, provides some interesting and challenging benchmarks for BO and establishes sensible baselines for neural network-based Bayesian optimization.

Three reviewers provided detailed feedback on the manuscript, engaged with the authors during the rebuttal and unanimously stated that the paper has become much stronger after the revisions. The reviewers concluded that the acceptance criteria for TMLR are met, and all reviewers recommend acceptance.

Since the reviewer argued that the numerical benchmarks could serve as a reference for future work, I would recommend adding a short table to the camera-ready version that summarizes the key features (such as dimensionality) of the benchmark problems to make this information more readily available to the readers.

Decision: Accept (poster)