k-Means Clustering

k-means is an unsupervised partitioning method. Given $n$ observations and a prespecified number of clusters $K$, it seeks a partition represented by $K$ centroids that minimizes a within-cluster squared-error objective under the standard Euclidean formulation taught here.

MacQueen introduced the name k-means for a process that partitions an $N$-dimensional population into $k$ sets on the basis of a sample, motivated by within-class variance about conditional means. That 1967 procedure is sequential: each new point is added to the current nearest mean, and that mean is then updated. This page does not teach MacQueen's sequential process as the computational algorithm. The iterative method taught below is the batch assignment and centroid-update cycle associated with Lloyd's least-squares quantization and with the modern multivariate statement of Lloyd's algorithm.

k-means does not use class labels when constructing the partition. It is not a classifier. Clusters are mathematical groups under the selected objective. They are not automatically known chemical classes.

Let the data matrix be $X\in {\mathbb{R}}^{n\times p}$, with rows as samples and columns as variables. Write $C_1,\dots ,C_K$ for the clusters and ${\mu }_1,\dots ,{\mu }_K$ for the corresponding centroids. The within-cluster sum-of-squares objective is

Arthur and Vassilvitskii write the same quantity as a potential $\varphi$ equal to the sum, over observations, of squared distance to the closest center. Current scikit-learn 1.6 documentation defines inertia_ as the sum of squared distances of samples to their closest cluster center, weighted by sample weights if provided. This page does not equate that quantity with SciPy cluster.vq distortion, which uses a different convention.

There is no one-step closed form that yields a globally optimal k-means partition in general. MacQueen noted that there is no feasible general method that always yields an optimal partition, and that the k-means procedure need not converge to one. The practical algorithm alternates assignment and centroid update.

With centroids held fixed, each observation is assigned to a nearest centroid under squared Euclidean distance:

Lloyd derived the corresponding rule for least-squares quantization: mass at a point is assigned to the quantum that minimizes squared error. Because squaring is strictly increasing on nonnegative distances, minimizing $\Vert x_i-{\mu }_j{\Vert }_2$ and minimizing $\Vert x_i-{\mu }_j{\Vert }_2^2$ select the same nearest centroid. The objective $J$ remains a sum of squared distances.

Exact ties are mathematically non-unique. Arthur and Vassilvitskii state that ties may be broken arbitrarily if the method is consistent. MacQueen assigned tied points to the set of lower index. The educational code below uses the lowest cluster index. That is an implementation convention, not a unique scientific rule.

With assignments held fixed, each non-empty cluster centroid is replaced by the arithmetic mean of its members:

Lloyd showed that, for a given partition, the least-squares representative of each set is its center of mass. MacQueen noted that the mean minimizes squared error. Arthur and Vassilvitskii record the identity that, for a set $S$ with center of mass $c(S)$ and any point $z$,

The right-hand side is nonnegative, so the arithmetic mean uniquely minimizes the within-cluster sum of squares for a fixed assignment. Empty clusters are omitted from this conceptual update. Software may apply a separate replacement rule; this page does not define one universal empty-cluster action.

Lloyd-style k-means alternates nearest-centroid assignment and mean updates. Lloyd's Method I imposed the nearest-set rule and the center-of-mass rule in turn and showed that the quantization noise does not increase. Arthur and Vassilvitskii state the multivariate cycle: assign each point to the nearest center, replace each center by the mean of its assigned points, and repeat until the centers no longer change. The potential $\varphi$ is monotonically decreasing, which prevents cycling through distinct clusterings.

That decrease is not a claim of global optimality. Lloyd showed that the first-order stationarity conditions are necessary at a minimum but not sufficient for a global minimum, and that several local minima can exist. sklearn 1.6 states that the algorithm converges, given enough time, but possibly to a local minimum that depends on initialization. MATLAB documents that kmeans can reach a local minimum that depends on the starting points.

Stopping rules are implementation-specific. Arthur's statement stops when centers no longer change. sklearn 1.6 tol is a relative tolerance on the Frobenius norm of the difference between successive cluster-center matrices. MATLAB repeats until assignments do not change or MaxIter is reached. This page does not present one software tolerance as the definition of k-means.

The local solution reached depends on the initial centroids. Multiple initializations can produce different values of $J$. Random initialization is not universally "bad," and libraries do not implement random identically. sklearn 1.6 init="random" chooses $K$ observations from the rows of $X$. MATLAB 'Start','sample' likewise selects $K$ observations at random.

k-means++ is an initialization method for the same k-means objective, not a different clustering criterion. Arthur and Vassilvitskii select the first center uniformly at random from the data. Let $D(x)$ be the Euclidean distance from observation $x$ to the nearest already chosen center. Each subsequent center $x^{\prime }$ is chosen with probability

