Bayesian model averaging (BMA) is the most common approach to linear weighting.

BMA weighs each model according to its a posteriori probability:

\(w_i = \mathrm{Pr}(M_i | \mathcal{D}) = \frac{p(M_i | \mathcal{D}) p(M_i)}{\sum_{j=1}^K p(M_j | \mathcal{D}) p(M_j)}\).

In the online setting, we obtain recursive weights via the prequential decomposition of the log-likelihood:

\(w_{t+1,k}^{\text{BMA}} = \frac{ w_{t, k} p_k(y_{t+1} | \mathbf{x}_{t+1}, \mathcal{D}_t)}{\sum_k w_{t, k} p_k(y_{t+1} | \mathbf{x}_{t+1}, \mathcal{D}_t)}\).

BMA is simple, and often effective. However, it does a poor job of model ensembling , and actually behaves like model selection.

We can show this explicitly in the online setting.

In fact, we can show this collapse behaves exactly like an online model selection algorithm!

One way to beat BMA is to optimize the weights themselves.

Bayesian stacking [Yao et al., BA '18] optimizes weights over a validation set:

\( w^* = \arg \max_{w \in \mathbb{S}^K} \sum_{n=1}^N \log \sum_{k=1}^K w_k p_k(y_n | \mathbf{x}_n, \mathcal{D}) \)

In the offline setting, this can be solved efficiently using leave-one-out cross-validation...

... but what if we want to solve for the weights online?

For the online scenario, we introduce a prequential version of the online stacking objective,

\({w}^* \coloneqq \arg \max_{w\in \mathbb{S}^K} \sum_{t=0}^{T} \log \sum_{k=1}^K w_k p_k(y_{t} | \mathbf{x}_{t}, \mathcal{D}_{1:t-1}).\)

This has been approached in time series settings before, but never for online learning.

We call the resulting problem and its analysis online Bayesian stacking (OBS).

The resulting optimization problem turns out to be exactly equivalent to a well-studied problem in finance and online convex optimization!

OPS algorithms come with robust regret guarantees, that guarantee strong performance in adversarial conditions. Let's look at one such algorithm.

One of the exciting things about this equivalence is that many algorithms exist to provide different bounds! Some variants include

• The Online Newton Step (ONS) algorithm, which provides logarithmic regret with bounded gradients.
• The dynamic ONS (D-ONS) algorithm, which adapts ONS to non-stationary environments and can beat the BCRP.
• Algorithms with data-dependent regret bounds.
• Algorithms with robust switching regret guarantees.

We begin with a famous toy example to illustrate the effectiveness of OBS.

We consider data from a linear model, where each model is missing input features.

Similar results hold in real, non-trivial Bayesian machine learning settings.

For brevity, here we consider two settings: online variational inference on MNIST, and an online forecasting task using a variety of GARCH models and sequential Monte Carlo.

Essentially all of these online Bayesian models amount to some sort of filtering algorithm.

Filtering algorithms use sequential Bayesian inference to update the posterior over some variable \(\theta\) as new data \(y_t\) is observed.

As OBS provides adversarial regret guarantees, it is particularly well-suited to classical multiple model filtering algorithms.

Indeed, O-BMA is the standard in bank of filter algorithms, e.g., in the GARCH example.

Ongoing Work (Email Me If Interested): OBS can be used to improve classical multiple model filtering, including static multiple model and interacting multiple model filters.

OBS was agnostic to how we learned the models in our ensemble -- but how can we train and ensemble simultaneously?

There are two primary challenges here:

• How do we obtain closed-form updates for online learning?
• How do we adapt to non-stationary environments?

These are both difficult to accomplish with GPs in their conventional formulation. We'll address both by reformulating using a linear-Gaussian approximation of GPs.

The point here is that, under a linear-Gaussian approximation to a GP, we can solve for the posterior exactly, via the Kalman filter.

The resulting models are generally known as basis expansions, which use a feature mapping \(\phi(x)\) to obtain non-linear features.

To obtain the online ensemble of basis expansions (OEBE), we can additionally combine with OBS or O-BMA.

One popular approximation for stationary GPs are Random Fourier Features (RFF) [Rahimi and Recht, NeurIPS '07].

These features are derived by a Monte Carlo approximation of the Wiener-Khinchin integral

\(\kappa(t, t') = \int_{\mathbb{R}^d} \exp(2\pi i s^\top (t - t')) \, dS(s),\)

where \(S(s)\) is a spectral density and \(i\) is the imaginary unit.

The resulting basis expansion is given by sampling \(F / 2\) spectral frequencies, \(s_1, \dots, s_{F/2}\), and computing the features

\(\phi(x) = \left[ \cos(2\pi s_1^\top x), \sin(2\pi s_1^\top x) \cdots \cos(2\pi s_{F/2}^\top x), \sin(2\pi s_{F/2}^\top x) \right].\)

Takeaway: The RFF provides an approximation that can be sequentially computed, and has computational cost independent of the data size.

RFF, while simple, suffers harshly from the curse of dimensionality.

Hilbert Space Gaussian Processes (HSGP) [Solin & Sarkka, S&C '20'] address this by instead deriving a quadrature rule over the Wiener-Khinchin integral.

For scalar GPs on \([-L, L]\), the corresponding basis function is

\( \phi_k(t) = S\left(\frac{k\pi}{2L}\right) \frac{\sin\left((k\pi(t+L) / 2L)\right)}{\sqrt{L}}, \quad k = 1, \dots, F. \)

where \(S(s)\) is a spectral density.

Corresponding approximations in higher dimensions are possible but expensive. We instead pursue additive HSGP approximations, which writes \(f(x) = \sum_{d=1}^D f_d(x_d)\).

Takeaway: There are many different basis expansions; these are also a promising place for ensembles.

From the filtering perspective, a random walk on parameters is an obvious choice.

In particular, with a Gaussian random walk, we maintain exact Kalman filter updates.

We adapt a hyperparameter \(\sigma_\text{rw}^2\), and the state-space model

\( \theta_{t+1}^{(m)} = \theta_t^{(m)} + \varepsilon_{\text{rw}}^{(m)}, \)

and observation model

\( y^{(m)} = \phi^{(m)}(x)^\top \theta_t^{(m)} + \varepsilon_{t}^{(m)}. \)

We call the resulting model a distributed online ensemble of basis expansions (DOEBE) [DW & Djuric, TMLR '24].

The explicit equivalence with Gaussian filtering opens many opportunities for extensions.

For example, we can use techniques from robust filtering theory for robust GPs.

Another example includes stronger non-linear filters through posterior linearization.

Through similar means, decentralized solutions are possible [Llorente*, DW*, & Djuric, ICASSP '25; Llorente, DW, et al., ICASSP '26].

These fall from using the information form of the Kalman filter, which makes all quantities additive.

Strong consensus algorithms are thus available, and provide decentralized inference that is asymptotically exact.

Our solutions to online inference have worked by translating to the language of dynamical systems.

This is quite powerful, for a few reasons:

• The resulting problem is naturally sequential, which is already close to online inference.
• The problem is well-studied, with many bespoke algorithms and approximations.

But how can non-experts exploit these powerful solutions easily?

Dynestyx gives a way to rapidly prototype for applications: let's find those applications!

With applications in-hand, understand how estimation methods actually fail, in practice, and when.

With our programmatic tools, learn how to do proper hybrid modeling, i.e., when we mix mechanistic understanding with data-driven components.