4.3 Gaussian mixtures and expectation-maximization

Motivation

A point between two groups need not have one overwhelmingly plausible membership. Hard assignment records only the selected group. A mixture model instead combines several probability densities and retains, for each observation, a probability of membership in every component. These probabilities become weights when estimating the components again.

The challenge is that the likelihood contains a logarithm of a sum. Even when each Gaussian is simple, the parameters of different components are coupled through that sum. Expectation-maximization, abbreviated EM, introduces responsibilities that separate the parameter update into weighted estimation problems. Its two steps optimize a lower bound that touches the observed log-likelihood at the current parameters.

We will derive those weights and updates before calculating a complete iteration. The distinction between a soft probability and a hard label is central: maximizing the density of one selected component is a different objective from maximizing a mixture density. Covariance normalization matters for the same reason. The denominator in a weighted covariance follows from the likelihood being optimized.

Geometric soft weights

Soft weights can also come from geometry. Let distinct sites $S=\{s_1,\ldots,s_K\}$ have old Voronoi cells $V_j$, with each distance tie assigned to one cell. Insert a query $x\notin S$ and let $V_x$ be its new Voronoi cell. Assume this cell is bounded and has finite, positive volume. Define

$$0<\operatorname{vol}(V_x)<\infty,\qquad w_j=\frac{\operatorname{vol}(V_x\cap V_j)} {\operatorname{vol}(V_x)},\qquad \sum_{j=1}^K w_j=1.$$

The old cells partition the new cell, so their overlap volumes add to its volume. Every weight is nonnegative; a site with no positive-volume overlap receives zero. The fractions measure how much of the query's new region previously belonged to each old site.

In one dimension, volume is length. For two sites and an inserted query,

$$S=\{0,4\},\quad x=1,\quad V_x=\left[\frac12,\frac52\right],\quad (w_1,w_2)=\left(\frac{3/2}{2},\frac{1/2}{2}\right) =\left(\frac34,\frac14\right).$$

The old boundary is $2$. Its left and right overlaps with the new interval have lengths $3/2$ and $1/2$, respectively, while the new interval has length $2$. These are geometric allocation weights. Gaussian responsibilities, derived next, are conditional probabilities under a fitted mixture model. Both constructions normalize across components, but their definitions use different information.

Definitions

Let $x_1,\ldots,x_n\in\mathbb R^d$ be independent observations, and fix a number $K$ of Gaussian components. Component $k$ has weight $\pi_k>0$, mean $\mu_k\in\mathbb R^d$, and positive-definite covariance $\Sigma_k$, with $\sum_{k=1}^K\pi_k=1$. Write $\theta$ for all these parameters. All logarithms below are natural logarithms.

Definition

Gaussian mixture and observed log-likelihood

The Gaussian density and mixture density are

$$\phi_d(x;\mu,\Sigma) =\frac{\exp\!\left(-\frac12(x-\mu)^\top\Sigma^{-1}(x-\mu)\right)} {(2\pi)^{d/2}\det(\Sigma)^{1/2}}, \qquad p_\theta(x)=\sum_{k=1}^K\pi_k\phi_d(x;\mu_k,\Sigma_k).$$

The observed log-likelihood is

$$\ell(\theta)=\sum_{i=1}^n\log p_\theta(x_i).$$

A latent component label $z_i$ takes value $k$ with probability $\pi_k$. The component density describes $x_i$ conditional on that label.

The mixture weight is part of the model, rather than a normalization applied after choosing a component. Its effect remains even when two component densities happen to have equal values at an observation. The sum over all possible labels is what makes the observed density.

If the labels were observed as well, the log-likelihood would be a sum of terms $\log\pi_{z_i}+\log\phi_d(x_i;\mu_{z_i},\Sigma_{z_i})$. Equivalently, a zero-or-one indicator for each component would select one term per observation. The unknown labels prevent us from grouping the observations that way. EM takes the conditional expectation of those indicators instead, producing the responsibilities that appear in $Q$. The fractional coefficients have a probabilistic origin: they are expectations of latent indicators under the current model.

This also distinguishes the three functions used in the derivation. The observed log-likelihood sums over possible labels inside each logarithm. The complete-data log-likelihood would use known labels. The function $Q$ averages that complete-data expression under the current posterior. Keeping track of which quantities are observed, latent, or held fixed explains why an apparently simple weighted Gaussian fit can improve a more complicated mixture likelihood.

