High-dimensional Additive Gaussian Processes under Monotonicity Constraints
Authors: Andrés F López-Lopera
Francois Bachoc
Olivier Roustant
Venue: NeurIPS 2022 Accept (2022)

Abstract: We introduce an additive Gaussian process (GP) framework accounting for monotonicity constraints and scalable to high dimensions. Our contributions are threefold. First, we show that our framework enables to satisfy the constraints everywhere in the input space. We also show that more general componentwise linear inequality constraints can be handled similarly, such as componentwise convexity. Second, we propose the additive MaxMod algorithm for sequential dimension reduction. By sequentially maximizing a squared-norm criterion, MaxMod identifies the active input dimensions and refines the most important ones. This criterion can be computed explicitly at a linear cost. Finally, we provide open-source codes for our full framework. We demonstrate the performance and scalability of the methodology in several synthetic examples with hundreds of dimensions under monotonicity constraints as well as on a real-world flood application.

High-dimensional Additive Gaussian Processes under Monotonicity Constraints
Andrés F. López-Lopera CERAMATHS, UPHF 59313 Valenciennes, France andres 1opez11perra@phh fr
François Bachoc IMT, UMR5219 CNRS Université Paul Sabatier 31062 Toulouse, France
Olivier Roustant IMT, UMR5219 CNRS INSA Toulouse 31077 Toulouse, France
Abstract
We introduce an additive Gaussian process framework accounting for monotonicity constraints and scalable to high dimensions. Our contributions are threefold. First, we show that our framework enables to satisfy the constraints everywhere the input space. We also show that more general componentwise linear inequality constraints can be handled similarly, such as componentwise convexity. Second, we propose the additive MaxMod algorithm for sequential dimension reduction. By sequentially maximizing squared-norn criterion, MaxMod identifies the active input dimensions and refines the most important ones. This criterion can be computed explicitly linear cost. Finally, we provide open-source codes for our full framework. We demonstrate the performance and scalability of the methodology several synthetic examples with hundreds of dimensions under monotonicity constraints as well as real-world flood application.
Introduction
The framework of additive and constrained Gaussian processes. In high dimension, many statistical regression models are based on additive structures of the form:
y( (..) 3/1 (211) +... ya(xa))
Although such structures may lead to more "rigid" models than non-additive ones, they result in simple frameworks that easily scale in high dimensions 2. Generalized additive models (GAMs) [I] and additive Gaussian processes (GPs) [ are the most common models in a wide range of applications. The latter can also be seen generalization of GAMs that allow uncertainty quantification. As shown in [3 additive GPs can significantly improve modeling efficiency and have major advantages for interpretability. Furthermore, in non-additive dimenioammniion GP models, adding inequality constraints (in particular monotonicity which will be the main focus of this paper) leads to more realistic uncertainty quantification learning from real data 10
Contributions. Our contributions are threefold. 1) We combine the additive and constrained frameworks propose additive constrained GP (cGP) prior. Our framework is based on fiite-dimensional representation involving one-dimensional knots for each active variable. The corresponding mode predictor can be computed and posterior realizations a scalable way to high dimension. 2) We suggest the additive MaxMod algorithm for sequential dimension reduction. At each step, MaxMod either adds new variable the model inserts knot for an active one. This choice is made by maximizing the squared-norm modification of the mode function, for which we supply exact expressions with linear complexity. 3) We provide open-source codes for our full framework. We demonstrate the performance and scalability of our methodology real-world application in dimension
be sampled, both in
with numerical examples in dimension up to 1000 as well 37. MaxMod identifies the most important input variables, with data size as low as n = 2d in dimension d. It also yields accurate models satisfying the constraints everywhere on the input space.
36th Conference on Neural Information Processing Systems (NeurIPS 2022).

---

