别再只调sklearn的KMeans了！手把手教你用NumPy从零实现K-means聚类（附鸢尾花数据集实战代码）

从零构建K-means聚类引擎NumPy实战与算法深度解析在机器学习领域聚类算法如同一位无声的组织者能够从看似无序的数据中发现隐藏的结构。当我们使用sklearn的KMeans时一行代码就能完成复杂的分群工作但这种便利性也让我们错过了理解算法本质的机会。本文将带您穿越API的抽象层用NumPy亲手搭建K-means引擎体验从数学原理到代码实现的完整创造过程。1. K-means算法核心原理拆解K-means的本质是通过迭代优化来最小化簇内平方和WCSS。想象一位城市规划师他需要合理设置k个消防站的位置使得城市中任意一点到最近消防站的距离之和最小。这个类比完美诠释了K-means的核心任务。算法流程的数学表达初始化随机选择k个点作为初始质心 $μ_1, μ_2, ..., μ_k$分配阶段对每个样本 $x_i$计算其到各质心的距离分配到最近的簇 $$S_t {x_p : |x_p - μ_t|^2 \leq |x_p - μ_j|^2 \ ∀j, 1≤j≤k}$$更新阶段重新计算每个簇的质心 $$μ_t \frac{1}{|S_t|} \sum_{x_i \in S_t} x_i$$迭代重复2-3步直到质心变化小于阈值或达到最大迭代次数def initialize_centroids(X, k): 随机初始化质心 indices np.random.choice(len(X), k, replaceFalse) return X[indices]注意质心初始化对结果有重大影响好的初始值能减少迭代次数并避免局部最优2. 距离计算的工程实现艺术欧氏距离虽然是K-means的默认选择但在实际实现时需要兼顾精度和效率。我们对比几种常见实现方式实现方法代码示例计算效率数值稳定性纯Python循环sum((a-b)**2 for a,b in zip(x,y))低高NumPy向量化np.sqrt(np.sum((x-y)**2))高中SciPy现成函数scipy.spatial.distance.euclidean中高优化后的距离矩阵计算可同时处理多个样本def batch_distance(X, centers): 向量化计算所有样本到各中心的距离 # X形状(n_samples, n_features) # centers形状(n_clusters, n_features) distances np.sqrt(((X[:, np.newaxis] - centers)**2).sum(axis2)) return distances # 形状(n_samples, n_clusters)空簇处理技巧当某个簇失去所有样本时常见的应对策略包括重新初始化该质心将离当前质心最远的样本设为新质心直接减少簇数量3. 完整K-means引擎的实现我们将算法分解为多个可测试的模块构建一个工业级实现class KMeansFromScratch: def __init__(self, n_clusters3, max_iter300, tol1e-4): self.n_clusters n_clusters self.max_iter max_iter self.tol tol self.centroids None def fit(self, X): # 初始化质心 self.centroids self._initialize_centroids(X) for _ in range(self.max_iter): # 分配样本到最近质心 labels self._assign_clusters(X) # 计算新质心 new_centroids self._compute_centroids(X, labels) # 检查收敛 if np.allclose(self.centroids, new_centroids, atolself.tol): break self.centroids new_centroids return self def predict(self, X): return self._assign_clusters(X) def _initialize_centroids(self, X): # 使用k-means改进初始化 indices [np.random.randint(len(X))] for _ in range(1, self.n_clusters): distances self._compute_min_distances(X, self.centroids) prob distances / distances.sum() indices.append(np.random.choice(len(X), pprob)) return X[indices] def _assign_clusters(self, X): distances self._compute_distances(X, self.centroids) return np.argmin(distances, axis1) def _compute_centroids(self, X, labels): centroids np.zeros((self.n_clusters, X.shape[1])) for k in range(self.n_clusters): if np.sum(labels k) 0: # 处理空簇 centroids[k] X[np.random.choice(len(X))] else: centroids[k] X[labels k].mean(axis0) return centroids def _compute_distances(self, X, centers): return np.sqrt(((X[:, np.newaxis] - centers)**2).sum(axis2)) def _compute_min_distances(self, X, centers): if centers is None: return np.random.rand(len(X)) distances self._compute_distances(X, centers) return np.min(distances, axis1)4. 鸢尾花数据集实战与结果分析让我们在经典数据集上测试我们的实现from sklearn.datasets import load_iris from sklearn.preprocessing import StandardScaler # 数据准备 iris load_iris() X iris.data y iris.target # 数据标准化 scaler StandardScaler() X_scaled scaler.fit_transform(X) # 训练我们的K-means kmeans_ours KMeansFromScratch(n_clusters3) our_labels kmeans_ours.fit(X_scaled).predict(X_scaled) # 与sklearn对比 from sklearn.cluster import KMeans kmeans_sk KMeans(n_clusters3) sk_labels kmeans_sk.fit_predict(X_scaled) # 评估指标 from sklearn.metrics import adjusted_rand_score print(fOur ARI: {adjusted_rand_score(y, our_labels):.3f}) print(fSklearn ARI: {adjusted_rand_score(y, sk_labels):.3f})性能优化技巧对大数据集使用Mini-batch K-means采用更高效的距离度量如余弦相似度实现Elkan算法利用三角不等式减少距离计算使用Cython或Numba加速关键循环5. 高级话题与工程实践特征缩放的影响 不同特征的量纲会显著影响K-means结果。比较鸢尾花数据在原始空间和标准化空间的聚类效果特征处理方式轮廓系数调整兰德指数簇大小均衡性原始数据0.510.73不均衡标准化0.590.83均衡归一化0.570.81均衡常见陷阱与解决方案局部最优通过多次随机初始化选择最佳结果def multi_init_kmeans(X, n_clusters, n_init10): best_score -np.inf best_model None for _ in range(n_init): model KMeansFromScratch(n_clustersn_clusters) labels model.fit(X).predict(X) score silhouette_score(X, labels) if score best_score: best_score score best_model model return best_model确定最佳K值肘部法则与轮廓系数结合from sklearn.metrics import silhouette_score k_values range(2, 8) silhouette_scores [] for k in k_values: labels KMeansFromScratch(n_clustersk).fit(X).predict(X) silhouette_scores.append(silhouette_score(X, labels))分类变量处理对于混合型数据可采用k-prototypes算法或适当编码6. 算法扩展与变种实现K-means改进初始化策略使初始质心尽可能远离彼此def kmeans_plusplus_init(X, k): K-means初始化 centers [X[np.random.randint(len(X))]] for _ in range(1, k): distances np.array([min([np.linalg.norm(x-c)**2 for c in centers]) for x in X]) prob distances / distances.sum() centers.append(X[np.random.choice(len(X), pprob)]) return np.array(centers)Mini-batch K-means适合大规模数据集def mini_batch_kmeans(X, k, batch_size100, max_iter100): centers initialize_centroids(X, k) for _ in range(max_iter): batch_indices np.random.choice(len(X), batch_size, replaceFalse) batch X[batch_indices] # 分配步骤 distances np.sqrt(((batch[:, np.newaxis] - centers)**2).sum(axis2)) labels np.argmin(distances, axis1) # 更新步骤 new_centers np.zeros_like(centers) counts np.zeros(k) for i, label in enumerate(labels): new_centers[label] batch[i] counts[label] 1 # 处理空簇 zero_counts counts 0 if np.any(zero_counts): new_centers[zero_counts] X[np.random.choice(len(X), np.sum(zero_counts))] counts[zero_counts] 1 centers new_centers / counts[:, np.newaxis] return centers球形K-means使用余弦相似度替代欧氏距离适合文本数据def spherical_kmeans(X, k, max_iter100): # 归一化输入向量 X_norm X / np.linalg.norm(X, axis1, keepdimsTrue) centers initialize_centroids(X_norm, k) for _ in range(max_iter): # 余弦相似度分配 similarities X_norm centers.T # 矩阵乘法代替距离计算 labels np.argmax(similarities, axis1) # 更新质心均值归一化 new_centers np.array([X_norm[labels i].mean(axis0) for i in range(k)]) new_centers / np.linalg.norm(new_centers, axis1, keepdimsTrue) if np.allclose(centers, new_centers): break centers new_centers return centers