Bayesian Transformed Gaussian Processes
Authors: Xinran Zhu
Leo Huang
Eric Hans Lee
Cameron Alexander Ibrahim
David Bindel
Venue: Accepted by TMLR ()

Abstract: The Bayesian transformed Gaussian (BTG) model, proposed by Kedem and Oliviera in 1997, was developed as a Bayesian approach to trans-Kriging in the spatial statistics community. In this paper, we revisit BTG in the context of modern Gaussian process literature by framing it as a fully Bayesian counterpart to the Warped Gaussian process that marginalizes out a joint prior over input warping and kernel hyperparameters.  As with any other fully Bayesian approach, this treatment introduces prohibitively expensive computational overhead; unsurprisingly, the BTG posterior predictive distribution, itself estimated through high-dimensional integration, must be inverted in order to perform model prediction.  To address these challenges, we introduce principled numerical techniques for computing with BTG efficiently using a combination of doubly sparse quadrature rules, tight quantile bounds, and rank-one matrix algebra to enable both fast model prediction and model selection. These efficient methods allow us to compute with higher-dimensional datasets and apply BTG with layered transformations that greatly improve its expressibility. We demonstrate that BTG achieves superior empirical performance over MLE-based models in the low-data regime ---situations in which MLE tends to overfit.  

Published in Transactions on Machine Learning Research (04/2023)
Bayesian Transformed Gaussian Processes
Xinran Zhu
m2584 &cornell eduu
Cornell University, Center for Applied Mathematics
Leo Huang'
ah839@cornell.ed edu
Cornell University, Department of Computer Science
Eric Hans Lee'
cric.lee com
SigOpt: An Intel Company
Cameron Ibrahim
Sarahim@udeer edu
University of Delowere, Department of Computer and Information Sciences
David Bindel
bindel@cornnlll.de edu
Cornell University, Deportment of Computer Science
Reviewed on OpenReview: https: // /pwp:///eeee.ee....mmmmm net/ form?? id= 4dCggggg)
Abstract
The Bayesian transformed Gaussian (BTG) model, proposed by Kedem and Oliviera in 1997, was developed as a Bayesian approach to trans-Kriging the spatial statistics com- munity. In this paper, we revisit BTG the context of modern Gaussian process litera- ture by framing as a fully Bayesian counterpart to the Warped Gaussian process that marginalizes out a joint prior over input warping and kernel hyperparameters. As with any other fully Bayesian approach, this treatment introduces prohibitively expensive com- putational overhead; unsurprisingly, the BTG posterior predictive distribution, itself esti- mated through high-dimensional integration, must be inverted in order to perform model prediction. To address these challenges, we introduce principled numerical techniques for computing with BTG efficiently using a combination of doubly sparse quadrature rules, tight quantile bounds, and rank-one matrix algebra to enable both fast model prodiction and model selection. These efficient methods allow us to compute with higher-dimennional datasets and apply BTG with layered transformations that greatly improve its expressibil- ity. We demonstrate that BTG achieves superior empirical performance over MLE-based models in the low-data regime situations in which MLE tends to overfit.
1 Introduction
Gaussian processes (GPs), also known as Kriging models in the spatial statistics community Cressie (1993) are a powerful probabilistic learning framework, including a marginal likelihood which represents the prob- ability of data given only GP hyperparameters. The marginal likelihood automatically balances model fit and complexity terms to favor simple models that explain the data well.
A GP assumes normally distributed observations. In practice, however, this condition is not always ade- quately met. The classic approach to moderate departures from normality is trans-Gaussian Kriging, which applies a normalizing nonlinear transformation to the data (Cressie 1993). This idea was reprised and expanded upon in the machine learning literature. One instance the warped GP (WGP), which maps the observation space to a latent space which the data is well- modeled by a GP and which learns GP hyperparameters through maximum likelihood estimation (Snelson al. 2004). The WGP paper employs
*These authors contributed equally this work.

---

Published in Transactions on Machine Learning Research (04/2023)
introduced compo- sitionally warped GPs (CWGP), which chain together sequence of parametric transformations with closed form inverses. Bayesian warped GPs further generalize WGPs by modelling the transformation as a GP [LazaroGrrddilee 2012). These are in turn generalized to Deep GPs by Damianou & Lawrence (2013). which stack GPs logether in the layers of a neural nelwork.
class of parametrized, hyperbolic tangent transformations. Later, Rios & Tobar 2019
Throughout this line of work, the GP transformation and kernel hyperparaameters are typically learned through joint maximum likelihood estimation (MLE). A known drawback of MLE is overconfidence in the data-sparse low-data regime, which may be exacerbated by warping (Chai & Garnett 2019 Bayesian approaches, on the other hand, offer a way to account for uncertainty in values of model parameters.
Bayesian trans-Kriging (Spöck et al. 2009) treats both transformation and kernel parameters in a Bayesian fashion. A prototypical Bayesian Mans-KKriiing model is the BTG model developed by Oliveira et al. (1997) The model places an uninformative prior on the precision hyperparamcter and analytically marginalizes it out to obtain a posterior distribution that is a mixture of Student's 's t-distributions. Then, it uses numerical integration scheme to marginalize out transformation and remaining kernel parameters. In this latter regard, BTG is consistent with other Bayesian methods in the GP literature, including those of Gibbs 19988 Adams ct al.
2009g; Lalchand & Rasmussen (2020).
While BTG was shown to have improved prediction accuracy and better uncertainty propagation, it comes with several computational challenges, which hinder its scalability and limit its competitiveness with the MLE approach.
First, the cost of numerical integration in BTG scales with the dimension of hyperparameter space, which can be large when transform and noise model parameters are both incorporated. Traditional methods such as Monte Carlo (MC) suffer from slow convergence. As such, we leverage sparse grid quadrature and quasi Monte Carlo (QMC), which have a higher degree of precision but require a sufficiently smooth integrand. Sccond, the postcrior mean predictor of BTG is not guaranteed to cxist, hence the neod to use the postcrior median predictor. The posterior median and credible intervals do not generally have closed forms, so one must resort to expensive numerical root-finding compute them. Finally, while fast vrossvaaliatton schemes www.oottt-roos--liiitiiiiioooooi (LOOCV) BTG, which incurs quartic
are known for vanilla GP models, cost naively, is less straightforward to perform because of an embedded generalized least squares problem.
In this paper, we reduce the overall computational cost of end-to-end BTG inference, including model prediction and selection. Our main contributions follow.
We revisit the Bayesian transformed Gaussian (BTG) model in the context of machine learning literature and GPs. We show that BTG fills an important gap in the GP research in which fully Bayesian treatment of input-warped GPs has not been thoroughly investigated.
We extend BTG to support multi-layered transformations.
We propose efficient methods for computing BTG predictive medians and predictive quantiles through a combination of doubly sparse quadrature and quantile bounds. We also propose fast LOOCV using rank-one matrix algebra.
We develop framework to control the tradeoff between speed and accuracy for BTG and analyze the error in sparsifying QMC and sparse grid quadrature rules.
We conduct comparisons between the Bayesian and MLE approaches and provide experimental results for BTG and WGP coupled with 1-layer and 2-layer transformations. We find evidence that BTG is well-suited for low-data regimes, where hyperparameters are under-specified by the data. In these regimes, our empirical testing suggests that BTG provides superior point estimation and uncertainty quantification.
We develop a modular Julia package for computing with transformed GPs (e.g., BTG and WGP) which exploits vectorized linear algebra operations and supports MLE and Bayesian inference
The code used run experiments this paper can be found https://github. com/xinranzhu//TC
2

---

Published in Transactions on Machine Learning Research (04/2023)
2 Background
In this section, we provide a brief overview of Gaussian processes and their transformed counterparts. This section is by no means comprehensive, and is meant to establish common notation and definitions for the reader's convenience. Note that we modified the standard GP treatment of Rasmussen & Williams (2008) to improve clarity when explaining the BTG model; notational differences aside, our formulations are equivalent.
2.1 Gaussian Process Regression
A GP / ~ G is a distribution over functions in RR where pu(x) is the expected value of [(x)) and F r"k(2,,') is the positive (semi)-definite covariance between f(x and f(z'). For later clarity, we separate the precision hyperparameter from lengthscales and other kernel hyperparameters (typically denoted by Unless otherwise specified, assuine linear mean field and the squared exponential kernel:
(a(x)) 3 m(xx)
m: Rd RR,
1 x x (B
ko(x,x') exp
2
Here m(x) is known function mapping a location to vector of covariates, 8 is vector of coefficients, and D is diagonal matrix of length scales determined by the parameter(s) 0. For any finite set of input locations, let:
= [1,,... xT,
X E Raxd
Mx = (m(x1). m(x))) fx
Mx E Raxp, fx E R",
f(x,))
where X is the matrix of observations locations, Mx is the matrix of covariates at X, and fx is the vector of observations. A GP has the property that any finite number of evaluations of f will have a joint Gaussian distribution: fx [3, T, NVMxx 1- -x),,), where (x-1 = - (2; x;) is the covariance matrix of fx. We assume Mx to be full rank.
Typically, the kernel hyperparameters T, and 0 are fit by minimizing the negative log likelihood:
- -log L(fx (x8,T,0)
/xx Mx + log Kx|-
This known as type maximum likelihood estimation (MLE) of the kernel hyperparameters. Having performed MLE to compute the optimal B, 7, and 9, the posterior predictive density of point æ is: Definition 1 (Posterior predictive density, GP model, type II MLE).
f(x) B, 0, Jx, X ~ N 80,3).
N(Ho, s0,8) is a normal distribution with the following mean and variance:
N0.3 33 33 + KXa X (fx Mx3),
$0,8 - K K Kxa),
where (Kx=); x). The Bayesian model selection approach forgoes MLE of the negative log like- lihood, and instead assumes a prior distribution over kernel hyperparameters p(3, 0). This inducos a posterior distribution over the kernel hyperparameters by noting that the posterior is proportional to the likelihood times the prior:
Standard GP literature includes additional hyperparameter which models the observed naise. omit in this section for brevity but treat in our experiments.
3

---

