Low-rank extended Kalman filtering for online learning of neural networks from streaming data

Suppose we observe a stream of labeled observations, ${\mathcal{D}}_{t}=\{({\bm{x}}_{t}^{n},{\bm{y}}_{t}^{n})\sim p_{t}({\bm{x}},{\bm{y}}):n=1{:}N_{t}\}$, where ${\bm{x}}_{t}^{n}\in\mathcal{X}=^{D}$, ${\bm{y}}_{t}^{n}\in\mathcal{Y}=^{C}$, and $N_{t}$ is the number of examples at step $t$. (In this paper, we assume $N_{t}=1$, since we are interested in rapid learning from individual data samples.) Our goal is to fit a prediction model ${\bm{y}}_{t}=h({\bm{x}}_{t},{\bm{\theta}})$ in an online fashion, where ${\bm{\theta}}\in^{P}$ are the parameters of the model. (We focus on the case where $h$ is a deep neural network (DNN), although in principle our methods can also be applied to other (differentiable) parametric models.) In particular, we want to recursively estimate the posterior over the parameters

without having to store all the past data. Here $p({\bm{\theta}}|{\mathcal{D}}_{1:t-1})$ is the posterior belief state from the previous step, and $p({\bm{y}}_{t}|{\bm{x}}_{t},{\bm{\theta}})$ is the likelihood function given by

For regression, we assume $h({\bm{x}}_{t},{\bm{\theta}})\in^{C}$ returns the mean of the output, and $\mathbf{R}_{t}=R\mathbf{I}_{C}$ is the observation covariance, which we view as a hyper-parameter. For classification, $h({\bm{x}}_{t},{\bm{\theta}})$ returns a $C$-dimensional vector of class probabilities, which is the mean parameter of the categorical distribution.

In many problem settings (e.g., recommender systems (Huang et al., 2015), robotics (Wołczyk et al., 2021; Lesort et al., 2020), and sensor networks (Ditzler et al., 2015)), the data distribution $p_{t}({\bm{x}},{\bm{y}})$ may change over time (Gomes et al., 2019). Hence we allow the model parameters ${\bm{\theta}}_{t}$ to change over time, according to a simple Gaussian dynamics model:¹¹1 We do not assume access to any information about if and when the distribution shifts (sometimes called a “task boundary”), since such information is not usually available. Furthermore, the shifts may be gradual, which makes the concept of task boundary ill-defined.

where we usually take $\mathbf{Q}_{t}=q\mathbf{I}$ and $\gamma_{t}=\gamma$, where $q\geq 0$ and $0\leq\gamma\leq 1$. Using $q>0$ injects some noise at each time step, and ensures that the model does not lose “plasticity”, so it can continue to adapt to changes (cf. Kurle et al., 2020; Ash & Adams, 2020; Dohare et al., 2021), and using $\gamma<1$ ensures the variance of the unconditional stochastic process does not blow up. If we set $q=0$ and $\gamma=1$, this corresponds to a deterministic model in which the parameters do not change, i.e.,

This is a useful special case for when we want to estimate the parameters from a stream of data coming from a static distribution. (In practice we find this approach can also work well for the non-stationary setting.)

Recursively computing eq. 1 corresponds to Bayesian inference (filtering) in a state space model, where the dynamics model in eq. 3 is linear Gaussian, but the observation model in eq. 2 is non-linear and possibly non-Gaussian. Many approximate algorithms have been proposed for this task (see e.g. Sarkka, 2013; Murphy, 2023), but in this paper, we focus on Gaussian approximations to the posterior, $q({\bm{\theta}}_{t}|{\mathcal{D}}_{1:t})=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t},\bm{\Sigma}_{t})$, since they strike a good balance between efficiency and expressivity. In particular, we build on the extended Kalman filter (EKF), which linearizes the observation model at each step, and then computes a closed form Gaussian update. The EKF has been used for online training of neural networks in many papers (see e.g., Singhal & Wu, 1989; Watanabe & Tzafestas, 1990; Puskorius & Feldkamp, 1991; Iiguni et al., 1992; Ruck et al., 1992; Haykin, 2001). It can be thought of as an approximate Bayesian inference method, or as a natural gradient method for MAP parameter estimation (Ollivier, 2018), which leverages the posterior covariance as a preconditioning matrix for fast Newton-like updates (Alessandri et al., 2007). The EKF was extended to exponential family likelihoods in (Ollivier, 2018; Tronarp et al., 2018), which is necessary when fitting classification models.