Definition

Responsibilities and effective component mass

At current parameters $\theta^{\mathrm{old}}$, Bayes' rule gives

$$r_{ik} =\Pr_{\theta^{\mathrm{old}}}(z_i=k\mid x_i) =\frac{\pi_k^{\mathrm{old}} \phi_d(x_i;\mu_k^{\mathrm{old}},\Sigma_k^{\mathrm{old}})} {\sum_{j=1}^K\pi_j^{\mathrm{old}} \phi_d(x_i;\mu_j^{\mathrm{old}},\Sigma_j^{\mathrm{old}})}, \qquad N_k=\sum_{i=1}^n r_{ik}.$$

Each row satisfies $\sum_k r_{ik}=1$. With positive mixture weights and positive-definite Gaussian covariances, every $r_{ik}$ is positive, so $N_k>0$. Also $\sum_kN_k=n$.

$N_k$ is the expected number of observations attributed to component $k$ under the current posterior labels. It can be fractional. A responsibility is conditional on the current model and the observation; it is not a known class label supplied with the data.

There are two different normalizations. The E-step normalizes across components for each individual observation, so its row sums are one. The mean update normalizes across observations for one component, so its denominator is the column sum $N_k$. Dividing a mean update by $n$ instead would treat responsibility assigned to other components as if it had contributed to this component's total mass. The two denominators correspond to different conditional questions.

A Gaussian density at a point is not itself a membership probability. It can exceed one and changes scale when measurement units change. Responsibilities are dimensionless probabilities obtained by comparing the weighted densities at the same observation. For general covariance, that comparison includes both the Mahalanobis distance in the exponent and the determinant in the normalization. A more spread-out component has a larger volume and a lower peak density, while $\pi_k$ sets its overall mixture weight. For equal weights and a common isotropic covariance, ranking component densities reduces to comparing Euclidean distances.

Theorem/Proposition

Theorem

The Gaussian M-step uses effective mass

Hold the responsibilities fixed and define the expected complete-data log-likelihood

$$Q(\theta\mid\theta^{\mathrm{old}}) =\sum_{i=1}^n\sum_{k=1}^K r_{ik}\left[\log\pi_k+\log\phi_d(x_i;\mu_k,\Sigma_k)\right].$$

When the weighted covariance matrices below are positive definite, its maximizers are

$$\pi_k^{\mathrm{new}}=\frac{N_k}{n},\qquad \mu_k^{\mathrm{new}}=\frac1{N_k}\sum_{i=1}^nr_{ik}x_i,$$$$\Sigma_k^{\mathrm{new}} =\frac1{N_k}\sum_{i=1}^nr_{ik} (x_i-\mu_k^{\mathrm{new}}) (x_i-\mu_k^{\mathrm{new}})^\top.$$

For a model with a prescribed common covariance $\sigma^2I_d$, $\sigma^2>0$, keep that covariance fixed and use the same weight and mean updates.

The weighted covariance is positive definite exactly when its weighted centred observations span all feature directions. Otherwise the displayed matrix is singular and is not a positive-definite Gaussian covariance. The fixed-variance model gives a useful setting in which every update above remains a valid model without estimating covariance matrices.

Theorem

A valid EM iteration does not decrease observed log-likelihood

Compute the exact posterior responsibilities at $\theta^{\mathrm{old}}$. Choose valid new parameters satisfying

$$Q(\theta^{\mathrm{new}}\mid\theta^{\mathrm{old}}) \ge Q(\theta^{\mathrm{old}}\mid\theta^{\mathrm{old}}).$$

Then

$$\ell(\theta^{\mathrm{new}})\ge\ell(\theta^{\mathrm{old}}).$$

An exact M-step maximizes $Q$; an update that only increases it is a generalized EM step. The conclusion concerns observed log-likelihood values. Global optimality and convergence of the parameter vectors are separate claims.

Proof sketch or proof idea

Deriving the weighted updates

For the weights, maximize $\sum_kN_k\log\pi_k$ subject to $\sum_k\pi_k=1$. The stationary equation for a multiplier $\lambda$ is $N_k/\pi_k=\lambda$. Summing over $k$ gives $\lambda=n$, hence $\pi_k=N_k/n$.

For each mean, the terms depending on $\mu_k$ are a negative weighted sum of squared Mahalanobis distances. Differentiating gives

