BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Authors: Toby Boyne
Jose Pablo Folch
Robert Matthew Lee
Behrang Shafei
Ruth Misener
Venue: ICML 2025 poster (2025)
Keywords: Bayesian optimization
tree kernels
Bayesian additive regression trees
Gaussian processes
Primary area: optimization->zeroorder_and_blackbox_optimization
TL;DR: A fully Bayesian tree kernel provides a richer posterior, providing better probabilistic models for mixed-feature Bayesian optimization.

Abstract: We perform Bayesian optimization using a Gaussian process perspective on Bayesian Additive Regression Trees (BART). Our BART Kernel (BARK) uses tree agreement to define a posterior over piecewise-constant functions, and we explore the space of tree kernels using a Markov chain Monte Carlo approach. Where BART only samples functions, the resulting BARK model obtains samples of Gaussian processes defining distributions over functions, which allow us to build acquisition functions for Bayesian optimization. Our tree-based approach enables global optimization over the surrogate, even for mixed-feature spaces. Moreover, where many previous tree-based kernels provide uncertainty quantification over function values, our sampling scheme captures uncertainty over the tree structure itself. Our experiments show the strong performance of BARK on both synthetic and applied benchmarks, due to the combination of our fully Bayesian surrogate and the optimization procedure.

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Toby Boyne' Jose Pablo Folch
Robert Matthew Lee? Behrang Shafei? Ruth Misener'
Abstract
such as the squared exponential and Matérn kemels are lim- ited to continuous inputs and assume function stationarity.
We perform Bayesian optimization using a Gaus- sian process perspective on Bayesian Additive Regression Trees (BART). Our BART Kernel (BARK) uses tree agreement to define a poste- rior over piecewise-constnt functions, and we explore the space of tree kernels using a Markov chain Monte Carlo approach. Where BART only samples functions, the resulting BARK model obtains samples of Gaussian processes defining distributions over functions, which allow us to build acquisition functions for Bayesian optimiza- tion. Our ree-based approach enables global optimization over the surrogate, even for mixed- feature spaces. Moreover, where many previous tree-based kernels provide uncertainty quantifica- tion over function values, our sampling scheme captures uncertainty over the tree structure itself. Our experiments show the strong performance of BARK on both synthetic and applied benchmarks, due to the combination of our fully Bayesian sur- rogate and the optimization procedure.
Tree regression provides a flexible model which is popular due to strong empirical performance, natural regularization and interprotability (Hill et al., 2020; Loh, 2011). Trees allow for mixed inputs and can model non-stationary func- tions, with hierarchical dependencies. A tree model defines sequence of domain splits. A single input will go through all the splits until a leaf node is reached which assigns value to the predictive model. A forest is an ensemble of trees, where each tree is regularized to be a weak learner to avoid overfitting. Despite their performance in regression tasks, tree models have limited application to BO since they lack explicit uncertainty quantification in their predictions, instead relying on empirical variance.
Given an ensemble of m trees, and two inputs x and x', we can define a function representing how often the trees agree on the inputs belonging to the same leaf. This notion of tree agreement can be used to construct a valid kernel for use in a GP (Balog et al., 2016; Davies & Ghahramani, 2014; Thebelt al., 2022a). This allows us to use tree models in applications which require well-calibrated estimates such Bayesian optimization. Tree kernels for model fitting and BO have previously been explored (Balog et al., 201 6; Thebelt et al., 2022a), but they have always used a two- step fitting process: first, a tree model is fitted on the data, and then the GP is trained using the tree kernel defined by the tree model. While this heuristic gives good results in a variety of benchmarks, it is not principled and hence potentially sub-optimal.
1. Introduction
Bayesian optimization (BO), which finds the optimal condi- tions of expensive-to--vaalute black-box functions (Frazier, 2018; Shahriari et al., 2015; Wang et al., 2023c), is lead- ing approach for data-driven design of experiments. BO applications range from molecule design (Griffiths et al., 2024) to battery and chemical reaction optimization (Folch et al., 2023). BO relies on a probabilistic surrogate model- ing our belief and uncertainty about the black-box function. One class of surrogates are Gaussian processes (GPs), flex- ible models giving well-calibrated uncertainty estimates (Rasmussen & Williams, 2005). The most important hyper- parameter choice when using GP is the kernel, which fully defines the model's covariance structure. Popular kernels
Since the tree structure has an infinite number of parameters, it is at risk of overfitting when optimizing the log-marginal likelihood (Rasmussen & Williams, 2005). To overcome this challenge, we take inspiration from Bayesian tree methods which impose a prior on the tree structure. We obtain a posterior over the tree structure, from which we sample using Markov chain Monte Carlo (MCMC).
We motivate using tree-kernels by exploiting the strong modeling power of tree-based approaches, which can be used for mixed input spaces. Moreover, we can formulate acquisition functions with these tree models as optimization problems that can be solved globally, where gradient- based methods often struggle with discrete inputs.
'Imperial College London (London, UK) "BASF SE (Lud- wigshafen, Germany). Correspondence to: Toby Boyne <t.boyne236 imperial.ac.uk>.
Proceedings of the 4nnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
MCMC
~ (a(
Figure I: BARK uses MCMC create samples of forests, each of which define a Gaussian process. We can then use all the GPs to obtain posterior samples of the process, to marginalize over the tree model uncertainty for the acquisition function.
Contributions. Our work achieves tree-based GP BO through improved Bayesian modeling. Our main contri- butions can be summarized as:
Linero (2017) shows that, under certain priors, BART con- verges to a GP with a Laplace kernel as the number trees increases. Linero (2017) also mentions that the Laplace kernel GP is sub-optimal in practice compared the tree model which has lower computational cost, and better generalization.
We propose fully Bayesian method for training tree Gaussian processes through Markov chain Monte Carlo, by exploiting connections between Bayesian tree models and Gaussian processes with tree kernels. We provide a brief summary of our approach in Fig.
Bayesian optimization over mixed feature spaces has been extensively researched in recent years. Techniques include using RKHS embeddings (Buathong et al., 2020), diffusion kernels (Deshwal et al., 2021), multi-armed bandit style op- timization (Ru et al., 2020), and plethora other methods SCarridoMeerchlle & Hernández-Lobato, 2020; Häse et al., 2021; Daxberger et al., 2021; Daulton et al., 2022; Deshwal et al., 2023; Wan et al., 2021) and benchmarks (Dreczkowski et al., 2023). This surge in interest also includes the pro- vision of software packages for mixed feature BO such as SMAC (Lindauer et al., 2022) and scikit-optimize (Head et al., 2018). More specifically, there also several ap- proaches using tree-based models for black-box optimiza- tion (Ammari et al., 2023; Mistry et al., 2021; Thebelt et al., 2021; 2022b). Jenatton et al. (2017) exploit a known tree structure for more efficient BO. Lei et al. (2021) use BART for BO, however, they do not fully exploit the implicit GP structure in the model as we do. Finally, Thebelt al. (2022a) show how to use GPs with tree kernels for BO with optimal acquisition function optimization, however, they fit the trees using a two-step procedure.
We provide computationally efficient training algo- rithm allowing for scalability and accelerating up the Markov chain sampling procedure.
We provide an empirical analysis of the benefits of our approach, showing its regression ability and achieving state-of-the-ar art results in variety of Bayesian opti- mization benchmarks.
2. Related Work
The relationship between Bayesian tree models and GPs has long history. Chipman et al. (1998; 2010) develop a Bayesian approach for fitting classification and regression trees (BART), which samples tree structure, and leaf values with Gaussian prior. Further work by several authors explores using GPs as tree leaf nodes (Gramacy & Lee, 2008; Maia et al., 2024; Wang al., 2023a), or within other partitions (Luo et al., 2021; 2023). These trees, or other spatial partitioning strategies, allow non-stationary models that may be effective for Bayesian optimization (Eriksson Jankowiak, 2021; Eriksson al., 2019; Papenmeier et al., 2022: Wang et al., 2020; Ziomek & Bou Ammar, 2023).
The rest of the paper will continue as follows: provides background on tree kernels and Bayesian opti- mization. Section 4 provides an overview the BART procedure, laying the groundwork for BARK. Sections and describe our method, highlighting the developments compared prior work and detailing the motivation behind these improvements, specifically in the context of BO. Sec- tion provides empirical results on suite of synthetic and applied BO benchmarks with mixed feature spaces. Sec- tion summarizes our work, addressing the strengths and weaknesses of the method.
Section 3
The use of tree models kernels has been investigated empirically (Thebelt et al., 2022a) and theoretically (Scor- net, 2016). Balog et al. (2016) consider tree kernels in the context of Mondrian processes and show that the Mondrian kernel approximates the Laplace kernel while providing a connection random forests.
2

---

BARK; A Fully Bayesian Tree Kernel for Black-box Optimization
3. Background
3.3. Bayesian Additive Regression Trees
3.1. Gaussian process Bayesian optimization
Bayesian Additive Regression Trees (BART) is a fully Bayesian tree model that has seen a lot of empirical success (Chipman et al., 2010). It is formulated as sum of trees:
Bayesian optimization (BO) maximizes a black-box func- tion, x, arg maxxex fx), through the sequential and possible noisy querying of the function. BO balances the exploration-exploitttinn trade-off via an acquisition fune- tion, The next query is X arg m axx; 0, D,), where Dj represents all the data gathered up to iteration i, and 0 are surrogate model parameters. The feature space X' contains all possible experiments, and may contain continu- ous, integer, or categorical features.
m
y(x) Σ Txx +, C~,(0,0)) (3)
where M, = [ueg.... "......ri. Mrs.] ~ NO, /m/m for some fixed parameter oi. T: represents the tree structure and M, the leaf values. Chipman et al. (2010) show how priors the tree structure enable Bayesian inference using Gibbs sampling (Gelfand, 2000). The model shares a close implicit relationship with the tree kernel. As noted by Linero (2017), the expectation of the BART posterior over the leaf values M (conditioned on a set of trees T) is:
The unknown function is typically modeled with Gaus- sian process (GP) prior. A GP is a probabilistic, nonpara- metric model over functions, defined by a mean and covari- ance function (Rasmussen & Williams, 2005). This prior, combined with a likelihood of data observations, induces a posterior over functions. The kernel parameters are often chosen by optimizing the marginal likelihood, however, it also possible to be fully Bayesian and perform inference over the kernel parameters. This cannot be done analyti- cally, and requires posterior samples of kernel parameters to estimate corresponding integrals. For example, we can marginalize over the kernel parameters when optimizing the acquisition function by taking S posterior samples of the parameters, { (Snoek et al., 2012):
EM [[(x) = Σ EM ((: T() MT
1=1
m
Φ (x: T.) T(x:
m
(4)
Eq. (4) is equivalent the kernel in Eq, (2), conditioned on of of. However, even with the same covariance structure as the tree kernel, all MCMC draws from BART result in deterministic functions, and require many more samples to approximate the distribution that the tree kernel explicitly de- fines. The effectiveness of BART in BO is therefore limited, since we do not have access to the uncertainty quantification required for building standard acquisition functions.
S. a(x; 0, D.)p(((D,)
X1+1 arg max XEX
s
Σ
axx; (()), Dj)
2 arg max. XEX S
8=1
3.2. Forest kernels
The forest kernel uses binary decision trees to define a dis- tribution over piecewise-constant functions. This kernel counts the proportion of its m trees that agree on the leaf node in which a pair of data falls into to define a covariance function. We denote the tth tree by Tú, which contains L, leaf nodes. We define the one-hot vector ((x; T:) where the non-zero entry corresponds to the leaf that contains the datapoint x. This enables the definition of the kernel:
4. The BART model
The BART procedure provides a basis for the BARK sam- pling mechanism. Here, we present an overview of the method for generating MCMC samples of tree functions- see Chipman al. (2010) for a comprehensive description.
The BART model consists of three key ingredients: the pri- ors placed tree functions, the likelihood of observations given a tree function, and the proposals used in sampling the posterior. The priors multiplied by the likelihood of observations induce a posterior p(T, ,M,,o(D) over the sum-of-tree model. This posterior is intractable, but MCMC sampling can obtain posterior samples.
m
1(xxx')
x')
(xx: T.) ((xx; TT),
m
where dg is a scale hyperparameter. Eq. (2) produces posi- tive semi-definite matrices, so is valid PSD kernel (Davies & Ghahramani, 2014). The data x and x' can belong to continuous, categorical, or mixed feature spaces, and the kernel is non-stationary. The forest used in the kernel can be trained independently of the GP (as in Davies & Ghahra- mani (2014), giving a supervised kernel), or fit jointly with the GP, which we propose in our method.
This MCMC algorithm is largely defined by the Metropolis- Hastings step (Hastings, 1970). Steps are taken in the hyper- parameter space according to a transition kernel, q(0 - 0*), which are rejected/accepted according to the acceptance probability a of the new state. By taking many steps, the samples converges to the posterior over 0. The acceptance
3

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
probability is decomposed into three terms,
9(0" 8) ( + 0")
p(y 0*) y 0)
p(0*)
a(0, 0*) - min
p(0)
transition
likelihood
prior
BART prior. The regularization prior defined in the BART model is prior on the depth of the nodes of a tree. The probability of any given node being a decision node is:
a(1 + dd) a E (0, 1), 3 [0, 00)
(5)
The number of trees m in the forest is fixed, e.g., m = 50 (Kapelner& & Bleich, 2016). To assign a decision rule to the node 77 in the tree, a sample is drawn from the decision rule prior. The feature on which node is split is drawn uni- formly, and the split value drawn uniformly from the set of unique values in the data that reaches node r7. The prior the observation noise, parameterized by (M, q), is given by of ~ InverseGamma(y/ 22, vt) For set of observations with variance σ², t chosen such that (r < d²) q-
Figure 2: Allowing the tree kernel to vary between samples produces a richer posterior distribution, whilst samples are still piecewise-constane (top) samples from a GP with a fixed tree kernel, and (bottom) samples from BARK.
forests (Breiman, 2001). Second, they maximize the model likelihood (ML-III) to obtain the noise and scale GP hyper- parameters. The tree model is fit using different objective than the GP, meaning the resulting tree structure may not necessarily maximize the likelihood of the data. Rather than fitting the tree and the GP in two steps, BARK jointly fits the tree kernel and the GP.
BART likelihood. The distribution of set of observations under the assumption of Gaussian noise, is
y/X, T, T,M, W
Σ M, (X)
Point estimates of tree posterior. Some existing tree kernel methods use one-step algorithm to fit the kernel and treat the noise hyperparameter fixed, for example (Cohen et al., 2022). However, existing methods only fit single tree model which a maximum a posteriori estimate. Since the tree kernel nonparametric, the GP effectively has infinite hyperparameters. The ML-II estimate tends overfit with many parameters (Rasmussen & Williams, 2005), which can be mitigated by a fully- Bayesian treatment of the kernel (Ober et al., 2021). Moreover, the posterior distribution over functions defined by the model is less rich than the full posterior over tree models, as the tree structure fixed for all functions sampled from the model, see Fig. 2.
The likelihood is conditioned on the sampled leaf values M. Since all other leaf values are fixed when sampling the jth tree, the authors compute the likelihood in terms of Rj- y Exe(mpy M. Ti)
BART proposals. To explore the posterior of trees. Chip- man et al. (2010) define set of proposals. Each proposal makes a small change to a tree's ' structure: growing leaf node by assigning a decision rule from the rule prior; prun- ing a pair of leaf nodes; changing the rule of singly internal node (a decision node where both children are leaves). We omit the swap rule in line with other BART implementations (Kapelner & Bleich, 2016; Linero& & Yang, 2018).
BARK therefore defines probabilistic model of the tree kernel, and performs MCMC explore the posterior of trees. We consider the tree model as nonparametric kernel. The BARK hyperparameters are the parameters defining the forest generating process, i.c. the parameters of the tree priors. These parameters are interpretable, and can be used with default values.
5. The BARK model
5.1. Motivating the Bayesian treatment of the kernel
Despite existing work on tree kernels for Gaussian pro- cesses, there is still limited literature on tree kernels in Bayesian optimization setting. We consider two key areas for improvement.
5.2. Differences from the BART model
Two-step fitting. Tree kernels such as those studied by Davies & Ghahramani (2014) and Thebelt et al. (2022a) fit the GP model in two steps. First, they fit tree indepen- dently, using gradient boosting (Friedman, 2001) or random
While we take inspiration from BART, relying upon its well established predictive performance, BARK has some key differences from the BART model. We explore these by

---

BARK; A Fully Bayesian Tree Kernel for Black-box Optimization
considering the three terms of the acceptance probability. BARK prior. In Chipman et al. (2010), the authors nor- malize observations then place a prior on the variance of the leaves as z6. where k is to be tuned. We remove the k parameter, instead opting to standardize the outputs and fix the scale of = 1, which has been shown to work well in Bayesian optimization (Hvarfner et al., 2024).
walk in an unconstrained space, using the softplus transform. Details for the BARK proposals are given in Appendix H.
5.3. Computational consideration
When we have sampled N data points, the matrix inversion required in the fitting procedure for BARK has ( com- plexity. This cubic order is unavoidable due inverting the covariance matrix when computing the log likelihood in the MCMC steps. However, we can reduce the cost in comput- ing this term for tree proposals by leveraging the low-rank nature of the contributions of individual trees. This is sim- ilar to other methods used for fast inverses of tree kernels (Balog al., 2016; Lee al., 2015), applied to sequential updates to trees in forest.
BARK adopts the same inverse gamma distribution for the noise, 2 o. as BART. Since the data is standardized, the noise
prior will be the same for any dataset.
BARK likelihood The marginal likelihood of the Gaussian process is computed from y distributed as
( - t=1
v/X,T,0,~> ~ N
Σ o (X: ((;;)) of
For dataset X [x()) ... ( we define & E (0, 1) (mxXXXL where Pjni = (m) Te)le and T. = max; L. This gives the following expression for the kernel:
Crucially, our kernel perspective means that we never sam- ple leaf values- the marginal likelihood integrates over all the values of the leaves. In this way, BARK samples distributions, which can then be used to build acquisition functions, motivating our use of the GP likelihood. In con- trast, BART samples the leaf values sequentially, hence uses the likelihood conditional on leaf values. In order to ap- proximate the posterior leaf distribution, BART requires a large number of function samples, which in turn means that it is infeasible to formulate optimization problem over the BART acquisition function. Both BART and BARK have the same probabilistic model: marginalizing the BART likelihood over leaf values yields the BARK likelihood, as in Eq. (4).
m
Ko
du m
3 + +
During the MCMC procedure, for a proposal where the 8th tree is updated, the difference in log- likelihoodd between the previous GP and the proposed GP can be computed:
p(y p 8*) p(y 0)
p = y [Ko- -1 Ko ")y
log
- -logg Kýl - Log Kol)
This requires computing Kg. and log Kg both of which are ((N3) operations. However, Kg and log Kgl will have already been computed in the previous iteration, and so we can rewrite this as matrix update,
BARK proposals. In BART, splitting points are typically uniformly sampled from the data seen by a given node (Chip- man et al., 2010; Kapelner & Bleich, 2016). This choice is made such that any leaf node will contain nonempty subset of the domain (no 'logically empty' leaves). How- ever, this approach performs poorly in the BO setting, since finding potential future queries requires good uncertainty estimation far from data (see Section 6). Instead, we sample the splitting rule uniformly from the domain. To avoid cre- ating logically empty nodes, each splitting rule is sampled from the subspace of the domain that reaches that node. For example, if x E [0, 1]. and the root splitting rule is æ > 0.5, then the splitting value of the left child is be sampled from [o, 0.5). For categorical features with S categories that reach node node, the splitting rule is sampled uniformly from P(S) {S, (), i.c., the power set of S excluding trivial splits. This enables greater sample efficiency than a one-hot encoding of categories. We investigate the impact of this modeling choice in Appendix B for toy regression examples.
o
Kg- m
D
+ m
+
jE(m
of 4 क T T
Kg m
T D.
m
Since each o, matrix has rank L,, and L, < N by the depth prior, an update to tree structure gives two low-rank updates. This can be exploited using the Aheemman-Morriso formula (Woodbury, 1950), and the matrix determinant lemma (Harville, 1997), to compute the matrix inverse of the updated kernel matrix in O(N L.) complexity:
(A= = AI UUU A UU 'TAA"
log/ A AU logl A| + log / + U A-
6. BARK for Bayesian optimization
BARK cannot sample directly from the posterior over o3 and so must sample the noise using MCMC. To enforce non-negativity, a proposal is generated by taking a Gaussian
While BARK can be used for regression, its strength lies in its performance in Bayesian optimization. In this section, we outline how the model has been designed for BO.
5

---

BARK; A Fully Bayesian Tree Kernel for Black-box Optimization
6.1. Optimizing the acquisition function
2
Typical BO approaches use gradient-based methods opti- mize the acquisition function (AF), such L-BFGS (Zhu et al., 1997). However, our piecewise constant prior has zero gradient almost everywhere. We turn to mixed integer linear programming to define the optimization problem. We use the framework provided by Thebelt et al. (2022a); Mišić (2020) to encode the tree structure using constraints on bi- nary variables. The full formulation is given in Appendix This has the additional benefit of global optimization, where gradient-based methods can be stuck in local optima.
CP prediction
1
2
avcn (xc
1
Integrated acquisition function. Each sample from the MCMC gives the parameters (T, ,y), defining a GP kernel and likelihood. function drawn from BARK is drawn from one of these GPs with equal probability; a BARK predictive distribution is mixture of Gaussians (Lalchand & Rasmussen, 2020). We therefore approximate the integrated acquisition function from Eq. (I).
-1
XXXXX
XXXXXX
0.0
0.5
1.0
20
Uniform split
Data split
AB Maximiser
Figure 3: GP predictions and UCB for two BARK GPs on minimization task: (blue) splits sampled uniformly, and (orange) splits sampled at datapoints. By sampling splits uniformly, the acquisition function (AF) maximizer exhibits better exploration.
We use the Upper Confidence Bound (UCB), which is a function of the mean, μI, and standard deviation, o, of the GP, OccB (x) = μ(x) + (x). The integrated UCB a convex function of fı and a, which leads to significant speed increase in solving the global optimization problem. After the optimization is solved and new datapoint is queried, we sample new tree kernels from the MCMC, using the final sample from each parallel chain of the previous BO iteration as an initialization.
the tree structure prior will concentrate on regions that are dense with observations, which leads to over-exploitation.
6.3. The Laplace approximation for regret bounds
The functions in the BARK prior are discontinuous and non- stationary. This violates the typical assumptions required for regret bounds in the Bayesian optimization literature (Srinivas et al., 2012). Previous work shows that, under some strong assumptions, the tree kernel can be approxi- mated as a Matern- kernel (Linero, 2017), which in turn admits sub-linear regret bounds of O(T where D is the dimension of the domain (Wang al., 2023b). Ap- pendix I provides further discussion on the quality of this approximation.
6.2. Uniform splitting rule
In Section 5, we describe the new prior over the splitting rules at nodes, where we sample the split values uniformly across the domain. While the original prior worked well in the high data regression setting, we show the importance of this change in the BO setting to encourage exploration.
A key benefit of the change is the better uncertainty quantifi- cation. When sampling splits only at datapoints, the uncer- tainty is constant in the region between datapoints. This no longer has the property that uncertainty increases as the dis- tance from observation increases, an essential ingredient in a well-performing BO surrogate. Moreover, the predictive uncertainty outside the range of observations constant, leading to poor extrapolation ability. This is demonstrated in Fig. 3; due to the poor uncertainty quantification when splits are sampled at datapoints, the UCB maximizer fails to properly explore the domain.
7. Experiments
This section describes experiments comparing the regression capabilities of BARK and BART, and evaluates the BARK model against selection of baselines on set of synthetic and applied Bayesian optimization benchmarks'.
Baselines. We compare BARK to several baselines. For each baseline, we use the upper confidence bound (UCB) AF BART uses the Chipman et al. (2010) surrogate, evaluated grid, and uses function samples to estimate the UCB (Wilson et al., 2018). GP-RBF uses the squared-expomential
further issue with sampling splits from the data is the ten- deney to cluster queries. Since BO aims find the optimum value of function, any exploitative function queries will cluster around known good values. This leads to clusters of observations near the optima. If the splitting rule prior assigns equal splitting probability to each datapoint, then
'The code these experiments is available at https: //github. thub. com wwwm//obyooyn///aboee
6

---

BARK; A Fully Bayesian Tree Kernel for Black-box Optimization
kernel, with the default priors from Hvarfner et al. (2024). For enumerable mixed feature spaces, we use linear combi- nation of sum and product kernels, as in Ru et al. (2020). To optimize over mixed feature spaces, we use an alternating approach to perform optimization, similar to (Wan et al., 2021), where we alternate between optimizing the continu- ous parameters with gradient-based methods, and exploring local perturbations in the discrete parameters. For domains where the number of combinations of categories is low, we also include enumeration approach, where each unique combination of discrete values is fixed and many parallel continuous optimizations are performed. LeafGP uses a two-step forest kernel, where the forest is fit using gradient boosted trees, then used as a kernel in a GP (Thebelt et al., 2022a). This model can be used in global optimization of the UCB acquisition function. SMAC uses an ensem- ble of trees, using sample variance to measure uncertainty (Lindauer et al., 2022). Entmoot defines mean function with tree surrogate, and uses distance-based uncertainty (Thebelt et al., 2021). Finally, BARKPrior samples tree kernels from the prior, without performing the MCMC steps to sample the posterior. Appendix D provides further details of the methods.
AutoMPG
Abalone
MLI
2
Sample
300
Sample
300
1.0
0.5
0.0
0
Lag, T
50
Lag, 50
Figure 4: Evolution of the marginal log likelihood (MLL) of the training data over MCMC samples, and the correspond- ing autocorrelation, p(T). Lines are plot for 20 MCMC chains, with a single random chain shown in solid colour.
entirely uncorrelated. We find that the samples become un- correlated after 50 samples, which suggests that the chosen thinning rate of 100 is sufficient.
7.1. Model fit comparison
7.2. Synthetic benchmarks
We perform regression with both BART and BARK. The purpose of this section is to demonstrate that we do not lose predictive power when using the kernel perspective of the BARK model. We show this in the tabular regression setting, where BART is already known to be strong. We use the default setting for the model hyperparameters: for discussion on the selection of these values, see Appendix G.
We perform Bayesian optimization on four mixed-space syn- thetic benchmarks, initialized with min min 30) datapoints sampled uniformly from the domain. For each method, we report the median and inter-quartile range across 20 runs.
TreeFunction and TreeFunctionCat are functions sampled from the BART prior, each with 10 continuous dimensions, and the latter with an additional 10 categorical dimensions. DiscreteAckley and DiscreteRosenbruck are partially dis- cretized functions from Dreczkowski et al. (2023); Blick et al. (2021). Appendix D has further benchmark details.
We use real datasets from the UCI Repository (Cortez & Silva, 2008; Nash et al., 1994; Quinlan, 1993; Yeh, 1998; Dua & Graff, 2017). We compare the predictive perfor- mances in Table 1. The performance of the two models is similar, suggesting that BARK can achieve similarly strong regression performance to BART. Our slight increase in NLPD performance may be due to our GP perspective, which effectively allows us to marginalize over the distri- bution of leaf values. Appendix F discusses the training time on these regression problems, and Table 7 provides comparison to the GP-RBF model.
We show the results in Fig. 5. A key takeaway here is that, despite TreeFunction and TreeFunctionCat being drawn from the BART prior, the BART baseline does not outper- form LeafGP and BARK. This shows the strength of the optimization formulation, and that a grid-based UCB evalu- ation is insufficient to find the global optimum. BARK both captures uncertainty over the tree structure, and optimizes the UCB efficiently, to achieve the strongest performance.
We further investigate the trajectory of the marginal log likelihood of the MCMC chain in Fig. 4. As the MCMC process makes local proposals, each sample is highly corre- lated with the previous sample, and the correlation decreases as the distance between samples (or lag) increases. This is quantified by the autocorrelation of the samples (X}s = p(T) = [x++++++ (X, / where is the lag.
DiscreteRosenbruck shows a failure case of BARK; the optimum lies very shallow region, which the piecewise- constant prior of BARK models worse than the smooth prior of the RBF, and Entmoot's distanced-based uncertainty quantification. These synthetic benchmarks tend to favor the smoothness assumptions of the RBF kernel over the sepwise-constand functions in the BARK posterior. We provide further results for continuous synthetic functions in Appendix E.2.
Since we are limited by the optimization procedure to a small number of samples, we want our samples to be almost
7

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Table 1: Benchmark performance of the BART model against BARK, given as mean (and standard deviation). The number
of training points used for each benchmark is given by n. The number of dimensions for each problem is given by D (continuous integer categorical). The better performing model bold for each metric.
Benchmark
D
NLPD /
MSE
BART
BARK
BART
BARK
Abalone
400
7+0+1
1.10 (0.06)
1.09 (0.06)
0.53 (0.06)
0.52 (0.06)
Auto MPG
100
4+3+0
0,46 (0.12)
0.45 (0.08)
0.14 (0.02)
0.15 (0.02)
Concrete
300
8+0+0
0.54 (0.02)
0.31 (0.06)
0.17 (0.01)
0.11 (0.01)
Student Performance
250
0+13+17
1.27 (0.09)
1.27 (0.06)
0.72 (0.08)
0.73 (0.07)
Hest abserved vallor
Trossunction
Troee ProePunntiinCaa
Best obaurved value
SVRBench
XGBooe
0.06
247
-10
246
0.05
-15
246
0.01
264
-20
50
100
50
100
100
Best ubserved value
DiscreteAckley
DiseroteRasembroct
CCOBench
Contral
0.20
16
3.5
3
Best observed
0.18
3.4
0.16
14
0.14
3.2
100
50
100
0
50 Iteration
100
Iteration
Iteration
Iteration
BARK
BART SMAC
BARKPrior
CP-RRP
BART
BARKFriar
LeutGP
Eatmoot
LeatGP
SMAC
Eutmoot
GP-RBF (esum)
Figure 5: Optimization of synthetic benchmarks. The shaded regions contain the 25th and 75th percentile of regret achieved across the 20 runs.
Figure 6: Optimization of applied benchmarks. XGBoostM- NIST allows enumeration GP-RBF, whereas other bench- marks require alternating optimization.
bination, GP-RBF is strong choice. However, in the other applied benchmarks where the alternating optimization ap- proach must be employed, BARK is highly competitive with GP-RBE
7.3. Applied benchmarks
We evaluate these methods on applied benchmarks. We use the hyperparameter optimization benchmarks SVRBench and XGBoostMNIST (Dreczkowski et al., 2023). We further evaluate using CCOBench (Dreifuerst et al., 2021), which optimizes the configuration of antennas maximize net- work coverage, and PestControl (Dreczkowski al., 2023), optimizing the choice of pesticide used at 25 different sta- tions. Fig. 6 gives the results.
Bandit optimization with material design. BART has been used to perform material selection using Bayesian optimization in Lei et al. (2021). Given dataset of 403 candidate compounds (specifically, class known as MAX phases), we maximize the bulk modulus and minimize the shear modulus (Talapatra et al., 2018). We present a compar- ison of BARK's performance on this task, alongside BART and GP-RBE, in Fig. 7.
BARK outperforms BARKPrior, which itself is a strong method. This shows that the global UCB optimization over tree kernels alone leads to good optimization performance, but proper inference leads to greater improvements. BARK also consistently outperforms LeafGP in both synthetic and applied problems, demonstrating the importance of captur- ing the uncertainty over tree structure.
Since this benchmark has (small) finite number of can- didates, it reveals that BART and BARK have similar BO performance when the acquisition function optimization is ablated. Furthermore, we see that the tree-based methods are better surrogates for this problem than the RBF product kernel, as BARK and BART identify the best candidate in fewer iterations than GP-RBE.
GP-RBF achieves the best final value on XGBoostMNIST.
When is feasible to enumerate over every categorical com-
8

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
MAX (balks modulus)
MAX (sheur modulus)
8. Conclusion
We present BARK, a fully Bayesian Gaussian process tree kemnll, that jointly trains the tree structure with the GP. This model explores similar posterior to the additive tree model BART, but provides explicit uncertainty quantification, mak- ing BARK strong surrogate for Bayesian optimization. We show competitive performance on BO, especially in com- binatorial settings, due to its strong surrogate and efficient optimizer, beating state-of sece--ff-heeatte the-art approaches.
Mean stupple zegret
Best ideut %
80
25
Iteration
Iteration
BARK
BART
GP-RUP
Impact statement
Figure 7: Optimization over compounds in the MAX dataset. We report the mean of the simple regret across 20 runs (lines), and the fractions of the runs that identify the com- pound with the best objective value (bars).
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Acknowledgments
Fitting
Optimization
We would like thank Alexander Larionov for support in implementing benchmarks and baseline methods. Funding for this work was provided by BASF SE, Ludwigshafen am Rhein, and EPSRC through the Modern Statistics and Statistical Machine Learning CDT for TB (EP/Y0348) 13/1) and JPF (EP/S023151//) TB also funded by the EPSRC IConIC Prosperity Partnership (EP/X025292/1) RM holds concurrent appointments
10a
Time ($)
10
10-3
50
100
100
Professor Imperial and as Amazon Scholar. This paper describes work performed at Imperial prior to joining Amazon and is not associated with Amazon.
Number iterations
BARK
BART
DARKPrioe
GP-REF GP-RBF (enum)
LeaUP
SMAC
Enimoot
Figure& 8: Time taken fit the surrogate and optimize the acquisition function for the XGBoostMNIST benchmark.
References
Ammari, B. L., Johnson, E. S., Stinchfield, G., Kim, T. Bynum, M., Hart, W. E., Pulsipher, J., and Laird, C. D. Linear model decision trees as surrogates optimization of engineering applications. Computers Chemical En- gineering, 108:10 October 2023. ISSN 0098-1354, doi: 10. 10.1016f.compphemmeee. 0000mmmppccmmmmm002200000433 2023..0347
7.4. Wall-time discussion
BARK's most significant limitation its wall-time perfor- mance. Fitting the model takes approximately 50 seconds, and the optimization of the acquisition function is limited to 100 seconds. This is due to the combination of the expensive MCMC procedure, and the large optimization formulation. We compare the fitting and optimization times against other methods in Fig. 8, and provide further selection of wall- time comparisons in Appendix F.2.
Balandat, M., Karrer, B., Jiang. D. R., Daulton, S., Letham, B., Wilson, A. G., and Bakshy, E. BoTorch Framework for Efficient Monte-Carle Bayesian Optimization. In Advances in Neural Information Processing Systems 33, 2020.
therefore recommend BARK for BO settings where the objective expensive to evaluate, taking least several minutes, or having large associated financial cost. Here, the benefit of an accurate model with strong uncertainty quantification is worth the cost paid in wall-time. Note that the PestControl and MAX (material design) benchmarks reflect examples of such black-box functions. For settings where experiments are cheap and/or function evaluations are quick, we recommend alternate methods.
Balog, M., Lakshminarayanan, B., Ghahramani, Z., Roy, M., and Teh, W. The Mondrian Kernel. 32nd Conference on Uncertainty in Artificial Intelligence (UAI), June 2016.
Blick, L.,, Guijt, A., Verwer, S., and de Weerdt, M. Black- box mixed-variable optimisation using surrogate model that satisfies integer constraints. Proceedings of the

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Genetic and Evolationary Computation Conference Com- panion, GECCO '21, pp. 1851-1859, New York, NY,
Deshwal, A., Ament, S., Balandat, M., Bakshy, E., Doppa,
J. R., and Eriksson, D. Bayesian Optimization over High- Dimensional Combinatorial Spaces via Dictionary-based Embeddings. In Ruiz, F., Dy, J., and van de Meent, J.-W. (eds.), Proceedings The 26th International Conference on Artificial Intelligence and Statistics, volume 206 of Proceedings of Machine Learning Research, pp. 7021- 7039. PMLR, April 2023.
USA, 2021. Association for Computing Machinery. ISBN
9781450383516. doi: 10.1 00114444444722334433333
Breiman, I.. Random Forests. Machine Learning, 45
(1):5-32, 2001. ISSN 0885-6125. doi: 10.1023/a: 1010933404324.
Buathong, P., Ginsbourger, D., and Krityakierne, T. Kernels over sets of finite sets using RKHS embeddings, with application to Bayesian (combinatorial) optimization. In AISTATS, pp. 2731-2741 PMLR, 2020.
Dreczkowski, K., Grosnit, A., and Ammar, H. Frame- work and Benchmarks for Combinatorial and Mixed- variable Bayesian Optimization. In l'hirty-seventh Confer- ence on Neural Information Processing Systems Datasets and Benchmarks Track, 2023.
Chipman, A., George, E. I., and McCulloch, R. E. Bayesian CART Model Search. Journal of the Ameri- can Statistical Association, 93(443):935-948, September 1998. ISSN 1537-274X. doi: 10. 10.1080001621459.1999 10473750.
Dreifuerst, R. M., Daulton, S., Qian, Y., Varkey, P., Balandat, M., Kasturia, S., Tomar, A., Yazdan, A., Ponnampalam, V., and Heath, R. W. Optimizing Coverage and Capacity in Cellular Networks using Machine Learning. In ICASSP 2021 2021 IEEE International Conference Acoustics, Speech and Signal Processing (ICASSP), pp. 8138-814 8142, 2021. doi: 10.1 10.1100/CAAS 94111555.
Chipman, H. A., George, E. I., and McCulloch, R. E. BART: Bayesian additive regression trees. The Annals of Applied Statistics, 4(1), March 2010. ISSN 1932-6157. doi: 24.111400-aoas2288
Dua, D. and Graff, C. UCI Machine Learning Repository, 2017. URL http: /archivee... ics edu/m1.
Cohen, M. K., Daulton, S., and Osborne, M. A. Log-linear- time Gaussian processes using binary tree kernels. In Proceedings of the 36th International Conference on Neural Information Processing Systems, NIPS '22, Red Hook, NY, USA, 2022. Curran Associates Inc. ISBN 9781713871088
Dürholt, P., Asche, S., Kleinekorte, J., Mancino-Ball, G., Schiller, B., Sung. S., Keupp. J., Osburg, A., Boyne, T., Misener, R., Eldred, R., Costa, W. S., Kappatou, C.,
Lee, R. M., Linzner, D., Walz, D., Wulkow, N., and Shafei, B. BoFire: Bayesian Optimization Framework Intended for Real Experiments, 2024.
Eriksson, D. and Jankowiak, M. High-dimensional Bayesian optimization with sparse axis-aligned subspaces. In UAI, volume 161, pp. 493-503, 2021.
Cortez, P. and Silva, A. M. G. Using data mining to predict secondary school student performance. In Pyschology, 2008. URL https:
semant org/CorpusID: 166211296 6621299.
Eriksson, D., Pearce, M., Gardner, J., Turner, R. D., and Poloczek, M. Scalable Global Optimization via Local Bayesian Optimization. In NeurlPS, volume 32, 2019.
Daulton, S., Wan, X., Eriksson, D., Balandat, M., Osborne, M. A., and Bakshy, E. Bayesian Optimization over Dis- crete and Mixed Spaces via Probabilistic Reparameteriza- tion. In Koyejo, S., Mohamed, S., Agarwal, A., Belgrave, D., Cho, K., and Oh, A. (eds.), Advances in Neural Infor- mation Processing Systems, volume 35, pp. 12760-12774. Curran Associates, Inc., 2022.
Folch, J. P., Lee, R. M., Shafei, B., Walz, D., Tsay, C., der Wilk, M., and Misener, R. Combining multi- fidelity modelling and asynchronous batch Bayesian Op- timization. Computers & Chemical Engineering, 172: 108194, 2023.
Frazier, P. I. A tutorial on Bayesian optimization. arXiv preprint arXiv:1807.0283 2018
Davies, A. and Ghahramani, Z. The Random Forest Kernel and other kernels for big data from random partitions, 2014.
Friedman, H. Greedy function approximation: gra- dient boosting machine. The Annals of Statistics, 29 (5), October 2001. ISSN 0090-5364, doi: 1214//s/ 1013203451.
Daxberger, E., Makarova, A., Turchetta, M., and Krause, A. Mixed-variable Bayesian optimization. In IJCAI, pp- 2633-2639, 2021.
Garrido-Merchán E. C. and Hermández-Lobato, D. Dealing with categorical and integer-valued variables in Bayesian optimization with Gaussian processes. Neurocomputing, 380:20-35, 2020.
Deshwal, A., Belakaria, S., and Doppa, J. R. Bayesian optimization over hybrid spaces. In ICML, pp. 2632- 2643. PMLR, 2021.
10

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Gelfand, A. E. Gibbs Sampling. Journal of the American Statistical Association, 2000. ISSN 1537-274X. doi: 10. 10.10800016214592200 10474335.
34th International Conference Machine Learning, vol- ume 70 of Proceedings of Machine Learning Research, pp. 1655-1664. PMLR, August 2017.
December
Kapelner, A. and Bleich, J. bartMachine: Machine Learning with Bayesian Additive Regression Trees. Journal of Statistical Software, 70(4), 2016. ISSN 1548-7660. doi: 10. (0.18637/ss..0700..0
Gramacy, R. B. and Lee, H. K. H. Bayesian Treed Gaussian Process Models With an Application to Computer Mod- eling. Journal of the American Statistical Association, (03(48)):11199111 September 2008. ISSN 1537-274X. doi: 10.1198001621455800000000000
Lalchand, V. and Rasmussen, C. E. Approximate Infer- ence for Fully Bayesian Gaussian Process Regression. Proceedings of The 2nd Symposium on Advances in Approximate Bayesian Inference, volume 118. pp. 1-12. PMLR, 2020.
Griffiths, R.-R., Klarner, L.., Moss, H., Ravuri, A., Truong, S., Du, Y., Stanton, S., Tom, G., Rankovic, B., Jamasb, A., et al. Gauche: A library for Gaussian processes in chemistry. Advances in Neural Information Processing Systems, 36, 2024.
Lee, D., Park, H., and Yoo, C. D. Face alignment using cascade Gaussian process regression trees. In CVPR, pp. 4204-4212.8 2015.
Harville, D. A. Matrix Algebra From Statistician 's Per- spective. Springer New York, 1997.
Lei, B., Kirk, T. Q., Bhattacharya, A., Pati, D., Qian, X., Ar- royave, R., and Mallick, B. K. Bayesian optimization with adaptive surrogate models for automated experimental de- sign. upj Computational Materials, 7(1), December 2021. ISSN 2057-3960. doi: 10. (..1038/s41524 0000:/4412---21-00000 021-00662-x.
Häse, F., Aldeghi, M., Hickman, R. J., Roch, L. M., and Aspuru-Guzik, A. Gryffin: An algorithm for Bayesian optimization of categorical variables informed by expert knowledge. Applied Physics Reviews, 8(3), 2021.
Hastings, W. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1):
Lindauer, M., Eggensperger, K., Feurer, M., Biedenkapp. A., Deng, D., Benjamins, C., Ruhkopf. T., Sass, R., and Hutter, F. SMAC3: A Versatile Bayesian Optimization Package for Hyperparameter Optimization. Journal of Machine Learning Research, 23(54):1-9, 2022.
97-109, 1970.
Head, T., MechCoder, Louppe, G., Iaroslav Shcherbatyi, Feharras, Zé Vinícius, Cmmalone, Schröder, C., Nel215, Campos. N., Young. T., Cereda, S., Fan, T., Rene-Rex, Kejia (KJ) Shi, Schwabedal, J., Carlosdanielcsantos, Hvass-Labs, Pak, M., SoManyUsernames Taken, Call- away, F., Estève, L., Besson, L., Cherti, M., Karlson Pfannschmidt, Linzberger, F., Cauet, C., Gut, A., Mueller, A., and Fabisch, A. scikii-ptimizzekcckittoptiee v0.5.2, 2018.
Linero, A. R. A review of tree-based Bayesian methods. Communications for Statistical Applications and Meth- ods, 24(6):543-559, November 2017. ISSN 2383-4757. doi: 10.29220ccam..011.24..0
Linero, A. R. and Yang, Y. Bayesian Regression Tree En- sembles that Adapt to Smoothness and Sparsity. Jour- nal of the Royal Statistical Society Series B: Statistical Methodology, 80(5): 1087-1110, September 2018. ISSN 1467-9868. doi: 10. 1/rssb. 12293.
Hill, J., Linero, A., and Murray, J. Bayesian Additive Regression Trees: A Review and Look Forward. An- nual Review of Statistics and Its Application, 7(1):251- 278, March 2020. ISSN 2326-831X.
doi: 1461
Loh, W. Classification and regression trees. WIREs Data Mining and Knowledge Discovery, 1(1):14-23, January 2011. ISSN 1942-4795. doi: 10. 10.1002/widm.8.
mmnurev-ttaistiis-011119000
Hvarfner, C., Hellsten, E. O., and Nardi, L. Vanilla Bayesian Optimization Performs Great in High Dimensions. In Salakhutdinov, R., Kolter, 7.., Heller, K., Weller, A., Oliver, N., Scarlett, J., and Berkenkamp. F. (eds.), Pro- ceedings of the 41 International Conference on Machine Learning, volume 235 of Proceedings of Machine Learn- ing Research, pp. 20793 20793-20817. PMLR, 21-27 Jul 2024.
Luo, Z. T., Sang, H., and Mallick, B. A Bayesian Contigu- Partitioning Method for Learning Clustered Latent Variables. Journal of Machine Learning Research, 22 (37):1-52, 2021.
Luo, Z. T., Sang, H., and Mallick, B. A Nonstationary Soft Partitioned Gaussian Process Model via Random Spanning Trees. Journal of the American Statistical As- sociation, 119(547):2105-211 2116, October 2023. ISSN 1537-274X. doi: 10. 159.00800162115922222222.2
Jenatton, R., Archambeau, C., González, J., and Seeger, M. Bayesian optimization with tree-structured dependencies. In Precup, D. and Teh, Y. W. (eds.), Proceedings of the

---

BARK; A Fully Bayesian Tree Kernel for Black-box Optimization
Maia, M., Murphy, K., and Parnell, A. C. GP-BART: A novel Bayesian additive regression trees approach using Gaussian processes. Computational Statistics Data Anal- ysis, 190: 107858, February 2024. ISSN 0167-9473. doi: (..0116jj.sda..2222 2tt.ddaa222..178555. 107858.
Shahriari, B., Swersky, K., Wang. Z., Adams, R. P., and De Freitas, N. Takinn the human out of the loop: A review of Bayesian optimization. Proceedings of the IEEE, 104(1):148-177) 2015.
Snoek, J., Larochelle, H., and Adams, R. P. Practical Bayesian Optimization of Machine Learning Algorithms. In Pereira, F., Burges, C., Bottou, L.., and Weinberger, K. (eds.), Advances in Neural Information Processing Systems, volume 25. Curran Associates, Inc., 2012.
Mistry, M., Letsios, D., Krennrich, G., Lee, R. M., and Mis- ener, R. Mixed-Integer Convex Nonlinear Optimization with Sradient-Boostter Trees Embedded. INFORMS Jour- nal on Compuring, 33(3):1003-11 1119 July 2021. ISSN 1526-5528. doi: 10. 118.2///0:.22220009999
Sobol, I. M. On the distribution of points in a cube and the approximate evaluation of integrals. Ussr Computational Mathematics and Mathematical Physics, 7:86-112, 1967.
Misić, V. V. Optimization of Tree Ensembles. Opera- fions Research, 68(5):1605-1622 September 2020. ISSN 1526-5463. doi: 10. 10.1287/0pr:221111922...
Srinivas, N., Krause, A., Kakade, S. M., and Seeger, M. W. Information-Thewne
Nash, W., Sellers, T., Talbot, S., Cawthorn, A., and Ford, W. The Population Biology of Abalone (Haliotis species) in Tasmania. I. Blacklip Abalone (H. rubra) from the North Coast and Islands of Bass Strait. Sea Fisheries Division, Technical Report No, 48, January 1994.
Regret Bounds for Gaussian Pro- cess Optimization in the Bandit Setting. IEEE Trans- actions on Information Theory, 58(5):3250-33) May 2012. ISSN 1557-9654, doi:
2011.2120233.2.
Surjanovic, S. and Bingham, D. Virtual Library of Simu- lation Experiments: Test Functions and Datasets. Re- trieved January 29, 2024, from http: / /www//wwwwww// sfu. jano, 2013.
Ober, S. W., Rasmussen, C. E., and van der Wilk, M. The promises and pitfalls of deep kernel learning. In UAI, pp- 1206-1216. PMLR, 2021.
Papenmeier, L., Nardi, L., and Poloczek, M. Increasing the Scope as You Learn: Adaptive Bayesian Optimization in Nested Subspaces. In NeurlPS, volume 35, pp. 11586- 11601, 2022.
Talapatra, A., Boluki, S., Duong, T., Qian, X., Dougherty, E., and Arróyave, R. Autonomous efficient experiment de- sign for materials discovery with Bayesian model averag- ing. Physical Review Materials, 2(11):11 113803,, November 2018. ISSN 2475-9953. doi: 10. 10.100//pyysreeaaatrraaa... 2.113803.
Petrillo, G. On the Gaussian process limit of Bayesian Additive Regression Trees, 2024.
Quinlan, J. R. Combining instance-based and model-based learning. In Proceedings of the Tenth International Conference on International Conference on Machine Learning, ICML'93, pp. 236-243, San Francisco, CA. USA, 1993. Morgan Kaufmann Publishers Inc. ISBN 1558603077.
Thebelt, A., Kronqvist, J., Mistry, M., Lee, R. M., Sudermann-Merx, N., and Misener, R. ENTMOOT: A framework for optimization over ensemble tree mod- cls. Computers & Chemical Engineering, 151:107343, August 2021. ISSN 0098-1354.
doi: 10.1016/j.
compchemeng.2021..o 2021.00333.
Quiroga, M., Garay, P. G., Alonso, J. M., Loyola, J. M., and Martin, O. A. Bayesian additive regression trees for probabilistic programming, 2022.
Thebelt, A., Tsay, C., Lee, R. M., Sudermann- Merx, N., Walz, D., Shafei, B., and Misener, R. Tree ensemble kernels for Bayesian optimization with known constraints over mixed-feature spaces. In Proceedings of the 36th In- ternational Conference on Neural Information Processing Systems, NIPS 22, Red Hook, NY, USA, 2022a. Curran Associates Inc. ISBN 9781 171117711888..
Rasmussen, C. F. and Williams, C. K. I. Gaussian Processes for Machine Learning. The MIT Press, 2005.
Ru, B., Alvi, A., Nguyen, V., Osborne, M. A., and Roberts, S. Bayesian optimisation over multiple continuous and categorical inputs. In III, H. D. and Singh, A. (eds.), Proceedings of the 37th International Conference on Ma- chine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 2276-8285. PMLR, July 2020.
Thebelt, A., Tsay, C., Lee, R. M., Sudermann-Mers N., Walz, D., Tranter, T., and Misener, R. Multi-objective constrained optimization for energy applications via tree ensembles. Applied Energy, 306: 18061, January 2022b. ISSN 0306-2619. doi: 10.1016/ 10.10166j.peeneryy2222 apenergy.2021.111661
Scornet, E. Random Forests and Kerel Methods. IEEE Transactions on Information Theory, (2(3)) March 2016. ISSN 1557-9654. doi: 10.1109/ti..016.
Wan, X., Nguyen, V., Ha, H., Ru, B., Lu, C., and Osboree, M. A. Think Global and Act Local: Bayesian Optimi- sation over High--Dienssional Categorical and Mixed
2514489.
12

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Search Spaces. International Conference on Machine Learning (ICML) 38, 2021.
Wang, L., Fonseca, R., and Tian, Y. Learning Search Space Partition for Black-box Optimization using Monte Carlo Tree Search. In NeurIPS, pp. 1511--1922, 2020.
Wang, M., He, J., and Hahn, P. R. Local Gaussian Process Extrapolation for BART Models with Applications to Causal Inference. Journal of Computational and Graph- ical Statistics, 33(2)7224-735, September 2023a. ISSN 1537-2715. doi: 10. 00000/006600000222222222
Wang, W., Zhang, X., and Zou, L. Regret Optimality of GP-UCB, 2023b.
Wang, X., Jin, Y., Schmitt, S., and Olhofer, M. Recent advances in Bayesian optimization. ACM Computing Surveys, 55(13s):1133 2023c.
Wilson, J., Hutter, F., and Deisenroth, M. Maximizing quisition functions for Bayesian optimization. Advances in neural information processing systems, 31, 2018.
Woodbury, M. A. Inverting modified matrices. In Mem- orandum Rept. 42, Statistical Research Group, pp. 4.
Princeton Univ., 1950.
Xie, Y., Zhang, S., Paulson, J., and Tsay, C. Global Op- timization of Gaussian Process Acquisition Functions Using Piecewise-Linear Kernel Approximation, 2024.
Yeh, I-C. Modeling of strength of high- performance concrete using artificial neural net- works. Cement and Concrete Research, 28 (12):1797-1808, 1998. ISSN 0008-8846.
doi:
http:://owoggoooooooooooooso4469 com/ URL
Zhu, C., Byrd, R. H., Lu, P., and Nocedal, J. Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound- constrained optimization. ACM Trans. Math. Softw., 23 (4):550-560, December 1997. ISSN 0098-3500. doi:
10.1145/279933.279966 URL https: http://oo..oo/ /rg/ 10. 10.1145/279232.27999
Ziomek, J. K. and Bou Ammar, H. Are Random Decom- positions all we need in High Dimensional Bayesian Op- timisation? ICML, volume 202, pp. 43347-43366 2023.
13

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Uljoozier Punction
GP Preficions
Otjectan Smation
UP Proficton
-0
0.00
GE-R0P
Otyistre fteetion
BAIK
Orgecther 20
(a) Continuous
(b) Discrete (ordered)
Oljectin Daskin
CP Proicio
RASS
Uijester Fumier
(c) Categorical (unordered)
Figure Regression on set ID toy functions, fit using a BARK GP and an RBF kernel. These functions highlight some notable differences in behavior between the two models.
A. Related work comparison
We provide summary of the existing literature in Appendix A, that highlights key similarities and differences between the various BO methods covered in the literature review.
B. Toy regression examples
As discussed in the main paper, our key motivation behind using tree-based kernels is the optimization of the acquisition function over mixed feature spaces. However, these models also model the underlying function differently, which we investigate in this section by considering some toy ID problems in Fig. 9.
First, the kernel is better able model non-stationary functions, where the lengthscale varies over the search space. For example, in Figs. 9a and 9b the very short lengthscale for small input values means that the RBF kernel learns short lengthscale over the whole space, leading to large uncertainty around x 0.5, unlike the BARK model which has much lower variance in this region.
Second, BARK has a prior belief about the correlation between categorical values. Since each split in the tree structure divides the categorical values into two subsets, the prior probability that two values will be placed in the same leaf non-zero. Conversely, the indicator kernel used in Ru al. (2020) explicitly assigns zero correlation between different categorical values. This results in increased variance for unseen categories, shown in Fig. 9c.
These properties are not necessarily 'better' different black-box functions require different modeling assumptions. For example, the MAX benchmark in Section 7.3 requires fewer BO iterations when modeled using trees. However, other benchmarks are better modeled with the assumption of smoothness provided by Euclidean GP kernels.
C. Notation
C.1. Table of symbols
We list selection of the key symbols used in this paper in Table
This method requires prior knowledge of structure the tree function.
14

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Table 2: Comparison of selection of relevant existing approaches for Bayesian optimization. For each method, describe the underlying surrogate used to model the objective, and the tool used to optimize the acquisition function. A 'mixture' kernel is the combination of a categorical kernel ky and a continuous kernel kg of the form k(z,z') (h, h²) (x, + kb, (h, h') + (x,x')
Trees/ Partitions
Mixed
Method
Reference
Surrogate
Optimization
GP
domain
BART
Chipman et al. (2010)
Bayesian sum-of-trees
Grid
X
ENTMOOT
Thebelt et al. (2021)
Gradient boosted trees
MIP
X
/
SMAC
Lindauer et al.
Random forest
Random local search
X
(2022)
LeafGP
Thebelt et al.
GP (Tree kernel)
MIP
(2022)
BARK
Ours
GP (Bayesian tree kernel)
MIP
Exact GP
Garrido-
GP (Matérn) with input rounding
Gradient descent
x
Merchán (2020)
CASMO- POLITAN
Wan et al. Local GP (Mixture) (2021)
Altemate gradient de- scent + discrete per-
turbation
CoCaBo
Ru et al. (2020)
GP (Mixture)
Multi-armed bandit +
X
gradient descent
RKHS
Em-
Buathong et al. (2020)
GP (Deep embeddings)
Genetic optimization
X
beddings
using derivatives
MiVaBo
Daxberger et al. (2021)
Linear model (with non- linear features)
Thompson sampling
Gryffin
Häse (2021)
al.
BNN with simplex projec- tion
Gradient descent
X
X
HyBO
Deshwal et al. (2021)
GP (Diffusion kernel)
CMA-ES: alternating continuous and dis- crete subspaces
X
PR
Daulton et al. (2022)
GP (Mixture)
Gradient descent
TuRBO
Eriksson et al. (2019)
Local GP (Matérn)
Thompson sampling
X
BAxUS
Papenmeier et al. (2022)
Local GP (Matérn, with em- beddings)
Thompson sampling
X
Tree tured
struc-
Jenatton al. (2017)
GP (Matérn)
shared
Path-based gradient descent
weights
15

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Table 3: Symbols used to define tree operations. Any symbol followed by a star, is defined with respect to proposal
Symbol
Meaning
m
Number of trees in forest
T
Tree structure
M
Vector of leaf values
T
Forest structure (set of trees, (TrH=²)
M
Forest Ieaf values (set of values, {M 1m
{ +
Proposal probability
Wo
Number of leaf nodes Number of singly-internal decision nodes (both children
W1
are leaves)
Node
dg
Depth of node 7
(a, (
Node depth prior parameters
(V,g) ()
Noise prior parameters
C.2. Model hyperparameters
For clarity, we provide an overview of the parameters used in the BARK model Appendix C.2.
Table 4: Hyperparameters and parameters of the BARK model. For fixed parameters, their value is given. Parameters that are sampled are given the label Sampled.
Symbol
Description
Value
Number of trees in ensemble Parameters of node depth prior
50
(a, 3)
(0.95, 2)
(v,
Parameters of noise prior
(3, 0.9)
Signal variance
Noise variance
Sampled Sampled
y
The ith tree structure
D. Experimental details
D.I. Method details
BARK. For both model fitting and BO, 16 samples for the BARK kernel. Increasing the number of samples would (marginally) improve the model fit, however this would increase the complexity of the optimization problem, meaning that the optimization would take longer to converge. BARK uses 1000 burn-in samples and 400 kernel samples, running in parallel with 4 chains. We use a thinning rate of 100 to obtain the 16 samples. We evaluate the effect of increasing the number samples in Appendix G.
In we only generate the burn-in samples in the first iteration. Since each iteration only adds a single datapoint, the posterior from the previous iteration will be close enough that of the current iteration that further mixing is not required.
For BARK (as well BART), use the default number of trees 50, suggested by Kapelner & Bleich (2016), and the default parametrization of the noise, (v, g) (3, 0.9).
BART. For model fitting and BO, BART takes 1000 samples to burn-in the MCMC chain, and 1000 posterior samples. This run in parallel for chains.
There no procedure in the literature for optimizing acquisition functions for BART, and the MIP formulation would be too
16

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
expensive for the 1000s of BART samples required. We therefore evaluate the surrogate grid of max 2 points using space-filling approach (Sobol, 1967). At each iteration, sample new grid to better explore the domain. This grid size was chosen to match the wall clock time of the BARK optimization (as shown in Fig. 11.
Each sample f(5) in BART is single, deterministic function draw, whereas a BARK sample f GP( .) is distribution over functions. Therefore, where BARK allows access to the exact UCB for each sample as detailed in Eq. (1), BART requires the sample-based estimation of Wilson et al. (2018), (C $ - / KV u4, where μ the empirical average of the observations.
We the PyMC-BART v0.8.2 implementation (Quiroga et al., 2022).
GP-RBE. We use the default priors from (Hvarfner et al., 2024), namely squared exponential kernel with log-normal lengthscale prior
p(c) = LV + log vD.
We use the BoFire v0.0.16 implementation (Dürholt et al., 2024), which in turn uses BoTorch v0.1 (Balandat et al., 2020). To optimize mixed space acquisition functions, we use the optimi ze_acqfff_ cqcffm_mm _mm xed aad_allleee alt.rnattatt....... in ng function provided by BoTorch.
LeafGP. We use the implementation provided by Thebelt al. (2022a).
SMAC. We the SMAC3 v.2.0 implementation (Lindauer al., 2022), using the ask-tell interface.
Entmoot. We use the BoFire v0.0.16 implementation (Dürholt et al., 2024).
D.2. Benchmark details
description of the benchmarks used in Bayesian optimization is given in Tables 5 and 6. For each BO benchmark, generate set of min(2D, (2)) initial points, and run the optimization loop for 100 iterations, reporting the minimum value observed to iteration t.
The experiments were run on a High Performance Computing cluster, equipped with AMD EPYC 7742 processors, with each core allocated 16GB of RAM. For models capable of multithreading, we use CPUs; otherwise, only 1.
Table Synthetic benchmarks for Bayesian optimization. For each benchmark, we give the number of dimensions D, and the bounds for the continuous (x), discrete ordinal (i), and categorical (h) features.
Benchmark
D
Features
TreeFunction
10
[0, 1]
TreeFunctionCat
20
x E [0, 1|
h E [1,..., 5
DiscreteAckley
13
x E [-1,1]3
i {-1,
DisereteRosenbroock
x [-5, 10j
ie (-5, 0,5, 1036
E. Additional results
E.1. Regression comparison with GP-RBF
In Section 7.1, we show that BART and BARK have similar modeling abilities, and that our kernel perspective on BART still leads to strong regression performance. For completeness, we also provide a comparison to the GP-RBF model, with linear combination of sum and product kernels, in Table
E.2. Continuous benchmarks
In the main experiments of the paper, we note the motivating use case for BARK is in optimizing over mixed feature spaces. However, we also provide two additional benchmarks here, for continuous synthetic problems. Specifically, we
17

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Table 6: Applied benchmarks for Bayesian optimization. The descriptions of SVRBench, XGBOostMNISST, and PestControl are from the source paper, Dreczkowski et al. (2023).
Benchmark
D
Features
SVRBench
53
21 is log 6, 22
log C, D3 log 7, are features to include.
21 E [-2,
22 E 2, 2
23 E 1
h E (exclude, include) 60
XGBoostMNIST
xi is log learning rate, æ2 min split loss, ca is min split loss, 14 reg lambda, 11 is max depth, h1 is booster, hg grow policy, h3 is objective
21 E 1-5,00 O]
42 E [0, 10]
23 E [0.001, 1]
I4 E [0,
11 E
h1 E (gbtree, dart}
h2 E (depthwise, lossguide]
ha E (multi:softmmxx, multi:softprob)
CCOBench
30 x are transmission powers, are downtilts.
x E 50
E (0.....5) 5515
PestControl
25
h are pesticide choices at each of 25 stages. h E {pestl, pest2, pest3, pest4, none] ]
MAX
28 ho, h1, h2 are the M, A, and X elements respectively. X(I,",,11 123 are chemical properties of the MAX compound. A(13,-..28) are nuisance features (uniform noise).
ho E Melements:
Melements
9
h1 E Aglements Aetements
= 12
h2 E Xelements
Xelements
= 2
(13,-r .28) E [-1,1]
use the Hartmann and Styblinksi-Tang problems from the SFU library (Surjanovic & Bingham, 2013). We show the BO performance in Fig. 10. The smooth functions more suited GP-RBF, however BARK outperforms all other tree-based methods, and still shows good performance in this unfavorable setting on the Styblinski-Tang benchmark.
F. Time comparison
F.1. Posterior sampling
BART and BARK both take significantly longer to train than other baselines, as they are MCMC-based methods. Given N observations, training BART has O(N) complexity, whereas BARK trains with ((N3) complexity (which is reduced to O(N2) when sampling tree structures due the low-rank update discussed in Section 5). We therefore expect BART to scale far more favorably with many observations. In Table 8. show the time taken to train the models on the regression problems from Section 7. The fitting time for BARK is reasonable in the low-data regime of BO, matching that of BART when only 40 training points are used. Since we apply our model in the BO setting, where the evaluation of the black-box function may take hours, or possibly even days, we argue that this time taken is not a significant limitation of the model. We also note that we use fewer samples for BARK (400 samples per parallel chain) compared to BART (1000 samples per chain).
18

---

Hartmann (6D) and BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Table Regression performance of the BART and GP-RBF models, given as mean (and standard deviation). The number of training points used for each benchmark is given by n. The number of dimensions for each problem is given by D (continuous integer categorical). The better performing model is in bold for each metric.
Benchmark
NLPD 4
MSE
BARK
GP-RBF
BARK
GP-RBF
Abalone
400
7+0+1
1.09 (0.06)
1.06 (0.05)
0.52 (0.06)
0.49 (0.05)
Auto MPG
100
443+0
0.45 (0.08)
0.50 (0.22)
0.15 (0.02)
0.15 (0.03)
Concrete Compressive Strength
300
8+0+0
0.31 (0.06)
0.39 (0.04)
0.11 (0.01)
0.15 (0.01)
Student Performance
250
0+13+17
1,27 (0.06)
1.43 (0.06)
0.73 (0.07)
1.00 (0.10)
Hartmann
Styblinsk Tang
200
Best olserved
-20
330
50
Iteration
Iteratinn
LeatGP
GP REF
BARKР
Figure 10: Bayesian optimization performance on benchmarks with a continuous ooly domain: Styblinski-Tang (10D), both from the SFU library of synthetic functions.
Table 8: Mean (and standard deviation) time to generate 1000 samples from the posterior for the BART and BARK models across 20 random data splits. We compare against BART using the PyMC-BART implementation.
Benchmark
Sampling time (seconds per 1000 samples)
BART
BARK
Student Performance
40
20.84 (5.86)
19.36 (7.80)
Auto MPG
100
18.09 (6.08)
71.04 (7.97)
Abalone
400
20.81 (4.98)
951.91 (14.28)
19

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
F.2. Optimization
We provide comparison of the time taken fit and optimize the models in the BO loop in Fig. 11.
BARK is slower than competing methods, taking approximately 50s to fit the model, and 100s to optimize. This is due to the combination of the expensive MCMC procedure, and the large optimization formulation. This is comparable to the fitting time for BART, and the time taken to evaluate BART on grid of 21 points. BARK is best applied in BO settings where the objective is expensive to evaluate (e.g. taking at least several minutes, or having a large associated financial cost). Note that the PestControl and MAX (material design) benchmarks reflect examples of such black-box functions. For settings where experiments are cheap and/or function evaluations are quick, we recommend alternate methods.
G. Hy Hyperparameter sensitivity
G.1. Node depth prior
Throughout our experiments, we use the default BART values of (a, B) = (0.95, for the node depth prior. These are the only hyperparameters for the tree generating process, so it is natural to ask how sensitive our method is to the selection of these values.
We evaluate the regression performance for grid of (o, 9) values in Fig. 12. We show that the default values have generally strong performance across both small- and large-data regimes. Furthermore, changing the hyperparameters do not have a significant impact on the observed NLPD all of the results are well within the standard deviation reported in Table 1. We therefore conclude that our method is largely insensitive to the hyperparameter choice, and use the default BART values.
We also note that, while a smaller value of 3 may lead to stronger regression in some cases, concentrating the prior on deeper trees will make the optimization problem more difficult. Maia et al. (2024) explore a similar sensitivity analysis for their BART-based model, and also opt for the default BART hyperparameters.
G.2. Number of samples
In all experiments with BARK, we use 16 samples from the posterior 4 samples collected from 4 parallel MCMC chains. This is a compromise between model fit, and optimization speed: increasing the number of samples improves the performance of the surrogate, the cost of adding more constraints to the optimization formulation, rendering the task of optimizing the acquisition function more expensive.
In Fig. 13, show the effect of increasing the number of samples regression performance. We observe that increasing the number of samples 16 gives marginal improvement on predictive performance, while rendering the optimization problem significantly more difficult.
H. Posterior sampling
Many of the derivations are included in Kapelner & Bleich (2016). We include of these terms here for completeness, however derivations are omitted where they presented in the literature. A table of symbols is provided in Appendix C.
The grow and prune proposals are the inverse of each other, and the change operation is its own inverse. As long as the probability of selecting each proposal is equal to the probability of its inverse proposal, the ratio of proposal probabilities (e.g., the probability of selecting the grow proposal) does not need to be included in the Metropolis Hastings ratio.
Transition ratio (prune and grow) Below is the ratio for the grow proposal. The probability of selecting the target node is we and the probability of selecting the target node for a prune (the inverse proposal) is 1. The ratio for the prune "["
proposal the inverse of this ratio.
(T* + + T) qT + T*)
Pr(rule) wi
(6)
Prior ratio (prune and grow) expression.
The grow prior ratio is given below, and the prune prior ratio is the inverse of this
20

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
PostControl Optimization
Fitting
Total
Time (s)
10a
50
Number of iterations
TreeFunctions Optimizzation
Fisting
Total
10a
Time (8)
50
100
50
100
50
100
Number iterations
XGBoostMNIST
Ficting
Optimization
Total
10a
Time (x)
10
50
Number of iterations
BARK
BART
BARKProor
GP-RH
(eaum))
Figure 11: Time taken to fit the surrogate, optimize the acquisition, and combined time for each method used in the Bayesian optimization experiments. We compare the time taken per iteration across three benchmarks. For the XGBoosIMNIST benchmark, GB-RBF is optimized by enumerating every discrete combination, leading to increased optimization times in line with BARK.
21

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
AutaMPG (u=20)
AutoMPG (n=100)
0.99
0.99
0.47
0.69
0.95
0.95
0.46
0.68
NLPD
NLPD
0.9
0.45
0.63
0.5
0.5
0.44
1.0
1.5
1.75
2.0
2.25
2.5
1.0
1.5
1.75
2.0
2.25
2.5
/
Student (u=40)
Student (n=250)
0.99
0.99
1.365
1.275
0.95
0.95
1.360
NLPP
1.244
0.9
0.9
1.355
1.273
0.5
1.0
1.75
2.25
1.5
2.0
2.5
B
Figure 12: Regression performance of BARK for four tasks, varying the node depth hyperparameters, averaged across 20 test/train splits. Larger, darker dots show stronger performance. The default value pair (0.95, 2.0) is circled. Note that the range of NLPDs for each task is small, and is far less than the standard deviations of performance given in Table 1.
(1 ((1 + d,) a) Pr(rule)
p(T*) p(T'
[2+d,1'
We note that the splitting rule probability terms in the transition ratio and the prior ratio cancel. This means that, as long as the same distribution is used for the transition kernel as the prior distribution on trees, the acceptance probability independent the choice of splitting rule distribution.
Transition ratio and prior ratio (change) As noted in Kapelner & Bleich (2016), the transition ratio and the prior ratio cancel in the expression for the acceptance probability. All that remains is the likelihood term.
Prior ratio (noise) This can be simply calculated using the ratio of the priors placed over the hyperparameter.
Transition ratio (noise) stated in Section we sample the noise hyperparameter transformed space, to ensure that the parameter remains positive. Specifically:
0 9 - (0-)
0r N0 N00000003 ;0,02) -² = 0(0))
22

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
AutoMPG 1n=100)
Student (n=250)
1.29
Number of chains
6
0.55
Number of chains
5
5
1.28
4
0.50
NILPD
NLPD
3
3
1.27
2
0.45
1.26
5 6 7
8
1
2
3
7
8
Samples per chain
Samples per chain
Figure 13: NLPD for varying the number of chains and number of samples per chain, averaged across 10 runs. Larger, darker dots show stronger performance. The BARK default (4, 4) circled.
where y: R- Rt. We use the softplus transform, g(0) = log(1 + exp(O)). Given this transform, we can then obtain the probability density function o²",
d
faz- (o²²) = =fa- (0=2))))00")
d0 9 (022)
log f-2--00")) (02*) = 20 g (0- -9 - (
log 1 - 1eexp
-exp(-o22) (-o2-)) + const
We therefore obtain the log-transition ratio,
L0 Soa
(1+
exp(o22) - exp(o>2 - + 0 o2
log
log
o²?
I. Regret bounds and kernel approximation
Bayesian optimization literature often provides bounds on the cumulative regret during the optimization process. this section, we diseuss why this is challenge for method, and provide some discussion how regret bound might be attained.
L11 General regret bounds: existing methods are not applicable
An initial exploration of the literature in Bayesian optimization and bandits gave us two promising avenues for proving regret bounds, however upon more detailed examination we found the problem to be far more difficult than we originally envisaged even under strong simplifying assumptions such as the tree structure being known.
The first attempt relates to the direct application of Theorem in Srinivas et al. (2012). Indeed, if we assume the tree structure was known, the kernel is well defined and the bounds would seemingly follow. However, the proof heavily relies on Lipschitz-continiitty the function, and the assumption of an RHKS norm implying differentiabilltty. our setting, this not necessarily the case near decision boundaries, there are discontinuities that cannot be bounded.
An additional avenue that could be explored relates to the bandit literature. Indeed, assume that the objective function, f. is sampled from tree-kernel GP. The complexity of f is bounded: this can be achieved by thresholding the prior probability of the tree structure. Observations of ff have Gaussian noise, i.e. y(x) f(x) + e. Even though we have continuous space, there a finite number of splits, so the resulting problem is trying to find the best subspace defining a multi-armed bandit problem. However, the number of resulting arms will be extremely large, therefore standard bounds for independent bandits under UCB will be far too loose.
23

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
We could further reformulate the problem in semi-bandit setting, where each of the trees in the forest bandit, and instead of observing individual rewards, observe the sum of the bandits in called "super-arms". This setting is explored in the combinatorial multi-armed bandit literature, where there is collection of M super-arms each with arms whose sum we can observe. However, in the literature M << but in our case we have far more splits than trees, so M >> which again results in very loose bounds. This approach would also be conditioned on a known tree structure, which would not match the setting where the tree structure is being explored.
I.2. Asymptotic bounds: kernel approximation fails
1.2.1. THE LAPLACE KERNEL AS AN APPROXIMATION
The connection between asymptotic tree kernels and the Matern- kernel (often referred to as the Laplace kernel) have been explored in the literature since Balog et al. (2016). This opens up the idea that Matern regret bounds (Wang et al., 2023b) can be combined with approximation guarantees (e.g. Xie et al. (2024)) to obtain asymptotic regret bounds (as the number of trees increases). However, as we will show, further inspection reveals the kernel approximation is too crude for the BART and BARK models to be meaningful.
Finstly, let begin by showcasing the scenario in which the kernel is good approximation, as it was originally explored in Linero (2017). We work in the domain X' = [0, 1] and in the following setting. Something that is important to note is that the following two assumptions do not hold for either BART or BARK models:
Assumptions.
1. For any leaf node v) E Lr, the depth of the node d,, ~ Poisson(A)
2. For any decision node 7 E Jr; the selected feature in Categorical ( b), and the decision rule C, u10 1).
Theorem. The limit of the kernel as m - co is given by
A k(x, x') = exp - 1x1xx(11)
(8)
Proof. We aim to find k(x, x'), which is equal to 1 - Pr(split E xx) of no splits in this interval is Pn (1 - 8)". For the leaf node that contains the point x, let dy ~ By assumption, the probability mass function of q is the Poisson distribution. Let x ~ x' denote the event that x and fall in the same leaf, i.e. Pr(x ~ x') = Ex(((x, TTx(( (xx', T)). Then:
For n splits. the probability
20 Pr(x~x))= = Σ JnPn
n=0
20
qn (1 8)"d
A" exp( -1)
(1 8)
n!
R=0
[A(1 - 5)}"
= exp(
X)
Σ
n!
n=0
= exp(
- exp (A(1 8))
- 1 x/11
= =xxp
where we use the Taylor expansion of exp(A(1 - -5) to obtain the approximation.
24

---

BARK; A Fully Bayesian Tree Kernel for Black-box Optimization
The t term in the Laplace approximation points to useful natural property of the BARK kernel: as the dimensionality of the problem increases, the effective lengthscale increases. This leads to similar behaviour Hvarfner et al. (2024), where we assume less complex models in higher dimension to mitigate the curse of dimensionality.
If the BARK kernel can indeed be approximated by the Laplace kernel, then we can use the result of Wang et al. (2023b) to show that the cumulative regret is bounded by Ö 3420) ), enjoying sub-linear regret. Even if the approximation is only correct up some error e, then could employ Xie et al. (2024) correct the regret bounds.
1.2.2. DROPPING SPLIT INDEPENDENCE AND POISSON ASSUMPTIONS
In the previous subsection, we assume that the splits chosen at each decision node are independent. This results in 'logically empty' leaf nodes that cannot be reached by any point in the domain. However, this assumption is necessary to obtain a stationary kernel.
We also note that the assumption of the Poisson distribution of node depths does not follow the BART (nor BARK) prior. The Poisson distribution is obtained by considering the successes of several independent identically distributed events, whereas the probability of split depends on the depth of the node. The Poisson distribution places too high a probability mass on large trees, which is in practice avoided by the BART depth prior. Moreover, the Poisson distribution is placed on the depth of a given node. However, for any distribution of node depths, datapoint is more likely to fall into a more shallow node, as the volume of the subset of the domain defined by that node will be larger. This assumption therefore places an even higher probability on a given datapoint falling into a deep node than is observed in practice.
This discrepancy was considered in Petrillo (2024), who show that the GP limit of BART follows a complicated kernel that can needs to be caleulated recursively. In this subsection, we alter the assumptions from Linero (2017) and show the approximation is too loose to be used for regret bounds. In particular we consider the following:
1. The splitting rules are sampled uniformly in the subspace defined by previous splits (i.e. there are no logically empty splits).
2. We use the true depth prior used by BART and BARK (i.e. we no longer assume the distribution to be Poisson).
We consider the 'chopping' process. The interval [0, 1] is repeatedly chopped, discarding the region to the right of the chop.
The location of each chop is uniformly distributed on the remaining interval: such that Ta [0, 1]. and the length of the remaining interval l = II t4. For a point æ E [0, 1], let the event Sy be the event that Ly E [0, æ] and la-1 & [0, x]- that is, the dth chop disconnects x from 0. This is equivalent the event that at a split of depth d at a decision node, the points 0 and x are placed in different leaf nodes.
Ta ~ u 1]
log ~ Exp(1)
log La
Σ log Ti ~ Gamma (d, 1)
Given the distribution of - log La, we can compute the cumulative distribution function of L as below:
(l. (d- - 1)! J.
1
log la
1d-le-tae dtt
1
(ta( (d - 1)!
(- log da)
25

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Pr(Sa) = Pr(La S and La-1 > x)
Pr(Sa)
Pr(Lo S La-1 > x)
Pr(Lo-1 2 x)
Pr(Ld S x (d-1)ffa- (ld-1 La-1 x) dl
{ x Sta (ld-1)
dld-1
l Pr(La-1 > x)
1
Pr(Sa)
( 2) dl
(d 2) J
(d - 1)! (- -logx)d
T
We now consider the probability of there being any split between 0 and æ. Let (rd) the prior probability of node being split, which depends only on the node depth. Then, the probability of any split ocurring between 0 and x is:
Pr(split E [0, æ]) = (1)((((()) + (2)() +..))
Σ
Pr(Sa) II (d)
1=1
This expression is the probability that the two points are separated, i.e., k(0, z). Unfortunately. from this expression we are able see that the approximation will quickly diverge away from the Laplace kernel we increase the preference for lew-hepth trees, be seen in Figure 14. Particularly, for the depth prior both BART and BARK use /3 = 2.0 resulting in vastly different kernel to the Laplace.
J. Optimizer details
Typical acquisition function optimization performed using gradient based approaches, such as the second-order method L-BFGS. However, functions sampled from the BARK kernel have zero gradient everywhere, and these cannot be used. Instead, we formulate the optimization problem mixed-integer program, optimizing directly over the tree structure using the formulation in Thebelt et al. (2022a). We use Gurobi 11 to optimize the acquisition function.
This approach comes with additional advantages:
We can naturally optimize over mixed spaces, as continuous and categorical features are treated equivalently by considering only the tree structure.
We can include nonlinear constraints.
In some cases, can provide guarantee of global optimality.
J.1. Optimizer hyperparameters
Table 9 contains list of the hyperparameters that have been set for the solver. We note that MIP solvers have hundreds of parameters, and it would be infeasible to tune them all.
We set the MIP gap to 10%, such that the best acquisition function value is guaranteed to be with in 10% of the global optimum. We found that in practice, the found point is generally much closer to the global optimum, however it takes more time to prove this optimality. We also note that it is not necessary to guarantee global optimality, since gradient- based methods do not guarantee global optimality either.
26

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
Cg014(011)
Cg (40, NompotyySpll
Poissont
0.2
1.0
1.0
>-rrePrioria. 3)
0.8
0.5
0.8
0.5
0.1
0.4
8 = 0.5
[3 1.0
B 2.0
Figure 14: Tree kernel approximation: we compare the effect of strengthening the prior to prefer smaller trees. Indeed, we observe that as increases, the kernel and the Laplace kernel become very distinct. We set IIn a). Plots k(0,.0 x:).
Parameter
Description
Value
NonConvex
0f the problem is convex in the con- tinuous variables; otherwise Gap between the upper and lower 10%
MIPGap
bound of the optimal value
Heuristics
Time spent using heuristics to find 20%
feasible points
Time Limit
Maximum time to find the optimum point
100s
Table 9: Some settings of the MIP solver parameters.
27

---

BARK: A Fully Bayesian Tree Kernel for Black-box Optimization
J.2. Optimization formulation
See Thebelt al. (2022a) for further discussion the formulation of the tree kernel optimization.
We encode the tree model below. This formulation is based on Thebelt et al. (2022a), with the only difference that the set T is the set of all trees across all MCMC samples. By using 5 samples from the tree kernel posterior, we increase the number of constraints in the optimization model by a factor of S.
Σ LEL
2,,1 = 1,
AA E T, (9a)
Σ
21,1 -
Σ W(s)
AA E T. ss E splits(t), (9b)
(Eleft(s)
jEC(s)
Σ
2r S
Σ
W(s)
3 A
T, 8A E splits(t), (9c)
(Eright(s)
jEC(+)
K
Σ f=1
Vijj 1,
Vi C, (9d)
M/ $ Vi,++1
A.W [K, 11, (9c) ViE [n] Vj E [Ki], (9f) VL VIET, VI L. (9g)
Vi E {0.1}.
2, 2 0,
The objective to be optimized is the UCB acquisition function. The key difference here is that we now have samples of the mean and standard deviation at each input x, which we sum over to obtain the integrated acquisition function. The equality of the variance relaxed to an inequality (S) to provide a second-order cone constraint.
8 Xib: Xih, X E arg max Σ (x)) + (x)))
(10a)
11)) (x) -x (Kxx000) + og 22(01 -1
y
(10b)
(n) (x) (xxxpt)) kxx (K (Kxx;e0 + t x -1
Kxx00)
(10c)
28


Reviews:
[{"review_id": "4xhSxI2tut", "reviewer": "Reviewer_zyg8", "summary": "The paper proposes a combination of forest kernel GPs and Bayesian tree models (BART) specifically tailored for the use in Bayesian Optimisation (BO). The main idea is to directly optimise the expected acquisition function values over the posterior distribution of the kernel parameters. This is done through posterior samples of those parameters obtained via a Metropolis Hastings Markov Chain Monte Carlo (MCMC) algorithm. The final optimisation of the mean acquisition function can apparently be done efficiently via a mixed integer programming (MIP) approach. Compared to a BART-BO approach where posterior sampling is performed to obtain individual trees, the proposed approach is said to lead to a much improved sample efficiency per MCMC sample.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their thorough review, and for identifying the strengths of our work.\n\n> The choice of UCI datasets\n\nSection 7.1 demonstrates that the BARK model performs similarly in regression to the BART model. BART is typically used in regression settings with tabular datasets, so the purpose using these UCI benchmarks is to be fair to BART. This section demonstrates that the improved performance of BARK in BO is due to the acquisition function optimization, and not any superior modeling capability.\n\nAgreed with the reviewer that we should state the purpose of  Section 7.1 more clearly: the revised paper will clearly signpost our objective of comparing regression tasks where BART is already known to be strong.\n\n> the convexity of the expected acquisition function\n\nThebelt et al. (NeurIPS, 2022) show the convexity of the nonlinear part of the acquisition function: see their development in their Section 4. The acquisition function (AF) itself is not convex in the input space, but only with respect to the mean and variance, which leads to improved performance with mixed-integer solvers.\n\n> for the comparison to BART I am wondering though whether the comparison does not have to be refined... the key claim of an improved efficiency seems to require some quantification\n\nMany BART samples are required to compute the predictive variance, e.g. experiments in Chipman et al. (2010) use at least 1000 function samples. We agree with the reviewer that, in regression, this lower sample efficiency may not be a disadvantage. However, the high number of samples means that it is infeasible to formulate an optimization problem over the BART AF, as the number of variables and constraints required to encode these tens of thousands of trees would be too great. The large number of samples therefore necessitates a grid-based approach to optimizing the AF. Our revised paper will clarify the relevance of the 'sample efficiency' claim in Section 5.2.\n\nFurthermore, the required density of the grid on which BART is evaluated increases exponentially with the problem dimension, and so it is not feasible to allow for a sufficiently dense grid in high dimensional problems. The grid sized used in our experiments was chosen to match the wall clock time of the BARK optimization (see [Rebuttal Figure 5](https://anonymous.4open.science/r/bark-rebuttal-8830/Fig\\_optimization\\_time.pdf), where evaluating BART at 2^14 gridpoints takes the same wall-clock time as the BARK optimization).\n\n> Could it be useful to have a hierarchical presentation of a general form of the model used in this and prior works\n\nThank you for the suggestion - we will add to the paper a summary of the existing literature (see [Rebuttal Table 3](https://anonymous.4open.science/r/bark-rebuttal-8830/Tab\\_literature\\_comparison.pdf)),  that highlights key similarities and differences between the various BO methods covered in the literature review.\n\n> no room to discuss [the MIP approach] in the main paper\n\nThebelt et al. (NeurIPS, 2022) develop a detailed discussion of the MIP. We have extended the optimization model, and will clarify these improvements in Appendix G. The contributions we will clarify in the revised version include: extending the MIP to include multiple tree-structure samples and linking the kernel samples in the (approximate) integrated AF.\n\n> see more compactly the various parameters\n\nThanks for the nice idea, which helps us improve clarity. We will add this summary of the BARK parameters (see [Rebuttal Table 1](https://anonymous.4open.science/r/bark-rebuttal-8830/Tab\\_BARK\\_parameters.pdf)) to the paper.", "claims_and_evidence": "The main claim is that the proposed method tends to outperforms similar method in their blackbox optimisation performance. The claim makes sense conceptually and experimental evidence is provided.", "methods_and_evaluation_criteria": "The optimisation benchmark problems are standard in the field. The choice of UCI datasets, which are not used for optimisation, to investigate the model fit is questionable though, as the main objective of the paper is to explain differences in optimisation performance. Here the paper would benefit from a more stringent formulation of objectives (hypotheses to be tested) of the experimental evaluation.", "theoretical_claims": "There is an interesting theoretical claim about the convexity of the expected acquisition function value at query points in terms of predictive mean and variance. This and other claims might be proven in the appendix, which I did not have time to check.", "experimental_designs_or_analyses": "Overall the experiments cover a reasonable range of test cases and the results look sensible. Specifically for the comparison to BART I am wondering though whether the comparison does not have to be refined. This is because we principally allow for large computation times of BO methods. Hence, the key claim of an improved efficiency seems to require some quantification in terms of wall clock computation time. Having less sample efficiency alone is not necessarily a big disadvantage if one can allow a lot of samples and the complexity per sample can differ.", "relation_to_broader_scientific_literature": "The paper makes an excellent job in surveying the key literature in an integrative and insightful way.", "essential_references_not_discussed": "None to my knowledge.", "other_strengths_and_weaknesses": "The paper makes a laudable approach to compactly yet accurately present the extensive amount of ideas, known and novel, that are required to understand the subject matter. Yet to my perception it is still hard to follow and to comprehensively see all the difference to alternative approaches. Could it be useful to have a hierarchical presentation of a general form of the model used in this and prior works, and from there develop their differences and commonalities?\n\nThe MIP approach for acquisition function optimisation seems to be very interesting in its own right. Unfortunately, there seems to be no room to discuss this in the main paper.", "other_comments_or_suggestions": "See other strength and weaknesses. It would really be useful to see more compactly the various parameters, especially the kernel parameters that are ultimately sampled.", "questions_for_authors": "None.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "SpXGiJr3yd", "reviewer": "Reviewer_PzbN", "summary": "The paper introduces a new tree-based surrogate model for use in Bayesian optimization. It extends of prior work in tree-based regression, particularly BART, to be more suitable for acquisition function optimization and Bayesian optimization. The tree model is fully Bayesian, including MCMC over tree structures. The acquisition function is an integer program. The method performs well on real mixed-variable problems.\n\n## update after rebuttal\nNo change", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We kindly thank the reviewer for their thoughtful review, and for recognizing the potential of our method.\n\n> lack of discussion of wall-time for BO\n\nAs mentioned in Appendix G.1, we limit each MIP optimization to 100 seconds. However, we agree that a more complete comparison is highly relevant. Please see requested results in [Rebuttal Figure 5](https://anonymous.4open.science/r/bark-rebuttal-8830/Fig\\_optimization\\_time.pdf) for a selection of 3 benchmarks. We will also include this figure and the following discussion in the main paper:\n\nBARK is slower than competing methods, taking ~50s to fit the model, and 100s to optimize. This is due to the combination of the expensive MCMC procedure, and the large optimization formulation. This is comparable to the fitting time for BART, and the time taken to evaluate BART on a grid of 2^14 points. BARK is best applied in BO settings where the objective is expensive to evaluate (e.g. taking at least several minutes, or having a large associated financial cost). Note that the PestControl and MAX (material design) benchmarks reflect examples of such black-box functions. For settings where experiments are cheap and/or function evaluations are quick, we recommend alternate methods.\n\n> Fig. 5 shows error bars but doesn't say what they are.\n\nThanks for noting our omission: the revised paper will state explicitly in the figure caption that the bars are the 25th and 75th percentile of regret achieved across the 20 runs.\n\n> this software will depend on Gurobi, which is very expensive for many people\n\nGurobi does provide a free academic license, but we recognize that requiring a licensed software will limit the reach of this work. Gurobi is a very strong optimizer for MIP problems, and open source solvers would take longer to maximize the AF. However, we would still be interested in providing such an interface, especially if there is community interest.", "claims_and_evidence": "The claims of the paper were well supported.", "methods_and_evaluation_criteria": "The method was well motivated and very clearly described. There was a strong set of ablations and results for understanding model performance. The choice of benchmark problems was reasonable.", "theoretical_claims": "No", "experimental_designs_or_analyses": "The empirical evaluation was typical for a BO paper like this.", "relation_to_broader_scientific_literature": "The framing with the broader scientific literature was well written and appropriate.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "I found the ideas in the paper to be an interesting, novel, and likely useful. The one major weakness in the paper is the complete lack of discussion of wall-time for BO. I expect this to be slow due to the combination of MCMC and integer programming (both of which tend to be slow). But I do not know if this is unusably slow or not, as the only wall-times given are for sampling, just for the regression problems and compared only to BART (Table 5 in the supplement). I expect to see total wall time for the full BO loop, compared with wall time for the other BO methods; and for that total wall time to be broken down into MCMC time vs. Gurobi time. I cannot recommend the paper very enthusiastically without that result as based on what I see now, the method may not be practical. It seems very promising though.", "other_comments_or_suggestions": "Fig. 5 shows error bars but doesn't say what they are.\n\nThe paper states that software the model will be released. It seems that this software will depend on Gurobi, which is very expensive for many people. If the authors would like to maximize the impact of their work, it would be worth the effort to introduce an interface for using a free/open source solver.", "questions_for_authors": "How does total BO wall time compare to all of the other baselines methods, and how much of that wall time is spent doing MCMC and how much is spent in Gurobi?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Njx3aaqHZk", "reviewer": "Reviewer_GKPb", "summary": "For Bayesian Optimization, the paper proposes to use tree-based functions for the basis for the Gaussian Process prior. The paper outlines the mechanics of this, including:\n  * Kernel definition\n  * Sampling MCMC method\n  * Acquisition definition\n\nThe paper then shows that the regression capabilities are reasonable, and experiments with BO over standard synthetic and applied blackbox-optimization benchmarks which mostly consist of mixed search spaces (i.e. containing categorical parameters). The conclusion is that BARK is fairly decent, sometimes optimal compared to previous standard GP baselines (LeafGP, SMAC, Entmoot, etc.)", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review, and for recognizing our contribution to mixed-space Bayesian optimization.\n\n> what exactly motivated the authors to consider that using trees would be good for BO\n\nWe agree with the reviewer that we should better motivate BARK to ensure practitioners use BARK appropriately. Our revised paper will include the requested figures and the following discussion.\n\nOur greatest motivation for using trees is for modeling - and specifically optimizing over - mixed feature spaces. Purely continuous spaces are well-modeled by standard Euclidean GP kernels. Purely categorical spaces can be addressed with multi-armed bandit approaches. However, many real-world problems have mixed domains. Trees have strong modeling performance in mixed spaces. Moreover, we can define a MIP to optimize the AF exactly and we can formulate this MIP so that it is effectively solved with off-the-shelf software. The resulting BO procedure does not require approximations during the optimization step.\n\nModeling with trees offers several additional benefits. First, we are able to model non-stationary functions, where the lengthscale changes throughout the domain. We demonstrate this behavior in [Rebuttal Figure 1](https://anonymous.4open.science/r/bark-rebuttal-8830/Fig\\_1D\\_toy\\_continuous.pdf). This is especially desirable in high dimensions, where a locally short lengthscale will not harm BARK's global uncertainty quantification.\n\nSecond, we assume prior correlation between values in categorical features due to the nature of splitting rules in trees.\nThe indicator kernel typically used for categorical features with GPs (Ru et al., 2020) assumes that each value is independent. BARK is able to capture correlations between values, as demonstrated in [Rebuttal Figure 3](https://anonymous.4open.science/r/bark-rebuttal-8830/Fig\\_1D\\_toy\\_categorical.pdf). Trees therefore have a more expressive prior on functions over categorical values.\n\nSome black-box functions are indeed better modeled by trees, as demonstrated in the MAX benchmark where a grid-based AF optimization is used for all methods. However, other benchmarks are better modeled with the assumptions of smoothness provided by Euclidean GP kernels.\n\n> Table 1... doesn't compare against other standard regression baselines\n\nSection 7.1 aims to show that BART and BARK have similar modeling abilities, and that our kernel perspective on BART still leads to strong regression performance. Since our only point is that BART and BARK are similar with respect to regression, we did not compare to a wide array of regression models. However, we agree that including RBF-GP would provide more context and the revised paper will include an extended version of the table: [Rebuttal Table 2](https://anonymous.4open.science/r/bark-rebuttal-8830/Tab\\_Regression\\_with\\_GP\\_RBF.pdf). Additionally, as noted in our reply to Reviewer zyg8, we will signpost the purpose of Table 1 more clearly.\n\n> I suspect that the tree-based kernel does indeed perform well on mixed search spaces, but the authors need to make this motivation more explicit. This can be boosted better with... benchmarking on continuous-only functions\n\nThanks for helping us make our motivation more explicit.\nIndeed, we do not expect BARK to outperform an RBF kernel in continuous-only BO. As explained by the reviewer, these continuous-only BO problems may be smooth, and therefore have a higher prior probability under an RBF kernel. The reviewer is correct that our original submission focuses on mixed spaces where we expect BARK to be more suitable. We also agree that adding continuous-only benchmarks offers a more complete view. We will provide two continuous-only benchmarks in the revised paper: Hartmann (6D) and Styblinski-Tang (10D), which we show in [Rebuttal Figure 4](https://anonymous.4open.science/r/bark-rebuttal-8830/Fig\\_continuous\\_benchmarks.pdf). Despite BARK's non-smooth assumption, we still observe reasonable performance from BARK in this setting.\n\n> how this new BARK-GP regresses on a 1D categorical / discrete function\n\nTo motivate using trees, our revised paper will include an example of regression on both a 1D discrete function and a 1D categorical function: [Rebuttal Figure 2](https://anonymous.4open.science/r/bark-rebuttal-8830/Fig\\_1D\\_toy\\_discrete.pdf) and [Rebuttal Figure 3](https://anonymous.4open.science/r/bark-rebuttal-8830/Fig\\_1D\\_toy\\_categorical.pdf). We will explain in the revised paper that BARK may not necessarily provide better regression for discrete, ordinal features (which are simply continuous features sampled on a grid) compared to Euclidean GPs, but that BARK provides a method of optimizing over such features without approximations.", "claims_and_evidence": "Generally yes, the paper did in fact satisfy its stated contributions in the intro, i.e.\n* Proposing tree-kernel based Bayesian regression method\n* Computationally efficient method for training this method (during the log-likelihood maximization phase) and sampling\n* Show that it performs reasonably on Bayesian optimization benchmarks.\n\nThe only contribution that's missing I feel, is motivating \"why propose trees in the first place at all\"? Please see \"Weaknesses\" section for more details. I suspect that BARK would not do as well on a wide variety of continuous-only space benchmarks, and that's fine - but the authors need to make it more explicit where BARK truly shines (mixed spaces).", "methods_and_evaluation_criteria": "Yes, paper uses standard benchmarking functions for BO, well known in literature Synthetic: (TreeFunction, Discrete-Ackley, Rosenbrocks) and Applied (PestControl, CCOBench, etc.). Its regression datasets (UCI) are also well-known.", "theoretical_claims": "Did not check carefully Section F (theoretical regret bounds).", "experimental_designs_or_analyses": "Yes, paper follows standard Bayesian Optimization evaluation procedures (e.g. running baselines, plotting best observed value over multiple seeded runs) - no red flags here.\n\nTable 1 compares two tree-based regression methods (BART and BARK) but doesn't compare against other standard regression baselines (e.g. Euclidean-GP) - why not? Adding these results establishes much better why tree-based BO works better than euclidean distance-based kernels.", "relation_to_broader_scientific_literature": "Mixed search spaces are in general have been fairly difficult to model with general euclidean-based GP kernels. This is because of having to deal with continuity issues / it's difficult to capture the shape of the objective very well with a smooth model designed originally for continuous spaces. \n\nSo this paper does make a solid contribution to this area of Bayesian Optimization.", "essential_references_not_discussed": "Not that I know of.", "other_strengths_and_weaknesses": "# Weaknesses\nWhile the paper discusses in detail the mechanics of tree-based regression modeling and how it would work in a GP setup, I'm still having a hard time understanding the fundamental conceptual benefits. I.e. what exactly motivated the authors to consider that using trees would be good for Bayesian Optimization (BO)? Is it some combination of:\n  * Tree-based functions are better at modeling categorical spaces (which from experience, classic GP-BO methods struggle with)?\n  * The structure of the tree-based function (i.e. being very piece-wise and not smooth) helps with modeling non-smooth functions better than using smooth functions as the basis (as a regular Euclidean GP would)?\n\nIf this is the case, I would strongly urge the authors to motivate these more, instead of jumping into the mechanics of the BARK method. Perhaps also provide figures on how this new BARK-GP regresses on a 1D categorical / discrete function better than a RBF-GP?\n\nEDIT: My weaknesses have been resolved. The authors have made their contributions over mixed spaces much more explicit.", "other_comments_or_suggestions": "For now, I will propose a weak reject - mainly because of the core weakness as I mentioned. But I am fully willing to upgrade my score once the authors respond to this issue.", "questions_for_authors": "Please see my overall weaknesses section.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The authors propose a Bayesian tree model estimated via MC that can be used as a kernel for a GP.  The AF is optimized via Gurobi.  They demonstrate that this kernel can yield improvements over a few baselines.  All reviewers provided at least an accept rating post-rebuttal. It was an interesting and gratifying paper to read. 

Two major weaknesses were pointed out by the reviewers: one was a question (PzbN) about the walltime (which had been evaluated and shared during the rebuttal; wall time is an order of magnitude slower (even w/ a limited AF optimization budget and the model itself would be prohibitively slow in high dimensions).  The other was a comment re: details of the AF optimization (GKPb), which while not flagged by any reviewers, is likely big issue in the empirical evaluation.

Here is my concern with the paper: In short, GP-RBF likely has much stronger empirical performance than reported if the authors were to use a more appropriate method for optimizing the AF.  GP-RBF has the best performance when the search space is enumerable (cf, XGBoost experiment).  When the cardinality becomes too large or the space is mixed, the authors use continuous relaxation for the discrete variables, which is well-documented to not work well.  A simple approach is to instead do alternating local search (see e.g., [2,4,5]), or alternatively gradient-based optimization via probabilistic reparameterization [3].  I expect this approach to be much more effective in the settings studied as well as higher-dimensional problems (studied in some of the papers below) which—as far as I understand—would be impractical for BARK.  It's also possible that instead of 1-hot encoding the categorical, the use of categorical kernels [1,2] or embeddings [6], could improve performance.

I expect that in the CR, the authors will utilize the alternating local search  instead of using the continuous relaxations and rounding.

The easiest way to accomplish this might be to just swap out BoFire for Ax, which should use the same priors and does alternating local search for mixed search spaces (via BoTorch's `optimize_acqf_mixed_alternating()`).  It wasn't clear from the text how the authors were modeling and optimizing the AF for ordered categorical or integer-valued inputs.  You may wish to double check that this information is being passed into BoFire or Ax as such so that they can be encoded directly into the RBF, rather than one-hot encoding.

As per the author-reviewer discussion, I also expect more wall-time discussion and experiments, including plots or tables that include overall walltimes for BARK and baselines.  Breaking these timings out by fitting and optimization time would be insightful, given reviewer feedback. It would be nice to see more benchmarks problems to illustrate scaling wrt the input dimensionality.  Some of the papers below could provide inspiration for such problems:

- [1] Ru et al. ICML 2020. propose the use of a categorical kernel (IIRC used in the two papers below) 
- [2] Wan et al. ICML 2021 use alternating local search, categorical kernels 
- [3] Daulton et al. NeurIPS 2022 optimize over ensembles of exact designs via gradient-based optimization 
- [4] Papenmeiher et al. NeurIPS 2023. demonstrate pathologies of common mixed-variable benchmark problems, and propose an improved embedding model for the discrete spaces 
- [5] Zhang et al.  Technometrics 2019. propose the LVGP, which uses a 1 or 2D embedding of categorical variables; this is fairly easy to implement in GPyTorch by adding an embedding layer.

Finally, none of the reviewers (or the meta-reviewer) verified the theoretical results in detail, so I encourage the authors to try to obtain some additional feedback on those before wrapping up the CR.
Decision: Accept (poster)