The main drawback of the EKF is that it takes $O(P^{3})$ time per step, where $P=|{\bm{\theta}}_{t}|$ is the number of parameters in the hidden state vector, because we need to invert the posterior covariance matrix. It is possible to derive diagonal approximations to the posterior covariance or precision, by either minimizing $D_{\mathbb{KL}}\left({p({\bm{\theta}}_{t}|{\mathcal{D}}_{1:t})}\mathrel{\|}{q({\bm{\theta}}_{t})}\right)$ or $D_{\mathbb{KL}}\left({q({\bm{\theta}}_{t})}\mathrel{\|}{p({\bm{\theta}}_{t}|{\mathcal{D}}_{1:t})}\right)$, as discussed in (Puskorius & Feldkamp, 1991; Chang et al., 2022; Jones et al., 2023). These methods take $O(P)$ time per step, but can be much less statistically efficient than full-covariance methods, since they ignore joint uncertainty between the parameters. This makes the method slower to learn, and slower to adapt to changes in the data distribution, as we show in section 4.

In this paper, we propose an efficient and deterministic method to recursively minimize $D_{\mathbb{KL}}\left({\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t},\bm{\Sigma}_{t})}\mathrel{\|}{p({\bm{\theta}}_{t}|{\mathcal{D}}_{1:t})}\right)$, where we assume that the precision matrix is diagonal plus low-rank, $\bm{\Sigma}_{t}^{-1}=\bm{\Upsilon}_{t}+\mathbf{W}_{t}\mathbf{W}_{t}^{\intercal}$, where $\bm{\Upsilon}_{t}$ is diagonal and $\mathbf{W}_{t}\in^{P\times L}$ for some memory limit $L$. The key insight is that, if we linearize the observation model at each step, as in the EKF, we can use the resulting gradient vector or Jacobian as “pseudo-observation(s)” that we append to $\mathbf{W}_{t-1}$, and then we can perform an efficient online SVD approximation to obtain $\mathbf{W}_{t}$. We therefore call our method LO-FI, which is short for low-rank extended Kalman filter. Our code is available at https://github.com/probml/rebayes.

We use the posterior approximation $p({\bm{\theta}}_{t}|{\mathcal{D}}_{1:t})$ in two ways. First, under Bayesian updating the covariance matrix $\bm{\Sigma}_{t}$ acts as a preconditioning matrix to yield a deterministic second-order Newton-like update for the posterior mean (MAP estimate). This update does not have any step-size hyperparameters, in contrast to SGD. Second, the posterior uncertainty in the parameters can be propagated into the uncertainty of the predictive distribution for observations, which is crucial for online decision-making tasks, such as active learning (Holzmüller et al., 2022), Bayesian optimization (Garnett, 2023), contextual bandits (Duran-Martin et al., 2022), and reinforcement learning (Khetarpal et al., 2022; Wang et al., 2021).

In summary, our main contribution is a novel algorithm for efficiently (and deterministically) recursively updating a diagonal plus low-rank (DLR) approximation to the precision matrix of a Gaussian posterior for a special kind of state space model, namely an SSM with an arbitrary non-linear (and possibly non-Gaussian) observation model, but with a simple linear Gaussian dynamics. This model family is ideally suited to online parameter learning for DNNs in potentially non-stationary environments (but the restricted form of the dynamics model excludes some other applications of SSMs). We show experimentally that our approach works better (in terms of accuracy for a given compute budget) than a variety of baseline algorithms — including online gradient descent, online Laplace, diagonal approximations to the EKF, and a stochastic DLR VI method called L-RVGA — on a variety of stationary and non-stationary classification and regression problems, as well as a simple contextual bandit problem.

Since exact Bayesian inference is intractable in our model family, it is natural to compute an approximate posterior at step $t$ using recursive variational inference (VI), in which the prior for step $t$ is the approximate posterior from step $t-1$ (Opper, 1998; Broderick et al., 2013). That is, at each step we minimize the ELBO (evidence lower bound), which is equal (up to a constant) to the reverse KL, given by

where $Z_{t}$ is a normalization constant and $q_{t}=\mathcal{N}({\bm{\theta}}_{t}|\bm{\mu}_{t},\bm{\Sigma}_{t})$ is the variational posterior which results from minimizing this expression. The main challenge is how to efficiently optimize this objective.