Range of applicability. The computational bottleneck of cGPs is sampling from the posterior distribution. Here, it boils down to sampling a constrained Gaussian vector with dimension equal to the number of knots. This is done with Hamiltonian Monte Carlo (HMC) [ which currently works with several hundreds of knots. Notice that MaxMod enables to minimize this number of knots. Our framework is illustrated with monotonicity constraints and can be directly applied to other componentwise constraints such as componentwise convexity. These constraints should be linear and such that satisfying them on the knots is equivalent to satisfying them everywhere (see Section 3.2
Related literature. Additive GPs have been considered in 3 12 14 to name a few. The benefit of considering inequality constraints in (non-additive) GPs is demonstrated in [5 [5] and in application fields such as nuclear safety [3), geostatistics [6. tree distributions 7 econometrics [16]
coastal flooding [17), and nuclear physics [18]. Note that there is bridge between cGPs and
(frequentist) splines-base regression with reproducing kernel Hilbert spaces (RKHS) 19 20
Our cGP model and MaxMod algorithm are extensions of the works in 5 and [1 (respectively)
to the additive case. Besides GPs, [22] provides a stochastic optimization framework with RKHS that yields regression with inequality constraints.
[23] suggests non-parametric model for non-negative functions; a problem related to that of monotonicity constraints. To the best of our knowledge, our framework is the first that enables to satisfy the constraints (especially monotonicity) everywhere and to scale to high dimension (up to one thousand our experiments). In particular, the aforementioned cGP works are not applicable in these dimensions.
Furthermore, [2 does not satisfy the constraints everywhere, but only approximately and in expec- tation, and 23 does not directly address monotonicity constraints. For instance, with increasing constraints in dimension one, our framework provides posterior sample functions h satisfying h(u) < h(v) for all S v. In contrast, the derivative of h, (w), would be guaranteed to be non-negative only at some points u [15]. Also, h((w) would only be guaranteed to be above some strictly negative threshold, in expectation in [22]
Note that for unconstrained GP regression, the likelihood and conditional distribution are explicit but their computation cost scales as O(n)), where is the size of the dataset. An active strand of literature then focuses on providing efficient numerical procedures [24, 25]]. Because of the constraints, our work is complementary to this strand of literature. Indeed, with cGPs, even if n small. the numerical challenge of sampling from the posterior distribution remains, discussed previously. The main indicator of complexity is thus the number of knots, rather than the dataset size.
The references 24 25
provide convergence guarantees. At this stage, we do not provide convergence guarantees in our setting. To reach that goal, a possibility would be to extend the convergence proof of MaxMod given in 21 from the non-additive to the additive case. Two comments are in order on this. First, the proof techniques in [2T] are different from these in [24] 25 (and also from standard asymptotic statistics settings, for instance 26)). Indeed, n is fixed and the number of knots increases in (21), while increases in (24] 25). Second, extending the proof of 22 to the additive case is not straightforward. since the various spaces of (non-additive) functions and their closures in [21] would have to be appropriately redefined for additive functions. Note also that (21) relies on a convergence proof in (27,, that would need to be extended from the non-additive to the additive setting.
describes the additive GP framework. Sections 3 and present our framework for additive cGPs and introduce the MaxMod algorithm. Section provides the numerical experiments. Appendix A Aaatherrs the technical proofs as well as additional details and experiments.
Paper organization. Section 2
2 Framework on additive Gaussian processes
In additive models, GP priors are placed here over the functions y1...., yd 11 & For i= 1, d, let (Yi(;;;;; v4 E [0, 1]} be a zero-mean GP with covariance function (or kernel) ky- Taking Yı,... , Yu as independent GPs, the process (x); x E [0, 1]d), with æ = (II,..., æd), that results from the addition of Y,...,, , Yd.
is also GP and its kernel is given by
k(x x²) = kn (x1, x4) + ... aa(a)
(2)
In regression tasks, we often train the GP Y to noisy data (æxs yn 1k)sses We denote x for K = 1....,m n and i = 1,..., d, the element corresponding to the i-th input of xg By considering additive
2

---

Figure 1: Additive GP predictions using (left) the unconstrained GP mean, (center) the cGP mode and (right) the cGP mean via HMC (see Section (.3). The constrained model accounts for both componentwise convexity and monotonicity conditions along Iı and 22, respectively.
Gaussian noises Ex. ~ N (0, with "i,...,, Y. then the conditional process Yx(x) := Y (Y + = "I...., (a,) + En = yy. is
En assumed to be independent and independent of
GP-distributed with mean function and covariance function given by 3
u(x) = T (x) + + yo:
cxe, x²) )- k( x²) k- (a) K + 1m1
where yn [8/1..... Un] k(x) Σi (x.) with (x,)) [ka xt )....
... h (x;,,;"")?
and K = -i 1=1 K; with (K,),.f = (xx 2 for S k, f n. The conditional mean and variance v(.) .) are used for predictions and prediction errors, respectively.
2D illustration. Figure shows the prediction an additive GP modeling the function (ar 22) H
4(z) - 0.5)² + 202 The GP is trained with a squared exponential (SE) kernel, k( x' - The covariance parameters 0- ((o (o) f2)), corresponding the variance and length-scale
E of -=,)² x/)2/2(f) and with x2):: (0.5, 0). (0.5, 0.5), (0.5, 1). (0, 0.5), (1,..
parameters (respectively), and the noise variance 7² are estimated via maximum likelihood (28). Although the resulting GP does preserve the additive condition, from Figure (left) we can observe that the quality of the prediction will depend on the availability of data. In our example, we can observe that the GP model does not properly capture the convexity condition along x1.
3 Contributions on additive Gaussian processes under inequality constraints
3.1 Finite-dimensional Gaussian process for fixed subdivisions
In order to satisfy the constraints everywhere (see the next subsection), we introduce the following finite-dimensional GP. For i= 1,..., we consider a one-dimensional subdivision S, (a finite subset of [0, 1] composed of knots) with at least knots at 0 and 1. Throughout Section S, is fixed, but its data-driven selection is studied Section If the number of knots of S; is m, then the total number of knots is given by m m. +... + my. We let S (Si...., Sa). The finite-dimensional
GP, denoted by Ys (,),in is written, for x [0, 1]d.
d m2
Ys(x) = Σ Yrs. (Ix) =
(3)
i=1
j=1
where Ex, (:( with .I...., Y. GPs as in Section 2 and with 0 = ((())
(S.)
(m,)
the knots in S.. We let (5.)
(0) = and (S.) (m:+1) 2 by convention, and Φijj = 0 (,15,)
"))
(8,) (9:)
[0, I is the hat basis function centered at the knot (())) of S, That is, for < < we let
u)
for u S L S v,
uu,o (t)
(w t)
for < t w,
(4)
10
for & & (u. w).
3

---

Y(2)
Yie)
0.50
0.25
0.00
01.00
03.25
1.00 03.00
0..50
0.50
0.75
1.00
Figure 2: (left) 1D illustrations of the hat basis functions in (4) and (center and right) two random realizations from standard GP (solid red line) and the corresponding finite dimensional GP (dashed blue line). For the latter, the location of the knots are plotted by vertical dashed lines.
Observe from 3 that, since Ei.j, for = L...., d and j = L,.... My are Gaussian distributed, then Yi,, is GP with kernel given by
ke( 21))
Dis (z. ) (1 (8))
(5)
j=l==
Moreover, Ys is a GP with kernel k(x x') - ka xt)). Figure shows random realizations from a standard GP and finite-dimensional GP using a SE kernel with (o (1...)) We
remark that consider asymetric basis functions here, which is novelty compared to [5] 6]. that will allow MaxMod (Section 4 to promote non-equispaced designs. The additive finiteddimensional kernel in (5 also seems not to appear the existing literature.
3.2 Satisfying inequality constraints
We consider the componentwise constraints Y,,,, E i I,..., d, where is mne-dimnsiioal function set. In line with [5, we assume that there are convex sets C; C IRm: such that
Yi,,, E E <> &: E: where &. = [S,,, Eimm Furthermore, for i .....,d. d, C, is itself composed by linear inequalities:
(6)
1....,
C, = f-ff CE Vn
m
100
ux
kJJ
j=1
(00 and (. represent the lower and upper where the X A.j encode the linear operations and the bounds. We write the constraints in l S A,c u.
Examples of such constraints are monotonicity and componentwise convexity. For instance, for the case where the function in ( is non-decreasing with respect to each input, then E, is the set of non-decreasing functions on [0, 1] and C, is given by
C, {c "" 2, m, c cj-1 0}.
Hence in this case, Ys is monotonic on the entire [0, 11 if and only if (by definition) each of its additive component Y,s, is monotonic on [o, 1). which happens and only if the m; constraints (8) satisfied. We discuss the example of componentwise convexity in Appendix A.
3.3 Constrained GP predictions
We let Σ; (S, S.) be the m; x m; covariance matrix of &. We consider æ1...., [0, 1] and write O; for the n x m; matrix with element (a, b) given by (x.( Then
Y, Ys(m)), Ys(s.) Σ E.E.
(9)
By considering noisy data (I,, Jr)Isn have the regression conditions Y, t. yn Es. where En ~ (0 r²I,) and is independent from E.....,e. Then given the observations and the constraints, the maximum à posteriori (MAP) estimate, also called the mode function, is given by
d m:
Es(x)
Σ
(x.).
(10)

---

The vector & = ET with E: ( ,,.,,...,,,...... is the mode of the Gaussian distribution N(p: 2e) of the values at the knots conditionally to the observations and truncated from the
constraints 1 S A S Ug,,., la S AaEa Ue in 7
argmin
(c Mc) D Mc).
(11)
4SA.
Above Me [ pedl T is the m 1 vector with block given by
Me = E [E + +
yn
and (E is the m x m matrix with block (i given by
+ Besy = -,= - E ( Σ E, + 1 In D 2
The expressions for fa. and Eg are obtained from the conditional formulas for Gaussian vectors 28 Note that the matrix Eg is not block-diagonal, that is cannot compute the posterior mode along each dimension separately. In Appendix A.2 based on the matrix inversion lemma [28 29 efficient
implementations to speed-up the computation of j., and Σg are given when m
Notice that the covariance parameters of the kernel (5) are estimated via maximum likelihood once the (unconstrained) additive GP framework is established. This is an intermediate step before solving the optimization problem in (I]). Hence, given N.. Ec. Al...., Ad, l...... Ld Ui....y Uds. and O, knowledge, the expression of the mode in [T] does not appear in the existing literature.
convex and can be solved via quadratic programming (5, 30). To the best of our
The cGP mode in (IO) can be used a point estimate of predictions. Since trajectories of .e...... E conditionally on the observations and constraints can be generated via HMC [ЛЛ, they can be used for uncertainty quantification. Furthermore, the mean of the HMC samples can be used for prediction purposes. Continuing the illustrative example in Figure 1 we see the improvement brought by the mode function and the mean of conditional simulations, compared to the unconstrained additive GP model.
Additive MaxMod algorithm
4.1 Squared-norm criterion
Consider an additive cGP model that uses only a subset J C (l....,d) of active variables, with cardinality [J]. We write its (vector of) subdivisions as S (S J). Its mode function Ys is defined similarly as in (I0), from RRIT to I., by, for x = (
P5(T) Σ> (z)))s(:)).
(12)
BEJ 1=1
Adding a new active variable i $ J to J, and allocating it the base (minimal) subdivision S,, (0. }} defines new mode function ss.-- RI55+1 + R. Adding knot / E [0, 1] to the subdivision S, for i* E J also defines new mode function Ys.r.. R 57 - R. Since a new variable or knot increases the computational cost of the model (optimization dimension for the mode computation and sampling dimension for generating conditional trajectories via HMC), it key to quantify its contribution benefit. We measure this benefit by the squared-none -orm modification of the cGP mode
15,1- Imanion (Ts(et) Ys,. (2)) dæ for J,
(13)
Is Jaanee
(Ps(x)) fs.), do for e
(14)
where in the first case we see Ds function of [1++1 variables that does not use variable it.
5

---

4.2
Analytic expressions of the squared-norm criterion
Adding a variable. For a new variable i* & J, the new mode function is
m.
(1))
ΣΣ SEJ
( Σ gng
where (Ers)s, and J), follow from the optimization problem in dimension Eie, my corresponding to (T0) and (II). We let Φ (u) = 11 and Φ5x [u) = for E [0, 1/. Note that even though, for i E J and j = L..... m, Ei. and Ei, correspond to the same basis function Φ..y. they not equal in general the new mode is reestimated with two more knots, which can modify the coefficients of all the knots. Next, we provide an analytic expression of the integral in (13)
Proposition 1 (See proof in Appendix A.3
We have
Is
Σ
m/
(SS)) 33'
Σ
Thes F
(S)
2
+
Th-
Σ IEJ j=1
nis (())) j
Si+ 2
12
IEJ
=1d=11 15-51121
NEJ
where Thi = Eis - Eng, Ti- E - E1,,1, Gi- = Ei,, + Eu F.S) = So
2 (t)d and
F/S)) = So (1)(1) (f)dt with explicit expressions in Lemma 1
Appendix
A.3
The matrices
(E)))) ( Imsss'mm. are 1-band and the computational cost is linear with respect to m ΣiEJ mi
Inserting knot variable. For a new added S... with i E J. the new mode function is
Ps.r = ΣΣ iEJ 1=
E. (
where m, = m, for + 19, m.- m.- + 1, Øi.j Di for + i", and di , obtained from
S.. U [t) as in Lemma As before, this follows from the optimization problem in dimension
Eier m, + 1 corresponding to (10) and (II. Next, we provide an analytic expression of (I4).
Proposition 2 (See proof in Appendix A. Forie E Miit), let S. S, Let S. = S U (t). (S,)) <.... < (( = 1. Let v E (l...., m. 13 be
Recall that the knots in Si. are written 0 =
(1) (1)
such that t (5.)) t (v41
(S..) Then, with linear computational cost with respect to m. Σier mi,
we have
m
The
i ΣΣ iEJ j=l
2
Is,5 = Σ Σ
Th
The Rh
(()) 333'
Σ
Thi
(3)) E(S)
+
Thi
E
(.))
NET
=-'=1
NEJ
1=1
13-1'1s1
where Thi = -G.j - Ev.j. E S.) and E $ are as in Proposition /
E.j = EoS for i Fi, = E.
jor S v. Err = 1-1 2 + 2, and
(S.)) t
(S..)
Er E, (S,,) ) (w+1) (v)
(v+1)
+ Err,41
(v)
((:+) (+:))
(v+1) (v)
4.3
The MaxMod algorithm
Algorithm]] [ [summrrizes the routine of MaxMod. When considering inserting knot i to an active variable i, we promote space filling by adding a reward of the form (t S,), where S, is the current subdivision and where Dtt, S,) is the smallest distance from ł to an element of S, When adding a new variable, we add a fixed reward A' Both A and A' are tuning parameters of the algorithm and allow to promote adding new variables over refining existing ones with new knots, or conversely. Step 3 in Algorithm Tis performed by a grid search, and involves multiple optimizations for computing the new modes in (1T). followed by applications of (14). Step 4 yields single

---

Algorithm 1 MaxMod for additive cGPs
Input parameters: A O, A' 0, d.
I: Initialize MaxMod with the dimension in which the mode MaxModdo maximizes the squared-norm. Sequential procedure: For m E N, m 2 0, do the following.
2: for i = I..... d do
3: if the variable is already active then compute the optimal position of the new knot is E [0,
that maximizes Is.i.t + (t, S,) over [0, 1], with /s,,, as in (14). Denote the resulting
mode as (-(1)
MaxModdmm+1 and the resulting value of Ts.i, I;-
4: else add two knots the boundaries of the selected new active dimension, i.e. (44,1; (2,2)
(0, 1), and denote the resulting mode as y-(1) M+x1 and the resulting value of (13) as I
5: Choose the optimal decision i" E {l,...,d} that maximizes the MaxMod criterion:
i" E arg.aX,e(L....d) (4 + AltE S.) + A'1gg)).
6: Update knots and active variables and set new mode to Y MaxModdmm+1 MaxMed,ma+11
(())
optimization for computing the new mode in followed by an application of (33 For each computation of a new mode, the covariance parameters of the kernels k; for i J. are estimated. For faster implementations, the covariance parameters can be fixed throughout each pass corresponding to fixed value of m (Steps 2 to 6) can be re-estimated every T values of m, for some period E N. Furthermore, steps and 4 can be parallelized and computed in different clusters
5 Numerical experiments
Implementations of the additive cGP framework are based on the R package 1 lineqCPR [1.. Experi- ments throughout this section are executed on an I 1tth Gen Intel(R) Core(TM) i5- 145G7 2.60GHz 1,50 GHz, 16 Gb RAM. Both R codes and notebooks to reproduce some of the numerical results are available in the Github repository: https:://// thub. com/ As training data, we use random Latin hypercube designs (LHDs). Note that using them is not necessary for the implementation of our methodology, but it is recommended when the user is able to design the experiments. Indeed, [32) shows that LH sampling reduces estimation variances compared to simple random sampling, especially for additive functions. Thus, using LHDs enables us to choose minimal design size 2d, corresponding to the number of parameters of the chosen additive GP kernel: an additive Matérn 5/2 kernel with one variance parameter of and one characteristic-lengtll parameter & per dimension (3). We denote O ((o7. (1). f1)...... (o La)). For other applications where n E [d 1, 2d], we may consider the kernel structure in 4 with only one global variance parameter, at the cost of more restricted cGP models. Finally, in the examples that require comparison with non-additive GP models, we use common GP settings: maximin LHDs and n = 10d. The larger value of accounts for the fact that additivity is not prior information of standard GPs.
We compare the quality of the GP predictions in terms of the @² criterion on unobserved data. is defined as Q² 1 SMSE, where SMSE is the standardized mean squared error (28). For noise-free data, Q² is equal to one if predictions are exactly equal to the test data and lower than one otherwise.
5.1 Additive Gaussian processes with monotonicity constraints
5.1.1 Monotonicity SD
We start by an example in small dimensions, in order to compare with non-additive constrained GPs which do not scale high dimensions. We thus consider the additive function given by
2
y(œ)= arctan(5x1) + arctan(2x2) I3 + +2x
(15)
expl- 10 1)
with æ E [0, 1]5. Observe that y is non-decreasing with respect to all its input variables. We evaluate y on a maximin LHD over [0, 1]5 at n 50 locations using (33). In this example, as explained in the introduction of the section, we have chosen maximin LHD rather than random LHD, and

---

Figure 3: Additive GP under monotonicity constraints in 5D with n = 10d. 2D projections of the true profiles and the constrained GP mean predictions are shown in the first and second row, respectively.
Table 1:
Results (mean + one standard deviation over 10 replicates) on the monotonic example in Section 5.1.2 with n 2d. Both computational cost and quality of the cGP predictions (mode and mean) are assessed. For the computation of the cGP mean, 10s (150) HMC samples are used.
CPU Time [s]
Qa [%]
d
cGP mode
cGP mean
GP mean
cGP mode
cGP mean
10
50
0,1 + 0.1
0.1
82.3 + 6.2
83.8 + 4.2
88.1 + 1.7
100
500
0.4 0.1
5.2
89.8 + 1.6
90.7 1.4
91.5 = 1.3
250
1250
4.2
132.3
91.7
92.9 0.6
93.4 0.6
500
2500
37.0 11.4
$156.9 40.5
92.5 0.6
0.5
194.3 0.5
1000
5000
262.4 35.8
110454.3 3399.3
92.6 0.3
94.6 0.2
195.1 0.2
10d rather than n 2d, because we also consider addiiivie GPs for which these settings are recommended. We fix 20 knots dimension, leading to total of m 100 knots.
Figure shows the additive cGP prediction under monotonicity constraints considering 10a HMC samples. Our framework leads to improvements on both CPU time and quality of predictions compared to the non-additive cGP (5). Due computational limitations, the non-additive cGP is performed with a reduced but tractable number of knots per dimension set to = 172 m4 = 4, my 2 and my 6, for total of 768 knots. With this setup, and considering factorial design with 11 test points per dimension, the additive cGP yields Q² = 99.8%, an absolute improvement of 1.3% compared to the non-additive cGP. In term the CPU times, the computation of the cGP mode and cGP mean using the additive framework are obtained in 6s and 0.8s (respectively). significant speed-up compared the non-additive cGP that required 28. 9s and 24.3 minutes.
5.1.2 Monotonicity in hundreds of dimensions
For testing the additive cGP in high dimensions, we consider the target function used in 21
arctan (5
v(x)
d +++1 2
(16)
with E [0, 1]d. This function exhibits decreasing growth rates the index i increases. For different values of d, we assess GPs with and without constraints. We fix 0 = (af 4 = (1, Although 0 can be estimated our focus here is assess both computational cost and quality of cGP predictors. The cGP mean is obtained by averaging 10 HMC samples. We set 5 knots per dimension.
Table | summarizes the CPU times and Q² values of the cGP predictors. The Q² criterion is computed considering 10s evaluations of y based on LHD fixed for the experiments. Results are shown over 10 replicates with different random LHDss From Table it can be observed that the cGPs lead

---

Table 2: Q² Performance of the MaxMod algorithm for the example in Section 5.2.1
with n 10D.
D
d
active dimensions
knots per dimension
Q² (YaaMoa) [%]]
Q² (P (MMaxMad) [%)]
2
2)
(4
99.5
99.8
10
(1,
(5,5,3)
97.8
99.8
(1,2,344,5)
(4, 4, 4,
91.4
99.8
(1,
(5,
99.7
99.8
20
3
(1,2,3)
(4, 4,
99.0
99.9
(1,2,3,4,5 4,
(5, 3, 3,
96.0
99.7
to prediction improvements, with Q² (median) increments between 1.7-5.8%. Although the cGP mean predictor provides Q values above 88%, it becomes expensive when 2 500. On the other hand, the cGP mode yields trade-off between computational cost and quality of prediction.
5.2 MaxMod algorithm
5.2.1 Dimension reduction illustration
We test the capability of MaxMod to account for dimension reduction considering the function
in (16). In addition to (I1,,..,,dd)) we include D d virtual variables, indexed (I FD). which will compose the subset of inactive dimensions. For any combination of D E (10, 20 and d E (,, 3,5), 5). apply Algorithm to find the true active dimensions. For each value of D, evaluate f at maximin LHD with = 10D. Similarly to Section [.1.1 where we compare to non-additive GPs, we have used the common GP settings rather than the ones used here for additive GPs (random LHD, 2d). We compare the accuracy of two modes:
MaxMod:: the mode resulting from the additive eGP and additive MaxMod.
M the mode resulting from the non-additive cGP in with equispaced one- dimensional knots but where the number of knots per dimension is the same as for MaxMoo-
Table 2 shows that MaxMod correctly identifies the dimensions that are actually active where the most variable ones have been refined with more knots. There, MaxMod is considered to converge if the squared-nont criterion is smaller than 5 10 10-4, In terms of the Q² criterion, YaaxMo leads more accurate results compared MaxMod, with Q² 99.7% in all the cases.
5.2.2 Real application: flood study of the Vienne river
The database contains flood study conducted by the French multinational electric utility company EDF in the Vienne river (34] is composed of N x 10 simulations of an output representing the water level and 37 inputs depending on: value of flow upstream, data on the geometry of the bed, and Strickler friction coefficients (see further details in Appendix A. This database was generated as random design, assuming that the inputs are independent and follow specific probability distributions. order to apply our framework, we use componentwise quantile transformations to make the inputs uniform. Such transformations do not modify the additivity and monotonicity properties, but may increase -ineaarrissss With expert knowledge, it possible to identify that H is decreasing along the first 24 input dimensions and increasing along dimension 37. The behavior along inputs 25 to is not clear. Although underlying functions modeling the flood phenomena may be additive, previous sensitivity analysis in (34), done with the full database, has shown that the additive assumption is realistic here. Morcover, it has shown that inputs 11, 35 and explain most of the variance. We aim at obtaining similar conclusions with fraction of the simulation budget, using the information of additivity and monotonocity constraints with respect to dimensions to 24 and 37.
In this experiment, MaxMod is applied for 10 replicates with random training sets and 10 iterations. The settings described at the beginning of Section 5 are followed much possible. the database is fixed, we cannot directly use random LHD of size size n (measured by the Fuclidean norm). The remaining data are used for testing the cGPs. Figure 4 shows the choice of MaxMod per iteration and for the 10 replicates (adding or refining variables).
2d. Thus, we pick the closest subset of

---

Input Dimension
16 11111222110011111122433332322221
1011 12 15111411100
reps
Max.od Iteration
3
5
10
1.32
0.97
Figure 4: Results on the flood application in (.2t.on with n 2d 74. The panels show: (top) the choice made by Max- Mod per iteration and (left) Q² boxplots per iteration of MaxMod. Results are com- puted over 10 random replicates. For the first panel, a bigger and darker square im- plies more repeated choice.
log Q9
1.00
0.90
0.70
0.80
0.40
0.60
0:22
2
MaxMod iteration
The first three iterations of the algorithm show that MaxMod activates mainly dimensions 11, 35 and 37. In the subsequent iterations, MaxMod either refines the these dimensions or activates dimensions 12 and 36. In terms of the quality of predictions (Figure 4), MaxMod leads to Q² > 0.92 after the first three iterations.
Similar results are obtained for n = 3d and n = 4d (see Appendix A.5.2) MaxMod indeed activates dimensions 11, 35 and 37 in the first three iterations and refines them mainly subsequent iterations. This allows to conclude that MaxMod correctly identifies the most relevant input dimensions and that accurate predictions are obtained once those dimensions are activated.
Acknowledgments and Disclosure of Funding
This research conducted within the frame of the ANR GAP Project (RR-21--44000007 and the consortium in Applied Mathematics CIROQUO, gathering partners in technological research and academia in the development of advanced methods for Computer Experiments 10.5:1 1tt.://o.eooo811111111
We thank Bertrand looss and EDF R&D/LNHE for providing the Mascaret test and Sébastien Petit who has performed the computations this model.
References
[1] Hastic and Tibshirani, "Generalized additive models," Statistical Science, vol. 1, no. 3, pp. 297-310, 08 1986.
[2] A. Buja, T. Hastie, and R. Tibshirani, "Linear smoothers and additive models," The Annals of Statistics, vol. 17, no. 2, pp. 453-510, 06 1989.
[3] N. Durrande, D. Ginsbourger, and O. Roustant, "Additive covariance kernels for high- dimensional Gaussian process modeling," Annales de la Faculté de Sciences de Toulouse, vol. 21, no. 3, pp. 481-499, 2012.
[4] D. K. Duvenaud, H. Nickisch, and E. Rasmussen, "Additive Gaussian processes," " in NewrIPS, 2011, pp. 226-234.

---

[5] A. F. López-Lopera, F. Bachoc, N. Durrande, and O. Roustant, "Finite-dimensional Gaus- sian approximation with linear inequality constraints," SIAM/ASA Journal on Uncertainty Quantification, vol. 6, no. 3, pp. 1224-1255, 2018.
[6] H. Maatouk and X. Bay, "Gaussian process emulators for computer experiments with inequality constraints," Mathematical Geosciences, vol. 49, no. 5, pp. 557-582, 2017.
[7] A. F. López-Lopper, S. John, and N. Durrande, "Gaussian process modulated Cox processes under linear inequality constraints," in AISTA7S, 16-18 Apr 2019, pp. 1997-2006.
[8] J. Rithimäki and A. Vehtari, "Gaussian processes with monotonicity information," in AISTATS, 2010, pp. 645-652.
[9] S. Golchi, D. R. Bingham, H. Chipman, and D. A. Campbell, "Monotone emulation of computer experiments," SIAM/ASA Journal on Uncertainty Quantification, vol. 3, no. 1, pp. 370-392, 2015.
[10] J. Wang, J. Cockayne, and C. J. Oates, "A role for symmetry in the Bayesian solution of differential equations," Bayesian Analysis, vol. 15, no. 4, pp. 1057 - 1085, 2020.
[11] A. Pakman and L. Paninski, "Exact Hamiltonian Monte Carlo for truncated multivariate Gaus- sians," " Journal of Computational and Graphical Statistics, vol. 23, no. 2, pp. 518-542, 2014.
[12] I. Delbridge, D. Bindel, and A. G. Wilson, "Randomly projected additive Gaussian processes for regression," in International Conference on Machine Learning, 2020, pp. 2453-2463.
[13] X. Deng, C. D. Lin, K.-W. Liu, and R. K. Rowe, "Additive Gaussian process for computer models with qualitative and quantitative factors," Technometrics, vol. 59, 3, pp. 283-292, 2017.
[14] G. Raskutti, M. J. Wainwright, and B. Yu, "Minimax-oppimal rates for sparse additive models over kernel classes via convex programming," Journal of Machine Learning Research, vol. 13, no. 13, pp. 389-427, 2012.
[15] Da Veiga and A. Marrel, "Gaussian process regression with linear inequality constraints." Reliability Engineering & System Safety, vol. 195, p. 106732, 2020.
[16] A. Cousin, H. Maatouk, and D. Rullière, "Kriging of financial term-structures," European Journal of Operational Research, vol. 255, no. 2, pp. 631-648 2016.
[17] A. F. López-Lopera, F. Bachoc, N. Durrande, J. Rohmer, D. Idier, and O. Roustant, "Approxi- mating Gaussian process emulators with linear inequality constraints and noisy observations via MC and MCMC." in Monte Carlo and Quasi-Monte Carlo Methods. Springer International Publishing, 2020, pp. 363-381.
[18] S. Zhou, P. Giulani, J. Piekarewicz, A. Bhattacharya, and D. Pati, "Reexamining the proton- radius problem using constrained Gaussian processes," Physical Review C, vol. 99, p. 055202, May 2019.
[19] X. Bay, Grammont, and H. Maatouk, "Generalization of the Kimeldorf Wahba correspon- dence for constrained interpolation," Electronic jourraal of statistics, vol. 10, no. 1, pp. 1580- 1595, May 2016.
[20] X. Bay and L. Grammont, "Constrained optimal smoothing and Bayesian estimation," arXiv preprint arXiv:2107.05273 2021.
[21] Bachoc, A. F. López-Loperal and O. Roustant, "Sequential construction and dimension reduction of Gaussian processes under inequality constraints," SIAM Journal Mathematics of Data Science, vol. 4, no. 2, pp. 772-800, 2022.
[22] A. Koppel, K. Zhang, H. Zhu, and T. Başar, "Projected stochastic primal-dual method for constrained online learning with kernels," IEEE Transactions on Signal Processing, vol. 67, no. 10, pp. 2528-2542, 2019.

---

[23] U. Marteau-Ferey, F. Bach, and A. Rudi, "Non-parametric models for non-negative functions," Advances neural information processing systems, vol. 33, pp. 12816-12824 826, 2020.
[24] D. Burt, C. E. Rasmussen, and M. Van Der Wilk, "Rates of convergence for sparse variational Gaussian process regression," in International Conference Machine Learning. PMLR, 2019, pp. 862-871.
[25] A. Koppel, H. Pradhan, and K. Rajawat, "Consistent online Gaussian process regression without the sample complexity bottleneck," Statistics and Computing, vol. 31, no. 6, pp. 1-18, 2021.
[26] A. W. Van der Vaart, Asymptotic statistics. Cambridge university press, 2000, vol. 3.
[27] X. Bay, L.. Grammont, and H. Maatouk, "A new method for interpolating in a convex subset of Hilbert space," Computational Optimization and Applications, vol. 68, no. 1, pp. 95-120, 2017.
[28] C. E. Rasmussen and C. K. I. Williams, Gaussian processes for machine learning (adaptive computation and machine learning). Cambridge, MA: The MIT Press, 2005.
[29] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes in C: The of scientific computing. Cambridge: Cambridge University Press, 1992.
[30] D. Goldfarb and A. Idnani, "Dual and primal-dual methods for solving strictly convex quadratic programs," in Numerical Analysis, 1982, pp. 226-239.
[31] A. F. López-Lopera, lineqGPR: Gaussian process regression models with linear inequality constraints, 2021, R package version 0.3.0. [Online]. Available: https
[32] M. Stein, 'Large sample properties of simulations using Latin hypercube sampling." Technomet- rics, vol. 29, no. 2, pp. 143-151, 1987.
[33] D. Dupuy, C. Helbert, and J. Franco, "DiceDesign and DiceFyal: Two packages for design and analysis of computer experiments," Journal Statistical Software, vol. 65, 11, pp. 1-38, 2015.
[34] S. Petit, F. Zaoui, A.L Popelin, C. Goeury, and N. Goutal, "Couplage entre indices à base de dérivées mode adjoint pour I' l'analyse de sensibilité globale. Application sur le code Mascaret," Sep. 2016, HAL e-prints.
[35] B. Iooss, S. D. Veiga, A. Janon, G. Pujol, with contributions from Baptiste Broto, K. Boumhaout, T. Delage, R. E. Amri, Fruth, L. Gilquin, Guillaume, M. Herin, M. Idrissi, L. Le Gratiet, P. Lemaitre, A. Marrel, A. Meynaoui, B. L. Nelson, F. Monari, R. Oomen, O. Rakovec, B. Ramos, O. Roustant, E. Song. J. Staum, R. Sueur, T. Touati, V. Verges, and F. Weber, sensitivity: Global Sensitivity Analysis of Model Outputs, 2021, R package version 1.27.0. [Online]. Available: https://CRAN.R projeet.org//pcccgg-sssnnni
[36] N. Goutal, J.-M. Lacombe, F. Zaoui, and K. T1-Kadi-Abbderrezzak "MASCARET: A 1-D open-source software for flow hydrodynamic and water quality in open channel networks," in River Flow, R. M. Munoz, Ed., 2012, pp. 1169-1174.
[37] N. Goutal and F. Maurel, "A finite volume solver for ID shallow-water equations applied to an actual river," International Journal for Numerical Methods in Fluids, vol. 38, pp. 1-29, 2002.
12

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] The outline in Section ("Paper organization") provides pointers to where the claimed contributions of the paper are provided.
(b) Did you describe the limitations of your work? [Yes] The paragraph "Range of applicability" in Section i dedicated to this.
(c) Did you discuss any potential negative societal impacts of your work? [N/A] After reading the Ethics Review Guidelines) we conclude that there is no negative societal impacts in submission. Part of our contributions are theoretical, and numerical experiments are proposed for general regression tasks.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes] See previous question.
2. If you including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] The assump- tions of our theoretical results are detailed Sections 3 a and Appendix A
(b) Did you include complete proofs of all theoretical results? [Yes] See Appendix 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] In Section S we included the URL of the repository containing R codes and notebooks of some numerical experiments.
(b) Did you specify all the training details (e.g.. data splits, hyperparameters, how they were chosen)? [Yes] The numerical setups to reproduce the experiments are discussed in each subsection of Section 5 Further details will be given in Jupyter notebooks.
(c) Did you report error bars (e.g.. with respect the random seed after running experi- ments multiple times)? [Yes] See Table[] and Figures 4 5
and
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] The computer specifications of the machine used in all the numerical experiments are provided in Section
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? Yes] We cited the main R
packages we used in our experiments (see references [1] 31 33 35)).
(b) Did you mention the license of the assets? [N/A] codes open-source. The flood database is private.
(c) Did you include any new assets either in the supplemental material or as URL? [Yes] In Scction we provided the URL of Github repository containing all the R codes.
(d) Did you discuss whether and how consent was obtained from people whose data you'l 're using/curating? [Yes] We thanked in the acknowledgments the persons and institutes that provided the flood database in Section 5.2.2
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] The flood database Section 5.2.2 does not have personally identifiable information or offensive content.
5. If you used crowdsoureing or conducted research with human subjects...
(a) Did you include the full text of instructions given participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to 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": "okrs-HKL2B9", "reviewer": "Reviewer_cQK3", "summary": "The authors consider the problem of regression tasks with an additive Gaussian process prior and linear inequality constraint. The authors propose a finite dimensional approximation to the Gaussian process as a linear combination of triangular basis functions with Gaussian weights. The weights are then estimated via solving a quadratic program (for MAP estimation) or approximately sampled with HMC to handle the inequality constraints. Additionally, the authors consider the problem of variable selection in this model, and propose a forward selection method based on the difference between the posterior mode with and without the inclusion of a particular variable (or knot). The authors demonstrate their method for inference and variable selection on synthetic data, that is generated to be additive and monotonic along certain dimensions. Finally, the authors apply the method to a flood dataset generated by a utility company via simulation.", "questions": "1. (This is a comment) I think it would be useful to include a figure showing a.) the hat basis functions, b.) a sample from a 1-dimensional Gaussian process prior with a given kernel (e.g. Matérn 5/2) as well as the sample (using the same realization of Gaussian noise) modified to the finite dimensional model (eqn 5) with a handful of knots.\n\n2. (Lines 96-97) You refer to a function as being monotonic on $[0,1]^d$. I don’t think this is a well-defined notion (since $[0,1]^d$ isn’t totally ordered). What do you mean by this? I assume you mean that restricted to each component it is monotonic, in which case the first equivalence you state is really just a definition. This should be made clear. \n\n3. Equation 8 seems more restrictive than the original formulation in which $\\mathcal{C}_i$ need only be convex. Are they somehow equivalent? If not, it would be best to initially introduce $\\mathcal{C}_i$ as a set defined by finitely many linear inequality to avoid confusion.\n\n4. Notationally, the first equality in (eqn 10) seems to suggest that the MAP estimate can be computed by computing a MAP estimate along each dimension separately and combining these. Is this what is intended? I don’t see why this would be true (and in equation 11 this does not seem to be the case). If it is, could you please give some explanation as to why? If not, consider changing the notation slightly.\n\n5. (Minor, line 116) “et” -> “and”\n\n6. (Minor) Is the information benefit directly tied to ideas from information theory? If not, it may be worth considering renaming.\n\n7. (Minor, Algorithm 1, 5) “D” -> “d”\n\n8. In the flood dataset experiments, you selected your training set via finding a subset as close as possible to a LHD (250-251). Would findings have differed if you had instead selected training data uniformly at random? Similarly, what if $n$ were somewhat larger than $2d$? In particular, is it important to the practical application of the method that the covariate design be something like a LHD? I didn’t see where an assumption like this was made in the discussion of the method, but it seems all of the experiments used this. \n\n9. How robust is the proposed MaxMod algorithm in the presence of model mis-specification? For example, if the function being modelled has non-additive components, can this lead to issues with variable selection?", "limitations": "There is some discussion of limitations in the paper. I think it would be nice to see a synthetic experiment where the assumptions of the model are not precisely true (e.g. the function used to generate the data has some, perhaps small, additive component). I expect most real-world data is not exactly additive, and it would be nice to see how robust the proposed method is to small mis-specifications. ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n- The scope of the paper is well-defined and of interest to those modeling with Gaussian processes. Incorporating prior knowledge such as inequality constraints into these models seems to be an important task. \n- The forward selection algorithm (MaxMod) for selecting which dimensions to include seems appealing for high-dimensional datasets, and the experiments seem to suggest that it can be effective. \n- The writing is generally quite clear.\n\n### Weaknesses\n- I am concerned that the contribution regarding extension of existing work to handle high-dimensional, additive models is over-claimed. It looks to me like the method considered is very similar to the references [5,6], applying the prior works approach to each dimension (separately).If there is a major distinction that I have missed between the approaches, I would appreciate if the authors could highlight this more explicitly. Otherwise, I think this contribution should be contextualized a bit more clearly relative to existing work, particularly in section 3 where it isn't entirely clear what is new and what is a based directly on existing methods.", "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": "DfktHssD1gT", "reviewer": "Reviewer_RpmA", "summary": "The authors propose an additive Gaussian process with support for monotonicity constraints and scalable to high dimensions. In particular, a sequential dimension reduction algorithm is presented which can identify the active input dimensions. The author also provided illustrative results on real-life data.\n\n", "questions": "1. Are the kernel parameters inferred jointly with other parameters in Eq. (11)?\n2. What are the desiderata for the basis functions Eq. (4)?\n", "limitations": "N/A", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "[Strengths]\n1. The paper is well-written. The proposed GP framework is an interesting generalization of additive GPs.\n2. The technical derivation looks solid, and the kernel of the proposed GP has a simple analytical form.\n3. The authors also validated the method on a Vienne river flood data where the proposed MaxMod algorithm correctly activated the relevant input variables.\n\n[Weaknesses]\n1. It's not clear how the kernel parameters are estimated. In Eq. (11), the objective is given by a constrained quadratic optimization. However, this objective could be nonconvex in terms of the kernel parameters. Are there efficient solvers?\n2. The proposed GP framework seems to rely on particular basis functions (Eq. (4)) which have domain [0, 1]. Does this domain impose limitations on the applicability of the GP framework? What are the desiderata for these basis functions?\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": "02RgJJ48ULT", "reviewer": "Reviewer_CWsg", "summary": "This work develops an algorithm for GP inference with constraints through a novel analytical update rule for an additive variant of GP inference. An active set technique for which model points to retain is developed based on Euclidean subspace projections. Numerical results illuminate the merits of the proposed approach.\n\n", "questions": "See response to previous field.", "limitations": "See response to previous field.", "rating": 4, "confidence": 5, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nThe algorithm appears novel to my knowledge, and addresses a fundamental and open problem in additive reparameterizations of Bayesian linear regression with Gaussian processes.\n\nThe 2d example, as well as the extended numerical evaluations, provide substantive evidence that the proposed technique works well In practice for imposing constraints in this setting.\n\n\nWeaknesses:\n\nThe authors should contrast their approach to imposing constraints through Bayesian regression with frequentist analogues. In particular, through a well-known link to kernel ridge regression, one can impose convex or linear constraints, which are very much related to the setting considered here. See, for instance:\n\nA. Koppel, K. Zhang, H. Zhu, and T. M. Basar. ”Projected Stochastic Primal-Dual Method\nfor Constrained Online Learning with Kernels” in IEEE Trans. Signal Processing, May. 2019.\n\nMarteau-Ferey, U., Bach, F., & Rudi, A. (2020). Non-parametric models for non-negative functions. Advances in neural information processing systems, 33, 12816-12826.\n\nAnd follow-on works.\n\nThe authors also consider a way to deal with complexity bottleneck of GP posterior inference through a sequential dimensionality reduction approach with the squared-norm. Substantial complexity reduction of computing an inference is achieved through this active set approach. What I wonder is how is this related or different from a large history on active set techniques for complexity reduction of GP inference. In particular, recent works have  established specific convergence guarantees as a function of model complexity. How does the proposed technique contrast with these performance certificates? Is the squared-norm approach sufficient due to the additive linear nature, but beyond this situation one should consider metrics in distribution space (such as Hellinger or Wasserstein?)\n\nMcAllester DA (1999) \"Pac-bayesian model averaging.\" In: Proceedings of the twelfth annual conference on Computational learning theory, pp 164–170\n\nBurt D, Rasmussen CE, Van Der Wilk M (2019) \"Rates of convergence for sparse variational gaussian process regression.\" In: International Conference on Machine Learning, pp 862–871\n\nA. Koppel, H. Pradhan, and K. Rajawat. “Consistent Online Gaussian Process Regression\nWithout the Sample Complexity Bottleneck,” in Statistics and Computing, Springer, Sept.\n2021\n\nSuch a contrast is not carefully addressed in this manuscript. Specifically, see Table 1 in the last of the preceding list of references.\n\nThe main theoretical results are of a nature of providing analytical expressions for the update rules. While this is credible and interesting, additional justification in the context of what would/would not guarantee constraint satisfaction is missing.\n\nThe numerical results do not really consider strong benchmarks for evaluating the tradeoffs in the constraint satisfaction/model fitness achieved by this method as compared with some alternatives. Moreover, it is unclear how does the active set approach appropriately capture the right statistical properties of the fully dense GP. \n\n\n\nOverall opinion: \n\nThe algorithm is novel, and its update expressions combined with the active set approach to dimensionality reduction, ensures its efficiency. However, these aspects are only addressed heuristically, and no performance certificates are provided. Moreover, numerically no strong baseline is considered to determine whether the proposed technique actually achieves state of the art performance. For these reasons, the work is below the bar.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "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