Unit 1: Introduction to Unsupervised Learning & K-Means

1. Foundations of the Unsupervised Learning Paradigm

1.1 Machine Learning Taxonomy and Unsupervised Learning Principles

Machine learning algorithms are traditionally classified into four primary paradigms based on the nature of the learning feedback and the presence or absence of supervisory ground-truth signals:

1. Supervised Learning: The learning algorithm is provided with a labeled training dataset $\mathcal{D} = \{(x^{(1)}, y^{(1)}), (x^{(2)}, y^{(2)}), \dots, (x^{(N)}, y^{(N)})\}$, where each instance $x^{(i)} \in \mathcal{X} \subseteq \mathbb{R}^d$ is paired with a continuous target or discrete category label $y^{(i)} \in \mathcal{Y}$. The objective is to induce a hypothesis function $f: \mathcal{X} \to \mathcal{Y}$ that minimizes an empirical risk objective $\frac{1}{N}\sum_{i=1}^N \mathcal{L}(f(x^{(i)}), y^{(i)})$ and generalizes accurately to unseen test distributions.
2. Unsupervised Learning: The algorithm receives an unannotated sample collection $\mathcal{D} = \{x^{(1)}, x^{(2)}, \dots, x^{(N)}\}$ consisting strictly of feature vectors without target labels $y^{(i)}$. Rather than predicting an external ground truth, the objective is to uncover underlying intrinsic structures, discover latent generative factors, estimate underlying probability density functions $P(X)$, identify topological manifolds, or partition the data into meaningful homogeneous subsets.
3. Semi-Supervised Learning: The algorithm operates over a small corpus of labeled instances $\mathcal{D}_L = \{(x^{(i)}, y^{(i)})\}_{i=1}^{N_L}$ augmented by a substantially larger volume of unannotated instances $\mathcal{D}_U = \{x^{(j)}\}_{j=1}^{N_U}$ (where $N_U \gg N_L$), leveraging structural smoothness and cluster assumptions over $\mathcal{D}_U$ to regularize decision boundaries.
4. Self-Supervised Learning: Supervisory signals are derived automatically from the intrinsic structure of unannotated data by formulating pretext tasks (e.g., masked token reconstruction, contrastive instance discrimination, next-frame video prediction), enabling representation learning without human annotations.

CODECopy Code
+-----------------------------------------------------------------------------------------+
|                               UNSUPERVISED LEARNING PARADIGM                            |
+-----------------------------------------------------------------------------------------+
|                               Input: Unlabeled Data X in R^d                            |
|                                                                                         |
|       +-----------------------+-------------------------+-----------------------+       |
|       |                       |                         |                       |       |
|       v                       v                         v                       v       |
|  [ Clustering ]     [ Dimensionality Reduction ]   [ Density Estimation ]  [ Anomaly Detect ]
|  Group data into    Project d dimensions into      Estimate underlying     Identify points  |
|  cohesive subsets   latent subspace k << d         probability P(X)        deviating from   |
|  (e.g., K-Means,    (e.g., PCA, t-SNE, UMAP,       (e.g., GMMs, KDE,       norm (e.g.,      |
|   DBSCAN, Hierarchical) Autoencoders)              Normalizing Flows)      iForest, LOF)   |
+-----------------------------------------------------------------------------------------+

1.2 Core Objectives of Unsupervised Learning

Unsupervised learning addresses distinct analytical challenges:

• Clustering: Discovering natural groupings within unlabeled data such that instances within the same group exhibit high internal similarity (cohesion) while instances across different groups exhibit high dissimilarity (separation).
• Dimensionality Reduction: Compressing high-dimensional feature spaces $\mathbb{R}^d$ into lower-dimensional manifolds $\mathbb{R}^k$ ($k \ll d$) while preserving maximal statistical variance (PCA), pairwise geometric distances (MDS), or local neighborhood topologies (t-SNE, UMAP).
• Density Estimation: Modeling the joint probability density function $P(x)$ underlying the observed samples, facilitating synthetic data sampling and probabilistic inference.
• Anomaly and Novelty Detection: Quantifying instance deviation from expected structural or density norms to flag outliers, sensor failures, and fraudulent transactions.
• Association Rule Mining: Uncovering transactional co-occurrence patterns (e.g., Apriori, FP-Growth) across multi-item relational databases.

