Title: CONTINUAL LEARNING WITH BAYESIAN NEURAL NETWORKS FOR NON-STATIONARY DATA

Abstract: This work addresses continual learning for non-stationary data, using Bayesian neural networks and memory-based online variational Bayes. We represent the posterior approximation of the network weights by a diagonal Gaussian distribution and a complementary memory of raw data. This raw data corresponds to likelihood terms that cannot be well approximated by the Gaussian. We introduce a novel method for sequentially updating both components of the posterior approximation. Furthermore, we propose Bayesian forgetting and a Gaussian diffusion process for adapting to non-stationary data. The experimental results show that our update method improves on existing approaches for streaming data. Additionally, the adaptation methods lead to better predictive performance for non-stationary data.

1. INTRODUCTION:

Continual learning (CL), also referred to as lifelong learning, is typically described informally by the following set of desiderata for computational systems: the system should (i) learn incrementally from a data stream, (ii) exhibit information transfer forward and backward in time, (iii) avoid catastrophic forgetting of previous data, and (iv) adapt to changes in the data distribution (Ring, 1997; Silver et al., 2013; Chen & Liu, 2016; Ruvolo & Eaton, 2013; Parisi et al., 2018) . The necessity to adapt to non-stationary data is often not reconcilable with the goal of preventing forgetting. This problem is also known as the stability-plasticity dilemma (Grossberg, 1987) .

The majority of current CL research is conducted in the context of online multi-task learning (Nguyen et al., 2018; Kirkpatrick et al., 2017; Schwarz et al., 2018; Rusu et al., 2016; Fernando et al., 2017) , where the main objective is to prevent catastrophic forgetting of previously learned tasks. This focus is reasonable since changes in the statistics of the data distribution are usually an artefact of learning different tasks sequentially. However, changes in the statistics of the data can also be real properties of the data-generating process. Examples include models of energy demand, climate analysis, financial market, or user-behavior analytics (Ditzler et al., 2015) . In such applications, the statistics of the current data distribution are of particular interest. Old data may be outdated and can even deteriorate learning if the drift in the data distribution is neglected. Consequently, CL systems for non-stationary data require adaptation methods, which deliberately forget outdated information.

In this work, we develop an approximate Bayesian approach for training Bayesian neural networks (BNN) (Hinton & van Camp, 1993; Graves, 2011; Blundell et al., 2015) incrementally with nonstationary streaming data. Similar to variational continual learning (VCL) (Nguyen et al., 2018) and the Virtual Vector Machine (VVM) (Minka et al., 2009) , we approximate the posterior using a Gaussian distribution and a complementary memory of previous data. Both components are updated sequentially, while adapting to changes in the data distribution. Our main contributions are as follows:

• We propose an online approximation consisting of a diagonal Gaussian distribution and a running memory, and we provide a novel sequential update method for both components. • We extend the online approximation by two alternative adaptation methods, thereby generalising online variational Bayes with Bayesian neural networks to non-stationary data.

We compare our sequential update method to VCL in the online-inference setting on several popular datasets, demonstrating that our approach is favorable. Furthermore, we validate our adaptation methods on several datasets with concept drift (Widmer & Kubat, 1996) , showing performance improvements compared to online variational Bayes without adaptation.

2. BACKGROUND: ONLINE INFERENCE:

Consider a stream of datasets {D t k } K k=1 , where t k are the time points at which datasets D t k are observed. For the moment, we assume that these datasets and the samples within are generated independently and identically distributed (i.i.d.) . Methods for non-i.i.d. data are considered in Sec. 4.

In the Bayesian approach to online learning, we want to infer the posterior distribution p(w|D t1:t k ) of our model parameters, with the restriction that the data is processed sequentially. 1 Using Bayes rule, a recursive posterior inference equation emerges naturally:

p(w|D t1:t k ) ∝ p(w|D t1:t k-1 ) p(D t k |w, D t1:t k-1 ) = p(w|D t1:t k-1 ) p(D t k |w),

(1) where the last step follows from the i.i.d. assumption of the data.

In this paper, we consider Gaussian and multinomial likelihoods, parametrised by a neural network with weights w and prior p(w|∅) = p 0 (w) = N (w; µ 0 , σ 0 ). Furthermore, we consider supervised learning, where

D t k = {d (n) t k } n = {(x (n) t k , y (n) t k )} n and p(d (n) t k |w) = p(y (n) t k | NN x (n) t k ; w) .

2.1. ONLINE VARIATIONAL BAYES:

Since exact Bayesian inference is intractable for non-trivial models, various approximations have been developed. Prominent examples include sequential Monte Carlo (Liu & Chen, 1998) , assumed density filtering (Maybeck, 1982) , and online variational Bayes (Opper, 1998; Ghahramani, 2000; Sato, 2001; Broderick et al., 2013) . Online variational Bayes (VB) approximates the posterior of Eq. ( 1) by a parametrised distribution q θt k (w) ≈ p(w|D t1:t k ) through a sequence of projections:

EQUATION

where Z t k is the normalisation constant. The above minimisation is equivalent to maximising the evidence lower bound (ELBO) L t k (θ; D t k ) = E q θ (w) log p(D t k |w) -KL q θ (w) || q θt k-1 (w) . In this work, we consider diagonal Gaussian posterior approximations q θt k (w) for the neural network weights, similar to Nguyen et al. (2018) .

2.2. ONLINE VARIATIONAL BAYES WITH MEMORY:

Online approximate Bayesian inference methods inevitably suffer from an information loss due to the posterior approximation at each time-step. An alternative approach to online learning is to store and update a representative dataset/generative model-and to use it as a memory-in order to improve inference (Robins, 1995; Lopez-Paz & Ranzato, 2017; Shin et al., 2017; Kamra et al., 2017) . Memory-based online learning has also been combined with online Bayesian inference methods (Minka et al., 2009; Nguyen et al., 2018) . A common property of these approaches is to represent the (current) posterior approximation by a product of two factors p(w|D t1:

EQUATION

and update them sequentially as new data

D t k is observed. The factor p(M t k |w) = M m p(m (m)

t k |w) is the likelihood of a set of M = |M| data points, which we refer to as running memory; and q θt k (w) is a Gaussian distribution, which summarises the rest of the data D1:

t k = D 1:t k \ M t k .

In case of VCL, the factors in Eq (3) are updated in two steps, which we refer to as (i) memory update and (ii) Gaussian update: (i) a new memory M t k ⊂ D t k ∪ M t k-1 is selected using heuristics such as random selection or the k-center method (a greedy algorithm that selects K data points based on geometric properties of

D t k ∪ M t k-1 .); (ii) the Gaussian distribution is updated with the remaining data Dt k = D t k ∪ M t k-1 \ M t k (using Eq. (2)) to obtain q θt k (w) ≈ q θt k-1 (w) p( Dt k |w).

Note that we cannot sample directly from the posterior approximation in Eq. ( 3) and thus we cannot easily evaluate quantities such as the posterior predictive distribution. VCL therefore performs a second projection

EQUATION

This distribution should not be confused with the recursively updated variational distribution (Eq. ( 2)).

3. IMPROVING MEMORY-BASED ONLINE VARIATIONAL BAYES:

In this section, we focus on two problems of existing approaches using online VB with a running memory: (i) the memory update does not take into account the approximation error or approximation capabilities of the variational distribution; (ii) the Gaussian update-performed by optimising the ELBO (Eq. ( 2)) only with data Dt k -can fail for streaming data. This is because VB yields poor posterior approximations if the dataset is too small or the neural network architecture has too much capacity (cf. Ghosh et al. (2018) , Fig. 1 ). In Secs. 3.2 and 3.3, we propose improvements to these two update methods. The mathematical background for our approach is provided in Sec. 3.1.

3.1. PROPERTIES OF THE GAUSSIAN VARIATIONAL APPROXIMATION:

There are two important properties of the Gaussian variational approximation that we will exploit later: (i) Gaussian approximate posterior distributions factorise into a product of Gaussian terms corresponding to the prior and each likelihood term; (ii) the ELBO can be written as the sum of the approximation's normalisation constant and a sum of residuals corresponding to these factors.

Let p 0 (w) = N (w; µ 0 , Σ 0 ) be a Gaussian prior and p(D|w) = n p(d (n) |w) be the likelihood of the observed data D. Furthermore, let q θ (w) = N (w; µ, Σ) denote the corresponding Gaussian variational approximation with θ = {µ, Σ}. Assume that µ and Σ are the optimal parameters corresponding to a (local) maximum of the ELBO L(µ, Σ; D). The optimality conditions ∂ µ L(µ, Σ; D) = 0 and ∂ Σ L(µ, Σ; D) = 0 can be rewritten as follows (Knowles & Minka, 2011; Opper & Archambeau, 2008; Cseke et al., 2013) (cf. App. C):

EQUATION

EQUATION

Since the sum of natural parameters corresponds to a product in distribution space, the above equations show that-at a local optimum-the approximation q θ (w) factorises in the same way as the posterior p(w|D). It can be written in the form q θ (w) = Z q -1 p 0 (w) n r (n) (w), where the factors r (n) (w) are Gaussian functions with natural parameters given by Eqs. (5a) and (5b), and where Z q = p 0 (w) n r (n) (w) dw is the normalisation constant. These Gaussian functions r (n) (w) each correspond to the contribution of the likelihood p(d (n) |w) to the posterior approximation q θ (w).

The resulting factorisation implies that the ELBO L(µ, Σ; D) can be written in the form (Opper & Winther, 2005) 

EQUATION

If the terms p(d (n) |w) were (diagonal) Gaussian in w, they would each cancel with the corresponding (diagonal) Gaussian term, leaving only log Z q . Intuitively, the residual terms in Eq. ( 6) can be used to quantify the quality of the Gaussian approximation.

3.2. MEMORY UPDATE:

The authors of VCL propose to use a memory to compensate the information loss resulting from the Gaussian approximation of the posterior distribution. However, their memory update is independent of the approximation error that is due to the chosen distributional family (diagonal Gaussian). An alternative memory update, which specifically targets the above mentioned information loss, has been introduced previously for VVM. Although the latter method was developed for expectation propagation in a (linear) logistic regression model-and is thus not directly applicable to online VBwe show that some of its properties can be transferred to the variational inference setting. The central idea is to replace the likelihood terms that can be well approximated by a Gaussian distribution by their Gaussian proxies p(d t k |w) ≈ r t k (w; d t k ) resulting in q θt k (w); and retain the data corresponding to the rest of the likelihood terms in the memory. To score a candidate memory, Minka et al. (2009) proposed to maximise the KL divergence between the model given in the form of Eq. ( 3) and a Gaussian posterior approximation, that is, maximise KL Z-1 t k q θt k (w) p(M|w) || qθt k (w) . However, this score function is intractable, because the expectation in the KL includes the likelihood p(M|w). In the following, we develop a tractable score function applicable to VB. Intuitively, we can use Eq. ( 6) to test how much L(µ, Σ; D) changes if we replace the exact likelihood terms (of all data which is not contained in the candidate memory) by their Gaussian approximations.

To achieve this, we need to find Gaussian approximations for every data point in the candidate memory. We first approximate the posterior distribution using both D t k and M t k-1 :

EQUATION

Next, we use Eqs. ( 5a) and (5b) to calculate the natural parameters of all Gaussian terms. In practice, we estimate the natural parameters using (unbiased) Monte-Carlo estimators for the expectations. We have now available the likelihood terms and their Gaussian approximations. This allows us to write

L(θ t k ; D t k ∪ M t k-1

) in the form of Eq. ( 6): 

L(θ t k ; D t k ∪ M t k-1 ) = log Z qt k + dt k ∈ Dt k ∪ Mt k-1 E qθ t k (w) log p(d t k |w) -log

EQUATION

and the corresponding memory update

M t k = argmax M S t k (M; D t k ∪ M t k-1

). Note that since all residual terms are computed independently, the update results in selecting the top M terms.

3.3. GAUSSIAN UPDATE:

The Gaussian update follows from the memory update presented in the previous section: once the memory M t k has been selected, we update the Gaussian distribution with the approximations corresponding to the rest of the data D t k ∪ M t k-1 \ M t k . We can update q θt k (w) in two equivalent ways:

q θ t k (w) = q θ t k-1 (w) d t k ∈ M t k rt k (w; dt k ), (9a) q θ t k (w) = qθ t k (w) / d t k ∈ M t k rt k (w; dt k ). (9b)

Note again that the natural parameters of r t k (w; d t k ) are estimated using Monte Carlo and the products in the above equations imply a summation of the natural parameters. In order to reduce the variance of this sum of estimators, we use Eq. (9a

) if |D t k | ≤ |M t k |, and Eq. (9b) if |D t k | > |M t k |.

Furthermore, we can compute the average bias from all natural parameter estimates (see App. C). We reduce the bias of our estimates by subtracting the average bias from all estimates. Note that a further option to update q t k (w) would be to use VB on the data

D t k ∪ M t k-1 \ M t k to compute the update q θt k (w) ≈ q θt k-1 (w) p(D t k ∪ M t k-1 \ M t k |w)

. The latter approach is numerically more stable but computationally more expensive. It also turned out that it is less favorable to the update using Eq. (9a) or Eq. (9b) in case of small datasets D t k , because VB applied to BNNs with small datasets often leads to a poor fit.

Previous work hypothesised that this problem is an artifact of the ELBO and not an optimisation problem (Trippe & Turner, 2018; Turner & M. Sahani, 2011) . We provide further evidence in Fig. 1 , where we infer the posterior of a Bayesian neural network with VB, using 70 and 100 data samples respectively and compare it to posterior inference with MCMC. In case of VB with 70 samples, the posterior approximation yields a model that is almost linear. These difficulties of posterior inference with variational Bayes are especially problematic in case of the streaming data setting, where the number of observations at each time-step is typically very small. The Gaussian update proposed above can alleaviate the problem of having to train BNNs with small datasets. Specifically, we have N t k + M instead of N t k data points to find a better optimum of the ELBO. (Figs. 1b, 1c) , the KL divergence of the ELBO is annealed from β = 0 to β = 1 over many iterations (450k annealing, 50k ELBO). Fig. 1d shows that the approximation trades off the expected log-likelihood for a better KL divergence as β is increased. With 70 data points, the annealed KL jumps to a significantly lower value, resulting in an almost linear decision boundary. By contrast, MCMC yields a much better predictive distribution for the same number of samples. Data is visualised in red and blue.

4. VARIATIONAL BAYES WITH MODEL ADAPTATION:

The incremental learning methods discussed so far assume i.i.d. data (cf. Sec. 2, Eq. ( 1)). This assumption can be reasonable even in scenarios with changing data distributions, e.g. when the data drift is an algorithmic artifact rather than a real phenomenon. For example, in online multi-task or curriculum learning we want to learn a model of all tasks, but we may choose to learn the tasks incrementally for various reasons (e.g. Nguyen et al., 2018; Kirkpatrick et al., 2017; Schwarz et al., 2018; Rusu et al., 2016; Fernando et al., 2017) . However, such approaches are not applicable for modeling non-stationary data: one of the properties of online VB is that the variance of the Gaussian posterior approximation shrinks at a rate of O(N ), where N is the total amount of data (e.g. Opper, 1998) . Consequently, learning comes to a halt as t → ∞. To overcome this issue, the model needs to be extended by a method that enables it to adapt to changes in the data distribution, e.g., by deliberately forgetting the belief inferred from previous data.

In the following, we describe two alternative methods for adapting to changing data. In Sec. 4.1, we impose Bayesian exponential forgetting, which forgets previous data exponentially by weighting the likelihood terms (or their approximations). In Sec. 4.2, we implement the adaptation through a diffusion process applied to the neural network parameters. Compared to the online learning scenario, we make the following assumptions: (i) we observe datasets D t k at potentially non-equidistant time steps t k ; (ii) data within D t k is assumed i.i.d., however, not between different datasets D t k and D t k+1 .

In both approaches, we realise adaptation by an additional forgetting step before observing the new data D t k+1 . We denote the distribution, which results from applying the forgetting step to the posterior approximation q θt k (w) p(M t k |w) by p t k+1 (w).

4.1. ADAPTATION WITH BAYESIAN FORGETTING:

Model adaptation through forgetting can be achieved by decaying the likelihood based on the temporal recency of the data (Graepel et al., 2010; Honkela & Valpola, 2003) . It has been explored previously as an alternative to filtering and is referred to as Bayesian exponential forgetting (Kulhavý & Zarrop, 1993) . This approach defines a forgetting operator that yields p(w t k+1 |D t1:t k ) directly. Here, we use a continuous-time version of this forgetting operation that can be formulated as

EQUATION

)

