Expectation Maximization and Mixture Models

Expectation Maximization

Suppose our model contains an observed variable $x$ and a latent variable $z$ that we do not observe directly. The joint distribution is:

p(x,z;\theta)

To find the probability of the observed data $x$ , we sum over all possible values of $z$ :

p(x;\theta) = \sum_z p(x,z;\theta)

For a continuous latent variable, the sum is replaced by an integral:

p(x;\theta) = \int p(x,z;\theta)\,dz

We use sums below for notational simplicity.

For a dataset of $n$ observations, the log-likelihood is:

\ell(\theta) = \sum_{i=1}^n \log p(x^{(i)};\theta)

= \sum_{i=1}^n \log \sum_{z^{(i)}} p(x^{(i)},z^{(i)};\theta)

The difficulty is the sum over the latent variable inside the logarithm:

\log \sum_z p(x,z;\theta)

This often makes the likelihood difficult to optimize directly. EM handles this by constructing a lower bound on the log-likelihood and optimizing that bound instead.

The Evidence Lower Bound

To construct this lower bound, introduce a distribution $Q(z)$ over the latent variable, where:

\sum_z Q(z)=1

and $Q(z)\geq0$ .

We assume $Q(z)>0$ wherever $p(x,z;\theta)>0$ , so the ratio below is well defined.

Jensen's Inequality

Jensen's inequality states that for a convex function $f$ :

\mathbb{E}[f(X)] \geq f(\mathbb{E}[X])

For a concave function, the inequality is reversed:

\mathbb{E}[f(X)] \leq f(\mathbb{E}[X])

Also, if $X$ is constant, then:

X=\mathbb{E}[X]

and therefore:

\mathbb{E}[f(X)] = f(\mathbb{E}[X])

Using Jensen's inequality:

\log p(x;\theta) = \log \sum_z Q(z) \frac{p(x,z;\theta)}{Q(z)}

\geq \sum_z Q(z) \log \frac{p(x,z;\theta)}{Q(z)}

The expression on the right is called the Evidence Lower Bound, or ELBO.

See derivation

\log p(x;\theta) = \log \sum_z p(x,z;\theta)

= \log \sum_z Q(z) \frac{p(x,z;\theta)}{Q(z)}

= \log \mathbb{E}_{z\sim Q} \left[ \frac{p(x,z;\theta)}{Q(z)} \right]

\geq \mathbb{E}_{z\sim Q} \left[ \log \frac{p(x,z;\theta)}{Q(z)} \right]

= \sum_z Q(z) \log \frac{p(x,z;\theta)}{Q(z)}

The bound becomes tight when the quantity inside the expectation is constant with respect to $z$ :

\frac{p(x,z;\theta)}{Q(z)} = c

Solving this condition for $Q$ gives:

Q(z) = p(z\mid x;\theta)

See derivation

Starting from:

\frac{p(x,z;\theta)}{Q(z)} = c

we get:

Q(z) = \frac{p(x,z;\theta)}{c}

Since:

\sum_z Q(z)=1

the right-hand side must also sum to $1$ :

\sum_z \frac{p(x,z;\theta)}{c} = 1

Therefore:

c = \sum_z p(x,z;\theta)

= p(x;\theta)

So:

Q(z) = \frac{p(x,z;\theta)} {p(x;\theta)}

= p(z\mid x;\theta)

We can verify that this makes the bound tight:

\sum_z Q(z) \log \frac{p(x,z;\theta)}{Q(z)}

= \sum_z p(z\mid x;\theta) \log \frac{ p(x,z;\theta) }{ p(z\mid x;\theta) }

= \sum_z p(z\mid x;\theta) \log \frac{ p(z\mid x;\theta)p(x;\theta) }{ p(z\mid x;\theta) }

= \sum_z p(z\mid x;\theta) \log p(x;\theta)

= \log p(x;\theta) \sum_z p(z\mid x;\theta)

= \log p(x;\theta)

Formally, we define:

\text{ELBO}(x;Q,\theta) = \mathbb{E}_{z\sim Q} \left[ \log \frac{p(x,z;\theta)}{Q(z)} \right]

