Streaming Bayesian inference: theoretical limits and mini-batch approximate message-passing

TL;DR

This paper analyzes the fundamental limits of mini-batch streaming algorithms in large-scale statistical models, introduces Mini-AMP, and demonstrates its effectiveness in real-world clustering tasks.

Contribution

It provides an information-theoretic analysis of mini-batch inference limits and introduces Mini-AMP, a new algorithm with proven optimality in certain settings.

Findings

Mini-AMP achieves near-optimal performance in theoretical models.

Mini-AMP is competitive with existing streaming algorithms in clustering.

The analysis reveals phase transitions depending on mini-batch size.

Abstract

In statistical learning for real-world large-scale data problems, one must often resort to "streaming" algorithms which operate sequentially on small batches of data. In this work, we present an analysis of the information-theoretic limits of mini-batch inference in the context of generalized linear models and low-rank matrix factorization. In a controlled Bayes-optimal setting, we characterize the optimal performance and phase transitions as a function of mini-batch size. We base part of our results on a detailed analysis of a mini-batch version of the approximate message-passing algorithm (Mini-AMP), which we introduce. Additionally, we show that this theoretical optimality carries over into real-data problems by illustrating that Mini-AMP is competitive with standard streaming algorithms for clustering.

Equations

$$P(\mathbf{x}|{\mathbf{y}}^{(k)},{\mathcal{D}}^{(k-1)})=\frac{P({\mathbf{y}}^{(k)}|\mathbf{x})P(\mathbf{x}|{\mathcal{D}}^{(k-1)})}{\int d\mathbf{x}P({\mathbf{y}}^{(k)}|\mathbf{x})P(\mathbf{x}|{\mathcal{D}}^{(k-1)})}.$$

$$P(\bm{x}|\mathcal{D}^{(k-1)})\approx Q^{(k-1)}(\bm{x})\!=\!\prod_{i=1}^{N}\!\mathcal{P}_{i}\!\big{[}P(\bm{x}|\bm{y}^{(k-1)},\mathcal{D}^{(k-2)})\big{]},$$

$$P(\mathbf{y},\mathbf{x}|\Phi )=\prod _{\mu =1}^{M}P(y_{\mu }|z_{\mu }\equiv {\mathbf{\Phi }}_{\mu }\cdot \mathbf{x})\prod _{i=1}^N{P}_X(x_i).$$

$$P(x_{i}|\Phi,\bm{y})\approx q(x_{i}|A,B_{i})=P_{X}(x_{i})\,e^{-\frac{1}{2}Ax_{i}^{2}+B_{i}x_{i}}\big{/}{Z(A,B_{i})}\,,$$

$${\mathbf{z}}^{(t)}=\mathbf{y}-\Phi {\hat{\mathbf{x}}}^{(t)}+{\alpha }^{-1}{\mathbf{z}}^{(t-1)}{A}^{(t-1)}{V}^{(t)},$$

$${\mathbf{B}}^{(t)}=A^{(t)}{\hat{\mathbf{x}}}^{(t)}+A^{(t)}{\alpha }^{-1}{\Phi }^T{\mathbf{z}}^{(t)},$$

$${\hat{x}}_i^{(t+1)}=\eta (A^{(t)},B_i^{(t)})\forall i,$$

$$A^{(t)}⇝{\mathcal{A}}^{(t)}=\frac{\alpha }{\Delta +{\mathcal{V}}^{(t)}},B_i^{(t)}⇝{\mathcal{B}}^{(t)}\sim \mathcal{N}\left({\mathcal{A}}^{(t)}x,\alpha \frac{{\Delta }_0+\mathcal{E}(t)}{(\Delta +\mathcal{V}(t){)}^2}\right).$$

$${\mathcal{A}}^{(t)}=\frac{\alpha }{\Delta +{\mathcal{E}}^{(t)}},{\mathcal{E}}^{(t)}=\mathbb{E}{\eta }^{\prime }({\mathcal{A}}^{(t-1)},{\mathcal{A}}^{(t-1)}x+\sqrt{{\mathcal{A}}^{(t-1)}}z),$$

$$i_{\text{RS}}(\mathcal{E})=\frac{\alpha }{2}\left[\frac{\mathcal{E}}{\Delta +\mathcal{E}}+\log \left(1+\frac{\mathcal{E}}{\Delta }\right)\right]-{\mathbb{E}}_{x,z}\left[\log {\mathbb{E}}_{\tilde{x}}{e}^{-\frac{{\left(\tilde{x}-\left(x+z\Sigma (\mathcal{E})\right)\right)}^2}{2{\Sigma }^2(\mathcal{E})}}\right]-\frac{1}{2},$$

$$P_{{\Lambda }_{k-1},{\mathbf{\Theta }}_{k-1}}^k(\mathbf{x})=P_X(\mathbf{x})\prod _{i=1}^N{e}^{-\frac{1}{2}{\Lambda }_{k-1}{x}_i^2+{\Theta }_{k-1,i}{x}_i},\text{where }{\Lambda }_{k-1}=\sum _{\ell =1}^{k-1}{A}_l,{\Theta }_{k-1,i}=\sum _{\ell =1}^{k-1}{B}_{l,i}.$$

$${\hat{x}}_{k,i}^{(t+1)}$$

$$V_k^{(t+1)}$$

$${\lambda }_k^{(t)}$$

$${\gamma }_k^{(t)}$$

$${\lambda }_k^{(t)}$$

$$i_{\mathrm{RS}}^b({\mathcal{E}}_k)=\frac{{\alpha }_b}{2}\left[\frac{{\mathcal{E}}_k}{\Delta +{\mathcal{E}}_k}+\log \left(1+\frac{{\mathcal{E}}_k}{\Delta }\right)\right]-{\mathbb{E}}_{x,z}\log {\mathbb{E}}_{\tilde{x}}{e}^{-\frac{{\left(\tilde{x}-\left(x+z{\Sigma }_b(\lambda ,{\mathcal{E}}_k)\right)\right)}^2}{2{\Sigma }_b^2(\lambda ,{\mathcal{E}}_k)}}-\frac{1+{\alpha }_b}{2}.$$

