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

从零构建K-means聚类引擎：NumPy实战与算法深度解析

在机器学习领域，聚类算法如同一位无声的组织者，能够从看似无序的数据中发现隐藏的结构。当我们使用sklearn的KMeans时，一行代码就能完成复杂的分群工作，但这种便利性也让我们错过了理解算法本质的机会。本文将带您穿越API的抽象层，用NumPy亲手搭建K-means引擎，体验从数学原理到代码实现的完整创造过程。

1. K-means算法核心原理拆解

K-means的本质是通过迭代优化来最小化簇内平方和（WCSS）。想象一位城市规划师，他需要合理设置k个消防站的位置，使得城市中任意一点到最近消防站的距离之和最小。这个类比完美诠释了K-means的核心任务。

算法流程的数学表达：

1. 初始化：随机选择k个点作为初始质心 $μ_1, μ_2, ..., μ_k$
2. 分配阶段：对每个样本 $x_i$，计算其到各质心的距离，分配到最近的簇 $$S_t = {x_p : |x_p - μ_t|^2 \leq |x_p - μ_j|^2 \ ∀j, 1≤j≤k}$$
3. 更新阶段：重新计算每个簇的质心 $$μ_t = \frac{1}{|S_t|} \sum_{x_i \in S_t} x_i$$
4. 迭代：重复2-3步直到质心变化小于阈值或达到最大迭代次数

def initialize_centroids(X, k): """随机初始化质心""" indices = np.random.choice(len(X), k, replace=False) return X[indices]

注意：质心初始化对结果有重大影响，好的初始值能减少迭代次数并避免局部最优

2. 距离计算的工程实现艺术

欧氏距离虽然是K-means的默认选择，但在实际实现时需要兼顾精度和效率。我们对比几种常见实现方式：

优化后的距离矩阵计算可同时处理多个样本：

def batch_distance(X, centers): """向量化计算所有样本到各中心的距离""" # X形状(n_samples, n_features) # centers形状(n_clusters, n_features) distances = np.sqrt(((X[:, np.newaxis] - centers)**2).sum(axis=2)) return distances # 形状(n_samples, n_clusters)

空簇处理技巧：当某个簇失去所有样本时，常见的应对策略包括：

• 重新初始化该质心
• 将离当前质心最远的样本设为新质心
• 直接减少簇数量

3. 完整K-means引擎的实现

我们将算法分解为多个可测试的模块，构建一个工业级实现：

class KMeansFromScratch: def __init__(self, n_clusters=3, max_iter=300, tol=1e-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, atol=self.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), p=prob)) return X[indices] def _assign_clusters(self, X): distances = self._compute_distances(X, self.centroids) return np.argmin(distances, axis=1) 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(axis=0) return centroids def _compute_distances(self, X, centers): return np.sqrt(((X[:, np.newaxis] - centers)**2).sum(axis=2)) 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, axis=1)

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_clusters=3) our_labels = kmeans_ours.fit(X_scaled).predict(X_scaled) # 与sklearn对比 from sklearn.cluster import KMeans kmeans_sk = KMeans(n_clusters=3) sk_labels = kmeans_sk.fit_predict(X_scaled) # 评估指标 from sklearn.metrics import adjusted_rand_score print(f"Our ARI: {adjusted_rand_score(y, our_labels):.3f}") print(f"Sklearn ARI: {adjusted_rand_score(y, sk_labels):.3f}")

性能优化技巧：

• 对大数据集使用Mini-batch K-means
• 采用更高效的距离度量如余弦相似度
• 实现Elkan算法利用三角不等式减少距离计算
• 使用Cython或Numba加速关键循环

5. 高级话题与工程实践

特征缩放的影响： 不同特征的量纲会显著影响K-means结果。比较鸢尾花数据在原始空间和标准化空间的聚类效果：

常见陷阱与解决方案：

1. 局部最优：通过多次随机初始化选择最佳结果

   def multi_init_kmeans(X, n_clusters, n_init=10): best_score = -np.inf best_model = None for _ in range(n_init): model = KMeansFromScratch(n_clusters=n_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

2. 确定最佳K值：肘部法则与轮廓系数结合

   from sklearn.metrics import silhouette_score k_values = range(2, 8) silhouette_scores = [] for k in k_values: labels = KMeansFromScratch(n_clusters=k).fit(X).predict(X) silhouette_scores.append(silhouette_score(X, labels))

3. 分类变量处理：对于混合型数据，可采用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), p=prob)]) return np.array(centers)

Mini-batch K-means：适合大规模数据集

def mini_batch_kmeans(X, k, batch_size=100, max_iter=100): centers = initialize_centroids(X, k) for _ in range(max_iter): batch_indices = np.random.choice(len(X), batch_size, replace=False) batch = X[batch_indices] # 分配步骤 distances = np.sqrt(((batch[:, np.newaxis] - centers)**2).sum(axis=2)) labels = np.argmin(distances, axis=1) # 更新步骤 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_iter=100): # 归一化输入向量 X_norm = X / np.linalg.norm(X, axis=1, keepdims=True) centers = initialize_centroids(X_norm, k) for _ in range(max_iter): # 余弦相似度分配 similarities = X_norm @ centers.T # 矩阵乘法代替距离计算 labels = np.argmax(similarities, axis=1) # 更新质心（均值+归一化） new_centers = np.array([X_norm[labels == i].mean(axis=0) for i in range(k)]) new_centers /= np.linalg.norm(new_centers, axis=1, keepdims=True) if np.allclose(centers, new_centers): break centers = new_centers return centers