Gaussian mixtures and the EM algorithm

The unsupervised half of the course has so far covered representation learning (PCA and kernel PCA, Parts 3 to 5), clustering with Lloyd's algorithm (Parts 6 and 7) and estimation by maximum likelihood and Bayesian methods (Part 8). This part continues with estimation, for data that is a little more realistic. The examples are one-dimensional because that is easiest to draw; everything carries over to higher dimensions.

Take real numbers on a line that fall into three separate groups. Estimation means proposing a generative story for the data and fitting its parameters. The values are real numbers, not 0s and 1s, so a Bernoulli model is out. The obvious choice is a single Gaussian. Maximum likelihood puts its mean at the sample mean, \(\hat\mu_{\text{ML}}\) (Part 8); take the variance as known, say 1.

Now look at the fitted density. Its peak sits in the middle of the data, and the two outer groups lie where the Gaussian has almost no probability. The model says data there should be rare, yet we saw plenty of it.

Maximum likelihood treats the assumed model as true and finds its best member. If no Gaussian can produce three separated groups, the best Gaussian is still a poor explanation, and maximum likelihood cannot do better than the family it searches.

What we want is a density with three modes, one high-density region for each group, so that the observed points come from places the model considers likely. A Gaussian has exactly one peak, so no single Gaussian will do. Each dense region on its own does look like a Gaussian, though. That suggests the new model: a mixture of Gaussians, in this example a mixture of three. The data clearly has a clustering behaviour, and the model will become a probabilistic way to describe it.

A generative story is the probabilistic mechanism we imagine produced the data. The stories so far had one step: toss a coin, or sample from a Gaussian. That is too simple for clustered data, and the mixture story has two steps. For now assume the number of mixtures \(K\) is given. "Mixture" and "cluster" mean the same thing here; mixture is the usual word.

1. Generating one data point \(x_i\)
2. Pick a mixture. Roll a die with \(K\) faces; face \(k\) comes up with probability \(\pi_k\). Call the face \(z_i\), the mixture that point \(i\) comes from.
3. Draw from that mixture. Every mixture has its own Gaussian with its own mean and variance. Draw \(x_i\) from the Gaussian of mixture \(z_i\).
4. To build a data set, repeat both steps for \(i = 1, \dots, n\). Each point gets its own roll of the die.

Standard name: the \(\pi_k\) are the mixing proportions. The lecture's example: with three mixtures and \(\pi_1 = 0.5\), \(\pi_2 = 0.25\), \(\pi_3 = 0.25\), half the points come from the first mixture and a quarter from each of the other two.

Step 2 is pictured as doors. If the die shows 5, walk through the fifth door; behind it sit a mean \(\mu_5\) and a variance \(\sigma_5^2\), and you sample a point from that Gaussian. The indicator \(z_i\) plays the same role as the cluster indicator \(z_i\) of K-means in Part 6, except that here it is drawn at random.

The data \(x_1, \dots, x_n\) are observed. The model also contains quantities we never see: the faces \(z_1, \dots, z_n\). Unobserved quantities that the observations depend on are called latent variables, and a model built on them is a latent variable model. The Gaussian mixture model (GMM) is one of the most commonly used. Its output depends on parameters we want to estimate and also on these hidden choices.

Which parameters pin the model down completely?

In Part 8 we estimated a single number, the bias of a coin or the mean of a Gaussian. Now there is a whole host of parameters, and looking at the data gives no obvious guess for the \(\pi\)'s, \(\mu\)'s and \(\sigma^2\)'s. Maximum likelihood is a principled method that gives an answer for any well-specified probabilistic model, so it is the next thing to try.

The likelihood depends on the unknowns \(\mu_1, \dots, \mu_K\), \(\sigma_1^2, \dots, \sigma_K^2\), \(\pi_1, \dots, \pi_K\) and on the observed data. As in Part 8, the data are held fixed and we maximize over the parameters. Write \(\theta\) for all the parameters together.

The i.i.d. assumption still holds. Every point goes through the same two steps (identically distributed), and the outcome for one point does not change the probabilities for another (independent). So the likelihood is a product over points of one density, the mixture density \(f(x_i;\theta)\):