$$I(Y,R;X)=H(Y,R)-H(Y|X)-H(R|X)={\mathbb{E}}_{\mathbf{y},\mathbf{r}}\log {\mathbb{E}}_{\mathbf{x}}{e}^{-\frac{\Vert \mathbf{y}-\Phi \mathbf{x}{\Vert }_2^2}{2\Delta }}{e}^{-\frac{\Vert \mathbf{r}-\mathbf{x}{\Vert }_2^2}{2{\lambda }^{-1}}}-\frac{1+{\alpha }_b}{2}.$$

$$\underset{N\to \infty }{\lim }\frac{I(Y,R;X)}{N}=i_{\text{RS}}^b({\mathcal{E}}_k)-{\int }_0^{\tau }{R}_{{\mathcal{E}}_k}(t)\mathrm{d}t+O(1),$$

$$\mathrm{MMSE}=\mathrm{argmin}{\mathrm{i}}_{\mathrm{RS}}^{\mathrm{b}}({\mathcal{E}}_{\mathrm{k}}).$$

$${\mathrm{MSE}}_x(\alpha )\sim e^{-\frac{1}{{\alpha }_{\text{b}}}\log (1-\frac{{\alpha }_{\text{b}}}{\rho })\alpha },$$

$$P(U,V,Y)=\prod _{ij}\mathcal{N}(Y_{ij};{\mathbf{U}}_i\cdot {\mathbf{V}}_j,\Delta )\prod _{i=1}^N\mathcal{N}({\mathbf{U}}_i;0,{\mathbb{I}}_R)\prod _{j=1}^M\frac{1}{R}\sum _{k=1}^R\delta ({\mathbf{V}}_j-\text{mathbbm}{e}_k),$$

$$P(\mathbf{y},\mathbf{x}|\Phi )=\prod _{\mu =1}^{M}P(y_{\mu }|z_{\mu }\equiv {\mathbf{\Phi }}_{\mu }\cdot \mathbf{x})\prod _{i=1}^N{P}_X(x_i).$$

$$q_x(x_i|A_i,B_i)=\frac{1}{Z_x(A_i,B_i)}{P}_X(x_i)e^{-\frac{1}{2}{A}_i{x}_i^2+B_i{x}_i},$$

$$q_z(z_{\mu }|y_{\mu },{\omega }_{\mu },V_{\mu })=\frac{1}{Z_z(y_{\mu },{\omega }_{\mu },V_{\mu })}P(y_{\mu }|z_{\mu })\frac{e^{-\frac{(z_{\mu }-{\omega }_{\mu })}{V_{\mu }}}}{\sqrt{2\pi V_{\mu }}}.$$

$${\mathbf{\omega }}^{(t)}=\Phi {\hat{\mathbf{x}}}^{(t)}-{\mathbf{V}}^{(t)}\circ {\mathbf{g}}^{(t-1)},$$

$$g_{\mu }^{(t)}=g_{\text{out}}(y_{\mu },{\omega }_{\mu }^{(t)},V_{\mu }^{(t)})\forall \mu ,$$

$${\mathbf{B}}^{(t)}={\Phi }^T{\mathbf{g}}^{(t)}+{\mathbf{A}}^{(t)}\circ {\hat{\mathbf{x}}}^{(t)},$$

$${\hat{x}}_i^{(t+1)}=\eta (A_i^{(t)},B_i^{(t)})\forall i,$$

$$V^{(t)}\approx \frac{1}{N}\sum _{i=1}^N{\sigma }^{(t)}(A^{(t)},B_i^{(t)}),A^{(t)}\approx \frac{\alpha }{\Delta +V^{(t)}},$$

$$KL[{\prod }_i{q}_i(x_i)\Vert P(\mathbf{x}|\Phi ,\mathbf{y})]=-{\mathbb{E}}_{\{q_i\}}\log P(\mathbf{y}|\Phi ,\mathbf{x})+\sum _{i=1}^N\mathrm{KL}[q_i(x_i)\Vert P_X(x_i)]-\log P(\mathbf{y}|\Phi )$$

$$\hat{x}_{i}^{(t+1)}=\eta\bigg{(}\frac{1}{\Delta}\sum_{\mu}\Phi_{\mu i}^{2},\frac{1}{\Delta}\sum_{\mu}\Phi_{\mu i}(y_{\mu}-{\textstyle\sum_{j\neq i}}\Phi_{\mu j}\hat{x}_{j}^{(t)})\bigg{)}\;\forall i,$$

$$\mathrm{KL}[P(\mathbf{x}|{\mathbf{\Phi }}_k,y_k)\Vert Q_k(\mathbf{x})].$$

Full text

1 Introduction

In current machine learning applications, one often faces the challenge of scale: massive data causes algorithms to explode in time and memory requirements. In such cases, when it is infeasible to process the full dataset simultaneously, one must resort to "online" or "streaming" methods which process only a small fraction of data points at a time — producing a step-by-step learning process. Such procedures are becoming more and more necessary to cope with massive datasets. For example, one can see the effectiveness of such approaches in deep learning via the stochastic gradient descent algorithm [1] or in statistical inference via the stochastic variational inference framework [2].

In this work, we treat streaming inference within a Bayesian framework where, as new data arrives, posterior beliefs are updated according to Bayes’ rule. One well known approach in this direction is assumed density filtering (ADF) [3, 4], which processes a single data point at a time, a procedure to which we refer to as fully online. A number of other works analyzed various related fully online algorithms [5, 6], especially in the statistical physics literature [7, 8, 9, 10, 11]. We are instead interested in the case where multiple samples – a mini-batch – arrive at once. Tuning the size of these mini-batches allows us to to explore the trade-off between the precision and efficiency.

Our motivation and setting are very much along the lines of streaming variational Bayes (VB) inference [12]. With respect to existing works, we bring three main contributions. (i) We introduce a streaming algorithm based on approximate message passing (AMP) [13, 14, 15] that we call Mini-AMP. As AMP treats some of the correlations which VB neglects, it is expected that AMP either outperforms or matches VB. (ii) Unlike other general-purpose algorithms for Bayesian inference, such as Gibbs sampling or VB, AMP possesses the state evolution method which asymptotically describes the performance of the algorithm for a class of generative models. We extend this state evolution analysis to Mini-AMP. (iii) For these generative models, we also analyze the optimal streaming procedure, within a class of procedures that retains only point-wise marginals from one step to another, and characterize regions of parameters where Mini-AMP reaches optimality.

