Gaussian Mixture Models and the EM Algorithm

k-means gives every point a hard label and assumes round clusters. But real data often has overlapping, elliptical groups, and we may want to say "this point is 70% likely from cluster A and 30% from cluster B". Gaussian Mixture Models (GMMs) provide exactly this, and fitting them introduces one of the most important algorithms in statistics and ML: Expectation–Maximisation (EM).

The model#

A GMM assumes each point is generated by first choosing a component $z \in \{1, \dots, K\}$ with probability $\pi_k$ (the mixing weights), then drawing $\mathbf{x}$ from that component's Gaussian:

p(\mathbf{x}) = \sum_{k=1}^{K}\pi_k\,\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k), \qquad \sum_k\pi_k = 1

The component label $z$ is a latent (hidden) variable. A GMM is both a clustering method and a flexible density estimator — with enough components, it can approximate any smooth density.

Why is fitting hard?#

The log-likelihood is

\ell(\boldsymbol{\theta}) = \sum_{i=1}^{n}\ln\sum_{k=1}^{K}\pi_k\,\mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)

The log of a sum does not decompose, so there is no closed-form maximiser. But if we knew each point's component, the problem would be trivial: fit each Gaussian to its own points. EM exploits this.

The EM algorithm#

E-step (Expectation). Given current parameters, compute the responsibility of component $k$ for point $i$ — the posterior probability that $i$ came from $k$:

\gamma_{ik} = \frac{\pi_k\,\mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_{j}\pi_j\,\mathcal{N}(\mathbf{x}_i \mid \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)}

M-step (Maximisation). Re-estimate parameters as responsibility-weighted statistics, with $N_k = \sum_i\gamma_{ik}$:

\pi_k = \frac{N_k}{n}, \qquad \boldsymbol{\mu}_k = \frac{1}{N_k}\sum_i\gamma_{ik}\mathbf{x}_i, \qquad \boldsymbol{\Sigma}_k = \frac{1}{N_k}\sum_i\gamma_{ik}(\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^\top

Repeat until the log-likelihood stops improving.

Why EM works#

For any distribution $q(z)$ over the latent variables, Jensen's inequality gives a lower bound on the log-likelihood — the Evidence Lower BOund (ELBO):

\ln p(\mathbf{x} \mid \boldsymbol{\theta}) \ge \sum_z q(z)\ln\frac{p(\mathbf{x}, z \mid \boldsymbol{\theta})}{q(z)} = \mathcal{L}(q, \boldsymbol{\theta})

The gap equals $D_{\text{KL}}(q(z)\,\|\,p(z \mid \mathbf{x}, \boldsymbol{\theta}))$.

The E-step sets $q$ to the exact posterior, closing the gap: the bound touches the likelihood.
The M-step maximises the bound over $\boldsymbol{\theta}$, which can only raise the likelihood.

Hence EM never decreases the likelihood. It converges to a local maximum (or saddle point), so initialisation matters — typically with k-means and several restarts. This ELBO derivation is also the foundation of variational autoencoders, where the E-step is replaced by a learned neural network.

Implementation from scratch#

python

import numpy as np
from scipy.stats import multivariate_normal

def gmm_em(X, K, iters=200, seed=0, tol=1e-6):
    n, d = X.shape
    rng = np.random.default_rng(seed)
    mu = X[rng.choice(n, K, replace=False)]
    Sigma = np.array([np.cov(X.T) + 1e-6 * np.eye(d) for _ in range(K)])
    pi = np.full(K, 1 / K)
    prev = -np.inf
    for _ in range(iters):
        # E-step (in log space for stability)
        logp = np.column_stack([np.log(pi[k]) + multivariate_normal.logpdf(X, mu[k], Sigma[k])
                                for k in range(K)])
        ll = np.logaddexp.reduce(logp, axis=1)
        gamma = np.exp(logp - ll[:, None])
        # M-step
        Nk = gamma.sum(0)
        pi = Nk / n
        mu = (gamma.T @ X) / Nk[:, None]
        for k in range(K):
            D = X - mu[k]
            Sigma[k] = (gamma[:, k, None] * D).T @ D / Nk[k] + 1e-6 * np.eye(d)
        if ll.sum() - prev < tol:
            break
        prev = ll.sum()
    return pi, mu, Sigma, gamma, ll.sum()

The small ridge 1e-6 * I prevents a degenerate solution: if a component collapses onto a single point, its variance shrinks to zero and the likelihood goes to infinity. This singularity is a real hazard of maximum-likelihood GMMs.

Covariance structures#

`covariance_type`	Parameters per component	Cluster shape
`spherical`	1 variance	Circles
`diag`	$d$ variances	Axis-aligned ellipses
`tied`	one shared full matrix	Same ellipse for all
`full`	$d(d+1)/2$	Any orientation and shape

Choosing the number of components#

Likelihood always increases with $K$, so use penalised criteria:

\text{BIC} = -2\ell(\hat{\boldsymbol{\theta}}) + p\ln n, \qquad \text{AIC} = -2\ell(\hat{\boldsymbol{\theta}}) + 2p

where $p$ is the number of free parameters. Choose the model with the lowest BIC. A Bayesian alternative, the variational Dirichlet-process GMM (BayesianGaussianMixture), switches off unnecessary components automatically.

python

from sklearn.mixture import GaussianMixture
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=800, centers=3, cluster_std=[1.0, 2.0, 0.6], random_state=2)
X = X @ np.array([[0.6, -0.6], [-0.4, 0.8]])            # stretch into ellipses
for K in range(1, 7):
    g = GaussianMixture(K, covariance_type="full", n_init=5, random_state=0).fit(X)
    print(f"K={K}: BIC={g.bic(X):9.1f}")

Uses of GMMs#

Soft clustering with probabilistic memberships.
Density estimation and anomaly detection (flag points with low $p(\mathbf{x})$).
Speaker identification — GMMs over acoustic features were the standard before deep learning.
Background subtraction in video.
Generating synthetic data by sampling from the fitted mixture.

#gmm#em algorithm#clustering#latent variables#density estimation

Written by

Janin A Apurba

B.Sc. in CSE, AUST · Advanced ICT Officer, CNRS-UNHCR. Teaching AI, ML and Deep Learning to the next generation of engineers and researchers.