尧图网站建设 
简介本资源是深圳大学计算机软件专业《最优化方法》课程配套实验材料面向机器学习初学者与高校算法实践者聚焦无监督学习核心任务——利用K-Means聚类实现MNIST手写数字图像的自动分组与结构发现。资源包共2个文件1个可运行Python脚本1份详细实验文档总大小873KB脚本完整实现数据加载、归一化预处理、肘部法则选K、聚类执行与结果可视化文档则系统梳理算法原理、关键参数设计依据、scikit-learn调用细节及聚类效果评估方法。已有3177人学习下载内容紧扣教学实验场景代码模块清晰、注释详尽附带可直接复现的绘图与分簇分析逻辑特别适合理解K-Means在真实图像数据上的行为边界与应用局限为后续监督分类打下扎实的无监督建模基础。1. K-Means 聚类真能分清手写数字——别急着调sklearn.cluster.KMeans先搞懂它在 MNIST 上为什么“看起来准、实际懵”你跑通了KMeans(n_clusters10).fit(X_train)聚类结果轮廓系数 0.12每个簇中心画出来像模像样甚至用匈牙利算法匹配标签后准确率冲到 53%——恭喜你成功复现了深大计软《最优化方法》实验1的经典翻车现场。这不是代码错了而是你还没意识到K-Means 在 MNIST 上不是分类器它是个“无监督的形状归组器”。它不认“3 是三”只认“这一坨像素和那一坨像素离得近”。当数字“1”和“7”在像素空间里挤在一起比如都写得细长倾斜K-Means 会毫不犹豫把它们划进同一个簇而“4”和“9”如果都带圆弧竖线也可能被强行拆散。这不是算法缺陷是欧氏距离 球形簇假设 无标签先验的必然结果。本实验的核心价值从来不是“用聚类替代分类”而是亲手推导目标函数、观察质心迭代如何卡在局部极小、验证维度灾难对距离失效的影响——这才是最优化方法课要锤你的点。适合人群刚学完梯度下降、拉格朗日乘子想把抽象公式落到像素矩阵上的计软本科生或想搞懂“为什么无监督评估指标总比有监督低一截”的算法初学者。别指望靠它上 Kaggle但如果你能手动实现 E-step/M-step 并画出每轮质心漂移轨迹你对“迭代优化”的理解就落地了。2. 从零推导 K-Means 目标函数为什么最小化平方误差等价于最大似然K-Means 表面是“找 10 个中心点让所有样本离最近中心最近”但它的数学根基藏在概率模型里。我们不直接讲 EM 算法而是从最优化视角拆解它本质是在求解一个带隐变量的非凸优化问题。先明确符号设训练集 $X {x^{(1)}, x^{(2)}, ..., x^{(N)}}$其中 $x^{(i)} \in \mathbb{R}^D$MNIST 是 784 维簇中心为 $\mu_1, \mu_2, ..., \mu_K \in \mathbb{R}^D$隐变量 $z^{(i)} \in {1,2,...,K}$ 表示第 $i$ 个样本所属簇。K-Means 的目标函数是$$ J(\mu, z) \sum_{i1}^{N} \sum_{k1}^{K} \mathbb{1}(z^{(i)} k) \cdot |x^{(i)} - \mu_k|^2 $$这个式子直白说就是对每个样本只算它到自己所属簇中心的距离平方再全部加起来。注意两点第一$z^{(i)}$ 是离散变量无法求导第二$J$ 关于 $\mu_k$ 是凸的但关于 $z^{(i)}$ 完全非凸。所以标准解法是坐标下降Coordinate Descent固定 $z$ 更新 $\mu$再固定 $\mu$ 更新 $z$交替进行直到收敛。这正是 K-Means 的两步迭代逻辑。2.1 E-step给定质心分配样本到最近簇硬分配这步不需要概率纯几何。对每个样本 $x^{(i)}$计算它到所有 $K$ 个质心的欧氏距离平方选最小的那个# 手动实现 E-step不用 sklearn import numpy as np def e_step(X, mu): X: (N, D) 样本矩阵 mu: (K, D) 质心矩阵 返回: (N,) 每个样本的簇索引数组 N, D X.shape K, _ mu.shape # 向量化计算所有距离平方避免 for 循环 # 利用 (a-b)^2 a^2 - 2ab b^2 展开 X_sq np.sum(X**2, axis1, keepdimsTrue) # (N, 1) mu_sq np.sum(mu**2, axis1, keepdimsTrue) # (K, 1) cross 2 * X mu.T # (N, K) # dist_sq[i, k] ||x_i - mu_k||^2 dist_sq X_sq - cross mu_sq.T # (N, K) return np.argmin(dist_sq, axis1) # (N,) # 验证用 MNIST 前 1000 张图测试 from sklearn.datasets import fetch_openml mnist fetch_openml(mnist_784, version1, as_frameFalse, parserauto) X, y mnist.data.astype(float32), mnist.target.astype(int) X_sample X[:1000] / 255.0 # 归一化到 [0,1] y_sample y[:1000] # 随机初始化 10 个质心均匀分布 np.random.seed(42) mu_init np.random.rand(10, 784) z_pred e_step(X_sample, mu_init) print(f初始分配簇大小分布 {np.bincount(z_pred, minlength10)}) # 输出类似[112 98 105 92 101 95 103 97 99 98] —— 基本均匀参数说明e_step中的dist_sq计算用了向量化技巧避免for i in range(N): for k in range(K):的 O(NKD) 复杂度。关键在X_sq - cross mu_sq.T这一行X_sq是每个样本自身模长平方mu_sq.T是每个质心模长平方转置cross是内积矩阵。这样一次广播运算就得到全部距离平方速度提升 100 倍以上。新手常在这里写嵌套循环跑 MNIST 直接卡死。2.2 M-step给定分配更新质心为簇内均值解析解这步是目标函数 $J$ 关于 $\mu_k$ 的最小化。对固定 $k$只看属于簇 $k$ 的样本集合 $C_k {i | z^{(i)} k}$则$$ \frac{\partial J}{\partial \mu_k} -2 \sum_{i \in C_k} (x^{(i)} - \mu_k) 0 \quad \Rightarrow \quad \mu_k \frac{1}{|C_k|} \sum_{i \in C_k} x^{(i)} $$看到没质心更新就是取簇内样本均值——这是唯一解析解不用梯度下降。这也是 K-Means 收敛快的原因M-step 是闭式解。def m_step(X, z, K): X: (N, D) 样本 z: (N,) 簇索引 K: 簇数 返回: (K, D) 新质心矩阵 N, D X.shape mu_new np.zeros((K, D)) for k in range(K): # 找出属于簇 k 的所有样本 mask (z k) if np.sum(mask) 0: # 防止空簇 mu_new[k] np.mean(X[mask], axis0) else: # 空簇处理重采样一个随机样本作为新质心常见策略 mu_new[k] X[np.random.randint(0, N)] return mu_new # 测试 M-step z_init e_step(X_sample, mu_init) mu_updated m_step(X_sample, z_init, K10) print(f更新后质心形状: {mu_updated.shape}) # (10, 784)关键细节m_step中必须处理空簇empty cluster。如果某轮 E-step 后某个簇没分到任何样本np.mean会报错或返回全零向量导致后续迭代崩溃。这里采用“重采样随机样本”策略也有用“距其他质心最远点”或“最大方差方向切分”的但对 MNIST 来说随机重采样最稳定。这是实操中第一个血泪经验永远检查np.bincount(z)是否有 0。3. 手动实现完整 K-Means 迭代监控损失下降、可视化质心漂移、验证收敛性现在把 E-step 和 M-step 串起来加上收敛判断和日志。重点不是“跑出结果”而是看见优化过程本身——这才是最优化方法课的灵魂。我们监控三个量目标函数值 $J$、质心移动距离、簇分配变化率。3.1 主循环带收敛判断与日志的迭代框架def kmeans_manual(X, K, max_iters100, tol1e-4, random_state42): 手动实现 K-Means 完整流程 X: (N, D) 归一化后的样本MNIST 已除以 255 K: 簇数对 MNIST 固定为 10 max_iters: 最大迭代次数 tol: 损失变化容忍阈值 返回: mu_history (list of (K,D)), z_history (list of (N,)), losses (list of float) np.random.seed(random_state) N, D X.shape # 初始化质心用 K-means 策略比随机好得多 mu np.zeros((K, D)) # 第一个质心随机选 mu[0] X[np.random.randint(0, N)] # 后续质心按距离平方概率选 for k in range(1, K): # 计算所有点到已选质心的最小距离平方 dist_sq_to_mu np.full(N, np.inf) for i in range(k): dist_sq np.sum((X - mu[i])**2, axis1) dist_sq_to_mu np.minimum(dist_sq_to_mu, dist_sq) # 概率正比于 dist_sq_to_mu probs dist_sq_to_mu / np.sum(dist_sq_to_mu) mu[k] X[np.random.choice(N, pprobs)] mu_history [mu.copy()] z_history [] losses [] for it in range(max_iters): # E-step z e_step(X, mu) z_history.append(z) # 计算当前损失 J loss 0.0 for k in range(K): mask (z k) if np.sum(mask) 0: loss np.sum((X[mask] - mu[k])**2) losses.append(loss) # M-step mu_new m_step(X, z, K) # 检查收敛质心移动距离 tol mu_shift np.max(np.sqrt(np.sum((mu_new - mu)**2, axis1))) mu_history.append(mu_new.copy()) mu mu_new # 打印进度 if it % 10 0 or it max_iters-1: print(fIter {it:3d} | Loss: {loss:.2e} | Max shift: {mu_shift:.4f}) if mu_shift tol: print(fConverged at iteration {it}) break return mu_history, z_history, losses # 在 MNIST 子集上运行 mu_hist, z_hist, losses kmeans_manual(X_sample, K10, max_iters50, tol1e-3)为什么用 K-means 初始化随机初始化可能导致质心全挤在数字“1”的区域其他数字如“8”、“6”永远分不到簇。K-means 通过距离加权采样强制质心分散通常减少 30% 迭代次数。dist_sq_to_mu的计算是核心对每个未选点算它到所有已选质心的最小距离平方再以此为权重抽样。这是工业级实现的标配不是炫技。3.2 可视化质心演化从噪声到“数字雏形”的 30 轮旅程K-Means 的质心不是静态图片是动态优化的产物。我们把每轮质心 reshape 成 28x28 并画出来import matplotlib.pyplot as plt def plot_centroids_evolution(mu_history, n_cols5, figsize(12, 8)): 绘制质心随迭代的变化 mu_history: list of (K, 784) 质心矩阵 n_iters len(mu_history) n_rows (n_iters n_cols - 1) // n_cols fig, axes plt.subplots(n_rows, n_cols, figsizefigsize) axes axes.flatten() if n_iters 1 else [axes] for i, mu in enumerate(mu_history): if i len(axes): break # 取第一个质心索引 0为例展示其变化 # 实际可循环画所有 10 个但这里简化 img mu[0].reshape(28, 28) axes[i].imshow(img, cmapgray) axes[i].set_title(fIter {i}) axes[i].axis(off) # 隐藏多余子图 for j in range(i1, len(axes)): axes[j].remove() plt.tight_layout() plt.show() # 运行绘图 plot_centroids_evolution(mu_hist)你会看到第 0 轮是随机噪声块第 5 轮开始出现模糊的竖线/横线第 15 轮能辨认出“1”、“0”的轮廓第 30 轮基本稳定但和真实数字仍有差距——因为 K-Means 学的是像素均值不是笔画结构。这就是最优化的真相它在约束下逼近局部最优而非生成完美图像。3.3 监控损失曲线与收敛诊断为什么有时损失不降反升画出losses曲线plt.figure(figsize(10, 4)) plt.plot(losses, b-o, markersize3) plt.xlabel(Iteration) plt.ylabel(Objective J (Sum of Squared Errors)) plt.title(K-Means Loss Curve on MNIST Subset) plt.grid(True) plt.show()正常情况是单调下降。但如果出现震荡或上升一定是代码 bug。常见原因e_step中距离计算错误比如忘了平方或用了曼哈顿距离m_step中空簇未处理导致mu_new[k]为全零下轮e_step计算距离时||x_i - 0||^2极大初始化质心超出数据范围如用了np.random.randn未归一化收敛性验证口诀K-Means 保证每轮 $J$ 不增但不保证全局最优。若losses曲线非单调立刻检查e_step的dist_sq公式——90% 的问题出在这里。4. 避坑MNIST 上 K-Means 的 4 个经典翻车现场与急救方案K-Means 在 MNIST 上不是不能跑而是处处是坑。这些坑不来自算法本身而来自数据特性与实现细节的碰撞。以下是我在深大计软助教三年批改上百份实验报告总结的最高频、最隐蔽、最致命的四个问题每条都附带现象、根因和一行修复代码。4.1 现象迭代 50 轮后某个簇的样本数为 0后续m_step报RuntimeWarning: Mean of empty slice原因E-step 分配时若某质心离所有样本都远可能没分到任何点M-step 中np.mean(X[mask])对空数组返回nan污染后续计算。解决在m_step中强制检查空簇并用安全策略重置质心。不要用np.nanmean它返回nan不解决问题# 错误示范引发连锁 nan # mu_new[k] np.mean(X[mask], axis0) # mask 全 False → nan # 正确修复在 m_step 函数内 if np.sum(mask) 0: # 方案1重采样随机样本推荐简单鲁棒 mu_new[k] X[np.random.randint(0, N)] # 方案2用所有样本的均值更稳定但可能偏移 # mu_new[k] np.mean(X, axis0)4.2 现象质心图像全是灰色块像素值集中在 0.4~0.6完全看不出数字形状原因MNIST 像素是uint80~255但你没归一化KMeans对量纲敏感未归一化的像素值0~255会让距离计算被高亮区域主导质心被拉向平均灰度。解决必须在输入前除以 255缩放到 [0,1]# 错误直接喂原始数据 # X_raw mnist.data # shape (70000, 784), dtype uint8 # 正确归一化是铁律 X_normalized mnist.data.astype(float32) / 255.0 # 验证print(X_normalized.min(), X_normalized.max()) → 应输出 0.0 1.04.3 现象e_step运行极慢10 秒CPU 占用 100%dist_sq计算卡死原因写了双重 for 循环计算距离复杂度 O(NKD)MNIST N70000, K10, D784 → 54.88 亿次运算。解决用向量化公式X_sq - 2*Xmu.T mu_sq.T一行替代循环# 错误慢如蜗牛 # dist_sq np.zeros((N, K)) # for i in range(N): # for k in range(K): # dist_sq[i,k] np.sum((X[i] - mu[k])**2) # 正确毫秒级 X_sq np.sum(X**2, axis1, keepdimsTrue) # (N, 1) mu_sq np.sum(mu**2, axis1, keepdimsTrue) # (K, 1) dist_sq X_sq - 2 * X mu.T mu_sq.T # (N, K) 广播4.4 现象用匈牙利算法匹配簇标签后准确率只有 10%~20%远低于预期原因你匹配的是“簇 ID”和“数字标签”但 K-Means 的簇 ID 是任意的0~9而数字标签 0~9 有语义。若簇 0 匹配到数字 5簇 1 匹配到数字 0直接按 ID 对齐当然错。必须用二分图匹配Hungarian Algorithm找最优映射。解决用scipy.optimize.linear_sum_assignment计算最佳匹配from scipy.optimize import linear_sum_assignment from sklearn.metrics import confusion_matrix def calculate_matching_accuracy(z_pred, y_true, K10): z_pred: (N,) 预测簇标签 y_true: (N,) 真实数字标签0~9 返回: 匹配后的准确率 # 构建混淆矩阵行簇ID列真实标签 cm confusion_matrix(y_true, z_pred, labelsrange(K)) # Hungarian 算法找最大匹配 row_ind, col_ind linear_sum_assignment(-cm) # 最大化故取负 # 计算匹配总正确数 total_correct cm[row_ind, col_ind].sum() return total_correct / len(y_true) acc calculate_matching_accuracy(z_hist[-1], y_sample) print(fMatching Accuracy: {acc:.3f}) # 正常应 0.5~0.6提示匹配准确率 50%~60% 是 MNIST 上 K-Means 的合理上限。别追求 95%——那属于有监督学习的领域。K-Means 的价值在于无标签下的结构发现不是替代分类器。5. 进阶验证用轮廓系数、Calinski-Harabasz 指标量化聚类质量并对比 K5 vs K10K-Means 的“好坏”不能只看准确率毕竟它没标签。我们要用无监督评估指标来回答当前 K10 是最优的吗质心真的分开了吗这里用两个黄金指标轮廓系数Silhouette Score和 Calinski-Harabasz 指数。5.1 轮廓系数衡量“簇内紧密度 vs 簇间分离度”对每个样本 $i$定义$a(i)$$i$ 到同簇其他点的平均距离簇内不相似度$b(i)$$i$ 到最近其他簇所有点的平均距离簇间不相似度轮廓值 $s(i) \frac{b(i)-a(i)}{\max(a(i),b(i))} \in [-1,1]$$s(i)$ 接近 1 表示样本聚得好接近 -1 表示分错簇。平均轮廓系数越高聚类越优。from sklearn.metrics import silhouette_score # 计算最终聚类的轮廓系数 silhouette_avg silhouette_score(X_sample, z_hist[-1]) print(fAverage Silhouette Score: {silhouette_avg:.3f}) # 典型值K10 时约 0.10~0.15K5 时可能升到 0.18因为簇更大更易分离为什么 MNIST 的轮廓系数这么低因为数字“4”和“9”、“3”和“8”在像素空间本就相似都含圆弧直线K-Means 无法用线性边界分开它们。低轮廓系数恰恰证明了数据本身的可分性有限不是算法不行。5.2 Calinski-Harabasz 指数簇间离散度 / 簇内离散度公式$CH \frac{Tr(B_k)}{Tr(W_k)} \times \frac{N-K}{K-1}$其中 $B_k$ 是簇间散度矩阵$W_k$ 是簇内散度矩阵。值越大越好。from sklearn.metrics import calinski_harabasz_score ch_score calinski_harabasz_score(X_sample, z_hist[-1]) print(fCalinski-Harabasz Score: {ch_score:.0f}) # K10 时典型值500~800K5 时可能达 1200簇少簇间距离更大5.3 系统性对比K 从 3 到 15 的指标曲线真正体现最优化思维的是做超参扫描。我们画出不同 K 下的指标K_range range(3, 16) sil_scores [] ch_scores [] for K_test in K_range: print(f\nTesting K{K_test}...) mu_hist_t, z_hist_t, _ kmeans_manual( X_sample, KK_test, max_iters30, tol1e-3, random_state42 ) z_final z_hist_t[-1] sil silhouette_score(X_sample, z_final) ch calinski_harabasz_score(X_sample, z_final) sil_scores.append(sil) ch_scores.append(ch) print(f Silhouette: {sil:.3f}, CH: {ch:.0f}) # 绘图 plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(K_range, sil_scores, bo-) plt.xlabel(Number of Clusters (K)) plt.ylabel(Silhouette Score) plt.title(Silhouette Score vs K) plt.grid(True) plt.subplot(1, 2, 2) plt.plot(K_range, ch_scores, ro-) plt.xlabel(Number of Clusters (K)) plt.ylabel(Calinski-Harabasz Score) plt.title(CH Score vs K) plt.grid(True) plt.tight_layout() plt.show()你会看到轮廓系数在 K5~7 达到峰值约 0.18之后缓慢下降CH 指数则随 K 增加持续上升因为分更多簇簇内更紧凑。这揭示了最优化的本质矛盾没有唯一“最优 K”只有根据目标权衡的选择。若你关心簇内一致性选 K6若你需精细区分如“手写体变体分析”可选 K12接受更低的轮廓系数。我带实验课时常让学生交两份报告一份用 K10满足实验要求一份用 K6指标最优。后者往往能画出更清晰的质心图像——因为 6 个簇天然对应“直线型1,7、圆弧型0,6,8,9、三角型3,4,5”等粗粒度结构。这比死磕 K10 更体现对数据的理解。希望帮到你。本文还有配套的精品资源点击获取 K-Means在MNIST上的原理、实现与避坑指南
简介本资源是深圳大学计算机软件专业《最优化方法》课程配套实验材料面向机器学习初学者与高校算法实践者聚焦无监督学习核心任务——利用K-Means聚类实现MNIST手写数字图像的自动分组与结构发现。资源包共2个文件1个可运行Python脚本1份详细实验文档总大小873KB脚本完整实现数据加载、归一化预处理、肘部法则选K、聚类执行与结果可视化文档则系统梳理算法原理、关键参数设计依据、scikit-learn调用细节及聚类效果评估方法。已有3177人学习下载内容紧扣教学实验场景代码模块清晰、注释详尽附带可直接复现的绘图与分簇分析逻辑特别适合理解K-Means在真实图像数据上的行为边界与应用局限为后续监督分类打下扎实的无监督建模基础。1. K-Means 聚类真能分清手写数字——别急着调sklearn.cluster.KMeans先搞懂它在 MNIST 上为什么“看起来准、实际懵”你跑通了KMeans(n_clusters10).fit(X_train)聚类结果轮廓系数 0.12每个簇中心画出来像模像样甚至用匈牙利算法匹配标签后准确率冲到 53%——恭喜你成功复现了深大计软《最优化方法》实验1的经典翻车现场。这不是代码错了而是你还没意识到K-Means 在 MNIST 上不是分类器它是个“无监督的形状归组器”。它不认“3 是三”只认“这一坨像素和那一坨像素离得近”。当数字“1”和“7”在像素空间里挤在一起比如都写得细长倾斜K-Means 会毫不犹豫把它们划进同一个簇而“4”和“9”如果都带圆弧竖线也可能被强行拆散。这不是算法缺陷是欧氏距离 球形簇假设 无标签先验的必然结果。本实验的核心价值从来不是“用聚类替代分类”而是亲手推导目标函数、观察质心迭代如何卡在局部极小、验证维度灾难对距离失效的影响——这才是最优化方法课要锤你的点。适合人群刚学完梯度下降、拉格朗日乘子想把抽象公式落到像素矩阵上的计软本科生或想搞懂“为什么无监督评估指标总比有监督低一截”的算法初学者。别指望靠它上 Kaggle但如果你能手动实现 E-step/M-step 并画出每轮质心漂移轨迹你对“迭代优化”的理解就落地了。2. 从零推导 K-Means 目标函数为什么最小化平方误差等价于最大似然K-Means 表面是“找 10 个中心点让所有样本离最近中心最近”但它的数学根基藏在概率模型里。我们不直接讲 EM 算法而是从最优化视角拆解它本质是在求解一个带隐变量的非凸优化问题。先明确符号设训练集 $X {x^{(1)}, x^{(2)}, ..., x^{(N)}}$其中 $x^{(i)} \in \mathbb{R}^D$MNIST 是 784 维簇中心为 $\mu_1, \mu_2, ..., \mu_K \in \mathbb{R}^D$隐变量 $z^{(i)} \in {1,2,...,K}$ 表示第 $i$ 个样本所属簇。K-Means 的目标函数是$$ J(\mu, z) \sum_{i1}^{N} \sum_{k1}^{K} \mathbb{1}(z^{(i)} k) \cdot |x^{(i)} - \mu_k|^2 $$这个式子直白说就是对每个样本只算它到自己所属簇中心的距离平方再全部加起来。注意两点第一$z^{(i)}$ 是离散变量无法求导第二$J$ 关于 $\mu_k$ 是凸的但关于 $z^{(i)}$ 完全非凸。所以标准解法是坐标下降Coordinate Descent固定 $z$ 更新 $\mu$再固定 $\mu$ 更新 $z$交替进行直到收敛。这正是 K-Means 的两步迭代逻辑。2.1 E-step给定质心分配样本到最近簇硬分配这步不需要概率纯几何。对每个样本 $x^{(i)}$计算它到所有 $K$ 个质心的欧氏距离平方选最小的那个# 手动实现 E-step不用 sklearn import numpy as np def e_step(X, mu): X: (N, D) 样本矩阵 mu: (K, D) 质心矩阵 返回: (N,) 每个样本的簇索引数组 N, D X.shape K, _ mu.shape # 向量化计算所有距离平方避免 for 循环 # 利用 (a-b)^2 a^2 - 2ab b^2 展开 X_sq np.sum(X**2, axis1, keepdimsTrue) # (N, 1) mu_sq np.sum(mu**2, axis1, keepdimsTrue) # (K, 1) cross 2 * X mu.T # (N, K) # dist_sq[i, k] ||x_i - mu_k||^2 dist_sq X_sq - cross mu_sq.T # (N, K) return np.argmin(dist_sq, axis1) # (N,) # 验证用 MNIST 前 1000 张图测试 from sklearn.datasets import fetch_openml mnist fetch_openml(mnist_784, version1, as_frameFalse, parserauto) X, y mnist.data.astype(float32), mnist.target.astype(int) X_sample X[:1000] / 255.0 # 归一化到 [0,1] y_sample y[:1000] # 随机初始化 10 个质心均匀分布 np.random.seed(42) mu_init np.random.rand(10, 784) z_pred e_step(X_sample, mu_init) print(f初始分配簇大小分布 {np.bincount(z_pred, minlength10)}) # 输出类似[112 98 105 92 101 95 103 97 99 98] —— 基本均匀参数说明e_step中的dist_sq计算用了向量化技巧避免for i in range(N): for k in range(K):的 O(NKD) 复杂度。关键在X_sq - cross mu_sq.T这一行X_sq是每个样本自身模长平方mu_sq.T是每个质心模长平方转置cross是内积矩阵。这样一次广播运算就得到全部距离平方速度提升 100 倍以上。新手常在这里写嵌套循环跑 MNIST 直接卡死。2.2 M-step给定分配更新质心为簇内均值解析解这步是目标函数 $J$ 关于 $\mu_k$ 的最小化。对固定 $k$只看属于簇 $k$ 的样本集合 $C_k {i | z^{(i)} k}$则$$ \frac{\partial J}{\partial \mu_k} -2 \sum_{i \in C_k} (x^{(i)} - \mu_k) 0 \quad \Rightarrow \quad \mu_k \frac{1}{|C_k|} \sum_{i \in C_k} x^{(i)} $$看到没质心更新就是取簇内样本均值——这是唯一解析解不用梯度下降。这也是 K-Means 收敛快的原因M-step 是闭式解。def m_step(X, z, K): X: (N, D) 样本 z: (N,) 簇索引 K: 簇数 返回: (K, D) 新质心矩阵 N, D X.shape mu_new np.zeros((K, D)) for k in range(K): # 找出属于簇 k 的所有样本 mask (z k) if np.sum(mask) 0: # 防止空簇 mu_new[k] np.mean(X[mask], axis0) else: # 空簇处理重采样一个随机样本作为新质心常见策略 mu_new[k] X[np.random.randint(0, N)] return mu_new # 测试 M-step z_init e_step(X_sample, mu_init) mu_updated m_step(X_sample, z_init, K10) print(f更新后质心形状: {mu_updated.shape}) # (10, 784)关键细节m_step中必须处理空簇empty cluster。如果某轮 E-step 后某个簇没分到任何样本np.mean会报错或返回全零向量导致后续迭代崩溃。这里采用“重采样随机样本”策略也有用“距其他质心最远点”或“最大方差方向切分”的但对 MNIST 来说随机重采样最稳定。这是实操中第一个血泪经验永远检查np.bincount(z)是否有 0。3. 手动实现完整 K-Means 迭代监控损失下降、可视化质心漂移、验证收敛性现在把 E-step 和 M-step 串起来加上收敛判断和日志。重点不是“跑出结果”而是看见优化过程本身——这才是最优化方法课的灵魂。我们监控三个量目标函数值 $J$、质心移动距离、簇分配变化率。3.1 主循环带收敛判断与日志的迭代框架def kmeans_manual(X, K, max_iters100, tol1e-4, random_state42): 手动实现 K-Means 完整流程 X: (N, D) 归一化后的样本MNIST 已除以 255 K: 簇数对 MNIST 固定为 10 max_iters: 最大迭代次数 tol: 损失变化容忍阈值 返回: mu_history (list of (K,D)), z_history (list of (N,)), losses (list of float) np.random.seed(random_state) N, D X.shape # 初始化质心用 K-means 策略比随机好得多 mu np.zeros((K, D)) # 第一个质心随机选 mu[0] X[np.random.randint(0, N)] # 后续质心按距离平方概率选 for k in range(1, K): # 计算所有点到已选质心的最小距离平方 dist_sq_to_mu np.full(N, np.inf) for i in range(k): dist_sq np.sum((X - mu[i])**2, axis1) dist_sq_to_mu np.minimum(dist_sq_to_mu, dist_sq) # 概率正比于 dist_sq_to_mu probs dist_sq_to_mu / np.sum(dist_sq_to_mu) mu[k] X[np.random.choice(N, pprobs)] mu_history [mu.copy()] z_history [] losses [] for it in range(max_iters): # E-step z e_step(X, mu) z_history.append(z) # 计算当前损失 J loss 0.0 for k in range(K): mask (z k) if np.sum(mask) 0: loss np.sum((X[mask] - mu[k])**2) losses.append(loss) # M-step mu_new m_step(X, z, K) # 检查收敛质心移动距离 tol mu_shift np.max(np.sqrt(np.sum((mu_new - mu)**2, axis1))) mu_history.append(mu_new.copy()) mu mu_new # 打印进度 if it % 10 0 or it max_iters-1: print(fIter {it:3d} | Loss: {loss:.2e} | Max shift: {mu_shift:.4f}) if mu_shift tol: print(fConverged at iteration {it}) break return mu_history, z_history, losses # 在 MNIST 子集上运行 mu_hist, z_hist, losses kmeans_manual(X_sample, K10, max_iters50, tol1e-3)为什么用 K-means 初始化随机初始化可能导致质心全挤在数字“1”的区域其他数字如“8”、“6”永远分不到簇。K-means 通过距离加权采样强制质心分散通常减少 30% 迭代次数。dist_sq_to_mu的计算是核心对每个未选点算它到所有已选质心的最小距离平方再以此为权重抽样。这是工业级实现的标配不是炫技。3.2 可视化质心演化从噪声到“数字雏形”的 30 轮旅程K-Means 的质心不是静态图片是动态优化的产物。我们把每轮质心 reshape 成 28x28 并画出来import matplotlib.pyplot as plt def plot_centroids_evolution(mu_history, n_cols5, figsize(12, 8)): 绘制质心随迭代的变化 mu_history: list of (K, 784) 质心矩阵 n_iters len(mu_history) n_rows (n_iters n_cols - 1) // n_cols fig, axes plt.subplots(n_rows, n_cols, figsizefigsize) axes axes.flatten() if n_iters 1 else [axes] for i, mu in enumerate(mu_history): if i len(axes): break # 取第一个质心索引 0为例展示其变化 # 实际可循环画所有 10 个但这里简化 img mu[0].reshape(28, 28) axes[i].imshow(img, cmapgray) axes[i].set_title(fIter {i}) axes[i].axis(off) # 隐藏多余子图 for j in range(i1, len(axes)): axes[j].remove() plt.tight_layout() plt.show() # 运行绘图 plot_centroids_evolution(mu_hist)你会看到第 0 轮是随机噪声块第 5 轮开始出现模糊的竖线/横线第 15 轮能辨认出“1”、“0”的轮廓第 30 轮基本稳定但和真实数字仍有差距——因为 K-Means 学的是像素均值不是笔画结构。这就是最优化的真相它在约束下逼近局部最优而非生成完美图像。3.3 监控损失曲线与收敛诊断为什么有时损失不降反升画出losses曲线plt.figure(figsize(10, 4)) plt.plot(losses, b-o, markersize3) plt.xlabel(Iteration) plt.ylabel(Objective J (Sum of Squared Errors)) plt.title(K-Means Loss Curve on MNIST Subset) plt.grid(True) plt.show()正常情况是单调下降。但如果出现震荡或上升一定是代码 bug。常见原因e_step中距离计算错误比如忘了平方或用了曼哈顿距离m_step中空簇未处理导致mu_new[k]为全零下轮e_step计算距离时||x_i - 0||^2极大初始化质心超出数据范围如用了np.random.randn未归一化收敛性验证口诀K-Means 保证每轮 $J$ 不增但不保证全局最优。若losses曲线非单调立刻检查e_step的dist_sq公式——90% 的问题出在这里。4. 避坑MNIST 上 K-Means 的 4 个经典翻车现场与急救方案K-Means 在 MNIST 上不是不能跑而是处处是坑。这些坑不来自算法本身而来自数据特性与实现细节的碰撞。以下是我在深大计软助教三年批改上百份实验报告总结的最高频、最隐蔽、最致命的四个问题每条都附带现象、根因和一行修复代码。4.1 现象迭代 50 轮后某个簇的样本数为 0后续m_step报RuntimeWarning: Mean of empty slice原因E-step 分配时若某质心离所有样本都远可能没分到任何点M-step 中np.mean(X[mask])对空数组返回nan污染后续计算。解决在m_step中强制检查空簇并用安全策略重置质心。不要用np.nanmean它返回nan不解决问题# 错误示范引发连锁 nan # mu_new[k] np.mean(X[mask], axis0) # mask 全 False → nan # 正确修复在 m_step 函数内 if np.sum(mask) 0: # 方案1重采样随机样本推荐简单鲁棒 mu_new[k] X[np.random.randint(0, N)] # 方案2用所有样本的均值更稳定但可能偏移 # mu_new[k] np.mean(X, axis0)4.2 现象质心图像全是灰色块像素值集中在 0.4~0.6完全看不出数字形状原因MNIST 像素是uint80~255但你没归一化KMeans对量纲敏感未归一化的像素值0~255会让距离计算被高亮区域主导质心被拉向平均灰度。解决必须在输入前除以 255缩放到 [0,1]# 错误直接喂原始数据 # X_raw mnist.data # shape (70000, 784), dtype uint8 # 正确归一化是铁律 X_normalized mnist.data.astype(float32) / 255.0 # 验证print(X_normalized.min(), X_normalized.max()) → 应输出 0.0 1.04.3 现象e_step运行极慢10 秒CPU 占用 100%dist_sq计算卡死原因写了双重 for 循环计算距离复杂度 O(NKD)MNIST N70000, K10, D784 → 54.88 亿次运算。解决用向量化公式X_sq - 2*Xmu.T mu_sq.T一行替代循环# 错误慢如蜗牛 # dist_sq np.zeros((N, K)) # for i in range(N): # for k in range(K): # dist_sq[i,k] np.sum((X[i] - mu[k])**2) # 正确毫秒级 X_sq np.sum(X**2, axis1, keepdimsTrue) # (N, 1) mu_sq np.sum(mu**2, axis1, keepdimsTrue) # (K, 1) dist_sq X_sq - 2 * X mu.T mu_sq.T # (N, K) 广播4.4 现象用匈牙利算法匹配簇标签后准确率只有 10%~20%远低于预期原因你匹配的是“簇 ID”和“数字标签”但 K-Means 的簇 ID 是任意的0~9而数字标签 0~9 有语义。若簇 0 匹配到数字 5簇 1 匹配到数字 0直接按 ID 对齐当然错。必须用二分图匹配Hungarian Algorithm找最优映射。解决用scipy.optimize.linear_sum_assignment计算最佳匹配from scipy.optimize import linear_sum_assignment from sklearn.metrics import confusion_matrix def calculate_matching_accuracy(z_pred, y_true, K10): z_pred: (N,) 预测簇标签 y_true: (N,) 真实数字标签0~9 返回: 匹配后的准确率 # 构建混淆矩阵行簇ID列真实标签 cm confusion_matrix(y_true, z_pred, labelsrange(K)) # Hungarian 算法找最大匹配 row_ind, col_ind linear_sum_assignment(-cm) # 最大化故取负 # 计算匹配总正确数 total_correct cm[row_ind, col_ind].sum() return total_correct / len(y_true) acc calculate_matching_accuracy(z_hist[-1], y_sample) print(fMatching Accuracy: {acc:.3f}) # 正常应 0.5~0.6提示匹配准确率 50%~60% 是 MNIST 上 K-Means 的合理上限。别追求 95%——那属于有监督学习的领域。K-Means 的价值在于无标签下的结构发现不是替代分类器。5. 进阶验证用轮廓系数、Calinski-Harabasz 指标量化聚类质量并对比 K5 vs K10K-Means 的“好坏”不能只看准确率毕竟它没标签。我们要用无监督评估指标来回答当前 K10 是最优的吗质心真的分开了吗这里用两个黄金指标轮廓系数Silhouette Score和 Calinski-Harabasz 指数。5.1 轮廓系数衡量“簇内紧密度 vs 簇间分离度”对每个样本 $i$定义$a(i)$$i$ 到同簇其他点的平均距离簇内不相似度$b(i)$$i$ 到最近其他簇所有点的平均距离簇间不相似度轮廓值 $s(i) \frac{b(i)-a(i)}{\max(a(i),b(i))} \in [-1,1]$$s(i)$ 接近 1 表示样本聚得好接近 -1 表示分错簇。平均轮廓系数越高聚类越优。from sklearn.metrics import silhouette_score # 计算最终聚类的轮廓系数 silhouette_avg silhouette_score(X_sample, z_hist[-1]) print(fAverage Silhouette Score: {silhouette_avg:.3f}) # 典型值K10 时约 0.10~0.15K5 时可能升到 0.18因为簇更大更易分离为什么 MNIST 的轮廓系数这么低因为数字“4”和“9”、“3”和“8”在像素空间本就相似都含圆弧直线K-Means 无法用线性边界分开它们。低轮廓系数恰恰证明了数据本身的可分性有限不是算法不行。5.2 Calinski-Harabasz 指数簇间离散度 / 簇内离散度公式$CH \frac{Tr(B_k)}{Tr(W_k)} \times \frac{N-K}{K-1}$其中 $B_k$ 是簇间散度矩阵$W_k$ 是簇内散度矩阵。值越大越好。from sklearn.metrics import calinski_harabasz_score ch_score calinski_harabasz_score(X_sample, z_hist[-1]) print(fCalinski-Harabasz Score: {ch_score:.0f}) # K10 时典型值500~800K5 时可能达 1200簇少簇间距离更大5.3 系统性对比K 从 3 到 15 的指标曲线真正体现最优化思维的是做超参扫描。我们画出不同 K 下的指标K_range range(3, 16) sil_scores [] ch_scores [] for K_test in K_range: print(f\nTesting K{K_test}...) mu_hist_t, z_hist_t, _ kmeans_manual( X_sample, KK_test, max_iters30, tol1e-3, random_state42 ) z_final z_hist_t[-1] sil silhouette_score(X_sample, z_final) ch calinski_harabasz_score(X_sample, z_final) sil_scores.append(sil) ch_scores.append(ch) print(f Silhouette: {sil:.3f}, CH: {ch:.0f}) # 绘图 plt.figure(figsize(12, 4)) plt.subplot(1, 2, 1) plt.plot(K_range, sil_scores, bo-) plt.xlabel(Number of Clusters (K)) plt.ylabel(Silhouette Score) plt.title(Silhouette Score vs K) plt.grid(True) plt.subplot(1, 2, 2) plt.plot(K_range, ch_scores, ro-) plt.xlabel(Number of Clusters (K)) plt.ylabel(Calinski-Harabasz Score) plt.title(CH Score vs K) plt.grid(True) plt.tight_layout() plt.show()你会看到轮廓系数在 K5~7 达到峰值约 0.18之后缓慢下降CH 指数则随 K 增加持续上升因为分更多簇簇内更紧凑。这揭示了最优化的本质矛盾没有唯一“最优 K”只有根据目标权衡的选择。若你关心簇内一致性选 K6若你需精细区分如“手写体变体分析”可选 K12接受更低的轮廓系数。我带实验课时常让学生交两份报告一份用 K10满足实验要求一份用 K6指标最优。后者往往能画出更清晰的质心图像——因为 6 个簇天然对应“直线型1,7、圆弧型0,6,8,9、三角型3,4,5”等粗粒度结构。这比死磕 K10 更体现对数据的理解。希望帮到你。本文还有配套的精品资源点击获取