2 Problem setting

Denoting the vector of $N$ values to be estimated by $\bm{x}$, the data presented at step $k$ by $\bm{y}^{(k)}$, and the collection of all previously presented data by $\mathcal{D}^{(k-1)}=\{\bm{y}^{(1)},\dots,\bm{y}^{(k-1)}\}$, the posterior distribution at step $k$ is given by

[TABLE]

In other words, with each presentation of new data, the prior distribution is updated with the posterior distribution derived from the previously presented data. Directly implementing this strategy is seldom feasible as the normalizing integral in (1) is intractable in general. Additionally, to keep the memory requirements small, we would like consider only the case where $O(N)$ parameters are passed from one step to the next. With this restriction, one cannot carry over high-order correlations from previous steps. Instead, following the strategy of [12], we resort to a factorized approximation of the "prior" term of (1),

[TABLE]

where $\mathcal{P}_{i}[\cdot]$ denotes the posterior marginals of parameter $x_{i}$ at a given step. Computing the marginals exactly is still computationally intractable for most models of interest. In the present work, this program is carried out with a scheme that is asymptotically exact for a class of generative models and that has the advantage of being amenable to a rigorous analysis.

We leverage the analysis of these models already conducted in the offline setting using concepts and techniques from statistical physics [16, 13, 14, 15, 17, 18, 19, 20] which have now been made almost entirely rigorous [21, 22, 23, 24, 25, 26]. We show in particular that – just as for the offline setting – phase transitions exist for mini-batch learning problems, and that their description provides information about the learning error that is achievable information-theoretically or computationally efficiently.

3 Generative models and offline learning

In our theoretical analysis, we consider inference in popular models with synthetic data generated from a given distribution, such as the perceptron with random patterns [27, 28], sparse linear regression with a random matrix (compressed sensing) [29] and clustering random mixtures of Gaussians [19]. For clarity, we restrict our presentation to the generalized linear models (GLMs), focusing on sparse linear estimation. Our results, however, can be extended straightforwardly to any problem where AMP can be applied. For offline GLMs, the joint distribution of the observation $\bm{y}\in\mathbb{R}^{M}$ and the unknown $\bm{x}\in\mathbb{R}^{N}$ is given by

[TABLE]

where $\bm{\Phi}_{\mu}$ is the $\mu$-th line of the $M\times N$ matrix $\Phi$. We consider the situation where $\Phi$ is a random matrix where each element is taken i.i.d. from ${\cal{N}}(0,1/N)$ and $\alpha=M/N$. Structured matrices have also been studied with AMP [30, 31]. Two situations of interest described by GLMs are (a) sparse linear regression (SLR) where the likelihood is Gaussian $P(y_{\mu}|z_{\mu})=\mathcal{N}(y_{\mu};z_{\mu},\Delta)$ and the parameters are sparse, for instance drawn from a Gauss-Bernoulli distribution $P_{X}(x_{i})=\rho\,\mathcal{N}(x_{i};0,1)+(1-\rho)\,\delta(x_{i})$, and (b) the probit regression problem P(y_{\mu}|z_{\mu})=\frac{1}{2}\operatorname{erfc}\big{(}-\frac{y_{\mu}z_{\mu}}{\sqrt{2\Delta}}\big{)} that reduces to the perceptron $P(y_{\mu}|z_{\mu})=\theta(y_{\mu}z_{\mu})$ when $\Delta\!\to\!0$ [32, 28, 15]. We first summarize the known relevant results for the fully offline learning problem, where one processes all data at once. Again, for clarity, we focus on the case of SLR.

The marginals estimated by AMP are given by [13, 14, 18]

[TABLE]

where $Z(A,B_{i})$ is a normalization factor. We shall denote the mean and variance of this distribution by $\eta(A,B)\equiv\frac{\partial}{\partial B}\log Z(A,B)$ and $\eta^{\prime}(A,B)\equiv\frac{\partial}{\partial B}\eta(A,B)$. The mean, in particular, provides an approximation to the minimum mean-squared error (MMSE) estimate of $\bm{x}$. The AMP iteration reads

[TABLE]

[TABLE]

One of the main strengths of AMP is that when the matrix $\Phi$ has i.i.d. elements, the ground truth parameters are generated i.i.d. from a distribution $P_{0}(x_{i})$, and $P(y_{\mu}|z_{\mu})={\cal N}(y_{\mu};z_{\mu},\Delta_{0})$, then the behavior and performance of AMP can be studied analytically in the large system limit ($N\rightarrow\infty$) using a technique called state evolution. This was proven by [21] who show that in the large $N$ limit, $A^{(t)}$ and $B_{i}^{(t)}$ converge in distribution such that, defining {\cal E}^{(t)}\equiv\mathbb{E}\,\big{(}\eta({\cal A}^{(t-1)},{\cal B}^{(t-1)})-x\big{)}^{2}, with $x\sim P_{0}(x)$, and ${\cal V}^{(t)}\equiv\mathbb{E}\,\eta^{\prime}({\cal A}^{(t-1)},{\cal B}^{(t-1)})$, one has

[TABLE]

The behavior of the algorithm is monitored by the computation of the scalar quantities ${\cal E}^{(t)}$ and ${\cal V}^{(t)}$.

The "Bayes-optimal" setting is defined as the case when the generative model is known and matches the terms in the posterior (1), i.e. when $P_{X}=P_{0}$, and $\Delta=\Delta_{0}$. One can show that in this case ${\cal E}^{(t)}={\cal V}^{(t)}$ (the so-called Nishimori property [15]), so that the state evolution further reduces to

[TABLE]

with $x\!\sim\!P_{X}(x)$, $z\!\sim\!\mathcal{N}(0,1)$, and ${\cal E}^{(t)}$ is the mean-squared error (MSE) achieved at iteration $t$.