$$\Sigma_k^{-1}\sum_i r_{ik}(x_i-\mu_k)=0.$$

Since $\Sigma_k$ is invertible and $N_k>0$, the mean must be the weighted average in the proposition. With a prescribed common isotropic covariance, this is the familiar weighted least-squares minimizer. Every observation contributes to each component in proportion to its responsibility, rather than only to the component with the largest one.

To derive the covariance update, let

$$S_k=\sum_i r_{ik}(x_i-\mu_k^{\mathrm{new}}) (x_i-\mu_k^{\mathrm{new}})^\top,\qquad \Lambda_k=\Sigma_k^{-1}.$$

The relevant part of $Q$ is

$$\frac{N_k}{2}\log\det\Lambda_k -\frac12\operatorname{tr}(\Lambda_kS_k).$$

Its stationary equation is $N_k\Lambda_k^{-1}=S_k$. For positive-definite $S_k$, concavity of the log determinant gives the maximizer $\Sigma_k=S_k/N_k$. This is a likelihood estimate, so the denominator is $N_k$, with no subtraction of one. Omitting that denominator would produce a weighted scatter matrix with the wrong scale.

The touching lower bound

For any positive probability row $q_i=(q_{i1},\ldots,q_{iK})$, Jensen's inequality gives

$$\log\sum_k\pi_k\phi_d(x_i;\mu_k,\Sigma_k) \ge\sum_kq_{ik}\log \frac{\pi_k\phi_d(x_i;\mu_k,\Sigma_k)}{q_{ik}}.$$

Define the summed right-hand side by

$$\mathcal F(q,\theta) =\sum_{i,k}q_{ik} \left[\log\pi_k+\log\phi_d(x_i;\mu_k,\Sigma_k)-\log q_{ik}\right].$$

The exact gap is

$$\ell(\theta)-\mathcal F(q,\theta) =\sum_i\operatorname{KL}\!\left(q_i\Vert r_i(\theta)\right)\ge0,$$

where $\operatorname{KL}(q\Vert r)=\sum_kq_k\log(q_k/r_k)$ and $r_i(\theta)$ denotes the posterior row at parameters $\theta$. The gap is zero when $q_i=r_i(\theta)$.

The E-step chooses $q=r(\theta^{\mathrm{old}})$, making the bound touch the current observed log-likelihood. Holding $q$ fixed, its entropy $-\sum_{i,k}q_{ik}\log q_{ik}$ is constant, so increasing $Q$ increases $\mathcal F$. Therefore

$$\ell(\theta^{\mathrm{new}}) \ge\mathcal F(r^{\mathrm{old}},\theta^{\mathrm{new}}) \ge\mathcal F(r^{\mathrm{old}},\theta^{\mathrm{old}}) =\ell(\theta^{\mathrm{old}}).$$

The two steps have different jobs: the E-step makes the bound exact at the current model, while the M-step improves the model under that fixed bound. Replacing posterior probabilities by hard labels removes the touching equality on which this argument depends.

The responsibilities remain fixed throughout one M-step. Every updated mean, covariance, and weight in that step uses the same posterior matrix calculated from the complete old parameter set. The covariance formula centres its residuals at the new weighted mean because it maximizes the same fixed-$Q$ problem. Once the new model has been formed, the next E-step recomputes the posterior from that model. This ordering makes the quantities in the chain of inequalities unambiguous: its middle comparison changes parameters while retaining one probability matrix.

The bound and KL identity also extend to probability rows containing zero entries, with $0\log0=0$ and zero-weighted terms interpreted by continuity. This lets us compare the tight posterior bound with the looser bound obtained from a hard assignment.

Worked examples

Worked example

1. Equal component densities retain unequal mixture weights

Consider a one-dimensional mixture with means $\mu_1=0$, $\mu_2=2$, common variance $\sigma^2=1$, and weights $(\pi_1,\pi_2)=(3/4,1/4)$. At $x=1$, the distances to the two means are equal, so the Gaussian densities are equal. Their common value cancels from Bayes' rule:

$$r_1=\frac{(3/4)\phi_1(1;0,1)} {(3/4)\phi_1(1;0,1)+(1/4)\phi_1(1;2,1)} =\frac34,\qquad r_2=\frac14.$$

The midpoint has equal compatibility with the component shapes, but the model gives the first component three times the prior weight. Normalizing the two densities without their mixture weights would instead give $1/2$ and $1/2$, which are posterior probabilities for a different mixture. Soft membership therefore involves both density and weight.

