arXiv CSopen access

Sequential Inference for Gaussian Processes: A Signal Processing Perspective

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

arXiv:2604.28163v1 [eess.SP] 30 Apr 2026

Sequential Inference for Gaussian Processes: A Signal Processing Perspective∗ Daniel Waxman1,2 , Fernando Llorente3 and Petar M. Djurić2 1 2

Basis Research Institute, e-mail: [email protected]

Department of Electrical and Computer Engineering, Stony Brook University, e-mail: [email protected]

Computing and Data Sciences Directorate, Brookhaven National Laboratory, e-mail: [email protected] Abstract: The proliferation of capable and efficient machine learning (ML) models marks one of the strongest methodological shifts in signal processing (SP) in its nearly 100-year history. ML models support the development of SP systems that represent complex, nonlinear relationships with high predictive accuracy. Adapting these models often requires sequential inference, which differs both theoretically and methodologically from the usual paradigm of ML, where data are often assumed independent and identically distributed. Gaussian processes (GPs) are a flexible yet principled framework for modeling random functions, and they have become increasingly relevant to SP as statistical and ML methods assume a more prominent role. We provide a self-contained, tutorialstyle overview of GPs, with a particular focus on recent methodological advances in sequential, incremental, or streaming inference. We introduce these techniques from a signal-processing perspective while bridging them to recent advances in ML. Many of the developments we survey have direct applications to state-space modeling, sequential regression and forecasting, anomaly detection in time series, sequential Bayesian optimization, adaptive and active sensing, and sequential detection and decision-making. By organizing these advances from a signal-processing perspective, we intend to equip practitioners with practical tools and a coherent roadmap for deploying sequential GP models in real-world systems.

∗

Accepted for publication in IEEE Signal Processing Magazine. This is the authors’ accepted manuscript version. © 2026 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works by sending a request to [email protected].

D. Waxman et al./Sequential Inference for Gaussian Processes

Contents 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Gaussian Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1 A Motivating Example: Multivariate Normal Distributions . . 6 2.2 Extending to Infinite Dimensions: Gaussian Processes . . . . . 7 2.3 Troubles in Sequential and Online Estimation . . . . . . . . . 13 3 Approximate Gaussian Process Inference . . . . . . . . . . . . . . . 14 3.1 Sparse Gaussian Process Approximations . . . . . . . . . . . . 14 3.2 Spectral Gaussian Process Approximations . . . . . . . . . . . 16 4 Sequential Estimation via Basis Expansions . . . . . . . . . . . . . 18 4.1 Exact Inference via Kalman Filtering . . . . . . . . . . . . . . 18 4.2 Adding Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.3 Approximate Inference for Non-Conjugate Likelihoods . . . . . 22 4.4 A Bank of Filters Approach . . . . . . . . . . . . . . . . . . . 23 4.5 Extensions of Basis Expansions GPs . . . . . . . . . . . . . . 24 5 Markovian GPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5.1 Gaussian Processes as State Space Models . . . . . . . . . . . 25 5.2 Spatiotemporal Regression via Markovian Gaussian Processes 28 5.3 Extensions of Markovian Gaussian Processes . . . . . . . . . . 29 6 Sequential Estimation in Sparse (Variational) Approximations . . . 30 6.1 Sequential Estimation in Sparse Approximations . . . . . . . . 31 6.2 Sequential Estimation in Sparse Variational Approximations . 32 7 Other Approaches to Sequential Estimation . . . . . . . . . . . . . 35 7.1 Recursive Kernel and Online Regression Methods . . . . . . . 36 7.2 Variational State–Space and Hybrid Models . . . . . . . . . . 36 7.3 Expectation Propagation and Message–Passing Frameworks . 37 7.4 Monte Carlo and Particle Approaches . . . . . . . . . . . . . . 38 8 Comparisons to Deep Learning Models . . . . . . . . . . . . . . . . 38 8.1 Conceptual Comparison: Function–Space versus Parameter–Space Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 8.2 Diffusion Processes and the Markovian View . . . . . . . . . . 39 8.3 Sequential Deep Learning Architectures . . . . . . . . . . . . . 40 8.4 Deep Kernels and Hierarchical Gaussian Processes . . . . . . . 41 8.5 Nonstationary and Adaptive Representations . . . . . . . . . . 41 8.6 Outlook: Complementarity and Integration . . . . . . . . . . . 42 9 Applications of Gaussian Processes to Sequential Signal Processing 43 9.1 Sequential Regression and Forecasting . . . . . . . . . . . . . 43

D. Waxman et al./Sequential Inference for Gaussian Processes

9.2 Distributed and Decentralized Online Inference . . . . . . . . 9.3 Robust Sequential Learning and Anomaly Detection . . . . . . 9.4 Sequential Bayesian Optimization . . . . . . . . . . . . . . . . 9.5 Adaptive and Active Sensing . . . . . . . . . . . . . . . . . . . 10 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43 44 45 45 45 47 48

1. Introduction A fundamental task of statistical signal processing (SP) is the estimation of a signal that is partially observed, embedded in noise, or both. This simplystated goal is rather difficult in practice and often relies on assumptions we make about the signal; to make these assumptions explicit, the signal is often represented as a state-space model, where an observed state (or signal) ft evolves in time in a way that depends only on ft−1 , emitting noisy observations yt that depend only on the current state ft . A common set of assumptions is that the state evolves as a linear function over time with a Gaussian random walk and that observations depend linearly on the state and are embedded in Gaussian noise. These assumptions are rather useful, and the estimation of ft can be optimally done via the Kalman filter. At the same time, the assumptions are also fairly restrictive. In this tutorial, we discuss Gaussian processes (GPs), a methodological tool with a number of desirable properties for the signal estimation task. First, GPs naturally cope with nonlinear signals and can handle a broad variety of estimation tasks as a result. They accomplish this by appealing to the more abstract property of correlation between samples of a signal and making assumptions about the structure of correlation, rather than on a specific parametric form of the model. Second, GPs naturally support Bayesian inference and provide quantitative uncertainty estimates in signal estimation and reconstruction. This is a critical feature in many SP-based applications, e.g., in healthcare. Finally, GPs operate over a continuum, i.e., they do not assume a discrete-time model of a signal. They, therefore, can be applied to a wide range of problems that may exhibit irregular sampling. For all their theoretical elegance, GPs are often marred by computational challenges. In particular, in their classical formulation that is often used in machine learning (ML) and statistics, inference of a GP using N samples of a signal requires the inversion of a dense N ×N matrix, which results in O(N 3 )

D. Waxman et al./Sequential Inference for Gaussian Processes

computational complexity. This is particularly problematic in the SP context, where inference of long time series is typical. Furthermore, the classical GP formulation provides no straightforward way to update our beliefs based on the first N samples after observing the (N + 1)-th sample without incurring at least O(N 2 ) operations. These computational limitations have limited the use of traditional GP inference in sequential or resource-constrained SP and hindered their widespread adoption. These computational challenges with scaling data are pernicious, and while particularly problematic in SP, they have been extensively studied in neighboring communities as well. Many algorithms have been proposed in the ML literature for scalable GP inference, with a modern trend toward formulations that are explicitly compatible with sequential and online operations. While these operations are fundamentally related to statistical SP and may be best understood in that setting, these innovations remain underutilized. One such approach relies on basis expansions, which approximate a GP prior using a finite set of features, chosen deterministically or randomly, and converge to the true GP prior as the number of features increases. This converts the complicated function-space inference of a GP into a parametric linear model. Examples include random Fourier features (RFF) and Hilbert–space GPs, both of which yield state–space structures that can be updated in real time through standard filtering equations (and, in particular, Kalman filtering, if the observations are embedded in Gaussian noise). Similarly, sparse and variational approximations summarize the influence of past observations through a compact set of inducing variables to allow tractable inference with constant memory and linear complexity. The Markovian formulation of a GP provides an alternative perspective by expressing many common GP priors as linear stochastic differential equations (SDEs). Although more mathematically elaborate, these models, too, allow exact inference through Kalman filtering and smoothing. Each of these approaches replaces the infinite–dimensional representation of the GP with an equivalent parametric model whose parameters evolve sequentially with data. From an SP viewpoint, these approximations return the GP framework to familiar territory. They recover the recursive structure of state–space models while retaining the interpretability and uncertainty quantification of Bayesian inference. The resulting algorithms extend classical tools such as recursive least squares and Wiener filtering into probabilistic forms that remain stable, adaptive, and computationally efficient. This synthesis also provides new

D. Waxman et al./Sequential Inference for Gaussian Processes

insight into the relationship between kernel methods and linear systems theory and reveals that many stationary kernels correspond to impulse responses of SDEs. The connection between continuous–time dynamics, spectral representations, and sequential estimation lies at the core of the modern SP interpretation of GPs. The methods discussed in this article show that GPs, when equipped with appropriate approximations, constitute practical and scalable tools for real–time inference. We emphasize their implementation as linear–Gaussian models, their extension to non–Gaussian likelihoods, and their integration with ensemble and distributed architectures. The resulting algorithms operate in a finite-dimensional parameter space while maintaining consistent probabilistic semantics. This parametric viewpoint also aligns GPs with recent developments in ML, where deep kernel learning, variational state–space models, and diffusion processes share a common goal of unifying functional priors with sequential estimation. The remainder of this article proceeds as follows. Section 2 introduces the foundation of GP regression and its finite–dimensional interpretation. Section 3 discusses various existing approaches in scalable GP inference, such as sparse GP regression, and regression via basis expansions. Section 4 presents the sequential version of the basis expansion case. Section 5 presents Markovian GPs and their equivalence to SDEs, linking GP inference to Kalman filtering and smoothing. Section 6 examines variational inference, non–conjugate likelihoods, and extensions to online ensembles and robust inference. Subsequent sections explore applications to forecasting, distributed estimation, and adaptive sensing, and conclude with a comparison to modern deep learning architectures. Our objective is to demonstrate that GPs, formulated through parametric and state–space approximations, provide computationally efficient and reliable tools for sequential SP. 2. Gaussian Processes In this section, we introduce GPs in their typical (non-sequential) formulation from the function-space perspective, where the model is parameterized by a covariance (kernel) function. This includes a motivating example of a joint Gaussian model of signal outputs, which generalizes to GP regression. We discuss several common kernels and how they can impact the properties of the resulting models. We illustrate that exact inference from this perspective is not scalable in online or sequential settings and requires simplifying

D. Waxman et al./Sequential Inference for Gaussian Processes

approximations in both sequential and non-sequential practices. 2.1. A Motivating Example: Multivariate Normal Distributions We begin by considering a simple, finite-dimensional example that will naturally extend to the infinite-dimensional case of GPs. Specifically, let us consider a one-dimensional, irregularly sampled, continuous-time signal, y(t). We observe samples at t ∈ {0.0, 1.0, 2.5}, with the goal of inferring the unobserved value y(2.0). For now, we assume that our observations are noiseless, but we will remove this assumption later. Our first example of GP inference will arise through building a finite-dimensional, multivariate normal distribution over the signal values y = [ y(0) y(1) y(2.0) y(2.5) ]⊤ , i.e., y ∼ N (µ, Σ).

(1)

This requires specifying the mean vector, µ, and covariance matrix, Σ. For simplicity, we will assume that µ = 0. The covariance matrix is more involved, as it tells us how samples are related to one another. One intuitive assumption decay with separation in t, for example, is that covariances ′ ′ 2 cov y(t), y(t ) = exp(−(t − t ) ). We may thus revise Eq. (1) as    1.000 0.368 0.018 0.002  0.368 1.000 0.368 0.105   (2) y ∼ N 0, 0.018 0.368 1.000 0.779 . 0.002 0.105 0.779 1.000 The core property that makes a Gaussian model, like Eq. (2), fundamental in GP inference is that marginal and conditional distributions remain analytically tractable once the joint distribution is specified. For clarity, we adopt the following notational conventions: we use ∗ to denote quantities associated with the unobserved y(2.0), and ⊙ to denote quantities corresponding to the observed signal. In particular, the observed signal is ⊤ y⊙ = y(0) y(1) y(2.5) , (3) which is described by the covariance matrix   1.000 0.368 0.002 Σ⊙⊙ = 0.368 1.000 0.105 . 0.002 0.105 1.000

(4)

D. Waxman et al./Sequential Inference for Gaussian Processes

According to the properties of the multivariate Gaussian distribution, this submatrix is obtained by deleting the rows and columns corresponding to y(2.0) from y. The covariance between the observed and unobserved signals, known as the cross-covariance, is therefore ⊤ (5) Σ⊙∗ = 0.018 0.368 0.779 . Following standard results of probability theory, we can derive exactly the conditional probability of the unobserved value y(2.0) given the observed samples, 2 p y(2.0) y(0), y(1), y(2.5) = N (µpost , σpost ), (6) 2 is its conwhere µpost is the conditional (posterior) mean of y(2.0) and σpost ditional variance. The mean µpost is a weighted average of the observed samples, ⊤ −1 µpost = Σ⊤ -0.104 0.328 0.744 y(0) y(1) y(2.5) . (7) ⊙∗ Σ⊙⊙ y⊙ ≈ | {z } weights

