BAST: Bayesian Additive Regression Spanning Trees for Complex Constrained Domain
Authors: Zhao Tang Luo
Huiyan Sang
Bani Mallick
Venue: NeurIPS 2021 Poster (2021)
Keywords: Bayesian nonparametric regression
Constrained domain
Ensemble learning
Manifold
Random spanning trees

Abstract: Nonparametric regression on complex domains has been a challenging task as most existing methods, such as ensemble models based on binary decision trees, are not designed to account for intrinsic geometries and domain boundaries. This article proposes a Bayesian additive regression spanning trees (BAST) model for nonparametric regression on manifolds, with an emphasis on complex constrained domains or irregularly shaped spaces embedded in Euclidean spaces. Our model is built upon a random spanning tree manifold partition model as each weak learner, which is capable of capturing any irregularly shaped spatially contiguous partitions while respecting intrinsic geometries and domain boundary constraints. Utilizing many nice properties of spanning tree structures, we design an efficient Bayesian inference algorithm. Equipped with a soft prediction scheme, BAST is demonstrated to significantly outperform other competing methods in simulation experiments and in an application to the chlorophyll data in Aral Sea, due to its strong local adaptivity to different levels of smoothness. 

BAST: Bayesian Additive Regression Spanning Trees for Complex Constrained Domain
Zhao Tang Luo Department of Statistics Texas A&M University zt ztluo@stat. tamu. edu
Huiyan Sang Department of Statistics Texas A&M University huiyan@stat. tamu. edu
Bani Mallick Department of Statistics
Texas A&M University bmallick@stat. tamu edu
Abstract
Nonparametric regression on complex domains has been challenging task most existing methods, such as ensemble models based on binary decision trees, are not designed to account for intrinsic geometries and domain boundaries. This article proposes a Bayesian additive regression spanning trees (BAST) model for nonparametric regression manifolds, with emphasis complex constrained domains or irregularly shaped spaces embedded Euclidean spaces. Our model is built upon random spanning tree manifold partition model as each weak learner, which is capable of capturing any irregularly shaped spatially contiguous partitions while respecting intrinsic geometries and domain boundary constraints. Utilizing many nice properties of spanning tree structures, we design an efficient Bayesian inference algorithm. Equipped with soft prediction scheme, BAST is demonstrated to significantly outperform other competing methods in simulation experiments and in application the chlorophyll data in Aral Sea, due its strong local adaptivity to different levels of smoothness.
Introduction
Over the past few decades, data collected from complex constrained domains have attracted much attention in machine learning and spatial statistics. Domains with non-trivial geometries, such as irregular boundaries, sharp concavities, and/or interior holes due geographic constraints (e.g., lakes and coasts), impose challenges on statistical modeling, as the Euclidean assumption underpinning many traditional statistical and machine learning methods no longer holds for data with intrinsic geometries.
In this paper, consider nonparametric regression problems with features lying constrained domains or, more generally, compact Riemannian manifolds. To be more specific, we model the response variable Y (s)) E IR location compact Riemannian manifold M
Y(s) = (s) + e(s). c(s) iid N(O,02).
for some unknown funetion M + R and noise variance o², In many applications, the true function f() may not be globally smooth but has discontinuities/abuupp changes across some narrow boundary regions in the domain. For example, housing prices can substantially different in two neighboring communities, and ocean chlorophyll data that are separated by narrow peninsula can exhibit distinct spatial patterns. is of great need develop new methodologies that can both respect intrinsic geometries of the domain and capture complicated local discontinuity patterns in the true function.
There growing literature on nonparametric regression and smoothing for complex domains. Spline 37 have been developed for data on constrained domains,
smoothing methods 33 22
41
43
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

