As one necessary step for a good modeling, this study was aimed at analyzing the sensitivity of GreenLab model for maize. When instantaneous value of biomass generation is considered as the output, the system tends to be linear, the level is above 94% in SRC(Standardized Regression coefficients)study. Conversion efficiency and characteristic surface are proved to be the most sensitive factors. In Sobol's measure, we excluded the two most sensitive factors in the analysis, then the system linearity tends to be weaker and we got the detailed sensitivity indexes for the other uncertain parameters, by which we get the clearer driven force of maize growth in different stages.

Sensitivity analysis of GreenLab model for maize

Qiongli WU * † ‡ , Paul-Henry Cournède * †

Introduction

Modeling represents a necessary tool for understanding plant growth and developing predictive tools for decision making. Good modeling practice requires that the modeler provides an evaluation of the confidence in the model. Sensitivity analysis (SA) is the study of how the variation (uncertainty) in the output of a mathematical model can be apportioned, qualitatively or quantitatively, to different sources of variation in the input factors of a model [START_REF] Saltelli | Global Sensitivity Analysis[END_REF].

So sensitivity Analysis (SA) has role of ordering by importance the strength and relevance of the inputs in determining the variations of the output variables of interest. Such information may provide some help for model assessment: SA can help the modeler to measure model adequacy (e.g.does the model fit observation) and relevance (e.g. is the model-based inference robust), to identify critical regions in the inputs space (e.g. which combination of factors corresponds to the highest risk), to detect interactions between factors, to establish priorities for research and experimentations and to simplify model structure [START_REF] Saltelli | Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models[END_REF].

If sensitivity analysis is quite usual in crop and plant growth models, it had long been restricted to local sensitivity analysis or to analysis of variance for linear models. An interesting exception is given by [START_REF] Ruget | Sensitivity analysis of a crop simulation model, stics, in order to choose the main parameters to be estimated[END_REF], for a variance-based analysis for the crop model STICS, with the objective of choosing the main parameters to be estimated. The analysis is first made module by module and then sensitivities of each module are compared for overall model outputs. The main factors addressed concern the interaction with the environment, which is of crucial interest.

Our objective in this paper is to study the interest of global sensitivity analysis and its last developments for the Green-Lab model (see [START_REF] Monod | Uncertainty and Sensitivity Analysis for Crop Models[END_REF], [START_REF] Cariboni | The role of sensitivity analysis in ecological modelling[END_REF], [START_REF] Yan | A dynamic, architectural plant model simulating resourcedependent growth[END_REF]), and more generally for a better understanding of source-sink dynamics and internal driving forces during plant growth. Potentially, this work should result in a better parameterization of the GreenLab model. As detailed by [START_REF] Guo | Parameter optimization and field validation of the functional-structural model greenlab for maize[END_REF], there are two groups of parameters in the model: the observed ones, that can be directly obtained by experimental observations, and the hidden ones, that must be estimated from experimental data by model inversion. For the observed parameters, we may need to set the level of accuracy of the experimental data and for those that mostly contribute to outputs' variabilities, more attention should be paid. Regarding hidden parameters, there is also a proper balance to find between the number of parameters used to describe the biophysical processes and the complexity of their estimation, which is always a bottleneck. By sensitivity analysis, we can rank the parameters by their significances to the system, or we can separate the parameters into different groups according to the interactions between the factors found by the sensitivity analysis. Then when we compute the estimation, according to the SA results, we can fix the least influential parameters, and we should pay more attention to those who play important roles in the outputs' variances.

In both cases, the sensitivity analysis may help to optimize the trade-off between experimental cost and accuracy. This is of crucial importance in the objective of developing a predictive capacity that scales from genotype to phenotype with the GreenLab model [START_REF] Letort | Quantitative genetics and functional-structural plant growth models: Simulation of quantitative trait loci detection for model parameters and application to potential yield optimization[END_REF].

This paper is organized as follows: In section 2, we present a general overview of the sensitivity analysis methods and what could be their insights regarding the specific issues raised by the GreenLab model. In section 3, we apply the methods to the GreenLab model for maize and analyze the result. In section 4, we summarize and discuss the methodology, as well as its perspective.

Sensitivity Analysis: an overview 2.1. Characteristics of SA method

A model is represented by a mapping f (a deterministic or stochastic function) which relates the inputs domain to the output space:

Y = f (X 1 , X 2 , ..., X k ) (1)
The input factors (X 1 , X 2 , ..., X k ) are supposed to be random variables described by identified probability distributions which reflect the uncertain knowledge of the system under analysis. Y is taken to be a scalar, i.e. even in the application we shall consider each output variable in turn [START_REF] Saltelli | Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models[END_REF]. In the analysis application of plant growth next in section 3, we consider the biomass for each kind of organ at one certain growth cycle as one output. So when the time goes, the output varies in turn, then we get one set of result for the analysis.

Brief overview of the available methods

The most common classifications of available SA methodologies distinguish between quantitative and qualitative methods and between local and global techniques.

Qualitative methods are aimed at screening, while quantitative techniques can be designed to give information on the amount of variance explained by each factor. In general, the choice of which kind of method to use is driven by cost, as local or qualitative methods are computationally less expensive [START_REF] Saltelli | Sensitivity analysis practices: Strategies for model-based inference[END_REF].

In local approaches (also known as one-at-a-time, OAT), the effect of the variation of a single factor is estimated by keeping all the others fixed at their nominal values. Yet they cannot include the effect of the shape of the density functions of the inputs, and they are not model independent.

Global approaches estimate the effect on the output of a factor when all the others are varying, enabling the identification of interactions in non-linear and/or nonadditive models. Generally, global approaches allow the use of model independent methods as they do not require assumptions of additivity or linearity. As a drawback, they are usually computationally expensive to estimate [START_REF] Cariboni | The role of sensitivity analysis in ecological modelling[END_REF].

The simplest and most intuitive way to obtain a local sensitivity index is to compute derivatives (see [START_REF] Grievank | Evaluating derivatives,Principles and techniques of algorithmic differentiation[END_REF], [START_REF] Tomovic | Modern analytical and computational methods in science and mathematics[END_REF], [START_REF] Varma | Parametric sensitivity in chemical systems[END_REF]). The sensitivity of the output Y to a perturbation of an input factor X i is estimated at a given value X * i , as

Y Xi = ∂Y ∂X i xi=x * i (2)
In situations where Y and X i have different range of uncertainty, a more balanced measure can be obtained normalizing the derivatives by the factors' standard deviations:

S σ Xi = σX i σY ∂Y ∂X i xi=x * i (3)
The estimation of these OAT methods can be easily implemented, but they are informative only if the model is linear/Quasi linear or if the range of uncertainty of the input factors is small [START_REF] Cariboni | The role of sensitivity analysis in ecological modelling[END_REF].

The Standardized Regression Coefficients (SRCs) can be viewed as the trade off between the local method and global method, concerning the advantages and shortcomings of the two: the accuracy of the analysis and the computing cost. It is based on the linear approximation of the model and Monte Carlo simulation. SRC method can demonstrate the shape of the probability distribution of every factor. The other important index produced by SRC is the model coefficient of determination, R 2 , which represents the fraction of the output variance explained by the regression model itself. In the condition of R 2 = 1, the system is linear and the SRCs can totally explain the variance of the output affected by each factor. Even when models are moderately non-linear (i.e.R 2 > 0.9), the SRCs index can still be valid from a qualitative way. It will only explain 90% of the model's variance. When R 2 gets small, then the SRCs could not be reliable sensitivity presentation.

The other important category of SA method is called 'Variance based' method. The basic concept for this kind of method is to decompose the output variance into the contributions imputable to each input factor. The most widely used ones are the FAST (Fourier Amplitude Sensitivity Test, see [START_REF] Cukier | Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients[END_REF], [START_REF] Cukier | Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients[END_REF], [START_REF] Koda | Automatic sensitivity analysis of kinetic mechanisms[END_REF]), and Sobol's method, see [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF]. FAST method decomposes the output variance V (Y ) by means of spectral analysis. Sobol's method bases on the same decomposition of variance, which is achieved by Monte Carlo methods in place of spectral analysis. What we used in this paper is Sobol's method.

Common Steps of all SA Methods

Here we give the readers the general steps of SA and specify them in our application to Greenlab model.

Specify the target function of interest.

It is easier to communicate the results of a sensitivity analysis when the target of interest has a direct relation to the problem tackled by the model. So to define the form of the output function that answers our questions would be the first thing before we do SA. In our analysis, the biomass generation and distribution of each kind of organ is the output we are interested in.

Assign a Probability