2 reduces the marginal variance Σ∗∗ by the weighted The posterior variance σpost cross-covariances −1 2 σpost = Σ∗∗ − Σ⊤ ⊙∗ Σ⊙⊙ Σ⊙∗ ≈ 0.3016. | {z }

(8)

weights

This process is illustrated in Fig. 1a. Through our simple probabilistic model, we not only impute y(2.0), but also quantify its uncertainty, which produces a posterior distribution that can be computed via simple linear algebra. In GP regression, the same principle extends to functions: the goal is to infer a posterior probability distribution over functions that, when covariances are specified appropriately, recovers the same posterior distribution p(y(2.0)|y⊙ ). GP regression refers to the Bayesian approach of assigning a GP prior over an unknown function and computing its posterior distribution after observing data. This enables both prediction at new inputs and rigorous uncertainty quantification. 2.2. Extending to Infinite Dimensions: Gaussian Processes GP regression is a natural extension of the procedure above for finite-dimensional random vectors to the case of random functions. To define a GP, we specify a mean function µ(t) and a covariance function κ(t, t′ ), which together

D. Waxman et al./Sequential Inference for Gaussian Processes

x(t)

2 0 2

0.0

0.5

1.0

1.5

2.0

2.5

1.5

2.0

2.5

(a)

x(t)

2 0 2

0.0

0.5

1.0

t (b)

Fig 1. From finite- to infinite-dimensional inference. (a) Inferring the unobserved signal at t = 2 using the correlations among a finite set of signal values. (b) Extending this idea to the function space: each gray curve represents a random function drawn from the GP posterior that gives rise to the same posterior distribution over y(2.0). GP regression generalizes finite-dimensional Gaussian inference to a continuous domain and provides a unified framework to specify, compute, and sample posterior distributions over functions.

D. Waxman et al./Sequential Inference for Gaussian Processes

determine the joint distribution p(y) for any finite collection of observed and unobserved samples y. The defining property of a GP is that this joint distribution, p(y), is always multivariate Gaussian. More generally, we may consider d-dimensional inputs, x ∈ Rd , provided that µ(x) and κ(x, x′ ) are specified accordingly. Through this construction, we obtain a random function f : Rd → R, which may be combined with a likelihood p(y | f (x)) to describe noisy observations. Formally, we write f ∼ GP(µ, κ), y | f, x ∼ p(y | f (x)).

(9) (10)

As the values of f are observed only through y, we will sometimes refer to f as the “latent function.” A common special case that makes inference analytically tractable arises when the likelihood p(y | f (x)) is Gaussian, i.e., y = f (x) + ε for some i.i.d. ε ∼ N (0, σε2 ). In this case, we may proceed using the same conditional normal identities as before to obtain analytical posterior distributions of the latent function values at any collection of test points, f ∗ = [f∗1 · · · , f∗P ], as f ∗ | X⊙ , X∗ , y⊙ ∼ N (µpost , Σpost ), where the posterior mean and covariance of f ∗ are given by −1 µpost = Σ∗⊙ Σ⊙⊙ + σε2 I y⊙ −1 Σpost = Σ∗∗ − Σ∗⊙ Σ⊙⊙ + σε2 I Σ⊙∗ ,

(11)

(12) (13)

where I is the identity matrix. For the resulting multivariate normal distribution to be mathematically well posed, κ(x, x′ ) must yield positive semidefinite covariance matrices, and no other mathematical restrictions are imposed. When the likelihood p(y | f (x)) is non-Gaussian, exact inference is no longer tractable, and we must rely on approximate inference, most commonly variational inference, which we will introduce shortly. As illustrated in Fig. 1b, GP regression generalizes the finite-dimensional Gaussian model of Fig. 1a to a posterior distribution over random functions. This means we obtain a stochastic process as our posterior, whereby any sample is a function f : Rd → R. Whereas Fig. 1a depicts inference at a single unobserved point based on correlations among a few samples, Fig. 1b shows multiple random functions (gray curves) drawn from the GP posterior. These

D. Waxman et al./Sequential Inference for Gaussian Processes

samples of the posterior function are also known as realizations of the GP. Each realization is consistent with the observed data and reflects the uncertainty specified by the covariance function. This visualization illustrates the key conceptual shift: rather than inferring a single value, GP regression infers a distribution over functions. Since a GP is just a random function for which any finite collection of points follows a multivariate Gaussian distribution, and since Eq. (11) is Gaussian, the posterior function f (x) is also a GP. Of the two functions that must be specified, the covariance function κ(x, x′ ) is the more informative and practically significant. While the mean function µ(x) provides a time-varying baseline, the covariance function is the component that represents the relationships between pairs of inputs x and x′ . As a result, and in connection to classical methods in ML, it is often called the kernel of the GP. The specific form of κ(x, x′ ) represents one of the most influential modeling choices, as it is primarily responsible for determining the shape of the posterior function and even influences the types of approximations or inference methods that can be employed later. There are many potential covariance functions, but several choices stand out as popular and useful. The most common class of kernels relates points according to their Euclidean distance. For example, the squared exponential (SE) kernel (sometimes called the radial basis function (RBF) kernel or exponentiated quadratic (EQ) kernel) assigns an exponentially decaying covariance between points x and x′ , ∥x − x′ ∥22 ′ 2 , (14) κSE (x, x ) = σf exp − 2ℓ2 where ∥x−x′ ∥22 = (x−x′ )⊤ (x−x′ ), and σf2 and ℓ are hyperparameters known as the process variance and length scale, respectively. The process variance determines the prior marginal variance of any given function evaluation, while the length scale controls how quickly correlations decay with distance, that is, how far apart two inputs can be before they lose relevance. As noted earlier, the kernel determines many properties of the resulting functions. One important property is their smoothness. For example, the SE kernel specifies a prior distribution over infinitely differentiable functions. In many real-world applications, however, this assumption can be overly restrictive, and the well-known Matérn family of kernels is used instead (Rasmussen and Williams, 2005, pp. 84–85). Each kernel in this family is parameterized by a smoothness parameter ν > 0, which controls the differentiability of

D. Waxman et al./Sequential Inference for Gaussian Processes

the resulting functions: larger values of ν correspond to smoother functions. The Matérn family is quite general, being defined for any ν > 0, though its form simplifies significantly for half-integer ν. For example, the Matérn-1/2 and Matérn-3/2 kernels correspond to certain non-differentiable and oncedifferentiable functions, respectively, and are given by ∥x − x′ ∥2 ′ 2 ; (15) κMat-1/2 (x, x ) = σf exp − ℓ ! ! √ √ ′ ′ 3∥x − x ∥ 3∥x − x ∥ 2 2 κMat-3/2 (x, x′ ) = σf2 1 + exp − , (16) ℓ ℓ p where ∥x − x′ ∥2 = (x − x′ )⊤ (x − x′ ). The Matérn-ν kernels will prove particularly useful in our later discussion of Markovian GPs. Each of the above kernels is an example of a stationary kernel, i.e., a kernel that depends only on the difference r = x − x′ . While we will later consider GPs with non-stationary kernels when discussing connections to deep learning, GP priors with stationary kernels remain highly versatile and unlock certain theoretical properties that are useful for approximation, e.g., a spectral decomposition. Building on the idea of spectral decompositions, one popular and expressive kernel is the spectral mixtures (SM) kernel (Wilson and Adams, 2013). The SM kernel represents the power spectral density (PSD) of the covariance function, that is, the Fourier transform of a stationary covariance, as a mixture of Gaussians, which allows it to approximate a broad class of stationary covariance functions. We will revisit the SM kernel in subsequent sections. The choice of kernel is crucial in GP regression, as it affects both the interpolation and extrapolation properties of the resulting model. For example, SE, Matérn-1/2, and Matérn-3/2 kernels each imply different smoothness properties: functions that are realizations from a GP with SE kernel are infinitely differentiable almost surely, whilst realizations from GPs with Matérn-1/2 and Matérn-3/2 kernels are almost surely non-differentiable and once-differentiable functions, respectively. In practical terms, this smoothness determines how rapidly the modeled function can vary between observations. The SE kernels produce very smooth, gently varying functions, while Matérn kernels allow for rougher or more abrupt changes, depending on the value of ν. More expressive kernels, such as the spectral mixture kernel, can also model quasi-periodic patterns that the SE kernels cannot. We illustrate several different choices of kernels in Fig. 2.

D. Waxman et al./Sequential Inference for Gaussian Processes

Matern-1/2

Matern-3/2

Spectral mixture

f(x)

2 0 2 0

6 0

f(x)

2 0 2

Fig 2. An illustration of prior (top) and posterior (bottom) samples from GPs with different kernels. The posterior plots include an extrapolation region (x ≥ 4) with no training data. While the SE and Matérn kernels yield smooth or locally correlated behavior that reverts quickly to the prior mean, the spectral mixture kernel captures the underlying quasiperiodic pattern and extrapolates more faithfully beyond the observed region.

D. Waxman et al./Sequential Inference for Gaussian Processes

Generally speaking, we may combine different covariance functions by either adding or multiplying them (Rasmussen and Williams, 2005, Sec. 4.2.4). In spatiotemporal settings, we may separate the input as x = [s t], where s denotes spatial coordinates and t denotes time. If the kernel is expressed as a product κ(x, x′ ) = κ(s, s′ )κ(t, t′ ), it is referred to as separable kernel. GPs are generally sensitive to the specific values of their hyperparameters. The predominant approach for selecting these hyperparameters is empirical Bayes, in which the marginal likelihood is maximized. Owing to the simple Gaussian form that the GP reduces to when conditioned on finite data, the resulting log marginal likelihood is simply log p(y⊙ ) = log N (y⊙ | 0, Σ⊙⊙ + σε2 I). This objective admits closed-form expressions for its gradients and is commonly optimized using non-linear conjugate gradients (as is done in the gpml toolbox for MATLAB (Rasmussen and Williams, 2005)), or by gradient descent and quasi-Newton methods (as is more commonly done for software in the Python ecosystem, e.g., GPyTorch (Gardner et al., 2018) or GPJax (Pinder and Dodd, 2022)). 2.3. Troubles in Sequential and Online Estimation GPs are now relatively mature tools within the ML and statistics literature and have been successfully applied to many challenging and practically important problems in SP. An important limitation that continues to restrict their broader use, however, is their scalability. While the goal of this article is to present scalable and sequential GP methods, it is worth examining the source of these challenges. The primary obstacle to scalability arises from the matrix inversion required in the posterior GP computations, Eqs. (12) and (13). Although explicit inversion may be avoided by solving equivalent linear systems, computing µpost and Σpost still incurs an O(N 3 ) computational cost and an O(N 2 ) memory requirement, where N is the number of data points. These quadratic and cubic scalings severely limit the applicability of exact GPs to large data sets. There were several early efforts to address these problems (or the related problems in kernel machines) in the SP literature. For example, kernel recursive least squares (Engel, Mannor and Meir, 2004) includes a recursive form of kernel regression, with a strategy for maintaining a sort of dictionary set of points. The kernel recursive least squares (KRLS) algorithm is interpreted under a Bayesian framework for GP regression by Van Vaerenbergh, Lázaro-Gredilla and Santamarı́a (2012), who also provide modifications for

D. Waxman et al./Sequential Inference for Gaussian Processes

time-varying environments. This connection between KRLS algorithms and sequential GP inference is made clear in the tutorial article of Pérez-Cruz et al. (2013). Unfortunately, apart from windowed approaches, there is typically no way to further reduce the computational cost of exact sequential GP inference beyond O(N 2 ) per update in the worst case. 3. Approximate Gaussian Process Inference As a result of the computational burden of exact GP inference via the kernel formulation, practical sequential inference requires the use of approximate methods. Over the past 15 years, a variety of powerful approximations have been developed that are far more amenable to online and sequential settings. In the following, we introduce two such families of approximations that will later be adapted for sequential inference: (i) sparse variational methods, which emulate the GP posterior using a set of “inducing points,” and (ii) basis expansion methods, which approximate the GP through a finite linear model. 3.1. Sparse Gaussian Process Approximations One form of approximation, likely the most common in ML, is a sparse GP. Whereas exact GP inference relies on all N training inputs when making predictions, sparse GP approximations seek to summarize the model using a smaller set of M inducing points, where M ≪ N . These inducing points generally consist of inducing point locations, xu = [xu1 · · · xuM ], and a corresponding multivariate distribution over the associated latent function values, u ∼ N (m, S).

(17)

