Improving the vector $\varepsilon$ acceleration for the EM algorithm using a re-starting procedure

Abstract

The expectation–maximization (EM) algorithm is a popular algorithm for finding maximum likelihood estimates from incomplete data. However, the EM algorithm converges slowly when the proportion of missing data is large. Although many acceleration algorithms have been proposed, they require complex calculations. Kuroda and Sakakihara (Comput Stat Data Anal 51:1549–1561, 2006) developed the $\varepsilon$-accelerated EM algorithm which only uses the sequence of estimates obtained by the EM algorithm to get an accelerated sequence for the EM sequence but does not change the original EM sequence. We find that the accelerated sequence often has larger values of the likelihood than the current estimate obtained by the EM algorithm. Thus, in this paper, we try to re-start the EM iterations using the accelerated sequence and then generate a new EM sequence that increases its speed of convergence. This algorithm has another advantage of simple implementation since it only uses the EM iterations and re-starts the iterations by an estimate with a larger likelihood. The re-starting algorithm called the $\varepsilon$R-accelerated EM algorithm can further improve the EM algorithm and the $\varepsilon$-accelerated EM algorithm in the sense of that it can reduces the number of iterations and computation time.

Appendices

Appendix 1: The vector $\varepsilon$ algorithm

Let $\theta ^{(t)}$ denote a vector of dimensionality $d$ that converges to a vector $\theta ^{(\infty )}$ as $t \rightarrow \infty$. Let the inverse $[\theta ]^{-1}$ of a vector $\theta$ be defined by

$$\begin{aligned}{}[\theta ]^{-1} = \frac{\theta }{\Vert \theta \Vert ^{2}}, \end{aligned}$$

where $|| \theta ||$ is the Euclidean norm of $\theta$.

In general, the vector $\varepsilon$ algorithm for a sequence $\{\theta ^{(t)}\}_{t \ge 0}$ starts with

$$\begin{aligned} \varepsilon ^{(t,-1)}=0, \qquad \varepsilon ^{(t,0)}=\theta ^{(t)}, \end{aligned}$$

and then generates a vector $\varepsilon ^{(t,k+1)}$ by

$$\begin{aligned} \varepsilon ^{(t,k+1)} = \varepsilon ^{(t+1,k-1)}+ \left[ \varepsilon ^{(t+1,k)}- \varepsilon ^{(t,k)} \right] ^{-1}, \qquad k=0,1,2,\ldots . \end{aligned}$$

(14)

For practical implementation, we apply the vector $\varepsilon$ algorithm for $k=1$ to accelerate the convergence of $\{\theta ^{(t)}\}_{t \ge 0}$. From Eq. (14), we have

$$\begin{aligned} \varepsilon ^{(t,2)}&= \varepsilon ^{(t+1,0)}+ \left[ \varepsilon ^{(t+1,1)}- \varepsilon ^{(t,1)} \right] ^{-1} \quad \text{ for } k=1, \\ \varepsilon ^{(t,1)}&= \varepsilon ^{(t+1,-1)}+ \left[ \varepsilon ^{(t+1,0)}- \varepsilon ^{(t,0)} \right] ^{-1} = \left[ \varepsilon ^{(t+1,0)}- \varepsilon ^{(t,0)} \right] ^{-1} \quad \text{ for } k=0. \end{aligned}$$

Then the vector $\varepsilon ^{(t,2)}$ becomes as follows:

$$\begin{aligned} \varepsilon ^{(t,2)}&= \varepsilon ^{(t+1,0)} + \left[ \left[ \varepsilon ^{(t,0)}-\varepsilon ^{(t+1,0)}\right] ^{-1} +\left[ \varepsilon ^{(t+2,0)}-\varepsilon ^{(t+1,0)} \right] ^{-1} \right] ^{-1} \\&= \theta ^{(t+1)} + \left[ \left[ \theta ^{(t)}-\theta ^{(t+1)} \right] ^{-1} +\left[ \theta ^{(t+2)}-\theta ^{(t+1)}\right] ^{-1} \right] ^{-1}. \end{aligned}$$

Appendix 2: Pseudo-code of the $\varepsilon$R-accelerated EM algorithm

Initialization

We set the initial value of the EM step $\theta _{0}$, the desired precision $\delta$, the threshold $\delta _{Re} (>\delta )$ and the size of decrement $10^{-k}$ and determine the maximum number of iterations ($itrmax$).

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 algorithm is that the change in $\theta '$ after an EM iteration $\tilde{\mathbf {g}}(\theta ') = M(\theta ')-\theta '$ 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 $\tilde{\mathbf {g}}(\theta ')$ is given by

$$\begin{aligned} \tilde{\mathbf {g}}(\theta ') = M(\theta ')-\theta ' \approx - \left. \frac{\partial ^{2} Q(\theta |\theta ')}{\partial \theta \partial \theta ^{\top }}^{-1} \right| _{\theta =\theta '} \left. \frac{\partial L_{o}(\theta )}{\partial \theta } \right| _{\theta =\theta '}. \end{aligned}$$

In the AEM algorithm, first the EM algorithm runs until the difference between $2 L_{o}(\theta ^{(t)})$ and $2 L_{o}(\theta ^{(t-1)})$ falls below one. Then the CG accelerator updates the EM estimate $\theta ^{(t)}$ for obtaining $\theta ^{(t+1)}$:

• 1. Set $t=0$ and $\mathbf {d}^{(t)}=\tilde{\mathbf {g}}(\theta ^{(t)})$.

• 2. Find $\alpha ^{(t)}$ to maximize $L_{o}(\theta ^{(t)}+\alpha \mathbf {d}^{(t)})$ using a line search algorithm.

• 3. Update $\theta ^{(t+1)} = \theta ^{(t)}+\alpha ^{(t)} \mathbf {d}^{(t)}$.

• 4. Compute

  $$\begin{aligned} \tilde{\mathbf {g}}(\theta ^{(t+1)})&= M(\theta ^{(t+1)})-\theta ^{(t+1)}, \\ \gamma ^{(t)}&= \frac{ \{\mathbf {g}(\theta ^{(t+1)})-\mathbf {g}(\theta ^{(t)})\}^{\top } \tilde{\mathbf {g}}(\theta ^{(t+1)}) }{ \{\mathbf {g}(\theta ^{(t+1)})-\mathbf {g}(\theta ^{(t)})\}^{\top } \mathbf {d}^{(t)}}, \end{aligned}$$

  where $\mathbf {g}(\theta )$ is the gradient of $L_{o}(\theta )$.

• 5. Update $\mathbf {d}^{(t+1)}=\tilde{\mathbf {g}}(\theta ^{(t+1)})-\gamma ^{(t)} \mathbf {d}^{(t)}$ and then $t=t+1$.

• 6. Repeat Steps 2 to 5 until $|| \tilde{\mathbf {g}}(\theta ^{(t)}) ||^{2} \le \delta$.