Unsupervised Learning

In supervised learning we have a dataset ${(x_i, y_i)}{i=1}^n$ and we want to learn $p(y|x)$ or a function $f: x \mapsto y$. In unsupervised learning we have only ${x_i}{i=1}^n$ — inputs with no labels — and we want to learn $p(x)$ itself, or some useful structure derived from it.

There are three canonical subproblems:

1. Clustering: Partition the data into $K$ groups such that points within a group are similar and points across groups are dissimilar. The challenge is that "similarity" and "the right $K$" must be defined without external ground truth.

2. Density estimation: Fit a model $\hat{p}(x)$ to the data. If you succeed, you can generate new samples, detect anomalies (low-$p$ points), and perform Bayesian inference.

3. Dimensionality reduction / representation learning: Find a mapping $z = f(x)$ where $z \in \mathbb{R}^d$ with $d \ll D$ (original dimension) that preserves as much structure as possible. This enables visualization, compression, and downstream supervised tasks.

These are not entirely separate — GMMs simultaneously cluster (which component did $x_i$ come from?) and do density estimation (what is $p(x)$?). PCA does dimensionality reduction but also implicitly defines a density model (the principal subspace captures most variance).

The central difficulty is the lack of a universal loss function. In supervised learning, prediction error on held-out data is a gold standard. In unsupervised learning, evaluation is harder: we use reconstruction error, held-out log-likelihood, silhouette scores, or simply downstream task performance.

K-means is the canonical clustering algorithm. Given $n$ points $x_1, \ldots, x_n \in \mathbb{R}^D$ and a desired number of clusters $K$, it finds cluster centers (centroids) $\mu_1, \ldots, \mu_K \in \mathbb{R}^D$ and an assignment $c_i \in {1, \ldots, K}$ for each point that minimizes the within-cluster sum of squares (WCSS):

$$J = \sum_{i=1}^n |x_i - \mu_{c_i}|^2$$

This is a hard combinatorial problem (NP-hard in general), but the Lloyd's algorithm gives a practical alternating minimization:

E-step (assign): Fix $\mu_1,\ldots,\mu_K$. Assign each point to its nearest centroid: $$c_i \leftarrow \arg\min_k |x_i - \mu_k|^2$$

M-step (update): Fix assignments. Set each centroid to the mean of its assigned points: $$\mu_k \leftarrow \frac{1}{|C_k|} \sum_{i: c_i = k} x_i$$

Repeat until convergence (assignments stop changing). This is guaranteed to converge because $J$ decreases monotonically at each step — the E-step is the unique minimizer over assignments given fixed centroids, and the M-step is the unique minimizer over centroids given fixed assignments. Since $J \geq 0$ and the number of possible assignments is finite, termination is guaranteed.

Analogy: Imagine placing $K$ magnets on a map and letting cities snap to their nearest magnet. Then move each magnet to the center of gravity of its city cluster. Repeat. The magnets settle.

Initialization matters enormously. A bad random initialization leads to local optima. The K-means++ initialization (Arthur & Vassilvitskii, 2007) fixes this: pick the first centroid uniformly at random, then pick each subsequent centroid with probability proportional to $D(x)^2$ (squared distance to the nearest already-chosen centroid). This gives an $O(\log K)$ approximation guarantee and dramatically better practical results. Always use K-means++ in practice.

Choosing K: Use the elbow method — plot $J$ vs. $K$ and look for a kink. More principled: silhouette score, which measures how tightly a point coheres to its own cluster vs. the next nearest cluster: $$s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))}$$ where $a(i)$ is mean intra-cluster distance and $b(i)$ is mean distance to the nearest other cluster. $s(i) \in [-1, 1]$; higher is better.

Limitations: K-means assumes spherical, equally-sized clusters (because it uses Euclidean distance to flat centroids). It fails on elongated clusters, clusters of very different density, or clusters with non-convex shapes. The Gaussian mixture model generalizes K-means to handle these cases.

import numpy as np