Published in Transactions on Machine Learning Research (04/2023)
p( 1xx, Xx & C(Jx (,8,7;;00(3,7,00
Given this posterior distribution, the posterior predictive density given by the following. Definition 2 (Posterior predictive densily, fully Bayesian GP).
p(f(z) fx, J Ja p(f(x) 18, T, 0, fx. p(8, fx.
We call this the fully Bayesian approach Gaussian process modeling. The well-known tradcoff between type II MLE and the fully Bayesian approach is that of tractability; the fully Bayesian approach requires computing integrals with no generic, analytic solution (often through MCMC) where type II MLE is much more straightforward.
2.2 Warped Gaussian Processes
While GPs powerful tools for modeling nonlinear functions, they make the fairly strong assumptions of
Gaussianity and homoscedasticityy WGPs (Snelson et al. (0044) challenge these assumptions by translorming
("warping") the observation space to a latent space, which itself is modeled by a GP. Given strictly increasing, differentiable parametric transformation 9A determined by a parameter X, WGPs model the composite function 9x of with a GP.
Definition 3 (Warped Goussion Process).
(gx o f) GP(US,
Let (9x(fx)); = (xx))) WGP jointly computes the parameters through MLE in the latent space, where the negative log likelihood
log ((xxfx) X, x
1 5 [1x(xx) - log Kx| log Jx-
Note that the difference between this equation and the negative log likelihood is the inclusion of J, which represents the Jacobian of the Iransformation:
Jx af)))) a2 (ʃ (xx))
WGPs predict the value of a point x by computing its posterior mean in the latent space and then inverting the transformation back the observation space: 9x (f())) Snelson al. (2004) uses the tanh transform family, whose members do not generally have closed form inverses; they must be computed numerically.
2.3 Other Transformed GP Models
There exist a variety of other transformed GP models, and this section we briefly summarize a few of them. This overview by no means comprehensive, but hope it succinctly communicates not only the breadth of existing work, but also the gaps that we believe BTG fills.
The compositionally warped GP (CWGP) (Rios & Tobar. 2019h considers a set of augmented warping functions used by the standard WGP, which possess analytic inverses (for computational simplicity) and are then composed logether for additional expressiveness.
The Bayesian warped GP (BWGP) (Lázaro-Gredilla 2012) takes a different approach transformations, placing GP prior on the warping function and variationally marginalizing out so that model prediction and selection does not require numerical integration.
4

---

Published in Transactions on Machine Learning Research (04/2023)
Table 1: table summarizing some of the different GP models in the literature today. The main differences among models are whether or not they use transformations (all do except for the standard GP), what sorts of transformations are used, and what type of inference method is used. Single denotes a single function used for transformations. Multiple refers to composition of such functions. At the bottom of the table is the topic of this paper, the Bayesian transformed GP (BTG).
Name
Transformations
Type
Inference
GP-MLE
No
N/A N/A
MLE
GP-Bayesian
No
Fully Bayesian
Warped GP
(WGP)
Yes
Single Analytic
MLE
Compositionally Warped GP
(CWGP)
Yes
Multiple Analytic
MLE
Bayesian Warped GP
(BWGP)
Yes
Single GP
Variational
Deep GP
(DGP)
Yes
Multiple GPs
Variational
Bayesian Transformed GP
(BTG)
Yes
Multiple Analytic
Fully Bayesian
2013 extends the "GP of GP idea by building network of GPs, and also uses variational inference to accelerate model prediction and selection by avoiding numerical integration. Indeed, a whole body of literature exists solely towards the application of variational inference to Gaussian process regression (Tran et al. 2015 Salimbeni & Deisenroth 2017 Jakkala 2021
Deep GPs (Damianou & Lawrence
Like all these other modols, BTG considers a family of warping functions single layer or multiple com- positions. However, instead of using lype II MLE (which can sometimes be overconfident) or variational inference (VI) (whose approximations often require significant data before showing their advantages). BTG opts for carefully chosen prior and explicit numerical integration during the model selection process for a fully Bayesian approach to transformed Gaussian process regression.
While we focus on improving the Bayesian treatment and scale BTG to higher dimensions using exact kernel inference, another line of work, orthogonal to ours, focuses on combining the warping setting with sparse GP techniques to scale to large datasets. Graßhooff et al. (2020) adopt the Structured Kernel Interpolation (SKI) method (iilson & Nickisch, 2015) to the WGP setting, which scales WGP up to large training data in low dimensions. Rossi al. (2021) considers the sparse GP setting with inducing points and focuses on a fully Bayesian mreatment of both the inducing points and model hyperparameters, based stochastic gradient Hamiltonian Monte Carlo.
3 The Bayesian Transformed Gaussian Process (BTG) Model
One might think of the Bayesian Transformed Gaussian (BTG) model (Oliveira et al. 1997 as fully Bayesian generalization of WGP. WGP uses MLE to learn optimal values for the transformation parameters X, mean function parameters 3, signal variance parameters т, and GP lengthscales 6.
Just like WGP, BTG models function f(x) as:
(gx o f) B 1,,T,,,00 OP(HS,T-11
Unlike WGP, BTG places joint prior p(B, r, 8, A) over all model parameters and performs fully Bayesian model selection by computing the usual marginal likelihood instead of performing the MLE. A key idea of BTG is to carefully select joint prior over X, 8, T, and 0 simplify the computation so that the BTG's posterior predictive density is mixture of t-distributions -hhis will be derived in later subsections.
Definition 4 (Posterior predictive density, BTG).
M Σ 1=1
f(x) Ix. ~
w. T Tm
5

---

Published in Transactions on Machine Learning Research (04/2023)
Gaussian Process Model
Bayesian Transformed Gaussian Process Model
3.5
3.0
3.0
2.
2.5
1.5
1.0
True Function
Training samples
GP Mean
GP Uncertainty
BTG Mean
BTG Uncertainty
Figure 1: A comparison of GP and BTG, both trained on 48 random samples from the rounded sine function with noise from (O0, 0.05). Left: the predictive mean and 95% equal tailed credible interval of Gaussian process using the squared exponential kernel. Right: the predictive median and 95% equal tailed credible interval of BTG using the squared exponential kernel. BTG better able to capture the irregularity in the training data despite a highly smooth kernel, largely due to it's Bayesian treatment of warping functions that map the data points to latent space.
where each T: is a different transformed distribution whose mean and variance we will derive in $.2 and $ for a comprehensive analysis, see Box Cox (1964). Tike the case of the WGP, this predictive distribution must then be inverted to perform prediction or uncertainty quantification.
BTG appears to possess the same advantages that any fully Bayesian model will possess over it's MLE counterpart: robustness to overconfidence. This results in a better fit, demonstrated in Figure 1 The underlying kernel is the highly smooth squared exponential, resulting in a GP that is unable to capture the underlying function, which not smooth at all. Due to a combination of input warping and Bayesian inference, BTG able to resolve the underlying datapoints much better than a GP despite using the same smooth kernel. In later sections, we explore the advantages of being Bayesian the low-data regime.
3.1 The Parameter Posterior
A key idea of the BTG model is that, conditioned on X, 0, and fx, the resulting WGP is a generalized lincar model (Oliveira et al. 1997) We estimate by Bo, the solution to the weighted least squares problem:
90.7 min [[x(fxx) Mx x
where qu is the residual norm. Explicity, we have
9x(xx)
00.A = (gx(xx) " (gx(xx) 3o.3.
Therefore, we have the marginal probability
p(fx x) / p( X)x
/ 1 Mx -
-1/2
1/2 2 J1-p/n
(nn-p//2
p()pppppp)()))))))p))))))))))))p da dx.
10,>
Further using Bayes theorem, we have
p00, fx, X) x 1
Mx
-1 (n-p)/2 ( J1-p/ (1
1/2 -1/2
40.2
(1)

---

Published in Transactions on Machine Learning Research (04/2023)
On the other hand, because appropriate values for T, and 0 depend nontrivially on X, BTG adopts the improper joint prior:
p(8,T,8,A) x ((0)p(A) (FJR(")
This effectively decouples the priors 0 and A; p(0) and p\)) independent and user-specified. BTG then selects the prior on (5, 7) to follow the conditional normal-inverso-game distribution:
8 81
T Mx
n2 - 2
(2)
+ =1,,0,1,,X~1 X~ r-
2 9x.0
Note that, due to the improper joint prior, the distribution 8 and depends the values of 0 and X.
Using 1 and the posterior distribution, up to some constant, is determined:
rBB fx. X) = p(3, 1,,,,0,,,,, 0,A, X)p(0, Xx;
The reason for the normal-inversed gamena prior becomes apparent; the marginal of normal-inverse- garmma is t-distribution, and once marginalize out equation 3 will be the result. We derive this in the following subsection.
the mixture of t-distributions seen in equation
3.2 The Predictive Density
BTG bases the prediction at æ on the Boyesian predictive density function (Aitchinson & Dunsmore
1975
p(f(z) fx, f
f p(f(x)
fx X) d3 d0 dr.
(4)
The Bayesian predictive density function is simply an integral of a product two probability densitics. In the previous section, we derived p(3,T, T, X, [x, so we know the second probability density. The first probability density that of a warped GP (Definition
9A f(x) B,,, T, X,0, fx, (mm3... Do). m3300 = 8T m(x) + K&_ Kx (ox fx Mxx3),
Dg - K K KxØ).
We consequently know everything needed to specify Equation Because adopted conditional normal- inverso-gamma distribution, once we marginalize out and t, the posterior of gx o fxx) is the i-distribution:
gA f(x) X, 9, fx. mm (g0,ACo,x)
(5)
where the mean largely resembles that of GP:
mx, K K (fx Mx {oom()))
and Cx,8 is the final Schur complement B(x)/[ki (x,
(x,x)) x]] of the bordered matrix:
0
Mf
m(x))
B(x) = Mx
Kx
Kx ko(x, x:)
m(x)
K
Therefore, by Bayes theorem, the marginal posterior of BTG is:
p((()) Xx Se.1 o(f(x)
0, A, fx. X)p(fx X. (0))(()(()))
(6)
Je p(0, X)x0))))(

---

Published in Transactions on Machine Learning Research (04/2023)
where the posterior being integrated is transformed f-distribution with the pdf
((((z) 0, X)
2 I( T 1/2 12 Co.x (f(x) mm,x)) "0ax) -ma,)]] (7)
Unlike WGP, BTG may not have first or second moments, because its marginal posterior may be for example, mixture of log-t distributions. this occurs, the probability density function (pdf) will not. have mean or variance. Therefore, BTG instead uses the median and credible intervals, computed by inverting its cumulative distribution function (edf).
3.3 The Approximate Predictive Density
For general nonlinear transformations, the posterior distribution of BTG (Equation 6 is intractable and therefore we approximate it using set of M quadrature nodes and weights ([0,, Ar), wi), yielding the mixture of transformed t-distributions
p(f(x) 1xx, = E wio(f(x) 0n A, A, fx. X),(fx A.)()))))) Enn w.p(O, A [x))(0)))(AA)
Here, (f(x)] Ai, X) is the pdf of the transformed I-distributions (Fquation 3.2 p(fx(O;,A) is the likelihood of data given hyperparameters, p(0;) and p(A;) are the user-specified hyperparameter priors, and w; is a quadrature weight at the quadrature point Az).
To simplify notation, we combine terms except for the transformed t-distribution pdf into weights i; simplify the predictive distribution into our desired mixture of transformed t-distributions:
w:p(Sx 0r, Ae (0)((())pAA))
W; = MM E wop(O;, As (x)r((:)p(A.))
(8)
M
p(f(x) | fx, s
o.)(f(z) 0, no Xx.
As mentioned earlier, p(f(z) X) is not guaranteed to have a mean, so we must. use the median predictor instead, We do by computing the quantile p-1(0.5) by numerical root-finding, where P is the edf of p(f(z) X), and therefore mixture of transformed t-distribution edfs.
4 Methodology
BTG regression via the median predictor (or any other quantile) of Equation challenging. The dimension- ality of the integral scales with the total number of hyperparameters, which contains the kernel lengthscale hyperparameters that grow with the ambient dimension of the data, as well as the transformation pa- rameters À that grow with the number of transformations used. Furthermore, its edf must be numerically inverted, requiring many such quadrature computations for even a single, point-wise regression task. This is further complicated by the difficulty in asscssing model fit through LOOCV, which must be repcated every quadrature node as well. As a result, naive implementation of BTG scales poorly.
In this section, discuss and propose efficient algorithms that make BTG model prediction and model validation far faster, and indeed, comparable to the speed of its MLE counterparts. First, we discuss our doubly sparse quadrature rules for computing the BTG predictive distribution ($4.1 and $4.2). We then provide quantile bounds that accelerate the root-finding convergence ($4.3). Next, we propose an O(n3) LOOCV algorithm for BTG using Cholesky downdates and rank-1 matrix algebra ($4.4), which greatly improves the naive Onn)) LOOCV cost. Finally we discuss the single and multi-layer nonlinear transformations used in experiments (64.5).
8

---

Published in Transactions on Machine Learning Research (04/2023)
4.1 Sparse Grid Quadrature
Sparse grid methods, Smolyak algorithms, are elfective for approximaling integrals of sulficient regularity in moderate to high dimensions. While the conventional Monte Carlo (MC) quadrature approach used by Oliveira et al. (1997) converges at the rate of O(1/VM), where M is the number of quadrature nodes, the approximation error of Smolyak's quadrature rule is O [M MM M((d-1)+)))) where d the dimensionality of the integral and is the integrand's regularity i.e., number of derivatives.
In this paper, use sparse grid rule detailed by Bungartz & Griebel 2004 and used for likelihood approximation by Heiss & Winschel (2008).
4.2 Quadrature Sparsification
Numerical integration schemes such sparse grids and QMC M fixed quadrature nodes, where M depends the dimensionality of the domain and fineness of the grid. In the Bayesian approach, expensive GP operations such as computing log determinant and solving a linear system are repeated across quadrature nodes, for a total time complexity of (MM)).
In practice many nodes are associated with negligible mixture weights, so their contribution to the posterior predictive distribution can effectively be ignored. We thus adaptively drop nodes when their associated weights fall below certain threshold. To do in principled way, we approximate the mixture with subset of dominant weights and then quantify the error in terms of the total mass of discarded weights.
We assume the posterior edf is the mixture of cdfs F(x) = E (xf( where each fi is a edf. Assume the weights }w are ordered by decreasing magnitude, Consider Fx(x), a truncated and re-scaled F(r). We first quantify the pointwise approximation error in Lemma 4.1 We then quantify the error in quantile computation: Propositions 4.1 4.2 show that the approximated quantile F (p) can be bounded by per- turbed true quantiles. Proposition 4.3 gives simple bound between Fr (p) and F-) within the region of interest, and applies to both QMC which uses positive weights and sparse grid quadraturo-whie uses positive and negative weights. See Appendix A.I for proofs.
Lemma 4.1. Let k be the smallest inleger such that Σ mizture
wi 2 -e. Then define the scaled, trancaled
(x(x)) = Σ wi wiJi(r),
c:=
Wi.
We have
F(())F((x))) < 2c.
Proposition 4.1 (Error Bound for Positive Weights). For any e (0, 1), let k be the smallest integer such that Σk Wi 2 1 c. Define the scaled, truncated mixture
Fa(z) = Σ wif((x), c:=
Let p E (0, 1) and assume that p 2c 1). Then the approsimate quantile (p) is bounded by perturbed true quantiles:
F-1(p- 2x) < (p) < F- (c)
Proposition 4.2 (Error Bound for Negative Weights). Let F(x) be defined as before, except each w; is longer required to be positive, Consider the split F(x) = FM (I) + (x) where Fxr(x) E wi wf((x)
AM and (x)(x) -M Ei=M4 wifi(x). Then for any x. we have Ras (Ir) [e_,etl, where the epsilons defoned the sum of positive (resp. negative) weights of (xr(z)
Σ [w]]- 50, 4
Σ
[w+]+ 0.
9

---

Published in Transactions on Machine Learning Research (04/2023)
Let 10 and assume p c-,p c E (0,1). Then the approzimate quantile F (p) is bounded by perturbed true quantiles:
(p+e) S FM (p) - (p++).
Proposition 4.3 (Error Bound at a quantile). Let F(x) be defined as before, E1, E (0, 1), and Fx(z) an approzimate to F(x) such that (p - <1) & (p) S (p- (2) for some E (0, 1). Assuming - E1.P 62 E (0, 1), we have the following error bound at a quantile,
dF-1
-F
FD)
max
dx (5)
where maxff (2)
By adaptively sparsifying our numerical quadrature schemes, we are able to discard a significant portion of summands the mixture F(x), which in turn, enables significant specdup of BTG model prodiction. Empirical results are shown in
4.3 Quantile Bounds
To compute posterior quantiles, we apply Brent's 's algorithun, a standard root-finding algorithm combining the secant and bisection methods, the edf defined in Equation [} Since Brent's algorithm is box-constrained, we use quantile bounds to narrow down the locations of the quantiles for p {0.025, 0.5, 0.975).
Let F(x) = EM uxfi(z). Then we have the following bounds for the quantile F-1(p)- See Appendix A.2 for proofs. Proposition 4.4 (Convex Hull). Let F(x) be defined as before with w; > 0 and E wi 1. Then min Ji (p) s F (p) I max fi (p).
Proposition 4.5 (Singular Weight). Let F(x) be defined as before with wi > and EM wi = 1. Let
W. Wy. Then
J (p- T.,) SF F (p)
min p PT
( w.)-
P- w. 20
When solving for = F-) we run Brent's algorithm using our quantile bounds the box ronstsaints. Furthermore, we adaptively set the termination conditions xtol and ftol to be on the same order of magnitude as the error in quadrature sparsification from $4.2 This greatly accelerates convergence in practice. A
performance comparison between quantile bounds outlined in this section can be found in $
4.4 Fast Leave- Satee-One-Qa Oite-Ot-Con Suut-Cross-Aae Crross-Validation (LOOCV)
LOOCV a standard measure model fit: in practice, it is most commonly used for tuning hyperparameters and for model selection. While fast LOOCV schemes are known for GP regression, it less straightforward to perform LOOCV on BTG. In particular, the computational difficulty lies in two LOOCV sub-problems: a generalized least squares problem and principle sub-matrix determinant computation. These correspond to the terms in the BTG likelihood function and the BTG conditional posterior in Equation K Being Bayesian about covariance and transform hyperparameters introduces additional layers of cost: LOOCV must be repeated al each quadrature node in the heerpaarameter space. This further motivates the need for an efficient algorithm.
For notational clarity, let (-i) denote the omission of the ith point. For a kernel matrix, this means deletion of the ith and column; for vector, this indicates the omission of the ith entry. We seek compute the
(-i) (C, ) -1/2 of the t-distributions (Equations 5 for each (-i) and standard deviation σg, mean ma
(-1) submodel, obtained by lcaving out the ith training point. Specifically, computing {9b Jim1 entails solving
the generalized least squares problems for 1....,
arg min ||( - Mx 1-0 (()0 x
||
((-1)
10

---

Published in Transactions on Machine Learning Research (04/2023)
Table 2: Elementary Transformations: analytic function forms and parameter constraints. Parameters are assumed to be in R unless stated otherwise.
Name
g(y)
Req.
Num. Params
Affine
+ by
b > 0
1
ArcSinh
+ b asinh
b d > 0
4
d
SinhArcSinh
sinh(b asinh(y
a))
> 0
2
if >
0
Box-Cox
2 0
log(y)
(-i) entails solves with KX-), which naively takes where Y go fx. In addition, computing C/,A)) and max O(n²) per sub-problem. Therefore, the BTG LOOCV proceedure naively takes Onn²) total time.
We develop an O(n²) fast LOOCV algorithm for BIG using three building blocks: fast determinant com- putations (Proposition 4.6), fast abridged linear system solves (Proposition 4.7) and fast rank-one O(p²) Cholesky down-dates (Proposition 4.8). We refer to (Stewart (1998) for the rank-1 Cholesky downdate algo- rithm. For algorithm details as well as proofs, see Appendix B The scaling behavior of our fast LOOCV algorithm shown in Figure 4
Proposition 4.6 (Determinant of Principal Minor).
det (2<(-) det(2)
Proposition 4.7 (Abridged Linear System). Let K E ux be of full rank, and let e, E IR" satisfy Ke = Then if r: = c//e/ K-le, et, we have:
) = (K(-1))-1 = - perrrK-'e
Proposition 4.8 (Rank one matrix downdate). If X puxm with m < n has full column rank and is a positive definite matriz in Raxn, then we have
E I-le,eTE- X,
x(- x((0) -X EΣ-
where ~(-)) E (1=-1))((-1) is the (i, i)th minor of E and e, is the ith canonical basis vector.
4.5 Transformations
The original BTG model Oliveira et al. (1997 uses the Box-Cox family of power transformations and places an uniform prior on A. Recent research has greatly expanded the set of flexible transformations available. Snelson et al. (2004) uses sum of tanh(-) transforms in the WGP model and Rios & Tobar (0019 composes various transformations to provide a flexible compositional framework in the CWGP model.
We apply BTG with more elementary transformations and compositions thereof, summarized in Table 2 As
show in $ these compositions have greater expressive power and generally outperform single transfor- mations, at the expense of greater computational overhead.
5
Experiments
This section contains:
Experimental details necessary to reproduce our experiments, including the datasets we used as well the decisions made regarding the BTG and GP models (such the kernel).
11

---

Published in Transactions on Machine Learning Research (04/2023)
Comparing likelihood functions when n=5 and when n=30
Likel reelioodd
Function Values
Likelihodd
02-
Function Values
02-
63-
63
60-
22
10.3
12.3
82
Function Domain
Figure 2: Plots of marginal log likelihoods for ShnhArreSin--ransfforneedd WGP and corresponding training data. In the first four columns, a, are transformation hyperparametors, and 0, o² are kernel hyperparame- ters; the optimal values of marginal log likelihoods are marked with a red star. The last column shows the underlying function (dashed curves) and the training data (red dots). The top row represents a data-sparse setting with training points, while the bottom row represents data-rich setting with 30 training points. We can that in the data-sparse setting, optimal marginal log likelihood values with respect to and o² are clearly not as well defined as the dala-rich setling.
small example highlighting that when data is sparse, the likelihoods for transformations and kernel hyperparameters be very flat, thus demonstrating the need for BTG in these scenarios.
Experiments validate the efficiency of our computational techniques.
Thorough regression experiments, which demonstrate BTG's strong empirical performance when compared to appropriately selected baselines. We also provide timing results of our fast BTG implementation, showing that BTG has a comparable computational cost to WGP.
We hope these experiments provide a thorough and nuanced picture that highlights when and why BTG can perform better than it's GP counterparts. All methods possess their stretgs and weaknesses, and BTG is not exception.
5.1
Motivation for the Bayesian Approach
5.2
Experiment Details
We provide the datasets use, the performance metrics measure and additional implementation details necessary to reproduce our experiments in this subsection.
Two synthetic datasets: IntSine and SixHumpCamel. The IntSine dataset, also used by Lázaro- Lázaro-Gredilaa (2012), sampled from a rounded 1-dimensional sine function with Gaussian noise given variance. The training set comprised of uniformly spacod samples on [- vr. The testing set consists of 400 uniformly spaced points on [- m). The SixHumpCamel function 2-dimensional benchmark optimization function usually evaluated on [-3,3] 3] [_ (Molga & Smutnicki, 2005). We shift its values to be strictly positive. The training set is comprised of 50 quasi-uniform samples, i.c., 2d Sobol sequence, from [-1, [- The testing set consists of 400 uniformly distributed points on the same domain.
Three real datasets: Abalone, WineQuality and Creep. Abalone an 8-dimensional dataset, for which the prediction task determine the age of an abalone using eight physical measurements (Dua Graff] 2017). The WineQuality dataset has 12-dimensional explanatory variables and relates the quality of wine to input attributes Cortez et al. 2009) The Creep dataset is 30-dimensional and relates the creep rupture
12

---

Published in Transactions on Machine Learning Research (04/2023)
Weight Concentration
Prediction Error
Prediction Time
102.
4.05
Percent (%)
4.00
Time (s)
1025
MSE
3.95
102.
3.90
0.7
101.
3.85
0.6
3.80
10
10
10
1025
101.9
10
10
1023
10th
10th
10th
1025
# Nodes
# Nodes
# Nodes
Figure 3: Comparison of sparse grid quadrature (blue) and QMC (red). Weights are ordered by decreasing magnitude. (Teft) Number of nodes versus percent of total mass they comprise (Middle) Number of nodes versus BTG prediction MSE (Right) Number of nodes versus prediction time. We can see that sparse grid quadrature has prediction error that decreases much faster than QMC.
stress (in MPa) for steel to chemical composition and other features (Cole al. 2000). To simulate data- sparse training scenarios, randomly select training samples of size 30, 200, and 100 from Abalone, WineQuality, and Creep, respectively, and test 500, 1000 and 1500 out-of-sample points.
5.2.1 Sparse Grids QMC
Our empirical experiments involve three loss functions to evaluate model perfor- mance: root mean squared error (RMSE), mean absolute error (MAE) and mean negative log predictive density (NLPD). Let (s*())) be true test labels and {f(a:)) be predictions, which are taken to be predictive medians in WGP and BTG. The predictive median is also used by Snelson et al. (2004). Let p be the predictive distribution from model. The RMSE, MAE and NLPD are defined as follows
Performance Metrics:
RMSE
1 (f)(s)) - f(mx))2
1=1
n
MAE Σ
(x;) - f((xs),
n i=1
1 Σ n
NLPD
log(p(f((xx))
Implementation: We run all experiments using our Julia software package, which supports a variety of models (WGP, CWGP and BTG) and allows for flexible treatment of hyperparameters. We also implement several single and composed transformations. For ML.E optimization, we the L-BFGS algorithm from the Julia Optim package (Mogensen & Riseth 2018
Kernel: We used the squared exponential (RBF) kernel for all experiments:
ka( x²) = 72 exp (- 1 --'x'13 x'I1b,>) 11-11=-' =)
o² 0gg
To model observation input noise for BTG, we add regularization term to make the analytical marginalization of mean and precision tractable. We also assume the constant covariate m(x) = 1, in the BTG model, and normalize observations to the unit interval. We assume the constant mean field for both BTG and WGP.
Model:
Our motivation for using a Bayesian approach like BTG rather than a MLE-based one like WGP is that the maximum value of the marginal log likelihood in a data-sparse setting not as well-defined as that in a data-rich setting. We demonstrate this motivation using 1D toy example. We examine the marginal log
13

---

Published in Transactions on Machine Learning Research (04/2023)
Table 3: Compute time (seconds) for computing predictive median and credible intervals using no quantile bound, convex hull quantile bound, and singular weight quantile bound. Results are averaged over 10 trials, We find that both bounds significantly reduce the total time used compute the predictive median and credible intervals -thouhh the convex hull bound performs better average than the singular weight bound, likely because the former is often a better bound than the latter.
Total (s)
Median (s) Credible Interval (s)
Baseline
13.0
3.24
7.87
Convex Hull
6.21
1.11
3.19
Singular Weight
11.0
2.52
6.54
likelihoods of transformation and kernel hyperparameters in the WGP model in Figure synthetic function, the 1D Levy function:
We consider a 1D
(x: + 3)
x -m111 4 1)
1
1++ (1+ sin"
(x 3) ))
f(x) = sin2
+
2
As a toy example, use 5 training data for a data-sparse setting and use 30 for data-rich setting, We use a SinhArcSinh transformation for both settings, which is parameterized by two parameters and (see Table 2. We thus have four total hyperparameters (a, b, 8, o). and we plot their marginal log likelihoods in Figure 2 for both the data-sparse and data-rich setting. We also identify the optimal hyperparameters (a*, ",,,b", 0", o²) using MLE, which we plot with red star.
We observe that in the data-sparse setting, the marginal log likelihoods tend to be quite flat, even log scale, which is especially true of I and o². This means that many possible hyperparameters explain the data, and consequently that an MLE-based model (which selects only one hyporparameter value) will likely be a poor fit. This also suggests that a fully Bayesian approach which is what BTG does could be more appropriate in the data-sparse setting than MLE estimation.
In the data-rich sctting, the distribution of and o² are tightly concentrated, and thus the MLE approach does make sense.
5.3 Scaling Experiments
We demonstrate the efficiency of the proposed computational techniques discussed in $ First, we illustrate our quadrature sparsification schome discussed in $4.2 on two different quadrature rules, and compare the performance under varying sparsification levels. Second, we evaluate the fast quantile bound for root-finding used BTG prediction proposed in $4.3 by comparing the total time cost of BTG with WGP. Last but not least, we verify the complexity improvement of the proposed fast cross validation.
We compare sparse grid and QMC quadrature rules under our quadrature sparsification framework. We use the SixHumpCamel dataset with 30 training data points and 100 testing data points toy problem. We train BTG with the composed transformation Affine-SinhAreSnntt The hyperparameter space is 7-dimensional. In our experiment, we begin with a handful of quadrature nodes, and gradually extend this set to the entire quadrature grid. As more and more quadrature nodes are incorporated into the approximate integral (the approximate postcrior distribution), we record the percentage of total weight that is used, and then compare how the prediction MSE and the prediction time cost vary. Intuitively, increasing quadrature nodes would reduce the prediction MSE and increase the prediction time cost. To better study our quadrature sparsification framework and balance the prediction accuracy and cost, we closely compare such phenomenon on the two quadrature rules allowed in BTG. Results are plotted in Figure K
From Figure 3 we observe that the sparse grid quadrature rule yields lower prediction MSE: QMC converges an MSE of 3.99, while sparse grid converges to an MSE of 3.80. We also observe that the sub-grids have similar weight concentration, which we showed was a proxy for quadrature approximation error in Therefore, as the mass of the dropped weights falls below 0.1, the error in the integration scheme
14

---

Published in Transactions on Machine Learning Research (04/2023)
10²
Naive LOOCV
Cholesky Downdate LOOCV
Ojn"
108
109
Time (s)
10a
Oin*)
10a
10
1.5
2.5
N
x10
Figure 4: Compare the LOOCV timing of computing posterior distribution parameters of n sub-moder with the proposed Cholesky downdate and without. It is verified that the naive way has Onn*) cost while our proposed Cholesky downdate reduce the cost, by factor of n, to O(n²).
increasingly be attributed to the error in the quadrature rule itself as opposed to sparsification. Since the error of the sparse grid rule decays faster than that of QMC, we expect sparse grid prediction error to also decay more quickly and this trend is supported by Figure We empirically find that the joint-likelihood function is sufficiently smooth, that sparse grids are effective. Finally, confirm that inference time scales linearly with the number of quadrature nodes.
5.3.1 Quantile Bounds Speedup
To assess the effectiveness of quantile bounds for root-finding. we record BTG prediction times using the convex hull bound and singular weight bound proposed $4.3 We set up a problem using the Levy1D dataset Molga & Smutnicki (2005) and 200 training points. The tolerance for the root-finding Brent's algorithm is set be 10-3 We record the total time cost of prediction, the time cost of computing medians and time cost of computing credible intervals. Intuitively, using
precise bound in root-finding algorithm would greatly improve convergence speed, and this is verified with results in Table 3
We find that the convex hull bound decreases the overall computational overhead by factor of at least two. The convex hull bound outperforms the singular weight bound for finding credible intervals and in overall time, but the singular weight bound was faster for finding the median many scenarios.
5.3.2 Fast Cross Validation
To assess our fast LOOCV proposed in $4.4 we infer the sub-model posterior distribution moments on toy problem in two different ways: with and without using Cholesky rank-one downdates Rx to compute the Cholesky factors for R We plot timing results Figure which confirm that our O(n²) method scales significantly better than the naive O(n²) method. We note that these timings are generated using exact GP inference, and that. more generally our fast cross validation is O(n) faster than standard
validation. We note that this speed-up applies not only to the typical O(n3) cost of dense GP inference, but also any other scalable GP method, such as variational GPs which would be Onn²) cost instead.
5.4 Regression and Uncertainty Quantification Experiments
Our efficient algorithms allow us to test BTG on set of synthotic and real-world problems. We consider datasets: IntSine, SixHumpCamel, Abalone, Wine, and Creep with dimension ranging from 1 30. The total dimension of the hyperparmeter space is further increased by transform parameters by much as 8.
15

---

Published in Transactions on Machine Learning Research (04/2023)
Table 4: comprehensive set RMSE and NLPD results for prediction using the IntSine, SixHumpCamel, Abalone, Wine and Creep datasets. We generally kept the training sets small, to highlight BTG's advantages in the data-sparse regime, though for the IntSine function there was good amount of data. We compare the following _ GP, DeepGP, WGP, CWGP, BTG, and use the following transformations transformations: I: Identity, A: ArcSinh, SA: SinhArcSinh, BC: BoxCox, L: affine. The best method is bolded; we find that BTG has lower RMSE and NLPD than the other models, though unsurprisingly, different transformations are more suitable for different datasets.
IntSine
Camel
Abalone
Wine
Creep
RMSE
NLPD
RMSE
NLPD
RMSE
NLPD
RMSE
NLPD
RMSE
NLPD
GP
0.227
-1.179
2.003
11.20
3.290
1.120
1.994
-0.074
37.88
-1.189
DeepGP
0.241
0.013
1.344
1.761
3.212
3.988
0.824
1.377
91.95
20.39
WGP-A
0.179
-1.392
2.055
44.10
3.097
1.105
0.811
-1.108
35.48
-1.421
WGP-SA
0.172
-1.693
2.012
45.10
2.992
4.879
0.809
-1.269
61.62
1.001
WGP-BC
0.184
-13771
1.964
44.70
2.826
0.837
1.045
-0.186
40.38
-1.357
CWGP-L-SA
0.172
-1.693
2.055
13.80
3.117
3.444
0.808
-1.415
39.71
-1.134
CWGP-A-BC
0.174
-1.457
1.960
37.90
3.088
1.151
0.808
-1.416
37.89
-1399
BTG-T
0.169
-1.429
1.820
2.342
2.804
-0.070
0.808
-1.356
38.11
-1.273
BTG-A
0.168
-1.44
1.827
2.331
2.890
-0.174
0.807
1.360
38.65
1.239
BTG-SA
0.165
-1759
1.801
0.907
2.791
-0.123
0.820
-0.56
91.69
-0.358
BTG-BC
0.170
-1439
1.675
0.659
3.225
0.110
0.808
1.358
39.10
-1.32
BTG-L-SA
0.145
-1.781
1.673
0.730
2.871
-0.023
0.809
-1.399
91.83
-0.374
BTG-A-BC
0.159
-15516
1.828
2.330
2.832
-0.299
0.802
-1.361
35.25
-1.433
We compare with standard GP model, WGP models and CWGP models using transformations and their compositions from Table 2
For each of these datasets, we perform both regression and uncertainty quantification (UQ) by computing the RMSE and NLPD, which we record in Table ! Our two metrics evaluate both the predictive mean and predictive variance, respectively: RMSE the root mean squared error metric evaluating predictive mean, and NLPD the mean negative log predictive density metric evaluating both the predictive mean and predictive variance (uncertainty). For space's sake, record additional performance metrics such MAE in the appendix (Tablee
We also compare end-to-end inference time cost in Table to give the reader a better idea of the computa- tional overhead associated with BTG. This compute time generally depends on the dimension of the integral, ie., the total number of hyperparameters that we must marginalize out the fully Bayesian approach, and thus the numbers vary by a bit depending on the datasct.
Our results largely favorable. We find that composed BTG models tend to outperform single- transformation BTG models, which themselves tend to outperform their GP, WGP, and CWGP counterparts. The DecpGP model performs poorly in most of these problems due to the lack of enough data, which would otherwise allow it to overcome the approximation error associated with variational inference. More generally, BTG outperforms other baselines on almost all problems for both regression and uncertainty quantification (UQ), and quite convincingly certain cases for example, BTG achieves an NLPD around two orders of magnitude better than existing models in the Camel dataset.
These results demonstrate the improved flexibility afforded by layered transformations, and is evidence of the superior performance possible with fully Bayesian approach small to medium sized datasets.
WWe omnt the tanh(-) transform used in the original WCP paper Suelson al. 2004 since Rios & Tobar 2019 shows that the elementary trausforms in Table are competitive with tanh(-). We tried compare against BWGP, but could find code release available. We suspect BWGP possess similar performance the DeepGP, both variational inference and highly flexible families of warping functions.
16

---

Published in Transactions on Machine Learning Research (04/2023)
Table 5: End-to-end regression times for SixHumpCamel, Abalone and Wine datasets, using the following models: GP, WGP, CWGP, BTG, and the following transformations: SA:SinAAreSinh, BC:BoxCox, Laffiin. We note that the timings for BTG depends on the number of quadrature nodes uses, but made sure that these settings consistent with those of Table 4 in which BTG outperforms GP, WGP, and CWGP.
SixHumpCamel
Abalone
Wine
| Dim
Time (min)
Dim
Time (min)
Dim
Time (min)
WGP-BC
5
0.68
10
1.28
1.52
WGP-SA
6
0.78
11
1.22
15
1.60
CWGP-L-SA
8
1.08
13
1.48
17
2.50
CWGP-A-BC
9
1.14
15
1.56
18
2.82
BTG-BC
4
1.10
9
1.02
13
1.40
BTG-SA
5
0.95
10
1.04
14
1.29
BTG-L-SA
1.74
12
1.11
1.79
BTG-A-BC
8
1.65
1.09
1.87
Additionally, due to efficient numerical methods, the end-to-end inference time of BTG is generally comparable to other baselines thee compute time of BTG is only bit higher than that of WGP.
6 Conclusion
In this paper, we revisited the Bayesian transformed Gaussian (BTG) model, which variant of Bayesian trans- Kriging. BTG can be seen fully Bayesian treatment of input-warped GPs, and we believe it fills important niche in the GP literature, which has yet not explored the combination of transformed GPs and Bayesian inference.
We have presented a set of efficient methods for efficiently computing with the BTG model. These methods combine sparse grid quadrature, quadrature sparsification, and tight. quantile bounds significantly reduces the expense of computing BTG's predictive median in certain cases rivaling even the speed of MLE without degrading prediction accuracy. Furthermore, we proposed fast. LOOCV algorithm for BTG for model selection and assessing model fit. An advantage of our framework is that it allows the practitioner to control the trade-off between the speed and accuracy of the Bayesian approach by modulating the sparsification of the grid and tolerance of the quantile-fnding routine.
Using these efficient methods, we were able to compute a set of experiments suggesting that BTG demon- strates superior performance over its comparable models low-data and modium-data regimes, where hy- perparameters likely under-specified.
In future work, would like to combine approach with approximate GP inference to further improve computational efficiency. Furthermore, though our sparse quadrature rules allow us to go to higher dimen- sions, we will still eventually hit a limit (depending on the problem, perhaps around 50); to go to dimensions in the hundreds or thousands, additional assumptions and methodologies will be needed. Lastly, we would like to apply BTG to Bayesian optimization, which is one application we believe to be particularly suitable for it, since we expect data be sparse.
References
Ryan Prescott Adams, Iain Murray, and David JC MacKay. Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities. Proceedings of the 26th Annual International Conference on Machine Learning, pp. 9-16, 2009.
J Aitchinson and JR Dunsmore. Statistical prediction analysis. Technical report, 1975.
17

---

Published in Transactions on Machine Learning Research (04/2023)
George EP Box and David R Cox. An analysis of transformations. Journal of the Royal Statistical Society. Series B (Melhodologienl); 26(2):211-2433 1964.
Hans-Joachim Bungartz and Michael Griebel. Sparse grids. Acta Numerica, 13:147-200 2004.
Henry R Chai and Roman Garnett. Improving quadrature for constrained integrands. In The 22nd Inter- national Conference on Artificial Intelligence and Statistics, pp. 2751-2755 PMLR, 2019.
D Cole, C Martin-Moran, AG Sheard, HKDH Bhadeshia, and DJC MacKay. Modelling creep rupture strength of ferritic steel welds. Science and Technology of Welding and Joining, 5(2):81 89, 2000.
Paulo Cortez, António Cerdeira, Fernando Almeida, Telmo Matos, and José Reis. Modeling wine preferences by data mining from physicochemical properties. Decision Support Systems, 47(4):544-553 2009.
Noel A. C. Cressie. Statislics for spolial dala. Wiley, 1993.
Andreas Damianou and Neil Lawrence. Deep Gaussian processes. In Artificial intelligence and statistics, pp. 207-215. PMLR, 2013.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017.
Mark N Gibbs. Bayesian Gaussian processes for regression and classification. PhD thesis, University of Cambridge, 1998.
Jan Graßhoff, Alexandra Jankowski, and Philipp Rostalski. Scalable gaussian process separation for kernels with non-stationary phase. In International Conference on Machine Learning, 3722-3731. PMLR, 2020.
Florian Heiss and Viktor Winschel. Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics, 144(1):62 - 80, 2008.
Kalvik Jakkala. Deep gaussian processes: A survey. ar preprint rrXiv:2106.12135 2021.
Vidhi Lalchand and Carl Edward Rasmussen. Approximate inference for fully Bayesian Gaussian process regression. In Symposium on Advances in Approzimate Bayesian Inference, pp. 12. PMLR, 2020.
Miguel Lázaro-Gredilla. Bayesian warped Gaussian processes. Advances in Newral Information Processing Systems, 25:1619-1622 1621, 2012.
Patrick Kofod Mogensen and Asbjørn Nilsen Riseth. Optim: A mathematical optimization package for Julia. Journal of Open Source Software, 3(24):615, 2018.
Marcin Molga and Czestaw Smutnicki. Test functions for optimization needs. Test functions for oplimization needs, 101:48, 2005.
Victor De Oliveira, Benjamin Kodem, and David A. Short. Bayesian prediction of transformed Gaussian random fields. Journal of the American Statistical Association, 922440):1422-1433, 1997.
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian processes for machine learning, MIT Press, 2008.
Gonzalo Rios and Felipe Tobar. Composiiionnally-warpp Gaussian processes. Neural Networks, 118:235 246, 2019.
Simone Rossi, Markus Heinonen, Edwin Bonilla, Zheyang Shen, and Maurizio Filippone. Sparse gaussian processes revisited: Bayesian approaches to aaducing-variable approximations. In International Conference on Arlificial Intelligence and Statistics, pp. 1837 1845. PMLR, 2021.
Hugh Salimbeni and Marc Deisenroth. Doubly stochastic variational inference for deep gaussian processes. Advances in neural information processing systems, 30, 2017.
18

---

Published in Transactions on Machine Learning Research (04/2023)
Edward Snelson, Zoubin Ghahramani, and Carl E. Rasmussen. Warped Gaussian processes. Advances in Newral Information Processing Systems, 16:337-344, 2004.
Gunter Spöck, Hannes Kazianka, and Jürgen Pilz. Bayesian trans- Gaussian kriging with log-log transformed skew data. In Interfacing Geostatisties and GTS, pp. 29 29-43. Springer, 2009.
Gilbert W. Stewart. Matriz algorithms,, volume I. SIAM, Society for industrial and applied mathematics, 1998.
Dustin Tran, Rajesh Ranganath, and David M Blei. The variational gaussian process. ar preprint wrXiv:1511.06..0 2015.
Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured gaussian processes (kiss- gp). In International conference machine learning, pp. 1775-1784. PMLR, 2015.
19

---

Published in Transactions on Machine Learning Research (04/2023)
A Methodology
A.1 Quadrature Sparsification
We assume the posterior cdf is the mixture of edfs F(x) - E wxf((z), 0 S f((x) S 1, and fi(a) are M
monotone increasing for = 1.....M. Assume the weights {w, (wx)i=1 are decreasingly ordered by magnitude from 1 to M. We consider the quantiles the approximant (z), a truncated and re-scaled F(x).
Lemma A.1. Define to be the smallest integer such that Σ wa l-c. Then define the scaled, trancated misture
k
Fx(x) = - Σ wif((z), c:=
Σ
i=1
We have
- 2e.
Proof. Let Rx(x)= = [Fa(x) F(x)). We have
&
M
(x) =
Σ
wif((r)
+
Σ wili(x)
E=1
1=++1
M
Σ We w)x((())
+ Σ wxf((x)
c
1 1-c Σ
wi)i(x) -
i=1
$ 2(1 c) S 2c.
Proposition A.1 (Error Bound for Positive Weights). For any e E (0, 1), let k be the smallest integer such
that Σ Wi 2 1 - e. Then define the scaled, truncated misture
Fx(z)
= Σ wr),(x), c:= Σ w.
1=1
1=1
Let p E 10 and assume that p + 2c (0,1). Then we have the bound:
F-1(p 2c) I Fr (p) (p+
Proof.
Let Fa(x")= = p. Then pp- Fx(x")) ' 2e, so
p 2- 2e < F(x") S 2
It follows that
(p- xxx < <F 2c).
Proposition A.2 (Error Bound for Negative Weights). Let F(x:) be defined as before, except each w, is longer required to be positive. Consider the split F(x) = Fas (x)) + Rx)(2)) where Fas((x) = EM' (.(x) and (xr(x) (x) = E Wxf((x). Then for any x, we have Ras (x) E [e_, where the epsilons are defined the sum of posilive (resp. negalive) weights of (xr(x)
M
M
M
[ws]- S 0, 6
[wi]+ 2 O.
$= M1+1
Then we bound Fp) as follows:
20

---

Published in Transactions on Machine Learning Research (04/2023)
Proof. Let Fas (x") p. Then we wrote F(x²) p+ + (xx. Since c- Rx S c+, it follows that
p S F(x²)) +p+
from which the result. follows.
Proposition A.3 (Error Bound at a quantile). Let F(x) be defined as before, 61, 62 E (0,1), and Fx( be an approximate F(x) such that (p - 61) (p) s (D - 62) for some p (0, 1). Assuming P- E1. E2 E (0, 1), we have the following error bound quantile,
dF-1
(D) - -F) S
max (e(p-c).pto)
dx (5)
where = m2)
Proof- We have
FF (p) - (p)
< max maxx(1F-1(2 (1) - (p) FF - -p- -FF
dF
<maffei, 62} max
(5)
(E(e=E1)ptes)
A.2 Quantile Bounds
Proposition A.4 (Convex Hull). Let F(x) be defined as before with w; > 0 and ΣM 1. Then
fin S (p) S max f. (p).
Proof. Assume for the sake of contradiction that F-1(pp) > max, f (p). Let l arg mag, fT (p). Using (p)) 2
for any # k. Therefore, we have the fact that ( 2 (p) for any j * k we have f, 1
p> F(FT(()) (p)
M
Σ
Wj wjfs(fk (p))
j=1
+
Uff (p))
jz
wap +
w/ pp = p:
which leads to contradiction. The lower bound is analogous.
j
Proposition A.5 (Singular Weight). Let F(x) be defined as before with wi 0 and MM w; 1. Let W; W. Then
max f- (p (p- W. S F-1(p) < min f f w.)
ptw:51
Proof. Assume for sake of contradiction that
F-)(p) > 1 w:)
21

---

Published in Transactions on Machine Learning Research (04/2023)
Then
p F (p+
(p+ wi))
N
Σ Wزlي (Ji (p)
wilp + w:) + Σ
زززر (J+ (p+
j&i
2 wi(p + w.)
However this implies that 0 > wi(w; p), which is contradiction because 2 and 2 (by the assumption that + W. < 1). The lower bound is analogous.
B Fast Cross Validation
In this section, we discuss results leading up to O LOOCV algorithms for BTG, which are given by Algorithm | and Algorithm 2 Naively, the BTG LOOCV procedure has O(n") time cost, due to the costs associated with solving generalized least squares problems related by single-point deletion and evaluating determinants of principle submatrices. We present relevant propositions used to solve these LOOCV sub- problems efficiently in $B.1 and derive our full algorithm in $B.2
Notation Let fx, Mx, Kx I, 00. 28.A. me.x and Ca.2 be defined same as in the paper. As before, we use the (-i) notiation to represent to omission of information from the ith data point. For the BTG LOOCV problem, we consider the n submodels [Model(=))p 3== trained on {x(-9) S " li=1 the location- covariate-label triples obtained by omitting data points one at time. We wish to efficiently compute the posterior predictive distributions of all n submodels indexed by i E (1,...n. n),
(f(x(-0))1{) 15 Σ w,L.JJp((,))())))
(9)
where LG - (fx((((((-0))) Aj. 15 and
Jj = -0((==0(00))) 10,,))))) for j E (1,2..... .M.
9 (()x fx) is the probability density function of the t-distribution Tn-p(mo,x. (go.ACo.x) and p(fx 0, A) is the likelihood of data given hyperparameters.
Recall that in Equation
Problem Formulation We have to efficiently compute the parameters that of the posterior mixture of t-distributions in Equation
c(-1) (Parammtters-01 (-0) : ( 2o.. Co,
i=1
For definitions of these quantities, refer to the main text. We instead emphasize here that solving for (--) entails solving perturbed gencralized least squares problems and that solving for m(-i) and C entail solving perturbed linear systems.
For the likelihood term in Equation
g we have
p(fx 0, A) Eo| 112M Mx] 1/2 11/2,111 (-(n-p)/2) 40,1
hence we are interesting in computing the following for i {1, .., n}:
(-t = = 111 1 10)
=
M M
(10)
22

---

Published in Transactions on Machine Learning Research (04/2023)
The perturbed least squares problems and linear systems be solved independently in O(n2) time, hence naive LOOCV procedure would take O(n²) time. However, using matrix decompositions, we can improve this to O(n²) total time.
Algorithms Algorithms and 2 are used for efficiently computing "TParameters-""): and (Det(=)) )im for fixed hyperparameters (8, AJ. The total time complexity is O(n)) because the dominant costs are precomputing Cholesky factorization for kernel matrix and repeating O operations across sub- models.
i=t
Algorithm T-Distributions of Sub-Modei.
Inputs Y = 9A fx, Mx, Kx:
Outputs: {m,_>) (- J
" p (C
f
Piml
Pre-compute R, Rx, and à; where R R = Kx, Rx Rx = MKK Mx, 2 Kx'Y
for i= 1...n do
4 = x xx K-tel
x (-- ---ooodattRRxx. 4x
(Proposition B.3
Y-/||R-Teiill
{(-1) - & -hi
1 (R-Te:) e.a)
(Proposition B.2
(/x0) - (nxx0) ((n
MMX ((())
x(1) y(-0) M1-00,((0 (()))
(-1) + ||r(-0)] |/
48,
Kx
1x - Kax (1'x -1 12x ()))) -I ((1) + (
c + B(x(--))) /Kk (z( ((( --(,)(1(1)) xx(-
0.2
end for
{mb,) in Jizi (ab,> (11) in |ل (0/-0 J M
return
Frozen Hyperparameters We remark that our LOOCV algorithm is possible because sparse grids and QMC are deterministic- since the underlying sampling grids in hyperparameter-spacc are frozen-in in con- trast to Monte Carlo (MC) methods, which are slochastic. Since we use fixed sparse grids, and we are in fact interested in evaluating the posterior distribution at fixed hyper-parameters (0,, X; ): If the sampling grid were not frozen across sub-models, our approach would not viable, because the sampled points in experparameter-spaane would be different for each sub-model. Likewise, in the MLE approach, hyperpa- rameters (0, M should theoretically be retrained on the submodels, hence cannot re-use computed values.
Algorithm 2 Fast Determinant Computation
Inputs Kx
Output (log K (-1) 10) X -
Precompute RD R Kx
Precompute log(IKxl)
5: for i 1...n n do
6: b, = e k
7: log | +- lol([Xxl) + log(b.) (Propsition B.1
end for return
B.1 Auxiliary Results
In this section, present linear algebra results used in the derivations of Algorithms and 2 in
B.2
23

---

Published in Transactions on Machine Learning Research (04/2023)
Proposition B.1 (Determinant of Principal Minor).
det (2(- del(E)) es)
Proposition B.2 (Abridged Linear System). Let E Thnxn be full rank, and let e, y E R satisfy Ke = y. Then r;;; = c//e K ei,,, we have:
d = (K(- rrKK
Lemma B.1 (Determinant of the Schur Complement of a Principal Minor). If X E Raxm with m < n has full column rank and RRxx a positive definite matria, then
det (-!-)))
det (Ii-1)) det (xx-
Σ
X
0
T
Σ
X
6
X
C2
Proof. Extend the Cholesky factorization R2
R11 R1
of E to obtain the LDL-decomposition
Σ X
x
R R
R1
R12
W
0
R
1
0
Ree
where Ryz Cholosky(RRT... 212) and R12 RTj X. Observe (((-)))((X111)) -1 (x(-) is Schur complement of W(-O), This implies that.
det (W(-1))
det (2(-1)) det (-(x(10))² ((((0)))(((-))) -1 -x((11))
By Proposition B.1
det (W(-0) det(W)), (W W-tes
Therefore
del ((x(-1))² (I(-1)) -1 -xx(-1))
1
dct(I((i)) det (W) (W)?? W -lee.
Lemma B.2 (Rank one downdate for bilinear forms). If R and E is a positive definite matriz in IR then
(x(- (D(())) (-0) = of 0 (x-1 Σ- -- D-1
Nー x
where E(-i) R(n=1)x(x-1) is the (i, i)th principal minor of E and x(-1) R. (n-1) results from deleting the ith enlry of x.
Proof. By Lemma B.1 we have
T (2(()) (-1)) = det ((x(-0)) ((01-0)((E(())) (I(- -1 (1(0))
det(I((1)) det(W)e W-te
24

---

Published in Transactions on Machine Learning Research (04/2023)
In this equation,
Σ
X
1
0 -1
Ru
R12 R2
W -
RºT
R,
R =
X
0
0
O
where Ru1 chol(2)), R12 Ríi X, and Ro22 term of eg. Since,
VIT
Using this decomposition, we may compute the
x e
R TT 6 =
RRT ei
Vx Σ x
we have
E (TT'T))(a tei
c @ e, = R RR-- RR-Te E
x 2-'2-x
xx I e,)
= "ー
x E
1
Lastly, we have
det(W) del(E) del del( 5-1 det (e)()
---x
e2t) 2 -i-te;
e'''i-te
These together imply that
-x(-0))))( (1(-1) = x-x--x
x TT
- T--1 I-1e;
Cj
e E, ---ee
as desired.
Proposition B.3 (Rank one matrix downdate). If X E RR with m has full column rank and Σ is positive definite matrix in Rax*, Let U, = E-'e, then
( -x" (X-1-
Viad o xx-les
X,
where Σ(-i) E ((-=1)x(x-1) is the (i, i)th principal minor of E and ((-i) E ((n- 1)xm results from deleting row from X.
E-leef It suffices to prove that
Proof. Let E:= E-1-
d
Vz,yeR
However, this follows from Lemma
B.2
because
((x + (-)) (x+t))) (x + v) TE(++ +w)
Expanding and canceling symmetric terms yields
(at (a (11-0)(((--1))) (1)--))-')((- + (((-0)((2((110) -1 (t))
x T Ey + +y Tix
implying the result.
25

---

Published in Transactions on Machine Learning Research (04/2023)
B.2 Algorithm Derivation
Recall the following definitions of elements in (TParameters(-")) ) from B
9 = min 11x (fx) Mx 6 Mx (KX)
(11)
B argming (fxx Mx (xx)* (Kx))
(12)
mA,0 = Kax x (fx( Mx + m)))
(13)
Cx =
m²] (Schur Complement)
(14)
We use the generalized least squares LOOCCV subroutine, outlined in Section 1.2.1 to compute % and
3 efficiently for all E n... We use Proposition B.2 to efficiently compute mx,0 and Cx whenever perturbed linear system arises. Generally, these routines involve precomputing a Cholesky decomposition and using it for back-substitution. These steps are enumerated in Algorithm [
The computation (Det(-0) Jm straightforward given Proposition B.1 and a Cholesky decomposition of the kernel matrix.
B.2.1 Generalized Least Squares
The gencralized least squares (GLS) LOOCV problem is that of solving the following set problems cffi- ciently:
n
arg min
- A(-0 D(-1) 1=1
ZERP
It is assumed that Σ E Raxx is positive definite, b E R", A E Raxp, Rank(A) = Rank(A(-1)) = p for some
< and for all i E m1.
We consider the normal equations for the ith subproblem:
arg min (A(-)) - -(-1) (Σ(-1)) (((-)) - (x(-)))-(-))) #ER
namely,
(A(-0)) ( ((-0)
(EX (xx)))-(---
(15)
We first show that Equation has a unique solution. By Proposition
B.3 we have
((--))((((1))) (((())))(((()))=)-)) - A
Σ-le;"
where u = AT The LHS a rank-1 downdate applicd to -T A. Morcover, the LHS is positive definite and hence invertible, because by assumption, Rank(A(-i) p, and E(-i is positive definite.
We find the solution to Equation 15 by first computing the Cholesky factorization of the LHS. Specifically, given a Cholesky factorization Rd R of - from the full problem, the Cholesky factorization of the subproblem can be computed by O(p²) Cholesky downdate R(-1) Downdate(R, u.je, Σ --'e,) such that
(R( T
" R(-))
(A(-1)) = -A
--0
We therefore can solve the normal equation 15
(-0) (a(-0)) a(-0) (A(-
where where y( We first evaluate (i) in (n time by making use of Proposition B.2 We then perform back-substitution
The cost of O(n2) is attained by evaluating terms from right to left.
using the cholesky factor -() in O(p)) time. The overall time complexity is thus O(n3).
26

---

Published in Transactions on Machine Learning Research (04/2023)
C Supplementary Numerical Results
In addition to the main results in Table Conclusions are similar as in Section 5
we provide results of the MAE metrics as well in Table 6
IntSine
Camel
Abalone
Wine
Creep
RMSE
MAE
RMSE
MAE
RMSE
MAE
RMSE
MAE
RMSE
MAE
GP
0.227
0.171
2.003
1.781
3.290
2.208
1.994
1.792
37.88
25.68
WGP-A
0.179
0.117
2.055
1.815
3.097
2.068
0.811
0.677
35.48
24.27
WGP-SA
0.172
0.109
2.012
1.788
2.992
2.022
0.809
0.670
61.62
45.85
WGP-BC
0.184
0.123
1.964
1.751
2.826
1.940
1.045
0.784
40.38
25.49
CWGP-L-SA
0.172
0.109
2.055
1.815
3.117
2.100
0.808
0.670
91.89
73.05
CWGP-A-BC
0.174
0.113
1.960
1.747
3.088
2.069
0.808
0.670
37.89
25.34
BTG-I
0.169
0.100
1.820
1.731
2.804
1.842
0.808
0.670
38.11
26.18
BTG-A
0.168
0.101
1.827
1.741
2.890
1.900
0.807
0.669
38.68
27.68
BTG-SA
0.165
0.104
1.801
1.796
2.791
1.822
0.820
0.696
91.69
74.04
BTG-BC
0.170
0.102
1.675
1.666
3.225
2.172
0.808
0.670
39.10
26.52
BTG-L-SA
0.145
0.082
1.673
1.658
2.871
1.870
0.809
0.670
91.83
73.05
BTG-A-BC
0.159
0.090
1.828
1.742
2.832
1.814
0.802
0.664
35.25
24.44
Table 6: A comprehensive set of RMSE and MAE results for prediction using the IntSine, SixHumpCamel, Abalone, Wine and Creep datasets. We generally kept the training sets small, to highlight BTG's advantages in the data-sparse regime, though for the IntSine function there was a good amount of data. We compare the following models: GP, WGP, CWGP, BTG, and use the following transformations transformations: T: Identity, A: AreSinh, SA: SinhAreSinh, BC: BoxCox, L: affine. The best method is bolded; find that BTG has lower RMSE and MAE than the other models, though unsurprisingly, different transformations are more suitable for different datasels.
27


Reviews:
[{"review_id": "IoyiQOvRbw", "reviewer": "Reviewer_5sB8", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- This paper presents relevant literature and content in a way that is easier to follow and sets up a clear motivation.\n- This work presents several theoretical and practical results that are interesting. E.g., analysis to quantify the quantile estimation error by sparsifying quadrature rules is interesting.  \n- Empirical analyses highlighting proposed algorithms' computational and speed aspects are convincing.  \n\nWeaknesses:\n- Time and accuracy compared with the MLE approach would be helpful. A top example demonstrates an 'overconfidence problem with the MLE approach, but an experimental setup involving a real dataset would add more value. \n- It seems worth mentioning about GraBhoff et al. 2020 and how the author's work compares to theirs. GraBhoff et al. 2020 adopt SKI to WGP setting and is fast low dimensional settings similar to the proposed work and sparse grids. \n- Simone et al. also propose a fully Bayesian sparse GPs approach that is relatively more scalable. It is worth mentioning why this approach would not apply to wrapped input settings considered in this paper. Or if it is, then it may be compared with them.\n- It'd be nice to mention when these approaches are more feasible and when they won't be possible. E.g., sparse grids are typically infeasible for high dimensions. Does the sparsification of quadrature allow it to go to higher dimensions? \n\nGraBhoff et al. 2020:- Jan GraBhoff, Alexandra Jankowski, Philipp Rostalski, \"Scalable {G}aussian Process Separation for Kernels with a Non-Stationary Phase\" Proceedings of the 37th International Conference on Machine Learning, PMLR 119:3722-3731, 2020.\n\nRossi, Simone, et al. \"Sparse Gaussian processes revisited: Bayesian approaches to inducing-variable approximations.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2021.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work revisits the BTG model, which has latent and observations spaces with a deterministic function mapping these two spaces. The central idea is to approximate the BTG posterior predictive density using a mixture of t-distributions. In doing so, several tricks are cleverly intermixed to address computational challenges in obtaining the GP posterior. The proposed work is evaluated in various relevant experimental settings to ascertain the advantages of introduced techniques, i.e., sparse grids, QMC, and fast way to perform LOOCV. ", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "One of the most critical changes is demonstrating that the proposed approach can indeed yield better uncertainty estimates on a real dataset. \n\nFor the remaining changes, please refer to the Weaknesses section for the requested change. Consider points with which you agree. \n\nMinor typo: 'quartic' in \"which incurs quartic cost naively.\" \n ", "broader_impact_concerns": "I don't see it as a concern. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "cnoF8iyGC4", "reviewer": "Reviewer_6uRk", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I have found the text well-written and relatively easy to follow. The motivation for the proposed methodology is clear, and the brief literature review is sufficient. I have not checked every single computation line-by-line, but the derivations are relatively standard. The numerical experiments are well executed and convincing. I think that it would be interesting to understand in more details the regime when the Bayesian treatment (BTG) is actually worth implementing when compared to the point-estimate approach (WGP/MLE). It is clear from the simulation (and the discussion in the text) that the data-scarce regime is such a scenarios -- what does happen with larger datasets? \n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The article proposes a computationally efficient approximation to the Bayesian Transformed Gaussian (BTG) model (i.e. composition of a GP with a deterministic function). The proposed method approximates the Bayesian posterior by a mixture of T-distributions, leverages sparse-quadrature methods (as well as a sparsification trick) instead of naive Monte-Carlo approximations, as well as clever numerical Linear-Algebra tricks for speeding-up computations. The proposed pipeline is tested on a range of test datasets and appears to provide competitive predictive performances without sacrificing scalability. \n\nBefore gathering a few comments, I would like to point out that I am **not** an expert in GP modeling.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "\nI have a few questions: these are not necessarily requests of changes, but questions that I have asked myself while reading the manuscript -- this also indicates that some of the points below are now crystal clear in the current version of the manuscript.\n\n1. The choice of prior distributions appears quite arbitrary (ie. chosen to ease computation, mainly). Is this correct? Does it impact the performance? If scalability were not an issue, are there better choices to consider?\n2. In Equation (7) for $\\phi(f|\\theta, \\lambda, \\ldots)$ , am I correct assuming that the integral with respect to $f(x)$ does not equal one, because of the term $|g'_{\\lambda}(f(x))|$ ?\n3. In Equation $(8)$ do the weights $\\widetilde{w}_i$ sum-up to one? I am finding this confusing.\n4. Why use both MSE and RMSE?\n5. I am not an expert in GP -- I could not help myself wonder whether choosing a more flexible family of kernel function would be enough to avoid using any \"wrapping\" function? In other words, would this lead to roughly equivalent predictive performances? All the experiments are performed with the extremely simple RBF kernel.\n6. Performing a simple Gaussian approximation (eg. Laplace approx) of the posterior is a possible approach in between MLE and the proposed approach -- how is this approach performing, numerically?\n7. One a possibly very simple example, this would be interesting, I think, to investigate WGP vs BTG as the number of points increases -- this would help to understand the regime when the \"full Bayesian treatment\" starts to not be necessary.\n\nI have enjoyed reading the text and I am looking forward to reading the authors' comments and clarifications.\n", "broader_impact_concerns": "NA.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "cQP62K9Lzt", "reviewer": "Reviewer_xyvo", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nThe paper is very well-written with clear & well-organized presentation. Despite being highly technical, the paper is actually easy to follow.\nThe theoretical results are interesting & refreshing to read.\n\nWeaknesses:\n\nDespite the above strengths, there are several aspects of this paper that are unclear or even debatable for me:\n\nFirst, I do not agree with the viewpoint that variational inference on deep GP is not considered a fully Bayesian treatment. To me, full Bayesian treatment means we have uncertainty calibration of the model parameters; and variational inference (VI) does provide that. In fact, exact full Bayesian treatment is rarely possible and mostly people would resort to approximation schemes which are either via VI or sampling, such as the MC approaches adopted by the authors. \n\nThus, in this line of thought, I would consider deep GP as a variant of multi-layer warped GP with full Bayesian treatment. As pointed out by the authors, deep GP provides a multi-layered transformation in which each layer is a random function distributed by a GP prior. On high level, deep GP and the revisited BTG basically aims to solve the same modeling problem with the same idea on provided multi-layered warping except that one is non-parametric while the other is parametric. \n\nGiven this, I believe the experiment needs comparison with deep GP baselines to demonstrate more convincingly the advantages of the revisited BTG over deep GP. \n\nSecond, regarding the improved scalability of BTG, I think the authors might have missed comparison with a set of potentially competitive baselines that involves the use of Gaussian processes. To elaborate, while warped Gaussian processes had been mostly demonstrated for full Gaussian processes, incorporating such warping idea to sparse Gaussian processes is straightforward and will likely lead to significant improvement in the compute cost. \n\nTo see this, note that the key to warped Gaussian processes is a simple application of the change-of-variable theorem to the marginal likelihood, which results in an augmentation of the original likelihood with an additional log det of the corresponding Jacobian. Most sparse GP (which scales linearly in the no. of observations) has an analytic marginal likelihood and hence, incorporating the warping idea is immediately possible. I suspect this will likely reduce the cost to O(Mn) where M is the no. of quadrature nodes.\n\nHere are a few samples in this direction:\n\nSparse online warped Gaussian process for wind power probabilistic forecasting (Applied Energy, 2013, pages 410-428)\nSparse Spectrum Warped Input Measures for Nonstationary Kernel Learning (NeurIPS-20)\n\nAlternatively, we can also pick any off-the-shelf sparse GPs [*] and apply the change-of-variable theorem on their marginal likelihood to derive fast variant of warped Gaussian processes. \n\n[*] https://www.jmlr.org/papers/volume6/quinonero-candela05a/quinonero-candela05a.pdf\n\nLast, the entire point of being fully Bayesian is to have good uncertainty calibration. But this aspect has not been demonstrated at all. Furthermore, the claim that the proposed method can scale well in high-dimensional data space is a not very well-demonstrated seeing that most experimental datasets have low dimensions.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper revisited the Bayesian Transformed Gaussian process (BTG) model, which is summarized below:\n\nThe original BTG (Kedem & Oliviera, 1997) imposes a Gaussian process (GP) prior on a parametric transformation of a (non-parametric) function of interest. This helps relaxing a strong assumption of Gaussian process which stipulates that the observation distribution is a multivariate normal distribution. This is similar to the idea of warped Gaussian processes (WGP) by (Snelson et al., 2004) but unlike WGP which optimizes the transformation/warped function via MLE; the BTG model further imposes a parametric prior on the warping function, which enables analytic derivation of the model evidence.\n\nHowever, the BTG model needs to use an artificial (and uninformative) prior to enable such analytic full Bayesian treatment, which is also costly in spaces of high-dimensional data. The BTG revisitation proposed in this paper therefore focuses on developing a more practical Bayesian treatment to reduce the overall computation cost of BTG while retaining its full Bayesian flavor. The BTG is also generalized to include multi-layered warping. Efficient methods for computing BTG predictive medians and quantiles as well as fast LOOCV with BTG are also proposed. The results are positively demonstrated on 6 regression datasets with input dim ranging from 5 to 18.\n", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. Comparison with Deep GP models as well as the warping variant of sparse Gaussian processes as elaborated above. \n2. Elaborate (at least empirically) on how the revisited BTG would have an advantage over deep GP as well as the (warped) sparse GP baselines in terms of uncertainty calibration.\n3. If possible, do consider providing extra experiments with higher-dimensional (e.g., d = 100) datasets.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "0sWDNyiBCe", "reviewer": "Reviewer_bVfk", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "*Strengths*\n- This paper is well written in general and easy to follow, although a few sections would require some extent of rewritten and clarification.\n- The introduced techniques to speed up BTG computation are sound and well explained. Using quantile bounds to speed up the search is a nice idea.\n- The theoretical analysis to quantify the quantile estimation error caused by sparsification is useful in practice.\n- The empirical analysis of each component of the algorithm is useful for readers to appreciate their contribution to the entire method.\n- The Julia package could be a useful tool for people to use more sophisticated and robust GP algorithms than a vanilla implementation.\n\n\n*Weaknesses*\n- Limited novelty. The main contribution is to speed up a BTG model proposed in 1997 using a series of techniques. Quadrature is a widely used method to speed Monte Carlo estimation. It’s hard to justify that revisiting BTG in the context of ML and GP literature as a contribution compared to the original paper that proposed BTG as a Bayesian version for Gaussian random fields.\n\n- Some claims about the efficiency of the proposed method are not well supported by the experiment. \n  - In Figure 3, to obtain a sufficient weight coverage, e.g. 0.9 or 0.95, the required number of nodes of QMC is similar to sparse grid quadrature. How do the authors choose the sparsification level in the experiments? Also, the prediction time is shown in the log-log scale. One should show it in the linear scale to tell whether it’s linear or estimate the exponent.\n  - In Figure 4, the authors claim “It is verified that the naive way has a O(n^4) cost while our … to O(n^3)” but the figure clearly shows that Naive LOOC scales better than n^4 while the prosed is worse than n^3. Can you estimate the exponent of the curves and present them in the paper?\n  - Please list the estimation error in Table 4 to check whether the differences are statistically significant.\n\n- Lack of comparison with BWGP\n\nThe closest method to the BTG is Bayesian WGP. The authors only compare their method with vanilla GP and WGP with MLE but what they really need to compare with is BWGP that also takes a Bayesian approach to estimate the warping parameters.\n\n- Not sufficient empirical verification.\n\nMany of the conclusions are drawn based on a single (synthetic) example, e.g. the smoothness of the $\\phi$ function in Figure 2 to justify the use of quadrature method; bounds speed up in Table 3; comparison of weight concentration in Figure 3. Also there are plenty of benchmarks for function regression, both synthetic and real ML applications. The empirical experiments on the five datasets in Table 4 (2 synthetic and 3 real) is probably not convincing enough.\n\n- There are several places that require some rewriting.\n  - The introduction of BTG in section 3 is not clear enough. A simple explanation for “trans-kriging model” would be very helpful to understand how the presentation of BTG here is different from that in Oliveira et al. (1997).\n  - Section 3.1, the notation n, p, and prior of parameters in the marginal probability are used without definition. Is the notation \\eta a typo? Statements like “BTG adopts the improper joint prior” and “BTG then adopts a conditional normal-inverse-gamma prior on (\\beta, \\tau)” are contradicting with each other.\n  - “input-waped” GPs has another meaning in the literature (Snoek et al., 2014) that warps the inputs x while this paper warps the outputs y.\n  - Figure 1 is not a good example to show the benefit of BTG compared to GP as “a fully Bayesian model” vs “its MLE counterpart”. To show the advantage against WGP, I’d expect a comparison with WGP instead of GP. Also, it looks like the length scale of Figure 1 is not well fitted. I guess the GP prediction will be better with a smaller length scale.\n  - In section 2.3 and Table 1, “the Bayesian Warped GP is not a fully Bayesian GP, and instead approximates the posterior distribution using variational inference”. I would not say that BWGP is not fully Bayesian while the proposed method is. Neither computes the exact posterior accurately but takes a different approximation, variational inference vs Monte Carlo.\n  - Conclusion: “we introduced the Bayesian transformed Gaussian (BTG) model”. I would not say so.\n\nRefs:\nSnoek, J., Swersky, K., Zemel, R. and Adams, R., 2014, June. Input warping for Bayesian optimization of non-stationary functions. In International Conference on Machine Learning (pp. 1674-1682). PMLR.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper revisits the Bayesian transformed Gaussian introduced in Oliveira et al. (1997) and presents it as a Bayesian counterpart of the Warp GP (with potentially multiple layers of transformation) by assigning a prior to the warping parameters. The main contribution of the paper is to introduce a few techniques to speed up the computational cost of BTG inference, prediction and selection, including:\n- Replace Monte Carlo estimation with sparse grid quadrature or Quasi-Monte Carlo. Additionally sparsify the samples to trade-off approximation error with speed\n- Use quantile bounds to speed up quantile estimation.\n- Fast implementation of Leave-One-Out-Cross-Validation that reduces the computation complexity from N^4 to N^3.\nThe paper also conducts empirical comparison of BTG with multiple transformation layers with GP and WGP methods, and provides a Julia package for transformed GPs.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please refer to Weaknesses section for requested change. The main points:\n- Comparison with BWGP\n- Experiments on additional benchmarks. For ablation studies, it’s fine to use one example for demonstration in the main text but it would be very useful to have additional examples in the appendix. For the main results in Table 4, I would expect more test functions.\n- Rewriting part of the text as pointed out above.", "broader_impact_concerns": "None. This paper introduces a few techniques to speed up the computation of an existing model.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This submission makes an important contribution to a remarkably flexible model class: the Bayesian transformed Gaussian (BTG) model. Although quite flexible and useful for modeling, the model presents numerous computational challenges in a practical setting. The authors outline best practices for working with this model -- introducing some useful computational "tricks" along the way -- resulting in an effective scheme for working with the BTG model in practice.

These algorithmic contributions are coupled with a series of well-designed empirical investigations offering additional insight into the BTG model.

I want to commend the authors in particular for their efforts in engaging the reviewers during the discussion period and incorporating numerous edits into their draft accordingly, which I believe strengthened the paper considerably.
Decision: Accept (poster)