Variational Inference over Nonstationary Data Streams for Exponential Family Models ^†

Abstract

In many modern data analysis problems, the available data is not static but, instead, comes in a streaming fashion. Performing Bayesian inference on a data stream is challenging for several reasons. First, it requires continuous model updating and the ability to handle a posterior distribution conditioned on an unbounded data set. Secondly, the underlying data distribution may drift from one time step to another, and the classic i.i.d. (independent and identically distributed), or data exchangeability assumption does not hold anymore. In this paper, we present an approximate Bayesian inference approach using variational methods that addresses these issues for conjugate exponential family models with latent variables. Our proposal makes use of a novel scheme based on hierarchical priors to explicitly model temporal changes of the model parameters. We show how this approach induces an exponential forgetting mechanism with adaptive forgetting rates. The method is able to capture the smoothness of the concept drift, ranging from no drift to abrupt drift. The proposed variational inference scheme maintains the computational efficiency of variational methods over conjugate models, which is critical in streaming settings. The approach is validated on four different domains (energy, finance, geolocation, and text) using four real-world data sets.

1. Introduction

One core problem in Bayesian statistics is the computation of posterior probability over the parameters (and the latent variables) of a model given a data set. In most relevant cases, the computation of this posterior probability is not feasible and requires the use of approximate inference algorithms [1,2]. Moreover, in many real-life settings, the data set is not static but arrives sequentially, in a streaming fashion. Computing the Bayesian posterior probability over the parameters (and the latent variables) of a model, in this case, is even more challenging. First, it requires the ability to handle a posterior distribution conditioned on an unbounded data set. Secondly, the posterior probability has to be updated frequently, at every time step, imposing restrictions over the speed of computation of the posterior probability. Finally, the classical i.i.d. (independent and identically distributed, or data exchangeability) assumption does not hold because the underlying data distribution may drift/change from one time step to the next. This last issue is known in the machine learning literature as learning from data streams that exhibit concept drift [3] in the sense that the underlying generative process may contain both gradual and abrupt changes.

In this paper, we look at the common setting in many real-world problems, where data arrives as a sequence of (potentially large) batches. Specifically, each batch is associated with a new time step and the data points within a batch are assumed to be i.i.d. However, the underlying distribution generating the batches may change from one time step to the next. For illustration, in Section 6, we will apply a probabilistic topic model [4] to identify the main topics of research papers published over a sequence of years. New papers are submitted yearly; hence, every year, the probabilistic model must be updated to infer the topics of that year’s papers. One approach could be to discard all data from previous years, but that would result in a loss of relevant data that could potentially have been used to draw more accurate inferences about the current year’s topic distribution (research topics typically span multiple years). On the other hand, we also have to take into account that the underlying distribution of the papers may change from one year to the next following the development of the research field. As another example domain, we will model the financial profile of customers asking for a loan/mortgage from a financial institution [5]. In this example, we update the financial profile of the customers on a monthly basis while simultaneously keeping in mind that the occurrence of exogenous events (like an economic crisis) may strongly affect the customers and therefore induce a drift in the underlying data distribution. Other situations similar to the ones above arise in many different real-life settings.

Standard temporal models [6] have the potential to capture the underlying temporal dynamics of the data set. However, for the problems addressed in this paper, temporal models are not readily applicable. The reason for this is two-fold. First, the set of objects differs from one time step to the next. For example, a new set of papers is submitted to a conference every year. Thus, temporal modeling needs to be applied at the level of the parameters in the model, not at the objects themselves. Secondly, drift (both gradual and abrupt) in the underlying distributions is not easily modelled with a stationary transition model as would typically be the case in temporal models. Instead, we will employ an implicit transition model to model the temporal dynamics, thus sidestepping the problem of specifying an explicit stationary transition distribution.

The approach presented in this paper is applicable to domains that, at each time point, can be described by a conjugate exponential family model. These models are widely used in probabilistic modeling [7,8] and include popular models like latent Dirichlet allocation (LDA) models [4] to uncover the hidden topics in a text corpora, a mixture of (multivariate) Gaussian models to discover hidden clusters in data [9], and probabilistic principal component analysis for revealing a low-dimensional representation of the data [10]. See, for instance, [9,11,12] for detailed reviews and applications of these models in data mining and machine learning settings.

Exact Bayesian inference is intractable for the model-class proposed in this paper, partly due to the implicit transition model and partly due to the high dimensionality of the parameters/latent variables in the model. We therefore resort to approximate Bayesian inference methods, and owing to the potential real-time requirements imposed by the streaming domains, we focus on variational methods (see [13] for an introduction). Markov Chain Monte Carlo methods (MCMC) [1,2], although widely used in Bayesian statistics, are not as computationally efficient as variational methods, especially in the presence of large datasets or complex models [13]. Moreover, variational inference can be performed very efficiently in conjugate exponential family models by exploiting properties of this model family.

The rest of the paper is organized as follows. Section 2 discusses related works, and Section 3 introduces preliminaries relevant for the remainder of the paper. In Section 4, we introduce an implicit transition model that works as an off-the-self parameter transition scheme to model both gradual and abrupt changes in the data distribution. Section 5 positions the transition model in a hierarchical Bayesian setting to model the rate of change of the data stream as an unobserved mechanism. Section 5 also includes the derivation of an ad hoc variational inference algorithm [13], based on a novel lower-bound of the data log-likelihood function for the proposed model family (i.e., a conjugate exponential model with an implicit transition scheme). As a consequence, the proposed approximate inference scheme is deterministic and scales to large data streams. Section 6 empirically evaluates the appropriateness of our approach using both synthetic and real-life data (covering energy, finance, geolocation, and text data), showing promising results. We conclude in Section 7, where the main conclusions and future works are discussed.

2. Related Work

Making inferences from nonstationary data streams has been extensively studied in statistics. Time series models [14] is a classic example of models that could be assumed relevant for this problem, but as we are not dealing with temporal data (where the the same set of objects are observed at each time step), this model class is not applicable in our situation. Work by [15,16] is closer to ours, but these contributions focus on simpler model classes (like generalized linear models without latent variables), and due to this simplicity, they do not consider variational inference.

Methods for change point detection [17,18] are also related to our situation. However, where these methods consider abrupt changes in the data stream, our method is able to deal with both abrupt and smooth changes in the distribution. Additionally, the main focus of change point detection methods is the detection of points where the changes occur, while our goal is to have a statistical model that accurately models the data at each time step by continuously adapting the model parameters.

Learning from nonstationary data streams has also been extensively studied in the machine learning literature, especially in the context of classification and clustering [3,19,20,21]. One of the main techniques employed in this area has been exponential forgetting [22,23]. Here, each new observation is initially assigned a weight (or “importance”) equal to 1, followed by an exponential decrease of the weight at each subsequent time step. In this way, older observations are less relevant than newer ones when the model is learned, thereby accounting for potential drift in the data stream. The main problem with this approach is to quantify the so-called exponential forgetting rate, which is usually determined by the analyst on a trial-and-error basis. In the present work, we provide a sound approach based on Bayesian methods to automatically adjust the exponential forgetting rate and show the advantages of this strategy empirically.

Bayesian modeling of nonstationary data streams for general probabilistic inference has not been extensively studied so far. An online variational inference method, which exponentially forgets old data, was proposed by Honkela and Valpola [24], but similarly to [22,23], this approach suffers from the problem of setting the exponential forgetting rate. Another proposal, called population variational Bayes (PVB), was introduced by [25]. It builds directly on the stochastic variational inference (SVI) algorithm [26]. SVI assumes the existence of a fixed data set observed in a sequential manner and, in particular, that this data set has a known finite size. This is unrealistic when modeling data streams. PVB addresses the problem by using the frequentist notion of a population distribution, $\mathbf{F}$, which is assumed to generate the data stream by repeatedly sampling M data points at a time; here, M parameterizes the size of the population and helps control the variance of the population posterior. By artificially having high variance in the posterior (i.e., by choosing a small value for M), PVB is able to accommodate drift in the data set. Unfortunately, M must be specified by the user, and no clear rule exists for specifying it. Furthermore, McInerney et al., [25] shows that the optimal value for M may differ between data stream. The streaming variational Bayes (SVB) algorithm by [27] also tries to address Bayesian inference in data streams. SVB builds on a Bayesian recursive updating approach but does not provide a mechanism for dealing with concept drift. In Section 6, we will show that our proposed method, which does not rely on hyperparameters that are hard to tune, outperforms these closely related approaches on several real-world data sets.

The so-called power prior approach [28] has also been studied in the context of data aggregation for Bayesian modeling. Power priors provide a sound mechanism for Bayesian updating in light of new data and introduces partial forgetting of older observations. The approach enjoys nice theoretical properties [29] but depends on a hyperparameter (set by the analyst) to control the forgetting rate. Our work builds partially on the power prior model but extends this model with a Bayesian model of the forgetting mechanism, thereby dispensing of the analyst-specified hyperparameter.

A modeling approach based on time-series models for concept drift using implicit transition models was pursued by [30,31]. Unfortunately, the implicit transition model also depends on a hyperparameter determining the forgetting factor, which has to be manually set. Our work also partially builds on this approach and establishes novel connections between the implicit transition models [30,31], power priors [28], and exponential forgetting [22,23].