def kmeans_plusplus_init(X, K, rng):
    n = X.shape[0]
    centers = [X[rng.integers(n)]]
    for _ in range(K - 1):
        dists = np.min([np.sum((X - c)**2, axis=1) for c in centers], axis=0)
        probs = dists / dists.sum()
        centers.append(X[rng.choice(n, p=probs)])
    return np.array(centers)

def kmeans(X, K, max_iter=300, seed=42):
    rng = np.random.default_rng(seed)
    mu = kmeans_plusplus_init(X, K, rng)          # (K, D)
    assignments = np.zeros(X.shape[0], dtype=int)
    for _ in range(max_iter):
        # E-step: assign each point to nearest centroid
        dists = np.linalg.norm(X[:, None, :] - mu[None, :, :], axis=2)  # (n, K)
        new_assignments = np.argmin(dists, axis=1)
        if np.all(new_assignments == assignments):
            break
        assignments = new_assignments
        # M-step: update centroids
        for k in range(K):
            members = X[assignments == k]
            if len(members) > 0:
                mu[k] = members.mean(axis=0)
    return assignments, mu

K-means requires you to commit to a fixed $K$ upfront. Hierarchical clustering avoids this by building a full tree (dendrogram) of cluster merges, which you can cut at any level to get any $K$.

Agglomerative (bottom-up) is by far the most common variant:

1. Start: each point is its own cluster. $n$ clusters total.
2. Repeat: merge the two closest clusters according to a linkage criterion until only 1 cluster remains.
3. The full history of merges is stored as a dendrogram — a binary tree where the height of each merge encodes the distance at which it happened.

The linkage criterion defines distance between clusters $A$ and $B$:

• Single linkage (nearest neighbor): $d(A,B) = \min_{a \in A, b \in B} d(a,b)$. Tends to produce long, stringy clusters (chaining effect).
• Complete linkage (farthest neighbor): $d(A,B) = \max_{a \in A, b \in B} d(a,b)$. Produces compact, roughly spherical clusters.
• Average linkage (UPGMA): $d(A,B) = \frac{1}{|A||B|}\sum_{a,b} d(a,b)$. A good compromise, widely used in bioinformatics.
• Ward's linkage: Merge the pair that minimizes the increase in total within-cluster variance. This is the most similar to K-means and usually gives the most interpretable results.

Reading the dendrogram: The y-axis is merge distance. A horizontal cut at height $h$ produces clusters: find all branches that cross the cut. This lets you explore cluster structure at multiple scales without rerunning the algorithm.

Complexity: Naive implementation is $O(n^3)$ time (recompute all pairwise distances after each merge). With priority queues and Lance-Williams update formulas, single/complete/average linkage can be done in $O(n^2)$ time and $O(n^2)$ space. Ward's is $O(n^2 \log n)$ using a heap. This makes hierarchical clustering unscalable beyond $\sim$10,000 points — use K-means or DBSCAN for large $n$.

Divisive (top-down) clustering starts with all points in one cluster and recursively splits. It is rarely used because choosing the best split is hard.

When to use: Hierarchical clustering shines when you genuinely do not know $K$, when you want to understand cluster structure at multiple scales (e.g., phylogenetic trees in biology), or when $n$ is small enough that $O(n^2)$ is acceptable. The dendrogram is often itself a scientific result.

from scipy.cluster.hierarchy import linkage, dendrogram, fcluster
import matplotlib.pyplot as plt
import numpy as np

# X: (n, D) data matrix
X = np.random.randn(30, 2)  # toy data

# Compute Ward linkage (returns condensed linkage matrix)
Z = linkage(X, method='ward')

# Plot dendrogram — the y-axis is merge distance
fig, ax = plt.subplots(figsize=(10, 4))
dendrogram(Z, ax=ax)
ax.set_title('Ward Hierarchical Clustering Dendrogram')

# Cut the tree at K=3 clusters
labels = fcluster(Z, t=3, criterion='maxclust')  # returns 1-indexed labels