= \sum_z Q(z) \log \frac{p(x,z;\theta)}{Q(z)}

Therefore:

\log p(x;\theta) \geq \text{ELBO}(x;Q,\theta)

This gives the two steps of EM.

In the E-step, we choose $Q$ to make the bound tight at the current parameters. In the M-step, we keep $Q$ fixed and choose new parameters that maximize the bound.

\text{Repeat until a stopping criterion is met \{}

\quad \text{(E-step) For each } i, \text{ set \{}

\quad\quad Q_i^{(t)}(z^{(i)}) \leftarrow p(z^{(i)}\mid x^{(i)};\theta^{(t)})

\quad \text{\}}

\quad \text{(M-step) Set \{}

\quad\quad \theta^{(t+1)} \leftarrow \arg\max_\theta \sum_{i=1}^n \text{ELBO} \left( x^{(i)}; Q_i^{(t)}, \theta \right)

\quad \text{\}}

\text{\}}

In the E-step, we compute the posterior distribution of the latent variable for each training example using the current parameters:

p(z^{(i)}=j\mid x^{(i)};\theta^{(t)}) = \frac{ p(x^{(i)}\mid z^{(i)}=j;\theta^{(t)}) p(z^{(i)}=j;\theta^{(t)}) }{ \sum_{l=1}^k p(x^{(i)}\mid z^{(i)}=l;\theta^{(t)}) p(z^{(i)}=l;\theta^{(t)}) }

The E-step choice makes the bound tight at $\theta^{(t)}$ .

In the M-step, we hold $Q^{(t)}$ fixed and choose $\theta^{(t+1)}$ to maximize that bound. The bound does not need to remain tight at the new parameter value.

This is exact EM when the E-step computes the exact posterior and the M-step globally maximizes the fixed- $Q^{(t)}$ objective.

Generalized EM may instead use any M-step that increases that objective. With an exact E-step, either version makes the observed-data log-likelihood nondecreasing, although neither guarantees a global optimum.

Monotonicity and Convergence

Each exact EM iteration guarantees that the observed-data log-likelihood does not decrease:

\ell(\theta^{(t+1)}) \geq \ell(\theta^{(t)})

This does not mean that EM finds the global maximum. It may converge to a local optimum or another stationary point.

For any parameter value $\theta$ and any distribution over the latent variables, the dataset log-likelihood:

\ell(\theta) = \sum_i \log p(x^{(i)};\theta)

is at least the sum of the per-example ELBOs.

The ELBO is tangent to the log-likelihood at theta-t; maximizing it reaches theta-(t+1), where the log-likelihood is higher still — For fixed $Q^{t}$ , the $\text{ELBO}$ touches $\ell$ at $\theta^{t}$ (A); the $M$ -step maximizes the bound over $\theta$ to obtain $\theta^{t+1}$ (B), and since $\ell \geq \text{ELBO}$ , $\ell(\theta^{t+1})$ (C) is at least as high.

First:

\ell(\theta^{(t+1)}) \geq \sum_{i=1}^n \text{ELBO} \left( x^{(i)}; Q_i^{(t)}, \theta^{(t+1)} \right)

In the M-step, $\theta^{(t+1)}$ is chosen to maximize the fixed- $Q^{(t)}$ ELBO.

Therefore:

\sum_{i=1}^n \text{ELBO} \left( x^{(i)}; Q_i^{(t)}, \theta^{(t+1)} \right)

\geq \sum_{i=1}^n \text{ELBO} \left( x^{(i)}; Q_i^{(t)}, \theta^{(t)} \right)

Finally, the E-step chooses $Q$ so that the bound is tight at the current parameter value:

\sum_{i=1}^n \text{ELBO} \left( x^{(i)}; Q_i^{(t)}, \theta^{(t)} \right) = \ell(\theta^{(t)})

Putting these together:

\ell(\theta^{(t+1)}) \geq \ell(\theta^{(t)})

Tightness at $\theta^{(t)}$ is what turns an increase in the fixed- $Q^{(t)}$ bound into a guaranteed non-decrease in the observed-data likelihood.

Additional assumptions are needed for stronger convergence claims.