where τ is a time-constant corresponding to the average of the time-lags ∆t k+1 = t k+1 -t k . The distribution defined in Eq. ( 10) can be formulated recursively (cf. App. F) as

EQUATION

This equation can be viewed as Bayes rule (Eq.( 1)) applied after the forgetting step. The first two terms of Eq. ( 11) can be identified as the forgetting operation, applied to the current posterior. In order to apply this operation to our posterior approximation q θt k (w) p(M t k |w), we modify it by an additional weighting factor for each likelihood term in the memory. Denoting the age of a memory item m by ∆t k (m), the forgetting operation for this new posterior approximation then results in

EQUATION

where ∆t k+1 (m) = ∆t k (m) + ∆t k+1 . As can be seen from Eq. ( 12), BF acts on both factors of the posterior approximation independently: in case of the memory, it re-weights the respective likelihood terms by updating ∆t k+1 (m). For the Gaussian term q θt k (w), BF leads to a weighted product with the prior distribution (i.e. the first two terms of Eq. ( 12)), resulting in a Gaussian with parameters

σ -2 t k+1 = 1 -(1 -) ∆t k+1/τ σ -2 0 + (1 -) ∆t k+1/τ σ -2 t k , σ -2 t k+1 µ t k+1 = 1 -(1 -) ∆t k+1/τ σ -2 0 µ 0 + (1 -) ∆t k+1/τ σ -2 t k µ t k .

For ∆t k+1 → ∞, the likelihood term in Eq. ( 12) converges to the uniform distribution and the Gaussian term reverts to the prior. We note, however, that while Eq. ( 11) is an exact recursive form of Eq. ( 10), the online VB approximation of Eq. ( 11) is not generally identical to the (offline) VB approximation of Eq. ( 10) due to its successive approximations. For tuning the hyperparameter , we note that the weighting of likelihood terms corresponds to an effective dataset size of 1 / • N (if all datasets are of equal size N ). In Fig. 2 , we also visualise the forgetting operation applied to the Gaussian part of the posterior approximation for multiple values of .

4.2. ADAPTATION WITH DIFFUSION PROCESSES:

Model adaptation can also be realised by using dynamic model parameters that evolve according to a stochastic process. In this case, adaptation is achieved by the stochastic transition p t k+1 ,t k (w |w) resulting in a prediction distribution

EQUATION

where we consider Gaussian transitions p t k+1 ,t k (w |w). However, this operation is generally not tractable for our posterior approximation q θt k (w) p(M t k |w). Moreover, the forgetting operation implied by the transition does not retain the product form as in the case of BF. For this reason, we consider only a Gaussian posterior approximation (without memory) for this approach, that is

p t k+1 (w ) = p t k+1 ,t k (w |w) q θt k (w) dw.

As mentioned in Sec. 4.1, BF yields the prior distribution for ∆t k+1 → ∞. This is a desirable property, since it corresponds to forgetting all information conveyed by the data. In case of a Gaussian prior, the only Gaussian process that fulfills this requirement is the Ornstein-Uhlenbeck (OU) process given by the stochastic differential equation

dw t = θ • (µ 0 -w t ) dt + σ 0 √ 2θ dW t ,

where θ is the stiffness parameter which controls the drift rate towards µ 0 . To decouple the adaptation parameter from the rate at which data is observed, we rescale the stiffness parameter as θ = a /τ. The resulting prediction distribution p t k+1 (w) = N µ t k+1 , σ 2 t k+1 is defined by the parameters

µ t k+1 = 1 -e -a ∆t k+1 τ µ 0 + e -a ∆t k+1 τ µ t k , σ 2 t k+1 = 1 -e -2a ∆t k+1 τ σ 2 0 + e -2a ∆t k+1 τ σ 2 t k .

An interesting observation is that both parameters evolve independently of each other. In contrast to BF, the mean and variance-instead of the natural parameters-follow an exponential decay. The hyperparameter a can be determined e.g. through the half-time of the exponential decay of the mean parameter, given as τ 1/2 = 1 /θ. We visualise the time evolution of the above parameters in Fig. 2 .

5. RELATED WORK:

There are many Bayesian approaches to online learning, which differ mostly in the approximation of the posterior distribution at each time-step. Sequential Monte Carlo (Liu & Chen, 1998) approximates the posterior by a set of particles. Assumed Density Filtering (ADF) (Maybeck, 1982) and Bayesian online learning (Opper, 1998) 2018) use a single (block-diagonal) posterior approximation, summarising all previous time-steps. The latter method is closer to Bayesian online inference, as it is an approximation of Eq. (1). Our work is based on online VB (Opper, 1998; Ghahramani, 2000; Sato, 2001; Broderick et al., 2013) , which approximates the posterior at every time-step by minimising the KL-divergence between a parametric (here Gaussian) and the true posterior distribution. In contrast to online VB, we approximate the posterior by a Gaussian distribution and a running memory.

Other approaches are based on various types of episodic memory, motivated by their empirical success in preventing catastrophic forgetting. The basic idea of rehearsal (Ratcliff, 1990) is to train on both the new data and a subset of previous data or pseudo samples (Robins, 1995; Shin et al., 2017; Kemker & Kanan, 2017) sampled from a generative model. The memory-based online inference methods most similar to our approach are VCL (Nguyen et al., 2018) and VVM (Minka et al., 2009) . Both methods use a Gaussian distribution and a running memory to approximate the posterior. VCL uses heuristics such as random selection or the k-center method to update the memory. However, both heuristics select the memory independently of the Gaussian approximation. By contrast, VVM updates the memory with data that cannot be well approximated by the Gaussian distribution. VVM uses expectation propagation for the posterior approximation in a logistic regression model and, therefore, it is not directly applicable to our work. We transferred the main idea of VVM to online VB and developed the corresponding memory update method. In our case, the memory is updated with data for which the ELBO changes most if the corresponding likelihood functions are approximated by a Gaussian. In contrast to these two approaches, we extend our model by an adaptation method that allows to cope with non-stationary data.

Many adaptation methods were developed in the context of concept drift, however, few of these approaches operate in the Bayesian framework. For example McInerney et al. (2015) treat the learning dynamics of their model as a non-stationary process that allows for adaptation. In contrast, our approach uses an evolving prior and a well defined forgetting mechanism that gives a better control over the learning process. A more closely related approach uses the extended Kalman filter to estimate the optimal parameters of a logistic regression classifier Su et al. (2008) . However, they consider a transition model which is equivalent to a Wiener process (in unit-time) and therefore does not revert to the prior. By contrast, our approach (Sec. 4.2) models the dynamics as a prior-reverting OU process. BF (Kulhavý & Zarrop, 1993) has been applied as an alternative to adaptation with an explicit transition model (e.g. Honkela & Valpola, 2003; Graepel et al., 2010) . Compared to previous work, we used a continuous-time version of BF and extended it to our posterior appoximation consisting of a Gaussian and a running memory. -3c ) of a BNN (gray) and running memory (rectangular shapes, size is proportional to the score), chosen by the memory update proposed in Sec. 3.2. Data from t k and t <k-1 is visualised as large circles and small dots, respectively. Fig. 3d shows the one-step-ahead (predictive) LML (divided by the number of samples) for data that will be selected for the memory and data that will be evicted. Data that will be selected in the memory tends to have a significantly lower predictive likelihood. 