That is a seeding probability, not a cluster-membership probability. After the $K$ centers are chosen, ordinary Lloyd iterations follow. The expected potential after k-means++ seeding satisfies $\mathbb{E}[\varphi ]\le 8(\ln K+2){\varphi }_{\mathrm{OPT}}$, which is an $O(\log K)$ expected approximation relative to the optimal k-means objective. It is not a guarantee of the globally optimal clustering.

sklearn 1.6 init="k-means++" implements greedy k-means++, which differs from vanilla k-means++ by making several trials at each sampling step and keeping the best candidate. MATLAB 'Start','plus' is documented as k-means++ seeding citing Arthur and Vassilvitskii. Those implementations are not claimed to be bit-identical. sklearn n_init and MATLAB Replicates rerun clustering from different seeds and retain a lowest-objective run. Current sklearn 1.6 n_init="auto" uses 10 runs for init="random" or a callable, and 1 run for init="k-means++" or an explicit center array. MATLAB Replicates defaults to 1. No fixed restart count is universally sufficient, and extra restarts do not guarantee the global optimum.

Standard squared-Euclidean k-means is a centroid definition, an SSE objective, a nearest-centroid assignment rule, and an alternating update. It is not a one-step formula for a globally optimal partition.

The procedure is Lloyd's batch k-means. Empty-cluster replacement is not inserted into the generic loop.

Ordinary k-means requires $K$ before fitting. The algorithm does not discover the number of clusters. Crase, Hall, and Thennadil note that clustering algorithms will divide data into any requested number of clusters, whether or not those groups exist naturally in the data.

Choosing $K$ is a separate exploratory or model-selection question. A plot of $J$ versus $K$ is sometimes examined as an elbow heuristic: additional clusters reduce within-cluster dispersion, and a kink is taken as a candidate $K$. That kink is not mathematically guaranteed and is not always clearly defined. Silhouette summaries and related indices belong to a future cluster-validation resource. They are not derived here, and a silhouette maximum is not a proof of the true $K$.

Crase, Hall, and Thennadil surveyed cluster analysis of IR and NIR spectra. k-means appeared in that literature. Euclidean distance was the most common similarity measure, and squared Euclidean distance was used where algorithms minimize a sum of squared errors, including k-means. IR and NIR data are typically high-dimensional with relatively few samples. k-means does not automatically fail in that geometry, and it is not a canonical solution for every spectral clustering problem. Preprocessing, feature extraction, and the selected representation determine the distances that enter $J$.

Because the objective uses squared Euclidean distances, variables with larger numerical dispersion can dominate the geometry. Autoscaling is not always required. Spectroscopic choices such as SNV, MSC, derivatives, baseline correction, normalization, scaling, or PCA change what differences count as large. See Normalization & Scaling. PCA scores may be supplied to k-means as an optional representation. Crase et al. report PCA used both for dimension reduction and for visualization in surveyed IR/NIR clustering studies. PCA is not required, and retaining components to a fixed explained-variance threshold is not a universal clustering rule.

sklearn 1.6 notes that inertia assumes clusters that are convex and isotropic, and that the criterion can respond poorly to elongated clusters or irregular manifolds. That is a statement about the SSE centroid geometry, not a theorem that clusters must be spheres, and not a claim that k-means cannot handle unequal cluster sizes in every data set. Squared residuals give unusual observations a large contribution to $J$, so extreme points can move centroids. Cluster labels such as 0, 1, 2 are arbitrary identifiers. Agreement with known classes after the fact is not proof that k-means recovered true chemistry.

Relative to HCA, k-means seeks a flat $K$-cluster partition under a centroid SSE objective, whereas agglomerative HCA builds a hierarchy under a chosen linkage before any cut. Neither is universally superior. k-means is not k-NN. In k-means, $K$ is the number of unsupervised clusters. In k-NN, $k$ is the number of labeled neighbors in a supervised predictor. LDA and PLS-DA use class labels; ordinary exploratory k-means does not. sklearn predict assigns new rows to the nearest fitted centroid. That is an extra assignment rule, not supervised class prediction.

The first Python listing implements Lloyd iterations from supplied initial centers: squared-Euclidean assignment, mean update, objective evaluation, and stopping when assignments stabilize. It does not implement k-means++. Empty clusters raise an error rather than inventing a replacement rule. Ties use the lowest cluster index. Production examples use sklearn KMeans with algorithm="lloyd". No scaler is applied automatically. SciPy vq.kmeans is not used.