1.3 Mathematical Formulation of the Clustering Problem

Let $\mathcal{X} = \{x^{(1)}, x^{(2)}, \dots, x^{(N)}\}$ represent a dataset of $N$ instances residing in a $d$-dimensional continuous vector space $\mathbb{R}^d$. The clustering problem seeks to partition $\mathcal{X}$ into $K$ non-empty, disjoint subsets (clusters) $\mathcal{C} = \{C_1, C_2, \dots, C_K\}$ satisfying:

1. Exhaustive Coverage: $\bigcup_{k=1}^K C_k = \mathcal{X}$ (every data point belongs to at least one cluster).
2. Mutual Exclusivity (Hard Clustering): $C_j \cap C_k = \emptyset$ for all $j \neq k$ (no data point belongs to more than one cluster).
3. Non-Triviality: $C_k \neq \emptyset$ for all $k \in \{1, \dots, K\}$.

Alternatively, in Soft (Fuzzy/Probabilistic) Clustering, each instance $x^{(i)}$ possesses a membership probability vector $w^{(i)} = [w_{i1}, w_{i2}, \dots, w_{iK}]^T \in [0, 1]^K$ such that $\sum_{k=1}^K w_{ik} = 1$, representing the degree of belonging to each respective cluster.

1.4 Taxonomy of Clustering Algorithms

Clustering methodologies differ fundamentally in how clusters are mathematically conceptualized:

• Partitional Clustering: Partitions data directly into $K$ non-overlapping subsets without hierarchical structure by optimizing an objective function (e.g., K-Means, K-Medoids / PAM, K-Modes).
• Hierarchical Clustering: Constructs nested tree structures (dendrograms) either agglomeratively (bottom-up merging of pairs) or divisively (top-down recursive splitting) based on linkage metrics (e.g., Single, Complete, Average, Ward's linkage).
• Density-Based Clustering: Defines clusters as continuous, dense spatial regions separated by low-density noise regions, capable of discovering arbitrary non-convex geometries (e.g., DBSCAN, OPTICS, HDBSCAN).
• Distribution-Based (Probabilistic) Clustering: Assumes data instances are generated from a mixture of underlying parametric probability distributions (e.g., Gaussian Mixture Models optimized via the Expectation-Maximization algorithm).
• Grid-Based Clustering: Quantizes continuous object space into a finite number of multi-dimensional cells forming a grid structure over which clustering operations execute independently of data size (e.g., STING, CLIQUE).
• Graph-Based and Spectral Clustering: Represents pairwise instance affinities as weighted graphs, performing dimensionality reduction via graph Laplacian graph spectra before partitional clustering.

2. K-Means Clustering: Mathematical Formulation & Lloyd's Algorithm

2.1 The Optimization Objective Function (Inertia / WCSS)

K-Means is a partitional clustering algorithm that assigns $N$ observations into $K$ clusters, where each cluster is parameterized by its mean vector (centroid) $\mu_k \in \mathbb{R}^d$.

Let $r_{ik} \in \{0, 1\}$ denote a binary indicator variable representing cluster membership:
$r_{ik} = \begin{cases} 1 & \text{if } x^{(i)} \in C_k \\ 0 & \text{otherwise} \end{cases}$
subject to the constraint $\sum_{k=1}^K r_{ik} = 1$ for all $i \in \{1, \dots, N\}$.

The objective of K-Means is to find cluster assignments $\{r_{ik}\}$ and centroid locations $\{\mu_k\}_{k=1}^K$ that minimize the Within-Cluster Sum of Squares (WCSS), also termed Inertia:
$J(r, \mu) = \sum_{i=1}^N \sum_{k=1}^K r_{ik} \|x^{(i)} - \mu_k\|_2^2 = \sum_{k=1}^K \sum_{i \in C_k} \sum_{j=1}^d (x_j^{(i)} - \mu_{kj})^2$
where $\|\cdot\|_2$ denotes the standard Euclidean norm ($L_2$ norm).

CODECopy Code
                      K-MEANS INERTIA (WCSS) GEOMETRY
                      
             Cluster C_1                      Cluster C_2
          x^(1)      x^(2)                x^(4)      x^(5)
             \        /                      \        /
              \  d_1 /                        \  d_4 /
               v    v                          v    v
               [ mu_1 ]                        [ mu_2 ]
               ^    ^                          ^    ^
              /  d_3 \                        /  d_6 \
             /        \                      /        \
          x^(3)      x^(N)                x^(6)      x^(7)
          
       WCSS = sum ||x^(i) - mu_1||^2  +  sum ||x^(j) - mu_2||^2

2.2 Coordinate Descent: Lloyd's Algorithm

Minimizing $J(r, \mu)$ simultaneously over discrete assignment variables $r$ and continuous coordinates $\mu$ is an NP-hard combinatorial optimization problem. Lloyd's algorithm solves this by decoupling the optimization via Block Coordinate Descent, alternating iteratively between two conditional minimization steps:

Step 1: Expectation / Cluster Assignment Step (E-Step)

Holding centroid positions $\{\mu_k\}_{k=1}^K$ fixed, minimize $J(r, \mu)$ with respect to cluster indicator variables $\{r_{ik}\}$. Because the inner sum over $k$ is unconstrained across different observations $i$, each instance $x^{(i)}$ is independently assigned to its closest centroid:
$r_{ik} = \begin{cases} 1 & \text{if } k = \arg\min_{j \in \{1, \dots, K\}} \|x^{(i)} - \mu_j\|_2^2 \\ 0 & \text{otherwise} \end{cases}$
Ties are broken arbitrarily. Geometrically, this step partitions the feature space into a Voronoi tessellation generated by the centroids $\{\mu_k\}$.

Step 2: Maximization / Centroid Update Step (M-Step)

Holding cluster assignments $\{r_{ik}\}$ fixed, minimize $J(r, \mu)$ with respect to centroid vectors $\mu_k$. Since $J$ is convex and unconstrained with respect to $\mu_k$, we compute the partial derivative with respect to $\mu_k$ and set it to zero:
$\frac{\partial J}{\partial \mu_k} = \frac{\partial}{\partial \mu_k} \sum_{i=1}^N r_{ik} \|x^{(i)} - \mu_k\|_2^2 = -2 \sum_{i=1}^N r_{ik} (x^{(i)} - \mu_k) = 0$
$\sum_{i=1}^N r_{ik} x^{(i)} - \mu_k \sum_{i=1}^N r_{ik} = 0 \implies \mu_k = \frac{\sum_{i=1}^N r_{ik} x^{(i)}}{\sum_{i=1}^N r_{ik}} = \frac{1}{|C_k|} \sum_{i \in C_k} x^{(i)}$
Thus, the optimal centroid $\mu_k$ is the empirical arithmetic mean of all data points currently assigned to cluster $C_k$.

CODECopy Code
+-----------------------------------------------------------------------------------------+
|                               LLOYD'S ALGORITHM WORKFLOW                                |
+-----------------------------------------------------------------------------------------+
|                                                                                         |
|      1. Initialize Centroids mu_1, mu_2, ..., mu_K in R^d                                |
|                                |                                                        |
|                                v                                                        |
|      +-------------------> [ Assignment Step ]                                          |
|      |                     Assign each x^(i) to nearest mu_k:                           |
|      |                     c^(i) = argmin_k ||x^(i) - mu_k||^2                          |
|      |                                |                                                 |
|      |                                v                                                 |
|      |                     [ Centroid Update Step ]                                     |
|      |                     Recalculate centroids as cluster means:                      |
|      |                     mu_k = (1 / |C_k|) * sum_{i in C_k} x^(i)                    |
|      |                                |                                                 |
|      |                                v                                                 |
|      |                     [ Convergence Check ]                                        |
|      |                     Did centroids shift by < tolerance epsilon,                  |
|      |                     or did assignments stabilize, or iter == max_iter?           |
|      |                                |                                                 |
|      |              No                |                Yes                              |
|      +--------------------------------+-----------------+                               |
|                                                         |                               |
|                                                         v                               |
|                                               [ Terminate & Output ]                    |
|                                               Final Clusters & Centroids                |
+-----------------------------------------------------------------------------------------+

2.3 Proof of Monotonic Convergence and Termination

1. Monotonic Non-Increasing Cost:
   • In Step 1, assigning each point to its closest centroid strictly decreases or maintains $J$: $J(r^{(t+1)}, \mu^{(t)}) \le J(r^{(t)}, \mu^{(t)})$.
   • In Step 2, setting each centroid to the cluster mean minimizes the sum of squared Euclidean distances, ensuring: $J(r^{(t+1)}, \mu^{(t+1)}) \le J(r^{(t+1)}, \mu^{(t)})$.
   • Therefore, at every iteration $t$:
     $J(r^{(t+1)}, \mu^{(t+1)}) \le J(r^{(t)}, \mu^{(t)})$
2. Finite Configuration Space:
   • For a dataset of $N$ points and $K$ clusters, there are at most $K^N$ distinct assignments of data points to clusters.
   • Because $J$ decreases strictly at every step unless assignments remain unchanged, no assignment partition can be revisited.
   • Consequently, Lloyd's algorithm is guaranteed to converge to a local minimum or stationary point in a finite number of iterations.
3. Non-Convexity and Local Minima:
   • While the algorithm is guaranteed to terminate, the objective function $J(r, \mu)$ is non-convex with numerous local minima, saddle points, and poor basin configurations.
   • The final partition reached depends heavily on initial centroid coordinates.
   • Finding the global minimum of WCSS is known to be NP-hard, even for $K=2$ in arbitrary dimensions, and for $d=2$ with arbitrary $K$.

3. Centroid Initialization Schemes & K-Means++

3.1 Standard Initialization Techniques and Their Vulnerabilities

Early formulations of K-Means relied on naive initialization heuristics:

• Forgy Method: Randomly selects $K$ distinct observations from dataset $\mathcal{X}$ to serve as initial centroids $\{\mu_1, \dots, \mu_K\}$.
• Random Partition Method: Randomly assigns each data point $x^{(i)}$ to one of the $K$ clusters with uniform probability, then computes the initial cluster means. This method places all initial centroids close to the global dataset center of mass, leading to poor partition initializations.

The Pathological Initialization Trap

When initial centroids are placed close together or clustered within a single dense cluster, K-Means can become trapped in suboptimal local minima. A cluster with high variance may be split into multiple sub-clusters, while two distinct distant clusters may be forced into a single centroid.

CODECopy Code
       PATHOLOGICAL INITIALIZATION                 DESIRED OPTIMAL INITIALIZATION
       
     True Cluster A     True Cluster B              True Cluster A     True Cluster B
      * * * * *          * * * * *                   * * * * *          * * * * *
      *  (X)(X) *        *         *                 *   (X)   *        *   (X)   *
      * * * * *          * * * * *                   * * * * *          * * * * *
     Two centroids       Zero centroids             One centroid       One centroid
     trapped in A!       assigned to B!

3.2 The K-Means++ Initialization Algorithm

To overcome initialization sensitivity, David Arthur and Sergei Vassilvitskii (2007) introduced K-Means++, which seeds initial centroids by selecting points spread out across the feature space using a distance-weighted probability distribution.

K-Means++ Algorithmic Protocol

1. First Centroid: Choose the first centroid $\mu_1$ uniformly at random from dataset $\mathcal{X}$.
2. Subsequent Centroids ($k = 2, \dots, K$):
   • For each data point $x^{(i)} \in \mathcal{X}$, compute its shortest squared Euclidean distance to any already selected centroid:
     $D(x^{(i)})^2 = \min_{j \in \{1, \dots, k-1\}} \|x^{(i)} - \mu_j\|_2^2$
   • Select the next centroid $\mu_k = x^{(m)}$ from $\mathcal{X}$ according to the categorical probability distribution:
     $P(x^{(i)}) = \frac{D(x^{(i)})^2}{\sum_{l=1}^N D(x^{(l)})^2}$
3. Repeat Step 2 until all $K$ centroids have been selected.
4. Proceed with standard Lloyd's coordinate descent iterations.

CODECopy Code
+-----------------------------------------------------------------------------------------+
|                              K-MEANS++ INITIALIZATION LOGIC                             |
+-----------------------------------------------------------------------------------------+
|                                                                                         |
|  1. Sample first centroid mu_1 ~ Uniform(X)                                             |
|                                                                                         |
|  2. For k = 2 to K:                                                                     |
|     a. For each x in X:                                                                 |
|        Calculate D(x) = min_{j < k} ||x - mu_j||                                        |
|                                                                                         |
|     b. Sample next centroid mu_k from distribution:                                     |
|        P(x) = D(x)^2 / sum_{x'} D(x')^2                                                 |
|                                                                                         |
|  Points far from existing centroids have high D(x)^2 -> high probability of selection.  |
|  Points close to existing centroids have low D(x)^2 -> low probability of selection.   |
+-----------------------------------------------------------------------------------------+

Theoretical Guarantee of K-Means++

Arthur and Vassilvitskii proved that K-Means++ initialization achieves an expected cost bound relative to the optimal clustering solution:
$\mathbb{E}[J_{\text{kmeans++}}] \le 8(\ln K + 2) J_{\text{opt}}$
This $O(\log K)$ competitive ratio ensures that K-Means++ is mathematically bounded from converging to catastrophically poor local minima, while simultaneously speeding up convergence by reducing the required number of coordinate descent iterations.

4. Determining the Optimal Number of Clusters ($K$)

Because K-Means requires the practitioner to specify $K$ a priori, heuristic diagnostic metrics are necessary to evaluate partitions across candidate values $K \in \{2, 3, \dots, K_{\max}\}$.

4.1 The Inertia Dilemma & The Elbow Method

As $K$ increases, the Within-Cluster Sum of Squares (Inertia) decreases monotonically:

• When $K = 1$, Inertia equals the total sum of squares around the global mean.
• When $K = N$, every data point is its own centroid, and Inertia reaches exactly zero ($J = 0$).

Therefore, minimizing Inertia directly cannot identify optimal $K$, as it trivially favors $K = N$.

The Elbow Heuristic

The Elbow Method plots Inertia $J(K)$ against candidate values of $K$:

• For $K < K^*$, adding a centroid captures major underlying cluster structure, producing a steep drop in Inertia ($\Delta J \gg 0$).
• For $K > K^*$, adding further centroids merely fragments existing clusters, producing diminishing marginal decreases in Inertia.
• The point of maximum curvature—resembling an "elbow"—indicates the parsimonious balance between model complexity and residual error.

CODECopy Code
     Inertia J(K)
      ^
      |  \
      |   \
      |    \   Underfitting Zone (Steep drop)
      |     \
      |      * <-- Optimal "Elbow" Point (K = 3)
      |       \
      |        \__________  Overfitting Zone (Diminishing returns)
      |                   \________
      +---------------------------------> K (Number of clusters)
         1   2   3   4   5   6   7

Mathematical Determination of the Elbow Knee

In complex datasets, the visual elbow may be subtle or ambiguous. Automated knee-detection algorithms (e.g., Satopaa et al.'s Kneedle algorithm) identify the point $K^*$ that maximizes the orthogonal distance between the discrete curve points $(K, J(K))$ and the secant line connecting $(K_{\min}, J(K_{\min}))$ to $(K_{\max}, J(K_{\max}))$.

4.2 Silhouette Analysis & Silhouette Coefficient

Silhouette Analysis provides an intrinsic validation metric that evaluates both cluster cohesion (tightness) and cluster separation (isolation) on an individual instance basis without requiring external ground truth.

Mathematical Formulation

For each observation $i \in \{1, \dots, N\}$ assigned to cluster $C_A$:

1. Mean Intra-Cluster Distance ($a(i)$):
   $a(i) = \frac{1}{|C_A| - 1} \sum_{j \in C_A, j \neq i} \|x^{(i)} - x^{(j)}\|_2$
   $a(i)$ quantifies how well instance $i$ is matched to its assigned cluster (lower indicates tighter cohesion).
2. Mean Nearest-Cluster Distance ($b(i)$):
   $b(i) = \min_{C_B \neq C_A} \frac{1}{|C_B|} \sum_{j \in C_B} \|x^{(i)} - x^{(j)}\|_2$
   $b(i)$ measures the distance between instance $i$ and its closest neighboring cluster (higher indicates greater separation).
3. Silhouette Coefficient ($s(i)$):
   $s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))}$
   Equivalently expressed as:
   $s(i) = \begin{cases} 1 - \frac{a(i)}{b(i)} & \text{if } a(i) < b(i) \\ 0 & \text{if } a(i) = b(i) \\ \frac{b(i)}{a(i)} - 1 & \text{if } a(i) > b(i) \end{cases}$