but most of them focus on domains in R². Sangalli et al. [6 generalized the spline methods to constrained regions in I For more general Riemannian manifolds, kemel based smoothing models, including kernel regressions [32 19 and local regressions
[ 12 were developed. Gaussian process (GP) regression is another popular tool for nonparametrie regression problems, and many works focus on developing valid covariance kernels on spheres [see (see[20]] 16 [5] among others]. More recently, a few practical GP models for constrained domains and Riemannian manifolds were studied in the literature (24 31
However, most of the aforementioned approaches assume globally smooth true functions and thus may not fully adapt to the ones with local discontinuities.
Ensemble tree models [5 6 have been widely used in traditional nonparametric regression problems. One prominent example is the Bayesian additive regression trees (BART) model [10], due to its versatility and capability of producing uncertainty measures. However, to our knowledge, ensemble methods have not been used for nonparametrie regression on complex constrained domains, and almost all ensemble tree methods rely on binary decision tree partition models as their ensemble members (weak learners). Nevertheless, binary trees may not be ideal to capture possibly highly irregular partitions complex domains as they can only make splits parallel to Euclidean coordinate axes. For instance, in the U-shape domain shown in Figure [c). a complicated and over-clustered binary treed partition is needed to approximate simple partition with three clusters (marked by different colors). Moreover, the rectangular partitions do not comply to irregular domain constraints, which may cause the called "leakage" problems complex domains. With similar partitioning idea, Menafoglio et al. [28) proposed a Voronoi tessellation based model to account for domain constraints, but their method imposes convexity restrictions partitions. Most recently, spanning treed partition models have been demonstrated as an effective tool for characterizing partitions with but their focus is on traditional two-dimensional Euclidean spaces and linear regression settings.
flexible shapes 23 39 26,
Our contribution in this paper is to propose novel Bayesian additive regression spanning trees (BAST) model for nonparametrie regressions on complex constrained domains with efficient Bayesian inference algorithms. The backbone of BAST is a new random spanning tree (RST) manifold partition model, which replaces binary decision trees in each weak learner. RST is capable of capturing irregularly shaped partitions with small number of spanning tree edge cuts while respecting intrinsic geometries and domain boundary constraints. Equipped with a soft prediction scheme, we show that BAST achieves superior prediction performance
other competing methods various tasks, thanks to its strong local adaptivity to different levels of smoothness.
The rest of the paper proceeds follows. In Section we present the RST partition models on mani- folds and develop a new Bayesian nonparametric regression model with RST ensembles. Section 3 discusses algorithms for Bayesian inference. In Section[4] we illustrate the model performance by simulation experiments and a real chlorophyll data set in Aral Sea. Section|]conclle with some discussions. Additional details Bayesian inference, hyperparameter selection, and sensitivity analysis are provided in Supplementary Materials. The code of BAST is available at
the paper
2 Bayesian Nonparametric Regressions with Additive Spanning Trees
2.1 A Spanning Treed Partition Model on Manifolds
Our novel nonparametric regression model is built upon an ensemble of partitions of observations. In this subsection, we introduce a stochastic partition model for data on compact Riemannian manifold via random spanning trees (RSTs), which will serve as a building block to develop the sum-of- sann-of-spannin--reee ttres model in Section 2.2
Let M be a d-dimensional compact Riemannian manifold that is known a priori, and S {s1 (S,....,S,, } C M be finite set of locations on M where the data are observed. We are interested in partitioning S into several disjoint subsets such that each subset consists of nearby locations where the data are relatively homogeneous that can be modeled separately. For spatial data, desired to impose contiguity constraints on partitions. Below, we introduce the notion of spatially contiguous partitions on a manifold based on a spatial graph whose edges encode the relationship of spatial adjacency or neighborhood. Let G = (S, E) be connected undirected graph with vertex set and edge set & that connects E to its "close neighbors". The construction of spatial graphs on a manifold will be discussed later in this subsection. We say T(S) (SI,... Skk J. where S, C S for
often

---

(a) 34
(B) 31
(c)
Figure 1: (a) A constrained Delaunay triangulation graph on a U-shaped domain. (b) A partition with three clusters obtained by removing the red edges in spanning tree. (c) A binary treed partition nested in the three-cluster partition (a, b).
i. .,.,, k, is a spatially contiguous partition of S with respect to G if U S = for all j * j'. and there exists connected subgraph Gj (S,, E,) of G for each j. We call each Sj a cluster, which consists of locations that are connected to each other and thus is spatially contiguous with respect the spatial graph. Henceforth, when there is risk of confusion, we will refer to spatially contiguous partitions simply as partitions. Figure Taa shows an example of a partition with three clusters.
For constrained domains in I such as the U-shaped domain in Figure spatial graphs can be constructed via constrained Delaunay triangulations [CDTs; 8) Specifically, let Go be a CDT mesh on S SUS where SB is a set locations on the domain boundaries. Then the induced subgraph of Go on S can be chosen desired. See Figure[[p for an example of G constructed via CDT. For a general manifold, motivated by the nice adaptive properties of nonparametrie regressions based on K nearest neighbor (K- (K-NN) graphs on manifolds 27.. one may construct G by a K-NN graph that connects E S to its K nearest neighbours with respect to geodesic distance.
spatial graph G. Edges longer than a certain threshold can be removed if
Given G, we model partitions on manifolds in a similar spirit as the spanning treed partition models developed for two-dimensional Euclidean spaces 23 39 26). Specifically, connected subgraph
T (S, &T) of G is called a spanning tree of G it has no cycle. well-known property of spanning trees is that if a set of & 1 edges in & is removed, we obtain a disconnected subgraph of T with k connected components, which naturally defines a partition I (S) with & clusters by letting S, be the vertex set of the jth component. In this case, say (S) is induced by T. See Figure [[b) for an example. This property implics that we can simplify a complicated graph partition modeling problem to modeling spanning trees well as the number and locations of removed edges.
Mathematically, conditional on and l we assume discrete uniform distribution on all possible partitions induced by T:
p (T()) k. 7) x (T((S is induced by T and has & clusters)
(2)
where (.)) is an indicator function.
Next, we consider a probabilistic model on the spanning tree space. Let We be the weight for an edge E & and w {we E E). We assume an iid uniform distribution the edge weights and let T be the resulting minimum spanning tree (MST), i.e., the spanning tree with minimum Eesey We:
T MST(W), We sid Unif (0, 1),
where MST(W) denotes an MST of G based on the edge weights w. Note that we are not assuming a discrete uniform distribution on the spanning tree space, with which spanning trees for Bayesian inference. As we will show in Section 3.I our model specification leads to an exact and fast spanning tree sampler, taking advantage of the Prim's algorithm for MST constructions.
challenging to sample
Finally, assume truncated Poisson distribution with mean parameter Å on the number of clusters:
Poisson(As) 1(1 $ k: s k),
where & is the pre-specified maximum number of clusters.
3

---

46
Pu tA
Ha
Pu thtss
Hi
++
Figure Demonstration of the partition obtained by adding two spanning treed partitions.
The following proposition states that the support RST is rich enough to accommodate all possible spatially contiguous partitions manifold with more than & clusters. The proof is postponed to Appendix C. Note that similar results do not hold for binary treed partition models. See Figure [Tc) for a counterexample where there does exist a rectangle containing all the blue points without including any green or red ones.
Proposition 1 Let (S) (S 5x) be arbitrary spatially contiguous partition. Then π (()) within the support of the partition model defined by (2) and if k.
2.2 A Sum-of-spanning-tree treess Regression Model
consider the nonparametric regression problem (1- Instead of assuming global continuity, we assume f(-) belongs broad class piecewise smooth functions. We propose model f(-) using summation of weak learners based the flexible RST partitions.
Given the data (Y(s
Given partition r(S) induced by T with & clusters and cluster-wise constants μ (. .p) R* we define mapping from S to R as
g(s T T) 1 e
where we write π (S) for conciseness. This piecewise constant function on S serves as a weak learner for f(-), which approximates f(-) by μu locally in Sy. Bayesian additive spanning trees (BAST) model is a summation of piecewise constant functions based on various spanning-treed partitions. Specifically, for pre-specified M N, BAST models f(-) as
f(s) Σ g(s| T Kims pm)) S,
(5)
where Tm = T. (S) S S } is partition with km clusters induced by a spanning tree Tm of G and 14m (pmml, )..nk..) Although in principle the spatial graphs for each weak learners need not be identical, focus the where they share common G for simplicity.
For S E S, the additive structure ( implics that f(s) equals the summation of the Mmy's 's corre- sponding the clusters from each weak learner that lies in. Figure 2 iillustrates summation of two spanning treed partitions. This together with the shrinkage priors to be discussed in Section 2.3 allows each weak learner to explain small amount of the variation in the response variable. The step function approximation also allows for capturing both smoothness and discontinuities/abuup changes model is more efficient than some existing ensemble binary tree and smoothing methods recovering irregularly shaped regions where discontinuities happen, thanks to the versatility of RST in capturing highly flexible cluster shapes for complex constrained domains.
in f(-) [3).. In particular,
2.3 Prior Regularization
Similar to BART [00, shrinkage priors play an important role in regularizing weak learners and preventing overfitting for BAST. In this subsection we discuss the specification of prior models, which admits the form
M II m==1
p ((Tmy Tmy kims Mm 1=
p(um Tmyy Tm, kim ( To,y Tm, kim) plo?).

---

(4)
Original Partition
(b)
Birth
[c)
Death M
[d)
Hyper
Figure 3: Partitions and spanning trees obtained after (b) a birth, (c) death, or (d) hyper move from the original partition and tree in (a).
We assign an iid RST prior for (Tm: Tm kim) given by 22, (3 and 4 Small values of A and in (4) are typically chosen to restrict the number clusters in each partition, leading to simpler piecewise constant structures in each weak learner which prevent overfitting and encourage better mixing in Markov chain Monte Carlo. Note that we regularize the number of clusters in more direct way than binary treed partition priors [9] I0] which implicitly penalize large numbers of clusters by increasing the prior probability that a node is terminal as the depth of the node.
Rescaling (s) into 0-5.5, 0.5], we opt to place shrinkage prior that concentrates around for Pm following 110. Conditional on (Tm: Tmy km). we choose conjugate prior for Jmm given by
Pm Tm, Tm km ~ Nam (0, (0,
(6)
independently for = 1,..., M, where ou is assumed to depend on the number of weak learners, specifically, 0μ (.5/(avM) with 0. Intuitively, when have larger number of weak learners, is desired to impose stronger shrinkage effects by choosing a smaller a₂ such that each learner is not too influential the overall fit. A default choice of can be 2, which assigns 0.95 prior probability for (s) that lies in 1-0.5.0.3) 0.5].
Finally, assign conjugate inverse-x2 prior on 02: 0² 02s We fix 3 and choose A, in a data-driven way such that the prior satisfies P(o² < x²²) = 0.90 similarly as in BART [10], where G² is the sample variance of Y
...., (s,)).
3 Bayesian Inference
3.1 Estimation
Inference of BAST is based on a tailored backfitting Markov chain Monte Carlo (MCMC) algorithm [8] in which we successively sample (x,, Ti, ki, (.),..... ( TM, (MM, and の² from their respective full conditionals. Our sampler for the full conditional (Tm; Tm, kim |-) where stands for all other parameters and the data Y, consists of two successive steps: (i) we first analytically integrate am out and sample from the collapsed conditional distribution for the partition parameters p(Tm, hi lm) and (ii) we then sample amm. from p(um Am. Tm, ki -). This design leads to better mixing and convergence performance of the sampler by avoiding the trans-dimensional problem for um 's. Thanks to the conjugate priors, sampling from pmu Tm Tmx Km, -) and p(o²l-) follows standard procedures, and we leave the details Appendix A.I.
Below, we focus on the sampling of the RST partition parameters. To draw samples of (Tm Ti km). one of the four moves - birth, death, change, and hyper - is performed with probabilities r& (kim). ra(kim), re (km ), and rh (k,,), respectively (26). The first three moves modify the partition by proposing a new partition induced by the current spanning tree, and the hyper move updates T, by sampling from its full conditional via an efficient sampling algorithm. Each move is detailed below
and some of them are visualized in Figure 3
In birth move, one of the clusters is split into two by randomly removing an edge in Tm that connects vertices belonging to the same cluster. Denoting the new partition by The the Metsopolls-Hastings (M-H) acceptance ratio is given by
ralkim + (b(km)
L C( ( Tms Kim + 1.-))
min
(km + +1 x
L (Y Tmy Tms kims -)
5

---

where L [(Y[mm Tm: km. -) is the integrated likelihood with pm marginalized out, whose closed form can be found in Appendix A.1. Opposite to the birth move, a death move randomly merges two adjacent clusters in Tm- Specifically, an edge in Tm that connects two distinct clusters in wm uniformly selected and then the two clusters are merged into one. The M-H ratio is analogous to (7) In a change move, death move is performed followed by birth move, so that the number of clusters is unchanged. This move is designed to encourage better mixing of the sampler.
Finally, a hyper move updates Tm using an exact sampler, which adaptively learns a spanning treed spatial order so that we can obtain a partition that is more compatible to the homogeneity pattern of data in subsequent MCMC iterations. Specifically, we sample the edge weight we of G from iid Unif 1/2, 1) if two endpoints of e are in different clusters under Ty. and otherwise from iid (0, 1/2). A new spanning tree is the MST generated by Prim's 's algorithm using the new edge weights. can be shown that the resulting MST induces the current partition 39 and an exact sample from its full conditional distribution 226
The overall computational complexity per MCMC iteration is
log n)), where Th is the probability that hyper step is selected which takes On log n1) using Prim's algorithm for CDT and KNNN graphs, and (n) is the computation complexity required when birth/death/change steps selected because closed form marginal likelihood without matrix inversion is available when calculating acceptance ratios. In practice, we suggest small value of Th such as 0.1 to reduce the computation and allow the algorithm to spend more iterations learning good partition compatible with the current tree. To further reduce computation complexity, we have done some preliminary exploration of using different but fixed spanning trees for each weak learner during MCMC (i.e., setting The = 0). As shown in Appendix B.1.3, this significantly speeds up the computation while the prediction performance remains comparable.
3.2 Prediction
The prediction at an unobserved location u E S involves two steps. First, in each weak learner, we randomly assign u to one of its nearby clusters subject to the manifold constraints, using a soft prediction scheme in a similar spirit to Linero and Yang [25] Second, the prediction is obtained by summing the constants corresponding to the clusters that u belongs to over all weak learners.
Specifically, to obtain cluster memberships for u, we define its neighbor set Nu C S as follows. For a constrained domain in IR .u is chosen as the vertices of the triangle containing u in the CDT mesh that belong to S (i.e., vertices on the domain boundary are excluded; see Appendix A.2 for detailed discussions). For general manifolds in higher dimensional spaces, Nu specified as the set of K nearest neighbors of u in S with respect to the geodesic distance.
Let 2m (v) (1...., km. be the cluster membership of generic location v E MM from the mth weak learner such that zm (s) j if s E St" Intuitively, u is expected to share the same cluster membership as one of its neighbors in S. and if is near the boundary of cluster in partition, it is more ideal to assign zu) probabilistically to reflect the partitioning uncertainty and adapt for smoother functions. This motivates us consider the following random assignment for 2m (u) is: given Nu and posterior sample of the partitions, 2m (u) is sampled independently such that IP {zm (u) = 2m (Nu.t)} ag, for l I,..., |Nul, where Nu.e the éth element in Nu and ag /alll satisfies L==1 ag- 1. One can specify as by setting ag 1/INN, for all E, or via inverse geodesic distance weighting such as positive power.
1/d2 (u, Nu,e), where dy (,, .) is the geodesic distance and is some
At the second step, we sum (m,s..( over m = 1,..., M to obtain a posterior predictive value of E(Y given samples of 2m (u)'s. A point predictor for (u) can then be taken as the mean of the posterior draws, which allows us average models with different RST partition structures.
Finally, we remark that the prediction algorithm highly parallelizable, as the predictive sampling for each RST partition is independent.

---

Table Prediction performance of BAST and its competing methods the U-shape domain example. Standard errors are given in parentheses.
BAST
BART
SFS
inGP
MSPE
0.189 (0.001)
1.541 (0.075)
0.418 (0.001)
0.814 (0.002) 0.610 (0.001)
0.1
MAPE
0.188 (0.001)
0.436 (0.010)
0.340 (0.001)
Mean CRPS
0.142 (0.001)
0.380 (0.009)
MSPE
0.464 (0.006)
1.704 (0.053)
0.680 (0.007)
1.057 (0.010)
0.5
МАPE
0.491 (0.004)
0.686 (0.008)
0.591 (0.004)
0.752 (0.004)
Mean CRPS
0.371 (0.003)
0,575 (0.007)
MSPE
1.283 (0.018)
2.650 (0.056)
1.491 (0.020)
1.823 (0.025)
MAPE
0.888 (0.007)
1,072 (0.008)
0.951 (0.007)
1.051 (0.008)
Mean CRPS
0.693 (0.006)
0.889 (0.008)
(8)
Truef
(b)
APE of BAST
(c)
APE of BART
APE oFSF
(e)
APE if
0
Figure 4: (a) True function in the training data. (b-e) APES of one test data set with 0.1. The
discontinuity boundaries are marked as red circles. Black squares in (c) indicate APE >> 4.10.
Experiments
4.1
U-shape Example
We first examine the BAST's performance of recovering piecewise smooth functions via some simulation experiments on a rotated U-shaped domain shown in Figure 4a). Our true function f()) is constructed based on the one in Ramsay [B3] denoted fa(-). We create discontinuities along a circle of radius 0.9 centered at the origin as follows. For locations inside the circle, we flip I'R by setting J -s-- For locations outside the circle, we f 2fn R for those in the lower arm of the domain, and f Jn for those in the upper arm, such that the jump in the lower arm has a larger magnitude. We uniformly generate out-of-sample locations for prediction. Figure 4[a) shows the true function in the training data. The responses are generated according level is replicated for 50 times.
500 locations in the domain as training data and 200
with different levels of noise o 0.1, 0.5, 1, and each noise
The spatial graph is constructed via CDT. We M 20 weak learners and set d 4 and = 10 to restrict the size of each partition. The prediction based on the CDT graph, and we use inverse distance weighting with moves are set as гb = Ta Tr = 0.3 and TΛ = 0.1, with adjustments for cases where km 1 or k. We compare BAST with BART [TO) and two other nonparametric regression methods for constrained domains, the soap film smoothing [SFS:43) and the sparse intrinsic Gaussian process (inGP) regression [II] For BART, we use the same number of weak learners as in BAST and its default settings. We run the MCMC for both BAST and BART for 20, 000 iterations, discarding the first half and retaining samples every 5 iterations. Our MCMC diagnostics suggest no convergence issue of BAST. We specify 32 equally spaced knots for SFS and set its basis dimension as 40. For sparse inGP, we use 24 equally spaced knots and simulate Brownian motions for 100, 000 times. The prediction performance over 200 out-of- testing locations is assessed by mean squared prediction errors (MSPEs) and mean absolute prediction errors (MAPEs). For the two Bayesian approaches BAST and BART, we also compare their probabilistic prediction performance gauged by continuous ranked probability scores (CRPSs) based on their posterior predictive distributions [see, c.g.J4]. For all the metrics, lower values indicate better performance.
to sample cluster memberships. The probabilities for MCMC
7

---

())
(b)
(c)
(d)
True
BAST
InGP
10
10
y
-
-10
11
-10
10
-10
10
-10
-5
10
-10
-
10
Figure 5: (a) True function and training locations (marked as black dots), (b-d) Predictive surfaces of BAST and its competing methods. All plots are viewed along the negative direction of the z-axis.
Table summarizess the average performance metrics over 50 replicates for each noise level. BAST outperforms all its competitors in terms of all the metrics under all the noise levels. This is because BAST partitions the training data in a way that adapts to both the domain constraints and the irregularly shaped discontinuity boundaries, thanks to the flexible RST partitions. In contrast, the binary treed partitions adapt to neither types of boundaries, and neither of SPS and inGP captures discontinuities in the true function. This is also evidenced by Figure&bbe)) bb-), where absolute prediction errors (APEs) of the test data in one replicate with o = 0.1 are shown. The APEs from BAST are small for most locations except for those near the discontinuity boundaries. SFS also has similar patterns; however, errors near the discontinuity boundaries are much higher due to the global smoothness assumption in SFS. The general APE pattern for inGP is similar SFS, except that it has larger errors, possibly due to the low-rank approximation of covariance functions. With the same number of weak learners, BART has larger errors near the upper discontinuity boundary, more rectangular partitions needed to well approximate irregular boundaries. Moreover, it also gives poor prediction at some locations near the domain boundary between the two arms, probably because binary treed partitions do not take into account the domain boundary when making axis parallel splits, and hence force some boundary locations in one arm share the same cluster memberships with those locations in the other arm. We have also experimented using more weak learners in BART, and the results in Appendix 1.1.1 suggest that BAST with fewer number of weak learners still outperforms BART. Computation time for each method is reported and compared in Appendix 1.1.3.
Hyperparameters of BAST can be tuned using standard vross-validation techniques. Our results in Appendix 1.2 show that the fine-tuned BAST with respect to M, le and the shrinkage parameter a for UN achieves better performance than the default version in Table but the performance of them is close to each other. We have also conducted additional sensitivity analyses to the hyperparameters M, k, and A in Appendix B.1.2 which suggests that the performance of BAST is in general robust to them.
4.2 Bitten Torus Example
To illustrate BAST for more general manifolds, consider a bitten torus example similar to Niu et al. B1]. A torus is a two-dimensional manifold embedded in R33 that is parameterized by (8, Φ). where 8 is the angle for the torus and é is the angle for the tube. Let R be the fixed distance from the center of the tube to the center of the torus, and r be the fixed radius of the tube. The Cartesian coordinate (Ir, y, on torus can be written as x (R r cos 0) cos c. y (R + r cos 0) sin and = = rin O. We create a bitten torus by setting Φ E [r/6,, 1.7x] and 0 E [0,2mT]
We consider piecewise smooth true function f(x, y, 2) defined the bitten torus. We divide the torus into three subregions corresponding to E [T/6, 3t/ /4) E (3T/4 /4), and 0 e (5n /4, 1.7mr. respectively. The true functions in the first and the third regions are the same as the one used in Niu et al. 31), while we set the one in the second region as the negative of the function 3T/4 /4 and 5 We generate responses at n = 500 random locations according to (I) with o = 0.1 as training data. The true function and the training locations are shown in Figure 5
in Niu et al. [31], such that there jumps along
We construct the spatial graph by a 10-NN graph based on the geodesic distance. Since the geodesic distance of a torus has no analytic form, we approximate it as in Isomap algorithm [40]. Specifically,
8

---

Table 2: Prediction performance of BAST and its competing methods in the bitten torus example. Standard errors are given in parentheses.
BAST
BART
inGP 2.283 (0.005)
MSPE
0.487 (0.002)
1.115 (0.041)
МАРE
0.307 (0.001)
0.406 (0.009)
1.159 (0.003)
Mean CRPS
0.225 (0.002)
0.355 (0.008)
we first construct weighted, Euclidean distance based nearest neighbor graph some fine grids and the training locations. Then we approximate the geodesic distance between two training locations by the length of the shortest path between them in the graph. For prediction at an unobserved location u, we use its 5 nearest neighbors in S based on the geodesic distance as its neighbor set Nu. We compare BAST with BART that uses Cartesian coordinates as features and sparse inGP that uses 24 equally spaced knots, as SFS is only applicable for domains in R Other model specifications are the same as those in Section 4.1
As in the previous experiment, we compare the prediction performance at 200 random out-of-sample locations. The average performance metrics over 50 replicates in Table Table suggest that BAST provides the most accurate prediction. We further compare the predictive surfaces of the three methods based on 2, 500 grid points in Figure(5[b-d). As expected, both BAST and BART capture the piecewise structure in the true function, whereas inGP does not. In the interior of each subregion, BAST and BART approximate the true smooth functions fairly accurately. The major difference between the two methods occurs near the discontinuity boundaries. In the surface from BART, some partition boundaries are parallel to the Euclidean coordinate axes such as those near x = 5 and y = 2 due to the use of binary trees, while this pattern does not appear in the one from BAST. Overall, the blue subregion recovered by BAST is more consistent with the truth. Notice that the estimated function at the discontinuity boundaries from BAST is smoother thanks to its soft prediction scheme. We have also experimented with different noise levels in Appendix B.2, and the findings are consistent.
4.3 Application to Chlorophyll Data
We apply BAST to analyze average remote sensed chlorophyll data in the Aral data over 1998-2002, which are available in the R package gamair [42). The chlorophyll measurements at 485 equally spaced locations are shown in Figure a(.. The southern part of the domain is separated by the isthmus of the peninsula near E9'E, and both shores of the peninsula have substantially different chlorophyll levels. It is thus desired take into account this geographical constraints when modeling the data. The goal of our analysis is to assess how well BAST captures the patterns of chlorophyll and predicts for unobserved locations in this complex spatial domain.
We follow Niu et al. [3 to rescale the domain and model the chlorophyll level as a function of the scaled longitude and latitude plus some Gaussian noise. We use same setup for BAST as in Section 4.1 except that we set M 30 and i to encourage smaller sizes of partitions as the number of weak learners is increased. For prediction, we sample the cluster membership of an out-oe-sample location u using equal sampling probabilities ag = 1/|Nul For BART, SFS, and sparse inGP, we also use the same settings as in the simulation experiments but with 30 trees for BART and 42 equally spaced knots for both SFS and inGP. The MCMC algorithms for BAST and BART
for 30, 000 iterations, keeping samples every 5 iterations from the second half.
We first compare the prediction performance of all the models via 10-fold cross-validation. In the end, the chlorophyll level at each observed location is predicted exactly once, and we compare MSPEs and MAPES based on all locations in the data set. Similarly, the mean CRPSs over all locations are compared between BAST and BART. Table 3 shows prediction performance metries for four models. BAST achieves the best performance among all the models.
Next, we turn to the predictive surfaces from each model, which are shown in Figure 6 b-e). All models capture the general patterns the data. The predictive surfaces from SFS and inGP are fairly smooth, while some sharp jumps can be observed in the one from BART. The surface from BAST somewhat in between, and preserves small-scale spatial dependence of the data, At the southern part of the eastern basin, BART identifies rectangular region with high chlorophyll level, while the

---

Table 3: Prediction performance BAST and its competing methods for the chlorophyll data.
BAST
BART
SFS
inGP
MSPE
2.346
2.933
2.894
3.191
MAPE
0.905
1.172
1.071
1,200
Mean CRPS
0.633
0.955
Deo
BAST
BART
SFS
InGP
10
10
10
Lottide
45.3
Scaled Let
Scaled Lat.
Scaled Lat
10
46.5
46.0
50
59
60
-15-10-5 0. Scaled Long.
-15-10-5 0. Scaled Long.
10
-15-10-5
15-10 -5 o. Scaled Long.
Longitude
Scaled Long.
Figure 6: (a) Observed chlorophyll data. (b-e) Predictive surfaces from BAST and its competitors.
corresponding region obtained from BAST has irregular shape, which is more consistent with the data and the results from SFS and inGP. This is due to the highly flexible RST partition model that can give irregularly shaped clusters.
Conclusion and Discussion
In this paper, we developed novel Bayesian nonparametric regression model known manifolds and complex constrained domains using additive RST partitions. The RST weak learner enjoys flexibly shaped partitions while respecting the intrinsic geometries and domain constraints. The additive piecewise constant structure further allows BAST to approximate piecewise smooth functions evidenced by our simulation studies and real data
with irregular boundaries of discontinuities. analysis. In the case where the manifold is unknown, one may estimate the geodesic distance from the data [see, e.g., 229 and references therein] to construct spatial graphs. We leave this scenario for future research.
Similar its binary treed counterpart BART, BAST is promising many other Bayesian hierarchical modeling settings, such as classification models with binary and multinomial responses [TO] 33), survival analysis 2 38), causal inference [17, and varying coefficient regressions [1] As BAST is built upon a spatial graph, it is an interesting direction to extend our methodology for classification and regression on general graphs and networks [e.g..B]. Finally, theoretical justifications are important but usually challenging for ensemble methods. For example, theoretical studies of BART have begun emerging only very recently [35] 34 Posterior concentration results of BAST for estimating the true function can be potentially established in similar manners as BART, non-trivial extensions are required to theoretically handle the complex spanning tree partition on manifolds and hence beyond the scope of this work.
prior models in
This work has foreseeable negative societal impacts, but users of BAST should be cautious of prediction outcomes based on ethically biased input data.
Acknowledgments and Disclosure of Funding
The research was partially supported by NSF DMS-1854655 NSF CCF-1934900 and NIH RO1A6000010. The authors thank the referees, the area chair, and Changwoo Lee from Texas A&M M University for valuable comments.
10

---

References
[1] Aswani, A., Bickel, and Tomlin, C. (2011). Regression manifolds: Estimation of the exterior derivative. The Annals of Statistics, 39(1):48-8) 811.
[2] Bonato, V., Baladandayuthapani, V., Broom, B. M., Sulman, E. P., Aldape, K. D., and Do, K.-A. (201 1). Bayesian ensemble methods for survival prediction in gene expression data. Bioinformatics, 27(3):359-367.
[3] Borovitskiy, V., Azangulov, I,, Teremin, A., Mostowsky, P., Deisenroth, M. P., and Durrande, N. (2021).
Matérn Gaussian processes on graphs. In International Conference on Artificial Intelligence and Statisties. PMLR.
[4] Borovitskiy, V. Teremin, A., Mostowsky, P., and Deisenroth, M. (2020). Matém Gaussian processes Riemannian manifolds. In Proceedings of Advances in Neural Information Processing Systems (NeurIPS).
[5] Breiman, L. (2001), Random forests. Machine learning, 45(1):5-32.
[6] Chen, T. and Guestrin, C. (2016). Xgboost: scalable tree boosting system. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 785-794.
[7] Cheng, M.-Y. and Wu, H.-T. (2013). Local linear regression on manifolds and its geometric interpretation. Journal of the American Statistical Association, 108(504):1421-444)
[8] Chew, L. P. (1989). Constrained Delaunay triangulations. Algorithmica, 4(1):97-108.
[9] Chipman, H. A., George, E. I., and McCulloch, R. E. (1998). Bayesian CART model search. Journal of the American Statistical Association, 93(443):955-948
[10] Chipman,
George, E. and McCulloch, E. (2010). BART: Bayesian additive regression trees. The Annals of Applied Statistics, 4(1 2:266-298.
[11] Deshpande, S. K., Bai, R., Balocchi, C., Starling, E., and Weiss, J. (2020). VCBART: Bayesian trees for varying coefficients. arXiv preprint 26X6v:2003.0641.
[12] Di Marzio, M., Panzera, A., and Taylor, C. C. (2014). Nonparametric regression for spherical data. Journal of the American Staristical Association, 109(566)748-765
[13] Dunson, D. B., Wu, H.-T., and Wu, N. (2020). Diffusion based Gaussian processes restricted domains. arxiv preprint 1xXi:20000.0222.
[14] Gneiting. and Raftery, A. (2007). Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):59-377,
[15] Guella, J. C., Menegatto, A., and Porcu, E. (2018). Strictly positive definite multivariate covariance functions on spheres, Journal of Multivariate Analysis, (66:150-155)
[16] Guinness, J. and Fuentes. M. (2016). Isotropic covariance functions on spheres: Some properties and modeling considerations. Journal of Multivariate Analysis, 143:143-152.
[17] Hahn, P. R., Murray, J. S., and Carvalho, C. M. (2020), Bayesian regression tree models for causal inference: Regularization, confounding, and heterogencous effects (with discussion). Bayesian Analysis, 15(3):965-10) 1056.
[18] Hastie, T. and Tibshirani, R. (2000). Bayesian backfitting (with comments and a rejoinder by the authors). Statistical Science, 15(3):196-22)
[19] Henry, G. and Rodriguez, D. (2009). Robust nonparametric regression on Riemannian manifolds. Journal of Nonparametric Statistics, (1(5)6111622)
[20] Jeong. J. and Jun, M. (2015). A class of Matém laké covariance functions for smooth processes sphere. Spatial Statistics, 11:1-18.
[21] Kpotufe, S. (2011). k-NN regression adapts to local intrinsic dimension. Shawe-Taylor, J., Zemel, R.. Bartlett, P., Pereira, F., and Weinberger, K. Q., editors, Advances in Neural Information Processing Systems, volume 24, Curran Associates, Inc.
[22] Lai, M.J. and Schumaker, L. L. (2007). Spline functions on triangulations, volume 110. Cambridge University Press.

---

[23] Li, F. and Sang, H. (2019). Spatial homogeneity pursuit of regression coefficients for large datasets. Journal of the American Statisticul Association, 1146577:1050-1062
[24] Lin, L., Mu, N., Cheung, P., and Dunson, D. (2019). Extrinsic Gaussian processes for regression and classification manifolds. Bayesian Analysis, (4(3):887-95)
[25] Linero, A. R. and Yang. Y. (2018). Bayesian regression tree ensembles that adapt to smoothness and sparsity. Journal of the Royal Statisticul Society: Series B (Statistical Methodology), 8065):1887-1100.
[26] Luo, Z. T., Sang, H., and Mallick, B. (2021). A Bayesian contiguous partitioning method for learning clustered latent variables, Journal of Machine Learning Research, 22(37):1-53 1-2..
[27] Madrid Padilla, O. H., Sharpnack, J., Chen, Y., and Witten, D. M. (2020). Adaptive nonparametric regression with the k-nearest neighbour fused lasso. Biometrika, 107(2)293-310.
[28] Menafoglio, A., Gaetani, G., and Secchi, P. (2018). Random domain decompositions for object-oriented kriging over complex domains. Stochastic environmental research and risk assessment, 32(12):3421 (3333)3321-3337.
[29] Meng, D., Leung, Y., Xu, Z., Fung, T., and Zhang, Q. (2008). Improving geodesic distance estimation based on locally linear assumption. Pattern Recognition Letters, 29(7):862-870.
[30] Murray, J. S. (2020). lon-linear Bayesian additive regression trees for multinomial logistic and count regression models. Journal of the American Statistical Association, (just- costtaccepted):11233
[31] Niu, M., Cheung, P., Lin, L.., Dai, 7., Lawrence, N., and Dunson, D. (2019). Intrinsic Gaussian processes complex constrained domains. Journal of the Royal Statistical Society: Series (Statistical Methodology). 81 63).606776
[32] Pelletier. B. (2006). Non-parametric regression estimation closed Riemannian manifolds. Journal of Nonparametric Statistics, 18(1):57-67.
[33] Ramsay, T. (2002). Spline smoothing over difficult regions. Journal of the Royal Statistical Society: Series B (Statistical Methodology,, 64(2)3077-319
[34] Ročková, V. and Saha, E. (2019). On theory for BART. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2839-2848. PMLR.
[35] Ročková, and van der Pas, S. (2020). Posterior concentration for Bayesian regression trees and forests. Annals of Statistics, (1((4)2208-21333
[36] Sangalli, L. M., Ramsay, J. O., and Ramsay, T. O. (2013). Spatial spline regression models. Journal of the Royal Statistical Society: Series B: Statistical Methodology, pages 681-703.
[37] Scott- Hayward, L. A. S., MacKenzie, M. L., Donovan, C. R., Walker, C., and Ashe, E. (2014), Complex region spatial smoother (CReSS). Journal of Computational and Graphical Statistics, 23(2):340-33) 3600.
[38] Sivaganesan, S., Müller, P., Huang, (2017). Subgroup finding via Bayesian additive regression trees. Statistics in medicine, 36( (5):23911124333
[39] Teixeira, L. V., Assunção, M., and Loschi, H. (2019). Bayesian space-time partitioning by sampling and pruning spanning trees. Journal of Machine Learning Research, 20:85-1. -.
[40] Tenenbaum, J. B., De Silva, V., and Langford, J. C. (2000). A global geometric framework for nonlinear dimensionality reduction. science, 290(5500)):211--2223
[41] Wang. H. and Ranalli, M. G. (2007). Low-rank smoothing splines complicated domains. Biometrics. 63(1):209-217,
[42] Wood, S. N. (2006). Generalized Additive Models: An Introduction with R. Chapman & Hall/CRC, Boca Raton, Florida, first edition. ISBN 1-58488-474-6.
[43] Wood, S. N., Bravington, M. V. and Hedley, S. L. (2008). Soap film smoothing. Journal of the Royal Statistical Society: Series B (Statistical Methodologyj. 70(5):931-955.
12

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] See Section 5
(c) Did you discuss any potential negative societal impacts of your work? [Yes] See
Section[5]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes]
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or URL)? [Yes] Code will be made publicly available upon the publication the paper.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes]
(c) Did you report error bars (e-g., with respect the random seed after running experi- ments multiple times)? [Yes] Standard errors are provided.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Appendix B.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you eite the creators? [Yes] See Appendix B.
(b) Did you mention the license of the assets? [Yes] See Appendix B.
(c) Did you include any new assets either in the supplemental material or as URL? [No]
(d) Did you discuss whether and how consent was obtained from people whose data you're 're using/curating? [N/A]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "jpOrqDU6Kds", "reviewer": "Reviewer_Kv2o", "summary": "This paper proposes a modification of Bayesian Additive Regression Trees tailored to complicated or constrained input space.\nThe basic idea, like BART, is to learn a collection of piecewise-constant weak learners to approximate an unknown regression function. Unlike BART, which parametrizes each weak learner using a binary regression trees with axis-aligned splits, the proposed BAST procedure parametrizes the weak learners using a spatially-contiguous partition of the input space that is index by a random spanning tree.\nThe proposed procedure works quite well on simulated and real world data. It represents a very interesting new direction in a growing literature around BART.", "questions": "", "limitations": "", "rating": 8, "confidence": 5, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality: The individual components of BAST are not especially novel; there are very active communities studying both BART and its extensions and spatially contiguous partitions using spanning trees . The main novelty, then, of the paper is the effective combination of these two ideas.\n\nQuality: On the whole, the ideas underlying BAST are sound and I commend the authors for a well-executed paper. A few minor questions about the method:\n\n- Sensitivity to $\\overline{k}$ & its interplay with $M$: Table S2 in the Supplementary Materials suggests that, at least for the two synthetic data experiments, BAST is not especially sensitive to the choice of the upper limit on number of clusters per partition. It would be helpful, however, to probe the sensitivity a bit more. E.g. how does the performance change as one increases or decreases $\\overline{k}$? I suspect that if one were to increase $M$, one could get away with a smaller choice of $\\overline{k}$ (and vice versa). A more comprehensive sensitivity analysis would be greatly appreciated. \n\n- Timing comparisons: I suspect that both BART and BAST are faster than the two other competing methods (SFS and inGP). Is this the case? \n\n- Complexity: What is the per-iteration time complexity? Will the worst-case complexity be O(M * n^2) as you may have to re-run Prim's Algorithm every time you propose a \"hyper\" move? How does this compare to the competing methods (BART, e.g., will have per-iteration time complexity O(M * n))\n\nClarity: The paper is well-written.\n\nSignificance: I think this paper represents a good step in pushing the boundaries of BART and moving away from strict axis-aligned splits. I anticipate it will be quite useful in spatial modeling. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "gzCrHNgEQu", "reviewer": "Reviewer_R75X", "summary": "Authors develop a novel Bayesian nonparametric regression model specially on known manifolds and complex constrained domains using additive RST partitions, which take advantage of flexibly shaped partitions. If the manifold is unknown, one should estimate it before applying the purposed method.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "- The main contribution of this paper is replacing the binary decision tree in BART model with random spanning tree, a flexibly shaped partitions. In this sense, the novelty is somehow limited per the standard of NeurIPS.\n\n- The experiments in the paper are not especially well-designed. For example, standard error of BAST goes up dramatically as $\\\\sigma$ increases, while BART told a different story. Authors should look into this or at least give some explanations.\n\n- The paper is well written with a clear mind. Readers outside this narrow area should find it readable.\n\n- Since the authors do not compare their method against some BART variants adapting to graph problem, it is hard for me to evaluate the significance of their work. In general, methods taking advantage of extra information (structure) should outperform those does not.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "There are a lot improvements can be made, a few are as followed:\n- Authors should find a way to fix the issues mentioned in the main review.\n- The number of tree should be larger, since it is essentially a truncation. Per the original paper, a recommended number is at least 50. However, the larger the better. Another way is to implement the Reversible-Jump MCMC.\n- Authors should examine whether BAST and BART converge or not after the purposed 20,000 iterations.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "KZQH7yjRMbt", "reviewer": "Reviewer_Tbpc", "summary": "This paper introduces a novel Bayesian model for nonparametric regression in complex spaces such as manifolds.  Let $S$ be the set of points $s$ at which the outcome $Y(s)$ is observed.  The basic idea is to take a known (or estimated) graph $G$ over $S$, and consider the class of piecewise constant functions over spanning trees of $G$.  The mean of the outcome $Y(s)$ at $s$ is then modeled as a finite sum of such functions $g_m$ for $m=1,\\ldots,M$.  An MH-within-Gibbs MCMC algorithm is proposed, in which the parameters of each $g_m$ function are updated in turn by making a move that preserves the full conditional distribution.  A well-chosen prior on trees facilitates exact sampling from the full conditional on spanning trees, using a previously known result from the literature. On simulated and real data, the method performs very favorably compared to competitors, including BART, SFS, and inGP.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I found the proposed model and inference technique to be quite interesting and thought-provoking.  The model appears to be novel, and while the individual pieces of the MCMC algorithm are not novel, the paper does a nice job of integrating these previous ideas.  The paper is well-written and clear, both in terms of exposition and mathematical presentation.  In terms of significance, if the good empirical results presented here generalize to other data sets, the method could be quite useful in many practical applications.\n\nMAIN COMMENTS\n1. Identifiability is not discussed in the article. Depending on how the model is used, this could be problematic.  For instance, if inferences on the individual component functions is desired, then there is a label-switching problem.  Also, lack of identifiability could affect MCMC mixing performance.\n2. Out-of-sample predictions (Sec 3.2): The rationale for the proposed method of making out-of-sample predictions was not clear to me.  It seems like a natural approach would be to consider a new model in which the new point $u$ is included in $S$, and look at the posterior predictive on $Y(u)$.  Does this coincide with the proposed approach?\n3. In Section 4.1, over what domain are the MSPEs and MAPEs computed?  For instance, are you averaging the prediction error over the whole domain $M$?  Or is this the error on the training set?  Or something else?\n4. It seems like the MCMC sampler might mix slowly since each component is updated individually, holding the others fixed.  This might make it difficult for the sampler to move around the parameter space.  I would recommend using standard diagnostic techniques to assess MCMC convergence.\n5. Computation time: How long does the BAST sampler take to run, compared to BART and the other competitors?\n6. On what types of problems does BAST not perform well?  I would imagine that high-dimensional spaces might be more difficult (for BAST compared to BART, for example), since the distances between points are less useful for defining local smoothness in a high-dimensional space.  Evaluating performance on a high-dimensional example would make the paper much more compelling.\n7. What about combining BAST with other types of functions, e.g., regression trees with axis-aligned splits like in BART?  A hybrid BAST-BART model might be interesting.\n\nMINOR COMMENTS\n- The paper assumes the data live on a compact Riemannian manifold. Is the compactness assumption needed?\n- Line 100: I think “cycle” is a more standard term than “circle”.\n- $M$ is used to denote both the manifold and the number of components in the model.\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"ethical_concerns": "No ethical concerns.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "-Qaq6-bGLfH", "reviewer": "Reviewer_kisJ", "summary": "This paper proposes a new Bayesian additive regression spanning trees (BAST) model for nonparametric regressions on manifolds. The model uses a new random spanning tree manifold partition model as weak learners, which allows it to capture complicated local discontinuity patterns in the true function while respecting intrinsic geometries of the domain. Building on nice properties of spanning tree structures, an efficient Bayesian inference algorithm is developed. Experiments on simulation data and the real chlorophyll data in Aral Sea show that BAST significantly outperforms other competing methods.", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "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": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is very well written and a pleasure to read. The Bayesian model and inference procedure in Sections 2-3 are clearly explained and technically sound. I found it is a solid piece of work as a whole and the BAST model is interesting itself and could be used in many applications.\n \nComments:\nLine 79: The word \"constrains\" should be spelled as \"constraints\".\nLine 120: Is would be useful to explain a bit more about the binary treed partition models, e.g., does it also need to specify the number k of clusters. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"ethical_concerns": "None.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This article introduces a novel Bayesian method for nonparametric regression on manifolds, referred to as Bayesian Additive Regression Spanning Trees (BAST). The basic idea is to build a sparse graph $\mathcal{G}$ on the observed covariate data points, and use a BART-like ensemble of tree-based models on $\mathcal{G}$.  The graph $\mathcal{G}$ captures the lower-dimensional geometry of the manifold on which the covariates live, and the tree-based models yield piecewise-constant functions on partitions of the data. A Gibbs sampler algorithm is provided for posterior inference. In particular, a well-chosen prior on spanning trees facilitates exact sampling from the full conditional on trees for each tree-based model component in the ensemble, using a previously known result from the literature. Experiments are performed on simulated and real data.

Generally speaking, the reviewers and I found the paper to be interesting and novel.  The paper is well-written and the method appears to provide a valuable contribution to this area of research.  The method provides a compelling improvement in performance relative to competing methods.

In my view, the main limitations are:

1) The computation time is quite costly compared to competitors.  For instance, on the U-shape example, BAST takes 651.49 seconds while BART takes only 15.83 seconds (see authors' reply to Reviewer Tbpc).  The authors mention that a more computationally efficient implementation of BAST is under investigation.

2) Performance is likely to degrade significantly as the dimensionality of the data grows; see point 5 by Reviewer Tbpc and the authors' reply.

3) The method for out-of-sample predictions seems *ad hoc* and could potentially be improved.  The method does not provide a coherent (i.e., projective) model for all points in the space, since it is defined only on the observed points due to reliance on the graph $\mathcal{G}$.  Thus, the model does not lead to a natural technique for making out-of-sample predictions, which is presumably why an *ad hoc* method was used.

The only score below the acceptance threshold was from Reviewer R75X, with a score of 5.  Although R75X has not responded further, I found the authors' reply to satisfactorily address the main criticisms and questions.
Decision: Accept (Poster)