Successive tangent lower bounds climbing the log-likelihood toward a local maximum, with shrinking steps — In this illustrative run, each fresh bound touches $\ell$ at the current $\theta^k$ and is maximized to obtain the next iterate. The likelihood climbs toward a local maximum, although EM's general guarantee is only non-decrease.

For an unconstrained Gaussian mixture model, the likelihood is not even bounded above: one component can collapse its covariance around a data point.

In practice, covariance constraints, regularization, or priors are used to prevent this degeneracy.

Other Interpretations

The evidence lower bound can also be written as:

\text{ELBO}(x;Q,\theta) = \mathbb{E}_{z\sim Q} \left[ \log p(x,z;\theta) \right] + H(Q)

where:

H(Q) = - \mathbb{E}_{z\sim Q} [\log Q(z)]

Factoring the joint distribution gives:

\text{ELBO}(x;Q,\theta) = \mathbb{E}_{z\sim Q} \left[ \log p(x\mid z;\theta) \right] - D_{KL} \left( Q\parallel p(z;\theta) \right)

If the prior $p(z;\theta)$ is fixed and has no parameters in $\theta$ , then the KL term is constant with respect to $\theta$ during the M-step.

This is not true in general. For example, the GMM mixing weights parameterize the prior over component assignments.

See derivation

For two distributions over a random variable $z$ :

D_{KL}(Q\parallel p) = \sum_z Q(z) \log \frac{Q(z)}{p(z)}

= \mathbb{E}_{z\sim Q} \left[ \log \frac{Q(z)}{p(z)} \right]

\text{ELBO}(x;Q,\theta) = \mathbb{E}_{z\sim Q} \left[ \log \frac{p(x,z;\theta)}{Q(z)} \right]

= \mathbb{E}_{z\sim Q} \left[ \log \frac{ p(x\mid z;\theta)p(z;\theta) }{ Q(z) } \right]

= \mathbb{E}_{z\sim Q} \left[ \log p(x\mid z;\theta) + \log p(z;\theta) - \log Q(z) \right]

= \mathbb{E}_{z\sim Q} \left[ \log p(x\mid z;\theta) - \log \frac{Q(z)}{p(z;\theta)} \right]

= \mathbb{E}_{z\sim Q} \left[ \log p(x\mid z;\theta) \right] - \mathbb{E}_{z\sim Q} \left[ \log \frac{Q(z)}{p(z;\theta)} \right]

= \mathbb{E}_{z\sim Q} \left[ \log p(x\mid z;\theta) \right] - D_{KL} \left( Q\parallel p(z;\theta) \right)

A second useful form relates the ELBO directly to the log-likelihood:

\text{ELBO}(x;Q,\theta) = \log p(x;\theta) - D_{KL} \left( Q \parallel p(z\mid x;\theta) \right)

This shows that for a fixed value of $\theta$ , the ELBO is maximized when:

Q(z) = p(z\mid x;\theta)

because this makes the KL divergence equal to $0$ .

This is exactly the choice made in the E-step.

Log-likelihood split into the ELBO plus the KL gap; setting Q to the posterior closes the gap so the bound is tight — At a fixed $\theta$ , $\log p(x;\theta)$ splits into the $\text{ELBO}$ and the gap $D_{KL}(Q\parallel p(z\mid x;\theta))$ . The E-step sets $Q=p(z\mid x;\theta)$ , closing the gap.

See derivation

\text{ELBO}(x;Q,\theta) = \mathbb{E}_{z\sim Q} \left[ \log \frac{p(x,z;\theta)}{Q(z)} \right]

= \mathbb{E}_{z\sim Q} \left[ \log \frac{ p(z\mid x;\theta)p(x;\theta) }{ Q(z) } \right]

= \mathbb{E}_{z\sim Q} \left[ \log p(z\mid x;\theta) + \log p(x;\theta) - \log Q(z) \right]

= \mathbb{E}_{z\sim Q} \left[ \log p(x;\theta) - \log \frac{Q(z)}{p(z\mid x;\theta)} \right]

= \mathbb{E}_{z\sim Q} \left[ \log p(x;\theta) \right] - \mathbb{E}_{z\sim Q} \left[ \log \frac{Q(z)}{p(z\mid x;\theta)} \right]

