k-Means Clustering: Initialization, k, and Assumptions

k-means seeks \(k\) centers that make the total squared distance from each observation to its assigned center small. It can partition a uniform cloud just as it partitions separated groups, so obtaining labels is not evidence that distinct populations exist. Its usefulness depends on the meaning of the distance and what the partition will be used for. Degenerate data, such as too few distinct observations, can also yield fewer than \(k\) occupied clusters.

Lloyd’s algorithm

The objective is \(\sum_i\|x_i-\mu_{c(i)}\|^2\), where \(x_i\) is observation \(i\), \(c(i)\) its assigned cluster, and \(\mu_j\) center \(j\). This sum is called inertia. The general global optimization problem is NP-hard. Lloyd’s algorithm alternates nearest-center assignment with replacing each nonempty cluster’s center by its mean. Each step minimizes the objective with the other part held fixed, so the objective cannot increase; it can stay unchanged. For points 0 and 2, center 1 gives squared error \(1+1=2\), compared with \(0+4=4\) at center 0.

In the code, broadcasting creates squared distances with shape (1500, 4): one row per observation and one column per center. argmin(1) chooses the closest center in each row. The teaching implementation retains an old center if its cluster becomes empty. It stops when assignments no longer change, and raises an error if the iteration budget is exhausted. The first library comparison uses exactly the same initial centers.

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import adjusted_rand_score

X, _ = make_blobs(n_samples=1500, centers=4, cluster_std=1.0, random_state=0)
rng = np.random.default_rng(0)
initial = X[rng.choice(len(X), 4, replace=False)].copy()
centers = initial.copy()
previous_labels = None
for it in range(1000):
    labels = ((X[:, None, :] - centers[None, :, :]) ** 2).sum(-1).argmin(1)
    centers = np.array([X[labels == j].mean(0) if (labels == j).any()
                        else centers[j] for j in range(4)])
    if previous_labels is not None and np.array_equal(labels, previous_labels):
        break
    previous_labels = labels.copy()
else:
    raise RuntimeError("assignment convergence not reached")

km = KMeans(n_clusters=4, init=initial, n_init=1, tol=0,
            max_iter=1000, algorithm="lloyd").fit(X)
best = KMeans(n_clusters=4, n_init=10, tol=0,
              max_iter=1000, random_state=0).fit(X)
print(f"from scratch: {it + 1} iterations, inertia {((X - centers[labels]) ** 2).sum():.4f}")
print(f"sklearn, same start: inertia {km.inertia_:.4f}")
print(f"same partition (ARI): {adjusted_rand_score(labels, km.labels_):.4f}")
print(f"sklearn, 10 k-means++ starts: inertia {best.inertia_:.4f}")
# from scratch: 20 iterations, inertia 2515.3439
# sklearn, same start: inertia 2515.3439
# same partition (ARI): 1.0000
# sklearn, 10 k-means++ starts: inertia 2515.3325

With the same start, both implementations give inertia 2515.3439 and adjusted Rand index (ARI) 1.0. ARI compares which pairs of observations share a cluster, accounting for chance agreement; cluster ID numbers need not match. This verifies the partition in this example, not equivalence for every input or stopping rule. Ten k-means++ starts find inertia 2515.3325, lower by about 0.0114. A different starting procedure can give a slightly different partition, and default tolerance-based stopping can also leave differences. Neither run establishes the global optimum.

Exact Lloyd updates reach a stable assignment in the usual nondegenerate case with consistent tie handling. The objective is bounded below by zero, and only finitely many assignments are possible; strict improvement prevents revisiting an assignment. Ties and empty clusters require a specified policy, and floating-point implementations also use tolerances and iteration limits. The resulting fixed point need not be globally best. Initial centers, stopping rules, and numerical details all matter.

Initialization and repeated starts

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