Another set of recent results [24, 25] allows for the exact computation of the Bayes-optimal MMSE and the mutual information between the observations and the unknown parameters. Given model (3) with Gaussian likelihood, the mutual information per variable is given by the minimum of the so-called replica mutual information: $\lim_{N\to\infty}I({\bm{X}},{\bm{Y}})={\rm min}~{}i_{\text{RS}}({\cal E})$ where, defining, $\Sigma^{-2}({\cal E})\equiv\frac{\alpha}{\Delta+{\cal E}}$,

[TABLE]

with $x\sim P_{X}(x),\tilde{x}\sim P_{X}(\tilde{x})$ and $z\sim\mathcal{N}(0,1)$. The MMSE is then given by $\arg\min~{}i_{\text{RS}}({\cal E})$.

Comparisons between the MMSE and the MSE provided by AMP after convergence are very instructive, as shown in [15]. Typically, for large enough noise, ${\cal E}^{(t\to\infty)}={\cal E}_{\rm AMP}=\rm{MMSE}$ and AMP achieves the Bayes-optimal result in polynomial time, thus justifying, a posteriori, the interest of such algorithms in this setting. In fact, since the fixed points of the state evolution are all extrema of the mutual information (10), it is useful to think of AMP as an algorithm that attempts to minimize (10). However, a computational phase transition can exist at low noise levels, where $i_{\text{RS}}({\cal E})$ has more than a single minimum. In this case, it may happen that AMP does not reach the global minimum, and therefore ${\cal E}_{{\rm AMP}}>{\rm MMSE}$. It is a remarkable open problem to determine whether finding the MMSE in this region is computationally tractable. The results we have just described are not merely restricted to SLR, but appear mutatis mutandis in various cases of low-rank matrix and tensor factorization [17, 19, 20, 33, 23, 26] and also partly in GLMs [14, 15] (in GLMs the replica mutual information is so far only conjectured).

4 Main results

4.1 Mini-AMP

Our first contribution is the Mini-AMP algorithm, which adapts AMP to the streaming setting. Again, we shall restrict the presentation to the linear regression case. The adaptation to other AMP algorithms is straightforward. We consider a dataset of $M$ samples with $N$ features each, which we split into $B$ mini-batches, each containing $M_{b}=M/B$ samples. We denote $\alpha=M/N$ and $\alpha_{b}=M_{b}/N$. Crucially, for each step, the posterior marginal given by AMP (4) is the prior multiplied by a quadratic form. Performing the program discussed in (2) is thus tractable as the ${\cal P}_{i}$ are given by a Gaussian distribution multiplied by the original prior. The only modification w.r.t. the offline AMP at each step is thus to update the prior by multiplying the former one by the exponential in (4). In other words, we use the following "effective" prior when processing the $k$-th mini-batch:

[TABLE]

In practice, the only change when moving from AMP to Mini-AMP is therefore the update of the arguments of the $\eta$ function. After $k$ mini-batches have been processed, one replaces (7) by

[TABLE]

The corresponding pseudo-code is given as Algorithm $1$. Each Mini-AMP iteration has a computational complexity proportional to $M_{b}\times N$. We note that, in the fully online scheme when $M_{b}=1$, Mini-AMP with a single iteration performed per sample gives the same as ADF [3, 11].

4.2 State evolution

Theorem 1** (State evolution of Mini-AMP).**

For a random matrix $\Phi$, where each element is taken i.i.d. from ${\cal{N}}(0,1/N)$, the MSE of Mini-AMP can be monitored asymptotically ($N\!\to\!\infty$ while $\alpha_{b}\!=\!O(1)$) by iterating the following state evolution equations,

[TABLE]

where $x\sim P_{0}(x)$, and $z\sim{\cal N}(0,1)$. For each $k=1,\dots,N_{\text{b}}$, these equations are iterated from $t=1,\dots,t_{\text{max}}$, at which point we assign $\lambda_{k+1}=\lambda_{k+1}^{(t_{\text{max}})}$ and $\gamma_{k+1}=\gamma_{k+1}^{(t_{\text{max}})}$. The MSE given by AMP after the $k$-th mini-batch has been processed is given by ${\cal E}_{k}$. In the Bayes-optimal case, in particular, one can further show that the state evolution reduces to

[TABLE]

Proof sketch.

We apply the proof of state evolution for AMP in [21] to each mini-batch step, each with its own denoiser function $\eta(.)$ Each step is an instance of AMP with a new, independent matrix, and an effective denoiser given by (12). Using (8), the statistics of the denoisers are known, and the application of the standard AMP state evolution leads to (13). The Bayes-optimal case (14) then follows by induction, as in Sec. V.A.2 of [34]. ∎

Note that the above Theorem holds for any value of $t_{\rm max}$. Hence, even a stochastic version of the Mini-AMP algorithm, where for every mini-batch one only performs a few iterations without waiting for convergence in order to further speed up the algorithm, is analyzable using the above state evolution.

4.3 Optimal MMSE and mutual information under mini-batch setting

Theorem 2** (Mutual information for each mini-batch).**

For a random matrix $\Phi$, where each element is taken i.i.d. from ${\cal{N}}(0,1/N)$, in the Bayes-optimal setting, assume one has been given, after $k-1$ mini-batches, a noisy version ${\bm{r}}$ of unknown signal ${\bm{x}}$ with i.i.d noise ${\cal N}(0,\lambda^{-1})$. Given a new mini-batch with $\alpha_{b}\!=\!O(1)$, the mutual information per variable between the couple $({\bm{r}},{\bm{y}})$ and the unknown ${\bm{x}}$ is asymptotically given by $i={\rm min}~{}{i}^{b}_{\text{RS}}({\cal E}_{k})$ where, defining $\Sigma_{b}^{-2}(\lambda,{\cal E}_{k})\equiv\lambda+\frac{\alpha_{b}}{\Delta+{\cal E}_{k}}$,

[TABLE]

Proof sketch.

The proof is a slight generalization of the Guerra construction in [24]. Using properties of the Shannon entropy, the mutual information can be written as

[TABLE]

Computation of this expectation is simplified by noticing that it appears as equation (41) in the Guerra construction of [24], where it was used as a proof method for the offline result by interpolating from a pure noisy Gaussian channel (at "time" $t\!=\!0$) to the actual linear channel (at "time" $t\!=\!1$). Authors of [24] denoted $\lambda(t)$ as the variance of the Gaussian channel and $\gamma(t)$ as the variance of the linear channel. Our computation corresponds instead to a "time" $0\leq\tau\leq 1$ where both channels are used. Using Sec. V of [24], with the change of notation $\gamma(\tau)\!\to\!\Delta$ and $\lambda(\tau)\!\to\!\lambda$, we reach