Many other contributions have proposed an ad hoc extension to specific statistical models in order to deal with nonstationary data streams [32,33], including dynamic extensions of latent Dirichlet allocation (LDA) models [34,35,36]. However, none are so far applicable to general conjugate exponential family models. Moreover, most of these contributions rely on complex and tailor-made inference mechanisms, instead of variational inference (which is well understood for conjugate exponential models).

3. Preliminaries

3.1. Conjugate Exponential Models

Let $\mathit{x}={\mathit{x}}_{1:N}$ denote a set of observed variables, $\mathit{z}={\mathit{z}}_{1:N}$ denote a set of latent variables, and $\mathbf{\beta }={\mathbf{\beta }}_{1:M}$ denote the parameters of the model. Let $\mathbf{\alpha }$ also denote a hyperparameter vector.

We assume that the joint distribution of our statistical model factorizes into a product of local terms and a global term:

$$p(\mathit{x},\mathit{z},\mathbf{\beta }|\mathbf{\alpha })=p\left(\mathbf{\beta }\right|\mathbf{\alpha })\prod _{n=1}^{N}p({\mathit{x}}_n,{\mathit{z}}_n|\mathbf{\beta }).$$

This kind of factorization is quite usual in many statistical models [26], for example, the case of the Bayesian mixture of Gaussian, where each latent variable is a categorical variable defining the component to which each data point belongs. We also consider the special case when there are no latent variables; this case would cover models like Bayesian linear regression. Figure 1a provides a graphical description of these model family in terms of probabilistic graphical models with plateau notation [11].

We further assume that the (conditional) distributions defining the statistical model belong to the conjugate exponential family [37]. This model family has been largely studied and covers a wide range of probability distributions such as multinomial, normal, gamma, Dirichlet, beta, etc. According to this assumption, the functional form of the conditional distribution of the local variables $({\mathit{x}}_n,{\mathit{z}}_n)$ given the parameters $\mathbf{\beta }$ has the well-known exponential family form [37]:

$$lnp({\mathit{x}}_n,{\mathit{z}}_n|\mathbf{\beta })=lnh({\mathit{x}}_n,{\mathit{z}}_n)+{\mathbf{\beta }}^{T}t({\mathit{x}}_n,{\mathit{z}}_n)-a_l\left(\mathbf{\beta }\right),$$

(1)

where the scalar functions $h(·)$ and $a_l(·)$ are the base measure and the log-normalizer, respectively, and the vector function $\mathit{t}(·)$ is the sufficient statistics vector. Moreover, the prior distribution $p\left(\mathbf{\beta }\right)$ also belongs to the exponential family and has the following structure:

$$lnp\left(\mathbf{\beta }\right)=lnh\left(\mathbf{\beta }\right)+{\mathbf{\alpha }}^{T}t\left(\mathbf{\beta }\right)-a_g\left(\mathbf{\alpha }\right)$$

(2)

where the sufficient statistics are $t\left(\mathbf{\beta }\right)=(\mathbf{\beta },-a_l\left(\mathbf{\beta }\right))$ and the hyperparameter $\mathbf{\alpha }$ has two components $\mathbf{\alpha }=({\mathbf{\alpha }}^{*},\xi )$, where the first component ${\mathbf{\alpha }}^{*}$ has the same dimension as $\mathbf{\beta }$ and encodes the prior belief about the distribution over $\mathbf{\beta }$ and the second component $\xi >0$ is a scalar and encodes the strength in our prior belief [38]. This second parameter is also known in the literature [39] as the equivalent sample size of the prior distribution ($ES{S}_{prior}$).

Our inference goal is to approximate the posterior distribution over the parameters and latent variables given the following observations:

$$p(\mathbf{\beta },\mathit{z}|\mathit{x})=\frac{p\left(\mathbf{\beta }\right){\prod }_{n}p({\mathit{x}}_n,{\mathit{z}}_n|\mathbf{\beta })}{{\int }_{\mathbf{\beta }}p\left(\mathbf{\beta }\right){\prod }_n{\int }_{{\mathit{z}}_n}p({\mathit{x}}_n,{\mathit{z}}_n|\mathbf{\beta })d{\mathit{z}}_{n}d\mathbf{\beta }}.$$

However, this posterior is intractable in general due to the integral in the denominator (i.e., the evidence integral), and we therefore have to resort to approximate inference algorithms such as variational inference. These (and other) inference algorithms benefit enormously when the complete conditional distribution of the parameters and the latent variables belongs to the exponential family [26]. Therefore, we make this extra assumption about the model class (Note however, that the presented approach also applies to the more general conjugate exponential family.), which states that the conditional distribution over $\mathbf{\beta }$ and $\mathit{z}$ given the rest of the variables has the same functional form as the priors:

$$\begin{array}{c}\hfill lnp\left(\mathbf{\beta }\right|\mathit{x},\mathit{z},\mathbf{\alpha })=lnh\left(\mathbf{\beta }\right)+{\eta }_g{(\mathit{x},\mathit{z},\mathbf{\alpha })}^{T}t\left(\mathbf{\beta }\right)-a_g(\mathit{x},\mathit{z},\mathbf{\alpha })\\ \hfill lnp\left({\mathit{z}}_n\right|{\mathit{x}}_n,\mathbf{\beta })=h\left({\mathit{z}}_n\right)+{\eta }_l{({\mathit{x}}_n,\mathbf{\beta })}^{T}t\left({\mathit{z}}_n\right)-a_l\left({\eta }_l({\mathit{x}}_n,\mathbf{\beta })\right),\end{array}$$

where the vector function ${\eta }_{·}(·)$ denotes the natural parameter vectors of the conditional probability distributions.

By Equations (1) and (2), the natural parameter vector of $p\left(\mathbf{\beta }\right|\mathit{x},\mathit{z},\mathbf{\alpha })$ can be expressed as

$${\eta }_g(\mathit{x},\mathit{z},\mathbf{\alpha })=\left({\mathbf{\alpha }}^{*}+\sum _{n=1}^{N}t({\mathit{x}}_n,{\mathit{z}}_n),\xi +N\right).$$

(3)

Therefore, computing the full posterior reduces to updating the natural parameters of the prior. Moreover, the equivalent sample size of the posterior ($ES{S}_{post}$) is equal to the equivalent sample size of the prior ($ES{S}_{prior}$) plus the size of the observations.

3.2. Variational Inference

Variational inference is a deterministic technique for finding a tractable posterior distribution, denoted by q, which approximates the true posterior, $p(\mathbf{\beta },\mathit{z}|\mathit{x})$, that is often intractable to compute. More specifically, by letting $\mathcal{Q}$ be a set of possible approximations of this posterior, variational inference solves the following optimization problem for any model in the conjugate exponential family:

$$\underset{q\left(\mathbf{\beta },\mathit{z}\right)\in \mathcal{Q}}{min}\mathit{KL}\left(q(\mathbf{\beta },\mathit{z})\left|\right|p(\mathbf{\beta },\mathit{z}|\mathit{x})\right),$$

(4)

where $\mathit{KL}$ denotes the Kullback–Leibler divergence between two probability distributions.

In the mean field variational approach, the approximation family $\mathcal{Q}$ is assumed to fully factorize. Following the notation of Hoffman et al. [26], we have that

$$q(\mathbf{\beta },\mathit{z}|\mathit{\lambda },\mathit{\varphi })=q\left(\mathbf{\beta }\right|\mathit{\lambda })\prod _{n=1}^{N}q\left({\mathit{z}}_n\right|{\mathit{\varphi }}_n).$$

Furthermore, each factor of the variational distribution is assumed to belong to the same family as the model’s complete conditionals:

$$\begin{array}{c}\hfill lnq\left(\mathbf{\beta }\right|\mathit{\lambda })=h\left(\mathbf{\beta }\right)+{\mathit{\lambda }}^{T}t\left(\mathbf{\beta }\right)-a_g\left(\mathit{\lambda }\right)\\ \hfill lnq\left({\mathit{z}}_n\right|{\mathit{\varphi }}_n)=h\left({\mathit{z}}_n\right)+{\mathit{\varphi }}_n^{T}t\left({\mathit{z}}_n\right)-a_l\left({\mathit{\varphi }}_n\right).\end{array}$$

It can be seen that $\mathit{\lambda }$ parameterizes the variational distribution of $\mathbf{\beta }$ while $\mathit{\varphi }$ plays the same role for the variational distribution of $\mathit{z}$.

To solve the minimization problem in Equation (4), the variational approach exploits the transformation:

$$lnp\left(\mathit{x}\right)=\mathcal{L}(\mathit{\lambda },\mathit{\varphi }|\mathit{x},\mathbf{\alpha })+\mathit{KL}\left(q(\mathbf{\beta },\mathit{z}|\mathit{\lambda },\mathit{\varphi }\left)\right|\left|p\right(\mathbf{\beta },\mathit{z}\left|\mathit{x}\right)\right),$$

(5)

where $\mathcal{L}(·|·)$ is a lower bound of $lnp\left(\mathit{x}\right)$ since $\mathit{KL}$ is nonnegative. This lower bound can be written as

$$\mathcal{L}(\mathit{\lambda },\mathit{\varphi }|\mathit{x},\mathbf{\alpha })={\mathbb{E}}_q[lnp\left(\mathit{x}\right|\mathit{z},\mathbf{\beta })]-{\mathbb{E}}_q\left[\mathit{KL}\left(q\left(\mathit{z}\right|\mathit{\varphi }\left)\right|\left|p\right(\mathit{z}\left|\mathbf{\beta }\right)\right)\right]-\mathit{KL}\left(q\left(\mathbf{\beta }\right|\mathit{\lambda }\left)\right|\left|p\right(\mathbf{\beta }\left|\mathbf{\alpha }\right)\right)].$$

