Streaming Bayesian inference: theoretical limits and
mini-batch approximate message-passing
Andre Manoel1 , Florent Krzakala2 , Eric W. Tramel3 , and Lenka Zdeborová4
1

arXiv:1706.00705v1 [stat.ML] 2 Jun 2017

2

Neurospin, CEA, Université Paris-Saclay
LPS ENS, CNRS, PSL, UPMC & Sorbonne Univ.
3
OWKIN
4
IPhT, CNRS, CEA, Université Paris-Saclay
June 5, 2017
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.

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–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–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
1

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 x, the data presented at step k by y (k) , and the collection
of all previously presented data by D(k−1) = {y (1) , . . . , y (k−1) }, the posterior distribution at step k is given
by
P (y (k) |x) P (x|D(k−1) )
P (x|y (k) , D(k−1) ) = R
.
(1)
dx P (y (k) |x) P (x|D(k−1) )
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),
P (x|D(k−1) ) ≈ Q(k−1) (x) =

N
Y


Pi P (x|y (k−1) , D(k−2) ) ,

(2)

i=1

where Pi [·] denotes the posterior marginals of parameter xi 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 [13–20] which have now been made almost entirely rigorous [21–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 y ∈ RM and the unknown x ∈ RN is given by
P (y, x|Φ) =

M
Y

P (yµ |zµ ≡ Φµ · x)

µ=1

N
Y

PX (xi ) .

(3)

i=1

where Φµ is the µ-th line of the M × N matrix Φ. We consider the situation where Φ is a random matrix
where each element is taken i.i.d. from N (0, 1/N ) and α = 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µ |zµ ) = N (yµ ; zµ , ∆) and the parameters are sparse, for instance drawn from a
Gauss-Bernoulli distribution
PX (xi ) = ρ N (xi ; 0, 1) + (1 − ρ) δ(xi ), and (b) the probit regression problem
y z 
P (yµ |zµ ) = 12 erfc − √µ2∆µ that reduces to the perceptron P (yµ |zµ ) = θ(yµ zµ ) when ∆ → 0 [15, 28, 32]. We
2

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]

2
1
P (xi |Φ, y) ≈ q(xi |A, Bi ) = PX (xi ) e− 2 Axi +Bi xi Z(A, Bi ) ,
(4)
where Z(A, Bi ) is a normalization factor. We shall denote the mean and variance of this distribution by
∂
∂
η(A, B) ≡ ∂B
log Z(A, B) and η 0 (A, B) ≡ ∂B
η(A, B). The mean, in particular, provides an approximation
to the minimum mean-squared error (MMSE) estimate of x. The AMP iteration reads
z (t) = y − Φx̂(t) + α−1 z (t−1) A(t−1) V (t) ,
B (t) = A(t) x̂(t) + A(t) α−1 ΦT z (t) ,
(t+1)

x̂i

A(t) =

(t)

= η(A(t) , Bi ) ∀i,

V (t+1)

α
,
∆ + V (t)
PN
(t)
= N1 i=1 η 0 (A(t) , Bi ).

(5)
(6)
(7)

One of the main strengths of AMP is that when the matrix Φ has i.i.d. elements, the ground truth
parameters are generated i.i.d. from a distribution P0 (xi ), and P (yµ |zµ ) = N (yµ ; zµ , ∆0 ), then the behavior and performance of AMP can be studied analytically in the large system limit (N → ∞) using a
technique called state evolution. This was proven by [21] who show that in the large N limit, A(t) and
2
(t)
Bi converge in distribution such that, defining E (t) ≡ E η(A(t−1) , B (t−1) ) − x , with x ∼ P0 (x), and
V (t) ≡ E η 0 (A(t−1) , B (t−1) ), one has


α
∆0 + E(t)
(t)
(t)
(t)
(t)
(t)
.
(8)
A
A =
,
Bi
B ∼ N A x, α
(∆ + V(t))2
∆ + V (t)
The behavior of the algorithm is monitored by the computation of the scalar quantities E (t) and 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 PX = P0 , and ∆ = ∆0 . One can show that in this case E (t) = V (t) (the
so-called Nishimori property [15]), so that the state evolution further reduces to
p
α
(t)
0
(t−1)
(t−1)
,
E
=
Eη
(A
,
A
x
+
A(t) =
A(t−1) z),
(9)
∆ + E (t)
with x ∼ PX (x), z ∼ N (0, 1), and 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: limN →∞ I(X, Y ) = min iRS (E) where, defining, Σ−2 (E) ≡ ∆+E





(x̃−(x+zΣ(E)))2
α
E
E
1
−
2 (E)
2Σ
iRS (E) =
+ log 1 +
− Ex,z log Ex̃ e
− ,
(10)
2 ∆+E
∆
2
with x ∼ PX (x), x̃ ∼ PX (x̃) and z ∼ N (0, 1). The MMSE is then given by arg min iRS (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, E (t→∞) = EAMP = 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 iRS (E) has more than a single minimum. In this case, it may
happen that AMP does not reach the global minimum, and therefore EAMP > 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, 23, 26, 33] and also partly in GLMs [14, 15] (in GLMs
the replica mutual information is so far only conjectured).
3

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 Mb = M/B samples. We denote α = M/N and αb = Mb /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 Pi 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:
PΛkk−1 ,Θk−1 (x) = PX (x)