= \log p(x;\theta) - D_{KL} \left( Q\parallel p(z\mid x;\theta) \right)

Gaussian Mixture Models

A Gaussian Mixture Model assumes that each observation comes from one of $k$ Gaussian components.

For each example $x^{(i)}$ , there is a latent variable:

z^{(i)} \in \{1,\ldots,k\}

indicating which Gaussian generated it.

First, a component is chosen according to:

p(z^{(i)}=j) = \phi_j

where:

\phi_j>0, \qquad \sum_{j=1}^k \phi_j = 1

Then the observation is drawn from that component:

x^{(i)} \mid z^{(i)}=j \sim \mathcal N(\mu_j,\Sigma_j)

where:

\mu_j\in\mathbb R^d

and:

\Sigma_j \in \mathbb R^{d\times d}

is a symmetric positive-definite covariance matrix.

Two panels over the same data: on the left only the unlabeled points are observed; on the right each point is colored by the latent z that chose one of two Gaussians, each drawn as covariance ellipses — All we observe is the unlabeled cloud $x^{(i)}$ (left). The model explains it with a hidden $z^{(i)}$ that chooses one of the Gaussians, which then generates the observation (right). If we knew the $z^{(i)}$ 's, fitting each Gaussian would be straightforward.

The component labels themselves are arbitrary: permuting all component labels gives the same mixture distribution.

In practice, covariance regularization is also commonly used to prevent a component from collapsing around a single training example.

The log-likelihood is:

\ell(\phi,\mu,\Sigma) = \sum_{i=1}^n \log p(x^{(i)};\phi,\mu,\Sigma)

= \sum_{i=1}^n \log \sum_{z^{(i)}=1}^k p(x^{(i)}\mid z^{(i)};\mu,\Sigma) p(z^{(i)};\phi)

The difficulty is that we do not know which component generated each observation.

If the $z^{(i)}$ 's were known, we could simply group the observations by component and estimate each Gaussian directly, similar to Gaussian Discriminant Analysis.

EM handles the missing assignments by alternating between estimating them and updating the Gaussian parameters.

Responsibilities in the E-Step

In the E-step, instead of assigning each point to exactly one component, EM gives every component a probability of having generated that point.

This probability is called a responsibility:

w_j^{(i)} = p \left( z^{(i)}=j \mid x^{(i)}; \phi,\mu,\Sigma \right)

Using Bayes' rule:

w_j^{(i)} = \frac{ p(x^{(i)}\mid z^{(i)}=j;\mu,\Sigma) p(z^{(i)}=j;\phi) }{ \sum_{l=1}^k p(x^{(i)}\mid z^{(i)}=l;\mu,\Sigma) p(z^{(i)}=l;\phi) }

The responsibilities for one example sum to $1$ :

\sum_{j=1}^k w_j^{(i)} = 1

After substituting the Gaussian densities, the responsibilities can also be written in a softmax-like form.

Top: two weighted Gaussian components and the mixture density they sum to. Bottom: the two posterior responsibilities as smooth curves crossing at one half, with sampled points sitting on them — The mixture $p(x)$ is the sum of the weighted components $\phi_j p(x\mid z=j)$ (top). The E-step assigns each point softly: the responsibilities $w_j^{(i)}=p(z=j\mid x^{(i)})$ move between $0$ and $1$ and sum to one (bottom).

See derivation

p(z^{(i)}=j\mid x^{(i)};\phi,\mu,\Sigma)

= \frac{ p(x^{(i)}\mid z^{(i)}=j;\mu,\Sigma) p(z^{(i)}=j;\phi) }{ \sum_{l=1}^k p(x^{(i)}\mid z^{(i)}=l;\mu,\Sigma) p(z^{(i)}=l;\phi) }

Substituting the Gaussian densities:

= \frac{ \frac{1}{(2\pi)^{d/2}|\Sigma_j|^{1/2}} \exp \left( -\frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j) \right) \phi_j }{ \sum_{l=1}^k \left[ \frac{1}{(2\pi)^{d/2}|\Sigma_l|^{1/2}} \exp \left( -\frac12 (x^{(i)}-\mu_l)^T \Sigma_l^{-1} (x^{(i)}-\mu_l) \right) \phi_l \right] }