6. EXPERIMENTS:

We validate our proposed inference methods in two stages. In Sec. 6.1, we compare our memory update and Gaussian update (Sec. 3) to existing memory-based online inference methods on several standard machine learning datasets. In Sec. 6.2, we evaluate our adaptation methods (Sec. 4) on commonly used datasets with concept drift (Widmer & Kubat, 1996) , where the conditional distribution of labels given the features changes over time (i.e. non-stationary data in the context of predictive models).

We found that training (variational) Bayesian neural networks on streaming data is challenging, specifically, our approach requires model parameters very close to a local optimum since Eqs. (5a) and (5b) hold only at local extrema of the ELBO. To overcome these difficulties, we use several methods to reduce the variance of the gradient estimates for learning: (i) we apply the local reparametrisation trick (Kingma et al., 2015) ; (ii) we use the Adam optimiser (Kingma & Ba, 2014) ; and (iii) we use multiple Monte Carlo samples to estimate the gradients (cf. Tab. 2 for details). Furthermore, we developed methods for determining hyperparameters of the Gaussian prior and the initialisation distribution of Bayesian neural networks. The idea is similar to the initialisation method proposed by Glorot & Yoshua Bengio (2010) and He et al. (2015) : we choose the prior and the posterior initialisation such that the mean and standard deviation of the activations in every layer are approximately zero and one, respectively. We refer to App. H and App. I for a derivation and further details.

We use the following metrics for evaluation: (i) the avg. test log-marginal likelihood (LML)

N -1 test n log Eq θ t k (w) p(d (n) test |w) , where d (n)

test is a sample from a heldout test dataset; (ii) the avg. onestep-ahead LML N -1

t k+1 n log Eq θ t k (w) p(d (n) t k+1 |w) , where d (n) t k+1 is data observed at time-step t k+1 .

Both metrics measure the predictive performance, however (i) can be used in the online setting, where the data is i.i.d.; and (ii) is typically used to evaluate models with non-stationary streaming data.

6.1. ONLINE LEARNING:

In this section, we evaluate our running memory (Sec. 3) in an online learning setting. To illustrate how our memory update works, we start our evaluation with a qualitative assessment: we train a model on 2-dimensional toy data (two-moons), where we can visualise the selected memory. The BNN has 2 layers with 16 units and tanh activations, and has a prior p 0 (w) = N (w; 0, 1) on all weights and biases. The memory-size is M = 30. The model observes 150 data samples at time-step t 0 and 15 samples at all consecutive time-steps. In Fig. 3 , we visualise the selected memory and the corresponding scores for time-steps t 0 , t 25 , and t 50 , respectively. We can make the empirical observation that our method favors data close to the decision boundary. Furthermore, in Fig. 3d , we visualise the one-step-ahead LML for data that will be selected and evicted (in the next time-step), respectively. This shows that our memory update tends to select data for which the model has a low predictive LML. These observations support our intuition that the memory is indeed complementary to the Gaussian approximation, selecting data for which the likelihood cannot be well approximated by a Gaussian function. In Fig. 9 of the supplementary material, we visualised the running memory for a model trained on MNIST, showing that the memory also accumulates diverse samples over time.

We evaluate our memory-based online inference method (Sec. 3) quantitatively on several standard machine learning datasets, including regression (UCI Boston, UCI Concrete, UCI Energy, UCI Yacht) and classification (MNIST, UCI Spam, UCI Wine) tasks. Here, we refer to our approach as Gaussian Residual Scoring (GRS). We compare GRS to the respective memory update and Gaussian update methods proposed in VCL (Nguyen et al., 2018) (cf. Sec. 2.2). Refer to App. B for an explanatory list of compared update methods. Online learning is performed by observing N t k samples per time-step (cf. Tab. 2 for the experiment setup and hyperparameters.). For evaluation, we use a random held-out test dataset (20% of the data). We perform each experiment with 16 different random data splits and random seeds for the model parameter initialisation. In Fig. 4 , we plot the test LML, averaged over the 16 runs, against the memory size, and the LML over all time-steps. In most cases, random selection and the k-center method start with a worse initial fit at t 0 . This is because these methods perform the initial Gaussian update by optimising the ELBO with N t0 -M samples at t 0 ; by contrast, GRS uses a Gaussian update that first optimises the ELBO with N t0 samples and subsequently discounts the contribution of the memory. In Tab. 1, we report the mean and std. deviation of the LML, where the mean and std. deviation are taken over the 16 independent runs, each averaged over the last 10% time-steps. The results demonstrate the superior predictive performance of our update methods. We also note that the experiments on the smaller datasets (cf. Tab. 2 in App. B) result in a high variance among the random data splits and random seeds. This is the case for all compared methods and it could not be remedied e.g. by using annealing or a different prior.

6.2. ADAPTATION:

In this section, we evaluate our adaptation methods (Sec. 4) in settings with concept drift. We begin with a simple logistic regression problem, where the data

D t k = {(x t k , y t k )} n , x t k ∈ R 2 , y t k ∈ R is sampled from x t k ∼ Uniform(-3, 3), y t k ∼ Bernoulli(σ(w t k x t k )

). The true model has two time-dependent parameters w 0 t k = 10 sin(α • t k ), w 1 t k = 10 cos(α • t k ), where α = 5 deg /sec and where we observe data at t k ∈ [0, 1, . . . , 720]. Fig. 5 shows the learned model parameters for standard online learning (without adaptation), OU process transitions, and Bayesian forgetting. If the time-dependence of the data is ignored (in case of online VB), the class labels are distributed with equal probability in the whole input space. Consequently, as t → ∞, the weights of the model without adaptation shrink to 0. By contrast, the posterior means of BF and the OU process follow a sinusoidal curve as the parameters of the true model. We also evaluate our adaptation methods quantitatively on 3 datasets with concept drift (Weather, Gas Sensor Array Drift, Covertype). We compare online VB (without adaptation), the Wiener process (a special case of the OU process), the OU process, and Bayesian Forgetting (with and without memory). All compared variants use the same model architecture and hyperparameters (cf. Tab. 2 in the supplementary material). We report the one-step-ahead LML, where the expectation is approximated with 500 Monte Carlo samples. Results are averaged over the last 50% time-steps, because we are interested in the continual learning performance, and the first few time-steps will be similar for most methods. We report the mean and std. deviation over 8 independent runs with different random seeds. In Fig. 6 (and Fig. 10 in the appendix), we plot the LML against 10 adaptation parameter values (of the respective adaptation method), where the value zero corresponds to online VB. The LML for BF with different memory sizes and a fixed forgetting rate = 0.11 is shown in Fig. 6 . As can be seen from the results, all adaptation methods significantly improve the performance compared to online VB. Interestingly, the Ornstein-Uhlenbeck process performs better than Bayesian Forgetting, however, using a running memory with Bayesian Forgetting closes the gap.

7. CONCLUSION:

In this work, we have addressed online inference for non-stationary streaming data using Bayesian neural networks. We have focused on posterior approximations consisting of a Gaussian distribution and a complementary running memory, and we have used variational Bayes to sequentially update the posteriors at each time-step. Existing methods update these two components without having an interaction between them, and they lack methods to adapt to non-stationary data. We have proposed a novel update method, which treats both components as complementary, and two novel adaptation methods (in the context of Bayesian neural networks with non-stationary data), which gradually revert to the prior distribution if no new data is observed.