for k, std in ((20, 0.4), (40, 0.3)):
    X, _ = make_blobs(n_samples=3000, centers=k, cluster_std=std,
                      center_box=(-30, 30), random_state=1)
    for init in ("random", "k-means++"):
        inertias = [KMeans(n_clusters=k, init=init, n_init=1,
                           random_state=s).fit(X).inertia_ for s in range(30)]
        print(f"k={k:3d}  {init:10s} best {min(inertias):9.1f}"
              f"  worst {max(inertias):9.1f}  mean {np.mean(inertias):9.1f}")
    print()
# k= 20  random     best    3110.2  worst   72845.4  mean   21009.2
# k= 20  k-means++  best     951.8  worst    1464.1  mean    1060.2
# 
# k= 40  random     best    2802.6  worst   28472.4  mean   10521.9
# k= 40  k-means++  best     521.4  worst     693.7  mean     558.1

For 20 generated groups, random starts have mean inertia 21,009 against 1,060 for k-means++, about 19.8 times as large. Comparing the worst runs gives 72,845 versus 1,464, about 49.8 times. Even the best random run, 3,110, is about 2.1 times the worst k-means++ run. These are comparisons over 30 seeds on this dataset; the table does not make a universal performance claim.

A random start can put several centers in one region while leaving another poorly represented. Subsequent local updates may split one group and merge others. Centers are not confined to regions containing observations: if a center receives points on both sides of a gap, their mean can lie in that gap. Poor results arise from the local assignments and updates settling into an unfavorable configuration, not from a rule forbidding movement across empty space.

In the original k-means++ scheme, the first center is sampled uniformly and later centers are sampled in proportion to squared distance from the nearest existing center. Distant, poorly represented regions then have a greater chance of receiving a center. Scikit-learn uses a greedy variant that tries several sampled candidates before choosing each new center. n_init specifies multiple complete runs for KMeans, keeping the lowest-inertia result. The benchmark sets it to 1 to isolate seed variability; the k-selection examples explicitly request 10 starts. Neither seeding nor restarts guarantee the best partition.

Choosing k

import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import silhouette_score

X, _ = make_blobs(n_samples=1500, centers=4, cluster_std=1.0, random_state=0)
print(f"{'k':>3} {'inertia':>12} {'silhouette':>12}")
for k in range(2, 10):
    m = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
    print(f"{k:3d} {m.inertia_:12.1f} {silhouette_score(X, m.labels_):12.4f}")
#   k      inertia   silhouette
#   2       7373.4       0.4619
#   3       4037.1       0.4796
#   4       2515.3       0.5086
#   5       2217.2       0.4366
#   6       1959.7       0.4258
#   7       1733.8       0.3947
#   8       1517.5       0.3412
#   9       1365.8       0.3437

The minimum achievable inertia cannot increase when another center is allowed: retaining the old centers remains an option. Independently fitted local solutions need not obey that ordering, although the printed values decrease here. An elbow is a change in the rate of improvement, not simply the largest drop. Here the reductions become much smaller after four centers. Choosing how sharp a bend counts as an elbow is still a judgment.

For a point, let \(a\) be its mean distance to other points in its assigned cluster and \(b\) the smallest mean distance to any other cluster. For \(\max(a,b)>0\), its silhouette is \((b-a)/\max(a,b)\). For example, \(a=1\) and \(b=3\) give \(2/3\). The average score lies between −1 and 1; larger values favor compact, separated groups under the chosen distance. Singleton clusters receive score zero in scikit-learn. The score peaks at four in this sweep, matching the number of generating groups, but that agreement is not guaranteed for other shapes, densities, or objectives. This calculation requires at least two occupied clusters and fewer clusters than observations.

Relate the partition to its intended use. Five customer segments and fifty impose different operational demands even if both have favorable internal scores. Resampling stability, comparison with a meaningful reference distribution, external information, or performance on a downstream task can provide different evidence. Choose the check that addresses the question rather than treating any one score as a test for the existence of clusters.

Which geometry the objective favors

Squared Euclidean distance favors compact groups around centers. The following are tendencies of the objective, not formal requirements that every cluster have the same shape, count, or variance.

