Gaussian Mixtures and Expectation-Maximization, machine learning lesson

A sample of measurements can mix two overlapping subpopulations, so that a value lying between them could plausibly come from either. A hard clustering such as k-means must still put that value in exactly one group. A Gaussian mixture model instead describes the whole distribution as a weighted sum of normal densities and gives every observation a probability of coming from each component. This lesson defines the model and its hidden component labels, computes those probabilities, carries out one exact iteration of the expectation-maximization algorithm on four constructed observations, and shows why the likelihood that the algorithm increases can grow without bound.

Key terms

Gaussian mixture model

Definition A Gaussian mixture model with $K$ components gives a feature vector $\mathbf{x}$ in $\mathbb{R}^p$ the density

\begin{aligned} &p(\mathbf{x}) \\[0.4em] &\quad = \sum_{k=1}^{K} \pi_k \, \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \Sigma_k), \end{aligned}

where the mixing weights satisfy $\pi_k \ge 0$ and $\pi_1 + \dots + \pi_K = 1$ , each mean $\boldsymbol{\mu}_k$ lies in $\mathbb{R}^p$ , each covariance matrix $\Sigma_k$ is symmetric positive definite, and

\begin{aligned} &\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \Sigma) \\[0.4em] &\quad = \frac{\exp\bigl(-\tfrac12\, q(\mathbf{x})\bigr)}{(2\pi)^{p/2} \det(\Sigma)^{1/2}} \end{aligned}

is the normal density with $q(\mathbf{x}) = (\mathbf{x} - \boldsymbol{\mu})^{\top} \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu})$ . In one dimension a component has a mean $\mu_k$ and a variance $\sigma_k^2 > 0$ , and

\begin{aligned} &\mathcal{N}(x \mid \mu, \sigma^2) \\[0.4em] &\quad = \frac{1}{\sqrt{2\pi\sigma^2}}\, e^{-(x - \mu)^2/(2\sigma^2)} . \end{aligned}

The parameters $\boldsymbol{\theta}$ are all weights, means and covariances, and the number of components $K$ is chosen by the analyst, for example by comparing log likelihoods on validation data. Exchanging the labels of two components together with their parameters gives the same density, so a component label has no meaning of its own.

In plain words The model blends $K$ bell-shaped densities, each with its own center, spread and share of the data. A blend of two components need not have two humps. With equal weights and a common variance $\sigma^2$ , two one-dimensional components give a single hump whenever their means are at most $2\sigma$ apart. A component is a piece of a density model, so it need not correspond to a real subgroup, and one skewed group can need several components.

Example Four constructed one-dimensional observations, without units, are $-2$ , $-1$ , $2$ and $3$ . A starting model has $K = 2$ , weights $\pi_1 = \pi_2 = \tfrac12$ , means $\mu_1 = -1$ and $\mu_2 = 2$ , and variances $\sigma_1^2 = \sigma_2^2 = 1$ . Its means are 3 apart, more than $2\sigma = 2$ , so its density has two humps. At $x = -1$ the component densities are about 0.398942 and 0.004432, and the mixture density is their average, about 0.201687. Swapping the two components gives the same density at every $x$ .

Latent variable

Definition A latent variable is a random variable of a model that is never observed. In a Gaussian mixture each observation $i$ has a latent component label $Z_i$ with values in $\{1, \dots, K\}$ , and the model generates the data in two stages,

\begin{gathered} P(Z_i = k) = \pi_k, \\[0.4em] \mathbf{X}_i \mid Z_i = k \sim \mathcal{N}(\boldsymbol{\mu}_k, \Sigma_k), \end{gathered}

with the pairs $(\mathbf{X}_i, Z_i)$ independent across $i$ . The joint density of an observation and its label is $\pi_k\, \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \Sigma_k)$ , and summing it over $k$ gives back the mixture density, so the mixture is the distribution of $\mathbf{X}_i$ once the label is ignored. The pairs $(\mathbf{x}_i, z_i)$ are called the complete data and the observations $\mathbf{x}_i$ alone the observed data.

In plain words Picture each observation as made in two moves. First a hidden component is drawn with probabilities $\pi_k$ , then a value is drawn from that component’s normal distribution, and only the value is recorded. If the labels were known, fitting would reduce to an ordinary proportion, mean and variance for each group. Not knowing them is what makes the problem harder. The hidden label is part of a model’s account of how the data could arise, not a measured fact.

