简介：本资源是深圳大学计算机软件专业《最优化方法》课程配套实验材料，面向机器学习初学者与高校算法实践者，聚焦无监督学习核心任务——利用K-Means聚类实现MNIST手写数字图像的自动分组与结构发现。资源包共2个文件（1个可运行Python脚本+1份详细实验文档），总大小873KB，脚本完整实现数据加载、归一化预处理、肘部法则选K、聚类执行与结果可视化；文档则系统梳理算法原理、关键参数设计依据、scikit-learn调用细节及聚类效果评估方法。已有3177人学习下载，内容紧扣教学实验场景，代码模块清晰、注释详尽，附带可直接复现的绘图与分簇分析逻辑，特别适合理解K-Means在真实图像数据上的行为边界与应用局限，为后续监督分类打下扎实的无监督建模基础。
1. K-Means 聚类真能分清手写数字？——别急着调sklearn.cluster.KMeans，先搞懂它在 MNIST 上为什么“看起来准、实际懵”
你跑通了KMeans(n_clusters=10).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_{i=1}^{N} \sum_{k=1}^{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, axis=1, keepdims=True) # (N, 1) mu_sq = np.sum(mu**2, axis=1, keepdims=True) # (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, axis=1) # (N,) # 验证：用 MNIST 前 1000 张图测试 from sklearn.datasets import fetch_openml mnist = fetch_openml('mnist_784', version=1, as_frame=False, parser='auto') 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, minlength=10)}") # 输出类似：[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], axis=0) 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, K=10) 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_iters=100, tol=1e-4, random_state=42): """ 手动实现 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, axis=1) 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, p=probs)] 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, axis=1))) mu_history.append(mu_new.copy()) mu = mu_new # 打印进度 if it % 10 == 0 or it == max_iters-1: print(f"Iter {it:3d} | Loss: {loss:.2e} | Max shift: {mu_shift:.4f}") if mu_shift < tol: print(f"Converged at iteration {it}") break return mu_history, z_history, losses # 在 MNIST 子集上运行 mu_hist, z_hist, losses = kmeans_manual(X_sample, K=10, max_iters=50, tol=1e-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_cols=5, 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, figsize=figsize) 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, cmap='gray') axes[i].set_title(f'Iter {i}') axes[i].axis('off') # 隐藏多余子图 for j in range(i+1, 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', markersize=3) 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], axis=0) # 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, axis=0)4.2 现象：质心图像全是灰色块（像素值集中在 0.4~0.6），完全看不出数字形状
- 原因：MNIST 像素是
uint8（0~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 N=70000, K=10, D=784 → 54.88 亿次运算。
- 解决：用向量化公式
X_sq - 2*X@mu.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, axis=1, keepdims=True) # (N, 1) mu_sq = np.sum(mu**2, axis=1, keepdims=True) # (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, K=10): """ z_pred: (N,) 预测簇标签 y_true: (N,) 真实数字标签（0~9） 返回: 匹配后的准确率 """ # 构建混淆矩阵：行=簇ID，列=真实标签 cm = confusion_matrix(y_true, z_pred, labels=range(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(f"Matching Accuracy: {acc:.3f}") # 正常应 0.5~0.6提示：匹配准确率 50%~60% 是 MNIST 上 K-Means 的合理上限。别追求 95%——那属于有监督学习的领域。K-Means 的价值在于无标签下的结构发现，不是替代分类器。
5. 进阶验证：用轮廓系数、Calinski-Harabasz 指标量化聚类质量，并对比 K=5 vs K=10
K-Means 的“好坏”不能只看准确率（毕竟它没标签）。我们要用无监督评估指标来回答：当前 K=10 是最优的吗？质心真的分开了吗？这里用两个黄金指标：轮廓系数（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(f"Average Silhouette Score: {silhouette_avg:.3f}") # 典型值：K=10 时约 0.10~0.15；K=5 时可能升到 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(f"Calinski-Harabasz Score: {ch_score:.0f}") # K=10 时典型值：500~800；K=5 时可能达 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, K=K_test, max_iters=30, tol=1e-3, random_state=42 ) 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()你会看到：轮廓系数在 K=5~7 达到峰值（约 0.18），之后缓慢下降；CH 指数则随 K 增加持续上升（因为分更多簇，簇内更紧凑）。这揭示了最优化的本质矛盾：没有唯一“最优 K”，只有根据目标权衡的选择。若你关心簇内一致性，选 K=6；若你需精细区分（如“手写体变体分析”），可选 K=12，接受更低的轮廓系数。
我带实验课时，常让学生交两份报告：一份用 K=10（满足实验要求），一份用 K=6（指标最优）。后者往往能画出更清晰的质心图像——因为 6 个簇天然对应“直线型（1,7）、圆弧型（0,6,8,9）、三角型（3,4,5）”等粗粒度结构。这比死磕 K=10 更体现对数据的理解。希望帮到你。
本文还有配套的精品资源，点击获取
