Gaussian Mixture Models and EM algorithm

Overview

Gaussian Mixture Models (GMM) are probabilistic models that assume data comes from a mixture of several Gaussian distributions, each representing a cluster. The Expectation-Maximization (EM) algorithm is an iterative method to find the parameters of these Gaussians when we don't know which data point belongs to which cluster.

Think of it like genres: a movie isn't purely "action" or "drama" — it's 60% action, 40% drama. GMMs capture this blend naturally.

The Mathematical Foundation

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

where:

• $K$ = number of mixture components (clusters)
• $\pi_k$ = mixing coefficient for component $k$ (satisfies $\sum_{k=1}^{K} \pi_k = 1$, $\pi_k \geq 0$)
• $\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$ = Gaussian distribution with mean $\boldsymbol{\mu}_k$ and covariance $\boldsymbol{\Sigma}_k$

WHY this form? We're saying: "Pick cluster $k$ with probability $\pi_k$, then generate a point from that cluster's Gaussian." The sum marginalizes over all possible cluster choices.

Deriving the Gaussian Component

The multivariate Gaussian distribution is:

$\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Sigma}) = \frac{1}{(2\pi)^{d/2} |\boldsymbol{\Sigma}|^{1/2}} \exp\left(-\frac{1}{2}(\mathbf{x} - \boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})\right)$

WHY this formula?

• The exponent $-\frac{1}{2}(\mathbf{x} - \boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})$ is the Mahalanobis distance squared — it measures how far $\mathbf{x}$ is from $\boldsymbol{\mu}$ accounting for correlations (captured by $\boldsymbol{\Sigma}$)
• The normalizing constant ensures $\int p(\mathbf{x}) d\mathbf{x} = 1$
• $(2\pi)^{d/2}$ comes from integrating the Gaussian over $d$ dimensions
• $|\boldsymbol{\Sigma}|^{1/2}$ accounts for the "volume" of the covariance ellipsoid

The EM Algorithm: From First Principles

THE PROBLEM: We have data $\mathbf{X} = \{\mathbf{x}_1, \ldots, \mathbf{x}_N\}$ but we don't know which cluster generated each point. We need to estimate $\boldsymbol{\theta} = \{\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k\}_{k=1}^{K}$.

WHY is this hard? If we knew cluster assignments, estimating parameters would be easy (just compute means and covariances). If we knew parameters, computing assignments would be easy (just use Bayes' rule). We know neither! This is a chicken-and-egg problem.

Introducing Latent Variables

Let $z_i \in \{1, \ldots, K\}$ be the hidden cluster label for $\mathbf{x}_i$. The complete-data log-likelihood would be:

$\log p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\theta}) = \sum_{i=1}^{N} \sum_{k=1}^{K} \mathbb{1}(z_i = k) \left[\log \pi_k + \log \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\right]$

WHY this form? The indicator $\mathbb{1}(z_i = k)$ picks out which cluster $i$ belongs to, and we log the probability of that choice ($\pi_k$) times the probability of $\mathbf{x}_i$ under that Gaussian.

But we don't observe $\mathbf{Z}$! So we maximize the incomplete-data log-likelihood:

$\log p(\mathbf{X} | \boldsymbol{\theta}) = \sum_{i=1}^{N} \log \left(\sum_{k=1}^{K} \pi_k \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\right)$

This is hard to maximize directly because the log of a sum has no closed form.

By Jensen's inequality, for any distribution $q(\mathbf{Z})$:

$\log p(\mathbf{X} | \boldsymbol{\theta}) = \log \sum_{\mathbf{Z}} p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\theta})$ $= \log \sum_{\mathbf{Z}} q(\mathbf{Z}) \frac{p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\theta})}{q(\mathbf{Z})}$ $\geq \sum_{\mathbf{Z}} q(\mathbf{Z}) \log \frac{p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\theta})}{q(\mathbf{Z})} = \mathcal{L}(q, \boldsymbol{\theta})$

WHY does this help? The lower bound $\mathcal{L}(q, \boldsymbol{\theta})$ is easier to maximize. We alternate:

1. E-step: Fix $\boldsymbol{\theta}$, choose $q$ to make the bound tight: $q(\mathbf{Z}) = p(\mathbf{Z} | \mathbf{X}, \boldsymbol{\theta}^{\text{old}})$
2. M-step: Fix $q$, maximize $\mathcal{L}$ with respect to $\boldsymbol{\theta}$ to get $\boldsymbol{\theta}^{\text{new}}$

E-Step: Computing Responsibilities

We need $p(z_i = k | \mathbf{x}_i, \boldsymbol{\theta})$, called the responsibility $\gamma_{ik}$.

By Bayes' rule:

$\gamma_{ik} = p(z_i = k | \mathbf{x}_i, \boldsymbol{\theta}) = \frac{p(z_i = k) p(\mathbf{x}_i | z_i = k, \boldsymbol{\theta})}{\sum_{j=1}^{K} p(z_i = j) p(\mathbf{x}_i | z_i = j, \boldsymbol{\theta})}$

WHY this step? Numerator: prior probability of cluster $k$ times likelihood of $\mathbf{x}_i$ under cluster $k$. Denominator: total probability over all clusters (normalization).

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

WHAT does $\gamma_{ik}$ mean? The probability that point $i$ was generated by cluster $k$. If $\gamma_{i3} = 0.8$, point $i$ is 80% likely from cluster 3.

M-Step: Updating Parameters

We maximize the expected complete-data log-likelihood:

$Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text{old}}) = \mathbb{E}_{\mathbf{Z} | \mathbf{X}, \boldsymbol{\theta}^{\text{old}}}[\log p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\theta})]$

Substituting responsibilities:

$Q = \sum_{i=1}^{N} \sum_{k=1}^{K} \gamma_{ik} \left[\log \pi_k + \log \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\right]$

Deriving $\boldsymbol{\mu}_k$ update:

$\frac{\partial Q}{\partial \boldsymbol{\mu}_k} = \sum_{i=1}^{N} \gamma_{ik} \frac{\partial}{\partial \boldsymbol{\mu}_k} \left[-\frac{1}{2}(\mathbf{x}_i - \boldsymbol{\mu}_k)^T \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k)\right]$

WHY this derivative? Only the quadratic term in the Gaussian depends on $\boldsymbol{\mu}_k$.

$= \sum_{i=1}^{N} \gamma_{ik} \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k)$

Setting to zero:

$\sum_{i=1}^{N} \gamma_{ik} \boldsymbol{\Sigma}_k^{-1} \mathbf{x}_i = \sum_{i=1}^{N} \gamma_{ik} \boldsymbol{\Sigma}_k^{-1} \boldsymbol{\mu}_k$

$\boldsymbol{\mu}_k = \frac{\sum_{i=1}^{N} \gamma_{ik} \mathbf{x}_i}{\sum_{i=1}^{N} \gamma_{ik}} = \frac{1}{N_k} \sum_{i=1}^{N} \gamma_{ik} \mathbf{x}_i$

where $N_k = \sum_{i=1}^{N} \gamma_{ik}$ is the effective number of points in cluster $k$.

WHAT does this mean? The new mean is the weighted average of all points, weighted by their responsibility to cluster $k$. Points strongly in cluster $k$ ($\gamma_{ik} \approx 1$) pull the mean toward them; points weakly in $k$ have little effect.

Deriving $\boldsymbol{\Sigma}_k$ update:

Similarly (taking derivative of log-determinant and quadratic terms):