K-means makes a hard assignment: every point belongs to exactly one cluster. But real data is messy — a point halfway between two clusters might genuinely belong to either. Gaussian Mixture Models (GMMs) replace hard assignments with soft, probabilistic ones.

A GMM with $K$ components models the data density as a weighted sum of Gaussians: $$p(x) = \sum_{k=1}^K \pi_k , \mathcal{N}(x \mid \mu_k, \Sigma_k)$$

where $\pi_k \geq 0$, $\sum_k \pi_k = 1$ are mixing weights, $\mu_k$ are component means, and $\Sigma_k$ are covariance matrices.

Generative story: To sample a point — first pick a component $z \sim \text{Categorical}(\pi_1, \ldots, \pi_K)$, then sample $x \sim \mathcal{N}(\mu_z, \Sigma_z)$. The variable $z$ is latent (hidden, unobserved).

K-means as a degenerate GMM: If all $\Sigma_k = \sigma^2 I$ and $\sigma \to 0$, the GMM posterior $p(z=k \mid x)$ collapses to a 1-hot vector (the nearest centroid wins), recovering K-means exactly.

Parameters: $\theta = {\pi_k, \mu_k, \Sigma_k}{k=1}^K$. We want the MLE: $$\hat{\theta} = \arg\max\theta \sum_{i=1}^n \log p(x_i; \theta) = \arg\max_\theta \sum_{i=1}^n \log \sum_{k=1}^K \pi_k \mathcal{N}(x_i \mid \mu_k, \Sigma_k)$$

The sum inside the log makes this intractable in closed form — the EM algorithm handles this.

Covariance structure choices (from most to least flexible):

• Full: $\Sigma_k$ is a general $D \times D$ PSD matrix ($O(D^2)$ parameters per component)
• Diagonal: $\Sigma_k = \text{diag}(\sigma_{k1}^2, \ldots, \sigma_{kD}^2)$ — axes-aligned ellipses
• Spherical: $\Sigma_k = \sigma_k^2 I$ — circles of varying size
• Tied: all components share one $\Sigma$ — like LDA

Model selection: Use the Bayesian Information Criterion $\text{BIC} = -2 \log \hat{L} + p \log n$ where $p$ is the number of free parameters. Pick $K$ minimizing BIC. BIC penalizes model complexity, preventing overfitting (a single Gaussian centered on each point perfectly memorizes the training set but generalizes poorly).

from sklearn.mixture import GaussianMixture
import numpy as np

# Fit GMMs with K=1..6 and pick K via BIC
X = np.vstack([
    np.random.randn(100, 2) + [0, 0],
    np.random.randn(80,  2) + [4, 0],
    np.random.randn(120, 2) + [2, 3],
])

best_bic, best_k, best_gmm = np.inf, None, None
for k in range(1, 7):
    gmm = GaussianMixture(n_components=k, covariance_type='full',
                          n_init=5, random_state=42)
    gmm.fit(X)
    bic = gmm.bic(X)
    if bic < best_bic:
        best_bic, best_k, best_gmm = bic, k, gmm

print(f"Best K={best_k}, BIC={best_bic:.1f}")
labels = best_gmm.predict(X)          # hard cluster assignment
probs  = best_gmm.predict_proba(X)    # soft posteriors r_ik

EM is one of the most important algorithms in statistics and machine learning. It is the standard tool for maximum likelihood estimation with latent variables.

The setup: We have observed data $X = {x_i}$, latent variables $Z = {z_i}$, and parameters $\theta$. The complete-data log-likelihood $\log p(X, Z; \theta)$ would be easy to maximize if we knew $Z$. But we only observe $X$, so we must maximize the incomplete-data log-likelihood $\log p(X; \theta) = \log \sum_Z p(X, Z; \theta)$, which has a sum inside the log.

The key insight — Jensen's inequality: For any distribution $q(Z)$, $$\log p(X; \theta) = \log \mathbb{E}{q} \left[\frac{p(X, Z; \theta)}{q(Z)}\right] \geq \mathbb{E}{q}\left[\log \frac{p(X, Z; \theta)}{q(Z)}\right] \equiv \mathcal{L}(q, \theta)$$ This lower bound $\mathcal{L}$ is called the ELBO (Evidence Lower BOund). The gap between $\log p(X;\theta)$ and $\mathcal{L}$ is exactly $\text{KL}(q(Z) | p(Z|X;\theta)) \geq 0$.

