Gaussian Mixture Models

A Gaussian Mixture Model is a probability density defined as a convex combination of $K$ multivariate Gaussian components. Formally, for a data point $\mathbf{x} \in \mathbb{R}^d$, ``` p(x) = Σ_{k=1}^{K} π_k · N(x | μ_k, Σ_k) ``` where $\pi_k \geq 0$ are the mixing weights satisfying $\sum_k \pi_k = 1$, and $N(\mathbf{x}\mid\boldsymbol\mu_k,\boldsymbol\Sigma_k)$ is the multivariate Gaussian density with mean $\boldsymbol\mu_k$ and covariance $\boldsymbol\Sigma_k$. This is strictly more expressive than any single Gaussian: with enough components a GMM can approximate essentially any smooth density arbitrarily well, which is why it is a workhorse of nonparametric-flavored Density Estimation Fundamentals even though each component is parametric. The generative story introduces a latent variable $z \in \{1, \dots, K\}$, a categorical variable indicating which component "produced" a given point, with $P(z=k) = \pi_k$. Conditional on $z = k$, the point is drawn from $N(\boldsymbol\mu_k, \boldsymbol\Sigma_k)$. Marginalizing $z$ out recovers the mixture density above — this is why GMMs are called latent-variable models: the joint distribution over $(\mathbf{x}, z)$ is simple and Gaussian-conditional, but the observed marginal over $\mathbf{x}$ alone is a complicated multimodal density. The key derived quantity is the responsibility, the posterior probability that component $k$ generated a given observed point $\mathbf{x}_n$: ``` γ_nk = P(z_n = k | x_n) = [π_k N(x_n | μ_k, Σ_k)] / [Σ_{j=1}^{K} π_j N(x_n | μ_j, Σ_j)] ``` Responsibilities are the soft analogue of the hard cluster assignments in K-Means: instead of assigning $\mathbf{x}_n$ to a single nearest centroid, we distribute a unit of "belongingness" across all $K$ components proportionally to how well each explains the point. This soft assignment is what lets GMMs express genuine uncertainty about cluster membership near decision boundaries — exactly the ambiguity that hard clustering erases.

Fitting a GMM means choosing $\{\pi_k, \boldsymbol\mu_k, \boldsymbol\Sigma_k\}_{k=1}^K$ to maximize the log-likelihood of the observed data, ``` log p(X | θ) = Σ_{n=1}^{N} log [ Σ_{k=1}^{K} π_k N(x_n | μ_k, Σ_k) ] ``` The sum inside the logarithm (over the unobserved $z_n$) makes this optimization non-convex and analytically intractable to maximize directly — unlike a single Gaussian, there is no closed-form MLE. The standard solution is Expectation-Maximization (EM), which iterates two steps until convergence.

E-step (Expectation). Fix the current parameters and compute the responsibilities $\gamma_{nk}$ for every point and every component using Bayes' rule, as shown above. This step "fills in" our best guess about the missing labels $z_n$, in the form of a full posterior distribution rather than a hard guess.

M-step (Maximization). Fix the responsibilities and re-estimate the parameters by weighted maximum likelihood, treating $\gamma_{nk}$ as a soft count of how many points "belong" to cluster $k$: ``` N_k = Σ_n γ_nk (effective size of cluster k) π_k^new = N_k / N μ_k^new = (1 / N_k) Σ_n γ_nk · x_n Σ_k^new = (1 / N_k) Σ_n γ_nk · (x_n − μ_k^new)(x_n − μ_k^new)^T ``` Each of these updates has a clean closed form because, given the responsibilities, maximizing the expected complete-data log-likelihood decouples into $K$ independent weighted Gaussian fits. These two steps are repeated until the log-likelihood stops improving by more than a small tolerance, or a maximum number of iterations is reached. ``` function FIT_GMM(X, K, max_iter, tol): N, d = shape(X)

# Initialization (e.g., k-means++ or random responsibilities) π, μ, Σ = INITIALIZE(X, K) log_likelihood_old = -infinity for iteration in 1 .. max_iter:

# ---- E-step ---- for n in 1 .. N: for k in 1 .. K: w[n][k] = π[k] * GAUSSIAN_PDF(X[n], μ[k], Σ[k]) total = sum_k w[n][k] for k in 1 .. K: γ[n][k] = w[n][k] / total # responsibilities

# ---- M-step ---- for k in 1 .. K: N_k = sum_n γ[n][k] π[k] = N_k / N μ[k] = (sum_n γ[n][k] X[n]) / N_k diff = X - μ[k] # broadcast, shape (N, d) Σ[k] = (γ[:,k] diff)^T @ diff / N_k Σ[k] += ε * I # regularize: avoid singular covariance

# ---- Convergence check ---- log_likelihood_new = sum_n log( sum_k π[k] * GAUSSIAN_PDF(X[n], μ[k], Σ[k]) ) if |log_likelihood_new - log_likelihood_old| < tol: break log_likelihood_old = log_likelihood_new return π, μ, Σ, γ ``` Note two implementation-critical details baked into the pseudocode: the small ridge term `ε·I` added to each covariance prevents singular or ill-conditioned matrices (a component collapsing onto a single point drives $\Sigma_k \to 0$ and the likelihood to infinity — a well-known pathology), and the log-likelihood, not the raw parameters, is the correct quantity to monitor for convergence, since parameter changes can be small while likelihood is still improving substantially, or vice versa.