CODECopy Code
                          SILHOUETTE GEOMETRY
                          
          Assigned Cluster C_A                 Neighboring Cluster C_B
        +-----------------------+            +-------------------------+
        |                       |            |                         |
        |      x^(j)            |            |       x^(k)             |
        |        \              |            |        /                |
        |       a(i)            |            |       /                 |
        |          \            |            |      /                  |
        |           x^(i)-------|------------|---->x   b(i)            |
        |          /            |            |                         |
        |      x^(l)            |            |                         |
        +-----------------------+            +-------------------------+
             Intra-Cluster                       Nearest-Cluster
               Distance                             Distance

Silhouette Metric Range and Interpretation

The silhouette coefficient $s(i)$ is bounded in the interval $[-1, 1]$:

• $s(i) \approx +1$: The sample is far from neighboring clusters ($b(i) \gg a(i)$) and tightly bound within its own cluster (excellent assignment).
• $s(i) \approx 0$: The sample lies directly along the decision boundary between two clusters ($a(i) \approx b(i)$).
• $s(i) < 0$: The sample is closer to points in a neighboring cluster than to points in its assigned cluster ($a(i) > b(i)$), signaling an erroneous assignment.

The Global Silhouette Score is the mean silhouette coefficient across all $N$ dataset instances:
$S(K) = \frac{1}{N} \sum_{i=1}^N s(i)$
The optimal cluster count $K^*$ is identified by maximizing the global silhouette score: $K^* = \arg\max_{K} S(K)$.