Example In the starting model, $P(Z = 1) = \tfrac12$ , and given $Z = 1$ the observation is drawn from $\mathcal{N}(-1, 1)$ . The observation $-1$ could also have come from component 2, since a normal density is positive everywhere. It is only much less likely. If the labels of the four observations were known to be 1, 1, 2 and 2, the group means would be $-1.5$ and $2.5$ , and the group variances, dividing by the group size, would be 0.25 each.

Responsibility

Definition Write $\mathcal{N}_k(\mathbf{x})$ for $\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \Sigma_k)$ . Under parameters $\boldsymbol{\theta}$ , the responsibility of component $k$ for observation $\mathbf{x}_i$ is the posterior probability of its label given the observation, obtained from Bayes’ rule with densities,

\begin{aligned} r_{ik} &= P(Z_i = k \mid \mathbf{x}_i) \\[0.4em] &= \frac{\pi_k\, \mathcal{N}_k(\mathbf{x}_i)}{\sum_{j=1}^{K} \pi_j\, \mathcal{N}_j(\mathbf{x}_i)} . \end{aligned}

Each $r_{ik}$ lies in $[0, 1]$ , and $r_{i1} + \dots + r_{iK} = 1$ . A soft assignment gives observation $i$ the whole vector $(r_{i1}, \dots, r_{iK})$ . A hard assignment keeps only a component with the largest responsibility, with a rule for ties. The values $\mathcal{N}_k(\mathbf{x}_i)$ are densities, not probabilities. They can exceed 1 and need not add to 1 over $k$ , and they become probabilities only after weighting by $\pi_k$ and normalizing. With equal weights and a common covariance $\sigma^2 I$ , the largest responsibility belongs to the nearest mean, which is the k-means assignment rule. In computation, work with $\log \pi_k + \log \mathcal{N}_k(\mathbf{x}_i)$ , subtract the largest of these values, exponentiate and normalize, which gives the same result without underflow.

In plain words A responsibility answers the question of how likely each component is to have produced this observation, given the current model. The weights play the part of prior probabilities and the densities say how well each component explains the value. A point near one mean and far from the others gets a share near 1, and a point between overlapping components gets sizable shares in several.

Example For $x = -1$ under the starting model, the densities are $1/\sqrt{2\pi} \approx 0.398942$ and $e^{-4.5}/\sqrt{2\pi} \approx 0.004432$ , and the equal weights cancel, so $r_{11} = 1/(1 + e^{-4.5}) \approx 0.989013$ . For any $x$ in this model,

\begin{aligned} &\log \frac{\mathcal{N}(x \mid -1, 1)}{\mathcal{N}(x \mid 2, 1)} \\[0.4em] &\ = \tfrac12 \bigl[(x - 2)^2 - (x + 1)^2\bigr] \\[0.4em] &\ = 1.5 - 3x , \end{aligned}

so the responsibility of component 1 is $1/(1 + e^{-(1.5 - 3x)})$ , the logistic sigmoid of $1.5 - 3x$ .

Observation	$r_{i1}$	$r_{i2}$
$-2$	$0.999447$	$0.000553$
$-1$	$0.989013$	$0.010987$
$2$	$0.010987$	$0.989013$
$3$	$0.000553$	$0.999447$

A hard assignment gives the labels 1, 1, 2 and 2 and discards the small chance of the other component. At $x = 0.5$ , halfway between the means, the two densities are equal, so with equal weights the responsibilities are $\tfrac12$ each and a hard assignment needs its tie rule. With weights $\tfrac34$ and $\tfrac14$ the responsibility of component 1 at $0.5$ becomes $\tfrac34$ , all of it from the prior.

E step and M step

Definition One iteration of the expectation-maximization algorithm for a Gaussian mixture starts from current parameters and has two steps. The E step computes every responsibility $r_{ik}$ from the current parameters. The M step holds those responsibilities fixed and sets

\begin{aligned} N_k &= \sum_{i=1}^{n} r_{ik}, \\[0.4em] \pi_k^{\text{new}} &= N_k / n, \\[0.4em] \boldsymbol{\mu}_k^{\text{new}} &= \frac{1}{N_k} \sum_{i=1}^{n} r_{ik}\, \mathbf{x}_i, \\[0.4em] \Sigma_k^{\text{new}} &= \frac{1}{N_k} \sum_{i=1}^{n} r_{ik}\, \mathbf{d}_{ik} \mathbf{d}_{ik}^{\top}, \end{aligned}