Recall that in exact GP inference, the joint distribution over y⊙ and y∗ is used for all predictions. In sparse GP approximations, we assume that the dependence within this joint distribution is mediated solely through the inducing variables. This leads to an approximate joint distribution q(y⊙ , y∗ ), defined as Z p(y⊙ , y∗ ) ≈ q(y⊙ , y∗ ) = p(y∗ | u)p(y⊙ |u)p(u) du. (18) u

There are many ways to choose the inducing points and their locations (Quinonero-Candela and Rasmussen, 2005), but the most widely used ap-

D. Waxman et al./Sequential Inference for Gaussian Processes

proach today is based on variational inference (Titsias, 2009), which places a distribution q(u) on the inducing points which is later optimized. Variational inference (VI) is a general framework for approximate Bayesian inference that consist in approximating the true posterior distribution by a tractable distribution from a simpler family. The best approximation is found by minimizing the Kullback-Leibler (KL) divergence, which is also intractable. In practice, the minimization of the KL divergence between the approximate and true posterior is equivalent to maximizing the so-called evidence lower bound (ELBO) on the marginal likelihood, which leads to a tractable objective for approximate Bayesian inference. Let p(f |y) and q(f ) be respectively the true posterior and its variational approximation, the KL divergence is Z q(f ) df (19) KL(q(f ) ∥ p(f |y)) = q(f ) log p(f |y) Z = q(f ) [log q(f ) − log p(y|f ) − log p(f ) + log p(y)] df (20) = log p(y) −Eq(f ) [log p(y|f )] + KL(q(f ), p(f )) . | {z }

(21)

−L(q)

Hence, so long as expectations w.r.t. q(f ) are tractable, VI gives us a tractable lower bound on the log-marginal likelihood, i.e., log p(y) ≥ L(q), such that maximizing L(q) is equivalent to minimizing KL(q(f ) ∥ p(f |y)). Note that L(q) contains a data-fit term and a regularization term, balancing the quality of the approximation with its complexity. In the context of sparse GPs, the ELBO allows us to optimize the locations of the inducing points and the hyperparameters simultaneously. In the sparse variational GP framework, we posit a variational distribution over the inducing variables, q(u) = N (m, S), and optimize the parameters m, S, and the kernel hyperparameters by maximizing the ELBO. Titsias (2009) showed that, when the variational parameters are optimized analytically (the so-called collapsed case), the bound simplifies to a closed-form objective that automatically regularizes the inducing points and avoids overfitting. Hensman, Fusi and Lawrence (2013) later generalized this formulation to the non-collapsed or stochastic variational GP (SVGP), in which the bound is optimized directly with respect to m and S and the hyperparameters using stochastic gradient descent and mini-batches of data. This stochastic formula-

D. Waxman et al./Sequential Inference for Gaussian Processes

tion makes GP training scalable to large data sets and facilitates incremental or mini-batch updates. 3.2. Spectral Gaussian Process Approximations Another approach to scalable GP regression, perhaps more reminiscent of traditional SP methods, is based on spectral approximation. Spectral methods aim to emulate the GP by matching its power spectral density (PSD). As previously mentioned in our discussion of the spectral mixture kernel, the PSD exists for GPs with stationary kernels, a result known as Bochner’s Theorem (Rasmussen and Williams, 2005, pp. 82). In this case, we may rewrite the kernel κ(x, x′ ) as κ(r). Subject to mild measurability constraints, the PSD of the kernel function is given by the Fourier transform of κ(r) (Rasmussen and Williams, 2005, pp. 82). For many common choices of kernel, the resulting PSD has a simple form and is easy to sample from. For example, the SE kernel results in a Gaussian PSD, up to a scaling factor, and Matérn-1/2 corresponds to a heavier-tailed Cauchy distribution. More expressive kernels, such as those from the MHM (Dowling, Sokól and Park, 2021) or SM families, are carefully constructed to have convenient spectral representations (in both of these cases, via mixture distributions). The fact that many common kernels induce GP priors with an easy-tosample PSD motivates a particular parametric approximation of the GP. In particular, the existence of the PSD implies a spectral representation of the kernel, Z ′ κ(x, x ) = exp(is⊤ (x − x′ ))S(s) ds, (22) RD

√ where S(s) is the PSD of the covariance function, and i = −1 denotes the imaginary unit. Since S(s) is easily sampled from, we create an approximate, parametric prior via direct Monte Carlo integration of Eq. (22). As the PSD of the kernel is given by its Fourier transform, the resulting approximate inference method (via naı̈ve Monte Carlo) is known as the random Fourier feature (RFF) approximation, introduced for frequentist kernel machines by Rahimi and Recht (2007). Application to GPs, from the perspective of Bayesian trigonometric regression, known as the sparse spectrum GP, was introduced by Lázaro-Gredilla et al. (2010). To avoid confusion with state-space GP formulations, as well as a version of the RFF-GP using opti-

D. Waxman et al./Sequential Inference for Gaussian Processes

mized frequencies, also introduced by Lázaro-Gredilla et al. (2010), we will refer to the resulting approximation as the “RFF-GP”. The RFF-GP makes only marginal refinements to a direct Monte Carlo approximation of Eq. (22) by exploiting the symmetric nature of the PSD of a stationary kernel. In the resulting approximation, F/2 spectral frequencies s1 , . . . , sF/2 are sampled from S(s), and a linear basis expansion established, r 2 [sin(x⊤ s1 ) cos(x⊤ s1 ) · · · sin(x⊤ sF/2 ) cos(x⊤ sF/2 )]⊤ ∈ RF . (23) φ(x) = F Thus, the GP prior f ∼ GP(0, κ) is reduced to a linear model, f (x) ≈ φ(x)⊤ θ,

i.i.d.

θ ∼ N (0, σf2 ).

(24)

The reduction to a linear model is particularly useful for sequential estimation, as it opens the door to the application of many well-established tools from classical SP. The use of a direct Monte Carlo approximation of the kernel is conceptually simple but can be inefficient, particularly when x is low-dimensional. For example, orthogonal random features (Yu et al., 2016) or quasi-Monte Carlo (Avron et al., 2016) can improve inference. An orthogonal approach is to abandon random features in favor of deterministic quadrature of Eq. (22). This is the fundamental approach of Hilbert space GPs (HSGPs) (Solin and Särkkä, 2020), so-called because of their use of Hilbert space methods from the analysis of partial differential equations (PDEs). As they correspond to certain quadrature rules, HSGPs are simplest to define over 1-dimensional inputs and are extended in a grid-like fashion. Additionally, they must be defined over some compact set for which the approximation is valid. The resulting approximation of the HSGP is another linear basis expansion (Solin and Särkkä, 2020) 1/2 kπ sin 2L (x + L) kπ √ , k = 1, . . . , F, (25) ϕk (x) = S 2L L where, we recall, S(s) denotes the PSD of the kernel. The HSGP thus inherits many convenient properties for sequential inference, similar to the RFF-GP. Using the grid-like extensions to multivariate inputs, HSGPs suffer from exponentially increasing costs. However, additive approximations (i.e., where f ∼ f1 + · · · + fD , where fd is a one-dimensional GP with xd as input) have

D. Waxman et al./Sequential Inference for Gaussian Processes

shown strong empirical performance (Solin and Särkkä, 2020; Waxman and Djurić, 2024a). The HSGP approximation is related to work on the stochastic PDE representation of GPs (Lindgren, Bolin and Rue, 2022), which efficiently solves GP regression problems for a variety of spatial fields. The key difference is that works in the SPDE approach to GPs typically consider a basis expansion around each spatial point s in a mesh, whilst HSGPs consider a global basis expansion. For more discussion, see Section 7.2 of (Lindgren, Bolin and Rue, 2022). 4. Sequential Estimation via Basis Expansions Finite basis-expansion approximations of GP priors are particularly powerful for sequential estimation, as they result in tractable finite-dimensional models for which a wide range of recursive estimation methods exist. In the conjugate prior case, these models reduce to linear-Gaussian dynamical systems that are directly amenable to Kalman filtering. In this section, we explore such methods and discuss how the connection to filtering theory extends naturally to non-conjugate likelihoods, dynamic environments, and online ensembling. 4.1. Exact Inference via Kalman Filtering The potential of the RFF-GP as a tool in sequential inference was first noted, to the best of our knowledge, by Gijsberts and Metta (2013), who derived incremental updates to the Cholesky decomposition of the covariance of θ in the RFF approximation. From an SP perspective, the existence of such a solution is intuitive. Indeed, the RFF approximation of a GP is simply a Bayesian linear model that naturally allows a trivial state-space representation, θ t = θ t−1 yt = φ(xt )⊤ θ t + εt , i.i.d.

t = 1, 2, . . .

(26)

Together with the prior θ 0 ∼ N (0, σf2 I), this state–space formulation yields recursive update rules that are formally identical to those of the Kalman filter. In this view, the RFF-based GP behaves as a Bayesian linear dynamical system with a static state vector θ and Gaussian observation

D. Waxman et al./Sequential Inference for Gaussian Processes

noise. Each new data point incrementally refines the posterior mean and covariance of θ and permits efficient online inference without recomputing the full covariance matrix. 4.2. Adding Dynamics In SP applications, it is often necessary to perform online estimation in dynamic environments, where the underlying response function f (xt ) may itself evolve over time, denoted as ft (xt ). When using basis expansions of GPs, several strategies can be employed to account for such time-varying behavior. The first approach, common in spatiotemporal regression problems, is to incorporate a time-dependent kernel. That is, instead of modeling a timeinvariant function f (x), we explicitly model a time-varying function f (x, t) (Cressie and Wikle, 2011). This formulation has the advantage of preserving the interpretability of the GP prior. Moreover, since it retains the static state–space structure of Eq. (26), the filtering solution obtained via the Kalman filter is the complete posterior distribution over θ; i.e., after a single forward pass over the data, we can predict f (x, t∗ ) at arbitrary t∗ using the corresponding p(θ | y1:T ). The main drawback of the time-dependent kernel approach to dynamic basis expansions is a gradual loss of spatial resolution as t increases. This occurs because the finite set of basis functions must simultaneously represent spatial structure and temporal evolution. As time progresses, the fixed F -dimensional coefficient vector θ must account for an increasing number of effective spatial frames, which causes its representational capacity to become more thinly allocated across time and results in coarser spatial detail. This effect, illustrated in Fig. 3, contrasts with the random-walk formulation introduced next, which maintains spatial fidelity by evolving θ t dynamically over time. An alternative approach is to introduce a nontrivial state–space model for the time-varying coefficients θ t . The simplest such model assumes that θ t follows a random walk, obtained by adding isotropic Gaussian noise to the parameter vector: θ t = θ t−1 + εrw,t , yt = φ(xt )⊤ θ t + εt , i.i.d.

t = 1, 2, . . .

(27)

2 where εrw,t ∼ N (0, σrw I). This formulation was adopted by Lu et al. (2023) in their so-called dynamic RFF-GP, and more generally for basis-expansion

D. Waxman et al./Sequential Inference for Gaussian Processes

t = 142

0.6

RMSE

t = 75

D-RFF-GP

t = 15

0.4 0.2

Dynamic RFF-GP

0.4 20

ST-RFF-GP

f(x * , t)

0.4

100 120 140

100

(b)

0.2

(a)

time step

0.0

0.8

ST RF-GP

0.2 0

0.6

0.5 0.0 0.5

time step

125

(c)

Fig 3. Illustrating the relative performance of the dynamic RFF-GP and spatiotemporal RFF-GP. (a) Snapshots of predictions at different points of online estimation of a spatiotemporal field. Note the decline in resolution of the ST-RFGP as time advances. (b) The error over time of each method, measured as RMSE of the mean. (c) Predictions at a particular point in the spatiotemporal field, x∗ = (0.4, 0.6).

models by Waxman and Djurić (2024a). In both cases, significant empirical benefits have been observed in time-varying scenarios for both RFF-based GPs and HSGPs. Under the augmented model Eq. (27), we obtain online updates once again via the Kalman filter: the predictive step propagates the posterior mean and covariance of θ t through the random-walk dynamics, θ̂ t|t−1 = θ̂ t−1|t−1 ,

2 Pt|t−1 = Pt−1|t−1 + σrw I,

(28)

and the measurement update incorporates the new observation yt through the feature vector φ(xt ): −1 (29) Kt = Pt|t−1 φ(xt ) φ(xt )⊤ Pt|t−1 φ(xt ) + σε2 , θ̂ t|t = θ̂ t|t−1 + Kt yt − φ(xt )⊤ θ̂ t|t−1 , Pt|t = (I − Kt φ(xt )⊤ )Pt|t−1 . (30) These recursive updates efficiently track the time-varying coefficients θ t and, by extension, the evolving function ft (x), thereby yielding an online Bayesian filtering solution for dynamic GP regression.

D. Waxman et al./Sequential Inference for Gaussian Processes