$\boldsymbol{\Sigma}_k = \frac{1}{N_k} \sum_{i=1}^{N} \gamma_{ik} (\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^T$

WHY? Weighted average of squared deviations from the mean — the definition of covariance, responsibility-weighted.

Deriving $\pi_k$ update:

We maximize $\sum_{k=1}^{K} N_k \log \pi_k$ subject to $\sum_{k=1}^{K} \pi_k = 1$ using a Lagrange multiplier $\lambda$:

$\mathcal{L} = \sum_{k=1}^{K} N_k \log \pi_k + \lambda\left(1 - \sum_{k=1}^{K} \pi_k\right)$

$\frac{\partial \mathcal{L}}{\partial \pi_k} = \frac{N_k}{\pi_k} - \lambda = 0 \implies \pi_k = \frac{N_k}{\lambda}$

Suming over $k$: $\sum_{k=1}^{K} \pi_k = \frac{\sum_{k=1}^{K} N_k}{\lambda} = \frac{N}{\lambda} = 1 \implies \lambda = N$

$\pi_k = \frac{N_k}{N} = \frac{1}{N}\sum_{i=1}^{N} \gamma_{ik}$

WHAT does this mean? The mixing coefficient is the fraction of all responsibility assigned to cluster $k$.

Initialize: Choose $K$, initialize $\pi_k$, $\boldsymbol{\mu}_k$, $\boldsymbol{\Sigma}_k$ randomly or via K-means

Repeat until convergence:

E-step: Compute responsibilities for all $i, k$: $\gamma_{ik} = \frac{\pi_k \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_{j=1}^{K} \pi_j \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)}$

M-step: Update parameters for all $k$: $N_k = \sum_{i=1}^{N} \gamma_{ik}$ $\boldsymbol{\mu}_k = \frac{1}{N_k} \sum_{i=1}^{N} \gamma_{ik} \mathbf{x}_i$ $\boldsymbol{\Sigma}_k = \frac{1}{N_k} \sum_{i=1}^{N} \gamma_{ik} (\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^T$ $\pi_k = \frac{N_k}{N}$

Convergence: Stop when log-likelihood change < threshold or max iterations reached.

Worked Examples

Data: $\mathbf{X} = \{1, 2, 9, 10\}$, $K=2$

Initialize:

• $\pi_1 = 0.5, \pi_2 = 0.5$
• $\mu_1 = 1.5, \mu_2 = 9.5$ (from data endpoints)
• $\sigma_1^2 = 1, \sigma_2^2 = 1$ (arbitrary small variance)

E-step iteration 1:

For $x_1 = 1$: $\mathcal{N}(1 | 1.5, 1) = \frac{1}{\sqrt{2\pi}} e^{-0.125} \approx 0.352$ $\mathcal{N}(1 | 9.5, 1) = \frac{1}{\sqrt{2\pi}} e^{-36.125} \approx 0$

WHY these values? Plug into the Gaussian formula. Point 1 is very close to $\mu_1 = 1.5$ (distance 0.5) but far from $\mu_2 = 9.5$ (distance 8.5).

$\gamma_{11} = \frac{0.5 \times 0.352}{0.5 \times 0.352 + 0.5 \times 0} \approx 1.0$

WHAT does this mean? Point 1 is essentially100% assigned to cluster 1.

Similarly:

• $\gamma_{21} \approx 1.0$ (point 2 → cluster 1)
• $\gamma_{31} \approx 0.0$ (point 9 → cluster 2)
• $\gamma_{41} \approx 0.0$ (point 10 → cluster 2)

M-step iteration 1:

$N_1 = \gamma_{11} + \gamma_{21} + \gamma_{31} + \gamma_{41} \approx 2.0$ $N_2 \approx 2.0$

WHY? Two points strongly in cluster 1, two in cluster 2.

$\mu_1 = \frac{1 \times 1 + 1 \times 2 + 0 \times 9 + 0 \times 10}{2} = \frac{3}{2} = 1.5$

WHY this step? Weighted average of all points, but only points 1 and 2 have non-zero weight for cluster 1.

$\mu_2 = \frac{0 \times 1 + 0 \times 2 + 1 \times 9 + 1 \times 10}{2} = 9.5$

$\sigma_1^2 = \frac{1 \times (1-1.5)^2 + 1 \times (2-1.5)^2}{2} = \frac{0.25 + 0.25}{2} = 0$

WHY? Variance is the weighted average of squared deviations. Both points are0.5 away from mean 1.5.

$\sigma_2^2 = \frac{1 \times (9-9.5)^2 + 1 \times (10-9.5)^2}{2} = 0.25$

$\pi_1 = \pi_2 = \frac{2}{4} = 0.5$

Result: After one iteration, the clusters have tightened ($\sigma^2$ decreased from 1 to 0.25), means are stable. The algorithm would converge quickly since the data is well-separated.

Data:

• Cluster 1: $(0,0), (1,1), (0.5, 0.5)$
• Cluster 2: $(4,4), (5,5), (4.5, 4.5)$
• Ambiguous point: $(2.5, 2.5)$

After convergence:

• $\boldsymbol{\mu}_1 \approx (0.5, 0.5)$
• $\boldsymbol{\mu}_2 \approx (4.5, 4.5)$
• For the ambiguous point $(2.5, 2.5)$: $\gamma_{71} \approx 0.5, \gamma_{72} \approx 0.5$

WHY 50/50? The point is equidistant from both cluster centers, so it has equal responsibility to both. This is the key advantage of GMM: soft assignments capture uncertainty.

If we forced hard assignment (like K-means), we'd arbitrarily put it in one cluster, losing information. GMM says "this point is genuinely ambiguous — it's half from each cluster."

For convergence checking, we compute:

$\log p(\mathbf{X} | \boldsymbol{\theta}) = \sum_{i=1}^{N} \log \left(\sum_{k=1}^{K} \pi_k \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\right)$

For point $\mathbf{x}_1 = (1, 1)$ with $K=2$:

$p(\mathbf{x}_1) = \pi_1 \mathcal{N}(\mathbf{x}_1 | \boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1) + \pi_2 \mathcal{N}(\mathbf{x}_1 | \boldsymbol{\mu}_2, \boldsymbol{\Sigma}_2)$

If $\pi_1=0.6$, $\pi_2=0.4$, and the Gaussians evaluate to $0.2$ and $0.001:

$p(\mathbf{x}_1) = 0.6 \times 0.2 + 0.4 \times 0.001 = 0.12 + 0.004 = 0.1204$

WHY sum? Total probability is the mixture of both components.

$\log p(\mathbf{x}_1) = \log(0.1204) \approx -2.116$

Sum over all points to get total log-likelihood. If it changes by less than $10^{-6}$ between iterations, we've converged.

Common Mistakes and Fixes

Wrong: Compute $\gamma_{ik} = \pi_k \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)$ and use directly.

