From K-means to GMM

K-means to GMM: From Distance Minimization to Probabilistic Models and the EM Derivation

Goal: starting from the geometric distance in K-means, derive step by step the Gaussian Mixture Model (GMM) and its EM algorithm, stringing everything together with detailed math:
distance → exponential-form likelihood → Gaussian distribution → mixture model → latent variables → joint distribution and marginalization → E-step posterior → M-step closed-form solution → EM monotonicity → K-means as a limiting case.

---

1. K-means: Distance-Based Hard Clustering

1.1 Model Setup

Dataset:

$$X = \{x_1,\dots,x_N\},\quad x_n\in\mathbb{R}^D.$$

Cluster centers (“prototypes”):

$$\{\mu_k\}_{k=1}^K,\quad \mu_k\in\mathbb{R}^D.$$

Hard assignment variables (the implicit cluster assignment):

$$r_{nk}\in\{0,1\},\quad \sum_{k=1}^K r_{nk} = 1,\quad \forall n.$$

1.2 The K-means Objective

K-means minimizes the within-cluster sum of squared distances:

$$J_{\text{kmeans}}(\{\mu_k\},\{r_{nk}\}) = \sum_{n=1}^N \sum_{k=1}^K r_{nk}\|x_n - \mu_k\|^2.$$

The algorithm alternates between two steps:

• Assignment step (analogous to the E-step): given the centers , assign each sample to its nearest cluster:

  $$r_{nk} = \begin{cases} 1, & k = \displaystyle\arg\min_j \|x_n - \mu_j\|^2,\\[4pt] 0, & \text{otherwise}. \end{cases}$$

• Update step (analogous to the M-step): given the assignments , update each center to the mean of the samples in its cluster:

  $$\mu_k = \frac{\sum_{n=1}^N r_{nk} x_n}{\sum_{n=1}^N r_{nk}}.$$

At this point K-means is purely a geometric optimization problem, with no explicit probabilistic structure yet.

---

2. From “Distance” to “Probability”: Exponential Form and the Gaussian Distribution

We want to formalize “closer distance → higher probability” as a probabilistic model.

2.1 The Energy View: Distance as “Energy”

Intuitive assumption:

The closer a point is to the center , the more likely it “belongs to this cluster”.

Construct an unnormalized likelihood/density:

$$\tilde{p}(x|\mu) = \exp\big(-\lambda \|x - \mu\|^2\big),\quad \lambda>0.$$

• Non-negative;
• Large distance ⇒ exponential decay;
• But it does not yet “integrate to 1”.

2.2 Normalization: The 1-D Case

Start with one dimension:

$$\tilde{p}(x|\mu) = \exp\big(-\lambda (x-\mu)^2\big).$$

Normalization constant:

$$Z = \int_{-\infty}^{\infty} \exp\big(-\lambda (x-\mu)^2\big)\,dx.$$

Shift the variable, :

$$Z = \int_{-\infty}^{\infty} \exp(-\lambda y^2)\,dy.$$

Use the Gaussian integral (a classic result):

$$\int_{-\infty}^{\infty} e^{-\lambda y^2} dy = \sqrt{\frac{\pi}{\lambda}}.$$

Therefore:

$$Z = \sqrt{\frac{\pi}{\lambda}}, \quad p(x|\mu) = \frac{1}{Z} e^{-\lambda (x-\mu)^2} = \sqrt{\frac{\lambda}{\pi}}\, e^{-\lambda (x-\mu)^2}.$$

Setting gives the standard 1-D Gaussian:

$$p(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).$$

2.3 The D-Dimensional Case

In -dimensional space:

$$\tilde{p}(x|\mu) = \exp(-\lambda\|x-\mu\|^2),\quad x,\mu\in\mathbb{R}^D.$$

Normalization constant:

$$Z = \int_{\mathbb{R}^D} \exp(-\lambda\|x-\mu\|^2) dx.$$

Shift again with :

$$Z = \int_{\mathbb{R}^D} \exp(-\lambda\|y\|^2) dy = \left(\int_{-\infty}^{\infty} e^{-\lambda y^2} dy\right)^D = \left(\sqrt{\frac{\pi}{\lambda}}\right)^D = \left(\frac{\pi}{\lambda}\right)^{\frac{D}{2}}.$$

Hence:

$$p(x|\mu) = \frac{1}{Z}\exp(-\lambda\|x-\mu\|^2) = \left(\frac{\lambda}{\pi}\right)^{\frac{D}{2}} \exp(-\lambda\|x-\mu\|^2).$$

