Expectation-Maximization (EM) Algorithm

There are times when we want to apply the ML or MAP estimation technique, but the data ${\bf x}$ is incomplete. This implies that the model consists of two parts: (a) the observed part: ${\bf x}$ and (b) the hidden part: ${\bf y}$. We can associate RVs $X$ and $Y$ with the observed and hidden parts, respectively. We can still apply ML or MAP estimation techniques if we assume a certain joint PDF $P (X,Y)$ between the observed and hidden RVs, and then marginalize over the hidden RVs $Y$. Marginalization of an RV $Y$ chosen from a set of RVs refers to the process of integration of the joint PDF over the values $y$ of the chosen RV. This is the key idea behind the EM algorithm.

Considering ML estimation, for example, we compute the optimal parameter as

$$\theta^{*} = \mathop{\mbox{argmax }}_{\theta} \Bigg( \log L (\theta \vert {\bf x}) \Bigg)$$

$$= \mathop{\mbox{argmax }}_{\theta} \Bigg( \log \bigg( \int_{\mathcal{S}_Y} P ({\bf x}, y \vert \theta) dy \bigg) \Bigg),$$ (37)

where $L (\cdot)$ is the likelihood function described previously in Section 2.3.1. This key idea is formalized in the expectation-maximization (EM) algorithm [43,104].

Herman O. Hartley [71] pioneered the research on the EM algorithm in the late 1950s. The first concrete mathematical foundation, however, was laid by Dempster, Laird, and Rubin [43] in the late 1970s. Neal and Hinton [111,112,108] presented the EM algorithm from a new perspective of lower-bound maximization. Over the years, the EM algorithm has found many applications in various domains and has become a powerful estimation tool [104,48].

The EM algorithm is an iterative optimization procedure. Starting with an initial parameter estimate $\theta^0$, it is guaranteed to converge to the local maximum of the likelihood function $L (\theta \vert {\bf x})$. The EM algorithm consists of two steps: (a) the E step or the expectation step and (b) the M step or the maximization step.

• 
  The E step constructs an optimal lower bound $B (\theta)$ to the log-likelihood function $\log L (\theta \vert {\bf x})$. This optimal lower bound is a function of $\theta$ that touches the log-likelihood function at the current parameter estimate $\theta^i$, i.e.,
  

  $$B (\theta^i) = \log L (\theta^i \vert {\bf x}),$$ (38)

  
  
  and never exceeds the objective function at any $\theta$, i.e.,
  

  $$\forall \theta \in (\-\infty,\infty): B (\theta) \le \log L (\theta \vert {\bf x}).$$ (39)

  
  
  Intuitively, maximizing this optimal lower bound $B (\theta)$ (in the M step) will surely take us closer to the maximum of the log-likelihood function $\log L (\theta \vert {\bf x})$, i.e., the ML estimate. We compute this optimal lower bound as follows [108,42,104].

  Let us rewrite the log-likelihood function as
  

  $$\log L (\theta \vert {\bf x}) = \log P ({\bf x} \vert \theta)$$

  $$= \log \int_{y} P ({\bf x}, y \vert \theta) dy$$

  $$= \log E_{f(Y)} \Bigg[ \frac { P ({\bf x}, Y \vert \theta) } { f (Y) } \Bigg],$$ (40)

  
  
  where $f (Y)$ is any arbitrary PDF. Applying Jensen's inequality [34], and using the concavity of the $\log (\cdot)$ function, gives:
  

  $$\log E_{f(Y)} \Bigg[ \frac { P ({\bf x}, Y \vert \theta) } { f (Y) } \Bigg] \ge E_{f(Y)} \Bigg[ \log \frac { P ({\bf x}, Y \vert \theta) } { f (Y) } \Bigg]$$

  $$\equiv B (\theta).$$ (41)

  
  

  Our goal is to try to find the particular PDF $f (Y)$ such that $B (\theta)$ is the optimal lower bound that touches the log-likelihood function at the current parameter estimate $\theta^i$. We can achieve this goal by solving the following constrained-optimization [137] problem:
  

  $$\mathrm {Maximize } B (\theta^i)$$

  $$\mathrm {with respect to } f (Y)$$

  $$\mathrm {under the constraint } \int_y f(y) dy = 1.$$ (42)

  
  
  Using the Lagrange-multiplier [137] approach, the objective function to be maximized is
  

  $$J (f(Y)) = B (\theta^i) + \lambda \Bigg( 1 - \int_y f(y) dy \Bigg).$$ (43)

  
  
  The derivative of the objective function $J (f(Y))$ with respect to $f(y)$ is
  

  $$\frac { \partial J } { \partial f(y) } = - \lambda + \int_{y} \log P ({\bf x}, y \vert \theta^i) dy - \bigg( 1 + \log f (y) \bigg).$$ (44)

  
  
  The derivative of the objective function $J (f(Y))$ with respect to $\lambda$ is
  

  $$\frac { \partial J } { \partial \lambda } = \int_y f(y) dy - 1.$$ (45)

  
  
  The objective function achieves its maximum value when both the aforementioned derivatives in (2.44) and (2.45) are zero. Using these conditions, and after some simplification, we get
  

  $$f (y) = \frac { P ({\bf x}, y \vert \theta^i) } { P ({\bf x} \vert \theta^i) }$$

  $$= P (y \vert {\bf x}, \theta^i).$$ (46)

  
  
  This gives our optimal lower bound as
  

  $$B (\theta) = \int_y P (y \vert {\bf x}, \theta^i) \log \frac { P ({\bf x}, y \vert \theta) } { P (y \vert {\bf x}, \theta^i) } dy$$ (47)

  
  
  We can confirm that $B (\theta^i)$ indeed equals $\Big( \log P ({\bf x} \vert \theta^i) \Big)$, which indicates that $B (\theta)$ touches the log-likelihood function $\Big( \log P ({\bf x} \vert \theta) \Big)$ at $\theta^i$ and is an optimal lower bound.
  