Future research could extend our work by drift detection methods and use them to infer the adaptation parameters. This work could also be extended by developing adaptation methods for gradual, abrupt, or recurring changes in the data distribution. Finally, we observed that variational Bayesian neural networks with a uni-modal approximate posterior often find poor local minima if the dataset is small and models are complex. This is especially challenging in scenarios with streaming data. While our Gaussian update alleviates this problem to a certain degree, further research in extending the approximation family beyond Gaussians could be beneficial. Progress in this direction would improve our proposed methods and allow to scale them to more complex models.

A FURTHER EXPERIMENTAL RESULTS:



A.1 MEMORY:

Here we provide additional experimental results for the memory update and Gaussian update from Sec. 3. We conducted experiments on 3 additional datasets (UCI Boston, UCI Yacht, UCI Red Wine). The influence of the memory size and the performance over time (for a specific memory size) are shown in Fig. 7 (corresponding to Fig. 4 in the main text). Furthermore, we test the memory update and Gaussian update of GRS separately on UCI Energy and UCI Concrete. For this purpose, i) we combine the k-center method with our Gaussian update from Sec. 3.3; and ii) we use our memory update from Sec. 3.2 and update the Gaussian distribution by optimizing Eq. ( 2) with D t k ∪ M t k-1 \M t k (re-fitting) . The results are shown in Fig. 8 . As can be seen, GRS performs better than one of the components used in combination with a baseline method. GRS with refit performs especially bad, similar to the baselines k-center and random. This is because refitting requires optimising the ELBO with a small dataset. As mentioned in Sec. 3.3 (cf. Fig. 1 ), Bayesian neural networks with VB perform bad on small datasets due to over-regularisation. Consequently, in case of refitting, a good memory update can lead to a worse overall performance due to a much worse Gaussian update. While this general issue with Bayesian neural networks (learned with VB) is beyond the scope of this work, it is an important future research direction. (bottom) . GRS denotes our approach (Sec. 3), GRS with k-center replaces our memory update by the k-center method; GRS with refit replaces our Gaussian update by the optimization of Eq. ( 2) with Dt k ∪ Mt k-1 \Mt k . Evaluated for several memory sizes (left), and evaluated over time (right) for 3 different memory sizes. Hyperparameters are chosen as in Sec. 6.1.

To better understand our memory update using the score function from Eq. ( 8), we visualise the running memory for a model trained on MNIST in Fig. 9 . We also evaluated our adaptation methods on 2 additional datasets (Gas Sensor Array Drift, Weather). In Fig. 10 , we visualise the influence of the adaptation parameter for these datasets. Note that the range of the adaptation parameters is on a much smaller range compared to the experiments on Covertype (Sec. 6.2). For larger values, the performance starts to degrade. Surprisingly, the memory degrades the performance in case of the Gas Sensor Array Drift dataset.

A.3 CATASTROPHIC FORGETTING WITH ONLINE VB AND BAYESIAN NEURAL NETWORKS:

Here we provide further experimental results for the behavior of online VB (Secs. 2.1, 3) in case of non-stationary data. For this purpose, we train Bayesian neural networks with different architectures on the toy classification problem with a rotating decision boundary from Sec. 6.2, however, with 150 data samples per time-step. In Fig. 11 , we visualise the training LML for different architectures, including a linear model. It can be seen that Bayesian neural networks with higher complexity (i.e. more layers or more units) drop slower in performance compared to the linear model. However, this is not a desired property for online VB, since exact online Bayesian inference would yield the same posterior distribution as offline Bayesian inference. In case of our toy classification data (where the time dependence is ignored), online inference should not be able to classify the data as t → ∞. Instead, this learning behavior shows that online VB with Gaussian approximate posterior distributions is prone to catastrophic forgetting. 

B EXPERIMENT SETUP:

The following is an explanatory list of the update methods used in Sec. 6.1:

• k-center (VCL): Uses the k-center method (Sec.2.2) for the memory update and Eq. ( 2) with (D t k ∪ M t k-1 ) \ M t k for the Gaussian update.

• random (VCL): Uses random selection (Sec.2.2) for the memory update and Eq. ( 2) with (D t k ∪ M t k-1 ) \ M t k for the Gaussian update.

• GRS (Gaussian Residual Scoring, ours): Uses Eq. ( 8) for the memory update (Sec. 3.2) and performs the Gaussian update by first using Eqs.

(2) with (D t k ∪ M t k-1 ) and subsequently using Eqs. (5a), (5b) for removing the local contributions of M t k (cf. Sec. 3.3).

Similarly, the following list summarises the adaptation methods used in Sec. 6.2:

• Wiener process: Posterior approximation consists of q θt k (w) only. Transition p(w t k+1 |w t k ) is given by a random walk. We used a diffusion that is proportional to the prior standard deviation in every neural network layer (cf. Sec. 4.2). No memory used.

• Ornstein-Uhlenbeck process: Posterior approximation consists of q θt k (w) only. Transition p(w t k+1 |w t k ) is given by the Ornstein-Uhlenbeck process (cf. Sec. 4.2). No memory used.

• Bayesian forgetting: Posterior approximation consists of q θt k (w) only. No state-space model assumption, instead uses Bayesian exponential forgetting (cf. Sec. 4.1).

• Bayesian forgetting with memory: Posterior approximation consists of q θt k (w) and M t . No state-space model assumption, instead uses Bayesian exponential forgetting (cf. Sec. 4.1).

In Tab. 2, we summarise experimental setup (hyperparameters) used for Secs. 6.1 and 6.2.

C FACTORISATION PROPERTY OF THE GAUSSIAN VARIATIONAL APPROXIMATION:

Here we derive the factorisation property of the Gaussian variational approximation distribution by expressing the natural parameters of the Gaussian approximation as a sum. This can be shown for the Gaussian approximation at a local optimum of the ELBO. For a Gaussian prior and posterior the ELBO is given as Hence, we obtain

L(µ * , Σ * ) = - 1 2 log |Σ 0 | -log |Σ * | -d + (µ * -µ 0 ) T Σ -1 0 (µ * -µ 0 ) + T r(Σ * Σ -1 0 ) + N n=1 E w∼q θ * (w) log p(d (n) |w) .

EQUATION

Similarly, we have ∂L(µ * ,Σ * )

EQUATION

Hence, we obtain

EQUATION

)

Next, we calculate the natural parameters from Eqs. ( 14), ( 16): n) .

EQUATION

Monte Carlo estimation: The natural parameters Λ (n) , η (n) can be estimated with Monte Carlo, by replacing the expectation with an empirical mean. Since the parameters Λ * and Λ 0 (and η * , η 0 respectively) are known, the total bias of the parameter estimates can be computed: n) .

EQUATION

We use this to reduce the bias for the individual terms:

Λ(n) = Λ (n) - 1 N Λ b , η(n) = η (n) - 1 N η b .

D ELBO IN RESIDUALS FORM:

Here we show how the ELBO can be written in the form of Eq. ( 6). Let us define the variational distribution in the factorised form q θ (w) = Z -1 q p(w) N n=1 r (n) (w) (cf. Sec.3.1). We can then write the ELBO as

L(µ, Σ; D) = E q θ (w) N n=1 log p(d (n) |w) + log p(w) -log q θ (w) = E q θ (w) N n=1 log p(d (n) |w) + log r (n) (w) -log r (n) (w) + log p(w) -log q θ (w) = E q θ (w) N n=1 log p(d (n) |w) -log r (n) (w) + log p(w) N n=1 r (n) (w) q θ (w) = E q θ (w) N n=1 log p(d (n) |w) -log r (n) (w) + log Z q • q θ (w) q θ (w) = log Z q + n E q θ (w) log p(d (n) |w) -log r (n) (w) .

E MEMORY UPDATE SCORE FUNCTION:

In Eq. 8, the expectation involving Gaussian terms can be calculated analytically:

E qθ t k (w) log r (m) t k (w) = qθt k (w) η (n) w - 1 2 Λ (n) w 2 dw = η (n) µ (n) - 1 2 Λ (n) (µ * ) 2 + Σ * = η (n) (Λ * ) -1 η * - 1 2 Λ (n) (Λ * ) -1 η * 2 - 1 2 Λ (n) (Λ * ) -1