E-step (Expectation): Set $q(Z) = p(Z|X; \theta^{\text{old}})$ — the posterior over latent variables given current parameters. This sets KL to zero, making $\mathcal{L} = \log p(X; \theta^{\text{old}})$. Equivalently, compute the responsibilities: $$r_{ik} = p(z_i = k \mid x_i; \theta^{\text{old}}) = \frac{\pi_k \mathcal{N}(x_i \mid \mu_k, \Sigma_k)}{\sum_{j} \pi_j \mathcal{N}(x_i \mid \mu_j, \Sigma_j)}$$

M-step (Maximization): Maximize $\mathcal{L}(q, \theta)$ over $\theta$ with $q$ fixed. This is maximizing the expected complete-data log-likelihood: $$\theta^{\text{new}} = \arg\max_\theta \sum_i \sum_k r_{ik} \log \left[ \pi_k \mathcal{N}(x_i \mid \mu_k, \Sigma_k) \right]$$ The closed-form updates for GMMs are (let $N_k = \sum_i r_{ik}$): $$\pi_k^{\text{new}} = \frac{N_k}{n}, \quad \mu_k^{\text{new}} = \frac{\sum_i r_{ik} x_i}{N_k}, \quad \Sigma_k^{\text{new}} = \frac{\sum_i r_{ik}(x_i - \mu_k^{\text{new}})(x_i - \mu_k^{\text{new}})^T}{N_k}$$

Why does EM converge? After each E-step, $\mathcal{L} = \log p(X; \theta^{\text{old}})$. After each M-step, $\theta^{\text{new}}$ increases $\mathcal{L}$ (or keeps it flat). Since $\mathcal{L} \leq \log p(X; \theta)$ and $\mathcal{L}(q^{\text{new}}, \theta^{\text{new}}) = \log p(X; \theta^{\text{new}})$ at the top of the next E-step, the marginal log-likelihood is non-decreasing. EM converges to a local maximum (not necessarily global).

EM is not just for GMMs: It applies to any model with latent variables — HMMs (Baum-Welch is EM), factor analysis, probabilistic PCA, topic models (LDA), and more. The pattern is always the same: introduce latent variables that make the complete-data likelihood easy, then alternate between inferring them and updating parameters.

Analogy: Imagine you are grading exams but you do not know which of two teaching assistants wrote each solution. In the E-step, you probabilistically attribute each solution to each TA based on current estimates of their writing style. In the M-step, you update your model of each TA's style using the probabilistic attributions. Repeat.

import numpy as np

def em_gmm(X, K, max_iter=100, tol=1e-6):
    n, D = X.shape
    # Initialize with K-means++ style random picks
    rng = np.random.default_rng(0)
    idx = rng.choice(n, K, replace=False)
    mu = X[idx].copy()                               # (K, D)
    Sigma = np.array([np.eye(D)] * K)                # (K, D, D)
    pi = np.ones(K) / K                              # (K,)

    log_lik_prev = -np.inf
    for iteration in range(max_iter):
        # ---- E-step: compute responsibilities r[i,k] ----
        log_r = np.zeros((n, K))
        for k in range(K):
            diff = X - mu[k]                         # (n, D)
            sign, logdet = np.linalg.slogdet(Sigma[k])
            inv_S = np.linalg.inv(Sigma[k])
            maha = np.einsum('ni,ij,nj->n', diff, inv_S, diff)  # (n,)
            log_r[:, k] = (np.log(pi[k])
                           - 0.5 * logdet
                           - 0.5 * maha
                           - 0.5 * D * np.log(2 * np.pi))
        # log-sum-exp for numerical stability
        log_r_max = log_r.max(axis=1, keepdims=True)
        log_norm = log_r_max + np.log(np.exp(log_r - log_r_max).sum(axis=1, keepdims=True))
        r = np.exp(log_r - log_norm)                 # (n, K)
        log_lik = log_norm.sum()

        # ---- M-step: update parameters ----
        Nk = r.sum(axis=0)                           # (K,)
        pi = Nk / n
        mu = (r.T @ X) / Nk[:, None]                # (K, D)
        for k in range(K):
            diff = X - mu[k]
            Sigma[k] = (r[:, k:k+1] * diff).T @ diff / Nk[k]
            Sigma[k] += 1e-6 * np.eye(D)            # regularize

        if abs(log_lik - log_lik_prev) < tol:
            break
        log_lik_prev = log_lik
    return pi, mu, Sigma, r

