K-Means Clustering

The Story That Explains K-Means

K-Means Clustering is an unsupervised machine learning algorithm that partitions a dataset into K distinct, non-overlapping clusters. Each data point belongs to the cluster whose centroid (centre point) is closest to it. The algorithm iterates until the cluster assignments stabilise.

The Algorithm — Step by Step

K-Means follows a beautifully simple loop. It takes five steps, and once you see them with a story, you will never forget them.

The Objective Function — What K-Means Minimises

K-Means is not guessing randomly — it is optimising a precise mathematical objective. Understanding this objective helps you understand both its power and its limitations.

K-Means++ — The Smarter Initialisation

Choosing K — The Hardest Part

K-Means requires you to specify K upfront. This is both its biggest weakness and the source of the most common beginner mistake (just picking K=3 because it sounds right). Here are the two most reliable methods.

When K-Means Works — and When It Fails

Feature Scaling — Non-Negotiable for K-Means

Python Code — K-Means from Scratch to Production

Let us build up from scratch. First, a minimal working example. Then the full production pattern with elbow method, silhouette scores, and visualisation.

8.1 — Minimal Example (Synthetic Data)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.preprocessing import StandardScaler

# 1. Generate synthetic data — 300 points in 4 natural clusters
X, y_true = make_blobs(n_samples=300, centers=4,
                        cluster_std=0.70, random_state=42)

# 2. Scale features — ALWAYS do this before K-Means
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 3. Fit K-Means with K=4 (K-Means++ initialisation is the default)
kmeans = KMeans(n_clusters=4, init='k-means++',
                n_init=10, random_state=42)
kmeans.fit(X_scaled)

# 4. Results
labels     = kmeans.labels_          # cluster ID for each point
centroids  = kmeans.cluster_centers_ # centroid coordinates
inertia    = kmeans.inertia_         # WCSS (lower = tighter clusters)
iterations = kmeans.n_iter_          # how many EM iterations ran

print(f"WCSS (Inertia):  {inertia:.2f}")
print(f"Iterations:      {iterations}")
print(f"Cluster counts:  {np.bincount(labels)}")

8.2 — Finding the Right K: Elbow + Silhouette

from sklearn.metrics import silhouette_score

wcss_list   = []
sil_list    = []
K_range     = range(2, 11)

for k in K_range:
    km = KMeans(n_clusters=k, init='k-means++',
                 n_init=10, random_state=42)
    km.fit(X_scaled)
    wcss_list.append(km.inertia_)
    sil_list.append(silhouette_score(X_scaled, km.labels_))

# Print both metrics side by side
print(f"{'K':>3}  {'WCSS':>10}  {'Silhouette':>12}")
print("-" * 30)
for k, w, s in zip(K_range, wcss_list, sil_list):
    print(f"{k:>3}  {w:>10.2f}  {s:>12.4f}")

best_k = K_range[sil_list.index(max(sil_list))]
print(f"\nBest K by Silhouette: {best_k}")

📌 Diagram — The K-Means Algorithm Visualised

The diagram below walks through four iterations of K-Means on a 2D dataset. Notice how the centroids (★) drift toward the true cluster centres with each step.

📋 Case Study 1 — Customer Segmentation for an E-Commerce Platform

Step 1 — Load & Explore the Data

import pandas as pd
import numpy as np
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score

# Simulate 50,000 customer records
np.random.seed(42)
n = 50000

df = pd.DataFrame({
    'customer_id':      range(1, n+1),
    'total_spend':      np.random.gamma(2, 500, n).round(2),
    'purchase_freq':    np.random.poisson(8, n),
    'avg_order_value':  np.random.gamma(3, 40, n).round(2),
})

print(df.describe().round(2))

# Select features and scale
features = ['total_spend', 'purchase_freq', 'avg_order_value']
X = df[features].values
scaler = StandardScaler()
X_sc   = scaler.fit_transform(X)