The expectation involving non-Gaussian terms (in Eq. 8) has no closed-form solution. We therefore estimate E qθ t k (w) log p(d

t k |w) using Monte-Carlo.

F BAYESIAN FORGETTING -RECURSIVE FORMULATION:

Here we show how Bayesian forgetting can be rearranged into a recursive formulation. We first bring this formula into a similar form as Eq. ( 1), extracting the most recent likelihood term:

p(w|D t1:t K+1 ) ∝ p 0 (w) • K+1 k=1 p(D t k |w) (1-) t K+1 -t k τ = p 0 (w) • K k=1 p(D t k |w) (1-) t K+1 -t k τ • p(D t K+1 |w).

The first two terms can be rewritten as

p 0 (w) • K k=1 p(D t k |w) (1-) t K+1 -t k τ = p 0 (w) • K k=1 p(D t k |w) (1-) t K+1 -t K +t K -t k τ = p 0 (w) • K k=1 p(D t k |w) (1-) t K -t k τ •(1-) t K+1 -t K τ = p 0 (w) • K k=1 p(D t k |w) (1-) t K -t k τ (1-) t K+1 -t K τ ∝ p 0 (w) • p(w|D t1:t K ) p 0 (w) (1-) t K+1 -t K τ = p 0 (w) 1-(1-) t K+1 -t K τ • p(w|D t1:t K ) (1-) t K+1 -t K τ

Hence, we have shown that the posterior can be expressed recursively as p(w|D t1:t k+1 ) ∝ p 0 (w) 1-(1-) ∆t k+1/τ p(w|D t1:t k ) (1-) ∆t k+1/τ p(D t k+1 |w).

The parameters of the Gaussian part q θt k (w) of the posterior approximation (after applying the forgetting operation) can be calculated easily from the above equation.

Natural parameters:

Λ t k+1 = Λ 0 • (1 -(1 -) t k+1 -t k τ ) + Λ t k • ((1 -) t k+1 -t k τ ) η t k+1 = η 0 • (1 -(1 -) t k+1 -t k τ ) + η t k • ((1 -) t k+1 -t k τ )

Covariance parameter:

Σ t k+1 = Λ 0 • (1 -(1 -) t k+1 -t k τ ) + Λ t k • ((1 -) t k+1 -t k τ ) -1 = Σ -1 0 • (1 -(1 -) t k+1 -t k τ ) + Σ -1 t k • ((1 -) t k+1 -t k τ ) -1

Location parameter:

µ t k+1 = Σ t k+1 • η t k+1 = Σ t k+1 η 0 • (1 -(1 -) t k+1 -t k τ ) + η t k • ((1 -) t k+1 -t k τ ) = Σ t k+1 Σ -1 0 µ 0 • (1 -(1 -) t k+1 -t k τ ) + Σ -1 t k µ t k • ((1 -) t k+1 -t k τ ) G PSEUDO-ALGORITHM

We provide the pseudo algorithm of GRS (Sec. 3) with Bayesian forgetting in Alg. 1. The computational complexity (at each of the K time-steps) is dominated by i) the minimisation of the KL divergence, ii) estimating the Gaussian factors, and iii) scoring the memory. The KL minimisation requires I t k sequential iterations with N t k + M data samples and K train Monte Carlo samples. The latter can both be processed in parallel on parallel hardware. The estimation of the Gaussian factors requires N t k + M sequential iterations with K term parallel Monte Carlo samples. The dominating computation for calculating the scores is the evaluation of the likelihood for N t k + M data samples and K train Monte Carlo samples, both of which can be processed in parallel. The highest scoring candidate memory is given by the top-M highest scoring data points, thus, the computational complexity of Eq. ( 8) is only linear in the number of samples.

Algorithm 1 Gaussian Residual Scoring with Bayesian forgetting. The function EstimateGaussian-Factors corresponds to Eqs. (5a), (5b) (cf. also App. C). The function ApplyForgetting corresponds to Eq. ( 12). Note that p t k (w) includes the adapted likelihood of the memory and all subsequent functions involving the memory use this adapted likelihood.

Inputs: p 0 (w),

q θt 0 (w), τ , K for k in [0...K] do t k = GetTimeStamp(k) ∆t k = t k -t k-1 D t k = GetData(t k ) if k == 0 then p t k (w) = p 0 (w) else p t k (w) = ApplyForgetting p 0 (w), q θt k-1 , M t k-1 , ∆t k {Sec. 4.1} end if qθt k (w) = argmin q θ KL q θ (w) || Z-1 t k p t k (w) p(D t k |w) {Sec. 3.2, Sec. 4.1} {r t k (w; d t k )} dt k ∈Dt k ∪Mt k-1 = EstimateGaussianFactors qθt k (w), D t k , M t k-1 {Sec. 3.1} M t k = argmax M S t k (M; D t k ∪ M t k-1 ) {Sec. 3.2} if |D t k | ≤ |M t k | then q θt k (w) = p t k (w) dt k ∈ Mt k r t k (w; d t k ) {Sec. 3.3} else q θt k (w) = qθt k (w) / dt k ∈ Mt k r t k (w; d t k ) {Sec. 3.3} end if end for H PRIOR PARAMETERS

Here we develop a heuristic to choose the initial prior p 0 (w). As this will not be specific to the online or continual setting, we drop the time index in this section, denoting the prior as p(w). Furthermore, we consider only Gaussian distributions with a diagonal covariance matrix. Assume that the data is standardised, that is, the first two moments are zero and one. A reasonable choice for the prior parameters is such that the first two moments of the prior-predictive distribution equals the first two moments of the data distribution. We go one step further and constrain the pre-activations of every neural network layer to have moments zero and one. Denote all weight matrices and weight biases by w = {W l } l ∪ {b l } l , and let x 0 denote the input data. Let us further denote the pre-activation (before non-linearity) of layer l and unit i as follows.

x i l = j W i,j l f l-1 (x j l-1 ) + b i l .

The constraints are then given as follows.

E w∼p(w) E x0∼p(D) x i l = 0, E w∼p(w) E x0∼p(D) (x i l ) 2 = 1.

The first constraint can be easily fulfilled by setting the prior mean to zero for all parameters.

µ i,j l = 0.

This follows immediately from W l ⊥ ⊥ f l-1 (x l-1 ) and the expectation of products of independent random variables. The second moment can then be calculated as follows. • E w∼p(w) Var x0∼p(D) f l-1 (x j l-1 ) =:

N l-1 j E w∼p(w) W i,j l 2 • c f l-1 = c f l-1 • N l-1 j Var w∼p(w) W i,j l = N l-1 • c f l-1 • Var w∼p(w) W i,j l .

Here we introduced c f l-1 to denote a correction factor for the non-linearity f l-1 . In case of the linear function, we will have c f l-1 = 1. For arbitrary non-linearities, we can estimate this factor numerically, assuming that the pre-activations are distributed according to N (0, 1).

c f l-1 = Var x j l-1 ∼N (0,1) f l-1 (x j l-1 )

This can be done beforehand and the factors for common activation functions can be stored in a lookup table. Finally, plugging in the constraint for the second moment in the above equation, we obtain the following prior variance.

EQUATION

I POSTERIOR INITIALISATION:

It is known that a proper initialisation of standard neural networks is crucial for the optimisation process Glorot & Yoshua Bengio (2010) ; He et al. (2015) . In Bayesian neural networks, the matter becomes even more complicated, since we have to deal additionally with the variance of the Monte Carlo estimator due to re-parametrisation. Analogous to the choice of prior parameters, we seek a posterior initialisation that yields the first two moments of zero and one. A naive attempt would be to initialise the posterior with the prior parameters. However, the significant noise in the Monte Carlo estimation typically leads to bad optimisation results and even numerical instabilities. We propose an initialisation method which fulfills our constraints but allows us determine the variance of the initial posterior with two hyperparameters α and β.