where $\mathbf{d}_{ik} = \mathbf{x}_i - \boldsymbol{\mu}_k^{\text{new}}$ . In one dimension the last line is a weighted variance with $(x_i - \mu_k^{\text{new}})^2$ in place of $\mathbf{d}_{ik} \mathbf{d}_{ik}^{\top}$ . These values maximize the expected complete-data log likelihood

\begin{aligned} Q(\boldsymbol{\theta}) &= \sum_{i=1}^{n} \sum_{k=1}^{K} r_{ik} \bigl[\log \pi_k \\[0.4em] &\qquad + \log \mathcal{N}_k(\mathbf{x}_i)\bigr] \end{aligned}

over all parameters, for the fixed responsibilities, provided every $N_k > 0$ and every $\Sigma_k^{\text{new}}$ is positive definite. Setting the derivative with respect to $\boldsymbol{\mu}_k$ to zero gives the weighted mean, and maximizing the weighted sum of $\log \pi_k$ under $\pi_1 + \dots + \pi_K = 1$ gives $N_k / n$ . The effective count $N_k$ need not be a whole number, and the covariance divides by $N_k$ , as maximum likelihood does, not by $N_k - 1$ . A component with $N_k = 0$ or a singular $\Sigma_k^{\text{new}}$ has no valid update and needs a declared rule.

In plain words The E step asks, for each observation, how much each component explains it. The M step refits every component as if each observation belonged to it with that fractional weight, giving a weighted proportion, a weighted mean and a weighted spread. With hard labels these would be the ordinary group proportions, means and variances.

Example From the responsibility table, $N_1 = 0.999447 + 0.989013 + 0.010987 + 0.000553 = 2$ , exactly, because the shares at $-2$ and $3$ and at $-1$ and $2$ mirror each other. So $N_2 = 2$ and both weights stay $\tfrac12$ . The weighted sum $\sum_i r_{i1} x_i$ is about $-2.964275$ , so $\mu_1^{\text{new}} \approx -1.482138$ . The data and the starting model are symmetric about $0.5$ , so $\mu_2^{\text{new}} = 1 - \mu_1^{\text{new}} \approx 2.482138$ . The weighted squared deviations of component 1 from its new mean are about 0.268033, 0.229903, 0.133220 and 0.011105, so

\begin{aligned} (\sigma_1^2)^{\text{new}} &\approx \frac{0.642261}{2} \\[0.4em] &\approx 0.321130 , \end{aligned}

and component 2 gets the same variance. The means moved outward toward $-1.5$ and $2.5$ , and the variances shrank from 1 to about 0.321.

Expectation-maximization algorithm

Definition The expectation-maximization (EM) algorithm is an iterative method for maximum likelihood estimation in a model with observed data $\mathbf{x}$ , latent variables $\mathbf{z}$ and parameters $\boldsymbol{\theta}$ . It seeks a large value of the log likelihood $\ell(\boldsymbol{\theta}) = \log p(\mathbf{x} \mid \boldsymbol{\theta})$ , where the density of the observed data has the latent variables summed out. For a Gaussian mixture with independent observations, $\ell(\boldsymbol{\theta}) = \sum_i \log p(\mathbf{x}_i)$ . Starting from $\boldsymbol{\theta}^{(0)}$ , iteration $t$ computes the posterior distribution of $\mathbf{z}$ under $\boldsymbol{\theta}^{(t)}$ in the E step, and in the M step chooses $\boldsymbol{\theta}^{(t+1)}$ to maximize the expected complete-data log likelihood under that posterior. It stops when the change in log likelihood or in the parameters falls below a stated tolerance, or after a stated maximum number of iterations. For every probability distribution $q$ on the latent values that is positive at each value, Jensen’s inequality for the concave logarithm gives

\begin{aligned} &\log p(\mathbf{x} \mid \boldsymbol{\theta}) \\[0.4em] &\ \ge \sum_{\mathbf{z}} q(\mathbf{z}) \log \frac{p(\mathbf{x}, \mathbf{z} \mid \boldsymbol{\theta})}{q(\mathbf{z})} , \end{aligned}