Time complexity. Each E-step evaluates the density of every point under every component: computing a single Gaussian log-density with a full $d \times d$ covariance costs $O(d^2)$ (dominated by the quadratic form $(\mathbf{x}-\boldsymbol\mu)^T \Sigma^{-1} (\mathbf{x}-\boldsymbol\mu)$ after a one-time $O(d^3)$ Cholesky factorization of $\Sigma_k$ per component per iteration). Summed over $N$ points and $K$ components, the E-step costs ``` O(N · K · d^2 + K · d^3) ``` The M-step recomputes weighted means ($O(NKd)$) and weighted covariances ($O(NKd^2)$), so it is dominated by the same $O(NKd^2)$ term. Over $T$ iterations until convergence, total cost is ``` O(T · (N K d^2 + K d^3)) ``` With diagonal covariances (a common practical restriction), the $d^2$ terms drop to $O(d)$, giving $O(TNKd)$ — the same order as K-Means' $O(TNKd)$, which makes sense since K-Means is a degenerate special case of a GMM (see below). Space complexity is $O(NK)$ for the responsibility matrix plus $O(Kd^2)$ for the covariances (or $O(Kd)$ if diagonal), plus $O(Nd)$ to store the data itself.

Correctness — why EM cannot decrease the likelihood. The proof rests on Jensen's inequality applied to a variational lower bound. For any distribution $q_n(z_n)$ over the latent label of point $n$, we can write ``` log p(x_n) = log Σ_k q_n(k) · [ π_k N(x_n|μ_k,Σ_k) / q_n(k) ] ≥ Σ_k q_n(k) · log [ π_k N(x_n|μ_k,Σ_k) / q_n(k) ] (Jensen, log concave) = L_n(q_n, θ) ``` Summing over $n$ gives a global lower bound $\mathcal{L}(q,\theta) \le \log p(X\mid\theta)$, with equality if and only if $q_n(k) = \gamma_{nk} = P(z_n = k \mid x_n, \theta)$ — exactly the E-step computation. So:

• The E-step sets $q_n = \gamma_n$, making the bound tight at the current $\theta$: $\mathcal{L}(q,\theta^{\text{old}}) = \log p(X\mid\theta^{\text{old}})$.
• The M-step maximizes $\mathcal{L}(q,\theta)$ over $\theta$ with $q$ held fixed — and this maximization has the closed form derived above, because with fixed responsibilities the expected complete-data log-likelihood decomposes into independent weighted Gaussian MLE problems, each solved exactly.

Chaining these two facts: ``` log p(X|θ^new) ≥ L(q, θ^new) ≥ L(q, θ^old) = log p(X|θ^old) ``` The first inequality is Jensen again (the bound never exceeds the true likelihood); the second holds because the M-step chose $\theta^{\text{new}}$ to maximize $\mathcal{L}$, so it cannot do worse than $\theta^{\text{old}}$. This chain proves the log-likelihood is monotonically non-decreasing across EM iterations, and since log-likelihood is bounded above (for non-degenerate covariance, ignoring the singular-covariance pathology), the sequence converges. This is the general argument underlying The EM Algorithm in any latent-variable model, not just Gaussian mixtures.