(6)

We introduce $\mathit{x}$ and $\mathbf{\alpha }$ in $ℒ$’s notation to make explicit the function’s dependency on $\mathit{x}$, the data sample, and $\mathbf{\alpha }$, the natural parameters of the prior over $\mathbf{\beta }$. As $lnp\left(\mathit{x}\right)$ is constant, and minimizing the $\mathit{KL}$ term is equivalent to maximizing the lower bound. Equation (6) shows the trade-off involved in the lower-bound. The first term measures the model’s fit to the data and favors variational posterior mass concentrated around the maximum likelihood estimate. The second and third terms are regularization terms and favor variational posteriors close to their respective prior distributions.

This lower bound can be maximized, for example, by a coordinate ascent method that iteratively updates each individual variational distribution while holding the others fixed. As shown in [26], these iterative updating equations have the following closed-form solutions:

$$\begin{array}{ccc}\hfill \mathit{\lambda }& =& \mathbf{\alpha }+\sum _{n=1}^N{\mathbb{E}}_{{\mathit{\varphi }}_n}\left[(t({\mathit{x}}_n,{\mathit{z}}_n),1)\right],\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\mathit{\varphi }}_n& =& {\mathbb{E}}_{\mathit{\lambda }}\left[{\eta }_l({\mathit{x}}_n,\mathbf{\beta })\right],\hfill \end{array}$$

(8)

where ${\mathbb{E}}_{\mathit{\lambda }}[·]$ and ${\mathbb{E}}_{{\mathit{\varphi }}_n}[·]$ denote the expected value according to $q\left(\mathbf{\beta }\right|\mathit{\lambda })$ and $q\left({\mathit{z}}_n\right|{\mathit{\varphi }}_n)$, respectively. If the number of data points is large, alternative scalable methods can also be used [26,40].

3.3. Variational Inference over Data Streams

As commented in the previous section, we envision a situation where the data stream is defined by a sequence of data batches generated at discrete points in time $\{{\mathit{x}}_1,\dots ,{\mathit{x}}_t\}$ and where each batch is composed by a set of data samples, ${\mathit{x}}_t={\mathit{x}}_{t,i=1:N_t}$. As new batches arrive, we want to update the posterior distribution over the parameters of the model. This can be addressed by applying a Bayesian recursive approach:

$$p\left(\mathbf{\beta }\right|{\mathit{x}}_1,\dots ,{\mathit{x}}_t)\propto p\left(\mathbf{\beta }\right|{\mathit{x}}_1,\dots ,{\mathit{x}}_{t-1})\prod _{i=1}^{N_t}\int p({\mathit{x}}_{t,i},{\mathit{z}}_{t,i}|\mathbf{\beta })d{\mathit{z}}_{t,i}.$$

Hence, updating the posterior at time t reduces to a problem of computing a posterior over $\mathbf{\beta }$ conditional on the data ${\mathit{x}}_t$ given that $\mathbf{\beta }$ gets a prior equal to $p\left(\mathbf{\beta }\right|{\mathit{x}}_1,\dots ,{\mathit{x}}_{t-1})$, i.e., the posterior in the previous time step.

When the above posterior is intractable to compute, we can use the streaming variational Bayes (SVB) algorithm [27]. This algorithm translates the above recursive updating approach to the variational settings described in the previous sections. Firstly, it approximates the posterior in the previous time step with a variational approximation, $p\left(\mathbf{\beta }\right|{\mathit{x}}_1,\dots ,{\mathit{x}}_{t-1})\approx q\left(\mathbf{\beta }\right|{\mathit{\lambda }}_{t-1})$, so that the new posterior is expressed as

$$p\left(\mathbf{\beta }\right|{\mathit{x}}_1,\dots ,{\mathit{x}}_t)\propto q\left(\mathbf{\beta }\right|{\mathit{\lambda }}_{t-1})\prod _{i=1}^{N_t}\int p({\mathit{x}}_{t,i},{\mathit{z}}_{t,i}|\mathbf{\beta })d{\mathit{z}}_{t,i}.$$

This posterior is still intractable due to integration. We therefore use variational inference to compute a new approximation $q\left(\mathbf{\beta }\right|{\mathit{\lambda }}_t)$ to the posterior at time t. The variational parameters are given as the solution to the optimization problem:

$$({\mathit{\lambda }}_t,{\mathit{\varphi }}_t)=arg\underset{{\mathit{\lambda }}_t,{\mathit{\varphi }}_t}{min}\mathcal{L}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1}),$$

which, similarly to the static case (see Equation (6)), can be solved by a coordinate ascent method that iteratively updates each individual variational distribution while holding the others fixed. The end result is defined by the following closed-form solutions:

$$\begin{array}{ccc}\hfill {\mathit{\lambda }}_t& =& {\mathit{\lambda }}_{t-1}+\sum _{i=1}^{N_t}{\mathbb{E}}_{{\mathit{\varphi }}_{t,i}}\left[(t({\mathit{x}}_{t,i},{\mathit{z}}_{t,i}),1)\right],\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\mathit{\varphi }}_{t,i}& =& {\mathbb{E}}_{\mathit{\lambda }}\left[{\eta }_l({\mathit{x}}_{t,i},\mathbf{\beta })\right],\hfill \end{array}$$

(10)

where ${\mathbb{E}}_{\mathit{\lambda }}[·]$ and ${\mathbb{E}}_{{\mathit{\varphi }}_{t,i}}[·]$ denote the expected value according to $q\left(\mathbf{\beta }\right|{\mathit{\lambda }}_t)$ and $q\left({\mathit{z}}_{i,t}\right|{\mathit{\varphi }}_{i,t})$, respectively.

3.4. Exponential Forgetting in Variational Inference

Exponential forgetting is a classic technique [19] that allows to gradually forget past data and to put more focus on more recent data samples when learning from a nonstationary data stream. This idea is usually implemented by exponentially down-weighting the loss function term associated with each data sample so that data samples closer in time have more impact on the model that older data samples.

In probabilistic terms, exponential forgetting is achieved by using a log-likelihood function of the following form:

$$lnp({\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_t|\mathbf{\beta })=\sum _{i=1}^t{\rho }^{t-i}lnp\left({\mathit{x}}_i\right|\mathbf{\beta })+cte,$$

where $\rho \in [0,1]$ is the exponential decay weight and $cte$ is a constant term. By using a small $\rho$ value, we aggressively forget old data samples. However, $\rho$ values close to 1 tend to mildly forget previous data samples.

Similarly, in Bayesian learning settings [24], we can use this scheme to compute the posterior:

$$\begin{array}{ccc}\hfill p\left(\mathbf{\beta }\right|{\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_t,\rho )& \propto & p({\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_t|\mathbf{\beta },\rho )p\left(\mathbf{\beta }\right)\hfill \\ & =& p\left({\mathit{x}}_t\right|\mathbf{\beta })p{\left({\mathit{x}}_{t-1}\right|\mathbf{\beta })}^{\rho }\cdots p{\left({\mathit{x}}_1\right|\mathbf{\beta })}^{{\rho }^{t-1}}p\left(\mathbf{\beta }\right).\hfill \end{array}$$

This scheme also applies to variational learning by considering this exponential down-weighted log-likelihood instead of the standard data log-likelihood, as used by [24]. Then, the lower bound function has the following form:

$${\mathcal{L}}_{\rho }(\mathit{\lambda },\mathit{\varphi }|\mathit{x},{\mathbf{\alpha }}_u)={\mathbb{E}}_q\left[\sum _{i=1}^t{\rho }^{t-i}lnp\left({\mathit{x}}_i\right|{\mathit{z}}_i,\mathbf{\beta })\right]-\mathit{KL}\left(q(\mathbf{\beta },\mathit{z}|\mathit{\lambda },\mathit{\varphi })\left|\right|p(\mathbf{\beta },\mathit{z}|{\mathbf{\alpha }}_u)\right),$$

(11)

where ${\mathbf{\alpha }}_u$ denotes the natural parameters of a non-informative prior, $p(\mathbf{\beta },\mathit{z}|{\mathbf{\alpha }}_u)$.

The updating equation of the coordinate gradient ascent algorithm described in Equation (7) can now be expressed as follows [41]:

$$\mathit{\lambda }={\mathbf{\alpha }}_u+\sum _{i=1}^t{\rho }^{t-i}{\mathbb{E}}_{{\mathit{\varphi }}_i}\left[(t({\mathit{x}}_i,{\mathit{z}}_i),1)\right].$$

(12)

The main point here is to highlight how, at the convergence point, the variational solution $\mathit{\lambda }$ exponentially down-weights the contribution (i.e., the expected sufficient statistics) of old data samples.

Exponential forgetting also addresses the problem of Bayesian learning over unbounded data streams. According to Equation (7), one of the components of the $\mathit{\lambda }$ parameter corresponds to the equivalent sample size of the variational posterior ($ES{S}_{post}$). In this case, this value can be computed as

$$ES{S}_{post,t}=ES{S}_{prior}+\sum _{i=1}^t{\rho }^{(i-1)}.$$

If $\rho <1$, then $ES{S}_{post,t}$ converges to a finite number:

$$\underset{t\to \infty }{lim}ES{S}_{post,t}=ES{S}_{prior}+\frac{1}{1-\rho },$$

(13)

avoiding the problem of having a degenerated Bayesian posterior distribution in the presence of an unbounded data stream. As noted in [30,42], this schema approximates a posterior distribution conditioned on the last $\frac{1}{1-\rho }$ data samples of the stream.

Stochastic varitional inference (SVI) [26] is a widely used variational learning algorithm for dealing with large data sets. Population variational Bayes [25] is a simple modification of SVI used when the total size of the data set is unknown. When these algorithms are applied in data streaming settings, they use a constant learning rate $\nu$ (It is usually set to small values like 0.1 or 0.01.), and the sequential updating equation of the global variational parameters $\mathit{\lambda }$ can be written as

$${\mathit{\lambda }}_t=(1-\nu ){\mathit{\lambda }}_{t-1}+\nu ({\mathbf{\alpha }}_u+S{\mathbb{E}}_{{\mathit{\varphi }}_i}\left[(t({\mathit{x}}_t,{\mathit{z}}_t),1)\right]),$$

(14)

where S is equal the total size of the data set. This size is equal to N in the case of SVI or to the size of the population M in the case of PVB. By expanding this equation, we find that

$${\mathit{\lambda }}_t=(1-{(1-\nu )}^t){\mathbf{\alpha }}_u+N\nu \sum _{i=1}^t{(1-\nu )}^{t-i}{\mathbb{E}}_{{\mathit{\varphi }}_i}\left[(t({\mathit{x}}_i,{\mathit{z}}_i),1)\right].$$

(15)

The above equation highlights that SVI and PVB also exponentially down-weight old data samples, with a forgetting rate of $\rho =1-\nu$ (compare the above equation with Equation (12)). Therefore, this is one of the mechanisms that these two methods use to adapt to drifts in the nonstationary data stream.

In the case of the PVB method, the parameter M helps to adapt to drifts in the data set through the effect it has on computing ${\mathit{\varphi }}_i$, as discussed by [25]. However, when the model does not contain local random variables, the variational parameters ${\mathit{\varphi }}_i$ do not exist and, then, the size of the population does not play any role in adapting to drifts in the data stream.

4. Implicit Transition Models

In order to extend the model in Figure 1a to data streams, we may introduce a transition model $p\left({\mathbf{\beta }}_t\right|{\mathbf{\beta }}_{t-1})$ to explicitly model the evolution of the parameters over time, enabling estimation of the predictive density at time t:

$$p\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})=\int p\left({\mathbf{\beta }}_t\right|{\mathbf{\beta }}_{t-1})p\left({\mathbf{\beta }}_{t-1}\right|{\mathit{x}}_{1:t-1})d{\mathbf{\beta }}_{t-1}.$$

(16)

However, this approach introduces two problems. First, in nonstationary domains, we may not have a single transition model or the transition model may be unknown. Secondly, if we seek to position the model within the conjugate exponential family in order to be able to compute the gradients of $\mathcal{L}$ in closed form, we need to ensure that the distribution family for ${\mathbf{\beta }}_t$ is its own conjugate distribution, thereby severely limiting the model’s expressive power (e.g., we cannot assign a Dirichlet distribution to ${\mathbf{\beta }}_t$).

Rather than explicitly modeling the evolution of the ${\mathbf{\beta }}_t$ parameters as in Equation (16), we instead follow a similar approach to Kárny [31] and Ozkan et al. [30] by defining the time evolution model implicitly by constraining the maximum $\mathit{KL}$ divergence over consecutive parameter distributions. Specifically, by defining

$$p_{\delta }\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})=\int \delta ({\mathbf{\beta }}_t-{\mathbf{\beta }}_{t-1})p\left({\mathbf{\beta }}_{t-1}\right|{\mathit{x}}_{1:t-1})d{\mathbf{\beta }}_{t-1}$$

(17)

one can restrict the space of possible distributions $p\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})$, supported by an unknown transition model, by the constraint

$$\mathit{KL}\left(p\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})\left|\right|p_{\delta }\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})\right)\le \kappa .$$

(18)

We then propose to approximate $p\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})$ by the distribution $\widehat{p}\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})$ having minimum Kullback–Leibler divergence w.r.t. (with respect to) a prior density $p_u\left({\mathbf{\beta }}_t\right)$. This approach ensures that we will not underestimate the uncertainty in the parameter distribution. The following result shows how this constraining optimization problem has an amenable solution.

Theorem 1

(Implicit Transition Models).The density $\widehat{p}\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1},{\rho }_t)$ which has minimum Kullback–Leibler divergence w.r.t. a prior density $p_u\left({\mathbf{\beta }}_t\right)$,

$$\widehat{p}\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1},{\rho }_t)=arg\underset{q}{min}\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right)\left|\right|p_u\left({\mathbf{\beta }}_t\right)\right)$$

and which satisfies the constrain imposed in Equation (18) takes the form

$$\widehat{p}\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1},{\rho }_t)\propto p_{\delta }{\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})}^{{\rho }_t}{p}_u{\left({\mathbf{\beta }}_t\right)}^{(1-{\rho }_t)},$$

(19)

where $0\le {\rho }_t\le 1$ is indirectly defined by Equation (18) and therefore depends on the user-defined parameter κ.

Proof. 

The proof can be found in Appendix A. □

This approach defines transition models without having to make explicit assumptions about the parametric family. In consequence, it provides a generic off-the-self mechanism for defining transition models. In subsequent subsections, we provide an intuitive interpretation of this approach by establishing a direct relationship with exponential forgetting and power prior [28] approaches, which has been previously use to deal with nonstationary data streams. We deviate from the original formulation given in Kárny [31] and Ozkan et al. [30], which proposes to approximate $p\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})$ by the distribution $\widehat{p}\left({\mathbf{\beta }}_t\right|{\mathit{x}}_{1:t-1})$ having maximum entropy under the constraint in Equation (18). However, in this case, they require that $p_u\left({\mathbf{\beta }}_t\right)$ is defined as an invariant measure (i.e., ${\mathbb{E}}_q[ln{p}_u\left({\mathbf{\beta }}_t\right)]$ has to be constant w.r.t. any density $q\left({\mathbf{\beta }}_t\right)$). Note that, in case we employ a prior $p_u\left({\mathbf{\beta }}_t\right)$ satisfying this property, Theorem 1 also chooses the maximum entropy distribution.

In our streaming data setting, we follow the assumed density filtering approach [43] and the SVB approach [27] and use the approximation $p\left({\mathbf{\beta }}_{t-1}\right|{\mathit{x}}_{1:t-1})\approx q\left({\mathbf{\beta }}_{t-1}\right|{\mathit{\lambda }}_{t-1}),$ where $q\left({\mathbf{\beta }}_{t-1}\right|{\mathit{\lambda }}_{t-1})$ is the variational distribution calculated in the previous time step. Using this approximation in Equations (16) and (17), we can express $p_{\delta }$ in terms of ${\mathit{\lambda }}_{t-1}$, in which case, Equation (19) becomes

$$\widehat{p}\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1},{\rho }_t)\propto p_{\delta }{\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1})}^{{\rho }_t}{p}_u{\left({\mathbf{\beta }}_t\right)}^{(1-{\rho }_t)},$$

(20)

which we use as the prior density for time step t. Now, if $p_u\left({\mathbf{\beta }}_t\right)$ belongs to the same family as $q\left({\mathbf{\beta }}_{t-1}\right|{\mathit{\lambda }}_{t-1})$, then $\widehat{p}\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1},{\rho }_t)$ will stay within the same family and have natural parameters ${\rho }_t{\mathit{\lambda }}_{t-1}+(1-{\rho }_t){\mathbf{\alpha }}_u$, where ${\mathbf{\alpha }}_u$ are the natural parameters of $p_u\left({\mathbf{\beta }}_t\right)$. Therefore, we can write

$$ln\widehat{p}\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1},{\rho }_t)=lnh\left({\mathbf{\beta }}_t\right)+({\rho }_t{\mathit{\lambda }}_{t-1}+(1-{\rho }_t){\mathbf{\alpha }}_u)\mathit{t}\left({\mathbf{\beta }}_t\right)-a_g({\rho }_t{\mathit{\lambda }}_{t-1}+(1-{\rho }_t){\mathbf{\alpha }}_u)$$

(21)

where $h\left({\mathbf{\beta }}_t\right)$ is the base measure, which does not depend on any parameter. Thus, under this approach, the transitioned posterior remains within the same exponential family, so we can enjoy the full flexibility of the conjugate exponential family (i.e., computing gradients of the $\mathcal{L}$ function in closed form), an option that would not be available if one were to explicitly specify a transition model as in Equation (16).

Therefore, at each time step, we simply have to solve the following variational problem, where only the prior changes with respect to the original SVB approach:

$$arg\underset{{\mathit{\lambda }}_t,{\mathit{\varphi }}_t}{max}\mathcal{L}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t|{\mathit{x}}_t,{\rho }_t{\mathit{\lambda }}_{t-1}+(1-{\rho }_t){\mathbf{\alpha }}_u).$$

(22)

We shall refer to the method outlined in this section as SVB with power priors (SVB-PP). The term power prior [28] will be explained in Section 4.2.

