从零实现K-Means聚类MNIST：无监督学习与向量化编程实战

简介本资源是深圳大学计算机科学与软件工程专业《最优化方法》课程的实验实践材料面向机器学习初学者与高校相关课程学习者聚焦无监督学习中K-Means聚类算法在手写数字识别任务中的原理理解与Python工程实现。资源包共2个文件1个可运行Python脚本 1份详细实验说明文档总大小873KB脚本完整实现MNIST数据加载、归一化预处理、肘部法则确定K值、scikit-learn K-Means建模、聚类结果可视化及与真实标签的匹配评估文档则系统梳理算法原理、关键步骤代码注释、实验现象分析与常见问题提示。已有3175人学习下载内容紧扣课程实验要求兼顾理论推导与动手实践特别适合用于课程作业参考、算法复现训练及无监督分类入门拓展。1. 项目概述从理论到实践的聚类探索最近在复盘一些经典的机器学习项目正好翻到之前带学生做的一个课程实验——用最基础的K-Means算法去处理大名鼎鼎的MNIST手写数字数据集。这个实验听起来有点“复古”毕竟现在动辄就是深度神经网络、Transformer。但恰恰是这种“复古”项目最能考验你对算法本质的理解和工程实现的基本功。K-Means作为一个无监督聚类算法要让它去完成一个有明确标签0-9数字的分类任务这中间的Gap怎么弥合参数怎么调效果怎么评估这里面每一步都藏着不少门道。这个实验不仅适合计软专业的学生巩固《最优化方法》课程里关于迭代优化、损失函数的概念对于任何想从“调包”迈向“懂原理、能实现、会优化”的Python开发者来说都是一次绝佳的练手机会。今天我就把这个实验的完整实现思路、代码细节以及我踩过的那些坑系统地梳理一遍。2. 核心思路与方案设计当聚类遇上分类2.1 任务本质无监督算法处理有监督数据首先我们必须清醒地认识到一点K-Means是一个纯粹的聚类算法而MNIST是一个带有明确标签的分类数据集。我们的任务本质上是“借用”聚类的结果来近似实现分类的目标。这中间存在一个根本性的矛盾聚类依据的是样本间的“相似度”通常是欧氏距离而分类的黄金标准是样本的“真实标签”。因此这个实验的成功与否不在于达到多高的分类准确率那应该用监督学习算法而在于我们能否通过K-Means挖掘出数据内在的结构并将聚类结果与真实标签进行合理的映射和评估。所以我们的核心思路可以拆解为三步特征准备将28x28的灰度图像784维向量作为K-Means的输入特征。聚类过程运行K-Means算法将6万张训练图像划分到K个簇中。标签匹配与评估为每一个聚类簇“分配”一个数字标签然后计算聚类“分类”的准确率。2.2 为什么选择K-Means和MNIST选择K-Means的原因算法透明原理直观它的核心——迭代优化簇中心和样本归属——完美体现了最优化方法中的坐标下降思想。每一步的“求平均”和“找最近”都清晰可循便于理解损失函数簇内误差平方和的下降过程。实现简单扩展性强自己从零实现一个可用的K-Means不超过50行代码但其中涉及矩阵运算、循环优化、收敛判断等多个编程要点。理解它之后对其变种如K-Means优化初始化、Mini-Batch K-Means等也能快速上手。高维数据考验784维的特征空间是检验算法实现效率和数值稳定性的试金石。选择MNIST的原因干净、标准图像大小统一、背景干净、数字居中省去了大量数据清洗和预处理的麻烦让我们能聚焦于算法本身。维度适中784维不算太低能体现算法处理规模数据的能力也不算太高在个人电脑上运行毫无压力。结果可解释我们可以直观地查看每个簇的中心点也是一张28x28的图像看看算法“学”到了什么样的“平均数字”这非常有趣。2.3 工具选型与环境搭建这个实验对环境要求极低核心就是Python和几个科学计算库。我强烈建议新手从头配置理解每个包的作用。# 1. 创建并激活一个干净的虚拟环境避免包冲突 python -m venv venv_kmeans # Windows: venv_kmeans\Scripts\activate # Mac/Linux: source venv_kmeans/bin/activate # 2. 安装核心依赖 pip install numpy matplotlib scikit-learnNumPy (1.24)所有数据运算的基石。我们的图像数据、簇中心、距离计算都将用NumPy数组来完成。务必熟悉它的广播机制和向量化操作这是提升代码效率的关键。Matplotlib (3.7)用于可视化。我们将用它来绘制簇中心图像、损失函数下降曲线直观感受算法过程。Scikit-learn (1.3)这里我们主要“借用”它的两个功能一是datasets模块中的fetch_openml来下载MNIST数据集比直接下载文件更稳定二是用它的metrics模块中的函数来评估聚类效果如adjusted_rand_score。注意我们不会用它的KMeans类来实现聚类那是“作弊”。我们只用它来验证我们自己的实现结果。注意关于MNIST数据下载直接使用torchvision.datasets.MNIST或在某些网络环境下可能会遇到404错误。一个更可靠的方法是使用sklearn.datasets.fetch_openml它从OpenML平台获取数据通常更稳定。from sklearn.datasets import fetch_openml mnist fetch_openml(mnist_784, version1, cacheTrue, as_frameFalse) X, y mnist[data], mnist[target].astype(int)3. K-Means算法核心原理与手撕实现3.1 算法步骤拆解与数学表达K-Means的目标是最小化簇内误差平方和也叫惯性Inertia $$J \sum_{i1}^{k} \sum_{x \in C_i} ||x - \mu_i||^2$$ 其中$k$是簇数$C_i$是第$i$个簇的样本集合$\mu_i$是第$i$个簇的中心。算法通过两步交替迭代来逼近最优解分配步骤E-step对于每一个样本$x^{(j)}$计算它到所有簇中心$\mu_i$的距离将其分配给距离最近的簇。 $$c^{(j)} \arg\min_{i} ||x^{(j)} - \mu_i||^2$$ 其中$c^{(j)}$是样本$x^{(j)}$被分配到的簇索引。更新步骤M-step对于每一个簇$i$重新计算其中心为该簇内所有样本的均值。 $$\mu_i \frac{1}{|C_i|} \sum_{x^{(j)} \in C_i} x^{(j)}$$重复以上两步直到簇中心的变化小于某个阈值或达到最大迭代次数。3.2 从零开始的Python实现理解了数学原理我们就可以动手实现了。关键在于利用NumPy进行向量化计算避免低效的Python循环。import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import fetch_openml from sklearn.metrics import confusion_matrix, accuracy_score import seaborn as sns class MyKMeans: def __init__(self, n_clusters10, max_iter300, tol1e-4, random_state42): 初始化K-Means参数 :param n_clusters: 簇的数量对于MNIST就是10 :param max_iter: 最大迭代次数 :param tol: 收敛阈值当中心点移动距离小于此值时停止迭代 :param random_state: 随机种子保证结果可复现 self.n_clusters n_clusters self.max_iter max_iter self.tol tol self.random_state random_state self.centers None # 簇中心 self.labels None # 每个样本所属簇的标签 self.inertia_ None # 最终的簇内误差平方和 def _init_centers(self, X): 初始化簇中心。这里采用最简单的随机选择样本点法。 np.random.seed(self.random_state) idx np.random.choice(len(X), self.n_clusters, replaceFalse) return X[idx].copy() # 必须使用copy避免后续修改影响原数据 def fit(self, X): 训练拟合模型 :param X: 形状为 (n_samples, n_features) 的NumPy数组 n_samples, n_features X.shape # 1. 初始化中心 self.centers self._init_centers(X) for i in range(self.max_iter): # 2. E-step: 计算所有样本到所有中心的距离并分配标签 # 利用广播机制一次性计算所有距离矩阵 (n_samples, n_clusters) distances np.sqrt(((X[:, np.newaxis, :] - self.centers[np.newaxis, :, :]) ** 2).sum(axis2)) new_labels np.argmin(distances, axis1) # 3. M-step: 更新每个簇的中心 new_centers np.zeros_like(self.centers) for k in range(self.n_clusters): # 找到属于当前簇k的所有样本 mask (new_labels k) if mask.any(): # 防止空簇 new_centers[k] X[mask].mean(axis0) else: # 如果出现空簇重新随机初始化该中心 new_centers[k] X[np.random.randint(0, n_samples)] # 4. 检查收敛条件中心点的平均变化是否小于阈值 center_shift np.sqrt(((new_centers - self.centers) ** 2).sum(axis1)).mean() self.centers new_centers self.labels new_labels if i % 10 0: print(fIteration {i}, center shift: {center_shift:.6f}) if center_shift self.tol: print(fConverged at iteration {i} with shift {center_shift:.6f}) break # 计算最终的惯性 self.inertia_ 0 for k in range(self.n_clusters): mask (self.labels k) if mask.any(): self.inertia_ ((X[mask] - self.centers[k]) ** 2).sum() return self def predict(self, X): 预测新样本所属的簇 distances np.sqrt(((X[:, np.newaxis, :] - self.centers[np.newaxis, :, :]) ** 2).sum(axis2)) return np.argmin(distances, axis1)关键实现细节解析距离计算的向量化X[:, np.newaxis, :] - self.centers[np.newaxis, :, :]这行代码利用了NumPy的广播机制一次性计算了所有样本与所有中心点的差值形成了一个形状为(n_samples, n_clusters, n_features)的三维张量。然后平方、求和、开方得到距离矩阵。这比写双层循环快几个数量级。空簇处理在更新簇中心时有可能某个簇没有分配到任何样本空簇。我们的策略是随机选择一个数据点作为该簇的新中心。这是一种简单的处理方式更优的方案是采用K-Means的初始化策略来从根本上降低空簇概率。收敛判断我们计算所有簇中心新旧位置之间的欧氏距离的平均值如果小于阈值tol则认为算法已收敛。3.3 进阶优化K-Means初始化随机初始化可能导致算法收敛到局部最优解且速度较慢。K-Means通过一种概率化的方法选择初始中心点使得它们彼此尽可能远离从而显著提升聚类效果和速度。def _init_centers_plusplus(self, X): K-Means 初始化 np.random.seed(self.random_state) n_samples, n_features X.shape centers np.zeros((self.n_clusters, n_features)) # 1. 随机选择第一个中心 first_idx np.random.randint(n_samples) centers[0] X[first_idx] # 2. 选择后续中心 for i in range(1, self.n_clusters): # 计算每个样本到已有最近中心的距离的平方 distances np.array([min([np.linalg.norm(x - c) ** 2 for c in centers[:i]]) for x in X]) # 将距离平方转换为概率分布 probs distances / distances.sum() # 根据概率分布随机选择下一个中心 cumulative_probs probs.cumsum() r np.random.rand() next_idx np.searchsorted(cumulative_probs, r) centers[i] X[next_idx] return centers将fit方法中的初始化改为self.centers self._init_centers_plusplus(X)你会发现算法收敛所需的迭代次数通常会减少且最终得到的惯性Inertia更小即聚类效果更好。4. MNIST数据加载、预处理与实验流程4.1 稳定获取MNIST数据如前所述我们使用sklearn的接口这是目前最稳定可靠的方式。def load_mnist_data(): 加载MNIST数据集 print(Loading MNIST data from OpenML...) # 使用fetch_openml指定数据名和版本并关闭DataFrame格式以直接获取数组 mnist fetch_openml(mnist_784, version1, cacheTrue, as_frameFalse, parserauto) X, y mnist[data], mnist[target] # 确保标签是整数类型 y y.astype(np.int32) # 数据已经是0-255的uint8我们将其转换为0-1之间的浮点数有利于距离计算和模型稳定性 X X.astype(np.float32) / 255.0 print(fData shape: {X.shape}, Labels shape: {y.shape}) print(fLabel range: {np.unique(y)}) return X, y X, y load_mnist_data()4.2 数据划分与实验设置通常我们会用全部6万张图像进行聚类以观察算法在完整数据集上的表现。但为了后续评估“分类”准确率我们需要一个映射关系。一种常见的做法是用全部数据X进行聚类得到聚类标签cluster_labels。为每个聚类簇分配一个数字标签。分配规则是统计一个簇内所有样本的真实标签将出现次数最多的真实标签作为该簇的“预测标签”。用这个映射关系将聚类标签转换为预测的数字标签再与真实标签y比较计算准确率。def assign_cluster_labels(cluster_labels, true_labels, n_clusters10): 为每个聚类簇分配一个数字标签。 策略对于每个簇找出其内部样本最多的真实标签作为该簇的标签。 assigned_labels np.zeros(n_clusters, dtypenp.int32) for k in range(n_clusters): # 找到属于簇k的样本索引 idx (cluster_labels k) if idx.any(): # 获取这些样本的真实标签 true_labels_in_cluster true_labels[idx] # 找出出现次数最多的标签 counts np.bincount(true_labels_in_cluster, minlength10) assigned_labels[k] np.argmax(counts) else: # 空簇分配一个默认标签如-1后续处理 assigned_labels[k] -1 return assigned_labels def evaluate_clustering(X, y, n_clusters10, use_plusplusTrue): 完整的训练与评估流程 # 1. 实例化并训练模型 kmeans MyKMeans(n_clustersn_clusters, max_iter300, tol1e-4, random_state42) if use_plusplus: kmeans._init_centers kmeans._init_centers_plusplus.__get__(kmeans, MyKMeans) print(Fitting K-Means...) kmeans.fit(X) print(fFinal inertia: {kmeans.inertia_:.2f}) # 2. 获取聚类标签并分配数字标签 cluster_labels kmeans.labels_ assigned_labels assign_cluster_labels(cluster_labels, y, n_clusters) # 3. 将聚类标签映射为预测的数字标签 # 注意assigned_labels可能为-1空簇需要处理 pred_labels np.array([assigned_labels[l] if assigned_labels[l] ! -1 else 0 for l in cluster_labels]) # 4. 计算准确率注意这只是一种近似评估因为聚类本身不是分类 accuracy accuracy_score(y, pred_labels) print(fClustering-to-Classification Accuracy: {accuracy:.4f}) return kmeans, assigned_labels, pred_labels, accuracy4.3 运行实验与初步结果现在让我们在MNIST训练集前60000张图这是约定俗成的划分虽然我们做无监督聚类但标签y用于评估上运行我们的模型。# 通常我们使用前60000张作为“训练集”进行聚类 X_train, y_train X[:60000], y[:60000] kmeans_model, assigned_labs, pred_labs, acc evaluate_clustering(X_train, y_train, n_clusters10, use_plusplusTrue)在我的实验中使用K-Means初始化经过大约70-100次迭代收敛最终得到的“分类”准确率大约在**53%-58%**之间。这个数字远低于监督学习模型如CNN可以达到99%以上但作为纯无监督聚类的结果它已经揭示出数据内在的一些结构算法成功地将大部分相同数字的图片聚集到了一起。5. 结果可视化与深度分析5.1 可视化簇中心——算法的“认知”查看K-Means学到的簇中心是理解其工作方式最直观的方法。每个中心是一个784维的向量我们可以将其重塑为28x28的图像。def plot_cluster_centers(centers, assigned_labels): 绘制10个簇的中心图像及其被分配的数字标签 fig, axes plt.subplots(2, 5, figsize(12, 5)) axes axes.ravel() for i in range(10): ax axes[i] center_img centers[i].reshape(28, 28) ax.imshow(center_img, cmapgray) ax.set_title(fCluster {i}\nAssigned: {assigned_labels[i]}) ax.axis(off) plt.tight_layout() plt.show() plot_cluster_centers(kmeans_model.centers, assigned_labs)你会看到10张模糊的“平均数字”图像。通常每个簇中心会对应一个主要的数字如一个清晰的“0”或“1”但有时也会出现“四不像”这可能是因为该簇混淆了多个书写风格相似的数字例如“4”和“9”“3”和“8”。5.2 混淆矩阵分析——错误聚在了哪里准确率只是一个数字混淆矩阵能告诉我们更详细的信息哪些数字容易被混淆。def plot_confusion_matrix(y_true, y_pred): 绘制混淆矩阵 cm confusion_matrix(y_true, y_pred) plt.figure(figsize(10, 8)) sns.heatmap(cm, annotTrue, fmtd, cmapBlues, cbarFalse, xticklabelsrange(10), yticklabelsrange(10)) plt.xlabel(Predicted Label (Cluster Assigned)) plt.ylabel(True Label) plt.title(Confusion Matrix for K-Means Clustering on MNIST) plt.show() plot_confusion_matrix(y_train, pred_labs)分析混淆矩阵你往往会发现对角线较亮说明大部分样本被正确聚类到对应的数字上。明显的错误模式例如数字“1”可能很少被分错因为它结构简单。而数字“8”可能被分到“3”、“5”或“9”的簇中因为它们在书写上都有环形或弯曲结构。数字“4”和“9”也经常互相混淆。5.3 惯性随迭代次数变化曲线绘制损失函数惯性随迭代次数下降的曲线可以观察算法的收敛过程。# 需要在MyKMeans的fit方法中记录每次迭代后的inertia # 修改fit方法添加 self.inertia_history [] # 在每次迭代结束后计算当前inertia并append到history中 def plot_inertia_history(inertia_history): plt.plot(range(1, len(inertia_history)1), inertia_history, markero) plt.xlabel(Iteration) plt.ylabel(Inertia (Within-cluster SSE)) plt.title(K-Means Inertia Convergence Curve) plt.grid(True) plt.show() # 假设kmeans_model有一个inertia_history属性 plot_inertia_history(kmeans_model.inertia_history)一个健康的曲线应该呈现快速下降后逐渐平缓的趋势。如果曲线震荡剧烈或下降缓慢可能是学习率这里没有显式学习率或初始化有问题。6. 性能瓶颈、常见问题与调优实战6.1 遇到的主要问题与解决方案问题算法运行慢根因在E-step计算距离时如果使用双重循环复杂度是O(n_samples * n_clusters * n_features)对于6万10784的数据量非常慢。解决如前所述必须使用向量化计算。我们的实现利用X[:, np.newaxis, :]进行广播将复杂度降为可高效执行的内存连续操作。这是NumPy编程的核心技巧。问题结果不稳定每次运行准确率差异大根因随机初始化敏感。不同的初始中心可能导致收敛到不同的局部最优解。解决使用K-Means初始化如前所述这是最有效的改进之一。多次运行取最优设置n_init参数像sklearn的KMeans一样用不同的随机种子运行多次选择惯性最小的那次作为最终模型。这需要我们在MyKMeans类外再包一层循环。问题出现空簇Empty Cluster根因在迭代过程中某个簇可能失去所有样本。解决我们在更新步骤中已经做了处理如果发现空簇就随机选择一个数据点作为其新中心。更优雅的方案是选择距离当前所有中心最远的点或者直接采用K-Means的初始化逻辑来重新选择。问题对噪声和异常值敏感根因K-Means使用欧氏距离和均值均值和平方误差对远离中心的点异常值非常敏感。解决对于MNIST这种干净数据问题不大。但在实际应用中可以考虑使用中位数代替均值K-Medoids算法或使用更鲁棒的距离度量。6.2 高级调优技巧特征工程原始的784维像素特征存在大量冗余和噪声。可以尝试降维使用PCA将维度从784降至50-100维既能去除噪声又能大幅提升计算速度有时甚至能提高聚类纯度。from sklearn.decomposition import PCA pca PCA(n_components100, whitenTrue, random_state42) X_pca pca.fit_transform(X_train) # 在X_pca上运行K-Means特征缩放虽然MNIST像素值已经在0-1之间但对于其他数据标准化StandardScaler通常是必须的因为K-Means基于欧氏距离量纲影响巨大。确定最佳K值对于MNIST我们知道数字是0-9所以K10。但在无监督任务中K是未知的。常用方法有肘部法则绘制不同K值对应的惯性曲线选择惯性下降速度突然变缓的点像肘部。轮廓系数计算所有样本的平均轮廓系数越接近1表示聚类效果越好。距离度量欧氏距离是默认选择但对于图像曼哈顿距离、余弦相似度有时可能更合适。修改距离计算部分即可尝试。6.3 与Sklearn的KMeans进行对比验证为了验证我们手写实现的正确性可以用sklearn的官方实现来跑一遍对比惯性值和聚类结果。from sklearn.cluster import KMeans as SKLearnKMeans sk_kmeans SKLearnKMeans(n_clusters10, initk-means, n_init1, max_iter300, tol1e-4, random_state42) sk_kmeans.fit(X_train) print(fSklearn Inertia: {sk_kmeans.inertia_:.2f}) print(fOur Inertia: {kmeans_model.inertia_:.2f}) # 比较两个模型聚类结果的一致性使用调整兰德指数ARI from sklearn.metrics import adjusted_rand_score ari adjusted_rand_score(kmeans_model.labels_, sk_kmeans.labels_) print(fAdjusted Rand Index between our and sklearn labels: {ari:.4f})如果我们的实现正确两个模型的惯性值应该非常接近且ARI指数会很高接近1说明两个模型将样本分到了几乎相同的簇中。7. 实验总结与延伸思考通过这个从零实现的K-Means聚类MNIST实验我们不仅巩固了最优化方法中迭代求解、损失函数最小化的思想更在代码层面实战了向量化编程、算法调试和结果评估的全流程。最终约55%的“分类”准确率虽然不高但它真实地反映了无监督学习在探索数据内在结构方面的能力与局限。几个关键的体会初始化决定上限K-Means带来的提升是实实在在的它用很小的计算开销换来了更稳定、更优的结果。这提醒我们在机器学习中一个好的起点往往事半功倍。评估指标需谨慎用分类准确率来评估聚类结果本身是一种妥协和近似。在真正的无监督场景中我们需要使用轮廓系数、Calinski-Harabasz指数等内部评估指标或者在有部分先验知识时使用调整互信息、调整兰德指数等外部指标。向量化是生命线在数据科学项目中放弃循环拥抱NumPy/SciPy的向量化和广播机制是写出高效代码的不二法门。这次的距离矩阵计算就是一个经典案例。这个项目还可以如何扩展挑战更高维度尝试在CIFAR-10彩色小图像或自己爬取的图片数据集上运行感受高维稀疏空间的“维度灾难”。实现其他聚类算法比如层次聚类、DBSCAN或高斯混合模型并在MNIST上对比它们的效果和特点。与降维技术结合先用t-SNE或UMAP将784维数据可视化到2维观察其自然分布再运行K-Means直观理解“簇”在低维空间中的形态。应用于半监督学习假设我们只有少量带标签数据可以先使用K-Means对大量无标签数据进行聚类然后基于聚类结果为无标签数据打上伪标签再用这些数据扩充训练集训练一个监督模型。这是聚类的一个非常实用的应用场景。本文还有配套的精品资源点击获取