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

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.

heading: 1 INTRODUCTION
text: 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 general-
ising 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. ∗Correspondence to richard.kurle@tum.de

heading: 2 BACKGROUND: ONLINE INFERENCE
text: Consider a stream of datasets {Dtk}Kk=1, where tk are the time points at which datasets Dtk 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|Dt1:tk) 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|Dt1:tk) ∝ p(w|Dt1:tk−1) p(Dtk |w,Dt1:tk−1) = p(w|Dt1:tk−1) p(Dtk |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|∅) = p0(w) = N (w;µ0, σ0). Furthermore, we consider supervised learning, where Dtk = {d (n) tk }n = {(x(n)tk ,y (n) tk )}n and p(d(n)tk |w) = p(y (n) tk |NN ( x (n) tk ;w) ) .

heading: 2.1 ONLINE VARIATIONAL BAYES
text: 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θtk (w) ≈ p(w|Dt1:tk) through a sequence of projections:
qθtk (w) = argmin qθ
KL [ qθ(w) ||Z−1tk qθtk−1 (w) p(Dtk |w) ] , (2)
where Ztk is the normalisation constant. The above minimisation is equivalent to maximising the evidence lower bound (ELBO) Ltk(θ;Dtk) = Eqθ(w) [ log p(Dtk |w) ] −KL [ qθ(w) || qθtk−1 (w) ] . In this work, we consider diagonal Gaussian posterior approximations qθtk (w) for the neural network weights, similar to Nguyen et al. (2018).