The big \(\prod\) is an ordinary product. It has nothing to do with the mixing proportions \(\pi_k\).

What is the density of one point? We never saw which face the die showed for \(x_i\), so we account for every possibility.

• "\(x_i\) came from mixture 1", "from mixture 2", …, "from mixture \(K\)" are mutually exclusive events, so their chances add up.
• Coming from mixture \(k\) needs two things to happen together: the die shows \(k\) (probability \(\pi_k\)), and mixture \(k\)'s Gaussian produces \(x_i\) (density \(\mathcal N(x_i;\mu_k,\sigma_k^2)\)). So these multiply.

The lecture's example has three mixtures with means \(-20\), \(0\) and \(15\), and an observed point \(x_1 = -15\). It is tempting to say \(x_1\) must have come from mixture 1, whose mean is closest. That does not follow. Mixture 2 could have been picked, in which case its Gaussian explains \(x_1\) with a smaller density. Mixture 3 gives an even smaller density, tiny but not zero, since a Gaussian is positive everywhere. None of the three can be dismissed.

Closeness to a mean is only half the story; the mixing proportions matter too. Suppose \(\pi_3 = 0.9\) and \(\pi_1 = \pi_2 = 0.05\). Then mixture 1 is picked only 5% of the time and mixture 3 90% of the time. The density of \(x_1\) weighs both effects: how often each mixture is chosen, and how well its Gaussian explains the point. That is exactly what each term \(\pi_k\,\mathcal N(x_1;\mu_k,\sigma_k^2)\) of the sum measures.

A product is harder to handle than a sum, so take logs as usual. The lecture labels the result ★:

In Part 8 the log did two useful jobs: it turned the product over points into a sum, and it cancelled the exponential in the Gaussian density. Here it only does the first. A sum now sits inside the log, and there is no nice rule for the log of a sum. (A product inside a log would split into a sum of logs.)

So how do we maximize ★ over all the \(\mu_k\), \(\sigma_k^2\) and \(\pi_k\)?

• No closed form. Setting the derivative with respect to \(\mu_k\) to zero does not give a formula like "\(\hat\mu_{\text{ML}}\) is the average of the data points". In general the equations cannot be solved analytically.
• A constraint. The \(\pi_k\) are not free variables: they must be nonnegative and sum to 1. Maximizing over them is a constrained optimization problem, which is harder still.
• General-purpose methods ignore the structure. One could write down the gradient and run gradient ascent, since this is a maximization (gradient methods appear in Part 10). It works for any function and would give some estimates. But it does not use what we know about this model: the data were generated in two well-defined steps, first a cluster indicator, then a Gaussian draw.

Goal: an efficient way to maximize ★ that exploits its two-step structure. The lecture first takes a short detour through convex functions.

Convexity could fill a course of its own, and convex functions appear all over machine learning. Only the bare essentials are needed here. Pick two points \(a\) and \(b\) and their midpoint \(\frac{a+b}{2}\). Compare two heights above the midpoint. One is the function itself, \(f\big(\frac{a+b}{2}\big)\). The other is the straight line joining \((a, f(a))\) to \((b, f(b))\), the linear interpolation between them, whose height there is \(\frac{f(a)+f(b)}{2}\), exactly halfway between \(f(a)\) and \(f(b)\).

The inequality must hold for every choice of \(a\) and \(b\), not just one pair you happened to draw. The names match the convex and concave mirrors of school physics.

• Both convex and concave. Both inequalities hold together only as an equality, \(f\big(\frac{a+b}{2}\big) = \frac{f(a)+f(b)}{2}\) for all \(a, b\). That is a straight line: interpolating linearly between two points of a line lands back on the line.
• Neither. A wavy function can have one pair of points whose chord lies above the curve, so it is not concave, and another pair whose chord lies below the curve, so it is not convex.

Write the midpoint as \(\frac12 a + \frac12 b\). For a convex function the definition reads \(f(\frac12 a + \frac12 b) \le \frac12 f(a) + \frac12 f(b)\). If this holds for every \(a\) and \(b\), one can argue that it holds for any weight \(\lambda\) between 0 and 1:

As \(\lambda\) goes from 1 to 0, the point \(\lambda a + (1-\lambda)b\) travels from \(a\) to \(b\), and at every stop along the way the function stays at or below the chord. For a concave function the inequality reverses.

The same idea extends to any number of points \(a_1, \dots, a_K\), with weights \(\lambda_1, \dots, \lambda_K\) that lie between 0 and 1 and sum to 1. Written for concave functions, the form used from here on:

In one dimension this adds little intuition, so the lecture pictures a concave function on the plane. Take three points \(a_1, a_2, a_3\). Varying the weights moves \(\lambda_1 a_1 + \lambda_2 a_2 + \lambda_3 a_3\) around the triangle they span, called their convex hull; for example \(0.1a_1 + 0.5a_2 + 0.4a_3\) is one point inside it. Over every point of the hull, the concave surface lies above the flat triangle joining the three function values. Jensen says this for any number of points.

Why convexity in a part about mixtures? The logarithm in the log-likelihood is a concave function, so Jensen's inequality applies to it. For positive \(a_k\) and valid weights:

Read it left to right. The left side has a sum inside a log, the very thing that blocked us in ★. The right side has the sum outside the log. Jensen lets us pull the sum out, at a price: we get an inequality instead of an equality. Lecture 25 shows how to live with that price.

The lecture leaves the concavity of the log as an exercise. One way, with the midpoint definition: for \(a, b > 0\),

and the last statement is always true. The first step is reversible because the log is increasing. (The middle inequality is the familiar fact that the arithmetic mean is at least the geometric mean.)

The video's YouTube title mentions supervised learning, but L25 is entirely about estimating the GMM's parameters.

Go back to ★. Inside each log is a plain sum of \(K\) numbers, but Jensen needs a weighted combination whose weights sum to 1. There are no such weights, so the lecture introduces them artificially, making the problem bigger on purpose. For every data point \(i\), add \(K\) new parameters \(\lambda_{i1}, \dots, \lambda_{iK}\) with

For each point they form a probability distribution over the \(K\) mixtures. Their meaning comes later; for now they are just extra parameters. Multiply and divide every term by \(\lambda_{ik}\), which changes nothing:

Each log now holds a \(\lambda\)-weighted combination of the numbers \(\pi_k\mathcal N(x_i;\mu_k,\sigma_k^2)/\lambda_{ik}\), and the log is concave. Jensen moves the sum outside:

Convention: a term with \(\lambda_{ik} = 0\) counts as 0, the limit of \(\lambda\log(c/\lambda)\) as \(\lambda \to 0\).

Why accept a different function with more parameters? L25 · 7:46

At first sight this is a bad trade. \(\ell\) is not the function we wanted to maximize, and it adds \(n \times K\) parameters. To be worth it, something must become very simple. Two things do.

• It is a lower bound for any valid λ. Hand me any \(\theta\). I can choose any \(\lambda\) I like (nonnegative, summing to 1 for each point), evaluate \(\ell(\theta,\lambda)\), and the number will never exceed the true log-likelihood \(\log L(\theta)\). That is Jensen, nothing more.
• Each half is easy when the other is fixed 10:10. With \(\lambda\) fixed, maximizing \(\ell\) over \(\theta\) is easy. With \(\theta\) fixed, maximizing over \(\lambda\) is easy. The next two sections check both claims by writing down the closed-form answers.

Key insight: one hard problem in θ becomes two easy problems, one in θ and one in λ.

Treat every \(\lambda_{ik}\) as a constant. The log of a product splits into a sum of logs, and the log cancels the exponential:

The last term involves only \(\lambda\), so it is a constant here and can be ignored. The first three terms contain the real parameters, and each appears in a simple way. Taking derivatives gives closed forms for all of them. The lecture states the results and recommends deriving them as an exercise; the steps are short.

