Bezier Gaussian Processes for Tall and Wide Data
Authors: Martin Jørgensen
Michael A Osborne
Venue: NeurIPS 2022 Accept (2022)
Keywords: Gaussian processes
high-dimensional
scalable
kernels

Abstract: Modern approximations to Gaussian processes are suitable for ``tall data'', with a cost that scales well in the number of observations, but under-performs on ``wide data'', scaling poorly in the number of input features. That is, as the number of input features grows, good predictive performance requires the number of summarising variables, and their associated cost, to grow rapidly. We introduce a kernel that allows the number of summarising variables to grow exponentially with the number of input features, but requires only linear cost in both number of observations and input features. This scaling is achieved through our introduction of the ``Bezier buttress'', which allows approximate inference without computing matrix inverses or determinants. We show that our kernel has close similarities to some of the most used kernels in Gaussian process regression, and empirically demonstrate the kernel's ability to scale to both tall and wide datasets.

Bézier Gaussian Processes for Tall and Wide Data
Martin Jørgensen Department of Engincering Science University of Oxford martinj@robots. ox. ac. uk
Michael A. Osborne Department of Engineering Science
University of Oxford mosbdrobots ox. ac uk
Abstract
Modern approximations to Gaussian processes are suitable for "tall data", with a cost that scales well in the number of observations. but under-performs data", scaling poorly in the number of input features. That is, as the number of input features grows, good predictive performance requires the number of summarising variables, and their associated cost, to grow rapidly. We introduce a kernel that allows the number of summarising variables grow exponentially with the number of input features, but requires only linear cost in both number of observations and input features. This scaling is achieved through our introduction of the Bézier buttress, which allows approximate inference without computing matrix inverses or determinants. We show that our kernel has close similarities to some the most used kernels in Gaussian process regression, and empirically demonstrate the kernel's ability to scale to both tall and wide datasets.
"wide
Gaussian processes (GPs) probabilistic approach to modelling functions that permit tractable Bayesian inference. They are, however, notorious for their poor scalability. In recent decades, this criticism has been challenged. Several approximate methods now allow GPs scale millions of data points, Yet, scalability in the number of data points is merely one challenge of big data. There are still problems associated with the inpul dimensionality one aspect of the famed curse of dimensionality. Burrt et al. (2200) analysed the most studied approximation, the so-called sparse inducing points methods, and showed it to be accurate for low dimensional inputs. Alarmingly, exponentially many inducing points are still needed in gigh-dimensional input spaces, that is, for problems with a large number features. As such, despite modern GP approximations scaling to tall data, they are still discounted when concerning wide data.
In response this, there exist GP approximations built simplices grid structtrsssurss in the input space [Wilson and Nickisch 2015 Gardner et al., 2018 Kapoor et al. [021]. These take advantage of attractive fast linear algebra, but are often limited by memory in higher dimensions. Their advantage is the ability to fill the input space with structured points, so all observations have close neighbour.
We propose
kernel for GP regression that requires neither matrix inversion determinant calculation GPs' two core computational sinners. Additionally, we cover the input spacd") with exponentially many points, but introduce approximation that grows only linearly in computational complexity. That is, our method scales linearly both the number of data points and the number of input dimensions, whilst being ppace-filling in the input domain.
GPs are indispensable to fields where uncertainty is driver in decision Docision-making mechanisms. Such fields include Bayesian optimisation, active learning and reinforcement learning. The critical decision trade-off. One ability useful such fields to assign high uncertainty unexplored regions, just as does an exact GP. We show that our proposed model also assigns high uncertainty to unexplored regions, suggesting our model as well-suited to decision- making problems.
mechanism the ppploaatio-exxplotattie
A limiting assumption of our kernel is its restriction to box nox-bounded domain in the input space.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

3
3
Control point
Points in hypercube (input space)
Figure I: Left: Illustration of 3d control points (orange) and the 2d grid spanning the input space [0, x [0, 1]. We see how the function (or surface) f gets impacted by the control points. Right:
How the points cover the hypercube in 3 dimensions. Here grid in 3d with orders 32)..).
Background
Bézier curves and surfaces are parametrised geometric objects that have found great usage in computer-aided aided design and robotics (Prautzsch et 2002]. The simplest Bézier curve the linear interpolation of two points Po and in R"; the Bézier curve of order
c(t) = (1 t)p. + (p1, LE [0.
Higher order curves are generalised in the following way: the order-v Bézier curve is defined
c(t) = Σ (t)t),,
E [0.
(2)
1=0
In Bézier terms, p; are referred to as control points. Notice an order-v Bézier curve has + control points. B, denotes the ith Bernstein polynomial of order v. They are defined as
v!
(x) =!(v - i)!
(1-t))
(3)
By going from curves to surfaces we wish to extend from the scalar t to spatial input æ E [0, 1] for d 1. Here, we can define Bézier d-surfaces
ca(x) = - ..
BB (F)) (Fa)Pe Σ B (01)
z (x1, T2........... ........,.y,, Fz) E [0, 1]d.
(4)
S1 1==0
12
1==0
Figure gives visual illustration of 2 2-dimensional surface embedded in R3. In the literature, it is difficult to find any studies of d-surfaces for > 2. This paper targets especially this high-dimensional input case. We restrict our output dimension tend to multidimensional outputs. The red points in Figure show how each control points has an asso- ciated location the input space. They are placed on a grid-like like structure, and the order of each dimen- sion determines how fine the mesh-grid of the hypercube is; i.c. how dense the input-space is filled.
the regression problem, but the methods naturally ex-
Gaussian processes (GPs) meticulously studied in probability and statistics (Williams and Rasmussen 2006). They provide a way to define a probability distribution over functions. This makes them useful, priors, to build structures for quantifying uncertainty in prediction. They defined as a probability measure over functions f X' + R, such that any collection elements following a joint Gaussian distribution. This distribution fully determined by a mean function m : X + IR and positive semi-definite kernel function A: X' X + R. They admit exact Bayesian inference. However, exact inference comes at a prohibitive worst-case computational cost of On""), where is the number of training points, due computing inverse and determinant of the kernel matrix.
(x), 322 ... x..) in X' have their associated output (J( ... f(xx)) ((x)).
Sparse Gaussian processes [Snelson and Ghahramani 2005 overcome this burden by conditioning inducing points, reducing complexity to O(nm²) where usually < The inducing points, denoted u, are then marginalised to obtain an approximate posterior of f. The variational posterior mean and variance, at a location x* are then given by
[] = k( (Z)k((Z, -) Mu.
(5)
Var(f(=*)) = k<==x=))) k(x" Z) Z. Z) (k(Z, Z) - kZ. 1k( x²),
(6)

---

under the assumption of a constant zero prior mean function. This assumption is easily relaxed if needed. Here Z denotes the inducing locations in the input space, i.e. f(Z) u NH Σ, ). Under further assumption of Gaussian observation noise e, i.e. y² = f((=') + 6, then Titsias 2009 showed the optimal μ and Σn are known analytically. a sought analogy to Figure 1, Z would he the red points and u would be the orange.
2 Bézier Gaussian Processes
Inspired by Bézier surfaces, we construct a Gaussian process f [0, 1j4 + R as
61
Vi
Vd
[(x)
Σ
B
(x1)....B B. (xa) Pi Piy....ka b1
=1 12=
where P Piy.ga....ds 12..1.2....ii ig..i..... are Gaussian variables and æ (x), 22,-... r.).
Here, Ty E [0, y lI...,,d. d. We write P. with capital letter to emphasise that is a random variable now. Further, we write it in boldface though we here only consider the scalar case, i.e. regression, but the multi-output case is not fundamentally different. is easy to verify that f satisfies the definition of GP since it, for any æ, is scaled sum of Gaussians. We assume that all P. i .......a are fully independent. With that assumption, we can make the following observation for the mean and kernel function
2/
p(x) =
Σ
B (2))
(xa) (n.)..i it....aaa and
(8)
12
k2=0
11
Vy
k(0, 2) = Cov (f(z),
Σ 11 +==0
B (21)
BB Er B (21). BBy (za). (9)
The Bernstein polynomials can approximate any continuous function given the order is large enough, thus they make good basis for GP regression (Hildebrandt and Schoenberg 1933). Naturally, selecting a prior over f comes down to selecting a prior over the random control points P. The most common prior mean in GP regression, the constant zero function, is then easily obtained by 0 01 = 0 for all i. By construction, the choice of Eig....fu needs consideration to yield a convenient prior over f. Mindlessly setting Ei.....is 1 would make Var (fx)) collapse quickly in the central region of the domain, especially as dimensions d grow. This, of course, gives much too narrow prior over f.
Figure 2 (middle) shows that in the central region the standard deviation of f is smaller due to the nature of the Bernstein polynomials. If we consider instead wwo-dimensional input the standard deviation would collapse even more, as we would then see the shrinking effect for both dimension and multiply them. We can, however, adjust for this.
We define the inverse squared Bernstein adjusted prior to counter this effect. all dimensions d, let
5y - A, A,'11,,41 where Ar,- (1/w)))
(10)
and By denotes the order of the dimension 7. Then setting E, Dij...he = II (ig)) ensures that Var (f()))
over the entire domain [0, 1jd. Eq. 10 solves a linear system, such that Var (f(i/e))) = 1, for = 0, O...., My. This means a prior hardly distinguishable from standard stationary ones such as the RBF kernel. Visual representation of this prior is shown in Figure[(rrghh This adjustment works up to By 25, after which negative values occur.
Summarising, we introduced kernel based on Bézier surfaces. An alternative viewpoint is that is polynomial GP, but with Bernstein basis rather than the canonical basis. We remark that f is defined outside the domain [0, 1]d, any intuition about the prior there is not considered, and we will not pursue investigating data points outside this domain. Of course, for practical purposes this domain generalises without loss of generality to any rectangular domain [a1, b1] x...x (aa, ba). For presentation purposes we keep writing [0, 1]d Next, we show how infer approximate posterior given data.
3

---

Bemslein
Unaduuted
Adjusted
1.0
1.0
1.0
0.
0.5
0.5
0.4
02
0.6 0.0
0.6
Figure 2: Left: The 21 Bernstein polynomials that make up the basis of a Bézier GP of order 20. We observe how they each impact 'local' region along the [0,1] domain. Middle: If not accounting for the Bernstein polynomials behaviour, the variance is of f is more narrow in the central region. Right: can enforce uniform variance the domain.
By adjusting, Eq. T0
2.1 Variational inference
Let P denote the set of all Pi Pi As the prior of f is fully determined by the random control points P. the posterior of f is determined by the posterior of these. As per above, set the prior
p(P)
II
11 +==0
P
II N(0, Ei........d.)........ 1==0===0 1==0
(11)
b==0
where Eg Di....ii II (1x)). We utilise variational inference to approximate the posterior of the control points, and hence f. This means we introduce variational control points. We assume they are fully independent (usually called the mean-field assumption), and have free parameters for the mean and variance, such that P......ka ~
Ei..fa.....da
Assume we have observed data D = {æy, 8y 1- The key quantity in variational inference is the Kullback-Leible divergence between the true posterior p(Ply) and the variational approximation which denote q(P). The smaller divergence, the better approximation of the true posterior. Without access to the true posterior, the quantity is not computable. However, it has been shown this divergence is equal to the slack in Jensen's inequality used of the log- marginal likelihood: log p(y).
logp(y) log / ply)Pip(P)dP 2 f log P)p()P)((P)
(P)dP
(12)
((P)
By [log(yy P)] KL. (I()))lp(())))
(13)
Knowing this, can approximate the true posterior with (P) by maximising Eq. 13 This the evidence lower bound, and it is maximised with respect the variational parameters 0 ig.iig.... and E, Eg........... This is fully analytical when the variational parameters and a Gaussian likelihood is assumed.
We assume our observation model is disturbed with additive Gaussian noise, which in other words means our likelihood is Gaussian
plu, P) N' (z)
(14)
for each j = 1...., and we assume they are independent conditioned on P. With these assumption the first term in Eq. (13) becomes
log(2r)) + log(o2)+ + (Ns - ((x))) Varg(e) (f(z))
Ep(P) [logp(y] P)]
2
o²
(15)
where [fx;)) and Var, (f(z)) are given Eq. variational parameters & da and Σ fa...kk used.
and
9 respectively, but with the