= \frac{ |\Sigma_j|^{-1/2} \exp \left( -\frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j) \right) \phi_j }{ \sum_{l=1}^k \left[ |\Sigma_l|^{-1/2} \exp \left( -\frac12 (x^{(i)}-\mu_l)^T \Sigma_l^{-1} (x^{(i)}-\mu_l) \right) \phi_l \right] }

Define the score:

\alpha_j(x^{(i)}) = \log\phi_j - \frac12 \log|\Sigma_j| - \frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j)

Then:

p(z^{(i)}=j\mid x^{(i)};\phi,\mu,\Sigma) = \frac{ \exp(\alpha_j(x^{(i)})) }{ \sum_{l=1}^k \exp(\alpha_l(x^{(i)})) }

To verify this:

\frac{ \exp(\alpha_j(x^{(i)})) }{ \sum_{l=1}^k \exp(\alpha_l(x^{(i)})) }

= \frac{ \exp \left( \log\phi_j - \frac12\log|\Sigma_j| - \frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j) \right) }{ \sum_{l=1}^k \exp \left( \log\phi_l - \frac12\log|\Sigma_l| - \frac12 (x^{(i)}-\mu_l)^T \Sigma_l^{-1} (x^{(i)}-\mu_l) \right) }

= \frac{ \phi_j \exp \left( -\frac12\log|\Sigma_j| - \frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j) \right) }{ \sum_{l=1}^k \phi_l \exp \left( -\frac12\log|\Sigma_l| - \frac12 (x^{(i)}-\mu_l)^T \Sigma_l^{-1} (x^{(i)}-\mu_l) \right) }

= \frac{ |\Sigma_j|^{-1/2} \exp \left( -\frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j) \right) \phi_j }{ \sum_{l=1}^k \left[ |\Sigma_l|^{-1/2} \exp \left( -\frac12 (x^{(i)}-\mu_l)^T \Sigma_l^{-1} (x^{(i)}-\mu_l) \right) \phi_l \right] }

Parameter Updates in the M-Step

In the M-step, the responsibilities $w_j^{(i)}$ are held fixed while we update the parameters of each component.

The resulting updates are responsibility-weighted versions of the usual sample estimates:

\phi_l = \frac1n \sum_{i=1}^n w_l^{(i)}

\mu_l = \frac{ \sum_{i=1}^n w_l^{(i)}x^{(i)} }{ \sum_{i=1}^n w_l^{(i)} }

\Sigma_l = \frac{ \sum_{i=1}^n w_l^{(i)} (x^{(i)}-\mu_l) (x^{(i)}-\mu_l)^T }{ \sum_{i=1}^n w_l^{(i)} }

So each component is fitted as if every observation belonged to it partially, with weight $w_l^{(i)}$ .

These updates come from maximizing the fixed-responsibility ELBO:

\sum_{i=1}^{n} \sum_{j=1}^{k} w_j^{(i)} \log \frac{ p(x^{(i)}\mid z^{(i)}=j;\mu,\Sigma) p(z^{(i)}=j;\phi) }{ w_j^{(i)} }

See derivation

Start with:

\sum_{i=1}^{n} \sum_{j=1}^{k} w_j^{(i)} \log \frac{ p(x^{(i)}\mid z^{(i)}=j;\mu,\Sigma) p(z^{(i)}=j;\phi) }{ w_j^{(i)} }

Substituting the Gaussian density:

= \sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \log \left( \frac{ \frac{1}{(2\pi)^{d/2}|\Sigma_j|^{1/2}} \exp \left( -\frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j) \right) \phi_j }{ w_j^{(i)} } \right)

= \sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \left[ \log \left( \frac{1}{(2\pi)^{d/2}|\Sigma_j|^{1/2}} \right) - \frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j) + \log\phi_j - \log w_j^{(i)} \right]

First, take the derivative with respect to $\mu_l$ :

\nabla_{\mu_l} \sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \left[ -\frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j) \right]