Setting , we get:

$$p(x|\mu,\sigma^2 I) = \frac{1}{(2\pi\sigma^2)^{D/2}} \exp\left(-\frac{1}{2\sigma^2}\|x-\mu\|^2\right),$$

i.e. the isotropic multivariate Gaussian with covariance .

Conclusion:
“Exponential of the squared distance + normalization” ⇒ Gaussian distribution.

---

3. From a Single Gaussian to the Gaussian Mixture Model (GMM)

Real data is often multimodal ⇒ superimpose several Gaussians.

3.1 Definition of the GMM

Suppose there are Gaussian components, each with parameters:

• Weight: , with
• Mean:
• Covariance matrix: , symmetric positive definite

Then the marginal distribution of the GMM is:

$$p(x|\theta) = \sum_{k=1}^K \pi_k \, \mathcal{N}(x|\mu_k,\Sigma_k),$$

where

$$\mathcal{N}(x|\mu_k,\Sigma_k) = \frac{1}{(2\pi)^{D/2}|\Sigma_k|^{1/2}} \exp\Big(-\frac{1}{2}(x-\mu_k)^\top\Sigma_k^{-1}(x-\mu_k)\Big),$$

and

$$\theta = \{\pi_k,\mu_k,\Sigma_k\}_{k=1}^K.$$

3.2 The Generative Process (Introducing Latent Variables)

To describe the process of “which Gaussian component a sample comes from”, introduce a latent variable :

$$z_n = (z_{n1},\dots,z_{nK})^\top,\quad z_{nk}\in\{0,1\},\quad \sum_{k=1}^K z_{nk}=1.$$

1. First pick a cluster (a discrete choice):

   $$p(z_{nk}=1) = \pi_k,$$

   or, in compact (one-hot) form:

   $$p(z_n) = \prod_{k=1}^K \pi_k^{z_{nk}}.$$

   Explanation: only one $ k^ z_{nk^}=1 $, so:

   $$\prod_{k} \pi_k^{z_{nk}} = \pi_{k^*}^1 \prod_{j\ne k^*}\pi_j^{0} = \pi_{k^*}.$$

2. Once the component is chosen, sample from the corresponding Gaussian:

   $$p(x_n|z_{nk}=1) = \mathcal{N}(x_n|\mu_k,\Sigma_k),$$

   or, in compact form:

   $$p(x_n|z_n) = \prod_{k=1}^K \mathcal{N}(x_n|\mu_k,\Sigma_k)^{z_{nk}}.$$

3.3 The Joint Distribution

By the chain rule:

$$p(x_n,z_n) = p(z_n)\,p(x_n|z_n) = \prod_{k=1}^K \left[ \pi_k \mathcal{N}(x_n|\mu_k,\Sigma_k) \right]^{z_{nk}}.$$

3.4 The Marginal Distribution : Summing over the Latent Variable

Since can only take the K one-hot values , we have:

$$\begin{aligned} p(x_n) &= \sum_{z_n} p(x_n,z_n) \\ &= \sum_{k=1}^K p(x_n,z_n=e_k) \\ &= \sum_{k=1}^K \pi_k \mathcal{N}(x_n|\mu_k,\Sigma_k), \end{aligned}$$

which brings us back to the mixture form.

3.5 The Posterior (Responsibility)

By Bayes’ rule:

$$\gamma_{nk} \triangleq p(z_{nk}=1|x_n) = \frac{p(z_{nk}=1)p(x_n|z_{nk}=1)} {\sum_{j=1}^K p(z_{nj}=1)p(x_n|z_{nj}=1)}.$$

Substituting the GMM parameters:

$$\gamma_{nk} = \frac{\pi_k\mathcal{N}(x_n|\mu_k,\Sigma_k)} {\sum_{j=1}^K \pi_j\mathcal{N}(x_n|\mu_j,\Sigma_j)}.$$

is the responsibility, which expresses “how much component k explains sample ”, and implements soft clustering.

---

4. Maximum Likelihood for the GMM and the EM Algorithm

4.1 Log-Likelihood and Latent Variables

The log-likelihood of the observed data:

$$\log p(X|\theta) = \sum_{n=1}^N \log p(x_n|\theta) = \sum_{n=1}^N \log\left( \sum_{k=1}^K \pi_k\mathcal{N}(x_n|\mu_k,\Sigma_k)\right).$$