PCA is the workhorse of dimensionality reduction. Given $n$ data points $x_1, \ldots, x_n \in \mathbb{R}^D$ (mean-centered: $\bar{x} = 0$), PCA finds the $d$-dimensional linear subspace that captures the most variance.

Formulation: Find unit vectors $u_1, \ldots, u_d$ (the principal components) such that the projected data $z_i = U^T x_i$ (where $U = [u_1 \cdots u_d]$) has maximum variance. The first PC $u_1$ solves: $$u_1 = \arg\max_{|u|=1} \frac{1}{n} \sum_i (u^T x_i)^2 = \arg\max_{|u|=1} u^T \Sigma u$$ where $\Sigma = \frac{1}{n} X^T X$ is the sample covariance matrix ($X$ is $n \times D$). By the Lagrangian, this is the eigenvector of $\Sigma$ with the largest eigenvalue $\lambda_1$. Each subsequent PC is the next eigenvector, constrained to be orthogonal to all previous ones.

The SVD connection: Let $X = U_s \Lambda V^T$ be the SVD of the (centered) data matrix (where $U_s \in \mathbb{R}^{n \times r}$, $\Lambda = \text{diag}(\sigma_1, \ldots)$, $V \in \mathbb{R}^{D \times r}$). Then the principal components are the right singular vectors (columns of $V$), and $\lambda_k = \sigma_k^2 / n$. In practice, always use SVD rather than computing $\Sigma$ explicitly (it is numerically more stable and works even when $D > n$).

Choosing $d$: Plot the explained variance ratio $\lambda_k / \sum_j \lambda_j$ and the cumulative sum. A common rule: keep enough PCs to explain 90-95% of variance. The scree plot (variance vs. PC index) often has a sharp elbow.

Reconstruction: Given codes $z_i = U^T x_i$, reconstruct $\hat{x}_i = U z_i$. The reconstruction error $\sum_i |x_i - \hat{x}i|^2 = \sum{k > d} \lambda_k$ is the sum of discarded eigenvalues.

Probabilistic PCA (PPCA): PCA has a probabilistic interpretation as a latent variable model: $x = Wz + \mu + \epsilon$ where $z \sim \mathcal{N}(0,I)$ and $\epsilon \sim \mathcal{N}(0, \sigma^2 I)$. As $\sigma \to 0$, the MLE of $W$ recovers classical PCA. PPCA enables missing data imputation via EM and gives a proper generative model.

Whitening: After projecting to $d$ PCs, you can whiten (spherize) by dividing each PC by $\sqrt{\lambda_k}$, yielding unit-variance coordinates. This is important as preprocessing for algorithms that assume isotropic input (e.g., ICA, some neural network initializations).

Limitations: PCA is a linear method. It will fail to discover curved manifolds, clusters arranged in concentric rings, or any nonlinear structure. For those cases, kernel PCA or t-SNE are needed.

import numpy as np

def pca(X, d):
    """PCA via SVD. X: (n, D). Returns codes Z (n, d), components V (D, d),
       explained variance ratios."""
    # 1. Center the data
    mu = X.mean(axis=0)
    Xc = X - mu

    # 2. Thin SVD — numerically stable, works for D > n
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)  # s: (min(n,D),)
    # Vt rows = principal components; s^2/n = eigenvalues

    # 3. Keep top-d components
    V = Vt[:d].T                        # (D, d)  — columns are PCs
    Z = Xc @ V                          # (n, d)  — low-dim codes

    # 4. Explained variance
    explained_var_ratio = (s**2) / (s**2).sum()

    return Z, V, explained_var_ratio[:d], mu

