Fast variable selection makes scalable Gaussian process BSS-ANOVA a speedy and accurate choice for tabular and time series regression
Authors: David Mebane
Kyle Hayes
Ali Baheri
Venue: NeurIPS 2022 Submitted (2022)
Keywords: scalable Gaussian process
timeseries
tabular data

Abstract: Many approaches for scalable GPs have focused on using a subset of data as inducing points. Another promising approach is the Karhunen-Loève (KL) decomposition, in which the GP kernel is represented by a set of basis functions which are the eigenfunctions of the kernel operator. Such kernels have the potential to be very fast, and do not depend on the selection of a reduced set of inducing points. However KL decompositions lead to high dimensionality, and variable selection thus becomes paramount. This paper reports a new method of forward variable selection, enabled by the ordered nature of the basis functions in the KL expansion of the Bayesian Smoothing Spline ANOVA kernel (BSS-ANOVA), coupled with fast Gibbs sampling in a fully Bayesian approach. It quickly and effectively limits the number of terms, yielding a method with competitive accuracies, training and inference times for tabular datasets of low feature set dimensionality. The new algorithm determines how high the orders of included terms should reach, balancing model fidelity with model complexity using $L^0$ penalties inherent in Bayesian and Akaike information criteria. The inference speed and accuracy makes the method especially useful for modeling dynamic systems, by modeling the derivative in a dynamic system as a static problem, then integrating the learned dynamics using a high-order scheme. The methods are demonstrated on two dynamic datasets: a 'Susceptible, Infected, Recovered' (SIR) toy problem, with the transmissibility used as forcing function, along with the experimental 'Cascaded Tanks' benchmark dataset. Comparisons on the static prediction of derivatives are made with a random forest (RF), a residual neural network (ResNet), and the Orthogonal Additive Kernel (OAK) inducing points scalable GP, while for the timeseries prediction comparisons are made with LSTM and GRU recurrent neural networks (RNNs). The GP outperforms the RF and ResNet on the static estimation, and is comparable to OAK. In dynamic systems modeling it outperforms both RNNs, while performing many orders of magnitude fewer calculations. For the SIR test, which involved prediction for a set of forcing functions qualitatively different from those appearing in the training set, BSS-ANOVA captured the correct dynamics while the neural networks failed to do so.

Fast variable selection makes Karhunen-Loève decomposed Gaussian process BSS-ANOVA a speedy and accurate choice for dynamic systems identification
Anonymous Author(s)
Affiliation
Address
email
Abstract
Many approaches for scalable GPs have focused on using a subset of data as inducing points. Another promising approach is the Karhunen-Loève (KL) decom- position, which the GP kernel is represented by a set of basis functions which are the eigenfunctions of the kernel operator. Such kernels have the potential to be very fast, and do not depend on the selection of reduced set of inducing points. However KL decompositions lead to high dimensionality, and variable selection thus becomes paramount. This paper reports a new method of forward variable selection, enabled by the ordered nature of the basis functions in the KI. expansion of the Bayesian Smoothing Spline ANOVA kernel (BSS-ANOVA), coupled with fast Gibbs sampling the number of terms, yielding a method with competitive accuracies, training and inference times for tabular datasets of low feature set dimensionality. The new gorithm determines how high the orders of included terms should reach, balancing model fidelity with model complexity using L penalties inherent in Bayesian and Akaike information criteria. The inference speed and accuracy makes the method especially useful for modeling dynamic systems, by modeling the derivative in a dynamic system as static problem, then integrating the learned dynamics using a high-order scheme. The methods are demonstrated on two dynamic datasets: 'Susceptible, Infected, Recovered' (SIR) toy problem, with the transmissibility used as forcing function, along with the experimental 'Cascaded Tanks' benchmark dataset. Comparisons
3
fully Bayesian approach. It quickly and effectively limits
10
11
12
13
14
16
18
20
the static prediction of derivatives are made with ran- dom forest (RF), a residual neural network (ResNet), and the Orthogonal Additive Kernel (OAK) inducing points scalable GP, while for the timeseries prediction com- parisons are made with LSTM and GRU recurrent neural networks (RNNs). The GP outperforms the RF and ResNet the static estimation, and comparable OAK. In dynamic systems modeling outperforms both RNNs, while performing many orders of magnitude fewer calculations. For the SIR test, which involved prediction for a set of forcing functions qualitatively different from those appearing the training set, BSS-ANOVA captured the correct dynamics while the neural networks failed to do so.
22
24
25
26
28
29
30
Karhunen-Loove decomposed Gaussian processes
31
32
1.1
Gaussian process fundamentals
33
Gaussian processes (GPs) are stochastic functions that are engines for nonparametric regression.
34
Initially developed for modeling and interpolation in geographic information systems datasets,
Submitted to 36th Conference on Neural Information Processing Systems (NeurIPS 2022). Do not distribute.

---

35 applications have multiplied across many fields of data science. A key advantage of the GP is its 36 broad, continuous nonparametric support and the frequent amenability of different GP kernels to
37 precise analysis.
GP is Gaussian in that it is covariance model linking pairs of points on functional draws. As such 39 a GP is completely described by a mean function (often zero in the prior) and covariance kernel. The
40
most famous and perhaps simplest of the covariance kernels the squared exponential:
x' C
<(a, x') = exp
(I)
41 where the sill 5² and range & parameters determine the scale and smoothness of the draws. In typical implementation modeling a static dataset Z, the statistical model
42
Z 6(x)s² {) +
43 with an observation error process, first used to infer the hyperparameters, after which predictions conditioned on the training dataset can be made. The draws on the squared exponential GP limiting case of the Matérn covariance family are infinitely differentiable.
44
45
From a practical standpoint the training of the above GP is O(N*). requiring a Cholesky decom- position of the full covariance matrix. This limits the use of the GP to moderately-sized datasets, generally of thousand instances or fewer.
46
48
1.2
Scalable Gaussian processes with inducing points
50
Liu et al. 20200) provide a thorough overview of efforts that aim to improve scalability while maintaining prediction accuracy using global kernel approximations derived in some sense from 52 set of M << N inducing points [Chalupka et al.] [2013] Caanoonero-Canded and Rasmussen, 2005
51
Deisenroth and Ng 2015 Rasmussen and Ghahramani 2001 Wang et al. 2022]. Generally the goal 54 is to approximate the full-rank kernel matrix with local approximations. Of particular note a (NN method that directly estimates the covariance with training and inference times that limits the increase
56 in M for large N developed by Wilson et al.) (2015). Some methods employ ANOVA decompositions
to the full kernel which break out contributions in terms of features and their combinations:
6(xxx' = Σ (1 2 +
Ka(2,, (Z), +...
(3)
=1)=
58 which presents opportunities for variable selection [Duvenaud et al. 2011]: of particular note is the recent Orthogonal Additive Kernel (OAK) which orthogonalizes the kernels in (3) in order to
59
60
minimize overlap between main effects and higher-order interactions [Lu al. 2022
61
1.3
Karhunen-Loève decomposition and BSS-ANOVA
Another approach to scalability in GPs that is distictive to the inducing points approach is the Karhunen-Loövv (KL) expansion, in which the kernel is expressed in terms of a sum over its
63
eigenfunctions:
8(z;3) ~ MVN(0, = Σ 3(d((x)
(4)
where
xx(x) ()uu(())
(5)
I k(xxx') )uil l''dd' = Ax((x))
(6)
B NO,
86 Such methods have the potential to fast: O(NP) in training and per point for inference, where P is the number of terms in the expansion. However such kernels have not been the subject of much
68 research in machine learning contexts generally. The main issues are tractable calculation the basis
functions } and dimensionality issues [Greengard and O'Neil 2021

---

70 In 2009 Reich et al. (009 introduced the Bayesian Smoothing Spline ANOVA (BSS-ANOVA)
71 kernel, which is subject first an ANOVA decomposition, followed by KL decomposition. The core of the BSS-ANOVA kernel is:
1 24
(2 Ki (z. x') Bi (())) + B4 x'
(8)
where B& is the htth Bernoulli polynomial, defined by the generating function
te et -
t i!
Σ (x)x) 1=
(9)
yielding
Bi (x) = x 1
(10)
2
B2(x) x² x ++
(11)
6
1
B1(x) x² 2x
(12)
30
75 This kernel is effectively sum of non-stationary quadratic response surface corresponding to the first two terms in 8 and a stationary deviation (the final term). As in (3) covariances for higher-order interactions are constructed with dyadic products of the main effect covariance:
76
77
<2[[2 mxl, (ayg xAl) = in (ry, )r (2k, a/k)
(13)
78 and on for higher-order interactions. Terms are then multiplied by scaling hyperparameters and added together to produce the full kernel:
79
प: + of Σ
K1,
02
K2,aj
(14)
The kernel so constructed is supported by second-order Sobolev space RRich et al. 2009 which is 81 very broad and dense set of continuous functions.
80
82 Building the kernel in this fashion effectively addresses the problem of generating the eigenfunctions from the KL decomposition: because all of the terms in (14) based on the generative kernel (8),
84 The KL decomposition of [4] will depend only on eigenfunctions of K1. Additionally all input features normalized to an [0, 1] interval (we restrict the discussion to continuous input features 86 for now), then it is only necessary to compute a single set of basis functions {Ø:}- The decomposed BSS-ANOVA GP is written:
(xx; 3) = 3o +
Baok (x:)
+
Σ
Σ
Bi (x:)) (x))dr +...
(15)
1=1
j=++1
k=
88
Given the assumption
0000 - -
... =.. 02,2
(16)
then the priors for the coefficients 3 are iid normal
B.
(17)
Following |Reich et al. 2009 generate the set {o,) by producing K1 for dense grid consisting
91 of 500 intervals on [0, 1]. eigendecompose and to cubic splines. Figure Tssows the first basis functions. These basis functions are nonparametric, pairwise orthogonal, and ordered: note the increase in frequency and decrease in amplitude as the orders increase.
88
Variable selection
94
95 It's clear from 115) that the number of terms in the expansion can increase rapidly, even for low- 96 dimensional input spaces. A key component of applying the GP to a modeling problem is thus the selection of terms. Effectively we seek to minimize the objective function
(()) Z 8(x; 3)]1 ((3)
(18)
where Ç is penalty function which leads to sufficiently sparse solution.

---

x10
0.5
0.2
0.02
0.01
M
0.01
W
M
0.1
0
0
0.01
-0.02
-0.5
0.02
0.01
0 0,5
1
O
0.5 O 0.5
0 0.5
0 0.5
1
0
0.5
Figure 1: The first six basis functions of the KL-decomposed BSS-ANOVA kernel. The basis is nonparametric, spectral, pairwise orthogonal and ordered.
99
2.1
Indicator variable methods
100
Reich et al. 2009) took a hierarchical Bayesian approach to the problem, estimating a separate
101
variance - for each term in the expansion. which turn expressed terms of an indicator variable with a Bernoulli prior. This approach, like other 'indicator variable' methods, accomplishes the variable selection and the training simultancously and comprchensively, at the cost of requiring
102
103
104
large number of variables in the prior model and computationally onerous Markov chain Monte Carlo (MCMC) sampling procedure.
105
106
Other sparse optimization methods such as ridge regression or LASSO share the limitation that many high-order terms must be included in the initial model before downselection occurs.
107
2.2
Forward variable selection
109
The ordered and orthogonal nature of the basis functions suggests a forward variable selection approach. Rewriting the model (15) for basis function set of maximum order q.
a-1 n ΣΣX
6(x: 3) = Bo
Ba
Σ
Bo +... (19)
111 then considering a model building procedure which increases q stepwise starting with 9 = reveals that each subsequent step adds main effect terms (each depending single input), (2q
Iwo-way interactions, and (5) 3(q 1)² + 3(q 1) 1] three-way interactions. As the model
113
114 order increases the 2 truncation error for the full kernel decreases as (for the case of single
input) (Greengad and O' Neilll 2021
1/2
[|n(r, x')
Se (x)) <
Σ
x3)
(20)
1=++1
Since the eigenvalues of the BSS-ANOVA kernel decomposition decrease quickly with increasing order, an approach to the optimization problem (18) focusing on low-order models will sacrifice little in the way of accuracy while realizing significant advantages in computing time.
117
119
The design and implementation of such an approach is the main contribution of this work. approaches the optimization of (18) with an iterative process, finding the most efficient truncation of the system while evaluating the cost function only for candidate models withffeer terms than the
121
optimum truncation. The method is fully Bayesian, with fast Gibbs sampling procedure at its core.
123
As such the form of the cost function is also Bayesian in nature, taking the form of the Bayesian or Akaike information criteria (BIC/AIC), which incorporate 7 penalties.
125
2.2.1
Gibbs sampling
126
Given a statistical model
(21) with white noise observation error, and given truncation the KL expansion (15). the model is linear in the coefficients and Gibbs sampling can be used to estimate parameters in fully Bayesian methodology.
z; 8, (x;; (3 +
128
129
the variance of the observation emrr o2, with inverse gamma prior o² TG(a, b), with a and the shape and scale parameters, respectively; and if 72 has inverse gamma prior 7² ~ 1G br).
130

---

then an iterative Gibbs sampler can be devised such that for fixed {o²,>²), 8 ~ (VN Σ), with
132
== (x"x + +1/-1)-1x²
(22)
133
Σ= 3=x2(""X1//1
(23)
where X E NxP is matrix constructed from the basis functions, whose rows correspond to instances and columns to terms in the expansion. For fixed {3, T }, 0² ~ IG(a", 6"), with
135
=-=++++ N/2 P/2
(24)
136
- 6' = b+ 2 (M 3) 3 -) - x
(25)
2
For fixed (3,02), 72 ~ IG b7), with
137
if = a, + -P-1)/2
(26)
138
6 = b, 20
B
(27)
139
This algorithm is implemented in the routine 'gibbs_Xie Xinss in the supplement.
2.2.2
Optimization algorithm
The algorithm constructs models with terms having up to three-way interactions. Terms added in stages labeled by an integer "index" that initializes at 1. At each stage, a series of substages cycle through all permutations of basis function orders that sum up that stage's index. Stage 1 adds only first order main effects. Stage 2 adds second order main effects and first order two way interactions corresponding di (z (2zy in two separate substages. The substages always occur such that terms involving lower-order basis functions (for example in the case of stage 2, this is the first order wvo-way interactions) come first. Each substage adds at once all combinations of inputs and all permutations among each combination, such that each substage adds (*) terms for two-way interactions and (3) terms for three-way way interactions. Then the sampler is called and the BIC or AIC calculated. Because there is not monotonic decrease increase pattern for the objective function, a "tolerance" setting controls how many substages the algorithm can iterate through without finding new minimum BIC or AIC before it terminates. The algorithm returns the optimum model.
142
144
146
147
148
149
150
152
153
This algorithm appears in the routine 'emulator Xin' in the supplement.
154
3
Experiments: Dynamic system identification
155
3.1
Procedure
156
BSS-ANOVA regression- as is the case for other GPs - is most effective for tabular datasets with continuous inputs and targets of moderate dimensionality. This suggests an application in dynamic systems identification. Indeed BSSVANOVA GPs have been utilized as components of other models ("intrusively") for this purpose in a number applications (Bhat al., 2017 Lei et al.) 2019 Ostace et al. 2020]. We demonstrate here that they may also be used directly to identity dynamics in more
157
158
159
180
161
general cases, without the aid of an accompanying model.
162
The procedure is a concurrent one, in that derivatives estimated from the datasets are modeled directly using BSS-ANOVA with forward variable selection, using the coneurrent values of the system states and other inputs; for example a two-state system is modeled using two separate GPs:
163
184
f 81 (x1, u)
(28)
i (2,1, u)
(29)
165
The identified system is then integrated to yield predictions with uncertainty.
186
The procedure was demonstrated on two nonlinear dynamic datasets: a synthetic dataset derived from the susceptible, infected, recovered model (SIR model) for infectious discase, and the 'Cascaded
5

---

Algorithm BSS-ANOVA forward variable selection algorithm
1: procedure FWDVARSELECTGE 2, Φ, tol, h)
D. is a vector of hyperparameters.
2: i
3: count=
4: while count tol do
5:
if ind is new then
6:
Find all combinations of integers that sum up to ind, ordering them by the maximum integer appearing in each combination, with the lowest maximum first.
7:
8:
1. Select the next combination in the and place the integers into vector with many elements as there are model inputs, buffering out with zeroes.
9:
10:
2. Produce matrix Ma the rows of which contain all permutations of that vector. & Each row corresponds to a term in the GP expansion. 3. Produce input matrix Xa where columns are model terms and rows are experiments,
11:
13:
for all terms appearing in Ma.
4. Recursively concatenate: X (X),, M [M; Ma]
15:
5. 3, BIC gibbs pioss_xi(XX,,,, Z, o. h)
if the BIC is a minimum for all models then
17:
save the model count
19:
else
count++
21:
if all combinations for ind have been utilized then
22:
ind++
23:
Return M. 3, BIC
Tanks' experimental benchmark dataset. In both cases comparisons were made to long short term memory (LSTM) and gated recurrent unit (GRU) neural netowrks for timeseries prediction. In the case of the cascaded tanks benchmark comparisons were made against random forest (RF), residual neural network (ResNet) and the state-of-the-ar OAK inducing points scalable GP [L et al. 2022] for the static derivative estimation problem.
169
171
173
3.2
Experimental benchmark: Cascaded tanks
174
The cascaded tanks nonlinear benchmark dataset is an experimental nonlinear dynamic system WWggren and Schoukens 2013]. The experiment consists of set of two tanks and reservoir of water. An upper tank is filled by pump from the reservoir. An outlet in the upper tank empties into the lower tank, which in turn empties through an outlet back into the reservoir. A signal sent the pump serves as the forcing function for the system, with the tank water level heights the two states of the system.
175
178
177
178
179
180
We first compared the performance of BSS-ANOVA with RF, ResNet and OAK static regressors. Derivatives were calculated via direct finite differences for the relatively noise-free dataset, yielding 10000 instances. Each method was trained on concurrent values of both states and the forcing function for each derivative. For the GP we used hyperparameters of
182
183
1000, b 1.001, a, 4 and b, 55 for hy and 69.1 for h2, with tolerances of 3 for h and 5 for h2. and the AIC discriminator.
184
185
Of 2000 draws the first 1000 were discarded. Only two-way interactions were required. For the RF
196
100 trees were used with a leaf size of 5. The ResNet had depth of (filter sizes ranging from 16 to 64) and in between each fully connected layer is batch normalization and relu layer. The mini batch size is 16, initial learn rate is 0.001, the data was shuffled every epoch for total of 30 epochs,
187
188
and the validation frequency was 1000. OAK was applied at a maximum dimension of 3 and with the default value of 200 inducing points. The 5-fold cross-validated results appear in Table T OAK performed best for both outputs, followed closely by BSS-ANOVA. Both GPs outperformed the RF and the ResNet by clear margins.
190
192
Timeseries predictions follow for the GP via a -odeerder Runge-Kutta integration routine. These were compared with LSTM and GRU recurrent neural networks (RNNs). For the LSTM there was one L.T layer and total of 128 hidden layers, the data was shuffled every epoch for a maximum of 125 epochs, verbose equal to 0, and the sequence was padded to the left. The GRU had one
195
195

---

Table 1: Cascaded tanks 5-fold cross validated accuracies: derivatives
Method
h1 (MAE/10
h2 (MAE/10
OAK
17:4.77
3612,4
BSS-ANOVA
1826.5
3943.6
ResNet
36+14
61+11
RF
30:9.44
49+:499
Table 2: Cascaded tanks 5-fold cross validated accuracies: timeseries
Method
h1 (MAE/MAPE)
l/2 (MAE/MAPE)
BSS-ANOVA
(1167+00.0382/4667-11
0.1577 5.99:
LSTM
0.2345 9.46:14.77 0.3243#0.1092 12.16:55.02
0.229610.0378 9.382 0.248110.0402 9.89-3..00
GRU
GRU layer and 150 total hidden layers, the data shuffled every epoch for total of 150 epochs, verbose was equal to zero and the sequence was padded to the left. The 5-fold cross-validated results (datapoints were not randomized before creating the folds
198
to preserve the timeseries order) appear in Table 2 BSS-ANOVA is most accurate, followed by the LSTM and the GRU. Figure shows the predictions of the GP and the LSTM for the upper tank for of the test folds. The GP predictions are superior near the sharp inflection and critical points where nonlinearities are strongest. provided the LSTT and GRU start-up
200
202
Note that the first 50 points of each test set, which set in the prediction phase, were removed from the calculation of error for both methods.
While it is reasonable to expect that OAK with 200 inducing points would outperform BSS-ANOVA in the time integration, it was not practical to make this comparison for reasons computing time. A comparison with reduced number of inducing points and increased time step in the integrator was made- results are discussed in section 3.4
206
208
209
3.3
Synthetic benchmark: Susceptible, infected, recovered model
The susceptible, infected, recovered model (SIR model) is common simulation for infectious disease. Though there are several versions, the simplest is three states, only two of which are independent. The system is written
211
$ = -BIS/N
(30)
BIS/NP -1
(31)
R=
(32)
6
GP Bounds
Twat Data
Test Data LSSTM Modee
Mean Prediction
2.4
2.6
2.8
2.45
2.8
Time 69)
x100
Time [s]
Figure 2: (a) BSS-ANOVA and (b) LSTM predictions vs. test set data for the water level height in tank of the cascaded tanks dataset. Shaded regions in (a) are 95% confidence bounds as estimated from a draw of 40 curves.

---

213
where S(t) is the susceptible population, [t) the infected, R(t) the recovered, (t)) is the transmis- sibility rate (which we utilize as forcing function), Y the recovery rate (which we leave fixed at
214
215
0.5) and Np is the total population. Because Np is fixed and S + 7 + R = Ny, only two states are independent, so the system dynamics can be captured by modeling only two of the three. We chose
216
217
(t) and R(t).
218
The training data consists of 58 curves. All curves in the training set have fixed B value ranging from 0.5 to 9, in six intervals of 1.7. For each value of B there are 8-10 siumulations corresponding to different initial conditions designed in such a way to provide coverage of the state space. (Exact
219
220
221
initial conditions used appear in the supplement.) Each simulation used Np = 1000.
The test data consists of 24 curves, each of which features temporally changing transmissibility BI). There are three initial Bo values: 1.35, 4.75 and 8.1.. For each starting point there are two types of transmissibility curves: a ramp and sinusoid. The Bg 1.35 and Bo 4.75 starting points have ramps with a positive slope of 1, while the By = 8.15 curves have a slope of -1. All ramps run from 1=0 to (=4, where they level off. The sinusoids have amplitudes between 0.5 and and period of 1.
223
225
228
227
Hyperparameters for BSS-ANOVA were: a a, 4 for both states, 6+_xt 8.95 and b,., 72.1, while 6, = 1.25 and bR = 20. 2000 draws were taken and the first 1000 discarded. The tolerance was 6. Hyperparameters for the LSTM and GRU were the same for the Cascaded Tanks.
229
231
partial display of the results are shown in Figures 3 for BSS-ANOVA and for for the GRU, which was the better performing of the two neural nets on this dataset. For the GP, the total test set MAF for R, corresponding to MAPEs of 8 8.99++4.92 for T and 2.80 for R. Statistics were not calculated for the GRU as it failed to replicate the dynamics
232
233
was 5. 5.2739++4.01. for 7 and 1 1.8345++11.33
234
235
in most test cases and was obviously inferior in a quantitative sense in every instance, as shown in Figure
500
1000
CP
400
GP Meen Model
Tout Duta
Recovend i-d
60
Tost Data
GP Mean Model
800
600
400
300
200
100
IGGP Beundh Test
10
15
-200
5
10
15
0
10
20
Time Weees)
Time (Werio))
Time (Weelel
(a)
(b)
(c)
1000
Boundo GP Meen Modell
Irheted indeoaae
2010
600
700-
200
cm
800
150
Paconend
700
100
HGP Dounce Text Date
600
HGG Boursts
Text Data
Model
3000
10
5
10
500
8
10
Time Wessi
Time (Wookk))
Time (Neoks)
(d)
(e)
(I)
Figure 3: BSS-ANOVA results for 3 curves in the test set: (a)-(b) sine wave transmissibility with low initial infections; (c)-(d) sine wave transmissibility with moderate initial infections; (e)-(f) ramp transmissibility. Shaded regions are 95% confidence bounds for the predictions as estimated from draw of 40 curves.
3.4 Training and inference times
238 Training and inference times for BSS-ANOVA were fast, with mean total train time of 6.3 seconds
239 for the cascaded tanks and 10.8 seconds for the SIR, with 8,000 and 20,000 training data points, respectively, on a 2019 6-core i7 processor with 16 GB of RAM. The routines were implemented in MATLAB, but not parallelized or optimized for speed. Models for h contain between 23 and 41
8

---

Intected holihaals
400
300
Tost Dusa
Recovred Indeedddddduu Indhicuale
Test Dora
300
300
800
Test -FU Mode
100
0
Time (Wesko))
Time (Weslo))
Time
(a)
(b)
(c)
400
1100
Reovened wrltouak
Infened trduidate
CIT MooM
1000
Recoend
1000
600
700
100
Test Duca
CFUU Mode!
10
30
5
0
5
Time
(Wookol
(Hoolol
(d)
(c)
(f)
Figure 4: GRU results for 3 curves in the test set: (a)-bb) sine wave transmissibility with low initial infections; (c)-(d) sine wave transmissibility with moderate initial infections; (e)-(f) ramp transmissibility
terms, while hy has between 38 and 57 terms. Prediction times for 2000 static points for the cascaded tanks averages 0.5437 s, and the time for evaluating integrals over the test set averages 20.22 For the SIR model the I model had 81 terms and the R model 9 terms, with mean integration time of 5.3 S. Analyses have shown that the rate limiting step in BSS-ANOVA build algorithms are the O(N native
242
245
construction of the X matrix from the inputs and basis functions. The neural networks MATLAB functions, parallelized and optimized for speed. Nonetheless train times were considerably longer, with mean train times of 130s for the ResNet and 175 and 123 s, respectively, for training the LSTM and GRU for the cascaded tanks. This is to be expected given that the number of weights in the neural nets are the order of 10a.
247
249
250
251
was not feasible to integrate OAK at the level of 200 inducing points to the same standard as that of BSS-ANOVA because of time considerations. A reduced set of 40 inducing points yielded accuracies in the static estimation problem that were approximately the same as BSS-ANOVA. A reduced time step (500 vs. 20,000 integration steps) brought the integration time down to 51 minutes for OAK, with MAE/MAPE of 1554//6. for h and 0.2378/9. II for h2. Reducing the integration step to the same level BSS-ANOVE (where could expect comparable integration accuracies) would require approximately 33 hours.
253
255
257
4
Limitations and future work
258
259
The two examples presented in this paper were the only two attempted for dynamic systems identi- fication. Other benchmark dynamic systems, especially those chaotic in nature, will be examined in future work. Desipite the advance in variable selection represented by this routine, datasets with higher dimensionalities in the feature space are more challenging and require additional methods for variable selection, which are already in development. More experiments and comparisons will be performed for dynamic systems as well, with larger datasets and more difficult identification problems. Like any GP BSSAANOVA is inaccurate in extrapolation: when test set inputs exceed bounds of the training set the resulting extrapolation sometimes causes instabilities causing the integration procedure to fail. These failures were eliminated by preventing extrapolation even in instances where the inputs exceeded the bounds, but more stable extrapolation strategies will possibly become necessary for longer prediction windows where extrapolation is unavoidable.
260
281
262
263
264
285
286
268

---

References
270
271
K. Sham Bhat, David S. Mebane, Priyadarshi Mahapatra, and Curtis B. Storlie. Upscaling uncertainty with dynamic discrepancy for mulli-scale carbon capture system. Journal of the American Statistical Association, (11520):1453-166) 2017.
273
274
Krzysztof Chalupka, Christopher KI Williams, and lain Murray. A framework for evaluating approximation methods for Gaussian process regression. Journal of Machine Learning Research, 14:333-350, 2013.
275
278
277
Mare Deisenroth and Jun Wei Ng. Distributed Gaussian processes. In International Conference on Machine Learning, pages 1481-1490. PMLR, 2015.
278
279
David Duvenaud, Hans Nickisch, and Carl Edward Rasmussen. Additive Gaussian processes. Advances in Neural Information Processing Systems 24 (NIPS 2011), 24, 2011.
Philip Greengard and Michael O'Neil. Efficient reduced-rank methods for Gaussian process with eigenfunction expansions. arXiv preprint arXiv:2108.05922 2021.
282
283
Yinkai Lei, Tian-Le Chen, David S. Mebane, and You-Hai Wen. Reduced-order model for mi- crostructure evolution prediction in the electrodes of solid oxide fuell cell with dynamic discrepancy reduced modeling, Journal of Power Sources, 416:37-44 49,, MAR 15 2019.
284
285
Haitao Liu, Yew-Soon Ong, Xiaobo Shen, and Jianfei Cai. When Gaussian process meets big data: review of scalable GPs, IEEE transactions on neural networks and learning systems, (11)) 4405-4423, 2020.
287
Xiaoyu Lu, Alexis Boukouvalas, and James Hansman. Additive Gaussian processes revisited.
290
International Conference on Machine Learning, pages 14358-14383, PMLR, 2022.
281
Anca Ostace, Keenan X. Kocan, David S. Mebane, J. Pieter Schmal, and Debangsu Bhattacharyya. Probabilistic model building with uncertainty quantification and propagation for dynamic fixed bed CO2 capture process. Energy & Fuels, (122):211-25332 FEB 2020.
293
294
Joaquin Ceinonero-Candell and Carl Edward Rasmussen. A unifying view of sparse approximate Gaussian process regression. The Journal of Machine Learning Rexearch, 6: 699999-9959, 2005.
295
296
Carl Rasmussen and Zoubin Ghahramani. Infinite mixtures Gaussian process experts. Advances in neural information processing systems, 14, 2001.
Brian J. Reich, Curtis B. Storlie, and Howard D. Bondell. Variable selection in Bayesian smoothing spline ANOVA models: Application to deterministic computer codes. Technometrics, 51)) 110-120, MAY 2009.
299
300
301
Liwei Wang, Suraj Yerramilli, Akshay Iyer, Daniel Apley, Ping Zhu, and Wei Chen. Scalable Gaussian processes for data-driven design using big data with categorical factors. Journal of Mechanical Design, 144(2), 2022.
302
303
304
Torbjorn Wigren and Johan Schoukens. Three free data sets for development and benchmarking in nonlinear system identification. In 2013 European Control Conference (ECC). pages 2933-2938. European Control Conference (ECC), ETH Zurich, Zurich, Switzerland, JUL 17-19, 2013.
306
307
Andrew Gordon Wilson, Christoph Dann, and Hannes Nickisch. Thoughts on massively scalable Gaussian processes. arXiv preprint arXiv: 2015.
308
Checklist
310
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
312
313
(b) Did you describe the limitations of your work? [Yes] see Section
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
10

---

315
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
316
2. If you including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A]
(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 ex- perimental results (either in the supplemental material or as a URL)? [Yes] See the supplement
322
323
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Sections 3.2 and 3.3
325
326
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes] See Sections 3..tiionss32aaann and 3.3
327
328
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See Section 3.4
330
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
331
(a) If your work uses existing assets, did you cite the creators? [Yes] See Section 3.2
332
(b) Did you mention the license of the assets? [N/A]
333
(c) Did you include any new assets either in the supplemental material or as URL? [Yes] See the supplement
334
335
(d) Did you discuss whether and how consent was obtained from people whose data you're 're using/curating? [N/A]
337
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
338
339
5. If you used crowdsourcing or conducted research with human subjects...
340
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you deseribe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
343
344
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
345


Reviews:
[{"review_id": "y23PVFjYMS9", "reviewer": "Reviewer_Yi9J", "summary": "The paper proposes a scalable GP based on Karhunen-Loeve expansions of kernels. The paper takes a fully Bayesian approach where Gibbs sampling is used to infer the kernel parameters, therefore, leading to forward variable selections.\n", "questions": "How to scale $\\mu$ and $\\Sigma$ in Eq. (10-11) which involve matrix inversion? If they are evaluated on 500 data points, such estimations are just local. \n", "limitations": "Besides the presentation issues, I think the paper does not provide an extensive comparison with existing models. For example,an obvious baseline is sparse inducing GPs. \n\nI believe the paper needs an additional revision to be ready for publication.\n", "rating": 3, "confidence": 4, "soundness": 1, "presentation": 1, "contribution": 1, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The variable selection using GP is an interesting problem.\n\nI think the weaknesses outweigh the strengths, include\n\n- The paper seems to be not well-structured, lacking a proper writing of what the problem that the paper tries to address here is.  It would be better to have an introduction section that would clarify those for us as readers.\n- The paper does not provide sufficient descriptions or analyses for the method, for example, lacking literature reviews and no comparison with alternative methods. \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": "W1-GD1btXlB", "reviewer": "Reviewer_abxx", "summary": "The authors describe a method for variable selection in a GP type model for large datasets called BSS-ANOVA. This (existing) approach uses a kernel that has a basis decomposition where the basis functions can be precomputed and the training complexity becomes O(NP) where P is the number of terms from the kernel decomposition used. The papers’ contribution is on selecting the components in this decomposition with non-zero coefficients using the natural order of complexity and a forward selection paradigm. The coefficients and an information criterion is calculated via Gibbs sampling.", "questions": "What is the worst-case complexity of this algorithm? I.e. BSS-ANOVA has complexity O(NP) if I pre-specify the number of terms P, but your algorithm appears to determine that number dynamically. Would it naturally stop at N^2, could it be unbounded or is guaranteed to be better than N^2?\n\nIs the Gibbs sampling algorithm described in Section 1.3 an original contribution of this paper or is it based on previous work? Does the crucial assumption (8) that enables the use of Gibbs sampling reduce the expressiveness of this model in a meaningful way or can it be shown that this assumption is wlog?\n\nPlease change the formula references to (N) to avoid confusion and clean up the citations ([Reich et al., 2009] vs Reich et al. [2009]).\n\nAre there any other variable selection methods for BSS-ANOVA in the literature? Is there a way to use e.g. L1 penalties to find solutions faster?", "limitations": "I can not see any negative societal impact of this work. The authors do address some limitations in that the paper underperforms on some datasets.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper describes a novel approach to fit a GP type model for large datasets and I think as such has a place in the literature.\n\nI would like the authors to be more clear about the impacts of the proposed Gibbs sampler, as described below, i.e. both with regards to the novelty of it and whether it changes the fitted model. My assumption is that this method is new, but if it already exists in the literature, I do not think the use of an information criterion for variable selection alone is contribution enough to warrant acceptance.", "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": "VDLdbBSuJc_", "reviewer": "Reviewer_jXjU", "summary": "Regarding the scalable Karhunen-Loeve decomposed kernel BSS-ANOVA for Gaussian process modeling, this paper proposed a BIC/AIC assisted forward variable selection method to have fast and accurate predictions for large tabular datasets. The performance of the proposed model was tested on two datasets.", "questions": "1) Additive Gaussian processes within the ANOVA formulation for high-dimensional inputs have been studied for a long history, see a recent paper\n[1] Lu X, Boukouvalas A, Hensman J. Additive Gaussian Processes Revisited. arXiv preprint arXiv:2206.09861, 2022.\nThese however have not been discussed and compared in this paper. \n2) Besides, the usage of BIC/AIC to perform truncation is not a new idea, what is the difference to previous works? Please highlight it.\n3) What benefits can we have when using the ANOVA type decomposition for GP? Reducing computational budget? Figuring out import variables? The authors should make it clear rather than just showcasing the accuracy of time series regression.\n4) The proposed model is suggested to be compared to the ANOVA-type GP competitors on more benchmarks to offer solid conclusions.\n", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Pros:\n1) A forward variable selection method via BIC/AIC is developed in the ANOVA regime for GP modeling;\n2) The new method has been applied to time series regression.\n\nCons:\n1) Benefits brought by forward variable selection are unclear;\n2) Differences to existing works have not been well discussed;\n3) Insufficient numerical experiments to showcase the superiority of proposed model.\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": "NdJeKCUoCFK", "reviewer": "Reviewer_j4ha", "summary": "The authors develop a Bayesian approach to  the Karhunen-Loéve (KL) decomposed kernel BSS-ANOVA", "questions": "See strengths and weaknesses", "limitations": "See strengths and weaknesses", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 1, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths and Weaknesses:\nKey weakness.  Most people in the GP community will not be familiar with the methods used.  So we need a _much_ clearer background, and perhaps more standard notation:\n* Why are we using curly $\\theta$ rather than $x$ for the input locations?\n* Define kth Bernoulli polynomial.\n* Tell us the shape of $\\mathcal{B}_1$ etc.\n* \"the full kernel for a system with n inputs is written $\\delta \\sim MVN(0, \\Gamma)$\" is a very strange statement.  The kernel is written $\\Gamma$, not $\\delta \\sim MVN(0, \\Gamma)$?  Where did $\\delta$ come from? Why are we using MVN rather than $\\mathcal{N}$            (especially given that it can cause confusion with a matrix-variate Normal distribution.\n* $\\Gamma = \\sigma_0^2 \\tau_0^2 + ...$ makes no sense.  $\\Gamma$ is a covariance matrix, while $\\sigma_0^2 \\tau_0^2$ are scalars?\n* $\\Gamma_{1, i}$ and $\\Gamma_{2, ij}$ aren't defined.\n* Where are the eigenvalues and eigenfunctions coming from?\n\nWithout these clarifications I do not believe most of the target audience will be able to understand the paper.", "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: Reject