N
Y

2

1

e− 2 Λk−1 xi +Θk−1,i xi ,

where Λk−1 =

i=1

k−1
X

Al , Θk−1,i =

`=1

k−1
X

Bl,i .

(11)

`=1

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


Pk−1
(t) Pk−1
(t)
(t+1)
+A
,
+B
x̂k,i = η
A
B
`
`,i
k
k,i ,
{z }
{z }
| `=1
| `=1
Θk−1,i

Λk−1

(t+1)
Vk

=

1
N

PN

i=1

η

0




Pk−1
(t)
(t) Pk−1
A` +Ak , `=1 B`,i +Bk,i .
| `=1
{z }
| {z }

(12)

Θk−1,i

Λk−1

The corresponding pseudo-code is given as Algorithm 1. Each Mini-AMP iteration has a computational
complexity proportional to Mb × N . We note that, in the fully online scheme when Mb = 1, Mini-AMP with
a single iteration performed per sample gives the same as ADF [3, 11].
Algorithm 1 Mini-AMP
1: initialize Λ0 = 0, Θ0,i = 0 ∀i
2: for k ← 1 to B do
(1)
3:
initialize zk,µ = 0 ∀µ
4:
5:
6:
7:
8:
9:
10:
11:
12:

(1)

(1)

initialize x̂k,i = η(Λk−1 , Θk−1,i ) ∀i, Vk
for t ← 1 to tmax do
(t)
compute zk using (5)
(t)
(t)
compute Ak , Bk using (6)
(t+1)
(t+1)
compute Vk
, x̂k
using (12)
end for
accumulate Λk ← Λk−1 + Ak
accumulate Θk ← Θk−1 + Bk
end for

4.2

=

1
N

PN

i=1

η 0 (Λk−1 , Θk−1,i )

State evolution

Theorem 1 (State evolution of Mini-AMP). For a random matrix Φ, where each element is taken i.i.d. from
N (0, 1/N ), the MSE of Mini-AMP can be monitored asymptotically (N → ∞ while αb = O(1)) by iterating the

4

following state evolution equations,
(t)
λk

= λk−1 +

αb
(t)

∆ + Vk

,

(t+1)
Vk

= Ex,z η

(t+1)
Ek



(t)

(t)
γk

= γk−1 + αb

∆0 + Ek

(t)

(∆ + Vk )2

,

0

(t)
(t)
λk , λk x

= Ex,z η

(t)
(t)
λk , λk x

q
(t) 
+ γk z ,
q
+

(t) 
γk z

−x

2

(13)
,

where x ∼ P0 (x), and z ∼ N (0, 1). For each k = 1, . . . , Nb , these equations are iterated from t = 1, . . . , tmax ,
(tmax )
(tmax )
at which point we assign λk+1 = λk+1
and γk+1 = γk+1
. The MSE given by AMP after the k-th mini-batch
has been processed is given by Ek . In the Bayes-optimal case, in particular, one can further show that the
state evolution reduces to
q

2
αb
(t+1)
(t)
(t) 
(t)
(t)
λk = λk−1 +
,
E
=
E
η
λ
,
λ
x
+
λ
z
−
x
.
x,z
(14)
k
k
k
k
(t)
∆ + Ek
Proof sketch. We apply the proof of state evolution for AMP in [21] to each mini-batch step, each with
its own denoiser function η(.) 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 tmax . 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 Φ, where each element is taken
i.i.d. from N (0, 1/N ), in the Bayes-optimal setting, assume one has been given, after k − 1 mini-batches,
a noisy version r of unknown signal x with i.i.d noise N (0, λ−1 ). Given a new mini-batch with αb = O(1),
the mutual information per variable between the couple (r, y) and the unknown x is asymptotically given by
αb
i = min ibRS (Ek ) where, defining Σ−2
b (λ, Ek ) ≡ λ + ∆+Ek ,
ibRS (Ek )




(x̃−(x+zΣb (λ,Ek )))2
−
αb
Ek
Ek
1 + αb
2Σ2 (λ,Ek )
b
=
− Ex,z log Ex̃ e
.
+ log 1 +
−
2 ∆ + Ek
∆
2

(15)

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
ky−Φxk2
2
2∆

kr−xk2
2
2λ−1

1 + αb
.
2
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
λ(t) as the variance of the Gaussian channel and γ(t) as the variance of the linear channel. Our computation
corresponds instead to a "time" 0 ≤ τ ≤ 1 where both channels are used. Using Sec. V of [24], with the
change of notation γ(τ ) → ∆ and λ(τ ) → λ, we reach
Z τ
I(Y, R; X)
lim
= ibRS (Ek ) −
REk (t) dt + O(1),
(16)
N →∞
N
0
I(Y, R; X) = H(Y, R) − H(Y |X) − H(R|X) = Ey,r log Ex e−

e−

−