Let us denote the mean and log-scale parameters of the approximate posterior as θ = {θ µ , θ log σ }.

We choose the following initialisation distributions.

q(w) = N (w; θ µ , e 2θ log σ ), where p(θ µ ) = N θ µ ; µ θµ , σ 2 θµ , and p(θ log σ ) = N θ log σ ; µ θ log σ , σ 2 θ log σ . Here and in the following, we dropped the time index for the approximate posterior, as well as the indices l, i, and j for the model parameters θ.

We follow a similar derivation as in Sec. H. As for the prior, the mean of the initialisation distribution will be zero for all parameters.

µ θµ = 0. For the second moment, the derivation is as follows. =:

N l-1 j E θ∼p(θ) θ 2 µ + e 2•θ log σ • c f l-1 = N l-1 • c f l-1 • E θ∼p(θ) θ 2 µ + e 2•θ log σ = N l-1 • c f l-1 • E θ∼p(θ) θ 2 µ + E θ∼p(θ) e 2•θ log σ = N l-1 • c f l-1 • E θ∼p(θ) θ µ 2 + Var θ∼p(θ) θ µ + e 2E[θ log σ ]+2Var[θ log σ ] = N l-1 • c f l-1 • µ 2 θµ + σ 2 θµ + e 2µ θ log σ +2σ 2 θ log σ = N l-1 • c f l-1 • σ 2 θµ + e 2µ θ log σ +2σ 2 θ log σ .

Hence, the second constraint is as follows.

1 N l-1 • c f l-1 = σ 2 θµ + e 2µ θ log σ +2σ 2 θ log σ .

In contrast to Sec. H, we are now under-constrained by 2 parameters. We therefore introduce two hyperparameters α and β. We first determine α := σ θ log σ , for which we generally choose small values α ≈ 0 (α = 0 corresponds to initialising all posterior variances in the given layer with the same value). The second hyperparameter β ∈ [0, 1] determines how much of the total variance is due to the variance of the location parameter and how much variance is due to the variance of the log-scale parameter. Inserting α and introducing β we obtain the following equations.

σ 2 θµ = β N l-1 • c f l-1

, and e 2µ θ log σ +2α 2 = 1 -β N l-1 • c f l-1 .

:

A strict definition of online learning requires single data samples at each time step instead of batches Dt k .

Figures and Tables:

Figure 1: Posterior predictive distribution in the xy-plane (grey) of a Bayesian neural network with 2 layers of 16 units, tanh activations, prior p0(w) = N (w; 0, 1), and Bernoulli likelihood. In case of variational Bayes (Figs. 1b, 1c), the KL divergence of the ELBO is annealed from β = 0 to β = 1 over many iterations (450k annealing, 50k ELBO). Fig.1dshows that the approximation trades off the expected log-likelihood for a better KL divergence as β is increased. With 70 data points, the annealed KL jumps to a significantly lower value, resulting in an almost linear decision boundary. By contrast, MCMC yields a much better predictive distribution for the same number of samples. Data is visualised in red and blue.

Figure 2: Time-evolution of distribution parameters of Bayesian Forgetting (top) and the Ornstein-Uhlenbeck process (bottom) for different adaptation parameter values. The initial distribution (at t = 0) can be seen as the approximate posterior at some time-step t k .

Figure 3: Two-moons dataset. Predictive distribution (Figs. 3a -3c) of a BNN (gray) and running memory

Figure 5: Mean and std deviation of the approximate posterior distributions of a logistic regression model over 720 time-steps. The model is trained on a toy classification problem with rotating class boundaries (cf. Sec. 6.2).Online VB (left) quickly converges to zero mean, whereas Bayesian forgetting and OU-process transitions lead to a sinusoidal curve as in the true model.

Figure 7: Average test LML on further datasets not included in the main text. Evaluated for several memory sizes (top), and evaluated over time (bottom) for a specific memory size. Cf. Sec. 6.1, App. A.1 for details.

Figure 8: Average test LML on UCI Energy (top) and UCI Concrete (bottom). GRS denotes our approach

Figure 9: Running memory at different time-steps on MNIST (cf. Sec. 6.1), with a memory size N = 100. The memory update tends to select non-typical data while showing diversity.

Figure 11: LML for toy classification problem with rotating class boundaries (cf.Sec. A)

At a local optimum, we have ∂L(µ * ,Σ * ) w∼q θ * (w) log p(d (n) |w) = Σ -1 0 (µ * -µ 0 ).

w∼p(w),x0∼p(D) (x i l ) 2 = E w∼p(w) E x0∼p(D) (x i l ) 2 = E w∼p(w) Var x0∼p(D) x i l + 0 = E w∼p(w) Var x0∼p(D)