Why it feels right: We computed the probability of $\mathbf{x}_i$ under cluster $k$, that should be the responsibility!

Why it's wrong: Without the denominator, $\sum_{k=1}^{K} \gamma_{ik} \neq 1$. Responsibilities are conditional probabilities and must sum to 1 for each point.

Fix: Always include the normalizing sum: $\gamma_{ik} = \frac{\pi_k \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)}{\sum_{j=1}^{K} \pi_j \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j)}$

Wrong: A cluster colapses to a single point, making $|\boldsymbol{\Sigma}_k| = 0$, which breaks the Gaussian evaluation (divide by zero).

Why it happens: If all points in a cluster have identical coordinates (or nearly identical due to numerical issues), the covariance becomes singular.

Why it's dangerous: The algorithm crashes or produces NaN values.

Fix: Add regularization by adding a small diagonal term: $\boldsymbol{\Sigma}_k = \frac{1}{N_k} \sum_{i=1}^{N} \gamma_{ik} (\mathbf{x}_i - \boldsymbol{\mu}_k)(\mathbf{x}_i - \boldsymbol{\mu}_k)^T + \epsilon\mathbf{I}$ where $\epsilon = 10^{-6}$ is a small constant. This ensures $\boldsymbol{\Sigma}_k$ is always positive definite.

Steel-man: "I don't want to add artificial variance — it's not in the data!" True, but with zero variance, the model is saying "this cluster generates only one exact point with infinite precision," which is physically unrealistic. The regularization represents measurement noise or model uncertainty.

Wrong: Initialize all $\boldsymbol{\mu}_k$ randomly from a uniform distribution without looking at the data.

Why it feels right: Random initialization works for gradient descent on neural networks, so it should work here too!