heading: 2.2 ONLINE VARIATIONAL BAYES WITH MEMORY
text: 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|Dt1:tk) ≈ qθtk (w) p(Mtk |w) (3) and update them sequentially as new dataDtk is observed. The factor p(Mtk |w) = ∏M m p(m (m) tk |w) is the likelihood of a set of M = |M| data points, which we refer to as running memory; and qθtk (w) is a Gaussian distribution, which summarises the rest of the data D̄1:tk = D1:tk\Mtk . 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 memoryMtk ⊂ Dtk ∪ Mtk−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 Dtk ∪ Mtk−1 .); (ii) the Gaussian distribution is updated with the remaining data D̄tk = Dtk ∪ Mtk−1 \ Mtk (using Eq. (2)) to obtain qθtk (w) ≈ qθtk−1 (w) p(D̄tk |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
q̃θtk (w) = argmin qθ
KL [ qθ(w) || Z̃−1tk qθtk (w) p(Mtk |w) ] . (4)
This distribution should not be confused with the recursively updated variational distribution (Eq. (2)). 1 A strict definition of online learning requires single data samples at each time step instead of batches Dtk .

heading: 3 IMPROVING MEMORY-BASED ONLINE VARIATIONAL BAYES
text: 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 D̄tk—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.

heading: 3.1 PROPERTIES OF THE GAUSSIAN VARIATIONAL APPROXIMATION
text: 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 p0(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):
Σ−1µ = Σ−10 µ0 + ∑ n ( ∂µEqθ(w) [ log p(d(n)|w) ] − 2∂ΣEqθ(w) [ log p(d(n)|w) ] µ ) , (5a)
Σ−1 = Σ−10 − 2 ∑ n ∂ΣEqθ(w) [ log p(d(n)|w) ] . (5b)
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) = Zq−1p0(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 Zq = ∫ p0(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) (c.f. App. D)
L(µ,Σ;D) = logZq + ∑ n Eqθ(w) [ log p(d(n)|w)− log r(n)(w) ] . (6)
If the terms p(d(n)|w) were (diagonal) Gaussian in w, they would each cancel with the corresponding (diagonal) Gaussian term, leaving only logZq. Intuitively, the residual terms in Eq. (6) can be used to quantify the quality of the Gaussian approximation.

heading: 3.2 MEMORY UPDATE
text: 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 VB— we 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(dtk |w) ≈ rtk(w;dtk) resulting in qθtk (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̃−1tk qθtk (w) p(M|w) || q̃θtk (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 Dtk andMtk−1 :
q̃θtk (w) = argmin qθ
KL [ qθ(w) || Z̃−1tk qθtk−1 (w) p(Dtk |w) p(Mtk−1 |w) ] . (7)
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(θtk ;Dtk ∪Mtk−1) in the form of Eq. (6):
L(θtk ; Dtk ∪ Mtk−1) = logZqtk + ∑
dtk ∈Dtk∪Mtk−1
Eq̃θtk (w) [ log p(dtk |w)− log rtk(w;dtk) ] ,
where dtk are the samples in Dtk ∪Mtk−1 and where rtk(w;dtk) are the Gaussian approximation of the corresponding likelihood terms. Note that rtk does not only depend on dtk , however, we omit the dependence on the remaining data for notational convenience.
If the likelihood p(dtk |w) is close to the Gaussian rtk(w;dtk) in expectation w.r.t. the approximate posterior qθtk (w), then its contribution to L(θtk ;Dtk ∪Mtk−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
Stk(M; Dtk ∪ Mtk−1) = ∑
dtk ∈M
Eq̃θtk (w) [ log p(dtk |w)− log rtk(w;dtk) ] , (8)
and the corresponding memory updateMtk = argmaxM Stk(M; Dtk ∪ Mtk−1). Note that since all residual terms are computed independently, the update results in selecting the top M terms.

heading: 3.3 GAUSSIAN UPDATE
text: The Gaussian update follows from the memory update presented in the previous section: once the memoryMtk has been selected, we update the Gaussian distribution with the approximations corresponding to the rest of the data Dtk ∪Mtk−1 \Mtk . We can update qθtk (w) in two equivalent ways:
qθtk (w) = qθtk−1 (w) ∏
dtk 6∈Mtk
rtk (w;dtk ), (9a) qθtk (w) = q̃θtk (w) / ∏
dtk ∈Mtk
rtk (w;dtk ). (9b)
Note again that the natural parameters of rtk(w;dtk) 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 |Dtk | ≤ |Mtk |, and Eq. (9b) if |Dtk | > |Mtk |. 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 qtk(w) would be to use VB on the data Dtk ∪ Mtk−1 \Mtk to compute the update qθtk (w) ≈ qθtk−1 (w) p(Dtk ∪ Mtk−1 \Mtk |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 Dtk , 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 Ntk +M instead of Ntk data points to find a better optimum of the ELBO.

heading: 4 VARIATIONAL BAYES WITH MODEL ADAPTATION
text: 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 Dtk at potentially non-equidistant time steps tk; (ii) data within Dtk is assumed i.i.d., however, not between different datasets Dtk and Dtk+1 . In both approaches, we realise adaptation by an additional forgetting step before observing the new data Dtk+1 . We denote the distribution, which results from applying the forgetting step to the posterior approximation qθtk (w) p(Mtk |w) by ptk+1(w).

heading: 4.1 ADAPTATION WITH BAYESIAN FORGETTING
text: 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(wtk+1 |Dt1:tk) directly. Here, we use a continuous-time version of this forgetting operation that can be formulated as
p(w|Dt1:tK ) ∝ p0(w) K∏ k=1 p(Dtk |w)(1− ) tK−tk τ , (10)
where τ is a time-constant corresponding to the average of the time-lags ∆tk+1 = tk+1 − tk. The distribution defined in Eq. (10) can be formulated recursively (cf. App. F) as
p(w|Dt1:tk+1) ∝ p0(w)1−(1− ) ∆tk+1/τ p(w|Dt1:tk)(1− ) ∆tk+1/τ
p(Dtk+1 |w). (11) 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θtk (w) p(Mtk |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 ∆tk(m), the forgetting operation for this new posterior approximation then results in
ptk+1(w) ∝ p0(w)1−(1− ) ∆tk+1/τ × [ qθtk (w) ∏ m∈Mtk p(m|w)(1− ) ∆tk(m)/τ ](1− )∆tk+1/τ
= [ p0(w) 1−(1− )∆tk+1/τ qθtk (w) (1− )∆tk+1/τ ] × ∏ m∈Mtk p(m|w)(1− ) ∆tk+1(m)/τ , (12)
where ∆tk+1(m) = ∆tk(m) + ∆tk+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 ∆tk+1(m). For the Gaussian term qθtk (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
σ−2tk+1 = ( 1− (1− )∆tk+1/τ ) σ−20 + (1− ) ∆tk+1/τσ−2tk ,
σ−2tk+1µtk+1 = ( 1− (1− )∆tk+1/τ ) σ−20 µ0 + (1− ) ∆tk+1/τσ−2tk µtk .
For ∆tk+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 .

heading: 4.2 ADAPTATION WITH DIFFUSION PROCESSES
text: 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 ptk+1,tk(w
′|w) resulting in a prediction distribution
ptk+1(w ′) = ∫ ptk+1,tk(w ′|w) p(w|Dt1:tk) dw, (13)
where we consider Gaussian transitions ptk+1,tk(w ′|w). However, this operation is generally not tractable for our posterior approximation qθtk (w) p(Mtk |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 ptk+1(w ′) = ∫ ptk+1,tk(w ′|w) qθtk (w) dw.
As mentioned in Sec. 4.1, BF yields the prior distribution for ∆tk+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 dwt = θ · (µ0 −wt) dt+ σ0 √ 2θ dWt, 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 ptk+1(w) = N ( µtk+1 , σ 2 tk+1 ) is defined by the parameters
µtk+1 = ( 1− e−a ∆tk+1 τ ) µ0 + e −a ∆tk+1 τ µtk ,
σ2tk+1 = ( 1− e−2a ∆tk+1 τ ) σ20 + e −2a ∆tk+1 τ σ2tk .
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.

heading: 5 RELATED WORK
text: 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) 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. (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.

heading: 6 EXPERIMENTS
text: 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−1test ∑ n logEq̃θtk (w) [ p(d (n) test|w) ] , where d(n)test is a sample from a heldout test dataset; (ii) the avg. one-
step-ahead LML N−1tk+1 ∑ n logEq̃θtk (w) [ p(d (n) tk+1 |w) ] , where d(n)tk+1 is data observed at time-step tk+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.

heading: 6.1 ONLINE LEARNING
text: 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 p0(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 t0 and 15 samples at all consecutive time-steps. In Fig. 3, we visualise the selected memory and the corresponding scores for time-steps t0, t25, and t50, 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 Ntk 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 t0. This is because these methods perform the initial Gaussian update by optimising the ELBO with Nt0 −M samples at t0; by contrast, GRS uses a Gaussian update that first optimises the ELBO with Nt0 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.

heading: 6.2 ADAPTATION
text: 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 Dtk = {(xtk ,ytk)}n, xtk ∈ R2, ytk ∈ R is sampled from xtk ∼ Uniform(−3, 3), ytk ∼ Bernoulli(σ(wtkxtk)). The true model has two time-dependent parameters w0tk = 10 sin(α · tk), w 1 tk
= 10 cos(α · tk), where α = 5 deg /sec and where we observe data at tk ∈ [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.

heading: 7 CONCLUSION
text: 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.

heading: 8 APPENDIX

heading: A FURTHER EXPERIMENTAL RESULTS
text: 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 Dtk ∪Mtk−1\Mtk (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.
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.

heading: B EXPERIMENT SETUP
text: 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 (Dtk ∪ Mtk−1) \Mtk for the Gaussian update.
• random (VCL): Uses random selection (Sec.2.2) for the memory update and Eq. (2) with (Dtk ∪ Mtk−1) \Mtk 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 (Dtk ∪ Mtk−1) and subsequently using Eqs. (5a), (5b) for removing the local contributions ofMtk (cf. Sec. 3.3).
Similarly, the following list summarises the adaptation methods used in Sec. 6.2:
• Wiener process: Posterior approximation consists of qθtk (w) only. Transition p(wtk+1 |wtk) 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θtk (w) only. Transition p(wtk+1 |wtk) is given by the Ornstein-Uhlenbeck process (cf. Sec. 4.2). No memory used.
• Bayesian forgetting: Posterior approximation consists of qθtk (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θtk (w) andMt. 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.

heading: C FACTORISATION PROPERTY OF THE GAUSSIAN VARIATIONAL APPROXIMATION
text: 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
L(µ∗,Σ∗) = −1 2
( log |Σ0| − log |Σ∗| − d+ (µ∗ − µ0)TΣ−10 (µ∗ − µ0) + Tr(Σ∗Σ −1 0 ) )
+ N∑ n=1 Ew∼qθ∗ (w) [ log p(d(n)|w) ] .
At a local optimum, we have ∂L(µ ∗,Σ∗)
∂µ∗ = 0, which yields
N∑ n=1 ∂ ∂µ∗ Ew∼qθ∗ (w) [ log p(d(n)|w) ] = Σ−10 (µ ∗ − µ0).
Hence, we obtain
µ∗ = µ0 + Σ0 N∑ n=1 ∂ ∂µ∗ Ew∼qθ∗ (w) [ log p(d(n)|w) ] . (14)
Similarly, we have ∂L(µ ∗,Σ∗)
∂Σ∗ = 0, which yields
N∑ n=1 ∂ ∂Σ∗ Ew∼qθ∗ (w) [ log p(d(n)|w) ] = −1 2 (Σ∗)−1 + 1 2 Σ−10 . (15)
Hence, we obtain
Σ∗ = ( Σ−10 − 2 N∑ n=1 ∂ ∂Σ∗ Ew∼qθ∗ (w) [ log p(d(n)|w) ])−1 . (16)
Next, we calculate the natural parameters from Eqs. (14), (16):
Λ∗ = Λ0 + N∑ n=1 −2 ∂ ∂Σ∗ Ew∼qθ∗ (w) [ log p(d(n)|w) ]
= Λ0 + N∑ n=1 Λ(n).
η∗ = Λ∗µ∗ = (
Λ0 + N∑ n=1 Λ(n) ) µ∗
= Λ0 ( µ0 + Σ0 N∑ n=1 ∂ ∂µ∗ Ew∼qθ∗ (w) [ log p(d(n)|w) ]) + N∑ n=1 Λ(n)µ∗
= Λ0µ0 + N∑ n=1 ( ∂ ∂µ∗ Ew∼qθ∗ (w) [ log p(d(n)|w) ] + Λ(n)µ∗ )
= η0 + N∑ n=1 η(n).
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:
Λb = (Λ ∗ − Λ0)− N∑ n=1 Λ(n), ηb = (η ∗ − η0)− N∑ n=1 η(n).
We use this to reduce the bias for the individual terms:
Λ̂(n) = Λ(n) − 1 N Λb, η̂ (n) = η(n) − 1 N ηb.

heading: D ELBO IN RESIDUALS FORM
text: 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−1q p(w) ∏N n=1 r
(n)(w) (cf. Sec.3.1). We can then write the ELBO as
L(µ,Σ;D) = Eqθ(w) [ N∑ n=1 log p(d(n)|w) + log p(w)− log qθ(w) ]
= Eqθ(w) [ N∑ n=1 ( log p(d(n)|w) + log r(n)(w)− log r(n)(w) ) + log p(w)− log qθ(w) ]
= Eqθ(w) [ N∑ n=1 ( log p(d(n)|w)− log r(n)(w) ) + log p(w) ∏N n=1 r (n)(w) qθ(w) ]
= Eqθ(w) [ N∑ n=1 ( log p(d(n)|w)− log r(n)(w) ) + log Zq · qθ(w) qθ(w) ] = logZq +
∑ n Eqθ(w) [ log p(d(n)|w)− log r(n)(w) ] .

heading: E MEMORY UPDATE SCORE FUNCTION
text: In Eq. 8, the expectation involving Gaussian terms can be calculated analytically: Eq̃θtk (w) [ log r (m) tk (w) ] = ∫ q̃θtk (w) ( η(n)w − 1 2 Λ(n)w2 ) 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 Eq̃θtk (w)
[ log p(d (n) tk |w) ] using Monte-Carlo.

heading: F BAYESIAN FORGETTING - RECURSIVE FORMULATION
text: 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|Dt1:tK+1) ∝ p0(w) · K+1∏ k=1 p(Dtk |w)(1− ) tK+1−tk τ
= p0(w) · K∏ k=1 p(Dtk |w)(1− ) tK+1−tk τ · p(DtK+1 |w).
The first two terms can be rewritten as
p0(w) · K∏ k=1 p(Dtk |w)(1− ) tK+1−tk τ = p0(w) · K∏ k=1 p(Dtk |w)(1− ) tK+1−tK+tK−tk τ
= p0(w) · K∏ k=1 p(Dtk |w)(1− ) tK−tk τ ·(1− ) tK+1−tK τ
= p0(w) · ( K∏ k=1 p(Dtk |w)(1− ) tK−tk τ )(1− ) tK+1−tKτ
∝ p0(w) · ( p(w|Dt1:tK ) p0(w) )(1− ) tK+1−tKτ = p0(w) 1−(1− ) tK+1−tK τ · p(w|Dt1:tK )(1− ) tK+1−tK τ
Hence, we have shown that the posterior can be expressed recursively as
p(w|Dt1:tk+1) ∝ p0(w)1−(1− ) ∆tk+1/τ p(w|Dt1:tk)(1− ) ∆tk+1/τ
p(Dtk+1 |w). The parameters of the Gaussian part qθtk (w) of the posterior approximation (after applying the forgetting operation) can be calculated easily from the above equation.
Natural parameters:
Λtk+1 = Λ0 · (1− (1− ) tk+1−tk τ ) + Λtk · ((1− ) tk+1−tk τ )
ηtk+1 = η0 · (1− (1− ) tk+1−tk τ ) + ηtk · ((1− ) tk+1−tk τ )
Covariance parameter: Σtk+1 = ( Λ0 · (1− (1− ) tk+1−tk τ ) + Λtk · ((1− ) tk+1−tk τ ) )−1
= ( Σ−10 · (1− (1− ) tk+1−tk τ ) + Σ−1tk · ((1− ) tk+1−tk τ ) )−1
Location parameter:
µtk+1 = Σtk+1 · ηtk+1
= Σtk+1 ( η0 · (1− (1− ) tk+1−tk τ ) + ηtk · ((1− ) tk+1−tk τ ) ) = Σtk+1 ( Σ−10 µ0 · (1− (1− ) tk+1−tk τ ) + Σ−1tk µtk · ((1− ) tk+1−tk τ )
)

heading: G PSEUDO-ALGORITHM
text: 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 Itk sequential iterations with Ntk +M data samples and Ktrain Monte Carlo samples. The latter can both be processed in parallel on parallel hardware. The estimation of the Gaussian factors requires Ntk +M sequential iterations with Kterm parallel Monte Carlo samples. The dominating computation for calculating the scores is the evaluation of the likelihood for Ntk +M data samples and Ktrain 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 EstimateGaussianFactors corresponds to Eqs. (5a), (5b) (cf. also App. C). The function ApplyForgetting corresponds to Eq. (12). Note that ptk(w) includes the adapted likelihood of the memory and all subsequent functions involving the memory use this adapted likelihood.
Inputs: p0(w), qθt0 (w), τ , K for k in [0...K] do tk = GetTimeStamp(k) ∆tk = tk − tk−1 Dtk = GetData(tk) if k == 0 then ptk(w) = p0(w)
else ptk(w) = ApplyForgetting ( p0(w), qθtk−1 ,Mtk−1 ,∆tk ) {Sec. 4.1} end if q̃θtk (w) = argminqθ KL [ qθ(w) || Z̃−1tk ptk(w) p(Dtk |w) ] {Sec. 3.2, Sec. 4.1}
{rtk(w;dtk)}dtk∈Dtk∪Mtk−1 = EstimateGaussianFactors ( q̃θtk (w),Dtk ,Mtk−1 ) {Sec. 3.1} Mtk = argmaxM Stk(M; Dtk ∪ Mtk−1) {Sec. 3.2} if |Dtk | ≤ |Mtk | then qθtk (w) = ptk(w) ∏ dtk 6∈Mtk
rtk(w;dtk) {Sec. 3.3} else qθtk (w) = q̃θtk (w) / ∏ dtk ∈Mtk
rtk(w;dtk) {Sec. 3.3} end if
end for

heading: H PRIOR PARAMETERS
text: Here we develop a heuristic to choose the initial prior p0(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 = {Wl}l ∪ {bl}l, and let x0 denote the input data. Let us further denote the pre-activation (before non-linearity) of layer l and unit i as follows.
xil = ∑ j Wi,jl fl−1(x j l−1) + b i l.
The constraints are then given as follows. Ew∼p(w) [ Ex0∼p(D) [ xil ]] = 0, Ew∼p(w) [ Ex0∼p(D) [ (xil) 2 ]] = 1.
The first constraint can be easily fulfilled by setting the prior mean to zero for all parameters.
µi,jl = 0.
This follows immediately from Wl ⊥ fl−1(xl−1) and the expectation of products of independent random variables. The second moment can then be calculated as follows.
Ew∼p(w),x0∼p(D) [ (xil) 2 ] = Ew∼p(w) [ Ex0∼p(D) [ (xil) 2 ]]
= Ew∼p(w) [ Varx0∼p(D) [ xil ] + 0 ]
= Ew∼p(w) [ Varx0∼p(D) [Nl−1∑ j Wi,jl fl−1(x j l−1) + b i l ]]
= Ew∼p(w) [ Nl−1∑ j ( Wi,jl )2 ·Varx0∼p(D) [ fl−1(x j l−1) ]]
= Nl−1∑ j Ew∼p(w) [( Wi,jl )2 ·Varx0∼p(D) [ fl−1(x j l−1) ]]
= Nl−1∑ j Ew∼p(w) [( Wi,jl )2] · Ew∼p(w) [ Varx0∼p(D) [ fl−1(x j l−1) ]]
=: Nl−1∑ j Ew∼p(w) [( Wi,jl )2] · cfl−1 = cfl−1 · Nl−1∑ j Varw∼p(w) [ Wi,jl
] = Nl−1 · cfl−1 ·Varw∼p(w) [ Wi,jl ] .
Here we introduced cfl−1 to denote a correction factor for the non-linearity fl−1. In case of the linear function, we will have cfl−1 = 1. For arbitrary non-linearities, we can estimate this factor numerically, assuming that the pre-activations are distributed according to N (0, 1).
cfl−1 = Varxjl−1∼N (0,1)
[ fl−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.(
σi,jl )2 = Varw∼p(w) [ Wi,jl ] =
1
Nl−1cfl−1 (17)

heading: I POSTERIOR INITIALISATION
text: 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; θµ, e2θ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.
Eθ∼p(θ),w∼q(w|θ),x0∼p(D) [ (xil) 2 ] = Eθ∼p(θ) [ Ew∼q(w|θ) [ Ex0∼p(D) [( xil )2]]]
= Eθ∼p(θ) [ Ew∼q(w|θ) [ Varx0∼p(D) [Nl−1∑ j Wi,jl · fl−1(x j l−1) ] + 0 ]]
= Eθ∼p(θ) [ Ew∼q(w|θ) [Nl−1∑ j (Wi,jl ) 2 ] Varx0∼p(D) [ xjl−1 ]]
= Eθ∼p(θ) [ Ew∼q(w|θ) [Nl−1∑ j (Wi,jl ) 2 ]] · Eθ∼p(θ) [ Ew∼q(w|θ) [ Varx0∼p(D) [ fl−1(x j l−1) ]]]
= Nl−1∑ j Eθ∼p(θ) [ Ew∼q(w|θ) [ (Wi,jl ) 2 ]] · Eθ∼p(θ) [ Ew∼q(w|θ) [ Varx0∼p(D) [ fl−1(x j l−1) ]]]
=: Nl−1∑ j Eθ∼p(θ) [ θ2µ + e 2·θlog σ ] · cfl−1
= Nl−1 · cfl−1 · Eθ∼p(θ) [ θ2µ + e 2·θlog σ ]
= Nl−1 · cfl−1 · ( Eθ∼p(θ) [ θ2µ ] + Eθ∼p(θ) [ e2·θlog σ ]) = Nl−1 · cfl−1 · ( Eθ∼p(θ) [ θµ ]2 + Varθ∼p(θ) [ θµ ] + e2E[θlog σ]+2Var[θlog σ ] )
= Nl−1 · cfl−1 · ( µ2θµ + σ 2 θµ + e 2µθlog σ+2σ 2 θlog σ ) = Nl−1 · cfl−1 · ( σ2θµ + e 2µθlog σ+2σ 2 θlog σ ) .
Hence, the second constraint is as follows.
1
Nl−1 · cfl−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θµ = β
Nl−1 · cfl−1 ,
and e2µθlog σ+2α 2 = 1− β
Nl−1 · cfl−1 .
Solving the last equation for µθlog σ , the result is as follows.
µθlog σ = 1
2 log
(1− β) · e−2α2
Nl−1 · cfl−1 We choose α = 0.001 and β = 0.999 in all experiments.
A note on the relation to initialisation methods for deterministic neural networks. Our result is similar to the initialisation methods from Glorot & Yoshua Bengio (2010) and He et al. (2015). The difference is in the correction factor cfl−1 . Whereas Glorot & Yoshua Bengio (2010) considers linear functions (or tanh in the linear regime), both methods base their derivation on the assumption that every data sample x0 is processed by a different, random neural network with independent weights, drawn from the initialisation distribution. The assumption is made explicit in (He et al., 2015) by the use of the variance of products of independent variables rule. We note that this assumption is false for both the initialisation of deterministic neural networks, as well as the graphical model assumption in Bayesian neural networks. Consequently, (He et al., 2015) obtains different correction factors (in their case for relu and leaky relu), taking into account the mean after the forward-pass through the non-linearity.

Review summaries:
This paper proposes a method for training Bayesian neural nets on a stream of non-stationary data. The authors propose a new way to approximate the variational distribution in the case where a memory (of past data) is used. The authors also use two forgetting mechanisms that allow the model to adapt to the changes in the data distribution.