从零实现K-Means++聚类算法：原理、Python实现与MNIST手写数字分类实战

简介本资源是深圳大学计算机软件专业《最优化方法》课程的实验1配套材料面向机器学习初学者与高校算法实践学生聚焦无监督学习中K-Means聚类原理及其在手写数字图像分类任务中的落地实现。资源包共2个文件1个Python源码文件1个Word实验文档总大小873KB结构精炼、即下即用Python脚本完整实现MNIST数据加载、归一化预处理、肘部法则选K、K-Means聚类训练、簇中心可视化及聚类结果与真实标签的匹配分析配套文档则系统梳理实验目标、算法推导、关键代码注释与结果评估逻辑。已有3175人学习下载内容紧扣课程教学要求兼顾理论深度与工程可复现性特别适合理解聚类本质、掌握scikit-learn实战技巧及完成课程实验报告。1. 项目缘起从“手写数字分类”到“无监督学习的敲门砖”最近在整理一些旧项目翻到了当年在深圳大学计算机与软件学院学习《最优化方法》这门课时做的实验报告。实验一的题目是“K-Means聚类之Python实现手写数字图像MNIST分类”现在回头看这个实验设计得相当巧妙。它没有一上来就让我们去调包跑一个分类模型而是用一个看似“简单”的聚类算法去处理一个经典的监督学习数据集这本身就充满了矛盾与挑战。很多刚接触机器学习的同学对MNIST的第一反应就是用卷积神经网络CNN去刷分追求99%以上的准确率。但这个实验反其道而行之它要求我们使用K-Means——一种典型的无监督学习算法去对已知标签的数据进行“分类”。这背后的教学意图非常明确第一让你亲手实现一个最基础、最核心的聚类算法理解其迭代优化的本质这正是《最优化方法》课程的核心第二让你深刻体会“聚类”与“分类”的根本区别理解无监督学习在面对有标签数据时的局限性及其独特价值第三通过Python从零实现打通“算法原理 - 数学推导 - 代码实现 - 结果分析”的完整链路这是夯实基础的关键一步。所以这篇内容不是一份简单的实验报告复现而是结合我后来多年的项目经验对这个实验进行一次深度复盘和扩展。我会带你从零开始用Python纯手工打造K-Means并把它应用到MNIST数据集上。我们会一起探讨为什么用聚类算法做分类听起来像个“歪招”K-Means在处理图像数据时会遇到哪些意想不到的坑如何评估一个聚类结果的好坏以及从这个简单的实验出发我们能窥见机器学习中哪些更深层的逻辑无论你是正在完成类似作业的学生还是想巩固基础的从业者相信这篇近万字的“超详细实验笔记”都能给你带来不少启发。2. K-Means核心原理拆解不止是“找中心点”在动手写代码之前我们必须把K-Means的“里子”和“面子”都搞清楚。很多人对它的理解停留在“随机选K个中心然后不断迭代更新”的层面这远远不够。我们需要从最优化理论的视角重新审视这个算法。2.1 问题形式化一个带约束的优化问题K-Means要解决的核心问题可以这样表述给定一个包含N个数据点的集合 ( X {x_1, x_2, ..., x_N} ) 其中每个 ( x_i ) 是一个D维向量。我们的目标是将这些数据点划分到K个簇cluster( C {C_1, C_2, ..., C_K} ) 中同时找到每个簇的中心 ( \mu {\mu_1, \mu_2, ..., \mu_K} ) 使得所有数据点到其所属簇中心的距离平方和最小。用数学公式表达这就是要最小化以下目标函数也称为畸变函数 Distortion Function [ J \sum_{i1}^{N} \sum_{k1}^{K} r_{ik} | x_i - \mu_k |^2 ] 其中( r_{ik} ) 是一个指示变量如果数据点 ( x_i ) 被分配到簇 ( k ) 则 ( r_{ik} 1 ) 否则为0。并且对于每个 ( i ) 有 ( \sum_{k1}^{K} r_{ik} 1 ) 即一个点只能属于一个簇。注意这里使用的是欧氏距离的平方。使用平方而非绝对距离在数学上会让求导更新中心点的过程变得非常简洁导数就是2倍的距离向量这也是算法高效的关键之一。但这也意味着K-Means对离群点Outliers非常敏感因为平方项会极大地放大远距离点的影响。2.2 迭代优化坐标下降法的经典案例直接同时优化所有 ( r_{ik} ) 和 ( \mu_k ) 是一个NP难问题。K-Means采用的是一种被称为坐标下降法Coordinate Descent的优化策略。它固定一组变量优化另一组变量交替进行分配步E-Step 期望步固定簇中心 ( \mu_k ) 优化分配变量 ( r_{ik} )。 这一步很简单对于每个数据点 ( x_i ) 将其分配给距离最近的簇中心对应的簇。 [ r_{ik} \begin{cases} 1 \text{if } k \arg\min_j | x_i - \mu_j |^2 \ 0 \text{otherwise} \end{cases} ] 这步操作将目标函数 ( J ) 关于 ( r_{ik} ) 的部分最小化了。更新步M-Step 最大化步固定分配变量 ( r_{ik} ) 优化簇中心 ( \mu_k )。 目标函数 ( J ) 对 ( \mu_k ) 求导并令导数为零 [ \frac{\partial J}{\partial \mu_k} -2 \sum_{i1}^{N} r_{ik} (x_i - \mu_k) 0 ] 解得 [ \mu_k \frac{\sum_{i1}^{N} r_{ik} x_i}{\sum_{i1}^{N} r_{ik}} ] 也就是说新的簇中心就是该簇内所有数据点的均值Mean。这就是“K-Means”名字的由来。这两步交替迭代直到簇中心的变化小于某个阈值或者分配不再发生变化。从优化角度看每一步都保证了目标函数 ( J ) 不会增加因此算法最终会收敛到一个局部最优解。2.3 关键特性与局限性为什么是“局部最优”理解K-Means的局限性和掌握其原理同等重要。对初始值敏感由于目标函数非凸算法收敛到的局部最优解严重依赖于初始簇中心的选择。不同的初始化可能导致完全不同的聚类结果和最终的目标函数值。这就是为什么实践中需要多次随机初始化并选择效果最好的一次。需要预先指定K你必须事先告诉算法要分成多少类。在实际的无监督场景中K往往是一个未知的超参数。假设簇为凸形且方差相近K-Means基于距离度量它隐含地假设每个簇是凸形的像一个球体或椭球体并且各个簇的方差大致相同。对于流形分布、环形分布或者大小差异悬殊的簇它的效果会很差。对离群点和噪声敏感如前所述平方误差项使得中心点的计算会被少数极端值拉偏。实操心得在实现时初始化策略是第一个需要精心设计的地方。完全随机选择数据点作为初始中心是最简单的方法但效果不稳定。更常用的方法是K-Means初始化它的核心思想是让初始的簇中心彼此尽可能远离。具体步骤是先随机选一个中心点然后对于每个非中心数据点计算其到已选中心点的最短距离 ( D(x) ) 接着按概率 ( \frac{D(x)^2}{\sum_x D(x)^2} ) 选取下一个中心点。重复此过程直到选满K个中心。这种方法能显著提高聚类质量并减少迭代次数我们后续的代码实现就会采用它。3. 实验环境搭建与MNIST数据预处理工欲善其事必先利其器。这个实验对环境的依赖非常轻量但数据预处理环节却藏着不少细节。3.1 极简Python环境配置你不需要复杂的深度学习框架。核心库就两个NumPy用于高效的矩阵和数学运算Matplotlib用于可视化结果。数据加载我们可以用scikit-learn或torchvision 这里为了通用性选择scikit-learn。# 使用pip安装所需库 pip install numpy matplotlib scikit-learn如果你使用torchvision下载MNIST有时会遇到网络问题导致404错误这也是相关热词中提到的问题。一个可靠的替代方案是使用scikit-learn自带的fetch_openml函数或者直接使用本地已下载的数据集文件。# 一个检查环境并导入库的示例代码块 import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import fetch_openml import warnings warnings.filterwarnings(ignore) # 忽略一些不重要的警告 print(fNumPy version: {np.__version__}) # 确保版本不要太旧即可3.2 MNIST数据加载与关键洞察MNIST数据集包含70000张28x28像素的灰度手写数字图像0-9通常分为60000张训练集和10000张测试集。每张图像被展平成一个784维的向量像素值在0-255之间。def load_mnist_data(): 加载MNIST数据集并返回展平后的像素数据和对应的真实标签。 使用scikit-learn的fetch_openml 它会自动下载数据到本地缓存。 print(正在加载MNIST数据集...) # data_home 可以指定缓存目录 mnist fetch_openml(mnist_784, version1, data_home./data, as_frameFalse) X, y mnist[data], mnist[target] # 将字符串标签转换为整数 y y.astype(np.uint8) print(f数据形状: X{X.shape}, y{y.shape}) print(f像素值范围: [{X.min()}, {X.max()}]) return X, y X_all, y_all load_mnist_data()加载后我们得到X_all的形状是 (70000, 784)y_all是 (70000,)。对于聚类实验我们通常不需要区分训练集和测试集但为了计算效率我们可以先使用一个子集比如前10000个样本。# 使用前N个样本进行实验加快速度 N_SAMPLES 10000 X X_all[:N_SAMPLES] y_true y_all[:N_SAMPLES] print(f使用前 {N_SAMPLES} 个样本进行实验。)3.3 数据预处理为什么以及怎么做原始像素值在0-255之间直接用于基于欧氏距离的K-Means会出问题。因为距离度量对特征的尺度非常敏感。假设有两个像素点一个在图像边缘值常为0一个在笔画中心值可能为255。如果不对数据进行缩放那么笔画中心的细微变化比如从250变到245对距离的贡献可能远大于边缘像素从0到50的变化这显然不合理。我们希望每个像素维度对距离的贡献是“公平”的。最常用的方法是标准化Standardization 也称为Z-score归一化。它将每个特征这里是每个像素点的分布转换为均值为0标准差为1的标准正态分布。def standardize_data(X): 对数据进行标准化处理 mean_val np.mean(X, axis0) std_val np.std(X, axis0) # 防止除零将标准差为0的特征该像素在所有图片中值相同保持不变 std_val[std_val 0] 1.0 X_std (X - mean_val) / std_val return X_std, mean_val, std_val X_std, mean_val, std_val standardize_data(X) print(f标准化后数据形状: {X_std.shape}) print(f标准化后均值近似: {np.mean(X_std, axis0)[:5]}...) print(f标准化后方差近似: {np.var(X_std, axis0)[:5]}...)踩坑记录这里有一个初学者极易忽略的细节。当我们后续用训练好的模型簇中心去预测新数据时必须用训练集计算得到的mean_val和std_val对新数据进行相同的标准化变换。绝对不能在新数据上重新计算均值和标准差否则模型所处的特征空间就变了预测结果将毫无意义。这个原则在几乎所有涉及数据预处理的机器学习流程中都适用。预处理完成后我们可以可视化几张图片感受一下数据。def plot_sample_images(X, y, n10): 随机展示n张图片及其标签 indices np.random.choice(len(X), n, replaceFalse) plt.figure(figsize(15, 2)) for i, idx in enumerate(indices): plt.subplot(1, n, i1) # 注意如果X是标准化后的需要反标准化才能正确显示或者直接用原始X img X[idx].reshape(28, 28) # 假设这里传入的是原始X或已反标准化的X plt.imshow(img, cmapgray) plt.title(fLabel: {y[idx]}) plt.axis(off) plt.tight_layout() plt.show() # 使用原始X或反标准化后的X来显示 X_display X # 或者 (X_std * std_val) mean_val 进行反标准化 plot_sample_images(X_display, y_true, n10)4. 从零实现K-Means算法理解了原理做好了数据准备现在进入核心环节手写K-Means。我们将实现包含K-Means初始化的完整算法。4.1 核心函数设计与实现我们将算法封装成一个类KMeansPlusPlus 这样便于保存模型状态簇中心、标签等。class KMeansPlusPlus: 手动实现的K-Means聚类算法包含K-Means初始化。 def __init__(self, n_clusters8, max_iter300, tol1e-4, random_stateNone): 初始化参数。 :param n_clusters: 簇的数量 (K) :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.centroids None # 簇中心 self.labels None # 每个样本所属簇的标签 self.inertia_ None # 最终的目标函数值距离平方和 def _initialize_centroids(self, X): K-Means 初始化中心点。 1. 随机选择一个点作为第一个中心。 2. 对于每个非中心点计算其到最近中心的距离D(x)。 3. 依据 D(x)^2 的概率分布选择下一个中心。 4. 重复步骤2-3直到选满K个中心。 n_samples, n_features X.shape centroids np.zeros((self.n_clusters, n_features)) rng np.random.RandomState(self.random_state) # 第一步随机选一个 first_idx rng.randint(n_samples) centroids[0] X[first_idx] for k in range(1, self.n_clusters): # 计算每个点到已选中心的最短距离的平方 distances np.array([np.min([np.linalg.norm(x - c) ** 2 for c in centroids[:k]]) for x in X]) # 转换为概率 probabilities distances / distances.sum() # 根据概率分布选择下一个中心点的索引 cumulative_prob probabilities.cumsum() r rng.rand() next_idx np.searchsorted(cumulative_prob, r) centroids[k] X[next_idx] return centroids def fit(self, X): 训练拟合模型。 :param X: 形状为 (n_samples, n_features) 的数组 :return: self n_samples, n_features X.shape # 1. 初始化中心点 self.centroids self._initialize_centroids(X) # 用于记录上一次迭代的标签用于判断收敛 prev_labels np.zeros(n_samples, dtypeint) for iteration in range(self.max_iter): # 2. E-Step: 分配样本到最近的中心 # 计算所有样本到所有中心的距离矩阵 (n_samples, n_clusters) distances np.zeros((n_samples, self.n_clusters)) for k in range(self.n_clusters): # 向量化计算距离比循环快很多 distances[:, k] np.linalg.norm(X - self.centroids[k], axis1) ** 2 self.labels np.argmin(distances, axis1) # 检查是否收敛标签不再变化 if np.array_equal(self.labels, prev_labels): print(f迭代 {iteration}: 标签未变化提前收敛。) break # 3. M-Step: 根据分配重新计算中心点 new_centroids np.zeros_like(self.centroids) for k in range(self.n_clusters): # 找到属于簇k的所有样本 members X[self.labels k] if len(members) 0: new_centroids[k] members.mean(axis0) else: # 如果一个簇没有分配到任何点重新随机初始化该中心避免空簇 new_centroids[k] X[np.random.randint(n_samples)] # 检查中心点变化是否小于容忍度 centroid_shift np.linalg.norm(new_centroids - self.centroids) if centroid_shift self.tol: print(f迭代 {iteration}: 中心点变化 {centroid_shift:.6f} tol 收敛。) break self.centroids new_centroids prev_labels self.labels.copy() if iteration % 20 0: # 计算当前目标函数值惯性 current_inertia np.sum([np.linalg.norm(X[self.labels k] - self.centroids[k]) ** 2 for k in range(self.n_clusters)]) print(f迭代 {iteration}, 惯性值: {current_inertia:.2f}) # 计算最终的惯性值 self.inertia_ np.sum([np.linalg.norm(X[self.labels k] - self.centroids[k]) ** 2 for k in range(self.n_clusters)]) print(f训练完成共迭代 {iteration1} 次最终惯性值: {self.inertia_:.2f}) return self def predict(self, X): 预测新样本的簇标签。 :param X: 形状为 (n_samples, n_features) 的数组 :return: 预测的簇标签数组 n_samples X.shape[0] distances np.zeros((n_samples, self.n_clusters)) for k in range(self.n_clusters): distances[:, k] np.linalg.norm(X - self.centroids[k], axis1) ** 2 return np.argmin(distances, axis1)4.2 代码实现的几个关键细节与避坑点距离计算的向量化在fit和predict函数中我们使用了np.linalg.norm(X - self.centroids[k], axis1)来计算一个样本矩阵与一个中心点之间的欧氏距离。axis1表示对每个样本计算范数距离。这种向量化操作比用for循环遍历每个样本要快几个数量级尤其是在处理像MNIST这样的大数据集时。空簇的处理在M-Step中有可能出现某个簇没有被分配到任何样本的情况特别是初始化不好或K值过大时。我们的代码中加入了检查if len(members) 0。如果为空我们选择随机选择一个数据点作为该簇的新中心。这是一种简单的处理策略。其他策略包括删除该簇减少K值、将离其最近的点分配给它、或者选择距离当前所有中心最远的点作为新中心。收敛条件我们设置了双重收敛判断一是中心点的移动距离小于阈值tol二是样本的簇分配不再发生变化。通常标签不变是一个更强的条件满足时一定收敛。在实际中由于浮点数计算精度有时中心点会微小震荡用tol判断可以避免无限循环。惯性值Inertiaself.inertia_存储了最终的目标函数值即所有样本到其所属簇中心的距离平方和。这个值越小说明聚类内部越紧凑。它是评估聚类效果和选择K值的重要内部指标。实操心得在实现K-Means初始化时计算每个点到已选中心的最短距离平方distances是性能瓶颈。上面的实现用了列表推导式对于大数据集可能较慢。一个更高效的优化是维护一个(n_samples,)的数组min_distances 每次新增一个中心后只更新那些到新中心距离更小的点的min_distances值而不是全部重新计算。这对于追求极致性能的场景是必要的优化但为了代码清晰上面的示例保持了直观的实现。5. 在MNIST上运行K-Means并分析结果现在让我们把实现的算法应用到处理好的MNIST数据上。我们已知手写数字有0-9共10个类别所以这里设置n_clusters10。5.1 训练模型与基础可视化# 设置随机种子保证结果可复现 RANDOM_STATE 42 # 实例化并训练模型 kmeans KMeansPlusPlus(n_clusters10, max_iter300, tol1e-4, random_stateRANDOM_STATE) kmeans.fit(X_std) # 使用标准化后的数据 # 获取预测标签 y_pred kmeans.labels print(f聚类结果标签示例前20个: {y_pred[:20]}) print(f真实标签示例前20个: {y_true[:20]})运行后你会看到迭代过程。由于K-Means的随机性每次运行的迭代次数和最终惯性值可能略有不同。接下来我们可视化一下学习到的簇中心。在图像聚类中簇中心可以看作是一个“平均数字”或“原型数字”。def plot_centroids(centroids, mean_val, std_val, n_cols5): 绘制簇中心图像。由于中心点是在标准化空间中学到的 需要反标准化到原始像素空间才能正确显示。 n_clusters centroids.shape[0] n_rows (n_clusters n_cols - 1) // n_cols # 计算需要多少行 plt.figure(figsize(2*n_cols, 2*n_rows)) # 反标准化centroids_original centroids * std_val mean_val centroids_original centroids * std_val mean_val # 将像素值裁剪到0-255的合理范围 centroids_original np.clip(centroids_original, 0, 255) for i in range(n_clusters): plt.subplot(n_rows, n_cols, i1) # 重塑为28x28图像 centroid_img centroids_original[i].reshape(28, 28) plt.imshow(centroid_img, cmapgray) plt.title(fCentroid {i}) plt.axis(off) plt.tight_layout() plt.show() plot_centroids(kmeans.centroids, mean_val, std_val, n_cols5)观察这些中心点图像你会发现一些有趣的现象有些中心点看起来像一个清晰的数字比如一个明显的“1”或“0”有些则像是几个数字的模糊混合体比如“4”和“9”可能被混在一起。这直观地反映了K-Means在MNIST上的局限性它基于像素距离的相似性进行划分但不同人写的“7”和“1”在像素空间上可能很接近而同一个“4”的不同写法可能相差很远。5.2 聚类效果评估当没有真实标签时怎么办在真正的无监督学习中我们是没有y_true的。那么如何评价聚类结果的好坏呢我们依赖内部评估指标和外部评估指标当有真实标签时用于验证。内部评估指标不依赖真实标签仅基于数据本身的紧凑性和分离性。惯性Inertia我们已经计算了越小越好。但它会随着K增大而单调减小所以不能单独用来确定K。轮廓系数Silhouette Score结合了簇内凝聚度和簇间分离度。对于样本 ( i ) 其轮廓系数 ( s(i) ) 计算如下( a(i) )样本 ( i ) 到同簇其他样本的平均距离簇内不相似度。( b(i) )样本 ( i ) 到其他所有簇中样本的平均距离的最小值簇间不相似度。( s(i) \frac{b(i) - a(i)}{\max(a(i), b(i))} ) ( s(i) ) 的取值范围为[-1, 1]。越接近1说明样本聚类越合理越接近-1说明可能被分错了簇接近0则说明样本在两个簇的边界上。所有样本轮廓系数的平均值即为整体轮廓系数。from sklearn.metrics import silhouette_score, silhouette_samples import matplotlib.cm as cm # 计算整体轮廓系数 silhouette_avg silhouette_score(X_std, y_pred) print(f整体轮廓系数: {silhouette_avg:.4f}) # 可以绘制轮廓分析图更详细地查看每个簇的情况 def plot_silhouette(X, y_pred, n_clusters): fig, (ax1, ax2) plt.subplots(1, 2, figsize(12, 5)) ax1.set_xlim([-0.1, 1]) # ax1.set_ylim([0, len(X) (n_clusters 1) * 10]) silhouette_vals silhouette_samples(X, y_pred) y_lower 10 for k in range(n_clusters): kth_cluster_silhouette_vals silhouette_vals[y_pred k] kth_cluster_silhouette_vals.sort() size_cluster_k kth_cluster_silhouette_vals.shape[0] y_upper y_lower size_cluster_k color cm.nipy_spectral(float(k) / n_clusters) ax1.fill_betweenx(np.arange(y_lower, y_upper), 0, kth_cluster_silhouette_vals, facecolorcolor, edgecolorcolor, alpha0.7) ax1.text(-0.05, y_lower 0.5 * size_cluster_k, str(k)) y_lower y_upper 10 ax1.set_xlabel(轮廓系数值) ax1.set_ylabel(簇标签) ax1.axvline(xsilhouette_avg, colorred, linestyle--) ax1.set_title(f轮廓分析图 (K{n_clusters}, avg{silhouette_avg:.3f})) # 在右侧子图展示簇中心图像反标准化后 # ... (此处省略中心点绘制代码可参考之前的plot_centroids函数) plt.tight_layout() plt.show() plot_silhouette(X_std, y_pred, 10)轮廓分析图可以直观显示每个簇的轮廓系数分布情况。红色虚线是整体平均值。理想情况下每个簇的“条形图”都应该较宽且大部分在平均线右侧。外部评估指标由于我们有真实标签y_true 可以用来验证聚类结果与真实分类的匹配程度。但注意聚类标签是算法自己赋予的0-9与真实数字标签0-9没有对应关系。我们需要找到一个最佳的映射关系将聚类标签映射到真实标签然后计算准确率等指标。这通常通过混淆矩阵和匈牙利算法或简单枚举来找到最佳匹配。from sklearn.metrics import confusion_matrix, accuracy_score import scipy.optimize def cluster_to_label_mapping(y_true, y_pred): 使用匈牙利算法找到聚类标签到真实标签的最佳一一映射。 目标是最大化映射后的整体准确率。 n_clusters len(np.unique(y_pred)) n_classes len(np.unique(y_true)) # 构建权重矩阵这里用负的计数因为匈牙利算法解决的是最小化问题 weight_matrix np.zeros((n_clusters, n_classes)) for i in range(n_clusters): for j in range(n_classes): weight_matrix[i, j] -np.sum((y_pred i) (y_true j)) # 使用线性求和分配匈牙利算法 row_ind, col_ind scipy.optimize.linear_sum_assignment(weight_matrix) # 创建映射字典 mapping {row_ind[i]: col_ind[i] for i in range(len(row_ind))} # 根据映射转换预测标签 y_pred_mapped np.array([mapping[label] for label in y_pred]) return y_pred_mapped, mapping y_pred_mapped, mapping_dict cluster_to_label_mapping(y_true, y_pred) print(f聚类标签到真实标签的映射: {mapping_dict}) accuracy accuracy_score(y_true, y_pred_mapped) print(f映射后的聚类准确率: {accuracy:.4f}) # 绘制混淆矩阵 cm confusion_matrix(y_true, y_pred_mapped) plt.figure(figsize(10, 8)) plt.imshow(cm, interpolationnearest, cmapplt.cm.Blues) plt.title(混淆矩阵 (聚类 vs 真实标签)) plt.colorbar() tick_marks np.arange(10) plt.xticks(tick_marks, range(10)) plt.yticks(tick_marks, range(10)) plt.ylabel(真实标签) plt.xlabel(映射后的聚类标签) # 在格子中填充数字 thresh cm.max() / 2. for i in range(cm.shape[0]): for j in range(cm.shape[1]): plt.text(j, i, format(cm[i, j], d), horizontalalignmentcenter, colorwhite if cm[i, j] thresh else black) plt.tight_layout() plt.show()通过混淆矩阵你可以清晰地看到哪些数字容易被K-Means混淆。例如“4”和“9”、“3”和“8”、“5”和“6”常常是“重灾区”。这是因为它们在像素空间上的形状确实有相似之处。5.3 如何确定最佳的K值肘部法则与轮廓分析在真实的无监督任务中K是未知的。如何选择两个最常用的方法是肘部法则和轮廓系数法。def find_optimal_k(X, k_rangerange(2, 16)): 通过肘部法则和轮廓系数寻找最佳K值。 inertias [] silhouette_scores [] for k in k_range: print(f正在训练 K{k}...) kmeans KMeansPlusPlus(n_clustersk, max_iter100, random_stateRANDOM_STATE) kmeans.fit(X) inertias.append(kmeans.inertia_) # 计算轮廓系数可能较慢对于大数据集可以抽样计算 if len(X) 5000: sample_idx np.random.choice(len(X), 5000, replaceFalse) X_sample X[sample_idx] labels_sample kmeans.predict(X_sample) score silhouette_score(X_sample, labels_sample) else: score silhouette_score(X, kmeans.labels_) silhouette_scores.append(score) print(f K{k}, Inertia{kmeans.inertia_:.2f}, Silhouette{score:.4f}) # 绘制肘部法则图 fig, (ax1, ax2) plt.subplots(1, 2, figsize(14, 5)) ax1.plot(k_range, inertias, bo-) ax1.set_xlabel(K值) ax1.set_ylabel(惯性 (Inertia)) ax1.set_title(肘部法则) ax1.grid(True) # 绘制轮廓系数图 ax2.plot(k_range, silhouette_scores, ro-) ax2.set_xlabel(K值) ax2.set_ylabel(轮廓系数 (Silhouette Score)) ax2.set_title(轮廓系数法) ax2.grid(True) plt.tight_layout() plt.show() # 找出轮廓系数最大的K best_k_by_silhouette k_range[np.argmax(silhouette_scores)] print(f根据轮廓系数建议的K值为: {best_k_by_silhouette}) return inertias, silhouette_scores # 注意在全量10000个样本上运行多个K值会很慢可以先用子样本测试 X_sample_for_k X_std[:2000] # 用2000个样本来寻找K inertias, scores find_optimal_k(X_sample_for_k, k_rangerange(2, 15))肘部法则观察惯性随K变化的曲线。随着K增大惯性会下降。曲线拐点像手肘一样对应的K值通常被认为是增加K的收益开始变小的点即一个合理的K值。轮廓系数法选择轮廓系数最大的K值。轮廓系数越高说明聚类效果越好。对于MNIST我们已知有10个数字所以K10应该是一个合理的选择。但通过上述分析你可能会发现轮廓系数在K10时未必是最高这可能因为K-Means本身对这类数据的假设并不完全成立。6. 深入探讨K-Means在图像聚类中的局限与改进思路通过前面的实验我们已经直观感受到了K-Means处理MNIST的“力不从心”。准确率可能只有50%-60%远低于有监督学习。我们来系统性地分析一下原因并探讨可能的改进方向。6.1 为什么K-Means对MNIST分类效果不佳特征表示问题K-Means直接在原始像素空间操作。784维的像素空间非常高维且稀疏而且欧氏距离在高维空间中会失效所谓的“维数灾难”。两个数字的相似性不能简单地用像素差的平方和来衡量。例如一个稍微倾斜的“1”和一个垂直的“1”在像素空间的距离可能很大但人眼看来它们是同一个数字。数据分布假设不匹配K-Means假设簇是球形的、各向同性的。但MNIST中数字的分布可能是复杂的流形。例如数字“2”的各种写法可能构成一个弯曲的“带”而不是一个紧致的球。距离度量单一欧氏距离对图像的平移、旋转、缩放等变化非常敏感。而人类识别数字对这些变化具有一定的不变性。类别不平衡与多模态同一个数字如“7”可能有多种完全不同的写法有横杠的和没有横杠的这在一个簇内形成了多模态分布而K-Means的一个簇中心只能代表一个“模式”。6.2 可能的改进方向与实践建议特征工程/降维主成分分析PCA在应用K-Means之前先使用PCA将784维数据降到较低维度如50维或100维。这可以去除噪声和冗余信息同时缓解维数灾难。在实践中这常常能提升聚类效果和速度。from sklearn.decomposition import PCA pca PCA(n_components50, random_stateRANDOM_STATE) X_pca pca.fit_transform(X_std) # 然后在 X_pca 上运行K-Means自编码器Autoencoder一种更强大的非线性降维方法可以学习到比PCA更有效的低维特征表示这些特征对数字的识别可能更有帮助。使用更先进的聚类算法谱聚类Spectral Clustering特别适合发现非凸形状的簇。它先对数据点构建相似度图然后对图进行切割。对于像MNIST这样可能存在复杂流形结构的数据谱聚类往往比K-Means效果好得多。DBSCAN基于密度的聚类不需要预先指定K能发现任意形状的簇并能识别噪声点。但对于MNIST这种密度差异可能不大的数据集参数调优会比较困难。改进的距离度量如果坚持使用K-Means可以考虑替换欧氏距离。例如余弦相似度对于文本或某些高维特征可能更有效因为它关注向量的方向而非绝对距离。但对于MNIST原始像素效果提升有限。集成与后处理多次运行K-Means不同初始化然后通过集成方法如投票来确定最终的簇标签可以提高稳定性。对聚类结果进行后处理例如将非常大的簇进行分裂或将非常接近的簇进行合并。一个简单的对比实验PCA K-Means# 使用PCA降维到50维后再进行K-Means聚类 print( 实验PCA降维后K-Means ) pca PCA(n_components50, random_stateRANDOM_STATE) X_pca pca.fit_transform(X_std) kmeans_pca KMeansPlusPlus(n_clusters10, random_stateRANDOM_STATE) kmeans_pca.fit(X_pca) y_pred_pca kmeans_pca.labels_ # 评估 y_pred_pca_mapped, _ cluster_to_label_mapping(y_true, y_pred_pca) accuracy_pca accuracy_score(y_true, y_pred_pca_mapped) silhouette_pca silhouette_score(X_pca, y_pred_pca) print(fPCAK-Means 准确率: {accuracy_pca:.4f}) print(fPCAK-Means 轮廓系数: {silhouette_pca:.4f}) print(f原始像素K-Means 准确率: {accuracy:.4f}) print(f原始像素K-Means 轮廓系数: {silhouette_avg:.4f})在我的多次实验中PCAK-Means的准确率通常能比原始像素的K-Means高出几个百分点这验证了特征预处理的重要性。7. 实验总结与延伸思考回顾整个实验我们从最优化理论的坐标下降法出发手动实现了K-Means算法并将其应用于经典的MNIST手写数字数据集。我们经历了数据加载、预处理、模型实现、训练、评估、可视化和分析改进的全流程。这个实验的价值远不止于完成一次聚类。它像一面镜子清晰地照出了无监督学习与有监督学习的根本差异也暴露了简单模型在复杂现实数据面前的局限性。K-Means在这里更像一个“探针”帮助我们理解MNIST数据在原始像素空间中的结构。几点重要的实操心得初始化决定下限数据决定上限K-Means初始化能极大提升算法的稳定性和结果质量这是必须掌握的技巧。但即使初始化再好如果数据本身的特征表示不合适如原始像素算法的上限也不会高。评估指标要组合使用惯性、轮廓系数、外部准确率如果有标签、以及最终的可视化中心点图像、混淆矩阵需要综合来看。单一指标可能会产生误导。聚类结果的可解释性在图像聚类中可视化簇中心是理解模型“学到了什么”的绝佳方式。如果中心点看起来是毫无意义的噪声那说明聚类很可能是失败的。K-Means是起点不是终点它是最基础、最直观的聚类算法理解它对于学习更复杂的模型如高斯混合模型GMM、谱聚类有莫大帮助。但也要清楚它的边界知道在什么情况下应该寻求更强大的工具。最后虽然在这个特定任务上K-Means的表现远不如有监督的深度学习模型但这个实验过程的训练价值丝毫不亚于训练一个高精度的CNN模型。它强迫你去思考数据的本质、距离的定义、优化的过程以及模型评估的复杂性。这些思考是构建扎实机器学习知识体系的基石。下次当你轻松调用sklearn.cluster.KMeans时希望你能想起手动实现时遇到的每一个细节和踩过的每一个坑这才是这个实验留给你的真正财富。本文还有配套的精品资源点击获取