Gaussian Mixture Models

Info

Gaussian Mixture Models (GMMs) are a probabilistic model that assumes all data points are generated from a mixture of several Gaussian distributions with unknown parameters. They are widely used for clustering and density estimation.

Mixture Density

A GMM models the data density as a weighted sum of $K$ Gaussian components:

$$p(x)=\sum _{k=1}^K{\pi }_k\mathcal{N}(x\mid {\mu }_k,{\Sigma }_k),\sum _{k=1}^K{\pi }_k=1,{\pi }_k\ge 0.$$

The parameters are the mixing weights ${\pi }_k$, means ${\mu }_k$, and covariances ${\Sigma }_k$. Unlike a single Gaussian, a mixture can represent multimodal, skewed, and elongated cluster shapes.

Latent Assignments

Introduce a latent one-hot variable $z\in \{1,\dots ,K\}$ naming which component generated $x$, with $p(z=k)={\pi }_k$. The responsibility is the posterior probability that component $k$ produced point $x_i$:

$${\gamma }_{ik}=p(z_i=k\mid x_i)=\frac{{\pi }_k\mathcal{N}(x_i\mid {\mu }_k,{\Sigma }_k)}{{\sum }_{j=1}^K{\pi }_j\mathcal{N}(x_i\mid {\mu }_j,{\Sigma }_j)}.$$

This is a soft assignment: each point is fractionally owned by every component.

EM for GMM

Because the latent $z$ makes the log-likelihood non-convex, we optimize with Expectation-Maximization, alternating two steps that never decrease the likelihood.

E step: given current parameters, compute responsibilities ${\gamma }_{ik}$ (formula above).

M step: given responsibilities, update parameters by weighted maximum likelihood. With $N_k={\sum }_i{\gamma }_{ik}$:

$${\mu }_k=\frac{1}{N_k}\sum _i{\gamma }_{ik}{x}_i,{\Sigma }_k=\frac{1}{N_k}\sum _i{\gamma }_{ik}(x_i-{\mu }_k)(x_i-{\mu }_k{)}^{\top },{\pi }_k=\frac{N_k}{N}.$$

Relation to K-Means Clustering

K-means is the limiting case of GMM where covariances are ${\sigma }^{2}I$ with $\sigma \to 0$ and assignments become hard (each point fully owned by its nearest center). GMM generalizes it with soft assignments and per-cluster shape via full covariances, so it handles elliptical, differently-sized, overlapping clusters that K-means cannot.

numpy code

import numpy as np
from scipy.stats import multivariate_normal
 
def em_gmm(X, K, iters=50):
    n, d = X.shape
    pi = np.ones(K) / K
    mu = X[np.random.choice(n, K, replace=False)]
    Sigma = np.stack([np.cov(X.T) for _ in range(K)])
 
    for _ in range(iters):
        # E step: responsibilities (n, K)
        r = np.stack([pi[k] * multivariate_normal.pdf(X, mu[k], Sigma[k])
                      for k in range(K)], axis=1)
        r /= r.sum(1, keepdims=True)
 
        # M step
        Nk = r.sum(0)
        pi = Nk / n
        mu = (r.T @ X) / Nk[:, None]
        for k in range(K):
            diff = X - mu[k]
            Sigma[k] = (r[:, k, None] * diff).T @ diff / Nk[k]
            Sigma[k] += 1e-6 * np.eye(d)   # regularize
    return pi, mu, Sigma, r

Pitfalls

• A component can collapse onto a single point, sending its covariance to zero and the likelihood to infinity; add a small ridge to ${\Sigma }_k$ as above.
• EM only finds a local optimum; use k-means++ initialization and several restarts.
• Choose $K$ with BIC / AIC rather than raw likelihood, which always improves with more components.

Related

• Trained via Expectation-Maximization.
• Soft-assignment generalization of K-Means Clustering.