BACK TO BYTELOGIC PLATFORMDOSSIER: CONCEPT // 001

MACHINE LEARNING · CLUSTERING•INTERMEDIATE•LLOYD'S ALGORITHM

K-Means Clustering

A visual and mathematical explanation of clustering through iterative expectation-maximization, Voronoi partitioning, and centroid displacement.

03 Derivation 05 Interactive Lab 06 NumPy Implementation

01 / THE FORMAL PROBLEM

Partitioning Continuous Feature Space

Given an unlabeled dataset of N observations X = {x_1, x_2, ..., x_N} where each observation x_i ∈ ℝ^d is a d-dimensional continuous vector, our goal is to partition the N observations into k non-empty, mutually disjoint subsets S = {S_1, S_2, ..., S_k} such that:

Disjoint Partition Requirement[PARTITION]

\bigcup_{j=1}^k S_j = X \quad \text{and} \quad S_a \cap S_b = \emptyset \quad \forall a \neq b

Finding the optimal partition that globally minimizes intra-cluster distance is an NP-hard combinatorial optimization problem even for k = 2 in general dimension. Lloyd's algorithm offers a deterministic heuristic that converges monotonically to a local optimum.

02 / GEOMETRIC INTUITION

Centers of Gravity and Gravitational Pull

Imagine placing $k$ anchors onto a plane covered with particles. Each particle feels an allegiance to whichever anchor is closest to it, forming distinct territories known as Voronoi cells.

Once every particle has declared allegiance to its nearest anchor, each anchor relocates to the exact center of gravity (the mathematical mean vector $\mu_j$) of all its devoted particles.

Because the anchors moved, the territories shift. Particles on the borders might now find a different anchor closer. We repeat this dance until no particle changes allegiance and the anchors cease to move.

03 / MATHEMATICAL DERIVATION

Objective Function & Coordinate Descent

The optimization objective minimizes the Within-Cluster Sum of Squares (WCSS), also referred to as inertia:

Within-Cluster Sum of Squares (WCSS)[EQ 01]

J(S, \mu) = \sum_{j=1}^k \sum_{x_i \in S_j} \|x_i - \mu_j\|^2

Where x_i is the observation vector, S_j is the j-th cluster subset, and \mu_j is the center coordinate vector of cluster j.

J:Total objective inertia loss

k:Number of clusters

S_j:Set of observations assigned to cluster j

x_i:Single d-dimensional data point

\mu_j:Centroid mean vector of cluster j

\|\cdot\|^2:Squared L2 Euclidean norm

The Coordinate Descent Decomposition

Lloyd's algorithm solves this non-convex problem by alternating minimization across two sets of variables: the discrete cluster assignments $S$ and the continuous centroid coordinates $\mu$.

Step 1: Voronoi Partition Assignment[EQ 02]

S_j^{(t)} = \left\{ x_i \in X : \|x_i - \mu_j^{(t)}\|^2 \le \|x_i - \mu_l^{(t)}\|^2 \quad \forall l = 1, \dots, k \right\}

Holding centroids fixed, we assign each point x_i to the closest centroid \mu_j, which strictly minimizes J with respect to S.

Step 2: Analytical Centroid Relocation[EQ 03]

\frac{\partial J}{\partial \mu_j} = -2 \sum_{x_i \in S_j} (x_i - \mu_j) = 0 \implies \mu_j^{(t+1)} = \frac{1}{|S_j^{(t)}|} \sum_{x_i \in S_j^{(t)}} x_i

Holding assignments fixed, taking the derivative with respect to \mu_j and setting to zero yields the arithmetic mean of the assigned points.

THEOREM: MONOTONIC CONVERGENCE

Because both Step 1 and Step 2 strictly decrease or preserve $J(S, \mu)$, and because the number of distinct partitions of $N$ points into $k$ subsets is finite (bounded by $k^N$), Lloyd's algorithm cannot cycle and must terminate at a local minimum in a finite number of iterations.

04 / ALGORITHM FLOW

The 5-Step Execution Cycle

01PHASE

INITIALIZE

Select k initial centroids (Random or K-Means++ D²)

→

02PHASE

DISTANCES

Compute pairwise L2 distances from all N points to k centroids

→

03PHASE

ASSIGN

Assign each point to its closest centroid (argmin Euclidean distance)

→

04PHASE

RECOMPUTE

Relocate each centroid to the mean of its assigned cluster

→

05PHASE

CONVERGE

Check if centroid shift ||Δμ|| < ε. If not, repeat from 02

05 / INTERACTIVE LABORATORY VISUAL

LIVE 2D CLUSTERING MANIFOLD

Observe Centroids Traversing the Manifold

Step through the alternating minimization phases. Watch the Voronoi partitioning update in real time as data points get captured by incoming centroids.

STATE:READY|ITER: 0

LOSS:0

K Clusters:

06 / IMPLEMENTATION FROM SCRATCH

Vectorized NumPy Lloyd Algorithm

Production-grade Python code implementing vectorized pairwise distance broadcasting: $(N, 1, d) - (1, k, d) \to (N, k, d)$ without external machine learning dependencies.