Density Function (pdf) to the selected factors. When this involves eliciting experts' opinion this is the most expensive and time consuming part of the analysis as it may concerning long-time experiment. Especially for the plant growth model, huge amount of practice is needed to get the pdfs of the parameters. Though at the beginning of modeling, the structure of the model is needed to be set for the experiment, the pdfs of the factors can also be derived from the experiential knowledge, which can to some sense reflect the natural regulation. Or we can consider the pdfs as one reasonable hypothesis that might be one of the uncertainties of the modeling, which we will take into account when we do the SA afterwards to see how much this hypothesis has affected the result. What we adopted here is uniform distribution for all the parameters. 2.3.3. Generate a matrix of inputs with that distribution(s) through an appropriate design. As in experimental design, a good design for numerical experiments should give a maximum of effects with a minimum of computed points. In our analysis in the environment of C++, we use the 'Mersenne Twister random number generator' to get sampling of the input data according to their pdfs.

2.3.4.

Evaluate the model and compute the distribution of the target function. This is the computer-time intensive step. The major index of the output (the mean value and the variance) is computed by Monte Carlo Method.

2.3.5. Select a method for assessing the influence or relative importance of each input factor on the target function. The factors affect the choosing of the method for SA concern many aspect: the purpose of the analysis, e.g. model simplification, factor prioritizing, uncertainty reduction, etc; the characteristic of the model: linear, nonlinear, quassi-linear, additive, etc; the computing cost the practitioner can afford...In [START_REF] Cariboni | The role of sensitivity analysis in ecological modelling[END_REF], the author has given out a good sketch of the various techniques available and their use as a function of computational cost of the model and dimensionality of the input space.

Case study for plant growth model 3.1. Model description

GreenLab is a functional-structural model that simulates plant development, growth and morphological plasticity. The model simulates individual organ production and expansion as a function of the growth cycle (GC), which corresponds to the phyllochron (thermal time in degree days between the appearances of two consecutive leaves on the main stem) before the end of plant organogenesis for maize [START_REF] Ma | Parameter optimization and field validation of the functional-structural model greenlab for maize at different population densities[END_REF].

Plant organogenesis is controlled by GC and plant morphogenesis depends on biomass production and allocation to expanding organs or competing sinks. Biomass production per plant is thus simplified according to the following mathematical equation:

q t (i) = E(i) • µ • Sp • [1 -exp(- λ Sp n(i) l=1 S l )] (4)
where E(i) is an environmental function at growth cycle i (generally related to the Photosynthetically Active Radiation), µ is a conversion efficiency, λ is analogous to the extinction coefficient of Beer-Lambert's Law, S p is a characteristic surface, S l is photosynthetic surface areas of all leaves. Here λ is set to 0.7.

Organs receive an incremental allocation of biomass that is proportional to their relative sink strengths. The relative sink strength for each type of organ is defined as a function of its age in terms of GCs: In the course of organ development, its relative sink strength is assumed to vary according to a beta function f o given by:

p o (j) = P o f o (j) (5
f o (j) = g o (j)/M o (1 ≤ j ≤ T o ) 0 (j ≥ T o ) g o (j) = (j + 0.5) αo (T o -j -0.5) βo (6)
M o = To j=1 g o (j)
The parameters α o and β o vary with organ type. This function is flexible to describe the shape of the sink variation and can be fitted to data by optimization. Only one parameter (B o ) was optimized to define the beta function for each organ type, and the values of two parameters α o and β o in eqn.( 6) are subsequently derived from B o using the constraints

α o + β o = 5 and B o = [α o /(α o + β o )].
At a given GC(i), d t (i) is the demand of all expanding organs :

d t (i) = o=b,s,e P o • i j=1 f o (j) + o=f,m P o • i-1 j=1 f o (j) (7)
Therefore, the instantaneous value of biomass got by each kind of organ at GC(i) is:

q t o (i) = min(To-1,i) k=max(0,i-Text) p o (k) q t (i) d t (i) (8)
where T ext stands for the GC at which organogenesis ceases for the whole plant. Note that for cob, the demand for biomass p f (k) is 0 before 15thGC, and tassel only exists for two GCs: the 21th and 22th. And the total time for the whole plant to be alive is 33GCs. What we are also interested in is the accumulated biomass value for each kind of organ, which is just to add the instantaneous value of the organs' corresponding GC, and we use the capital letter Q to denote this class of values, in order to distinguish them with the instantaneous ones:

Q t o (i) = i n=0 q o (n) (9)
Note that there is the similar exception for cob and tassel as mentioned above. We assume each uncertain input parameter has a uniform distribution which relies on the empirical knowledge, and we use the data from [START_REF] Guo | Parameter optimization and field validation of the functional-structural model greenlab for maize[END_REF], [START_REF] Ma | Parameter optimization and field validation of the functional-structural model greenlab for maize at different population densities[END_REF], [START_REF] Ma | Calibration of greenlab model for maize with sparse experimental data[END_REF] to set the mean value and variance of all the parameters, which are listed in Table .1

Local sensitivity analysis and its normalized version

The time evolution of q t o (i) for the mean values of uncertain input factors is given in Fig. 1. We took q t (i) for example to do the local sensitivity analysis and its normalized version presented in eqn.( 2) and eqn.(3) respectively, the result is showed in Fig. 2.

We used numerical simulation to get the derivatives of output q t (i). There are two ways to compute the standard deviations σY : the first one is Monte Carlo simulation, which we used and can be more reliable to calibrate the stochastic characteristics of the output variances but with heavy computing cost. To avoid this shortcoming, the other possibility is to use the following approximation:

(σY ) 2 ∼ = i (σX i ) 2 ∂Y ∂X i 2 (10)
Such an approximation is justified when the system is highly linear, which we will evaluate in the following section. For Fig. 2, the sensitivity is better calibrated by the normalized version, because we took into account the differences in the ranges of variations of all the uncertain input factors. As for Sp, we can not even pick it out from other factors in the pure local analysis, but with the normalized one, the great contribution from Sp is better demonstrated from about 12th GC to the end of the plant's growth. And the basic knowledge we got from this local analysis is that the factor µ contributes the most to the variance of biomass generation from the beginning to the end.

Note that the point at which we calculate the derivatives is important, it could be not reliable for the system which has the parameter-to-output curve with many apices, then the derivative will be either too big or too small, and the value at these special points can not stand for the whole character of sensitivity for this parameter. In this case, the result could be misleading, and we may miss the important details of the system. So we will introduce applications of other methods next.

Standardized Regression coefficients

To check whether the growth model for maize is linear in its input factors corresponding to the output we are interested in, the Standardized Regression coefficients (SRCs) can be used (see section 2.2). Indicating with the parameters listed in table.1 as X i , still with q t , q t o , Q t , Q t o the output variables as Y respectively, the following regressions are fit:

Y = b 0 + i b X i X i (11)
To change this dimensioned coefficients values to the standardized ones, the following comes the SRCs:

γ X i = σX i σY b X i (12)
Note that if the factors are independent and the true model (for example, the plant growth model for maize) is linear, then the following equation is hold (see [START_REF] Draper | Applied Regression Analysis[END_REF]):

i (γ X i ) 2 = 1 (13)
The sum of (γ X i ) 2 is therefore can be considered as a index of the model linearity. This sum is described also by the called 'model coefficient of determination', which is computed as:

R 2 Y = N d=1 (Y * -µ Y ) 2 (Y -µ Y ) 2 (14)
where µ Y is the mean value of Y in the Monte Carlo simulation (200000 runs for our simulation), Y * is the fitted value in eqn. [START_REF] Tomovic | Modern analytical and computational methods in science and mathematics[END_REF], Y is the actual value in each run, and N is the total number of runs.

Still we take the result with output q t as example, which is shown in Fig. 3. Note that the SRCs result is similar to the squared normalized local measures, which is not a coincidence, for the linear system, the two should be equal to each other. To affirm this linearity conclusion, the two methods of linearity computing are showed in Fig. 4.

Both of the linearity indexes show that even lower at the beginning and drop twice from around 5th to 15th GC, the whole linearity of the system is high (all above 0.94) when the output of q t is considered. So the result of SRCs and squared normalized local measures tend to be equal in our case. For there is not uniform standard to say the system could be considered linear when the linearity index is above one certain value for different systems (in [START_REF] Cariboni | The role of sensitivity analysis in ecological modelling[END_REF], the standard is set at 0.8), especially for the result in Fig. 4 with waves. So to be cautious, we consider the smooth part which maintain a relatively higher value as linear in our simplification, for maize, it is from 21th GC to the end. In this way, the regression analysis could be used as a preliminary step in SA to save the computing cost before we really do some complicate analysis. And if the linearity index is low, for instance below 90%, we must find one appropriate global analysis result for the analysis, which we will present next.

Variance decomposition-based sensitivity measure: Sobol's method

The basic idea of Sobol's method (see [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF]) is to decompose the function of eqn.(1) into terms of increasing dimensionality:

f (X 1 , ..., X k ) = f 0 + k i=1 f i (X i ) (15)
+ 1≤i<l≤k f il (X i , X l ) + • • • + f 1,2,...,k (X 1 , ..., X k )
Figure 3. The SRCs for output of q t . The vertical axes denotes for the SRCs depending on the curves explaining different parameters. If the input factors are mutually independent then there exists a unique decomposition of eqn. [START_REF] Koda | Automatic sensitivity analysis of kinetic mechanisms[END_REF], such that all the summands are mutually orthogonal. The variance of the output variable Y can be therefore decomposed into:

V (Y ) = k i=1 V i + 1≤i<l≤k V il + • • • + V 1,2,...,k (16)
Where V i , V il , V 1,2,...,k denote the variance of f i , f il , f 1,2,...,k respectively. In this approach the first-order sensitivity index for factor X i defined in eqn.( 1) is given by:

S i = V (E(Y |X i )) V (Y ) = V i (E -i (Y |X i )) V (Y ) (17)
where E and V indicate, respectively, the mean an variance operators and -i indicate all factors but i. The inner expectation is taken at a generic point in the space of variable X i , while the outer variance is over all possible values of this generic point. The higher order sensitivity indexes S i1,...,is are given by:

S i1,...,is = V i1,...,is V (Y )
for s > 1, eqn.( 16) can be rewritten in terms of sensitivity indexes as:

1 = k i=1 S i + 1≤i<l≤k S il + • • • + S 1,2,...,k (18)
In [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF], a computationally efficient design is discussed, which we used in our simulation processing.

As the result we got here with the same parameter space as in section 3.2 and section 3.3 with output q t is similar to the two for the high linearity of the system we mentioned in section 3.3. And we have already known that the parameters Sp and µ take as high as 75%∼98% of the summation for first-order index and the whole first-order sensitivity index contribute to 93%∼99% of the systems' variance in the former analysis. So as one alternative parameter space, we fixed Sp and µ at their mean value, and applied Sobol's method again, to our expectation, since the majority part of the uncertainty is excluded by us and the detailed sensitivity of the other uncertainty input should be clearer. And the result proved our deducing. The other important conclusion is that after we fixed the two most sensitive parameters, the system becomes less linear. Fig. 5 and Fig. 6 are showed to demonstrate these points. Note that without Sp and µ, α s and α b play the most important role in the system at the first 21 GCs, and after which α f takes the role instead, but with less amount than α s , until the 30th GC, the sensitivity of P f rises dramatically while all the other parameters dropping. This trend can be explained in the botany world: at the beginning of the maize's growth, before the fruit appears, the factors that control the competition of biomass acquisition must play an important role, here they are α s and α b for internodes' and blades' growth. Afterwards, the fruit (cob) starts appearing, then most of the biomass produced goes to the fruit, during which α f has a strong influence. But its value is generally Figure 6. The three kinds of linearity indexes with parameter space of Sp and µ fixed not as high as α s , because in this time period, the other organs still play relatively important roles. The biomass absorbed by the fruit is mainly produced by the blades and simultaneously, new blades and internodes still keep appearing while some of the old ones keep expanding and some die. So the whole variance of the system output tends to be shared more uniformly during this stage between most parameters. After the 31th GC, all other organs stop growing except the cobs, which still accumulate biomass, and the sensitivity coefficient of P f increases dramatically.

In Fig. 6, we show the two indexes which give the linearity level of the system, along with the summation of all S i . We noticed that the trends which they express agree to each other. So the summation of S i also has the characteristic of linearity index. Actually it can be explained from the definition of S i : the more the sum of the first-order indexes is close to 1 in eqn. [START_REF] Ma | Calibration of greenlab model for maize with sparse experimental data[END_REF], the more mutually independent the factors are, and the more linear the system is. In the new parameter space without Sp and µ, the system linearity is weaker, especially from 10th GC to 23th GC, which is the most non-linear stage of the growth, corresponding to abrupt changes in the allocation dynamics due to cob appearance.

Discussion

As in GreenLab model we need valid tools to locate the sources of uncertainty and evaluate them in a quantitive way. Sensitivity analysis is an efficient tool for this purpose and was used to assess the significance of GreenLab model parameters, based on the methods described in [START_REF] Cariboni | The role of sensitivity analysis in ecological modelling[END_REF], [START_REF] Saltelli | Sensitivity analysis practices: Strategies for model-based inference[END_REF].

For the simple ranking of the whole pack of parameters, we can get that the most important parameters for plant growth model of maize are µ and Sp, which can take responsible for almost 75%∼98% of the first-order sensitivity index. Coefficient of sink strength P o generally does not contribute so much to the system's variance, except for P f at the last 2 GCs of the growth, which means in the parameter estimation for modeling afterward, we can fix this group of parameters to simplify the computing.

If we fix the most important factors that we have already found, we can get the detailed sensitivity characters of the other factors, but notice that this kind of method relies on the fact that the second order sensitivity coefficients of the other parameters to the ones fixed should be relatively small, which means the co-relationship between the fixed factors and the kept variables should be slim, so that the disappearance of the variance for the most important ones will not imply missing the corresponding second order intercourse variances. In our case, this condition is confirmed by the linearity analysis of the system. Linearity analysis can also be used to choose the best SA method, with proper computing cost.

Though the distribution of the parameters is hard to get, if the experiment allows, we can get a much more specified analysis result when the distribution information can be as plentiful as possible. In fact, we have tried both Gaussian and uniform distribution for the parameters in our analysis, the results do not get so much different. Maybe it just needs an approximate range that is easier for us to get from the experiment than to get the accurate distribution of the parameters, see [START_REF] Helton | Survey of sampling based methods for uncertainty and sensitivity analysis[END_REF] and [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF] for the sampling method and proper method to get experimental data for this analysis.

In future works, we should account for the interactions with the environmental factors, then to assess the most important environment factors at different growth stage, according to the stress levels. As for this, the method we use should be changed for large scale sensitivity analysis, see [START_REF] Cacuci | Sensitivity and Uncertainty Analysis: Theory[END_REF], [START_REF] Cacuci | Sensitivity And Uncertainty Analysis: Applications to Large-Scale Systems[END_REF], [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF], in which the most popular method is Morris's [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF].

  )where o denotes organ type (b: leaf blade; s: sheath; e: internode; f : cob; m: tassel). P o is the sink strength associated to organ type o. For leaf blade, P b = 1 is set as a normalized reference. The relative sink strength for the first six short internodes is K e P e , with K e an empirical coefficient. f o (j) is an organ type-specific function of sink variation. A normalization constraintTo j=1 f o (j) = 1 (5A)is set, with T o being the maximum expansion duration for organ o that depends on the organ position.

Figure 1 .

 1 Figure 1. Biomass distributed to different organs at each GC

Figure 2

 2 Figure2. a. Pure local sensitivity of q t with respect to the factors listed in table.1, absolute value. b. Normalized local sensitivity of q t with respect to the factors listed in table.1, absolute value. Y in the vertical axes label indicates q t , X in the vertical axes label indicates the parameters, depending on the curves.

Figure 4 .

 4 Figure 4.(γ X i ) 2 and R 2 Y for output of q t

Figure 5 .

 5 Figure 5. Sobol's method for SA with parameter space of Sp and µ fixed

Table 1 .

 1 . Distribution characteristics for uncertain parameters

	Parameters	Mean value	Variance
	Ps	0.93	0.0537
	Pe	1.63	0.0941
	P f	162.18	9.3635
	Pm	1.33	0.0768
	B b	0.53	0.0306
	Bs	0.7	0.0404
	Be	0.88	0.0508
	B f	0.58	0.0335
	Sp	0.23	0.0133
	µ	0.0046	0.0003


Olivier Cappé

Kerrie Mengersen

Mike Titterington

Christian P Robert

Online Expectation-Maximisation *

Before entering into any more details about the methodological aspects, let's discuss the motivations behind the association of the two phrases "online (estimation)" and "Expectation-Maximisation (algorithm)" as well as their pertinence in the context of mixtures and more general models involving latent variables.

The adjective online refers to the idea of computing estimates of model parameters on-thefly, without storing the data and by continuously updating the estimates as more observations become available. In the machine learning literature, the phrase online learning has been mostly used recently to refer to a specific way of analysing the performance of algorithms that incorporate observations progressively [START_REF] Césa-Bianchi | Prediction, learning, and games[END_REF]. We dot not refer here to this approach and will only consider the more traditional setup in which the objective is to estimate fixed parameters of a statistical model and the performance is quantified by the proximity between the estimates and the parameter to be estimated. In signal processing and control, the sort of algorithms considered in the following is often referred to as adaptive or recursive [START_REF] Ljung | Theory and Practice of Recursive Identification[END_REF]Söderström, 1983, Benveniste et al., 1990). The word recursive is so ubiquitous in computer science that its use may be somewhat ambiguous and is not recommended. The term adaptive may refer to the type of algorithms considered in this chapter but is also often used in contexts where the focus is on the ability to track slow drifts or abrupt changes in the model parameters, which will not be our primary concerns.

