😊 Springer

1 Introduction
The expectation–maximization (EM) algorithm formulated by Dempster et al. ( 1977 ) is a general and popular algorithm for finding maximum likelihood estimates (MLEs) from incomplete data due to stability in convergence, simplicity in implementation and applicability in practice. However, the drawback of the EM algorithm is that its convergence is linear and very slow when the proportion of missing data is high. In order to circumvent the problem of slow convergence of the EM algorithm, various acceleration algorithms incorporating optimization methods with faster con- vergence rate have been proposed. The optimization methods include the multivariate Aitken method of Louis ( 1982 ) and Laird et al. ( 1987 ), the conjugate gradient method of Jamshidian and Jennrich ( 1993 ) and the quasi-Newton method of Lange ( 1995 ) and Jamshidian and Jennrich ( 1997 ). However, they require the matrix computation such as matrix inversion or evaluation of Hessian and Jacobian matrices and a line search for step length optimization. Therefore their acceleration algorithms tend to lack one or more of the nice properties of the EM algorithm, although they converge faster than the EM algorithm. Kuroda and Sakakihara ( 2006 ) developed the ε-accelerated EM algorithm for accel- erating the convergence of the sequence of EM iterations using the vector ε algorithm of Wynn ( 1962 ). The algorithm consists two steps: The first step is the expectation and maximization steps (the EM step) of the EM algorithm and the second step is the acceleration step using the vector ε algorithm. The vector ε algorithm is a fairly simple computational procedure and its implementation is very easy. Moreover, its computational cost is much less than that of any optimization method. The merit of the ε -accelerated EM algorithm is that it requires only the sequence of EM iterations for acceleration and maintains the nice properties of the EM algorithm. In the numer- ical experiments, Kuroda and Sakakihara ( 2006 ) demonstrated that the ε -accelerated EM algorithm significantly accelerates the convergence of the sequence of EM itera- tions. Wang et al. (2008) provided theorems concerning convergence and acceleration of the ε -accelerated EM algorithm. In order to further reduce the number of iterations and computation time, we improve the ε -accelerated EM algorithm using a re-starting procedure. The re-starting proce- dure embedding in the acceleration step finds an initial value for re-starting the EM step such that a newly generated sequence of EM iterations from the value moves quickly into a neighborhood of a stationary point. When applying the ε -accelerated EM algo- rithm to the newly generated sequence, its speed of convergence can be increased. Therefore the use of the re-starting procedure makes the ε -accelerated EM algorithm converge faster. We refer the ε -accelerated EM algorithm with a re-starting procedure to the ε R-accelerated EM algorithm. The paper is organized as follows. Section 2 describes the ε -accelerated EM algo- rithm. In Sect. 3 , we provide the ε R-accelerated EM algorithm and show some theo- retical results concerning its acceleration. Section 4 presents numerical experiments to illustrate the behavior of convergence of the ε R-accelerated EM algorithm. In Sect. 5 , we present our concluding remarks.
2 The ε-accelerated EM algorithm
Let y be the incompletely observed data in a sample space Ω y and x be the complete data augmented from y in a sample space Ω x . Assume that there exists some function h(x) = y relating x to y. Let f( |θ) denote a density function with an unknown d-dimensional parameter vector θ = (θ 1 , . . ,θ d ) T in a parameter space Θ and L (θ) = log ƒ ( · | θ ) be the log-likelihood function of θ . We denote the log-likelihood function for observed data y by L o (θ) = log f (y|θ), and the log-likelihood function for complete data x by L c (θ) = log f (x |θ). Denote the conditional expectation of L c (θ) given y and θ ′ by
The EM algorithm iteratively finds the sequence of EM estimates by
at the t th iteration for t = 0, 1, . . . . First we describe the EM algorithm for an initial value θ (0) ∈ Θ as follows: –  E-step: Calculate the expectation
– M-step: Find
The E- and M-steps are repeated to produce the sequence of estimates {θ ( ' ) t ≥ 0. We denoted by
The EM algorithm has the property
and {θ(t)}t>0 converges to a stationary point θ* . Next we briefly introduce the vector ε algorithm presented by Wynn ( 1962 ) to accelerate the convergence of a slowly convergent vector sequence. The algorithm is very effective for linearly convergent sequences. Let {θ ( t ) } t>0 be a linearly convergent vector sequence generated by an iterative computational procedure and converge to a stationary point θ* as t → ∞. Then the vector ε algorithm generates the accelerated sequence {θ(t)}t>0 for {θ(t)}t>0 by
where [θ]− 1 = θ/||θ||² and ||θ|| is the Euclidean norm of θ, see "Appendix 1" for details. Below we describe the ε -accelerated EM algorithm proposed by Kuroda and Sakak- ihara ( 2006 ). Given an initial value θ (0) ∈ Θ , the ε -accelerated EM algorithm repeats the following steps: –  EM step: Find
− ε-acceleration: Use the EM sequence (θ( t +1) , θ( t ) , θ( t −1) ) to generate the accel- erated sequence from until
where δ is a desired accuracy. The ε -accelerated EM algorithm generates two sequences: One is the EM sequence {θ ( t )} t>0 at the EM step, and the other is the accelerated sequence {θ ( t )} t>0 at the ε -acceleration step. The accelerated sequence does not make any effect on the EM sequence. Wang et al. ( 2008 ) showed that { θ ( t ) } r ≥ 0 converges to the same stationary point of {θ(t)}t>0 and it converges faster than {θ(t)}t>0. The ε-accelerated EM algo- rithm only uses the EM sequence, and thus it maintains the stability and simplicity of the EM algorithm.
3 Improvement of the ε-accelerated EM algorithm by using a re-starting
procedure: the εR-accelerated EM algorithm
The ε -accelerated EM algorithm generates two parallel sequences, the accelerated sequence { ∂( ′ ) t ≥ 0 and the EM sequence { θ( ′ ) t ≥ 0. But at the ε -acceleration step, ∂(t−1) may make the next EM estimate M ( δ(t−1) ) have a larger likelihood value than the current EM estimate θ ( ¹+¹ ) , that is, L o ( M ( θ ( − ¹ )) > L o ( θ ( ¹+¹ ) . Thus, when this occurs, we re-start the EM iterations with the initial value ∂( t − 1 ), stop the original EM sequence, and get ớ( t ) from ( ð(t − 1), M( ð(t − 1))). Notice that at the re-starting point, we still generate the EM sequence using three estimates obtained from the same initial value ∂( t − 1 ) . That is, we keep to always apply the ε -acceleration to a sequence obtained by the EM mapping M () from the same initial value. This re-starting algorithm proposed here is called the ε R-accelerated EM algorithm. By our experiments, the re-starting procedure is performed almost every time only by the re-starting condition L o ( M ( θ ( '− ' )) > L o ( θ ( t + ' ), and it inefficiently takes much computation time. Thus we add one more condition for re-starting ‖ ∂ ( t − 1 ) __ ∂(t-2) || 2 ≤ δ Re (> δ), and we reset δ Re = δ Re /10 k at each re-starting, where k is an integer, such as 1. By this condition, we control the re-starting frequency. For example, let δ = 10− 12 for stopping condition, and initialize δ R e = 1 and k = 1. Then the re- starting procedure is performed at most 12 times. The conditions for re-starting are summarized as follows: (i) L o (M(θ(t-l)) > L o (θ(t+l)), and (ii) ||θ(t-1) _ θ(t-2)||2 < δRe Condition (i) means that the likelihood can be increased by the re-starting. Condition (ii) is used to reduce the frequency of re-starting. This is the key idea of the re-starting procedure. The ε R-accelerated EM algorithm repeats the following steps: –  EM step: Find
− ε-acceleration step: Use (θ( t +1) , θ( t −1) ) to generate the accelerated sequence from – re-starting step: If Lo(M(d(t-l))) > Lo(d(t+l)) and ||d(t-l) - d(t-l)||2 < dRe, then set
update
and reset δ Re = δ Re /10 k . Set t = t + 1. Repeat the above steps until
The initial value δ R e and the size of decrement 10− k are related to the improvement of the computational efficiency of the ε -accelerated EM algorithm. When setting δ Re  = 1 and the decrement of 10 0 , the algorithm may re-start in every iteration after several iterations. Then the computation time for obtaining the re-starting EM sequence is twice that for doing {θ(t)} t>0 in each iteration. In the case of setting the larger size of decrement such as 10− 8 or 10− 10 for δ = 10− 12 , the re-starting is performed a few times. In both cases, the computation time of the ε R-accelerated EM algorithm may take longer than or the same as that of the ε -accelerated EM algorithm. Thus when the re-starting step effectively finds initial values, the ε R-accelerated EM algorithm greatly reduces the number of iterations and computation time for convergence. The advantage of the ε R-accelerated EM algorithm over the ε -accelerated EM algorithm is that it re-starts the iterations of the EM algorithm at a better current estimate and also keeps that the likelihood increases in the iterations. We give the pseudo-code of the ε R-accelerated EM algorithm in "Appendix 2". We describe the sequences obtained by two algorithms at iterations: The sequences of the ε -accelerated EM algorithm –EM
{θ(t)}t>0 : θ(0) θ(1) θ(2)     ... ーε-acceleration
The sequences of the ε R-accelerated EM algorithm —re-started EM
—re-started ε-acceleration
When the re-starting procedure is performed at the ( t + 1)th iteration, we obtain a new EM sequence { δ ( s )} s ≥ t +1 starting from θ ( t +1 ) = M ( θ ( t −1 ), and then the ε R-accelerated EM algorithm generates { δ ( s ) } s ≥ t +1 from the re-started EM sequence {∂(s) , ∂(s+l) , ∂(s+2)} . Note that ∂(t) is calculated using {∂(t) , ∂(t+l) , ∂(t+2)} . Now we discuss and compare the convergence of these sequences. We have from Meng and Rubin ( 1994 )
where
is the d × d Jacobian matrix for the mapping M(θ) = (M 1 (θ), . . . , M d (θ)) T . We assume that t is sufficiently large and θ ( t ) is in a neighborhood of θ * . Then the EM algorithm is essentially a linear iteration θ(t+1) = θ* = DM(θ*)(θ(t) = θ*) (4) with iteration matrix DM(θ*). The largest eigenvalue of DM(θ*) governs the rate of convergence of the EM algorithm and the linear iteration ( 4 ) has the same rate of convergence as the EM algorithm. Schafer ( 1997 , p. 59) provides the approximation
where λ is the largest eigenvalue of DM (θ*). By Eq. (5), we have
and therefore obtain
Under some conditions (Little and Rubin 1987 ), the estimate of θ has asymptotically the normal distribution with mean θ * . Then the closer to θ * the estimate is, the larger the density is. By L o (M(θ(t-1)) > L o (θ(t+1)), we have
The mapping of  θ ( ¹ ) is the same as the mapping of θ ( ¹ ) , and they are different only by the parameter values. Thus we have
Similar to Eq. ( 6 ), we obtain
Thus we have from Eqs. ( 6 ) and ( 7 )
When inequality ( 8 ) holds, we obtain the following result of the speed of convergence of {θ (t)}t>0 and {θ (t)}t>0. Theorem 1  Let { δ ( t )} t ≥ 0 be the sequence generated by the εR-accelerated EM algo- rithm and { θ ( ′ ) } r ≥ 0 be the sequence by the ε -accelerated EM algorithm. Then { θ ( ′ ) r ≥ 0 converges to θ* more quickly than {θ(t)}t>0. Proof  First we provide the result that {θ ( t ) t ≥ 0 and {θ (t ) t ≥ 0 converge to θ * faster than {θ(t)} t>0 and {θ(t)} , respectively. Wang et al. ( 2008 ) showed the following lemma.
🔥 Springer
Lemma 1
By Lemma 1 , we have || ∂( t − l ) − θ*|| = σ(|| θ( t + l ) − θ*||). Similarly applying the ε-acceleration to the re-starting sequence { δ ( ′ ) t ≥ 0, we also have || δ ( t − ¹ ) − θ * || = o(|| ğ ( t +1) − θ*||). That is, the ε -acceleration can speed up the convergence of the sequence {θ (t)} t ≥0. Next we show that the speed of convergence of { δ( t )} t ≥ 0 is faster than that of {6(')}r>0. We set Δθ(t-1) = θ(t) = θ(t-1) and η(t) = [[Δθ(t)]-1 = [Δθ(t-1)]-1]-1. Then we have from Eq. ( 3 )
Since
for sufficiently large t and θ ( t ) in a neighborhood of θ *, we have
and
We obtain from Wang et al. ( 2008 ) and the above equations
Then we have
and
😊 Springer Equation ( 9 ) becomes
In a similar way, we obtain We have from Eqs. ( 6 ) and ( 7 ) and inequality ( 8 ) By the above two equations, we get
and thus obtain
which completes the proof of the theorem. □
4 Numerical experiments
In this section, we provide numerical experiments using linear models and Poisson and normal mixture models. Then we investigate how much faster the εR-accelerated EM algorithm converges than the EM and ε -accelerated EM algorithms and compare the performance of the ε R-accelerated EM algorithm with that of the AEM algorithm pro- posed by Jamshidian and Jennrich ( 1993 ). The AEM algorithm is a conjugate gradient acceleration of the EM algorithm and is described in "Appendix 3". As the line search algorithm in the AEM algorithm, we use the gold-section method (R code is given in Jones et al. ( 2009 )). All computations are performed with the statistical package R (R Development Core Team 2013 ) executing on Intel Core 2 Duo 2.4 GHz with 4 GB of memory. The CPU times (in seconds) are measured by the function proc.time. 1 For all experiments, we set δ = 10− 12 for convergence of the algorithms, δ Re = 1 and k = 1 for the re-starting condition. 4.1 Linear models Consider the following linear model
where y is an n × 1 vector of n observations, X is an n × ( p + 1) known design matrix, β is a vector of ρ + 1 fixed parameters, and e is an n × 1 random vector of errors with e j ～ N (0, σ 2 ) for j = 1, . . . , n :
Thus we have
The analysis of variance and multiple linear regression belong to this model. We consider the case that y is missing at random. Then y is partitioned into the observed part y o and the missing part y m . Then the E-step for the ( t + 1)th estimate of y m is given by and the M-step for finding β( t +1) is given by where y(r+l) = (yo, y 4.1.1 Two-way analysis of variance We consider the two-way analysis of variance (two-way ANOVA) with two factors A and B having levels J and K , respectively. Then the two-way ANOVA model can be written as
where the terms λ A and λ β correspond to main effects of A and B . For the parameters, we impose the sum-to-zero constraints such that
In this numerical experiment, y consists of A and B with J = K = 4 levels and is generated from the normal distribution N (0, 100). We set the proportion of missing values of y to 50 %, and then find β = (λ 0 , {λ f } 1 ≤ j ≤ 4 , {λ β } 1 ≤ k ≤ 4 ) and σ 2 . The procedure is replicated 1,000 times. In order to evaluate the effect of the re-starting procedure, we examine the con- vergence behavior of the EM, ε -accelerated EM and ε R-accelerated EM algorithms measured by log 10 ‖ β( t ) − β( ∞ )‖ 2 (or log 10 ‖ β( t ) − β( ∞ )‖ 2 ), where β( ∞ ) is the MLE founded by the algorithm. Figure 1 illustrates the traces of these algorithms till the convergence of the ε R-accelerated EM algorithm attains. The figure shows that the sequence generated by the sR-accelerated EM algorithm converges after 37 iterations with δ = 10−12, while the sequence of the ε-accelerated EM algorithm matches four- digit precision to the MLE. For the experiment, the EM algorithm converges after 217
iterations and the ε -accelerated EM algorithm does after 86 iterations. We can see that the re-starting procedure works effectively to reduce the number of iterations. Table 1 presents the summary statistics of the numbers of iterations and CPU times of these algorithms. The ε R-accelerated EM algorithm requires about a quarter of the number of iterations and about one third of CPU time of the EM algorithm. The AEM algorithm greatly reduces the number of iterations but increases CPU time comparing with the other two acceleration algorithms.
4.1.2 Multiple linear regression
For the explanatory variables that are all continuous, X has the first column (1, . . . , 1) T corresponding to β 0 and the observed values of explanatory variables as other columns:
In this numerical experiment, we generate ( y, X ) for n = 100 and p = 30 from the multivariate normal distribution N(0, Σ ρ ), where Σ ρ is a ρ × ρ covariance matrix and is randomly chosen using a Wishart random number generator. The procedure is replicated 1,000 times for every algorithm. We show the results in Table 2 . The ε -accelerated EM algorithm converges about 1.7 times faster than the EM algorithm in both of the number of iterations and CPU time. The ε R-accelerated EM algorithm furthermore speeds up the convergence of the EM algorithm. The algorithm requires the number of iterations about 3.5 times smaller and CPU time about 2.8 times shorter than those of the EM algorithm. The AEM algorithm reduces the number of iterations as well as the εR-accelerated EM algorithm but not CPU time. The line search algorithm for finding a suitable α for each iteration may take a longer computation time.
CPU time We compare the performance of the ε R-accelerated EM algorithm with that of the ε -accelerated EM algorithm. Figure 2 presents the scatter plots of the ε R-accelerated EM algorithm by the ε-accelerated EM algorithm for the number of iterations and CPU time. We can see from the figure that the s R-accelerated EM algorithm converges in a smaller number of iterations than the ε -accelerated EM algorithm and well accelerates the convergence of {β ( t ) } t>0 when the ε -accelerated EM algorithm requires more iterations and longer CPU time. 4.2 Linear mixed model The linear mixed model is given by
for i  = 1, . . . , m , where y i is an n i × 1 observed vector, X i and Z i are known n i × p and n i × q design matrices corresponding to the p × 1 fixed effects vector β to be estimated and the q × 1 random effects vector u i . We assume that e j and u j are independent of each other, and e is distributed N (0, σ 2 R i ) and u i is N (0, D), where R i is a known n i × n i matrix and D is an unknown q × q positive definite covariance matrix. We treat u = ( u 1 , . . . , u m ) as missing data. The log-likelihood of θ = ( β, σ 2 , D ) given y = (y1, . . . , ym) is
where n =       m i=1 n i and V i = Z i DZ i + σ 2 R i . We denote the complete data vector m) T , where w i = (y T , u T ) T . The log-likelihood function of θ by w for w is
where
At the ( t + 1)th iteration, the E-step computes
and
and the M-step updates the estimate of θ by
The data shown in Table 3 are the average daily gain of two pigs of each litter in pounds (Snedecor and Cochran 1967 ). The experiment was designed so that each sire ( i  = 1, . . . , 5) is mated to a random group of dams ( j  = 1, 2) and each mating producing a litter in which two pigs are chosen ( k = 1, 2). The gain in weight of those two pigs is the criterion. The considered model is the linear mixed model with one random effect
where y i j is the observed average gain of two pigs by day in pounds produced by the i th sire and j th dam, β = (β 0 , β 1 , . . . , β 4 ) T is a sire effect, u ij is a random effect associ- ated with the i th sire and the j th dam and e ij is a random term. We assume that e ij and 😊 Springer
Table 3   Average daily gain of two pigs of each litter in pounds (Snedecor and Cochran 1967)
u ij are distributed N (0, σ ) and N (0, σ 2 ), respectively. Thus the parameters to be estimated are θ = (β,σ 2 ,σ 2 ). We start with σ 2 (0) and σ 2 (0) from U (0.001,1.5) and with β (0) fixed to the value of mean vector (2.57,0.098, − 0.040,0.063, −0.010). Figure 3 presents the traces of the log-likelihoods of the EM, ε -accelerated EM and ε R-accelerated EM algorithms till the converge of the ε R-accelerated EM algo- rithm attains. The figure indicates that the ε R-accelerated EM and ε -accelerated EM algorithms increase the log-likelihood L o (θ) much more than the EM algorithm. It is obvious the faster convergence of these acceleration algorithms over the EM algo- rithm. Table 4 is the results on 1,000 replications. The ε R-accelerated EM algorithm is about 5 times faster than the EM algorithm and also is about 3 times faster than the ε -accelerated EM algorithm in the number of iterations. We see that the ε R-accelerated EM algorithm can greatly improve the speed of convergence of the ε -accelerated EM algorithm.
4.3 Mixture models
Mixture models become increasingly popular due to the modeling flexibility and are one of the most interesting application areas of the EM algorithm. McLachlan and Peel (2000) provided a comprehensive book of the theory and applications of mixture models. We consider two-component Poisson and normal mixture models. A two- component mixture model for density of an observation y has the form
where λ is an unknown mixing proportion (0 < λ < 1), and f 1 (y |θ 1 ) and f 2 (y |θ 2 ) are the component density functions with parameters θ 1 and θ 2 , respectively. The log-likelihood function of θ = (λ,θ 1 ,θ 2 ) given observations y = (y 1 , . . . ,θ n ) is
In the setting of the EM algorithm, we regard y as incomplete data and introduce latent variables z = (z 1 , . . . , z n ), where z i  = (z i 1 , z i 2 ) is a binary vector defined as z ik = 1 if observation y i arises from the k th component of the mixture model, and z ik = 0 otherwise. Then the complete data x are given by x = (y, z) T . The density of x can be written as
where p ( z i 1 = 1 λ ) = λ and p ( z i 2 = 1 λ ) = 1 − λ . The log-likelihood function of θ given x is
Then, the E-step computes the conditional expectations w( y i , θ ′ ) of z i 1 given θ ′ and y i for i = 1, . . . , n and the M-step finds the MLEs of θ = (λ, θ 1 , θ 2 ). 4.3.1  Poisson mixture model We consider a mixture model of two Poisson distributions with θ = (λ, μ ι , μ 2 ) given by
The log-likelihood function of θ given observed data y = ( γ 1 , . . , Yn ) is
Let c j be the frequency of Y = y j for j = 0, 1, 2, . . . . Then Eq. ( 13 ) can be rewritten as
Then the EM estimates for the ( t + 1)th iteration are given by

🔥 Springer

where
In this numerical experiment, y is generated from the mixture of two Poisson distri- butions with θ = ( λ , μ 1 , μ 2 ) = (0.4, 3, 5) and is given in Table 5 . Then we generate 1,000 initial values of λ from the uniform distribution U (0.05, 0.95), and those of μ 1 and μ 2 from U (1, 15) under restriction μ (0) < μ 2 0) We compare the MLEs θ MLE from the EM and three acceleration algorithms with the true parameter values θ true = (0.4, 3, 5). The mean values of θ MLE from the EM algorithms are (0.268, 3.05, 4.78). Each acceleration algorithm can also find the same values as those from the EM algorithm. The standard errors of θMLE  from each algorithm are less than 10− 4 . Thus it can be seen that the MLEs of ( μ 1 , μ 2 ) are closed to the true values. The values of L o ( θ MLE ) from the algorithms are −1091.30 and are larger than L o (θ true) = −1092.74 that is the value of the log-likelihood of θ true We see from the mean values in Table 6 that the ε -accelerated EM algorithm reduces 1/1.6 of the number of iterations and 1/1.3 of CPU time of the EM algorithm. The AEM algorithm converges in a substantially smaller number of iterations than those of the other acceleration algorithms but takes longer CPU time. We note that the ε R- accelerated EM algorithm considerably converges faster than the EM algorithm in both of the number of iterations and CPU time. As shown in the mean values of the table, the algorithm requires about 1/10 of the number of iterations and about 1/7 of CPU time of the EM algorithm. Figure 4 is the scatter plots of the ε -accelerated EM and ε R-accelerated EM algorithms by the EM algorithm for the number of iterations and CPU time. The figure illustrates that the ε -accelerated EM algorithm increases linearly with the number of iterations and CPU time as the EM algorithm takes a larger number of iterations, while there is little variation in those for the εR-accelerated EM algorithm.

4.3.2 Univariate normal mixture model
Next we consider a mixture model of two univariate normal distributions. Let Φ k (y | μ k , σ ζ ) be the k th normal density function with a mean μ k and a variance σ ζ The density of an observation y is
where θ = (λ, μ 1 , μ 2 , σ 2 ). Then the E-step for the ( t + 1)th iteration compute the conditional expectations w(yi,θ) = λφι(yi|μι,σ 2 )/f(yi|θ) (i = 1, . . . ,π) and the M-step estimates θ(ï+ï) by We generate y of size n = 1,000 from the mixture of two univariate normal distribu- tions with θ = ( λ, μ 1 , μ 2 , σ 2 , σ 2 ) = (0.6, −1, 3, 8, 6). In setting 1,000 initial values of θ , λ ( ∅ ) is generated from the beta distribution Be (100λ , 100(1 − λ )), μ N ( μ k , 3) and σ 2(0) is from U (σ ξ /2, 2σ ξ ) for k = 1, 2. We also impose μ
🔥 Springer

We obtain from 1,000 simulation runs of the EM and its acceleration algorithms that the mean values of θ MLE are (0.599, −0.43, 2.59.8.96, 7.27), and the standard error for the MLE of μ is 0.001 and the errors of other parameters are 0.01. We see that the MLEs of ( λ , μ 2 ) are closed to the true values and, for other parameters, the errors between the MLEs and the true values are about ±1. We obtain L o (θ true) = −2599.45 and L o ( θ MLE ) = −2596.74 for all θ MLE  from the algorithms. Table 7 reports the results of the number of iterations and CPU time. The EM algo- rithm converges very slowly and takes more than 20000 iterations for convergence. The number of iterations and CPU time of the ε -accelerated EM algorithm reduces about one third of those of the EM algorithm. The AEM algorithm takes much fewer iterations than the ε -accelerated EM and ε R-accelerated EM algorithms and its CPU time is shorter than that of the ε -accelerated EM algorithm. The ε R-accelerated EM algorithm converges 30 times faster than the EM algorithm for both the number of iterations and CPU time. In terms of the performance of acceleration of the EM algo- rithm, the ε R-accelerated EM algorithm is clearly the best. Figure 5 shows the scatter plots of the ε -accelerated EM and ε R-accelerated EM algorithms by the EM algorithm for the number of iterations and CPU time. We obtain from the figure the same result as in the experiment of the Poisson mixture model. They show that the re-starting procedure can work effectively to reduce the number of iterations and CPU time for slowly convergent EM sequences.
4.3.3 Bivariate normal mixture model
Finally we consider a mixture model of two bivariate normal distributions. Let Φ k ( y | μ k , Σ k ) be the kth bivariate normal density function with a mean vector μ k and a covariance matrix Σ k . For the two-component mixture model, the density of an observation y is given by
where θ = (λ, μ 1 , μ 2 , Σ 1 , Σ 2 ). Then we update the EM estimates for the ( t + 1)th iteration by
where w(yi, θ) = λφι(yi |μι, Σι)/f(yi |θ).
😊 Springer
In this numerical experiment, we generate y of size n = 1,000 from the mixture of two bivariate normal distributions with
In choosing 1,000 initial values of θ , λ ( 0 ) is generated from Be (100 λ , 100 (1 − λ )), μ (0) is from the bivariate normal distribution N 2 (μ k , diag(3)) and Σ (0) is from the Wishart distribution W 2 ( Σ k ) for k = 1, 2. Table 8 shows the mean values of the MLEs from the EM algorithm and their standard errors. The values from the ε -accelerated and ε R-accelerated EM algorithms are also same in the table. We have the similar results for the errors between the MLEs and the true values but the larger standard errors in the experiments of the Poisson and univariate mixture models. The values of L o ( θ MLE ) from the algorithms are −5350.29, while L o (θ t rue) = −5352.81. Table 9 indicates that the EM algorithm requires many iterations and takes long computation time for each convergence, and thus it is very valuable to accelerate the convergence of the algorithm. The ε -accelerated EM algorithm takes a half of the number of iterations and CPU time of the EM algorithm. We see that the ε R- accelerated EM algorithm greatly reduce the number of iterations and CPU time and
The values in parenthesis are the standard errors of the MLEs of the EM algorithm
its speed of convergence is 12 times faster than that of the EM algorithm. The results illustrate that the ε R-accelerated EM algorithm can greatly improve the computational efficiency of the EM algorithm more than the ε -accelerated EM algorithm. Table 10 gives the total number of iterations and total CPU time for each algorithm. The EM algorithm requires 65 × 10 5 iterations and its computation time is 2.8 h, while the ε R-accelerated EM only takes 5.3 × 10 5 iterations and 14min and reduces 1/12 of those of the EM algorithm. The choice of initial values is of great importance for finding the highest likelihood in mixture models due to the local convergence of the EM algorithm. Several different initial values are employed to ensure that the global maximum is obtained, see Biernacki et al. ( 2003 ) and Karlis and Xekalaki (2003 ). The results show the possibility that the ε R-accelerated EM algorithm greatly contributes to shorten the total number of iterations and the CPU time required.
5 Concluding remarks
In this paper, we provided the ε R-accelerated EM algorithm for increasing the speed of convergence of the ε -accelerated EM algorithm by embedding the re-starting pro- cedure in the ε -acceleration step. The re-starting step is a very simple rule, and more- over is performed with increasing a little bit of computation time for each iteration. When the re-starting step finds θ ( t − 1 ) satisfying that L o ( M ( θ (t − 1 ))) > L o ( θ ( t +1 )) and ‖ ∂( t − 1 ) −  ∂( t − 2 ) ‖ 2 is smaller than a threshold, we re-start the EM iterations using M ( θ ( t − 1 ). We showed under some conditions that the ε R-accelerated EM algorithm accelerates the convergence of the EM sequence more than the ε -accelerated EM algorithm. Numerical experiments demonstrate that the ε R-accelerated EM algorithm gener- ates a faster convergent sequence than the ε -accelerated EM algorithm. We see that the speed of convergence of the ε R-accelerated EM algorithm is much faster than that of the ε -accelerated EM algorithm when the EM algorithm requires many iterations for convergence. Then the re-starting procedure can work effectively to reduce greatly the number of iterations and computation time of the ε -accelerated EM algorithm. We also compared the performance of the ε R-accelerated EM and AEM algorithms. Although the ε R-accelerated EM algorithm requires a larger number of iterations than the AEM algorithm, its CPU time is much shorter than that of the AEM algorithm. On the other hand, the AEM algorithm needs to evaluate the gradient of likelihood and find α using a line search algorithm, while the ε R-accelerated EM algorithm does not. Thus the ε R-accelerated EM algorithm has more advantageous than the AEM algorithm in terms of the computational efficiency. Mixture models are increasingly interest and popularity with the numerous devel- opments and the frequent applications, see Lee and Scott ( 2012 ), Pynea et al. ( 2009 ), Lee et al. ( 2011 ) and Lin ( 2009 ). The EM algorithm is largely used for the maxi- mum likelihood estimation of mixture models but its convergence tends to be slow. Improvement of convergence of the EM algorithm is an important topic. The results from the experiments show that the ε R-accelerated EM algorithm is useful due to its fast speed of convergence. Acknowledgments The authors would like to thank the editor and two referees for their valuable com- ments and helpful suggestions. This research is supported by the Japan Society for the Promotion of Science (JSPS), Grant-in-Aid for Scientific Research (C), No. 24500353.
Appendix 1: The vector ε algorithm
Let θ ( t ) denote a vector of dimensionality d that converges to a vector θ ( ∞ ) as t → ∞ . Let the inverse [θ]− 1 of a vector θ be defined by
where ||θ|| is the Euclidean norm of θ . In general, the vector ε algorithm for a sequence {θ (t)} t ≥ 0 starts with
and then generates a vector ε ( t,k+1) by
For practical implementation, we apply the vector ε algorithm for k = 1 to accelerate the convergence of {θ ( t )} t ≥ 0. From Eq. ( 14 ), we have
Then the vector ε ( t, 2 ) becomes as follows:
Appendix 2: Pseudo-code of the εR-accelerated EM algorithm Initialization We set the initial value of the EM step θ 0 , the desired precision δ , the threshold δ Re (> δ) and the size of decrement 10− k and determine the maximum number of iterations ( itrmax ). Iterations θ 1 ← M (θ 0) θοια ← θι itr ← 0 repeat itr ← itr + 1

🔥 Springer
if ‖θnew — θold‖² < δ or itr > itrmax then Termination of iterations end if θtmp ← M (θnew) if L o (θτmp) > L o (θ₂) then θշ ← θιmp θ」 ← θnew δRe ← δRe × 10− k end if end if θ old  ← θ new θο ← θι θ」 ← θշ end repeat
Appendix 3: The AEM algorithm of Jamshidian and Jennrich (1993)
We briefly introduce the AEM algorithm, which applies the generalized conjugate gradient (CG) algorithm to accelerate the EM algorithm. The idea of the AEM algo- rithm is that the change in θ ′ after an EM iteration ğ(θ ′ ) = M(θ ′ ) − θ ′ can be viewed approximately as a generalized gradient. Thus the AEM algorithm treats the EM step as a generalized gradient and uses the generalized CG algorithm as an EM accelerator. The generalized gradient g(θ ′ ) is given by
In the AEM algorithm, first the EM algorithm runs until the difference between 2L o (θ ( t ) and 2L o (θ (t − 1) ) falls below one. Then the CG accelerator updates the EM estimate θ ( t ) for obtaining θ ( t +1). 1. Set t = 0 and d(t) = g(θ(t)). 2. Find α ( t ) to maximize L o (θ ( t ) + α d ( t ) using a line search algorithm. 3. Update θ(r+1) = θ(r) + α(r)d(r). 4. Compute
where g(θ) is the gradient of Lo(θ) . 5. Update d ( i + 1 ) = g(θ ( i + 1 ) = y( i )d( i ) and then t = t + 1. 6. Repeat Steps 2 to 5 until ||g(θ ( ′ )||² ≤ δ . 🔥 Springer