t = 150

t = 400

t = 700

t = 950

f(x, t)

2 1 0 0.5

1.0

0.0

0.5

1.0

0.0

0.5

1.0

0.0

0.5

1.0

trajectory (PCA)

1 0.0

Dynamic (time-colored)

Static

True

Fig 4. Illustration of the adaptive behavior of an RFF-GP with time-varying parameters. Top: True time-varying function and corresponding model predictions with (colored) and without (gray) random-walk dynamics. Recent data points are shown as darker dots, with older points omitted for clarity. Bottom: Evolution of the parameter vector θ t , visualized through a two-dimensional PCA projection (PC1 vs. PC2) that reveals the gradual adaptation of the model over time.

Compared with the spatiotemporal kernel approach, the random-walk state–space model has the advantage of retaining the entire finite-dimensional approximation for a single, time-varying spatial frame. As a result, it tends to produce crisper spatial reconstructions, albeit at the cost of requiring O(T ) storage and an additional backward pass over the data to obtain smoothing estimates. The potential effectiveness of the random-walk formulation of dynamic RF GPs is illustrated in Fig. 4. More generally, the state–space evolution of θ t can be extended to include control inputs and non-isotropic process noise to obtain a variety of time-varying formulations. This leads to parameterized dynamics of the form (cf. Llorente and Djurić (2024)) θ t = at θ t−1 + ut + εθ,t , yt = φ(xt )⊤ θ t + εt ,

t = 1, 2, . . .

(31)

where at is a scalar autoregressive coefficient, ut is a control input, and

D. Waxman et al./Sequential Inference for Gaussian Processes

i.i.d.

εθ,t ∼ N (0, C) is GP noise with covariance C. The random-walk model 2 I. of Eq. (27) is recovered as the special case of at = 1, ut = 0, and C = σrw Alternative parameterizations of these quantities result in a variety of dynamic formulations for time-varying functions discussed in the literature. One example is back-to-prior (B2P) forgetting, first introduced for recursive kernel-based learning (Van Vaerenbergh, Lázaro-Gredilla and Santamarı́a, 2012). B2P forgetting blends the current posterior knowledge with the prior, controlled by a forgetting factor λ ∈ [0, 1], such that λ = 1 corresponds to fully trusting the current posterior (no forgetting), and λ = 0 corresponds to reverting to the prior at each time step. B2P forgetting can be recovered from the general state–space formulation in Eq. (31) by setting √ ut = 0, C = (1 − λ)σf2 I. (32) at = λ, This choice introduces controlled diffusion around the prior and allows the model to adapt gradually to nonstationary or time-varying functions while it retains stability and bounded uncertainty. 4.3. Approximate Inference for Non-Conjugate Likelihoods The reduction of the GP model to a linear–Gaussian state–space form has proven extremely convenient for online and sequential inference. This perspective is equally useful in the case of non-conjugate likelihoods, such as those arising in GP classification. In such cases, we obtain a modified state–space model of the form θ t = at θ t−1 + ut + εθ,t , ft = φ(xt )⊤ θ t , yt | ft ∼ p(yt | ft ),

(33) t = 1, 2, . . .

where, for example, the likelihood p(yt | ft ) may be softmax-based in GP classification or Poisson for count data. For non-conjugate models, a range of filtering-based approximations can be employed. These include simple approaches such as the Laplace approximation (Kass, Tierney and Kadane, 1991) and more advanced methods such as posterior linearization (Tronarp, Garcı́a-Fernández and Särkkä, 2018). A comparison of these techniques is shown in Fig. 5. Empirically, iterated posterior linearization filters achieve the best performance but at a substantially

D. Waxman et al./Sequential Inference for Gaussian Processes IPLF (converged) class 0 class 1

1.5

1.0

0.5

0.0

1.0

0.5

0.0

0.5

1.0

(a)

(b)

Online predictive log-likelihood

class 0 class 1

1.5

1.0 x2

PLF (1 sweep)

class 0 class 1

0.10 predictive log-likelihood

Laplace (RFF-GP classification) 1.5

0.15 0.20 0.25 0.30 0.35

Laplace PLF IPLF

0.40

(c)

100

150

train size n

200

250

(d)

Fig 5. Comparison of inference strategies for non-conjugate likelihoods in an RFF-GP classification task. (a)–(c) Posterior estimates obtained using the Laplace approximation, the posterior linearization filter (PLF), and the iterated posterior linearization filter (IPLF), respectively. (d) Predictive log-likelihoods over time, illustrating the improved accuracy of the IPLF relative to the other two methods, at the cost of higher computational complexity.

higher computational cost than both the standard posterior linearization filter and the Laplace approximation; they may also be more complex to implement. 4.4. A Bank of Filters Approach A persistent challenge in online learning lies in the optimization of hyperparameters. GPs, in particular, are highly sensitive to their kernel hyperparameters, and their marginal likelihood is generally non-convex and expensive to optimize. Moreover, efficient online algorithms for hyperparameter adaptation remain limited. As a result, approaches such as the ISSGP of Gijsberts and Metta (2013) often rely on pre-trained hyperparameters or costly periodic retraining. An alternative strategy is to maintain a collection of diverse models in the form of a kernel dictionary and perform online ensembling. This approach transforms the problem of hyperparameter adaptation into a more tractable online model-selection task, for which several effective algorithms exist. The first work to introduce ensembles of online RFF-GPs was that of Lu, Karanikolas and Giannakis (2022), which employed online Bayesian model averaging (O-BMA). O-BMA is an online and exact variant of the well-known Bayesian model averaging (BMA) framework (Hoeting et al., 1999), in which models are weighted according to their marginal likelihood. Specifically, BMA combines K probabilistic models given data D, denoted pk (yt | xt , D), for

D. Waxman et al./Sequential Inference for Gaussian Processes

k = 1, . . . , K, into a linear mixture, pBMA (yt | xt ) =

K X

wt,k pk (yt | xt , D),

(34)

k=1

where the weights wt,k are proportional to the model evidences, pk (D) wt,k = PK . j=1 pj (D)

(35)

In the online setting, D corresponds to the available observations (x1:t−1 , y1:t−1 ). The weights satisfy a recursive relationship that allows for exact recursive updates, wt+1,k ∝ wt,k pk (yt | x1:t , y1:t−1 ). (36) This property arises from the recursive factorization of the marginal likelihood, pk (D1:t ) = pk (D1:t−1 ) pk (yt | D1:t−1 ), which allows each model’s weight to be updated incrementally using only its previous value and the likelihood of the latest observation. From a Kalman filtering perspective, this corresponds directly to the wellknown bank of filters approach (Bar-Shalom, Li and Kirubarajan, 2004, Sec. 11.6) for online model combination. This connection suggests several potential improvements to the algorithms proposed in the literature, which have typically implemented naı̈ve O-BMA or simple Markov-switching weight updates. By contrast, the SP literature has demonstrated the advantages of more flexible switching and weighting mechanisms for combining filter outputs (El-Laham et al., 2021). A more recent alternative to O-BMA is online Bayesian stacking (OBS) (Waxman, Llorente and Djurić, 2026), which computes the weights wt,k by solving an online optimization problem. Unlike marginal-likelihood maximization, which is non-convex and computationally demanding, OBS seeks the set of weights that maximizes the marginal likelihood of the combined mixture, a convex objective that can be optimized efficiently using established online methods. Waxman, Llorente and Djurić (2026) demonstrate that OBS achieves favorable performance across several settings, including those involving the online optimization algorithms discussed above. 4.5. Extensions of Basis Expansions GPs The reduction of GPs to linear filters through basis-expansion approximations has proven highly useful for a broad range of SP tasks. In particular,

D. Waxman et al./Sequential Inference for Gaussian Processes

it enables the application of many established techniques from linear filtering theory, extending GPs to new settings. Examples include decentralized inference via consensus algorithms and the information form of the Kalman filter (Llorente, Waxman and Djurić, 2025), as well as robust and decentralized inference (Llorente et al., 2026), which leverages existing results from robust filtering theory (Chang and Li, 2017). The linear formulation of the resulting GP is also advantageous for analytical purposes. For instance, accessible bounds on the RFF-GP have led to new theoretical results in Bayesian optimization (Lu et al., 2023) and have been applied in online conformal inference (Xu, Lu and Giannakis, 2025). 5. Markovian GPs In the previous section, we showed that sequential estimation in GPs can be achieved by applying exact Kalman filtering to an approximate linear model. This approach relied on the assumption of a stationary kernel, which permits spectral representations of the GP. Remarkably, for many one-dimensional GPs with stationary kernels, referred to here as Markovian GPs, the process can be represented exactly as an SDE, which allows exact inference in linear time through Kalman filter and smoother formulations. 5.1. Gaussian Processes as State Space Models The central idea of the Markovian GP approach to sequential inference is to represent the GP as a continuous–discrete linear dynamical system. One of the simplest models is the Ornstein–Uhlenbeck (OU) process. This process is described by an SDE over a latent function f (t), which resembles an ordinary differential equation but is driven by a noise process W (t). In the OU case, the deterministic part of the SDE defines a mean-reverting dynamic: the function f (t) is continuously pulled toward zero at a rate λ, while the stochastic term introduces random fluctuations through the standard Wiener process W (t). The parameter λ controls the strength of mean reversion, and dW (t) provides a continuous-time analog of white Gaussian noise with variance scaled to be q. Noisy measurements ytn of the latent function are observed at discrete time points tn . Together, the OU process with discrete observations defines

D. Waxman et al./Sequential Inference for Gaussian Processes

f(t)

Online filter early Online filter later RTS smoother mid RTS smoother final

Fig 6. Equivalence of an OU–SDE prior and a GP with a Matérn– 12 kernel. We show filtering and smoothing solutions of the OU model alongside the reference GP posterior; illustrated are filtering solutions at an “early” observation and a “later” observation, and the smoothing solution at a “mid”-time solution and the “final” solution.

the continuous–discrete state–space model df (t) = −λf (t) dt + ytn = f (tn ) + εtn ,

√

q dW (t),

(37)

where q and λ are hyperparameters, and εtn represents white Gaussian observation noise. It can be shown that this process corresponds exactly to a GP prior with a Matérn– 21 kernel when λ = ℓ−1 and σf2 = q/(2λ), given a prior f (0) ∼ N (0, σf2 ). Using this representation yields a surprisingly efficient algorithm for GP regression via Kalman filtering. Given data D = {(tn , yn )}N n=1 , we first discretize the OU state–space model over irregular time steps ∆tn = tn − tn−1 , ftn = e−λ∆tn ftn−1 + ηtn ytn = f (tn ) + εtn , n = 1, 2, . . . where ηtn ∼ N 0, Qtn with Qtn =

q 1 − e−2λ∆tn = σf2 1 − e−2λ∆tn , 2λ

(38)

(39)

and εn ∼ N (0, σε2 ). We then apply Kalman filtering and, if desired, Rauch–Tung–Striebel (RTS) smoothing to obtain exact GP inference. For a Markovian GP with state dimension d, the total cost is O(N d3 ). Thus, for fixed (small) d, the complexity is linear in N . This SDE–GP duality naturally raises a question: can other GP priors be represented as linear time-invariant (LTI) SDEs to enable linear-in-N exact

D. Waxman et al./Sequential Inference for Gaussian Processes

inference? Although the OU process alone produces rough sample paths, the answer is yes by augmenting the state with derivatives. As a first family, half-integer Matérn kernels admit exact linear-SDE representations. For Matérn–ν with half-integer ν, we define the augmented state to include f and its first ν − 21 derivatives, ⊤ f̃ (t) = f (t) f (1) (t) · · · f (ν−1/2) (t) . The driving noise acts only on the highest derivative, with lower orders coupled by a polynomial in λ. (See, e.g., Särkkä and Solin (2019).) For the Matérn– 32 kernel, one convenient continuous–discrete form is (Särkkä and Solin, 2019, Ex. 12.7) df̃ (t) = F f̃ (t) dt + L

√

q dW (t),

(40)

ytn = H f̃ (tn ) + εtn , with

0 1 F= , −λ2 −2λ

0 L= , 1

H= 1 0 ,

q = 4λ3 σf2 ,

and εtn ∼ N (0, σε2 ). Kalman filtering/smoothing on this two-dimensional state recovers the exact GP posterior for Matérn– 23 . A natural question is: how expressive are the GPs that can be represented in this way? Perhaps surprisingly, they are arbitrarily expressive within the class of stationary GPs. Recent work (Dowling, Sokól and Park, 2021; Loper et al., 2021) shows that a GP with any stationary kernel can be approximated arbitrarily well by a simple modification of the Matérn kernel. Following Dowling, Sokól and Park (2021), we refer to these as mixture Hida–Matérn (MHM) kernels, although they are also known as latent exponentially generated processes in Loper et al. (2021). To construct the MHM family, we begin with the Hida–Matérn (HM) kernel, which combines a phase shift b with a Matérn–ν kernel: κHM(b,ν) (t, t′ ) = cos b |t − t′ | κMat-ν (t, t′ ), (41) where all hyperparameters are left implicit. We refer to b as a phase shift because it corresponds to a shift by b in the PSD of the kernel relative to the underlying Matérn kernel. The MHM family is defined as linear mixtures of