• Means. \(\dfrac{\partial\ell}{\partial\mu_k} = \sum_i \lambda_{ik}\,\dfrac{x_i-\mu_k}{\sigma_k^2} = 0\) gives \(\mu_k = \sum_i\lambda_{ik}x_i \big/ \sum_i\lambda_{ik}\).
• Variances. Write \(v = \sigma_k^2\), so \(\log(\sqrt{2\pi}\,\sigma_k) = \frac12\log 2\pi + \frac12\log v\). Then \(\dfrac{\partial\ell}{\partial v} = \sum_i \lambda_{ik}\Big(\dfrac{(x_i-\mu_k)^2}{2v^2} - \dfrac{1}{2v}\Big) = 0\) gives \(v = \sum_i\lambda_{ik}(x_i-\mu_k)^2 \big/ \sum_i\lambda_{ik}\).
• Mixing proportions. Only \(\sum_i\sum_k\lambda_{ik}\log\pi_k\) involves \(\pi\), and the \(\pi_k\) must be nonnegative and sum to 1. The lecture solves this constrained problem with the method of Lagrange multipliers (constrained optimization comes back in Part 16). With \(N_k = \sum_i\lambda_{ik}\), set the derivative of \(\sum_k N_k\log\pi_k + \nu\big(1 - \sum_k\pi_k\big)\) to zero: \(N_k/\pi_k = \nu\), so \(\pi_k = N_k/\nu\). Summing over \(k\) fixes \(\nu = \sum_k N_k = \sum_i\sum_k\lambda_{ik} = n\).

Reading λ as a clustering L25 · 16:22

Think of \(\lambda_{ik}\) as the probability that point \(i\) belongs to cluster \(k\).

• Hard clustering. Put each point entirely into one cluster: \(\lambda_{ik} = 1\) for its cluster and 0 for the rest. Then \(\hat\mu_k\) is the mean of the points assigned to cluster \(k\), \(\hat\sigma_k^2\) their sample variance, and \(\hat\pi_k\) the fraction of points assigned to \(k\).
• Soft clustering. With \(\lambda_{ik}\) between 0 and 1, every point tells each cluster how much it belongs there, and is weighted by that importance. The estimates become weighted means and variances, and \(\hat\pi_k\) is the average chance that a point belongs to cluster \(k\).

Now hold all the \(\mu_k\), \(\sigma_k^2\) and \(\pi_k\) fixed. Write \(a_{ik} = \pi_k\,\mathcal N(x_i;\mu_k,\sigma_k^2)\), the weighted Gaussian density, which is a constant now. Splitting the log of the ratio leaves

No term mixes the \(\lambda\)'s of two different points. The function is a sum of separate pieces, one per data point, that share no variables, so maximizing each piece on its own maximizes the whole. Fix any point \(i\) and solve

This is a constrained problem again, and again not a hard one; Lagrange multipliers solve it analytically. The lecture skips the step, which goes like this: the derivative of \(\sum_k(\lambda_{ik}\log a_{ik} - \lambda_{ik}\log\lambda_{ik}) + \nu_i\big(1-\sum_k\lambda_{ik}\big)\) with respect to \(\lambda_{ik}\) is \(\log a_{ik} - \log\lambda_{ik} - 1 - \nu_i\). Setting it to zero makes \(\lambda_{ik}\) proportional to \(a_{ik}\), and the constraint fixes the scale.

It is Bayes' theorem L25 · 27:35

The formula looks complicated, but its meaning is plain. \(\lambda_{ik}\) was meant to be the chance that point \(i\) comes from cluster \(k\). If all the parameters are known and we are shown \(x_i\), that chance is \(P(z_i = k \mid x_i)\), and Bayes' theorem gives

Here \(P(x_i \mid z_i = k)\) is the density of \(x_i\) under mixture \(k\), \(P(z_i = k) = \pi_k\) is the chance from step 1 of the story, and the denominator adds up every way of generating \(x_i\), which is \(P(x_i)\). So maximizing the modified likelihood over \(\lambda\) gives exactly what Bayes' theorem gives. Standard name: \(\lambda_{ik}\) is often called the responsibility of mixture \(k\) for point \(i\).

Both halves are now simple formulas. Fix λ: weighted mean, weighted variance and weighted fraction. Fix θ: the posterior P(zᵢ = k | xᵢ) by Bayes' theorem.

Going back and forth between the two easy maximizations gives an iterative algorithm. The superscript \(t\) counts iterations.