---

300,0.2
W2.11
W3,1
3,3332
Figure 3: A Bézier buttress visualised. Here, the input dimension hence 3 layers. Each layer has 4 nodes, since each dimension has order 3. There exist 43 paths from source (left square) to sink (right square), each path represents one unique control point. We can sum over all control points by sequential matrix multiplication from source to sink.
The second term in Eq.
13 enjoys the independence of control points to split into sums
P1
P2
KL (L())))((())
KL (
(16)
02=0
11
P2
Vd
Σ Eiy
1 Dis 1 1
En Eny
Σ
+ log
(17)
E Eiy
When inspecting the evidence lower bound, Eq. 13 we see it has data term forcing control points to fit the data, and a KL-term to make control points revert to the prior. Knowing how control points are allocated in the input domain, we expect control points in regions of no data revert to the prior. This is similar to what stationary kemells do in said regions. We verify visually in Section 4
All together, Bézier GPs can be adjusted have priors similar stationary GPs, and have analogous posterior behaviour, which is favourable to many practitioners. But Bézier GPs scale. None of the terms in the evidence lower bound require matrix inversions or determinants. It is simple to mini-batch batch over the data points, utilising stochastic variational inference [Hoffman et 2013). to scale to large n. However, nearly all terms require evaluations of huge sums the input dimension is high. The next section is aimed at this problem.
2.2 Scalability with the Bézier buttress
Until this point, we have omitted addressing the number random control points needed for Bézier GPs. Let us denote this number 7. It can quickly be checked that - = II (My + 1). This implies that to evaluate f we must sum over 7 summands, which as d increases, quickly becomes computationally cumbersome. increases exponentially with d. It is justifiable view the random control points as inducing points; after all, they are Gaussian variables in the output space, Thus, it would be extremely valuable manage exponentially many of them.
To overcome this, we introduce the Bézier buttres.
We assume parameters of the random control
points, say 9, can parametrise U4..4......4 tion is the key of the Bézier buttress. Figure provides visualisation. It visualises a source-sink graph, where each unique path from source to sink represents one unique control point with above parametrisation. The cyan highlighted path represents the 91,2,3 Wi,1 21,2,2 where we multiply the values along the path from source sink. Notice last edges have value 1.
Wiy-1.h.77+ where Wo, := W,,,1- This assump-
In the Bézier buttress there are layers, one for each input dimension, and V4 + nodes each layer y l,...,dd d. Borrowing from neural network terminology, a forward -aaardd--sssssae is a sequential series of matrix multiplications which are element-wil wise warped with non- non-linearities, such as tanh or ReLU.
A buttress is architectural structure that provides support building.
5

---

we let our sequence of matrices be w1, W......, wa where wy is the matrix with entries {wi,k. hikk Let 1, denote the vector of size with in all entries, then fixing 'the input" to 1 I and the last matrix to 1 a forward pass is
Vd
1 I 1 w2 W2-
wal
I by P2
(18)
f1=0 12=0
id=0
We see this forward pass computes sum over all control points means. Naturally, we construct Bézier buttress for summing over variances too, which must be restricted to positive weights. was these sums that formed our bottleneck, in this parametrisation it is just sequence of matrix products.
What about computing f It comes down to a use of nonnlincaritties in the buttress. Multiplying element-wit wise the Bernstein polynomials as seen in Figure 3 (visualised only on 3rd layer), a forward pass computes either [z), or Var(f()) if using squared Bernstein polynomials. Each control point is then exactly multiplied by the correct polynomials from its way from source to sink. Notice, the 'input' is fixed to 1 in the source, but the observed appearing via the Bernstein polynomials along the way. We can write this too sequence of matrix products
11 12
Md Σ
E[f(x)]
Σ
B (21) -. -B (xa)
- ........... ....Bx. ..WaBr.uuat (19)
12=0
where Bz, diagonal matrix with the vy + Bernstein polynomials on its diagonal. For the variance of f(x there exists similar expression, of course with squared polynomials, and with the positive weights associated with the variance Bézier buttress.
On this inspection, all terms needed to compute Eq. (15) are available. All terms for the KL divergences are algebraic manipulations of Eq. (T8) these are explicit in the supplementary material. The key takeaway for Bézier buttress is that parametrise each random control point a product of weights. This can be seen as an amortisation, and we do inference on these weights rather than the control points themselves. Hence, no matrix inversions are needed to backpropagate through our objective, and a forward pass is only a sequence of matrix products.
2.3 Marginalising matrix commutativity
Matrix multiplication is not commutative. This implies the ordering of the matrices in Eq. 18 matter, which again implies how we order the input dimensions in the Bézier buttress is of importance. This the price for the computational benefit this parametrisation gives. An ordering of the input dimensions somewhat unnatural for spatial regression, so we present a way to overcome this in approximate Bayesian manner. Define f = fi + f2 + +.... + fr. where each of these individual fe. k (l,....r), are Bézier GPs with a random permutation of the ordering in the associated Bézier buttress. In other words, we let f be ensemble of Bézier GPs to (approximately) marginalise over all possible orderings. The number of all possible orderings quickly becomes too large for which to account; in practice, we set in feasible region, say 20, which gives satisfactory empirical performance.
Another lens on this that each control point a sum of r control points each control point's standard deviation is scaled by r-1, to obtain the same prior so far discussed. Oddly, we have then circumvented the problem of too many control points by introducing even more control points.
That is, each control point mean parametrise i i2....... = -h II Wi, (o where l denotes random permutation of (1.....d) d). We remind again that similar expression exist for control point variances, restricted to positive weights.
As remarked, inference comes down to a forward pass in the Bézier buttress, a sequence of d matrix multiplications. Assume all dimension are of order v, then the computational complexity of one forward pass is O (d(v 1)2). Now we need forward passes to marginalise the ordering, and forward passes, one for each observation, leaving final complexity of O(rdv²)). Linear in n and d.
3 Related Work
Variational inference in GPs was initially considered by Csató al. [1999] and Gibbs and MacKay 2000]. In recent times, the focus has shifted focus to scalable solutions to accommodate the big

---

data era. In this respect, Titsias 2009] took a variational approach to the inducing points methods [Quinonero-Came Candela and Rasmussen 2005]; later Hensman et al. [0033 further enhanced scalability to allow for mini- mini-batching and wider class of likelihoods. Still, the need for more inducing points is of importance, especially as the number of input features grows.
The response to this has mostly revolved around exploiting some structure of the kernel. Wilson
and Nickisch [2015] and Wu al. [202I] exploit specific structure that allow for fast linear algebra
methods; similar lar to our method inducing locations tend to lie on grid. These grids expand fast as the input dimension increases, as also pointed out earlier in the article. Kapoor et al. [2021] remedy this
by instead of rectangular grid, they consider the permutohedral lattice. such that each observation
only embeds to d + neighbours instead of 2", as in (iilson and Nickisch 2015).
Another approach to allowing for more inducing points is incorporating nearest neighbour search in the approximation [Kim al. 2005 "guyen--Tuong al.. 2008 Tran et al.] [2021] introduced saarse-withnn-sparss where they have many inducing points in memory, but each time search for the k-nearest ones to any observation. They discard the remaining ones as they have little to no influence
on the observation. Wu et al. [0022] made a variational pendant to this method.
Lastly, when dealing with nigh-diimensionll inputs it is worth mentioning feature learning. That is, learning more low low-dimensional features where GPs have better performance. The success of deep models has been transported to GPs by Damiaaou and Lawrence (2013) and later scaled to large datasets in [Salimbeni and Deisenroth 2017]. Another approach is Deep Kernel Learning [Wilson] et al. 2016 [Bradshaw al. 2017], where feature extraction happens inside the kernel function; lately Ober et al. [2021] has investigated the limitations and benefits of these models.
We have treated structured control points as our version of inducing points; and by parametrising them with a Bézier buttress, we limit the expansion of grids to linear growth in parameters. We are not the first to consider the Bernstein polynomials
basis for learning functions. Petrone [1999b] used it to model kernel estimate probability density functions, and several follow up works Petrone 1999a Petrone and Wasserman, 2002). Hug et al.) [2020) recently introduced Bézier GPs, but with
a focus on time series [Hug et al. [0022]. J. Our emphasis has been on spatial input; even the Bézier
surface literature contain close to nothing more than 2-dimensional surfaces.
Evaluation
We split our evaluation into four parts. First, we visually inspect the posterior toy dataset to show how the control points behave, and indicate that there indeed is a stationary- some standard UCI
dimensional
like behaviour
the domain of the hypercube. Next, Benchmark datasets, to gives insight into when Bézier GPs are applicable. After that, we switch to tall and wide data large both in the input dimension and in number of data points. These experiments give certainty that the method delivers on its key promise: scalability. Lastly, we turn our eyes to the method itself and investigate how performance influenced by the ordering of dimensions.
test empirically
Care is needed in optimising a Bézier GP not all parameters are born equal. We split optimisation into two phases. First, we optimise all variational parameters, keeping the likelihood variance o² fixed as T with T being the number of control points. After this initial phase, we optimise o² with all variational parameters fixed. We let both phases run for 10000 iterations with a mini-batch size of 500, for all datasets. Both phases use the Adam optimiser (Kingma and Ba 2015). the first phase with learning rate 0.001, and the second with learning rate 0.01. If not following such a bespoke training scheme, we see a tendency for the posterior to revert to the prior, because the KL-term becomes too dominating initially. This training scheme is designed for the Gaussian likelihood, but we wish to emphasise that, in principle, the loss function is accurate for any choice of likelihood.
4.1 One dimensional visual inspection
We hypothesised the objective function, Eq. 16, would ensure, within the hypercube domain, that f reverts to its prior in regions where data is scarce. To verify this we construct small dataset to inspect. We generate one-dimensional inputs uniformly in the regions [0, 0.33) and [0.66, 1]; we sample 20 observation in each region. The responsive variable is generated as v(²) = sin (16r).). According to the hypothesis, I should in the region [0.33, 0.66), tend towards zero in mean, and increase its variation here. We use BézierGP of order 20 to model the observations, since they are highly
7

---

Oiavevation
Prior control point
Proteeiar control point
off
10
0.10
0.05
0.00
0.2
0.5
Figure 4: Left: Posterior distribution of f. Middle: Posterior and prior distribution of control points.
Right: Posterior variance function over the domain. Variance increases in scarce data regions.
power
boston
bike
Regadrected
9566
45770
758 8
100
17379
48537
1030 8
165699 18
10
37
Tewilog-Baellbove
100
6:01
-2 939307
-16450
0.07
-2
0.05
06:0R 0.08
-3.1819
0.9167
-2.7400
228479
2.8199-0000
0.13
-2.4002
0.04
07715 4 M0S
-0000 0.08
0,5077 M/M
SimplexCiP
3,4416
3.2745
NA
1,0502
0.19
410008
0.02
-2 9015
0.01
-2.9855
0.01
-0607
0:25
16 179
186
0.472
0.14
0.580 + 0.06
-3. -3.0112
0.62
0-9634
BezierTPP
10
2.7796
0.04
2.9267
0.00
0.7163
4 0.46
197
140.5
0,7000
0.15
D.6799 4 0.06
4,5613 à: 1.20
0,5965 d: 0,62
30
-2. 7301
0.03
-2,8002
0.01
-00501 0.02
1550.2
52.7
0.7475
0.32
01009 + 0.06
-7.7176 + 2.70
0.5055 :: 0,06
Test RMSI
SGPR
3.7U81
3
9115
4 =
0.00
1
011
012
03389 0.02
00303= 0.02
0.1279 #
01212
5,3512
SimptexCip
3.1147
0.13
0:257 0.07
5.4457 4.75
0,2256
3,9050 0.15
LG020 a 0.04
0.4348 1 D0D5
47007 = 0.0L
0.1474
0.02
D.1019 2 D.B
4,8127 2 0.56
0,2030 :: 0,55
BesterCIP
10
3 3.9675 0:15
4.5122
a 0.04
0.4157
2 0.13
5,0712 1.76
0.1300 0.01
0.1735 0.22
4,7917 2: 1154
1028.5 4c3
20
3,7427
0.18
4.5533
= 0.04
0.5575
8 0.11
2.00
01821
0.01
90.05 + 318.4
4.822 8 0.55
Tell + Se12
Table I: Results on eight standard UCI Benchmark datasets. We list test-set log-likehood and RMSE,
both are averages over 20 train/test splits. On top are test log-likelihoods (higher is better), and bottom is test RMSE (lower is better). NA indicates Cholesky error.
non-linear. Figure 4 shows the posterior distribution of f to the left: we observe f tends towards the prior in the middle region. The middle plot illustrates the distribution, both prior and posterior, of the 21 control points. There is a clear tendency for the central- points to align the posterior and posterior, enforcing this behaviour in f. The non-equal priors are due to the inverse- inverse-squared Bernstein adjusted prior which ensures uniform variation in f over the domain, see Figure 2 right in Figure 4 shows the behaviour foundational to practitioners of Bayesian optimisation and active learning etc., the variance increase away from data regions.
The plot to the
4.2 UCI Benchmark
We evaluate on eight small to mid-size real world datasets commonly used to benchmark regression [Herandezz-Loatto and Adams 2015). We split each dataset into train/test-spli with the ratio 90/ 10. We do this over 20 random splits and report test set RMSE and log log-likelihood average and standard deviation over splits. We choose baselines to be SGPR, following the method from Titsias [2009]; we do both for 100 and 500 inducing variables. SimplexGP is another baseline, they suggest their approximation is beneficial for semi- semi-high dimensional inputs (between and 20) [Kapoor al., 2021], hence they are an obvious baseline. SimplexGP usually use validation set to choose the final model. This due to a highly noisy optimisation scheme using a high error-tolerance (1.0) for conjugate gradients. We remedy this by setting the error-tolerance to (0.01), which harms scalability, but we can omit using a validation set for better comparability. Wang et al. (0119] recommend this error-tolerance, but remark it is more stable for RMSE than for log-likelihood. For BézierGP, we fix the number of permutations to r = 20, and vary the order in v = 5, 10, 20. The order is identical over input dimensions. The inputs are pre- pre-processed such that the training set is contained in [0, 1]dd
Table contains the results of this experiment. We make the following observations about our presented BézierGP. On keggdirected and elevators there are test points outside the defined domain on some splits, which cause the RMSE be extreme, but the likelihood more forgiving. This highlights the constraint of our model: it needs a box- bou-deedde domain to be priori known. Had we standardised such that both test and train data were in the hypercube, BézierGP (v 20) would have a average test RMSE of 0.0937. We could not reproduce results from Kapoor et al. 2021 on keggdirected, the optimisation was too noisy and with no use of validation. On concrete, boston and energy overfitting tendencies. Even though BézierGP is the optimal choice on energy
8

---

RUSE
RUSSE
400
5.40
SVISP/00
0.25
-330
0.08
3..
20.0
0.20
0.50
35-
#0012
-.0
#000P123
400RR512
0.25
3.4
60000
0.10
50092112
0.26
4.5
4500
DocoPP
+830
0.05
Figure 5: Test negative log-likel (lower is better) and RMSE (lower is better) for large scale regression. Colours indicate the origin of numbers: red from Salimbeni and Deisenroth (2017). purple from Wang et al. 2019], and cyan from Kapoor et al. 2021 Orange is ours. We observe our BézierGP is highly competitive on large datasets.
there is a mismatch between train and test error. On concrete this shows in better test RMSE, than the baselines, but the variance is overfitted yielding non-optimal likelihood. We conjecture this happens because the n/ /drratii is low; which makes it more likely to overfit the control points especially for higher orders v. Knowing these model fallacies, we observe that BézierGP outperforms on baselines on multiple datasets, most notably the 17-dimensional bike dataset.
4.3 Large scale regression
Figure 5 shows the results of regression tasks in regimes of high dimensions and one in high number of observations. Here, we follow exactly the experimental setup of either Salimbeni and Deisenroth [2017) or Waag et al. [2019]. If the latter, we use the validation-spli they use as training data. Our optimisation scheme for BezierGP is consistent with above, except for slice, where the first training phase runs for 30000 iterations. We discard test points that not in the [0, 1]d domain in no situation did this remove more than 0.001% % of the test set. The number after DGP, denotes the number of hidden layers in Deep GP (Salimbeni and Deisenroth (017]. after SGPR and SVGP it denotes the number of inducing points. SVGP refers to the method from Hensman et al. 2013 After B, it denotes the order used in BézierGP.
On year, we observe our (non-deep) model is on-par with 2-layered Deep GPs, and closer to 3 in RMSE. The highest dimensional dataset, slice, sees in the low nıd-ratio again, and
power (24)
protein (302880)
bike (3e14)
RMSE
LL
RMSE
L.L
RMSE
LL
r=1
4.1476
-2.8418
4.8418
-2.9962
0.1477
0.4813
we are again faced with a too flexible model.
3.8028
-2.7549
-2.9080
0.7275
This is why we report results for orders 3 and 5, rather than 20, since the overfitting kicks
r- 20
4.2514
2.8826
4:3713
2.8043
0.0853
0.5454
30
3.8257
-2.7551
4.3377
-2.8677
0.0775
0.8372
3.6612
2.7100
2,8888
0.5024
in. Even for these small orders the test log- likelihood has high variance and under-performs compared RMSE. With respect RMSE it top-performer signalling again it is overfitting the variance. On the remaining two datasets BézierGP is best-performing among baselines.
r= 50
3.7841
-2.7503
4.2614
-286666
0.0718
0.9009
Table 2: Performance in test RMSE and log- likelihood against the number of permutations, on three datasets. In parenthesis shows the num- ber of possible permutations of input dimensions. Higher dimensional datasets show improving per- formance with increasing r.
4.4 Influence of number of permutations
All experiments so far used r = 20; that is, 20
random permutations of the ordering of dimension used in the Bézier buttress. For a problem with input dimension d, there exist d! possible permutations. Table shows results with varying for each dataset, the results are over the same train/test split (0.9/ 0.1)). We fixed // = 20. Up to some for the two highest dimensional datasets, protein and bike,
noise in the optimisation phase, performance improves with higher r. Bike has over 50% reduction in RMSE from r to r = 50.
Table 2 emphasises the results we have presented are not optimised over hyperparameter r and v. They also illustrate an interesting direction of future research: optimising these hyperparameters. We chose permutations by random sampling, but choosing them in principled deliberate manner could yield good performance with a computationally manageable r. This result indicates, at least, protein and bike would see increased performance in Table I from better (or just more) permutations.