4.1. Exponential Forgetting as Implicit Transition Models

In this section, we show that the exponential forgetting mechanism used in variational inference, as described in Section 3.4, is an implicit transition model with constant forgetting rate $\rho$.

The updating equation detailed in Equation (7) to optimize the lower-bound function described in Equation (6) can be easily adapted to optimize the lower-bound associated to the implicit transition models given in Equation (22). This new updating equation for implicit transition models can be expressed as follows:

$${\mathit{\lambda }}_t={\mathbb{E}}_{{\mathit{\varphi }}_t}\left[(t({\mathit{x}}_t,{\mathit{z}}_t),1)\right]+\rho {\mathit{\lambda }}_{t-1}+(1-\rho ){\mathbf{\alpha }}_u.$$

(23)

Expanding the above equation, we have

$${\mathit{\lambda }}_t=\sum _{i=1}^t{\rho }^{t-i}{\mathbb{E}}_{{\mathit{\varphi }}_i}\left[(t({\mathit{x}}_i,{\mathit{z}}_i),1)\right]+{\mathbf{\alpha }}_u,$$

(24)

which exactly matches exponentially down-weighting scheme of old data samples given in Equation (12). Therefore, it is clear that the classic technique of exponential forgetting, which was usually supported by heuristic arguments, has a sound interpretation in terms of implicit transition models.

4.2. Power Priors as Implicit Transition Models

Power priors [28] is a widely used class of informative priors for dealing with situations in which historical data are available. Let ${\mathit{x}}_0$ denote a previously obtained data set, and let ${\mathit{x}}_1$ be our current data set. According to the power priors scheme [28], the posterior probability over the model parameters should be computed as

$$p\left(\mathbf{\beta }\right|{\mathit{x}}_1,{\mathit{x}}_0,\rho )\propto p\left({\mathit{x}}_1\right|\mathbf{\beta })p{\left({\mathit{x}}_0\right|\mathbf{\beta })}^{\rho }p\left(\mathbf{\beta }\right),$$

(25)

where $\rho \in [0,1]$ is a scalar parameter down-weighting the likelihood of historical data relative to the likelihood of the current data.

As stated in the following lemma, power priors can also be interpreted as implicit transition models.

Lemma 1.

The Bayesian updating scheme described by Figure 1b and Equation (19) but with ${\rho }_t$ fixed to a constant value is equivalent to the recursive application of the Bayesian updating scheme of power priors given in Equation (25).

Proof. 

The proof can be found in Appendix A. □

This connection allows us to introduce well-known results of power priors [29], finding that

$$\begin{array}{ccc}\hfill p\left(\mathbf{\beta }\right|{\mathit{x}}_1,{\mathit{x}}_0,\rho )& =& arg\underset{r\in \mathcal{P}}{min}\left\{\rho \mathit{KL}\left(r\left|\right|p\left(\mathbf{\beta }\right|{\mathit{x}}_1,{\mathit{x}}_0,\rho =1)\right)\right.\hfill \\ & & +\left.(1-\rho )\mathit{KL}\left(r\left|\right|p\left(\mathbf{\beta }\right|{\mathit{x}}_1,{\mathit{x}}_0,\rho =0)\right)\right\}\hfill \end{array}$$

where $\mathcal{P}$ denotes the set of all possible densities over $\mathbf{\beta }$. Citing [29], “power priors minimize the convex combination of KL (Kullback–Leibler) divergences between two extremes: one in which no historical data is used and the other in which the historical data and current data are given equal weight.”

5. Hierarchical Power Priors

In the approach taken by Ozkan et al. [30] (and, by extension, SVB-PP), the forgetting factor ${\rho }_t$ is user-defined. In this paper, we instead pursue a (hierarchical) Bayesian approach and introduce a prior distribution over ${\rho }_t$, allowing the distribution over ${\rho }_t$ (and thereby the forgetting mechanism) to adapt to the data stream.

5.1. A Hierarchical Prior over the Forgetting Rate $\rho$

In this section, we extend the model in Figure 1a to also account for the dynamics of the data stream being modeled. We shall assume here that only the parameters $\mathbf{\beta }$ in Figure 1a are time-varying, which we will indicate with the subscript t, i.e., ${\mathbf{\beta }}_t$. The resulting model can be illustrated as in Figure 1b. We shall refer to models of this type as hierarchical power prior (HPP) models.

We will show in Section 5.3 that the exponential and normal distributions, both of which truncated to the interval [0,1], are valid alternatives as prior distributions, $p\left({\rho }_t\right|\gamma )$. The densities of these distributions have the following forms:

$$p\left({\rho }_t\right|\gamma )=\frac{\gamma exp(-\gamma {\rho }_t)}{1-exp(-\gamma )},0\le {\rho }_t\le 1,$$

(26)

$$p\left({\rho }_t\right|\mu ,\sigma )=\frac{exp\left(-{({\rho }_t-\mu )}^2/\left(2{\sigma }^2\right)\right)}{\sqrt{2\pi {\sigma }^2}\left(\Phi \left(\frac{1-\mu }{\sigma }\right)-\Phi \left(\frac{-\mu }{\sigma }\right)\right)},0\le {\rho }_t\le 1,$$

(27)

where $\Phi$ represents the standard normal cumulative distribution function, $\mu \in \mathbb{R}$ (can be outside the interval $[0,1]$), and $\sigma >0$. Since the natural parameters of the normal distribution are $\left(\mu /{\sigma }^2,-1/\left(2{\sigma }^2\right)\right)$, it is sometimes convenient to parameterize in terms of the precision $\eta =1/{\sigma }^2$ (the reciprocal of the variance). Using the precision, the natural parameters become $\left(\mu \eta ,-\eta /2\right)$. Notice that precision $\eta$ appears in both components.

Figure 2 plots examples of both densities for different values of the parameters. The truncated exponential allows to model a uniform prior and priors that either favors ${\rho }_t$ values close to 1 (i.e., non-forgetting past data) or ${\rho }_t$ values close to 0 (i.e., forgetting past data). The truncated normal distribution when using a mean equal to 0.5 tends to favor non-extreme ${\rho }_t$ values (i.e., partial forgetting of past data), where the variance parameter defines the strength of this belief.

For later use, we also detail here the equation for computing the expected value of ${\rho }_t$ for both distributions:

$$\mathbb{E}\left[{\rho }_t\right|\gamma ]=\frac{1}{(1-e^{-\gamma })}-\frac{1}{\gamma },$$

(28)

$${\mathbb{E}}_q\left[{\rho }_t\right|\mu ,\sigma ]=\mu +\sigma \frac{\mathit{\varphi }\left(\frac{-\mu }{\sigma }\right)-\mathit{\varphi }\left(\frac{1-\mu }{\sigma }\right)}{\Phi \left(\frac{1-\mu }{\sigma }\right)-\Phi \left(\frac{-\mu }{\sigma }\right)},$$

(29)

where $\gamma$ is the mean parameter parameter of the truncated exponential, $\mu$ and $\sigma$ are the parameters of the truncated normal distribution in $[0,1]$, and $\mathit{\varphi }$ and $\Phi$ are respectively the probability density function and the cumulative distribution function of the standard normal distribution.

5.2. The Double Lower Bound

For updating the model distributions, we pursue a variational approach where we seek to maximize the evidence lower bound $\mathcal{L}$ in Equation (5) for time step t. However, since the model in Figure 1b does not define a conjugate exponential distribution due to the introduction of $p\left({\rho }_t\right|\gamma )$, we cannot maximize $\mathcal{L}$ directly. Instead, we will derive a (double) lower bound $\widehat{\mathcal{L}}$ (with $\widehat{\mathcal{L}}\le \mathcal{L}$) and use this lower bound as a proxy for updating the rules of the variational posteriors.

First, by instantiating the lower bound ${\mathcal{L}}_{HPP}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t,{\mathit{\omega }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1})$ in Equation (5) for the HPP model, we obtain

$$\begin{array}{cc}\hfill {\mathcal{L}}_{HPP}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t,{\mathit{\omega }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1})& ={\mathbb{E}}_q[lnp\left({\mathit{x}}_t\right|{\mathit{Z}}_t,{\mathbf{\beta }}_t)]-{\mathbb{E}}_q\left[\mathit{KL}\left(q\left({\mathit{Z}}_t\right|{\mathit{\varphi }}_t)\left|\right|p\left({\mathit{Z}}_t\right|{\mathbf{\beta }}_t)\right)\right]\hfill \\ \hfill & -{\mathbb{E}}_q\left[\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|\widehat{p}\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1},{\rho }_t)\right)\right]\hfill \\ \hfill & -\mathit{KL}\left(q\left({\rho }_t\right|{\mathit{\omega }}_t)\left|\right|p\left({\rho }_t\right|\gamma )\right)\hfill \end{array}$$

(30)

where ${\mathit{\omega }}_t$ is the variational parameter of the variational distribution of ${\rho }_t$. For ease of presentation, we shall sometimes drop from ${\mathcal{L}}_{HPP}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t,{\mathit{\omega }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1})$ the subscript and the explicit specification of the parameters when this is otherwise clear from the context.

We now define the double lower bound ${\widehat{\mathcal{L}}}_{HPP}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t,{\mathit{\omega }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1})$ as