Silhouette Diagram Interpretation

A silhouette plot visualizes sorted silhouette coefficients for each cluster:

• Cluster Balance: The relative vertical thickness of each cluster silhouette reflects the number of samples in that cluster.
• Cluster Quality: If a cluster's average silhouette score falls below the global average vertical dashed line, or contains substantial negative tails, that cluster is poorly separated or artificial.

5. Geometric Assumptions, Failure Modes, and Limitations

While computationally fast and intuitive, K-Means imposes strong structural assumptions on data distribution that can result in poor clustering when violated.

CODECopy Code
+-----------------------------------------------------------------------------------------+
|                                K-MEANS FAILURE MODES                                    |
+-----------------------------------------------------------------------------------------+
|                                                                                         |
|  1. Non-Convex / Ring Manifolds         2. Unequal Cluster Sizes / Densities            |
|       * * * * * * *                          * * * * * * * * *                          |
|     *   +-------+   *                       *      Dense      *         *  Sparse *     |
|    *   | Centroid|   *                      *    (X)  Cluster *        *    (X)    *    |
|    *   |  falls  |   *                      * * * * * * * * *           * * * * *       |
|     *  | in gap! |  *                       Large dense cluster        Centroid pulled  |
|       * * * * * * *                         split arbitrarily!         toward dense!    |
|                                                                                         |
|  3. Elongated / Anisotropic Clusters    4. Extreme Outlier Sensitivity                  |
|     * * * * * * * * * * * *                  * * * * *                                  |
|               (X)                           *   mu_1  *                * [ Outlier ]    |
|     * * * * * * * * * * * *                  * * * * *                      |           |
|     Cut across natural axis                  Centroid dragged toward outlier ------>    |
+-----------------------------------------------------------------------------------------+