Worked example

2. One complete E-step and M-step with fixed variance

Take observations $x_1=0$, $x_2=2$, initial means $(0,2)$, equal mixture weights, and prescribed common variance $1$. Set

$$a=\frac1{1+e^{-2}}\approx0.880797078.$$

At each observation, the squared distance to the other mean exceeds the squared distance to its matching mean by $4$, giving the density ratio $e^{-2}$. Thus

$$R=\begin{bmatrix}a&1-a\\1-a&a\end{bmatrix}, \qquad N_1=N_2=1.$$

The M-step gives

$$\pi^{\mathrm{new}}=(1/2,1/2),\qquad \mu^{\mathrm{new}}=(2(1-a),2a) \approx(0.238405844,1.761594156).$$

Although each observation initially lies at one component mean, its positive responsibility for the other component pulls both updated means inward. The common variance stays exactly $1$. Using the full Gaussian density, the observed log-likelihood changes from $-2.970315405$ to $-2.886562696$, an increase of about $0.083752710$.

This NumPy code performs the same single iteration. Logarithmic normalization computes the posterior from all components; it does not select a component before forming the denominator.

import numpy as np

X = np.array([0.0, 2.0])
means = np.array([0.0, 2.0])
weights = np.array([0.5, 0.5])
variance = 1.0

def log_components(means, weights):
    return (np.log(weights)[None, :]
            - 0.5 * np.log(2.0 * np.pi * variance)
            - 0.5 * (X[:, None] - means[None, :]) ** 2 / variance)

joint = log_components(means, weights)
normalizer = np.logaddexp.reduce(joint, axis=1)
responsibilities = np.exp(joint - normalizer[:, None])
mass = responsibilities.sum(axis=0)
new_weights = mass / len(X)
new_means = (responsibilities.T @ X) / mass
before = normalizer.sum()
after = np.logaddexp.reduce(
    log_components(new_means, new_weights), axis=1
).sum()
assert np.allclose(responsibilities.sum(axis=1), 1.0)
assert after >= before
print("responsibilities =", np.round(responsibilities, 6).tolist())
print("mass =", np.round(mass, 6).tolist())
print("weights =", np.round(new_weights, 6).tolist())
print("means =", np.round(new_means, 6).tolist())
print(f"log_likelihood: {before:.6f} -> {after:.6f}")

The output is

responsibilities = [[0.880797, 0.119203], [0.119203, 0.880797]]
mass = [1.0, 1.0]
weights = [0.5, 0.5]
means = [0.238406, 1.761594]
log_likelihood: -2.970315 -> -2.886563

Worked example

3. Covariance collapse can make likelihood unbounded

Take observations $0$ and $3$. Give two components weight $1/2$ each, means $0$ and $3$, and variances $\varepsilon^2$ and $1$, respectively. For every $\varepsilon>0$, both Gaussian components are valid. Their mixture densities satisfy

$$p_\theta(0)\ge\frac1{2\sqrt{2\pi}\varepsilon}, \qquad p_\theta(3)\ge\frac1{2\sqrt{2\pi}}.$$

Consequently,

$$\ell(\theta)\ge-\log\varepsilon-\log(8\pi) \longrightarrow+\infty \qquad(\varepsilon\downarrow0).$$

Shrinking a component around one observation creates an arbitrarily tall density there, while the other component continues to explain the other observation. This unrestricted likelihood has no finite maximum. The issue is the statistical objective itself, not an instruction to evaluate a Gaussian with a singular covariance. Prescribing a positive common variance, as in Example 2, specifies a different model in which this collapse is impossible.

Reading the progress of EM

In the prescribed isotropic-variance model, each component density is at most $(2\pi\sigma^2)^{-d/2}$. A weighted average of those densities has the same upper bound, so

$$\ell(\theta)\le-\frac{nd}{2}\log(2\pi\sigma^2).$$

The nondecreasing likelihood values therefore have a finite limit in this model. This statement concerns a real sequence of objective values. It does not identify the best possible mixture or establish a limit for every parameter coordinate. Component-label permutations also describe the same mixture density.

Initialization determines which posterior weights are calculated first, and those weights determine the first weighted fitting problem. Two different initial models can therefore follow different sequences of touching bounds. Inspecting the likelihood is a check of optimization progress within the chosen model; it does not choose the number of components or compare different covariance assumptions automatically.

