Recurrent Memory for Online Interdomain Gaussian Processes
Authors: Wenlong Chen
Naoki Kiyohara
Harrison Bo Hua Zhu
Jacob Curran-Sebastian
Samir Bhatt
Yingzhen Li
Venue: NeurIPS 2025 poster ()
Keywords: Gaussian Processes
Online Learning
Continual Learning
Sequence Model
HiPPO
Primary area: probabilistic_methods

Abstract: We propose a novel online Gaussian process (GP) model that is capable of capturing long-term memory in sequential data in an online learning setting. Our model, Online HiPPO Sparse Variational Gaussian Process (OHSVGP), leverages the HiPPO (High-order Polynomial Projection Operators) framework, which is popularized in the RNN domain due to its long-range memory modeling capabilities. We interpret the HiPPO time-varying orthogonal projections as inducing variables with time-dependent orthogonal polynomial basis functions, which allows the SVGP inducing points to memorize the process history. We show that the HiPPO framework fits naturally into the interdomain GP framework and demonstrate that the kernel matrices can also be updated online in a recurrence form based on the ODE evolution of HiPPO. We evaluate OHSVGP with online prediction for 1D time series, continual learning in discriminative GP model for data with multidimensional inputs, and deep generative modeling with sparse Gaussian process variational autoencoder, showing that it outperforms existing online GP methods in terms of predictive performance, long-term memory preservation, and computational efficiency.