One common approach is to assume the variational family consists of a diagonal Gaussian. By linearizing the likelihood, we can solve the VI objective in closed form, as shown in (Chang et al., 2022); this is called the “variational diagonal EKF” (VD-EKF). They also propose a diagonal approximation which minimizes the forwards KL, $D_{\mathbb{KL}}\left({p({\bm{\theta}}_{t}|{\mathcal{D}}_{1:t})}\mathrel{\|}{q({\bm{\theta}}_{t})}\right)$, and show that this is equivalent to the “fully decoupled EKF” (FD-EKF) method of (Puskorius & Feldkamp, 1991). Both of these methods are fully deterministic, which avoids the high variance that often plagues stochastic VI methods (Wu et al., 2019; Haußmann et al., 2020).

It is also possible to derive diagonal approximations without linearizing the observation model. In (Kurle et al., 2020; Zeno et al., 2018) they propose a diagonal approximation to minimize the reverse KL, $D_{\mathbb{KL}}\left({q({\bm{\theta}}_{t})}\mathrel{\|}{p({\bm{\theta}}_{t}|{\mathcal{D}}_{1:t})}\right)$; this requires a Monte Carlo approximation to the ELBO. In (Ghosh et al., 2016; Wagner et al., 2022), they propose a diagonal approximation to minimize the forwards KL, $D_{\mathbb{KL}}\left({p({\bm{\theta}}_{t}|{\mathcal{D}}_{1:t})}\mathrel{\|}{q({\bm{\theta}}_{t})}\right)$; this requires approximating the first and second moments of the hidden units at every layer of the model using numerical integration.

(Farquhar et al., 2020) claims that, if one makes the model deep enough, one can get good performance using a diagonal approximation; however, this has not been our experience. This motivates the need to go beyond a diagonal approximation.

One approach is to combine diagonal Gaussian approximations with memory buffers, such as the variational continual learning methd of (Nguyen et al., 2018) and other works (see e.g., (Kurle et al., 2020; Khan & Swaroop, 2021)). However, we seek to find a richer approximation to the posterior that does not rely on memory buffers, which can be problematic in the non-stationary setting.

(Zeno et al., 2021) proposes the FOO-VB method, which uses a Kronecker block structured approximation to the posterior covariance. However, this method requires 2 SVD decompositions of the Kronecker factors for every layer of the model, in addition to a large number of Monte Carlo samples, at each time step. In (Ong et al., 2018) they compute a diagonal plus low-rank (DLR) approximation to the posterior covariance matrix using stochastic gradient applied to the ELBO. In (Tomczak et al., 2020) they develop a version of the local reparameterization trick for the DLR posterior covariance, to reduce the variance of the stochastic gradient estimate.

In this paper we use a diagonal plus low-rank (DLR) approximation to the posterior precision. The same form of approximation has been used in several prior papers. In (Mishkin et al., 2018) they propose a technique called “SLANG” (stochastic low-rank approximate natural-gradient), which uses a stochastic estimate of the natural gradient of the ELBO to update the posterior precision, combined with a randomized eigenvalue solver to compute a DLR approximation. Their NGD approximation enables the variational updates to be calculated solely from the loss gradients, whereas our approach requires the network Jacobian. On the other hand, our EKF approach allows the posterior precision and the DLR approximation to be efficiently computed in closed form.

In (Lambert et al., 2021a), they propose a technique called “L-RVGA” (low-rank recursive variational Gaussian approximation), which uses stochastic EM to optimize the ELBO using a DLR approximation to the posterior precision. Their method is a one-pass online method, like ours, and also avoids the need to tune the learning rate. However, it is much slower, since it involves generating multiple samples from the posterior and multiple iterations of the EM algorithm (see fig. 7 for an experimental comparison of running time).

The GGT method of (Agarwal et al., 2019) also computes a DLR approximation to the posterior precision, which they use as a preconditioner for computing the MAP estimate. However, they bound the rank by simply using the most recent $L$ observations, whereas LO-FI uses SVD to combine the past data in a more efficient way.

The ORFit method of (Min et al., 2022) is also an online low-rank MAP estimation method. They use orthogonal projection to efficiently compute a low rank representation of the precision at each step. However, it is restricted to regression problems with 1d, noiseless outputs (i.e., they assume the likelihood has the degenerate form $p(y_{t}|{\bm{x}}_{t},{\bm{\theta}}_{t})=\mathcal{N}(h({\bm{x}}_{t},{\bm{\theta}}_{t}),0)$.)