D. Waxman et al./Sequential Inference for Gaussian Processes

HM kernels, which provide a Matérn-based analog of the spectral mixture kernel (Wilson and Adams, 2013), with the intuitive goal to perform density estimation of the PSD. The resulting state-space formulation is quite similar to the Matérn state space above. By the linearity of GPs, the sum of kernels corresponds to the sum of independent GPs. Thus, for a simple two-component MHM kernel, κ(t, t′ ) = w1 cos(b1 |t − t′ |) κMat-ν1 (t, t′ ) + w2 cos(b2 |t − t′ |) κMat-ν2 (t, t′ ), (42) the corresponding state–space representation is # " √ F1 0 f̃ 1 (t) q1 L 1 df̃ (t) = dt + √ dW(t), q2 L 2 0 F2 f̃ 2 (t) √ √ w1 H1 w2 H2 f̃ (tn ) + εtn . ytn =

(43)

Further generalizations follow in the same way and yield block-diagonal covariance matrices that can be computed efficiently and exploited for more efficient filtering and smoothing. 5.2. Spatiotemporal Regression via Markovian Gaussian Processes Although Markovian GPs were originally formulated for temporal processes, they can be naturally extended to spatiotemporal settings. The mathematical strategy remains similar, exploiting the spectral representation of the spatiotemporal kernel to derive an appropriate stochastic partial differential equation (SPDE) (Sarkka and Hartikainen, 2012). This approach shares deep connections with more general work on the SPDE formalism of GPs (Lindgren, Bolin and Rue, 2022), with the infinite-dimensional filtering presented by Sarkka and Hartikainen (2012) being an efficient implementation for simple spatial fields (cf. Section 4.3 of Lindgren, Bolin and Rue (2022)). This general formulation, however, is mathematically involved and typically requires advanced approximate inference methods for SPDEs. In this tutorial, we focus instead on a simplified and practically useful case that involves separable kernels. Recall that a spatiotemporal kernel κ((x, t), (x′ , t′ )) is separable if it can be expressed as the product of a spatial kernel and a temporal kernel, κs (x, x′ ) κt (t, t′ ). This assumption is convenient because it preserves the Markov property: spatial correlations act only

D. Waxman et al./Sequential Inference for Gaussian Processes

through a filtered and augmented spatial state f̃ s . The remaining task is to define f̃ s and specify how spatial dependencies are incorporated into the overall model. In the general case, f̃ s is a functional, that is, an infinite-dimensional spatial state evolving over time, on which infinite-dimensional Kalman filtering and smoothing can, in principle, be applied (Sarkka and Hartikainen, 2012). While this formulation is theoretically elegant and occasionally practical, we again adopt a simpler approach. Our final simplification for spatiotemporal GPs is to discretize the spatial domain into a finite set of locations x1 , . . . , xNs . The resulting state is then defined by Ns copies of the temporal kernel state f̃ , i.e., i⊤ h ⊤ ⊤ (44) f̃ S = f̃ 1 · · · f̃ Ns . The associated SDE has dynamic and feedback matrices formed as the blockdiagonal concatenation of F and L. Spatial dependencies enter through the diffusion term, whose covariance takes the Kronecker-structured form Qt ⊗ ΣSS , where ⊗ denotes the Kronecker product, and ΣSS is the spatial covariance matrix induced by the spatial kernel. Assuming that the temporal kernel induces a dt -dimensional state space (using dt to differentiate between dimensions induced by the temporal kernel and spatial kernel) , the resulting spatiotemporal model has a (dt × Ns )dimensional joint state. Inference through the Kalman filter and smoother then requires O(d3t Ns3 T ) computations, which yield quadratic savings in T compared to standard GP inference. In the following, we discuss several extensions that further reduce the cubic scaling with respect to the number of spatial locations Ns . 5.3. Extensions of Markovian Gaussian Processes As with basis-expansion approximations, much of the recent progress in extending Markovian GPs arises from their convenient formulation as linear–Gaussian systems, which allows the direct use of tools from SP. A primary challenge in this setting involves non-Gaussian likelihoods. Analogous to the basis-expansion case (cf. Fig. 5), standard nonlinear filtering and smoothing methods can be applied, such as posterior linearization (Garcı́aFernández, Tronarp and Särkkä, 2019). The widespread use of GP classification has even inspired the development of new approximate inference methods

D. Waxman et al./Sequential Inference for Gaussian Processes

for state-space models, such as state-space expectation propagation (Wilkinson et al., 2020). Variational inference methods are also effective and represent the state of the art in large-scale spatiotemporal applications (Hamelijnck et al., 2021) and in modeling count data, including neuroscience applications (Dowling, Zhao and Park, 2023). Another major challenge for Markovian GPs, especially in spatiotemporal settings, is scalability. Although Markovian formulations avoid the cubic scaling in time typical of standard GPs, they still exhibit O(Ns3 ) complexity in the number of spatial locations, which can become prohibitive. One approach to avoid this is the spatiotemporal sparse variational GP, which introduces a grid of Ms fixed spatial inducing points and achieves O(Ms2 Ns T ) training and inference costs (Hamelijnck et al., 2021). While such models are not sequential in the strict sense of being recursive or online, they still perform inference sequentially in time and thus retain favorable scaling properties. Markovian GPs have also motivated new computation-aware filtering and smoothing algorithms (Pförtner et al., 2025), which leverage probabilistic linear algebra for scalable inference. Finally, efforts to make Markovian GPs robust have benefited directly from the robust filtering literature. Simple approaches include rejection filtering (Bock et al., 2022; Waxman and Djurić, 2024b), while more recent methods employ generalized Bayesian inference to achieve efficient and outlierrobust spatiotemporal GPs (Laplante et al., 2025). We will further discuss robustness, and its implications for SP, in the following section. 6. Sequential Estimation in Sparse (Variational) Approximations Exact GP inference remains computationally expensive for large datasets. Sparse and variational approximations overcome this limitation by representing the latent function through a compact set of inducing variables that summarize the information from past data. In sequential settings, these lowdimensional representations can be updated recursively to allow for efficient online inference without revisiting earlier observations. When data arrive sequentially, these sparse representations can be updated online rather than recomputed from scratch. This section introduces the main ideas of sequential inference within sparse and sparse-variational GPs. We show how the algebra of Gaussian conditioning allows recursive updates of the inducing-point statistics, discuss connections to classical filtering concepts, and describe modern variational formulations that preserve scalability in both

D. Waxman et al./Sequential Inference for Gaussian Processes

time and data size. 6.1. Sequential Estimation in Sparse Approximations The classical sparse GP approximation assumes that the latent function values f = [f (x1 ), f (x2 ), . . . , f (xN )]⊤

(45)

depend on a much smaller set of inducing variables u = [f (xu1 ), f (xu2 ), . . . , f (xuM )]⊤ ,

M ≪ N,

located at inducing inputs xu = {xu1 , . . . , xuM }. The joint prior is 0 Σf f Σf u p(f , u) = N , , 0 Σuf Σuu

(46)

(47)

where Σuu = [κ(xui , xuj )]i,j is the covariance among inducing inputs, Σf u = [κ(xi , xuj )]i,j is the cross-covariance between data and inducing inputs, and Σf f is the full covariance of the observations. Conditioning on u yields p(f | u) = N (Σf u Σ−1 uu u, Qf f ),

(48)

where Qf f = Σf f − Σf u Σ−1 uu Σuf represents the residual covariance. All dependence on the data enters only through u, which becomes the effective state of the approximation. This greatly reduces the computational burden of prediction, as we require the inverse of an M × M matrix, rather than an N × N matrix, where M ≪ N . When data arrive sequentially, the posterior over u can be updated recursively. Suppose at time t − 1 we have pt−1 (u) = N (mt−1 , St−1 ).

(49)

The new observation (xt , yt ) introduces a likelihood p(yt | ft ) = N (yt ; ft , σε2 ),

(50)

−1 and the marginal ft is linearly related to u through ft = k⊤ t Σuu u + εt , where kt = [κ(xt , xu1 ), . . . , κ(xt , xuM )]⊤ . Using Gaussian conditioning identities,

D. Waxman et al./Sequential Inference for Gaussian Processes

the posterior mean and covariance of u follow 2 −1 Kt = St−1 ht (h⊤ t St−1 ht + σε ) ,

(51)

mt = mt−1 + Kt (yt − h⊤ t mt−1 ),

(52)

St = St−1 − Kt h⊤ t St−1 ,

(53)

where ht = Σ−1 uu kt . These updates are mathematically identical to those of a Kalman filter with state u. Each step requires O(M 2 ) computation, which enables online inference without the need to revisit previous observations. This correspondence with the Kalman filter provides an intuitive interpretation: the inducing-point vector u serves as the hidden state, and each new observation performs a measurement update on that state. The matrices Σuu and kt act as prior and observation covariances, respectively. Under this viewpoint, the sparse GP operates as a linear–Gaussian dynamical system with a trivial transition model ut = ut−1 . This perspective highlights both the strengths and limitations of the approach. Updates are exact within the approximation, but the inducing inputs locations xu remain fixed. They do not adapt to new data unless explicitly optimized. Early placement of xu therefore heavily influences long-term performance. Adaptive selection or relocation of inducing inputs can improve accuracy, though at the cost of additional optimization. Sparse approximations differ from subset-of-data (SoD) or subset-of-regressor (SoR) approaches. In subset methods, only a portion of past observations contributes to future predictions, whereas in sparse GPs the inducing inputs interpolate the full posterior. Sequential updates thus preserve the information content of all past data rather than discarding it, which makes sparse GPs well suited for long-term estimation tasks in slowly varying systems. 6.2. Sequential Estimation in Sparse Variational Approximations The probabilistic formulation above assumes exact Gaussian identities. In practice, the marginal likelihood p(y) is rarely tractable, especially for nonGaussian likelihoods p(yt | ft ). Variational inference introduces a principled way to approximate this posterior while maintaining scalability. In the variational sparse GP (VSGP), we approximate the true posterior p(u | y) by a Gaussian variational distribution q(u) = N (m, S).

(54)

D. Waxman et al./Sequential Inference for Gaussian Processes

True

Posterior mean

Data (current/prev)

Inducing mean (current/prev)

Inducing ±2σ (current/prev)

Fig 7. An illustration of the streaming VSGP (via the implementation provided by Chang et al. (2023)): as new data arrive (black) the GP approximation adapts online, distributing the inducing locations (light and darker yellow) along the interval in order to produce an accurate approximation.

The evidence lower bound (ELBO) is log p(y) ≥

T X

Eq(ft ) [log p(yt | ft )] − KL q(u) ∥ p(u) ,

(55)

t=1

where q(ft ) denotes the marginal of ft under q(u). Optimizing this bound yields both the posterior parameters (m, S) and the kernel hyperparameters. Sequential inference follows by decomposing the ELBO over time. When a new observation (xt , yt ) arrives, we add a new term to the ELBO and update (m, S) using Kalman filter-like recursion (Schürch et al., 2020), i.e., −1 ⊤ −1 S−1 t = St−1 + Ht Rt Ht ,

(56)

−1 ⊤ −1 S−1 t mt = St−1 mt−1 + Ht Rt yt ,

(57)

where Ht and Rt are determined by the local likelihood model. These relations resemble recursive least squares but arise from the variational objective and allow online updates of (m, S) without revisiting past data. In the streaming VSGP (Bui, Nguyen and Turner, 2017), the current posterior qt−1 (u) acts as the prior for the next step, and the new posterior qt (u) minimizes p(yt | u) qt−1 (u) , (58) KL(qt (u) ∥ p(u | D1:t )) ≈ KL qt (u) Zt

D. Waxman et al./Sequential Inference for Gaussian Processes