with equality when $q$ is the posterior under $\boldsymbol{\theta}$ . Let $B_t(\boldsymbol{\theta})$ be the right side with $q$ equal to the posterior under $\boldsymbol{\theta}^{(t)}$ . The M step maximizes $B_t$ , because $B_t$ differs from the expected complete-data log likelihood by a term that does not depend on $\boldsymbol{\theta}$ . With exact steps,

\begin{aligned} &\ell(\boldsymbol{\theta}^{(t+1)}) \ge B_t(\boldsymbol{\theta}^{(t+1)}) \\[0.4em] &\quad \ge B_t(\boldsymbol{\theta}^{(t)}) \\[0.4em] &\quad = \ell(\boldsymbol{\theta}^{(t)}) . \end{aligned}

So the log likelihood never decreases, and when it is bounded above its values converge. The limit can be a local maximum or another stationary point rather than a global maximum, and different starting values can lead to different fits.

In plain words EM alternates between guessing the hidden labels softly from the current parameters and refitting the parameters to those soft guesses. Each guess builds a lower bound on the log likelihood that touches it at the current parameters, and the refit climbs that bound, so the log likelihood cannot fall. That makes EM stable, but a likelihood that has stopped rising only shows that the algorithm reached a stationary point. Running several starting values and comparing the results is standard practice.

Example With natural logarithms, the log likelihood of the four observations is about $-7.425142$ under the starting model and $-5.735514$ after the iteration above. A second iteration brings the means to about $-1.5$ and $2.5$ and the variances to about 0.25, with log likelihood about $-5.675754$ , and later iterations change these values only in far decimal places. This limit is a local maximum. It is not the largest possible value of the likelihood, because, as the next card shows, this likelihood has no largest value at all. Starting instead from the means 2 and $-1$ reaches the same fit with the two labels exchanged.

Covariance collapse

Definition Covariance collapse is the failure in which a component’s covariance matrix approaches a singular matrix while the component concentrates on very few observations. In one dimension its variance tends to 0 while its mean sits on one observation. In $p$ dimensions it can happen whenever a component’s weight concentrates on $p$ or fewer points, since such points lie in a lower-dimensional flat. For a mixture with $K \ge 2$ , set $\mu_1 = x_1$ with $\pi_1 > 0$ , and hold the other components fixed with $\pi_2 > 0$ . As $\sigma_1^2 \to 0$ ,

\begin{aligned} &\pi_1\, \mathcal{N}(x_1 \mid x_1, \sigma_1^2) \\[0.4em] &\quad = \frac{\pi_1}{\sqrt{2\pi \sigma_1^2}} \to \infty , \end{aligned}

while the mixture density at every other observation stays above the positive value $\pi_2\, \mathcal{N}_2(x_i)$ . So the log likelihood grows without bound, an unconstrained maximum likelihood estimate does not exist, and EM can drift toward such a degenerate fit. The usual safeguards are a floor $\sigma_k^2 \ge c$ with a stated $c > 0$ , or a small positive amount added to the diagonal of every covariance matrix, a prior on the covariances, several starting values, a rule for components whose effective count $N_k$ falls near 1, and a comparison of fits by log likelihood on validation data. With a floor, every one-dimensional component density is at most $1/\sqrt{2\pi c}$ , so the likelihood is bounded. A fit with a safeguard maximizes a constrained or penalized objective, and a report should say which.

In plain words The likelihood rewards a component for placing high density where observations are. A component that wraps itself around a single observation can make that density as large as it likes, so the likelihood can grow without limit while the model becomes useless for new data. A record-high likelihood is then a symptom, not a better model. A tiny component that claims one observation with almost no spread is nearly always this failure rather than a discovered group.

Example For the four observations, give each component weight $\tfrac12$ , put component 2 at mean $0.5$ with variance 4, and put component 1 on the observation $-2$ with variance $s^2$ . The log likelihood is about $-5.264$ for $s^2 = 10^{-4}$ , about $-0.661$ for $s^2 = 10^{-8}$ and about $8.549$ for $s^2 = 10^{-16}$ . Dividing $s^2$ by 100 raises it by about $\ln 10 \approx 2.303$ , and already at $s^2 = 10^{-4}$ it exceeds the $-5.676$ of the sensible fit from the previous card. With a floor $c = 0.01$ on the variances, every component density is at most $1/\sqrt{0.02\pi} \approx 3.989$ , so the log likelihood of four observations cannot exceed $4 \ln 3.989 \approx 5.535$ .