• Elongation. Nearest-center regions are convex Voronoi cells, not necessarily spheres. A long group may be split if that reduces distortion enough.
• Unequal counts and spreads. A group with many observations or large spread can gain more from an extra center than a small tight group. K-means does not enforce equal cluster sizes, and sufficiently separated unequal groups can still be recovered.
• Concentric components. Two components with the same mean but different variances are distinguished by their density profiles. Nearest-centroid assignment alone does not model those profiles.

Euclidean distance gives equal coefficients to squared coordinate differences; it does not explicitly weight features by their variance. Multiplying one feature’s units by 1,000 multiplies its contribution by a million. Standardization can remove that unit effect, but equalizing scales is a modeling choice, not a mandatory step for data already expressed in meaningful units. Squared distance also gives outliers substantial influence. Fit any learned scaling or feature transformation on training data before a held-out assessment.

A full-covariance Gaussian mixture can represent elongated groups with different spreads and mixing proportions, although fitting remains sensitive to initialization. Other clustering methods use different geometry or density criteria and have their own limitations. Lloyd assignment costs \(O(nkd)\) per iteration for \(n\) points and \(d\) features. Full-covariance mixture calculations typically include \(O(nkd^2)\) work and covariance factorization; diagonal or spherical mixtures have different costs. Compare the actual covariance model rather than assigning one cost to every mixture.

Exercises

1. Partitioning a uniform cloud. Run k-means on a uniform square with no separated generating components. Inspect inertia and silhouette across candidate values of \(k\).

Can either curve favor a partition even though the generating distribution contains no distinct groups?

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

rng = np.random.default_rng(0)
X = rng.uniform(0, 1, (2000, 2))               # uniform noise, no clusters

print(f"{'k':>3} {'inertia':>10} {'drop':>8} {'silhouette':>12}")
prev = None
for k in range(2, 9):
    m = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
    drop = "" if prev is None else f"{prev - m.inertia_:8.1f}"
    print(f"{k:3d} {m.inertia_:10.2f} {drop:>8} {silhouette_score(X, m.labels_):12.4f}")
    prev = m.inertia_
#   k    inertia     drop   silhouette
#   2     206.53                0.3563
#   3     130.18     76.4       0.3808
#   4      82.89     47.3       0.4080
#   5      70.27     12.6       0.3866
#   6      59.20     11.1       0.3789
#   7      49.61      9.6       0.3769
#   8      41.94      7.7       0.3822

2. Different group sizes and spreads. Fit k-means and Gaussian mixtures to one large diffuse group and two small tight ones. Compare ARI, sizes, and cross-tabulations against the known generating labels.

A column with 240 points might contain both small groups. Use the cross-tabulation to tell it apart from a recovered component.

import numpy as np
from sklearn.cluster import KMeans
from sklearn.mixture import GaussianMixture
from sklearn.metrics import adjusted_rand_score
from sklearn.metrics.cluster import contingency_matrix

rng = np.random.default_rng(0)
big = rng.normal([0, 0], 2.0, size=(1500, 2))
small_a = rng.normal([9, 2], 0.35, size=(120, 2))
small_b = rng.normal([9, -2], 0.35, size=(120, 2))
X = np.vstack([big, small_a, small_b])
y = np.repeat([0, 1, 2], [1500, 120, 120])

km = KMeans(n_clusters=3, n_init=10, random_state=0).fit(X)
print(f"k-means ARI {adjusted_rand_score(y, km.labels_):.4f} sizes {np.bincount(km.labels_)}")
print("k-means: true rows vs fitted columns")
print(contingency_matrix(y, km.labels_))
true_centers = np.array([X[y == j].mean(0) for j in range(3)])
print(f"inertia: k-means {km.inertia_:.1f}, generating groups {((X - true_centers[y])**2).sum():.1f}")