With equal weights and a common isotropic variance, responsibilities are normalized exponentials of negative squared distances. Smaller variance makes them more concentrated around a uniquely closest mean. For every positive variance the posterior remains soft, whereas a hard nearest-centre update changes the optimization problem. This connection explains the resemblance between weighted mean updates and centre-based clustering without identifying their objectives.

Common mistakes

Common mistake

Replacing an E-step by a hard argmax

Starting from the updated means in Example 2, a hard assignment sends $0$ to the first component and $2$ to the second. Replacing each mean by its assigned observation returns the means to $(0,2)$. With the same fixed variance and equal weights, observed log-likelihood falls from about $-2.886563$ to $-2.970315$. The hard assignment is not the posterior distribution that makes the EM lower bound touch the current likelihood.

Common mistake

Using scatter as covariance or omitting mixture weights

The Gaussian covariance update divides weighted scatter by $N_k$. The E-step normalizes $\pi_k\phi_d(x_i;\mu_k,\Sigma_k)$ across components. Both factors arise from the mixture likelihood. Changing either normalization changes the fitting equations and breaks the stated EM derivation.

Summary

A Gaussian mixture sums weighted component densities. Its E-step uses normalized posterior responsibilities, and its M-step uses their effective masses to update weights, means, and any estimated covariances. The Jensen lower bound is exact at the current posterior and increases when the M-step increases $Q$, proving nondecreasing observed log-likelihood. Positive weights, positive-definite covariances, and positive effective masses keep the formulas within the stated model. The choice between estimated covariance and prescribed variance is part of the objective, and unrestricted covariance can yield unbounded likelihood through component collapse.

Exercises

Checkpoint

1. At $x=1$, two components have means $0$ and $2$, common variance $1$, and weights $(3/4,1/4)$. Find both responsibilities. What information is lost if only the largest responsibility is retained?

Separate equal component densities from unequal mixture weights.

Checkpoint

2. For observations $(0,2)$, suppose fixed responsibilities are $R=\begin{bmatrix}0.9&0.1\\0.3&0.7\end{bmatrix}$. With prescribed common variance, find the M-step weights and means.

Compute each $N_k$ before dividing a weighted sum.

Checkpoint

3. For observations $(0,2,4)$, one component has responsibilities $(4/5,1/2,1/5)$. Find its effective mass, updated mean, and one-dimensional covariance estimate when covariance is estimated.

The covariance is weighted squared deviation from the updated mean, divided by effective mass.

Checkpoint

4. Apply a hard assignment and mean update to the updated model in Example 2. Compare the resulting observed log-likelihood with its current value. Which equality in the EM proof no longer holds?

The midpoint between the current means is still $1$.

Solutions

Solution · 1. The posterior preserves both fractions

The responsibilities are $(3/4,1/4)$ because the common density cancels. Keeping only the first component records a single chosen label and discards the remaining posterior mass $1/4$.

Solution · 2. Fractional mass updates the mixture weights too

The masses are $(N_1,N_2)=(6/5,4/5)$. Thus $\pi^{\mathrm{new}}=(3/5,2/5)$ and $\mu_1^{\mathrm{new}}=(0.9\cdot0+0.3\cdot2)/(6/5)=1/2$, while $\mu_2^{\mathrm{new}}=(0.1\cdot0+0.7\cdot2)/(4/5)=7/4$.

Solution · 3. Normalize the weighted scatter by effective mass

The mass is $N=3/2$, and the mean is $((4/5)0+(1/2)2+(1/5)4)/(3/2)=6/5$. The weighted squared-deviation sum is $(4/5)(6/5)^2+(1/2)(4/5)^2+(1/5)(14/5)^2=76/25$. The covariance estimate is $(76/25)/(3/2)=152/75$.

Solution · 4. The lower bound no longer touches the current likelihood

The hard assignments are one observation per component, returning the means to $(0,2)$ and the weights to $(1/2,1/2)$. Observed log-likelihood decreases from $-2.886562696$ to $-2.970315405$. For the hard probability rows $q$, $\mathcal F(q,\theta^{\mathrm{old}})=\ell(\theta^{\mathrm{old}})$ no longer holds: those rows differ from the strictly positive posterior. Improving that hard-assignment bound cannot use the EM chain of inequalities to guarantee improvement of the observed likelihood.

Continue with least squares and residuals or the course reading list.