Check yourself

A two-component mixture has weights 0.2 and 0.8. At an observation $x$ , the two component densities are 2.0 and 0.1. Find both responsibilities and the hard assignment. A classmate says the value 2.0 shows an error, since a probability cannot exceed 1. Is that right?

Show the answer

The weighted densities are $0.2 \times 2.0 = 0.4$ and $0.8 \times 0.1 = 0.08$ , with total $0.48$ , so

\begin{aligned} r_1 &= \frac{0.4}{0.48} = \frac{5}{6} \approx 0.833, \\[0.4em] r_2 &= \frac{0.08}{0.48} = \frac{1}{6} \approx 0.167 . \end{aligned}

The shares add to 1 and the hard assignment is component 1. The classmate is wrong. The value 2.0 is a density, a probability per unit of $x$ , not a probability. A normal density reaches $1/(\sigma\sqrt{2\pi})$ at its mean, which is above 2 whenever $\sigma < 1/(2\sqrt{2\pi}) \approx 0.1995$ . For example, with $\sigma = 0.15$ the peak is about 2.660. Probabilities come from areas under the density, and the responsibilities, which are probabilities, stay in $[0, 1]$ .

Four observations $0$ , $2$ , $4$ and $10$ have responsibilities $1$ , $0.75$ , $0.25$ and $0$ for component 1 of a two-component mixture. Carry out the M step for component 1, then compare its mean and variance with those from a hard assignment to the component with the larger responsibility.

Show the answer

The effective count is $N_1 = 1 + 0.75 + 0.25 + 0 = 2$ , so $\pi_1 = 2/4 = 0.5$ . The weighted mean is

\mu_1 = \frac{0 + 1.5 + 1 + 0}{2} = 1.25 .

The squared deviations from $1.25$ are $1.5625$ , $0.5625$ , $7.5625$ and $76.5625$ , and weighting them by the responsibilities gives $1.5625$ , $0.421875$ , $1.890625$ and $0$ . Their sum is $3.875$ , so $\sigma_1^2 = 3.875/2 = 1.9375$ .

A hard assignment gives component 1 the observations 0 and 2, whose responsibilities are at least 0.75, and gives 4 and 10 to component 2. The hard mean is 1 and the hard variance, dividing by 2, is 1. The soft update lets the observation 4 pull the mean up with a quarter of its weight and widen the spread, and it counts 2 as only three quarters of a member. For component 2 the soft update gives $N_2 = 2$ , mean $6.75$ and variance $10.9375$ , against the hard mean 7 and variance 9.

A three-component mixture is fitted to 200 one-dimensional training observations from several random starts. One run reports a component with weight 0.005, a mean equal to the observation 7.31 and variance $10^{-10}$ , and its training log likelihood is far above that of every other run. What has happened, and what should be done?

Show the answer

The effective count of that component is $0.005 \times 200 = 1$ , so it describes a single observation. Its standard deviation is $10^{-5}$ , and its density at its own mean is

\frac{1}{\sqrt{2\pi \times 10^{-10}}} \approx 39{,}894 .

Weighted by 0.005 this is about 199.5, so the mixture density at 7.31 is at least 199.5, and that single observation contributes at least $\ln 199.5 \approx 5.30$ to the log likelihood. Each further halving of the variance would multiply the spike’s height by $\sqrt{2}$ and add about $\ln \sqrt{2} \approx 0.347$ , so this direction has no maximum to find. This is covariance collapse. The high training likelihood is an artifact, not evidence of a real group of identical values, and on new observations the spike contributes almost nothing, since a new value almost never falls within a few multiples of $10^{-5}$ of 7.31.

The run should be discarded or refitted with a safeguard, such as a stated floor on the variances or a small amount added to each variance, and the remaining runs and choices of $K$ compared by log likelihood on validation data rather than on the training data. If collapse keeps recurring, the model may simply have more components than the data support.

Gaussian Mixtures and Expectation-Maximization

Key terms

Gaussian mixture model

Latent variable

Responsibility

E step and M step

Expectation-maximization algorithm

Covariance collapse

Check yourself

How this lesson connects

Builds on

Used later in