[TABLE]

where $0\leq{\tau}\leq 1$ and $R_{{\cal E}_{b}}(t)$ a non-negative function called the reminder. The validity of the mutual information formula in the offline situation [24, 25] implies that the integral of the reminder in $[0,1]$ is zero when ${\cal E}^{*}=\rm{argmin}\,~{}{i}^{b}_{\text{RS}}({\cal E}_{b})$. Since $R(t)$ is non-negative, this implies that it is zero almost everywhere, thus $\int_{0}^{\tau}R_{{\cal E}^{*}_{b}}(t)~{}~{}{\rm d}t=0$. ∎

Theorem 3** (MMSE for each mini-batch).**

With the same hypothesis of Theorem 2, the MMSE when one has access to a noisy estimate with i.i.d. noise ${\cal N}(0,\lambda^{-1})$ and the data from the mini-batch at step $k$, is

[TABLE]

Proof sketch.

The proof follows again directly from generic results on the Guerra interpolation in [24] and the so-called I-MMSE formula ${\rm d}i(\Delta)/{\rm d}\Delta^{-1}=\alpha_{b}y_{\rm MMSE}/2$ and y-MMSE formula $y_{\rm MMSE}={\rm MMSE}/(1+\Delta\,{\rm MMSE})$ linking the mutual information and the MMSE. ∎

One can show through explicit computation that the extrema of $i^{b}_{\rm RS}({\cal E})$ correspond – just as in the offline case – to the fixed points of the state evolution.

Using these results, we can analyze, both algorithmically and information theoretically, the mini-batch program (2). Indeed, the new information on the parameters $\bm{x}$ passed from mini-batch $k-1$ to $k$ contained in (11) is simply a (Gaussian) noisy version of $\bm{x}$ with inverse variance $\lambda_{k-1}$. This is true for the AMP estimate (see (11)) and, in the large $N$ limit, for the exact marginalized posterior distribution as well (see e.g. [35, 36, 26]). The optimal MSE at each mini-batch is thus given by the recursive application of Theorem 3, where at each mini-batch $k=1,\ldots,B$ we minimize (15) using $\lambda_{k-1}=\lambda_{k-2}+{\alpha_{b}}/{(\Delta+\rm{MMSE}_{k-1})}$.

Now we can compare the MSE reached by the Mini-AMP algorithms to the MMSE. If, for each mini-batch, the MMSE is reached by the state evolution of the Mini-AMP algorithm starting from the previously reached MSE, then we have the remarkable result that Mini-AMP performs a Bayes-optimal and computationally efficient implementation of the mini-batch program (2). Otherwise, the Mini-AMP is suboptimal. It remains an open question whether in that case any polynomial algorithm can improve upon the MSE reached by Mini-AMP.

All our results can be directly generalized to the case of AMP for matrix or tensor factorization, as derived and proven in [37, 26, 20]. They can also be adapted to the case of GLMs with non-linear output channels. However, in this setting, the formula for the mutual information has not been yet proven rigorously.

5 Performance and phase transitions on GLMs

5.1 Optimality & efficiency trade-offs with Mini-AMP

We now illustrate the above results on some examples. In Figure 1, we consider the SLR model and the perceptron with binary $\pm 1$ parameters, both with random matrices $\Phi\in\mathbb{R}^{M\times N}$, $\Phi_{\mu i}\sim\mathcal{N}(0,1/N)$. Our analysis quantifies the loss coming from using mini-batches with respect to a fully offline implementation. In the limit of small mini-batch $\alpha_{\text{b}}\to 0$, we recover the results of the ADF algorithm which performs fully online learning, processing one sample at a time [9, 11]. This suggests that the state evolution accurately describes the behavior of Mini-AMP beyond the theoretical assumption of $\alpha_{b}\!=\!O(1)$, even for mini-batches as small as a single sample.

The effect of the mini-batch sizes varies greatly with the problem. For the perceptron with $\pm 1$ weights, a zero error is eventually obtained after a sufficient number of mini-batches have been processed. Moreover, the dependence on the mini-batch size is mild: while the offline scheme achieves zero error at $\alpha\approx 1.5$ [27, 28], the fully online does it at $\alpha\approx 4.4$ [9], that is, going from offline to a fully online scheme costs only about three times more data points. The behavior of the Mini-AMP for SLR shows instead rather drastic changes with the mini-batch size. The MSE decays smoothly when the mini-batch size is small. However, as we increase it, a sudden decay occurs after a few mini-batches have been processed. For the noiseless case ($\Delta=0$), the study of the state evolution shows that the asymptotic (in $\alpha$) MSE is given by

[TABLE]

if $\alpha_{\text{b}}\leq\rho$, and by 0 otherwise. These results provide a basis for an optimal choice of mini-batch size. Given the drastic change in behavior past a certain mini-batch size, one concludes that small investments in memory might be worthwhile, since they can lead to large gains in performance.

Finally, we have compared the Mini-AMP scheme with the streaming VB approach [12] using the mean-field algorithm described in [38] for SLR. While the mean-field approach is found to give results comparable to AMP in the offline case in [38], we see here that the results are considerably worse in the streaming problem. In fact, as shown in Figure 1 (center), mean-field can give worse performance than the fully online ADF, even when processing rather large mini-batches.

5.2 Phase transitions

It turns out that, just as for the offline setting, there are phase transitions appearing for mini-batch learning, in terms of the learning error that is achievable information-theoretically (MMSE) or computationally efficiently (by Mini-AMP). These can be understood by an analysis of the function $i^{b}_{\rm RS}$, since the minimum of $i^{b}_{\rm RS}$ gives the MMSE, and since AMP is effectively trying to minimize $i^{b}_{\rm RS}$ starting from the MSE reached at the previous mini-batch steps.

