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

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?

• Can't we just iteratively sample given the current estimates and maximize to get a new estimate? Sorry for the new post, I can't post comments yet Jul 23, 2023 at 0:04

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.

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.