where 0 ≤ τ ≤ 1 and REb (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
E ∗ = argmin ibRS (Eb ). Since R(t) is non-negative, this implies that it is zero almost everywhere, thus
Rτ
REb∗ (t) dt = 0.
0
5

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 N (0, λ−1 ) and the data from the mini-batch at step k, is
MMSE = argmin ibRS (Ek ) .

(17)

Proof sketch. The proof follows again directly from generic results on the Guerra interpolation in [24] and the
so-called I-MMSE formula di(∆)/d∆−1 = αb yMMSE /2 and y-MMSE formula yMMSE = MMSE/(1 + ∆ MMSE)
linking the mutual information and the MMSE.
One can show through explicit computation that the extrema of ibRS (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 x passed from mini-batch k − 1 to k contained
in (11) is simply a (Gaussian) noisy version of x with inverse variance λ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.
[26, 35, 36]). The optimal MSE at each mini-batch is thus given by the recursive application of Theorem 3,
where at each mini-batch k = 1, . . . , B we minimize (15) using λk−1 = λk−2 + αb /(∆ + MMSEk−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 [20, 26, 37]. 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
5.1

Performance and phase transitions on GLMs
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 ±1 parameters, both with random matrices Φ ∈ RM ×N , Φµi ∼ 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 αb → 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 α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 ±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 α ≈ 1.5 [27, 28], the
fully online does it at α ≈ 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 (∆ = 0), the
study of the state evolution shows that the asymptotic (in α) MSE is given by
MSEx (α) ∼ e

− α1 log(1−
b

αb
ρ

)α

,

(18)

if αb ≤ ρ, 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.

6

10−11
−14

10−17
10−20

10−2

fully online
αb = 0.25
αb = 0.35

0.5

offline

0.4

10−4
10−6

0.3
0.2
0.1

10−8
0 0.49 2
4
6
8
10
total amount of data α = M/N

fully online
αb = 0.01
αb = 0.625
αb = 1.25
offline

offline

10−8

10

100

MSE(x)

MSE(x)

10−5

fully online
αb = 0.05
αb = 0.25
αb = 0.3
αb = 0.35

error(x)

10

−2

0
1
2
3
4
5
total amount of data α = M/N

0.0

0
1
2
3
4
5
total amount of data α = M/N

Figure 1: Accuracy of Mini-AMP inference as a function of the total amount of presented data. Left: SLR with
∆ = 0 and sparsity ρ = 0.3. Center: SLR with ∆ = 10−8 and ρ = 0.3. Right: Perceptron learning with Rademacher
(±1) parameters/synapses. For different mini-batch sizes (colors), we show both the state evolution predictions
for Mini-AMP (solid lines), the predicted MMSE (dashed lines, only center pannel), and empirical experiments for
Mini-AMP (•) and streaming VB [12, 38] (×) averaged over 10 realizations of size N = 2000. We also show the results
for both fully offline (solid black line) and fully online (dash-dot black line) algorithms. Even for moderate N the
state evolution is found to almost perfectly describe Mini-AMP’s behavior. For the parameters of the center plot we
observed that for αb . 0.33 Mini-AMP is asymptotically optimal.
Finally, we have compared the Mini-AMP scheme with the streaming VB approach [12] using the meanfield 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 ibRS , since the minimum of ibRS
gives the MMSE, and since AMP is effectively trying to minimize ibRS 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 ibRS (E) for the streaming SLR problem as a function of the MSE 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 α ≥ αoffline or α ≤ ρ [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.

7

100

che
s

6

bat

4

10−8

0.3
0.2

10−12

0.1

10−16

0.0

of

0 0

2

10
8

10−4
0.4

no.

−12 −10 −8
−6 −4 −2
log10 MSE

relative batch size αb

mutual information

0.5

0

1
2
3
4
5
total amount of data α = M/N

< 10−20

Figure 2:

Phase transitions in streaming SLR. Left: evolution of the mutual information in the streaming SLR
problem as each mini-batch is processed. Parameters are set to ρ = 0.3, ∆ = 10−8 and αb = 0.35. Right: MSE of
Mini-AMP for different mini-batch sizes. Mini-AMP achieves the MMSE for αb < ρ. For αb > ρ the MMSE is zero
after processing a single batch, while for batch sizes between ρ < αb < 0.49 the Mini-AMP is suboptimal unless a
sufficient number of mini-batches is processed.

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
P (U, V, Y ) =

Y
ij

N (Yij ; Ui · Vj , ∆)

N
Y
i=1

N (Ui ; 0, IR )

M
R
Y
1 X
δ(Vj − ek ),
R
j=1

(19)

k=1

where Ui and Vj 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, 19, 23, 26, 37, 40]. 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 ≤ 4. If the number of samples is such that α = M/N < αc = R2 ∆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 αb < α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
αb > α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].
8

15

0.5

MSE(U )

MSE(U )

offline

0.0

0

1

2

3

4

error(V )

error(V )

0.4
0.2

5
0

1

2

3

0.5

0.6
0.4
0.2
0

1
2
3
total amount of data α = M/N

R2 ∆2 = 0.25

4

mini-AMP
mini-batch K-means

10

0

0.8

0.0

newsgroups

mnist

αb = 0.275
αb = 0.35

1.0

0.4

0.0

0

1

2

3

4

5

0.7
0.6

0.3
0.2
0.1

0
1
2
3
amount of data α = M/N

0.5
0.4

0 1 2 3 4 5
amount of data α = M/N

Figure 3: Clustering with Mini-AMP on synthetic Gaussian mixture data (left) and real-world data (right). Left:
mean-squared error in U (centroids) and 0-1 loss in V (labels) using different batch sizes. Solid lines give state
evolution and symbols give averages over 100 instances of size N = 1000. A transition at αc = R2 ∆2 = 0.25 prevents
Mini-AMP from giving non-zero overlap when αb < αc . Parameters are set to R = 5, ∆ = 0.1. Right: clustering
on MNIST and the 20 newsgroups dataset using Mini-AMP for model (19), with the prior on U replaced with a
non-negative Gaussian of mean zero and variance 0.1, and the noise variance ∆ estimated from the data. On MNIST,
digits of size N = 784 were clustered on K = 3 classes (0, 1 and 2), whereas for the 20 newsgroups dataset, frequency
statistics of N = 1000 words and K = 3 top-level hierarchies (comp, rec, and sci) were used. Batch sizes were set
so that αb = 0.05. Blue/red circles give the cumulative performance of Mini-AMP and mini-batch K-means [39],
respectively, averaged over a 100 different orders of presentation. For the newsgroups dataset, a 2nd pass over the
data was performed so that all labels could be recomputed with accurate estimation of the centroids; results of 1st
and 2nd passes are shown in light/dark lines respectively.

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.

References
[1] L. Bottou. Large-scale machine learning with stochastic gradient descent. In Proc. COMPSTAT, pages
177–186, 2010.
[2] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. J. Machine
Learning Research, 14(1):1303–1347, 2013.
[3] M. Opper. A Bayesian approach to online learning. In D. Saad, editor, On-line learning in Neural
Networks, pages 363–378. Cambridge University Press, 1998.
9

[4] T. P. Minka. Expectation propagation for approximate Bayesian inference. In Proc. Conf. on Uncertainty
in Artificial Intelligence, pages 362–369, 2001.
[5] I. Mitliagkas, C. Caramanis, and P. Jain. Memory limited, streaming pca. In Adv. in Neural Info. Proc.
Systems, 2013.
[6] C. Wang and Y. M. Lu. Online learning for sparse PCA in high dimensions: Exact dynamics and phase
transitions. In Proc. IEEE Info. Theory Workshop, pages 186–190, 2016.
[7] O. Kinouchi and N. Caticha. Optimal generalization in perceptrons. Journal of Physics A, 25(23):6243,
1992.
[8] M. Biehl and P. Riegler. On-line learning with a perceptron. EPL, 28(7):525, 1994.
[9] S. Solla and O. Winther. Optimal perceptron learning: as online Bayesian approach. In D. Saad, editor,
On-line Learning in Neural Networks, pages 379–398. Cambridge, 1998.
[10] D. Saad. On-line learning in neural networks, volume 17. Cambridge, 1999.
[11] P. V. Rossi, Y. Kabashima, and J. Inoue. Bayesian online compressed sensing. Phys. Rev. E, 94(2):
022137, 2016.
[12] T. Broderick, N. Boyd, A. Wibisono, A. C. Wilson, and M. I. Jordan. Streaming variational Bayes. In
Adv. in Neural Info. Proc. Systems, pages 1727–1735, 2013.
[13] D. L. Donoho, A. Maleki, and A. Montanari. Message-passing algorithms for compressed sensing. Proc.
Nat. Acad. Sci., 106(45):18914–18919, 2009.
[14] S. Rangan. Generalized approximate message passing for estimation with random linear mixing. In Proc.
IEEE Int. Symp. on Info. Theory, pages 2168–2172, 2011.
[15] L. Zdeborová and F. Krzakala. Statistical physics of inference: Thresholds and algorithms. Advances in
Physics, 65(5):453–552, 2016.
[16] M. Mézard and A. Montanari. Information, Physics, and Computation. Oxford, 2009.
[17] S. Rangan and A. K. Fletcher. Iterative estimation of constrained rank-one matrices in noise. In IEEE
Int. Symp. on Info. Theory, pages 1246–1250, 2012.
[18] F. Krzakala, M. Mézard, F. Sausset, Y. Sun, and L. Zdeborová. Probabilistic reconstruction in compressed
sensing: algorithms, phase diagrams, and threshold achieving matrices. J. of Statistical Mechanics:
Theory and Experiment, (08):P08009, 2012.
[19] T. Lesieur, C. De Bacco, J. Banks, F. Krzakala, C. Moore, and L. Zdeborová. Phase transitions and
optimal algorithms in high-dimensional Gaussian mixture clustering. Proc. Allerton Conf. on Comm.,
Control, and Computing, pages 601–608, 2016.
[20] T. Lesieur, L. Miolane, M. Lelarge, F. Krzakala, and L. Zdeborová. Statistical and computational phase
transitions in spiked tensor estimation. arXiv Preprint []:1701.08010, 2017.
[21] M. Bayati and A. Montanari. The Dynamics of Message Passing on Dense Graphs, with Applications to
Compressed Sensing. IEEE Trans. on Inf. Th., 57(2):764–785, February 2011.
[22] M. Bayati, M. Lelarge, and A. Montanari. Universality in polytope phase transitions and message
passing algorithms. Annals of Applied Probability, 25(2):753–822, 2015.
[23] J. Barbier, M. Dia, N. Macris, F. Krzakala, T. Lesieur, and L. Zdeborová. Mutual information for
symmetric rank-one matrix estimation: A proof of the replica formula. In Adv. in Neural Info. Proc.
Systems, 2016.
10

[24] J. Barbier, N. Macris, M. Dia, and F. Krzakala. Mutual information and optimality of approximate
message-passing in random linear estimation. arXiv Preprint []:1701.05823, 2017.
[25] G. Reeves and H. D. Pfister. The replica-symmetric prediction for compressed sensing with Gaussian
matrices is exact. In IEEE Int. Symp. on Info. Theory, pages 665–669, 2016.
[26] L. Miolane. Fundamental limits of low-rank matrix estimation. arXiv Preprint [math.PR]:1702.00473,
2017.
[27] E. Gardner and B. Derrida. Three unfinished works on the optimal storage capacity of networks. J. of
Phys. A: Mathematical and General, 22(12):1983, 1989.
[28] G. Györgyi. First-order transition to perfect generalization in a neural network with binary synapses.
Phys. Rev. A, 41:7097–7100, 1990.
[29] E. J. Candès, J. K. Romberg, and T. Tao. Stable signal recovery from incomplete and inaccurate
measurements. Communications on Pure and Applied Mathematics, 59(8):1207–1223, 2006.
[30] S. Rangan, P. Schniter, and A. Fletcher. Vector approximate message passing. arXiv Preprint
[]:1601.03082, 2016.
[31] B. Çakmak, O. Winther, and B. H. Fleury. S-amp: Approximate message passing for general matrix
ensembles. In IEEE Info. Theory Workshop, pages 192–196. IEEE, 2014.
[32] E. Gardner and B. Derrida. Optimal storage properties of neural network models. J. of Phys. A:
Mathematical and General, 21(1):271, 1988.
[33] F. Krzakala, J. Xu, and L. Zdeborová. Mutual information in rank-one matrix estimation. In IEEE
Info. Theory Workshop, pages 71–75, 2016.
[34] Y. Kabashima, F. Krzakala, M. Mézard, A. Sakata, and L. Zdeborová. Phase transitions and sample
complexity in bayes-optimal matrix factorization. IEEE Trans. on Info. Theory, 62:4228–4265, 2016.
[35] M. Talagrand. Spin glasses: a challenge for mathematicians: cavity and mean field models, volume 46.
Springer Science & Business Media, 2003.
[36] M. Lelarge and L. Miolane. Fundamental limits of symmetric low-rank matrix estimation. arXiv Preprint
[math.PR]:1611.03888, 2016.
[37] T. Lesieur, F. Krzakala, and L. Zdeborová. Constrained low-rank matrix estimation: Phase transitions,
approximate message passing and applications. arXiv Preprint []:1701.00858, 2017.
[38] F. Krzakala, A. Manoel, E. W. Tramel, and L. Zdeborová. Variational free energies for compressed
sensing. In IEEE Int. Symp. on Info. Theory, pages 1499–1503, 2014.
[39] D. Sculley. Web-scale k-means clustering. In Proc. Int. Conf. on World Wide Web, pages 1177–1178,
2010.
[40] R. Matsushita and T. Tanaka. Low-rank matrix reconstruction and clustering via approximate message
passing. In Adv. in Neural Info. Proc. Systems, 2013.
[41] T. Lesieur, F. Krzakala, and L. Zdeborová. MMSE of probabilistic low-rank matrix estimation: Universality with respect to the output channel. In Proc. Allerton Conf. on Communication, Control, and
Computing, pages 680–687, 2015.
[42] D. Arthur and S. Vassilvitskii. k-means++: The advantages of careful seeding. In Proceedings of the
18th Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1027–1035. Society for Industrial
and Applied Mathematics, 2007.
11

[43] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer,
R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.
Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825–2830, 2011.

12

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 Λ, Θ variables and replacing the η(A, B) function with (12).

A.1

Generalized linear models

Denote by y ∈ RM the response variable, by Φ ∈ RM ×N the design matrix, and by x ∈ RN the parameter
vector that we want to estimate. Then our generative model reads
P (y, x|Φ) =

M
Y

P (yµ |zµ ≡ Φµ · x)

µ=1

N
Y

PX (xi ).

(20)

i=1

The GAMP algorithm provides the following approximation to the marginals of x
qx (xi |Ai , Bi ) =

2
1
1
PX (xi ) e− 2 Ai xi +Bi xi ,
Zx (Ai , Bi )

(21)

and, to the marginals of z
−

(zµ −ωµ )

Vµ
e
1
P (yµ |zµ ) p
qz (zµ |yµ , ωµ , Vµ ) =
.
Zz (yµ , ωµ , Vµ )
2πVµ

(22)

The parameters A, B, ω and V are determined by iterating the GAMP equations. We denote the mean
∂η
∂
and variance of qx (A, B) by η(A, B) = ∂B
log Zx (A, B) and η 0 (A, B) = ∂B
(A, B) respectively; moreover, we
∂
define gout (y, ω, V ) = ∂ω
log Zz (y, ω, V ). The GAMP equations then read [14, 15]
ω (t) = Φx̂(t) − V (t) ◦ g (t−1) ,

V (t) = (Φ ◦ Φ) σ̂ (t) ,

gµ(t) = gout (yµ , ωµ(t) , Vµ(t) ) ∀µ,

∂ω gµ(t) = ∂ω gout (yµ , ωµ(t) , Vµ(t) ) ∀µ,

B (t) = ΦT g (t) + A(t) ◦ x̂(t) ,

A(t) = −(Φ ◦ Φ)T ∂ω g (t) ,

(t+1)

x̂i

(t)

(t)

= η(Ai , Bi ) ∀i,

(t+1)

σ̂i

(t)

(23)

(t)

= η 0 (Ai , Bi ) ∀i.

y−ω
Note that, for a Gaussian likelihood P (y|z) = N (y; z, ∆), gout (y, ω, V ) = ∆+V
. Then, by defining z̃ (t) =
(t)
(t)
(t)
y − ω and replacing V
and A by its averages, we get via the central limit theorem

V (t) ≈

N
1 X (t) (t) (t)
σ (A , Bi ),
N i=1

A(t) ≈

α
,
∆ + V (t)

(24)

where we have assumed the Φµ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 qi (xi ) which minimize
KL[

Q

i qi (xi )kP (x|Φ, y)]

= −E{qi } log P (y|Φ, x) +

N
X
i=1

13

KL[qi (xi )kPX (xi )] − log P (y|Φ)

(25)

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

P
1
1 X
(t+1)
(t)
2
x̂i
=η
Φ ,
Φµi (yµ − j6=i Φµj x̂j ) ∀i,
∆ µ µi ∆ µ

(26)

where η is defined as before. A closer inspection shows that the same equations are obtained by setting
V (t) = 0 in (23).
As shown in [38], this iteration leads to good results when performed sequentially, if the noise ∆ 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 (x|Φk , yk ) ∝
P (yk |Φk , x) Qk−1 (x), by the distribution Qk (x) that minimizes
KL[P (x|Φk , yk )kQk (x)].

(27)

Note this is the direct KL divergence, and not the reverse one KL[Q(x)kP (x|Φ, y)] that is minimized in
QN
Variational Bayes. In particular if Qk (x) = i=1 qk (xi ), then minimizing this KL divergence leads to the
following integral

Z Y
qk (xi ) = qk−1 (xi )
dxj qk−1 (xj ) P (yk |Φk , x) ∀i,
(28)
j6=i

which is tractable since we are processing a single sample yk at a time, i.e. since the likelihood consists of a
single factor. These are actually the exact marginals of P (x|Φk , yk ), 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 −V (t) g (t−1) (analogously, if a single
iteration is performed).
Algorithm 2 Assumed Density Filtering for GLMs [11]
1: initialize Λ0,i = 0, Θ0,i = 0 ∀i
2: for k ← 1 to M do
3:
compute ωk = Φk · x̂k−1 , Vk = (Φk ◦ Φk ) · σ̂k−1
4:
compute gk , ∂ω gk following (23)
5:
compute Ak , Bk following (23)
6:
compute σ̂k , x̂k following (12)
7:
accumulate Λk ← Λk−1 + Ak
8:
accumulate Θk ← Θk−1 + Bk
9: end for

A.2

Low-rank matrix factorization

Denote by Y ∈ RN ×M the matrix we want to factorize, and by U ∈ RN ×R , V ∈ RM ×R the matrices which
product approximates Y . The generative model then reads
P (Y, U, V ) =

Y

P (Yij |Wij ≡ Ui · Vj )

ij

N
Y
i=1

14

PU (Ui )

M
Y
j=1

PV (Vj ),

(29)

where Ui and Vj denote the i-th and j-th rows of U and V respectively. The algorithm provides the following
approximation to the marginal of Ui
qU (Ui |AU , BU,i ) =

T
T
1
1
PU (Ui ) e− 2 Ui AU Ui +BU,i Ui .
ZU (AU , BU,i )

(30)

and qV is analogously defined as the marginal of Vj . As in the previous case, AU ∈ RR×R and BU,i ∈ RR are
to be determined by iterating a set of equations. We denote by BU the N × R matrix which rows are given
0
by BU,i , i = 1, . . . , N . The functions ηU (A, B) = EqU U = ∇B log ZU (A, B) and ηU
(A, B) = ∇B ηU (A, B)
0
give the mean and covariance of qU , and ηV and ηV , defined analogously, the mean and covariance of qV .
In order to write the AMP equations, we first introduce
1 ∂ ln P (Yij |w = 0)
Jij = √
,
∂w
N

2
1
∂ ln P (y|w = 0)
.
β=
EP (y|w=0)
N
∂w

(31)

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

(t)

Ûi

(t)

T

(t)

(t)

BU = J V̂ (t) − βΣV Û (t−1) ,

AU = β V̂ (t) V̂ (t) ,

(t)

ΣU =

(t)

(t)
βΣU V̂ (t) ,

(t)
AV

= η(AU , BU,i ) ∀i,

N
X

(t)

(t)

η 0 (AU , BU,i ),

i=1
(t)
BV

T

= J Û

(t+1)
V̂j

=

(t)

−

(t)
(t)
η(AV , BV,j )

(t+1)
ΣV

∀j,

(32)

T

= β Û (t) Û (t) ,
=

M
X

(t)

(t)

η 0 (AV , BV,j ).

j=1

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 ηU (AU , BU ) by
!
k−1
k−1
X
X
(t)
(t)
BU,` +BU,k .
ηU
AU,` +AU,k ,
(33)
`=1

`=1

| {z }

| {z }
Θk,i

Λk−1

We assume that U is fixed and Vk changes for each batch k; thus, the calls to ηV do not change.

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 PX (x) are
( (t)
m̂ = −α Ey,z,w ∂w g(y, w, ρ − m(t) ),
m(t+1) = Ex,b xη(m̂(t) , b),

15

(34)

where we denote ρ = Ex2 , and the averages are taken with respect to P (x, b) = PX (x)N (b; m̂x, m̂) and
P (y, ω, z) = P (y|z)N (z; ω, ρ − m)N (ω; 0, m). The MSE at each step is obtained from E (t) = ρ − m(t) . For a
α
Gaussian likelihood P (y|z) = N (y; z, ∆), m̂(t) = ∆+E
(t) and we recover (9).
The fixed points of the state evolution extremize the so-called replica free energy
φ(m, m̂) =

1
mm̂ − Eb,x log Zx (m̂, b) − αEy,ω,z log Zz (y, ω, ρ − m)
2

(35)

which gives the large system limit of the Bethe free energy (extremized by AMP). The mutual information
(10) differs from φ(m) = extrm̂ φ(m, m̂) by a constant – more specifically, by the entropy of P (y|z), iRS (m) =
φ(m) − αH[P (y|z)]. For a Gaussian likelihood, H[P (y|z)] = 12 log(2πe∆).
Pk−1 (t )
(t)
(t)
In order to adapt this to the streaming setting, we introduce λk = `=1 m̂` max + m̂k and iterate
instead, for each mini-batch k
( (t)
(t)
λk = λk−1 − αb Ey,z,w ∂w g(y, w, ρ − mk ),
(36)
(t+1)
(t)
mk
= Ex,bk xη(λk , b),
with the averages now computed over P (x, bk ) = PX (x)N (bk ; λk x, λk ). These equations should be iterated
(tmax )
for t = 1, . . . , tmax , at which point we assign λk = λk−1
. The MSE on x after mini-batch k is processed is
(t

)

then given by Ek = ρ − mk max .
Note that in the small batch size limit (αb → 0), the equation for λ becomes an ODE
dλ
= −Ey,z,w ∂w g(y, w, ρ − m(λ)),
dα

(37)

which describes the performance of the ADF algorithm [3, 11].
The free energy is also easily rewritten
φk (mk , m̂k ; λk−1 ) =

1
mk m̂k − Eb,x log Zx (λk−1 + m̂k , b) − αb Ey,ω,z log Zz (y, ω, ρ − mk ),
2

(38)

or analogously, by working with λk = λk−1 + m̂k instead of m̂k
φk (mk , λk ; λk−1 ) =

1
mk (λk − λk−1 ) − Eb,x log Zx (λk , b) − αb Ey,ω,z log Zz (y, ω, ρ − mk ).
2

(39)

from which it is clear that the extrema of φk are given by the fixed points of (36).
B.1.1

Asymptotic behavior

Equations (36) can be put in the following form
(
Ek = ε(λk ),
λk = λk−1 + αb δ(Ek ),

(40)

where ε(λk ) and δ(Ek ) are functions that depend on the prior/channel respectively. Assuming ε is invertible,
we rewrite this system of equations as a function of Ek only
ε−1 (Ek ) = ε−1 (Ek−1 ) + αb δ(Ek ),

(41)

and then solve this equation for Ek ; that gives us a recurrence relation which is unsolvable in most cases. We
use instead asymptotic forms for ε and δ, obtained in the λ → ∞, E → 0 limit. For the Bernoulli-Gaussian
prior P0 (xi ) = ρ N (xi ; 0, 1) + (1 − ρ) δ(xi ), we have
ε(λ) ∼
16

ρ
,
λ

(42)

while a Gaussian likelihood gives, in the ∆ → 0 limit, δ(E) =
ρ
αb
ρ
≈
+
⇒ Ek ≈
Ek
Ek−1
Ek



1
E.

Thus, for SLR

αb
1−
ρ

k
E0 ,

(43)
α

leading to (18). In the αb → 0 limit we recover the expression obtained by [11], MSE(α) ∼ e− ρ .

B.2

Low-rank matrix factorization
(t)

(t)

Also for low-rank models the large N limit can be analyzed by taking into account that AU,V and BU,V
converge in distribution to [41]
(t)

(t)

BU ∼ N (βMV U , βMV ),

(t)

(t)

BV ∼ N (βMU V , βMU ),

AU = βMV ,
AV = βMU ,

(t)

(t)

(t)

(t)

(44)

where MU , MV ∈ RR×R are the overlap matrices between the ground truth and the estimate at time t, that
is
(t)
(t)
MU = EU ,BU U T ηU (βMV , BU ),
(45)
(t+1)
(t)
MV
= α EV ,BV V T ηV (βMU , BV ).
While computing the expectations might become unfeasible for R > 1, an ansatz of the following form can
often be used [19, 41]
MU,V = aU,V IR + bU,V JR ,
(46)
with JR denoting the R × 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 MV by
(t+1)

λV,k

(t)

= λV,k−1 + αb EV ,BV V T ηV (βMU,k , BV ).

(47)

Note that since U is fixed, the equation for MU does not change.
The replica free energy reads, in the offline case
φ(MU , MV ) =

β
Tr MU MVT − EU ,BU log ZU (βMV , BU ) − EV ,BV log ZV (βMU , BV ),
2

(48)

and we adapt it to the online case by taking into account that λV is being incremented
β
Tr MU,k (λV,k − λV,k−1 ) −
2
log ZU (βλV,k , BU ) − EV ,BV log ZV (βMU,k , BV ).

φ(MU,k , λV,k ; λV,k−1 ) =
EU ,BU

C

(49)

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 N1 kx̂(t) − x̂(t−1) k1 < 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 tmax instead. In Figure
4, we investigate the performance of Mini-AMP under the same settings of Figure 1 (center), for different
values of tmax . We observe that the performance deteriorates if convergence is not reached.

17

tmax
tmax
tmax
tmax
tmax

100

MSE(x)

10−2

=1
=5
= 10
= 50
= 1000

10−4
10−6
10−8
0
1
2
3
4
5
total amount of data α = M/N

Figure 4: Performance of the Mini-AMP algorithm for different values of tmax , under the same settings of
Figure 1 (center), and αb = 0.35. Solid lines give state evolution, and symbols results of empirical experiments
averaged over 10 realizations of size N = 2000. The performance deteriorates if the algorithm is not iterated
until convergence.

D

Experiments on real-world data

For the experiments with real data, we have used the following model
P (U, V, Y ) =

Y

N (Yij ; Ui · Vj , ∆)

ij

N
Y

N≥0 (Ui ; 0, σ 2 IR )

i=1

where
N≥0 (x; µ, σ 2 IR ) =

M
R
Y
1 X
δ(Vj − ek ),
R
j=1

(50)

k=1

R
Y
1
2
N
(x;
µ,
σ
I
)
θ(xk )
R
Z(µ, σ 2 )

(51)

k=1

is a truncated normal distribution supported on the positive quadrant of a R-dimensional space, and Z(µ, σ 2 )
ensures proper normalization.
Note that evaluating the η(A, B) function in this case is not trivial, since it depends on the following
integral
Z
R
Y
1
− 12 xT Ax+B T x
dx
e
P0 (xk )
ZU (A, B) =
Z(µ, σ 2 )
k=1
(52)

Z Y
P
R
− 21 Akk x2k + Bk +
Ak` x` xk
`6
=
k
∝
dxk P0 (xk ) e
k=1

where, in this case, P0 (xk ) = N (xk ; µk , σ 2 ) θ(xk ). We proceed by performing a mean-field approximation.
We first define
Z
2
1
∂
η̃(A, B) =
log dx P0 (x) e− 2 Ax +Bx .
(53)
∂B
for scalar A and B. Then, for each i = 1, . . . , N , we iterate
Ûik = η̃ Akk , Bik −


1X
Ak` Ûi`
2

(54)

`6=k

sequentially in k = 1, . . . , R until convergence is reached, at which point we use the values obtained for

18

assigning Û (t) in AMP. The variances are computed from
σ̃ik = η̃ 0 Akk , Bik −


1X
Ak` Ûi`
2

(55)

`6=k

and the covariance matrix used in AMP is obtained as a function of these variances

N
X




σ̃ik ,
if k = `,


i=1
σk` =
N

X

1


σ̃ik σ̃i` ,
otherwise,

 − 2 Ak`
i=1
where in order to assign the off-diagonal terms we have used a linear response approximation, σk` =
We proceed by detailing other aspects of the experiments

(56)

PN

∂ Ûik
i=1 ∂Bi` .

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.
P
(t)
Stopping criterion For each batch, Mini-AMP was iterated either for 50 steps or until N1R ik |Ûik −
P
(t−1)
(t)
(t−1)
Ûik | + M1R jk |V̂jk − V̂jk | < 10−7 .
Noise learning We do not assign a fixed value for ∆, but instead update it after each mini-batch is processed
using a simple learning rule
X
ˆk = 1
∆
(Yk,ij − Ûk,i · Vk,j )2 .
(57)
N M ij
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.

19