Let us illustrate the reason behind the sharp phenomenon in the behavior of AMP in Fig.1. We show, in Fig. 2 (left), an example of the function $i^{b}_{\rm RS}({\cal E})$ for the streaming SLR problem as a function of the MSE ${\cal E}$ as each mini-batch is being processed. Initially, it presents a “good” and a “bad” minimum, at small and large MSEs respectively. In the very first batch, AMP reaches the bad minimum. As more batches are processed, the good minimum becomes global, but AMP is yet not able to reach it, and keeps returning the bad one instead. This indicates a computational phase transition, and we expect that other algorithms will, as AMP, fail to deliver the MMSE in polynomial time when this happens. Eventually, the good minimum becomes unique, at which point AMP is able to reach it, thus yielding the sudden decay observed in Figure 1.

Consider now the Bayes-optimal "streaming-MMSE" given by the global minimum of the mutual information at each step, regardless of whether AMP achieves it. In the offline noiseless case, the MMSE is achieved by AMP only if the processed batch has size $\alpha\geq\alpha_{\rm offline}$ or $\alpha\leq\rho$ [18]. In the streaming case, we also observe that Mini-AMP reaches the streaming-MMSE if the mini-batch size is sufficiently small or sufficiently large. In Figure 2 (right) we compare MMSE to the MSE reached by Mini-AMP, with a region between the full and dashed line where the algorithm is sub-optimal.

6 Mini-AMP for matrix factorization problems

We now consider the case of low-rank matrix factorization, and in particular clustering using the Gaussian mixture model (GMM) with $R$ clusters. For such problems, the generative model reads

[TABLE]

where $\bm{U}_{i}$ and $\bm{V}_{j}$ give the $i$-th row of $U$ and $j$-th row of $V$ respectively. Each of the $R$ columns of $U$ describe the mean of a $N$-variate i.i.d. Gaussian, and $V$ has the role of picking one of these Gaussians. Finally, each column of $Y$ is given by the chosen column of $U$ plus Gaussian noise. In clustering, these are the data points, and the objective is to figure out the position of the centroids as well as the label assignment, given by the columns of $U$ and the rows of $V$ respectively. In the streaming setting the columns of the matrix $Y$ are arriving in mini-batches. The offline AMP algorithm, its state evolution, and corresponding proofs are known for matrix factorization from [17, 40, 19, 37, 23, 26]. The Mini-AMP is obtained by adjusting the update of the estimators using (12).

In GMM clustering with prior on U having zero mean, there is an interesting "undetectability" phase transition for $R\leq 4$. If the number of samples is such that $\alpha=M/N<\alpha_{c}=R^{2}\Delta^{2}$, then the Bayes optimal posterior asymptotically does not contain any information about the ground truth parameters [19]. This transition survives even when $R>4$, in the sense that AMP and other tractable algorithms are unable to find any information on the ground truth parameters.

In the streaming problem, this undetectability implies that for mini-batches of relative size $\alpha_{\text{b}}<\alpha_{\text{c}}$, Mini-AMP does not improve the error of the random estimator, no matter the number of mini-batches presented. In particular, the fully online algorithm does not provide any useful output in this scenario. If $\alpha_{\text{b}}>\alpha_{\text{c}}$, on the other hand, an accurate reconstruction of the unknown values becomes possible. We illustrate the MSE as a function of the mini-batch size in Figure 3.

While we have presented Mini-AMP as a means for a theoretical analysis, it can be applied to real data, performing concrete learning tasks. To illustrate its efficacy, we have considered the classical problem of unsupervised clustering using the GMM. In Figure 3, Mini-AMP is shown to obtain better performance for real data clustering than mini-batch K-means, a state-of-the-art algorithm for streaming clustering [39].

7 Conclusion

Let us conclude by stating that the Mini-AMP algorithm can be applied to any problem for which the streaming can be defined and for which offline AMP exists. Therefore, we expect that this novel development will improve the usefulness of AMP algorithms in more practical situations.

Acknowledgments

This work has been supported by the ERC under the European Union’s FP7 Grant Agreement 307087-SPARCS. AM thanks Paulo V. Rossi and Thibault Lesieur for insightful discussions.

Appendix A AMP equations for different classes of models

We present here the AMP equations for different models. As before, adapting them to the streaming setting is done by introducing $\bm{\Lambda}$, $\bm{\Theta}$ variables and replacing the $\eta(A,B)$ function with (12).

A.1 Generalized linear models

Denote by $\bm{y}\in\mathbb{R}^{M}$ the response variable, by $\Phi\in\mathbb{R}^{M\times N}$ the design matrix, and by $\bm{x}\in\mathbb{R}^{N}$ the parameter vector that we want to estimate. Then our generative model reads

[TABLE]

The GAMP algorithm provides the following approximation to the marginals of $\bm{x}$

[TABLE]

and, to the marginals of $\bm{z}$

[TABLE]

The parameters $\bm{A}$, $\bm{B}$, $\bm{\omega}$ and $\bm{V}$ are determined by iterating the GAMP equations. We denote the mean and variance of $q_{x}(A,B)$ by $\eta(A,B)=\frac{\partial}{\partial B}\log Z_{x}(A,B)$ and $\eta^{\prime}(A,B)=\frac{\partial\eta}{\partial B}(A,B)$ respectively; moreover, we define $g_{\text{out}}(y,\omega,V)=\frac{\partial}{\partial\omega}\log Z_{z}(y,\omega,V)$. The GAMP equations then read [14, 15]

[TABLE]

Note that, for a Gaussian likelihood $P(y|z)=\mathcal{N}(y;z,\Delta)$, $g_{\text{out}}(y,\omega,V)=\frac{y-\omega}{\Delta+V}$. Then, by defining $\bm{\tilde{z}}^{(t)}=\bm{y}-\bm{\omega}^{(t)}$ and replacing $\bm{V}^{(t)}$ and $\bm{A}^{(t)}$ by its averages, we get via the central limit theorem

[TABLE]

where we have assumed the $\Phi_{\mu i}$ are i.i.d. and have zero mean and variance $1/N$. From these (5)-(7) follow through.

A.1.1 Variational Bayes

We compare AMP equations to the Variational Bayes (VB) ones, which we use with the Streaming Variational Bayes scheme. For simplicity we restrict ourselves to the Gaussian case. As usual, VB is derived by determining the $q_{i}(x_{i})$ which minimize