1. EM for a Gaussian mixture
2. Initialize \(\theta^0\): some means \(\mu_1, \dots, \mu_K\), variances \(\sigma_1^2, \dots, \sigma_K^2\) and proportions \(\pi_1, \dots, \pi_K\).
3. E step. \(\lambda^{t+1} = \arg\max_\lambda \ell(\theta^t, \lambda)\): for every \(i\) and \(k\), \(\lambda_{ik}^{t+1} = \pi_k\mathcal N(x_i;\mu_k,\sigma_k^2) \big/ \sum_l \pi_l\mathcal N(x_i;\mu_l,\sigma_l^2)\), with the parameters taken from \(\theta^t\).
4. M step. \(\theta^{t+1} = \arg\max_\theta \ell(\theta, \lambda^{t+1})\): the weighted means, variances and proportions, with weights \(\lambda^{t+1}\).
5. Repeat the E and M steps until \(\lVert\theta^{t+1} - \theta^t\rVert \le \epsilon\), a tolerance such as \(10^{-3}\) or \(10^{-2}\): stop once the estimates barely change.

By introducing the \(\lambda\)'s and using Jensen, the problem splits into two parts, each easy when the other is held fixed, and we go back and forth between them. This is an instance of a famous algorithm, the EM algorithm. E stands for expectation and M for maximization. The first step is called the expectation step because the modified log-likelihood, and the \(\lambda^{t+1}\) it produces, can be written as an expectation of some quantity. (For instance, \(\sum_k\lambda_{ik}\log\frac{a_{ik}}{\lambda_{ik}}\) is an average of \(\log\frac{a_{ik}}{\lambda_{ik}}\) with weights \(\lambda_{i1},\dots,\lambda_{iK}\).) The course does not develop the general version.

The algorithm was developed in the 1970s. It was known in other forms before, but Dempster's paper set it down and gave it the name EM (standard reference: Dempster, Laird and Rubin, 1977), and it has been widely used ever since. The lecture also says one can argue that it converges, without giving a proof.

Why the log-likelihood never decreases

This standard argument is not in the lecture, but it shows what EM achieves and explains the widget below. It needs one more fact: after an E step, the bound is tight. Put \(\lambda_{ik} = a_{ik}/\sum_l a_{il}\) into \(\ell\). Then \(a_{ik}/\lambda_{ik} = \sum_l a_{il} = f(x_i;\theta)\), the mixture density of \(x_i\), for every \(k\). So

since the \(\lambda_{ik}^{t+1}\) sum to 1 over \(k\). Now chain the three facts:

The M step can only raise \(\ell\), because \(\theta^{t+1}\) maximizes it and \(\theta^t\) was one of the candidates. So every iteration keeps or raises the true log-likelihood.

EM is closely related to Lloyd's algorithm (Part 6). The \(\lambda_{ik}\) it ends with can be read as the chance that point \(i\) belongs to cluster \(k\), so EM produces a soft clustering where Lloyd's produces a hard one. At the end you have the estimates of \(\mu\), \(\sigma^2\) and \(\pi\), and also the optimal \(\lambda\)'s, which you can use as a clustering. Lloyd's algorithm also alternated two steps, and they line up with EM's:

Estimating variances matters. In higher dimensions they become covariance matrices, so each cluster gets its own shape. If the points of one cluster spread along one direction and those of a neighbouring cluster along another, Lloyd's algorithm may do poorly where the two overlap. EM can pick up that each cluster has its own direction of spread, simply because it estimates more parameters.

To turn EM's output into a hard clustering, run the algorithm, then put each point \(i\) in the cluster with the largest of \(\lambda_{i1}, \lambda_{i2}, \dots, \lambda_{iK}\). Because the variances shape these clusters, they need not be Voronoi regions 8:55. Especially where clusters overlap, the lecture argues EM assigns points better than Lloyd's, even after this conversion to hard clusters.

How does EM relate to the problem we started with, maximizing \(\log L(\theta)\)? We wrote down the log-likelihood and then worked entirely with the modified one instead. Picture the parameter space as a horizontal axis with \(\log L\) plotted above it. For simple models such as a single Gaussian or a Bernoulli, this curve is a nice concave function and easy to maximize. For a Gaussian mixture it is a complicated landscape, and we want the \(\theta^*\) at its highest point.

