Online EM Algorithm for Latent Data Models Olivier Capp´ e & Eric Moulines LTCI, TELECOM ParisTech, CNRS. 46 rue Barrault, 75013 Paris, France. Abstract In this contribution, we propose a generic online (also sometimes called adaptive or recursive) version of the Expectation-Maximisation (EM) algorithm applicable to latent variable models of independent observations. Compared to the algorithm of Titterington (1984), this approach is more directly connected to the usual EM algorithm and does not rely on integration with respect to the complete data distribution. The resulting algorithm is usually simpler and is shown to achieve convergence to the stationary points of the Kullback-Leibler divergence between the marginal distribution of the observation and the model distribution at the optimal rate, i.e., that of the hal-00201327, version 1 - 27 Dec 2007 maximum likelihood estimator. In addition, the proposed approach is also suitable for conditional (or regression) models, as illustrated in the case of the mixture of linear regressions model. Keywords: Latent data models, Expectation-Maximisation, adaptive algorithms, online estima- tion, stochastic approximation, Polyak-Ruppert averaging, mixture of regressions. 1 Introduction The EM (Expectation-Maximisation) algorithm (Dempster et al., 1977) is a popular tool for maximum- likelihood (or maximum a posteriori) estimation. The common strand to problems where this ap- proach is applicable is a notion of incomplete data, which includes the conventional sense of missing data but is much broader than that. The EM algorithm demonstrates its strength in situations where some hypothetical experiments yields complete data that are related to the parameters more conveniently than the measurements are. Problems where the EM algorithm has proven to be useful include, among many others, mixture of densities (Titterington et al., 1985), censored data models (Tanner, 1993), etc. The EM algorithm has several appealing properties. Because it relies on com- plete data computations, it is generally simple to implement: at each iteration, (i) the so-called E-step only involves taking expectation over the conditional distribution of the latent data given the obser- vations and (ii) the M-step is analogous to complete data weighted maximum-likelihood estimation. Moreover, (iii) the EM algorithm naturally is an ascent algorithm, in the sense that it increases the (observed) likelihood at each iteration. Finally under some mild additional conditions, (iv) the EM algorithm may be shown to converge to a stationary point (i.e., a point where the gradient vanishes) of the log-likelihood (Wu, 1983). Note that convergence to the maximum likelihood estimator cannot in general be guaranteed due to possible presence of multiple stationary points. When processing large data sets or data streams however, the EM algorithm becomes impractical due to the requirement that the whole data be available at each iteration of the algorithm. For this reason, there has been a strong interest for online variants of the EM which make it possible to estimate the parameters of a latent data model without storing the data. In this work, we consider online algorithms for latent data models with independent observations. The dominant approach (see also Section 2.2 below) to online EM-like estimation follows the method proposed by Titterington (1984) which consists in using a stochastic approximation algorithm, where the parameters are updated after each new observation using the gradient of the incomplete data likelihood weighted by the complete 1 data Fisher information matrix. This approach has been used, with some variations, in many different applications (see, e.g., Chung and B¨ ohme, 2005; Liu et al., 2006); a proof of convergence was given by Wang and Zhao (2006). In this contribution, we propose a new online EM algorithm that sticks more closely to the principles of the original (batch-mode) EM algorithm. In particular, each iteration of the proposed algorithm is decomposed into two steps, were the first one is a stochastic approximation version of the E-step aimed at incorporating the information brought by the newly available observation, and, the second step consists in the maximisation program that appears in the M-step of the traditional EM algorithm. In addition, the proposed algorithm does not rely on the complete data information matrix, which has two important consequences: firstly, from a practical point of view, the evaluation and inversion of the information matrix is no longer required, secondly, the convergence of the procedure does not rely on the implicit assumption that the model is well-specified, that is, that the data under consideration is actually generated by the model, for some unknown value of the parameter. As a consequence, and in contrast to previous work, we provide an analysis of the proposed algorithm also for the case where the observations are not assumed to follow the fitted statistical model. This consideration is particularly relevant in the case of conditional missing data models, a simple case of which is used as an illustration of the proposed online EM algorithm. Finally, it is shown that, with the additional use of Polyak-Ruppert averaging, the proposed approach converges to the stationary points of the limiting normalised log-likelihood criterion (i.e., the Kullback-Leibler divergence between the marginal density of the observations and the model pdf) at a rate which is optimal. The paper is organised as follows: In Section 2, we review the basics of the EM and associated algorithms and introduce the proposed approach. The connections with other existing methods are discussed at the end of Section 2.3 and a simple example of application is described in Section 2.4. Convergence results are stated in Section 3, first in term of consistency (Section 3.1) and then of con- vergence rate (Section 3.2), with the corresponding proofs given in Appendix A. Finally in Section 4, the performance of this approach is illustrated in the context of mixture of linear regressions. 2 Algorithm Derivation 2.1 EM Basics In this section, we review the key properties of the EM algorithm as introduced by Dempster et al. (1977). The latent variable statistical model postulates the existence of a non-observable or latent variable X distributed under f (x; ?) where {f (x; ?), ? ? ?} denotes a parametric family of probability density functions indexed by a parameter value ? ? ? ? Rd? . The observation Y is then viewed as a deterministic function of X. This latent variable mechanism provides a unified framework for situa- tions which includes missing data, censored observations, noisily observed data, . . . (Dempster et al., 1977). We will denote by g(y; ?) the (observed) likelihood function induced by the latent data model. In addition, the notations E? [·] and E? [·|Y ] will be used to denote, respectively, the expectation and conditional expectation under the model parameterised by ?. Likewise, let ? denote the probability density function of the observation Y , where we stress again that we do not restrict ourselves to the case where ?(·) = g(y; ? ? ), for an unknown value ? ? of the parameter. The notations P? and E? will be used to denote probability and expectation under the actual distribution of the observation. def Given n independent observations Y1:n = (Y1 , . . . , Yn ), the maximum likelihood estimator is defined as ??n def = argmax??? n?1 log g(Y1:n ; ?), where n def log g(Y1:n ; ?) = log g(Yi ; ?) . (1) i=1 2 At this point of the exposition, the fact that we maximise the normalised log-likelihood rather than the log-likelihood is introduced only to ease the transition to the online setting, in which n increases when new observations become available. The EM algorithm is an iterative optimisation algorithm that maximises the above (normalised) log-likelihood function despite the possibly complicated form of g resulting from the latent data model. Traditionally, each EM iteration is decomposed in two steps. The E-step consists in evaluating the conditional expectation n def?1 ?k (Y1:n ; ?) = n Q? ?k log f (Xi ; ?) Yi E? (2) i=1 where ??k is the current estimate of ?, after k iterations of the algorithm. In the M-step, the value of ? maximising Q? ? ?k (Y1:n ; ?) is found. This yields the new parameter estimate ?k+1 . This two step process is repeated until convergence. The essence of the EM algorithm is that increasing Q? ?k (Y1:n ; ?) forces an increase of the log-likelihood log g(Y1:n ; ?) (Dempster et al., 1977). For m : ? ? R a differentiable function, denote by ?? m = (?m/??1 , . . . , ?m/??d? )T the gradient of m. If m is twice differentiable, we denote by ?2 ? m the Hessian which is a d? × d? matrix whose 2 ?2m components are given by [?? m]i,j = ??i ??j , 1 ? i, j ? d? . Following Lange (1995), if the M-step of the EM algorithm is replaced by a Newton update, one obtains, assuming some regularity, the following recursion n ?1 ? ?k + ?k+1 J(Y1:n ; ? ?k+1 = ? ?k ) ? ?k ?? log f (Xi ; ?k ) Yi , E? (3) i=1 where ?k+1 is a step size (?k+1 = 1 correspond to the actual Newton update) and J(Y1:n ; ?) = n?1 n 2 i=1 J(Yi ; ?) with J(y; ?) = ?E? ?? log f (X; ?) Y = y . Note that due to the so-called Fisher?s identity, the gradient term indeed coincides with the (observed data) score function as E? ?? log f (X; ?) Y = ?? log g(Y ; ?) . (4) The algorithm in (3) can be shown to be locally equivalent to the EM algorithm at convergence (Lange, 1995). In practise, it is often required to adjust the step-size ?k+1 using line searches to ensure global convergence of the algorithm. In addition, J(Y1:n ; ?) is not necessarily a positive definite matrix or could be badly conditioned; therefore, some adjustment of the weight matrix J(Y1:n ; ?) may be necessary to avoid numerical problems. A well-known modification of the Newton recursion consists in replacing J(Y1:n ; ?) in (3) by the Fisher Information Matrix (FIM) associated to a complete observation, def I(?) = ?E? ?2 ? log f (X; ?) . (5) Under the mild assumption that the complete data model is regular, I(?) is guaranteed to be positive definite. This modified recursion, which is more robust, may be seen as an approximation of the scoring method (McLachlan and Krishnan, 1997), where the complete data FIM is used in place of the actual (observed) FIM def Iobs (?) = ?E? ?2 ? log g(Y ; ?) , (6) despite the fact that, in general, Iobs (?) and I(?) are different. I(?) usually also differs from J(Y1:n ; ?), as J(Y1:n ; ?) converges, as n tends to infinity, to def I? (?) = ?E? E? ?2 ? log f (X; ?) Y , (7) which doesn?t correspond to a Fisher information matrix in the complete data model, except when ? coincides with f (·; ?). In the particular case, however, where the complete data model belongs to a 3 canonical (or naturally parameterised) exponential family of distributions, J(Y1:n ; ?) coincides with I(?) and thus does not depend on ? anymore. Hence, except in some specific cases or if one assumes that the model is well-specified (i.e., ? = g(·; ? ? )), the convergence behaviour of the recursion in (3) will be different when I(?) is used instead of J(Y1:n ; ?). 2.2 Stochastic Gradient EM Algorithms Being able to perform online estimation means that the data must be run through only once, which is obviously not possible with the vanilla EM algorithm. To overcome this difficulty we consider in the sequel online algorithms which produce, at a fixed computational cost, an updated parameter estimate ??n for each new observation Yn . Note that in the online setting, the iteration index (which was previously denoted by k) is identical to the observation index n and we will use the latter when describing the algorithms. To our best knowledge, the first online parameter estimation procedure for latent data models is due to Titterington (1984) who proposed to use a stochastic approximation version of the modified gradient recursion: ?n+1 = ? ? ?n + ?n+1 I ?1 (? ?n )?? log g(Yn+1 ; ? ?n ) , (8) where {?n } is a decreasing sequence of positive step sizes. One may also consider using a stochastic approximation version of the original (Newton) recursion in (3): ?n+1 = ? ? ?n + ?n+1 I ?1 (? ?n )?? log g(Yn+1 ; ? ?n ) . (9) ? Note that (9) does not correspond to a practical algorithm as I? (? ?n ) is usually unknown, although it can be estimated, for instance, by recursively averaging over the values of J(Yn ; ? ?n ). As discussed ? above however, this algorithm may be less robust than (8) because J(Yn ; ?n ) is (usually) not guar- anteed to be positive definite. In the following, we will refer to (8) as Titterington?s online algorithm and to (9) as the online gradient algorithm (in reference to the title of the paper by Lange, 1995). Note that both of these algorithms are based on the stochastic gradient approach and bear very little resemblance with the original EM algorithm. 2.3 The Proposed Online EM Algorithm We now consider an online approach which is more directly related to the principle underpinning the EM algorithm. The basic idea is to replace the expectation step by a stochastic approximation step, while keeping the maximisation step unchanged. More precisely, at iteration n, consider the function Q ? n (?) + ?n+1 E ? log f (Xn+1 ; ?) Yn+1 ? Q ? n+1 (?) = Q ? n (?) , (10) ?n and set ??n+1 as the maximum of the function ? ? Q ? n+1 (?) over the feasible set ?. One important advantage of (10) compared to (8) is that it automatically satisfies the parameter constraints without requiring any further modification. In addition, (10) does not explicitly require the inversion of a (d? × d? ) matrix. For further comparisons between both approaches, both practical and in terms of rate of convergence, we refer to the example of Section 2.4 and to the analysis of Section 3. Of course, this algorithm is of practical interest only if it is possible to compute and maximise ? Qn (?) efficiently. In the following we focus on the case where the complete data likelihood belongs to an exponential family satisfying the following assumptions. Let ·, · denotes the scalar product between two vectors of Rd and | · | the associated norm. Assumption 1. The complete data likelihood is of the form f (x; ?) = h(x) exp {??(?) + S(x), ?(?) } . (11) 4 The function def ¯(y; ?) = E? S(X) Y = y , s (12) is well defined for all (y, ?) ? Y × ?. There exists a convex open subset S, which is such that ? for all s ? S, (y, ?) ? Y × ? and ? ? [0, 1), (1 ? ?)s + ?¯ s(y; ?) ? S. def ? for any s ? S, the function ? ? ?(s; ?) = ??(?) + s, ?(?) has a unique global maximum over ? denoted ?¯ (s),i.e. ¯ (s) def ? = argmaxs?S ?(s; ?) . (13) Assumption 1 implies that the evaluation of E? log f (X; ?) Y , and hence the E-step of the EM algorithm, reduces to the computation of the expected value E? S(X) Y of the complete data sufficient statistic S(X). Indeed, the EM reestimation functional Q?? (Y1:n ; ?) is then defined by n Q?? (Y1:n ; ?) = ? n?1 ¯(Yi ; ? ? ); ? s . i=1 The (k + 1)-th iteration of the (batch mode) EM algorithm may thus be expressed as n ?k+1 = ? ? ¯ n?1 s ?k ) ¯(Yi ; ? , (14) i=1 where the M-step corresponds to the application of the function ?.¯ Note that the construction of the set S in Assumption 1 reflects the fact that in most applications of EM, the M-step is unambiguous only when a sufficient number of observations have been gathered. This point will be illustrated in the example to be considered in Section 4 below. Assumption 1 takes care of this issue in the case of the online EM algorithm. As an additional comment about Assumption 1, note that we do not require that ? be a one to one mapping and hence the complete data model may also correspond to a curved exponential family, where typically ? is of much lower dimension than ?(?) (see, for instance, Chung and B¨ ohme (2005); Capp´ e et al. (2006) for an example involving Gaussian densities with structured covariance matrices). In this setting, the proposed online EM algorithm takes the following form s ?n+1 = s ?n + ?n+1 (¯ ?n ) ? s s(Yn+1 ; ? ?n ) , ?n+1 = ? ? ¯ (? sn+1 ) . (15) Algorithm of that kind have a rather long history in the machine learning community. The idea of sequentially updating the vector of sufficient statistics has apparently been first proposed by Nowlan (1991), using a fixed step size (or learning rate) ?n = ? (see also Jordan and Jacobs, 1994). The online EM algorithm (15) is also closely related to the ?incremental? version of the EM algorithm derived by Neal and Hinton (1999). The incremental setting is more general than the recursive setting considered here, because the observations are not necessarily processed sequentially in time and might be used several times. The incremental EM algorithm of Neal and Hinton (1999) defines the (k+1)-th parameter estimate as ? ? min(k+1,n) ?k+1 = ? ? ¯ ?[min(k + 1, n)]?1 ?k+1,i ? , s (16) i=1 where s ?k+1,i = s ?k,i if i = Ik+1 and s ?k+1,Ik+1 = s ?k ). The index Ik+1 is typically chosen as ¯(YIk+1 ; ? k + 1 while k ? n ? 1 and runs through the data set, that is, Ik ? {1, . . . , n}, in a fixed or pseudo 5 random scanning order for subsequent iterations. When used in batch mode (that is when k > n) it is seen that it mostly differs from the traditional EM strategy in (14) by the fact that the parameters are updated after each computation of the conditional expectation of the complete data sufficient statistic corresponding to one observation. When used in online mode (k ? n), the algorithm of Neal and Hinton (1999) coincides with the proposed online EM algorithm with a step-size of ?k = 1/k (see Section 3 for further discussion of this particular choice of step sizes). A specific instance of the proposed online EM algorithm has been derived by Sato and Ishii (2000) for maximum likelihood estimation in the so-called normalised Gaussian network; this algorithm was later extended by Sato (2000) to a canonical exponential family (?(?) = ? in (11)) and a sketch of the proof of convergence, based on stochastic approximation results, was given. The online EM algorithm defined in (15) may be seen as a generalisation of this scheme. 2.4 An Example: Poisson Mixture Before analysing the convergence of the above algorithm, we first consider a simple example of ap- plication borrowed from Liu et al. (2006): consider the case of a mixture of m Poisson distributions m ?y j g(y; ?) = ?j e??j , for y = 0, 1, 2, . . . , (17) y! j=1 where the unknown parameters ? = (?1 , . . . , ?m , ?1 , . . . , ?m ) satisfies the constraints ?j > 0, m i=1 ?j = 1 and ?j > 0. In the mixture problem, the incompleteness is caused by the ignorance of the com- ponent of the mixture. Let W be a random variable taking value in {1, . . . , m} with probabilities {?1 , . . . , ?m }. The random variable W is called the regime or state and is not observable. The probability density defined in (17) corresponds to the assumption that Y is distributed, given that W = j, according to a Poisson random variable with parameter ?j . Note that in this case, as in all examples which involve the simpler missing data mechanism rather than the general latent data model introduced in Section 2.1, the complete data X simply consists of the couple (Y, W ) and hence conditional expectations of X given Y really boils down to expectations of W given Y . For the Poisson mixture, the complete data log-likelihood is given by m m log f (y, w; ?) = ? log(y!) + [log(?j ) ? ?j ] ?w,j + log(?j )y?w,j , (18) j=1 j=1 where ?i,l is the Kronecker delta symbol: ?i,l = 1 if i = l and ?i,l = 0 otherwise. The complete data likelihood may be rewritten as in (11) with h(y, z) = ? log(y!), S(y, w) = (S1 (y, w), . . . , Sm (y, w)) and ?(?) = (?1 (?), . . . , ?m (?)), where def ?w,j def log(?j ) ? ?j Sj (y, w) = , and ?j (?) = . y?w,j log(?j ) In this case, the conditional expectation of the complete data sufficient statistics is fully determined by the posterior probabilities of the mixture components defined by def ? j ?y je ??j w ¯j (y; ?) = P? [W = j|Y = y] = m y ??l , for i = 1, . . . , m . l=1 ?l ?l e The (n + 1)-th step of the online EM algorithm consists in computing, for j = 1, . . . , m, w ?n ) ¯j (Yn+1 ; ? s ?j,n+1 = s ?j,n + ?n+1 w ?n )Yn+1 ? s ¯j (Yn+1 ; ? ?j,n , ? ? j,n+1 = s ?j,n+1 (1) , ? j,n+1 = s ? ?j,n+1 (2) . (19) s ?j,n+1 (1) 6 Comparing with the generic update equations in (15), one recognises the stochastic approximation ¯ version of the E-step, in the first line of (19), followed by the application of ?. To compare with Titterington?s online algorithm in (8), one first need to evaluate the complete m Fisher information matrix I(?). To deal with the equality constraint j=1 ?j = 1, only the first m ? 1 weights are used as parameters and the remaining one is represented as ?m = 1 ? m?1 j=1 ?j as in Titterington (1984). The complete data Fisher information matrix defined in (5) is then given by ?1 1 T ?1 ?1 ?m m?1 1m?1 ? diag(?1 , . . . , ?m?1 ) 0(m?1)×m I(?) = , 0m×(m?1) diag(?1 /?1 , . . . , ?m /?m ) where the superscript T denotes transposition, 1 and 0 respectively denote a vector of ones and a matrix of zeros, whose dimensions are specified as subscript. Upon inverting I(?), the following expression for the (n + 1)-th step of Titterington?s online algorithm is obtained : ? ? j,n + ?n+1 w ? j,n+1 = ? ?n ) ? ? ¯j (Yn+1 , ? ? j,n , ?n ) ? ? j,n + ?n+1 w ? j,n+1 = ? ¯j (Yn+1 ; ? ? j,n Yn+1 ? ? . (20) ? ? j,n To make the connection more explicit with the update of the online EM algorithm, note that due to the fact that, in this simple case, there is an identification between some components of the vector ¯ (sj,n ) = ?j ), it is possible to rewrite (19) in of sufficient statistics and the weight parameters (i.e., ? terms of the latter only: ? j,n + ?n+1 w ? j,n+1 = ? ? ?n ) ? ? ¯j (Yn+1 ; ? ? j,n , ? j,n ? ? ? j,n + ?n+1 w ?n )Yn+1 ? ? ¯j (Yn+1 ; ? ? j,n ? ? j,n ? j,n+1 = ? . ? ? j,n+1 In this example, the two algorithms differ only in the way the intensities of the Poisson components are updated. Whereas the online EM algorithm in (19) does ensure that all parameter constraints are satisfied, it may happen, in contrast, that (20) produces negatives values for the intensities. 2.5 Extensions As previously mentioned, Neal and Hinton (1999) advocate the use of online algorithms also in the case of batch training with large sample sizes. The online algorithm then operates by repeatedly scan- ning through the available sample. In our setting, this use of the proposed online EM algorithm may be analysed by letting ? denote the empirical measure associated with the fixed sample X1 , . . . , Xn . The results to follow thus also apply in this context, at least when the data scanning order is random. In semi-parametric regression models, each observation Y comes with a vector of covariates Z whose distribution is usually unspecified and treated as a nuisance parameter. To handle latent data versions of regression models (mixture of regressions, mixture of experts, etc.?see Gr¨ un and Leisch, 2007; Jordan and Jacobs, 1994, as well as the example of Section 4) in our framework, one only needs to assume that the model consists of a parametric family {f (x|z; ?), ? ? ?)} of conditional pdfs. In this setting however, it is not possible anymore to compute expectations under the complete data distribution and the model can never be well-specified, as the distribution of Z is left unspecified. Thus Titterington?s algorithm in (8) does not directly apply in this setting. In contrast, the proposed algorithm straightforwardly extends to this case by considering covariate-dependent expectations of the sufficient statistics of f (x|z; ?), of the form s ¯(y, z; ?) = E? [S(X)|Y = y, Z = z], instead of (12). For notational simplicity, we do not explicitly consider the case where there are covariates but the following results also apply in this case, with suitable notation changes. 7 3 Convergence Issues 3.1 Consistency In this section, we establish the convergence of the proposed algorithm towards the set of stationary points of the Kullback-Leibler divergence between the actual observation density and the model likelihood. These results are the analogous of those given by Wang and Zhao (2006) for Titterington?s online algorithm, with a somewhat broader scope since we do not assume that the model is well- specified. The proofs corresponding to this section are given in Appendix A. As a preliminary, we state below a set of additional assumptions which are requested, in addition to those listed in Assumption 1, for the following results to hold. Assumption 2. The parameter space ? is a convex open subset of Rd? and ? and ? in (11) are twice continuously differentiable on ?. The function s ? ?¯ (s), defined in (13), is continuously differentiable on S, ? For some p > 2, and all compact subsets K ? S, ¯ (s)) p sup E? s ¯(Y ; ? <?. s?K To analyse the recursion (15), the first step consists in expressing it as a standard Robbins-Monro stochastic approximation procedure operating on the complete data sufficient statistics: s ?n+1 = s ?n + ?n+1 (h(? sn ) + ?n+1 ) , (21) where h : S ? Rd? is the so-called mean field given by def ¯ (s)) ? s , h(s) = E? s ¯(Y ; ? (22) and {?n }n?1 is a sequence of random variables representing stochastic perturbations defined by def ¯ (? ¯ (? ?n+1 = s ¯(Yn+1 ; ? sn )) ? E s ¯(Yn+1 ; ? sn )) Fn , (23) s0 , {Y1 }n where Fn is the ?-field generated by (? i=1 ). The aim of the Robbins-Monro procedure (21) is to solve the equation h(s) = 0. As a preliminary step, we first characterise the set of roots of the mean field h. The following proposition shows that, if s? belongs to def ? = {s ? S : h(s) = 0} , (24) then ?? = ¯ (s? ) ? is a stationary point of the Kullback-Leibler divergence between ? and g? , def ?(Y ) K (? g? ) = E? log . (25) g(Y ; ?) ¯ (s? ) is Proposition 3. Under Assumptions 1?2, if s? ? S is a root of h, i.e., h(s? ) = 0, then ? ? = ? a stationary point of the function ? ? K (? g? ), i.e., ?? K (? g? )|?=?? = 0. Conversely, if ? ? is a stationary point of ? ? K (? g? ), then s? = E? [¯ s(Y ; ? ? )] is a root of h. We then show that the function w : S ? [0, ?) defined by def w(s) = K ? g?(s) ¯ , (26) is a Lyapunov function for the mean field h and the set ?, i.e. for any s ? S, ?s w(s), h(s) ? 0 and ?s w(s), h(s) = 0 if and only if h(s) = 0. The existence of a Lyapunov function is a standard argument to prove the global asymptotic stability of the solutions of the Robbins-Monro procedure. This property can be seen as an analog of the monotonicity property of the EM algorithm: each unperturbed iteration s¯k+1 = s ¯k + ?k+1 h(¯ sk ) decreases the Kullback-Leibler divergence to the target distribution ?, provided that ?k+1 is small enough. 8 Proposition 4. Under Assumptions 1?2, ? w is continuously differentiable on S, ? for any compact subset K ? S \ ?, sup ?s w(s), h(s) < 0 . s?K Using this result, we may now prove the convergence of the sequence {? sk }. Denote by L = {? ? ? : ?? K (? g? ) = 0} the set of stationary points of the Kullback-Leibler divergence, and, for x ? Rm and A ? Rm , let d(x, A) = inf{y ? A, |x ? y|}. Theorem 5. Assume 1?2 and that, in addition, ? ? 2 1. 0 < ?i < 1, i=1 ?i = ? and i=1 ?i < ?, ?0 ? S and with probability one, lim sup |? 2. s sn , S c ) > 0. sn | < ? and lim inf d(? 3. The set w(?) is nowhere dense. Then, limk?? d(? ?k , L) = 0, with probability one. sk , ?) = 0 and limk?? d(? The first condition of Theorem 5 is standard for decreasing step-size stochastic approximation procedures (Kushner and Yin, 1997). It is satisfied for instance by setting ?i = ?0 i?? , with ? ? (1/2, 1]. The additional requirements that ?i be less than 1 and s ?0 be chosen is S is just meant to ensure that the whole sequence {?sk } stays in S (see Assumption 1). The rest of the second assumption of Theorem 5 correspond to a stability assumption which is not trivial. In general settings, the stability of the algorithm can be enforced by truncating the algorithm updates, either on a fixed set (see, e.g., Kushner and Yin, 2003, chapter 2) or on an expanding sequence of sets (see, e.g., Chen, 2002, chapter 2, or Andrieu et al., 2005). We do not explicitly carry out these constructions here to keep the exposition concise. 3.2 Rate of Convergence In this section, we show that when approaching convergence, the online EM algorithm is comparable to the online gradient algorithm in (9). The existence of such links is hardly surprising, in view of the discussions in Section 4 of Titterington (1984) and Section 3 of Lange (1995), and may be seen as a counterpart, for stochastic approximation, of the asymptotic equivalence of the gradient EM algorithm of Lange (1995) and the EM algorithm. To evidence these relations, we first express the online EM algorithm as a stochastic approximation procedure on ?. ?n }n?0 given by (15) Proposition 6. Under the assumptions of Theorem 5, the online EM sequence {? follows the recursion ?n+1 = ? ? ?n + ?n+1 I ?1 (? ?n )?? log g(Yn+1 ; ? ?n ) + ?n+1 ?n+1 (27) ? where limn?? ?n = 0 a.s. and I? (?) is defined in (7). Hence, the online EM algorithm is equivalent, when approaching convergence, to the online gradi- ent algorithm defined in (9) which coincides with Titterington?s online algorithm with I? (??n ) substi- ?n ). It is remarkable that the online EM algorithm can achieve a convergence performance tuted for I(? similar to that of the online gradient algorithm without explicit matrix approximation nor inversion. Although the recursion (27) will not lead to asymptotic efficiency, we can, under appropriate ?1/2 additional conditions, guarantee ?n -consistency and asymptotic normality. We use the weak con- vergence result presented in Pelletier, 1998, Theorem 1. 9 Theorem 7. Under the assumptions of Theorem 5, let ? ? be a (possibly local) minimum of the Kullback-Leibler divergence: ? ? K (? g? ). Denote by def H(? ? ) = I? ?1 ? (? ) ??2 ? K (? g? )|?=? ? , def ?(? ? ) = I? (? ) ? ?? log g ?? {?? log g ?? }T ?1 ? ?1 ? I? (? ) . Then, 1. H(? ? ) is a stable matrix whose eigenvalues have their real part upper bounded by ??(? ? ), where ?(? ? ) > 0. 2. Let ?n = ?0 n?? , where ?0 may be chosen freely in (0, 1) when ? ? (1/2, 1) but must satisfy ?0 > ?1/2 ? ?n = ? ? }, the sequence ?n ?(? ? ) when ? = 1; then, on the event ?(? ? ) = {limn?? ? ?n ? ? ? converges in distribution to a zero mean Gaussian distribution with covariance ?(? ? ), where ?(? ? ) is the solution of the Lyapunov equation (H(? ? ) + ?Id) ?(? ? ) + ?(? ? ) H T (? ? ) + ?Id = ??(? ? ) , (28) ?1 where ? = 0 if ? ? (1/2, 1) and ? = ?0 if ? = 1. Solving (28) is easy for a well-specified model, that is, when ? = g?? , as the FIM Iobs (? ? ) associated to the (observed) data model then satisfies Iobs (? ? ) = ?E?? ?2 ?1 ? ? log g? ? (Y ) = I? (? ) T = ??2 ? K (g? ? g? )|?=? ? = ? ?? log g ? ? {?? log g ? ? } . ?1 ? When ? = 0, the solution of the Lyapunov equation is given by ?(? ? ) = Iobs (? )/2, the asymptotic covariance matrix is equal to one-half the inverse of the FIM. When ? = 0, the Lyapunov equation cannot be solved in explicitly, except when the parameter is scalar (the result then coincides with (Titterington, 1984, Theorem 1). Note that using ?n = ?0 n?? with ? = 1 provides the optimal ? convergence rate of 1/ n but only at the price of a constraint on the scale ?0 , which is usually impossible to check in practice. On the other hand, using ? ? (1/2, 1) results in a slower convergence rate but without constraint on the scale ?0 of the step-size (except for the fact that it has to be smaller than 1). To circumvent this difficulty, we recommend to use the so-called Polyak-Ruppert averaging tech- nique (Polyak, 1990; Ruppert, 1988) as a post-processing step. Following Polyak (1990) ?see also Polyak and Juditsky (1992); Mokkadem and Pelletier (2006)?, if ? = ?0 n?? , with ? ? (1/2, 1), then the running average n ?n def ? = n?1 ?j , ? n?1 (29) j=n0 ? converges at rate 1/ n, for all values of ?0 . Furthermore, on the event ?(? ? ) defined in Theorem 7 ? ? above, n(?n ? ? ? ) is asymptotically normal, with asymptotic covariance matrix ?(? ? ) = H ?1 (? ? )?(? ? )H ?1 (? ? ) = ?1 ?1 ??2 ? K (? g? )|?=? ? ? ?? log g ?? {?? log g ?? }T ??2 ? K (? g? )|?=? ? , (30) which is known to be optimal (Kushner and Yin, 1997). If ? = g?? , the previous result shows that ?n is an asymptotically efficient sequence of estimates of ? ? , i.e. the asymptotic the averaged sequence ? ? ? covariance of n(?n ? ? ? ) is equal to the inverse of the (observed data) FIM Iobs (? ? ). 10 4 Application to Mixtures of Gaussian Regressions To illustrate the performance of the proposed method, we consider a regression model which, as discussed in Section 2.5, correspond to a case where the complete data FIM is not available. In contrast, we illustrate below the fact that the proposed algorithm, without explicitly requesting the determination of a weighting matrix does provide asymptotically efficient parameter estimates when Polyak-Ruppert averaging is used. The model we consider is a finite mixture of Gaussian linear regressions, where the complete data consists of the response variable R, here assumed to be scalar for simplicity, the dz -dimensional vector Z that contains the explanatory variables, and, W which corresponds, as in the example of Section 2.4, to a latent mixture indicator taking its value in the finite set {1, . . . , m}. We assume that given W = j and Z, R is distributed as a N (?j T Z, ? 2 ) Gaussian variable, while W is independent j of Z and such that P? (W = j) = ?j . Thus the parameters of the model are the mixture weights ?j and the regression vectors ?j and variances ?j 2 , for j = 1, . . . , m. As is usually the case in conditional regression models, we specify only the part of the complete data likelihood that depends on the parameters, without explicitly modelling the marginal distribution of the vector of regressors Z. In terms of our general notations, the complete data X is the triple (R, Z, W ), the observed data is the couple (R, Z) and the model is not well-specified, in the sense that the distribution of the observation (R, Z) is not fully determined by the model. We refer to Hurn et al. (2003) or Gr¨ un and Leisch (2007) and references therein for more information on mixture of regression models and their practical use. In the mixture of Gaussian regressions model, the part of the complete data log-likelihood that depends on the parameters may be written as 2 ? ? ?? m ? r ? ? Tz ? ? 1? 2 j ? log f (r, w, z; ?) = log(?j ) ? ?log ?j + ?w,j , (31) ? ? 2 ?j2 ? ? j=1 ? ? where ? denotes, as before, the Kronecker delta. To put (31) in the form given in (11), one need to define the statistics S = (S1,j , S2,j , S3,j , S4,j )1?j?m where S1,j (r, w, z) = ?w,j (scalar) , S2,j (r, w, z) = ?w,j rz (dz × 1) , T S3,j (r, w, z) = ?w,j rzz (dz × dz ) , S4,j (r, w, z) = ?w,j r 2 (scalar) . (32) As in the simple Poisson mixture example of Section 2.4, the E-step statistics only depend on the conditional expectation of the indicator variable W through s ¯1,j (r, z; ?) = w ¯j (r, z; ?) , s ¯2,j (r, z; ?) = w ¯j (r, z; ?)rz , s ¯j (r, z; ?)rzz T , ¯3,j (r, z; ?) = w s ¯j (r, z; ?)r 2 , ¯4,j (r, z; ?) = w (33) where T z)2 (r??j ?j ?j exp ? 1 2 2 ?j def w ¯j (r, z; ?) = E? [W = j|R = r, Z = z] = T z)2 . m ?l 1 (r??l l=1 ?l exp ? 2 2 ?l 11 Finally, it is easily checked that that the M-step is equivalent to an application of the function ¯: s ? ? ? ¯j (s), ? ¯ j (s), ? ¯j (s) 1?j?m where ? ¯ j (s) = s1,j , ¯j (s) = s?1 s2,j , ? 3,j ? ¯j2 ?T s2,j /s1,j . (s) = s4,j ? ? (34) j In this example, the role played by the set S in Assumption 1 is important: In order to apply (34), its required that the scalars s1,j belong to the open set (0, 1) and that the (dz +1)-dimensional matrices block-defined by s s2,j Mj = 3,j T , s2,j s4,j be positive definite, since ? 2 (s) is, up to normalisation by s , the Schur complement of M . These ¯j 1,j j constraints, for j = 1, . . . , m define the set S which is indeed open and convex. The function s ¯ defined in (33) however does never produce values of s which are in S. In particular, s ¯3,j (r, z; ?) is a rank one matrix which is not invertible (unless dz = 1). Hence the importance of using an initialisation s0 which is chosen in S. For the simulations below, we took care of this issue by inhibiting the parameter re-estimation step in (34) for the first twenty observations of each run. In other words, the first twenty observations are used only to build a up a value of s ?20 , using the first line of (15), which is in S with great probability. For illustration purpose, we consider a variation of a simple simulation example used in the flexmix R package (Leisch, 2004), where m = 2, ?1 = ?2 = 0.5, and 5U + V (when W = 1) R= , 15 + 10 U ? U 2 + V (when W = 2) where U ? Unif(0, 10) and V ? N (0, 92 ). In order to balance the asymptotic variances of the regres- sion parameters (see below) we used Z T = (1, U, U 2 /10) as the vector of regressors, hence the actual value of the regression parameter is ?1T = (0, 5, 0) for the first component and ? T = (15, 10, ?10). The 2 corresponding data is shown in Figure 1 where the points corresponding to both classes are plotted differently for illustration purpose, despite the fact that only unsupervised estimation is considered here. The labelling is indeed rather ambiguous in this case as the posterior probability of belonging to one of the two classes is between 0.25 and 0.75 for about 40% of the points. Clearly, the mixture of regressions model is such that the associated complete data sufficient likelihood has the form given in Assumption 1, where the marginal density of the explanatory variables in Z appears in the term h(x) since it does not depend on the parameters. Hence the previous theory applies straightforwardly and the online EM algorithm may be used to maximise the conditional likelihood function of the responses R given the regressors Z. However, the explicit evaluation of the complete data FIM I(?) defined in (5) is not an option here because the model does not specify the marginal distribution of Z. Titterington?s online algorithm may thus not be used directly. Applying the recursion in (8) without a weighting matrix is not recommended here as the regression parameters are greatly correlated due to the non-orthogonality of the regressors. In order to determine a suitable weighting matrix, one can use Fisher?s relation (4) which gives, for the regression parameters, ??j log g(r|z; ?) = E? ??j log f (R, W |Z; ?) R = r, Z = z; ? T Z)Z (R ? ?j T z)z (r ? ?j = E? ?W,j 2 R = r, Z = z; ? = w ¯j (r, z; ?) 2 . ?j ?j 12 80 70 60 50 40 Z 30 20 10 0 ?10 ?20 0 1 2 3 4 5 6 7 8 9 10 U Figure 1: 500 points drawn from the model: circles, points drawn from the first class, and, crosses, points drawn from the second class (the algorithms discussed here ignore the class labels). 40 ?2(1) 20 0 20 15 ?2(2) 10 5 0 0 ?5 ? (3) ?10 2 ?15 ?20 EM5 OL1 OL06 OL06a Figure 2: Box-and-whisker plots of the three components of ?2 (from top to bottom) estimated from 500 independent runs of length n = 100 for, EM5: five iterations of the batch EM algorithm, OL1: online EM algorithm with ?i = 1/i, OL06: online EM algorithm with ?i = 1/i0.6 , OL06a: online EM algorithm with ?i = 1/i0.6 and averaging started from the 50th iteration. The horizontal dashed line corresponds to the actual parameter value and the interval in bold at the right of each plot to the interquartile range deduced from the asymptotic normal approximation. 13 Hence, if we assume that the model is well-specified, the (observed) FIM Iobs (?) may be approximated, near convergence, by computing empirical averages of the form n T 1/n ??j log g(Ri |Zi ; ?) ??j log g(Ri |Zi ; ?) . i=1 As the online EM algorithm does not require such computations, this estimate has been used only to determine the FIM at the actual parameter value for comparison purpose. It is easily checked that due to the linearity of the model and the fact that both components have equal weights and variance, the covariance matrices for ?1 and ?2 are the same. The numerical approximation determined from a million simulated observations yields asymptotic standard deviations of (47.8, 22.1, 21.1) for the coordinates of ?j , with an associated correlation matrix of ? ? 1 ?0.87 0.75 ??0.87 1 ?0.97? . 0.75 ?0.97 1 As noted above, the coordinates of the regression vector are very correlated which would make the unweighted parameter-space stochastic approximation algorithm (i.e., (8) with an identity matrix ?n )) very inefficient. instead of I ?1 (? For run-lengths of n = 100 and n = 10, 000 observations, we illustrate the performance of the following four algorithmic options : EM5 Five iterations of the batch EM algorithm, using the whole data. OL1 The online EM algorithm with step size ?n = 1/n. OL06 The online EM algorithm with step size ?n = n?0.6 . OL06a The online EM algorithm with step size ?n = n?0.6 and averaging started from n0 = n/2 according to (29). Note that whereas OL1 and OL06a have the same computational complexity (as the averaging post- processing has a negligible cost), EM5 is significantly more costly requiring five times as much E-step computations; it is also non-recursive. All algorithm are started from the same point and run for 500 independent simulated replicas of the data. The results (for ?2 ) are summarised as box-and- whisker plots in Figure 2, for n = 100, and Figure 3 for n = 10, 000. Comparing both figures, one observes that OL06a is the only approach which appears to be consistent with a variance compatible with the asymptotic interquartile range shown on the right of each plot. EM5 (five iterations of batch EM) is clearly the method which has the less variability but Figure 3 suggests that it is not ? 1/ n consistent, which was indeed confirmed using longer runs not shown here. This observation confirms the fact that the use of the standard EM algorithm is prohibitive for large sample sizes as the fact that the batch size is larger cannot be compensated by the use of fewer EM iterations (see also Neal and Hinton, 1999). The online EM with step size ?n = 1/n (OL1) presents a bias which becomes very significant when n increases. According to Theorem 7, this problem could be avoided (asymptotically) by choosing a sufficiently small value of ?0 . For fixed n however, lowering ?0 can only reduce the perceived speed of convergence, which is already very slow, as illustrated by Figure 4. In contrast, the online EM algorithm with Polyak-Ruppert averaging (OL06a) appears to be very efficient: averaging significantly reduces the variability of the OL06 estimate, reducing it to a level which is consistent with the asymptotic interquartile range, while maintaining a systematic bias which vanishes as n increase, as expected. 14 17 16 ?2(1) 15 14 13 11 ?2(2) 10.5 10 9.5 9 ?9 ?9.5 ? (3) ?10 2 ?10.5 ?11 EM5 OL1 OL06 OL06a Figure 3: Same plots as in Figure 2 for signals of length n = 10, 000 (OL06a uses averaging started from the 50, 000th iteration). 30 ?2(1) 20 10 0 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 20 15 ?2(2) 10 5 0 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 ?5 ?2(3) ?10 ?15 ?20 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 samples Figure 4: Example of parameters trajectories for the three components of ?2 (from top to bottom) for a signal of length n = 3, 000: OL1, dashed line, OL06 doted line, OL06a solid line (with averaging started after 1, 000 iterations. 15 5 Conclusion Compared to other alternatives, the main advantages of the proposed approach to online parameter estimation in latent data models are its analogy with the standard batch EM algorithm, which makes the online algorithm easy to implement, and its provably optimal convergence behaviour. In addition, the combination of a slow parameter decrease (?n = n?0.5+? being a typical choice) with Polyak-Ruppert averaging appears to be very robust and, in many cases, the method can be used without tuning algorithmic parameters. ¯ (s) has to be explicit, which, for instance, would not A limitation is the fact that the function ? be the case for mixture of regression models with generalised link functions. Another extension of interest concerns non independent models and in particular hidden Markov models or Markov random fields. References ´ Moulines, and P. Priouret. Stability of stochastic approximation under verifiable C. Andrieu, E. conditions. SIAM J. Control Optim., 44(1):283?312 (electronic), 2005. O. Capp´ e, M. Charbit, and E. Moulines. Recursive EM algorithm with applications to DOA estima- tion. In Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., Toulouse, France, may 2006. H.-F. Chen. Stochastic approximation and its applications, volume 64 of Nonconvex Optimization and its Applications. Kluwer Academic Publishers, Dordrecht, 2002. P.-J. Chung and F. B. B¨ ohme. Recursive EM and SAGE-inspired algorithms with application to DOA estimation. IEEE Trans. Signal Process., 53(8):2664?2677, 2005. A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. J. Roy. Statist. Soc. Ser. B, 39(1):1?38 (with discussion), 1977. B. Gr¨ un and F. Leisch. Fitting finite mixtures of generalized linear regressions in R. Comput. Statist. Data Anal., 51(11):5247?5252, July 2007. M. Hurn, A. Justel, and C. P. Robert. Estimating mixtures of regressions. Comput. Statist. Data Anal., 12(1):55?79, 2003. ISSN 1061-8600. M. Jordan and R. Jacobs. Hierarchical mixtures of experts and the EM algorithm. Neural computa- tion, 6:181?214, 1994. H. J. Kushner and G. G. Yin. Stochastic Approximation Algorithms and Applications. Springer, 1997. H. J. Kushner and G. G. Yin. Stochastic Approximation and Recursive Algorithms and Applications, volume 35. Springer, New York, 2nd edition, 2003. K. Lange. A gradient algorithm locally equivalent to the EM algorithm. J. Roy. Statist. Soc. Ser. B, 57(2):425?437, 1995. F. Leisch. FlexMix: A general framework for finite mixture models and latent class regression in R. Journal of Statistical Software, 11(8):1?18, 2004. Z. Liu, J. Almhana, V. Choulakian, and R. McGorman. Online EM algorithm for mixture with application to internet traffic modeling. Comput. Statist. Data Anal., 50(4):1052?1071, 2006. G. McLachlan and T. Krishnan. The EM Algorithm and Extensions. Wiley, 1997. 16 A. Mokkadem and M. Pelletier. Convergence rate and averaging of nonlinear two-time-scale stochastic approximation algorithms. Ann. Appl. Probab., 16(3):1671?1702, 2006. ISSN 1050-5164. R. M. Neal and G. E. Hinton. A view of the EM algorithm that justifies incremental, sparse, and other variants. In M. I. Jordan, editor, Learning in graphical models, pages 355?368. MIT Press, Cambridge, MA, USA, 1999. S. Nowlan. Soft competitive adaptation: neural network learning algorithms based on fitting statistical mixtures. PhD thesis, School of Computer Science, Carnegie Mellon University, 1991. M. Pelletier. Weak convergence rates for stochastic approximation with application to multiple targets and simulated annealing. Ann. Appl. Probab., 8(1):10?44, 1998. ISSN 1050-5164. B. T. Polyak. A new method of stochastic approximation type. Autom. Remote Control, 51:98?107, 1990. B. T. Polyak and A. B. Juditsky. Acceleration of stochastic approximation by averaging. SIAM J. Control Optim., 30(4):838?855, 1992. D. Ruppert. Efficient estimation from a slowly convergent Robbins-Monro process. Technical Report 781, Cornell University, School of Operations Research and Industrial Engineering, 1988. M. Sato. Convergence of on-line EM algorithm. In prooceedings of the International Conference on Neural Information Processing, volume 1, pages 476?481, 2000. M. Sato and S. Ishii. On-line EM algorithm for the normalized Gaussian network. Neural Computa- tion, 12:407?432, 2000. M. A. Tanner. Tools for Statistical Inference. Springer, New York, 2nd edition, 1993. D. M. Titterington. Recursive parameter estimation using incomplete data. J. Roy. Statist. Soc. Ser. B, 46(2):257?267, 1984. D. M. Titterington, A. F. M. Smith, and U. E. Makov. Statistical Analysis of Finite Mixture Distri- butions. Wiley, Chichester, 1985. S. Wang and Y. Zhao. Almost sure convergence of Titterington?s recursive estimator for mixture models. Statist. Probab. Lett., 76:2001?2006, 2006. C. F. J. Wu. On the convergence properties of the EM algorithm. Ann. Statist., 11:95?103, 1983. A Proofs For m = (m1 , . . . , mr )T a differentiable function from ? to Rr , we define by ?? mT the d? × r matrix def whose columns are the gradients, ?? mT = [?? m1 , . . . , ?? mr ]. Thus the symbol ?? denotes either a vector or a matrix, depending on whether the function to which it is applied is scalar or vector-valued. With this convention, the matrix ?? mT is the transpose of the usual Jacobian matrix. Proof of Proposition 3. Let s? be a root of h, and set ? ? = ? ¯ (s? ). For any s ? S, the function ¯ (s). Therefore, ? ? ?(s; ?) has a unique stationary point at ? ??? ?(? ? ) + ?? ?T (? ? )s? = 0 . (35) 17 Note that, from Fisher?s identity, ?? log g(y; ?) = E? [?? log f (X; ?)| Y = y] = ??? ?(?) + ?? ?T (?)¯ s(y; ?) . (36) As ?? K (? g? ) = E? [?? log g(Y ; ?)], the latter relation implies ?? K (? g? ) = ??? ?(?) + ?? ?T (?)E? [¯ s(Y ; ?)] . (37) Since h(s? ) = 0, s? = E? [¯ s(Y ; ? ? )], and thus (35) and (37) imply ?? K (? g? )|?=?? = ??? ?(? ? ) + ?? ?T (? ? )s? = 0 , establishing the first assertion. Conversely, suppose that ?? K (? g? )|?=?? = 0 and set s? = E? [¯ s(Y ; ? ? )]. By (37), ??? ?(? ? ) + ?? ?T (? ? )s? = 0 . ¯ (s? ), Under Assumption 1, the function ? ? ??(?) + ?(?), s? has a unique stationary point at ? ? which is a maximum. Hence, ? = ?¯ (s ) establishing the second assertion. ? Proof of Proposition 4. Using (37) and the chain rule of differentiation, ¯T (s) ??? ?(? ?s w(s) = ??s ? ¯ (s)) + ?? ?T (? ¯ (s))E? s ¯ (s)) ¯(Y ; ? . (38) ¯ (s) is in ? and is the maximum of ? ? ?(s; ?), thus, For any s ? S, ? ?? ?(s; ?)|?=?(s) ¯ ¯ (s)) + ?? ?T (? = 0 = ??? ?(? ¯ (s))s . (39) Plugging this relation into (38) and using the definition (22), we obtain ¯T (s) ?? ?T (? ?s w(s) = ??s ? ¯ (s))h(s) . (40) def ¯ (s)) where ?(s, ?) = ?? ?(s; ?), we obtain By differentiating the function s ? ?(s; ? ¯ (s)) = ?s ?T (s; ? ?s ?T (s; ? ¯ (s)) + ?s ? ¯T (s) ?? ?T (s, ? ¯ (s)) . T Since ?s ?T (s; ?) = ?s [?? ?(s; ?)]T = ?? ?T (?) and ?? ?T (s, ?) = ?? 2 ?(s, ?), the latter relation may be alternatively written as T T ?s ?? ?(s; ?)|?=?(s) ¯ ¯ (s)) = ?? ?T (? ¯T (s) ?2 ?(s; ?) + ?s ? ¯ . ? ?=?(s) Because the function s ? ?? ?(s; ?)|?=?(s) is identically equal to 0, ?s [ ?? ?(s; ?)|?=?(s) T ¯ ¯ ] = 0, and thus, ¯ (s)) = ? ?2 ?(s; ?) ¯ ?? ?T (? ¯T (s) T . ?s ? (41) ? ?=?(s) Plugging this relation into (40) shows that ¯T (s) ?2 ?(s; ?) ¯T (s) T ?s w(s), h(s) = hT (s)?s ? ? ¯ ?=?(s) ?s ? h(s) . (42) Assumptions 1?2 imply that, for any s ? S, the matrix ?2 ? ?(s; ?) ¯ ?=?(s) is non-negative definite. Hence, for any s ? S, ?s w(s), h(s) ? 0 with equality if and only if ?s ? ¯ (s) , h(s) = 0. Assump- 2 tion 1 guarantees that the matrix ?? ?(?) is non-singular and thus (41) shows that for any s ? S, ?s ?¯ (s) is non-singular. Therefore, the condition ?s w(s), h(s) = 0 implies that h(s) = 0. The proof follows upon noting that s ? ?s w(s), h(s) is continuous. 18 Proof of Theorem 5. Under the stated assumption, for any ? > 0, there exists a compact K ? S and n, such that, P k?n {? sn ? K} ? 1 ? ?. Therefore, for any ? > 0, k P sup ?i ?i ? ? k?n i=n ? ? ? ? k ? P ?sup ?i ?i ½K (? si ) ? ?, {? si ? K}? + P ? {? si ? / K}? , k?n i=m i?n i?k k ? ? + P sup ?i ?i ½K (? si ) ? ? . k?n i=m i=n ?i ?i ½K (? k Note that Mn,k = si ) is a L2 -martingale, and that its angle-bracket is bounded by n ? 2 sup E s?K ? [|¯ s ((; Y ) , ¯ ? (s))|2 ] < ?. Using the Doob martingale inequality, we conclude that i=1 i 1 n P(sup |Mn,k | ? ?) ? M ?2 2 ?i sup E? [|¯ ¯ (s))|2 ] , s((; Y )1 , ? k?n i=1 s?K and the proof is concluded by applying Theorem 2.3 of Andrieu et al. (2005). Before proving Proposition 6, we need the following simple stability result. Lemma 8. Let p ? 1. Assume that for any compact subset K ? S, sups?K E? |¯ ¯ (s))|p < ? for s(Y ; ? ?n exists and belongs to S. Then, the sequence {¯ some p > 0 and that P? -a.s. lim s s(Yn+1 ; ?¯ (? sn ))}n?0 is bounded in probability, i.e. lim lim sup P |¯ ¯ (? s(Yn+1 ; ? sn ) |)| ? M = 0 . M ?? n?? Proof. Let K be a compact subset of S. We may decompose P |¯ ¯ (? s(Yn+1 ; ? sn ) |)| ? M as follows P |¯ ¯ (? s(Yn+1 ; ? sn ) |)| ? M ? P(? sn ? K) + P |¯ ¯ (? s(Yn+1 ; ? sn ) |)| ? M, s ?n ? K , ?p ¯ p sn ? K) + M sup E? |S(Y1 , ? (s))| P (? ? P(? sn ? K) , s?K which implies that lim sup P |¯ ¯ (? s(Yn+1 ; ? sn ) |)| ? M ? P( lim s ¯ (s))|p . ?n ? K) + M ?p sup E? |S(Y1 , ? n?? n?? s?K Proof of Proposition 6. A Taylor expansion with integral remainder shows that ?n+1 = ? ? ¯ s ?n + ?n+1 s ?n ) ? s ¯(Yn+1 ; ? ?n =? ¯T (? ?n + ?n+1 ?s ? sn ) s ?n ) ? s ¯(Yn+1 ; ? ?n + ?n+1 rn+1 , (43) where rn+1 is given by T def ?n ) ? s rn+1 = s ¯(Yn+1 ; ? ?n 1 × ¯T s ?s ? ?n + t?n+1 s ?n ) ? s ¯(Yn+1 ; ? ?n ¯T (? ? ?s ? sn ) dt . (44) 0 19 We will first show that limn?? rn = 0 a.s. Lemma 8 shows that {¯ ?n )}n?0 is bounded in s(Yn+1 ; ? probability, which we denote by s ? ¯(Yn+1 ; ?n ) = OP (1). Under the assumption of Theorem 5, s ?n = OP (1), which implies that s ?n ) ? s ¯(Yn+1 ; ? ?n = OP (1). Choose ? > 0 and then a compact subset K and a constant M large enough such that sn ? K) + lim sup P lim sup P (? s ?n ) ? s ¯(Yn+1 ; ? ?n ? M ? ? . (45) n?? n?? ¯ (·) is assumed continuous, it is uniformly continuous over every compact subset, Since the function ?s ? def i.e., there exists a constant ?0 such that K?0 = {s ? S, d(s, K) ? ?0 } ? S and sup ¯ (s + h) ? ?s ? |?s ? ¯ (s) | ? ? . (46) |h|??0 ,s?K Since limn?? ?n = 0 and (¯ ?n )?? s(Yn+1 ; ? sn ) is bounded in probability, ?n+1 (¯ ?n )?? s(Yn+1 ; ? sn ) converges to zero in probability, which we denote by ?n+1 (¯ ? s(Yn+1 ; ?n ) ? s ?n ) = oP (1). For ?0 > 0 satisfying (46), this implies in particular that lim P ?n+1 s ?n ) ? s ¯(Yn+1 ; ? ?n ? ?0 = 0 . n?? Therefore, lim sup P (|rn+1 | ? M ?) ? lim sup P ?n+1 s ?n ) ? s ¯(Yn+1 ; ? ?n ? ?0 n?? n?? sn ? K) + lim sup P + lim sup P (? s ?n ) ? s ¯(Yn+1 ; ? ?n ? M ? ? , n?? n?? showing that rn = oP (1). We now proceed with the main term (43). Eqs. (36) and (39) show that ¯ (s)) = ?? ?T ? ?? log g(y; ? ¯ (s) s ¯ (s)) ? s ¯(y; ? . (47) In addition, (41) shows that T ?1 ¯T (s) = ?? ?T (? ?s ? ¯ (s)) ? ?2 ? ?(s; ?) ¯ ?=?(s) . (48) Therefore, by combining (47) and (48), ¯T (? T ¯ (? ?1 ? ?n ) ?s ? sn ) s ¯(Yn+1 ; ? sn )) ? s ?n = I? (?n )?? log g(Yn+1 ; ? ?1 2 ?1 ? ?n ) . + ? ?? ?(? sn ; ?) ?n ?=? ? I? (?n ) ?? log g(Yn+1 ; ? By definition, the complete data FIM may be expressed as 2 I(?) = ?? ? (s; ?) s=E? [¯ s(Y ;?)] , We will show that: ?2 ? ?(? sn ; ?) ?n ?=? ? ?2 ? ?(s; ?) (s,?)=(E[¯ ?n )|Fn ],? s(Yn+1 ;? ?n ) = oP (1) . (49) Since the function (s, ?) ? ?2 ? ?(s; ?) is continuous, there exists ?1 > 0 such that, K?1 = {s ? S, d(s, K) ? ?1 } ? S and ?2 ¯ 2 ¯ sup ? ?(s + h; ? (s)) ? ?? ?(s; ? (s)) ? ? , |h|??1 ,s?K 20 where the set K is defined in (45). Under the stated assumption limn?? d(? sn , L) = 0 a.s., which implies that lim h(?sn ) = lim {E[¯ ¯ (? s(Yn+1 ; ? sn ))|Fn ] ? s ?n } = 0 , a.s. n?? n?? By combining the latter two equations, we therefore obtain lim sup P ?2 ? ?(? sn ; ?) ?n ?=? ? ?2 ? ?(s; ?) (s,?)=(E[¯ ?n )|Fn ],? s(Yn+1 ;? ?n ) ?? ? n?? sn ? K) + lim sup P (|h(? lim sup P (? sn )| ? ?0 ) ? ? , n?? n?? showing (49). Proof of Theorem 7. We use the definition of the recursive EM sequence given by Proposition 6. def ?1 (?)? K (? g ). The Jaco- The mean field associated to this sequence is given by h? (?) = ?I? ? ? bian of this vector field at ? ? is equal to H(? ? ). Since ? ? is a (possibly local) minimum of the Kullback-Leibler divergence, the matrix ?2 ? K (? g? ) ?=? ? is positive definite. Because the two ma- ?1 (? ? ) ?2 K (? g ) ?1/2 ?1/2 trices I? ? ? ?=? ? and I? (? ? ) ?2 ? K (? g? ) ?=? ? I? (? ? ) have the same eigenvalues, counting multiplicities, the eigenvalues of the matrix H(? ? ) = I? ?1 (? ? ) ?2 K (? g ) ? ? ?=? ? are all real and strictly positive. This shows the first assertion of the theorem; the proof of the second assertion follows directly from (Pelletier, 1998, Theorem 1). Acknowledgement The authors thank Maurice Charbit and Christophe Andrieu for precious feedback during the early stages of this work. 21