模式识别课程实验包：K-means、GMM与感知机实战MNIST

简介这份资源是哈尔滨工业大学模式识别课程的实验配套代码包面向正在学习机器学习与模式识别的高校学生及自学者帮助其通过动手实践掌握均值聚类、高斯混合模型与感知机等经典算法。包内共7个文件以6个Python脚本和1份Markdown说明为主压缩包约12KB体量轻便便于快速运行与阅读。内容围绕四个实验展开K-means用于手写字体特征的初步聚类GMM对复杂笔画分布进行概率建模感知机结合LMSE实现线性分类并最终在MNIST数据集上完成十分类识别挑战。代码结构清晰覆盖从无监督到有监督的完整流程读者可借此理解算法原理、参数估计与分类边界优化并直接复用脚本进行调试与扩展。目前已有167人学习适合作为课程作业参考或模式识别入门练手项目。1. 从一份课程实验包说起K-means、GMM、感知机怎么串起手写数字识别如果你正在补模式识别的课设或者想找一个能把无监督聚类、概率生成模型和线性分类器串起来练手的小项目这份课程实验包值得拆开看看。它把四个实验放在同一个目录里实验1用K-means做聚类实验2用GMM做概率建模实验3用感知机加LMSE做线性分类实验4直接上MNIST做多分类考试。每个实验都有独立的main.py实验2还额外留了main_back2.py和back.py两个备份版本说明作者在调参和改结构时踩过反复回滚的坑。这份资源最适合两类人一类是刚学完模式识别理论、需要把公式落成代码的学生另一类是已经工作、想快速回顾经典算法实现细节的工程师。它不依赖深度学习框架核心逻辑基本靠NumPy手写所以你能清楚看到质心怎么更新、协方差矩阵怎么估计、权重怎么迭代。代价是代码风格偏课程作业缺少工程封装直接拿去做生产级任务需要自己补数据管道和评估模块。下面按“先跑通、再调参、后避坑”的顺序拆。2. 实验1与实验3K-means聚类和感知机LMSE的代码骨架2.1 K-means的初始化陷阱与迭代终止条件实验1的main.py通常包含三个核心函数距离计算、质心更新、类别分配。课程作业里最常见的写法是随机从样本中选K个点作为初始质心然后反复执行“分配—更新”直到质心不再移动或达到最大迭代次数。这里第一个翻车点就是初始化如果随机选到的质心全部落在同一个密集区域K-means会收敛到一个局部最优簇的划分严重偏斜。我一般会建议把初始化改成K-means或者至少跑多次随机初始化取组内平方误差和最小的那次。下面是一个可抄的骨架保留了课程作业的NumPy风格同时补上了多次重启import numpy as np def kmeans(X, K, max_iter100, n_init10, tol1e-4): best_labels, best_centers, best_inertia None, None, np.inf for _ in range(n_init): # K-means 初始化第一个质心随机选后续按距离平方概率选 centers [X[np.random.randint(len(X))]] for _ in range(1, K): dist_sq np.array([min([np.sum((x - c) ** 2) for c in centers]) for x in X]) prob dist_sq / dist_sq.sum() centers.append(X[np.random.choice(len(X), pprob)]) centers np.array(centers) for _ in range(max_iter): # 分配步骤计算每个样本到各质心的距离取最近 distances np.linalg.norm(X[:, None] - centers[None, :], axis2) labels np.argmin(distances, axis1) # 更新步骤按类别求均值 new_centers np.array([X[labels k].mean(axis0) if np.any(labels k) else centers[k] for k in range(K)]) if np.linalg.norm(new_centers - centers) tol: break centers new_centers inertia sum(np.min(distances, axis1) ** 2) if inertia best_inertia: best_inertia, best_labels, best_centers inertia, labels, centers return best_labels, best_centers, best_inertia逻辑说明n_init10表示跑10次不同初始化取组内平方误差和最小的结果。tol1e-4是质心移动的容忍度小于这个值就认为收敛。参数上K值需要根据手写数字的类别数来定MNIST是10类但K-means是无监督的你可以先设K10看聚类结果和真实标签的混淆矩阵再调整。注意np.any(labels k)这个判断如果某个簇没有分配到任何样本质心保持不变否则会报空数组求均值的错。2.2 感知机LMSE的权重更新与学习率选择实验3的main.py做的是感知机加最小均方误差LMSE训练。感知机的经典更新规则是对误分类样本权重加上学习率 * 样本 * 标签。LMSE则把目标改成最小化均方误差权重更新变成学习率 * (标签 - 预测) * 样本。课程作业里通常用批量梯度下降或随机梯度下降手写数字二分类时一般把某个数字当作正类其余当作负类。import numpy as np def train_lmse(X, y, lr0.01, epochs200): # X: (N, D), y: (N,) 取值为 1 / -1 N, D X.shape w np.zeros(D) b 0.0 losses [] for epoch in range(epochs): # 批量梯度下降 y_pred X w b error y - y_pred grad_w -2 * X.T error / N grad_b -2 * error.sum() / N w - lr * grad_w b - lr * grad_b loss np.mean(error ** 2) losses.append(loss) return w, b, losses逻辑说明lr是学习率课程作业里常见取值是0.01或0.001太大容易震荡太小收敛慢。epochs是遍历整个数据集的次数200轮通常够用但如果你把学习率调得很小可能需要500轮以上。y必须转成1/-1如果原始标签是0/1要先用y 2 * y - 1映射。注意LMSE对异常值敏感如果某个样本的特征尺度特别大梯度会被它主导所以跑之前最好做标准化。2.3 两个实验的串联思路实验1和实验3可以串起来用先用K-means对训练样本做聚类把每个样本到各质心的距离作为新特征再喂给感知机做分类。这种“无监督特征提取有监督分类”的流程在早期模式识别里很常见能帮你理解为什么聚类结果可以作为分类的预处理。代价是K-means的簇标签没有语义你需要自己建立簇编号到真实类别的映射通常用匈牙利算法或者简单的混淆矩阵匹配。3. 实验2GMM参数估计与EM算法的落地细节3.1 协方差矩阵的正则化与数值稳定实验2的main.py和main_back2.py是GMM的核心。GMM用多个高斯分布的加权和来建模数据分布参数包括每个分量的均值、协方差矩阵和权重。EM算法分两步E步计算每个样本属于每个分量的后验概率M步用这些概率加权更新均值、协方差和权重。课程作业里最容易翻车的地方是协方差矩阵奇异——当某个分量只分配到极少样本或者样本在某个维度上方差接近零协方差矩阵不可逆行列式接近零对数似然直接变成-inf或NaN。常见做法是在协方差矩阵的对角线上加一个小的正则项比如1e-6 * np.eye(D)。另外初始化权重时避免让某个分量一开始就拿到零权重可以统一设为1/K。下面是一个简化但可运行的GMM实现import numpy as np def gmm_em(X, K, max_iter100, tol1e-4, reg1e-6): N, D X.shape # 初始化均值从样本中随机选协方差设为单位阵权重均匀 means X[np.random.choice(N, K, replaceFalse)] covs np.array([np.eye(D) for _ in range(K)]) weights np.ones(K) / K log_likelihoods [] for iteration in range(max_iter): # E步计算后验概率 gamma gamma np.zeros((N, K)) for k in range(K): diff X - means[k] cov_reg covs[k] reg * np.eye(D) inv_cov np.linalg.inv(cov_reg) det_cov np.linalg.det(cov_reg) coef 1.0 / ((2 * np.pi) ** (D / 2) * np.sqrt(det_cov)) exponent -0.5 * np.sum(diff inv_cov * diff, axis1) gamma[:, k] weights[k] * coef * np.exp(exponent) gamma_sum gamma.sum(axis1, keepdimsTrue) gamma_sum[gamma_sum 0] 1e-10 # 防止除零 gamma / gamma_sum # M步更新均值、协方差、权重 Nk gamma.sum(axis0) for k in range(K): means[k] (gamma[:, k] X) / Nk[k] diff X - means[k] covs[k] (gamma[:, k] * diff.T) diff / Nk[k] reg * np.eye(D) weights[k] Nk[k] / N # 计算对数似然用于判断收敛 log_likelihood np.sum(np.log(gamma_sum)) log_likelihoods.append(log_likelihood) if iteration 0 and abs(log_likelihood - log_likelihoods[-2]) tol: break return means, covs, weights, log_likelihoods逻辑说明reg1e-6是协方差正则项防止矩阵奇异。gamma是后验概率矩阵形状(N, K)。E步里用对数域计算更稳定但课程作业里直接算概率再归一化也能跑。M步更新协方差时用了(gamma[:, k] * diff.T) diff这是加权协方差的向量化写法。log_likelihood用np.log(gamma_sum)求和注意gamma_sum在归一化前是未归一化的概率取对数后求和就是对数似然。如果发现对数似然震荡不收敛通常是K设得太大或者数据维度太高可以先用PCA降到2维可视化看看。3.2 分量数K的选择与BIC准则GMM的K值不能像K-means那样靠肘部法随便看因为GMM的似然函数随K增大单调递增K越多似然越高但会过拟合。课程作业里一般直接指定K10对应MNIST的10类但如果你想选得更合理可以用BIC贝叶斯信息准则BIC -2 * log_likelihood p * log(N)其中p是参数个数GMM的参数个数是K * (D D*(D1)/2 1) - 1。BIC越小越好。我一般会跑K从2到15画一条BIC曲线选拐点附近的K。3.3 备份文件main_back2.py和back.py的启示实验2目录里多了两个备份文件这其实是个信号GMM的调参过程很容易让人反复回滚。常见的情况是改了初始化方式后对数似然反而下降或者加了正则项后收敛变慢。我的习惯是每次改参数前先复制一份main.py为main_back.py在备份里保留上一版能跑通的参数新版本跑不通就回退。另外main_back2.py可能对应另一种初始化策略比如用K-means的聚类结果来初始化GMM的均值这样EM算法收敛更快也不容易陷入糟糕的局部最优。4. 实验4MNIST多分类的评估与常见翻车点4.1 数据加载与标签对齐实验4的main.py做MNIST分类考试通常会调用前面实验的模块或者直接用一个多层感知机、SVM做对比。MNIST的原始文件是IDX格式不是图片文件需要自己解析。课程作业里常见做法是用sklearn.datasets.fetch_openml(mnist_784)或者手动读二进制文件。手动读的话注意大端序前4个字节是魔数接着4个字节是图像数量再4个字节是行数再4个字节是列数然后才是像素数据。import numpy as np def load_mnist_images(filename): with open(filename, rb) as f: magic int.from_bytes(f.read(4), big) num_images int.from_bytes(f.read(4), big) rows int.from_bytes(f.read(4), big) cols int.from_bytes(f.read(4), big) images np.frombuffer(f.read(), dtypenp.uint8) images images.reshape(num_images, rows * cols) return images def load_mnist_labels(filename): with open(filename, rb) as f: magic int.from_bytes(f.read(4), big) num_labels int.from_bytes(f.read(4), big) labels np.frombuffer(f.read(), dtypenp.uint8) return labels逻辑说明int.from_bytes(..., big)按大端序读整数这是IDX格式的规定。np.frombuffer直接把二进制缓冲区转成数组比逐字节读快很多。标签文件只有魔数和数量两个头字段后面全是标签。注意图像数据要除以255做归一化否则梯度下降会震荡。4.2 评估指标准确率之外的混淆矩阵课程作业里通常只看准确率但MNIST的10类分类里某些数字容易混淆比如4和9、3和5、7和1。只看准确率会掩盖这些问题。我一般会额外算混淆矩阵看看哪两类之间误判最多。如果发现4和9混淆严重可以针对这两类单独训练一个二分类器做后处理。另外训练集和测试集的划分要固定随机种子否则每次跑出来的准确率波动几个百分点没法对比不同算法的优劣。4.3 多分类策略一对多与一对一感知机本身是二分类器做MNIST的10类分类需要扩展。常见做法是一对多训练10个二分类器每个把某个数字当作正类其余当作负类预测时取置信度最高的那个。另一种是一对一训练10*9/245个二分类器每个只区分两个数字预测时用投票法。一对多的训练样本不平衡正类1份负类9份一对一每个分类器的样本更均衡但数量多。课程作业里一般用一对多因为代码简单但如果你发现某个数字的召回率特别低可以换成一对一试试。5. 避坑与排查课程实验里最容易翻车的五件事5.1 现象K-means跑出来的簇全是空的或者某个簇只有一个样本原因初始化质心太集中或者K值设得比实际类别数大很多。解决改用K-means初始化或者先跑PCA降维再聚类。如果某个簇只有一个样本检查是不是有离群点把质心拉过去了可以在距离计算时用中位数代替均值或者先做异常值剔除。5.2 现象GMM的对数似然出现NaN程序直接崩原因协方差矩阵奇异行列式为0或负数取对数报错。解决在协方差矩阵上加正则项reg * np.eye(D)reg取1e-6到1e-4之间。另外检查是不是某个分量的权重变成了0导致Nk[k]为0除零产生NaN。可以在M步更新前判断Nk[k] 1e-10就跳过该分量的更新。5.3 现象感知机LMSE的损失不下降或者震荡原因学习率太大或者特征没有标准化。解决先把学习率降到0.001或0.0001试如果还是震荡检查特征尺度。MNIST像素值在0到255之间直接拿去训练梯度会爆炸必须除以255。另外如果标签没有转成1/-1LMSE的误差计算会完全错误。5.4 现象MNIST分类准确率只有10%左右相当于随机猜原因标签和图像没有对齐或者数据加载时字节序读反了。解决先打印前10个样本的标签和图像均值看看标签是不是0到9均匀分布图像均值是不是在0.1到0.9之间。如果标签全是0或者图像全是噪声检查int.from_bytes的字节序参数是不是big。另外如果用了sklearn的fetch_openml注意它返回的标签可能是字符串需要转成整数。5.5 现象实验2的备份文件越改越乱最后不知道哪个版本能跑原因没有版本管理习惯直接在主文件上改改坏了就复制一份加后缀。解决每次改参数前先提交一次git或者至少把能跑通的版本重命名为main_stable.py。备份文件不要用back.py、back2.py这种无意义命名改成main_kmeans_init.py、main_random_init.py一眼能看出区别。我一般会在文件头写一行注释记录当前参数和对应的对数似然值回滚时不用重新跑。6. 进阶技巧用K-means初始化GMM把收敛速度提上来GMM的EM算法对初始均值很敏感随机选初始均值经常收敛到差的局部最优。一个实用的技巧是先用K-means跑一遍把K-means的质心作为GMM的初始均值协方差统一设为单位阵权重设为均匀。这样EM算法通常能在20轮内收敛而且对数似然比随机初始化高。下面是对实验2的改造def init_gmm_with_kmeans(X, K): # 先用K-means得到质心和硬分配 labels, centers, _ kmeans(X, K, n_init5) means centers covs np.array([np.eye(X.shape[1]) for _ in range(K)]) weights np.array([np.sum(labels k) / len(X) for k in range(K)]) return means, covs, weights逻辑说明kmeans函数复用第2章的实现n_init5跑5次取最优。weights用每个簇的样本占比初始化比均匀权重更贴近数据分布。把返回的means, covs, weights直接传给gmm_em作为初始值EM算法的迭代次数通常能从100轮降到20轮左右。注意K-means的硬分配和GMM的软分配不同但作为初始化足够好。另一个技巧是监控对数似然的变化曲线。如果曲线在某个值附近震荡说明K设得太大或者正则项太小。我一般会把每次迭代的对数似然存下来画一条曲线如果看到明显的平台期就提前停止。还有如果发现某个分量的权重越来越小趋近于0说明这个分量是多余的可以把它去掉重新跑或者增大正则项让它和其他分量合并。从那以后我每次跑GMM之前都会先用K-means初始化并且把对数似然曲线画出来看一眼。这个习惯帮我省了很多调参时间也避免了对数似然突然崩掉的玄学问题。希望帮到你。本文还有配套的精品资源点击获取