where Zt is a normalization constant. This formulation enables efficient onepass updates of the variational parameters. In the simpler scenario, the inducing inputs xu remain fixed (along with the hyperparameters) and maintain constant memory even as the data stream grows. Bui, Nguyen and Turner (2017) formalize the intuition in Eq. (58) but consider the more challenging problem of online updating for the inducing inputs and hyperparameters. Denoting with θ the set of hyperparameters, they consider the intermediate t |f ) qt−1 (f ) , hence target density to be instead pb(f |Dt−1 , yt , θt−1 , θt ) ∝ p(f |θt ) p(yp(f |θt−1 ) effectively accounting for the new θt in the qt (f ), with possibly different inducing locations xu,t ̸= xu,t−1 . An instance of the streaming VSGP by Bui, Nguyen and Turner (2017) is shown in Figure 7. Chang et al. (2023) propose a memory-based dual-SVGP that augments the inducing-point posterior summary with a small, actively selected buffer of past data, helping sequential conditioning and hyperparameter/inducing-input adaptation track the offline ELBO more closely. The streaming VSGP also connects to natural-gradient descent in the space of variational parameters (Hamelijnck et al., 2021; Khan and Lin, 2017). Because q(u) is Gaussian, its natural parameters are (S−1 m, S−1 ), and small updates along these directions correspond to moment matching. This relationship links online variational inference to recursive information filters in classical SP (Opper and Archambeau, 2009). Non-Gaussian likelihoods and Robustness The variational formulation easily extends to non-Gaussian likelihoods. The expectation Eq(ft ) [log p(yt | ft )] occasionally has strong approximations available (e.g., in the case of logit or probit models), or numerically otherwise. Sequential updates depend only on the local gradients of this expectation. Robust likelihoods, such as Student-t or Huber forms, reduce sensitivity to outliers without altering the recursive structure. A further advantage of the variational approach is its modularity. It can combine with structured likelihoods, multi-output GPs, or hierarchical priors without changing the recursive pattern. This flexibility has made online variational GPs the preferred method in complex spatiotemporal models that demand both scalability and adaptability. The per-step complexity of sequential variational updates is O(M 3 ) due to matrix inversions in Σuu or S. Because M is fixed and typically small, the total cost grows linearly with time. Mini-batches of size B reduce this to O(BM 2 ) per update. Only the parameters (m, S) and the kernel matrices

D. Waxman et al./Sequential Inference for Gaussian Processes

associated with xu are stored, so memory remains constant, unlike the O(N 2 ) storage requirement of exact GPs. The algebraic similarity between online variational GPs and Kalman filtering is not accidental. Both maintain a Gaussian belief that is updated through the integration of new information. The difference lies in interpretation: Kalman filtering assumes a known generative model, whereas variational inference optimizes a lower bound that approximates it. The update of (m, S) parallels the correction step of an information filter, with S−1 acting as the information matrix. This analogy helps build intuition: we may view the variational updates as approximate filtering, where the innovation term is replaced by a variational expectation. This viewpoint unifies probabilistic inference and adaptive SP within a common mathematical framework. A few practical comments are in order. In sequential applications, numerical conditioning is very important. Since updates of S involve inversion and accumulation, round-off errors can degrade positive definiteness. Cholesky-based or square-root parameterizations maintain stability. Furthermore, adaptive placement of inducing inputs can improve accuracy in nonstationary domains. The initialization of xu also affects convergence: k-means clustering or Latin hypercube sampling over the input domain provide effective starting points. 7. Other Approaches to Sequential Estimation While the previous sections focused on sparse, variational, and Markovian formulations, several alternative approaches to sequential estimation have been developed in the GP literature and in related signal-processing frameworks. These methods attempt to reconcile GP expressiveness with the computational and memory constraints of real-time applications. They include recursive kernel methods, variational state-space formulations, online expectation propagation, and combinations of GPs with neural or particle-based estimators. Although conceptually diverse, they share a common goal: to maintain a compact belief state that summarizes past information while adapting efficiently to new data.

D. Waxman et al./Sequential Inference for Gaussian Processes

7.1. Recursive Kernel and Online Regression Methods Early work on sequential GP inference originated in the signal-processing community through kernel adaptive filters (Engel, Mannor and Meir, 2004). The kernel recursive least-squares (KRLS) algorithm and its variants (Van Vaerenbergh, Lázaro-Gredilla and Santamarı́a, 2012) replace explicit GP conditioning with recursive updates of kernel coefficients. Each new observation introduced a basis function, and sparsification criteria determined whether it should be retained in the “dictionary.” The approximate posterior mean took the form of a weighted sum of these dictionary elements, with the corresponding covariance tracked implicitly through the weight update rule. Although KRLS does not maintain a full probabilistic representation, its updates are algebraically equivalent to those of a GP with a dynamically pruned set of inducing points. Parallel developments in the signal-processing community established a probabilistic foundation for these methods (Pérez-Cruz et al., 2013), leading to the kernel Kalman filter and GP adaptive filters that explicitly propagated the mean and covariance of the function values . These methods provided the first bridge between adaptive filtering and GP regression, and they remain competitive in real-time nonlinear tracking and denoising tasks. A practical distinction between kernel adaptive filters and modern sparse GPs lies in their treatment of uncertainty. Kernel filters rely on deterministic or heuristic forgetting factors, whereas variational or state-space GPs retain explicit probabilistic uncertainty through Gaussian posteriors. Nevertheless, the recursive kernel viewpoint remains influential, particularly in systems that require rapid adaptation with limited computation, such as streaming sensor networks or embedded platforms. 7.2. Variational State–Space and Hybrid Models A second class of methods merges GP priors with explicit state-space formulations where the transition or observation functions are learned via neural networks. In these hybrid models, latent dynamics evolve according to a GP prior, but inference proceeds through variational filtering rather than exact Kalman updates. The variational Kalman GP and related frameworks (Dowling, Zhao and Park, 2023; Hamelijnck et al., 2021) maintain a structured Gaussian posterior over both latent states and inducing variables. Sequential updates of the variational parameters proceed by minimizing the local

D. Waxman et al./Sequential Inference for Gaussian Processes

Kullback–Leibler divergence between the approximate and true filtering distributions. This approach inherits the expressive power of non-Markovian kernels while preserving linear-time updates in the temporal dimension. A practical advantage of these hybrid models is their ability to incorporate non-Gaussian likelihoods, latent regime switches, and heteroscedastic noise without resorting to particle approximations. They have become standard tools in modern spatiotemporal modeling, particularly in applications involving neural activity (Dowling, Zhao and Park, 2023) and environmental monitoring (Hamelijnck et al., 2021), among others. The combination of GP priors with amortized variational inference, implemented through recurrent neural networks, further extends their reach to high-dimensional or irregularly sampled data streams (Fortuin et al., 2020). 7.3. Expectation Propagation and Message–Passing Frameworks Another family of sequential GP algorithms derives from expectation propagation (EP). In the EP formulation, each data point contributes a site approximation that updates the global posterior through moment matching. Sequential EP updates can be written in recursive form, qt (f ) ∝ qt−1 (f )

p̃t (ft ) , Zt

(59)

where p̃t (ft ) denotes the local Gaussian site and Zt a normalization factor. When expressed in natural parameters, this recursion is identical to an information-form Kalman filter with adaptive site precision. The state-space expectation propagation (SSEP) method (Wilkinson et al., 2020) generalizes this idea to temporal GPs and achieves linear-time exactness for Markovian kernels and accurate approximations for more general covariances. Because each update corresponds to a local moment-matching step, SSEP preserves the calibration of posterior variances, a property often lost in purely variational schemes. Message-passing perspectives also provide a unifying view of GP inference. From this perspective, the predictive mean and covariance evolve through forward and backward messages that obey the same algebra as Kalman filtering but operate in function space. This view clarifies relationships among variational methods and EP, and it supports extensions to hierarchical or deep GP architectures.

D. Waxman et al./Sequential Inference for Gaussian Processes

7.4. Monte Carlo and Particle Approaches In strongly nonlinear or non-Gaussian settings, sequential estimation may rely on Monte Carlo techniques. GP state-space models — state-space models where the transition and/or observation models are modeled by a GP — are a particularly challenging setting, necessitating Monte Carlo methods. In this setting, GP particle filters approximate the joint distribution over latent states and function values by a set of weighted particles (Ko and Fox, 2009). Each particle carries its own GP belief, which evolves through resampling and weight updates based on predictive likelihoods. Although computationally expensive, these methods enable fully nonlinear and non-Gaussian inference with GP priors and have been applied to robotics and control problems where model uncertainty must be tracked explicitly (Ko and Fox, 2009). Recent work has further incorporated alternative filtering methods into the GP SSM problem. Namely, Lin et al. (2024) use the ensemble Kalman filter for efficient variational inference in GP SSM models. 8. Comparisons to Deep Learning Models GP models and deep learning represent two complementary paradigms for modeling sequential data. Both attempt to model dependencies across time, uncertainty in prediction, and generalization to unseen conditions. Yet their philosophical foundations differ. GPs permit two complementary views. In the function–space view, the kernel defines a prior on functions, and Bayes’ rule yields a posterior over functions given observed data. In the parameter–space view, finite–basis approximations such as RFF or HSGP represent the GP through a fixed set of coefficients that evolve sequentially with constant memory. Related low–rank functional approximations, such as variational inducing–point GPs, maintain the same efficiency by summarizing past information in a finite set of inducing variables. On the other hand, deep learning constructs deterministic parametric mappings optimized for predictive accuracy. This section contrasts the sequential GP frameworks developed earlier with several major classes of deep learning methods, including recurrent and convolutional architectures, diffusion processes, and deep kernel models. We emphasize similarities in structure, differences in inductive bias, and recent unifying viewpoints that link probabilistic and neural approaches.

D. Waxman et al./Sequential Inference for Gaussian Processes

8.1. Conceptual Comparison: Function–Space versus Parameter–Space Inference Exact GPs define inference directly in function space, where the kernel specifies a prior over functions, and conditioning on observations produces a posterior distribution over all admissible functions. We discussed sequentially tractable approximations such as RFF, HSGP, and variational inducingpoint models, which retain the same probabilistic interpretation but express the posterior through a finite set of parameters or inducing variables. Thus, these models allow for constant-memory recursive updates with calibrated uncertainty. Deep neural networks (DNNs), by contrast, define explicit parametric maps fθ (xt ) with parameters θ updated in parameter space. Learning in this setting typically involves stochastic gradient descent, which yields a single point estimate of θ rather than a full distribution. Uncertainty estimates are therefore approximate and are usually obtained through ensembles, dropout sampling, or stochastic variational layers. The contrast between fully functional inference and its finite-dimensional approximations determines how models balance data efficiency and robustness. Exact function-space inference relies on strong priors and remains sampleefficient, whereas parameter-space formulations such as RFF, HSGP, or variational inducing-point GPs trade some prior fidelity for scalability and numerical stability. Furthermore, GPs excel in low–data regimes, where priors constrain functional variability. By contrast, DNNs tend to require larger data volumes to learn similar inductive biases but achieve greater flexibility in high–dimensional domains. Sequential GP algorithms occupy a middle ground between analytical tractability and dynamic adaptability. They maintain explicit uncertainty through recursive updates while scaling linearly in time. 8.2. Diffusion Processes and the Markovian View A particularly interesting point of comparison arises with diffusion–based models, which have become a dominant class in modern deep learning (Yang et al., 2023). Diffusion models define a forward stochastic process that gradually corrupts data with Gaussian noise and a reverse process that learns to denoise and recover the original distribution. Both directions form Markov processes in continuous time and thus bear formal similarity to the Marko-

D. Waxman et al./Sequential Inference for Gaussian Processes

vian GP models discussed in Section V. Each defines an SDE over latent trajectories, but their objectives differ: diffusion models learn (intractable) reverse dynamics by optimizing an ELBO term on the data likelihood, whereas Markovian GPs use known forward and reverse dynamics of a simpler, linear diffusion model to perform exact Bayesian inference. Recent work has begun to connect these two frameworks explicitly. Verma, Adam and Solin (2024) build upon the work of Archambeau et al. (2007) by examining the relationship between general diffusion models (now commonplace in ML) and the LTI SDEs (described by Markovian GPs) through a variational framework. Efficient variational algorithms are developed, and they compare favorably to other variational algorithms on SDEs (such as neural SDEs). The comparison between diffusion processes and Markovian GPs is particularly interesting in that it contrasts tractable inference under a defined stochastic-process prior with more data-driven, generative deep learning models with similar mathematical foundations. 8.3. Sequential Deep Learning Architectures Sequential data modeling in deep learning primarily relies on recurrent neural networks (RNNs), long short–term memory networks (LSTMs), temporal convolutional networks (TCNs), and, more recently, transformers. These architectures propagate information through hidden states that are updated deterministically by learned parameters. The hidden state of an RNN plays a role analogous to the latent state in a GP state–space model, but its update follows a nonlinear map rather than a probabilistic transition. This distinction affects both interpretability and uncertainty quantification. For instance, the basis expansion GP defines a linear Gaussian SSM, θ t+1 = θ t + εt ,

yt = ht θ t + ηt ,

(60)

where εt and ηt are Gaussian noise terms. The resulting filtering solution are nonlinear updates of ht θt and yt . An RNN replaces these stochastic transitions with a nonlinear function ft+1 = ϕ(W ft + U yt ), where ϕ is a learned function and ft are a set of state variables. While this deterministic formulation achieves strong empirical performance, it lacks the uncertainty calibration inherent in GP updates. Hybrid designs, such as variational RNNs and deep state–space models, reintroduce latent stochasticity and thus approximate the recursive inference mechanisms of sequential GPs.

D. Waxman et al./Sequential Inference for Gaussian Processes

Transformers extend this idea through self–attention, which models all pairwise dependencies between time steps (Vaswani et al., 2017). From a GP perspective, self–attention can be interpreted as defining a learned, data–dependent kernel. Indeed, several recent analyses interpret attention weights as low–rank approximations to a GP covariance matrix (e.g., Bui et al. (2025)). Sequential GPs with adaptive kernels, therefore, share structural parallels with transformer models but maintain an explicit probabilistic interpretation. 8.4. Deep Kernels and Hierarchical Gaussian Processes Another natural point of contact between deep learning and GPs lies in kernel design. Traditional GPs rely on hand–crafted kernels such as the squared–exponential or Matérn family, which impose fixed notions of similarity in input space. Deep kernel learning (DKL) (Wilson et al., 2016) replaces this assumption by parameterizing the kernel through a neural feature map, κ(x, x′ ) = κ0 (ψθ (x), ψθ (x′ )) ,

(61)

where ψθ is a deep network and κ0 is a base kernel such as the RBF. This composition allows the GP to operate in a learned latent representation that captures the hierarchical structure in the data while retaining Bayesian inference in the output space. Sequential DKL variants integrate temporal recurrence into the embedding process, exploiting the sequential nature of data in the embedding prior (Al-Shedivat et al., 2017). An alternative approach are deep Gaussian processes (DGPs) (Damianou and Lawrence, 2013), which stack multiple GPs hierarchically, with each layter defining a random function that feeds into the next. The resulting process represents a nonparametric analog of a DNN, with explicit uncertainty propagation across layers. Although exact inference in DGPs is intractable, variational approximations yield tractable sequential updates analogous to those described in Section VI. From a signal–processing perspective, these models resemble cascaded dynamical systems whose intermediate outputs remain probabilistically defined. Versions of the DGP based on the Markovian GP are proposed in Zhao, Emzir and Särkkä (2021). 8.5. Nonstationary and Adaptive Representations The intersection of deep learning and GP modeling has also led to new perspectives on nonstationarity. Titsias et al. (2024) provide an online neural

D. Waxman et al./Sequential Inference for Gaussian Processes

network learning solution based on a method akin to a combination of DKL and the Kalman-filter-based form of linear regression. This is equivalent to taking neural network-based embeddings and a linear kernel in the basis expansion GPs previously outlined. This approach connects neural adaptive filtering to filtering-based Bayesian models. Using the same tricks as the non-stationary GP regression discussed earlier, the model adapts dynamically to nonstationary environments while maintaining analytic tractability in its posterior updates. These developments suggest a continuum between purely kernel–based and purely neural approaches. At one end, classical GPs provide closed–form sequential inference with strong uncertainty calibration but limited representational flexibility. At the other end, deep networks provide extreme flexibility but weak uncertainty control. Hybrid methods, including DKL, DGP, and neural–SDE formulations, offer a balance by embedding learnable representations within probabilistic filters. 8.6. Outlook: Complementarity and Integration From the viewpoint of SP, sequential GPs and deep learning differ not as competing paradigms but as complementary tools. GPs serve as interpretable, uncertainty–aware priors that can guide learning and quantify confidence, whereas deep networks supply expressive, data–driven parameterizations of complex nonlinear dependencies. Recent work increasingly combines the two, with neural networks that amortize kernel computation, approximate marginal likelihoods, or perform variational updates. Conversely, GP perspectives have improved the calibration and sample efficiency of deep generative and sequential models. In practice, hybrid architectures that integrate GP priors with deep representations already achieve state–of–the–art performance in time–series forecasting, robotics, and spatiotemporal modeling. They preserve the recursive update structure central to SP while drawing on the representational power of deep learning. This convergence underscores a broader trend: probabilistic and neural methods are not distinct domains but parts of a unified framework for sequential inference under uncertainty.

D. Waxman et al./Sequential Inference for Gaussian Processes

9. Applications of Gaussian Processes to Sequential Signal Processing The methods discussed in the preceding sections—basis expansions, Markovian representations, and variational approximations—enable GP models to operate efficiently in sequential environments. In this section, we explore how these formulations translate into practical SP applications. We organize the discussion around several representative domains that highlight different aspects of GP inference: sequential regression and forecasting, distributed and decentralized inference, robust learning and anomaly detection, sequential Bayesian optimization, and adaptive sensing. Across these areas, GPs serve as a unifying probabilistic framework for uncertainty–aware estimation and decision making. 9.1. Sequential Regression and Forecasting Sequential regression remains one of the canonical applications of GPs in SP. In this setting, the objective is to infer a latent, possibly time–varying signal from an incoming data stream while maintaining calibrated uncertainty estimates. We saw that when combined with RFF or Hilbert–space approximations, GP regression reduces to a linear–Gaussian model that allows recursive updates through the Kalman filter. This formulation allows real–time prediction, filtering, and smoothing with scalable memory usage and linear time complexity in the number of observations. Recent works demonstrate the flexibility of this approach. Lu, Karanikolas and Giannakis (2022) proposed incremental ensemble GPs that adapt model weights online through Bayesian model averaging. Waxman and Djurić (2024a) extended this concept to more general basis expansions. Other approaches focus on variational learning; Chang et al. (2023) build upon the dual GP of Adam et al. (2021) to include memory of past points and derive efficient sequential updates. 9.2. Distributed and Decentralized Online Inference In multi–agent and networked systems, sequential GPs provide a principled mechanism for decentralized learning. Each node maintains a local GP model that encodes partial information about the global process, while message

D. Waxman et al./Sequential Inference for Gaussian Processes

passing or consensus updates ensure global coherence. Because the GP posterior is Gaussian, local summaries can be transmitted as information quantities — these are equivalent to transmitting means and covariances, but decompose additively with more data — which allows for communicationefficient fusion of information. Llorente, Waxman and Djurić (2025) recently introduced decentralized online ensembles of GPs, which employ the information form of the Kalman filter to propagate local sufficient statistics. The framework extends classical results from distributed SP to probabilistic nonparametric settings and allows agents to maintain uncertainty–aware beliefs without centralized coordination. Related work on robust decentralized GPs (Llorente et al., 2026) incorporates outlier resistance and adaptive weighting and draws from robust filtering theory (Chang and Li, 2017). Such approaches illustrate how sequential GP inference can support scalable, fault–tolerant estimation in sensor networks and multi–robot systems. 9.3. Robust Sequential Learning and Anomaly Detection Sequential GP models are also effective for detecting anomalies or regime changes in streaming data. Because GPs maintain predictive variances that reflect uncertainty about future observations, deviations from expected predictive likelihoods naturally indicate potential outliers or distributional shifts. Markovian and basis–expansion GPs both accommodate such analysis by tracking evolving latent states and adapting their uncertainty in real time. Bock et al. (2022) used state–space GPs to perform online anomaly detection and robust inference in time series, leveraging predictive residuals as test statistics for unexpected events. Waxman and Djurić (2024b) proposed a GP–based streaming algorithm for regime switching and outlier prediction, which use Markovian GPs for efficient updates. More recently, Laplante et al. (2025) introduced robust spatiotemporal GPs based on generalized Bayesian inference, which avoid the influence of heavy–tailed noise and unmodeled dynamics. Laplante et al. explicitly note that the sequential formulation helps resolve open problems in robust GP estimation, by allowing more informed weighting of new observations. These developments connect GP–based filtering to long–standing interests in the SP community, namely resilience, robustness, and principled uncertainty quantification.

D. Waxman et al./Sequential Inference for Gaussian Processes

9.4. Sequential Bayesian Optimization GPs are widely used in Bayesian optimization due to their ability to quantify uncertainty and guide exploration. In sequential settings, these methods adapt the acquisition function as new data arrive, often under real–time or resource–constrained conditions. Sequential GP optimization has proven valuable in adaptive control, communication systems, and autonomous experimentation. Lu et al. (2023) proposed surrogate–model ensembles that update their posterior in real time through RFF–based approximations and achieve scalable optimization beyond a single GP. Maddox, Stanton and Wilson (2021) introduced conditioning techniques that stabilize GP posteriors for streaming optimization tasks while maintaining uncertainty calibration. By combining approximate inference with efficient online acquisition strategies, these methods preserve the statistical rigor of Bayesian optimization while remaining computationally feasible for embedded or interactive systems. 9.5. Adaptive and Active Sensing Sequential GP inference naturally lends itself to adaptive sensing, where the goal is to select measurements that maximize information gain or minimize uncertainty about an evolving process. In this context, the GP posterior defines an acquisition criterion that guides sensor placement or measurement timing. Because GP updates are recursive, sensor trajectories or sampling schedules can be optimized in real time. Ziatdinov et al. (2022) demonstrated this idea for atomic–scale imaging, where Bayesian GP frameworks adaptively guide scanning probe measurements to capture dynamic physical phenomena. In spatiotemporal SP, similar strategies enable mobile sensor networks to track dynamic fields while minimizing energy consumption and communication costs. The probabilistic nature of GPs ensures that sensing decisions are informed by both estimated signal structure and model uncertainty, and that they align naturally with long-standing principles of active and adaptive filtering. 10. Conclusion GPs offer a principled probabilistic framework that unites many classical ideas of SP under a single mathematical formalism. Throughout this article, we have presented how sequential and online inference in GPs connects

D. Waxman et al./Sequential Inference for Gaussian Processes

directly to established concepts such as recursive least squares, Kalman filtering, and SDEs. By introducing finite basis expansions, Markovian representations, and variational approximations, we have shown that the apparent computational limitations of GPs can be overcome while preserving their ability to quantify uncertainty. These formulations restore the recursive structure fundamental to SP and make GPs suitable for real-time and large-scale applications. The equivalence between linear-Gaussian state–space models and GP priors has clarified the role of stochastic dynamics in kernel design. Spectral and Hilbert-space approximations translate nonparametric GP priors into compact linear models that evolve sequentially through Kalman updates. Markovian GPs extend this correspondence to continuous time and yield exact inference in linear time for a wide class of stationary kernels. Sparse and variational formulations summarize the information contained in past observations through inducing variables or basis coefficients, which ensures constant memory and linear update complexity. Together, these developments establish a bridge between Bayesian nonparametrics and classical estimation theory. The application domains discussed illustrate how sequential GPs have matured from theoretical constructs into practical tools for modern SP. In regression and forecasting, they provide uncertainty-aware alternatives to deep recurrent models. In distributed inference, they enable probabilistic cooperation across agents through consensus updates in information form. In robust learning and anomaly detection, they produce calibrated likelihoods that expose distributional shifts or sensor faults. In sequential Bayesian optimization and adaptive sensing, they guide exploration and data acquisition under resource constraints. Across these tasks, the GP framework provides both predictive performance and interpretability, which are two qualities seldom achieved simultaneously in contemporary ML. Despite their progress, important research challenges remain. Current sequential GP algorithms still depend on stationary kernels or a limited number of basis functions, which restricts their expressiveness in high-dimensional and strongly nonstationary systems. More flexible kernel parameterizations, possibly driven by neural feature maps, promise to extend their applicability while maintaining the probabilistic interpretation that distinguishes GPs from purely data-driven models. Another open direction concerns scalability in multi-agent networks, where thousands of sensors or robots must coordinate their inference under communication constraints. Consensus-based GP

D. Waxman et al./Sequential Inference for Gaussian Processes

filtering and distributed variational inference may offer a path forward. Robustness also remains an active frontier, particularly under heavy-tailed or adversarial noise; generalized Bayesian inference and α-divergence formulations suggest principled solutions that retain analytic structure. Several emerging areas invite further exploration. In physical modeling, hybrid GP–physics approaches link latent stochastic dynamics with partial differential equations and could redefine data assimilation in geoscience, fluid mechanics, and biomedicine. In signal intelligence and communications, GPs may provide a foundation for probabilistic waveform synthesis and adaptive spectrum management, where uncertainty quantification directly influences control decisions. In neuroscience, GP state–space models already enable interpretable decoding of neural activity and could evolve into real-time brain–computer interfaces through scalable sequential inference. In robotics, Markovian GPs integrated with diffusion or neural-SDE representations may offer a unified framework for perception, prediction, and control under uncertainty. From a broader perspective, GPs exemplify the convergence of statistical learning and SP. They formalize estimation as inference in function space while maintaining the recursive and modular character that defines signal-processing methodology. As computing architectures continue to favor streaming and decentralized operation, the relevance of sequential GP formulations will only increase. Future progress will likely emerge through cross-fertilization with deep generative modeling, reinforcement learning, and information-theoretic control, where uncertainty must be represented, propagated, and acted upon in real time. In summary, GPs have evolved into a versatile paradigm that extends the classical theory of signals and systems to the probabilistic domain. Their sequential formulations inherit the rigor of Bayesian inference and the efficiency of recursive estimation, which provides a coherent path from theory to implementation. Continued integration of these models with advances in optimization, numerical linear algebra, and neural representation learning promises to expand their role as central tools for signal processing in the coming decade. Acknowledgment The authors would like to acknowledge the support of the National Science Foundation under Award 2212506.

D. Waxman et al./Sequential Inference for Gaussian Processes

References Adam, V., Chang, P., Khan, M. E. E. and Solin, A. (2021). Dual parameterization of sparse variational Gaussian processes. Advances in Neural Information Processing Systems 34 11474–11486. Al-Shedivat, M., Wilson, A. G., Saatchi, Y., Hu, Z. and Xing, E. P. (2017). Learning scalable deep kernels with recurrent structure. Journal of Machine Learning Research 18 1–37. Archambeau, C., Cornford, D., Opper, M. and Shawe-Taylor, J. (2007). Gaussian process approximations of stochastic differential equations. In Gaussian Processes in Practice 1–16. PMLR. Avron, H., Sindhwani, V., Yang, J. and Mahoney, M. W. (2016). Quasi-Monte Carlo feature maps for shift-invariant kernels. Journal of Machine Learning Research 17 1–38. Bar-Shalom, Y., Li, X. R. and Kirubarajan, T. (2004). Estimation with Applications to Tracking and Navigation: Theory, Algorithms, and Software. Wiley. Bock, C., Aubet, F.-X., Gasthaus, J., Kan, A., Chen, M. and Callot, L. (2022). Online time series anomaly detection with state space Gaussian processes. arXiv preprint arXiv:2201.06763. Bui, T. D., Nguyen, C. and Turner, R. E. (2017). Streaming sparse Gaussian process approximations. Advances in Neural Information Processing Systems 30. Bui, L. M., Huu, T. T., Dinh, D., Nguyen, T. M. and Hoang, T. N. (2025). Revisiting Kernel Attention with Correlated Gaussian Process Representation. In The 40th Conference on Uncertainty in Artificial Intelligence (UAI). Chang, L. and Li, K. (2017). Unified form for the robust Gaussian information filtering based on M-estimate. IEEE Signal Processing Letters 24 412–416. Chang, P. E., Verma, P., John, S., Solin, A. and Khan, M. E. (2023). Memory-based dual Gaussian processes for sequential learning. In International Conference on Machine Learning 4035–4054. PMLR. Cressie, N. and Wikle, C. K. (2011). Statistics for spatio-temporal data. John Wiley & Sons. Damianou, A. and Lawrence, N. D. (2013). Deep gaussian processes. In Artificial intelligence and statistics 207–215. PMLR. Dowling, M., Sokól, P. and Park, I. M. (2021). Hida-Matérn Kernel.

D. Waxman et al./Sequential Inference for Gaussian Processes

arXiv preprint arXiv:2107.07098. Dowling, M., Zhao, Y. and Park, I. M. (2023). Linear time GPs for inferring latent trajectories from neural spike trains. In Proceedings of the 40th International Conference on Machine Learning 8424–8448. El-Laham, Y., Yang, L., Djurić, P. M. and Bugallo, M. F. (2021). Particle filtering under general regime switching. In 2020 28th European Signal Processing Conference (EUSIPCO) 2378–2382. IEEE. Engel, Y., Mannor, S. and Meir, R. (2004). The kernel recursive leastsquares algorithm. IEEE Transactions on Signal Processing 52 2275–2285. Fortuin, V., Baranchuk, D., Rätsch, G. and Mandt, S. (2020). GPVAE: Deep probabilistic time series imputation. In International conference on artificial intelligence and statistics 1651–1661. PMLR. Garcı́a-Fernández, Á. F., Tronarp, F. and Särkkä, S. (2019). Gaussian process classication using posterior linearization. IEEE Signal Processing Letters 26 735–739. Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D. and Wilson, A. G. (2018). gPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. Advances in Neural Information Processing Systems 31. Gijsberts, A. and Metta, G. (2013). Real-time model learning using incremental sparse spectrum Gaussian process regression. Neural Networks 41 59–69. Hamelijnck, O., Wilkinson, W., Loppi, N., Solin, A. and Damoulas, T. (2021). Spatio-temporal variational Gaussian processes. Advances in Neural Information Processing Systems 34 23621–23633. Hensman, J., Fusi, N. and Lawrence, N. D. (2013). Gaussian processes for Big data. In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence 282–290. Hoeting, J. A., Madigan, D., Raftery, A. E. and Volinsky, C. T. (1999). Bayesian model averaging: a tutorial (with comments by M. Clyde, David Draper and EI George, and a rejoinder by the authors). Statistical Science 14 382–417. Kass, R. E., Tierney, L. and Kadane, J. B. (1991). Laplace’s method in Bayesian analysis. Contemporary Mathematics 115 89–99. Khan, M. and Lin, W. (2017). Conjugate-computation variational inference: Converting variational inference in non-conjugate models to inferences in conjugate models. In Artificial Intelligence and Statistics 878–887. PMLR.

D. Waxman et al./Sequential Inference for Gaussian Processes

Ko, J. and Fox, D. (2009). GP-BayesFilters: Bayesian filtering using Gaussian process prediction and observation models. Autonomous Robots 27 75–90. Laplante, W., Altamirano, M., Duncan, A. B., Knoblauch, J. and Briol, F.-X. (2025). Robust and Conjugate Spatio-Temporal Gaussian Processes. In Forty-second International Conference on Machine Learning. PMLR. Lázaro-Gredilla, M., Quinonero-Candela, J., Rasmussen, C. E. and Figueiras-Vidal, A. R. (2010). Sparse spectrum Gaussian process regression. The Journal of Machine Learning Research 11 1865–1881. Lin, Z., Sun, Y., Yin, F. and Thiéry, A. H. (2024). Ensemble Kalman filtering meets Gaussian process SSM for non-mean-field and online inference. IEEE Transactions on Signal Processing 72 4286–4301. Lindgren, F., Bolin, D. and Rue, H. (2022). The SPDE approach for Gaussian and non-Gaussian fields: 10 years and still running. Spatial Statistics 50 100599. Llorente, F. and Djurić, P. M. (2024). Dynamic random feature Gaussian processes for Bayesian optimization of time-varying functions. In ICASSP 2024-2024 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) 9756–9760. IEEE. Llorente, F., Waxman, D. and Djurić, P. M. (2025). Decentralized Online Ensembles of Gaussian Processes for Multi-Agent Systems. In 2025 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). Llorente, F., Waxman, D., Jantre, S., Urban, N. M. and Minkoff, S. E. (2026). Robust, Online, and Adaptive Decentralized Gaussian Processes. In ICASSP 2026-2026 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) 1–5. IEEE. Loper, J., Blei, D., Cunningham, J. P. and Paninski, L. (2021). A general linear-time inference method for Gaussian Processes on one dimension. Journal of Machine Learning Research 22 1–36. Lu, Q., Karanikolas, G. V. and Giannakis, G. B. (2022). Incremental ensemble Gaussian processes. IEEE Transactions on Pattern Analysis and Machine Intelligence 45 1876–1893. Lu, Q., Polyzos, K. D., Li, B. and Giannakis, G. B. (2023). Surrogate modeling for Bayesian optimization beyond a single Gaussian process. IEEE Transactions on Pattern Analysis and Machine Intelligence 45 11283–11296.