Traditional applications of online algorithms involve situations in which the data cannot be stored, due to its volume and rate of sampling as in real-time signal processing or stream mining. The wide availability of very large datasets involving thousands or millions of examples is also at the origin of the current renewed interest in online algorithms. In this context, online algorithms are often more efficient -i.e., converging faster towards the target parameter value-and need fewer computer resource, in terms of memory or disk access, than their batch counterparts [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. In this chapter, we are interested in both contexts: when the online algorithm is used to process on-the-fly a potentially unlimited amount of data or when it is applied to a fixed but large dataset. We will refer to the latter context as the batch estimation mode.

Our main interest is maximum likelihood estimation and although we may consider adding a penalty term (i.e., Maximum A Posteriori estimation), we will not consider "fully Bayesian" methods which aim at sequentially simulating from the parameter posterior. The main motivation for this restriction is to stick to computationally simple iterations which is an essential requirement of successful online methods. In particular, when online algorithms are used for batch estimation, it is required that each parameter update can be carried out very efficiently for the method to be computationally competitive with traditional batch estimation algorithms. Fully Bayesian approaches -see, e.g., [START_REF] Chopin | A sequential particle filter method for static models[END_REF]-typically require Monte Carlo simulations even in simple models and raise some challenging stability issues when used on very long data records [START_REF] Kantas | An overview of sequential Monte Carlo methods for parameter estimation in general state-space models[END_REF].

This quest for simplicity of each of the parameter update is also the reason for focussing on the EM (Expectation-Maximisation) algorithm. Ever since its introduction by [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], the EM algorithm has been criticised for its often sub-optimal convergence behaviour and many variants have been proposed by, among others, [START_REF] Lange | A gradient algorithm locally equivalent to the EM algorithm[END_REF], [START_REF] Meng | The EM algorithm -an old folk song sung to a fast new tune[END_REF]. This being said, thirty years after the seminal paper by Dempster and his coauthors, the EM algorithm still is, by far, the most widely used inference tool for latent variable models due to its numerical stability and ease of implementation. Our main point here is not to argue that the EM algorithm is always preferable to other options. But the EM algorithm which does not rely on fine numerical tunings involving, for instance, line-searches, re-projections or pre-conditioning is a perfect candidate for developing online versions with very simple updates. We hope to convince the reader in the rest of this chapter that the online version of EM that is described here shares many of the attractive properties of the original proposal of [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] and provides an easy to implement and robust solution for online estimation in latent variable models.

Quite obviously, guaranteeing the strict likelihood ascent property of the original EM algorithm is hardly feasible in an online context. On the other hand, a remarkable property of the online EM algorithm is that it can reach asymptotic Fisher efficiency by converging towards the actual parameter value at a rate which is equivalent to that of the Maximum Likelihood Estimator (MLE). Hence, when the number of observations is sufficiently large, the online EM algorithm does becomes highly competitive and it not necessary to consider potentially faster-converging alternatives. When used for batch estimation, i.e., on a fixed large dataset, the situation is more contrasted but we will nonetheless show that the online algorithm does converge towards the maximum likelihood parameter estimate corresponding to the whole data. To achieve this result, one typically needs to scan the data set repeatedly. In terms of computing time, the advantages of using an online algorithm in this situation will typically depend on the size of the problem. For long data records however, this approach is certainly more recommendable than the use of the traditional batch EM algorithm and preferable to other alternatives considered in the literature.

The rest of this chapter is organised as follows. The first section is devoted to the modelling assumptions that are adopted throughout the chapter. In Section 2, we consider the large-sample behaviour of the traditional EM algorithm, insisting on the concept of the limiting EM recursion which is instrumental in the design of the online algorithm. The various aspects of the online EM algorithm are then examined in Sections 3 and 4.

Although the chapter is intended to be self-contained, we will nonetheless assume that the reader is familiar with the fundamental concepts of classical statistical inference and, in particular, with Fisher information, exponential families of distributions and maximum likelihood estimation, at the level of [START_REF] Lehmann | Theory of Point Estimation[END_REF]Casella (1998), Bickel andDoksum (2000) or equivalent texts.

Model and Assumptions

We assume that we are able to observe an independent sequence of identically distributed data (Y t ) t≥1 , with marginal distribution π. An important remark is that we do not necessarily assume that π corresponds to a distribution that is reachable by the statistical model that is used to fit the data. As discussed in Section 3.4 below, this distinction is important to analyse the use of the online algorithm for batch maximum-likelihood estimation.

The statistical models that we consider are of the missing-data type, with an unobservable random variable X t associated to each observation Y t . The latent variable X t may be continuous or vector-valued and we will not be restricting ourselves to finite mixture models. Following the terminology introduced by [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], we will refer to the pair (X t , Y t ) as the complete data. The likelihood function f θ⋆ is thus defined as the marginal

f θ (y t ) = p θ (x t , y t )dx t ,
where θ ∈ Θ is the parameter of interest to be estimated. If the actual marginal distribution of the data belongs to the family of the model distributions, i.e., if π = f θ⋆ for some parameter value θ ⋆ , the model is said to be well-specified; but, as mentioned above, we dot not restrict ourselves to this case. In the following, the statistical model (p θ ) θ∈Θ is assumed to verify the following key requirements.

Assumption 1. Modelling Assumptions

(i) The model belongs to a (curved) exponential family

p θ (x t , y t ) = h(x t , y t ) exp ( s(x t , y t ), ψ(θ) -A(θ)) , (1)
where s(x t , y t ) is a vector of complete-data sufficient statistics belonging to a convex set S, •, • denotes the dot product and A is the log-partition function.

(ii) The complete-data maximum-likelihood is explicit, in the sense that the function θ defined by

θ : S → Θ S → θ(S) = arg max θ∈Θ S, ψ(θ) -A(θ)
is available in closed-form.

Assumption 1 defines the context where the EM algorithm may be used directly (see, in particular, the discussion of [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. Note that (1) is not restricted to the specific case of exponential family distributions in canonical (or natural) parameterisation. The latter correspond to the situation where ψ is the identity function, which is particular in that p θ is then log-concave with the complete-data Fisher information matrix I p (θ) being given by the Hessian ∇ 2 θ A(θ) of the log-partition function. Of course, if ψ is an invertible function, one could use the reparameterisation η = ψ(θ) to recover the natural parameterisation but it is important to recognise that for many models of interest, ψ is a function that maps low-dimensional parameters to higher-dimensional statistics. To illustrate this situation, we will use the following simple running example.

Example 2 (Probabilistic PCA Model). Consider the probabilistic Principal Component Analysis (PCA) model of [START_REF] Tipping | Probabilistic principal component analysis[END_REF]. The model postulates that a d-dimensional observation vector Y t can be represented as

Y t = uX t + λ 1/2 N t , (2)
where N t is a centred unit-covariance d-dimensional multivariate Gaussian vector, while the latent variable X t also is such a vector but of much lower-dimension. Hence, u is a d × r matrix with r ≪ d.

Eq. ( 2) is thus fully equivalent to assuming that Y t is a centred d-dimensional Gaussian variable with a structured covariance matrix given by Σ(θ) = uu ′ + λI d , where the prime denotes transposition and I d is the d-dimensional identity matrix. Clearly, there are in this model many ways to estimate u and λ that do not rely on the probabilistic model of (2); the standard PCA being probably the most well-known. [START_REF] Tipping | Probabilistic principal component analysis[END_REF] discuss several reasons for using the probabilistic approach that include the use of priors on the parameters, the ability to deal with missing or censored coordinates of the observations but also the access to quantitative diagnostics based on the likelihood that are helpful, in particular, for determining the number of relevant factors.

To cast the model of (2) in the form given by (1), the complete-data model

  X t . . . Y t   ∼ N       0 . . . 0   ,     I r . . . u ′ . . . . . . . . . . . . . . . . u . . . uu ′ + λI d        
must be reparameterised by the precision matrix Σ -1 (θ), yielding

p θ (x t , y t ) = (2π) -d/2 exp trace Σ -1 (θ)s(x t , y t ) - 1 2 log |Σ(θ)| ,
where s(x t , y t ) is the rank one matrix

s(x t , y t ) = - 1 2   X t . . . Y t     X t . . . Y t   ′ .
Hence in the probabilistic PCA model, the ψ function in this case maps the pair (u, λ) to the (d + 1)-dimensional symmetric positive definite matrix Σ -1 (θ). Yet, Assumption 1-(ii) holds in this case and the EM algorithm can be used -see the general formulas given in [START_REF] Tipping | Probabilistic principal component analysis[END_REF] as well as the details of the particular case where r = 1 below.

In the following, we will more specifically look at the particular case where r = 1 and u is thus a d-dimensional vector. This very simple case also provides a nice and concise illustration of more complex situations, such as the Direction Of Arrival (DOA) model considered by [START_REF] Cappé | Recursive EM algorithm with applications to DOA estimation[END_REF]. The case of of a single factor PCA is also interesting as most required calculations can be done explicitly. In particular, it is possible to provide the following expression for Σ -1 (θ) using the block matrix inversion and Sherman-Morrison formulas:

Σ -1 (θ) =     1 + u 2 /λ . . . -u ′ /λ . . . . . . . . . . . . . . . . . . . . . . -u/λ . . . λ -1 I d     .
The above expression shows that one may redefine the sufficient statistics s(x t , y t ) as consisting solely of the scalars

s (0) (y t ) = y t 2 , s (2) (x t ) = x 2
t and of the d-dimensional vector s (1) (x t , y t ) = x t y t . The complete-data log-likelihood is then given by

log p θ (x t , y t ) = C st - 1 2 nd log λ + s (0) (y t ) -2u ′ s (1) (x t , y t ) + s (2) (x t ) u 2 λ , (3)
ignoring constant terms that do not depend on the parameters u and λ.

Note that developing an online algorithm for (2) in this particular case is equivalent to recursively estimating the largest eigenvalue of a covariance matrix from a series of multivariate Gaussian observations. We return to this example shortly below.

The EM Algorithm and the Limiting EM Recursion

In this section, we first review core elements regarding the EM algorithm that will be needed in the following. Next, we introduce the key concept of the limiting EM recursion which corresponds to the limiting deterministic iterative algorithm obtained when the number of observations grows to infinity. This limiting EM recursion is important both to understand the behaviour of classic batch EM when used with many observations and for motivating the form of the online EM algorithm.

The Batch EM Algorithm

In light of Assumption 1 and in particular of the assumed form of the likelihood in (1), the classic EM algorithm of [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] takes the following form.

Algorithm 3 (Batch EM Algorithm). Given n observations, Y 1 , . . . , Y n and an initial parameter guess θ 0 , do, for k ≥ 1,

E-Step S n,k = 1 n n t=1 E θ k-1 [ s(X t , Y t )| Y t ] , (4)

M-

Step

θ k = θ (S n,k ) . (5)
Returning to the single factor PCA model of Example 2, it is easy to check from the expression of the complete-data log-likelihood in (3) and the definition of the sufficient statistics that the E-step reduces to the computation of

s (0) (Y t ) = Y t 2 , E θ s (1) (X t , Y t ) Y t = Y t Y ′ t u λ + u 2 , E θ s (2) (X t ) Y t = λ λ + u 2 + (Y ′ t u) 2 (λ + u 2 ) 2 , (6)
with the corresponding empirical averages

S (0) n,k = 1 n n t=1 s (0) (Y t ) , S (1)
n,k = 1 n n t=1 E θ k-1 s (1) (X t , Y t ) Y t , S (2)
n,k = 1 n n t=1 E θ k-1 s (2) (X t ) Y t .
The M-step equations which define the function θ are given by

u k = θ(u) (S n,k ) = S (1) n,k /S (2) n,k , λ k = θ(λ) (S n,k ) = 1 d S (0) n,k -S (1) n,k 2 /S (2) n,k . (7)

The Limiting EM Recursion

Returning to the general case, a very important remark is that Algorithm 3 can be fully reparameterised in the domain of sufficient statistics, reducing to the recursion

S n,k = 1 n n t=1 Eθ (S n,k-1 ) [ s(X t , Y t )| Y t ]
with the convention that S n,0 is such that θ 0 = θ(S n,0 ). Clearly, if an uniform (wrt. S ∈ S) law of large numbers holds for the empirical averages of Eθ

(S) [ s(X t , Y t )| Y t ],
the EM update tends, as the number of observations n tends to infinity, to the deterministic mapping T S on S defined by

T S (S) = E π Eθ (S) [ s(X 1 , Y 1 )| Y 1 ] . (8)
Hence, the sequence of EM iterates (S n,k ) k≥1 converges to the sequence (T k S (S 0 )) k≥1 , which is deterministic except for the choice of S 0 . We refer to the limiting mapping T S defined in (8), as the limiting EM recursion. Of course, this mapping on S also induces a mapping T Θ on Θ by considering the values of θ associated to values of S by the function θ. This second mapping is defined as

T Θ (θ) = θ {E π (E θ [ s(X 1 , Y 1 )| Y 1 ])} . (9)
Using exactly the same arguments as those of [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] for the classic EM algorithm, it is straightforward to check that under suitable regularity assumptions, T Θ is such that

1. The Kullback-Leibler divergence D(π|f θ ) = log π(x) f θ (x) π(x)dx is a Lyapunov function for the mapping T Θ , that is, D π|f T Θ (θ) ≤ D (π|f θ ) .
2. The set of fixed points of T Θ , i.e., such that T Θ (θ) = θ, is given by

{θ : ∇ θ D(π|f θ ) = 0} ,
where ∇ θ denotes the gradient.

Obviously, ( 9) is not directly exploitable in a statistical context as it involves integrating under the unknown distribution of the observations. This limiting EM recursion can however be used in the context of adaptive Monte Carlo methods [START_REF] Cappé | Adaptive importance sampling in general mixture classes[END_REF] and is known in machine learning as part of the information bottleneck framework [START_REF] Slonim | Maximum likelihood and the information bottleneck[END_REF].

Limitations of Batch EM for Long Data Records

But the main interest of ( 9) is to provide a clear understanding of the behaviour of the EM algorithm when used with very long data records, justifying much of the intuition of [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. Note that a large part of the post 1980's literature on the EM algorithm focusses on accelerating convergence towards the MLE for a fixed data size. Here, we consider the related but very different issue of understanding the behaviour of the EM algorithm when the data size increases.

Figure 1 displays the results obtained with the batch EM algorithm for the single component PCA model estimated from data simulated under the model with u being a 20-dimensional vector of unit norm and λ = 5. Both u and λ are treated as unknown parameters but only the estimated squared norm of u is displayed on Figure 1, as a function of the number of EM iterations and for two different data sizes. Note that in this very simple model, the observation likelihood being given the N (0, uu ′ + λI d ) density, it is straightforward to check that the Fisher information for the parameter u 2 equals {2(λ + u 2 ) 2 } -1 , which has been used to represent the asymptotic confidence interval in grey. The width of this interval is meant to be directly comparable to that of the whiskers in the box-and-whisker plots. Boxplots are used to summarise the results obtained from one thousand independent runs of the method.

Comparing the top and bottom plots clearly shows that when the number of iterations increases, the trajectories of the EM algorithm converge to a fixed deterministic trajectory which is that of the limiting EM recursion (T k Θ (θ 0 )) k≥1 . Of course, this trajectory depends on the choice of the initial point θ 0 which was fixed throughout this experiment. It is also observed that the monotone increase in the likelihood guaranteed by the EM algorithm does not necessarily imply a monotone convergence of all parameters towards the MLE. Hence, if the number of EM iterations is kept fixed, the estimates returned by the batch EM algorithm with the larger number of observations (bottom plot) are not significantly improved despite the ten-fold increase in computation time (the E step computations must be done for all observations and hence the complexity of the E-step scales proportionally to the number n of observations). From a statistical perspective, the situation is even less satisfying as estimation errors that were not statistically significant for 2,000 observations can become significant when the number of observations increases. In the upper plot, interrupting the EM algorithm after 3 or 4 iterations does produce estimates that are statistically acceptable (comparable to the exact MLE) but in the lower plot, about 20 iterations are needed to achieve the desired precision. As suggested by [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF], this paradoxical situation can to some extent be avoided by updating the parameter (that is, applying the M-step) more often, without waiting for a complete scan of the data record. 3 Online Expectation-Maximisation

The Algorithm

The limiting-EM argument developed in the following section shows that when the number of observations tends to infinity, the EM algorithm is trying to locate the fixed points of the mapping T S defined in (8), that is, the roots of the equation

E π Eθ (S) [ s(X 1 , Y 1 )| Y 1 ] -S = 0 . ( 10
)
Although we cannot compute the required expectation wrt. the unknown distribution π of the data, each new observation provides us with an unbiased noisy observation of this quantity through 10) is thus an instance of the most basic case where the Stochastic Approximation (or Robbins-Monro) method can be used. The literature on Stochastic Approximation is huge but we recommend the textbooks by [START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF], [START_REF] Kushner | Stochastic Approximation and Recursive Algorithms and Applications[END_REF] for more details and examples as well as the review paper by [START_REF] Lai | Stochastic approximation[END_REF] for a historical perspective. The standard stochastic approximation approach to approximate the solution of ( 10) is to compute

Eθ (S) [ s(X n , Y n )| Y n ]. Solving (
S n = S n-1 + γ n Eθ (S n-1 ) [ s(X n , Y n )| Y n ] -S n-1 , (11)
for n ≥ 1, S 0 being arbitrary and (γ n ) n≥1 denoting a sequence of positive stepsizes that decrease to zero. This equation, rewritten in an equivalent form below, is the main ingredient of the online EM algorithm.

Algorithm 4 (Online EM Algorithm). Given S 0 , θ 0 and a sequence of stepsizes (γ n ) n≥1 , do, for n ≥ 1,

Stochastic E-Step S n = (1 -γ n )S n-1 + γ n E θ n-1 [ s(X n , Y n )| Y n ] , (12)

M-

Step

θ n = θ(S n ) . (13)
Rewriting ( 10) under the form displayed in ( 12) is very enlightening as it shows that the new statistic S n is obtained as a convex combination of the previous statistic S n-1 and of an update that depends on the new observation Y n . In particular it shows that the stepsizes γ n have a natural scale as their highest admissible value is equal to one. This means that one can safely take γ 1 = 1 and that only the rate at which the stepsize decreases needs to be selected carefully (see below). It is also observed that if γ 1 is set to one, the initial value of S 0 is never used and it suffices to select the initial parameter guess θ 0 ; this is the approach used in the following simulations.

The only adjustment to Algorithm 4 that is necessary in practice is to omit the M-step of (13) for the first observations. It typically takes a few observations the complete-data maximum likelihood solution is well-defined and the parameter update should be inhibited during this early phase of training [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF]. In the simulations presented below in Sections 3.2 and 3.4, the M-step was omitted for the first five observations only but in most complex scenarios a longer initial parameter freezing phase may be necessary.

When γ 1 < 1, it is tempting to interpret the value of S 0 as being associated to a prior on the parameters. Indeed, the choice of a conjugate prior for θ in the exponential family defined by (1) does result in a complete-data Maximum A Posteriori (MAP) estimate of θ given by θ(S n,k +S 0 /n) (instead of θ(S n,k ) for the MLE), where S 0 is the hyper-parameter of the prior [START_REF] Robert | The Bayesian Choice[END_REF]. However, it is easily checked that the influence of S 0 in S n decreases as n k=1 (1γ k ), which for the suitable stepsize decrease schemes (see beginning of Section 3.2 below), decreases faster than 1/n. Hence, the value of S 0 has a rather limited impact on the convergence of the online EM algorithm. To achieve MAP estimation (assuming a conjugate prior on θ) it is thus recommend to replace (13) by θ n = θ(S n + S 0 /n), where S 0 is the hyperparameter of the prior on θ.

A last remark of importance is that Algorithm 4 can most naturally be interpreted as a stochastic approximation recursion on the sufficient statistics rather than on the parameters. There does not exists a similar algorithm that operates directly on the parameters because unbiased approximations of the rhs. of (9) based on the observations Y t are not easily available. As we will see below, Algorithm 4 is asymptotically equivalent to a gradient recursion on θ n that involves an additional matrix weighting which is not necessary in (12).

Convergence Properties

Under the assumption that the stepsize sequence satisfies n γ n = ∞, n γ 2 n < ∞ and other regularity hypotheses that are omitted here (see [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF] for details) the following properties characterise the asymptotic behaviour of the online EM algorithm.

(i) The estimates θ n converge to the set of roots of the equation ∇ θ D(π|f θ ) = 0.

(ii) The algorithm is asymptotically equivalent to a gradient algorithm

θ n = θ n-1 + γ n J -1 (θ n-1 )∇ θ log f θ n-1 (Y n ) , (14)
where

J(θ) = -E π E θ ∇ 2 θ log p θ (X 1 , Y 1 ) Y 1 .
(iii) For a well specified model (i.e., if π = f θ⋆ ) and under Polyak-Ruppert averaging, θ n is Fisher efficient: sequences that do converge to θ ⋆ are such that √ n(θ nθ ⋆ ) converges in distribution to a centred multivariate Gaussian variable with covariance matrix I π (θ ⋆ ), where

I π (θ) = -E π [∇ 2 θ log f θ (Y 1
)] is the Fisher information matrix corresponding to the observed data.

Polyak-Ruppert averaging refers to a postprocessing step which simply consists in replacing the estimated parameter values θ n produced by the algorithm by their average

θn = 1 n -n 0 n t=n 0 +1 θ n ,
where n 0 is a positive index at which averaging is started [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF]Juditsky, 1992, Ruppert, 1988). Regarding the statements (i) and (ii) above, it is important to understand that the limiting estimating equation ∇ θ D(π|f θ ) = 0 may have multiple solutions, even in well-specified models. In practice, the most important factor that influences the convergence to one of the stationary points of the Kullback-Leibler divergence D(π|f θ ) rather than the other is the choice of the initial value θ 0 . An additional important remark about (i)-(iii) is the fact that the asymptotic equivalent gradient algorithm in ( 14) is not a practical algorithm as the matrix J(θ) depends on π and hence cannot be computed. Note also that J(θ) is (in general) neither equal to the complete-data information matrix I p (θ) nor to the actual Fisher information in the observed model I π (θ). The form of J(θ) as well as its role to approximate the convergence behaviour of the EM algorithm follows the idea of [START_REF] Lange | A gradient algorithm locally equivalent to the EM algorithm[END_REF].

From our experience, it is generally sufficient to consider stepsize sequences of the form γ n = 1/n α where the useful range of values for α is from 0.6 to 0.9. The most robust setting is obtained when taking α close to 0.5 and using Polyak-Ruppert averaging. The latter however requires to chose an index n 0 that is sufficiently large and, hence, some idea of the convergence time is necessary. To illustrate these observations, Figure 2 displays the results of online EM estimation for the single component PCA model, exactly in the same conditions as those considered previously for batch EM estimation in Section 2.3.

From a computational point of view, the main difference between the online EM algorithm and the batch EM algorithm of Section 2.1 is that the online algorithm performs the M-step update in (7) after each observation, according to (13), while the batch EM algorithm only applies the M-step update after a complete scan of all available observations. Both algorithms however require the computation of the E-step statistics following (6) for each observation. In batch EM, these local E-step computation are accumulated, following (4), while the online algorithm recursively averages these according to (12). Hence, as the computational complexity of the E and M steps are, in this case, comparable, the computational complexity of the online estimation is equivalent to that of one or two batch EM iterations. With this in mind, it is obvious that the results of Figure 2 compare very favourably to those of Figure 1 with an estimation performance that is compatible with the statistical uncertainty for observation lengths of 2,000 and larger (last two boxplots on the right in each display). ||u|| 2 α = 0.9

Figure 2: Online EM estimates of u 2 for various data sizes (200, 2,000 and 20,000 observations, from left to right) and algorithm settings (α = 0.9, α = 0.6 and α = 0.6 with Polyak-Ruppert averaging, from top to bottom). The box-and-whisker plots (outliers plotting suppressed) are computed from 1,000 independent replications of the simulated data. The grey regions corresponds to ±2 interquartile range (approx. 99.3% coverage) under the asymptotic Gaussian approximation of the MLE.

Regarding the choice of the stepsize decrease exponent α, it is observed that while the choice of α = 0.6 (middle plot) does result in more variability than the choice of α = 0.9 (top plot), especially for smaller observation sizes, the long-run performance is somewhat better with the former. In particular, when Polyak-Ruppert averaging is used (bottom plot), the performance for the longest data size (20,000 observations) is clearly compatible with the claim that online EM is Fisher-efficient in this case. A practical concern associated with averaging is the choice of the initial instant n 0 where averaging starts. In the case of Figure 2, we choose n 0 to be equal to half the length of each data record, hence averaging is used only on the second half of the data. While it produces refined estimates for the longer data sizes, one can observe that the performance is rather degraded for the smallest observation size (200 observations) due to the fact that the algorithm is still very far from having converged after just 100 observations. Hence, averaging is efficient but does require to chose a value of n 0 that is sufficiently large so as to avoid introducing a bias due to the lack of convergence.

In our experience, the fact that choices of α close to 0.5 are more reliable than values closer to the upper limit of 1 is a very constant observation. It may come to some surprise for readers familiar with the gradient descent algorithm as used in numerical optimisation, which shares some similarity with (12). In this case, it is known that the optimal stepsize choice is of the form γ n = a(n + b) -1 for a broad class of functions [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF]. However, the situation here is very different as we do not observe exact gradient or expectations but only noisy version of Figure 3: Four superimposed trajectories of the estimate of u 1 (first component of u) for various algorithm settings (α = 0.9, α = 0.6 and α = 0.6 with Polyak-Ruppert averaging, from top to bottom). The actual value of u 1 is equal to zero. them1 . Figure 3 shows that while the trajectory of the parameter estimates appears to be much rougher and variable with α = 0.6 than with α = 0.9, the bias caused by the initialisation is also forgotten much more rapidly. It is also observed that the use of averaging (bottom display) makes it possible to achieve the best of both worlds (rapid forgetting of the initial condition and smooth trajectories). [START_REF] Liang | Online EM for unsupervised models[END_REF] considered the performance of the online EM algorithm for large-scale natural language processing applications. This domain of application is characterised by the use of very large-dimensional models, most often related to the multinomial distribution, involving tens of thousands of different words and tens to hundreds of semantic tags. As a consequence, each observation, be it a sentence or a whole document, is poorly informative about the model parameters (typically a given text contains only a very limited portion of the whole available vocabulary). In this context, [START_REF] Liang | Online EM for unsupervised models[END_REF] found that the algorithm was highly competitive with other approaches but only when combined with mini-batch blocking: rather than applying Algorithm 4 at the observation scale, the algorithm is used on mini-batches consisting of m consecutive observations (Y m(k-1)+1 , Y mk+2 . . . , Y mk ) k≥1 . For the models and data considered by [START_REF] Liang | Online EM for unsupervised models[END_REF], values of m of up to a few thousands yielded optimal performance. More generally, mini-batch blocking can be useful in dealing with mixture-like models with rarely active components.

Application to Finite Mixtures

Although we considered so far only the simple case of Example 2 which allows for the computation of the Fisher information matrix and hence for quantitative assessment of the asymptotic performance, the online EM algorithm is easy to implement in models involving finite mixture of distributions.

Figure 4 displays the Bayesian network representation corresponding to a mixture model: for each observation Y t there is an unobservable mixture indicator or allocation variable X t that takes its value in the set {1, . . . , m}. A typical parameterisation for this model is to have a separate sets of parameters ω and β for, respectively, the parameters of the prior on X t and of the conditional distribution of Y t given X t . Usually, ω = (ω (1) , . . . , ω (m) ) is chosen to be the collection of component frequencies, ω (i) = P θ (X t = i), and hence ω is constrained to the probability simplex (ω (i) ≥ 0 and m i=1 ω (i) = 1). The observation pdf is most often parameterised as

f θ (y t ) = m i=1 ω (i) g β (i) (y t ) ,
where (g λ ) λ∈Λ is a parametric family of densities and β = (β (1) , . . . , β (m) ) are the component-specific parameters. We assume that g λ (y t ) has an exponential family representation similar to that of (1) with sufficient statistic s(y t ) and maximum likelihood function λ : S → Λ, which is such that λ(

1 n n t=1 s(Y t )) = arg max λ∈Λ 1 n n t=1 log g λ (Y t ).
It is then easily checked that the complete-data likelihood p θ belongs to an exponential family with sufficient statistics

s (ω,i) (X t ) = ½{X t = i} , s (β,i) (X t , Y t ) = ½{X t = i}s(Y t ) , for i = 1, . . . , m.
And the function θ(S) can then be decomposed as (ω,i) , for i = 1, . . . , m.. Hence the online EM algorithm takes the following specific form.

ω (i) = S (ω,i) , β (i) = λ S (λ,i) /S
Algorithm 5 (Online EM Algorithm for Finite Mixtures). Given S 0 , θ 0 and a sequence of stepsizes

(γ n ) n≥1 , do, for n ≥ 1, Stochastic E-Step Compute P θ n-1 (X t = i|Y t ) = ω (i) n-1 g β (i) n-1 (Y t ) m j=1 ω (j) n-1 g β (j) n-1 (Y t ) , and
S (ω,i) n = (1 -γ n )S (ω,i) n-1 + γ n P θ n-1 (X t = i|Y t ) , S (β,i) n = (1 -γ n )S (β,i) n-1 + γ n s(Y t )P θ n-1 (X t = i|Y t ) , (15)
for i = 1, . . . , m.

π being the empirical distribution πN (dy) = 1 N N t=1 δ Yt (dy) associated to the observed sample2 . Hence, if the online EM algorithm is applied by randomly drawing (with replacement) subsequent "pseudo-observations" from the finite set {Y 1 , . . . , Y N }, it converges to points such that

∇ θ D(π N |f θ ) = - 1 N ∇ θ N t=1 log f θ (Y t ) = 0 ,
that is, stationary points of the log-likelihood of the observations Y 1 , . . . , Y N . Property (ii) also provides an asymptotic equivalent of the online EM update, where the index n in ( 14) should be understood as the number of online EM steps rather than the number of actual observation, which is here fixed to N . In practice, it doesn't appear that drawing randomly the pseudo-observations make any real difference when the observations Y 1 , . . . , Y N are themselves already independent, except for very short data records. Hence, it is more convenient to scan the observation systematically in tours of length N in which each observation is visited in a predetermined order. At the end of each tour, k = n/N is equal to the number of batch tours completed since the start of the online EM algorithm.

To compare the numerical efficiency of this approach with that of the usual batch EM it is interesting to bring together two results. For online EM, based on ( 14), it is possible to show that √ n(θ nθ ∞ ) converges in distribution to a multivariate Gaussian distribution, where θ ∞ denotes the limit of θ n3 . In contrast, the batch EM algorithm achieves so-called linear convergence which means that for sufficiently large k's, there exists ρ < 1 such that θ kθ ∞ ≤ ρ k , where θ k denotes the parameter estimated after k batch EM iterations [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Lange | A gradient algorithm locally equivalent to the EM algorithm[END_REF]. In terms of computing effort, the number k of batch EM iterations is comparable to the number k = n/N of batch tours in the online EM algorithm. Hence the previous theoretical results suggest that

• If the number of available observations N is small, batch EM can be way faster than the online EM algorithm, especially if one wants to obtain a very accurate numerical approximation of the MLE. Note that from a statistical viewpoint, this may be unnecessary as the MLE itself is only a proxy for the actual parameter value, with an error that is of order 1/ √ N in regular statistical models.

• When N increases, the online EM algorithm becomes preferable and, indeed, arbitrary so if N is sufficient large. Recall in particular from Section 3.2 that when N increases, the online EM estimate obtained after a single batch tour is asymptotically equivalent to the MLE whereas the estimate obtained after a single batch EM iteration converges to the deterministic limit T Θ (θ 0 ).

In their pioneering work on this topic, [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] suggested an algorithm called incremental EM as an alternative to the batch EM algorithm. The incremental EM algorithm turns out to be exactly equivalent to Algorithm 4 used with γ n = 1/n up to the end of the batch tour only. After this initial batch scan, the incremental EM proceeds a bit differently by replacing one by one the previously computed values of

E θ t-1 [ s(X t , Y t )| Y t ] with E θ t-1 +(k-1)N [ s(X t , Y t )| Y t ]
when processing the observation at position t for the k-th time. This incremental EM algorithm is indeed more efficient than batch EM, although they do not necessarily have the same complexitya point that will be further discussed in the next section. For large values of N however, incremental EM becomes impractical (due to the use of a storage space that increases proportionally to N ) and less recommendable than the online EM algorithm as shown by the following example (see also the experiments of [START_REF] Liang | Online EM for unsupervised models[END_REF] for similar conclusions). Figure 5 displays the normalised log-likelihood values corresponding to the three estimation algorithms (batch EM, incremental EM and online EM) used to estimate the parameters of a mixture of two Poisson distributions (see Example 6 for implementation details)4 . All data is simulated from a mixture model with parameters ω (1) = 0.8, β (1) = 1 and β (2) = 3. In this setting, where the sample size is fixed, it is more difficult to come up with a faithful illustration of the merits of each approach as the convergence behaviour of the algorithms depend very much on the data record and on the initialisation. In Figures 5 and6, the two data sets were kept fixed throughout the simulations but the initialisation of both Poisson means was randomly chosen from the interval [0.5,5]. This randomisation avoids focussing on particular algorithm trajectories and gives a good idea of the general situation, although some variations can still be observed when varying the observation records.

Figure 5 corresponds to the case of an observation sequence of length N = 100. In this case it is observed that the performance of incremental EM dominates that of the other two algorithms, while the online EM algorithm only becomes preferable to batch EM after the second batch tour. For the data record of length N = 1,000 (Figure 6), online EM now dominates the other two algorithms, with incremental still being preferable to batch EM. Notice that it is also the case after the first batch tour, illustrating our claim that the choice of γ n = n -0.6 used here for the online EM algorithm is indeed preferable to the value of γ n = n -1 , which coincides with the update used by the incremental EM algorithm during the first batch tour. Finally, one can observe on Figure 6 that even after five EM iterations there are a few starting values for which batch EM or incremental EM gets stuck in regions of low normalised log-likelihood (around -1.6). These indeed correspond to local maxima of the log-likelihood and some trajectories of batch EM converge to these regions, depending on the value of initialisation. The online EM with the above choice of step-size appears to be less affected by this issue, although we have seen that in general only convergence to a stationary point of the log-likelihood is guaranteed.

Discussion

We conclude this chapter by a discussion of the online EM algorithm. First, the approach presented here is not the only option for online parameter estimation in latent variable models. One of the earliest and most widely used (see, e.g., [START_REF] Liu | Online EM algorithm for mixture with application to internet traffic modeling[END_REF] algorithm is that of [START_REF] Titterington | Recursive parameter estimation using incomplete data[END_REF] consisting of the following gradient update:

θ n = θ n-1 + γ n I -1 p (θ n-1 )∇ θ log f θ n-1 (Y n ) , (17)
where the matrix I p refers to the complete-data Fisher information matrix. For complete-data models in natural parameterisation -with ψ ≡ 1 in (1), I p (θ) coincides with ∇ 2 θ A(θ) and, as it does not depend on X t or Y t , is also equal to the matrix J(θ) that appears in (14). Thus, Titterington's algorithm is in this case asymptotically equivalent to online EM. In other cases, where I p and J differs, the recursion of ( 17) converges at the same rate as the online EM algorithm but is not Fisher-efficient. Another difference is the way Algorithm 4 deals with parameter constraints. Assumption 1 implies that the M-step update, taking into account the possible parameter constraints, is explicit. Hence, θ n = θ(S n ) does satisfy the parameter constraint by definition of the function θ. In the case of Example 6 for instance, the proposed update does guarantee that the mixture weight vector ω stays in the probability simplex. This is not the case with the update of (17) which requires reparameterisation or reprojection to handle possible parameter constraints.

As discussed in Section 3.4, the online EM algorithm is inspired by the work of [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] but distinct from their incremental EM approach. To the best of our knowledge, the online EM algorithm was first proposed by [START_REF] Sato | Convergence of on-line EM algorithm[END_REF] and [START_REF] Sato | On-line EM algorithm for the normalized Gaussian network[END_REF] who described the algorithm and provided some analysis of convergence in the case of exponential families in natural parameterisation and for the case of mixtures of Gaussians.

In Section 3.4, we have seen that the online EM algorithm is preferable to batch EM, in terms of computational effort, when the observation size is sufficiently large. To operate in batch mode, the online EM needs to be used repeatedly by scanning the data record several times so as to converge toward the maximum likelihood estimate. It is important to note however that one iteration of batch EM operating on N observations requires N individual E-step computations but only one application of the M-step update; whereas online EM applies the M-step update at each observation and, hence, requires N M-step updates per complete batch scan. Thus, the comparison between both approaches also depends on the respective numerical complexities of the E and M steps. The online approach is more attractive for models in which the M-step update is relatively simple. The availability of parallel computing resources would be most favourable to batch EM for which the E-step computations pertaining to each observation may be done independently. In contrast, in the online approach the computations are necessarily sequential as the parameter is updated when processing each observation.

As indicated in the beginning of this chapter, the main strength of the online EM algorithm is its simplicity and ease of implementation. As discussed above, this is particularly true in the presence of constraints on the parameter values. We have also seen in Section 3.2, that Algorithm 4 is very robust with respect to the choice of the stepsize γ n . In particular, the absolute scale of the stepsize is fixed due to the use of a convex combination in (12) and it is generally sufficient to consider stepsizes of the form γ n = n -α . We have shown that values of α of the order of 0.6 (i.e., closer to the lower limit of 0.5 than to the upper one of 1) yield more robust convergence. In addition, Polyak-Ruppert averaging can be used to smooth the parameter estimates and reach the optimal asymptotic rate of convergence that makes the online estimate equivalent to the actual MLE.

As illustrated by Figure 2, the online EM algorithm is not optimal for short data records of, say, less than 100 to 1,000 observations and in this case performing batch mode estimation by repeatedly scanning the data is recommended (see Figure 5 for the typical improvement to expect from this procedure). The main limitation of the online EM algorithm is to require that the M-step update θ be available in closed-form. In particular, it would be interesting to extend the approach to cases where θ needs to be determined numerically, thus making it possible to handle mixtures of Generalised Linear Models for instance (and not only mixtures of linear regressions as in [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF].

We finally mention two more directions in which recent works have proposed to extend the online EM framework. The first one concerns the case of non-independent observations and, in particular, of observations that follows a Hidden Markov Model (HMM) with Markov dependencies between successive states. [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] and Cappé (2009) have proposed an algorithm for HMMs that appears to be very reminiscent of Algorithm 4, although not di-rectly interpretable as a stochastic approximation recursion on the expected sufficient statistics. The other topic of importance is motivated by the many cases where the E-step computation of E θ n-1 [ s(X n , Y n )| Y n ] is not feasible. This typically occurs when the hidden variable X n is a continuous variable. For such cases, a promising solution consists in approximating the E-step using some form of Monte Carlo simulations (see, e.g., Cappé, 2009[START_REF] Del Moral | Forward smoothing using sequential Monte Carlo[END_REF] for methods that use sequential Monte Carlo approaches). However much conceptual and theoretical work remains to be done, as the consistency result summarised in Section 3.2 only extends straightforwardly to the -rather limited-case of independent Monte Carlo draws X (θ n-1 ) n from the conditionals p θ n-1 (x n |Y n ). In that case, s(X (θ n-1 ) n , Y n ) still provides an unbiased estimate of the limiting EM mapping in θ n-1 and the theory is very similar. In the more general case where Markov chain or sequential Monte Carlo simulations are used to produce the draws X (θ n-1 ) n , the convergence of the online estimation procedure needs to be investigated with care.

Figure 1 :

 1 Figure1: Convergence of batch EM estimates of u 2 as a function of the number of EM iterations for 2,000 (top) and 20,000 (bottom) observations. The box-and-whisker plots (outliers plotting suppressed) are computed from 1,000 independent replications of the simulated data. The grey region corresponds to ±2 interquartile range (approx. 99.3% coverage) under the asymptotic Gaussian approximation of the MLE.

Figure 5 :

 5 Figure5: Normalised log-likelihood of the estimates obtained with, from top to bottom, batch EM, incremental EM and online EM as a function of the number of batch tours (or iterations, for batch EM). The data is of length N = 100 and the box an whiskers plots summarise the results of 500 independent runs of the algorithms started from randomised starting points θ 0 .

Figure 6 :

 6 Figure6: Same display as in Figure5for a data record of length N = 1,000.

For a complete-data exponential family model in natural parameterisation, it is easily checked that E θ [ s(Xn, Yn)| Yn] = E θ [s(Xn, Yn)] + ∇ θ log f θ (Yn) and hence that the recursion in the space of sufficient statistics is indeed very close to a gradient ascent algorithm. However, we only have access to ∇ θ log f θ (Yn) which is a noisy version of the gradient of the actual limiting objective function D(π|f θ ) that is minimised.

Where δu(dy) denotes the Dirac measure localised in u.

Strictly speaking, this has been shown only for random batch scans.

Obviously, the fact that only log-likelihoods normalised by the length of the data are plotted hides some important aspects of the problem, in particular the lack of identifiability caused by the unknown labelling of the mixture components.

M-

Step

Example 6 (Online EM for Mixtures of Poisson). We consider a simplistic instance of Algorithm 5 corresponding to the mixture of Poisson distribution (see also Section 2.4 of [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF]. In the case of the Poisson distribution, g λ (y t ) = 1 yt! e -λ λ y t , the sufficient statistic reduces to s(y t ) = y t and the MLE function λ also is the identity λ(S) = S. Hence, the online EM recursion for this case simply consists in instantiating (15)-( 16) as

and

Use for Batch Maximum-Likelihood Estimation

An interesting issue that deserves some more comments is the use of the online EM algorithm for batch estimation from a fixed data record Y 1 , . . . , Y N . In this case, the objective is to save on computational effort compared to the use of the batch EM algorithm. The analysis of the convergence behaviour of online EM in this context is made easy by the following observation: Properties (i) and (ii) stated at the beginning of Section 3.2 do not rely on the assumption that the model is well-specified (i.e., that π = f θ⋆ ) and can thus be applied with