= \sum_{i=1}^n w_l^{(i)} \nabla_{\mu_l} \left[ -\frac12 (x^{(i)}-\mu_l)^T \Sigma_l^{-1} (x^{(i)}-\mu_l) \right]

Expanding the quadratic term:

= \sum_{i=1}^n w_l^{(i)} \nabla_{\mu_l} \left[ -\frac12 \left( x^{(i)T}\Sigma_l^{-1}x^{(i)} - \mu_l^T\Sigma_l^{-1}x^{(i)} - x^{(i)T}\Sigma_l^{-1}\mu_l + \mu_l^T\Sigma_l^{-1}\mu_l \right) \right]

= \sum_{i=1}^n w_l^{(i)} \left[ -\frac12 \left( -\Sigma_l^{-1}x^{(i)} - \Sigma_l^{-T}x^{(i)} + 2\Sigma_l^{-1}\mu_l \right) \right]

= \sum_{i=1}^n w_l^{(i)} \left[ -\frac12 \left( -\Sigma_l^{-1}x^{(i)} - \Sigma_l^{-1}x^{(i)} + 2\Sigma_l^{-1}\mu_l \right) \right]

Therefore:

= \sum_{i=1}^n w_l^{(i)} \left( \Sigma_l^{-1}x^{(i)} - \Sigma_l^{-1}\mu_l \right)

Setting this equal to $0$ :

\sum_{i=1}^n w_l^{(i)} \Sigma_l^{-1}x^{(i)} = \sum_{i=1}^n w_l^{(i)} \Sigma_l^{-1}\mu_l

\Sigma_l^{-1} \left[ \sum_{i=1}^n w_l^{(i)}x^{(i)} \right] = \Sigma_l^{-1} \left[ \sum_{i=1}^n w_l^{(i)} \right] \mu_l

Canceling $\Sigma_l^{-1}$ :

\sum_{i=1}^n w_l^{(i)}x^{(i)} = \left[ \sum_{i=1}^n w_l^{(i)} \right] \mu_l

Therefore:

\mu_l = \frac{ \sum_{i=1}^n w_l^{(i)}x^{(i)} }{ \sum_{i=1}^n w_l^{(i)} }

Now take the derivative with respect to $\phi_l$ .

Only the terms:

\sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \log\phi_j

depend on $\phi$ .

Because:

\sum_{j=1}^k \phi_j = 1

we use a Lagrange multiplier:

\mathcal L(\phi) = \sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \log\phi_j + \beta \left( \sum_{j=1}^k \phi_j - 1 \right)

Taking the derivative:

\frac{\partial}{\partial\phi_l} \mathcal L(\phi) = \sum_{i=1}^n \frac{w_l^{(i)}}{\phi_l} + \beta

Setting this equal to $0$ :

\sum_{i=1}^n \frac{w_l^{(i)}}{\phi_l} + \beta = 0

\phi_l = \frac{ \sum_{i=1}^n w_l^{(i)} }{ -\beta }

Using the constraint:

\sum_{j=1}^k \phi_j = 1

gives:

\sum_{j=1}^k \frac{ \sum_{i=1}^n w_j^{(i)} }{ -\beta } = 1

\sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} = -\beta

\sum_{i=1}^n 1 = -\beta

Therefore:

-\beta=n

and:

\phi_l = \frac1n \sum_{i=1}^n w_l^{(i)}

We use the following matrix derivatives:

\nabla_A|A| = |A|(A^{-1})^T

and:

\nabla_A \left( x^TA^{-1}x \right) = -A^{-T}xx^TA^{-T}

When $A$ is symmetric:

A^{-T}=A^{-1}

Finally, take the derivative with respect to $\Sigma_l$ :

\nabla_{\Sigma_l} \sum_{i=1}^n \sum_{j=1}^k w_j^{(i)} \left[ \log \frac{1}{(2\pi)^{d/2}|\Sigma_j|^{1/2}} - \frac12 (x^{(i)}-\mu_j)^T \Sigma_j^{-1} (x^{(i)}-\mu_j) \right]