D. Waxman et al./Sequential Inference for Gaussian Processes

Maddox, W. J., Stanton, S. and Wilson, A. G. (2021). Conditioning sparse variational Gaussian processes for online decision-making. Advances in Neural Information Processing Systems 34 6365–6379. Opper, M. and Archambeau, C. (2009). The variational Gaussian approximation revisited. Neural computation 21 786–792. Pérez-Cruz, F., Van Vaerenbergh, S., Murillo-Fuentes, J. J., Lázaro-Gredilla, M. and Santamaria, I. (2013). Gaussian processes for nonlinear signal processing: An overview of recent advances. IEEE Signal Processing Magazine 30 40–50. Pförtner, M., Wenger, J., Cockayne, J. and Hennig, P. (2025). Computation-Aware Kalman Filtering and Smoothing. In International Conference on Artificial Intelligence and Statistics 2071–2079. PMLR. Pinder, T. and Dodd, D. (2022). GPJax: A Gaussian process framework in JAX. Journal of Open Source Software 7 4455. Quinonero-Candela, J. and Rasmussen, C. E. (2005). A unifying view of sparse approximate Gaussian process regression. Journal of machine learning research 6 1939–1959. Rahimi, A. and Recht, B. (2007). Random features for large-scale kernel machines. Advances in Neural Information Processing Systems 20. Rasmussen, C. E. and Williams, C. K. I. (2005). Gaussian Processes for Machine Learning. The MIT Press. Sarkka, S. and Hartikainen, J. (2012). Infinite-dimensional Kalman filtering approach to spatio-temporal Gaussian process regression. In Artificial Intelligence and Statistics 993–1001. PMLR. Särkkä, S. and Solin, A. (2019). Applied Stochastic Differential Equations 10. Cambridge University Press. Schürch, M., Azzimonti, D., Benavoli, A. and Zaffalon, M. (2020). Recursive estimation for sparse Gaussian process regression. Automatica 120 109127. Solin, A. and Särkkä, S. (2020). Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing 30 419–446. Titsias, M. (2009). Variational learning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics 567–574. PMLR. Titsias, M., Galashov, A., Rannen-Triki, A., Pascanu, R., Teh, Y. W. and Bornschein, J. (2024). Kalman Filter for Online Classification of Non-Stationary Data. In The Twelfth International Conference on Learning Representations.