def reconstruct(Z, V, mu):
    """Reconstruct from PCA codes back to original space."""
    return Z @ V.T + mu

# Example usage
X = np.random.randn(200, 50)           # 200 samples, 50 dims
Z, V, evr, mu = pca(X, d=10)          # compress to 10 dims
print(f"Explained variance (10 PCs): {evr.sum()*100:.1f}%")
Xhat = reconstruct(Z, V, mu)           # (200, 50) reconstruction

t-SNE (van der Maaten & Hinton, 2008) is the dominant method for visualizing high-dimensional data in 2D or 3D. It produces stunning, interpretable visualizations where clusters are clearly separated — but it must be understood carefully, as it is frequently misinterpreted.

Core idea: Preserve local neighborhood structure. Nearby points in high-D space should remain nearby in low-D; distant points' relative positions are not preserved. The method:

Step 1 — High-dimensional similarities: Define a symmetric probability distribution over pairs of points in the high-dimensional space. For each pair $(i, j)$, define the similarity as a Gaussian kernel: $$p_{j|i} = \frac{\exp(-|x_i - x_j|^2 / 2\sigma_i^2)}{\sum_{k \neq i} \exp(-|x_i - x_k|^2 / 2\sigma_i^2)}, \quad p_{ij} = \frac{p_{j|i} + p_{i|j}}{2n}$$ The bandwidth $\sigma_i$ is set per-point to achieve a target perplexity $\text{Perp}(P_i) = 2^{H(P_i)}$ where $H$ is the Shannon entropy. Perplexity loosely controls the number of effective neighbors (typically 5–50).

Step 2 — Low-dimensional similarities: In the map space $y_i \in \mathbb{R}^2$, use a Student-t distribution with 1 degree of freedom (Cauchy kernel) rather than a Gaussian: $$q_{ij} = \frac{(1 + |y_i - y_j|^2)^{-1}}{\sum_{k \neq l}(1 + |y_k - y_l|^2)^{-1}}$$

Why the t-distribution? The "crowding problem": a $D$-sphere has volume $\propto r^D$, but a 2D circle has volume $\propto r^2$. In high dimensions, the moderately-nearby neighbors must be mapped to a 2D circle that is too small — they crowd together. Using a heavy-tailed t-distribution in low-D gives extra room for moderately distant points, resolving this.

Step 3 — Minimize KL divergence: Minimize $\text{KL}(P | Q) = \sum_{i \neq j} p_{ij} \log \frac{p_{ij}}{q_{ij}}$ over the map points $y_i$ using gradient descent with momentum: $$\frac{\partial \text{KL}}{\partial y_i} = 4 \sum_j (p_{ij} - q_{ij})(y_i - y_j)(1 + |y_i - y_j|^2)^{-1}$$

Intuitively: when $p_{ij} > q_{ij}$, nearby high-D neighbors are too far in 2D — they attract. When $p_{ij} < q_{ij}$, distant high-D points are too close in 2D — they repel.

Critical warnings for correct interpretation (following Wattenberg et al., Distill 2016):

1. Cluster sizes are meaningless — t-SNE normalizes $Q$, so relative cluster densities are not preserved.
2. Distances between clusters are meaningless — only local structure is trustworthy.
3. Perplexity matters a lot — run with multiple values (5, 30, 50); different perplexities reveal different scales of structure.
4. Random initializations can look very different — run multiple times; trust features that are stable.
5. t-SNE is not for inference — it is a visualization tool, not a generative model or a metric space.

Computational cost: Naive t-SNE is $O(n^2)$. The Barnes-Hut t-SNE approximation (the default in scikit-learn) reduces this to $O(n \log n)$ by treating far-away points as clusters via a quad-tree. For very large datasets ($n > 100\text{k}$), use UMAP (Uniform Manifold Approximation and Projection), which preserves more global structure and is $O(n)$.