The online Laplace method of (Ritter et al., 2018; Daxberger et al., 2021) also computes a Gaussian approximation to the posterior, but makes different approximations. In particular, for “task” $t$, it computes the MAP estimate ${\bm{\theta}}_{t}=\operatornamewithlimits{argmax}_{{\bm{\theta}}}\log p({\mathcal{D}}_{t}|{\bm{\theta}})+\log\mathcal{N}({\bm{\theta}}|\bm{\mu}_{t-1},\bm{\Sigma}_{t-1})$, where $\bm{\Sigma}_{t-1}=\bm{\Lambda}_{t-1}^{-1}$ is the approximate posterior covariance from the previous task. (This optimization problem is solved using SGD applied to a replay buffer.) This precision matrix is usually approximated as a block diagonal matrix, with one block per layer, and the terms within each block may be additionally approximated by a Kronecker product form, as in KFAC (Martens & Grosse, 2015). By contrast, LO-FI computes a posterior, not just a point estimate, and approximates the precision as diagonal plus low rank. In the appendix, we show experimentally that LO-FI outperforms online Laplace in terms of NLPD on various classification and regression tasks.

It is possible to go beyond Gaussian approximations by using particle filtering (see e.g., (Yang et al., 2023)). However, we focus on faster deterministic inference methods, since speed is important for many real time online decision making tasks (Ghunaim et al., 2023).

There are many papers on continual learning, which is related to online learning. However the CL literature usually assumes the task boundaries, corresponding to times when the distribution shifts, are given to the learner (see e.g., (Delange et al., 2021; De Lange & Tuytelaars, 2021; Wang et al., 2022; Mai et al., 2022; Mundt et al., 2023; Wang et al., 2023).) By contrast, we are interested in the continual learning setting where the distribution may change at unknown times, in a continuous or discontinuous manner (c.f., (Gama et al., 2013)); this is sometimes called the “task agnostic” or “streaming” setting. Furthermore, our goal is accurate forecasting of the future (which can be approximated by our estimate of the “current” distribution), so we are less concerned with performance on “past” distributions that the agent may not encounter again; thus “catastrophic forgetting” (see e.g., (Parisi et al., 2019)) is not a focus of this work (c.f., (Dohare et al., 2021)).

In LO-FI, we approximate the belief state by a Gaussian, $p({\bm{\theta}}_{t}|{\mathcal{D}}_{1:t})=\mathcal{N}(\bm{\mu}_{t},\bm{\Sigma}_{t})$, where the posterior precision is diagonal plus low rank, i.e., it has the form $\bm{\Sigma}_{t}^{-1}=\bm{\Upsilon}_{t}+\mathbf{W}_{t}\mathbf{W}_{t}^{\mkern-1.5mu\mathsf{T}}$, where $\bm{\Upsilon}_{t}$ is diagonal and $\mathbf{W}_{t}$ is a $P\times L$ matrix. We denote this class of models by $\text{DLR}(L)$, where $L$ is the rank. Below we show how to efficiently update this belief state in a recursive (online) fashion. This has two main steps — predict (see algorithm 2) and update (see algorithm 3) — which are called repeatedly, as shown in algorithm 1. The predict step takes $O(PL^{2}+L^{3})$ time, and the update step takes $O(P(L+C)^{2})$ time, where $C$ is the number of outputs.

In this section, we report experimental results on various classification and regression datasets. using the following approximate inference techniques: LO-FI (this paper); FDEKF (fully decoupled diagonal EKF) (Puskorius & Feldkamp, 2003); VDEKF (variational diagonal EKF) (Chang et al., 2022); SGD-RB (stochastic gradient descent with FIFO replay buffer), with memory buffer of size $B$, using either sgd or adam as the optimizer; online gradient descent (OGD), which corresponds to SGD-RB with $B=1$; the LRVGA method of (Lambert et al., 2021a) (for the NLPD results in section D.1); and the online Laplace approximation of (Ritter et al., 2018) (for the NLPD results in section D.3 and section D.6). For additional results, see appendix D. For the source code to reproduce these results, see https://github.com/probml/rebayes.

We have presented an efficient new method of fitting neural networks online to streaming datasets, using a diagonal plus low-rank Gaussian approximation. In the future, we are interested in developing online methods for estimating the hyper-parameters, perhaps by extending the variational Bayes approach of (Huang et al., 2020; de Vilmarest & Wintenberger, 2021), or the gradient based method of (Greenberg et al., 2021). We would also like to further explore the predictive uncertainty created by our posterior approximation, to see if it can be used for sequential decision making tasks, such as Bayesian optimization or active learning. This may require the use of (online) deep Bayesian ensembles, to capture functional as well as parametric uncertainty.