06 /K-Means from Scratchpython

GitHub

01import numpy as np
02
03class KMeans:
04    """
05    K-Means Clustering via Lloyd's Algorithm.
06    Vectorized implementation in pure NumPy.
07    """
08    def __init__(self, k=3, max_iter=300, tol=1e-4, init='kmeans++'):
09        self.k = k
10        self.max_iter = max_iter
11        self.tol = tol
12        self.init = init
13        self.centroids = None
14        self.inertia_ = None
15
16    def _init_centroids(self, X):
17        n_samples = X.shape[0]
18        if self.init == 'random':
19            indices = np.random.choice(n_samples, self.k, replace=False)
20            return X[indices].copy()
21        
22        # K-Means++ D^2 Initialization
23        centroids = [X[np.random.choice(n_samples)]]
24        for _ in range(1, self.k):
25            # Compute distance to closest centroid
26            dists = np.min([np.sum((X - c)**2, axis=1) for c in centroids], axis=0)
27            probs = dists / np.sum(dists)
28            next_centroid_idx = np.random.choice(n_samples, p=probs)
29            centroids.append(X[next_centroid_idx])
30        return np.array(centroids)
31
32    def fit(self, X):
33        X = np.asarray(X, dtype=np.float64)
34        n_samples = X.shape[0]
35        self.centroids = self._init_centroids(X)
36
37        for iteration in range(self.max_iter):
38            # Step A: Vectorized pairwise Euclidean distance calculation
39            # Shape: (n_samples, k)
40            distances = np.linalg.norm(X[:, np.newaxis, :] - self.centroids[np.newaxis, :, :], axis=2)
41            
42            # Step B: Expectation step - assign to closest centroid
43            labels = np.argmin(distances, axis=1)
44            
45            # Step C: Maximization step - update centroids to cluster means
46            new_centroids = np.zeros_like(self.centroids)
47            for j in range(self.k):
48                mask = (labels == j)
49                if np.any(mask):
50                    new_centroids[j] = np.mean(X[mask], axis=0)
51                else:
52                    # Handle empty cluster: re-initialize with random sample
53                    new_centroids[j] = X[np.random.choice(n_samples)]
54
55            # Check convergence via tolerance
56            shift = np.linalg.norm(new_centroids - self.centroids)
57            self.centroids = new_centroids
58            if shift < self.tol:
59                break
60
61        # Compute final inertia (Within-Cluster Sum of Squares)
62        final_dists = np.min([np.sum((X - c)**2, axis=1) for c in self.centroids], axis=0)
63        self.inertia_ = np.sum(final_dists)
64        return self

07 / EMPIRICAL EXPERIMENT

Random Initialization vs. K-Means++

Standard Lloyd initialization randomly selects $k$ observations uniformly, which frequently places two centroids within the same true cluster. Arthur & Vassilvitskii (2007) introduced K-Means++, choosing subsequent centroids with probability proportional to their squared distance $D(x)^2$ from already chosen centroids, guaranteeing an $O(\log k)$ competitive ratio.

METHOD	APPROXIMATION RATIO	AVG RUNTIME ITERATIONS	LOCAL MINIMA SENSITIVITY
K-Means++ (Arthur & Vassilvitskii)	O(log k) Guaranteed	7.8 steps	Very Low (Optimal dispersion)
Uniform Random Initialization	Unbounded (arbitrarily bad)	18.4 steps	High (Frequent local traps)

08 / LIMITATIONS & PATHOLOGIES

When K-Means Fails

Non-Spherical Geometry

Because distance is measured with isotropic Euclidean norms, K-Means assumes convex, spherical clusters. It completely fails on concentric circles, crescent moons, or manifold ribbons (DBSCAN or Spectral Clustering are required).

Varying Cluster Densities & Sizes

If one cluster contains 10,000 points and a neighboring cluster contains 100 points, K-Means will split the large cluster in half and merge the small cluster into the neighbor to minimize squared distance.

Scale Sensitivity

Features with large numerical variances dominate the squared distance computation. Features must strictly be standardized ($\mu = 0, \sigma = 1$) prior to clustering.

Sensitivity to Extreme Outliers

Because distances are squared in the objective function, a single rogue point far from the origin will drag a centroid away from legitimate data (K-Medoids / PAM provides an L1 robust alternative).

09 / CONNECTED KNOWLEDGE GRAPH

Concepts Orbiting K-Means

In ByteLogic, no concept lives in isolation. Explore the theoretical connections from K-Means to general latent variable models:

Generalization

Gaussian Mixture Models (GMM)

Soft probabilistic cluster assignments using Expectation-Maximization with full covariance matrices.

Explore Concept

Continuous Relaxation

Principal Component Analysis (PCA)

Ding & He (2004) proved that PCA subspaces are continuous relaxations of K-Means cluster indicator vectors.

Explore Concept

Geometric Dual

Voronoi Tessellation & Delaunay

Partitioning space into polygonal convex regions closest to a discrete set of seed sites.

Explore Concept