$$\begin{array}{cc}\hfill {\widehat{\mathcal{L}}}_{HPP}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t,{\mathit{\omega }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1})& ={\mathbb{E}}_q[lnp\left({\mathit{x}}_t\right|{\mathit{Z}}_t,{\mathbf{\beta }}_t)]-{\mathbb{E}}_q\left[\mathit{KL}\left(q\left({\mathit{Z}}_t\right|{\mathit{\varphi }}_t)\left|\right|p\left({\mathit{Z}}_t\right|{\mathbf{\beta }}_t)\right)\right]\hfill \\ \hfill & -{\mathbb{E}}_q\left[{\rho }_t\right]\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p({\mathbf{\beta }}_t|{\mathit{\lambda }}_{t-1}\right))\hfill \\ \hfill & -(1-{\mathbb{E}}_q\left[{\rho }_t\right])\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p\left({\mathbf{\beta }}_t\right|{\mathbf{\alpha }}_u)\right)\hfill \\ \hfill & -\mathit{KL}\left(q\left({\rho }_t\right|{\mathit{\omega }}_t)\left|\right|p\left({\rho }_t\right|\gamma )\right)\hfill \end{array}$$

(31)

which, according to Theorem 2, provides a lower bound for $\mathcal{L}$.

Theorem 2.

Let ${\widehat{\mathcal{L}}}_{HPP}$ be as defined in Equation (31). Then,

$${\widehat{\mathcal{L}}}_{HPP}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t,{\mathit{\omega }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1})\le {\mathcal{L}}_{HPP}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t,{\mathit{\omega }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1}).$$

Proof. 

The proof can be found in Appendix A. □

Even though Equation (31) defines an alternative objective function, when we compare this double lower bound with Equation (6), we can observe that the double lower bound still has the intuitive interpretation of the standard lower bound in terms of data fitting and Kullback–Leibler (KL) regularization. The only difference is that the KL regularization term associated to $q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)$ appears now as a convex combination of two KL terms, one regularizing with respect to $p\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1})$ and the other regularizing with respect to $p\left({\mathbf{\beta }}_t\right|{\mathbf{\alpha }}_u)$, whith ${\mathbb{E}}_q\left[{\rho }_t\right]$, acting as a combination factor.

Rather than seeking to maximize $\mathcal{L}$, we will instead maximize $\widehat{\mathcal{L}}$; see Equation (31). Thus, maximizing $\widehat{\mathcal{L}}$ w.r.t. the variational parameters ${\mathit{\lambda }}_t$ and $\mathit{\varphi }$ also maximizes $\mathcal{L}$. By the same observation, we also have that the (natural) gradients are consistent relative to the two bounds, as stated by the next corollary.

Corollary 1.

$$\begin{array}{cc}\hfill {\nabla }_{{\mathit{\lambda }}_t}^{nat}{\widehat{\mathcal{L}}}_{HPP}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t,{\mathit{\omega }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1})& ={\nabla }_{{\mathit{\lambda }}_t}^{nat}{\mathcal{L}}_{HPP}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t,{\mathit{\omega }}_t|{\mathit{x}}_t,{\mathit{\lambda }}_{t-1})\hfill \\ & ={\nabla }_{{\mathit{\lambda }}_t}^{nat}\mathcal{L}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t|{\mathit{x}}_t,{\mathbb{E}}_q\left[{\rho }_t\right]{\mathit{\lambda }}_{t-1}+(1-\mathbb{E}\left[{\rho }_t\right]){\mathbf{\alpha }}_u)\hfill \end{array}$$

The same result holds for ${\mathit{\varphi }}_t$.

Proof. 

The proof can be found in Appendix A. □

Thus, updating the variational parameters ${\mathit{\lambda }}_t$ and ${\mathit{\varphi }}_t$ in HPP models can be done in the same way as for regular conjugate exponential models of the form in Figure 1. A pseudo-code description of the updating process can be found in Algorithm 1 when ${\rho }_t$ is assumed to follow a truncated exponential distribution.

Algorithm 1 Streaming variational Bayes (SVB) with Hierarchical Power Priors (SVB-HPP).

Input: A data batch xt, the variational posterior in previous time step λt−1.

Output: (λt, ϕt, ωt), a new update of the variational posterior.

1: λt ← λt−1.

2: ${\mathbb{E}}_q\left[{\rho }_t\right]\leftarrow 0.5$.

3: Randomly initialize ϕt.

4: repeat

5:  (λt, ϕt) = arg minλt,ϕt $\mathcal{L}({\mathit{\lambda }}_t,{\mathit{\varphi }}_t|{\mathit{x}}_t,{\mathbb{E}}_q\left[{\rho }_t\right]{\mathit{\lambda }}_{t-1}+(1-\mathbb{E}\left[{\rho }_t\right]){\mathbf{\alpha }}_u)$

6:  ${\mathit{\omega }}_t=\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p_u\left({\mathbf{\beta }}_t\right)\right)-\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p_{\delta }\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1})\right)+\gamma$

7:  Update ${\mathbb{E}}_q\left[{\rho }_t\right]$ according to Equation (28) or Equation (29).

8: until convergence

9: return (λt,ϕt,ωt)

In order to update ${\mathit{\omega }}_t$, we rely on $\widehat{\mathcal{L}}$, which we can maximize using the natural gradient w.r.t. ${\mathit{\omega }}_t$ [44] and which can be calculated in closed form for a restricted distribution family for ${\rho }_t$, as stated in the following result.

Lemma 2.

Assuming that the first component of the sufficient statistics function for ${\rho }_t$ is the identity function, i.e., ${\mathit{t}}_1\left({\rho }_t\right)={\rho }_t$, it holds that

$$\begin{array}{cc}\hfill \frac{{\partial }^{nat}\widehat{\mathcal{L}}}{\partial {\mathit{\omega }}_{t,1}}=& \mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p_u\left({\mathbf{\beta }}_t\right)\right)-\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p_{\delta }\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1})\right)+{\gamma }_1-{\mathit{\omega }}_{t,1},\hfill \\ \hfill \frac{{\partial }^{nat}\widehat{\mathcal{L}}}{\partial {\mathit{\omega }}_{t,k}}=& {\gamma }_k-{\mathit{\omega }}_{t,k}(k\ne 1).\hfill \end{array}$$

(32)

Proof. 

The proof can be found in Appendix A. □

From the above lemma, we can easily deduce that the truncated exponential distributions, for which the sufficient statistics are $t\left({\rho }_t\right)={\rho }_t$, and the truncated normal distribution, for which the sufficient statistics are $t\left({\rho }_t\right)={({\rho }_t,{\rho }_t^2)}^T$, satisfy the criteria to be considered as hierarchical priors for ${\rho }_t$.

The problem of the above result is that, for $k\ne 1$, the optimal ${\mathit{\omega }}_{t,k}$ is just equal to the the prior value, i.e., ${\mathit{\omega }}_{t,k}={\gamma }_k$. In the case of the truncated normal, which has a two-dimensional natural parameter vector, it would imply that the variance of the posterior $q\left({\rho }_t\right|{\mathit{\omega }}_t)$, denoted by ${\sigma }_q^2$, will be equal to the variance of prior, denoted by ${\sigma }_p^2$, which has to be set manually (We dropped the t-index in ${\sigma }_q^2$ for simplicity.). To address the issue of having to manually fix the variance of the truncated normal prior, ${\sigma }_p^2$, we employ an empirical Bayes approach and consider ${\sigma }_p^2$ as another free parameter of the double lower bound that we want to optimize. Therefore, we need to compute the gradient of the double lower bound w.r.t. this parameter:

$$\begin{array}{ccc}\hfill \frac{\partial \widehat{\mathcal{L}}}{\partial {\sigma }_p^2}& =& \frac{\partial {\gamma }_1}{\partial {\sigma }_p^2}\frac{\partial \widehat{\mathcal{L}}}{\partial {\gamma }_1}+\frac{\partial {\gamma }_2}{\partial {\sigma }_p^2}\frac{\partial \widehat{\mathcal{L}}}{\partial {\gamma }_2}\hfill \\ & =& -\frac{{\mu }_p}{{\sigma }_p^4}(\mathbb{E}\left[{\rho }_t\right|{\mu }_q,{\sigma }_q^2]-\mathbb{E}\left[{\rho }_t\right|{\mu }_p,{\sigma }_p^2])+\frac{1}{2{\sigma }_p^4}(\mathbb{E}\left[{\rho }_t^2\right|{\mu }_q,{\sigma }_q^2]-\mathbb{E}\left[{\rho }_t^2\right|{\mu }_p,{\sigma }_p^2]),\hfill \end{array}$$

where $\gamma =({\mu }_p/{\sigma }_p^2,-1/\left(2{\sigma }_p^2\right))$ is the natural parameter vector of the truncated normal prior, and ${\mu }_p$ and ${\mu }_q$ denote the mean of the truncated normal prior and posterior over ${\rho }_t$, respectively. We set ${\mu }_p$ to 0.5, trying to define a non-informative and symmetric prior.

Note that, in this case, we have a plain gradient instead of a natural gradient w.r.t. ${\sigma }_p^2$. Also, note that there is no closed-form solution for the stationary point of ${\sigma }_p^2$. To optimize along this direction, we use a simple gradient ascent approach with backtracking line-search to set the learning rate.

5.3. Towards a Measure of Concept Drift

Observe that the form of the natural gradient of ${\mathit{\omega }}_t$ given in Lemma 2 has an intuitive semantic interpretation in terms of measure of concept drift. If we follow a coordinate ascent algorithm, at every iteration, we should set