[TABLE]

If done by means of a fixed-point iteration, this minimization leads to [38]

[TABLE]

where $\eta$ is defined as before. A closer inspection shows that the same equations are obtained by setting $\bm{V}^{(t)}=0$ in (23).

As shown in [38], this iteration leads to good results when performed sequentially, if the noise $\Delta$ is not fixed but learned. We employ this same strategy in our experiments.

A.1.2 Assumed Density Filtering

The assumed density filtering (ADF) algorithm [3, 4] replaces the posterior at each step $k$, $P(\bm{x}|\bm{\Phi}_{k},y_{k})\propto P(y_{k}|\bm{\Phi}_{k},\bm{x})\,Q_{k-1}(\bm{x})$, by the distribution $Q_{k}(\bm{x})$ that minimizes

[TABLE]

Note this is the direct KL divergence, and not the reverse one $\operatorname{KL}[Q(\bm{x})\|P(\bm{x}|\Phi,\bm{y})]$ that is minimized in Variational Bayes. In particular if $Q_{k}(\bm{x})=\prod_{i=1}^{N}q_{k}(x_{i})$, then minimizing this KL divergence leads to the following integral

[TABLE]

which is tractable since we are processing a single sample $y_{k}$ at a time, i.e. since the likelihood consists of a single factor. These are actually the exact marginals of $P(\bm{x}|\bm{\Phi}_{k},y_{k})$, and also the equations given by the belief propagation (BP) algorithm in the single sample limit.

The ADF equations for GLMs can be derived by using (28) together with the central limit theorem [11]. Because BP gives ADF in the single sample limit, AMP (which is based on BP) gives the equations derived by [11] when $M=1$, if one additionally neglects the correction term $-\bm{V}^{(t)}\bm{g}^{(t-1)}$ (analogously, if a single iteration is performed).

A.2 Low-rank matrix factorization

Denote by $Y\in\mathbb{R}^{N\times M}$ the matrix we want to factorize, and by $U\in\mathbb{R}^{N\times R}$, $V\in\mathbb{R}^{M\times R}$ the matrices which product approximates $Y$. The generative model then reads

[TABLE]

where $\bm{U}_{i}$ and $\bm{V}_{j}$ denote the $i$-th and $j$-th rows of $U$ and $V$ respectively. The algorithm provides the following approximation to the marginal of $\bm{U}_{i}$

[TABLE]

and $q_{V}$ is analogously defined as the marginal of $\bm{V}_{j}$. As in the previous case, $A_{U}\in\mathbb{R}^{R\times R}$ and $\bm{B}_{U,i}\in\mathbb{R}^{R}$ are to be determined by iterating a set of equations. We denote by $B_{U}$ the $N\times R$ matrix which rows are given by $\bm{B}_{U,i}$, $i=1,\dots,N$. The functions $\eta_{U}(A,\bm{B})=\mathbb{E}_{q_{U}}\,\bm{U}=\nabla_{\bm{B}}\log Z_{U}(A,\bm{B})$ and $\eta^{\prime}_{U}(A,\bm{B})=\nabla_{\bm{B}}\eta_{U}(A,\bm{B})$ give the mean and covariance of $q_{U}$, and $\eta_{V}$ and $\eta^{\prime}_{V}$, defined analogously, the mean and covariance of $q_{V}$.

In order to write the AMP equations, we first introduce

[TABLE]

so as to define an effective Gaussian channel [41]. The equations to be iterated are then

[TABLE]

In order to adapt this algorithm to the online setting, we repeat procedure (12) and, as the $k$-th batch is processed, replace calls to $\eta_{U}(A_{U},\bm{B}_{U})$ by

[TABLE]

We assume that $U$ is fixed and $V_{k}$ changes for each batch $k$; thus, the calls to $\eta_{V}$ do not change.

Appendix B State evolution and asymptotic limits

Through the state evolution equations, we analyze the behaviour and performance of the algorithms described in the previous section. We restrict ourselves to the Bayes-optimal case (i.e. the Nishimori line), where the generative model is known. The strategy we use to go from the offline to the streaming setting is easily adapted to the non-optimal case.

B.1 Generalized linear models

The state evolution equations for a GLM with likelihood $P(y|z)$ and prior $P_{X}(x)$ are

[TABLE]

where we denote $\rho=\mathbb{E}x^{2}$, and the averages are taken with respect to $P(x,b)=P_{X}(x)\mathcal{N}(b;\hat{m}x,\hat{m})$ and $P(y,\omega,z)=P(y|z)\mathcal{N}(z;\omega,\rho-m)\mathcal{N}(\omega;0,m)$. The MSE at each step is obtained from ${\cal E}^{(t)}=\rho-m^{(t)}$. For a Gaussian likelihood $P(y|z)=\mathcal{N}(y;z,\Delta)$, $\hat{m}^{(t)}=\frac{\alpha}{\Delta+{\cal E}^{(t)}}$ and we recover (9).

The fixed points of the state evolution extremize the so-called replica free energy

[TABLE]

which gives the large system limit of the Bethe free energy (extremized by AMP). The mutual information (10) differs from $\phi(m)=\operatorname*{extr}_{\hat{m}}\,\phi(m,\hat{m})$ by a constant – more specifically, by the entropy of $P(y|z)$, $i_{\text{RS}}(m)=\phi(m)-\alpha H[P(y|z)]$. For a Gaussian likelihood, $H[P(y|z)]=\frac{1}{2}\log(2\pi e\Delta)$.

In order to adapt this to the streaming setting, we introduce $\lambda_{k}^{(t)}=\sum_{\ell=1}^{k-1}\hat{m}_{\ell}^{(t_{\text{max}})}+\hat{m}_{k}^{(t)}$ and iterate instead, for each mini-batch $k$

[TABLE]

with the averages now computed over $P(x,b_{k})=P_{X}(x)\mathcal{N}(b_{k};\lambda_{k}x,\lambda_{k})$. These equations should be iterated for $t=1,\dots,t_{\text{max}}$, at which point we assign $\lambda_{k}=\lambda_{k-1}^{(t_{\text{max}})}$. The MSE on $\bm{x}$ after mini-batch $k$ is processed is then given by ${\cal E}_{k}=\rho-m_{k}^{(t_{\text{max}})}$.

