Why use the EM Algorithm and not just marginalise the complete likelihood?

Question

On the wikipedia article for Expectation-Maximization it states

Given the statistical model which generates a set $\mathbf{X}$ of observed data, a set of unobserved latent data or missing values $\mathbf{Z}$, and a vector of unknown parameters $\boldsymbol\theta$, along with a likelihood function $L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z}) = p(\mathbf{X}, \mathbf{Z}\mid\boldsymbol\theta)$, the maximum likelihood estimate (MLE) of the unknown parameters is determined by maximizing the marginal likelihood of the observed data $$L(\boldsymbol\theta; \mathbf{X}) = p(\mathbf{X}\mid\boldsymbol\theta) = \int p(\mathbf{X},\mathbf{Z} \mid \boldsymbol\theta) \, d\mathbf{Z} = \int p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol\theta) p(\mathbf{Z} \mid \boldsymbol\theta) \, d\mathbf{Z} $$ However, this quantity is often intractable since $\mathbf{Z}$ is unobserved and the distribution of $\mathbf{Z}$ is unknown before attaining $\boldsymbol\theta$.

I don't understand this last sentence. Surely $\mathbf{Z}$ being unobserved is why we integrate over $\mathbf{Z}$, and we set values of $\boldsymbol{\theta}$ during the optimisation procedure which then allows us to compute $p(\mathbf{Z} \mid \boldsymbol{\theta})$?

I initially thought it was intractable because there are too many settings of $\mathbf{Z}$ to consider, but in the expectation step of the EM algorithm we calculate \begin{align} Q(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{(t)}) &= \mathbb{E}_{\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{(t)}}\left[\log L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z})\right]\\ &= \int_Z p(\mathbf{Z} \mid \mathbf{X}, \boldsymbol{\theta}^{(t)})\log L(\boldsymbol\theta; \mathbf{X}, \mathbf{Z})d\mathbf{Z} \end{align} which doesn't look any easier than the first integral. If we can approximate this expectation using some method like monte-carlo, why not just approximate the first integral instead?

Answer 1

So I've discussed this with a colleague.

Consider the marginalisation $$p(\mathbf{X} \mid \boldsymbol{\theta}) = \int p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol{\theta})p(\mathbf{Z} \mid \boldsymbol{\theta}) d\mathbf{Z}.$$

This can be rewritten as the expectation $$\mathbb{E}_{\mathbf{Z} \mid \boldsymbol{\theta}}\left[ p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol{\theta}) \right].$$

If this can be calculated exactly, we're fine. However, if the expectation is intractable, we need to compute a numerical approximation. Hence to maximise w.r.t. $\boldsymbol{\theta}$ we will need to use some iterative procedure like gradient-ascent.

Suppose we have the current estimate $\boldsymbol{\hat{\theta}^{(t)}}$ for $\boldsymbol{\theta}$. If we now try and approximate the expectation by sampling from $\mathbf{Z}$, we get $$\frac{1}{n}\sum_{i=1}^n p(\mathbf{X} \mid \mathbf{Z}_i, \boldsymbol{\theta})$$ but the $\mathbf{Z}_i$ were sampled from the distribution $\mathbf{Z} \mid \boldsymbol{\hat{\theta}^{(t)}}$, so this actually approximates the expectation $$\mathbb{E}_{\mathbf{Z} \mid \boldsymbol{\hat{\theta}^{(t)}}}\left[ p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol{\theta}) \right] \approx \frac{1}{n}\sum_{i=1}^n p(\mathbf{X} \mid \mathbf{Z}_i, \boldsymbol{\theta})$$ which is not the same as the original expectation we wanted. Hence if we maximised this new expectation w.r.t. $\boldsymbol{\theta}$ we won't be doing Maximum Likelihood Estimation.

The EM algorithm therefore takes a different approach.

Answer 2

This is a few years old but I don't believe the current answer really covers it.

What the question boils down is whether $\mathbb{E}_{p(z\mid x;\theta')}[\log p(x, z\,;\, \theta')]$ (in EM lower-bound) is tractable even when $\log p(x ; \theta)=\log(\mathbb{E}_{p(z;\theta)}[ p(x \mid z\,;\,\theta)])$ (standard marginal likelihood) is not. The answer is often Yes. In your question you had left the logarithm out of the marginal likelihood, but that is a key observation; in general $\int \log p(x, z) dz$ may be tractable while $\log (\int p(x, z) dz)=\log p(x)$ is not -- one way to think about this is that sum of logs is much easier than log of sums. The point of EM is to lower-bound this intractable marginal likelihood into a lower-bound that has, possibly, tractable parts.

The EM algorithm is especially popular with models where $p(x, z)$ is in the exponential family, but the marginal $p(x)$ is an intractable non-exponential family distribution. This situation includes household names such as mixture of Gaussians, hidden Markov models, Gaussian state-space models etc. To see the tractability of the lower-bound term in exponential families, consider $\log p(x, z)$ in a general exponential family: \begin{align} \log p(x, z; \theta) = \langle \eta(\theta), t(x, z)\rangle - \log Z(\theta) + constant \end{align}
and therefore: \begin{align} \mathbb{E}_{p(z\mid x;\theta)}[\log p(x, z)] = \langle \eta(\theta), \mathbb{E}_{p(z\mid x;\theta)}[t(x, z)]\rangle - \log Z(\theta) + constant \end{align}
all we need are the expected sufficient statistics under the posterior which is, usually, tractable with conjugate exponential families. In contrast, as mentioned above $\log p(x ; \theta)$ is typically not tractable and nor even exponential family, hence the motivation for the whole EM algorithm.

---

Addition: the reason (related to the other answer) learning parameters by sampling is tricky:

$\nabla_{\theta}\log \mathbb{E}_{p(z;\theta)}[p(x \mid z; \theta)]=\frac{\nabla_{\theta}\mathbb{E}_{p(z;\theta)}[p(x \mid z; \theta)]}{\mathbb{E}_{p(z;\theta)}[p(x \mid z; \theta)]}$

ignoring the challenge of sampling the denominator, taking gradients of the numerator alone is tricky already as it doesn't only involve gradients $\nabla_{\theta}p(x \mid z;\theta)$, but instead you have to take gradients through the monte-carlo sample: \begin{align} \nabla_{\theta}\mathbb{E}_{p(z;\theta)}[p(x \mid z; \theta)]&=\int p(x \mid z; \theta)\nabla_{\theta}p(z;\theta)dz + \int p(z;\theta)\nabla_{\theta}p(x \mid z; \theta)dz \\ &= \mathbb{E}_{p(z;\theta)}[p(x \mid z; \theta)\nabla_{\theta}\log p(z;\theta)]+\mathbb{E}_{p(z;\theta)}[\nabla_{\theta} p(x \mid z; \theta)] \\ &\approx \frac{1}{M}\sum_{i=1}^M p(x \mid z_i; \theta)\nabla_{\theta}\log p(z_i;\theta) + \frac{1}{M}\sum_{i=1}^M \nabla_{\theta} p(x \mid z_i; \theta) \end{align} Usually these gradients have a very large variance (there are some tricks to alleviate such issues, for example the reparametrization trick). However, I believe if this was easy to get these gradients accurately I would think this approach could work, which would contradict the earlier answer.