Problem: there is a sum inside the log, , which makes it hard to differentiate with respect to and maximize directly.

Introduce the latent variables and define the complete-data log-likelihood:

$$\log p(X,Z|\theta) = \sum_{n=1}^N\sum_{k=1}^K z_{nk}\big[ \log\pi_k + \log\mathcal{N}(x_n|\mu_k,\Sigma_k) \big].$$

The idea of EM: use the “soft label” $\mathbb{E}[z{nk}] z{nk} $.

4.2 Definition of the Q Function

At iteration , with the current parameters known, define:

$$Q(\theta,\theta^{(t)}) = \mathbb{E}_{Z|X,\theta^{(t)}}[\log p(X,Z|\theta)].$$

Since

$$\mathbb{E}_{Z|X,\theta^{(t)}}[z_{nk}] = \gamma_{nk}^{(t)},$$

we get:

$$Q(\theta,\theta^{(t)}) = \sum_{n=1}^N \sum_{k=1}^K \gamma_{nk}^{(t)} \Big( \log\pi_k + \log\mathcal{N}(x_n|\mu_k,\Sigma_k) \Big).$$

4.3 E-Step: Compute the Responsibilities

Given the current parameters :

$$\gamma_{nk}^{(t)} = p(z_{nk}=1|x_n,\theta^{(t)}) = \frac{\pi_k^{(t)} \mathcal{N}(x_n|\mu_k^{(t)},\Sigma_k^{(t)})} {\sum_{j=1}^K \pi_j^{(t)} \mathcal{N}(x_n|\mu_j^{(t)},\Sigma_j^{(t)})}.$$

4.4 M-Step: Maximize Q to Obtain Closed-Form Updates

Let

$$N_k = \sum_{n=1}^N \gamma_{nk}^{(t)}.$$

4.4.1 Updating the Mixing Coefficients

Take the part that depends on :

$$Q_\pi = \sum_{n,k} \gamma_{nk}\log\pi_k.$$

Constraint:

$$\sum_k \pi_k = 1.$$

Form the Lagrangian:

$$\mathcal{L}(\pi,\lambda) = \sum_{k=1}^K \big(\sum_{n=1}^N \gamma_{nk}\big)\log\pi_k + \lambda\left(\sum_{k=1}^K \pi_k -1\right).$$

Differentiate with respect to :

$$\frac{\partial\mathcal{L}}{\partial \pi_k} = \frac{\sum_n\gamma_{nk}}{\pi_k} + \lambda = 0 \quad\Rightarrow\quad \pi_k = -\frac{\sum_n\gamma_{nk}}{\lambda}.$$

Apply the constraint:

$$\sum_k \pi_k = 1 \Rightarrow -\frac{1}{\lambda}\sum_{k}\sum_n\gamma_{nk} = 1.$$

Also, since for every n, $\sum_k\gamma_{nk}=1$, we have:

$$\sum_k\sum_n\gamma_{nk} = N \Rightarrow -\frac{N}{\lambda} = 1 \Rightarrow \lambda = -N.$$

Hence:

$$\pi_k^{(t+1)} = \frac{\sum_n\gamma_{nk}}{N} = \frac{N_k}{N}.$$

4.4.2 Updating the Means

Take only the part of Q that depends on :

$$Q_{\mu,\Sigma} = \sum_{n,k} \gamma_{nk} \log\mathcal{N}(x_n|\mu_k,\Sigma_k).$$

Expand the Gaussian log:

$$\log\mathcal{N}(x_n|\mu_k,\Sigma_k) = C - \frac{1}{2}(x_n-\mu_k)^\top\Sigma_k^{-1}(x_n-\mu_k),$$

where does not depend on .

Therefore

$$Q_{\mu_k} = -\frac{1}{2}\sum_n \gamma_{nk} (x_n-\mu_k)^\top\Sigma_k^{-1}(x_n-\mu_k).$$

Differentiate with respect to :

$$\frac{\partial}{\partial\mu_k} (x_n-\mu_k)^\top\Sigma_k^{-1}(x_n-\mu_k) = -2\Sigma_k^{-1}(x_n-\mu_k)$$

So

$$\frac{\partial Q_{\mu_k}}{\partial\mu_k} = -\frac{1}{2} \sum_n \gamma_{nk}\big(-2\Sigma_k^{-1}(x_n-\mu_k)\big) = \sum_n \gamma_{nk} \Sigma_k^{-1}(x_n-\mu_k).$$