D. Waxman et al./Sequential Inference for Gaussian Processes

Tronarp, F., Garcı́a-Fernández, Á. F. and Särkkä, S. (2018). Iterative Filtering and Smoothing in Nonlinear and Non-Gaussian Systems Using Conditional Moments. IEEE Signal Processing Letters 25 408–412. Van Vaerenbergh, S., Lázaro-Gredilla, M. and Santamarı́a, I. (2012). Kernel recursive least-squares tracker for time-varying regression. IEEE Transactions on Neural Networks and Learning Systems 23 1313– 1326. Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L. and Polosukhin, I. (2017). Attention Is All You Need. In Advances in Neural Information Processing Systems. Verma, P., Adam, V. and Solin, A. (2024). Variational Gaussian process diffusion processes. In International Conference on Artificial Intelligence and Statistics 1909–1917. PMLR. Waxman, D. and Djurić, P. M. (2024a). Dynamic Online Ensembles of Basis Expansions. Transactions on Machine Learning Research. Waxman, D. and Djurić, P. M. (2024b). A Gaussian Process-based Streaming Algorithm for Prediction of Time Series With Regimes and Outliers. In 2024 27th International Conference on Information Fusion (FUSION). Waxman, D., Llorente, F. and Djurić, P. (2026). Bayesian Ensembling: Insights from Online Optimization and Empirical Bayes. Transactions on Machine Learning Research. Wilkinson, W., Chang, P., Andersen, M. and Solin, A. (2020). State space expectation propagation: Efficient inference schemes for temporal Gaussian processes. In International Conference on Machine Learning 10270–10281. PMLR. Wilson, A. and Adams, R. (2013). Gaussian process kernels for pattern discovery and extrapolation. In International Conference on Machine Learning 1067–1075. PMLR. Wilson, A. G., Hu, Z., Salakhutdinov, R. and Xing, E. P. (2016). Deep kernel learning. In Artificial Intelligence and Statistics 370–378. PMLR. Xu, J., Lu, Q. and Giannakis, G. B. (2025). Online scalable Gaussian processes with conformal prediction for guaranteed coverage. In ICASSP 2025-2025 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) 1–5. IEEE. Yang, L., Zhang, Z., Song, Y., Hong, S., Xu, R., Zhao, Y., Zhang, W., Cui, B. and Yang, M.-H. (2023). Diffusion models: A com-

D. Waxman et al./Sequential Inference for Gaussian Processes

prehensive survey of methods and applications. ACM Computing Surveys 56 1–39. Yu, F. X. X., Suresh, A. T., Choromanski, K. M., HoltmannRice, D. N. and Kumar, S. (2016). Orthogonal random features. Advances in Neural Information Processing Systems 29. Zhao, Z., Emzir, M. and Särkkä, S. (2021). Deep state-space Gaussian processes. Statistics and Computing 31 75. Ziatdinov, M., Liu, Y., Kelley, K., Vasudevan, R. and Kalinin, S. V. (2022). Bayesian active learning for scanning probe microscopy: From Gaussian processes to hypothesis learning. ACS Nano 16 13492–13512.

Record · ID 146914 · SHA-256 cbefe27b98bded5a
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.