Gaussian Mixture Models
1 Introduction
In this section, we introduce Gaussian Mixture Models, or GMMs. Like K-means, fitting a GMM is an unsupervised method that can be used to find clusters in unlabelled data. Unlike K-means, which did not explicitly have anything to say about the likelihood of data or its method of generation, a GMM is a generative model that assumes data was generated from a collection of \(K\) Gaussians in the data space.
2 Understanding GMMs
A GMM represents a distribution as a mixture of Gaussians \[p({\bf x}) = \sum_{k=1}^K \pi_k \, \mathcal{N}({\bf x} | \boldsymbol\mu_k,\boldsymbol\Sigma_k)\] where:
- \(\textbf{x}\) is a vector representing one data point,
- \(K\) is the number of Gaussians (this is a hyperparameter, just like in K-means),
- \(\boldsymbol\mu_k\) and \(\boldsymbol\Sigma_k\) are the mean and covariance of the \(k^\text{th}\) Gaussian (introduced here),
- and \(\pi_k\) are the mixing coefficients.
For the mixing coefficients, \(\pi_k\),
\[\sum_{k=1}^K \pi_k = 1\quad \text{ and }\quad \pi_k\geq 0 \quad \forall k\]
We can interpret these mixing coefficients as the prior probability that any data point comes from a particular Gaussian. Or in other words, if we were to generate a random data point using the GMM, \(\pi_k\) gives the probability that it would be generated from Gaussian \(k\).
The function \(p({\bf x})\) defines a probability density over \({\bf x}\) (i.e. gives us the likelihood of our data point). This means that we can use maximum likelihood to find all the parameters for our GMM above that match the data density of \({\bf x}\) as closely as possible.
For example, we can take the leaf datasets from the previous section on K-means and instead fit GMMs to them. These are given in the figures below, where the individual Gaussians are indicated by their centers (outlined in black) and coloured contours at 0.5 and 1 standard deviations. (All of the contour lines you see in this section will be at 0.5 and 1 standard deviations.)
Leaf Data - 3 Cluster GMM
A Troublesome Dataset for K-means, but not for a GMM
Take a moment to compare these to their K-means equivalents (here, and here respectively). Notice that a GMM is better able to fit non-circular clusters, particularly in the second figure. In fact, GMMs are universal approximators of density, which is to say that a GMM can approximate any data generating function to an arbitrary level of precision.
2.1 Generating Data from a GMM
Before we move on to learning parameter values, let’s discuss how we can use a GMM that has already been fit to generate data. This can be a useful thing to do in data augmentation (i.e. generating synthetic data for training other models), and is also useful for developing intuition about how a GMM works.
Our task is to take a trained GMM and generate a completely new set of data points that looks like the original dataset the GMM was trained on (i.e. it should appear as if our new dataset was generated by the same underlying process as the original). We can generate each new data point via a two-step process:
Choose a “cluster” (i.e., a Gaussian) from which to sample \({\bf x}^{(i)}\). \[z^{(i)} \sim \textrm{Categorical}(\boldsymbol\pi)\]
Sample the observable features \({\bf x}^{(i)}\) from the chosen Gaussian. \[{\bf x}^{(i)} \sim \mathcal{N}(\boldsymbol\mu_{z^{(i)}}, \boldsymbol\Sigma_{z^{(i)}})\]
In other words, first choose one of the individual Gaussians with probability given by the mixing coefficients, and then use the parameters of that Gaussian to place the point in the data space. So, the point will be placed with highest probability near the Gaussian’s center, and diminishing probability as distance from the center increases.
In the figure below, we have copied the GMM from Figure 1 above. Try generating some new data, and observe how it has the same shape as the data in the original figure, but is not actually made up of the same individual data points.
Generating Data from a GMM
2.2 Fitting GMMs (a First Attempt)
In the figure above, you may have noticed that if you generated enough points, there was some intermingling between the clusters. Because we were implementing the data generation process, we knew which Gaussian each point actually came from, but if you removed the coloured labels, would you be able to tell? Can you say which Gaussian each of the following points was generated from?
3 Clusters - Clustering A
It certainly looks as if Point A could have been generated by the brown circle Gaussian, given how close it is to that Gaussian’s center, and how far it is from the others. On the other hand, Point B looks as if it were about equally likely to have been generated by the brown circle or orange square Gaussians. If we were fitting a GMM and these points were part of the unlabelled dataset, we could never know to which Gaussian they truly “belonged”. (And, strictly speaking, real data is generated by some natural, probably non-Gaussian process that we only approximate with Gaussians.) So, we say that the cluster assignments, \(z^{(i)}\), are latent variables, meaning that we cannot observe them.
If we knew the cluster assignment of each data point, then finding the parameters of our Gaussians would reduce to just Gaussian Discriminant Analysis, as discussed in the last chapter. However, because the cluster assignments are latent, we are forced to consider the contribution of each Gaussian to the probability of every data point. This is why, when we introduced the likelihood of a single data point, we summed over all Gaussians:
\[p({\bf x}) = \sum_{k=1}^K \pi_k \, \mathcal{N}({\bf x} | \boldsymbol\mu_k,\boldsymbol\Sigma_k)\]
Then, the likelihood of our full dataset, \({\bf X}= \{{\bf x}^{(i)}\}_{i=1}^{N}\), is given by multiplying together the likelihood of each individual point:
\[p({\bf X}) = \prod_{i=1}^{N}\sum_{k=1}^K \pi_k \, \mathcal{N}({\bf x^{(i)}} | \boldsymbol\mu_k,\boldsymbol\Sigma_k)\]
We can then use our trick of taking the log-likelihood to simplify our expression for maximum likelihood:
\[\log p({\bf X}) = \sum_{i=1}^N \log \left( \sum_{k=1}^K \pi_k \, \mathcal{N}({\bf x}^{(i)} ; \boldsymbol\mu_k,\boldsymbol\Sigma_k)\right)\]
…and we are stuck. When we take partial derivatives of this function, instead of getting neatly separated equations for each parameter, we find that the parameters are all depend on each other. In general, there are no closed-form solutions to this problem.
While gradient descent is technically possible, it has a number of issues that we state here without proof:
Non-convexity (due to permutation symmetry).
We need to enforce non-negativity constraint on \(\pi_k\).
We need the covariance matrices \(\boldsymbol\Sigma_k\) to be positive semi-definite.
Derivatives with respect to \(\boldsymbol\Sigma_k\) are expensive/complicated.
Fortunately, as with K-means, there is an algorithm we can use that alternates back and forth between two simpler problems.
3 Expectation-Maximization
We have already observed that if we somehow knew the latent cluster assignments, then the parameters \(\pi_k\), \(\boldsymbol\mu_k\) and \(\boldsymbol\Sigma_k\) (which we collectively refer to as \(\theta\)) could be fit using GDA. Conversely, if we knew the values of the parameters, \(\theta\), it would be straightforward to find the most likely cluster assignments with Bayes’ Rule, just as we did for inference after GDA:
\[P(z=k | {\bf x}) = \frac{P(z=k) \, p({\bf x} | z=k)}{\sum_\ell P(z=\ell) \, p({\bf x} | z=\ell)}\]
If we wanted a hard cluster assignment, we could assign each point to the cluster for which \(P(z=k | {\bf x})\) was the greatest. However, since we already have the probabilities above that neatly coincide with the expected contribution of each Gaussian, it is much more common for GMMs to use soft responsibilities rather than hard cluster assignments, just as we did for Soft K-means.
We now find ourselves in a similar situation to the one we had for K-means, and a similar solution presents itself. If we fix \(\theta\), we can calculate a set of responsibilities for each data point. If we then fix those responsibilities, we can recalculate more accurate values for \(\theta\). Alternating between these steps, we gradually improve our clustering via block coordinate descent. This is the Expectation-Maximization (EM) algorithm.
3.1 The Expectation Step
So named because we are taking the expectation of the latent cluster assignments, the Expectation Step fixes the parameter values \(\theta\) and produces new responsibilities for each data point, \(r_k^{(i)}\):
\[r_k^{(i)}= P(z^{(i)}=k | {\bf x}^{(i)}; \theta) = \frac{P(z=k) \, p({\bf x}^{(i)} | z=k)}{\sum_\ell P(z=\ell) \, p({\bf x}^{(i)} | z=\ell)} = \frac{\pi_k \mathcal{N}({\bf x}^{(i)} ; \boldsymbol\mu_k,\boldsymbol\Sigma_k)}{\sum_\ell \pi_\ell \mathcal{N}({\bf x}^{(i)} ; \boldsymbol\mu_\ell,\boldsymbol\Sigma_\ell)}\]
Question: For a practice problem that is easy to do by hand, let’s consider some one-dimensional data. Suppose we have the following measurements for the temperature in Toronto this month: \[x^{(1)} = 10,\ x^{(2)} = 12,\ x^{(3)} = 5,\ x^{(4)} = 4,\ x^{(5)} = 9\] We wish to fit a mixture of \(K=2\) Gaussians to this dataset. Suppose we initialize \(\mu_1 = 2\), \(\mu_2 = 7\), \(\sigma_1 = \sigma_2 = 1\), and \(\pi_1 = \pi_2 = 0.5\). What is the value of the responsibility vector \({\bf r}^{(1)}\) for the first example? Recall that the probability density function of the univariate Gaussian is \[p(x) = \mathcal{N}(x ; \mu, \sigma) = \frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{1}{2} \left( \frac{x - \mu}{\sigma} \right)^2}\]
Answer:
To get started, use Bayes’ Rule: \[\begin{align*} r_k^{(1)} &= P(z^{(1)}=k | x^{(1)}; \theta) &= \frac{p(x^{(1)} | z^{(1)}=k)\, \pi_k}{\sum_{\ell=1}^2 p(x^{(1)} | z^{(1)}=\ell)\, \pi_\ell} \end{align*}\]
\[\begin{align*} p(x^{(1)}=10 | z^{(1)}=1)\, \pi_1 &= \mathcal{N}(10 ; \mu=2, \sigma=1)\, \pi_1 &\approx 5.05 \times 10^{-15} \\ p(x^{(1)}=10 | z^{(1)}=2)\, \pi_2 &= \mathcal{N}(10 ; \mu=7, \sigma=1)\, \pi_2 &\approx 0.00443 \end{align*}\]
So, since the second term is much, much larger than the first, \(r_1^{(1)} \approx 0\) and \(r_2^{(1)} \approx 1\), so \({\bf r}^{(1)} \approx \begin{bmatrix}0 \\ 1\end{bmatrix}\).
3.2 The Maximization Step
In the Maximization Step, we fix the responsibilities for each data point and optimize the parameters \(\theta\) to maximize the data likelihood. To do this, we carry over the equations for GDA found by maximum likelihood in the last chapter, and modify them slightly to account for soft responsibilities. (Instead of each point counting as “1 point” in exactly one class, each point counts as a fraction of a point for each cluster, given by that cluster’s responsibility for it.)
\[\begin{eqnarray*} \pi_k &=& \frac{1}{N} \sum_{i=1}^N r_k^{(i)} \\ \boldsymbol\mu_k &=& \frac{\sum_{i=1}^N r_k^{(i)}\cdot {\bf x}^{(i)}}{\sum_{i=1}^N r_k^{(i)}} \\ \boldsymbol\Sigma_k &=& \frac{1}{\sum_{i=1}^N r_k^{(i)}} \sum_{i=1}^N r_k^{(i)} ({\bf x}^{(i)} -\boldsymbol\mu_k)({\bf x}^{(i)} -\boldsymbol\mu_k)^\top \end{eqnarray*}\]
Question: Continuing the temperature example above, suppose that after an E-step we have \[\begin{align*} {\bf r}^{(1)} \approx \begin{bmatrix}0 \\ 1\end{bmatrix}, \quad {\bf r}^{(2)} \approx \begin{bmatrix}0 \\ 1\end{bmatrix}, \quad {\bf r}^{(3)} \approx \begin{bmatrix}0.08 \\ 0.92\end{bmatrix}, \quad {\bf r}^{(4)} \approx \begin{bmatrix}0.92 \\ 0.08\end{bmatrix}, \quad {\bf r}^{(5)} \approx \begin{bmatrix}0 \\ 1\end{bmatrix} \end{align*}\] Re-estimate \(\mu_1\), \(\mu_2\), \(\sigma_1\), \(\sigma_2\), \(\pi_1\), and \(\pi_2\).
Answer:
Apply the M-step updates with \(N=5\): \[\begin{align*} \pi_1 &= \frac{1}{5} \sum_{i=1}^5 r_1^{(i)} = \frac{1}{5}(0 + 0 + 0.08 + 0.92 + 0) = 0.2 \\ \pi_2 &= \frac{1}{5} \sum_{i=1}^5 r_2^{(i)} = 1 - \pi_1 = 0.8 \\ \mu_1 &= \frac{\sum_{i=1}^5 r_1^{(i)}\cdot x^{(i)}}{\sum_{i=1}^5 r_1^{(i)}} = \frac{0.08 \cdot 5 + 0.92 \cdot 4}{0.08 + 0.92} = 4.08 \\ \mu_2 &= \frac{\sum_{i=1}^5 r_2^{(i)}\cdot x^{(i)}}{\sum_{i=1}^5 r_2^{(i)}} = \frac{10 + 12 + 0.92 \cdot 5 + 0.08 \cdot 4 + 9}{4} = 8.98 \\ \sigma_1^2 &= \frac{1}{\sum_{i=1}^5 r_1^{(i)}} \sum_{i=1}^5 r_1^{(i)} (x^{(i)} - \mu_1)^2 \approx 0.074 \implies \sigma_1 \approx 0.27 \\ \sigma_2^2 &= \frac{1}{\sum_{i=1}^5 r_2^{(i)}} \sum_{i=1}^5 r_2^{(i)} (x^{(i)} - \mu_2)^2 \approx 6.68 \implies \sigma_2 \approx 2.58 \end{align*}\]
3.3 Expectation-Maximization - Putting It All Together
The Expectation-Maximization (EM) algorithm is as follows:
Initialize the values of \(\theta\). This may be done randomly or, for example, by spreading the initial clusters evenly throughout the data space. (Strictly speaking, we could instead choose an initial set of responsibilities for each data point; it would just have the effect of swapping the order of steps 2 and 3.)
Perform the Expectation Step described above, treating the parameters \(\theta\) as fixed.
Perform the Maximization Step described above, treating the responsibilities as fixed.
Repeat steps 2 and 3 until convergence. We call one pass through both the E and M Steps an iteration. Since we are using soft responsibilities, we need to modify our definition of convergence slightly. Choose some tolerance value, \(\epsilon\). Stop when no responsibility changes by more than \(\epsilon\). (There are many possibilities for determining how to stop, including setting a limit on the number of iterations.)
You will probably have noticed that the EM algorithm has many parallels with K-means. They are very closely related algorithms, where a GMM has the ability to fit more varied distributions at the cost of more complexity and computation.
Using the figures below, you can step through the EM algorithm with up to \(K=5\) for each of the two datasets we began this section with. You should see that, compared to K-means, the GMM is much more easily able to capture clusters that are not circular. And, because each Gaussian has its own covariance parameter, a GMM is also better able to capture clusters that are of different sizes from each other. Note that, because we are maximizing log-likelihood rather than minimizing cost, we expect the line graphs to increase in value with each iteration, rather than decrease.
GMM, One Step at a Time
GMM, One Step at a Time, on a Tricky Dataset
4 Summary
In this section, we have introduced the Gaussian Mixture Model and the Expectation-Maximization algorithm for optimizing its parameters. More complex and more powerful than K-means, a GMM is better suited to capturing clusters of varying sizes and irregular (non-spherical) shapes.
Despite this chapter on clustering being a different kind of machine learning from the supervised learning topics we saw in previous chapters, we can make several connections back to the Fundamental Ideas of machine learning. We have seen how, via block coordinate descent, clustering tasks can be framed as an optimization problem (Idea #1: Learning is Optimization). We have also seen that both K-means and GMMs are explicitly geometrical in nature, dealing with distances in the data space (Idea #4: ML Describes Geometric Processes), and that a GMM is a generative, probabilistic model (Idea #5: ML Demands a Probabilistic Lens).
In the next chapter, we will discuss a completely different kind of unsupervised learning: Principal Component Analysis (PCA). Like clustering, PCA gives us insight into the structure of unlabelled data, and in addition, it finally gives us a tool to fight back against the curse of dimensionality.