---

Discussion
We introduced the Bézier Gaussian Process GP, with polynomial kernel the Bernstein basis, that scales to large number of observations, and remains space-filling for high number of input features (limited to a box-bounded domain). We illustrated that, with slight adjustments, the prior and posterior have similar behaviour 'usual' stationary kernels. We presented the Bézier buttress, weighted graph to which we amortise the inference, rather than inferring the control points themselves. The Bézier buttress allows GP inference without any matrix inversion. Using the Bézier buttress, we inferred 6385 control points for the high- nigh-dimensional slice dataset.
We highlighted weaknesses of the proposed model: most crucially the tendency of overfitting when the n/ n/d-ratio is low. The results demonstrate scalability in both n and d. but does not solve the short, but wide problem. The paper did not optimise over the hyperparameters of the proposed kernel, namely v and r, but it showeased briefly that doing so might enhance BézierGPs empirically; especially smart selection of the permutations is an interesting direction for future research. We speculate that optimising over orders, on a validation set, would alleviate some of the overfitting issues.
Acknowledgments and Disclosure of Funding
MJ is supported by the Carlsberg Foundation.
References
J. Bradshaw, A. G. d. G. Matthews, and Z. Ghahramani. Adversarial examples, uncertainty, and trans- fer testing robustness in Gaussian process hybrid deep networks. arXiv preprint (rXXiv:707.0044/: 2017.
D. R. Burt, C. E. Rasmussen, and M. van der Wilk. Convergence of sparse variational inference in Gaussian processes regression. Journal of Machine Learning Research (JMLR), 21(131): 2020.
L. Csató, E. Fokoué, M. Opper, B. Schouty, and O. Winther. Efficient approaches gaussian process classification. Advances in neural information processing systems, 12, 1999.
A. Damianou and N. D. Lawrence. Deep Gaussian processes. In Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS), 2013.
J. Gardner, G. Pleiss, R. Wu, K. Weinberger, and A. Wilson. Product kernel interpolation for scalable gaussian processes. In A. Storkey and F. Perez-Cruz, editors, Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 1407-1416. PMLR, 09-11 Apr 2018. URL https:// proceedings mlr
hat
M. N. Gibbs and D. J. MacKay. Variational gaussian process classifiers. IEEE Transactions Neural Networks, I166): 1458-1464, 2000.
J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In Proceedings of the 29th Conference on Uncertainty in Artifical Intelligence (UAI), 2013.
Hernandez-Looatto and R. Adams. Probabilistic backpropagation for scalable learning of bayesian neural networks. In F. Bach and D. Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 1861-1869 Lille, France, 07-09 Jul 2015. PMLR. URL https:: mww..ng..........mmmmm..mm.mmm.m miw.press/ v37 +37/hernandez: hww.aaaanndee-oootoo.......... www.loco......tcoooooooo
T. Hildebrandt and I. Schoenberg, On linear functional operations and the moment problem for finite interval in one or several dimensions. Annals of Mathematics, pages 317-328, 1933.
M. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. Journal of Machine Learning Research, 2013.
10