5.1 Assumption of Spherical, Isotropic Clusters

K-Means implicitly assumes that:

1. Clusters are spherical (isotropic): Feature variance is identical in all coordinate directions.
2. Clusters have equal variance: The spatial spread $\sigma_k^2$ is uniform across all clusters $k \in \{1, \dots, K\}$.
3. Clusters have balanced sample sizes: All clusters contain roughly comparable observation counts.

This behavior stems directly from using the squared Euclidean distance metric $\|x - \mu\|_2^2$. In probabilistic terms, K-Means is mathematically equivalent to a Gaussian Mixture Model (GMM) where all components share an identical covariance matrix constrained to the spherical form $\Sigma_k = \sigma^2 I$.

5.2 Failure on Non-Convex and Arbitrary Geometries

Because K-Means partitions feature space via linear Hyperplane Voronoi boundaries, it cannot capture clusters with non-convex topologies, such as:

• Concentric circles or nested rings.
• Interlocking crescents (two-moons dataset).
• Elongated, serpentine, or sinusoidal manifolds.
  For non-convex geometries, density-based algorithms (DBSCAN) or spectral clustering are required.

5.3 Sensitivity to Feature Scaling

Because Euclidean distance computes squared coordinate differences $(x_j - \mu_j)^2$, features with larger numerical magnitudes dominate the distance calculations:

• If feature $x_1$ spans $[0, 100000]$ (e.g., Annual Income) and feature $x_2$ spans $[1, 5]$ (e.g., Credit Rating), distance calculations are dominated by $x_1$, rendering $x_2$ statistically irrelevant.
• Mandatory Preprocessing: Features must be standardized using StandardScaler ($z = \frac{x - \mu}{\sigma}$) or scaled using RobustScaler prior to fitting K-Means.

5.4 Sensitivity to Extreme Outliers

Because K-Means squares Euclidean distances, anomalous observations located far from the main data manifold exert large leverage on centroid calculations during the M-step. A single distant outlier can pull a centroid away from the true cluster center, sometimes isolating the outlier in its own degenerate cluster.

5.5 Handling Empty Clusters During Execution

During Lloyd's iterations, an initialization configuration may result in zero samples being assigned to a centroid $\mu_k$ ($|C_k| = 0$). Dividing by $|C_k| = 0$ during the centroid update produces an arithmetic exception. Standard scikit-learn implementations handle this by re-seeding the empty centroid at the observation currently exhibiting the largest squared distance $D(x)^2$ from its assigned centroid.

6. Scalability, Variations, and Alternative Formulations

6.1 Mini-Batch K-Means

To cluster web-scale datasets that exceed available RAM, D. Sculley (2010) proposed Mini-Batch K-Means. Rather than computing cluster assignments and updating centroids across the entire dataset $N$ in each iteration, Mini-Batch K-Means uses small, randomly sampled subsets (mini-batches) of size $B \ll N$ (typically $B \in [100, 1000]$).

Online Centroid Update Mechanism

For each mini-batch $\mathcal{B} \subset \mathcal{X}$ of size $B$:

1. Assign each instance $x^{(i)} \in \mathcal{B}$ to its nearest existing centroid: $c^{(i)} = \arg\min_k \|x^{(i)} - \mu_k\|_2^2$.
2. For each assigned instance $x$, update the associated centroid $\mu_k$ using a running stochastic gradient step:
   $v_k \leftarrow v_k + 1$
   $\eta = \frac{1}{v_k}$
   $\mu_k \leftarrow (1 - \eta)\mu_k + \eta x$
   where $v_k$ is a per-centroid counter tracking the cumulative number of samples assigned to centroid $k$ across all prior batches.

• Advantages: Drastically reduces disk I/O and RAM overhead, achieving an order-of-magnitude reduction in computation time with minimal degradation in clustering quality (typically within 1–3% of standard K-Means inertia).

6.2 K-Medoids (Partitioning Around Medoids - PAM)