• 
  The M step performs the maximization of the function $B (\theta)$ with respect to the variable $\theta$.
  

  $$\mathop{\mbox{argmax }}_{\theta} B (\theta) = \mathop{\mbox{argmax }}_{\theta} \int_y P (y \vert {\bf x}, \thet... ...\log \frac { P ({\bf x}, y \vert \theta) } { P (y \vert {\bf x}, \theta^i) } dy$$

  $$= \mathop{\mbox{argmax }}_{\theta} \int_y P (y \vert {\bf x}, \theta^i) \log P ({\bf x}, y \vert \theta) dy$$

  $$= \mathop{\mbox{argmax }}_{\theta} \int_y P (y \vert {\bf x}, \theta^i) \log P ({\bf x}, y \vert \theta) dy$$

  $$= \mathop{\mbox{argmax }}_{\theta} Q (\theta)$$ (48)

  
  
  where the Q function is
  

  $$Q (\theta) = E_{P (Y \vert {\bf x}, \theta^i)} \Big[ \log P ({\bf x}, Y \vert \theta) \Big]$$ (49)

  $$= \int_{\mathcal{S}_Y} P (y \vert {\bf x}, \theta^i) \log P ({\bf x}, y \vert \theta) dy.$$ (50)

  
  
  Observe that $Q (\theta)$ also depends on the current parameter estimate $\theta^i$ that is considered a constant. The M step assigns the new parameter estimate $\theta^{i+1}$ as the one that maximizes $Q (\theta)$, i.e.,
  

  $$\theta^{i+1} = \mathop{\mbox{argmax }}_{\theta} Q (\theta).$$ (51)

  
  

The iterations proceed until convergence to a local maximum of $L (\theta \vert {\bf x})$. Actually, the M step need not find that $\theta^{i+1}$ corresponding to the maximum value of the Q function, but rather it is sufficient to find any $\theta^{i+1}$ such that

$$Q (\theta^{i+1}) \ge Q (\theta^i).$$ (52)

This modified strategy is referred to as the generalized-EM (GEM) algorithm and is also guaranteed to converge [43].