Why it's wrong: EM is sensitive to initialization. If two means start very close together, they may stay close (both trying to model the same cluster) while another cluster is missed entirely.

Fix: Use K-means++ for initialization:

1. Run K-means to convergence (fast, simple)
2. Use K-means cluster centers as $\boldsymbol{\mu}_k$
3. Use K-means cluster covariances as $\boldsymbol{\Sigma}_k$
4. Use K-means cluster sizes to initialize $\pi_k$

Why this works: K-means gives a good rough clustering, and EM refines it with soft assignments.

Wrong: Run EM once, get result, assume it's the best possible clustering.

Why it feels right: EM is guaranteed to increase the likelihood each iteration, so it must converge to the global maximum!

Why it's wrong: EM is guaranteed to converge to a local maximum, not necessarily the global one. Different initializations lead to different local maxima.

Fix: Run EM multiple times (5-10) with different random initializations, pick the result with highest final log-likelihood.

Steel-man: "But that's computationally expensive!" Yes, but it's necessary. The log-likelihood landscape is non-convex with many local maxima. Think of it like training neural networks: you try different random seeds and pick the best.

Choosing the Number of Clusters K

Unlike K-means, GMM is probabilistic, so we can use model selection criteria:

Akaike Information Criterion (AIC): $\text{AIC} = 2k - 2\ln(\hat{L})$

where $k$ is the number of parameters, $\hat{L}$ is the maximum likelihood.

Bayesian Information Criterion (BIC): $\text{BIC} = k \ln(N) - 2\ln(\hat{L})$

WHY BIC over AIC? BIC penalizes model complexity more heavily (factor of $\ln(N)$ vs. 2). For GM, BIC often performs better at avoiding overfitting.

HOW to use: Fit GMMs with $K = 1, 2, 3, \ldots$ and plot BIC. Choose $K$ with the lowest BIC.

Parameter count for GMM:

• $K-1$ mixing coefficients (constraint: sum to 1)
• $K \times d$ means
• $K \times d(d+1)/2$ covariance parameters (symmetric matrices)

Total: $k = (K-1) + Kd + K \frac{d(d+1)}{2}$

Connections

• K-Means Clustering: GM is a soft version of K-means; K-means is GM with spherical covariances and hard assignments
• Maximum Likelihood Estimation: EM maximizes the likelihood of observed data under the mixture model
• Bayes Theorem: Responsibilities are computed via Bayes' rule
• Expectation Maximization Algorithm: GM is the most common application of the general EM framework
• Principal Component Analysis: GMM covariances capture data spread like PCA, but per-cluster
• Hidden Markov Models: HMMs useEM with similar E-step (forward-backward) and M-step updates
• Variational Inference: Modern Bayesian GMs use variational EM (more robust, automatic $K$ selection)
• Anomaly Detection: Points with low probability under all mixture components are outliers
• Density Estimation: GMs model the full probability density, useful for generative modeling

Or think: Evaluate, then Modify. You evaluate where points belong, then modify the cluster parameters accordingly.

You start by guessing: "I think the smaller marbles were red and the bigger ones were blue."

Then you look at each marble and say: "This is medium-sized... it's probably 60% likely to be red and 40% likely to be blue." You're not sure, so you hedge your bets.

Now you update your guess: "Okay, if I treat this marble as 60% red, then the average size of red marbles is probably this much..." You recalculate the average size of reds and blues based on these partial assignments.

Then you go back and look at each marble again with your new averages: "Hmm, with this new average, maybe this marble is actually 70% red, not 60%."

You keep going back and forth — updating your guess about which marbles are which color, then updating your guess about what the average red and blue sizes are — until your guesses stop changing. That's the EM algorithm! You're alternating between estimating assignments (E-step) and updating parameters (M-step).

GMs do this with data points and clusters instead of marbles and colors.

Flashcards

#flashcards/ai-ml

What does GMM stand for and what is its key advantage over K-means? :: Gaussian Mixture Model. Key advantage: soft assignments — each point has a probability of belonging to each cluster, capturing uncertainty, rather than hard assignment to one cluster.

What

Concept Map

Hinglish (regional understanding)