EM never climbs \(\log L\) directly. At \(\theta^t\), the E step picks \(\lambda^{t+1}\). As a function of \(\theta\), the bound \(\ell(\theta, \lambda^{t+1})\) is a nice function that, by Jensen, stays below \(\log L\) everywhere. The M step moves to its top, \(\theta^{t+1}\). The next E step builds a new bound there, its top gives \(\theta^{t+2}\), and so on. The guesses keep getting better, but they can settle on a peak that is not the highest one.

Where EM ends up therefore depends on where it starts. A very bad initialization may only reach a peak much lower than the best one. The "poor start" in the EM widget shows this: it converges, but to a log-likelihood about 8 below what the Lloyd's start reaches on the same data.

Since initialization decides which maximum EM reaches, the usual practice is to start EM from a K-means solution.

1. Building θ⁰ from Lloyd's
2. Run Lloyd's algorithm with the same \(K\) to get a hard clustering (Part 6; Part 7 covers its own initialization).
3. For each cluster, compute the sample mean and the sample variance of its points: these are \(\mu_k\) and \(\sigma_k^2\) in \(\theta^0\).
4. Set \(\pi_k\) to the fraction of data points in cluster \(k\).
5. Run EM from this \(\theta^0\).

These are the M-step formulas of Lecture 25 applied to hard (0 or 1) weights, so \(\theta^0\) is what one M step would produce from the K-means assignment.

The thread of these five lectures: data with cluster structure called for a latent variable model whose latent variable is a cluster indicator. Maximizing its log-likelihood was hard. Jensen's inequality introduced new variables so that the problem could be solved by alternate maximization: fix one set of variables and maximize over the other, then swap. Initialize with Lloyd's, run EM, and use the estimates however the later task needs, for instance converting the \(\lambda\)'s into a clustering.

For estimation this is a very powerful technique. In practice EM converges fast and typically gives very good parameter estimates. The principle reaches far beyond Gaussian mixtures: it applies to any reasonable latent variable model whose log-likelihood has a similar structure that can be separated. The course writes down EM only for the Gaussian mixture.

This ends unsupervised learning in the course: representation learning (PCA and kernel PCA), clustering (K-means) and estimation (maximum likelihood, Bayesian estimation, and mixture models fitted by EM). Part 10 begins a different paradigm, supervised learning.

Gaussian mixture model (GMM)

A density \(\sum_k \pi_k\mathcal N(x;\mu_k,\sigma_k^2)\): pick a component with probability \(\pi_k\), then draw from its Gaussian.

Mixing proportions \(\pi_k\)

The probabilities of the \(K\) faces of the die; nonnegative and summing to 1.

Generative story

The probabilistic mechanism assumed to have produced the data; for a GMM it has two steps.

Latent variable

An unobserved quantity the data depend on; here the component label \(z_i\) of each point.

Multimodal density

A density with several peaks. A single Gaussian has only one, so it cannot fit clustered data.

Convex function

\(f(\frac{a+b}{2}) \le \frac{f(a)+f(b)}{2}\) for all \(a, b\): the curve lies below every chord.

Concave function

The reverse inequality: the curve lies above every chord. The logarithm is concave.

Jensen's inequality

For concave \(f\) and weights \(\lambda_k \ge 0\) summing to 1: \(f(\sum_k\lambda_k a_k) \ge \sum_k\lambda_k f(a_k)\).

Modified log-likelihood \(\ell(\theta,\lambda)\)

Jensen's lower bound on \(\log L(\theta)\), valid for every valid \(\lambda\) and easy to maximize in each argument.

Responsibility \(\lambda_{ik}\)

The weight of point \(i\) on mixture \(k\); after an E step it equals \(P(z_i = k \mid x_i)\).

Soft clustering

Each point gets a probability for every cluster instead of a single label; take the largest to get a hard clustering.

EM algorithm

Alternate the E step (best \(\lambda\) for the current \(\theta\)) and the M step (best \(\theta\) for that \(\lambda\)) until \(\theta\) stops changing.

Local maximum

A peak of \(\log L\) that need not be the highest; EM converges to one, chosen by its initialization.