$${\mathit{\omega }}_t=\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p_u\left({\mathbf{\beta }}_t\right)\right)-\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p_{\delta }\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1})\right)+\gamma .$$

(33)

Specifically, using the constant $\gamma$ as a threshold, we see that, if the uniform prior $p_u\left({\mathbf{\beta }}_t\right)$ is closer to the variational posterior at time t, in terms of KL divergence than the variational posterior at the previous time step (i.e., $\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p_u\left({\mathbf{\beta }}_t\right)\right)+\gamma <\mathit{KL}\left(q\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_t)\left|\right|p_{\delta }\left({\mathbf{\beta }}_t\right|{\mathit{\lambda }}_{t-1})\right)$), then we will get a negative value for ${\mathit{\omega }}_t$.

This in turn implies that ${\mathbb{E}}_q\left[\rho \right]<0.5$, according to Equations (28) and (29) (plotted in Figure 3), which means that we have a higher degree of forgetting for past data. If ${\mathit{\omega }}_t>0$, then ${\mathbb{E}}_q\left[\rho \right]>0.5$ and less past data is forgotten. Figure 3 (left) graphically illustrates this trade-off, while Figure 4 shows a particular example of this situation using Gaussian posterior distributions.

The difference between the use of a truncated normal over a truncated exponential is that, with the former, the relation between ${\mathit{\omega }}_t$ and ${\mathbb{E}}_q\left[{\rho }_t\right]$ can be tuned by a change in the precision of the truncated normal prior, as it is graphically illustrated in Figure 3 (right). By using a prior with higher precision, we impose a stronger belief about the fact that ${\rho }_t$ values are close to neither 1 nor 0. In this way, the approach has the possibility to enforce smooth drift regimes.

5.4. The Multiple Hierarchical Power Prior Model

In this section, we propose a modification of our HPP model to deal with complex concept drift patterns which involve only a part of the model. For example, let us consider the application of an LDA model [4] for tracking over time the evolution of topics in a text corpus. Under these settings, a drift could eventually affect only a subset of the topics. Using our current approach, we might detect this drift and forget part of the data to adapt to the new situation and to learn the new topics. However, if some topics have not changed, we are losing information that could provide better estimations for these topics.

We propose an immediate extension of HPP, which include multiple power priors ${\rho }_t^{\left(i\right)}$, one for each parameter ${\mathbf{\beta }}_i$. In this model, the ${\rho }_t^{\left(i\right)}$s are pair-wise independent. The latter ensures that optimizing $\widehat{\mathcal{L}}$ can be performed as above, since the variational distribution for each ${\rho }_t^{\left(i\right)}$ can be updated independently of the other variational distributions over ${\rho }_t^{\left(j\right)}$ for $j\ne i$. This extended model allows local model substructures to have different forgetting mechanisms, thereby extending the expressiveness of the model. We shall refer to this extended model as a multiple hierarchical power prior (MHPP) model.

6. Results

In this section we will evaluate the following methods:

• Streaming variational Bayes (SVB) as described in Section 3.3.
• Four versions of Population Variational Bayes (PVB) (We do not compare with SVI because SVI is a special case of PVB when M is equal the total size of the stream.) resulting from combining the values of the population size parameter $M=1000$ (Section 6.1) and $M=10,000$ (Section 6.2) with the learning rate values $\nu =0.1$ and $\nu =0.01$. In the four cases, the mini-batch size was set to 1000. Note that, however, for the LDA case, we set $M=1000$ rather than $M=10,000$ in Section 6.2 and use a mini-batch size of 1000 instead of $10,000$.
• Two versions of the SVB method with power priors (SVB-PP) or fixed exponential forgetting (as described in Section 4.1) with $\rho =0.9$ or $\rho =0.99$.
• Three versions of our method based on the SVB method with adaptive exponential forgetting using hierarchical power priors (as described in Section 5):
  • SVB-HPP-Exp using a single shared $\rho$ with a truncated exponential distribution as a prior over $\rho$ with $\gamma =0.1$ (i.e., close to uniform).
  • SVB-MHPP-Exp using separate ${\rho }^{\left(i\right)}$ for each parameter (as described in Section 5.4) with truncated exponential distributions as priors over each ${\rho }^{\left(i\right)}$ with $\gamma =0.1$.
  • SVB-MHPP-Norm using also separate ${\rho }^{\left(i\right)}$ for each parameters but with truncated normal distributions as priors over each ${\rho }^{\left(i\right)}$. In this case, we use ${\mu }_p=0.5$ and learn the variance ${\sigma }_p^2$ using the empirical Bayes approach described at the end of Section 5.2.

The underlying variational inference method used in the experiments is the variational message passing (VMP) algorithm [41] for all models; VMP was terminated after 100 iterations or if the relative increase in the lower bound fell below $0.01\%$ ($0.0001\%$ for LDA). We considered non-informative priors, i.e., flat Gaussians, flat Gamma, or uniform Dirichlet (full details in Appendix B). For the LDA model [4], we use standard priors for this model which include a Dirichlet prior over topics with $\alpha =\frac{1}{\left|V\right|}$, where $\left|V\right|$ denotes the size of the vocabulary and another Dirichlet prior over topics assignments with $\alpha =0.1$. Variational parameters were randomly initialized using the same seed for all methods.

6.1. Evaluation Using an Artificial Data Set

First, we illustrate the behavior of the different approaches in a controlled experimental setting: We produced an artificial data stream by generating 100 samples (i.e., $|{\mathit{x}}_t|=100$) from a Binomial distribution at each time step. We artificially introduced concept drift by changing the parameter p of the Binomial distribution: $p=0.2$ for the first 30 time steps, then $p=0.5$ for the following 30 time steps, and finally $p=0.8$ for the last 40 time steps. The data stream was modeled using a beta-binomial model.

Parameter estimation:Figure 5 shows the evolution of ${\mathbb{E}}_q\left[{\mathbf{\beta }}_t\right]$ for the different methods. We recognize that SVB simply generates a running average of the data, as it is not able to adapt to the concept drift. The results from PVB strongly depend on the learning rate $\nu$, where the highest learning rate, which results in the more aggressive forgetting, works better in this example. Recall, however, that $\nu$ has to be hand-tuned to achieve optimal performance. As expected, the choice of the size of the population M for SVB does not have an impact because the present model has no local hidden variables (cf. Section 4.1). SVB-PP yields results almost identical to PVB when $\rho$ matches the learning rate of PVB (i.e., $\rho =1-\nu$). Finally, SVB-HPP provides the best results, almost mirroring the true model.

Equivalent sample size of the posterior ($ES{S}_{post}$):Figure 6 (left) gives the evolution of the equivalent sample size of the posterior, $ES{S}_{post,t}$, for the different methods (For this model, $ES{S}_{post,t}$ is simply computed by summing the components of ${\lambda }_t$ defining the beta posterior.). $ES{S}_{post,t}$ of PVB is always given by the constant M. For SVB, $ES{S}_{post,t}$ monotonically increases as more data is used, while SVB-PP exhibits convergence to the limiting value computed in Equation (13). A different behaviour is observed for SVB-HPP: It is automatically adjusted. Notice that the values for this model are to be read off the alternative y-axis. We can detect the concept drift by identifying where ESS rapidly declines.

Evolution of expected forgetting factor: In Figure 6 (right), the series denoted “$E\left[\rho \right]-100$” shows the evolution of ${\mathbb{E}}_q\left[{\rho }_t\right]$ for the artificial data set. Notice how the model clearly identifies abrupt concept drift at time steps $t=30$ and $t=60$. The series denoted “$E\left[\rho \right]-1000$” illustrates the evolution of the parameter when we increase the batch size to 1000 samples. We recognize a more confident assessment about the absence of concept drift as more data is made available.

6.2. Evaluation Using Real Data Sets

For this evaluation, we consider four real data sets from four different domains:

Electricity market [45]: The data set describes the electricity market of two Australian states. It contains 45,312 instances of 6 attributes, including a class label comparing change of the electricity price related to a moving average of the last 24 h. Each instance in the data set represents 30 min of trading; during our analysis, we created batches such that ${\mathit{x}}_t$ contains all information associated with month t.

The data is analyzed using a Bayesian linear regression model. The binary class label is assumed to follow a Gaussian distribution in order to fit within the conjugate model class. Similarly, the marginal densities of the predictive attributes are also assumed to be Gaussian. The regression coefficients are given Gaussian prior distributions, and the variance is given a Gamma prior. Note that the overall distribution does not fall inside the conditional conjugate exponential family [26]; hence, we do not apply SVI (and PVB) in this setting.

GPS [46,47,48]: This data set contains $17,621$ GPS trajectories (time-stamped x and y coordinates), totalling more than 4.5 million observations. To reduce the data size, we kept only one out of every ten measurements. We grouped the data so that ${\mathit{x}}_t$ contains all data collected during hour t of the day, giving a total of 24 batches of this stream.

Here, we employ a model with one independent Gaussian mixture model per day of the week, each mixture with 5 components. This enables us to track changes in the users’ profiles across hours of the day and to monitor how the changes are affected by the day of the week.

Finance [5]: The data contains monthly aggregated information about the financial profile of around $50,000$ customers over 62 (nonconsecutive) months. Three attributes were extracted per customer in addition to a class-label indicating whether the customer will default within the next 24 months.