θ∼p(θ),w∼q(w|θ),x0∼p(D) (x i l ) 2 = E θ∼p(θ) E w∼q(w|θ) E x0∼p(D) x i l E θ∼p(θ) E w∼q(w|θ) Var x0∼p(D) E θ∼p(θ) E w∼q(w|θ) N l-1 j (W i,j l ) 2 Var x0∼p(D) x j l-1 = E θ∼p(θ) E w∼q(w|θ) N l-1 j (W i,j l ) 2 • E θ∼p(θ) E w∼q(w|θ) Var x0∼p(D) f l-1 (x j lθ∼p(θ) E w∼q(w|θ) (W i,j l ) 2 • E θ∼p(θ) E w∼q(w|θ) Var x0∼p(D) f l-1 (x j l-1 )

r t k (w; d t k ) , where d t k are the samples in D t k ∪ M t k-1 and where r t k (w; d t k ) are the Gaussian approximation of the corresponding likelihood terms. Note that r t k does not only depend on d t k , however, we omit the dependence on the remaining data for notational convenience. If the likelihood p(d t k |w) is close to the Gaussian r t k (w; d t k ) in expectation w.r.t. the approximate posterior qθt k (w), then its contribution to L(θ t k ; D t k ∪ M t k-1) is small. Similarly, likelihood terms that cannot be well approximated by the respective Gaussian have a large contribution, and, hence, the corresponding data should be kept in the memory. For this reason, we propose the score function

Table content: <table/>

are deterministic posterior approximations based on moment matching. Other deterministic approaches are based on Laplace's approximation (MacKay, 1992): Kirkpatrick et al. (2017) use multiple diagonal Gaussian posterior approximations of previous time-steps to regularise future tasks; Ritter et al. (

Table content: <table/>



Table content: <table><tr><td/><td>GRS (ours)</td><td>k-center</td><td>random</td></tr><tr><td colspan="4">Concrete -0.779 ± 0.039 -0.798 ± 0.039 -0.800 ± 0.039</td></tr><tr><td>Boston</td><td colspan="3">-0.619 ± 0.111 -0.638 ± 0.093 -0.664 ± 0.156</td></tr><tr><td>Energy</td><td>0.365 ± 0.440</td><td colspan="2">-0.119 ± 0.128 -0.078 ± 0.087</td></tr><tr><td>Yacht</td><td>1.925 ± 0.229</td><td>1.658 ± 0.291</td><td>1.648 ± 0.254</td></tr><tr><td>Spam</td><td colspan="3">-0.216 ± 0.016 -0.219 ± 0.015 -0.217 ± 0.016</td></tr><tr><td>Wine</td><td colspan="3">-1.165 ± 0.056 -1.212 ± 0.059 -1.194 ± 0.070</td></tr><tr><td>MNIST</td><td colspan="3">-0.148 ± 0.005 -0.158 ± 0.005 -0.153 ± 0.005</td></tr></table>

Average test LML, evaluated for several memory sizes (top), and evaluated over time (bottom) for a specific memory size (cf. corresponding legend). Cf. Sec. 6.1 for details and App. A for further results.

Table content: <table><tr><td/><td/><td/><td/><td/><td/><td/><td/><td>UCI Energy</td><td/><td/><td/><td/><td/><td/><td>UCI Spam</td><td/><td/><td/><td>MNIST</td></tr><tr><td/><td/><td/><td/><td/><td>log-expected-likelihood</td><td>1 0.2 0.3</td><td/><td/><td colspan="2">GRS kcenter random</td><td>log-expected-likelihood</td><td>0.219 0.218 0.217 0.216</td><td/><td>GRS kcenter random</td><td/><td/><td/><td>log-expected-likelihood</td><td>0.165 0.160 0.155 0.150</td><td>GRS kcenter random</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td>0</td><td>10</td><td>20 memory size 30</td><td>40</td><td>50</td><td/><td>0</td><td>50</td><td>100 150 200 250 memory size</td></tr><tr><td>log-expected-likelihood</td><td>1.0 0.9 0.8</td><td>0</td><td>20 UCI Concrete 40 t</td><td>60 GRS 15 kcenter 15 random 15 no memory</td><td>log-expected-likelihood</td><td>0.4 0.2 0.0 0.2 0.4</td><td>0</td><td>10 GRS 15 kcenter 15 20 UCI Energy 30 t random 15 no memory</td><td>40</td><td>50</td><td>log-expected-likelihood</td><td>0.35 0.30 0.25</td><td>0</td><td colspan="4">25 50 75 100 125 t UCI Spam GRS 30 kcenter 30 random 30 no memory</td><td>log-expected-likelihood</td><td>0.275 0.250 0.225 0.200 0.175 0.150</td><td>0</td><td>t 20 40 60 80 100 MNIST GRS 150 kcenter 150 random 150 no memory</td></tr><tr><td colspan="3">Figure 4:</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr></table>

Experiment setup for online experiments. N t0 is the number of observed samples at the first time-step and N t 1:k is the dataset size of all other time-steps. M refers to the number of samples in the memory. We evaluated 5 different memory sizes for each experiment in Sec. 3 and 10 sizes for experiments in Sec. 4, where the sizes are equally spaced in the given range. K train and K term is the number of MC samples used for training and for estimating the Gaussian terms respectively. I t0 and I t 1:K refer to the number of iterations for the first time-step and all subsequent time-steps, respectively. The architecture denotes the number of units for each hidden layer (e.g.[16, 16]  denotes 2 hidden layers with 16 units each).

Table content: <table><tr><td/><td>Nt 0</td><td colspan="2">Nt 1:k M</td><td colspan="3">Ktrain Kterm It 0</td><td colspan="2">It 1:K Architecture</td></tr><tr><td>Concrete</td><td>100</td><td>10</td><td>[5 . . . 25]</td><td>1k</td><td>50k</td><td colspan="2">50k 10k</td><td>[16, 16]</td></tr><tr><td>Boston</td><td>100</td><td>10</td><td>[5 . . . 25]</td><td>1k</td><td>50k</td><td colspan="2">50k 10k</td><td>[16, 16]</td></tr><tr><td>Energy</td><td>100</td><td>10</td><td>[5 . . . 25]</td><td>1k</td><td>50k</td><td colspan="2">50k 10k</td><td>[16, 16]</td></tr><tr><td>Yacht</td><td>100</td><td>10</td><td>[5 . . . 25]</td><td>1k</td><td>50k</td><td colspan="2">50k 10k</td><td>[16, 16]</td></tr><tr><td>Spam</td><td>250</td><td>25</td><td>[10 . . . 50]</td><td>400</td><td>50k</td><td colspan="2">50k 10k</td><td>[32, 32]</td></tr><tr><td>Wine</td><td>250</td><td>25</td><td>[10 . . . 50]</td><td>400</td><td>50k</td><td colspan="2">50k 10k</td><td>[32, 32]</td></tr><tr><td>MNIST</td><td colspan="2">2.5k 250</td><td colspan="2">[50 . . . 250] 40</td><td>50k</td><td colspan="2">50k 10k</td><td>[64, 64]</td></tr><tr><td>GasSensorArray</td><td>500</td><td>50</td><td>[5 . . . 50]</td><td>200</td><td colspan="3">100k 50k 10k</td><td>[8, 8]</td></tr><tr><td>Weather</td><td>1k</td><td>100</td><td>[5 . . . 50]</td><td>100</td><td colspan="3">100k 50k 10k</td><td>[16, 16]</td></tr><tr><td>Covertype</td><td>1k</td><td>1k</td><td colspan="2">[10 . . . 200] 400</td><td colspan="3">100k 50k 10k</td><td>[32, 32]</td></tr></table>

Peer review (iclr2020):
Contributions: The paper extends variational continual learning with memory. In this setting the posterior distribution of the model parameters is approximated using a mean-field Gaussian approximation and a small set of datapoints is kept in a memory to combat catastrophic forgetting. - The paper proposes a new rule for updating the memory and the Gaussian approximation after examining each batch of points. For the memory update, instead of the previously used k-means or random sampling, it examines the factor *r(d|w)* that each datapoint contributes and selects the most significant ones. - A second contribution is the use of two adaptation methods in the online learning setting in case there is a distribution shift in the data: Bayesian forgetting and the Ornstein-Uhlenbeck process. (1) The first contribution, the new memory update rule and the Gaussian update, is novel to my knowledge. The idea is to calculate the Gaussian factor that each datapoint contributes to the posterior (the update corresponding to each datapoint in the message-passing interpretation) and the select the points that contribute the most. Effectively, it calculates the change in the ELBO if each candidate datapoint was moved to the memory and selects the set that minimizes the ELBO of the remaining points. This is a sensible heuristic because we want to keep datapoints in the memory that contribute the most to the posterior approximation. If a datapoint contributes nothing to the posterior then it shouldn*t be kept in memory. The idea clever and it seems to work well in practice. Questions: - The Gaussian that each datapoint contributes is an approximation in the sense that if we take that datapoint and its contribution away, the remaining approximate posterior might be suboptimal. This has the consequence that S_tk (eq. 8) is only an approximation to what the optimal ELBO would be without the datapoints selected for the memory. Did the authors conduct experiments or have thoughts on how well S_tk approximates the ELBO of the variational approximation of the remaining points? -Regarding the experiments, there is a comparison to k-centre and random selection, both of which are proposed in VCL. I am very surprised to see in the experiments that both of these are outperformed by *no-memory*. For example on *energy*, the baselines at t50 are outperformed by the non memory model at t0. Can this all be contributed to the baselines training their variational approximation without the memory at t0? - It is unclear whether the gains can be contributed to the new Gaussian update rule or the new memory update rule. To see the contribution of each, there should be an experiment where the proposed Gaussian update rule is used along with k-centre for the memory update. The fact the the no-memory model performed to close to GRS suggests that k-centre with the new Gaussian update rule would be even closer to it. (2) The second contribution of the paper is the use of Bayesian forgetting and OU to use to deal with the data distribution shift. It shows how these two approaches can be used along with the proposed Gaussian and memory updates. The paper is already very long, so there likely isn*t any space for fleshing out this section, but it would be nice to have the experimental results included in the main paper, because most of the experiments are left to the appendix. Perhaps it would be a good idea to focus only on Bayesian forgetting (or the OU process) and try to shorten the section a bit. In terms of writing and clarity, I have no complaints. The paper is well written and easy to understand. Minor: - 3.2 *If the likelihood term p(dtk|w) is well approximated by r(w;dtk)* - I find this sentence a bit confusing. What does it mean for p(dtk|w) to be well approximated by r(w;dtk)? - 3.3 *In order to reduce the variance of the Monte-Carlo approximation*. What is the MC approximation made here? - The method never actually uses the assumption that the datapoints are sampled i.i.d.. The algorithm should still work if the datapoints are not examined in a random order, e.g. consider seeing all the images with label *0* and then all the images with label *1* etc.. This would likely degrade the performance significantly but the method should still work. Overall assessment: Pros: I like the ideas in the paper and they are presented well. Cons: The paper is a bit too long. The experiments could more thoroughly investigate the source of the gains.