# Evaluate K=2 to K=8
results = []
for k in range(2, 9):
    km = KMeans(n_clusters=k, init='k-means++',
                 n_init=10, random_state=42)
    km.fit(X_sc)
    sil = silhouette_score(X_sc, km.labels_, sample_size=5000)
    results.append({'K': k, 'WCSS': round(km.inertia_, 1), 'Silhouette': round(sil, 4)})

print(pd.DataFrame(results).to_string(index=False))

# Fit final model
final_km = KMeans(n_clusters=4, init='k-means++',
                   n_init=20, random_state=42)
df['segment'] = final_km.fit_predict(X_sc)

# Profile each segment (un-scaled for interpretability)
profile = df.groupby('segment')[features].mean().round(2)
profile['count']   = df.groupby('segment').size()
profile['% share'] = (profile['count'] / n * 100).round(1)
print(profile)

import numpy as np
from sklearn.cluster import KMeans
from PIL import Image
import matplotlib.pyplot as plt

# Load image
img   = np.array(Image.open('beach_sunset.jpg'))
h, w, c = img.shape
print(f"Image shape: {h}×{w} pixels, {c} colour channels")
print(f"Total unique colours (theoretical): {256**3:,}")
print(f"Total unique colours (actual):       {len(set(map(tuple, img.reshape(-1,3)))):,}")

# Reshape: from (H, W, 3) to (H*W, 3)
# Each row is now one pixel with [R, G, B]
pixels = img.reshape(-1, 3).astype('float32')
print(f"Pixel matrix shape: {pixels.shape}")

# Fit K-Means: find 64 representative colours
# Note: RGB values 0–255 are already on the same scale, so no scaling needed
K = 64
km = KMeans(n_clusters=K, init='k-means++',
             n_init=3, max_iter=100,
             random_state=42, verbose=0)

# Use a random sample for speed, then predict all pixels
sample_idx = np.random.choice(len(pixels), size=50000, replace=False)
km.fit(pixels[sample_idx])

# Assign every pixel to nearest centroid colour
labels = km.predict(pixels)

# Reconstruct: replace each pixel with its centroid colour
palette     = km.cluster_centers_.astype('uint8')  # 64×3 colour palette
compressed  = palette[labels].reshape(h, w, 3)      # back to (H, W, 3)

print(f"Original unique colours:    412,881")
print(f"Compressed unique colours:  {K}")
print(f"Original storage (raw):     {h*w*3:,} bytes")
print(f"Compressed storage:         {h*w*np.ceil(np.log2(K))/8:,.0f} bytes")
print(f"Compression ratio:          {h*w*3 / (h*w*np.ceil(np.log2(K))/8):.1f}×")

# Compare visual quality at different K values
for k_val in [4, 8, 16, 32, 64, 128]:
    km_temp = KMeans(n_clusters=k_val, init='k-means++',
                      n_init=3, max_iter=100, random_state=42)
    km_temp.fit(pixels[sample_idx])
    lbl_temp  = km_temp.predict(pixels)
    pal_temp  = km_temp.cluster_centers_.astype('uint8')
    bits_px   = np.ceil(np.log2(k_val))
    bytes_tot = int(h * w * bits_px / 8)
    ratio     = h * w * 3 / bytes_tot
    print(f"K={k_val:4d} | Bits/pixel: {bits_px:.0f} | "
          f"Size: {bytes_tot:>9,} bytes | Ratio: {ratio:.1f}×")

from sklearn.cluster import MiniBatchKMeans
import time

# Simulate 1 million records
X_big = np.random.randn(1_000_000, 5)

# Standard K-Means (sample 50k for timing fairness)
t0 = time.time()
km_std = KMeans(n_clusters=8, n_init=3, random_state=42)
km_std.fit(X_big[:50000])
print(f"Standard K-Means (50k rows):    {time.time()-t0:.2f}s  WCSS={km_std.inertia_:.1f}")

# Mini-Batch K-Means (full 1M rows)
t1 = time.time()
km_mb = MiniBatchKMeans(n_clusters=8, batch_size=2048,
                           n_init=3, random_state=42)
km_mb.fit(X_big)
print(f"Mini-Batch K-Means (1M rows):   {time.time()-t1:.2f}s  WCSS={km_mb.inertia_:.1f}")