Recurrent Memory for Online Interdomain Gaussian Processes
Wenlong Chen' 1+, Naoki Kiyohara' 12,,, Harrison Bo Hua Zhu
Jacob Curran-Sebastian' Samir Bhatt Yingzhen Li', Imperial College London "Canon Inc. University of Copenhagen
wenlong chen21@enieeee chen21@imperial. ac. n. a.kkyyohara2@ddmeriaal ac. uk harrison zhuōsund, ku. dk yingzhen. www.mmeriiia.....mmmmmm........m ac.
Abstract
We propose a novel online Gaussian process (GP) model that capable of capturing long-term memory in sequential data in online learning setting. Our model, On- line HiPPO Sparse Variational Gaussian Process (OHSVGP), leverages the HiРPO (High-order Polynomial Projection Operators) framework, which is popularized in the RNN domain due to its long-range memory modeling capabilities. We interpret the HiPPO time-varying orthogonal projections inducing variables with time- dependent orthogonal polynomial basis functions, which allows the SVGP inducing variables to memorize the process history. We show that the HiPPO framework fits naturally into the interdomain GP framework and demonstrate that the kernel matri- ces can also be updated online in recurrence form based on the ODE evolution of НоРРO. We evaluate OHSVGP with online prediction for ID time series, continual learning in discriminative GP model for data with multidimensional inputs, and deep generative modeling with sparse Gaussian process variational autoencoder, showing that it outperforms existing online GP methods terms of predictive performance, long-term memory preservation, and computational efficiency.
Introduction
Gaussian processes (GPs) popular choices for modeling time series due their functional expressiveness and uncertainty quantification abilities [Roberts et al.) 2013 Fortuin et al. 2020]
However, GPs are computationally expensive and memory intensive, with cubic and quadratic complexities, respectively. In online regression settings, such as weather modeling, the number of time steps can be very large, quickly making GPs infeasible. Although variational approximations, such as utilizing sparse inducing points (SGPR [Tissia)] 2009]: SVGP [ensmma et al. 2013 2015a)
and Markovian GPs (Särkkä and Solin 2019 Wilkinson et 2022]. have been proposed to address the computational complexity. it would still be prohibitive to re-fit the GP model from scratch every time new data arrives. BBui et al. [2077] proposed an online sparse variational GP (OSVGP) learning method that sequentially updates the GP posterior distribution only based on the newly arrived data. However, indicated in their paper, their models may not maintain the memory of the previous data, as the inducing points will inevitably shift as new data arrive. This is a major drawback, as their models may not model long-term memory unless using a growing number of inducing points.
In deep learning, as an alternative to Transformers VVsswani et al. 2017]. significant works on state space models (SSMs) have been proposed model long-term memory in sequential data. Originally proposed to instill long-term memory in recurrent neural networks, the HiPPO (High-order Polynomial Projectioo Operators) framework [Gu et al. 2020 provides mathematical foundations for
* "Equal contribution. Source Code: http://wppppgggtttt..../////ssaa
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

compressing continuous-time signals into memory states through orthogonal polynomial projections. HiPPO is computationally efficient and exhibits strong performance in long-range memory tasks, and forms the basis for the state-of-thee SSMs, e.g., structured state space sequential (S4) model Gu
et al.
2022 and Mamba [Gu and Dao
2023 Dao and Gu 2024].
Inspired by HiPPO, we propose Online HiPPO SVGP (OHSVGP), by applying the HiPPO framework SVGP in order to leverage the long-range memory modeling capabilities. Our method interprets the
HiPPO time-varying orthogonal projections as inducing variables of an interdomain SVGP Lázaro-
Gredilla and Figueiras- Vidal 2009 Leibfried et al. 2020, Van der Wilk et al. 2020]. where the basis functions are time- dependenn orthogonal polynomials. We show that we are able to significantly
resolve the memory-loss issue in OSVGP, thereby opening up the possibility of applying GPs to long-term online learning tasks. In summary, our contributions include:
(Section We demonstrate that HiPPO integrates into the interdomain GPs by interpreting the HiPPO projections as inducing variables with time-dependent orthogonal polynomial basis functions. This allows the inducing variables compress historical data, capturing long-term information.
(Section 3.2 & [5.]] We show that the kernel matrices can leverage the efficient ODE evolution of the HiPPO framework, bringing extra layer of computational efficiency to OHSVGP.
(Section We demonstrate OHSVGP a variety of online/continual learning tasks includ- ing time series prediction, continual learning on UCI benchmarks, and continual learning in Gaussian process variational autoencoder, showing that it outperforms other online sparse GP baselines in terms of predictive performance, long-term memory preservation, and computational efficiency.
2 Background
In this section, we provide a brief overview of GPs, inducing point methods, online learning with GPs, and Gaussian process variational autoencoders. In addition, we review the HiPPO method, which is the basis of our proposed method.
2.1 Gaussian processes
Let X' be the input space. For time series data, X' = [0, 20), the set of non-negative real numbers. Gaussian process (GP) f GP k) is defined with covariance function k X x X - R. It has the property that for any finite set of input points X = [r1....,I,]] the random vector Kr),,,), where Krr is the kernel matrix with entries 1K(X. X)ly = Kwle = k(,,, 2. For notational convenience and different sets of input points X, and X2, we denote the kernel matrix as Kr or k(X,, X2). The computational and memory complexities of obtaining the GP posterior on X, scale cubically and quadratically respectively. according n2 = | Given responses y and inputs X, a probabilistic model can be defined as y: ~ p(1. fxx)) with GP prior f GP(O, k). where p(y: ((x:) is the likelihood distribution.
f f(X) [f(z))...,, fT
However, for non-conjugate likelihoods, the posterior distribution is intractable, and approximate inference methods are required, such as, but not limited to, variational inference [Titsias] 2009
Hensman et al. 2013 2015a| and Markov chain Monte Carlo (MCMC) (Hensman et al. 2015b
2.2
Variational inference and interdomain Gaussian processes
To address the intractability and cubic complexity of GPs, Sparse Variational Gaussian Processes (SVGP: (Titsias) 2009 Hensman et al., 2013 2015b) cast the problem as an optimization prob- lem. By introducing M inducing points Z E x that correspond to M inducing variables [f(z),.... (z (TT) the variational distribution qf,, u) is defined as q(f, u) := p(f (u)a(u) (u). where (u) is the variational distribution of the inducing variables with parameters 6. Then, the evidence lower bound (ELBO) is defined as
p(y) 2 Ein Eq(ra) (log p(u [f] KL [g (u)]lp(u) =: Lo.
(1)
where v(/)) = fr(s u)go (u)du is the posterior distribution of fa = f(xz). Typical choices for the variational distribution are go (u) Nu M,, Su), where m. and Sn are the free-form mean

---

and covariance of the inducing variables, and yields the posterior distribution:
- - K Ktf: f Kun - S.
(2)
When the likelihood conjugate Gaussian, the ELBO can be optimized in closed form and m, and S can be obtained in closed form (SGPR: Titsias (0099). In addition to setting the in- ducing variables as the function values, interdomain GPs Lazaro-Grediee 2009) propose to generalize the inducing variables to um f(z) (z)dr, where @m (x) are basis functions, allow for further flexibility. This yields [K5u]m = k((x,, x' D0 te')d'' and [Kuulam [] k(x, x²' ) (2) (m) (2' lz'd'''''.
and Figueiras- Vidal
We see that the interdomain SVGP bypasses the selection of the inducing points E RM and reformulates it with the selection of the basis functions Φ.. The basis functions dictate the structure of the kernel matrices, which in turn modulate the function space of the GP approximation. In contrast, SVGP relies on the inducing points Z. which can shift locations according to the training data. Some examples of basis functions include Fourier basis functions [Hensman et al., 2018) and the Dirac delta function 8z,.. the latter recovering the standard SVGP inducing variables.
2.3 Online Gaussian processes
In this paper, we focus on online learning with GPs, where data arrives sequentially in batches (X., Yt, (Xrz1 .t.)....eee etc. For example, in the time series prediction setting, the data arrives in in- tervals of (1), (2),6))... etc. The online GP learning problem is to sequentially update the GP pos- terior distribution data arrives. Suppose that we have already obtained Pt, (ylf)),, (flu,, )g,, (ur,) of the likelihood and variational approximation (with inducing points Z, ). from the first data batch (Xe,, ).,) - Online SVGP (OSVGP; [Bui et al. 2017) utilizes the online learning ELBO
Σ R2g
E (x,)
log (vi fi)l KL (āt, (u,,) Pt, (ur))
KI. (i, qt, (u,,) KI. (4t, (uz²) Ptg (urz))
where yi E Yt for i == 1.... "t, and ttg (u,,) )= == Ptg (u,, [u, )t2 (Utg )du; Unfortunately, with more and more tasks, OSVGP may not capture the long-term memory in the data since data arrives, it is not guaranteed that the inducing points after optimization can sufficiently cover all the previous tasks' data domains.
2.4
Gaussian process variational autoencoders
Gaussian processes can be embedded within a variational autoencoder ( (VAE; (Kingma and Welling
2014) framework, giving rise to the Gaussian process variational autoencoder (GPVAF; Casale et al.
Fortuin et al. 2020 Ashman et al.] [2020] Jazbec et al.) 202T] Zhhu et al. 2023 For sparse GPs with inducing variables, Jazbec et al. [[021] introduced the SVGPVAE, which combines the sparse variational GP (SVGP) with the VAE formulation. The likelihood ply (P())) is parameterized by decoder network ype, which takes GP latent draws input, together with the variational inducing posterior go (u] y). This posterior, go ((y)), is parameterized by the encoder network d. Finally, the latent GP f is typically modeled as multi-on aulti-outpuu GP with independent components. GPVAES have been shown to successfully model high-dimension
2018
time series such as weather data and videos [Zhu et al., 2023 Fortuin al.) 2020). In this work, we consider the SVGPVAE model defined in Jazbec et al. [02 for one set of our experiments, and the detailed specification of the model and training objective can found Appendix D.2
2.5 HÎPPO: recurrent memory with optimal polynomial projections
The HiPPO framework Gu et al. 2020 provides mathematical foundations for compressing continuous-tim signals into jections. Given a time series y(!), HiPPO maintains memory state c(l) E R M that optimally approximates the historical signal (r()))=< The framework consists of time-dependent measure (x) over (-oo, |] that defines input importance, along with normalized polynomial basis func- tions {gh (1) (x) } M- that orthonormal under (x)), satisfying f" gm (x)g
memory states through optimal polynomial pro-
(2) (!) (6) (t)
3

---

8mn. The historical signal is encoded through projection coefficients given by (t) = L-x v(z)gn (z)) This yields the approximation y(a) EN (t)))n(())) for x (-o0, t), minimizing the Le-error S" [lv(x) -x. on (z))l (x)dx. Differentiating c(t)
[co(t)...,,cM ( induces a linear ordinary differential equation a: c(t) = (t) + (t)) with matrices (t). (t) encoding measure-basis dynamics. Discretization yields the recurrence of the form C₁ = A, + Bey enabling online updates. The structured state space sequential (S4) model [Gu et al. 2022 extends HiPPO with trainable parameters and convolutional kernels, while Mamba [Gu and Dao 2023, Dao and Gu, 2024] introduces hardware-aware selective state mechanisms, both leveraging HiPPO for efficient long-range memory modeling. HiPPO supports various measure-base configurations [Gu et al. 2020 2023). A canonical instantiation, HiPPO-LegS, uses a uniform measure w (x) = } (0.4 with scaled Legendre polynomials adapted to [o, d).
In (a) = (2m + 1)//2 Pm (2r- -1) This uniform measure encourages HiPPO-Leeg to keep the whole past in memory.
3 Interdomain inducing point Gaussian processes with HiPPO
We bridge the HiPPO framework with interdomain Gaussian processes by interpreting HiPPO's state vector defined by time-varying orthogonal projections as interdomain inducing points. This enables adaptive compression of the history of a GP while preserving long-term memory.
3.1 HiPPO as interdomain inducing variables
Recall that in interdomain setting in Section 2.2 inducing variables are defined through integral transform against a set of basis functions. Let f ~ GP(O, k), and consider time-dependent basis functions 2 (x) = 9nd (x)) where g. are the orthogonal functions of HiPPO and (() is the associated measure. We define the corresponding interdomain inducing variables as (( - f (t) (x)dz, which not a standard random variable as in Section 2.2 Rather, it is random functions (i.e. stochastic processes) functions. These inducing variables adapt in time, capturing long-range historical information in compact form via HiPPO's principled polynomial projections.
time (um Um um due time-de dependent basis
3.2 Adapting the kernel matrices over time
When new observations arrive at later times in a streaming scenario, we must adapt both the prior cross-covariance Kf, and the prior covariance of the inducing variables K In particular, the basis functions in our HiPPO construction evolve with time, so the corresponding kernel quantities also require updates. Below, we describe how compute and update these matrices at time fz given their values at time t, For clarity, we first discuss Kfn, then K.n-
Prior cross-covariance K Recall that for single input In. the prior cross-covariance with the m-th inducing variable is (0) = f (In, om (x)dæ. We can compute the temporal k(t) evolution of fu in manner consistent with the HiPPO approach, leveraging the same parameters (t) and (t). Specifically,
K (1)
A())
K (0)
+ BB())
n,
where K is the n-th of K (t) fu The matrices A(t) and B(t) depend the specific choice of the HiPPO measure and basis functions. In our experiments, we employ HiPPO-LegS, whose explicit matrix forms are provided in Appendix A One then discretizes in (e.g. using an Euler method or a bilinear transform) to obtain recurrence update rule.
n.:
Prior inducing covariance (t)
The ém-th th element of the prior covariance matrix for the inducing ond (x') 'zdd'' Since k(æ, x') depends on both = If (ax, x') (0 (x)
(t)
variables given by K
fm

---

( (1
6P (4)
GP
fi
data
65% (2
New data.
69%00 (2)
OHSVGP
a faiteh tei
0,60 - faut faa
*04")
Conft. giwitt
ats"")
Recon, from gix
s14"))
Recon
Recon, from glell)
1/2
Figure 1: Online HiPPO Sparse Variational Gaussian Process (OHSVGP) toy time series with tasks. Here x is used to denote arbitrary time index. (a) Time-dependent basis functions with end time index l and 12. (b) Evolution of optimal approximate posterior of inducing variables (mean +22 marginal standard deviation). (c), (d), (e) illustrate predictive mean +-2 standard deviation of posterior online GP, OHSVGP and finite basis reconstruction of posterior OHSVGP, respectively.
and x', a recurrence update rule based on the original HiPPO formulation, which is designed for single integral, can not be obtained directly for KM Fortunately, for stationary kernels, Bochner Theorem [Rudin, 1994] can be applied factorize the double integrals into two separate single integrals, which gives rise to Random Fourier Features (RFF) approximation [Rahimi and Recht 2007]: for stationary kernel k< x²) A([x= x') I), RFF approximates it follows: k(x, x') z N Em= [cos (w,.)) cos (W,x²) + sin (w,,.r) sin (w,,²)), where wn, ~ p(w) is the spectral density of the kernel. Substituting this into the double integral factorizes the dependency on x and x', reducing [Kūü) to addition of products of one-dimensional integrals. Each integral, with the form of either cos(sws (.az) (2) (0 (xr) or (in(wax))) (xx dxz, again corresponds to a HiPPO-ODE in time. By sampling multiple random features, updating them recurrently to time i, and averaging,
obtain RFF approximation of k<(t) In addition, more advanced Fourier feature approximation techniques (e.g.. [Ton et al., 2018)) ) can be leveraged for non-stationary kernels. The details of the ODE for recurrent updates of the RFF samples appear in Appendix B. Alternatively, one may differentiate k directly with respect to t. This yields a matrix ODE of the form different from the original HiPPO formulation. For details, see Appendix [B.] Empirically, vanilla implementation of this approach shows numerical unstability. Hence, we conduct our experiments based on RFF approximation.
Sequential variational updates. Having obtained (f2) (1(ta) at new time t2 > f1, we perform fu variational updates following the online GP framework described in Section 2.3 This ensures the posterior at time to remains consistent with both the new data and the previous posterior at time t1 based on K fu K(tt)) Overall, this procedure endows interdomain HiPPO-based GPs with the ability to capture long-term memory online. By viewing the induced kernel transforms as ODEs in time, we efficiently preserve the memory of past observations while adapting our variational posterior in online fashion. Figure Pp illustrates the evolution of the optimal posterior q(u(z)) time increases on a toy online time series regression problem with two tasks, where x determines the end of the recurrent update for the prior cross and inducing covariance matrices (evolved up K(#t) and K(t), respectively). Furthermore, when I > t1, we will update (((())) online with the two data points from the second task by optimizing the online ELBO (Eq. 33). which gives the discrete jump at x = f1. Figure TI shows the posterior OHSVGP compared with the fit of the gold-standard
online GP in Figure 1 Notably, if f ~ qe (f), then q(u!t)) 1 f f(x) (æ)dr (detailed derivation in Appendix C Therefore, our framework also provides finite basis approximation of the posterior (t)
OHSVGP as a byproduct: / = ΣM- um gi (z)), ( ~ q(u()) (()) ). Figure le plots the finite basis
approximation//rcoossttttooo and it close to the posterior OHSVGP for this simple example.
5

---

3.3 Extending OHSVGP to multidimensional input
For multidimensionall input data, suppose there is a time order for the first batch of training points with inputs {ath [=:)), )== such that x;; appears after ) if > j. and we further assume 2; appears at time index i (i.e., x(i = x ()), where At a user-specified constant step size. In this case, we can again obtaín interdomain prior covariance matrices via HiPPO recurrence. For example, forward Euler method applied to the ODE in Eq. for K yields
(5) The equation above can be viewed as a discretization (with step size At) of an ODE solving path integrals of the form Io At k [an
K fu
((+1))4)) m, [I + At
fu (150) In. + AtB(i (1) (x
xx
(1) ar(s) om (s)ds. The i-th training input x²" is assumed to be
z (1) == (iLt) and thus the path integral is approximately solved with discretized recurrence based on the training inputs corresponding to {a(istt) We continue the recurrence for the second task with ordered training inputs {>/ { |n by assigning time index (N1 + i) to its i-th instance. and keep the recurrence until we learn all the tasks continually. In practice, one may use a multiple of Át as the step size to accelerate the recurrence, e.g., instead of using all the training inputs, one can compute the recurrence based on {æ1, T3, 4s, ...) only by using step size 2 When there no natural time order for training instances in each task, such as in standard continual learning applications, we need to sort the instances with some criterion to create pseudo time order to fit OHSVGP, similar to the practice of applying SSMs to non-sequence data modalities, e.g., SSMs, when applied to vision tasks, assign order to patches in an image for recurrence update of the memory [Zhu et al. 2024]. In our experiments, we show that the performance of VHPV when applied to continual learning, depends the sorting criterion used.
4 Related work
Online sparse GPs. Previous works mainly focus on reducing the sparse approximation error with different approximate inference techniques, such as variational inference [Bui et al. 2017
Maddox al., 202T), expectation propagation [Csató and Opper 2002. Bui et al. 2017] Laplace approximation [Mddox et al., [222]] J. and approximation enhanced with replay buffer [Chang et al.)
2023). The orthogonal research problem of online update of inducing points remains relatively underexplored, and Hivoted-Choleske
(Burt et al. 2019) as deployed in Maddox et al. [2021] Chang al. 2023] is one of the most effective approaches for online update inducing points up to date. We tackle this problem by taking advantage of the long-term memory capability of HiPPO to design an interdomain inducing variable based method and the associated recurrence based online update rules. Notably, our HiPPO inducing variables in principle are compatible with all the aforementioned approximate inference frameworks since only the way of computing prior covariance matrices will be different from standard online sparse GPs.
Inferdomain GPs. To our knowledge, OHSVGP the first interdomain GP method in the context of online learning. Previous interdomain GPs typically construct inducing variables via integration based on predefined measure (c.g.. uniform measure over a fixed interval [Hensman et al. 2018) or a fixed Gaussian measure [Lázaro-Gredilla and Figueiras- Vidal 2009 to prevent diverging covariances, and this predefined measure may not cover all regions where the time indices from future tasks are, making them unsuitable for online learning. In contrast, OHSVGP bypasses this limitation by utilizing adaptive basis functions constructed based on time-dependent measure which keeps extending to the new time region as more tasks arrive.
Markovian GPs. Markovian GPs (ääkkkä and Solin 2019 Wilkinson et al.) 2021 have similar recurrence structure during inference and training due to their state space SDE representation. However, Markovian GPs are popularized due their O(x) computational complexity and is not explicitly designed for online learning.
Experiments
Applications & datasets.
We evaluate OHSVGP against baselines in the following tasks.

---

Time series prediction. We consider regression benchmarks, Solar Irradiance Lean 2004 and Audio Signal (Bui and Turner 2014] produced from the TIMIT database Gaaifolo et al.
[993]. We preprocess the two datasets following similar procedures described in Gal
and Turner 2015] and BBui et al. [017]. respectively (the train-test split different due to random splitting). In addition, we consider daily death-count time series from Santa Catarina State, Southern Brazil spanning the March 2020 to February 2021 COVID-19 pandemic, obtained from Hawryluk et al. [2021]. We construct online learning tasks by splitting each dataset into IO (5 for COVID) sequential partitions with an equal number of training instances.
Continual learning. We consider continual learning on two UCI datasets with multi-dim inputs, Skillcraft [Blair et al., 2013] and Powerplant (Teekii and Kaya 2014], using the same data preprocessing procedure as in Stanton et al. [021]] We construct two types of continual learning problems by first sorting the data points based on either the values in their first dimension or their 1.2 distance from the origin, and then splitting the sorted datasets into 10 sequential tasks with an equal number of training instances.
High dimensional time series prediction. We evaluate GPVAEs on hourly climate data from ERA5 [Copernieus Climate Change Service, Climate Data Store 2023 Hersbach et al., 2023]. comprising 17 variables across randomly scattered locations around the UK from January 2020 onward. The dataset is split into 10 sequential tasks of 186 hourly time steps each.
Baseline. We compare OHSVGP with OSVGP [Bui et al., 2017 and OVC (Online Variational Conditioning: [Maddox et al.) 202T)) D. At the beginning of each task, OSVGP initialize the induing points by sampling from the old inducing points and the new data points, while OVC initializes them via pivoted-Cholesky [Burt et al. 2019] and we consider both fixing the initialized inducing points as in Chang et al. [023] (OVC) or keep training them as in Maddox et al. [2021] (OVC-optZ). For time series regression with Gaussian likelihood, we consider OHSGR and OSGPR (OHSVGP and OSVGP based on closed form ELBO), and we further consider OVFF (OSGPR based on variational Fourier feature (VFF), an interdomain inducing point approach from Hensman al. 2018). )
Hyperparameters. Within each set of experiments, all the models are trained using Adam [Kingma and Ba. 2015) with the same learning rate and number of iterations. For OHSVGP, we construct inducing variables based on HiPPO-LegS [Gu et al. 2020] (see Appendix F:4 for visualizations of using other HiPPO variants) and use 1000 RFF samples. We use ARD-RBF kernel, except for OVFF tailored specifically to Matérn kernels, where we use Matérn- kernel instead. Similar to Maddox al. [022]], we do not observe performance gain by keeping updating kernel hyperparameters online, and we include results with trainable kernel byperparameters in Appendix F2 for time series regression, but the performance becomes unstable when number of tasks is large. Thus, we either only train the kernel hyperparameters during the initial task and keep them fixed thereafter (Section 5.3) or obtain them from a full GP model trained over the initial task. It is also worth noting that OVFF requires computing covariances as integrals over a predefined interval covering the whole range of the time indices from all tasks (including unobserved ones), which impractical in real online learning scenarios. For our experiments, we set the two edges of this interval to be the minimum and maximum time index among the data points from all the tasks, respectively.
Evaluations & metrics. We report results in Negative Log Predictive Density (NLPD) in the main text, and Root Mean Squared Error (RMSE) in Appendix (expected calibration error (ECE; [Guo] al. 20177) for COVID data instead), which shows consistent conclusions as NLPD. We report the mean and the associated 95% confidence interval obtained from 5 (3 for experiments on ERAS) independent runs.
5.1 Online time series prediction
Time series regression. Figure D NLPD (over the past tasks) of different methods during online learning through the 10 tasks for Solar Irradiance and Audio dataset. Overall, OHSGPR consistently achieves the best performance with OVC performing competitively, especially as we learn more and more tasks, suggesting OHSGPR effectively preserves long-term memory through its HiPPO-based memory mechanism. OSGPR shows catastrophic forgetting starting around task
7

---

050M
0017
(a) Solar, M=50
(b) Solar, M=150
(c) Audio, M=100
(d) Audio, M=200
Figure 2: Test set NLPD over the learned tasks number of learned tasks for Solar Irradiance and Audio signal prediction dataset.
Trann
Test
Splits
(a) OSGPR
(b) OVC
(c) OVFF
(d) OHSGPR
Figure 3: Predictive mean +2 standard deviation of OSGPR, OVC, OVFF, and OHSGPR after task 10 of the Solar dataset. M 50 inducing variables are used.
especially when the number of inducing points M is small. Although OVC-optZ also initializes inducing points with pivoted-Cholesky as OVC, with further optimization, its performance starts to degrade starting from task for the audio dataset when M 100, which suggests the online ELBO objective cannot guarantee optimal online update of inducing points that preserve memory. OVFF tends perform well at the later stage. However, during the first few tasks, underfits the data significantly compared with other methods since its inducing variables are computed via integration over a predefined interval capturing the region of all the tasks, which unnecessarily long and suboptimal for learning at the early stage.
In Figure Bl compare the final predictive distributions for different methods after finishing online learning 10 tasks of Solar Irradiance. The inducing points Z for OSGPR tend to move to the regions where the later tasks live after online training, and the prediction of OSGPR in the initial regions without sufficient inducing points becomes close to the uninformative prior GP. In contrast, OHSGPR maintains consistent performance across both early and recent time periods.
Infectious disease modeling We replace the Gaussian likelihood with a non-conjugate Negative Binomial likelihood to capture the over-dispersion in COVID-19 death counts. All methods use M E {15, 30) inducing points and are trained for 5000 iterations per task with learning rate of 0.01. Figure 4 reports the change of NLPD through online learning for the first four out of five tasks. The wide metric variance reflects the noisy nature of death-count data as it is difficult to accurately track down COVID-19 death counts. OHSVGP achieves the best performance overall while OSVGP forgets Task with small M.
Runtime comparison. Table 1 shows the accumulated wall-clock runtime for different meth- ods to learn all the tasks. Unlike OSVGP and OVC-optZ, which must iteratively optimize induc- ing points (for which train e.g., 1000 iterations for time-series regression tasks), OHSVGP, OVFF (both based on interdomain inducing points), and OVC (based on one-time pivoted-Cholesky update of inducing points for each task) bypass this cumbersome optimization. In particular,
OHSVGP recurrently evolves K and Kuu for each new task. For regression problems where closed- form posterior can be obtained, OHSGPR requires no training at all. As result, OHSGPR, OVC and OVFF run significantly faster, adapt- ing to all tasks within a couple of sec- onds for Solar Irradiance and Audio data. For COVID data, even when free-form variational parameters of
Table 1: Wall-clock accumulated runtime for learning all the tasks on a single NVIDIA RTX3090 GPU in seconds, of dif- ferent models for time series prediction experiments.
Solar Irradiance
Audio Data
COVID
M
M
M
Method
50
150
100
200
15
30
OSGPR/OSVGP
140
149
144
199
525
530
OVC
0.450
0.620
0.558
0.863
345
360
0.327
0.295
0.356
OHSGPRJOHSVGP
0.297
0.394
0.392
0.655
370
380

---

OSVGP
OVC optz
CHSVGP
CVC
OHEVGP
Task
Took
Task
Teck
2,4
3.4
2.4
3.1
2.1
3.2
3.2
3.0
3.2
2.0
3,0
Tasks
Tasks
Tasks
Tasks
Tasks
(a) M=15
(b) M=30
Figure 4: Test set NLPD on COVID dataset right after learning Task and after learning all the tasks.
OSVEP
ovc
CHSVEPK
CHSSVPPP
OSVGP
OHSVEP-
CHS4GRP
Tosk
NLPO
1.3
1.5
2.5
2.5
1.2-
1.0
1.0
10
10
10
10
10
10
Tasks
Tasks
Tasks
Tasks
Tisks
Tasks
Tasks
(a) Skillcraft (Ist dimension)
(b) Skillcraft (L.2)
0SVGP
DVC Task: 2
DRSVCP-K
CHSVGG--
DSVGP
ove Task
OHEVEE--
OHAVEPO
Task 1
Task
Task
Tb.sk
Taskc
Task
4.0
0.7
0.0
2.0
0.2
0.0
0.5
01
0.0
0.0
10
10
Tasks
Tasks
Tasks
Tasks
Tasks
Tasks
Tasks
Tasks
(c) Powerplant (1st dimension)
(d) Powerplant (1.2)
Figure 5: Test set NLPD after continually learning Task i and after learning all the tasks for 1,
2, 4, 8. Tasks are created by splitting Powerplant and Skillcraft datasets with inputs sorted either according to the 1st input dimension or L2 distance to the origin).
inducing variables are learned using uncollapsed ELBO, OHSVGP and OVC are still significantly faster than OSVGP since no gradient computation is required for the inducing points.
5.2 Continual learning on UCI datasets
We use 256 inducing variables for all methods, and for each task, we train each method for 2000 iterations with a learning rate of 0.005. We only consider OVC here since initial trials show OVC-opiZ give worse results these two datasets. As described in Section [.3] within each task, OHSVGP requires sorting the data points to compute prior covariance matrices via recurrence. We consider two sorting criteria. The first one, which we call OHSVGP-o, uses the oracle order compatible with how the tasks are created (c.g., sort with L2-distance to the origin if the tasks are initially splitted based on it). In real- world problems, we typically do not have the information on how the distribution shifts from task to task. Hence, we also consider OHSVGP-K, which uses a heuristic sorting method based on kernel similarity: we select the i-th point in task j to be (5) arg max k(x,, x (i)
for i 1, and the first point in first task is set to be zz(1) compares the two variants of OHSVGP with OSVGP and OVC. Overall, OS achieves the worst performance and is again prone to forgetting the older tasks, especially in Figure 5c OVC performs
arg max (1)) k(x, 0). Figure
decently for Skillcraft but it also demonstrates catastrophic forgetting in Figure 5c While OHS VGP-k achieves similar performance other methods across all 4 scenarios, suggesting the importance of a sensible sorting method when applying OHSVGP for continual learning. Here, we only report the results for Task 2, 4, and for concise presentation, and in Appendix [.I] we include the complete results for all the tasks (the overall conclusion the same). In Appendix F3 we further visualize how different sorting methods impact OHSVGP's performance in continual learning with 2D continual classification problem.
OSVGP on Skillcraft, OHSVGP-0 consistently outperforms the
5.3 Continual learning for high dimensional time series prediction
All models share a two-layer MLP encoder-decoderr 20-dimensional latent space, and a multi-output GP with independent components; we use M E (50, 100) and train each task for 20 epochs with learning rate 0.005 on single NVIDIA A6000 GPU. The continual learning in SVGPVAE is achieved

---

OVC-per-dim Task
GHSVGP Task
OSVGP
OVC-per-dim Task
OHSvGP
Task 1
Taskk
Task
Task 2
Task 8
10-
1
10
16
10
10
10
10
10
Tasks
Tasks
Tasks
Tasks
Tasks
Tasks
Tasks
Tasks
(a) M=50
(b) M=100
Figure 6: Test set NLPD after continually learning Task i and after learning all the tasks for i = 1, 2. 4, 8, on ERA5 dataset.
by further imposing Elastic Weight Consolidation (EWC; KKirkpatrick et al. 2017 loss the encoder and decoder, which yields the vanilla baseline, Online SVGPVAE (OS VGP). Since EWC alone leaves inducing locations non-regularized, a principled online placement rule for the inducing points will improve the model. Thus, we further consider OVC-SVGPVAE (OVC) which adjusts inducing points online via Swvoted-Choleske and OVC-SVGPVAR per dimension (OVC-per-ding) which makes OVC more flexible by allocating a separate set of M inducing points to every latent dimension. Our method, Online HiPPO SVGPVAE (OHSVGP), replaces standard inducing points in SVGPVAE with HiPPO inducing variables and updates them online via recurrence. Figure plots the change of NLPD during continual learning for Task 1, 2, and 8 (full results in Appendix [.II. The performance of OHSVGP remains stable throughout, while the other methods all demonstrate obvious catastrophic forgetting shown by the large gaps between performances after learning current task i and after learning final task 10. Two factors plausibly explain the gap: first, standard inducing points cannot adequately cover the long time axis, whereas OHSVGP ties its inducing variables to basis functions rather than time locations; second, the added secoder-decode complexity makes optimization harder for models that must reuse limited inducing set. Increasing M narrows the gap but scales at (M3) computational and (M2) memory cost respectively, underscoring OHSVGP's superior efficiency.
Conclusion
We introduce OHSVGP, a novel online Gaussian process model that leverages the HiPPO framework for robust long-range memory in online/continual learning. By interpreting HiPPO's time-varying orthogonal projections as adaptive interdomain GP basis functions, we leverage SSM for improved online GP. This connection allows OHSVGP to harness HiPPO's efficient ODE-based recurrent updates while preserving GP-based uncertainty- aware prediction. Empirical results on a suite of online and continual learning tasks show that OHSVGP outperforms existing online sparse GP methods, especially in scenarios requiring long- tengttermmo memory. Moreover, its recurrence- based covariance updates yield far lower computational overhead than OSVGP's sequential inducing point optimization. This efficient streaming capability and preservation of historical information make OHSVGP well-suited for real-worid applications demanding both speed and accuracy.
Broader impact. This paper presents work whose goal is to advance machine learning research. There may exist potential societal consequences of our work, however, none of which we feel must be specifically highlighted here.
Acknowledgments and Disclosure of Funding
Samir Bhatt acknowledges funding from the MRC Centre for Global Infectious Disease Analysis (reference MR/X020258/1), funded by the UK Medical Research Council (MRC). This UK funded award carried out in the frame of the Global Health EDCTP3 Joint Undertaking. Samir Bhatt acknowledges support from the Danish National Research Foundation via a chair grant (DNRF160) which also supports Jacob Eurran-Sebastian. Samir Bhatt acknowledges support from The Eric and Wendy Schmidt Fund For Strategic Innovation via the Schmidt Polymath Award (G-22-63345) which also supports Harrison Bo Hua Zhu. Samir Bhatt acknowledges support from the Novo Nordisk Foundation via The Novo Nordisk Young Investigator Award (NNF200C0059309).
10

---

References
Matthew Ashman, Jonathan So, Will Tebbutt, Vincent Fortuin, Michael Pearce, and Richard E Turner.
Sparse Gaussian process variational autoencoders. arXiv preprint (tt0:20000/01777 2020.
Mark Blair, Joe Thompson, Andrew Henrey, and Bill Chen. SkillCraft Master Table Dataset. UCI Machine Learning Repository, 2013. DOI: https://bi.ogg/1122444////1111
Thang D. Bui and Richard E. Turner. Tree-structured Gaussian process approximations, In Advances in Neural Information Processing Systems, 2014.
Thang D. Bui, Cuong V. Nguyen, and Richard E. Turner. Streaming sparse Gaussian process approximations. In Advances Neural Information Processing Systems, 2017.
David R. Burt, Carl E. Rasmussen, and Mark van der Wilk. Rates of convergence for sparse variational Gaussian process regression. In International Conference Machine Learning (ICML), 2019.
Francesco Paolo Casale, Adrian Dalca, Luca Saglietti, Jennifer Listgarten, and Nicolo Fusi. Gaussian process prior variational autoencoders. Advances in neural information processing systems, 31, 2018.
Paul E. Chang, Prakhar Verma, S.T. John, Arno Solin, and Mohammad Emtiyaz Khan. Memory- based dual Gaussian processes for sequential learning. In International Conference on Machine Learning. 2023.
Copernicus Climate Change Service, Climate Data Store. Era5 hourly data on single levels from 1940 to present, 2023. URL https///:/// org/ 10. 24381/cd2 cds. adbb2d477
Accessed: DD-MMM-
YYYY.
Lehel Csató and Manfred Opper. Sparse on-line Gaussian processes. Neural Computation, 14(3): 641-668, 2002.
Tri Dao and Albert Gu. Transformers are SSMs: Generalized models and efficient algorithms through structured state space duality. In International Conference on Machine Learning (ICML), 2024.
Seth Flaxman, Swapnil Mishra, Axel Gandy, Juliette T Unwin, Thomas A Mellan, Helen Coupland, Charles Whittaker, Harrison Zhu, Tresnia Berah, Jeffrey W Eaton, et al. Estimating the effects of non-pharmaceutical interventions on covid-19 19 in europe. Nature, 584(7820):257-61,, 2020.
Vincent Fortuin, Dmitry Baranchuk, Gunnar Rätsch, and Stephan Mandt. GP-:AE: Deep probabilistic time series imputation. In International Conference on Artificial Intelligecce and Statistics, pages 1651-1661. PMLR, 2020.
Yarin Gal and Richard E. Turner. Improving the Gaussian process sparse spectrum approximation by representing uncertainty in frequency inputs. In International Conference on Machine Learning (ICML), 2015.
Yaroslav Ganin, Evgeniya Ustinova, Hana Ajakan, Pascal Germain, Hugo Larochelle, François Laviolette, Mario March, and Victor Lempitsky. Domain-adversarial training neural networks. Journal of Machine Learning Research, 17(59):1-35, 2016. URL http://jmlr: org/ papp//aprs/// v17/15- 239. html
J. Garifolo, L. Lamel, W. Fisher, J. Fiscus, D. Pallett, N. Dahlgren, and V. Zue. TIMIT acoustic- phonetic continuous speech corpus LDC93S1. In Philadelphia: Linguistic Data Consortium, 1993.
Albert Gu and Tri Dao. Mamba: Linear-time sequence modeling with selective state spaces, arXiv preprint arXiv:2312.00052,, 2023.
Albert Gu, Tri Dao, Stefano Ermon, Atri Rudra, and Christopher Ré. HiPPO: Recurrent memory with optimal polynomial projections. In Advances in Neural Information Processing Systems, 2020.
Albert Gu, Karan Goel, and Christopher Ré. Efficiently modeling long sequences with structured state spaces, In The International Conference on Learning Representations (ICLR), 2022.
11

---

Albert Gu, Isys Johnson, Aman Timalsina, Atri Rudra, and Christopher Re. How to train your HIPPO: State space models with generalized orthogonal basis projections. In Internationaal Conference on Learning Representations, 2023. URL https::/ http:///....ppoeeeeee/e//////tta ret//foouummxd=x211000008
Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q. Weinberger. On calibration of modern neural networks. In International Conference on Machine Learning, 2017.
Iwona Hawryluk, Henrique Hoeltgebaum, Swapnil Mishra, Xenia Miscouridou, Ricardo P Schneken- berg, Charles Whittaker, Michaela Vollmer, Seth Flaxman, Samir Bhatt, and Thomas A Mellan. Gaussian process noweasting: application to covid-19 mortality reporting. In Uncertainty in Artificial Intelligence, pages 1258-1268 PMLR, 2021.
James Hensman, Nicolò Fusi, and Neil Lawrence. Gaussian processes for big data. In Proceedings of the Twetty-Ninth Conference on Uncertainty in Artificial Intelligence, UAI' 13, pages 282-290, Arlington, Virginia, USA, 2013. AUAI Press.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. In Artificial Intelligence and Statistics, pages 351-360. PMLR, 2015a.
James Hensman, Alexander G Matthews, Maurizio Filippone, and Zoubin Ghahramani. MCMC for variationally sparse Gaussian processes. Advances in Neural Information Processing Systems, 28, 2015b.
James Hensman, Nicolas Durrande, and Arno Solin. Variational Fourier features for Gaussian processes. Journal of Machine Learning Research, 18(151):1-52, 2018.
H. Hersbach, B. Bell, P. Berrisford, G. Biavati, A. Horányi, J. Muñoz Sabater, J. Nicolas, C. Peubey, R. Radu, I. Rozum, D. Schepers, A. Simmons, C. Soci, D. Dee, and J.-N. Thépaut. Era5 hourly data single levels from 1940 to present, 2023, Accessed: DD-MMM-YYYY.
Roger A. Horn and Charles R. Johnson. Topics in Matrix Analysis. Cambridge University Press, 1991.
Metod Jazbec, Matt Ashman, Vincent Fortuin, Michael Pearce, Stephan Mandt, and Gunnar Rätsch. Scalable Gaussian process variational autoencoders. In International Conference on Artificial Intelligence and Statistics, pages 351 1-3519. PMLR, 2021.
Sanyam Kapoor, Theofanis Karaletsos, and Thang D. Bui. Variational auto-regressive Gaussian processes for continual learning. In International Conference on Machine Learning. 2021.
D. P. Kingma and M. Welling. Auto-encoding variational Bayes. In International Conference on Learning Representations, 2014.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International
Conference on Learning Representations, 2015.
James Kirkpatrick, Razvan Pascanu, Neil Rabinowitz, Joel Veness, Guillaume Desjardins, Andrei A Rusu, Kieran Milan, John Quan, Tiago Ramalho, Agnieszka Grabska- Barnkkaa-Baiiiaaaa..a et al. Overcoming catastrophic forgetting in neural networks. Proceedings of the national academy of sciences, 114 (13):3521-355 3526,, 2017.
Judith Lean. Solar irradiance reconstruction. In Data contribution series # 2004-035, IGBP PAGES/World Data Center for Paleoclimatology NOAA/NGDC Paleoclimatology Program, Boul- der, CO, USA, 2004.
Felix Leibfried, Vincent Dutordoir, ST John, and Nicolas Durrande. A tutorial on sparse Gaussian processes and variational inference. arXiv preprint arXi::2012.33962 2020.
Miguel
Vigdal. and Anibal Figueiras- Vidal. Inter-domain Gaussian processes for sparse inference using inducing features. In Advances in Neural Information Processing Systems, 2009.
Wesley J. Maddox, Samuel Stanton, and Andrew Gordon Wilson. Conditioning sparse variational Gaussian processes for online decision-making In Advances Neural Information Processing Systems, 2021.
12

---

Mélodie Monod, Alexandra Blenkinsop, Xiaoyue Xi, Daniel Hebert, Sivan Bershan, Simon Tietze, Mare Baguelin. Valerie C Bradley, Yu Chen, Helen Coupland, al. Age groups that sustain resurging covid-19 19 epidemics in the united states. Science, 37166336)) 2021.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems, 2007.
Stephen Roberts, Michael Osborne, Mark Ebden, Steven Reece, Neale Gibson, and Suzanne Aigrain. Gaussian processes for time-series modelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1984):20110555 2013.
W. Rudin. Fourier analysis on groups. Wiley Classics Library. Wiley-Interscience New York, reprint edition, 1994.
Simo Särkkä and Arno Solin. Applied stochastic differential equations, volume 10. Cambridge University Press, 2019.
Samuel Stanton, Wesley J. Maddox, Ian Delbridge, and Andrew Gordon Wilson. Kernel interpolation for scalable online Gaussian processes. In International Conference on Artificial Intelligence and Staristics. PMLR, 2021.
Pnar Tfekci and Heysem Kaya. Combined Cycle Power Plant. UCI Machine Learning Repository, 2014. DOL: https://di..or//00 http://di.orr/1..2445///// 20t/1000443//00002..
Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial intelligence and statistics, pages 567-574. PMLR, 2009.
Jean-Francois Ton, Seth Flaxman, Dino Sejdinovic, and Samir Bhatt. Spatial mapping with gaussian processes and nonstationary fourier features. Journal of Spatial Statisties, 28:59-78, 2018.
H Juliette T Unwin, Swapnil Mishra, Valerie C Bradley, Axel Gandy, Thomas A Mellan, Helen Coupland, Jonathan Ish-Horowicz, Michaela AC Vollmer, Charles Whittaker, Sarah L Filippi. et al. State-level tracking of covid-19 in the united states. Nature communications, 11(1):6189, 2020.
Mark Van der Wilk, Vincent Dutordoir, ST John, Artem Artemev, Vincent Adam, and James Hensman. A framework for interdomain and multioutput Gaussian processes. arXiv preprint arXiv:2003.01 115, 2020.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in Neural Information Processing Systems, 2017.
William Wilkinson, Amno Solin, and Vincent Adam. Sparse algorithms for Markovian Gaussian processes. In International Conference Artificial Intelligence and Statistics, pages 1747-1755. PMLR, 2021.
Harrison Zhu, Carles Balsells Rodas, and Yingzhen Li. Markovian Gaussian process variational autoencoders. In International Conference on Machine Learning, 2023.
Lianghui Zhu, Bencheng Liao, Qian Zhang, Xinlong Wang. Wenyu Liu, and Xinggang Wang. Vision Mamba: Efficient visual representation learning with bidirectional state space model. In International Conference on Machine Learning (ICML), 2024.
13

---

NeurIPS Paper Checklist
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: We accurately describe the paper's 's contributions and scope, and include the necessary motivations and backgrounds required to understand contributions. The claims are verified with our experiments.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
We discuss the effect of a sub-optimal sorting method when applying our
method to continual learning in experiments (Section 3.3 5.2 & Appendix F3).
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "L_imitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and prool's
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
14

---

Answer: [NA]
Justification:
No theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification:
We include the detailed dataset composition and hyperparameter information.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
specific model and empirical evaluation, may
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with open-source dataset instructions for how construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
15

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
We have released our code at HLLPOSVGP/tree///g
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co
ooliiicy) for more details.
While encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https //nips wwtpppplccgguiieeeoccceococoecoe http//uddees//eee (eooodeSumisisioonooeo
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Question:
Answer: [Yes]
We fully describe the training and test details, including data splits, random seeds, hyperparameter tuning, optimizer type and any necessary information needed for reproducibility in our experiments section.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification; runs.
We include the error bars of the performance metries over several random seed
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
16

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
If error bars reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: We indicate which GPUs were used to run the experiments. We also include wallclock time information (T.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't itdn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics http:///neaarrp..cccpppp//e/aaaa
Answer: [Yes]
Justification: In this paper, we introduce work designed to push the boundaries of machine learning. While our methods could carry ethical implications, we do not believe any require specific discussion at this point in the submission process.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: We include statement for broader impacts the end of the paper.
Guidelines:
17

---

The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g.. gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: This paper poses such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: Whenever we use any assets, we always cite the original asset source.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
18

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, datasets.paperssiii www.ppppptwwhhooooo thcode www.coo//aaa//o candsatasss
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We fully describe our contributions and any assets used in the paper. We will also be releasing code after the publishing of the paper.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/coodemmodll as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and sereenshots, if applicable, as well details about compensation (if any)?
Question:
Answer: [NA]
Justification: This paper does not involve crowdsourcing research with human subjects.
Guidelines:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
19

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if' applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: LLM was only used for editing the writing, helping with plotting and assisted coding.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy otttop////eerrrrsscc//////oeerrp / for what should or should not be described.
20

---

A HiPPO-LegS matrices
Here we provide the explicit form of matrices used in our implementation of HiPPO-LegS Gu et al. 2020]. For given time i, the measure w.tt (x) [. (x) and basis functions Øm
(x) 9m (z) <2m-1 P... (2r 1) [0.4] (a) are used, where (.) is the m-th Legen- dre polynomial and (0.0 (2:) is the indicator function on the interval [0, 1]. These basis functions are orthonormal, i.e.,
(1)
L 9m (z)gm (x) dr 8mn
(6)
Following Gu et al. 2020 the HiPPO-LegS framework maintains coefficient vector c(t) RMx1 that evolves according to the ODE:
(1)) (1))(()) + (())((1)
dt
where f(t) is the input signal at time t. The matrices A()) E RMXM and B(t) E RM are given by:
(2n 1)(2
if k
[t) = 7A, [A]nk
n
if = k
(8)
0
if < k
and
B
+1
[())
These matrices govern the evolution of the basis function coefficients over time, where the factor 1/t reflects the time-dependent scaling of the basis functions to the adaptive interval [0, f]. When discretized, this ODE yields the recurrence update used in our implementation.
B Computing prior covariance of the inducing variables K(())
We provide the detailed derivation for the following two approaches when the inducing functions are defined via HiPPO-LegS. Recall that the fm-th element of the prior covariance matrix for the inducing variables is given by
K ( 1eleem Em = ff k(x, x' (1) (2 d'dddd''
(10)
where (x) (x) (x) are the mareingrying basis functions under the HiPPO-LegS frame- work.
B.1 RFF approximation
Since k(z, x') depends on both x and x', a recurrence update rule based on the original HiPPO formulation, which is designed for single integral, not be obtained directly for K'ft Fortunately, for stationary kernels, Bochner Theorem [Rudin] 1994 can be applied to factorize the above double integrals into separate single integrals, which gives rise to Random Fourier Features (RFF) approximation [Rahimi and Recht 2007]: for a stationary kernel k(æ, x') = k(x x')), RFF approximates follows:
k(x, x Ep(w) cos(wrr) cos(wz') + sin(wx) sin (mxa)
(11)
Jcos (w,,rr cos (wnx') + sin (wn,x) sin (w,
21

---

where Wh p(w) is the spectral density of the kernel. Substituting this into the double integral (Eq. 10p factorizes the dependency on z and x', reducing [Kulken to addition of products of one- dimensional integrals. Each integral based on Monte Carlo sample w has the form of either
T = cos(wxz) & (x) da or
Zu,t = sin(wx)00 (x) dx,
(12)
which corresponds to a standard projection coefficient the HiPPO framework. We further stack these integrals based on M basis functions and define
(0)
(0) W,M
200
w.
2 ur
T
(13)
Collecting N Monte Carlo samples Jn form the feature matrix
Z
Z
Zay
Z
(t)
ZuN
(14)
and the RFF approximation of the covariance
1 (()) (2(0) T
KO
(15)
Since ( and Z are standard HiPPO projection coefficient, their computation is governed by the HiPPO ODE evolution as before
d z (+) A() 2 (t,)
Z A(t 2 B(t)
(16)
dt
with (t) = cos( () and W, (t) = (w,t) and these ODEs can be solved parallel across dif- ferent Monte Carlo samples. In summary, the procedure involves sampling multiple random features, updating them recurrently to time t, and averaging across samples to obtain RFF approximation of Kuu (4)
For non-stationary kemels,, more advanced Fourier feature approximation techniques (c.g., Ton et al. 2018)) can be applied.
B.2 Direct ODE evolution
Differentiating [ K
f.m. with respect to gives
d K (1)) k.m
② (0) (x)o (x') dodr'.
k( x')
oi 01.
(17)
di
Applying the product rule:
dt K k. Jf k(x,z') ② ot (1) (z')drd'' ff k(2, 20 (z) ② al (e')dadr' (18)
(2'
St
In HiPPO-LegS, each 2² (x) obeys an ODE governed by lower-order scaled Legendre polynomials on [0, f] and Dirac delta boundary term at t. Concretely,
/2 +1
σ1 22 (z)
11 (x) V2f - 16 (x)
/22+1
(19)
1 (z) ...1. 5((x))
28
30
3012
where (x(x) is the Dirac delta (see Appendix D.3 in Gu et al. 2020 for details).
Substituting this expression into the integrals yields the boundary terms of the form
kktt, 1 (a') dx', along with lower-order terms involving { Summarizing in matrix form leads to
(t)
(1) K }. etc.
d (!) A(t) K(A(T)] dt
B(t) + B(O)).
(20)
22

---

where the Im--th entry of K RAXX is [Kile, ( E RRXM M is the same lower-triangular
matrix from the HiPPO-LegS framework defined in Eq. and Bt) -Mxx is built from the boundary contributions as
((t) c(t)
(21)
where 1M E 1I is row vector of ones of size M and c(t) E RMx1 is the coefficient vector with each element being
cx(() - f l: (1, of (z)dz. dr.
(22)
After discretizing in & (e.g. an Euler scheme), one repeatedly updates (9 and the boundary vector c(t) over time.
B.2.1 Efficient computation of (()
Computing (t)) directly at each time step requires evaluating M integrals, which can be computa- tionally intensive, especially when t changes incrementally and we need to update the matrix B()) repeatedly.
To overcome this inefficiency, we propose an approach that leverages the HiPPO framework to compute Btt) recursively as evolves. This method utilizes the properties of stationary kernels and the structure of the Legendre polynomials to enable efficient updates.
Leveraging Stationary Kernels Assuming that the kernel k(x,, t) is stationary, it depends only on the difference d |æc tl so k( t) = k(d). In our context, since we integrate over [tam t] with x I f, we have = x 2 Therefore, we can express k(x, t) as a function of d over the interval [0, lstan]:
(xx =)(d)) with dE [0, smar))
(23)
Our goal is to approximate k(d) over the interval [0, tstar] using the orthonormal Legendre basis functions scaled to this interval. Specifically, we can represent k(d) as
M Em (t))t)g! 9!00 (d),
k(d)
(24)
where ( (d) are the Legendre polynomials rescaled the interval [0, faan).
To efficiently compute the coefficients ém (t), we utilize the HiPPO-LegS framework, which provides a method for recursively updating the coefficients of a function projected onto orthogonal basis the interval expands. In our case, increases, the interval [fam t) over which k(d) is defined also expands, and we can update 7m (f) recursively.
Recursive Computation via HiPPO-LegS
Discretizing time with step size At and indexing to Esttu + kAt, the update rule using the Euler method is:
[
1 Bk((2), K
A
(25)
where d₁ = [Co(t)), (t4)..... CM [T and A E RMM and B E RM are again matrices defined by the HiPPO-LegS operator as Eq. 8 andd
Accounting for Variable Transformation and Parity The change of variables from to d= t- x introduces reflection in the function domain. Since the Legendre polynomials have definite parity, specifically,
(m((-) = (-1)" (2),
(26)
need to adjust the coefficients accordingly when considering the reflected function.
As result of this reflection, when projecting k(d) onto the Legendre basis, the coefficients 7. (1) computed via the HiPPO-LegS updates will correspond to a reflected version of the function. To
23

---

account for this, we apply a parity correction to the coefficients. Specifically, the corrected coefficients (t) are related (t) by sign change determined by the degree mi
(4) (-1)" (m))
(27)
This parity correction ensures that the computed coefficients properly represent the function over the interval [fstatt t] without the effect of the reflection.
By computing c(t) recursively as t evolves, we can efficiently update Bt) = c(t) (I...., at each time step without the need to evaluate the integrals directly. This approach significantly reduces the computational burden associated with updating (t) and allows for efficient computation K via the ODE.
B.2.2 Unstability of directly evolving Ki as ODE.
Empirically, we find that the direct ODE approach is less stable compared with RFF approach. Intuitively, it can be seen from the difference in the forms of their evolutions, especially in the first term. In RFF approach, the first term of the evolution of Fourier feature is of the form A (1) Z which includes evolving vectors with the operator C1 X + At) In direct ODE approach, the first term in the direct evolution of (5) is of the form (() Kin + K (1) which requires the Lyapunov operator L2 X - Att) + XA(t) The critical difference is that Cg has cigenvalues {} + } (where A and X, are eigenvalues of (t)) (Horn and Johnson, 1991]. while L has eigenvalues {Ar]. Since HiPPO-LegS uses a lower-triangular (?)) with negative diagonal entries, the eigenvalues all negative X < 0. Hence, the eigenvalues of the Lyapunov operator L2 approximately as twice negative as the eigenvalues of L), leading to a stiff ODE system with poorer numerical conditioning.
C Finite basis approximation of posterior OHSVGP
Here, we show that (u(()) is the distribution of HiPPO coefficients of the posterior OHSVGP (f).. From Eq,
2 the posterior of the function values evaluated at arbitrary indices X is (fx)
(fxx; K (t)
K iu m Kixfx
((t)
fxu K ((t))-
k(t)
S " ( ((()) ufx)
Ixu
in
Based on this, compute the mean of the m-th HiPPO coefficient for qe(f) as follows,
Eaih V (1x))(((
= Eais) [f() (x)
( (((1) KuuK
K(D- 1ol (z)dz
m!!)
(19(2')
/ ((zz,a))
d. (1'k(()-11
600 (x)dz m (1)
m
(28)
(1))
(2')
pom (x:)
f M(z,z')
dx dx
k(()-1mm m
oom (x))
(())
((00~mm)
uu
uu
24

---

which is exactly the variational mean of g(ul)). Similarly, the covariance between the 1-th and the m-th HiPPO coefficient for (f( can be computed as
Eyicn I// (w))))(a))t (f (x' ydd' )] m
mu
= =ʃ PaiS (xx))(('))) (1)))) (o' d'ddd'dd' [m() (m(]
leu = =/ / (M(2,x') K! K K S /us usz of (x) dd'd''
+ 1 Eaii)) [f())] Eace (x) (x' d')dddda m L
( J.u Kin-l (z)) (
[K] - s (ʃ
[(t)]
di (x' (x') [
K (t)
us de'
im
(
k K
K
K
(()
[ un
uu
uu
im
1,
a.
(((()
Jim
(29)
which is exactly the variational covariance between ui (1)
and in (u
(0)
Hence, if f ~ qe (f), then ) d f(x) (t) (æ)dæ, which implies that we can approximate the posterior OHSVGP with finite basis: ʃ = ΣM witl ( (xr), will ~ q(u ).
D Additional experimental details
D.1 Infectious disease modeling
For a sanity check, we also fit a weekly renewal-equation model Flaxman et al. 2020 Unwin et al. 2020 Monod al. 202T]. trained offline on the full history. This is well-established infectious disease model that has been widely utilized by scientists during the COVID-19 pandemic, and we include its results in section F.T (denoted as AR(2) Renewal in the legend). Interestingly, OHSVGP achieves better predictive performance than this traditional infectious disease model. Although this may be partly due to the strong inductive biases of the renewal equations, it nevertheless highlights OHSVGP's suitability for long infectious-diseane time series.
The details of renewal-equation model are follows. Let Yt.o be the number of deaths on day t = 1,..., '' in state a = 1,..., A (in our experiments, a = since we only fit on state), The probabilistic model is:
(expected deaths),
Pe, Σ my IFR(t -
(30)
(expected infections),
My,a
Ridd
Σ
mi, SI(T -i)
(31)
i=1
No -Cr. RT,
(Adjusted reproductive number),
padj
(32)
Na
Ría 3.3 x (T)))
(33)
(cumulative infections), Cr = NumDeathsInit,
(34)
where is the sigmoid function, N₂ is the population of country a, NumDeathslnit, is the initial number of deaths in country a, the IFR the infection fatality rate and SI the serial interval. The choice of (T) requires stochastic process and here choose to model it weekly random
25

---

effect day т using AR(2) process as in Monod et al. [2021], Unwin et al. (2020). Typically, Markov chain Monte Carlo (MCMC) methods are used infer the posterior distribution of fa and the other parameters, so it is less scalable than the online sparse GP based models shown in main text, especially when the number of states and time steps increase.
D.2 SVGPVAE model details
With SVGPVAE, we utilize Jazbec et al. [2021] and notation from Zhu et al. (2023), and we have the following encoder-decoder model:
utilizeJazbec
p(y1:)) = p(fer) II p(y)
(35)
with likelihood p(y, [f) (N(y))lp(f)) and decoder network & Re - Rdd, The encoder
s 2 2R yields (y)-2, v.2 o(y.). fr follows an L-dimensional multi-output GP and its approximate posterior given by:
L
(furr)
= II (qu'), ((u')=((ulmm) A'), 1=1
(36)
st Kuu + (iadiva(Vier)) Khu m' = Kyu (S') da
(37)
A' (S') "
(38)
where p(fi.T] ) the prior conditional distribution.
Following Jazbec et al. [2021], the objective function is defined as:
LSvGPVAE (0)
[Eyt log p(y:(f) log.V(ffl5.. v() Σ ch. 1=1
(39)
where C22 is the "Hensman" ELBO described in Equation 7 of Jazbec et al.
2021 ). Since the
variational parameters m' and S and the likelihood are all amortized by neural networks, we further add EWC (Elastic Weight Consolidation; (Kirkpatrick et al.) (017]) regularization for both encoder and decoder networks to the loss above for continual learning.
E Algorithmic Breakdown of OHSVGP
Algorithm The HIPPO-SVGP ELBO for single task of data. Differences with SVGP blue,
Require: x= [I1,.... ".1,...r,, F, = t1} (training time steps up to time t1).
[V1..... Bn) (training targets),
Z E IR Mxl (inducing points) A(t) E MRM B(t) PRMM Mxl (HIPPO matrices),
mu E Mx1 S,, E 18 M M (variational params)
1: K = (XX, Z), K = kZZ, Z). K& KG, from HIPPO ODEs evolved from 0 to the final time step with HIPPO matrices A(t), B(t)
2: p(zz)) K K, (Ktt))'mm,, mm
> Variational Posterior Mean
3: (x = -K Kh
Sal(K(i))) -k u/
Variational Posterior Variance
4: Evarexp Ei= (. log p(y: fi]
closed form quadrature/MC
5: KL t- KL(N(m,.. (.))
(t, K.,)))
6: return Evarexp KL.
26

---

Algorithm 2 The OHSVGP ELBO on the second task after learning the first task. Differences with OSVGP in blue.
Require:
X' {t1 < 21..... z = t2) (training time steps up to time t2),
(V'...., [.......V.... (training targets),
Zu E R Mxl (frozen and learned inducing points from task 1) with inducing variables
ur fZz,),
Zeg gMxx (new inducing points for task 2) with inducing variables ur f(Zi2).
A(t) E RMXM B(t) E R g (HIPPO matrices),
E ERAX1 Sun, E RRMX (learned variational params from task 1),
mui
Muzz
R Mx1
Sung ER MxM (new variational params for task 2),
K ue
1: Ke'nez = k(X', Z12), Ku."", (Z1 Ztz), Kin, evolved from 0 to the final time step /2,
Kla uuyuu, u., evolved from to the final time step t2 with HIPPO matrices A(t), B(t)
2: Mez (x,)) = K fluz
(K uiz
) Muaz
D Variational Posterior Mean
3: ( = =K)) P.st
KK
(Kis ury
)
K ur ura
S Su,,,(KK
K niy
Variational
f/ne
Posterior Variance
4: vvarexp + Σi=1 (ws). (ax)) logp(1/ f()]
closed form or quadrature/MC
5: m Kl ui, Uez Heg Uez Mus,
D Mean dts (u Pta (u (u MMg)) (Uga du
6: St.t. = K u. uqg 122
K's (Kig uzz )- "ungne + K 1222 K 102 Sue, (Ki, uz
K
b Covariance of qt, (uG,) = Ptg (ur, (uee )4t )dus,
7: KL t KL(N(mu,, Su.,) [l K. Ku.jue!! u.s
8:
TorrectionTerm(1n (2) + KL(N(mn Scta)llV Kú, un
KKLWN(m ( Su,,)
9: return Evarexp KL + CorrectionTerm((t.
Additional results
F1 Full results including RMSE and ECE
We also include the full results of NLPD and RMSE (ECE for experiments on COVID due to non-Gaussian likelihood) containing evaluation of Task after learning tasks j = i, + ... 10 (5 for experiments on COVID) for all i.
We define the Root Mean Squared Error (RMSE) and Negative Log Predictive Density (NLPD) as the following:
N (ys - -x)² NLPD = N
N log (ugi(ys wi).
RMSE
Σ
(40)
N
(41)
The Expected Calibration Error Guo et al. 2017 ; ECE) is defined as
N
ECE
1ly [4 1-cz (2x), á 1+0, (I:))) ck.
(42)
K
N
27

---

where
K = 10, Ch E {0.05, 0.15, ..., 0.95)
N is the number of test points (Ii Vi),
(x.) = the empirical p-quantile of the predictive samples { (1) J S = 100
1, if the argument true,
1() 0, otherwise.
Figure shows the RMSE results for time series regression experiments (results of NLPD are in Figure in the main text). Figure 8 and show NLPD and ECE results on COVID dataset with an additional sanity check baseline described Appendix D. Figure T0 show full results of RMSE and NLPD for UCI datasets, and Figure 14 and 15 show full results of RMSE and NLPD for ERA5 dataset.
6.30
0.10
0.00
0.31
(a) Solar, M=50
(b) Solar, M=150
(c) Audio, M=100
(d) Audio, M=200
Figure 7: Test set RMSE over the learned tasks vs. number of learned tasks for Solar Irradiance and Audio signal prediction dataset.
2R Renewal
OSVGP
OVC
OVC-optZ
OHSVGP
3.55
322
51
3.0
1.45
(a) M=15
AR(2) Renewal
OSVGP
OVC
OVC.
OHSVGP
Task
lask
Tack
Task
Taak
4.2
3.5
3.54
3.2
32
3.52
222
2.7
3.5
2.2
T55K3
Tayk3
T5515
Tisks
Th51s
(b) M=30
Figure 8: Test set NLPD per task after continually learning each task for all the 5 tasks on COVID dataset.
28

---

AR(2) Renewal
OSVGP
OVC
OVC-optZ
OHSVGP
Tusk
Task
0335
6:22
0152
0:35
0.75
0:00
0.1
0.03
0.00
0.352
0335
Tasks
Ticks
Tasks
Dakk
T2cks
(a) M=15
AR(2) Renewal
OSVGP
ovc
GVC-aptZ
OHSVGP
Task
0.32
0115
6.25
0.150
5.22
6.90
0.15
0.1
8.30
6:375
1.08
0.352
0.04
5.05
Tesls
Teaks
Tasts
Tenlss
(b) M=30
Figure 9: Test set ECE per task after continually learning each task for all the 5 tasks on COVID dataset.
29

---

Tak
Taee2
Tase3
Taxt
0.35
14
0.5
13
0.75
12
0.65
0.55
30
18466
1.55
1.30
0.05
12
0.35
1.5
0.35
0.70
0.30
0.75
0.70
Tases
Tasks
Tasks
Task
0.35
0.750
0.90
0.755
0.35
0.30
0875
OSVGP
0.75
0.690
OVC
OHSVGP-K
0.85
0.000
OHSVGP-0
0.755
(a) Skilleraft (1st dimension)
Desc
Tack
0.75
0.70
20932
0.85
0.50
0.7
0.6
0.5
0.50
baks
0.45
0.35
0.60
0.55
0.55
0.75
1.05
0.59
0.45
0.59
0.35
Tadks
Tosks
365
1.075
1.925
OSVGP
0.975
0.223
OHSVGP-a
0.575
Tasks
(b) Skillcraft (L2)
Figure 10: Test set RMSE per task after continually learning each task for all the 10 tasks on UCI Skillcraft dataset.
30

---

Task
Task 2
Task
Task
18
19
MF7
1.7
1.6
14
1.2
1.2
1.0
1.0
1.4
1.3
Tasch
Tascs
Tascs
Taaca
1.0
1.9
1.7
2.0
1.7
15
1.6
1.8
1.4
1.5
12
1.
1.3
15
1.0
0.9
1.0
T655 10
1.6
1.4
OSVGP
1.2
OVC
OHSVGP-K
OHSVGP-0
0.9
(a) Skillcraft (1st dimension)
Task
Task
Task
Task
1.8
2.0
1.6
1.4
1.2
3.0
1.0
0.0
1.0
0.0
3.0
1.9
1.9
1.7
2.0
1.2
1.4
1.5
1.2
3.0
1.3
1.0
12
1.0
Tack
Task:
1.6
3.7
1.44
3.5
1.5
1.0
OVC
OSVGP
OHSVGP-K
1.4
OHSVGP-0
1.3
(b) Skillcraft (1.2)
Figure 11: Test set NLPD per task after continually learning each task for all the 10 tasks on UCI Skillcraft dataset.
31

---

Teat
Tasik
0.30
0.850
0.5
1992
2.400
2.350
0.30
0.4
3.325
2.300
1adk
1335
0.20
5.27
0.25
0.26
0.27
0.26
0.24
0.24
0.25
0.24
0.22
0.22
0.22
0.23
0.22
Tasks
Tasks
0.31
0.285
0.30
0.255
OSVGP
0.230
0.245
0.28
OHSVGP-k
0.235
0.27
OHSVGP-o
0.225
(a) Powerplant (Ist dimension)
Taak
Tasc
Taat
Task
0.26
0.240
0.250
2.22
0.22
2.225
0.240
0.23
To.k55
bac
Tou6
0.33
0.270
0.32
0.263
0.235
0.253
0.230
0.250
2.260
0.29
0.225
0.28
0.245
0.250
0.25
0.210
Taxkc
0.275
2305
0.265
2300
OSVGP
260
2275
OVC
2250
OHSVGP-k
0.55
2225
OHSVGP-o
0.250
2200
2175
125
135
(b) Powerplant (1.2)
Figure 12: Test set RMSE per task after continually learning each task for all the 10 tasks on UCI Powerplant dataset.
32

---

Task
Task3
1.4
0.4
3.0
1.2
3.5
2.0
3.5
0.3
0.2
0.5
0.2
0.3
06.3
Tenos
Tocks
Junce
Toace
0.150
0.25
0.15
0.125
0.10
0.075
0.23
0.05
0.023
0.08
0.00
000
3.023
0.05
0.00
-000
Ticl5
basks
bsss
Tosks
0.225
0.10
0.250
0.225
MA9
0.02
0.150
OHSVGP-K
0.00
0.125
OHSVGP-0
0.100
124
(a) Powerplant (Ist dimension)
Task: 1
0.03
Task
Task33
0.57
Thsk
0.02
0.10
0.0
0.04
0.00
0.86
0.54
0.01
0.02
0.02
0.02
9.00
0.64
0.04
-0.00
0.08
0.06
Thsk 5
Tasc
Task 7
0.12
0.09
0.10
0.00
0.88
0.20
-0.4
0.02
0.05
0.00
0.03
0.10
Tasks
15505
T3545
-0.01
-0.02
OSVGP
0.00
DVC
DHSVGP-K
OHSVGP-0
0.04
0.06
0.07
13535
(b) Powerplant (1.2)
Figure 13: Test set NLPD per task after continually learning each task for all the 10 tasks on UCI Powerplant dataset.
33

---

Tesk
1.1
0.9
1.0
10
0.7
0.7
0.3
0.2
0.4
0.2
0.2
Tencs
1.1
1.4
1.0
1.0
1.2
0.9
12
1.0
0.0
0.5
0.6
0.4
0.3
0.3
bash 30
1.1
OSVGP
0.7
OVC
OVC-per-dlim
0.6
0.6
0.5
0.5
OHSVGP
0.4
0.4
(a) M-50
Tait
Tast
0.85
Task:44
0.7
0.6
0.7
0.55
0.5
0.45
0.5
0.35
0.3
0.25
0.20
0.1
Taxcs
Taxcs
TaxEs
0.05
0.70
0.6
0.0
0.55
0.6
0.50
0.40
0.45
0.4
0.40
0.3
0.30
0.35
bah.10
0.41
0.500
0.475
0.430
OSVGP
OVC
0.30
OVC-per-dlm
0.38
0.350
OHSVGP
0.24
0.325
(b) M=100
Figure 14: Test set RMSE per task after continually learning each task for all the 10 tasks ERA5 dataset.
34

---

Task
Taak
Task
NUD
Thrles
Tayls
11-kk5
70506
Dy5k7
Thck
22.5
1.0
12
17.5
12.5
10
5.0
0.0
Tasks
Tisks
13
Tock 10
14
OSVGP
OVC
OVC-per-dim
OHSVGP
Tisks
3,555
(a) M=50
2
Tolk
Teak 4
2,5
2.
NUD
1.5
0.5
0.3
0.0
Tal
18
4.3
2.5
3.0
1.5
15
0.5
15
055
2.5
1.5
1.0
1.4
OSVGP
1.2
1.4
OVC
1.2
OVC-per-dim
1.0
OHSVGP
0.8
Tascs
(b) M=100
Figure 15: Test set NLPD per task after continually learning each task for all the 10 tasks on ERA5 dataset.
35

---

F.2 Results for time series regression with trainable kernel hyperparameters
Figure 16 and show RMSE and NLPD results for time series regression experiments based on trainable kernel hyperparameters (i.e.. keep optimizing kernel hyperparameters online in all the tasks). Notice that OVC only compatible with fixed kemel [Maddox al., 202T]. so we don't consider here. Compared with the results based on fixed kernel in the main text, here all methods show less stable performance. Previous works either find a well-performed fixed kernel [Maddox et al. 2021]
or scale the KI. terms in the online ELBO objective with positive factor requiring careful funing to
mitigate the unstable online optimization of kernel hyperparameters (Stanton al. 2021 Kapoor
al. 2021
oscFs
0.75
ove
100
0.00
0.-
(a) Solar, M=50
(b) Solar, M=150
Audio, M=100
(d) Audio, M=200
Figure 16: Test set RMSE over the learned tasks vs. number of learned tasks for Solar Irradiance and Audio signal prediction dataset (keep updating kernel hyperparameters).
0507A Deff 0-60P
17
-15
38
(a) Solar, M=50
(b) Solar, M=150
(c) Audio, M=100
(d) Audio, M=200
Figure 17: Test set NLPD over the learned tasks vs. number of learned tasks for Solar Irradiance and Audio signal prediction dataset (keep updating kernel hyperparameters).
F.3
Visualization of impacts of sorting criterion for OHSVGP in continual learning
Here, we consider fitting OHSVGPs with 20 inducing variables for continual binary classification problem the Two-moon dataset [Ganin et al. 20166. The data is splitted into three task and we use a Bernoulli likelihood to model binary labels. We consider three different sorting criteria:
Random, denoted as OHSVGP-rand. The order of data points in each task is obtained via random permutation.
Kernel similarity maximization, denoted OHSVGP-k max. We select the i-th point in task j to be x (g) = arg max k( 221 for i 1, and the first point in first task
is set to be x0 = arg max k(x; 0). The intuition is that the signals to memorize, when computing the prior covariance matrices, tend be more smooth if the consecutive x's are close to each other.
Kernel similarity minimization, denoted as OHSVGP-k-min. We select the i-th point in task j to be 2 arg min k(x x (1) for > 1, and the first point in first task is (1) = arg minge x k(æ; 0). In this case, we deliberately make it difficult to memorize the signals the recurrent computation for the prior covariance matrices.
(j)
set to be xi
Figure show the decision boundaries after each task for different OHSVGPs based on different sorting entteri. We also include the decision boundaries of an OSVGP model for reference. Both OHSVGP-k-man and OHSVGP-ram return decision boundaries achieving 100% accuracy, while OHSVGP-k-min show catastrophic forgetting after Task 3, which suggests OHSVGP requires sensible sorting criterion to perform well in continual learning tasks.
36

---

100
10
100
0.5
05
0.5
0.0
0.0
0.0
-0.5
-0.5
-0.5
(a) OSVGP, after Task
(b) OSVGP, after Task 2
(c) OSVGP, after Task 3
1.0
10:
0.5
0.5
0.0
0.0
0.0
-0.5
-0.
-0.
2
(d) OHSVGP-rand, after Task
(e) DHSVGP-rand, after Task 2
(I) OHSVGP-rand, after Task 3
0.5
0.0
0.0
0.5
0.5
0.5
(g) CHSVGP--kmaat after Task
(h) OHSVGP- k-rxx,---ma,,,, after Task 2
(i) OHSVGH --op--k-ma m-xxx after Task
1.0
1.0
(i) OHSVGP-k-min, after Task
(k) OHSVGP-k-min, after Task 2 (I) OHSVGP-K-min, after Task 3
Figure 18: Decision boundaries of OSVGP, and OHSVGP models with different sorting criteria after each task (3 in total) on the Two-moon dataset. For OSVGP, we visualize the inducing points with red color.
37

---

F.4 Comparison of meassree variants
Figure shows the results OHSGPR applied to a toy time-series regression dataset, where the data is split chronologically into three equal segments, used Tasks 1-3 in an online learning setup. The figure compares the effect of several variants of the HiPPO operators [Gu al., 2020 2023] when used for OHSGPR. Subfigures (a-c) correspond to HiPPO-LegS as used in all of our main experiments. Subfigures (d-f) apply HiPPO-LegT, (g-i) apply HiPPO-L.gT,, based on the Laguerre polynomial basis, and (j-l) apply HiPPO-Fou'T based on Fourier basis functions. While OHSGPR-LegS successfully memorizes all the past tasks, OHSGPR-LegT, OHSGPR-LagT and OHSGPR-Foul I all demonstrate catastrophic forgetting to certain degree since instead of the uniform measure over the past (as used in HiPPO-LegS), they are based on measures which place more mass over the recent history. LegT and FouT use fixed-length sliding window measure, while LagT uses exponentially decaying measure, which assigns more importance to recent history.
0-
(a) OHSGPR-LegS, after Task
(b) OHSGPR-LegS after Task
(c) OHSGPR-Legpe after Task 3
M
(d) OHSGPR LegT, after Task 1
(e) OHSGPR-LegT, after Task
(1 OHSGPR LegT, after Task
(g) OHSGPR LagT, after Task 1
(h) OHSGPR LagT, after Task
OHSGPR LagT, after Task
(i) OHSGPR-FouT, after Task
(k) OHSGPR-FouT, after Task
(I) OHSGPR-FouTT after Task
Figure 19: Comparison of OHSGPR based on different HiPPO variants on toy online regression dataset.
38


Reviews:
[{"review_id": "RGIUf5qD8o", "reviewer": "Reviewer_CC7b", "summary": "The paper proposes a new GP model based on interdomain GP and HiPPO. Inspired by the usage of HiPPO in approximating time series in RNN, the authors adopt the bases in HiPPO for online GP. They show how the corresponding kernel matrices can be updated with differential equations (which is later approximated by difference equations). In the experiments, the proposed algorithm is evaluated in several problems (including online prediction, continual learning, and sparse GP VAE).", "questions": "1. Can you formally explain why using HiPPO can address the forgetting issue?\n\n2. What is the benefit of using time-varying bases here, when compared with stationary bases?\n\n3. In sec 3, in the integral with basis functions (e.g., line 133, line 146, etc.) what's the domain the integral is defined? Previously in in sec 2.5, the integral is from `-\\inf` to `t`, is it also the case here? or from `-\\inf` to `\\inf`?\n\n4. Can you explain why the matrix ODE version leads to numerical instability?\n\n5. From the current paper, I still don't fully follow how the mult-dimensional input case is derived. How is the integral originally defined on the real axis extended to here?", "limitations": "No. There's no explicit sections or paragraphs on limitation, although some limitations are mentioned in the writing. I suggest to make them more obvious.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths\n\n1. The combination of GP and HiPPO to my understanding is novel \n2. The overall structure and the logic is clear (though some details are missing see Questions).\n3. The authors conducted diverse experiments to validate the proposed approach. \n\nWeakness. \n\n1. It is unclear from the theoretical side the benefits of using basis functions in HiPPO in the online setting. The current paper merely draws an analogy to the time series approximation in RNN. \n\n2. The paper is motivated by updating inducing points may lead to forgetting. But it doesn't give a formal theoretical argument why using the HiPPO bases can avoid this issue.\n\n2. It is unclear whether the proposed algorithm is scalable to high-dim, complex problems with large data. The experimental results are mainly low-dimensional. However, I acknowledge this is a quite prevalent limitation of GPs in general.", "was_revised": true, "final_justification": "The rebuttal addressed my concern especially about the adaptivity.", "rebuttal": "### Global Response:\n\nWe sincerely thank all reviewers for the constructive feedbacks and the overall positive reviews. We are encouraged that reviewers appreciated our **novel (CC7b, GFdF, Aht3)** and **technically sound (D83n, GFdF)** method that addresses the **critical (D83n)** catastrophic forgetting problem in online GPs. In general, the reviewers are happy with the **clarity (D83n, Aht3)** of the manuscript and our **comprehensive (CC7b, HjVk) experimental evaluation** of OHSVGP, which shows **strong performance (D83n, HjVk, GFdF, Aht3)** compared with baselines.\n\n---\n### Individual Response to Reviewer CC7b:\n\nBelow we address the questions from Reviewer CC7b.\n\n> Q1. Benefit of HiPPO basis in online setting and why it can address the forgetting issue.\n\nOne of the key benefits of HiPPO basis is its adaptive nature, and to our knowledge, OHSVGP is the first work that considers adaptive basis for interdomain GP. Our motivation for this adaptive basis is exactly trying to address the forgetting issue: the basis function solely focus on history, and we discuss it in details in our response to Q3. \n\nMoreover, in conventional OSVGP, when new task data arrives, updating the inducing points by optimizing online ELBO cannot gaurantee there are sufficient inducing points left in the previous task regions as shown in our experiments, which makes the prediction in these past regions close to the uninformative prior since each inducing point can only capture local behavior of GP approximation. In contrast, OHSVGP defines inducing variables as time-dependent orthogonal polynomial basis functions which by design can capture the global functional behavior in all past data regions. In particular, the basis of HiPPO-LegS is constructed based on uniform measure over the whole history which imposes strong inductive bias in the GP approximation to ensure the posterior predictive is equally well in all the past regions (see Appendix E.4 for the comparison between OHSVGP based on HiPPO-LegS and other HiPPO variants). \n\nFrom a theoretical persepective, given regularity assumption of the function to memorize, its approximation based on orthogonal polynomials basis in HiPPO frameworks are theoretically shown to be asymptotically exact in L2 sense with $M=\\infty$ (Proposition 6 in [1]). Although it is a result in the limit, this asymptotically perfect memory property does not hold for any arbitrary basis functions (i.e., the memory constructed based on many other basis functions are gauranteed to be lossy even in the limit). Empirically, HiPPO based RNNs and the following SSMs (S4, Mamba) inspired by it have shown profound success in sequence modeling tasks requiring long-term memory. \n\n[1] \"HiPPO: Recurrent Memory with Optimal Polynomial Projections\", NeurIPS 2020\n\n> Q2. Scalability to high-dimensional problems.\n\nWe conduct the experiments with relevant baselines on commonly used datasets in the context of GPs. Based on the overall positive results, we believe our method is one of the state-of-the-art GP based methods in online/continual learning. Indeed, GPs may have issues when scaling up to high-dimensional inputs. This is an interesting and open research question of generic GP or kernel methods orthogonal to the contribution of this work, which is about building reliable long-term memory.\n\n> Q3. Benefit of time-varing basis against stationary ones.\n\nThe \"time-varying\" nature of HiPPO bases is a principled way to handle the expanding time domain in online learning. Specifically, at time $t$, the orthogonal polynomial basis functions are defined on $[0,t]$. As new data arrives, this interval expands to $[0,t+Δt]$ to cover the new time region where the new data live, and the basis functions adapt accordingly. This is fundamentally different from stationary bases on a fixed interval. The key benefits are:\n* Adaptive coverage: The bases always span exactly the observed history, no more, no less. An unnecessarily large coverage including regions where no data exists may lead to underfitting in the historical region (see the comparison between OVFF and OHSVGP in Sec 5.1), and a small coverage excluding some regions with data will make the prediction in those regions close to the uninformative prior. The issue of suboptimal coverage has also been discussed in previous works such as [2].\n* No need for pre-specification: Unlike stationary bases, we don't need to know the final time horizon covering all the future tasks in advance, which is in general impossible for real world online learning problems.\n\n[2] \"Variational Fourier Features for Gaussian Processe\", JMLR 2018\n\n\n> Q4. Domain of the integral in Sec 3.\n\nThe domain of the integral for generic HiPPO is from $-\\infty$ to $t$ as in Sec 2.5. For HiPPO-LegS, since the measure it uses is uniform over [0, t], $\\omega^{(t)}(x)=\\frac{1}{t}\\mathbb{1}_{[0,t]}(x)$, the domain of the integral essentially becomes from $0$ to $t$ by noticing that the time varing basis is defined as $\\phi^{(t)}(x)=g^{(t)}(x) \\omega^{(t)}(x)$.\n\n> Q5. Why the matrix ODE version of $K_{uu}^{(t)}$ leads to numerical instability.\n\nThe instability of the direct ODE approach compared with RFF approach can be seen from the difference in the forms of their evolutions (especially the first term):\n\n* RFF approach: The evolution of Fourier fearure is of the form $\\frac{d}{d t} \\mathbf{Z}_w^{(t)}=A(t) \\mathbf{Z}_w^{(t)}+\\cdots$, which includes evolving vectors with the operator $\\mathcal{L}_1: X \\mapsto A(t) X$.\n\n* Direct ODE approach: The direct evolution of $K_{uu}^{(t)}$ is of the form $\\frac{d}{dt} K_{uu}^{(t)} = [A(t) K_{u u}^{(t)}+K_{u u}^{(t)} A(t)^{\\top}] + \\cdots$, which requires the Lyapunov operator $\\mathcal{L}_2: X \\mapsto A(t) X+X A(t)^{\\top}$.\n\nThe critical difference is that $\\mathcal{L}_2$ has eigenvalues $\\lambda_i+\\lambda_j$ (where $\\lambda_i, \\lambda_j$ are eigenvalues of $A(t)$), while $\\mathcal{L}_1$ has eigenvalues $\\lambda_i$. Since HiPPO-LegS uses a lower-triangular $A(t)$ with negative diagonal entries, the eigenvalues $\\lambda_i<0$. The Lyapunov operator can thus have eigenvalues $\\lambda_i+\\lambda_j$ that are approximately twice as negative, leading to a stiff ODE system with poorer numerical conditioning.\n\nWhile we find simple numerical methods, such as Euler method may be unstable to solve this ODE, we will explore fast and numerically stable solvers for this matrix ODE in future work.\n\n> Q6. How is the mult-dimensional input case is derived\n\nGiven an ordered training batch $\\\\{x_1, \\cdots, x_N\\\\}$, Eq. 5 can be viewed as a discretization (with step size $\\Delta t$) of an ODE solving a path integral $\\int k(x_n, x(s)) \\phi^{(t)}(s) ds$. The $i$-th training input $x_i$ is assumed to be $x_i:=x(i \\Delta t)$ and thus the path integral is approximately solved with discretized recurrence based on the training inputs corresponding to $\\\\{x(i\\Delta t)\\\\}_{i=1}^N$. Notice that different orders of the training batch implicitly define different path integrals to approximate. While they all define valid interdomain inducing variables, the order does matter for the performance of OHSVGP in multi-dimensional input continual learning as demonstrated in our UCI experiments and the 2D-input visualization example in Appendix E.3."}, {"review_id": "b4fD3imB2C", "reviewer": "Reviewer_D83n", "summary": "The paper proposes Online HiPPO Sparse Variational Gaussian Process (OHSVGP), a novel method for online Gaussian process that preserves long-term memory through recurrent updates. The key idea is to interpret the HiPPO framework—originally developed for state-space models—as defining interdomain inducing variables via time-varying orthogonal projections. This construction allows the model to capture history compactly and update kernel matrices via ODEs, avoiding expensive optimization over inducing points. The method is applied to time series prediction, continual learning, showing better accuracy and computational efficiency compared to existing approaches.", "questions": "* Do the authors expect this approach to work for non-stationary kernels beyond the RFF approximation used here? Would something like non-stationary Fourier features (e.g., [43]) be feasible?\n* Kernel hyperparameters are fixed after the first task in most experiments. While this is understandable for stability, it limits the model’s flexibility in adapting to data shifts. Is there a principled way to integrate kernel hyperparameter updates into the recurrence framework? \n* Can you comment on the sensitivity of the method to the number of inducing variables $M$?", "limitations": "Nothing major to flag", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "This is a well-written and technically sound paper addressing a critical issue in online GP learning. The integration of HiPPO into interdomain SVGPs is both elegant and well-motivated, and the empirical results are strong. The ODE-based kernel updates combined with RFF approximations are especially good from a scalability standpoint.\n\nOverall, the method is original, effective, and cleanly presented.\n\nSome comments:\n\n* The idea of viewing the HiPPO projection coefficients as time-evolving inducing variables is elegant. However, it assumes a time-ordered input sequence. In Section 3.3 (and 5.2), the authors mention a pseudo-time ordering heuristic for unordered data. It would be good to clarify how robust the method is to different sorting heuristics, including the overhead for sorting.\n* The HiPPO basis function evolution introduces a dependency on accurate numerical integration over time. The recurrence is discretized (Euler or bilinear), but no error analysis or stability bounds are provided. This is important since the quality of the memory encoding depends on these dynamics.\n* In Figure 6, increasing the number of inducing variables from $M = 50$ to $M = 100$ yields no clear performance improvement and occasionally degrades performance. This contradicts the common intuition that more inducing variables should improve approximation quality. An ablation or sensitivity analysis would help clarify the cause.", "was_revised": true, "final_justification": "I maintain a positive overall opinion. The rebuttal addressed the main concerns I raised. The main unresolved aspect is the lack of formal stability/error analysis for the HiPPO recurrence discretization, but given the strong empirical results and sound methodological design I confirm the borderline-accept score.", "rebuttal": "### Global Response:\n\nWe sincerely thank all reviewers for the constructive feedbacks and the overall positive reviews. We are encouraged that reviewers appreciated our **novel (CC7b, GFdF, Aht3)** and **technically sound (D83n, GFdF)** method that addresses the **critical (D83n)** catastrophic forgetting problem in online GPs. In general, the reviewers are happy with the **clarity (D83n, Aht3)** of the manuscript and our **comprehensive (CC7b, HjVk) experimental evaluation** of OHSVGP, which shows **strong performance (D83n, HjVk, GFdF, Aht3)** compared with baselines.\n\n---\n\n### Individual Response to Reviewer D83n:\n\nBelow we address the questions from Reviewer D83n.\n\n> Q1. How robust OHSVGP is to sorting methods and its computational overhead.\n\nThe sorting is conducted only once every time when a new task arrives and its cost is negligible compared with the training loops. The performance of OHSVGP does depend on the sorting method and we include a 2D-input binary continual classification example in Appendix E.3 for visualization of the impacts of different sorting methods (there is a typo in the manuscript - the results of OHSVGP-k-max and OHSVGP-k-min there should be swapped. See our response to Q2 from Reviewer Aht3 for details). In that example, we found that OHSVGP is robust as long as we do not use an advesarially constructed sorting strategy. In our experiments, we found OHSVGP-k, which heuristically sorts the data points based on the kernel distances among them, indeed achieves worse performance, than OHSVGP-o, based on oracle ordering. However, it is still decently robust in the sense that it achieves competitive or better performance than the OSVGP baseline as shown in Figure 5.\n\n> Q2. Discretization error of HiPPO recurrence.\n\nWhile we do not have a formal stability bound yet, empirically we observe that the specification of the discretization step size depends on how smooth the time series is (e.g., for smooth time series (implying larger length-scale), the step size can be relataively large). In our experiments, we found that setting the step size to be within $2 \\times \\text{length-scale}$ is able to gaurantee reliable performance. \n\n> Q3. In Figure 6, increasing $M$ does not improve performance.\n\nIncreasing $M$ from 50 to 100 does improve the performance in Figure 6. **The scales of the y-axis in Figure 6(a) and Figure 6(b) are different.** We choose to maintain different scales in y-axis because otherwise the results in Figure 6(b) will be too squeezy.\n\n> Q4. Do the authors expect this approach to work for non-stationary kernels.\n\nYes, as long as there is an available feature approximation for a chosen kernel (e.g., [1] provides an example of feature approximation for non-stationary kernel).\n\n[1] \"Spatial Mapping with Gaussian Processes and Nonstationary Fourier Features\", Journal of Spatial Statistics 2018.\n\n> Q5. Fixed kernel hyperparamters may limit adaptation.\n\nThe instability of online kernel hyperparameters update is a common problem in online SVGPs. Previous works either fix the kernel hyperparameters to ensure stable performance [2] or rely on additional tricks, such as generalized variational inference [3,4] or replay buffer that keeps a subset of past data for retraining to make the online variational inference closer to batch variational inference [5]. Here, we adopt the simpler strategy of fixing kernel hyperparameters since it already gives us good results for the benchmarks considered in this work and the long-term memory capability of OHSVGP, which is our main contribution, has been sufficiently demonstrated. Moreover, the aforementioned tricks of mitigating the unstable online kernel hyperparameters updates can all be easily integrated into OHSVGP since they are completely orthogonal and compatible to our key contribution.\n\n[2] \"Conditioning Sparse Variational Gaussian Processes for Online Decision-making\", NeurIPS 2021\n\n[3] \"Kernel Interpolation for Scalable Online Gaussian Processes\", AISTATS 2021\n\n[4] \"Variational Auto-regressive Gaussian Processes for Continual Learning\", ICML 2021\n\n[5] \"Memory-based Dual Gaussian Processes for Sequential Learning\", ICML 2023\n\n> Q6. Sensitivity of the method to the number of inducing variables $M$.\n\nSimilar to standard SVGP, typically performance improves when $M$ increases (possibly with diminishing return though). In terms of memory preservation, we found that catastrophic forgetting in OSVGP tends to be more severe when $M$ is small (e.g., Figure 2(a) vs Figure 2(b)). This is possibly because when $M$ is small, the location of each inducing point needs to be representative, while when $M$ is relatively large, the inducing points are more likely to sufficiently cover the historical regions even though some of them are placed at suboptimal locations."}, {"review_id": "GxMyxIWObX", "reviewer": "Reviewer_HjVk", "summary": "This paper incorporates the Hippo framework within the interdomain GP for developing an online novel GP method that can memory the history information.", "questions": "Some minor comments: \n\nIn the background section, the introduction of GPVAE appears somewhat unnecessary, as it is not essential to the main narrative. Moreover, the summary at the beginning of the section does not mention GPVAE, which makes its inclusion feel somewhat disconnected.", "limitations": "See the weaknesses and questions part.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. Establishes a connection between the Hippo framework and inducing points from the perspective of interdomain GPs.\n\n2. Leverages RFF to decouple the two inputs, enabling the ODE-based evolution of Kuu.\n\n3. Provides comprehensive simulations that demonstrate not only the model’s superior performance on various tasks but also its computational efficiency.\n\n4. The paper is clearly written and easy to follow.\n\nWeaknesses: \nOverall, the paper is well-executed. While the connection between the Hippo framework and interdomain GP is clearly presented, it may be seen as a natural extension rather than a fundamentally novel contribution, as it is quite straightforward.", "was_revised": true, "final_justification": "The author addressed my concern and I will keep my positive score.", "rebuttal": "### Global Response:\n\nWe sincerely thank all reviewers for the constructive feedbacks and the overall positive reviews. We are encouraged that reviewers appreciated our **novel (CC7b, GFdF, Aht3)** and **technically sound (D83n, GFdF)** method that addresses the **critical (D83n)** catastrophic forgetting problem in online GPs. In general, the reviewers are happy with the **clarity (D83n, Aht3)** of the manuscript and our **comprehensive (CC7b, HjVk) experimental evaluation** of OHSVGP, which shows **strong performance (D83n, HjVk, GFdF, Aht3)** compared with baselines.\n\n---\n\n### Individual Response to Reviewer HjVk:\n\nBelow we address the questions from Reviewer HjVk.\n\n> Q1. Clarification of the novelty and contribution.\n\nOHSVGP is indeed built on two well-established frameworks, HiPPO and online SVGP. However, the identification of the connection between the two frameworks from two different domains (deep sequence modeling and probabilistic machine learning) is non-trivial and requires in-depth understanding of both areas. The novelties and contributions of many interdomain GP methods are about finding good basis functions that enable desired properties in the approximate GP (e.g., reduced computational complexity [1,2]). In our case, we leverage the long-term memory preservation property of HiPPO framework to build an online GP free of severe catastrophic forgetting. To our knowledge, this is the first work that considers interdomain GP in online setup and OHSVGP is the first interdomain GP method that considers time-varing basis function which keeps memorizing the whole history and also adapting to the new task. As a result, it is naturally suitble for online learning.\n\n[1] \"Variational Fourier Features for Gaussian Processe\", JMLR 2018\n\n[2] \"Sparse Gaussian Processes with Spherical Harmonic Features\", ICML 2020\n\n> Q2. The introduction of GPVAE in the Background appears somewhat unnecessary.\n\nWe consider continual learning in GPVAE to show that OHSVGP can be scaled up to higher dimensional problems, which will be of interest to the wider machine learning community. Hence, we also give a short description of GPVAE in the Background section and include its detailed model specification in Appendix D.2 for readers that are not familiar with it. It is indeed not a core part of our algorithm, but rather an add-on. We will mention GPVAE in the summary at the beginning of the Background section and add a pointer to Appendix D.2 in future revision of the manuscript."}, {"review_id": "MavRAvHQVv", "reviewer": "Reviewer_GFdF", "summary": "This paper proposes OHSVGP, an online Gaussian Process model that integrates the HiPPO framework into sparse variational GPs (SVGP) to address forgetting in continual learning. This allows recurrent online updates of kernel matrices via HiPPO's ODE structure, enabling long-term memory compression. Results show good predictive performance, memory retention, and computational efficiency compared to OSVGP, OVC, and OVFF baselines.", "questions": "What is the convergence behavior and computational cost of using RFF to approximate Kuu w.r.t. the number of RFF samples? \nFurthermore, have you conducted a time complexity analysis to facilitate a general comparison with other baseline models?\nSince your work incorporates both RFF-based and inducing point-based approximations, have you considered comparing it with RFF-based online Gaussian Processes (SSM)?", "limitations": "Yes, and see the weakness.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths: Good work to bridge HiPPO (SSM-based memory) with interdomain GPs for online learning. The reformulation of HiPPO projections as time-dependent inducing variables is novel and theoretically grounded. \n\nWeaknesses: Performance drops significantly with heuristic sorting (OHSVGP-k vs. OHSVGP-o). This reliance on \"oracle\" task ordering limits applicability to real-world continual learning.", "was_revised": false, "final_justification": "", "rebuttal": "### Global Response:\n\nWe sincerely thank all reviewers for the constructive feedbacks and the overall positive reviews. We are encouraged that reviewers appreciated our **novel (CC7b, GFdF, Aht3)** and **technically sound (D83n, GFdF)** method that addresses the **critical (D83n)** catastrophic forgetting problem in online GPs. In general, the reviewers are happy with the **clarity (D83n, Aht3)** of the manuscript and our **comprehensive (CC7b, HjVk) experimental evaluation** of OHSVGP, which shows **strong performance (D83n, HjVk, GFdF, Aht3)** compared with baselines.\n\n---\n\n### Individual Response to Reviewer GFdF:\n\nBelow we address the questions from Reviewer GFdF.\n\n> Q1. Performance drops with OHSVGP-k (heuristic sorting) in UCI experiments.\n\nFor continual learning problems where the order of data points within each task is unknown, the sorting method will impact the performance of OHSVGP and what sorting method may be preferable remains an open question. OHSVGP-k, which sorts the data points based on the kernel distances among them, indeed achieves worse performance, than OHSVGP-o, based on oracle ordering. However, OHSVGP-k is still decently robust in the sense that it achieves competitive or better results than the OSVGP baseline as shown in Figure 5.\n\nThis sorting step also can not be avoided in deep SSMs, such as Mamba (advanced extension of HiPPO), when applied to data modalities without pre-defined ordering (e.g., Vision [1], where Mamba approach requires the user to specify an ordering of the patches split from an image). Empirically, some heuristic sorting strategies can still enable good performance of deep SSMs on these unordered data.\n\n[1] \"Vision Mamba: Efficient visual representation learning with bidirectional state space model\", ICML 2024\n\n> Q2. RFF approximation quality and computational cost.\n\nWe expect there is a tradeoff between computational cost due to RFF sample size and the performance. However, with the current number of RFF samples considered in our experiments, we are already able to both obtain better results compared with other baselines and also reduce the training time. \n\nOHSVGP and OVC are almost equally fast and both of them are significantly more efficient than OSVGP. We report the time metric for time series regression experiments in Table 1. For continual learning on UCI datasets, the training time of OHSVGP and OVC are about half compared with OSVGP. The time complexity of OHSVGP and OVC are the same after obtaining the kernel matrices: $O(M^3 + N^2M)$, where M is the number of inducing points and N is the number of data points in the current task. Both OVC and OHSVGP update kernel matrices just once before entering the training loops every time when a new task arrives. OHSVGP further requires an additional step to sort the data points in the task for continual learning. However, all these one-time costs are negligible compared with the training loops for the variational parameters, and inducing locations (for OSVGP only).\n\nAll the baselines considered in the experiments are not based on RFF but based on exact kernel matrices, so even RFF may introduce some approximation bias, the benefit of OHSVGP in terms of long-term memory completely outweights it. We can turn the baselines into RFF-based online Gaussian processes by also approximating kernel evaluations with RFF, but we believe it is not likely to improve their performances. We are not aware of online GP methods tailored to RFF approximations and we would be interested if you can provide some references."}, {"review_id": "Yw0a2PfYCl", "reviewer": "Reviewer_Aht3", "summary": "The authors propose the OHSVGP model, which marries the HiPPO and the interdomain GP frameworks to allow online learning with long-range memory support. The paper also indicate how the kernel matrices can be updated recurrently by handling the ODE involved in the learning procedure. The performed experiments, in different scenarios, show promising results, compared to other methods from the literature.", "questions": "- Section 5.2 - Line 319: The heuristic sorting method for OHSVGP-k considers the minimization of k(x, x_{i-1}^{(j)}). Wouldn't it result in the \"farthest\" point, since the kernel is inversely proportional to the distance between the inputs? It seems that if the intent is to find the closest point, one should maximize the kernel value between the points. Appendix E.3 also covers this point.\n\n- How exactly the \"High dimensional time series prediction\" task was conducted with GPVAEs in the experiments? The description in Section 5.3 focus on the hyperparameters of the model components, but does not detail the overall setup.\n\n- Line 99: Wrong reference for the original VAE paper.", "limitations": "The authors state the limitations of their work.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 4, "originality": 4, "strengths_and_weaknesses": "The paper is well organized and written, with all the necessary background properly summarized for the reader. In that regard, both the related work and the baseline comparisons are also well covered. The experimental section is detailed and transparent in terms of settings and hyperparameters choices. I also praise the authors for the provided comprehensive supplementary material, which contains several interesting additional descriptions and experiments.\n\nThe main weak point might be the quantity of steps and decisions involved in the training of the proposed approach, with some particularities for each experimental scenario. A complete algorithm, summarizing all the calculations and methodological decisions in each case, would be great for better reader understanding.", "was_revised": true, "final_justification": "The authors addressed my points and I maintain my favorable score.", "rebuttal": "### Global Response:\n\nWe sincerely thank all reviewers for the constructive feedbacks and the overall positive reviews. We are encouraged that reviewers appreciated our **novel (CC7b, GFdF, Aht3)** and **technically sound (D83n, GFdF)** method that addresses the **critical (D83n)** catastrophic forgetting problem in online GPs. In general, the reviewers are happy with the **clarity (D83n, Aht3)** of the manuscript and our **comprehensive (CC7b, HjVk) experimental evaluation** of OHSVGP, which shows **strong performance (D83n, HjVk, GFdF, Aht3)** compared with baselines.\n\n---\n### Individual Response to Reviewer Aht3:\n\nBelow we address the questions from Reviewer Aht3.\n\n>Q1. A complete algorithm, summarizing all the calculations and methodological decisions in each case.\n \nWe agree that a complete algorithm will improve the manuscript and we will include it in future revision of the paper. For now, we give a brief sketch of the algorithm in the cases with conjugate Gaussian likelihood (time series regression) and with non-conjugate likelihoods (all the other experiments):\n\n* Case 1 (Gaussian likelihood). Assume we have the old kernel matrices ($K_{fu}^{old}$ and $K_{uu}^{old}$) and optimal variational parameters ($m_{old}$ and $S_{old}$) computed from the last task.\n    * Update $K_{fu}^{old}$ and $K_{uu}^{old}$ to $K_{fu}^{new}$ and $K_{uu}^{new}$, respectively, based on HiPPO recurrences to cover the time region of the current task (Eq. 4 and Appendix B).\n    * Since Gaussian likelihood is conjugate to GP prior, the optimal $m_{new}$ and $S_{new}$ are analytically tractable and can directly be computed based on $m_{old}$, $S_{old}$, and these kernel matrices without training (detailed formulas can be found in Appendix B in [1]). \n    * Construct the new posterior predictive based on $m_{new}$, $S_{new}$, and kernel matrices.\n\n* Case 2 (non-conjugate likelihood). Assume we have the old kernel matrices ($K_{fu}^{old}$ and $K_{uu}^{old}$) and variational parameters ($m_{old}$ and $S_{old}$) obtained from the last task.\n    * Update $K_{fu}^{old}$ and $K_{uu}^{old}$ to $K_{fu}^{new}$ and $K_{uu}^{new}$, respectively, based on HiPPO recurrences to cover the time region of the current task (Eq. 4 and Appendix B).\n    * Construct the online ELBO (Eq. 3), based on $m_{old}$, $S_{old}$, and kernel matrices, as training objective for $m_{new}$ and $S_{new}$, and train $m_{new}$ and $S_{new}$ by gradient-based optimization of the online ELBO.\n    * Construct the new posterior predictive based on $m_{new}$, $S_{new}$, and kernel matrices.\n\n[1] \"Streaming Sparse Gaussian Process Approximations\", NeurIPS 2017\n\n> Q2. Clarification of OHSVGP-k.\n\nYou are right. **This is a typo - OHSVGP-k actually considers $ \\argmax_x k(x, x_{i-1})$ (the results of OHSVGP-k-min and OHSVGP-k-max in Appendix E.3 should also be swapped).** We double-checked that our implementation of OHSVGP-k in facts considers $\\argmin_x \\sum_{d=1}^D \\frac{(x^{(d)} - x_{i-1}^{(d)})^2}{l_d^2}$ (where $D$ is the dimension of $x$ and $l_d$'s are per-diemension length-scales), which is equivalent to the maximization of the kernel values. Thank you for pointing it out and we will fix this typo in our next revision of the manuscript.\n\n> Q3. Background and experimental Setups of GPVAEs.\n\nIn GPVAE experiment, each data point is a multi-output/high-dimensional time series and GPVAE uses amortized inference based on an encoder to output an approximate GP posterior for each time series. We include the model specification of GPVAE used in this work in Appendix D.2 and we will include a pointer to it in the main text in our next revision of the manuscript. The dataset of time series is splitted into multiple tasks and the continual learning is achieved by (i) imposing an additional EWC loss to the encoder and decoder networks and (ii) updating the kernel matrices via HiPPO recurrence to keep adapting to the new time region associated with the time series in the new task.\n\n> Q4. Wrong reference of VAE paper.\n\nThank you for pointing out the typo. We will fix it in our next revision of the manuscript."}]

Author rebuttal:

Meta-review: The authors propose OHSVGP, a new online Gaussian Process model that leverages HiPPO to efficiently retain past information and update as new data arrives. Reviewers found the work to be both novel and elegant—particularly the use of HiPPO projection coefficients as time-evolving inducing variables—and effective based on the experimental results. The paper was also considered well presented, though I encourage the authors to further refine the presentation based on the reviewers’ feedback and discussions.
Decision: Accept (poster)