Set the derivative to 0:

$$\sum_n \gamma_{nk} \Sigma_k^{-1}(x_n-\mu_k) = 0 \Rightarrow \sum_n \gamma_{nk} (x_n-\mu_k) = 0.$$

Expanding:

$$\sum_n \gamma_{nk} x_n = \mu_k \sum_n \gamma_{nk} = \mu_k N_k,$$

so:

$$\mu_k^{(t+1)} = \frac{1}{N_k}\sum_{n=1}^N\gamma_{nk}x_n.$$

4.4.3 Updating the Covariances

Again start from:

$$Q_{\Sigma_k} = \sum_n \gamma_{nk}\Big( -\frac{1}{2}\log|\Sigma_k| -\frac{1}{2}(x_n-\mu_k)^\top\Sigma_k^{-1}(x_n-\mu_k) \Big)$$

and differentiate with respect to .

Using the matrix-derivative formulas (for symmetric matrices):

we obtain:

$$\frac{\partial Q_{\Sigma_k}}{\partial\Sigma_k} = -\frac{1}{2}\sum_n\gamma_{nk}\left[ \Sigma_k^{-1} - \Sigma_k^{-1}(x_n-\mu_k)(x_n-\mu_k)^\top\Sigma_k^{-1} \right].$$

Set the derivative to 0:

$$\sum_n\gamma_{nk}\left[ \Sigma_k^{-1}(x_n-\mu_k)(x_n-\mu_k)^\top\Sigma_k^{-1} -\Sigma_k^{-1} \right]=0.$$

Multiply by on the left and on the right:

$$\sum_n\gamma_{nk}\left[ (x_n-\mu_k)(x_n-\mu_k)^\top - \Sigma_k \right]=0.$$

→

$$\sum_n\gamma_{nk}(x_n-\mu_k)(x_n-\mu_k)^\top = \sum_n\gamma_{nk}\Sigma_k = N_k \Sigma_k.$$

Therefore:

$$\Sigma_k^{(t+1)} = \frac{1}{N_k}\sum_{n=1}^N \gamma_{nk} (x_n-\mu_k^{(t+1)})(x_n-\mu_k^{(t+1)})^\top.$$

---

5. EM Monotonicity: Why Does Each Step Never Decrease the Log-Likelihood?

5.1 Constructing the Variational Lower Bound (ELBO)

For any distribution :

$$\log p(X|\theta) = \log \sum_Z p(X,Z|\theta) = \log \sum_Z q(Z)\frac{p(X,Z|\theta)}{q(Z)}.$$

By Jensen’s inequality:

$$\log p(X|\theta) \ge \sum_Z q(Z)\log\frac{p(X,Z|\theta)}{q(Z)} \triangleq \mathcal{F}(q,\theta).$$

Define:

$$\mathcal{F}(q,\theta) = \mathbb{E}_q[\log p(X,Z|\theta)] + H(q),$$

where is the entropy.

5.2 KL Decomposition: The Relation between the Log-Likelihood and the Lower Bound

One can show:

$$\log p(X|\theta) = \mathcal{F}(q,\theta) + \mathrm{KL}\big(q(Z)\,\|\,p(Z|X,\theta)\big),$$

where the KL divergence is always non-negative.

Derivation:

$$\begin{aligned} \mathcal{F}(q,\theta) &= \sum_Z q(Z)\log\frac{p(X,Z|\theta)}{q(Z)} \\ &= \sum_Z q(Z)\log\frac{p(Z|X,\theta)p(X|\theta)}{q(Z)} \\ &= \log p(X|\theta) + \sum_Z q(Z)\log\frac{p(Z|X,\theta)}{q(Z)} \\ &= \log p(X|\theta) - \mathrm{KL}\big(q(Z)\,\|\,p(Z|X,\theta)\big). \end{aligned}$$

Rearranging gives the result.

5.3 E-Step: Fix , Find the Optimal

With fixed, is a constant.

We want to maximize ,
which is equivalent to minimizing the KL:

$$\mathrm{KL}(q(Z)\,\|\,p(Z|X,\theta^{(t)})).$$

The minimum is 0, attained if and only if:

$$q^{(t)}(Z) = p(Z|X,\theta^{(t)}).$$

At that point:

$$\mathcal{F}(q^{(t)},\theta^{(t)}) = \log p(X|\theta^{(t)}).$$

That is: the E-step “pulls” the lower bound up to the current log-likelihood, making it tangent to it.

5.4 M-Step: Fix , Update

In the M-step, by definition:

$$\theta^{(t+1)} = \arg\max_\theta \mathcal{F}(q^{(t)},\theta),$$

so:

$$\mathcal{F}(q^{(t)},\theta^{(t+1)}) \ge \mathcal{F}(q^{(t)},\theta^{(t)}) = \log p(X|\theta^{(t)}).$$

On the other hand, for the new parameters the decomposition still holds:

$$\log p(X|\theta^{(t+1)}) = \mathcal{F}(q^{(t)},\theta^{(t+1)}) + \mathrm{KL}\big(q^{(t)}\,\|\,p(Z|X,\theta^{(t+1)})\big) \ge \mathcal{F}(q^{(t)},\theta^{(t+1)}).$$

Combining the two:

$$\log p(X|\theta^{(t+1)}) \ge \mathcal{F}(q^{(t)},\theta^{(t+1)}) \ge \mathcal{F}(q^{(t)},\theta^{(t)}) = \log p(X|\theta^{(t)}).$$

Therefore:

$$\boxed{ \log p(X|\theta^{(t+1)}) \ge \log p(X|\theta^{(t)}), }$$

i.e. no EM iteration ever decreases the log-likelihood.

5.5 The Relation between the Q Function and

After the E-step chooses :

$$\mathcal{F}(q^{(t)},\theta) = \mathbb{E}_{Z|X,\theta^{(t)}}[\log p(X,Z|\theta)] + H(q^{(t)}),$$

where does not depend on .

Therefore, maximizing in the M-step is equivalent to maximizing:

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

which is exactly the Q function used in the M-step derivation above.

---

6. K-means as a Limiting Case of GMM-EM

Simplify the GMM with the following constraints:

1. All components share the same spherical covariance: $$\Sigma_k = \sigma^2 I,\quad \forall k.$$
2. Uniform weights: $$\pi_k = \frac{1}{K},\quad \forall k.$$

Then the responsibilities:

$$\gamma_{nk} = \frac{\exp\left(-\frac{1}{2\sigma^2}\|x_n-\mu_k\|^2\right)} {\sum_{j=1}^K \exp\left(-\frac{1}{2\sigma^2}\|x_n-\mu_j\|^2\right)}$$

are a “softmax over squared distances”.

As :

• The exponential term of the nearest cluster dominates all the others by far;
• The softmax degenerates into the one-hot argmin assignment.

That is:

$$\gamma_{nk} \to r_{nk} = \begin{cases} 1, & k = \arg\min_j \|x_n - \mu_j\|^2,\\ 0, & \text{otherwise}. \end{cases}$$

At this point:

• E-step ⇒ hard assignment to the nearest center (the K-means assignment step);
• M-step: $$\mu_k = \frac{\sum_n \gamma_{nk}x_n}{\sum_n \gamma_{nk}} \to \frac{\sum_n r_{nk}x_n}{\sum_n r_{nk}},$$ i.e. the K-means update.

Therefore:

K-means = the Hard-EM special case of GMM-EM in the limit of “covariance → 0, uniform weights”.

---

7. Summary

1. K-means: minimizes the within-cluster sum of squared distances;

   • a geometric model with no probabilistic structure;
   • hard (0/1) assignments.

2. From distance to Gaussian:

   • exponential energy model ;
   • normalization ⇒ Gaussian distribution.

3. GMM: a weighted sum of several Gaussians;

   • introduce the one-hot latent variable to build the joint ;
   • marginalize ⇒ the mixture form ;
   • the posterior ⇒ soft responsibilities.

4. EM algorithm:

   • E-step: compute $\gamma_{nk}=p(z_{nk}=1|x_n,\theta^{(t)})$
   • M-step: maximize Q to get closed-form updates for .

5. EM monotonicity:

   • construct the variational lower bound and use $$\log p(X|\theta) = \mathcal{F}(q,\theta) + \mathrm{KL}(q\|p(Z|X,\theta))$$
   • the E-step makes KL=0, so the lower bound equals the current log-likelihood;
   • the M-step raises the lower bound ⇒ the log-likelihood does not decrease.

6. The relation between K-means and GMM:

   • GMM is the probabilistic generalization of K-means;
   • K-means is the Hard-EM limit of GMM-EM with shared spherical covariance and .

---

Translated from the Chinese original.