Honest caveats. EM converges to a local maximum (or saddle point) of the likelihood, not necessarily the global one — the objective is non-convex because of label-switching symmetry and multiple basins of attraction, so different initializations can converge to different, non-equivalent solutions; practitioners run multiple random restarts and keep the best. The likelihood is also unbounded above in the strict sense (a component's covariance can shrink around a single data point, driving likelihood to $+\infty$), so regularization (a minimum covariance eigenvalue, a Bayesian prior giving a MAP rather than MLE objective) is not optional in practice. Finally, choosing $K$ is a model-selection problem outside EM itself, typically handled with penalized criteria like BIC or AIC, or with a Bayesian nonparametric extension (a Dirichlet process mixture).

Consider a 1-D toy dataset (to make every number checkable by hand) with 6 points: ``` x = [1.0, 1.5, 2.0, 8.0, 8.5, 9.0] ``` Initialize $K=2$ components with $\pi_1=\pi_2=0.5$, $\mu_1 = 2$, $\mu_2 = 8$, and equal variances $\sigma_1^2 = \sigma_2^2 = 4$.

E-step (iteration 1). For each point compute the unnormalized weight $\pi_k \cdot N(x\mid\mu_k,\sigma_k^2)$ under each component, then normalize. Using $\sigma^2=4$ for both, points near 2 get overwhelmingly assigned to component 1 and points near 8 to component 2:

(The tiny cross-component values, e.g. $N(8\mid 2,4) \approx 4\times10^{-5}$, are negligible next to $N(8\mid8,4)\approx0.1995$, which is why the table rounds them to 0.000 — this well-separated example makes responsibilities nearly hard, which is instructive: it shows GMM soft assignments reduce to K-Means-like hard assignments when clusters are cleanly separated, and only differ meaningfully near overlaps.)

M-step (iteration 1). With responsibilities essentially binary here, the updates match a weighted average within each hardened group: ``` N_1 = 3, π_1 = 3/6 = 0.5 μ_1 = (1.0 + 1.5 + 2.0)/3 = 1.5 σ_1^2 = [(1.0-1.5)^2 + (1.5-1.5)^2 + (2.0-1.5)^2] / 3 = 0.5/3 ≈ 0.1667 N_2 = 3, π_2 = 3/6 = 0.5 μ_2 = (8.0 + 8.5 + 9.0)/3 = 8.5 σ_2^2 = [(8.0-8.5)^2 + (8.5-8.5)^2 + (9.0-8.5)^2] / 3 ≈ 0.1667 ```

Convergence. Because the two groups are so well separated, this problem essentially converges after this single M-step: recomputing responsibilities with the new, tighter Gaussians ($\mu_1=1.5,\sigma_1^2\approx0.167$ and $\mu_2=8.5,\sigma_2^2\approx0.167$) still assigns each point with responsibility indistinguishable from 1 or 0 to its natural group, so subsequent E/M steps leave the parameters essentially unchanged and the log-likelihood plateaus. In a real dataset with overlapping clusters, this convergence would take many more iterations and responsibilities would settle at genuinely fractional values (e.g., $\gamma_{n1}=0.7,\gamma_{n2}=0.3$) for points near the boundary between components — the exact behavior that distinguishes GMM soft clustering from hard K-Means partitions.

The mixture-of-Gaussians idea and the general EM framework that fits it emerged from a lineage connecting 19th-century statistics to modern machine learning.

Gaussian Mixture Models sit at the intersection of density estimation and clustering, and this dual identity is why they appear across so many production systems. In speech and speaker recognition, GMMs historically modeled the distribution of acoustic feature vectors (before being supplanted, but not entirely replaced, by neural acoustic models); GMM-based Universal Background Models are still used in speaker verification pipelines. In computer vision, per-pixel GMMs are a classic technique for adaptive background subtraction in surveillance video, modeling each pixel's color distribution as a small mixture that adapts over time. In genomics and bioinformatics, GMMs cluster gene expression profiles or flow-cytometry measurements where cell populations naturally form multiple modes. GMMs are also a standard baseline and building block for Anomaly Detection Fundamentals: once a GMM is fit to "normal" data, a new point's mixture density $p(\mathbf{x})$ directly gives an anomaly score — points in low-density regions between or outside components are flagged as outliers, with a natural probabilistic threshold rather than an arbitrary distance cutoff. Finance uses this for fraud and risk scoring, and manufacturing uses it for defect detection from sensor telemetry. Finally, GMMs are the pedagogical entry point into the broader The EM Algorithm framework, whose descendants — hidden Markov models, latent Dirichlet allocation, mixture-of-experts networks, and Gaussian mixture priors inside variational autoencoders — power speech recognition, topic modeling, and modern generative models alike.

Ask the tutor: Can you show how the EM updates change if we constrain all components to share a single covariance matrix?