MATLAB examples use official kmeans with 'Distance','sqeuclidean'. That default makes each centroid the mean of its members. Cityblock, cosine, and correlation options exist in MATLAB and are not ordinary squared-Euclidean k-means. 'Start','plus' is MATLAB's k-means++-style seeding. Replicates in the listing is an explicit example count, not a universal default. MATLAB may enable an additional online-update phase; current documentation defaults OnlinePhase to off. sklearn and MATLAB stopping rules are not claimed to be identical. Compare partitions only after allowing for label permutation, and only under matched squared-Euclidean settings.

import numpy as np  def _validate_matrix(X):    X = np.asarray(X, dtype=float)    if X.ndim != 2:        raise ValueError("X must be a 2D array of samples by variables.")    if min(X.shape) == 0:        raise ValueError("X must have at least one row and one column.")    if not np.all(np.isfinite(X)):        raise ValueError("X must contain only finite values.")    return X  def cluster_centroid(X, members):    X = _validate_matrix(X)    idx = np.asarray(list(members), dtype=int)    if idx.size == 0:        raise ValueError("Cluster member sets must be non-empty.")    return X[idx].mean(axis=0)  def within_cluster_sse(X, members):    X = _validate_matrix(X)    idx = np.asarray(list(members), dtype=int)    centroid = cluster_centroid(X, idx)    return float(np.sum((X[idx] - centroid) ** 2))  def kmeans_objective(X, labels, centers):    X = _validate_matrix(X)    labels = np.asarray(labels, dtype=int)    centers = np.asarray(centers, dtype=float)    if labels.shape[0] != X.shape[0]:        raise ValueError("labels must have one entry per sample.")    if centers.ndim != 2 or centers.shape[1] != X.shape[1]:        raise ValueError("centers must have shape (K, p).")    J = 0.0    for j in range(centers.shape[0]):        members = np.flatnonzero(labels == j)        if members.size == 0:            continue        J += float(np.sum((X[members] - centers[j]) ** 2))    return J  def assign_nearest_centroid(X, centers):    X = _validate_matrix(X)    centers = np.asarray(centers, dtype=float)    if centers.ndim != 2 or centers.shape[1] != X.shape[1]:        raise ValueError("centers must have shape (K, p).")    if centers.shape[0] < 1:        raise ValueError("At least one centroid is required.")    n = X.shape[0]    labels = np.empty(n, dtype=int)    for i in range(n):        d2 = np.sum((centers - X[i]) ** 2, axis=1)        labels[i] = int(np.argmin(d2))    return labels  def update_centroids(X, labels, n_clusters):    X = _validate_matrix(X)    labels = np.asarray(labels, dtype=int)    n_clusters = int(n_clusters)    if n_clusters < 1:        raise ValueError("n_clusters must be at least 1.")    centers = np.zeros((n_clusters, X.shape[1]), dtype=float)    for j in range(n_clusters):        members = np.flatnonzero(labels == j)        if members.size == 0:            raise ValueError(                "Empty cluster encountered. This educational implementation "                "does not define a replacement rule."            )        centers[j] = X[members].mean(axis=0)    return centers  def lloyd_kmeans(X, init_centers, max_iter=100):    X = _validate_matrix(X)    centers = np.asarray(init_centers, dtype=float).copy()    if centers.ndim != 2 or centers.shape[1] != X.shape[1]:        raise ValueError("init_centers must have shape (K, p).")    max_iter = int(max_iter)    if max_iter < 1:        raise ValueError("max_iter must be at least 1.")    n_clusters = centers.shape[0]    history = []    n_iter = 0    labels = assign_nearest_centroid(X, centers)    for _ in range(max_iter):        n_iter += 1        centers = update_centroids(X, labels, n_clusters)        history.append(kmeans_objective(X, labels, centers))        new_labels = assign_nearest_centroid(X, centers)        if np.array_equal(new_labels, labels):            break        labels = new_labels    return labels, centers, np.asarray(history, dtype=float), n_iter 

1. 1.

   MacQueen, J. (1967). Some Methods for Classification and Analysis of Multivariate Observations. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1, 281-297.

2. 2.

   Lloyd, S. P. (1982). Least Squares Quantization in PCM. IEEE Transactions on Information Theory, 28(2), 129-137.

   doi:10.1109/TIT.1982.1056489
3. 3.

   Arthur, D., & Vassilvitskii, S. (2007). k-means++: The Advantages of Careful Seeding. Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, 1027-1035.

4. 4.

   Crase, S., Hall, B., & Thennadil, S. N. (2021). Cluster Analysis for IR and NIR Spectroscopy: Current Practices to Future Perspectives. Computers, Materials & Continua, 69(2), 1945-1965.

   doi:10.32604/cmc.2021.018517
5. 5.

   scikit-learn Developers (n.d.). KMeans and the clustering user guide. scikit-learn 1.6 documentation.

6. 6.

   The MathWorks, Inc. (n.d.). kmeans. MATLAB documentation.