We fit a naïve Bayes model to this data set, where the distribution at the leaf nodes is a 5-component mixture of Gaussians. The distribution over the mixture node is shared by all attributes but not between the two classes of customers.

NIPS [36]: This data set consists of the abstracts of published papers in the NIPS (Neural Information Processing Systems) conference between 1987 and 2015 (5804 documents in total). The data was preprocessed by choosing the most relevant individual terms across the whole dataset. This was done by ordering the words (11,463 in total) by their importance in the dataset, using the TF-IDF (term frequency-inverse document frequency) metric. The top 10 words after this filtering were “policy", “image", “kernel", “network", “neurons", “training", “graph", “images", “matrix", and “tree", while the last 5 words in the ranking were “ralf", “ciated", “havior", “references", and “abstract". Only the top 100 words were kept, according to this criterion. In this way, we removed words that were not significant to tracking the concept drift in this data set. The documents were grouped by year, yielding a total of 29 batches of documents of different sizes. An LDA model with ten topics was employed to analyze the vocabulary and to detect changes in the evolution of the major topics of the papers of this conference every year. Note that the temporal extension of this model involves dealing with dynamics at the Dirichlet distributions over the topics. As commented in Section 2, there have been many previous approaches trying to deal with this problem [34,35,36], but none of them were applicable to general conjugate exponential family models and, in general, rely on much more complex inference schemes.

To evaluate the different methods discussed, we use the test marginal log-likelihood (TMLL). Specifically, each data batch is randomly split into a train data set, ${\mathit{x}}_t$, and a test data set, ${\tilde{\mathit{x}}}_t$, containing two thirds and one third of the data batch, respectively. Then, TMLL_t is computed as ${\mathrm{TMLL}}_t=\frac{1}{|{\tilde{\mathit{x}}}_t|}\int p({\tilde{\mathit{x}}}_t,{\mathit{z}}_t|{\mathbf{\beta }}_t)p\left({\mathbf{\beta }}_t\right|{\mathit{x}}_t)d{\mathit{z}}_{t}d{\mathbf{\beta }}_t$ (For LDA, $|{\tilde{\mathit{x}}}_t|$ refers to the number of words in the test set; we then compute the so-called per-word perplexity.).

A detailed description of all models, including their structure and their variational families, is given in Appendix B.

6.3. Discussion of Results

In this first part, we will highlight how the basic versions of SVB-HPP and SVB-MHPP outperform the rest of the approaches in most of the cases. Figure 7 shows for each method the difference between its TMLL_t and that obtained by SVB (which is considered the baseline method). To improve readability, we only plot the results of the best performing method inside each group of methods. Figure 8 shows the development of ${\mathbb{E}}_q\left[{\rho }_t\right]$ over time for SVB-HPP-Exp, SVB-MHPP-Exp, and SVB-MHPP-Norm. For SVB-HPP-TExp, we only have one ${\rho }_t$-parameter, and its value is given by the solid line. SVB-MHPP utilizes one ${\rho }^{\left(i\right)}$ for each variational parameter. (The numbers of variational parameters are 14, 78, 33, and 10 for the electricity, GPS, financial, and NIPS models, respectively.) In this case, we plot ${\mathbb{E}}_q\left[{\rho }_t^{\left(i\right)}\right]$ at each point in time to indicate the variability between the different estimates throughout the series. We also report the average of the ${\mathbb{E}}_q\left[{\rho }_t^{\left(i\right)}\right]$ values at every time step. Finally, we compute the aggregated test marginal log-likelihood measure ${\sum }_{t=1}^T{\mathrm{TMLL}}_t$ for each method and report these values in

Table 1. Aggregated test marginal log-likelihood: maximum values for each data set are boldfaced.

Data Set	SVB	PVB				SVB-PP		SVB-HPP	SVB-MHPP
		(1)	(2)	(3)	(4)	$\mathit{\rho }=0.9$	$\mathit{\rho }=0.99$	Exp	Exp	Norm
Electricity	−44.91	−51.01	−52.19	−51.11	−61.70	−43.92	−44.80	−40.05	−40.02	−39.91
GPS	−1.98	−2.10	−2.77	−1.97	−4.49	−1.94	−1.97	−1.97	−1.86	−1.86
Finance	−19.84	−22.29	−22.57	−20.40	−20.73	−19.05	−19.78	−19.83	−19.83	−19.82
NIPS	−4.07	−4.04 *	−4.21 *	−4.01	−4.12	−4.02	−4.06	−4.01	−4.00	−4.00

Population variational Bayes (PVB) parameters: (1) $M=10k,\nu =0.1$; (2) $M=10k,\nu =0.01$; (3) $M=|{\mathit{x}}_t|,\nu =0.1$; and (4) $M=|{\mathit{x}}_t|,\nu =0.01$. (*) For NIPS, $M=1k$ was used in (1) and (2).

.

For the electricity data set, we can see that the two proposed methods (SVB-HPP and SVB-MHPP) perform the best, as shown in Table 1. According to Figure 7 (electricity plot), all models are comparable during the first nine months, which is a period where our models detected no or very limited concept drift. However, after this period, both SVB-HPP and SVB-MHPP detected substantial drift and were able to adapt better than the other methods, which appeared unable to adjust to the complex concept drift structure in the latter part of the data. SVB-HPP and SVB-MHPP continued to behave at a similar level, mainly because when drift happened, it typically involved a high proportion of the parameters of the model (see the electricity plot in Figure 8).

For the GPS data set, we can observe how the SVB-MHPP perfoms quite well (see Table 1), particularly towards the end of the series (see the GPS plot in Figure 7). We can see in Figure 8 (GPS plot), that a significant proportion of the model parameters drifted (i.e., ${\mathbb{E}}_q\left[{\rho }_t^{\left(i\right)}\right]\le 0.05$) at all times, while another proportion of the parameters showed a quite stable behavior ($\rho$-values above $0.9$). This complex pattern is not properly captured by SVB-HPP, which ends up assuming no concept drift. PVB with $M=|{\mathit{x}}_t|$ and $\nu =0.1$ does well here, but it strongly depends on the hyperparameters, as in fact any other hyperparameter combination yields poorer results (see Table 1).

The financial data set shows a different behavior. As can be seen in Figure 7 (Finance plot), during the first months, no major differences among the methods were found. However, after month 30, SVB-PP with $\rho =0.9$ is superior. Looking at the $E\left[{\rho }_t^{\left(i\right)}\right]$-values of SVB-MHPP in Figure 8 (Finance plot), we observe that there is remarkable concept drift in some of the parameters over the first few months. However, only a few parameters exhibit noteworthy drift after the first third of the sequence. Apparently, the simple SVB-PP approach has the upper hand when drift is constant and fairly limited, at least when the optimal forgetting factor $\rho$ has been identified.

In the case of the NIPS data set, we see again that HPP approaches capture drift in the data. As can be seen in Figure 7 (NIPS plot), SVB-HPP hardly detects any drift in the first 20 years and performs quite similarly to SVB (i.e., relative performance close to zero). However, in the last 10 years, SVB-HPP clearly outperforms SVB because it detects two strong drifts at years 23 and 29. Therefore, at these time steps, the whole LDA model is almost reestimated from scratch. In this case, SVB-MHPP is able to capture more fine-grained drifts in the data, as can been deduced from the $E\left[{\rho }_t^{\left(i\right)}\right]$-values of SVB-MHPP in Figure 8 (NIPS plot). Mainly, it detects changes in some topics while other topics remain constant over time. This allows SVB-MHPP to outperform SVB-HPP during some periods.

We have also observed that there are no major differences between SVB-MHPP-Exp and SVB-MHPP-Norm. Therefore, it seems that the inclusion of alternative priors does not have a big impact on the performance, at least if the variance parameter of the truncated normal is fixed automatically using an empirical Bayes approach. However, this is something that may require further investigations.

We conclude this section by highlighting that the performances of SVB-PP and PVB strongly depends on the hyperparameters of the model, cf. Table 1. As an example, consider SVB-PP and how its performance varies by changing the $\rho$ parameter. Similarly, PVB’s performance is sensitive both to $\nu$ (see in particular the results for the GPS data) and M (financial data). These hyperparameters are hard to fix, as their optimal values depend on data characteristics (see McInerney et al. [25] and Broderick et al. [27] for similar conclusions). We, therefore, believe that the fully Bayesian formulation is an important strong point of our approach.

7. Conclusions

In this paper, we have introduced a novel Bayesian approach for learning general latent variable models from nonstationary data streams. For this purpose, we introduce implicit transition models as a general method for transitioning the parameters of a latent variable model. We also show that previous approaches like exponential forgetting and power priors can be seen as specific cases of this general transition model. However, these approaches are only able to model slowly changing data streams. Our approach is able to handle both abrupt and gradual drifts in the data stream by explicitly modeling the rate of change of the data stream. For this purpose, we introduce a novel hierarchical prior which allows the model to adapt to the different drifts that one can encounter in a data stream. We then develop an efficient variational inference scheme that optimizes a novel lower bound of the likelihood function.

As future work, we aim to provide a sound approach to semantically characterize concept drift by inspecting the $\mathbb{E}\left[{\rho }_t^{\left(i\right)}\right]$ values provided by SVB-MHPP and to investigate the effects in variational approximation introduced by the use of the double lower bound approximation.