---

R. Hug, W. Hübner, and M. Arens. Introducing probabilistic bézier curves for n-step sequence prediction. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 10162-10169, 2020.
R. Hug, S. Becker, W. Hübner, M. Arens, and J. Beyerer. Bézier curve gaussian processes. arXiv preprint arXiv:2205.0/77 2022.
Kapoor, M. Finzi, K. A. Wang, and A. G. G. Wilson. Skiing on simplices: Kernel interpolation on the permutohedral lattice for scalable gaussian processes. In International Conference on Machine Learning, pages 5279-5289. PMLR, 2021.
H.-M. Kim, B. K. Mallick, and C. C. Holmes. Analyzing nonstationary spatial data using piecewise gaussian processes. Journal of the American Statistical Association, 100(470);653-666 668, 2005.
D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In 3rd International Conference on Learning Representations (ICLR), 2015.
D. Nguyen-Tuooe J. Peters, and M. Seeger. Local gaussian process regression for real time online model learning. Advances in neural information processing systems, 21, 2008.
S. W. Ober, C. E. Rasmussen, and M. der Wilk. The promises and pitfalls of deep kemnel learning. In Proceedings of the 37th Conference on Uncertainty in Artifical Intelligence (UAI), 2021.
S. Petrone. Bayesian density estimation using bernstein polynomials. Canadian Journal of Statistics, 27(1):105-126, 1999a.
S. Petrone. Random bernstein polynomials. Scandinavian Journal of Statistics, 26(3):373-393, 1999b.
S. Petrone and L. Wasserman. Consistency of bernstein polynomial posteriors. Journal of the Royal Statistical Society: Series (Statistical Methodology). 64(1):79-100 1000 2002.
H. Prautzsch, W. Boehm, and M. Paluszny. Bézier and B-spline techniques, volume 6. Springer, 2002.
J. Quinonero-Candele and C. E. Rasmussen. A unifying view of sparse approximate gaussian process regression. The Journal of Machine Learning Research, 6:1939-1959, 2005.
H. Salimbeni and M. Deisenroth. Doubly stochastic variational inference for deep gaussian processes. Advances in neural information processing systems, 30, 2017.
E. Snelson and Z. Ghahramani. Sparse gaussian processes using pseudo-inputs. Advances in neural information processing systems, 18, 2005.
M. Titsias. Variational learning of inducing variables in sparse gaussian processes. In Proceedings of the 12th International Conference on Artificial Intelligence and Statistics (AISTATS), 2009.
G.-L. Tran, D. Milios, P. Michiardi, and M. Filippone. Sparse within sparse gaussian processes using neighbor information. In International Conference on Machine Learning, pages 10369-10378. PMLR, 2021.
K. Wang, G. Pleiss, J. Gardner, S. Tyree, K. Q. Weinberger, and A. G. Wilson. Exact gaussian processes on million data points. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances Neural Information Processing Systems, volume 32, Curran Associates, Inc., 2019, URL http:/::///popeeepeeeeeee htwpeeeeiiinnnneniiiiiieeeeeeeee neurips cc pc/paperr009///ile 0Le86988065500d5516666633333 Paper pdf
C. K. Williams and C. E. Rasmussen. Gaussian processes for machine learning. MIT Press Cambridge, MA, 2006.
A. Wilson and H. Nickisch. Kernel interpolation for scalable structured gaussian processes (kiss-gp).
In International conference on machine learning, pages 1775-1784. PMLR, 2015.
A. G. Wilson, Z. Hu, R. Salakhutdinov, and P. Xing. Deep kernel learning. In Proceedings of the 19th International Conference on Artificial Intelligence and Statistics (AISTATS), 2016.
11

