Marc Lambert
email: marc.lambert@inria.fr

Silvère Bonnabel
email: silvere.bonnabel@mines-paristech.fr

Dga / Catod

Francis Bach
email: francis.bach@inria.fr

The limited-memory recursive variational Gaussian approximation (L-RVGA)

published or not. The documents may come

The limited-memory recursive variational Gaussian approximation (L-RVGA)

Introduction

In machine learning and statistics, Bayesian inference aims at estimating the full posterior distribution of a parameter of interest θ, to better track the uncertainty of the learning process. Exact Bayesian inference is generally not tractable, and approximations are required. Variational approximation (Hinton andvan Camp, 1993, Jordan et al., 1998) consists in rewriting the estimation problem as an approximate optimization problem over a restricted class of posterior distributions, such as Gaussian distributions, as advocated in the present paper.

More precisely, assume we want to approximate the posterior distribution p(θ|Y N ) of a Bayesian parameter θ given N observations Y N = y 1 , . . . , y N . The variational Gaussian approximation (Barber and Bishop, 1998b[START_REF] Opper | The variational Gaussian approximation revisited[END_REF][START_REF] Challis | Gaussian Kullback-Leibler approximate inference[END_REF] consists in minimizing the Kullback-Leibler divergence between this unknown posterior and a restricted class of 1 parameterized distributions, namely a Gaussian distribution q(θ|µ, P ) ∼ N (θ|µ, P ), leading to the following optimization problem: min µ,P KL (q(θ|µ, P ) p(θ|Y N )) = q(θ|µ, P ) log q(θ|µ, P ) p(θ|Y N ) dθ.

(1)

In the recursive variational approach, the observations at each time step consist only of the last sample y t , and the previous Gaussian estimate q t-1 serves as a prior that sums up information obtained so far. The posterior is then approximated by p(θ|y t ) ∝ p(y t |θ)q t-1 (θ). This leads to a recursive variational approximation problem that writes at each step as follows: q t (θ) = arg min µ,P KL (N (θ|µ, P ) c t p(y t |θ)q t-1 (θ)) ,

(2) where c t is a normalization constant that disappears when we consider the derivatives with respect to the variational parameters. It can be shown that the necessary stationary conditions for (2) yield the following updates, which are implicit-see [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF]:

µ t = µ t-1 + P t-1 E θ∼N (µt,Pt) [∇ θ log p(y t |θ)] (3)
P -1 t = P -1 t-1 -E θ∼N (µt,Pt) [∇ 2 θ log p(y t |θ)]. (4)
We recognize a second-order online algorithm doing an adaptive gradient descent on the stochastic loss function t (θ) =log p(y t |θ, x t ). Although this algorithm seamlessly scales with the number of observations N , it has quadratic cost in the parameter dimension d, hindering its use in high-dimensional problems. In this paper, we consider a factor analysis (FA) structure to approximate the precision matrix P -1 at each time step with a "low rank + diagonal" matrix, that is, P -1 ≈ W W T + Ψ where W ∈ M(d × p) is a rank p matrix with p d and Ψ ∈ M d (R) is a diagonal matrix. This choice of decomposition makes it possible to manipulate matrices of size d × d while only storing d + d × p parameters in memory.

We reformulate the considered factorization in a recursive way such that the implicit equation (4) becomes:

W t-1 W T t-1 + Ψ t-1 -E θ∼N (µt,(WtW T t +Ψt) -1 ) [∇ 2 θ log p(y t |θ)] ≈ FA W t W T t + Ψ t . (5)
In the particular case of a linear regression model y t = x T t θ + w t with output noise w t ∼ N (0, 1) and with x t the input associated to the observation y t , this recursive update becomes explicit:

W t-1 W T t-1 + Ψ t-1 + x t x T t ≈ FA W t W T t + Ψ t . (6)
We recognize a classical factor analysis problem for the scaled empirical covariance of the inputs N t=1 x t x T t , but done recursively. We show that the variational parameters W t and Ψ t can be computed through a recursive variant of the EM algorithm that is free of step tuning and provides a good approximation of the precision matrix after only one pass through the data.

We then generalize this recursive EM approach in the nonlinear case. We have to assume additional approximations in this case: the Hessian term in ( 5) is approximated with an explicit outer product and the expectations are approximated by sampling. To this aim we propose a novel sampling method to sample from N (µ, (W W T + Ψ) -1 ) without storing a d × d matrix, resulting in updates which scale linearly with d in terms of both computation and storage costs.

The paper is organized as follows: in Section 2 we discuss how our approach differs from several related studies on large-scale variational inference. In Section 3, we model our problem as a two-stage variational approximation scheme and show how the second stage can be recast as an expectationmaximization (EM) problem. Considering first the linear case, we reformulate our solution as a recursive EM algorithm with computation and storage cost being linear in d. In Section 4, we extend our approach to the nonlinear case using several approximations when updating the covariance matrix. We also introduce our novel method for sampling from a Gaussian whose dispersion is encoded by a precision matrix structured with a factor analysis model.

In Section 5 our algorithm is evaluated on synthetic data for large-scale matrix approximation and linear and logistic regression.

Related work

Variational inference aims to approximate the distribution of latent parameters (Hinton andvan Camp, 1993, Jordan et al., 1998) in Bayesian machine learning, and has received a lot of attention in recent years for its ability to estimate a distribution without computing its normalization constant. Although the approximation may provide biased results, it converges much faster than the gold standard, Markov Chain Monte Carlo or MCMC [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF], which is eventually accurate but can be slow. Moreover, the model used in variational methods can be better specified using latent parameters (Minh-Ngoc et al., 2016, Linda andDavid, 2013).

To tackle large-scale datasets, stochastic solvers are used to compute the unknown parameters leading to general frameworks such as the black box variational inference algorithm [START_REF] Ranganath | Black box variational inference[END_REF]. In this article, we consider variational Gaussian approximation (VGA) where the variational parameters boil down to a mean vector and a covariance matrix. An approximation of the covariance matrix is needed to tackle high-dimensional problems in VGA. [START_REF] Challis | Gaussian Kullback-Leibler approximate inference[END_REF] have shown that, for a generalized linear model (GLM), the KL divergence is convex in the square root of the covariance matrix and have proposed different types of approximations for this form. Other works consider the factor analysis structure in the context of variational inference. [START_REF] Ghahramani | Variational inference for bayesian mixtures of factor analysers[END_REF] used a mean field variational Bayes to obtain closed-form updates for the factor analysis parameters of each Gaussian of a mixture model. Their model is more general than ours since they use a mixture of Gaussian but all observations are required at each iteration. [START_REF] Victor | Gaussian variational approximation with a factor covariance structure[END_REF] estimated the factor analysis parameters of a Gaussian using a stochastic gradient descent which allows the observation to be processed sequentially. We also consider a factor analysis structure in this article but we do not directly optimize the parameters to avoid using stochastic gradient descent which can be difficult to tune. [START_REF] Mishkin | Slang: Fast structured covariance approximations for bayesian deep learning with natural gradient[END_REF] solved the variational problem in high dimension using a singular value decomposition (SVD) to compute a factor analysis structure for the precision matrix at each step of a natural gradient descent. The two approximation schemes, Gaussian approximation, and factor analysis approximation are performed sequentially. Our contribution extends this approach and differs from it in several respects. First, we do not consider a natural gradient descent algorithm on the variational parameters but we use instead the RVGA scheme which requires no step tuning [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF]. Second, we propose a parameter-free recursive EM algorithm to compute the factor analysis precision matrix approximation without explicitly using an SVD. SVD may not be efficient in high dimensions and the decoupled approach proposed in [START_REF] Mishkin | Slang: Fast structured covariance approximations for bayesian deep learning with natural gradient[END_REF] does not allow approximating accurately the correlation between the low-rank part and the diagonal part of the factor analysis structure. Moreover, we show that our approach has a nice interpretation as a two-stage variational inference that performs two "I-projections" [START_REF] Csiszár | Information theory and statistics: A tutorial[END_REF] at each step: the first one on the space of Gaussian distributions, the second one on the space of Gaussian distributions with a structured precision matrix constrained to be "diagonal + low-rank".

Beyond variational inference, our work is also related to second-order stochastic optimization where the Hessian is approximated online to fit in the memory like in Adagrad [START_REF] Duchi | Adaptive subgradient methods for online learning and stochastic optimization[END_REF] or TONGA [START_REF] Roux | Topmoumoute online natural gradient algorithm[END_REF]. Interesting connections may be done also with large-scale Kalman filtering like the SEEK filter [START_REF] Dinh | A singular evolutive extended Kalman filter for data assimilation in oceanography[END_REF] which approximates the covariance matrix with a singular value decomposition or the ensemble filter [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics[END_REF] which approximates it with sampling.

Two-stage variational approximation

To introduce a variational loss for factor analysis, we recast the original recursive problem (2) as a two-stage variational problem as follows:

Two-stage variational approximation q * t (θ) = N (θ|µ * t , P * t ) = arg min µ,P KL(N (θ|µ, P ) c t p(y t |θ)q t-1 (θ)) (7) q t (θ) = N (θ|µ * t , Pt ) = arg min W,Ψ KL(N (θ|µ * t , (W W T + Ψ) -1 ) q * t (θ)). (8)
The first variational problem is the projection of the posterior distribution onto the space of Gaussian distributions whereas the second variational problem is the projection of a Gaussian onto the space of structured Gaussian having a low rank plus diagonal structure for the precision matrix. This factorization aims to exploit the parsimony of the covariance matrix to limit the memory requirement.

The first projection was the object of our previous work [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF] and leads to the RVGA update giving the optimal solution for the mean (3) and covariance (4), whereas the second one can be addressed through an expectation-maximization algorithm as we show now.

Low rank + diagonal approximation via EM

We first show that the divergence in (8) may be algebraically related to a maximum likelihood problem. Indeed, since the means are the same on both sides of the KL divergence, we may ignore them and rewrite the update as follows:

KL(N (θ|µ * t , (W W T + Ψ) -1 ) N (θ|µ * t , P * t )) (9) = KL(N (θ|0, P * t -1 ) N (θ|0, W W T + Ψ)) (10) = 1 2 Tr((W W T + Ψ) -1 P * t -1 ) + 1 2 log det(W W T + Ψ) - 1 2 log det(P * -1 ) - d 2 . ( 11
)
Let us consider K samples v 1 , . . . , v K supposed to be centered and such that their empirical covariance is

S K = 1 K K i=1 v i v T i = P * t -1
. The log-likelihood on these samples is:

max W,Ψ log N (v 1 , . . . , v K |0, W W T + Ψ) (12) = max W,Ψ - K 2 Tr((W W T + Ψ) -1 1 K K i=1 v i v T i ) - K 2 log det(W W T + Ψ) - K 2 d log(2π), (13)
where we have used the relation

K i=1 v T i (W W T + Ψ) -1 v i = Tr((W W T + Ψ) -1 K i=1 v i v T i ).
We see that minimizing the divergence ( 11) is equivalent to maximizing the log-likelihood (13).

The variational parameters can be obtained by zeroing the derivative of the maximum likelihood (ML) and using a singular value decomposition to find the solution. We will rather consider an expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] which better scales to highdimensional problems and which is guaranteed to increase the (total) likelihood at each step (see [START_REF] Donald | EM algorithms for ML factor analysis[END_REF]. It can be shown that the EM approach computes implicitly a singular value decomposition to find the latent factors but in a memory-efficient way. The connection between the fixed-point ML, the fixed-point EM approach, and the associated eigenvalue algorithms is discussed in more detail in Appendix B.

To apply the EM algorithm we need to introduce latent variables z 1 , . . . , z K such that our samples take the form v i = W z i + ε i where ε i ∼ N (0, Ψ) and i = 1, . . . , K.

At the expectation step (E-step), the latent variables z i are estimated conditionally on the observation and the current variational parameters W and Ψ using the conditioning formula:

p(z i |v i , W, Ψ) = N (M -1 W T Ψ -1 v i , M -1 ) where M = I p + W T Ψ -1 W. ( 14
)
At the maximization step (M-step), the variational parameters are adjusted to maximize the expected total likelihood defined by:

E[log p(v 1 , z 1 , . . . , v K , z K |W, Ψ)] = E[ K i=1 log p(v i |z i , W, Ψ) + K i=1 log p(z i )] (15) = E[- K 2 (v i -W z i ) T Ψ -1 (v i -W z i ) - K 2 log det Ψ - K 2 log(2π) + K i=1 log p(z i )], (16)
where the expectations are taken over the conditional distribution z i ∼ p(z i |v i , W, Ψ), where parameters W, Ψ are there fixed to their current value. After some calculations recapped in Appendix A.1, the EM algorithm for the factor analysis problem yields the following updates:

E-Step: E[z i |v i ] = M -1 W T Ψ -1 v i E[z i z T i |v i ] = M -1 + E[z i |v i ]E[z i |v i ] T . (17)
M-Step:

W (n) = K i=1 v i E[z i |v i ] T ( K i=1 E[z i z T i |v i ]) -1 Ψ (n) = diag 1 K K i=1 v i v T i -W (n) 1 K K t=1 E[z i |v i ]v T i , (18)
where the (n) stands for "new".

The expectation and maximization steps can be computed through a single update leading to the following fixed-point scheme (see Appendix A.1 for further details):

W (n) = S K Ψ -1 W (I p + M -1 W T Ψ -1 S K Ψ -1 W ) -1 , (19)
where

M = I p + W T Ψ -1 W Ψ (n) = diag(S K -W (n) M -1 W T Ψ -1 S K ) (20) W = W (n) , Ψ = Ψ (n) .
This update no longer depends on the samples v i nor on the latent parameters z i explicitly. It only makes use of

S K = 1 K K i=1 v i v T i = P * t -1 . The matrix P * t -1
is the output of the R-VGA update (4) and depends on the previous estimate P t-1 -1 , supposed already in a factor analysis form, such that we have the relation:

P * t -1 = W t-1 W T t-1 + Ψ t-1 -E θ∼N (µt,Pt) [∇ 2 θ log p(y t |θ)]. (21)
Substituting this expression into our fixed point EM algorithm, we obtain a closed form for the second-stage variational update. But this update involves the d × d matrix ∇ 2 θ log p(y t |θ), which does not hold in memory in large-scale problems. Moreover, this matrix is not well defined since it depends on an expectation under unknown parameters E θ∼N (µt,Pt) .

We will first consider a linear model where the Hessian matrix and expectations disappear and updates reduce to simple recursive equations leading to Algorithm 1 which is memory-efficient. We will then show in Section 4 that this same algorithm can be used in the nonlinear case if we approximate the Hessian matrix and the expectations in a suitable way.

The linear case: a recursive EM algorithm for factor analysis

Let us consider for now the simpler case of a linear model y t = x T t θ + w t , where w t ∼ N (0, σ 2 w ). The variance σ 2 w is a scaling factor we will set equal to 1. The RVGA updates (3)-( 4) is then equivalent to the following explicit updates, see [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF]:

µ t = µ t-1 + P t x t (y t -x T t µ t-1 ) (22)
P -1 t = P -1 t-1 + x t x T t . (23)
If we apply our two-stage variational approximation in this linear case, the first stage approximation (7) is exactly solved by the update above and the second stage approximation (8) addresses limitedmemory requirements. Using the factor analysis approximation in a recursive way, we obtain:

µ t = µ t-1 + (W t W T t + Ψ t ) -1 x t (y t -x T t µ t-1 ), (24)
W t W T t + Ψ t ≈ F A W t-1 W T t-1 + Ψ t-1 + x t x T t . (25)
Using the Woodbury formula 24) can be rewritten in a limited-memory fashion involving only operations linear in d:

(W t W T t + Ψ t ) -1 x t = Ψ -1 t (x t -W t M -1 t (W T t Ψ -1 t x t )) with M t = I p + W T t Ψ -1 t W t , (
(24) ⇔ µ t = µ t-1 + Ψ -1 t (x t -W t (I p + W T t Ψ -1 t W t ) -1 (W T t Ψ -1 t x t ))(y t -x T t µ t-1 ). ( 26
)
To address (25), we can use our fixed-point equation replacing the matrix S K in (20) by

W t-1 W T t-1 + Ψ t-1 + x t x T t .
This defines a recursive EM algorithm, namely Algorithm 1, which consists in successively performing a few cycles on the fixed-point equation ( 20), where we expand and rearrange the terms to avoid any operations requiring storing or multiplying d × d matrices. The Algorithm 1 is given in a more general setting where we solve a factor analysis approximation of the form:

W t W T t + Ψ t ≈ F A α t (W t-1 W T t-1 + Ψ t-1 ) + β t x t x T t , (27)
with α t , β t scalar constants being equal to 1 in this section. We will show in the next Section that we can use this same algorithm in the nonlinear case (or for the factorization of a covariance matrix) by using different values for the parameters α t , β t and replacing the entries x t by rectangular matrixes.

Using the recursive EM in this way gives a factor analysis decomposition of the precision matrix P * t -1 , which corresponds to an increasing amount of information P -1 0 + t i=1 x i x T i as more inputs are processed, where P 0 corresponds to the covariance of the prior.

This prior matrix P 0 can not be stored in memory and the algorithm 1 contains an initialization procedure which computes W 0 and Ψ 0 such that W 0 W T 0 + Ψ 0 ≈ P -1 0 , see Appendix D for more details.

Algorithm 1: Recursive expectation-maximization algorithm for solving factor analysis approximation (27) Result: W, Ψ Given N inputs x 1 , . . . , x N in high dimension d; Given a latent dimension p d; Given a prior covariance P 0 ; -Initialization (consistent with P 0 , see D)-

Initialize W ∈ M(d × p) and Ψ ∈ M d (R) diagonal: Ψ = Ψ 0 ; W = W 0 ; -Update- for t ← 1 to N do Access an input x t ; for k ← 1 to nbInnerLoop do M = I p + W T Ψ -1 W ; V = β t x t (x T t Ψ -1 W ) + α t (W t-1 (W T t-1 Ψ -1 W ) + Ψ t-1 Ψ -1 W ) ; W (n) = V (I p + M -1 W T Ψ -1 V ) -1 ; Ψ (n) = β t x t * x t + α t (W t-1 * W t-1 + Ψ t-1 ) -W (n) M -1 * V ; W = W (n) ; Ψ = Ψ (n) ; end W t-1 = W ; Ψ t-1 = Ψ;
end /* From 1 to 3 inner loops (nbInnerLoop) may be sufficient to make the algorithm converge. The initialization derived in Section D compute W 0 and Ψ 0 such that W 0 W T 0 + Ψ 0 ≈ P -1 0 . The operation " * " which appears at the last line aims to compute X * Y = diag(XY T ) in a memory-efficient way. It is applied on two matrices X and Y of the same size d × p. If p = 1, this operator matches the component-wise operator :

X * Y = x y. If p > 1, it writes X * Y = p i=1 x[., i] y[., i].
The diagonal matrix Ψ is stored and used as a vector: any multiplication of the diagonal matrix with a vector may be computed faster as an element-wise product: Ψ -1 W = W/ψ where ψ = diag(Ψ) and /ψ operates on columns of W and W T Ψ -1 gives as well W T /ψ T where /ψ T operates on the rows of W T . */

4 The limited memory recursive variational Gaussian approximation (L-RVGA)

In this Section, we consider the extension of our algorithm when the observations depend nonlinearly on the latent parameter θ. We first extend the previous linear case to generalized linear models in Section 4.1. Section 4.2 considers the wholly general case, which necessitates additional approximations.

L-RVGA for generalized linear models

The latter approach may be extended whenever the observation y t follows an exponential family distribution, which may be advantageous in the context of classification problems. An exponential family distribution [START_REF] Pitcher | Review: O. Barndorff-Nielsen, Information and exponential families in statistical theory[END_REF] takes the form:

p(y t ) = c(y t ) exp(η T y t -F (η)), ( 28
)
where F is the log partition function and c(y) is a normalization function. It can represent a large family of distributions like the Gaussian, multinomial, or Dirichlet distribution. The natural parameter η is related to the expectation parameter m = E[y] through the link function g such that η = g(m). In the machine learning framework, the natural parameter is also related to the input x and the hidden latent parameter θ through a linear function η = θ T x, this model is called the generalized linear model with a canonical natural parameter. The advantage of this model is that the Hessian term takes the form of a (generalized

) outer product -∇ 2 θ log p(y t |θ) = x t ∇ 2 η F (η(θ))x T t .
In this case the RVGA updates (3)-( 4) become:

P -1 t = P -1 t-1 + γ t x t x T t (29) µ t = µ t-1 + ξ t P t-1 x t , (30)
where we have introduced the scalar weighting parameter γ t and the scalar error ξ t which depend on θ as follows

γ t = E θ∼N (µt,Pt) [Cov(y t |θ)], (31)
ξ t = y t -E θ∼N (µt,Pt) [m(y|θ)], (32)
and where we have used the properties of exponential family [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF]. The factorized approximation of this update can be performed through Algorithm 1, since the factor analysis approximation 29) fits into the problem ( 27), letting α t = 1, β t = γ t . The computation of the scalar terms γ t and ξ t is nontrivial and requires resorting to approximations, though. Those approximations may rely on the tools developed in the next subsection, devoted to the general case. However, in the case of logistic regression with m(y|θ) = 1 1+exp(-θ T x) and Cov(y|θ) = m(y|θ)(1m(y|θ)) these parameters were shown in our prior work on RVGA Lambert et al. (2021) to be easily obtainable numerically and will serve as a baseline for our experimentation in Section 5.3.

∇ η F (η(θ)) = m(y|θ) and ∇ 2 η F (η(θ)) = Cov(y|θ), see
W t-1 W T t-1 +Ψ t-1 +γ t x t x T t ≈ FA W t W T t + Ψ t of (

Limited memory RVGA (L-RVGA): nonlinear model

As before, we suppose that the observation y t follows the exponential family distribution (28) but with a nonlinear dependence of the form η = h(θ, x) with h is a nonlinear function. In our setting, we need to approximate the matrix P -1 t recursively, as follows

P -1 t = W t-1 W T t-1 + Ψ t-1 -E θ∼N (µt,Pt) [∇ 2 θ log p(y t |θ, x t )] (33) ≈ W t-1 W T t-1 + Ψ t-1 + X t X T t (34) ≈ FA W t W T t + Ψ t ,
where X t is a rectangular matrix supposed to fit into memory which will be defined more precisely in Section 4.2.4. Using this approximation, we can run the recursive EM Algorithm 1 in a memoryefficient way, replacing the input x t by the matrix X t . The mean is computed as before:

µ t = µ t-1 -(W t-1 W T t-1 -Ψ t-1 ) -1 E θ∼N (µt,Pt) [∇ θ log p(y t |θ, x t )], (35)
where we can use the Woodbury formula again to compute efficiently the weighting term

(W t-1 W T t-1 - Ψ t-1 ) -1 .
A practical implementation of these updates raises the following difficulties:

1. Computing P -1 t and µ t in ( 33)-( 35) involves computing an expectation over a distribution parameterized by µ t and P t , i.e., the scheme is implicit.

2. As no closed form is generally available for the computation of the expectations E θ∼N (µ,P ) , one may resort to Monte-Carlo sampling. However, in this paper we maintain a limited memory approximation of the inverse of the covariance matrix, that is, P -1 = W W T + Ψ and we need to sample from a Gaussian N (µ, P ) without storing or inverting a d × d matrix.

3. The Hessian matrix ∇ 2 θ log p(y t |θ, x t ) (or its averaged value) must be computed and stored with linear cost in the dimension d.

In the remainder of the section, we address all these points.

Using extra-gradients for the implicit scheme

The more direct way to address Point 1 above is to open the loop and to replace the expectations under N (µ t , P t ) with expectations under N (µ t-1 , P t-1 ). However, experiments we conducted showed that this naive scheme can lead to instability. The importance of managing the implicit scheme was the object of our previous work [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF], where we developed closed-form formulas to solve the implicit scheme in the linear and logistic regression case. In the general case, one may resort to extra-gradient, i.e., to update first the covariance and then the mean, and to iterate twice:

Iterated RVGA Pt -1 = P -1 t-1 -E θ∼N (µ t-1 ,P t-1 ) [∇ 2 θ log p(y t |θ)] μt = µ t-1 + Pt E θ∼N (µ t-1 ,P t-1 ) [∇ θ log p(y t |θ)] P t -1 = P -1 t-1 -E θ∼N (μ t , Pt ) [∇ 2 θ log p(y t |θ)] µ t = µ t-1 = P t E θ∼N (μ t , Pt ) [∇ θ log p(y t |θ)]. (36)
It turns out that this iterated scheme is equivalent to the "Mirror Prox" algorithm [START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF], a.k.a., extra-gradient, with a unit step size and applied to the function f (µ, P ) = E θ∼N (µ,P ) [log p(y|θ)]. Mirror Prox is known to help convergence on a large set of problems (convex optimization, variational inequalities) and we have experimentally observed that this iterated scheme reduces significantly the bias of our estimator (see Appendix G for further details).

However, when we combine extra-gradient with factor analysis, the extra covariance update can make the Mirror Prox scheme unstable. We have observed it is then preferable to skip the extra covariance update, that is, the third line of the iterated scheme above (see Appendix G). This choice has been made in Section 5.4 dedicated to experiments in the general case.

Gaussian sampling from the precision matrix

To address Point 2 above, we need to sample efficiently from E θ∼N (µ,(W W T +Ψ) -1 ) without storing a d × d matrix. This problem was already addressed by [START_REF] Mishkin | Slang: Fast structured covariance approximations for bayesian deep learning with natural gradient[END_REF] and Ambikasaran et al. ( 2014) using a square-root form and two Cholesky decompositions on the latent space. However, we propose here a faster method inspired by the ensemble Kalman filter [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics[END_REF] which does not require computing any Cholesky decomposition. Our method is close to the one developed in [START_REF] Orieux | Sampling high-dimensional gaussian distributions for general linear inverse problems[END_REF] for a sparse precision matrix but is more suitable for factor analysis.

Proposition 1. Let us define the quantities:

M = I p + W T Ψ -1 W and L = Ψ -1 W M -1 . Draw x ∼ N (0, Ψ -1
) and ∼ N (0, I p ) independently, and define

x + = x + L( -LW T x) = (I d -LW T )x + L .
We have then x + ∼ N (0, P ), with

P = (W W T + Ψ) -1 = Ψ -1 -Ψ -1 W (I p + W T Ψ -1 W ) -1 W T Ψ -1 .
Proof. We have obviously E(x + ) = 0. Moreover using the independence of the variables

E x + (x + ) T = (I d -LW T )Ψ -1 (I d -LW T ) T + LI p L T (37) = Ψ -1 -LW T Ψ -1 -Ψ -1 W L T + L[W T Ψ -1 W + I p ]L T (38) = Ψ -1 -Ψ -1 W (I p + W T Ψ -1 W ) -1 W T Ψ -1 , (39)
since the three rightmost terms of ( 38) are all equal to ±Ψ -1 W (

I p + W T Ψ -1 W ) -1 W T Ψ -1 .
This suggests that starting from a decomposition of the form W W T + Ψ of the precision matrix P -1 , we know how to sample from a law N (0, P ). We need to draw x 1 , . . . , x K from N (0, Ψ -1 ) and 1 , . . . , K from N (0, I p ), and to let x

+ i = x i + L( i -LW T x i ) = (I d -LW T )x i + L i , for 1 ≤ i ≤ K.
All operations involved are linear in d, so that drawing K samples is of order O(Kd). As soon as K is kept moderate with respect to d, the complexity is linear in d.

Dealing with the d × d Hessian matrix

In this paragraph, we address Point 3 of the list above, where we propose an approximation of the Hessian based on an outer product. The Gauss-Newton approximation, which consists in approximating the Hessian with the outer product of the gradients, called also empirical Fisher, does not capture well second-order information (see for instance [START_REF] Kunstner | Limitations of the empirical Fisher approximation for Natural gradient descent[END_REF]). As we have supposed the probability distribution belongs to an exponential family, we should rather consider the Generalized Gauss-Newton (GGN) approximation [START_REF] Martens | New insights and perspectives on the natural gradient method[END_REF] which better approximates the local curvature [START_REF] Kunstner | Limitations of the empirical Fisher approximation for Natural gradient descent[END_REF]. The GGN approximation exploits the structure of the exponential family p(y t |η = h(θ, x)) as a composition of functions involving the natural parameter η and the nonlinear model h. If the Hessian oflog p with respect to θ can be complicated, the Hessian with respect to η is simply Cov(y|θ). We can use this property through the chain rule to generalize the Gauss-Newton approximation. Our expected Hessian term can be written as follows:

-E θ [∇ 2 θ log p(y t |θ)] ≈ E θ [ ∂h ∂θ Cov(y|θ) ∂h ∂θ T ] = E θ [F(θ)], (40)
where

F(θ) = E y∼p(y|θ) [-∇ 2 θ log p(y|θ)
] is the Fisher matrix. Under the GGN approximation, the covariance update no longer depends on the labels y t and only depends on the inputs x t , as in the linear case. The derivation is detailed in Appendix F. We can now combine this approximation with the approximation of the expectation and the implicit scheme to implement our final update.

Final algorithm

We now define the form of matrix X t involved in the outer product in (34). Assuming that we make the scheme explicit using extra-gradients (as described in Section 4.2.1) and that we generate K samples to approximate the expectation with Ensemble sampling described in Section 4.2.2, the GGN approximation becomes:

P -1 t ≈ GGN W t-1 W T t-1 + Ψ t-1 + E θ∼N (µt,Pt) [ ∂h ∂θ Cov(y|θ) ∂h ∂θ T ] (41)
≈ Extragrad W t-1 W T t-1 + Ψ t-1 + E θ∼N (µ t-1 ,(W t-1 W T t-1 +Ψ t-1 ) -1 ) [ ∂h ∂θ Cov(y|θ) ∂h ∂θ T ] (42)
≈ Sampling W t-1 W T t-1 + Ψ t-1 + 1 K K i=1 ∂h ∂θ (θ i )Cov(y|θ i ) 1/2 Cov(y|θ i ) 1/2 ∂h ∂θ T (θ i ) (43) = W t-1 W T t-1 + Ψ t-1 + 1 K K i=1 c i c T i where c i = ∂h ∂θ (θ i )Cov(y|θ i ) 1/2 (44) = W t-1 W T t-1 + Ψ t-1 + X t X T t where X t = 1 √ K c 1 • • • c K (45) ≈ FA W t W T t + Ψ t , (46)
where the last line refers to the FA approximation developed in Section 3.1 applied to the mini-batch matrix X t of size d × K. As long as the number of samples K d, the memory cost is kept linear in d. Regarding the other approximations, the approximation of the expectation with sampling seems not very sensitive and we recommend using a limited number of samples to reduce the memory cost. Indeed, we have observed few samples are sufficient to obtain a good approximation of the logistic regression problem considered in Section 5.4.

Experiments

Experiments have been performed on synthetic data for different classes of problems, from simpler linear problems to more difficult nonlinear problems. The first class of problems addressed in Section 5.1 deals with the approximation of large-scale covariance matrices that are too large to hold in memory. We show our approach requires far less memory than the batch EM algorithm which requires storing a d × d matrix in memory and competes with the online EM algorithm [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF]. These results are directly applied to linear regression in Section 5.2 where we propose computationally cheap updates for the linear Kalman filter in high dimensional spaces. In Section 5.3, we consider the logistic regression problem for which we can solve the recursive variational scheme in a closed form at each step [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF]. Finally, we turn to more general cases in Section 5.4 where we assess our approximation method for L-RVGA updates in the general case. The logistic regression problem, where we know how to compute the expectation analytically, offers a baseline for comparison purposes.

We observe the mirror-prox method combined with ensemble sampling reaches the baseline results in terms of KL divergence. This is promising regarding the generalization of our algorithms to arbitrary nonlinear problems.

Limited-memory approximation of large-scale covariance matrices

In this section, we assess the limited-memory recursive EM algorithm on large-scale empirical covariance matrices. In this setting, we do not update a precision matrix but a covariance matrix with a moving average, this amounts to using Algorithm 1 with α t = (t -1)/t and β t = 1/t (see Appendix C for details). To assess the method, we generate the data from an actual low-rank covariance plus diagonal matrix S = Diag(ψ) + W W T for which the factor analysis approximation can be exact. We construct such a matrix with random parameters where W is of size d × p and ψ is a vector of size d. We then generate N samples from the zero-mean Gaussian distribution equipped with this covariance matrix S. To make the experiments more appealing, we have also assessed our algorithm on the NIPS Madelon dataset and the Breast Cancer real dataset. We have kept only the input data, using the LIBSVM library [START_REF] Chang | Libsvm: A library for support vector machines[END_REF], and normalized them using the same process as for the synthetic dataset.

We run the batch EM algorithm on this set of N samples and our limited memory EM algorithm recursively on each sample. The recursive EM generally converges to the result given by the batch version after only one pass through the data as shown in Figure 1. Our recursive EM is also compared to the online EM algorithm [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF] which is also memory efficient. The derivation of the online EM for factor analysis and the implementation detailed have been moved to Appendix E. While our method is parameter-free it yields results similar to the online EM on the synthetic dataset and better results on the real datasets as shown in upper plots of Figure 1.

The advantage of online or recursive EM versions for factor analysis is that they maintain a memory cost linear in d, which is an important feature in high dimension. To illustrate this feature, we performed large-scale covariance matrix approximation in dimension 1 million. We see in Table 1 that the proposed algorithm may be performed using a regular laptop, contrary to the batch EM.

Algorithm

Nb. iterations Time per iteration Memory cost Batch EM (estimated) The lower plots show additional results on two covariance matrices computed from the Madelon and Breast Cancer datasets. We show in dash green line our recursive EM algorithm, which is parameter-free, and in the red plain line, the online EM tuned with a learning rate γ = 1 t 0.6 (see Appendix E for implementation details). We show also the batch EM, in the solid blue line. The online and batch methods should be compared with respect to the number of observations processed. Batch estimation requires processing all the data several times, although more efficiently since 5 steps are enough to converge. The online variants use only a single passe. With 3 inner loops in dimension d = 100 (upper left plot), the recursive EM is clearly superior to the online EM. In the higher dimension (upper right plot), we achieve a similar loss to the online EM using only 1 inner loop in the recursive EM. Moreover, the recursive EM gives better results than the other algorithms on real datasets as illustrated in the lower left plot and lower right plot.

Fixed ( N ) - 8000 GB Online EM-1 inner loop N 1 s 1.6 GB Recursive EM-1 inner loop N 2 s 0.8 GB

Application to Bayesian linear regression

We apply in this section the previous covariance approximation to derive a limited memory version of the linear regression problem as described in Section 3.2. In the linear general case, we have an analytical form of the solution given by the linear Kalman filter with a static state which is our baseline.

The linear Kalman filter is a parameter-free online algorithm that gives the exact solution of a linear regression problem after only one pass on the dataset, that is, it recursively computes exactly the posterior given all the data it has considered so far. However, it requires estimating online the d × d covariance matrix of the inputs, which is intractable in high dimension. This matrix is approximated with the recursive EM algorithm in this experiment to provide large-scale linear regression.

The inputs x t are generated using the following Gaussian distribution:

x t ∼ N (0, C) with C = M T Diag(1, 1/2 c , • • • , 1/d c )M, ( 47
)
where M is an orthogonal rotation matrix and c a coefficient driving the condition number (c = 1 by default). Since the matrix C is rotated, it is not directly available in a factor analysis form. The inputs are normalized on average. The outputs y t are generated with a normal noise:

y t = x T t θ * + ε t with ε t ∼ N (0, 1), (48)
We consider the Bayesian setting where θ is supposed to be a Gaussian random value and we want to estimate the parameters of its distribution N (µ, P ) where P -1 = W W T + Ψ, given all the examples (x t , y t ) seen so far. We suppose that the prior follows a centered isotropic Gaussian distribution θ 0 ∼ N (0, σ 2 0 I d ) and compute the initial values of the factor analysis parameters W 0 and Ψ 0 as detailed in Appendix D. Our recursive scheme is sensitive to the parameter σ 0 : a high value of σ 0 (flat prior) may lead to bad conditioning at the first steps and make our recursive algorithm diverge. A low value for σ 0 (strong prior) speeds up the convergence of the algorithm. This value defines the shape of the posterior and the difficulty of the estimation problem and can be used to test the robustness of our algorithm in the following experiments. In this Section, we consider a prior deviation σ 0 = 1.

When d = p we have checked that our algorithm matches the linear Kalman filter's estimates. To assess how the factor analysis approximation degrades the results, we have done experiments for different values of the latent space dimension p. We see in Figure 2 that the divergence decreases globally for all latent dimensions even if for lower values of p the divergence may temporarily increase. As expected the convergence is faster for higher values of p.

Table 2 details the memory cost for higher dimensions.

Remark 1 (Uncertainty quantification and calibration). Using a lower value for p may lead to a poor estimation of the uncertainty. In particular, the estimated covariance may not be consistent with the empirical covariance. Variational inference is known to underestimate the target distribution but in critical applications, we may often prefer to overestimate the distribution. A common approach in Kalman filtering is to add a process noise to recover consistency. In our setting, it will be equivalent to adding a covariance matrix of noise Q t in the L-RVGA update (23). The covariance update will then rewrite P t|t-1 = P t-1 + Q t and P -1

t = P -1 t|t-1 + x t x T t .
To formulate this update in terms of the factor analysis parameters W t and ψ t , we can assume Q t = η t P t-1 which corresponds to a fading memory filter. The parameter(s) of the covariance of noise may be learned online using a variational approach or an adaptive filter that adapts the covariance based on the prediction error (see for instance Zhang et al. ( 2020)). Table 2: Memory requirements test for linear regression with d = 10 6 and variable p: The L-RVGA has been executed on a laptop (Intel core i7 at 2.3 GHz on CPU) in dimension one million, the memory cost for the full Kalman has been estimated as it scales quadratically in d and did not fit into the memory, so that execution time is not available for it. A reduced number of samples N has been considered in this experiment since they do not influence the memory cost. When p is low, the cost in memory as well as the time per iteration (with 3 inner loops for the recursive EM) is drastically reduced.

Application to Logistic regression

We now apply the method to derive a limited-memory version of a logistic regression problem which is a particular case of the generalized linear model described in Section 4.1. In this model, we seek to learn the parameter θ encoding the hyperplane from N examples (x t , y t ). We consider the Bayesian setting where θ ∼ N (µ, P ) and the observations are generated from y t = σ(µ T x t ), where σ denotes here the logistic function. The estimation of µ t and P t are given by the RVGA updates for general linear model ( 29)-( 30). These updates involve the expectation of the gradient and the Hessian of the logistic loss which can be approximated as follow:

E θ [∇ θ t (θ)] = -y t x t + E θ [σ(x T t θ)]x t ≈ -y t x t + σ(k t x T t µ t ) (49) E θ [∇ 2 θ t (θ)] = E θ [σ(x T t θ)(1 -σ(x T t θ))]x t x T t ≈ k t σ (k t x T t µ t )x t x T t ( 50
)
where

k t = β x T t P t x t + β 2 and β = 8 π . ( 51
)
These equations were derived in our previous work (see [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF], Section 4.1) and come from the approximation of the logistic function σ with the inverse probit function φ (Barber and Bishop, 1998a)

:σ(x) ≈ φ( 1 β x) = 1 2 (1 + erf ( x √ 2β
)). We can then rewrite ( 29)-( 30) as:

µ t = µ t-1 + P t-1 x t (y t -σ(k t x T t µ t )) (52)
P -1 t = P -1 t-1 + k t σ (k t x T t µ t )x t x T t . ( 53
)
The scheme is now implicit owing only to the two scalar parameters ν = x T t P t x t and α = x T t µ t which can efficiently be computed by solving a scalar fixed point equation with a Newton solver, see Lambert et al. 2021, Section 4.2. We can in turn apply our recursive EM algorithm on the input X t := x t k t σ (k t x T t µ t ) for highdimensional logistic regression. To assess the performance, we generate N synthetic pairs (x t , y t ) where the inputs are generated from the same Gaussian distribution as in the linear regression case. The inputs are normalized in mean and the prior is set to σ 0 = 4 which corresponds to a sharp posterior distribution. In Figure 3 we plot the KL divergence between our current Gaussian estimate N (µ t , P t ) and the true posterior, that is KL (N (µ t , P t )||p(θ|y 1 , x 1 , . . . y N , x N )). In logistic regression, we have a simple expression for the posterior at each θ, up to a normalizing constant. As the former left KL divergence is an expectation under N (µ t , P t ), it may be approximated via Monte-Carlo sampling. This offers a way to perform relative comparisons between the algorithms in terms of divergence with respect to the true (unnormalized) posterior. For more details see [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF] We compare the results for different values of the parameter p with the Laplace approximation which gives a batch approximation of µ and P . The L-RVGA converges to the batch Laplace approximation and may even yield lower divergence. This is because the covariance given by the Laplace approximation spills out the true posterior to regions of very low probability whereas the L-RVGA avoids them. Contrary to the linear case, even with p = 1, we obtain very good results. The memory requirements in dimension one million yield identical costs to those reported in Table 2 for linear regression.

General nonlinear method evaluation

We address in this Section the general nonlinear L-RVGA, based on the approximations detailed in Section 4.2. To assess the method, we choose to apply the general nonlinear method to logistic regression, again, as we have closed-form expressions that may serve as ground truth.

As the generalized Gauss-Newton approximation, presented in Section (4.2.3), is always exact in the logistic case, only the other approximations will be evaluated: sampling to approximate the expectation, extra-gradient to approximate the implicit scheme, and factor analysis to approximate the covariance matrix.

As logistic regression is based on a generalized linear model, µ t and P t are given by the RVGA updates for generalized linear models ( 29)-( 30). These updates involve the expectation of the The estimated mean is roughly aligned with the maximum a posteriori (MAP). Only one inner loop proves sufficient to ensure convergence. The prior is set to be P 0 = σ 0 I d with σ 0 = 4 and the algorithm uses it for initialization. The KL divergence is computed with sampling and is unnormalized.

gradient and the Hessian of the logistic loss:

E θ [∇ θ t (θ)] = -y t x t + E θ [σ(x T t θ)]x t (54) E θ [∇ 2 θ t (θ)] = E θ [σ(x T t θ)(1 -σ(x T t θ))]x t x T t . (55)
Those expectations could be easily computed using equations ( 49)-( 50). For comparison purposes, we may approximate those expectations using the general sampling procedure of Section 4.2.2, that is,

E θ [σ(x T t θ)] ≈ 1 K K i=1 σ(x T t θ i ) (56) E θ [σ(x T t θ)(1 -σ(x T t θ))] ≈ 1 K K i=1 σ(x T t θ i )(1 -σ(x T t θ i )),
where the samples θ i are drawn from the Gaussian distribution N (µ, (Ψ + W W T ) -1 ) exactly using ensemble sampling as described in Section 4.2.2. Table 3 shows that the method approximates well the closed-form expression used in Section 5.3. We then combined sampling with the extra-gradient method, where we skipped the extra covariance update, as described in Section 4.2.1 since we observed it helped convergence. We found few samples may be sufficient to provide good convergence of the algorithm in terms of KL divergence, as shown in Figure 4.

Method

Tr(W W T + Ψ) ] is performed with 10 samples. The closedform expression with the inverse probit approximation is considered as the baseline. We also show how the trace of the matrix (W W T +Ψ) -1 is approximated. The Cholesky sampling uses the square root decomposition of (W W T + Ψ) -1 . The ensemble sampling is our method described in Section 4.2.2. It achieves similar results as Cholesky sampling which is not memory-limited since it uses the full matrix and inverts it. The inputs have been normalized on average. We compare the L-RVGA where the expectations are computed analytically (baseline) and the L-RVGA where the expectations are approximated with sampling for 1, 10, and 100 samples. All the algorithms use a Mirror prox (where we have skipped the extra covariance update) and use the same factor analysis approximation p = 10. We see K = 10 samples are sufficient even in the more difficult case to converge to the batch Laplace KL.

In the general nonlinear case, the L-RVGA algorithm can be advantageously used when the observations are in large dimensions and also in large numbers. A batch approach cannot be used on those problems where one can only afford multiple passes, or even a single pass, over the data set. Standard variational Bayesian methods based on a mean-field approach are not suitable for this class of problem where all observations are required at each iteration. [START_REF] Opper | A Bayesian Approach to On-Line Learning[END_REF] consider an online variational method using only one observation at each update. However, they used the right-side KL (moment matching) and lose the guarantee to increase an (evidence) lower bound.

Approaches based on the stochastic maximization of the evidence lower bound Ong et al. ( 2018), [START_REF] Mishkin | Slang: Fast structured covariance approximations for bayesian deep learning with natural gradient[END_REF] tackle the factor analysis problem in an online way using an adaptive step descent in contrast to the L-RVGA which doesn't use a step parameter. These algorithms may be more adaptable to the smoothness of the problem using step-tuning. We anticipate that adding process noise and learning it online can provide adaptability, in particular, it can compensate for all the approximation errors we added in the process.

Sources:

The sources of the code are available on Github on the following repository: https: //github.com/marc-h-lambert/L-RVGA.

Conclusion

We have developed a new second-order algorithm, called L-RVGA, for online variational inference which scales to both large data sets and high dimensions. This algorithm is based on a twostage variational problem that combines a variational Gaussian approximation followed by a factor analysis approximation of the inverse of the covariance matrix. L-RVGA is able to estimate the mean and the covariance of the distribution of the latent parameters in a memory-efficient and parameter-free way and with only one pass through the data. We have tested it on linear and logistic regression problems and shown how to extend it to more general nonlinear problems with extra-gradients, memory-efficient sampling, and the generalized Gauss-Newton approximation. To build our generic algorithm we have introduced two new tools: a recursive EM algorithm for factor analysis which is parameter-free and faster than the online EM in this context; and a sampler for Gaussian distribution with a structured precision matrix.

Beyond variational inference, we anticipate the L-RVGA may prove useful for stochastic optimization keeping only the mean and using the covariance to build an adaptive learning step. This version could be compared to state-of-the-art algorithms in limited memory optimization such as Adagrad or even L-BFGS. We will investigate this direction in future work and will extend our algorithm using adaptive filtering techniques to further increase its performance. We believe that a deeper connection between the communities of stochastic optimization and Kalman filtering may bring new ideas to tackle nonlinear stochastic problems. Sam Roweis. EM algorithms for PCA and SPCA. Advances in Neural Information Processing Systems, 10, 1998.

Michael E. Tipping and Christopher M. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society., 61:611-622, 1999.

Lingyi Zhang, David Sidoti, Adam Bienkowski, Krishna Pattipati, Yaakov bar shalom, and David Kleinman. On the identification of noise covariances and adaptive kalman filtering: A new look at a 50 year-old problem. 05 2020.

A The recursive EM

A.1 Derivation of fixed point equation for factor analysis with EM

In this Section, we show that finding the factor analysis parameter with EM is equivalent to iterating through the fixed point equation ( 20). In the context of probabilistic principal component analysis Tipping and Bishop (1999) have highlighted (in their Appendix B) how the EM can be rewritten as a fixed point equation. The following proof is an extension of this result to the factor analysis case.

In factor analysis, we want to approximate the empirical covariance matrix 1 K K i=1 v i v T i with a "diagonal + low rank" structure Ψ + W W T . In the classical EM approach, we introduce latent variables z i such that our samples can be rewritten

v i = W z i + ε i where ε i ∼ N (0, Ψ). We note V = v 1 • • • v K the sample matrix of size d × K such that 1 K K i=1 v i v T i = 1 K V V T and Z = z 1 • • • z K the latent matrix of size p × K.
At the expectation step (E-step), the latent variables z i are estimated conditionally on the observation and the current parameter estimates using the conditioning formula:

p(z i |v i , W, Ψ) = N (M -1 W T Ψ -1 v i , M -1 ) where M = I p + W T Ψ -1 W. ( 57
)
At the maximization step (M-step), the latent variables are assumed fixed and the parameters are adjusted to maximize the total likelihood defined by:

L(W, Ψ) = log p(V, Z|W, Ψ) = K i=1 log p(v i |z i , W, Ψ) + K i=1 log p(z i ) (58) = - K 2 (v i -W z i ) T Ψ -1 (v i -W z i ) - K 2 log det Ψ - K 2 log(2π) + K i=1 log p(z i ). ( 59
)
The derivative of the expectation of L with respect to W gives:

∂ ∂W E Z∼p(Z|V ) [L(W, Ψ)] = - K i=1 Ψ -1 W E[z i z T i |v i ] - K i=1 Ψ -1 v i E[z i |v i ] T = 0, (60)
and we obtain the optimum

W (n) = K i=1 v i E[z i |v i ] T ( K i=1 E[z t z T t |v i ]) -1
, where the (n) stands for "new".

Taking the derivative of the expectation of L with respect to Ψ -1 and keeping only the diagonal terms yields:

Ψ (n) = diag( 1 K K i=1 v i v T i + 1 K K i=1 W (n) E[z i z T i |v i ]W (n) T - 2 K K i=1 W (n) E[z i |v i ]v T i ) (61) = diag( 1 K K i=1 v i T i -W (n) 1 K K i=1 E[z i |v i ]v T i ), (62)
where we have replaced W (n) T by its optimal value in equation ( 61).

And finally, the EM updates give:

(63)

E-Step: E[z i |v i ] = M -1 W T Ψ -1 v i E[z i z T i |v i ] = M -1 + E[z i |v i ]E[z i |v i ] T . (64)
M-Step:

W (n) = K i=1 v i E[z i |v i ] T ( K i=1 E[z i z T i |v i ]) -1 Ψ (n) = diag( 1 K K i=1 v i v T i -W (n) 1 K K i=1 E[z i |v i ]v T i ). (65)
The EM updates ( 64) and ( 65) can be rewritten in batch form using the matrix notation V and Z:

(66) E-Step:

E[Z|V ] = M -1 W T Ψ -1 V K i=1 E[z i z T i |v i ] := E[ZZ T |V ] = KM -1 + E[Z|V ]E[Z|V ] T = K(M -1 + M -1 W T Ψ -1 S K Ψ -1 W M -1 ). (67)
M-Step:

W (n) = V E[Z|V ] T E[ZZ T |V ] -1 Ψ (n) = diag( 1 K V V T -W (n) 1 K E[Z|V ]V T ). (68)
And finally the E-step and M-step can be fused to form a fixed-point equation:

W (n) = V E[Z|V ] T E[ZZ T |V ] -1 = S K Ψ -1 W M -1 (M -1 + M -1 W T Ψ -1 S K Ψ -1 W M -1 ) -1 = S K Ψ -1 W (I p + M -1 W T Ψ -1 S K Ψ -1 W ) -1 Ψ (n) = diag( 1 K VV T -W (n) 1 K E[Z|V]V T ) diag(S K -W (n) M -1 W T Ψ -1 S K ).
which is the fixed point equation given in (20).

B The fixed point EM is equivalent to the MLE fixed point

In this Section we show new results concerning the equivalence of the fixed points for the marginal likelihood (MLE algorithm) and the total likelihood (EM algorithm) for the particular case of factor analysis. We consider here the batch factor analysis problem. This result may not be applied to our RVGA algorithm but to each inner loop which is guaranteed to increase the one sample likelihood.

The marginal likelihood and the total likelihood are both related as follows:

max W,Ψ log p(v 1 , . . . , v N |W, Ψ) where S N = 1 N N i=1 v i v T i = P -1 (MLE) . (69)
≤ max

W,Ψ E z [log p(v 1 , . . . , v N , z|W, Ψ)] with p(v | z) ∼ N (W z, Ψ), p(z) ∼ N (0, I p ) (EM) . ( 70
)
The advantage of the EM approach is that the total likelihood appearing in ( 70) is guaranteed to increase at each step, i.e., the algorithm is stable. The EM algorithm may not always converge to the maximum of the marginal likelihood appearing in (69) but converges to a stationary point which turns out to be also a stationary point for the maximum likelihood. This fact was already highlighted by [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] and [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF] but we specify the result in the case of factor analysis in an algebraic way. We first write the maximum likelihood as a fixed point, then show an equivalence between the two fixed points. This result is finally used to show an equivalence between different eigenvalues decomposition algorithms.

B.1 Factor analysis with maximum likelihood (MLE) as a fixed point.

In this section, we write the maximum likelihood over the parameters of the factor analysis in the form of a fixed-point equation. This is a well-known result derived from chapter 21.2 of the book of [START_REF] Barber | Bayesian Reasoning and Machine Learning[END_REF] but we want to use this result to make a connection with the fixed point obtained with EM. We want to approximate a matrix S with a matrix W W T +Ψ by maximizing the following likelihood (where C is a constant):

max W,Ψ L(W, Ψ) = max W,Ψ - N 2 Tr((W W T + Ψ) -1 S N ) - N 2 log det(W W T + Ψ) + C. (71)
A necessary condition on the optimal solution is to zeroing the gradients:

∂L(W, Ψ) ∂W = (W W T + Ψ) -1 W T -(W W T + Ψ) -1 S N (W W T + Ψ) -1 W = 0 (72) ∂L(W, Ψ) ∂Ψ = (W W T + Ψ) -1 -(W W T + Ψ) -1 S N (W W T + Ψ) -1 = 0, (73)
leading to the fixed-point equations:

Fixed-point equations for the MLE (74)

W (n) = S N (W W T + Ψ) -1 W (75) Ψ = diag(S N -W (n) W (n) T ). (76)
B.2 Equivalence of fixed points and convergence

The following proposition shows the algebraic equivalence of the fixed points.

Proposition 1. The factor analysis parameters which maximize the marginal likelihood, that is, the solution to (69) satisfy the following fixed point equation:

W (n) = S N (W W T + Ψ) -1 W (77) Ψ = diag(S N -W (n) W (n) T ). ( 78
)
This fixed point equation is equivalent to the EM fixed-point equation:

W (n) = S N Ψ -1 W (I p + M -1 W T Ψ -1 S N Ψ -1 W ) -1 (79)
where:

M = I p + W T Ψ -1 W (80) Ψ = diag(S N -W (n) M -1 W T Ψ -1 S N ). (81)
As a consequence, the EM algorithm converges to a stationary point which is also a stationary point for the likelihood.

Proof. The proof for the equivalence of the fixed points is straightforward, the update for W can be rewritten as:

W = S N (W W T + Ψ) -1 W (82) = S N Ψ -1 W (I p + W T Ψ -1 W ) -1 (From the Woodbury formula) (83) = S N Ψ -1 W M -1 (84) = S N Ψ -1 W (I p + (S N Ψ -1 W M -1 ) T Ψ -1 W ) -1 (85) = S N Ψ -1 W (I p + M -1 W T Ψ -1 S N Ψ -1 W ) -1 , (86)
where we have replaced the term W T in (83) by its development in (84).

The update for Ψ can be rewritten as:

Ψ = diag(S N -WW T ) (87) = diag(S N -W (S N Ψ -1 W M -1 ) T ) (88)
where we have replaced the term W T in (87) by its development in (84), leading to

Ψ = diag(S N -W M -1 W T Ψ -1 S N ) (89) = diag(S N -W n M -1 W T Ψ -1 S N ). (90)

B.3 Relation to singular value decomposition

The fixed point equation from the maximum likelihood is solved using a singular value decomposition (SVD) of S N (see [START_REF] Barber | Bayesian Reasoning and Machine Learning[END_REF], chapter 21.2). The equivalence with the EM fixed point suggests that the EM make implicitly an SVD decomposition. It can be shown it is the case if we consider an asymptotically Probabilistic Principal Component analysis form (PPCA), ie Ψ = σI where we let tends the parameter σ to 0. This result was shown by Roweis (1998) for the fixed point EM with PPCA and is related to the MLE fixed point in the following Corollary.

Corollary 1 (of Prop. 1).

The factor analysis MLE fixed-point equation for PPCA:

W (n) = S N (W W T + σ 2 I d ) -1 W = S N W (σ 2 I p + W T W ) -1 (91)
converge for σ → 0 to a fixed point :

W (n) = S N W (W T W ) -1 , (92)
which corresponds to the power method:

w ← S N w ||w|| 2
(given here in vectorial form).

The factor analysis EM fixed-point for PPCA:

W n = S N W (σ 2 I p + (σ 2 I p + W T W ) -1 W T S N W ) -1 (94)
converge for σ → 0 to a fixed point :

W (n) = S N W (W T S N W ) -1 W T W, (95)
which is equivalent to the EM-PCA method (Roweis, 1998):

X = (W T W ) -1 W T Y (96) W (n) = Y X T (XX T ) -1 where Y is defined such that S N = Y Y T . ( 97
)
Proof. The proof is direct.

C Recursive factor analysis of the covariance matrix

The retained factorization W t W T t + Ψ t of the precision matrix P -1 t associated with the Bayesian parameter θ t is not conventional in factor analysis, which is usually applied to the empirical covariance matrix S t = 1 t t i=1 x i x T i of the inputs x t . In the linear case, both can be related as follows S t := 1 t P -1 t for t > 0, and at t = 0 we let S 0 = P -1 0 . S t can be expressed in a recursive way as:

S t = 1 t P -1 t = t -1 t P -1 t-1 t -1 + 1 t x t x T t = t -1 t S t-1 + 1 t x t x T t . (98)
The online EM algorithm [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF] considers the following fixed point equation with respect to the sufficient statistics S:

S = E vt E zt∼p(z|vt,θ * (S)) [S(z t , v t )] = T (S) (100)
where θ * (S) = arg max

θ∈Θ N t=1 < S(z t , v t ), φ(θ) > -F (θ), (101)
and solve T (S) -S using a stochastic root solver based on the Robbins-Monro algorithm [START_REF] Robbins | A stochastic approximation method[END_REF]. The sufficient statistics S and the optimal parameter θ are update online with an adaptive step γ t at each new incoming observations v t as follows:

S t = (1 -γ t )S t-1 + γ t E zt∼p(z|vt,θ t-1 ) [S(z t , v t )] (E-step) (102) θ t = arg max θ∈Θ < S t , φ(θ) > -F (θ) (M-step). (103)
In the case of factor analysis, the joint distribution is:

log p(v t , z t ) = - 1 2 Tr[Σ -1 S(v t , z t )] - 1 2 log det Σ + c, (104)
where the joint covariance is :

Σ = W W T + Ψ W W T I p , (105)
and the sufficient statistics are :

S(v t , z t ) = v t z t v t z t T = v t v T t v t z T t z t v T t z t z T t . (106)
The expectation of sufficient statistics gives

E zt∼p(z|vt,θ t-1 ) [S(v t , z t )] = v t v T t v t E[z T t |v t ] E[z t |v t ]v T t E[z t z T t |v t ] (107)
Rather than update the full matrix S(v t , z t ) we will update the blocks:

S 1.t = v t v T t , S 2.t = E[z t |v t ]v T t and S 3.t = E[z t z T t |v t ],
which are necessary to compute the M-step. Finally, using the same notation as in the previous Section, the online EM updates for factor analysis become:

Online EM E-step S 1.t = (1 -γ t )S 1.t-1 + γ t v t v T t S 2.t = (1 -γ t )S 2.t-1 + γ t E[z t |v t ]v T t = (1 -γ t )µ t-1 + γ t M -1 t-1 W T t-1 Ψ -1 t-1 v t v T t S 3.t = (1 -γ t )S 3.t-1 + γ t E[z t z T t |v t ] = (1 -γ t )S 3.t-1 + γ t (M -1 t-1 + M -1 t-1 W T t-1 Ψ -1 t-1 v t v T t Ψ -1 t-1 W t-1 M -1 t-1 ) M-step W t = S T 2.t S -1 3.t Ψ t = diag(S 1.t ) -diag(W t S 2.t ). (108)
To develop a limited memory version of this algorithm, we store only the diagonal of the high dimensional squared matrix S 1.t and update it as follows:

diag(S 1.t ) = (1 -γ t )diag(S 1.t-1 ) + γ t u t * u t , (109)
where x * y is a component wise operation giving x * y = diag(xy T ) for two vectors.

To choose the step size γ t we must satisfy the Robins Monro rules:

N t=1 γ t = ∞ N t=1 γ 2 t < ∞. (110)
We consider the following step recommended by [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF]:

γ 0 = 1 (111) γ t = 1 t 0.6 ∀t > 0. (112)
Finally, in post-processing, we use a Polyak-Ruppert halfway averaging to improve the convergence as recommended by [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF]:

∀t > N/2 st t = t -N/2 > 0 do : Wt = t -1 t W t-1 + 1 t W t (113) Ψt = t -1 t Ψ t-1 + 1 t Ψ t . (114)
F The outer product approximation in the general case.

We show in this Section the relation between the generalized Gauss-Newton approximation and the Fisher matrix. The proof proposed here comes from Ollivier (2018)[Appendix A]:

Proof. Let's p be an exponential family such that:

p(y|θ) = m(y t ) exp(η T T (y) -A(η)), (115)
where T (y) is the sufficient statistics, A is the log partition function which satisfies ∇ 2 η A(η) = Cov(y|θ) and ∇ η A(η) = m = E[(y|θ)] and finally η is the natural parameter which depends on θ through a function h, ie η = h(θ). Using twice the chain rules, the second derivative of the negative likelihood of p writes :

- ∂ 2 ln p(y|θ) ∂θ 2 = - ∂h ∂θ ∂ 2 ln p(y|θ) ∂h 2 ∂h ∂θ T - ∂ 2 h ∂θ 2 ∂ ln p(y|θ) ∂h (116) = ∂h ∂θ Cov(y|θ) ∂h ∂θ T - ∂ 2 h ∂θ 2 (T (y) -m). (117)
Taking the expectation under y on both sides, we obtain directly the relation:

E y [- ∂ 2 ln p(y/η) ∂θ 2 ] = F(θ) = ∂η ∂θ Cov(y|θ) ∂η ∂θ T . (118)
And finally:

E θ [F(θ)] = E θ [ ∂h ∂θ Cov(y|θ) ∂h ∂θ T ]. (119)
Now for a generalized linear model such that h(θ) = x T t θ the GGN approximation is exact since

∂ 2 h
∂θ 2 = 0, this completes the proof.

G Mirror prox

In this section, we show that the iterated scheme defined equation ( 36) is equivalent to a Mirror prox update [START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF]. We recall first the connection between the recursive variational scheme and the Mirror descent using the results initially derived by Khan et al. (2018) for the batch variational approach and extended by [START_REF] Lambert | The recursive variational gaussian approximation (r-vga)[END_REF] for the recursive variational approach.

Considering an exponential family q η of natural parameter η, mean parameter m and a strictly convex log partition function F such that q η (θ) = h(θ) exp(< η, θ > -F (η)), the recursive variational approximation problem between a target distribution q η and the one-sample posterior p(θ|y t ) ∝ p(y t |θ)q η t-1 (θ) writes: arg min ηt KL(q ηt (θ)|p(θ|y t ))

(120)

=arg min ηt E qη t [-log p(y t |θ)] + B F (η t-1 , η t ), (121)
where B F is the Bregman divergence associated with the strictly convex log partition function F . The critical point must satisfy:

∇ ηt E qη t [-log p(y t |θ)] + (η tη t-1 )∇ 2 F (η t ) = 0, (122) which gives the following implicit fixed point equation on the natural parameter: This implicit scheme can be approximated using the Mirror prox algorithm [START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF] with a step size one:

η t = η t-1 + (∇ 2 F (η t )) -1 ∇ η E qη [log p(y t |θ)](η t ) (123) = η t-1 + ∇ m E qm
Mirror prox ηt = ∇F * ( mt ) = ∇F * (m t-1 ) -∇ m f (m t-1 ) mt = ∇F (η t )

η t = ∇F * (m t ) = ∇F * (m t-1 ) -∇ m f ( mt ) m t = ∇F (η t ).
In fact, we have rather computed here a stochastic version of Mirror prox. The stochastic version of Mirror prox may inherit the good convergence properties of the original batch version if the function f is convex and the step is adaptive [START_REF] Bach | A universal algorithm for variational inequalities adaptive to smoothness and noise[END_REF]. Here the setting is different: our function f is not jointly convex in µ and P and we consider a constant step of size one. However, we have found the stochastic version behaves empirically well.

In the case where q is a multivariate Gaussian distribution, the mean and natural parameters are given by η = ∇F * (m) = η 1 = P -1 µ η 2 = -1 2 P -1 and m = ∇F (η) = m 1 = µ m 2 = P + µµ T .

30

The gradient with respect to the mean parameters m 1 , m 2 can be expressed as the gradient with respect to the source parameters µ, P using the chain rule:

∂f
If the expectations are replaced with a linearization around the last estimated, these updates are also equivalent to the extended iterated Kalman filter scheme [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF].

Applying the mirror-prox scheme to our logistic regression problem 5.3 without factor analysis, we see that the Gaussian well approximates the logistic posterior in figure 5. However, when we combine mirror-prox with factor analysis, the extra covariance update can make the mirror-prox scheme unstable. We have observed it is then preferable to skip the extra covariance update 139, i.e. using: P-1 t = P -1 t-1 -E θ∼N (µ t-1 ,P The contour lines of the true posterior are displayed in green. The batch Laplace ellipsoid is shown in a dashed line. We compare the mirror-prox updates 137-140 (right column) with other variants of the updates: explicit using only 137-138 (second column), implicit using 52-53 (third column).

The natural gradient (left column) corresponds to a variant where the models are linearized. The mirror-prox schemes clearly better approximate the Bayesian logistic posterior.

KLKLFigure 1 :

 1 Figure1: Factorization of a synthetic large scale covariance matrix of dimension d = 100 (upper left plot) and d = 1000 (upper right plot). The lower plots show additional results on two covariance matrices computed from the Madelon and Breast Cancer datasets. We show in dash green line our recursive EM algorithm, which is parameter-free, and in the red plain line, the online EM tuned with a learning rate γ = 1 t 0.6 (see Appendix E for implementation details). We show also the batch EM, in the solid blue line. The online and batch methods should be compared with respect to the number of observations processed. Batch estimation requires processing all the data several times, although more efficiently since 5 steps are enough to converge. The online variants use only a single passe. With 3 inner loops in dimension d = 100 (upper left plot), the recursive EM is clearly superior to the online EM. In the higher dimension (upper right plot), we achieve a similar loss to the online EM using only 1 inner loop in the recursive EM. Moreover, the recursive EM gives better results than the other algorithms on real datasets as illustrated in the lower left plot and lower right plot.

Figure 2 :

 2 Figure2: Sensitivity to the latent dimension p. We show the error as the KL divergence between the true posterior distribution and the L-RVGA posterior. The full Kalman error converges exactly to 0, as the Kalman filter exactly finds the posterior. When p is very low, oscillations occur.

Figure 3 :

 3 Figure3: Sensitivity to the latent dimension p. The L-RVGA algorithms may provide a lower divergence than the batch Laplace approximation in the sense of KL divergence to the true posterior. The estimated mean is roughly aligned with the maximum a posteriori (MAP). Only one inner loop proves sufficient to ensure convergence. The prior is set to be P 0 = σ 0 I d with σ 0 = 4 and the algorithm uses it for initialization. The KL divergence is computed with sampling and is unnormalized.

Figure 4 :

 4 Figure4: Approximation of the logistic posterior for different deviations on the priors σ 0 = 1, 2 and 3. The inputs have been normalized on average. We compare the L-RVGA where the expectations are computed analytically (baseline) and the L-RVGA where the expectations are approximated with sampling for 1, 10, and 100 samples. All the algorithms use a Mirror prox (where we have skipped the extra covariance update) and use the same factor analysis approximation p = 10. We see K = 10 samples are sufficient even in the more difficult case to converge to the batch Laplace KL.

  [log p(y t |θ)](m t ). (124) If we consider the function f (m) = E qm [-log p(y t |θ)] and use the relation η = ∇F * (m), this update can be interpreted as an implicit version of a Mirror descent on f with step size one: ∇F * (m t ) = ∇F * (m t-1 ) -∇ m f (m t ) (dual descent with step 1) (125) m t = ∇F (η t ). ( projection) (126)

Figure 5 :

 5 Figure 5: Gaussian approximation of the Bayesian logistic posterior with a sharp prior σ 0 = 10 for different algorithms. The confidence ellipsoids of the Gaussians at the final time are shown in red.The contour lines of the true posterior are displayed in green. The batch Laplace ellipsoid is shown in a dashed line. We compare the mirror-prox updates 137-140 (right column) with other variants of the updates: explicit using only 137-138 (second column), implicit using 52-53 (third column). The natural gradient (left column) corresponds to a variant where the models are linearized. The mirror-prox schemes clearly better approximate the Bayesian logistic posterior.

Table 1 :

 1 Memory test for large scale matrix factorization with d = 10 6 , and p = 100: Test dedicated for memory requirement, the recursive and online EM have been executed on a laptop (Intel core i7 at 2.3 GHz on CPU) in dimension one million, the memory cost for the batch EM have been estimated as it scales quadratically in d and did not fit the memory, so no time is available for it. A reduced number of samples N have been considered in this experiment since they do not influence the memory cost. Recursive EM requires more operations by iteration than online EM but consumes less memory than online EM and much less than batch EM.

Table 3 :

 3 Test for the ensemble sampling approximation in dimension d = 10000 and p = 10: Approximation of the expectation E θ [σ(x T t θ)

	-1 E θ [σ(x T t θ)]

  ∇F * (m t ) = ∇F * (m t-1 ) -∇ m f (m t-1 ) (129) ⇐⇒ η t = η t-1 -∇ m f (m t-1 ),(130)become, if we write f (m(µ, P )) = E θ∼N (µ,P ) [log p(y t |θ)]:,P t-1 -2 ∂f ∂P | µ t-1 ,P t-1 µ t-1 ∂f ∂P | µ t-1 ,P t-1Rearranging terms and applying two times the update as in Mirror prox gives the iterated scheme defined in equation (36): -E θ∼N (µ t-1 ,P t-1 ) [∇ 2 θ log p(y t |θ)](137)μt = µ t-1 + Pt E θ∼N (µ t-1 ,P t-1 ) [∇ θ log p(y t |θ)] -E θ∼N (μ t , Pt ) [∇ 2 θ log p(y t |θ)](139) µ t = µ t-1 + P t E θ∼N (μ t , Pt ) [∇ θ log p(y t |θ)].

		∂m 1	=	∂f ∂µ	∂µ ∂m 1	+	∂f ∂P	∂P ∂m 1	=	∂f ∂µ	-2	∂f ∂P	µ	(127)
		∂f ∂m 2	=	∂f ∂µ	∂µ ∂m 2	+	∂f ∂P	∂P ∂m 2	=	∂f ∂P		(128)
	A step of the Mirror prox update:
		P -1 t µ t -1 2 P -1 t	=		P -1 t-1 µ t-1 -1 t-1 2 P -1		+	∂f ∂µ | µ t-1 (131)
	⇐⇒	P -1 t µ t -1 2 P -1 t	=		(P -1 t-1 -2 ∂f ∂P | µ t-1 ,P t-1 µ t-1 -1 2 P -1 t-1	+	∂f ∂µ | µ t-1 ,P t-1 ∂f ∂P | µ t-1 ,P t-1	(132)
	⇐⇒ ⇐⇒	P -1 t µ t -1 2 P -1 t P -1 t µ t -1 2 P -1 t	= =		P -1 t µ t-1 2 P -1 t-1 -1 P -1 t µ t-1 -1 2 P -1 t-1								(133)
														(135)
				∇ P N (θ|µ, P ) =	1 2	∇ 2 θ N (θ|µ, P ).	(136)
		P-1 t = P -1 t-1 (138)
		P t	-1 = P -1 t-1

+ ∂f ∂µ | µ t-1 ,P t-1 ∂f ∂P | µ t-1 ,P t-1 + -E θ∼N (µ t-1 ,P t-1 ) [∇ θ log p(y t |θ)] 1 2 E θ∼N (µ t-1 ,P t-1 ) [∇ 2 θ log p(y t |θ)] .

(134)

The last derivation (134) comes from the Bonnet & Price formulas (Lin et al., 2019):

∇ µ N (θ|µ, P ) = -∇ θ N (θ|µ, P )

  t-1 ) [∇ 2 θ log p(y t |θ)] (141) μt = µ t-1 + Pt E θ∼N (µ t-1 ,P t-1 ) [∇ θ log p(y t |θ)] (142) µ t = µ t-1 + Pt E θ∼N (μ t , Pt ) [∇ θ log p(y t |θ)].

Acknowledgements

This work was funded by the French Defence procurement agency (DGA) and by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001(PRAIRIE 3IA Institute). We also acknowledge support from the European Research Council (grant SEQUOIA 724063).

If we let S 0 be null, we obtain exactly 1 t t i=1 x i x T i , otherwise, we obtain a regularized version of the empirical covariance. The corresponding recursive factor analysis form is:

which fits into the problem (27). It may thus be addressed through Algorithm 1 with α t = (t -1)/t and β t = 1/t. We may guess the value of S 0 to initialize the procedure and make it more stable.

Observing that Tr 1

we see that the trace of the unknown covariance must match the expectation of the square norm of inputs. This expectation may be estimated on a batch of size M . We use this property to compute S 0 as follows: get a batch formed by the first M data inputs x 1 , . . . , x M , and let S 0 = 1

. This is equivalent to normalizing the inputs in mean and set σ 0 = √ d.

D Initialization and prior information

Given a prior covariance P 0 on θ, we want to initialize W 0 and Ψ 0 such that W 0 W T 0 +Ψ 0 ≈ P -1 0 . We must suppose P 0 is sparse enough to fit in memory and has an inverse which can be approximated by W 0 W T 0 + Ψ 0 for example using the Woodbury formulas. For the experiments, we consider an isotropic initial covariance P 0 = σ 2 0 I d . We then compute W 0 and Ψ 0 such that

I d and W 0 = 0 d×p would make the algorithm run into problems as W 0 = 0 d×p is a stationary point of the fixed-point equation ( 20).

We use the following rule Ψ 0 = ψ 0 I d where ψ 0 > 0 is a scalar and generate W 0 as a d × p matrix whose columns are the vectors u 1 0 , u 2 0 , . . . , u p 0 independently drawn from an isotropic Gaussian distribution in R d and which have been normalized so that ∀k : ||u k 0 || = w 0 . We then let:

T so that we have:

Remark 2. This initialization can be extended to the case where P 0 is a diagonal matrix P 0 = Diag(σ 2 1 , . . . , σ 2 d ) where σ 2 i represents now a variance on the i th coordinate.

E Derivation of the online EM algorithm for factor analysis

In this section, we derive the online EM algorithm [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF] for the factor analysis problem to compare it to our Recursive EM algorithm. We use the same notation as in the previous section where v 1 , . . . , v N are our N observations in dimension d and z 1 , . . . , z N are our latent variables in dimension p. Moreover, we suppose that each couple of variables (z t , v t ) are independent and belong to an exponential family given by log p(z t , v t |θ) =< S(z t , v t ), φ(θ) > -F (θ), where F is the log partition function, φ is a function which map the natural parameter and S(z t , v t ) are the sufficient statistics.