from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler
import numpy as np

def visualize_tsne(X, labels=None, perplexity=30, n_iter=1000):
    """
    Best practice t-SNE pipeline:
    1. PCA to 50 dims first (denoises, speeds up t-SNE)
    2. Standardize
    3. Run t-SNE
    """
    from sklearn.decomposition import PCA
    # Step 1: reduce to 50 dims with PCA (standard practice for high-D data)
    d_init = min(50, X.shape[1])
    X_pca = PCA(n_components=d_init, random_state=42).fit_transform(X)

    # Step 2: standardize
    X_scaled = StandardScaler().fit_transform(X_pca)

    # Step 3: t-SNE
    tsne = TSNE(
        n_components=2,
        perplexity=perplexity,
        n_iter=n_iter,
        init='pca',          # PCA init is more stable than random
        learning_rate='auto',
        random_state=42
    )
    Y = tsne.fit_transform(X_scaled)   # (n, 2) map coordinates
    return Y

# NOTE: interpret only local cluster membership, NOT cluster size/distance

The absence of labels makes unsupervised evaluation genuinely hard. Here is a systematic toolkit:

Internal indices (no ground truth needed):

• Inertia (WCSS): Sum of squared distances to centroids. Always decreases with $K$ — use only relative to identify the elbow, never as an absolute measure.
• Silhouette score: $s(i) = (b(i) - a(i)) / \max(a(i), b(i))$ where $a(i)$ is mean intra-cluster distance and $b(i)$ is mean distance to the nearest other cluster. Range $[-1, 1]$; higher is better. Works for any clustering, not just K-means.
• Davies-Bouldin index: Average of max cluster-pair ratios of scatter to separation. Lower is better.
• Calinski-Harabasz index: Ratio of between-cluster to within-cluster dispersion. Higher is better.

External indices (requires ground truth, used for benchmarking):

• Adjusted Rand Index (ARI): Measures agreement between predicted and true partitions, adjusted for chance. Range $[-1, 1]$; $1$ = perfect.
• Normalized Mutual Information (NMI): $\text{NMI}(U, V) = 2I(U;V) / (H(U) + H(V))$. Range $[0, 1]$.

Likelihood-based (for probabilistic models like GMMs):

• Held-out log-likelihood: Fit on training set, evaluate $\sum_i \log p(x_i^{\text{test}})$ on test set. The gold standard for generative models.
• BIC/AIC: Penalized training log-likelihood. BIC penalizes more aggressively and is more consistent for model selection.

For dimensionality reduction:

• Reconstruction error: $\sum_i |x_i - \hat{x}_i|^2$ (PCA, autoencoders)
• Downstream task performance: Does the embedding help a supervised classifier? Often the most meaningful measure in practice.
• Trustworthiness and Continuity (Venna & Kaski 2006): Formalize the notion that local neighborhoods are preserved in both directions.

The fundamental tension: A model that memorizes the data (one cluster per point) achieves perfect internal metrics but is useless. Model selection must balance fit and complexity — this is the deep reason BIC and held-out likelihood are preferred over training metrics.

from sklearn.metrics import silhouette_score, adjusted_rand_score
from sklearn.cluster import KMeans
import numpy as np

def select_k_kmeans(X, k_range=range(2, 11), n_init=10):
    """Select K for K-means via silhouette score."""
    scores = {}
    for k in k_range:
        km = KMeans(n_clusters=k, n_init=n_init, random_state=42).fit(X)
        sil = silhouette_score(X, km.labels_)
        inertia = km.inertia_
        scores[k] = {'silhouette': sil, 'inertia': inertia}
        print(f"K={k}: silhouette={sil:.3f}, inertia={inertia:.1f}")
    best_k = max(scores, key=lambda k: scores[k]['silhouette'])
    print(f"\nBest K by silhouette: {best_k}")
    return best_k, scores