1 Introduction
The EM algorithm has been a general and popular algorithm for finding maximum likelihood estimates (MLEs) from incomplete data since Dempster, Laird and Ru- bin [DLR77]. However, it is often criticized that the convergence of the EM algorithm is slow when the proportion of missing data is high. In order to improve the speed of convergence of the EM algorithm, various algorithms incorporating optimization algorithms with faster convergence speed have been proposed, see Krishnan and McLachlan [MK97]. Applying the accelerator based on the Newton-type algorithms, it requires the computation of a matrix inversion at each iteration. Then its compu- tation is likely to become rapidly complicated as the number of parameters increases, and is also expected numerical instabilities. Therefore the Newton-type accelerators for the EM algorithm are lost the attractive features of the EM algorithm, such as its stability, flexibility and simplicity. Kuroda and Sakakihara [KS05] proposed the ε-accelerated EM algorithm that is the fairly simple computational procedure and speeds up the convergence of EM sequences using the vector ε algorithm of Wynn [Wy62], and demonstrated that the algorithm produces a sequence to converge to the MLEs much faster than the EM sequence in the numerical experiments. In this paper, we apply the ε-accelerated EM algorithm to the maximum likeli- hood estimation for log-linear models with missing data. Moreover we propose the acceleration of the EM algorithm using the Aitken δ 2 method of Aitken [A26] to ac- celerate the convergence of scalar sequences of iterates. For some log-linear models, there are no closed-form to compute the MLEs. The ECM algorithm of Meng and Rubin [MR93] is also used to the models. Then the ε- and Aiken δ 2 -accelerators are applicable to the ECM algorithms. In Section 2 we introduce log-linear models for contingency tables, and, in Section 3, show the EM and ECM algorithms for the log-linear models with missing data. In Section 4, we provide the the ε- and Aitken δ 2 -accelerations for the EM and ECM algorithms, and give the fundamental theories for the convergence of these acceleration algorithms.
2 Log-linear models with missing data
Let Xv = (X 1 , . . , X k ) be a k-dimensional discrete random vector and V = {1,...,k} be the index set of Xv. We also denote a finite set of values of Xi as Ω i , i ∈ V , and the space of possible values of X V as the Cartesian product ΩV = Π iEV Ω i. Associated with each i E V, we shall have a random variable X i taking values in a sample space Ω . For a subset A ⊂ V , we write X A for { X i | i ∈ A } and Ω = Π ic A Ω . Let pv(xv) denote the cell probability of Xv  = xv  E  Ωv  and let θ  = {pv(xv)|xv E Ωv} be the set of cell probabilities. The marginal probability of X Λ = x Λ ∈ Ω Λ for A ⊂ V can be also calculated by
where the symbol "\" denotes the operator of a difference set. In this paper, we consider the class of hierarchical log-linear models for contin- gency tables in the case of multinomial sampling. An hierarchical log-linear model is specified by a generating class E = {e 1, . . , e 1, . . , e 1, . . , e 1, . . , e 1, . . , e 1, . . , e 1, . . , e 1, . . , e 1, . , e 1, . , e 1, in minimal interaction terms, and e ∈ E is called a generator. For example, considen the log-linear model with the generating class E = {{1,2}, {1,3}, {2,3}}. Then the log-linear model is represented by
+ u12(x {1,2}) + u12(x {1,3}) + u23(x {2,3}).
In the cases with missing data, the completely classified data are represented in the full contingency tables and the partially classified data are represented in supplemental tables. These tables are identified with an observing pattern T = t s } where each configuration t is a set of observed variables for a sup- plemental table. The cell counts of the contingency table classified by observing variables Xt are also denoted as {nt(xt)|xt E Ωt}. Assume that the missing data are missing at random in the sense that Rubin [R76].
3 The EM and ECM algorithms for log-linear models
with missing data
We show the EM algorithm to find the MLE of the parameter vector θ for the log- linear model with a generating class E and an observing pattern T = {V, t 1 , . . , t s}. or the case that the MLE of θ is given by closed-forms, the EM algorithm repeats wing two steps till the convergence to the desired accuracy is obtained: tep: Calculate the expected marginal cell counts {ñ e ( x e )} for each generator
where C = (c 1 , . . , c M –1) is the set of separators of E and ν(c) is an index that counts the number of times of each separator c in C . It is known that certain log-linear models do not have closed-form MLEs, such as the no three-interaction model for a three-way contingency table given by equation (1). Then the ECM algorithm of Meng and Rubin [MR93] is applicable to the models. The ECM algorithm finds the MLEs by replacing a complicated M-step of the EM algorithm with several computationally simpler CM-step: CM-step: Obtain θ (r+1) by performing the Iterative Proportional Fitting proce- dure of Bishop, Fienberg and Holland [BFH77] as follows: pv (x v ) (r+(i-1)/M) ñ c i (x c i ) (r) pv ( x ) ( r + i M )  = pc (x c i ) (r + (i − 1) / M) N for i = 1, 2, . . , M .
4 Acceleration algorithms
In order to accelerate the convergence of the EM and ECM algorithms, we incorpo- rate the ε- and Aitken δ²-accelerators into them. The first accelerator is the vector ε algorithm of Wynn [Wy62] and accelerates the convergence of vector sequences. The second accelerator is the Aitken δ 2 method of Aitken [A26] and accelerates the convergence of scalar sequences.
4.1 The ε acceleration
Let θ be a d-dimensional vector. Assume that a sequence {θ(r)}r>o converges to a vector θ* as r → ∞. Define the inverse [ x ] − ¹ of a vector x by
in which < x, x > is the scalar product of x by itself. In general, the vector ε algorithm for a sequence {θ(r)},r>0 starts with
and then generates a vector ε(r,k+1) by
For the case of k + 1 = 2l + 2, we have the iteration form such as
from equation (2), see Brezinski and Zaglia [BZ91]. For practical implementation, we apply the case of l = 0 to equation (??). Then, from initial conditions ε(r,o) = θ(r) and ε (r,−2) = ∞ oo of Brezinski and Zaglia [BZ91], the iteration becomes as follows:
To accelerate the convergence of the EM and ECM algorithms, we apply the ε -acceleration process to the sequence {θ(t)} t>o generated by the M- or CM-steps: ● ε-acceleration: Generate a vector θ(r) by from equation (1) and check the convergence to a desired accuracy. Then the acceleration algorithm via the vector ε algorithm does not improve the E- and M-steps or the CM-step but accelerates the convergence of the sequence {θ(r)}r>0 using the ε-acceleration process.
4.2 Convergence of the ε-acceleration process
The EM algorithm is a first-order successive substitution method and then implicitly defines a map θ → M(θ) from the parameter space to itself such that
Expanding M(θ(r)) in a first term of a Taylor series of about θ(r-1), the following approximation holds:
where J(r-1) is the Jacobian matrix for the mapping M(0) evaluated at 0(r-1). Suppose θ(r) converges to a stationary point θ*. For r large enough, we have J(r) = J* and thus
where J* is J evaluated at a stationary point θ*. Then J* determines the rate of convergence of the EM algorithm near θ*. In fact, equation (4) becomes
because the largest eigenvalue λ of J* dominates the convergence, see Louis [L82]. Kuroda and Sakakihara [KS05] provided the convergence of the ε-acceleration process. Theorem 1. The sequence { θ(r )} r>o generated by the ε-acceleration process for the EM sequence converges to the stationary point θ* of the EM sequence. Like the EM algorithm, equation (5) also holds for each iteration of the ECM algorithm when r tends to be large. Thus Theorem 1 is valid for ECM sequences. Corollary 1. The sequence { θ( r ) } r>o generated by the ε-acceleration process for the ECM sequence converges to the stationary point θ° of the ECM sequence. For the convergence and acceleration properties of the vector ε algorithm, Brezin- ski and Zaglia [BZ91] has completely described only for special classes of vector sequences. Thus the mathematical consideration of the ε-acceleration process is less trivial.
4.3 The Aitken δ²-acceleration
Let a scalar sequence {φ{r}}r>>o converges to a limit φ *. For a scalar sequence {φ{r}}r>o, the Aitken δ² algorithm generates a sequence {φ{r}}r>o by
Note that the ε algorithm for scalar sequences is identical to the Aitken δ 2 algorithm. For the sake of simplicity, we re-denote θ = {pv(xv)|xv E Sv} as θ = (θ 1 , θ 2 , . . , θ d ) where θ d = 1 − ∑ j −1 θ j . We can find an alternative parameterization φ = φ(θ) such as
Then the MLE of φ(θ) is φ(θ) that is the function evaluated at t the MLE θ of θ, because φ(θ) is a one-one function of θ. For the parameterization of φ = φ(θ) in equation (7), the property of the invariance of MLEs holds. In order to accelerate the convergence of the EM or ECM sequence {θ(r)}r>o, the Aitken δ²-acceleration process generates a sequence {φ ( x )} r>o: ● Aitken δ 2 -acceleration: Generate a vector φ (r) = (φ i r) i =1,...,a from
and check the convergence to a desired accuracy. After the convergence of the Aitken δ 2 acceleration, we compute θ from the limit ù* of { ì(r } r >o using

4.4 Convergence of the Aitken δ²-acceleration process
Next we give the convergence properties of the Aitken δ 2 acceleration for the EM and ECM algorithms. For a scalar sequence {φ{τ}}r≥o generated by equation (6), Traub [T64] provided the following lemma: Lemma 1. If {φ(r)}r>0 converges to a stationary point φ* as r → ∞, then {φ(r)}r>0 generated by equation (6) converges to the same stationary point φ*. Then we have the following lemma: Lemma 2. If {<p(r)}r>o is a convergent vector sequence to the vector <p*, then ||ç(r) — ç * ||ç  converges  to  zero as r  → ∞, moreover  |ç {r}  — ç !  | converges  to  zero as r → ∞ , where ||ϕ||∞ = maxi{|φi|}. From Lemmas 1 and 2, we obtain the result: Theorem 2. If ϕ (r) = (φ (r)) i=1,...,a is a convergent vector sequence to the vector ϕ* = (Φ i )i=1,...,a and the vector sequence φ(r) = (φ(r))i=1,...,a is generated by equa- tion (6), then φ(r) → φ* as r → ∞. Furthermore, from the invariance of MLEs of ϕ = ϕ(θ) in (7) and Theorem 2, we can give the following result: Theorem 3. If the EM sequence {θ(τ)} r>0 converges to the MLE δ, then θ obtained by equation (??) is fi. To compare the speed of convergence of the Aitken δ 2 acceleration for the EM al- gorithm with that of the EM algorithm, we provide the following notion of Brezinski and Zaglia [BZ91]. r>o be a scalar sequence obtained by applying an extrapo- r>0. Assume that lim t →∞∞ φ(r) = lim t →∞∞ φ(r) = φ*. If say that the sequence {φ(r)}r≥0 converges to φ* faster than {φ(r)}r≥0 or the ion method accelerates the convergence of {φ(r)}r≥0. proved that Aitken δ 2 method accelerates the convergence of >o in the sense that
cts of Theorem 3 and equation (8) give us the validity of the Aitken δ 2 acceleration for the EM algorithm.
nowledgement
The authors would like to thank a referee for valuable comments. This is orted by the Japan Society for the Promotion of Science (JSPS), Grant-in-Aid Young Scientists, No 16700264, and the Wesco Scientific Promotion Foundation.