---

L. Wu, A. Miller. L. Anderson, G. Pleiss, D. Blei, and J. Cunningham. Hierarchical inducing point gaussian process for inter-domian observations. In International Conference on Artificial Intelligence and Statistics, pages 2926-2934. PMLR, 2021.
L. Wu, G. Pleiss, and J. Cunningham. Variational nearest neighbor gaussian processes. arXiv preprint arXiv:2202.01694, 2022.
12

---

Checklist
The checklist follows the references. Please read the checklist guidelines carefully for information on how to answer these questions. For each question, change the default [TODO] to [Yes], [No], , or [N/A]. You are strongly encouraged include justification your answer, either by referencing the appropriate section of your paper or providing a brief inline description. For example:
Did you include the license to the code and datasets? [Yes] See Section ??.
Did you include the license the code and datasets? [No] The code and the data are proprictary.
Did you include the license to the code and datasets? [N/A]
Please do not modify the questions and only use the provided macros for your answers. Note that the Checklist section does not count towards the page limit. In your paper, please delete this instructions block and only keep the Checklist section heading above along with the questions/answers below.
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [No] The work does not have any foresecable negative societal impact.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL)? [Yes]
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes]
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes]
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [No] Our method ran on a laptop for all but experiment. GPUs used for some baselines.
4. If you using cxisting assets (e.g., code, data, models) or curating/releusing assets...
(a) If your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license of the assets? [No]
(c) Did you include any new assets either in the supplemental material or URL? [No] Code is in supplements.
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] Only open source data used, personal information.
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "pOjDvRqbqz", "reviewer": "Reviewer_qTRQ", "summary": "The paper proposes a new Bézier buttress-based kernel that allows for scaling GP regression both in number of inputs (tall data) and in number of input features (wide data). The scalability is shown both theoretically and empirically through evaluation on popular datasets.\n", "questions": "I have a several questions to the authors.\n\n1. Am I correct in understanding that computing full covariance would be just rendered as another type of \"non-linearity\" in the Bézier buttress, that is, multiplying be a cross-product of Bernstein polynomials instead of the squares of them?\n2. Have you tried to make even larger orders $\\nu$ work? Does it make sense to try larger orders of the Bézier?\n3. Have you tried to faithfully compute the true posterior (on a small example perhaps) and see how good the proposed variational approximation is?\n4. How easy is it to introduce non-Gaussian likelihood?\n5. In your opinion, does it make sense to vary the order of each dimension (that is, non-equal $\\nu_\\gamma$)? It would make sense as not all dimensions are equal. That would just render a Bézier buttress with different layer widths.\n6. A major thing is an assumption of the line 142 which makes Bézier buttress possible. This introduces inter-relations between variational parameters. Do you think this is not very restrictive of an assumption? Again, comparing to the true posterior would probably be good enough of an answer.\n7. What happens when jointly optimizing likelihood noise and variational parameters? Have you observed a worsening of the performance?\n", "limitations": "The authors address the limitations and constraint of their proposed method. The distribution of the data must be supported on an a-priori known compact cube. Another constraint is a tendency to overfit when $n/d$ ratio is small. Orders $\\nu$ larger than 25 do not work. \n\n", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "*Originality*. The work is original. While there are other papers leveraging Bézier curves for Gaussian processes, the paper provides a way to parametrize Bézier control points to allow for a scalable training and inference (termed Bézier buttress).\n\n*Quality*. The research is sound. The statements the authors make are supported by evidence. The authors make a detailed analysis of empirical results, addressing both advantages and limitations of the suggested method. They provide a review of the related work, and highlight areas of possible further research.\n\nThat said, a more detailed overview of performance against order of the Bézier would be nice. \nAnother thing beneficial for the paper would be reports on the wall time consumed and computing power used in the experiments. \n\n*Clarity*. The work is written in a clear language, the presentation is on-point and laid out in a way that facilitates reading. The notation might be cumbersome sometimes.\n\nOne thing to pay attention to is in the supplementary material, line 5: the definition of the variance buttress weights. I think there is a forgotten dependence of $v$ on $\\gamma$ and also $\\zeta_{\\gamma_i}$ doesn't make sense, I guess it should be $\\zeta_\\gamma(i_\\gamma)$.\nThere is also a typo on line 142: \"can be parametrized\".\n\n*Significance*. Since applying Gaussian processes to high-dimensional inputs is known to be hard, the proposed method that deals with both \"wide\" (high-dimensional inputs) and \"tall\" (many inputs) data is a significant contribution to the world of Gaussian process machine learning. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "dSkx2heLThq", "reviewer": "Reviewer_mrad", "summary": "This work introduces a type of kernel for Gaussian process regression that requires only a linear computational cost in both number of observations and input features. The scaling is achieved by means of the introduction of the Bézier buttress structure, which allows approximate inference without computing matrix inverses or determinants.", "questions": "---Specific comments---\n\nlines 13 to 16:\nIt might be useful to include references for Gaussian process model and the approximate methods for their scalability.\n\nline 33:\nInclude references for Bayesian optimisation, active learning and reinforcement learning.\n\nline 63:\nAs per Titsias' work (2009), I would suggest to use inducing variables when referring to $\\boldsymbol{u}$, and inducing inputs/locations or points when referring to $\\mathbf{Z}$.\n\nline 65:\nbold for $\\boldsymbol{\\Sigma_u}$?\n\nline 72:\nWhy was the lower case variable $\\mathbf{p}_{i_1,...,i_d}$  from Eq. (4)\n\nchanged to upper case variable $\\mathbf{P}_{i_1,...,i_d}$?\n\nis it a matrix now? What is the dimensionality of $\\mathbf{P}_{i_1,...,i_d}$? Is it the same of Equation (1) where it was defined $\\mathbf{p}_i \\in \\mathbb{R}^D$?\nIs it just a scalar? if so, then most of the bold notations are problematic.\n\nline 73:\nIt is necessary to introduce the mean $\\boldsymbol{\\vartheta}_{i_1,...,i_d}$ \n\nand covariance $\\boldsymbol{\\Sigma}_{i_1,...,i_d}$ of the Gaussian variables. \n\nSince they are defined in bold it means they are not scalars, doesn't it?\n\nline 81:\nTypo? shouldn't it be $\\mathbf{0}$ since the mean is a vector?\nAlso why bold \"$i$\" in \"for all $\\mathbf{i}$\"? What is $\\mathbf{i}$?\n\nline 82:\nwhy is the covariance matrix $\\boldsymbol{\\Sigma}_{i_1,...,i_d}$ equal to a scalar? \n\nI mean in $\\boldsymbol{\\Sigma}_{i_1,...,i_d}=1$.\n\nline 90: What is $1_{\\nu_{\\gamma}+1}$? The term was not introduced. \n\nShould it be $(A_{\\gamma})_{i,j}$? It is not clear how both equations in (10) are related. \n\nline 91:\nHow is $\\zeta_{\\gamma}(i_{\\gamma})$ related to equation (10)? It is not clear.\n\nline 95:\nThe statement: \"This adjustment works up to $\\nu_{\\gamma}= 25$, after which negative values occur.\" is an empirical fact? or is there any mathematical explanation for such a fact? Does it happen regardless of the dimensionality?\n\nline 104 (Figure 2):\nTypo, get rid of \"is\" in \"the variance is of $f$ is more...\"\n\n---The comments below will appear separated a bit strange, but it is due to a problem in the reviewing platform---\n\nline 108 :\nIt is a bit confusing to write, $\\mathbf{P}_{i_1,...,i_d}$ following a Gaussian with\n\nthe mean $\\boldsymbol{\\hat{\\vartheta}}_{i_1,...,i_d}$ \n\nand covariance $\\boldsymbol{\\hat{\\Sigma}}_{i_1,...,i_d}$\n\nsince it was already defined for the prior, $\\mathbf{P}_{i_1,...,i_d}$ (after Eq. (7)) to follow a Gaussian with\n\nthe mean $\\boldsymbol{\\vartheta}_{i_1,...,i_d}$ \n\nand covariance $\\boldsymbol{\\Sigma}_{i_1,...,i_d}$.\n\nActually for the \"variational control points\" it should be something like $\\mathbf{P}_{i_1,...,i_d}|\\mathbf{y}$ following a Gaussian with\n\nthe mean $\\boldsymbol{\\hat{\\vartheta}}_{i_1,...,i_d}$ \n\nand covariance $\\boldsymbol{\\hat{\\Sigma}}_{i_1,...,i_d}$,\n\ngiven that it is the (approximate) posterior. Nonetheless, I suggest to define in terms of each variational distribution $q(\\mathbf{P}_{i_1,...,i_d})$ equal to  a Gaussian with\n\nthe mean $\\boldsymbol{\\hat{\\vartheta}}_{i_1,...,i_d}$ \n\nand covariance $\\boldsymbol{\\hat{\\Sigma}}_{i_1,...,i_d}$,\n\nthis in order to ease the understanding of the following equations (15) and (16).  \n\nSpecially to differentiate the prior $p(\\mathbf{P}_{i_1,...,i_d})$ \n\nfrom the approximate posterior $q(\\mathbf{P}_{i_1,...,i_d})$.\n\nline 111:\nHow is $q(\\mathbf{P})$ related to each $q(\\mathbf{P}_{i_1,...,i_d})$? It need to be defined.\n\nlines 110, 113 and 120:\nthe notation for the variable $y$ should be bold, I mean a vector or set that stacks all the $y_j$, say $\\mathbf{y}$={$y_j$}$^n_{j=0}$. \n\nline 116:\nWhy is it mentioned a variable $\\sigma^2$ if the likelihood has not even been introduced?\n\nline 120 Eq. (15):\nThere should be big brackets containing all the terms of the summatory.\n\nline 127:\nIt might be better to make Figures appear after they are mentioned in the text.  \n\nline 151:\nplease introduce the variables $1_{\\nu_{1}+1}^{\\top}$ and $1_{\\nu_{d}+1}^{\\top}$. \n\nline 225:\ncorrect wording in \"..., to gives insight into when...\"\n\nline 243:\nTypo? include \"be\" in \"..., $f$ should in the region...\"\n\nline 256:\nintroduce acronym SGPR.\n\nlines 267 to 269:\nIt might be useful to know what are the results for keggdirected and elevators, but guaranteeing that the test points are not outside of the domain for the splits.\n\nline 274:\nThe authors explain that for the concrete dataset the variance overfits yielding a non-optimal likelihood, and conjecture this might be due to $n/d$-radio being low. Nonetheless, the $n/d$-radio for energy is not far from concrete's, and we see from the Table 1 that the energy results were much better than SGPR. How did the authors come up with such a conjecture?  \n\nline 325: please upper case names in the references like: Gaussian, Bayesian, Bernstein,\n\n---Other Questions---\n\nI do not remember having seeing it in the paper, how are the predictions of the model performed? How does the computational complexity of the predictions differ from the classical SGPR?\n\nAn interesting property of a Sparse GP is that when the number of inducing points approaches to the number of data observations, we would expect to obtain almost the same results of a full (non-sparse) GP. Is this fact tested in the Bézier GP approach, I mean that when the number of control points approaches to the number of observations we would almost end up with the original Non-sparse GP?\n\nAlso, if we set the inducing points as exactly the same to the original input observations, we would expect to recover the original full (non-sparse) GP, is it possible to think about the Bézier GP approach as able to recover the full GP where the input data observations are the same control points?\n\nWould it be feasible to think about a Deep Bézier Gaussian Process? What could it be the challenges and limitations?\n\nThe construction of the GP allow the use of different orders $\\nu_d$ per dimension, when might that case of having different orders per dimension be useful or applicable?", "limitations": "The authors were critic about the limitations of their work.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Originality: The work seems original and broadly different to the classical sparse Gaussian process approaches.\n\nQuality: The experiments provided in the paper can be useful for practitioners in the field of Gaussian process models and can open an alternative path for researching about the scalability of Gaussian process models. \n\nClarity: The work is generally well written and not difficult to follow, though it presents some problems in the mathematical notation.\n\nSignificance (importance): Nowadays the application of Gaussian process models in the context of large datasets with a wide number of features in the input space is still limited due to scalability problem. Therefore, alternative GP approaches that aim to enhance the applicability of GP in such contexts of large and wide data are of great value for the research community.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "bQCNTI2b9oR", "reviewer": "Reviewer_W66W", "summary": "The paper presents a Gaussian process method focusing on data in which the number of samples and the number of dimensions are both high. The GP prior is built upon Bezier surfaces and is formulated as a combination of Gaussian random variables and coefficients taken from Bernstein polynomials. The paper suggested adjusted priors and scalable construction via Bezier buttress. Experiments are performed on various data sets. ", "questions": "- How to perform the matrix inverse of $k(Z,Z)$? Is it done analytically by using inverse squared Bernstein.  If so, I would suggest to put more explanation or reference for inverse squared Bernstein and Bernstein polynomials which I feel is missing in this paper. (I’m not familiar with the literature here).\n \n- The construction of $f(x)$ somewhat resembles deep structures (of course with no activations), the comparison can be done extensively with other models like deep Gaussian process or deep kernel learning. However, the paper reports just one result with DGP for the “year” data set. Are the results for the remaining data sets (buzz, house electric, slice) available?\n\nMinor points:\n\t-The optimization is somehow pragmatic. Since the paper suggests KL terms may be big at initial steps, I think KL annealing techniques can be used here. \n", "limitations": "Please see above.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The problem of handling high-dim GP is interesting , the approach is novel to tackle the problem in GP research.\n- The paper gives a good analysis for the method and provides a scalable solution for GP.\n\nWeakness:\n- I have one concern of how well Bezier d-surfaces can approximate data. What is the error as the number of dimensions increases? If there is any reference for this, it would be nice to add the clarification.\n- It might be out of scope in this paper, there is no theoretical justification for that the proposed approach can address effectively high-dim data .This may be related to the previous point and the way the paper suggests constructing Bezier buttress. The empirical results (first plot of Fig 5.) may. The main concern is that we do not know whether the ability to handle high-dim data comes from deep networks or from the proposed approach.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}, {"review_id": "BUSfYTrUfp", "reviewer": "Reviewer_rsfL", "summary": "The manuscript proposes a methodology to extend the Gaussian process framework towards datasets that are both high-dimensional (wide) and have a large sample size (tall) using multivariate Bézier curves. In particular, the authors propose to perform variational inference over a multivariate linear model over Bernstein polynomials having exponentially many control points as parameters. To render the model scalable, the authors propose two approximations: They use a factorisation assumption to reduce the number of control point parameters and they neglect the order of the summation over the dimensions. Experiments on several UCI benchmarks conclude the paper.", "questions": "* Computations with Bernstein polynomials usually require special care regarding numerical stability i.e. De Casteljau's algorithm. Do you encounter issues in that respect?\n* The \"adjustment\" of the Bernstein polynomials might be numerically challenging as shown by Figure 2 right Panel. Is this correct? What is the condition of the linear system in Equation (10)?\n* How would experimental results change without that correction?\n* I'm missing an analysis and more detailed evaluation on the effects of the two approximations made to render the model scalable: the factorisation structure and the order of summation.\n* In particular, how is the marginal likelihood bound affected? Is it still a bound or a mere approximation that cannot be trusted as a comparative measure in the experiments? Columns 2,4,6 of Table 2 suggests exactly that behavior. Then, one could only compare RMSE and not probabilistic model properties. More specifically, the comparisons in Figure 5, in the left panels and Table upper part would be questionnable in their meaning.\n* What would need to change if the likelihood was non-Gaussian? I'm assuming that the same construction can be used.", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "1) Clarity\nThe paper is mostly well written and the technical content is rather well accessible.\n\n2) Originality\nThe proposed polynomial covariance structure using Bernstein polynomials seems to not have been explored before.\n\n3) Significance\nThe method is (currently) specific to large scale regression and could be useful for GP regression. However, in the current state, the underlying approximation assumption need to be better analysed before the model can be used as a black box.\n\n4) Reproducibility\nUnfortunately, only some vague outline of the code is contained in the appendix suggesting that there might be an implementation that might be shared at some point. The relevant classes i.e. BezierButtress are not contained hence this cannot be judged. Assuming that the code will be shared, the experiments should be well reproducible. If not, this should be rather difficult.\n\n5) Empirical analysis\nThe manuscript contains experiments on several standard regression datasets regarding marginal likelihood and prediction error.\na large number of numerical results gathered in a big table regarding marginal likelihood, prediction error and runtime. Some more results on more subtle aspects of the covariance function itsel would have been much appreciated.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept