Dynamic Online Ensembles of Basis Expansions
Authors: Daniel Waxman
Petar Djuric
Venue: Accepted by TMLR ()

Abstract: Practical Bayesian learning often requires (1) online inference, (2) dynamic models, and (3) ensembling over multiple different models. Recent advances have shown how to use random feature approximations to achieve scalable, online ensembling of Gaussian processes with desirable theoretical properties and fruitful applications. One key to these methods' success is the inclusion of a random walk on the model parameters, which makes models dynamic. We show that these methods can be generalized easily to any basis expansion model and that using alternative basis expansions, such as Hilbert space Gaussian processes, often results in better performance. To simplify the process of choosing a specific basis expansion, our method's generality also allows the ensembling of several entirely different models, for example, a Gaussian process and polynomial regression. Finally, we propose a novel method to ensemble static and dynamic models together.

Published in Transactions on Machine Learning Research (05/2024)
Dynamic Online Ensembles of Basis Expansions
Daniel Waxman
Department of Electrical and Computer Engineering Stony Brook University
Petar M. Djurić
petor. detaordjjuricssttnnnoook..k
Deportment Electrical and Computer Engineering Stony Brook Universiby
Reviewed on OpenReview: https: // openreview. net/
met/ forum? id- a==aVVii aVOzWH1Nc5
Abstract
Practical Bayesian learning often requires (1) online inference, (2) dynamic models, and (3) ensembling over multiple different models. Recent advances have shown how to use random feature approximations to achieve scalable, online ensembling of Gaussian processes with desirable theoretical properties and fruitful applications. One key to these methods' success is the inclusion of a random walk on the model parameters, which makes models dynamic. We show that these methods can be generalized casily to any basis expansion model and that using allernative basis expansions, such as Hilbert space Gaussian processes, often results in better performance. To simplify the process choosing a specific basis expansion, our mothod's generality also allows the ensembling of several entirely different models, for example, Gaussian process and polynomial regression. Finally, we propose a novel method to ensemble static and dynamic models together.
1 Introduction
Many machine learning applications require real-time, online processing of data, a setting that often neces- sitates significant modifications to standard methods. Online adaptations of many different methods have
been derived, including kernel machines (Kivinen et al. 2004 Lu et al. 2016), (kernel) least-squares (Kay
1993 Engel et al., 2004). and Gaussian processes (Gijsberts & Metta 2013 Stanton ct al. 2021). ). Online
learning has also been studied extensively from the optimization perspective (Oraoonal 2019
Hazan 2022
Online learning may be additionally complicated by the need for model selection, as it rarely clear at the
beginning of the learning process which model will perform the best. One option in this case is train several models in parallel and ensemble them. In Bayesian context, Bayesian model averaging (BMA) (Hoeting et al. 1999) has long been used to ensemble online models, and works by weighting each "expert" model by its evidence.
More recently, Lu et al. 2022) showed how to adapt BMA to online Gaussian processes (GPs) in a technique they call incremental ensembles of GPs (TE-GPs) (sec also Lu et al. (2020)). GPs are a flexible, non- parametric tool of Bayesian machine learning that enjoys universal approximation properties and principled uncertainty estimates (Rasmussen & Williams) 2005). Through random Fourier feature (RFF) approx- imation to Gaussian processes (Rahimi & Recht 2007 Lázaro-Credilla et al.) 2010), online learning can be performed (Gijsberts & Metta 2013), with closed-form Bayesian model averaging updates and tractable regret. analysis (Lu et. al.) (022)..
illusraate the benefits of allowing random walks on model parameters, which they call dynamic IE-GP G-GPss (DIE-GPs) This can dramatically increase performance when the learning task changes slightly over time. For example, when tested the SARCOS dataset
In addition to an online ensemble GPs, Lu al. 2022

---

Published in Transactions on Machine Learning Research (05/2024)
Rasmuusse & Williams 2005). which consists of learning the inverse dynamics of robotic arm, DIE-GPs significantly outperform TF-GPs LLu et al. 2022
Taking more global view, the ensembling of random feature GPs, as introduced by IE-GPs, has proven ftexible and effective; extensions to the framework include Gaussian process state-space models (Liu et al. 2022), deep Gaussian processes (uuu et al. 2023a); ; and graph learning (Polyzos al. 2021 p. Together with
its extensions, DIE-GPs have been successfully used for Bayesian optimization (Lu al. 2023 and causal inference (Liu et al. 2023b
However, reliance on the RFF approximation means that TE-GPs inherit the weaknesses of random feature GPs. Namely, the RFF approximation is a direct Monte Carlo approximation of the Wiener-Khinchin integral and therefore suffers harshly from the curse of dimensionality. We show in our results that on several real-world datasets, (D) provide similar or worse performance than simpler models, including online Bayesian lincar regression and onc-layer RBF networks.
In this paper, we introduce online ensembles of basis expansions (OEBEs); generalization of IE-GPs that break their reliance RFF GPs and improve performance over several real datasets. In particular, we make the following contributions:
We note that the derivation of DIE-GPs does not rely at all on the RFF approximation other than being a linear basis expansion. In fact, the same derivations and code can be reused to ensemble arbitrary Bayesian lincar models with any design matrix. This allows not only the ensembling of the same lype of model logether, but also the ensembling of many different basis expansions (e.g. B-splines, one-layer RBF networks, and more).
2. We argue that, if GP regression is of interest, GP with generalized additive model (GAM) structure is often more appropriate. To this end, we use GAM Hilbert space Gaussian processes (HSGPs) (Solin & Särkkä, 2020), which can be vicwed as a quadrature rule of the same integral that the RFF approximation largets by direct Monte Carlo. Besides theoretical arguments, empirical results suggest HSGPs converge the true approximated GP more quickly (in terms of the number of basis functions) than RFF GPs (Riutort-Mayol et al. 2023 analysis.
We provide a similar empirical
3. We introduce novel method to combine static and dynamic models, enabling the use of principled posteriors of static methods when appropriate and expanding expressiveness of dynamic methods otherwise. We show that this method is necessary by providing a constructive example on real data where the naive approach to ensembling static and dynamic methods fails.
We provide Jax/Objax (Bradbury et al. 2018 Objax Developers 2020 code at httpss://ww: github. com
that only requires the user to specify the design matrix, with several choices alrcady implemented.
The rest of this paper is structured as follows: in Section 2 we review background concepts in linear basis expansions, GP regression, spectral approximations of GPs, and BMA. These concepts are operationalized in Section 3 whore we introduce the OEBEs and several extensions, including applications to non-Gaussian likelihoods, and offer some brief theoretical remarks. We provide additional practical comments regarding the construction of OEBEs, including discussion on the contents of an ensemble and how to ensemble static and dynamic models in Section 4 The proposed models are tested empirically in Section 5 Finally, we offer some concluding remarks and future directions in Section
2 Background
In this section, we begin by reviewing Bayesian linear regression and linear basis expansions (Section 2.1 We then discuss GP regression (Section 2.2 These topics are then connected via the sparse spectrum and Hilbert space approximations, which utilize lincar basis expansions to approximate GPs (Section 2.3 Finally, we discuss BMA (Section 2.4), which can be used to combine distributions from several different models.
2

---

Published in Transactions on Machine Learning Research (05/2024)
2.1 Linear Basis Expansions
A common setting for Bayesian linear regression is with an iid. Gaussian prior on D-dimensional weight vector and Li.d. Gaussian noise, i.c.,
== x
(1)
0 ~ N(0,0,00īD) e N,
(2)
where is a x N matrix [x] ... xx] of imput data, y is an N-dimensional output vector, M stands for a Gaussian distribution, is an N-dimensional Gaussian vector, ID and In identity matrices of dimensions D and N, respectively, and of and o² are variances corresponding to 0 and respectively. If an intercept term is desired, X can be augmented with row of ones and 0 taken to be (D + 1)-dimensional.
By swapping each X, with its mapping under a function Φ: RD - RF, the Bayesian linear regression of Eq. (T] becomes linear basis expansion, with Eqs. to being replaced by
y + T
(4)
0 N(0:0,
(5)
E ~ N,
(6)
where & is an F' N matrix, [o(xi) ((xx)]. The functions ex(-) = proje (Ø(-) are known basis functions. We refer to + as the "design matrix" of X with its rows representing the basis functions.
2.1.1 Predictive Distribution and Online Estimation
is conjugate and Gaussian, which means that efficient, closed- form posterior updates and predictive distributions are available (sce, c.g., (Bishop) 2007 Chapter 3)). In particular, let μ4 and Eg denote the mean and covariance function of the posterior p x, = =N Eg).
The model described by Eqs.
Then the posterior predictive distribution of N. test output/inpat points, and X., respectively, is given by another Gaussian, i.e.
Kyyly,X,.,,, = -N(yylmy, Ey.),
(7)
pyy = - 140
(8)
Ey. D Eo IN
(9)
where T. the design matrix of X..
Using the predictive distribution in Eq. online regression can be performed by what is often called "predicting then correcting" in the filtering literature (see, e.g., (Särkkä & Svensson 2023)). To set up notation, let X1æt and Yıa denote the first t data points, Φ denote the design matrix of X1s, and no,1, Eo. denote the parameters of the posterior p(0)xi:e; Wia) (OOOpo.. Ee.). Then given Xx+1 and Bt41, the posterior can be updated by the following two steps (ääkkkä & Svensson 2023):
(Predict) Obtain 1y+1 and oyt1 from Eqs. and 9
(Corvect) Update the posterior to pl (1 (1.+1) -N(0ppo.41 Eo,+1), where
Eg. (0.10(x1110011 Mytti)
p4.11 poot +
(10)
03__ Eg Ea (xxx))x+))xxxx Eg.
Eg.41 Σo,-
(11)
y++1
3

---

Published in Transactions on Machine Learning Research (05/2024)
2.1.2 Fitting Hyperparameters with Empirical Bayes
The approach of Bayesian lincar basis expansions is attractive in that it is simple, has closed-form posterior updates, and as we will see, can be quile effective. However, there are several hyperparaneters to the model, namely the parameters of the priors og and o². It is also typically the case that ((-) depends on additional hyperparameters, which refer to as up.
A conventional approach to address the lack of knowledge of these parameters is to maximize the marginal likelihood; we can immediately identify this as the prior predictive, and using Eq. calculate it as of
ol. x)) = (y((o, $ T + INo?) This technique is known by several names, including "empir-
ical Bayes," "type-II maximum likelihood estimation," and "evidence maximization" " (Theodoridis 2020 p. 611). Empirical Bayes has a long history of being effectively used for Bayesian linear basis expansions and Gaussian process regression [Rasmussen & Williams 2005 Bishop 2007).
Another approach is rooted in the assumption of conjugate priors. For example, closed-form posterior updates are still possible if the joint prior of the linear parameters A and the unknown variance of follows a multivariate normal inverse Gamma distribution. In this case, it can be demonstrated that the marginal predictive distribution of y. (i.e., the predictive distribution after integrating out the unknowns O and o²) is a Student's t-distribution (Zellner 1996). This analysis has previously been used for a generalizations of IE-GPs in Liu et al. (2023a).
2.1.3 Examples of Basis Functions
Many examples of basis functions exist, and we only provide a few here.
The simplest possible basis function is the identity function Φ(x) = x, in which case we recover ordinary linear regression. Slightly more sophisticated is the polynomial regression, where for scalar æ, ( = * In either case, an intercept term ag can be added by appending a row of ones to O. Multivariable polynomial regression is formulated similarly with polynomials over Xp, but becomes unwieldy with growing dimension, the number of terms grows exponentially.
From polynomial regression come splines, which are piecewise polynomials between locations x XX called knots dd Boor 1978). From the perspective of a generalized linear model, a design matrix can be calculated for so-called basis splines (B-splines) Ide Boor 1972 however, like polynomials, generalizing splines to high-dimensional settings is expensive and not straightforward Due to their added complexity, B-splines are not included in the experiments of the current work, but we note they form an interesting topic for research.
Another common example is the set of (Gaussian) radial basis functions (RBFs). They depend on location parameters μn and length scale parameter E, i.e.,
& (x:x = exp ( (- (x 11x) [xx) - 1&5))
2
The locations Me may be optimized taken to be fixed by some initialization (for example, via a K-means clustering of the data).
In what follows, we will introduce GPs and two basis expansions designed to approximate them.
2.2 Gaussian Process Regression
Gaussian processos (GPs) are an infinite-dimensional generalization of the lincar basis expansion model given above, in the sense that there may be an infinite number of basis functions dx(-) (Rasmussen & Williams, 2005). Inference is then made possible by the kernel trick [Rasmussen & Williams 2005 pp. 12), where basis functions are instead described by a kernel function x(:,)) which summarizes the pairwise relationships of inputs. Indeed, each of the linear basis expansion models explored above is special case of Gaussian
For numerical stability, when F N. one should use the equivalent expression given in Bishop 2007 of F and N, caleulations should be handled in the log domain using Cholesky decompositions.
16.1 Regardless
4

---

Published in Transactions on Machine Learning Research (05/2024)
process and can be written with the kernel trick using the notion of the equivalent kernel Bishop 2007 pp. 159).
While this description of a GP is correct, it is more conventional to describe GPs as a prior distribution over functions. To this end, we must specify a mean function m(-) in addition to the kernel k(-;-); then GP regression (GPR) consists of estimating f()) in the problem
y= f(x) +c,
(12)
f(x) ~ (P(m x(x,x))), rid N(c,0,0)))
(13)
(14)
2.2.1 Predictive Distribution
Estimating the value of f on test points in closed form is simple using the conjugate model of Eqs.
12
to [4]. As before, let X. be a set of test points and s(X.) be lest values to be inferred. Additionally, let K denote the matrix whose entries are given by (X(X1,)))1 = : Then the predictive distribution is again given by a Gaussian distribution,
p(f(X.)) X, X.) (1), E,.),
(15)
1y mX.) K (X,X) ( (X)),
(16)
Σ,. = K(X, x.) K(X,, X K(X,X) + o.1 -1 K(,,X,).
(17)
In contrast. to the basis expansion case, the values f(X.) are noiseless evaluations of the random function f- From Eq. [2), it is easy to see that the predictive distribution of the noisy quantity y. can be obtained by simply adding IN, of to the predictive covariance Eq. (15)
These expressions are simple, but Eqs. 16 and 17 require solving the linear systems
[KXX X) + o - m (X) and (K(X,X)+ K(,,
By using the Cholesky decomposition of [KXXXX) + +o2 and solving the associated triangular systems, the calculation of Eqs. (16) and 177 can be made numerically stable and memory efficient relative to the calculation of inverses, but incurs 3(N3) time complexity (Rasmussen & Williams 2005 App. A). The poor computational scaling with respect to N of solving these linear systems is a major obstacle to implementing GPs in practice and will motivate our discussion of scalable approximations. But first, we will specify the mean and covariance functions that are most commonly used.
2.2.2 The Choice of Mean and Kernel
Without loss of generality, the mean function m(-) can be chosen be identically zero. The kernel function x(;-) can in principle be any positive semi-dcfinite function, but few choices of kernel dominate the literalure. The squared exponential (SE) kernel and the Matérn-3/2 kernels are likely the most common
Rasmussen & Williams 2005 Chapter 4):
KSE(X, x') = oxeep ( (-1xx (11x-xx11
282
3/ (x, x') 0² 1+ /3Ixx 1x11
x'1l
(v3 x x'I
exp
where l and of are the length scale and process scale hyperparameters, and are both positive. The process scale of is directly related prior variance oe of the linear basis expansion case, and in subsequent sections, will refer to of as dg when convenient.
We can additionally allow for the length scale to vary in each dimension; this modification to kernels is called automotic relevance determination (ARD). By expressing
D
(Ii,d - r2a)²
2
(x,,
Σ
=|[x1 xellaingg((h....f)))
B
d=1
5

---

Published in Transactions on Machine Learning Research (05/2024)
the SE-ARD and Matérn-8/2-ARD kernels can be expressed follows:
RSE-ARD (x,x) of MMatém-3/22ARDD ( x = = /3RARD) exp ( VBAARD)
exp exp 0.5r AND),
To fully specify a kernel function, optimal values of o, and & must also be determined. In this regard, it is common to employ the empirical Bayes approach discussed in Section 2.1.2 Again similar to the linear basis expansion case, the marginal likelihood is casily seen to be ((yoo, (( X) +
2.3 Sparse Spectrum and Hilbert Space Approximations
While GPs are flexible tool that makes minimal assumptions on the data-geneee sitaageeeraaingg process, calculating the predictive distribution requires the (symbolic) inversion of an N x N matrix, alternatively, the solving of the corresponding linear system. This is ( operation, and quickly becomes unwieldy as N grows, A typical way to deal with this is the inclusion of inducing points, a set of M < N points that are chosen to summarize the data basis expansion via the spectra of the GPs, as is done in the random Fourier feature (RFF) (Lázaro-Gredilla
Caanonerr-Canndlle & Rasmussen 2005). An alternative approach is to take a finite
et al. 2010 and Hilbert space GP (HSGP) (Solin Särkkä 2020) approximations.
2.3.1 Random Fourier Feature Gaussian Processes
The main idea of RFF GPs is to take Monte Carlo approximation of the Wiener-Khinchin integral, which describe below. The kernel of GP is said to be stationary if depends only on the difference xx'. In this case, Bochner's theorem guarantees spectral representation of the kernel in the form of the Wiener- Khinchin integral (Rasmussen & Williams 2005 pp. 82):
x(x,x) = / exp(2Tis' (x - x')) du(s)
(18)
for some positive finite measure μ(()).
In the case that p(-) admits a density S(s) w.r.t. the Lebesgue measure on RD, it known as the power spectral density (PSD) and the Wiener-Khinchin theorem says given by the Fourier transform of g(x -x) x') (Rasmussen & Williams 2005 pp. 82). Then we can sample Sp....; S(s)
and approximate Eq. 18] by
6(x, x)x
Σ exp Smppsmis (x-
(19)
By noting the symmetry of the PSD, and to ensure vallev kernels, we in practice sample F/2 spectral points Si...., Sp/2 and use both sn and Sm Eq. (19) Seizaro-Grrdille et al. 2010 This results in a linear basis expansion model la Eq. 44 with the basis Tunction
2 2 [sin (1) cos(x Ts) ... sin(x s5/2) cos(x 8F/22)1 cos(xTsF/2)) e F... V
(xx)
(20)
2.3.2 Hilbert Space Gaussian Processes
Hilbert space GPs (Solin & Särkkä 2020 take different approach, which can be viewed as quadrature
rule on the same integral of Eq. (18). Their derivation is rooted in targeting the power spectral density and approximating it via Hilbert space methods from partial differential equations and can be found in Solin &
Särkkä 2020 We begin by specifying v(-) in the scalar case, 1
For the Hilbert space approximation be taken, a compact subset of Ω R must be specified where the approximation is valid. Riutort- SiuuorttMMayol al. (0023) explore the case where S2 = [-L, L, L}, with emaxp |z], representing the maximum absolute value takes over the training data D scaled by value <>> 1.

---

Published in Transactions on Machine Learning Research (05/2024)
The resulting approximation is given by the basis functions (Solin & Särkkä 2020
A:T
sin
KIT x + L))
2n (x(x)) S
A == R,..,,F F
(21)
2L
VL
where we recall S(s) to be the PSD of the kernel.
The multivariable case written similarly (Solin & Särkkä 2020
However, the approximation suffers from exponential growth, in the sense that using Pi...., : FD basis functions for each inpul dimension results in a F1- F2..... FD dimensional embedding. As a result, the approximation is typically more efficient than the exact GP only when D 3 (Riutort-Mayol et al. 2023
In low dimensions (D S 3), the HSGPs empirically converge to the true GP much more quickly than RFF GPs (Solin & Särkkä 2020 Riutort-Mayol et al.) 2023). While HSGPs are not tractable in high dimensions, we note that RFF GPs also suffer harshly from the curse of dimensionality, which is typical in Monte Carlo we argue why better performance in the univariate case sufficiently attractive.
sampling. In Section 4.3
2.4 Bayesian Model Averaging
Bayesian model averaging long-standing technique to combine Bayesian models (Hoeting et al. 1999) With sequence of models M,,..., MM, BMA computes the distribution p(y(D) by process that simulates the marginalization of the choice of model, i.c.
(lD Σ Pi(MMilippllD,MM,n
(22)
The resulting distribution can be understood as a mixture model with mixands p(y MM)) and weights Wm PP)(M.DD).
BMA is simple and asymptotically optimal when the true data-generating process is one of the models Mi...., (Hoeting et al. 1999). This is on account of "collapsing" in the limit of infinite data, in the sense that Pr(MM D) approaches for the Mm, which minimizes a statistical distance to the true data-generating distribution YYa et al. 2018
3 Dynamic Online Ensembles of Basis Expansions
In this section. we re-introduce the work of Lu el al. (022 with no explicit reference to GPs. This comprises an online learning and ensembling step (Section 3.1 as well as an adjustment for problems where dynamies (Section 3.2) or switching (Section 3.3) are present. Finally, we discuss adaptations to non- non-Ganssian likelihoods (Section 3.4 and theoretical regret analysis (Section 3.5
3.1 Online Ensembles Basis Expansions
The Online Ensemble of Basis Expansions (OEBE) consists of learners with basis expansions (f(m)()) and parameters ((m) for m = ,,...,M, M, and a meta-learner based on BMA. The meta-learner initializes weights we to a uniform prior over models, and then learning begins.
A ffter the retrieval of new data pair (x++1. M++1), each learner predicts p(yc+1 (xa Vias x++1) via Eqs. The BMA weights are updated follows:
(m) (m)) (mm) X1: V1a, X141).
w wi
(23)
Note that, if weight (m) (numerically) collapses to O at time i, will remain at

---

Published in Transactions on Machine Learning Research (05/2024)
Next, posteriors p( p(0 X V1:+1) are calculated via Eqs. (10 and Finally, we note that the minimum mean squared error (MMSE) estimate of yt+1 is given by
Pyr+1 Σ m
W "211"
(m) (m)
(24)
and the variance of this estimate by
2 go+1
Σ ((m) (m) (Myee (mm))) )²)
(25)
The OEBE contains as special case the IE-GP if (.(m)()) is given by Eq. (20 for every m. However, the functions ((()(() be of any type of basis expansion, including different types of basis expansions for different values of m.
Note that while each model is linear in ((m), each ((m) corresponds to a different, (nonlinear) basis expansion (.)) Therefore each ((m) represents a different quantity, perhaps even defined on a different space. particular, even with the MMSE estimate Eq. 24 the resulting estimate a mixture of several different nonlinear functions of X,.
We remark that with the MMSE estimate, the resulting prediction is lincar in the "stacked" vector 0 (0(2) ... ((m)), but this does not reflect how each model is trained. For example, an ensemble consisting of the basis expansions
e (x)
and ( [x²]]
is not equivalent to the quadratic basis expansion model. Depending on the basis expansions used, this may result in less expressive model (as in the example above), but it also lessens the computational burden of empirical Bayes significantly, which scales cubically with F.
3.2 Adding a Random Walk on the Parameters
If the dynamics of a problem are expected to change in the dataset, we can add random walk on the parameters 0. For the RFF GP case, this was referred to as a dynamic JE-GP (DIE-GP), so we follow in calling the generalization a dynamic OEBE (DOEBE). To differentiate OEBE from DOEBE, we will often refer OEBE models as "static," juxtaposed to "dynamic."
Incorporating dynamics can be done simply by injecting noise term (m) with variance (m) on the
posterior (mm) prior to prediction.
p(0 yia) = N(((m)) (m)
Mo E + (m) I).
(26)
Note that this slightly differs from the formulation in Lu et al. (2022), which does not allow (m)² to vary with We somewhat generalize to allow OEBE and DOEBE models to be ensembled together, in
Section 4.2 Additionally, note that OEBE the special case where noise degenerate, i.e., when ( 0.
From signal processing perspective, each model in the DOEBE can be viewed as the Kalman filter formu- lation for a linear basis expansion with drift, with Eq. 26 being a closed form solution of the Chapman- 2023 pp. 34). Accordingly, the mth learner corresponds to a
Kolmogorov equation (Särkkä & Svensson
state space model
(fm) (m)
0 0!m) + Erw.t
(m)
1+1
y o (m) (xx) 0 (m)
(m)
+
IIn principle, we could record the entire Gaussian mixture model the estimate. Because of the rapid collapse of BMA single mixand, daes not much change results, and simplifies computations.
8

---

Published in Transactions on Machine Learning Research (05/2024)
i1+11 Qw,
SDOEBE
((m)
ml))
(m)
DOEBE
rw..
Ve ((m) (xx)
To.
+ Ee
p(ye+1 X1:t YI:t, xx+1) Em m0 (m p(m) (ye+1 X1st VIst, X
Vias Wi+1 o wi+11 (m)) (ve+1 Y x++1)
OFBE
Figure 1: The hierarchy of the proposed models. All models have the same distribution for ys conditioned on 0,, while DOEBE adds a random walk to parameters, and SDOEBE adds a random walk BMA weights.
(m) where e(m) iid ~ (( 0, orw (m)2 more complex, i.c.,
I). The state-space formulation shows how this random walk may be even
0 0 =0 (m + 1/22..
where L/2 is the Cholesky factor of positive semi-definite matrix and E, is white Gaussian noise. While this more general form may be useful some applications, particularly when domain expertise is available, we choose to focus on the special case L L=a(m) (m) I in the current work for simplicity.
From this state-space perspective, we can also explore a formal justification for placing a random walk on the parameters. Namely, if the estimation task is expected to change slightly over time, it is useful to model this change drift in the model parameters. From statistical perspective, the posterior distribution of Eq. (26) merely appears as a "more conservative" estimate of the posterior, but this is beneficial in the recursive setting: since we use the posterior at time as the prior for time t t- 1, a more conservative posterior at time corresponds exactly to a weaker prior at time + 1. This allows for better adaptation to new data by "forgetting" old data.
As the weight update step remains unchanged, the DOEBE inherits the weight collapse property of the OEBE. Likewise, the DOEBE contains the DIE-GP as a special case when the noise hyperparameters oiw are all the same and every basis function Φf (-) is given by Eq. (20
3.3 Adding a Random Walk to the BMA Weights
Allowing for a random walk on the BMA weights is also a simple addition, which becomes important for ensembling dynamic and static models. Analogous the Gaussian random walk on 0, dynamics can be introduced to w, via discrete-time Markov chain, specified by stochastic matrix Q:
M
Wr+1
Qwi
Σ
m wi
(m²)
(27)
m²=1
We follow Lu et al. (2022) in calling this a switching (dynamic) OEBE (S(D)OEBE).
The term "switching" is suggestive of applications where regime switching (i.e., sudden changes from generating model to another) may be present. Additionally, Lu et al. (2022) show that adding switching allows for regret analysis in adversarial settings.
However, the more important, property introduced by Eq.. (27 to the current work is that the revival of collapsed weights is possible [u et al. 2022). In particular, note that even if (mm)) has (mumerically) collapsed to 0, it can be revived so long as (mmmw(!") (mm)) is positive for some m'. This allows us to preserve the full ensemble structure after the numerical collapse of the BMA weights, which often occurs in the first few thousand samples.
The OEBE, DOEBE, and SDOEBE models, as well as the relationships between them, are summarized in Fig. 1

---

Published in Transactions on Machine Learning Research (05/2024)
3.4 Inference With Non- Con-Gaus Likelihoods
One seeming weakness of the DOEBE approach is the necessity of Gaussian likelihoods. However, by casting the DOEBE linear Gaussian state-space model, we are able to use several approximations for arbitrary non- non-Gaussian likelihoods.
The simplest such approximation is the Laplace approximation (Kass et al. 1991), which is used for TF- GPs. The Laplace approximation takes a second-order Taylor series approximation of the log-posterior log l(g! [X11, centered around the MAP estimate of 0, which in turn is the solution to the optimization problem
((, EMMP = arg max p(0 (T) X1a.N.a).
(28)
The resulting approximation is a Gaussian distribution whose covariance is given by the inverse of the Hessian matrix for log.(()."") X1:1 Y1:), i.c.,
1 = ((m) log p(0 (m) (1:t,
(m) O.tpaae
06m))
gim)
(29)
By noting the linearity of the Laplacian and factorizing the posterior log p(0, of Eq. (29 can be done recursively using only the Hessian with respect to the likelihood, i.e.
X Via), the computation
Σ (m) ( Pot:Laplaee)
(m) I:Lapace
V (fm) log p(u ((m). xx) 00m)
(30)
EMAP
In this work, our principal concern (Gaussian) regression, so we do not undertake extensive experiments for classification. Nevertheless, we provide a simple experiment using the Laplace approximation in Appendix A where the proposed method is tested on toy dataset. We show performance comparable other recent online GP approximations, and demonstrate that dynamic models can outperform static models when the distribution of data changes gradually over time.
If classification is of particular concern, several other approximation methods are available that one could experiment with. For example, particle filtering (Djurić et al., 2003) or posterior linearization filters (Tonnrp cct al. 2018) offer other methods for inference in state-space models with non-Gaussian mcasurement models. We believe the adaptation of DOEBE to these non-Gaussian settings is an interesting avenue for future work, and speculate that they may provide better results than the Laplace approximation.
3.5 Theoretical Aspects
In Lu ct al. 20222. several online regret bounds were derived. The first regret bound is similar to that of Kakade & Ng (2004) and bounds the cunulative loss between the IE-GP and any individual expert of the ensemble. This result does not directly depend on the use of RFF GPs in particular and is easily adapted to the OEBE, with minor additional assumption on each ((m))
For this, we will define the loss function (xf(x,) yr) as the negative log- ligeeiill of data yr using estimator (x,). We will define the loss of expert at time + ((m) as the negative predictive log- lideliodd of 3441 under the posterior predictive of estimator m at time t. Likewise, 4411 is the negative predictive log- likelihood of yt+1 under the posterior predictive of the OEBE estimator at time t. We are then able to bound the gap between the cumulative loss of any expert m with fixed parameters !!m) and the OEBE estimator.
Theorem (Online Regret Bound for OEBE). Let the negative log likelikeddoood L(:iy) be C2 with bounded d8
second derivative, i.e., dz, C(zrizz)) & for all and some constant Let us consider an OEBE with prior
(((())) = N (glm), 00 0.00 (m))
(m)² IFM)
where F... dim (()) Purthermore, assume that for any x, the norm of (fim) is bounded by unity.
10

---

Published in Transactions on Machine Learning Research (05/2024)
Then the curallatve negotive log-likelihood incurred by the OEBE is bounded with respect to any single expert m with fized parameters m by
T
T Σ
| Fm
T TCT glm) Fm
Σ
4
(x,)) T _o ;2r) 35)) 20
log
1+
+ log M.
(31)
2
Proof. The proof follows immediately from the proof of Lemma 1 in Lu et al. (2022, Sec. 9.1), which makes no reference to RFF GPs besides thoir being a (bounded) basis expansion. For notational consistency and completeness, we provide a proof in Appendix B
The authors of Lu et al. 2022) then show regret bounds with respect to any learner in the RKHS of any expert. by essentially using the uniform convergence of RFF GPs (Rahimi & Recht 2007). A higher probability bound could also be derived with the classical results of Sutherland & Schneider (2015) Similar uniform convergence results exist for HSGPs (Solin & Särkkä 2020), which we posit can be used to derive similar regret bounds with respect to the RKHS of each expert. However, we do not discuss them here, our main focus is the ensembling of several different types of basis expansions.
The switching regret bounds in Lemma 2 and Lemma 3 of Lu et al. (022) also do not depend on the RFF approximation and can be adapted to SOEBE. The resulting bounds require more bookkeeping on account of the possibly varying dimension of each ((m), but the proofs again remain essentially the same as for IE-GPs. We also omit derivation of these bounds to not distract from our main interest, which more practical in nature.
4 Choosing Models for Ensembling
In this section, we make general comments on how an ensemble can be constructed. This includes how to sample kernel hyperparameters (Section 4.1), how handle problems where it is unelear a change in the prediction problem will occur (Section 4.2 and why one may want include generalized additive models in the ensemble (Section 4.3
4.1 Sampling Hyperparameters from the Marginal Likelihood
For IE-GPs, some initial data Do is set aside to determine og and of via empirical Bayes (Section 2.1.2 There is no guidance given on how obtain kernel hyperparameters. To allow for expressive basis functions allow all hyperparameters to be optimized with respect to the depending on many hyperparameters, marginal likelihood.
Clearly, there is no diversity among the same type of model in the ensemble if we simply take the estimates via empirical Bayes. We provide a two-pronged approach to promote ensemble diversity, targeting both the multimodality and probabilistic nature of the marginal likelihood.
The first strategy hinges on the observation that marginal likelihoods are often multimodal (Rasmussen & Williams 2005, pp. 115). To leverage this fact, optimize with several initial sets of hyperparameters. Tn practice, find that initializing models with length scales to la s(maxp, Iz mammi xa) for E 1, works well 5
Second, we note that the marginal likelihood is indeed a valid probability distribution and propose taking a Laplace approximation of the marginal likelihood at cach local maximum. In particular, we sample M- sets of hyperparameters from the Laplace approximation, in addition to keeping the empirical Bayes estimate, for a total of M estimates at each mode the variance terms whilst allowing for legitimate ensemble. Additionally, if the Hossian can be computed quickly (e.g., using automatic differentiation when the number of hyperparameters is small), this sampling
Fl This maintains the same spirit of using empirical Bayes for
TThe 1 case used initialization the code of Lizaro-Credilla al. 2010 4 Foor optimization, each constrained hyperparameter transtormed an unconstrained one. The sampling occurs in the unconstrained space.
11

---

Published in Transactions on Machine Learning Research (05/2024)
step is also quick to compute. We observe that this provides better accuracy, significant at the p < 0.05 level according to a Wilcoxon signed-rank test (Appendix C
It is tempting to retrain the hyperparameters periodically as an additional method of dealing with a changing task. While this may be worthwhile in some applications, it has the unfortunate consequence of changing the computational complexity of estimation over time. particular, the ensemble is retrained periodically, the online algorithm will drastically increase its complexity periodically, which is often undesirable in an online algorithm. Continual learning of hyperparameters, for example with online stochastic gradient descent also adds unfortunate complexities by changing the measurement model of the state-space without changing the current posterior. In our experiments, we therefore fix all hyperparancters after the initialization above is performed.
4.2 Ensembles of Ensembles for Dynamic and Static Models
Dynamic models are expected to perform significantly better than static models if a progressive change in the online estimation problem occurs. However, if the data consists of i.i.d. draws from a GP, it clear that dynamic models will overestimate the predictive variance (see Eq. 26 Therefore, it is desirable to be able to ensemble static and dynamic models.
(m²) is zero for some models and nonzero for A straightforward solution is to use a DOEBE model where orw others. However, recall from Section 2.4 that a central property of BMA is that it collapses in the limit of infinite data to the model "closest" to the generating model. In practice, this collapse happens rapidly, often within a few thousand samples in our experiments, and weights may become exactly zero due to numerical underflow or implementation. This highly problematic, this collapse can (numerically) occur before it becomes clear that dynamic models are necessary.
To theoretically avoid the collapse of BMA weights, any matrix Q which describes an ergodic Markov chain is sufficient | However, the collapse of BMA weights to zero is numerical property rather than a theoretical property in the first place. Moreover, we can incorporate the intuition that the methods that perform best statically are likely to also perform best dynamically, and vice versa.
We propose an SDOEBE with particular form of Q which incorporates both these points. In particular, we specify an ensemble of 2M models, the first M dynamic and the second M static, and we build an SDOEBE matrix Q with the following block structure:
(1 - 8
81M
Q
SIM
[M
Because of the block structure of Q we call this model Ensemble of DOEBEs (E-DOEBEs). It can be interpreted as switching between dynamic and static models, where the probability of switching is given by 8.
The E-DOEBE unfortunately introduces the hyperparameter 8, which needs to be selected, in addition to the (m) for cach m. These hyperparameters are not casily amendable to the empirical Bayes approach
values of a(m) either, since their entire purpose is to improve performance on timescales larger than the pretraining period.
In the case of ó, we performed an experiment in Appendix D where the value of ő is varied from 0.01 to 0.2; we found that the effects are relatively small, and a value Meanwhile, for of the block structure of the E-DOEBE can be repeated to ensemble several values of (lt) oh For example, to ensemble R different models with of E { orw }, we would use the RM x RM matrix with (1 - -R6) as on the block-diagonals and 8IM elsewhere.
0.05 works well on all tested datasets.
While the E-DOEBE a special case of the SDOEBE, is novel formulation and consistently outperforms any given DOEBE model in our experiments.
Thiis follows immediately from the definition of ergodicity and that the Gaussian distribution's support of B.
12

---

Published in Transactions on Machine Learning Research (05/2024)
4.3 Why Additive Models May be Preferable
In the case of either RFF GPs HSGPs, performance degrades rapidly for fixed feature dimension F and increasing input dimension D. However, GPs with "additive" structure have been widely explored (Duvenaud ct al. 2011). The simplest case is when the GP is given the structure of a generalized additive model (GAM) (astte & Tibshirani 1990p, known as a GAM-GP:
f(x) ~ SP(1. 100~sp((p
C.. Exaltassta) = f(x) = fa(ea). E..Enn d=1
D
(32)
d=1
where fa(aad) GP (O,,) GPs with respective kernel function na(;,))
In this case, the number of basis functions in the HSGP approximation grows linearly with D, and likewise, so does the volume integrated over in the RFF approximation. This can be seen from the hand side of Eq. (32). as we can simply take D independent approximations. Therefore, while the GAM-GP structure may be restrictive, we can obtain better-quality approximation as the input dimension of the data D grows for fixed number of features F. This is useful because, in medium-to-hige dimensions the approximation error stemming from high-dimensional quadrature may dominate the approximation error from assuming additive structure.
For large datasets, Hensman et al. 2018) show how a variational version of RFF GPs with additive structure can be effective and efficient, and Solin & Särkkä (2020) show the same for HSGPs. [Debrrigee et al. (2000) explore GAM-GPs and a generalization using random projections and note their "surprisingly" good performance on large datasets and in high dimensions.
To test the effectiveness of additive HSGPs compared RFF GPs, we experimented with the "large" datasets used in Delbridge et al. (2020) and found that DOEBEs consisting of additive HSGPs may often outperform DOEBEs consisting of RFF GPs more details are available in Appendix F
5 Experiments
We provide three different experiments in the main text, with additional experiments in the appendices. In the first experiment (Section 5.1], we compare ensembles of several different basis expansions, showing that the best-performing model varies widely. In the second experiment (Section 5.2), we show how model collapse between static and dynamic models occur, and how the model introduced in Section 4.2 alleviates it. Finally, show that E-DOEBE can effectively ensemble methods that are static and dynamic, and of different basis expansions (Section 5.3).
The metrics we use are the normalized mean square error (nMSE), and the predictive log-likelihood (PLL). The nMSE is defined to be the MSE of ye with the predictive mean, divided by the variance of yur. In particular, at time t, the nMSE
E (My- yr) MMEE 6. Var(y1:))
The predictive log-likclihood (PLL), which is the average value of 1 ie., PLL, Σ/=1 1
Across all experiments, use several publicly available datasets, ranging in both size and the number of features. A table of dataset statistics is available in Table [ Friedman # and Friedman #2 are synthetic datasets designed to be highly nonlinear, and notably, are i.idd The Elevators dataset is related to controlling the elevators on an aircraft. The SARCOS dataset uses simulations of a robotic arm, and Kuka #1 is a similar real dataset from physical experiments. CaData comprised of California housing data, and the task of CPU Small to predict a type of CPU usage from system properties.
All hyperparameter optimization was done on the first 1,000 samples of each dataset; since already assume access, each dataset was additionally standardized in both and y using the statistics of the first 1,000 samples. We follow Lu et al. 2022) in setting a weight to 0 when it falls below the threshold of 10 10 16
Tn this context, by "high-dimensional" " we mean D 20. This similar its usage in, e.g. Binois & Wycolf 2022
13

---

Published in Transactions on Machine Learning Research (05/2024)
Datasct Name
Number of Samples
Dimonsionality
Friedman #1
Friedman 1991
40,000
10
Friedman #2
Friedman 1991
40,000
Elevators (Torgo)
16,599
17
SARCOS [Rasmussen Wil iiams
2005
44,484
21
Kuka #
MMiee et al. 2014
197,920
21
CaData
Pace & Barry 1997
20,640
CPU Small
Delve Developers)
1996
8,192
12
Table 1: Dataset statistics, including the number of samples, the number of features, and the original source. In addition to the original sources above, several of these datasets were curated by the UCT Machine Learning Repository (Kelly et al.) or LibSVM (Chang & Lin
2011.
5.1 Comparing Different Basis Expansions
To show that having a wide varicty of basis expansion models available is useful, we test several model types on each dataset of Table D Additionally, we lest static and dynamic versions of models to compare performance.
Models used for comparison include additive HSGP model (D)OE-HSGP), an RFF GP (())EE-FFF), ensemble of quadratic, cubic, and quartic polynomials with additive structure (D))EE-Pol]], linear regression (D)OE-Linear), and a one-layer RBF network (D))E-RBF). Apart from additional hyperparameter tuning in an ARD kernel, the (D)OE-RFF model is identical to the (D)IE-GP of Lu et al. 20022)
For RFF GPs, 50 Fourier features were used (so = x 50), and for HSGPs, [100/ DJ features were used for each dimension (so F < 100). In either case, an SE-ARD kernel was used. For RBF networks, 100 locations were initialized via K-means and then optimized with empirical Bayes, well ARD length scales. For all models except RBF networks, ensembles were created by the process proposed in Section 4.1 for RBF networks, the computation of the Hessian was too expensive, and parameters were randomly perturbed by white Gaussian noise with variance 10-3 instead.
For dynamic models, of was set to 10-3 following Lu et al. 2022 The initial values of oj and of were 1.0 and 0.25, respectively. Optimization was performed with Adam (Kingma & Ba 2014).
Results Results of the average nMSE and PLL can be found in Fig. 2 nMSE and PLL a function of samples seen in Appendix F We find that the best-performing class of models differs dramatically from dataset to dataset. In particular, in terms of both nMSE and PLL, HSGPs, RFF GPs, and RBF networks all attain the best performance on at least one dataset. This reinforces the idea that ensembling with respect to several different models is beneficial, as there not a single method consistently outperforming the others.
We provide full plots of (online)
Additionally, expected, dynamic models can dramatically outperform static models in specific cases (for example, on SARCOS and Kuka #1), but attain worse PLL on datasets where the data is reasonably i.i.d. (for example, Friedman #1).
As expected, when an additive structure is a reasonable approximation, additive HSGP methods outperform RFF GPs, for example on Kuka #1 and CaData. The RFF GP approximation rarely performs particularly poorly, making it consistently "good" estimator, and achieves the highest PLL on Friedman #2, SARCOS, and CPU Small. However, it is also at times outperformed by simpler methods, such as the RBF network, highlighting the potential benefits of using diverse basis expansions.
Key Takeaways Key takeaways from this experiment include: (1) neither dynamic nor static methods are strictly superior across all sottings, (2) no single basis expansion superior across all datasets, and (3) RFF GPs consistently provide good performance, but this performance can often be improved upon by using other basis expansions.
14

---

Published in Transactions on Machine Learning Research (05/2024)
DOE-HSGP OE-HSGP
DOE-RFF OE-RFF
DOE-Poly OE-Poly
DOE-Linear OE-Linear
DOE-RBF
OE-RBF
100
10
OMSE
10
10-3
2
PLL
0
-2
Friedman #1
Friedman #2
Elevators
SARCOS
Kuka #1
CaData
CPU Small
Figure 2: Results of Experiment 1. Pictured are the nMSE (lower is better) and PLL (higher is better) with error bars denoting one standard deviation over 10 random trials. The best-performing method on each dataset and metric highlighted with bold edges to preserve readability, the nMSE axis was bounded +1.0 and the PLL axis is bound at +2.5, even if points extend past this.
5.2 The Necessity of Ensembles of Dynamic Ensembles
In this experiment, we show that the E-DOEBE, model proposed in Section 4.2 can indeed prevent premature collapse of BMA weights. While this premature collapse of BMA weights does not seem typical in real datasets, it is not hard to illustrate its possibilily, even on real datasets with high-performing methods.
As a constructive example, we can create an ensemble of additive HSGPs on the Kuka #1 dataset, where dynamic models performed significantly botter in Section 5.1 In particular, we created an ensemble of two additive HSGPss where the first model is dynamic (in particular, drw (1)2 = 10-3) and the second model is
.2)² 0). The ensemble hyperparameters were fit with empirical Bayes, with the initial length scale values being the vector ones. Then, the resulting ensemble was trained online as a DOEBE and as an E-DOEBE, with 5 10-2 Note that in this carefully controlled setting, each basis expansion entirely deterministic given the hyperparameters, the results are purely deterministic and cannot be attributed to poor random seeds.
static (i.e., orw
The resulting weights, pictured in Fig. 3 demonstrate that premature collapse of BMA weights can be a problem. Numerically, the log-likelihood of the E-DOEBE model is dramatically better than that of the DOEBE model (Table 2), showing this collapse can be catastrophic.
This issue can be partially averted by climinating the threshold of 10-16 when ensembling. Indeod, this example, the weights reach minimum of approximately 10-72 However, with any finite precision arithmetic, there is always the potential for this type of collapse to occur due to numerical underflow. It is
15

---

Published in Transactions on Machine Learning Research (05/2024)
trivial to construct such examples by generating the first Ni samples with oew! (m) 0 until weight collapse occurs, and the rest of the dataset with orw >
Key Takeaway The key takeaway this experiment is that an ensemble dynamic and static models can catastrophically collapse- even when the discrepancy performance along the entire dataset is large and that the E-DOEBE approach proposed in Section 4.2 can avoid this collapse.
Weight of Dynamic Model on Kuka #1
Weight of Dynamic Model
1.0
DOEBE E-DOEBE
0.8
0.4
0.2
10"
10a
102
10a
10*
10
Number of Data Points Seen
Figure 3: Results of Experiment 2. Pictured is the 10-sample moving average of the BMA weight corre- sponding the dynamic additive HSGP when trained DOEBE and E-DOEBE.
Method DOEBE E-DOEBE
Predictive Log-
-403.41
0.55
Table 2: Predictive log-likelihood of DOEBE and E-DOEBE models in Experiment 2 (higher is better).
5.3 E-DOEBE Outperforms Other Methods
The ultimate goal of the E-DOEBE model is to combine static and dynamic models of several different types. To do so, we repeat the experiments of Section 5.1 while comparing to an E-DOEBE model. We restrict our allention to static and dynamic versions of the three best-performing families of models in Experiment ((D)OE-HGGP, (D)OE-RFF, and (D)OE-RBF), and an E-DOEBE ensemble containing all them. The E-DOEBE model is created with = 10-2, which was not tunod. Results are presented in Fig. Full (online) nMSE and PLL plots from each experiment can be found in Appendix
As desired, the E-DOEBE model effectively ensemble dynamic and static models of different basis expansions. Across all experiments, the E-DOEBE model performs the best in terms of PLL, and is the best in terms of NMSE for all but one dataset (Friedman #2).
Key Takeaway The E-DOEBE can effectively ensemble several different ensembles of high-performing basis expansions, resulting in consistently better performance than any single method.
5.4 Additional Experiments
We undertook several additional experiments which are omitted from the main text but included in appen- dices. The first such experiment is simple streaming classification task using the Laplace approximation outlined in Section 3.4 (Appendix A). We also perform an experiment validating that the sampling tech- nique of Section 4.1 can be beneficial (Appendix C,, and a small ablation study on the parameter of the F-DOEBE (Appendix D). We also compared the (D)OE-RFF and (D)OE-HSGP methods on additional datasets appearing the GP literature (Appendix E Finally, we provide some comments and experiments regarding a variant of (D)OF-RFF where the frequencies are optimized (Appendix H
16

---

Published in Transactions on Machine Learning Research (05/2024)
DOE-H5GP
OE-HSGP
DOE-RFF
OE-RFF
DOE-RBF
OE-RBF
E-DOEBE
100
10 10-1 1
nMSE
10-2
10-3
2
PLL
0
-2
Friedman #1
Friedman #2
Elevators
SARCOS
Kuka #1
CaData
CPU Small
Figure 4: Results of Experiment 3. Pictured are the nMSE (lower is better) and PLL (higher is better) with error bars denoting one standard deviation over 10 random trials. The best-performing method each dataset and metric is highlighted with bold edges to preserve readability, the nMSE axis was bounded at +1.0 and the PL. axis is bound at +2.5, even if points extend past this.
6 Conclusion & Future Directions
In this paper, showed that recent progress in online prediction using RFF GPs can be expanded to arbitrary linear basis expansions. This included several basis expansions that outperform RFF GPs on real and synthetic datasets. We show how, in simple framework, different linear basis expansions can be ensembled together, increasing ensemble diversity. While several popular choices of basis expansions were used, it would be interesting expand the tests even further, particularly with splines.
We also showed that premature collapsing of BMA weights can be an issue in online ensembling. We introduced the E-DOEBE model, which alleviates this issue, and showed its effectiveness. However, this meta- ensembling may be seen adding complex bandage to BMA rather than addressing the underlying issues. Further research may explore how to incorporate other Bayesian ensembling methods, such as Bayesian (hierarchical) stacking (Yao et al. 2018 2022).
While we provide guidance on how to initialize ensembles given a set of basis expansions, deciding which basis expansions to use an important open topic. A naive approach would be to expand on the existing use of the marginal likelihood for model selection, but this may be "unsafe" when using different basis expansions [Llorente et al. 2023) and therefore requires care. We additionally provided several ideas for inference with ann-Gasssian likelihoods, for example, for classification tasks. Determining which, if any, of these tasks is superior to the Laplace approximation is another interesting topic for future study.
Finally, could beneficial to modify or add new basis expansions in the online setting, Indeed, recent progress GPs has worked towards selecting and adapting kernels online great benefit (Duvenaud al., 2013 Shin et al.) 2023). such techniques could be adapted to DOEBE, it could eliminate the pre-training period, and allow for adapting the domain of approximations when data arrives.
17

---

Published in Transactions on Machine Learning Research (05/2024)
Acknowledgments
This work was supported by the National Science Foundation under Award 2212506.
References
Mickael Binois and Nathan WycolI. A survey on high-dimensional Gaussian process modeling with applica- tion to Bayesian optimization. ACM Transactions on Evolutionary Learning and Optimization, 2(2):1-26, 2022.
Christopher M. Bishop. Pattern Recognition and Machine Learning. Information Science and Statistics. Springer, 2007. ISBN 9780387310732.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Nilne, and Qiao Zhang. JAX: compos- able transformations of Python-NumPy programs, 2018. URL http://github. com www. com/google/jar
Thang Bui, Cuong Nguyen, and Richard Turner. Streaming sparse Gaussian process approximations. Advances Neural Information Processing Systems, 30, 2017.
Chil-Chung Chang and Chih-Jen Lin. LIBSVM: library for support vector machines. ACM Transaclions on Intelligent Systems and Technology (TIST), 2(3):1 27, 2011.
Carl de Boor. On calculating with B-splines. Journal of Approzimation Theory, 6(1):50 62, 1972.
Carl de Boor A Practical Guide to Splines, volume 27 Applied Mathematical Sciences. Springer New York, NY, 1978. ISBN 9780387953663.
Ian Delbridge, David Bindel, and Andrew Gordon Wilson. Randomly projected additive Gaussian processes for regression. In International Conference on Machine Learning, pp. 2453 2463. PMLR, 2020.
Delve Developers. Delve datasets. http::///w.....ttooooo.o... http::/w...rrrtoo..eeeaaaaaaaaaa
html
1996.
T. Diethe. 13 benchmark datasets derived from the UCT, DELVE and STATLOG repositories. https:
//githu.. http://dietheegguaarrrnee
datasets/ 2015.
Petar M Djurić, Jayesh H Kotecha, Jianqui Zhang, Yufei Huang, Tadesse Ghirmai, Mónica F Bugallo, and Joaquin Miguez. Particle filtering. IEEE Signal Processing Magazine, 20(5):19 38, 2003.
David Duvenaud, James Lloyd, Roger Grosse, Joshua Tenenbaum, and Ghahramani Zoubin. Structure discovery in nonparametric regression through compositional kernel search. In International Conference on Machine Learning, pp. 1166-1174. PMLR, 2013.
David K Duvenaud, Hannes Nickisch, and Carl Rasmussen. Additive Gaussian processes. Advances in. Neural Information Processing Systems, 24, 2011.
Yaakov Engel, Shie Mannor, and Ron Meir. The kernel recursive least-squares algorithm. IEEE Transactions on Signal Processing, 52(8):2275-22) 2285, 2004.
Jerome H Friedman. Multivariate adaptive regression splines. The Annols of Statistics, 19(1):1-6) 1991.
Arjan Gijsberts and Giorgio Metta. Real-time model learning using incremental sparse spectrum Gaussian process regression. Neural Networks, 41:59 69, 2013.
T.J. Hastic and R.J. Tibshirani. Generalized Additive Models. Chapman & Hall/ Monographs on Statistics & Applied Probability. Taylor Francis, 1990. ISBN 9780412343902.
E. Hazan. Introduction to Online Conver Optimization. Adaptive Computation and Machine Learning Series. MIT Press, second edition, 2022. ISBN 9780262046985..
18

---

Published in Transactions on Machine Learning Research (05/2024)
James Hensman, Nicolas Durrande, and Arno Solin. Variational Fourier features for Gaussian processes. The Journal of Machine Learning Research, 18(151):1-52, 2018.
Jennifer Hoeting, David Madigan, Adrian Raftery, and Chris Volinsky. Bayesian model averaging: tutorial (with comments by M. Clyde, David Draper and EI George, and rejoinder by the authors). Statistical Science, 14(4):382-411 1999.
Sham M Kakade and Andrew Ng. Online bounds for Bayesian algorithms. Advances in Neural Information Processing Systems, 17, 2004.
Robert E Kass, Luke Tierney, and Joseph B Kadane. Laplace's method in Bayesian analysis. Contemporory Mathemalics, 115:89 99, 1991.
Steven M. Kay. Fandamentals of Statistioul Signal Processing: Estimalion Theory. Fundamentals of Statis- tical Signal Processing, Prentice- Hall, 1993.
Markelle Kelly, Rachel Longjohn, and Kolby Nottingham. The UCI machine learning repository. https:
//archive. ics. w.c.suci.com
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. ar preprint
arXiv:1412 6980, 2014.
Jyrki Kivinen, Alexander J Smola, and Robert C Williamson. Online learning with kernels. IEEE Transac- lions on Signal Processing, 52(8))2165 2176, 2004.
Miguel Lázaro-Gredillla Joaquin Quinonero-Candela, Carl Edward Rasmussen, and Aníbal R Figueiras- Vidal. Sparse spectrum Gaussian process regression. The Journal of Machine Learning Research, 11: 1865-1881, 2010.
Yuhao Liu, Marzich Ajirak, and Petar M Djurié. Inference with deep Gaussian process state space models. In 2022 30th European Signal Processing Conference (EUSIPCO), pp. 792 796. IEEE, 2022.
Yuhao Liu, Marzieh Ajirak, and Petar M Djurić. Sequential estimation of Gaussian process-based deep stale-space models. IEEE Transactions on Signal Processing, 71:2968 2980, 2023a.
Yuhao Liu, Chen Cui, Daniel Waxman, Kurt Butler, and Petar M. Djurić. Detecting confounders in mul- tivariate time series using strength of causation. In 2023 31st Europeon Signal Processing Conference (EUSIPCO), pp. 1400-1404, 2023b.
Fernando Llorente, Luca Martino, David Delgado, and Javier Topez-Santiago. Marginal likelihood compu- tation for model selection and hypothesis testing: extensive review. SIAM Review, 65(1):3 58, 2023.
Jing Lu, Steven CH Hoi, Jialei Wang, Peilin Zhao, and Zhi- Yong Liu, Large scale online kernel learning. The Journal of Machine Learning Research, 17(47):1, 2016.
Qin Lu, Georgios Karanikolas, Yanning Shen, and Georgios B Giannakis. Ensemble Gaussian processes with spectral features for online interactive learning with scalability, In International Conference on Artificial Intelligence and Statistics, pp. 1910-1920. PMLR, 2020.
Qin Lu, Georgios V Karanikolas, and Georgios B Giannakis. Incremental ensemble Gaussian processes. IEEE Transactions on Pattern Analysis and Machine Intelligence, 45(2):1876 1893, 2022.
Qin Lu, Konstantinos D. Polyzos, Bingcong Li, and Georgios B. Giannakis. Surrogate modeling for Bayesian optimization beyond single Gaussian process. IEEE Tronsactions on Pattern Analysis and Machine Intelligence, 45(9):11283 -12966, 2023.
Franziska Meier, Philipp Hennig, and Stefan Schaal. Incremental local Gaussian regression. Advonces in Newral Information Processing Systems, 27, 2014.
Objax Developers. Objax, 2020. URL https: httpss///ggthbbb.ccoggglgl//////
19

---

Published in Transactions on Machine Learning Research (05/2024)
Francesco Orabona. A modern introduction to online learning. arXiv preprint arXit:1912.1991: 2019.
R. Kelley Pace and Ronald Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33(3): 291-297, 1997.
Konstantinos D Polyzos, Qin Lu, and Georgios Giannakis. Ensemble Gaussian processes for online learning over graphs with adaptivity and scalability. IEEE Transactions on Signal Processing, 70:17-30, 2021.
Joaquin Cainonero-Candela and Carl Edward Rasmussen. A unifying view of sparse approximate Gaussian process regression. The Journal of Machine Learning Research, 6:1939-1999 1959, 2005.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in Newral Information Processing Systems, 20, 2007.
Carl Edward Rasmussen and Christopher K. Williams. Gaussion Processes for Machine Learning. MIT Press, 2005.
Gunnar Rätsch, Takashi Onoda, and K-R Müller. Soft margins for AdaBoost. Machine Learning, 42:287-322 320, 2001.
Gabriel Riutort-Mayol, Paul-Christian Bürkner, Michael R Andersen, Arno Solin, and Aki Vehtari. Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. Statistics and Computing, 33(1):17, 2023.
Simo Särkkä and Lennart Svensson. Boyesion Fillering and Smoolhing, volune 17 of Institute of Mathemal- ical Statistics Tertbooks. Cambridge University Press, second edition, 2023. ISBN 9781108912303.
Eura Shin, Predrag Klasnja, Susan Murphy, and Finale Doshi-Velez. Online model selection by learning how compositional kernels evolve. Transactions Machine Learning Research, 2023. URL https:
//openreview.nnt//oooueee wtt.:nteooooommmmmdddmaappa
Arno Solin and Simo Särkkä. Hilbert space methods for reduced-rank Gaussian process regression. Statistics
and Computing, 30(2):419-446, 2020.
Samuel Stanton, Wesley Maddox, lan Delbridge, and Andrew Gordon Wilson. Kernel interpolation for scalable online Gaussian processes, In International Conference on Artificial Intelligence and Statistics, pp. 3133 3141. PMLR, 2021.
Danica J Sutherland and Jeff Schneider. On the error of random Fourier features. In Proceedings of the Thirty-First Conference on Uncertainty in Artificial Intelligence, pp. 362-871. 2015.
Sergios Theodoridis. Machine Learning: A Bayesian and Optimization Perspective. Academic Press, second edition, 2020. ISBN 9780128188033.
Luš Torgo Regression datasets. URL Last Accessed March 19, 2024.
Filip Tronarp, Ángel Garcia- and Simo Särkkä. Iterative filtering and smoothing nonlincar and non-Gaussian systems using conditional moments. IEEE Signal Processing Letters, 25(3):408 412, 2018.
Andrew Wilson and Ryan Adams. Gaussian process kernels for pattern discovery and extrapolation. International Conference on Machine Learning, pp. 1067-1075. PMLR, 2013.
Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman. Using stacking to average Bayesian predictive distributions (with discussion). Bayesian Analysis, 13(3):917 11007, 2018.
Yuling Yao, Gregor Piră, Aki Vehtari, and Andrew Gelman. Bayesian hierarchical stacking: Some models are (somewhere) useful. Bayesian Analysis, 1744):1043-1077 2022.
A. Zellner. An Introduction to Bagesian Inference in Fconometrics. John Wiley & Sons, 1996.
20

---

Published in Transactions on Machine Learning Research (05/2024)
A
Classification Experiment
In the main body of the text, we focused on the empirical performance of DEEEE on Gaussian regression problems. While there are several possible improvements to classification (or other non-Gaussian) problems, as outlined in Section 3.4 we perform a simple experiment on the common toy dataset "banana" " (Rätsch
let al. 2001
The banana dataset is notably not exchangeable, with data being loaded from left to right. Based on model knowledge and the experimental results in Section 5.1 we would expect dynamic models to perform better on the non-exchangeable dataset, and static models to perform better on a shuffled (i.e., ii.d.) version. We first experiment with the shuffled version of the dataset, and then with the unshuflled version. Models in each experiment are identical to those of Section 5.3 except using logistic likelihood via the Laplace approximation.
A.1 Shuffled Data
By shuffling the data, we expect to be i.i.d., in which case we expect static models to perform well and the classification error to be steadily decreasing. We show the classification error as function of the number of data seen and an example of the predictive distribution for the OE-RFF, both pictured in Fig. 5
The data is obviously non-additive, so the HSGPs do not perform particularly well. The dataset is quite
simple and, consistent with our regression experiments on simple i.i.d. datasets, the RBF networks perform quite well. Notably, the E-DOEBE consistently performs as well or better than the best-performing model. Finally, we note that the results with RFF models are comparable to those reported in Lu et al. (0022) for IE-GPs.
DE GP Prosictive (Shuffied)
Bariana (Shuffled)
Ener
2-
-0.6
0..4
Number Data
x10
DCE rog
10+ bOCAEr
0C45522
ECCCES
DCE te+
1.0
Figure 5: (Left) Classification errors for the shuffled banana dataset. Lines indicate the mean over 5 trials, with shaded regions denoting +1 standard deviation. (Right) The predictive distribution of an OE-RFF ensemble after training on the shuffled banana dataset.
A.2 Unshuffled Data
We also experimented with the case where the data is not shuffled. In this case, the data is presented from left-to-rihh,, as visualized in Fig. The resulting classification errors are pictured in Fig. 7 and are comparable several other works performing streaming GP classification (Stanton et al., 2021 Bui el al. 2017), where classification errors are or around 10%.
Again consistent with results for regression tasks, the dynamic models perform better in this non- exchangeable setting. Intuitively, the RBF models perform quite poorly as the locations μ4 are only trained
7Tee original link to the dataset longer available, but copy available online at Diethe 2015
21

---

Published in Transactions on Machine Learning Research (05/2024)
the first 1,000 points, which are not representative of the entire dataset. The RFF models still perform well, and interestingly, the DOF-HSGP model performs competitively well. This may be because "local" additivity is more realistic condition.
2
10
10
10
10
Figure 6: The predictive distribution of an OF-RFF ensemble after observing the first. 25%, 50%, 75%, and 100% of the data, respectively.
Danana (Not Shuffledi
DCE-HSGP D64-500
OC-FF DOE-BI
CERBD F-00Fa6
DCe RT?
Figure 7: Classification errors for the unshuffled banana dataset. Lines indicate the mean over 10 trials, with shaded regions denoting #1 standard deviation.
22

---

Published in Transactions on Machine Learning Research (05/2024)
B Proof of Theorem 1
In this appendix, we provide the proof for the online regret analysis bound given in Theorem | The proof here is essentially the same as the proof in Lu et al. (2022) Section 9) (and in turn, the proof of Theorem 2.2 Kakade Ng (2004)). ) It. provided for notational consistency and completeness.
Proof of Theorem 1 The proof consists of three steps: (1) bounding the regret of the OEBE to any given expert, (2) bounding the regret of any expert with respect to variational distribution, and (3) bounding the regret of the variational distribution with respect to fixed strategy.
Step 1 For step (1), note that by definition of the BMA weights, we have that
(m)
exp(-(+)
W1 mm)
exp(-1,
Since the initial weights are we (m)
= 1/M, this implies
(m))
woforr
W
Σ 4
(m)
exp
(m)
M Mwi (m)
T
Together, these bound the regret of the OEBE with respect to expert m,
T
T
Σ
4
Σ
(m) = log M + log wi [m) < log M,
(33)
T=1
T=1
since log mm() by construction.
Step For stop (2), we introduce a variational distribution (0(m)) = (m((m))
with varia-
tional parameter g². We then bound the regret of expert m with respect to the expected loss under g(-),
Σ (mm) E.([Coio)
where Lgrm) (() .x.T)- By Bayes's rule, we get that ET=1 ((m) marginal probability of yıat given x1a under model m), and therefore
log.(vi:r( (i.e., the
T
Eg( / ((((m) log p(y1.r(0(m) XL.T d0(m)
(m)
p(ynT) X1.T)
since logp(yir(xx:r) does not depend on O. By the definition of conditional densities, we rewrite p (((m) x1T) p(0(m)) (1.T,x1T))
p())T(xIIT)
p(0(m)
We will then multiply by 1 = (m)) (0((0l)))g())))))))))))))))))))) (m(()()) and rearrange, to get
Σ (m)) E. = / (m)) ((0(m)) dd(m) / (01m)) ((0(m) BI.r, UI.T,X..T)
(0(m))
d0 (m).
p(((m))
Recognizing each term as Kullback Liebler (KL) divergence, this is clearly bounded by DEL (g(0(m)) p(O(m))) by non-negativity of the KL divergence, resulting in the bound
(m) [o DKL. (p(O(m)) || p(((()))
112 + F2
F. (1 2oo E
Fm log (m))
(34)
20 (m)
23

---

Published in Transactions on Machine Learning Research (05/2024)
Step 3 Finally, for (3), we will bound the difference
E mpm)
By taking Taylor series expansion of L(,, Vr) centered at 2t,= = (,()((2)) (0. we find
d
(2t 2T,a)
EylCorm)] o Eg
((g(2))iV+)) d2?
2
for 2, = ((())())))) (x,) T ((m) and some function g(-). By invoking our assumption about boundedness of the second derivative of ((:;y+), we see
d ((g(z));))) d=
(2r 2T,+)
(2r 27,x
2
Eg
S En
c
2
2
and by invoking our assumption that (1es(m)) (x)]12 S 1, see
(2r Zr,0) 121
[gg Lo/m: E.g
(35)
With each step complete, must combine them to finish proof. By combining Eqs. we find
33 to 35
Teg² 2
|0 (m) 12
ng²
Fm (1 2 2og E)
(( T :.r) S
+ Fm log og +
(36)
20 (m) 2
0
After minimizing with respect to the variational parameter ç² and simplifying, the theorem is proven.
24

---

Published in Transactions on Machine Learning Research (05/2024)
C Experiments with Sampling Hyperparameters for Ensembles
To show the benefits of using a Laplace approximation to sample hyperparameters, we repeat the experiments of Section 5.I on the Elevators, SARCOS, and CaData datasets using just DOE-HSGP and DOE-RFF models, creating two copies: one using the type-IT MLE hyperparameters, and the other using the Taplace approximation of the marginal likelihood. The experiment consisted of 100 trials where 10 samples were drawn from each Laplace approximation. For RFF GPs, the frequencies of the Monte Carlo approximation were fixed to be the same for each trial. The method with better performance one-sided Wilcoxon signed-rank test in the relevant direction. Results be found in Table 3
determined via the
Predictive Log-Likelihood
Normalized Mean Square Error
Method
Elevators
SARCOS
CaData
Elevators
SARCOS
CaData
DOE-HSGF-MLE
-0.753 = 0.000
0.421 = 0.000
0.081 0,000
0.221 + 0,000
0.017 = 0.000
0,055 + 0,000
DOE- NOE-HSGP-Sampp
0.778 + 0.003
0.466
0.010
0.120 0.010
0.219 + 0.001
0.018 + 0.000
0.052 0.001
DOE-RFF-MLE
-0.010 0.007
0.756
0.213 0.009
0.178 + 0.003
= 0.001
0.0400 0.000
DOF.RFF-Sample
0.639 + 0.007
0.766 + 0.019
0.247 0.009
0.177 + 0.004
0.018 + 0.001
0.040 +: 0.002
Table Predictive likelihood (higher better) and normalized MSE (lower better) of type-II MLE initialization, plus/minus minus one standard deviation over 100 trials. Bolded entries denote superior performance significant the p 0.05 level according to one-sided Wilcoxon rank-sum test.
and
Overall, the benefit of sampling from the posterior on these datasets appears be small but significant. In particular, in terms of predictive log-likelihood (which DOEBE targets), the Laplace approximation initialization never performs worse.
25

---

Published in Transactions on Machine Learning Research (05/2024)
D Sensitivity of E-DOEBE With Respect to 5
In this section, we test the performance of the E-DOEBE model used in Section 5.3 across different values of d E [0.01, 0.05, 0.1, 0.2]. We ran the E-DOEBE model on each datasel once, and recorded the PLL. Results can be found in Fig. S We conclude that a value of ő = 0.05 typically performs well, but performance is reasonably similar for all 0.01 < & 0.2.
= 0.01
6= 0.05
6= 0.1
6 0.2
10-1
10-3
2
PLL
0
-2
Friedman #1
Friedman #2
Elevators
SARCOS
Kuka #1
CaData
CPU Small
Figure 8: Results of the 8 ablation experiment. Pictured are the nMSE (lower is better) and PLL (higher is better) over one trial. The best-performing mothod on cach datasct and metric is highlighted with bold edges.
26

---

Published in Transactions on Machine Learning Research (05/2024)
E Further Experiments on Large Datasets
To support our assertion that additive HSGPs may be superior to RFF GPs in medium-high dimension, we
perform experiments on the UCI datasets (Kelly et al.) used in Delbridge et al. (2020) A summary of each
dataset used can be found in Table 4
Dataset Name
Number of Samples
Dimensionality d
autos servo machine yacht autompg
159
25
167
209
7
308
6
392
7
housing
506
13
stock
536
11
energy
768
concrete
1,030
airfoil
1,503
5
gas
2,565
128
skillcraft
3,338
19
sml
4,137
26
pol
15,000
26
bike
17,379
17
kin40k
40,000
Table 4: Dataset statistics, including the number of samples and the number features for datasets used Delbridge al. (2020) All datasets are available on the UCI Machine Learning Repository [Kelly et al.).
the HD)O and (D)OF-RFF models of Section 5.1 with the following modifications: the number of Fourier features was set to 250 (so F = 250), and [500/D] features were used for each dimension (so F x 500). Because the datasets are often much shorter, empirical Bayes estimates are taken over the first 100 points only. Dimensions of the input vectors which were identically zero for all samples were removed. Results can be found in Fig. g
In particular,
Overall, we find that HDJOE-HSGPs outperform (D)OE-RFFs on several lasks, with varying input dimension D and length N. This includes several datasets which would typically be considered high- eigh-dimensional in the GP literature, for example, the gas dataset (D 128) and the sml datasct (D 26).
27

---

Published in Transactions on Machine Learning Research (05/2024)
DOE-HSGP
DSGH
DOE-RFF
OE-RFF
100
10-1
nMSE
10-2
10-3
2
PLL
0
-2
autos
servo
machine
yacht
autompg
housing
stock
energy
concrete
airfoil
gas
skillcraft
sml
po/
elevators bike
kin40k
Figure 9: Results of additional experiments on UCT datasets. Pictured are the nMSE (lower is better) and PLL (higher is better) with error bars denoting one standard deviation over 10 random trials. The best- performing method on each dataset and metric is highlighted with bold edges to preserve readability, the nMSE axis was bounded at +1.0 and the PLL axis is bound at +2.5, even if points extend past this.
28

---

Published in Transactions on Machine Learning Research (05/2024)
F Results by Dataset for Experiment 1
Here, we provide the full nMSE and PLL plots from the experiment in Section 5.1
Friedman
Elevators
0.35
0.175
0.150
0.25
0.125
0.25
0.20
0.35
0.15
0.10
0.10
0.025
0.06
0.00
0.000
0.00
0.00
0.25
0.75
1.00
Number of Data Points Seen
x10
Number of Data Points Seen
X10*
Number Data Points Seen
x10*
SARCOS
Kuica #1
Data
0.6
0.4
0.5
0.3
0.6
MASE
OMSE
0.1
0.0
0.0
0.0
Mumber of Data Points Sean
x10*
Number of Data Points Seon
Number of Data Points Sean x10*
CPU Small
0.9
DOE-HSGP OE-HSGP
DOE-Poly OE-Poly
OE-Linear
DOE-RBF
DOE-RFF
DOE-Linear
OE-RBF
OE-RFF
0.2
Number of Data Paints Seen
Figure 10: nMSE plots of each method tested in the experiment of Section 5.1 Lines indicate the mean over 10 trials, with shaded regions denoting =1 standard deviation.
29

---

Published in Transactions on Machine Learning Research (05/2024)
Friedman
Friedman #2
Elevators
-0.6
1.5
-0.7
Average PLL
1.0
1.5
Avarage
0.5
0.0
-05
-1.0
-0.5
-1.0
-1.1
Number of Data Points Seen
Number Data Points Seen
Number Data Points Seen x10*
SARCOS
Kulka #1
CaData
0.50
0.75
-10
Averege PLL
-1000
Averace PLL
-20
-1500
1.25
-2000
-1.50
-30
-2000
-1.75
3000
-2.00
-3500
0.0
0.5
1.0
Number of Deta Points Sean
x10
Number of Data Points Seen X10
Number of Data Points Seen x10
CPU Small
0.4
DOE-HSGP
DOE-Poly
OE-Linear
Average
0.2
OE-H5GP
OE-Poly
DOE-RBF
DOE-RFF
DOE-Linear
OE-RBF
-0.4
OE-RFF
-0.8
1.0
Number
x100
Figure 11: PLL plots of each mothod tested in the experiment of Scction 5.1 Lines indicate the mean over 10 trials, with shaded regions denoting =1 standard deviation.
30

---

Published in Transactions on Machine Learning Research (05/2024)
G Results by Dataset for Experiment 3
Here, we provide the full nMSE plots from the experiment in Section 5.3
Friedman
Friedman
Elevators
0.175
0.25
0.12
0.150
0.30
0.125
0.08
0.30
0.060
0.02
0.025
0.00
0.000
0.00
0.00
0.25
0.75
1.00
125
1.50
Number of Data Points Seon
x10
Number of Data Points Seen
X10*
Number Data Points Seen
x10*
SARCOS
Kuka #1
Ca
1.0
0.40
0.5
0.35
0.8
0.30
0.25
MMSE
0.3
0.20
0.35
0.30
0.1
0.05
0.00
0.0
0.0
Number of Data Points Sean
x10
Number of Data Points Seon x10*
Number of Data Points Sean x10*
CPU Small
0.06
0.07
0.06
0.05
DOE-HSGP OE-HSGP DOE-RFF
OE-RFF DOE-RBF
OE-RBF E-DOEBE
0.04
0.08
0.01
0.00
Number Data Points Seen
Figure 12: nMSE plots of each mothod testod in the experiment of Section 5.3 Lines indicate the mcan over 10 trials, with shaded regions denoting =1 standard deviation.
31

---

Published in Transactions on Machine Learning Research (05/2024)
Friedman
Friedman #2
Elevators
-0.5
2.0
2.5
1.5
Average PLL
1.5
Avarage PLL
-0.7
0.5
0.9
0.5
-0.5
Number of Data Points Seen
Number Data Points Seen
Number Data Points Seen x10*
SARCOS
Kulka #1
CaData
-0.75
-10
Average PLL
Averace PLL
1.00
-1500
-2000
-1.50
-30
-2500
-1.75
-2.00
0.0
0.5
1.0
Number of Deta Points Sean
x10
Number of Data Points Seen X10
Number of Data Points Seen x10
CPU Small
0.5
0.4
Average
DOE-H5GP OE-HSGP DOE-RFF
OE-RFF DOE-RBF
OE-RBF E-DOEBE
0.1
0.0
Number
x10
Figure 13: PLL plots of each mothod tested in the experiment of Scction 5.3 Lines indicate the mean over 10 trials, with shaded regions denoting =1 standard deviation.
32

---

Published in Transactions on Machine Learning Research (05/2024)
H Optimizing Fourier Features
Our framework and implementation make it simple to additionally optimize the Fourier features of the RFF model, which we will refer to as the optimized Pourier features (OFF) model. Such an optimization is not novel and was indeed explored the original paper introducing RFF for GPs [Lazaro-Gredila et al. 2010). The idea of optimizing spectral points also forms the basis of the popular spectral mixture kernel (Wilson & Adams 2013), which uses a Gaussian mixture in the spectral domain instead of particle representation.
In some sense, optimizing frequencies defeats the main motivation of RFF GPs; we the results converge to any particular GP pre-specified kernel. Indeed, one can view the resulting model as a one hidden-laylrr neural network, with activation function ΦRFF(-); likewise, we expect the results to inherit both benefits of neural networks (e.g. flcxibility) and drawbacks (o.g. harsh overfitting).
longer claim that
o0, since
longer approximating Eq. (19) for
We additionally note that the resulting model bears significant resemblance to the Gaussian RBF basis expansion model. As will become evident in our results, there is likewise much similarity in their performance: "simple" datascts, optimizing frequencies provides large benefits. On more complex dynamic datasets, it is instead prone to overfitting.
For experiments, we repeat Section 5.1 using the DOE-RFF model, and a version with optimized frequen- cics (DOE-OFF). Hyperparameters for both were sampled with the "Gaussian" " strategy, as the number of hyperparameters too large the optimized version for the Laplace approximation. Results can be found in Fig. 14
33

---

Published in Transactions on Machine Learning Research (05/2024)
DOE-RFF
DOE-OFF
friedman2
elevators
0.05
0.007
0.006
0.005
mAse
0.004
OMSE
0.02
0.002
0.01
0.001
0.06
0.00
0.000
0.00
0.00
0.25
0.50
0.75
1.00
125
1.50
Number of Data Points Seeen
x10*
Number of Data Points Seen
X10*
Number Data Points Seen
x10*
kukal
cadata
SARCOS
0.35
0.04
0.0177
0.0150
0.03
0,0125
ANSR
0.01.00
0.02
0.0075
0.10
0.0050
0.0025
0.05
0.0000
0.00
0.00
1.
2
Number of Data Po intss Seen
K20f
Number of Deta Points Seon
X10*
Number of Data Points Soon x10
cpusmall
10'
10-1
0.06
10-4
0.05
0.04
0.02
0.01
Mlddmnt
Medmonn2
elsvators
SARCOS
koko1
Caidao
qussmall
0.00
Number Data Points Seen
Figure 14: (Continuous Plots) nMSE plots of random and optimized Fourier features. Lines indicate the mean over 10 trials, with shaded regions denoting #1 standard deviation. (Bar Plots) Summaries of the nMSE (lower better) and PLL (higher is better) with error bars denoting one standard deviation over 10 random trials. The best pess-performin method on each dataset and metric is highlighted with bold edges preserve readability, the nMSE axis bounded +1.0 and the PLL axis is bound at +2.5, even if points extend past this.
34


Reviews:
[{"review_id": "DhBTbNt5vX", "reviewer": "Reviewer_AyHk", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. The authors provide a comprehensive discussion of Bayesian inference based on basis expansions. Quite a few different methods in the literature are properly introduced and compared, and the paper is easy to follow.\n2. The method is proposed following clear motivation: avoiding numerical collapse in BMA due to underflow, and Fig. 3 supports why this issue could be problematic.\n3. The method is simple and effective according to the experiments.\n\nWeaknesses & questions:\n1. What is the benefit of adding a Random Walk on the parameters? The dynamic model in the form of Eq. 26 seems to just make the posterior distribution more conservative. Is there any principled justification for this? \n2. How to select $\\sigma_{rw}^{(m)}$, and $\\delta$ in $Q$? Are they hard to tune?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors provide a thorough survey of Bayesian methods based on basis expansions including discussions of popular basis choices, online inference, and the ensemble of models using Bayesian model averaging (BMA), and empirical analyses are conducted to support their claims. They further introduce a method to combine static and dynamic models within the BMA framework. The proposed method is compared with other baselines through experiments.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See the weaknesses & questions above.", "broader_impact_concerns": "I have no concerns about the broader impacts and the ethical implications.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "WA2FtQHTSe", "reviewer": "Reviewer_j5fF", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: \n\n-\tThe ideas and their ideological histories are clearly presented and well explained. I really appreciate the longer length so that each concept can be clearly introduced.\n\n-\tThe writing and summary are nicely done and the statement of contributions at the beginning summarizes the paper nicely.\n\n-\tI also like that the code was provided and in a pretty readable state.\n\nWeaknesses:\n\n-\tOverall, there’s not a huge amount of novelty especially compared to Lu et al, ’22. The new method seems to be enabling just other basis expansions (many of which are severely dimensionality limited). The experimental results don’t currently suggest that these other options are a very good idea as RFFs require less tuning and as the authors state “it rarely performs particularly poorly”.\n\n       o\tFurthermore, this comes at a cost of not being able to handle non-Gaussian likelihoods very well which the IE-GP can handle.\n\n-\tTo distinguish between static and dynamic models, there’s no real experiments (e.g. Bayesian optimization and active learning) which limits the “practical” applicability in some sense (more on this in the suggestions).\n\n\n-\tThe hyperparameters seem to be fixed after the first 1000 data points. Including a periodic retraining or continuous optimization steps ought to help make even static models even more “dynamic”.\n\n(review edited for viewing clarity)", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper uses a hierarchical ensemble of basis expanded linear models for online learning. The weights in the ensemble are given by the previous steps’ predictive log likelihood, while collapse of the ensemble to focus on a single model is prevented by using a structured matrix. Furthermore, noise is added to each model’s weights to ensure that it can update effectively to new data. Overall, these additions make each model a type of linear state space model, while the Bayesian model averaging piece is also a state space model in some sense. \n\nThe main novelty of the work seems to be that it expands Lu et al, ‘22’s efforts from random fourier features to other Bayesian linear models.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Comments / corrections:\n\n-\tFigure 3 / Table 2: It may be better presentation wise to inline both (for the table, use the tabular environment). As currently placed together, it looks like the captions are running on to each other.\n\n-\tFigure 3: show +- 2sd of mean for these trials.\n\n-\tFigure 2: please don’t trim the NLL for sarcos / kuka #1. I personally tend to think that each dataset should be its own small plot for this reason; however, I understand that this is more time consuming. This would enable different scales for each dataset.\n\n-\tFigures – remove “Results of Experiment []” from captioning.\n\n-\tFootnote 5: missing period at end.\n\n-\tFigure 6 – I believe what is marked as “S-DOEBE” should be “E-DOEBE”.\n\n-\tPresentation wise, I have a bit of trouble comparing the different online learning experiments, which all seem to be mixed up with a comparison of different basis expansions (quite unimportant).\n\n      o\tI would suggest that the authors pick one basis expansion they like and instead compare it against an explicitly defined IEGP with a similar enough kernel from Liu et al, ’22 as well as other online learning techniques using either GPs or linear models (e.g. Stanton et al, ’21, Yadav e al, ’21, Zhang et al, ’23). \n\n-\tIt’s not especially clear to me why these ensembles should be better than a single Gaussian process / RFF, especially as the authors use the MMSE estimate with a single Gaussian distribution. So, by linearity, shouldn’t this mean that there’s really only one set of weights \\theta_t? \n\n      o\tI think I understand that there’s differences in the additive setting when there are various different inputs for each model (each dimension). However, if the inputs are the same, shouldn’t there really only be one global set of weights?\n\nExperiments (to me, these are of high importance):\n\n-\tI’d suggest that the authors focus primarily on the additive Hilbert space features and conduct experiments on high dimensional datasets, e.g. the ones from Delbridge et al, ’20.\n\n-\tI’d suggest that the authors try some sort of online classification problems like in some of the other online GP literature. A useful first approximation may be to turn classification into a regression problem via reformulating the likelihood like was done in Stanton et al, ’21. \n\n        o\t Comparison could then be done to the IE-GP on the datasets from Liu et al, ’21. \n\n-\tBasic baselines that would be nice to include for the learning curves are the performance of non-online models trained on various amounts of data.\n\n-\tTo alleviate my concern that the ensembles are not necessarily better than a single model, it would be good to study M = 1, 2, 3, …, 5 (or so) on the kuka dataset. \n\nCode: \n\nI took a quick look at the code and it is nicely documented with good comments. A few small comments here:\n\n-\tL72 of src/DOEBE/sdoebe.py: should ` final_log_w, log_ws` be the same vector?\n\n-\tL118 of src/DOEBE/doebe.py: why is the svd of the pseudo-inverse necessary? Shouldn’t it just be something like \n\n```python\nhess_matrix = jnp.block(hess)\nu, s, _ = jnp.linalg.svd(hess_matrix)\ninv_hess_factor = u @ jnp.diag(jnp.sqrt(1.0 / s))\n```\n\nModulo potentially needing to clip s at zero…\n\n\nReferences:\n\nLiu et al, 2022, Incremental Ensemble Gaussian Processes.\n\nStanton et al, 2021, Kernel Interpolation for Scalable Gaussian Processes.\n\nYadav et al, 2021, Faster Kernel Interpolation for Gaussian Processes.\n\nZhang et al, 2023, Sequential Gaussian Processes for Online Learning of Nonstationary Functions.", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "bzNCeUb5iy", "reviewer": "Reviewer_jUqT", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**strengths:**\n\nSeveral approaches to GPs are describeded (infinite dimensional, RFF, HSGP, ...), as are approaches to hyperparameter fitting (point estimate, sampling from Laplace approx of marignal likelihood) and Bayesian model averaging. The text is of appropriate length and detail.\n\nAll of these techniques are evaluated empirically in a suite of appropriate empirical benchmarks. \n\nThe authors do a good job in acknowledging limitations of their study and potential future work, e.g. \"We introduced the E-DOEBE model, which alleviates this issue, and showed its effectiveness. However, this meta-ensembling may be seen as adding a complex bandage to BMA rather than addressing the underlying issues\".\n\n**weaknesses:**\n\nThe results are entirely empirical. While this is okay, it would be nice if the authors could at least mention some of the existing theoretical results so that they can relate their empirical findings to the theory (see requested changes below).\n\nThere are some mostly minor errors (see requested changes below).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Random (Fourier) features are a widely used method for approximating Gaussian process inference, and have recently been used to ensemble GPs in an online manner. The authors show how such methods can be used and extended in the setting where any basis expansion model is used (e.g. Hilbert space GP). Empirically, they test various features and basis expansions and find that different settings and datasets result in different rankings of performance. Bayesian model averaging is discussed and investigated. Ensembling of dynamic models and static models is also investigated.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- There is an error in equation (15)-(17). See Equation (2.19) of Rasmussen and Williams book. The authors provided the posterior predictive distribution of the noise free prediction f(X*), but they wrote that this was the posterior predictive distribution of the noiseless distribution y*. I believe what they really want to write is the noise free predictive distribution f(X*). This is in contrast with the earlier finite dimensional setting, where predict and correct updates require the predictive distribution with noise for y*.\n- The authors incorrectly state that inference requires inversion of an $N \\times N$ matrix. What is really required is to \\emph{solve} a linear system, which also scales cubically in $N$, but requires less memory than inversion, has a smaller constant, and is more numerically stable. Again in section 2.3.\n- In section 2.3.1, you could mention any of the classical variance estimates of RFF. E.g. \"On the Error of Random Fourier Features\", Sutherland and Schneider. I am not sure if you can also provide variance estimates for the more recent Hilbert Space GP. If these are available, it would be great to compare and contrast, since the claim empirically seems to be that sometimes RFF is better than HSGP but not consistently, at least in your current context.\n- The last part of the Figure 3 caption does not make any sense and contains some weird whitespace.\n\nMinor:\n- First page, BMA (Bayesian Machine Averaging) is used without introducing the acronym.\n- Second paragraph of 2.1.2. Citation for  “empirical Bayes,” “type-II maximum likelihood estimation,” and “evidence maximization” seems oddly recent. Maybe also Rasmussen and Williams?\n- Third paragraph of 2.1.2. Is this an instance of conjugate prior updates? Perhaps that could be worth mentioning.\n- No date on the Torgo dataset in table 1. Is this known?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: As mentioned above, a few points in the presentation could be clarified, hence the minor revision, but otherwise the paper seems to be well-written and technically solid enough to warrant acceptance.
Decision: Accept (poster)