Note that in the small batch size limit ($\alpha_{b}\to 0$), the equation for $\lambda$ becomes an ODE

[TABLE]

which describes the performance of the ADF algorithm [3, 11].

The free energy is also easily rewritten

[TABLE]

or analogously, by working with $\lambda_{k}=\lambda_{k-1}+\hat{m}_{k}$ instead of $\hat{m}_{k}$

[TABLE]

from which it is clear that the extrema of $\phi_{k}$ are given by the fixed points of (36).

B.1.1 Asymptotic behavior

Equations (36) can be put in the following form

[TABLE]

where $\varepsilon(\lambda_{k})$ and $\delta({\cal E}_{k})$ are functions that depend on the prior/channel respectively. Assuming $\varepsilon$ is invertible, we rewrite this system of equations as a function of ${\cal E}_{k}$ only

[TABLE]

and then solve this equation for ${\cal E}_{k}$; that gives us a recurrence relation which is unsolvable in most cases. We use instead asymptotic forms for $\varepsilon$ and $\delta$, obtained in the $\lambda\to\infty$, ${\cal E}\to 0$ limit. For the Bernoulli-Gaussian prior $P_{0}(x_{i})=\rho\,\mathcal{N}(x_{i};0,1)+(1-\rho)\,\delta(x_{i})$, we have

[TABLE]

while a Gaussian likelihood gives, in the $\Delta\to 0$ limit, $\delta({\cal E})=\frac{1}{{\cal E}}$. Thus, for SLR

[TABLE]

leading to (18). In the $\alpha_{b}\to 0$ limit we recover the expression obtained by [11], $\operatorname{MSE}(\alpha)\sim e^{-\frac{\alpha}{\rho}}$.

B.2 Low-rank matrix factorization

Also for low-rank models the large $N$ limit can be analyzed by taking into account that $A_{U,V}^{(t)}$ and $B_{U,V}^{(t)}$ converge in distribution to [41]

[TABLE]

where $M_{U},M_{V}\in\mathbb{R}^{R\times R}$ are the overlap matrices between the ground truth and the estimate at time $t$, that is

[TABLE]

While computing the expectations might become unfeasible for $R>1$, an ansatz of the following form can often be used [41, 19]

[TABLE]

with $\mathbb{J}_{R}$ denoting the $R\times R$ matrix of ones. This significantly simplifies the iteration above.

Again we adapt this to the online case by incrementing the matrices obtained as each batch is processed, that is, we replace $M_{V}$ by

[TABLE]

Note that since $U$ is fixed, the equation for $M_{U}$ does not change.

The replica free energy reads, in the offline case

[TABLE]

and we adapt it to the online case by taking into account that $\lambda_{V}$ is being incremented

[TABLE]

Appendix C Performance for different number of iterations

For our experiments in Figure 1, Mini-AMP has been iterated until, for each block, convergence is achieved – that is, until $\frac{1}{N}\|\hat{\bm{x}}^{(t)}-\hat{\bm{x}}^{(t-1)}\|_{1}<10^{-13}$. Remarkably, our framework allows us to study the performance of the algorithm even if we do not iterate it until convergence, but only for a few steps $t_{\text{max}}$ instead. In Figure 4, we investigate the performance of Mini-AMP under the same settings of Figure 1 (center), for different values of $t_{\text{max}}$. We observe that the performance deteriorates if convergence is not reached.

Appendix D Experiments on real-world data

For the experiments with real data, we have used the following model

[TABLE]

where

[TABLE]

is a truncated normal distribution supported on the positive quadrant of a $R$-dimensional space, and $\mathcal{Z}(\bm{\mu},\sigma^{2})$ ensures proper normalization.

Note that evaluating the $\eta(A,\bm{B})$ function in this case is not trivial, since it depends on the following integral

[TABLE]

where, in this case, $P_{0}(x_{k})=\mathcal{N}(x_{k};\mu_{k},\sigma^{2})\,\theta(x_{k})$. We proceed by performing a mean-field approximation. We first define

[TABLE]

for scalar $A$ and $B$. Then, for each $i=1,\dots,N$, we iterate

[TABLE]

sequentially in $k=1,\dots,R$ until convergence is reached, at which point we use the values obtained for assigning $\bm{\hat{U}}^{(t)}$ in AMP. The variances are computed from

[TABLE]

and the covariance matrix used in AMP is obtained as a function of these variances

[TABLE]

where in order to assign the off-diagonal terms we have used a linear response approximation, $\sigma_{k\ell}=\sum_{i=1}^{N}\frac{\partial\hat{U}_{ik}}{\partial B_{i\ell}}$.

We proceed by detailing other aspects of the experiments

Initialization

At the first few mini-batches (usually the first five), we reinitialize the position of the centroids. We use the same strategy as the k-means++ algorithm [42]: the first centroid is picked at random from the data points, and the next ones are sampled so as to have them far apart from each other. The labels are initialized according to the closest centroid.

Stopping criterion

For each batch, Mini-AMP was iterated either for 50 steps or until $\frac{1}{NR}\sum_{ik}|\hat{U}_{ik}^{(t)}-\hat{U}_{ik}^{(t-1)}|+\frac{1}{MR}\sum_{jk}|\hat{V}_{jk}^{(t)}-\hat{V}_{jk}^{(t-1)}|<10^{-7}$.

Noise learning

We do not assign a fixed value for $\Delta$, but instead update it after each mini-batch is processed using a simple learning rule

[TABLE]

Preprocessing

For MNIST, we work with all samples of digits 0, 1 and 2. They are rescaled so that the pixel intensities are between 0 and 1. For the 20 newsgroups dataset, we build Term Frequency Inverse Document Frequency (TF-IDF) features for 3 top-level hierarchies (comp, rec and sci), and use the 1000 most frequent words; we rescale each feature vector so that its maximum is equal to 1.

Mini-batch K-means

We use the mini-batch K-means [39] implementation available on scikit-learn [43]. Default parameters were used, apart from the centroids initialization, which was set to random normal variables of zero mean and variance $10^{-3}$ – this seemed to improve the algorithm performance with respect to the standard choices.