Abstract

1. Introduction
The EM algorithm of Dempster et al. (1977) is a very general and popular iterative computational algorithm that is used to find maximum likelihood estimates (MLEs) from incomplete data and is widely used to perform statistical analysis with missing data, bécause of its stability, flexibility and simplicity. However, a common criticism is that the convergence of the EM algorithm is slow, since its speed of convergence is linear. In particular, the convergence of the EM algorithm is quite slow when the proportion of missing data is high. In order to speed up the convergence of the EM algorithm, various acceleration algorithms have been proposed. Louis (1982) suggested the EM algorithm incorporating the multivariate version of Aitken's acceleration method. Jamshidian and Jennrich (1993) proposed an acceleration algorithm based on conjugate gradients. Lange (1995) used a quasi-Newton algorithm to accelerate the EM algorithm. Using the Newton-type accelerators requires the computation of a matrix inversion at each iteration. Then its computation is likely to become rapidly complicated as the number of parameters increases, and numerical instabilities are also expected. Therefore, these accelerations do not have some of the attractive features of the EM algorithm, such as stability, flexibility and simplicity. Some extensions of the EM algorithm have been proposed. Meng and Rubin (1993) presented the expectation conditional maximization (ECM) algorithm that improves the maximization process in the M-step. Liu and Rubin (1994) suggested the expectation conditional maximization either (ECME) algorithm that is an extension of the ECM algorithm. McLachlan and Krishnan (1997) provided a comprehensive account of the EM algorithm. Geng et al. (2000) presented the partial imputation EM algorithm that imputes a part of missing data such that the observed data and imputed data form a monotone data pattern. In this paper, we propose the " ε-accelerated EM algorithm" that accelerates the convergence of the EM sequence using the vector ε algorithm of Wynn (1961). Section 2 describes the vector ε algorithm. In Section 3, we presented the ε-accelerated EM algorithm. Then we give the key results that the sequence generated by the ε-accelerated EM algorithm converges to the stationary point of the EM sequence and, for scalar sequences, it converges faster than the EM sequence. Section 5 examines the performance and properties of the ε-accelerated EM algorithm using numerical experiments.
2. The vector ε algorithm
The ε algorithm of Wynn (1961) is utilized to accelerate the convergence of,a slowly convergent scalar sequence. The algorithm has also been extended to vector sequences by Wynn (1962): It is known that the algorithm is very effective for linearly converging sequences. In this section, we illustrate the vector ε algorithm to obtain the proposed acceleration method for the EM algorithm. Let θ (י) denote a vector of dimension d that converges to a vector β − as t → ∞ . Let the inverse [ x ] − 1 of a vector x be defined by
in which 〈 x , x 〉 is the scalar product of x by itself. In general, the vector ε algorithm for a sequence
and then generates a vector ε ( t,k+1) by
For the case of k + 1 = 2r + 2, we have, the iteration form
from Eq. (1), see Brezinski and Zaglia (1991). For practical implementation, we apply the case of r = 0 to Eq. (2):
Then, from the initial conditions ε ( ι , 0 = θ ( ι ) and ε ( − − 2 ) = ∞ of Brezinski and Zaglia (1991), the iteration (3) becomes as follows:
Note that, at each iteration, the vector ε algorithm requires only O( d ) arithmetic operations while the Newton–Raphson and quasi-Newton algorithms are achieved at O( d 3 ) and O( d 2 ) and thus the computational cost is likely to become more expensive as d becomes large.
3. The ε-accelerated EM algorithm
Let y be observed data with a sample space Ω y and x be complete data augmented by y with a sample space Ω χ . We assume that there exists a many-to-one mapping from Ω x to Ω χ . Instead of observing x in Ω χ , we observe y = γ( x ) in Ω γ . Let f ( x | θ ) denote the probability density function of x and g ( y | θ ) denote the probability density function of y, where θ is an unknown parameter vector with a parameter space Θ . Define the conditional expectation of the log-likelihood function log f ( x | θ ) given y and θ  as
The EM algorithm chooses
at each iteration t = 1, 2, Given an initial value θ (0) ∈ Θ , the ε -accelerated EM algorithm incorporating the ε -acceleration process performs the following steps: E-step: Calculate M-step : Choose θ(t) such that for all θ ∈ Θ. ε-acceleration: Calculate from Eq. (4) and check the convergence using where ‖ x ‖∞ = max i {| x i |} for a vector x = { x i } and δ is a desired accuracy. Note that the ε -accelerated EM algorithm does not improve the E- and M-steps but it accelerates the convergence of the EM sequence using the ε-acceleration process. 4.  Convergence of the ε-accelerated EM algorithm The EM algorithm is a first-order successive substitution method and implicitly defines a map θ → M (θ) from Θ to Θ such that Suppose that θ ( 1 ) converges to a stationary point θ * and that M (θ) is differentiable at θ *. Expanding M (θ ( 1 ) in a Taylor series about θ *, Meng and Rubin (1994) gave the following equation:
where J* is the Jacobian matrix for M(θ) evaluated at θ*. For sufficiently large t , the distance between θ( 1 ) and θ* tends to be small, and then Eq. (5) becomes
where λ is the largest eigenvalue of J* (Schafer, 1997). For the parameter vector θ , the iterative sequence
where c is some constant, 0 < c < 1. The EM sequence converges linearly and λ corresponds to c . Now we can show the convergence of the ε-acceleration process. generated by the ε-accelerated EM algorithm converges to the stationary point θ* of the EM sequence. Proof.  See Appendix B. Next we investigate the speed of convergence of the ε -accelerated EM algorithm for a scalar sequence t ≥ 0 To compare the speed of convergence of the ε -accelerated EM algorithm with that of the EM algorithm, we provide the following notion (Brezinski and Zaglia, 1991).
t ≥ 0.
} θ (τ) t ≥ 0
and is identical to the Aitken δ² method of Aitken (1926). Traub (1964) proved that the Aitken δ² method accelerates the convergence of linear convergent sequences in the sense that
From the fact that the EM sequences converge linearly, we can see that the ε -accelerated EM algorithm accelerates the convergence of the EM sequence. Brezinski and Zaglia (1991) have completely described the convergence and acceleration properties of the vector ε algorithm, only for special classes of vector sequences. Thus, the mathematical consideration of the ε -accelerated EM algorithm is not trivial. Louis (1982) gave an acceleration formula that resembles Eq. (4) based on the multivariate version of the Aitken's acceleration method as follows:
where I d denotes the d × d identity matrix and J ( t − 1 ) is the Jacobian matrix for M ( θ ) evaluated at θ ( t − 1 ) . In order to J(t−1))−1 , Louis suggested making use of the equation
where After performing the E- and M-steps, Aitken's acceleration using Eq. (9) generates ğ ( 1 ) and then uses it in the next E- and M-steps. Meilijson (1989) pointed out that Louis's acceleration method using Eq. (9) is essentially equivalent to the Newton–Raphson algorithm in the neighborhood of MLEs. The ε-accelerated EM algorithm always uses acceleration (4), which does not depend on statistical models. By contrast, the Newton-type accelerations for the EM algorithm are needed to derive the acceleration formula for every statistical model. Thus, the ε -accelerated EM algorithm accelerates the convergence without affecting its simplicity and stability.
5. Numerical experiments
In this section, we provide numerical experiments to investigate how much faster the ε -accelerated EM algorithm converges than the EM algorithm. Note that the speed of convergence of the ε -accelerated EM algorithm does not depend on statistical models but on the EM sequence. To examine the performance of the ε-accelerated EM algorithm, we apply the algorithm to data sets for which the speed of convergence of the EM algorithm is quite slow. Example 1 ( Contingency Tables with partially classified observations ). Consider a 2 × 2 contingency table with completely and partially classified observations. Let X and Y be dichotomous variables and θ = be a set of joint probabilities of X and Y. Denote the cross-classified data of X and Y as n χy = {n χy (i, j)}i, j=1,2, and the partially classified data of X as nx = {nx(i)}i=1.2 and Y as ny = {ny(j)}j=1.2 . Assume that the data sets have a multinomial distribution with an unknown parameter θ . The data sets are shown in Table 1. For these data patterns, the convergence of the EM algorithm is quite slow, because its convergence is deeply associated with the proportion of missing data. In Table 2, we summarize the number of iterations for the EM and ε -accelerated EM algorithms for each δ = 10− 5–10− 8 and the data sets (a)–(e), and plot these results in Fig. 1. Table 3 gives the estimates of θ by the EM and ε -accelerated EM algorithms. As shown in Table 2 and Fig. 1, the EM algorithm increases linearly with the number of iterations as the data set changes from (a) to (e), while there is little variation in the number of iterations for the ε -accelerated EM algorithm and its convergence is significantly faster. For example, for δ = 10− 6 and the data set (d), the ε -accelerated EM algorithm takes only 41 iterations using 43 EM iterations to obtain the final values, while the EM algorithm takes 476 iterations Table 1
to achieve the same values. For other data sets, the EM algorithm requires the number of iterations roughly 3–10 times greater than that of the ε -accelerated EM algorithm. Table 3 demonstrates that the estimates by the EM algorithm for δ = 10− 5 match only three digits to the MLEs that are the estimates with δ = 10− 7 and 10− 8 . Then the ε-accelerated EM algorithm finds estimates that are identical to the MLEs. Next we investigate the speed of convergence of the ε-accelerated EM algorithm. In these numerical experiments, the ( i , j )th component-wise speeds of convergence of the EM and ε -accelerated EM algorithms are assessed as
If there exist MLEs of θ , θMLE i, j =1,2, and the ε -accelerated EM algorithm accelerates the convergence of the sequence { θ(t) > 0 , then it holds that Rij = 0 for all i and j . Fig. 2 shows the traces of converges to fMLE of the data sets (a)–(e) for each iteration with δ = 10− 6 , and illustrates that the sequence faster than   θ(t) t ≥ 0
μ = ( μ 1 , μ 2 ) and Σ = ( σ 11 , σ 22 , σ 12 ). The observed data are shown in Table 4, where the observed data from No. 1 to 3 are complete, and the data from No. 4, 5 and No. 6, 7 are incomplete with missing X 2 and X 1 , respectively. The EM algorithm fills in missing data of the incompletely observed data in the E-step and estimates μ and Σ in the M-step. For the data set (a) and δ = 10− 6, after 317 iterations of the EM algorithm, the MLEs of μ and Σ are obtained as follows: μ = (1.3005, 1.4163),    Σ = (0.2371, 4.9603, − 1.0478). The ε -accelerated EM algorithm finds the same values after 172 iterations using 174 EM iterations. Using both algorithms with the data set (b) and δ = 10− 6 , we obtain the same MLEs of μ and Σ : Σ = (70.1051, 79869.7113, 2182.2234). μ = (78.3977, 2247.1084), In this case the ε -accelerated EM algorithm takes 133 iterations, while the EM algorithm requires 312 iterations.
From these numerical experiments, we can also see that the sequence generated by the ε -accelerated EM algorithm converges to the MLEs faster than the EM sequence, as for the results of Example 1.
6. Concluding remarks
In this paper, we proposed the ε -accelerated EM algorithm that speeds up the convergence of the EM sequence incorporating the vector ε algorithm. The vector ε algorithm in itself is a fairly simple computational procedure and its computational cost is less expensive than that of Newton-type algorithms. The ε -accelerated EM algorithm does not require computation of the matrix inversion at each iteration and thus is numerically stable. Therefore, the ε -accelerated EM algorithm is an extension algorithm within the framework of the EM algorithm without affecting its simplicity, stability and flexibility. Theoretically we showed that the ε-accelerated EM algorithm is guaranteed to converge to the stationary point of the EM sequence and that it speeds up the convergence of the scalar EM sequence. Numerical experiments demonstrate that the ε-accelerated EM algorithm produces a sequence that converges to MLEs much faster than the sequence generated by the EM algorithm. Hence, the ε-accelerated EM algorithm finds sufficiently accurate estimates using a smaller number of EM iterations. In the future we intend to evaluate theoretically the speed of convergence of the ε -accelerated EM algorithm for vector cases. The authors would like to thank the editor and two referees whose valuable comments and kindly suggestions that led to an improvement of this paper, and Zhi Geng and Mingfeng Wang for helpful discussion on the revision of this paper. This research is supported by the Japan Society for the Promotion of Science (JSPS), Grant-in-Aid for Young Scientists, no. 16700264, and the Wesco Scientific Promotion Foundation. We here give the following equations used in the proof of Theorem 1. 1558 M. Kuroda, M. Sakakihara / Computational Statistics & Data Analysis 51 (2006) 1549–1561
and then
Thus, we can obtain
M. Kuroda, M. Sakakihara / Computational Statistics & Data Analysis 51 (2006) 1549–1561 because lim t →∞ θ( t ) = lim t →∞ θ( t +1) = θ* . Then, from Eq. (13), we have
and, from Eqs. (12) and (13), we also have
Thus, for t → ∞, we can obtain
Next we consider the case that some element of θ ( ) has converged, that is, θ ( ) − θ ∗ = 0 for some i ∈ {1, 2, . . . , d } and all t that are larger than some constant M . Let η ( ′ ) denote the ith element of η ( ′ ) . Then, for all t > M + 2, we obtain = 0 from Eq. (14). Since lim t →∞ θ ( t + 1) = θ * , we have
References Aitken, A.C., 1926. On Bernoulli's numerical solution of algebraic equations. Proc. Roy. Soc. Edinburg 46, 289–305. Brezinski, C., Zaglia, M.R., 1991. Extrapolation Methods: Theory and Practice. Elsevier Science Ltd. North-Holland, Amsterdam. Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm. J. Roy. Statist. Soc. Ser. B 39, 1–22. Geng, Z., Wan, K., Tao, F., 2000. Mixed graphical models with missing data and the partial imputation EM algorithm. Scand. J. Statist. 27, 433–444. Jamshidian, M., Jennrich, R.I., 1993. Conjugate gradient acceleration of the EM algorithm. J. Amer. Statist. Assoc. 88, 221–228. Lange, K., 1995. Acceleration of the EM algorithm by using Quasi-Newton methods. J. Roy. Statist. Soc. Ser. B 59, 569–587. Liu, C., Rubin, D.B., 1994. The ECME algorithm: a simple extension of EM and ECM with faster monotone convergence. Biometrika 81, 633–648. M. Kuroda, M. Sakakihara / Computational Statistics & Data Analysis 51 (2006) 1549–1561 Louis, T.A., 1982. Finding the observed information matrix when using the EM algorithm. J. Roy. Statist. Soc. Ser. B 44, 226–233. McLachlan, G.J., Krishnan, T., 1997. The EM Algorithm and Extensions. Wiley, New York. Meilijson, I., 1989. A fast improvement to the EM algorithm on its own terms. J. Roy. Statist. Soc. Ser. B 51, 127–138. Meng, X.L., Rubin, D.B., 1993. Maximum likelihood estimation via the ECM algorithm: a general framework. Biometrika 80, 267–278. Meng, X.L., Rubin, D.B., 1994. On the global and componentwise rates of convergence of the EM algorithm. Linear Algebra Appl. 199, 413–425. Schafer, J.L., 1997. Analysis of Incomplete Multivariate Data. Chapman & Hall, London. Traub, J.F., 1964. Iterative Methods for the Solution of Equations. Prentice-Hall Inc., Englewood Cliffs, NJ. Wynn, P., 1961. The epsilon algorithm and operational formulas of numerical analysis. Math. Comp. 15, 151–158. Wynn, P., 1962. Acceleration techniques for iterated vector and matrix problems. Math. Comp. 16, 301–322.