In standard K-Means, the centroid $\mu_k$ is a virtual point calculated as the arithmetic mean, which may not correspond to any actual observation and is sensitive to outliers.

• K-Medoids restricts each cluster center to be an actual observation from the dataset (termed a medoid):
  $m_k \in \mathcal{X} \cap C_k$
• Objective Function: Minimizes the sum of dissimilarities using arbitrary distance metrics (e.g., Manhattan $L_1$, Cosine, or Jaccard distances):
  $J_{\text{medoids}} = \sum_{k=1}^K \sum_{x \in C_k} D(x, m_k)$
• Robustness: Because it uses medoids and $L_1$-style dissimilarities rather than squared Euclidean distances, K-Medoids is less sensitive to extreme outliers.

6.3 K-Modes and K-Prototypes for Categorical & Mixed Data

• K-Modes (Huang, 1998): Replaces the Euclidean distance metric with a simple matching dissimilarity measure:
  $D_{\text{cat}}(x, y) = \sum_{j=1}^d \delta(x_j, y_j) \quad \text{where } \delta(a, b) = \begin{cases} 0 & \text{if } a = b \\ 1 & \text{if } a \neq b \end{cases}$
  Centroids are represented by cluster modes (the most frequent categorical attribute values) rather than means.
• K-Prototypes: Combines K-Means and K-Modes using a weighted objective function:
  $D_{\text{proto}}(x, y) = \sum_{j \in \text{num}} (x_j - y_j)^2 + \gamma \sum_{j \in \text{cat}} \delta(x_j, y_j)$
  where $\gamma$ is a balance hyperparameter penalizing categorical mismatches.

7. Scikit-Learn Implementation Architecture

7.1 Production Clustering Pipeline

The following script demonstrates standard preprocessing, K-Means++ clustering, silhouette diagnostics, and elbow evaluation using Scikit-Learn:

pythonCopy Code
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans, MiniBatchKMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score, silhouette_samples
from sklearn.datasets import make_blobs

# 1. Generate Synthetic Data
X, y_true = make_blobs(
    n_samples=2000, 
    n_features=6, 
    centers=4, 
    cluster_std=1.2, 
    random_state=42
)

# 2. Mandatory Feature Scaling
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 3. Model Instantiation & Fitting
kmeans = KMeans(
    n_clusters=4,
    init='k-means++',        # K-Means++ seeding protocol
    n_init=10,               # 10 independent restarts with distinct seeds
    max_iter=300,            # Maximum coordinate descent iterations per run
    tol=1e-4,                # Convergence threshold for centroid shift
    algorithm='lloyd',       # Standard Lloyd coordinate descent
    random_state=42
)

labels = kmeans.fit_predict(X_scaled)
centroids = kmeans.cluster_centers_
inertia = kmeans.inertia_

print(f"Final WCSS (Inertia): {inertia:.2f}")

# 4. Silhouette Validation Metric
sil_score = silhouette_score(X_scaled, labels)
print(f"Mean Silhouette Coefficient: {sil_score:.4f}")

# 5. Mini-Batch K-Means for Streaming / Scalability
minibatch_kmeans = MiniBatchKMeans(
    n_clusters=4,
    init='k-means++',
    batch_size=256,
    max_iter=100,
    random_state=42
)
mb_labels = minibatch_kmeans.fit_predict(X_scaled)
print(f"Mini-Batch WCSS (Inertia): {minibatch_kmeans.inertia_:.2f}")

7.2 Parameter Reference Guide

• n_clusters (default=8): The number of clusters $K$ and centroids to generate.
• init ('k-means++', 'random', or array): Seeding method. 'k-means++' is default and recommended.
• n_init (default=10): Number of times the algorithm runs with distinct centroid seeds. The run with the lowest WCSS/Inertia is chosen.
• max_iter (default=300): Maximum iterations of the Lloyd coordinate descent per run.
• tol (default=1e-4): Relative tolerance with respect to Frobenius norm of difference in centroid positions to declare convergence.
• algorithm ('lloyd', 'elkan'): 'lloyd' uses standard coordinate descent; 'elkan' uses triangular inequality bounds to accelerate convergence on dense datasets.