for name, model in (
    ("default 1", GaussianMixture(3, n_init=1, random_state=0)),
    ("default 10", GaussianMixture(3, n_init=10, random_state=0)),
    ("kmeans++ 10", GaussianMixture(3, n_init=10, init_params="k-means++", random_state=0)),
    ("known start", GaussianMixture(3, n_init=1, random_state=0,
        means_init=np.array([[0., 0.], [9., 2.], [9., -2.]]),
        weights_init=np.array([1500, 120, 120]) / 1740,
        precisions_init=np.array([np.eye(2)/4, np.eye(2)/0.35**2, np.eye(2)/0.35**2])))):
    fitted = model.fit_predict(X)
    print(f"{name:12s} ARI {adjusted_rand_score(y, fitted):.4f}"
          f" mean loglik {model.score(X):.4f} sizes {np.bincount(fitted)}")
    if name in ("kmeans++ 10", "known start"):
        print(contingency_matrix(y, fitted))
# k-means ARI 0.3088 sizes [761 245 734]
# k-means: true rows vs fitted columns
# [[761   5 734]
#  [  0 120   0]
#  [  0 120   0]]
# inertia: k-means 8955.6, generating groups 11925.8
# default 1    ARI 0.3223 mean loglik -4.3806 sizes [644 856 240]
# default 10   ARI 0.3201 mean loglik -4.3806 sizes [842 240 658]
# kmeans++ 10  ARI 0.5668 mean loglik -4.3795 sizes [1272  240  228]
# [[1272    0  228]
#  [   0  120    0]
#  [   0  120    0]]
# known start  ARI 1.0000 mean loglik -4.2426 sizes [1500  120  120]
# [[1500    0    0]
#  [   0  120    0]
#  [   0    0  120]]

3. Full-data and incremental fitting. Compare one full k-means fit with one pass of partial_fit on two million generated rows. Evaluate both final models on exactly the same rows, regenerated in batches. Report objective, fitting time, and input batch size separately.

Distinguish the largest input batch from measured peak process memory. Which work is excluded from the fit timer?

import time, gc
import numpy as np
from sklearn.cluster import KMeans, MiniBatchKMeans
from threadpoolctl import threadpool_limits

n, d, k, chunk_size = 2_000_000, 20, 25, 50_000
true_centers = np.random.default_rng(0).uniform(-10, 10, (k, d))

def chunks():
    rng = np.random.default_rng(1)
    for start in range(0, n, chunk_size):
        size = min(chunk_size, n - start)
        labels = rng.integers(k, size=size)
        yield true_centers[labels] + rng.normal(0, 1.5, (size, d))

with threadpool_limits(limits=1):
    X = np.empty((n, d))
    for j, block in enumerate(chunks()):
        X[j * chunk_size:j * chunk_size + len(block)] = block
        del block
    t0 = time.perf_counter()
    full = KMeans(n_clusters=k, n_init=1, algorithm="lloyd", random_state=0).fit(X)
    t_full = time.perf_counter() - t0
    del X
    gc.collect()

    mb = MiniBatchKMeans(n_clusters=k, batch_size=chunk_size, n_init=1, random_state=0)
    t_stream = 0.0
    for block in chunks():
        t0 = time.perf_counter()
        mb.partial_fit(block)
        t_stream += time.perf_counter() - t0
        del block

    full_sse = stream_sse = 0.0
    for block in chunks():
        full_sse -= full.score(block)
        stream_sse -= mb.score(block)
        del block

print(f"full fit {t_full:.2f} s")
print(f"partial_fit calls {t_stream:.2f} s")
print(f"full inertia {full_sse:.1f}")
print(f"stream inertia {stream_sse:.1f}")
print(f"largest input batch: full {n:,} rows, stream {chunk_size:,} rows")
print(f"feature array bytes: full {n*d*8/2**20:.2f} MiB, chunk {chunk_size*d*8/2**20:.2f} MiB")
# full fit 5.65 s # varies by machine
# partial_fit calls 0.67 s # varies by machine
# full inertia 109942554.4
# stream inertia 109941009.9
# largest input batch: full 2,000,000 rows, stream 50,000 rows
# feature array bytes: full 305.18 MiB, chunk 7.63 MiB

References

• Arthur, D., & Vassilvitskii, S. (2007). k-means++: The Advantages of Careful Seeding. SODA.
• Sculley, D. (2010). Web-Scale K-Means Clustering. WWW.