K-Means Clustering

Explore unsupervised learning with K-means, centroid-based clustering, and elbow method

beginner45 min
On this page

The Big Picture

Unlike supervised learning (where we have labels), clustering is unsupervised: we have just features and want to discover natural groups in the data. KMeans is the simplest and most popular clustering algorithm. It partitions data into K clusters by minimizing within-cluster variance — think of it as "placing K centers so that points cluster tightly around them."

From Raw Data → KMeans Clusters
RAW DATA No labels, just points KMeans with K=3 3 CLUSTERS C1 C2 C3 Each point assigned to nearest center
Key insight: KMeans minimizes the within-cluster sum of squares (WCSS) — the sum of distances from each point to its cluster center. This makes clusters as "tight" as possible.

When is clustering useful? When you want to discover structure without pre-existing labels. Examples include customer segmentation, document clustering, image compression, and anomaly detection. It's called "unsupervised" because the algorithm learns the clusters from data alone.

AdTech Example: Audience Segmentation for Bing Ads

At Microsoft Bing Ads, you have millions of users with behavioral data. You want to segment them into distinct audience groups to tailor ad strategies per group. Unlike supervised learning (where you'd predict a specific outcome), you don't have pre-labeled "audience segments" — you discover them via clustering.

User Features → KMeans Clusters → Personalized Campaigns
USER FEATURES • Search frequency • Avg time per session • Device preference • Category affinity • Purchase history • Click-through rate • Time zone • Competitor clicks • Seasonal patterns 10M users, 9 features KMEANS K=5 clusters Learn 5 centers Assign each user to nearest center Output: cluster_id ∈ {1..5} per user CAMPAIGN STRATEGIES Cluster 1: Power Users → Premium, high-bid ads Cluster 2: Mobile First → Mobile-optimized ads Cluster 3: Price Sensitive → Deals & discount ads Clusters 4 & 5: ... → Tailored bidding & creative
AdTech benefit: Once you have clusters, you can bid differently per cluster. Power users → higher bids. Price-sensitive users → lower bids or discount-focused ads. This audience segmentation is a core strategy in ad platforms.

How KMeans Works — Step-by-Step Algorithm

KMeans is simple: initialize K centers, then repeatedly (1) assign points to nearest center, (2) recompute centers. Repeat until convergence.

KMeans Algorithm Flow
STEP 1: Initialize K Centers Randomly select K points from the dataset (or use KMeans++ for better initialization). centers ← random_k_points(data) STEP 2: Assign Points to Nearest Center (E-Step) For each point, compute distance to all K centers. Assign to the closest one. cluster[i] ← argmin_k distance(point[i], center[k]) STEP 3: Recompute Centers (M-Step) For each cluster, compute the mean of all points in that cluster as the new center. center[k] ← mean(points in cluster k) STEP 4: Check Convergence If centers haven't moved (or moved less than ε), STOP — clustering is done. Else go back to STEP 2 and repeat. if max(distance(old_center, new_center)) < ε: done
Why two steps? E-step (Expectation) assigns points. M-step (Maximization) optimizes centers. This is the Expectation-Maximization (EM) framework in action. Each iteration decreases the total within-cluster variance, so the algorithm is guaranteed to converge.
KMeans Convergence: Visual Example (K=2)
Iteration 0 Random centers Iteration 1 Centers recomputed Iteration 2 CONVERGED! Centers stable

Distance Measures (Euclidean, Manhattan, Cosine, Mahalanobis)

The choice of distance metric fundamentally changes KMeans. Different metrics suit different data types and problem domains.

1. Euclidean Distance (L2)

Euclidean Distance: Straight-Line Distance in 2D/3D Space
EUCLIDEAN (L2) d(p,q) = √(Σᵢ (pᵢ − qᵢ)²) Most common; assumes continuous numeric data Example: 2D Space x y p = (2,3) q = (5,6) Δx=3 Δy=3 d = √(3²+3²) ≈ 4.24

2. Manhattan Distance (L1)

Manhattan Distance: City Block Distance
MANHATTAN (L1) d(p,q) = Σᵢ |pᵢ − qᵢ| Sum of absolute differences; robust to outliers Example: Grid Movement x y p = (2,3) q = (5,6) 3 steps 3 steps d = |3| + |3| = 6

3. Cosine Similarity / Distance

Cosine Distance: Angle-Based, Ignores Magnitude
COSINE DISTANCE cos_sim(p,q) = (p·q) / (||p|| ||q||) d(p,q) = 1 − cos_sim(p,q) Good for text, sparse vectors; ignores magnitude Example: Two Vectors p q θ (angle) Small θ → vectors similar cos_sim → 1, distance → 0

4. Mahalanobis Distance

Accounts for correlation between features and feature variance. More sophisticated than Euclidean but computationally expensive.

Mahalanobis Distance: Accounting for Covariance
MAHALANOBIS DISTANCE d(p,q) = √((p−q)ᵀ Σ⁻¹ (p−q)) where Σ = covariance matrix of data • Handles correlated features (e.g., height & weight) • More expensive to compute; requires inverting covariance matrix • Rarely used in practice for KMeans due to cost; better for anomaly detection

When to Use Which Distance Metric?

Metric Best For Characteristics Cost
Euclidean Continuous numeric data (user age, revenue, time spent). Most common default. Geometric distance; assumes isotropic (equal variance) clusters O(d) per pair
Manhattan Data on grid/network (lat/lon, chess moves). Robust to outliers. Sum of absolute differences; less sensitive to extreme values O(d) per pair
Cosine Text, TF-IDF vectors, sparse high-dimensional data (document clustering) Angle-based; ignores magnitude. Symmetric. O(d) per pair
Mahalanobis When features are correlated; when you want to account for variance Statistically principled; handles covariance O(d³) initial + O(d²) per pair
AdTech note: For user feature clustering, Euclidean distance is standard. Features like CTR, spend, session time are roughly Gaussian and independent. If you have text features (e.g., search queries), use cosine distance after vectorization (TF-IDF or embeddings).

Handling Mixed Data Types (K-Prototypes, Gower Distance)

Real-world data often has both numeric (age, revenue) and categorical (device type, region) features. Standard KMeans (Euclidean) can't handle categorical data well. Here are three strategies.

Strategy 1: One-Hot Encoding + Euclidean

Convert Categorical to Binary via One-Hot Encoding
ORIGINAL DATA User 1: age = 28 (numeric) device = "mobile" (categorical) region = "US_West" (categorical) User 2: age = 35 (numeric) device = "desktop" (categorical) region = "EU" (categorical) ONE-HOT ENCODED User 1: [28, 1, 0, 1, 0] ↑age ↑mobile ↑desktop ↑US_West ↑EU • device="mobile" → [1, 0] • device="desktop" → [0, 1] • region="US_West" → [1, 0] • region="EU" → [0, 1] Now use standard Euclidean KMeans!
Pros & Cons:
✓ Simple, uses standard KMeans.
✓ Works well when categorical features have few categories.
✗ Creates sparse, high-dimensional vectors (curse of dimensionality).
✗ Treats all categorical levels equally in distance (no semantic similarity).

Strategy 2: K-Prototypes (Numeric + Categorical)

K-Prototypes extends KMeans to handle both data types simultaneously by using different distance functions per type.

K-Prototypes Algorithm: Mixed Distance Metric
DISTANCE IN K-PROTOTYPES d(x,c) = Σⱼ(xⱼ − cⱼ)² + λ·Σₖ δ(xₖ,cₖ) Mix numeric (squared difference) + categorical (Hamming distance) with weight λ NUMERIC PART Features: age, ctr, spend, session_time d_numeric = √(Σ (xᵢ − cᵢ)²) Standard Euclidean distance on continuous features CATEGORICAL PART Features: device, region, gender d_cat = Σ δ(xₖ,cₖ) where δ = 0 if values match, 1 if different (Hamming distance) λ controls influence
K-Prototypes advantage: No need to one-hot encode. Prototypes (centers) store actual values: age=31.2, device="mobile", region="US_West". Much more interpretable than binary vectors!

Strategy 3: Gower Distance (Unified Metric)

Gower distance normalizes all features to [0,1] range, then computes a weighted average of per-feature distances. Works for any data type.

Gower Distance: Normalized per-Feature, Averaged
GOWER DISTANCE FORMULA d_gower(x,y) = (1/p) · Σⱼ dⱼ(xⱼ,yⱼ) where for each feature j: Numeric: dⱼ = |xⱼ − yⱼ| / range(j) (Normalize by max−min of that feature) Categorical: dⱼ = 0 if xⱼ = yⱼ, else 1 (Match or mismatch) Result: distance in [0,1], comparable across feature types
Gower vs K-Prototypes:
Gower: More flexible, works with any feature type. Distance between points, symmetric.
K-Prototypes: Simpler to implement, faster, prototypes are interpretable feature vectors.

Comparison: Which Strategy?

Strategy When to Use Pros Cons
One-Hot + Euclidean Few categorical features with few categories. Simple baselines. Standard, well-understood. Most tools support it. High dimensionality. Loses categorical semantics. Doesn't scale with many categories.
K-Prototypes Mixed numeric + categorical data. Want interpretable prototypes. Handles both types natively. Prototypes are human-readable. Less common; fewer libraries. Need to tune λ parameter.
Gower Distance Many feature types. Need unified distance metric. Robust approach. Principled. Handles ordinal, binary, continuous, categorical. Normalized [0,1]. More complex. Slower. Still need to pair with some clustering algorithm.

Choosing K — Elbow Method, Silhouette Score, Gap Statistic

The number of clusters K is not known a priori in unsupervised learning. Here are three common methods to determine it.

1. Elbow Method

Compute WCSS (within-cluster sum of squares) for different K values. Plot WCSS vs K. Look for the "elbow" — the point where adding more clusters no longer significantly reduces WCSS.

Elbow Method: WCSS Curve
WCSS K (clusters) 8000 0 1 2 3 4 5 6 7 ← Elbow at K=4 Big drop before, small drop after Steep drop Marginal gain
Elbow method pros: Simple, fast, intuitive.
Cons: Elbows aren't always clear (especially in noisy data). Subjective.

2. Silhouette Score

For each point, measure how similar it is to points in its own cluster vs. other clusters. Silhouette score ranges from -1 (bad, in wrong cluster) to +1 (good, far from other clusters).

Silhouette Coefficient: Cohesion vs. Separation
SILHOUETTE COEFFICIENT s(i) = (b(i) − a(i)) / max(a(i), b(i)) a(i) = mean dist to points in own cluster | b(i) = min mean dist to other cluster Sample Silhouettes (K=3) Point in Cluster 1: 0.78 Point on boundary: 0.25 Misclassified: -0.12 Avg Silhouette vs K K=2 K=3 K=4 K=5 K=6 Pick K=3 (highest average)
Silhouette pros: Quantitative, considers both cohesion & separation.
Cons: Computationally expensive O(n²) for large datasets.

3. Gap Statistic

Compares WCSS of real data to WCSS of random uniform data. A large gap means the real clusters are much tighter than random noise, indicating good clustering.

Gap Statistic: Real Data vs. Random Reference
GAP STATISTIC Gap(k) = E[log WCSS_random(k)] − log WCSS_real(k) Choose K where Gap(k) ≥ Gap(k+1) − s(k+1) (most significant gap) Log WCSS K Random ref Real data Gap Largest gap suggests K=3 or 4

Which Method to Use?

Method Speed Reliability When to Use
Elbow Fast O(nk) Moderate (subjective) Quick exploration. Good starting point. Use when you want simplicity.
Silhouette Slow O(n²) Good (quantitative) When you want quantitative quality metric. Works well for small-medium datasets.
Gap Statistic Very slow (needs random samples) Very good (statistically principled) When you need statistical rigor. Slower but most robust. For research/critical applications.
AdTech recommendation: Start with Elbow method to get a rough sense. Then compute Silhouette scores for those candidate K values. If dataset is large, use Silhouette on a sample. Domain knowledge (e.g., "we want ~5 audience segments") should also inform K.

Initialization Strategies — Random vs KMeans++

KMeans is sensitive to initial center placement. Bad initialization can trap the algorithm in a poor local optimum. KMeans++ initialization dramatically improves convergence quality.

Random Initialization (Naive)

Random Init: Can Get Stuck in Poor Local Optima
BAD RANDOM INIT C1 C2 C3 Centers not near actual clusters GOOD KMEANS++ INIT C1 C2 C3 Centers near actual clusters

KMeans++ Initialization

Smart initialization: place first center randomly, then place each subsequent center far from existing centers. This spreads centers out, avoiding poor local optima.

KMeans++ Algorithm: Spread Out Centers
STEP 1: Pick First Center Randomly center₁ ← random_point(data) Each point has equal probability (1/n). STEP 2: Pick 2nd Center (Far from 1st) For each point, compute d(point, center₁) Pick next center with probability ∝ d² STEP 3: Pick 3rd, 4th, ... K-th Centers (Iteratively Far) For each remaining point, compute d(point, nearest_existing_center) Pick next center with probability ∝ d² (far from all existing centers) Why d² Weighting? Squaring distance makes far points much more likely. Point at distance 3 is 9x more likely than point at distance 1. Result: centers spread across data, covering multiple clusters from the start.
Impact of KMeans++: Same algorithm (E-step, M-step), but much better starting points.
• Better final WCSS (tighter clusters).
• Faster convergence (fewer iterations).
• Less dependent on random seed (more consistent).
Cost: O(nk) additional initialization (negligible for training).

Recommendation

Always use KMeans++ for production. It's the default in scikit-learn, nearly cost-free, and dramatically improves robustness. If you have domain knowledge about centers (e.g., "we know 5 cities matter"), you can seed with those instead.

Convergence, Local Minima & Practical Tips

Convergence Guarantee

KMeans is guaranteed to converge (reach a fixed point where assignments don't change). However, it may converge to a local optimum, not the global optimum. This is why good initialization (KMeans++) matters.

Local vs. Global Optima in KMeans
WCSS (cost) Global Min (best clustering) Local Min 1 (suboptimal) Local Min 2 (suboptimal) Depends on where KMeans starts! KMeans++ reduces this risk.
Practical mitigation:
1. Run multiple times: KMeans 10+ times with different random seeds, pick the best (lowest WCSS).
2. Use KMeans++: Reduces local optima issues significantly.
3. Set convergence tolerance: Don't require perfect convergence; stop when centers move < ε.

Practical Tuning Tips

Convergence Settings

max_iter: 300 (default)
tol: 1e-4 (center movement threshold)
n_init: 10 (number of runs)
Run 10× with different seeds, keep best.

For Large Datasets

Mini-batch KMeans: Subset data per iteration (faster, lower WCSS quality)
Sampling: Cluster sample, apply to full data
Sparse matrices: If data is sparse, use sparse distance metrics

When KMeans Struggles

Challenge Why Solution
Unequal cluster sizes KMeans minimizes WCSS, not cluster size. Tends to split large clusters into many small ones. Use DBSCAN or hierarchical clustering. Or apply K-means to balanced subsamples.
Non-convex clusters KMeans assumes spherical (convex) clusters. Crescent or ring shapes confuse it. Use DBSCAN, spectral clustering, or density-based methods.
High dimensions Distance metrics become less meaningful ("curse of dimensionality"). All points equidistant. Reduce dimensionality first (PCA, UMAP). Use cosine distance instead of Euclidean.
Outliers KMeans (Euclidean distance) sensitive to outliers. One outlier can drag a center away. Preprocess: remove/cap outliers. Or use robust distance (Mahalanobis, Manhattan).

Python Implementation from Scratch

Here's a complete, class-based KMeans implementation using only numpy and pandas. Includes KMeans++, multiple distance metrics, and elbow method.

import numpy as np
import pandas as pd
from typing import Tuple, Optional

class KMeansClustering:
    """
    KMeans clustering algorithm with KMeans++ initialization,
    multiple distance metrics, and helper methods.
    """

    def __init__(
        self,
        n_clusters: int = 3,
        max_iter: int = 300,
        tol: float = 1e-4,
        random_state: int = 42,
        metric: str = 'euclidean',
        init: str = 'kmeans++'
    ):
        """
        Parameters:
        -----------
        n_clusters : int
            Number of clusters (K)
        max_iter : int
            Maximum iterations before stopping
        tol : float
            Tolerance for convergence (change in center positions)
        random_state : int
            Random seed for reproducibility
        metric : str
            Distance metric: 'euclidean', 'manhattan', 'cosine'
        init : str
            Initialization method: 'random', 'kmeans++'
        """
        self.n_clusters = n_clusters
        self.max_iter = max_iter
        self.tol = tol
        self.random_state = random_state
        self.metric = metric
        self.init = init

        self.centers = None
        self.labels = None
        self.wcss = None
        self.wcss_history = []

        np.random.seed(random_state)

    def _euclidean_distance(self, X: np.ndarray, center: np.ndarray) -> np.ndarray:
        """Compute Euclidean distance from points to a center."""
        return np.sqrt(np.sum((X - center) ** 2, axis=1))

    def _manhattan_distance(self, X: np.ndarray, center: np.ndarray) -> np.ndarray:
        """Compute Manhattan distance from points to a center."""
        return np.sum(np.abs(X - center), axis=1)

    def _cosine_distance(self, X: np.ndarray, center: np.ndarray) -> np.ndarray:
        """Compute cosine distance from points to a center."""
        # Cosine similarity: (A·B) / (||A|| ||B||)
        dot_product = np.dot(X, center)
        norm_X = np.linalg.norm(X, axis=1)
        norm_center = np.linalg.norm(center)

        # Avoid division by zero
        with np.errstate(divide='ignore', invalid='ignore'):
            similarity = dot_product / (norm_X * norm_center + 1e-8)
            similarity = np.nan_to_num(similarity)

        return 1 - similarity  # Convert similarity to distance

    def _compute_distance_matrix(self, X: np.ndarray) -> np.ndarray:
        """
        Compute distance from each point to each center.
        Returns: (n_samples, n_clusters) array
        """
        distances = np.zeros((X.shape[0], self.n_clusters))

        for k in range(self.n_clusters):
            if self.metric == 'euclidean':
                distances[:, k] = self._euclidean_distance(X, self.centers[k])
            elif self.metric == 'manhattan':
                distances[:, k] = self._manhattan_distance(X, self.centers[k])
            elif self.metric == 'cosine':
                distances[:, k] = self._cosine_distance(X, self.centers[k])
            else:
                raise ValueError(f"Unknown metric: {self.metric}")

        return distances

    def _initialize_centers(self, X: np.ndarray):
        """Initialize cluster centers using specified method."""
        if self.init == 'random':
            indices = np.random.choice(X.shape[0], self.n_clusters, replace=False)
            self.centers = X[indices].copy()

        elif self.init == 'kmeans++':
            # KMeans++ initialization
            self.centers = np.zeros((self.n_clusters, X.shape[1]))

            # Pick first center randomly
            first_idx = np.random.randint(X.shape[0])
            self.centers[0] = X[first_idx]

            # Pick remaining centers
            for k in range(1, self.n_clusters):
                # Compute distances to nearest center
                distances = self._compute_distance_matrix(X)
                min_distances = np.min(distances, axis=1)

                # Pick next center with probability ∝ distance²
                probabilities = min_distances ** 2
                probabilities /= np.sum(probabilities)

                cumulative = np.cumsum(probabilities)
                r = np.random.rand()
                next_idx = np.searchsorted(cumulative, r)
                self.centers[k] = X[next_idx]
        else:
            raise ValueError(f"Unknown init method: {self.init}")

    def fit(self, X: np.ndarray) -> 'KMeansClustering':
        """
        Fit the KMeans model to data.

        Parameters:
        -----------
        X : np.ndarray
            Data array of shape (n_samples, n_features)

        Returns:
        --------
        self : KMeansClustering
        """
        X = np.asarray(X, dtype=np.float64)

        # Initialize centers
        self._initialize_centers(X)

        self.wcss_history = []

        # Main EM loop
        for iteration in range(self.max_iter):
            # E-step: Assign points to nearest center
            distances = self._compute_distance_matrix(X)
            self.labels = np.argmin(distances, axis=1)

            # Compute WCSS for this iteration
            wcss = np.sum(np.min(distances, axis=1) ** 2)
            self.wcss_history.append(wcss)

            # M-step: Recompute centers
            old_centers = self.centers.copy()

            for k in range(self.n_clusters):
                mask = self.labels == k
                if np.sum(mask) > 0:
                    self.centers[k] = np.mean(X[mask], axis=0)
                # If cluster is empty, keep old center

            # Check convergence
            center_shift = np.max(np.linalg.norm(self.centers - old_centers, axis=1))

            if center_shift < self.tol:
                print(f"Converged at iteration {iteration}")
                break

        # Final assignment
        distances = self._compute_distance_matrix(X)
        self.labels = np.argmin(distances, axis=1)
        self.wcss = np.sum(np.min(distances, axis=1) ** 2)

        return self

    def predict(self, X: np.ndarray) -> np.ndarray:
        """
        Predict cluster labels for new data.

        Parameters:
        -----------
        X : np.ndarray
            Data array of shape (n_samples, n_features)

        Returns:
        --------
        labels : np.ndarray
            Cluster labels for each point
        """
        if self.centers is None:
            raise ValueError("Model not fitted yet. Call fit() first.")

        X = np.asarray(X, dtype=np.float64)
        distances = self._compute_distance_matrix(X)
        return np.argmin(distances, axis=1)

    def fit_predict(self, X: np.ndarray) -> np.ndarray:
        """Fit model and return cluster labels."""
        self.fit(X)
        return self.labels


def elbow_method(X: np.ndarray, K_range: range = range(1, 11)) -> Tuple[list, list]:
    """
    Compute WCSS for different K values (for elbow plot).

    Parameters:
    -----------
    X : np.ndarray
        Data array
    K_range : range
        Range of K values to test

    Returns:
    --------
    K_values : list
        K values tested
    wcss_values : list
        WCSS for each K
    """
    wcss_values = []
    K_values = list(K_range)

    for k in K_values:
        kmeans = KMeansClustering(n_clusters=k, init='kmeans++')
        kmeans.fit(X)
        wcss_values.append(kmeans.wcss)

    return K_values, wcss_values


def silhouette_coefficient(X: np.ndarray, labels: np.ndarray) -> float:
    """
    Compute average silhouette coefficient.

    Parameters:
    -----------
    X : np.ndarray
        Data array
    labels : np.ndarray
        Cluster labels

    Returns:
    --------
    avg_silhouette : float
        Average silhouette score in [-1, 1]
    """
    n = X.shape[0]
    silhouettes = []

    for i in range(n):
        # a(i): mean distance to points in same cluster
        same_cluster = X[labels == labels[i]]
        if len(same_cluster) > 1:
            a_i = np.mean(np.linalg.norm(X[i] - same_cluster, axis=1))
        else:
            a_i = 0

        # b(i): min mean distance to other clusters
        b_i = np.inf
        for c in np.unique(labels):
            if c != labels[i]:
                other_cluster = X[labels == c]
                b_c = np.mean(np.linalg.norm(X[i] - other_cluster, axis=1))
                b_i = min(b_i, b_c)

        # Silhouette coefficient
        if max(a_i, b_i) == 0:
            s_i = 0
        else:
            s_i = (b_i - a_i) / max(a_i, b_i)

        silhouettes.append(s_i)

    return np.mean(silhouettes)


# ==============================================================================
# EXAMPLE: AdTech User Segmentation
# ==============================================================================

if __name__ == "__main__":
    # Create synthetic user data
    np.random.seed(42)

    # Cluster 1: Power users (high spend, high CTR)
    power_users = np.random.normal(loc=[50, 0.08, 1000], scale=[5, 0.01, 200], size=(100, 3))

    # Cluster 2: Casual users (medium spend, low CTR)
    casual_users = np.random.normal(loc=[25, 0.02, 200], scale=[5, 0.005, 50], size=(100, 3))

    # Cluster 3: New users (low spend, medium CTR)
    new_users = np.random.normal(loc=[10, 0.03, 50], scale=[3, 0.01, 30], size=(100, 3))

    X = np.vstack([power_users, casual_users, new_users])

    print("=" * 70)
    print("ADTECH USER SEGMENTATION WITH KMEANS")
    print("=" * 70)
    print(f"Data shape: {X.shape}")
    print(f"Features: [session_time (min), ctr, spend ($)]")

    # 1. ELBOW METHOD
    print("\n1. ELBOW METHOD")
    print("-" * 70)
    K_values, wcss_values = elbow_method(X, range(1, 8))
    for k, wcss in zip(K_values, wcss_values):
        print(f"  K={k}: WCSS={wcss:.2f}")

    # 2. FIT WITH K=3 (TRUE K)
    print("\n2. FITTING WITH K=3")
    print("-" * 70)
    kmeans = KMeansClustering(n_clusters=3, init='kmeans++', metric='euclidean')
    kmeans.fit(X)

    print(f"Final WCSS: {kmeans.wcss:.2f}")
    print(f"Convergence: {len(kmeans.wcss_history)} iterations")

    print("\nCluster Centers:")
    for k, center in enumerate(kmeans.centers):
        print(f"  Cluster {k}: session_time={center[0]:.1f}min, ctr={center[1]:.3f}, spend=${center[2]:.0f}")

    print("\nCluster Assignments:")
    for k in range(3):
        count = np.sum(kmeans.labels == k)
        print(f"  Cluster {k}: {count} users")

    # 3. SILHOUETTE SCORE
    print("\n3. SILHOUETTE SCORE")
    print("-" * 70)
    sil_score = silhouette_coefficient(X, kmeans.labels)
    print(f"Average Silhouette: {sil_score:.3f}")

    # Compare with K=2 and K=4
    for k in [2, 4]:
        kmeans_k = KMeansClustering(n_clusters=k, init='kmeans++')
        kmeans_k.fit(X)
        sil = silhouette_coefficient(X, kmeans_k.labels)
        print(f"K={k}: Silhouette={sil:.3f}")

    print("\n" + "=" * 70)
    print("✓ Clustering complete! Use Cluster IDs for audience segmentation.")
    print("=" * 70)

Key Features

Class Structure

KMeansClustering: Full implementation
• fit(X): Train the model
• predict(X): Label new data
• fit_predict(X): Both
• Support for Euclidean, Manhattan, Cosine distances

Helper Functions

elbow_method(): Compute WCSS for K range
silhouette_coefficient(): Quantify cluster quality
Both use only numpy, no external deps

Usage Example

# Basic usage
kmeans = KMeansClustering(n_clusters=5, init='kmeans++')
kmeans.fit(user_data)
cluster_ids = kmeans.labels

# Find optimal K
K_values, wcss = elbow_method(user_data, range(1, 11))
# Plot K_values vs wcss, find elbow

# Evaluate cluster quality
silhouette = silhouette_coefficient(user_data, cluster_ids)
if silhouette > 0.5:
    print("Good clustering!")

# Predict on new data
new_users = get_new_user_data()
new_cluster_ids = kmeans.predict(new_users)
Production notes: This implementation is educational and works well for datasets up to ~1M rows. For larger datasets, use scikit-learn's KMeans (C-optimized) or MiniBatchKMeans. For streaming/online scenarios, use incremental variants. Always normalize features first if they have different scales!
End of K-Means Clustering