= \sum_{i=1}^n w_l^{(i)} \nabla_{\Sigma_l} \left[ -\log|\Sigma_l|^{1/2} - \frac12 (x^{(i)}-\mu_l)^T \Sigma_l^{-1} (x^{(i)}-\mu_l) \right]

We evaluate the two derivatives separately.

First:

\nabla_{\Sigma_l} \log|\Sigma_l|^{1/2}

= \frac{1}{|\Sigma_l|^{1/2}} \nabla_{\Sigma_l} |\Sigma_l|^{1/2}

= \frac{1}{|\Sigma_l|^{1/2}} \left( \frac12 |\Sigma_l|^{-1/2} \right) \nabla_{\Sigma_l} |\Sigma_l|

= \frac12 \left( \frac{1}{|\Sigma_l|} |\Sigma_l| \right) \Sigma_l^{-1}

= \frac12 \Sigma_l^{-1}

For the quadratic term:

\nabla_{\Sigma_l} \left[ (x^{(i)}-\mu_l)^T \Sigma_l^{-1} (x^{(i)}-\mu_l) \right]

= - \Sigma_l^{-1} (x^{(i)}-\mu_l) (x^{(i)}-\mu_l)^T \Sigma_l^{-1}

Putting these together:

\sum_{i=1}^n w_l^{(i)} \left[ \frac12 \left( -\Sigma_l^{-1} + \Sigma_l^{-1} (x^{(i)}-\mu_l) (x^{(i)}-\mu_l)^T \Sigma_l^{-1} \right) \right]

Setting this equal to $0$ gives:

\sum_{i=1}^n w_l^{(i)} \Sigma_l^{-1} = \sum_{i=1}^n w_l^{(i)} \Sigma_l^{-1} (x^{(i)}-\mu_l) (x^{(i)}-\mu_l)^T \Sigma_l^{-1}

\left( \sum_{i=1}^n w_l^{(i)} \right) \Sigma_l^{-1} = \Sigma_l^{-1} \left[ \sum_{i=1}^n w_l^{(i)} (x^{(i)}-\mu_l) (x^{(i)}-\mu_l)^T \right] \Sigma_l^{-1}

Multiplying through by $\Sigma_l$ gives:

\Sigma_l \left( \sum_{i=1}^n w_l^{(i)} \right) = \sum_{i=1}^n w_l^{(i)} (x^{(i)}-\mu_l) (x^{(i)}-\mu_l)^T

Therefore:

\Sigma_l = \frac{ \sum_{i=1}^n w_l^{(i)} (x^{(i)}-\mu_l) (x^{(i)}-\mu_l)^T }{ \sum_{i=1}^n w_l^{(i)} }

\text{Repeat until convergence \{}

\quad \text{(E-step) For each }i\text{ and }j,\text{ set \{}

\quad\quad w_j^{(i)} \leftarrow p \left( z^{(i)}=j \mid x^{(i)}; \phi,\mu,\Sigma \right)

\quad \text{\}}

\quad \text{(M-step) For each }j,\text{ set \{}

\quad\quad \phi_j \leftarrow \frac1n \sum_{i=1}^n w_j^{(i)}

\quad\quad \mu_j \leftarrow \frac{ \sum_{i=1}^n w_j^{(i)}x^{(i)} }{ \sum_{i=1}^n w_j^{(i)} }

\quad\quad \Sigma_j \leftarrow \frac{ \sum_{i=1}^n w_j^{(i)} (x^{(i)}-\mu_j) (x^{(i)}-\mu_j)^T }{ \sum_{i=1}^n w_j^{(i)} }

\quad \text{\}}

\text{\}}

It is useful to define the effective number of examples assigned to component $j$ :

N_j = \sum_{i=1}^n w_j^{(i)}

Because the assignments are soft, $N_j$ can be fractional.

The M-step then sets:

\phi_j = \frac{N_j}{n}

and computes responsibility-weighted means and covariances.

The responsibilities stay fixed throughout the M-step.

The E-step makes the ELBO equal to $\log p(x;\theta^{(t)})$ at the old parameters. After updating to $\theta^{(t+1)}$ , the bound is generally no longer tight, although the observed-data log-likelihood cannot decrease under exact EM.