简介这份代码练习包面向正在学习周志华《机器学习》的读者适合希望把教材公式落到代码层面的初学者与课程实验者用于巩固经典算法的实现思路、训练流程与评估方式。压缩包共16个文件以15个Python脚本为主另附1份README说明整体约20KB体量轻便便于按章节逐个运行和修改。内容围绕线性模型、决策树、集成学习、支持向量机以及偏差—方差分析等主题展开包含逻辑回归、留一法与交叉验证评估、LDA、Iterated Bagging、MultiBoosting、简化SVM与SVR对比等练习可帮助读者对照教材理解算法细节、观察参数与泛化表现。已有510人学习适合作为课程配套练习和自学复现的参考。1. 从一行for循环开始的机器学习练习第一次把自己手写的感知机训练过程打印成日志看着损失从 0.69 一路抖到 0.02那种感觉和调通一个线上接口完全不是一回事——调接口你知道成功就在眼前训练模型你只能靠数字判断它有没有在学。这个练习包做的就是这件事它把一本经典机器学习教材里每一章的算法用 NumPy 从零写一遍线性回归、对数几率回归、决策树、支持向量机、AdaBoost、k-means、PCA 一个不落每个文件都能单独跑出数值和图表。它适合两类人。一类是刚啃完理论、怀疑自己对公式的理解有偏差的跑一遍就知道梯度推得对不对另一类是工作几年、调包调出手感但心里没底的想看看那些fit()背后到底在算什么。这个包不负责教你调参它负责让你把拉格朗日乘子、信息增益、EM 的 E 步真正写成能跑的循环。2. 练习包的目录结构与技术选型动手之前先摸清骨架拿到一个算法练习包最忌讳的操作是直接python xxx.py跑一遍报错了再从头翻。正确顺序是先看目录——目录结构决定了这份代码是能复用的工程还是一次性脚本两者的阅读方式完全不同。2.1 从章节目录看代码组织一份可预测的文件地图这类配套代码的组织方式高度趋同基本是一章一个目录、一个算法一个文件。整理出来大致长这样ml-practice/ ├── ch03_linear_model/ │ ├── linear_regression.py # 最小二乘 梯度下降 │ ├── logistic_regression.py # 对数几率回归牛顿法 │ └── lda.py # 线性判别分析 ├── ch04_decision_tree/ │ └── tree.py # ID3/C4.5 的划分与递归建树 ├── ch05_neural_network/ │ └── bp.py # 单隐层 BP误差逆传播 ├── ch06_svm/ │ └── smo.py # 简化版 SMO 求解对偶问题 ├── ch07_bayes/ │ └── naive_bayes.py # 朴素贝叶斯 拉普拉斯平滑 ├── ch08_ensemble/ │ ├── adaboost.py # 加权错误率与样本权重更新 │ └── bagging.py # 自助采样 基学习器投票 ├── ch09_cluster/ │ ├── kmeans.py │ ├── gmm.py # EM 求解高斯混合 │ └── dbscan.py ├── ch10_dim_reduce/ │ ├── pca.py # 协方差矩阵特征分解 / SVD │ └── lda.py # 有监督的线性判别降维 ├── data/ │ └── watermelon.csv # 教材自带的小样本数据 └── utils/ ├── metrics.py # 准确率、混淆矩阵、轮廓系数 └── plot.py # 决策边界、聚类散点、收敛曲线看清这个骨架之后你会发现一个规律每个算法文件内部基本都是三段式——顶部构造或加载数据、中间是算法本体一个类或一组函数、底部if __name__ __main__:里调用并画图。这个结构看着土但对学习型代码是最优解单文件可独立运行改坏了不影响别人也方便你只关心某一个算法时把它单独拎出来。utils/的存在是个信号。如果这个目录里的metrics.py被多个章节 import说明作者已经做过一轮解耦你可以直接在这些工具函数上加自己的评估指标如果每个算法文件里都各自写了一份准确率计算那说明代码是各写各的二次开发时第一件事就是把重复逻辑抽出来。我一般拿到这类包会先做一次grep -rn def --include*.py把所有函数名列出来扫一眼。命名风格能告诉你作者的抽象粒度如果全是train_xxx/predict_xxx这种成对命名说明已经有一层朴素接口约定如果全是main、run、test1那就只能一个文件一个文件读。2.2 技术选型NumPy 手写为主sklearn 只做对照这份代码的依赖清单非常短短到有点反直觉依赖用途是否必需numpy矩阵运算、linalg求逆与特征分解必需matplotlib决策边界、收敛曲线、聚类散点必需不需要绘图时可换 Agg 后端scipyspecial.expit之外的优化工具、距离矩阵可选pandas读取 csv、看数据结构可选np.loadtxt能完全替代scikit-learn只用来对拍结果不参与实现可选但强烈建议装核心算法一律用 NumPy 手写这是这份练习包最大的价值点。sklearn 在这里的角色是裁判而不是选手你写完一个决策树用同一个数据集跑一遍sklearn.tree.DecisionTreeClassifier如果准确率差出十几个点那一定是你的划分准则或者递归终止条件写错了。这种对拍是排查算法实现错误最快的手段比盯着公式看半小时有效得多。有一点要提前说明不同年份整理出来的练习代码NumPy 版本差异会带来 API 变化。老代码里常见的np.float、np.int、np.bool在 NumPy 1.24 之后已经彻底移除np.matrix也早就不推荐用了。如果你拿到的是几年前整理的版本第一件事就是全局搜一下这几个词替换成内置的float、int、bool以及用np.ndarray配合运算符代替np.matrix的*乘法。2.3 环境准备把 Python 版本和依赖钉死不要直接在系统 Python 里装依赖。这类练习包往往对 NumPy 版本有隐性要求比如某些写法依赖 1.x 的行为装到全局环境里后面做别的项目时容易翻车。用一个干净的虚拟环境# 创建并激活虚拟环境Python 3.9~3.11 兼容性最好 python -m venv .venv source .venv/bin/activate # Windows 下是 .venv\Scripts\activate # 装最小依赖集先不锁版本跑通再冻结 pip install numpy matplotlib scipy pandas scikit-learn # 跑通之后把当前版本固化下来方便复现 pip freeze requirements.txt装完之后先做一次冒烟测试确认这几个库能正常 import 且版本匹配import numpy as np import matplotlib matplotlib.use(Agg) # 无显示器环境服务器/容器必须设置否则 import pyplot 会卡住或报错 import matplotlib.pyplot as plt print(np.__version__, matplotlib.__version__) print(np.linalg.inv(np.array([[2.0, 1.0], [1.0, 3.0]]))) # 验证 linalg 可用matplotlib.use(Agg)这一行值得单独拎出来说。很多人在服务器上跑练习代码卡在import matplotlib.pyplot as plt十几秒甚至直接报no display name其实不是代码问题是默认后端在找图形界面。设置成Agg之后只出文件不弹窗plt.savefig()照常工作。如果要在 Jupyter 里看图再单独用%matplotlib inline两者不冲突。环境准备好之后建议按从线性模型到集成学习的顺序跑而不是按目录顺序。线性回归和对数几率回归是最短的完整闭环构造数据 → 定义损失 → 优化 → 评估 → 画图跑通这两个后面所有算法的阅读模式都是一样的。3. 核心算法模块拆解线性模型、对数几率回归与决策树这三块是整份练习代码的地基。后面的 SVM 对偶问题、AdaBoost 的指数损失、神经网络的链式求导全都能在这三章里找到原型。3.1 线性回归闭式解和梯度下降两条路都得会线性回归的代码一般会写两遍一遍用正规方程闭式解一遍用梯度下降。不是为了炫技是因为这两条路对应的数学直觉完全不同闭式解告诉你最优解长什么样梯度下降告诉你怎么一步步走到那。import numpy as np def linear_regression_closed_form(X, y): 最小二乘闭式解w (X^T X)^{-1} X^T y X: (n_samples, n_features1)第一列已补 1 y: (n_samples,) # 用 pinv 而不是 inv特征共线时 inv 会炸pinv 给最小范数解 w np.linalg.pinv(X.T X) X.T y return w def linear_regression_gd(X, y, lr0.01, epochs1000, tol1e-8): 批量梯度下降版返回权重和每轮损失方便画收敛曲线 n, d X.shape w np.zeros(d) losses [] for _ in range(epochs): y_pred X w err y_pred - y loss (err err) / (2 * n) # 均方误差除以 2 是为了求导后系数变干净 losses.append(loss) grad X.T err / n # dL/dw X^T(Xw - y)/n w - lr * grad if np.linalg.norm(grad) tol: # 梯度足够小就提前退出 break return w, losses参数上最容易出问题的是lr。对未归一化的特征0.01 这个学习率可能直接发散成nan因为损失对某个量纲大的特征求导后数值巨大一步就跨过了谷底。我一般会在gd函数内部先做一次特征标准化减均值除标准差训练完再把权重按同样比例还原回去这样学习率的取值就稳定在 0.01~0.1 这个区间不用一个数据集试一次。np.linalg.pinv代替inv是个血泪经验。教材里的数据往往只有两个特征、十几条样本X^T X很容易接近奇异比如两个特征线性相关或者样本数少于特征数用inv会抛LinAlgError: Singular matrix用pinv走 SVD 求伪逆得到的是最小范数解不会崩。两种写法在良态数据上结果完全一致你没有理由不用更稳的那个。收敛曲线一定要画。只看最终 loss 的数字你没法区分已经收敛和在震荡画出来一目了然——正常情况下应该是单调下降然后变平如果出现锯齿说明学习率偏大如果一千轮还在缓慢下降说明轮数不够或者学习率偏小。3.2 对数几率回归为什么这里用牛顿法而不是梯度下降对数几率回归名字叫回归干的是分类的活的代码通常会给出两种优化方式教材里重点讲的是牛顿法。原因很实际这个问题的目标函数是二阶可导的凸函数牛顿法利用二阶信息收敛速度是二次的几轮就能收敛梯度下降只有一阶信息在数据条件数不好的时候慢得让人怀疑人生。import numpy as np def sigmoid(z): # 数值稳定的写法避免 exp 溢出z 很大时 exp(-z) 上溢 return np.where(z 0, 1.0 / (1.0 np.exp(-z)), np.exp(z) / (1.0 np.exp(z))) def logistic_newton(X, y, epochs20, tol1e-8): 牛顿法求解对数几率回归 X: (n, d) 已补 1 列 y: (n,) 取值 0/1 n, d X.shape beta np.zeros(d) # 待求参数 (b; w) for i in range(epochs): p sigmoid(X beta) # 预测为正类的概率 grad X.T (p - y) / n # 一阶导 # 二阶导海森矩阵X^T diag(p(1-p)) X W np.diag(p * (1 - p)) H X.T W X / n # 加一个极小的对角项防止海森矩阵病态时求逆失败 step np.linalg.solve(H 1e-8 * np.eye(d), grad) beta - step if np.linalg.norm(step) tol: # 参数几乎不再变化就停 break return beta两个细节值得掰开说。第一是sigmoid里的分支写法直接写1 / (1 np.exp(-z))当z -800时np.exp(800)会溢出成inf接着整条损失链变成nan。分成正负两支处理后exp的输入永远是负数不会上溢。这个坑在小数据集上几乎遇不到一旦换成真实数据、特征没归一化第一次迭代就能炸。第二是海森矩阵那一行。教材上的形式是X^T R X其中R是对角阵对角元是p_i(1-p_i)。用np.diag显式构造一个n×n的矩阵再相乘在样本量大时会吃掉大量内存。规范的写法是直接用广播计算(X * (p * (1 - p))[:, None]).T X避免构造稠密对角阵。练习代码为了可读性通常选前者你自己扩展的时候记得换掉。牛顿法的停止条件用的是参数变化量||step||而不是梯度范数。这在接近最优解时更灵敏因为梯度可能已经很小但参数还在缓慢漂移。3.3 决策树递归划分与剪枝的落点决策树的代码是整个练习包里最容易写错、也最容易写对之后最有成就感的一块。核心就三件事选哪个特征划分、怎么递归、什么时候停。import numpy as np from collections import Counter def entropy(y): 信息熵用 log2 保持和教材一致单位 bit probs np.array(list(Counter(y).values())) / len(y) return -np.sum(probs * np.log2(probs 1e-12)) def info_gain(X, y, feat_idx): 按第 feat_idx 个特征做二分划分的信息增益 base entropy(y) vals np.unique(X[:, feat_idx]) if len(vals) 2: return 0.0, None # 连续特征取相邻取值中点作为候选切分点遍历找最优 best_gain, best_thr -1, None for i in range(len(vals) - 1): thr (vals[i] vals[i 1]) / 2 left, right y[X[:, feat_idx] thr], y[X[:, feat_idx] thr] if len(left) 0 or len(right) 0: continue cond len(left) / len(y) * entropy(left) len(right) / len(y) * entropy(right) gain base - cond if gain best_gain: best_gain, best_thr gain, thr return best_gain, best_thr def build_tree(X, y, depth0, max_depth5, min_samples2): 递归建树返回嵌套字典表示的树 # 终止条件一全部同类 if len(set(y)) 1: return {leaf: True, label: y[0]} # 终止条件二达到深度上限或样本太少预剪枝 if depth max_depth or len(y) min_samples: return {leaf: True, label: Counter(y).most_common(1)[0][0]} best (-1, None, None) for f in range(X.shape[1]): gain, thr info_gain(X, y, f) if gain best[0]: best (gain, f, thr) # 找不到有效划分增益为 0 或所有特征取值唯一退化成叶节点 if best[1] is None or best[0] 1e-9: return {leaf: True, label: Counter(y).most_common(1)[0][0]} _, f, thr best mask X[:, f] thr return { leaf: False, feat: f, thr: thr, left: build_tree(X[mask], y[mask], depth 1, max_depth, min_samples), right: build_tree(X[~mask], y[~mask], depth 1, max_depth, min_samples), }这段代码里有三个点是你在别人的实现里也一定会撞见的。entropy里的 1e-12是为了防止log2(0)。理论上当某一类占比为 0 时该类别根本不进Counter但浮点误差会让概率算出极小的正数log之后变成大负数最后熵可能出现-0.0000。加上这个 epsilon结果稳定在 0 附近。连续特征的处理用的是排序后取相邻中点这是 C4.5 的做法代价是每个特征要扫O(n)个候选切分点每算一次增益都是O(n)整体O(n^2)。教材数据集只有十几条样本随便跑换成几百条以上最好先对特征做一次分箱把候选点降到十几二十个。预剪枝用max_depth和min_samples两个参数控制写在递归入口。后剪枝先长满再自底向上砍在练习代码里通常不实现因为它需要一个独立的验证集来评估剪之前和剪之后哪个泛化更好工程量比建树本身还大。真想练手的话我的做法是留出 20% 数据当验证集写一个prune(node, X_val, y_val)递归函数用剪枝后的验证准确率是否下降作为判据——这个改动只涉及一个函数但能让你彻底理解为什么教材里说预剪枝有欠拟合风险。4. 核心算法模块拆解SVM、集成学习、聚类与降维过了线性模型这道坎后面的算法在代码结构上是同一套模板一个fit循环 一个predict分支。差别在于循环里更新什么和用什么指标判断收敛。4.1 支持向量机简化版 SMO 的两个循环怎么读SVM 的代码是整份练习包里最吓人的翻开就是两个嵌套循环加一堆if看着像天书。但它的本质很朴素每次挑两个拉格朗日乘子固定其他所有乘子把这两个调到一个让目标函数最小的位置反复直到满足 KKT 条件。import numpy as np def smo_simple(X, y, C1.0, tol1e-3, max_iter40): 简化版 SMO X: (n, d)y: (n,) 取值 -1/1 返回 alpha, b以及支持向量的索引 n, d X.shape alpha np.zeros(n) b 0.0 K X X.T # 线性核先算好避免循环里重复计算 it 0 while it max_iter: alpha_pairs_changed 0 for i in range(n): # 第 i 个样本的预测值sum_j alpha_j y_j K[j,i] b Ei (alpha * y) K[:, i] b - y[i] # KKT 违反判定乘子在边界内却偏离间隔或在边界上却没贴合 if (y[i] * Ei -tol and alpha[i] C) or (y[i] * Ei tol and alpha[i] 0): j i while j i: # 随便挑一个不等于 i 的 j j np.random.randint(n) Ej (alpha * y) K[:, j] b - y[j] ai_old, aj_old alpha[i], alpha[j] # 根据 y_i 与 y_j 是否同号计算 alpha_j 的上下界 L 和 H if y[i] ! y[j]: L, H max(0, aj_old - ai_old), min(C, C aj_old - ai_old) else: L, H max(0, ai_old aj_old - C), min(C, ai_old aj_old) if L H: continue # 目标函数对 alpha_j 的二阶导二阶差分形式 eta 2 * K[i, j] - K[i, i] - K[j, j] if eta 0: # 核函数不满足 Mercer 条件时的兜底 continue aj_new aj_old - y[j] * (Ei - Ej) / eta aj_new np.clip(aj_new, L, H) # 裁剪回合法区间 if abs(aj_new - aj_old) 1e-5: continue ai_new ai_old y[i] * y[j] * (aj_old - aj_new) # 由约束 sum alpha*y0 推出 alpha[i], alpha[j] ai_new, aj_new # 重新计算偏置 b两个乘子都在界内时取平均否则取其中一个 b1 b - Ei - y[i] * (ai_new - ai_old) * K[i, i] - y[j] * (aj_new - aj_old) * K[j, i] b2 b - Ej - y[i] * (ai_new - ai_old) * K[i, j] - y[j] * (aj_new - aj_old) * K[j, j] if 0 ai_new C: b b1 elif 0 aj_new C: b b2 else: b (b1 b2) / 2 alpha_pairs_changed 1 # 一整轮没有任何乘子被更新说明已经收敛 it it 1 if alpha_pairs_changed 0 else 0 sv np.where(alpha 1e-6)[0] return alpha, b, sv读这段代码的关键是抓住不变量sum(alpha_i * y_i) 0和0 alpha_i C必须始终成立。所以aj_new算出来之后必须立刻裁剪到[L, H]区间再反推ai_new。如果顺序反了或者裁剪时忘了考虑y_i y_j的符号跑出来的alpha会越界最后预测全是一类——这个现象很典型后面排查章节还会再提。eta 0这个判断容易被忽略。线性核下eta恒为负永远不会触发一旦你换成rbf核但核矩阵写得不对比如忘了距离平方eta就可能非负此时除法直接让alpha变成inf。加这一行兜底的成本极低。tol控制 KKT 违反的容忍度C控制软间隔的惩罚力度。数据集线性可分时把C调大接近硬间隔数据有噪声时把C调小允许更多样本落在间隔内。练习时我建议先把C设成 1.0 跑一遍画出决策边界然后改成 0.1 和 10 各跑一遍看支持向量数量怎么变——支持向量越多说明模型对噪声越宽容这个直观感受比看公式管用。4.2 集成学习AdaBoost 的权重更新与 Bagging 的自助采样集成学习这两段代码放在一个目录里特别合适因为它们刚好代表了两种相反的思路AdaBoost 是串行的每个基学习器盯着前一个犯的错改Bagging 是并行的每个基学习器看一份随机抽出来的数据。import numpy as np def adaboost_stump(X, y, n_estimators10): AdaBoost 决策树桩单层决策树 y 取值 -1/1 n len(y) w np.full(n, 1.0 / n) # 初始样本权重均匀 models [] # 保存 (特征, 阈值, 方向, alpha) for _ in range(n_estimators): best {err: np.inf} # 树桩就是找一个特征一个阈值把样本切成两半 for f in range(X.shape[1]): for thr in np.unique(X[:, f]): for direction in (-1, 1): pred np.ones(n) pred[direction * X[:, f] direction * thr] -1 err w[pred ! y].sum() # 加权错误率 if err best[err]: best {err: err, feat: f, thr: thr, dir: direction} # 错误率 0.5 说明这个树桩还不如瞎猜直接终止 if best[err] 0.5 - 1e-9: break eps max(best[err], 1e-10) # 防止 log(0) alpha 0.5 * np.log((1 - eps) / eps) # 基学习器权重 pred np.ones(n) pred[best[dir] * X[:, best[feat]] best[dir] * best[thr]] -1 w w * np.exp(-alpha * y * pred) # 分错的样本权重上升 w w / w.sum() # 必须重新归一化成概率分布 models.append({**best, alpha: alpha}) return models def predict_adaboost(models, X): 加权投票sign(sum alpha_m * h_m(x)) agg np.zeros(len(X)) for m in models: pred np.ones(len(X)) pred[m[dir] * X[:, m[feat]] m[dir] * m[thr]] -1 agg m[alpha] * pred return np.sign(agg)权重更新那一行w w * exp(-alpha * y * pred)是 AdaBoost 的核心。当预测正确时y*pred 1指数项小于 1权重被压低预测错误时y*pred -1指数项大于 1权重被抬高。归一化那一行绝对不能省——省了之后权重总和会指数级衰减几轮之后就全部下溢成 0err计算出来永远是 0模型直接在第一个基学习器上完美收敛你还会以为是代码写对了。err 0.5的提前终止也值得留意。在二分类问题里基学习器的错误率超过一半意味着它比随机猜还差此时alpha会变成负数加权投票的逻辑就反过来了。教材的写法是直接跳出循环更严谨的做法是把这类基学习器整个丢弃然后继续找下一个。Bagging 那边的代码短得多核心只有一句自助采样def bagging_fit(X, y, base_learner, n_estimators10, seed42): 自助采样有放回地抽 n 条约 63.2% 的样本会被抽到至少一次 rng np.random.RandomState(seed) models [] for _ in range(n_estimators): idx rng.choice(len(X), sizelen(X), replaceTrue) # 有放回 models.append(base_learner().fit(X[idx], y[idx])) return modelsreplaceTrue是 Bagging 和普通交叉验证的分水岭。每轮抽样大约有1 - (1-1/n)^n ≈ 63.2%的样本被抽中剩下 36.8% 的袋外样本天然构成一个验证集可以直接拿来估泛化误差不需要额外划分。这个性质在练习代码里往往没写出来但它正是 Bagging 比单模型稳的数学来源。想验证的话在bagging_fit里顺便把没被抽中的索引存下来用这些样本算一遍准确率你会发现它和留出法估出来的数字非常接近。4.3 聚类k-means 的初始化陷阱与 GMM 的 EM 两步k-means 的代码是全篇最短的二十行能写完但它的失败模式也是最多的。import numpy as np def kmeans(X, k, max_iter100, seed42): rng np.random.RandomState(seed) # k-means 初始化第一个中心随机之后按距离平方成正比的概率选 centers [X[rng.randint(len(X))]] for _ in range(1, k): d2 np.min([np.sum((X - c) ** 2, axis1) for c in centers], axis0) probs d2 / d2.sum() centers.append(X[rng.choice(len(X), pprobs)]) centers np.array(centers) for _ in range(max_iter): # E 步每个样本归到最近的中心 dist np.linalg.norm(X[:, None, :] - centers[None, :, :], axis2) labels dist.argmin(axis1) # M 步中心更新为簇内均值 new_centers np.array([ X[labels j].mean(axis0) if np.any(labels j) else centers[j] for j in range(k) ]) if np.allclose(new_centers, centers, atol1e-6): break centers new_centers return labels, centers两个地方必须防。一是空簇如果某个中心初始化得离所有样本都很远E 步之后没有任何样本归它M 步算mean会得到nan然后整个矩阵污染。上面代码里用if np.any(labels j)兜底保留原中心至少不会崩。二是初始化敏感纯随机挑 k 个样本点当初始中心跑十次能得到十种不同的聚类结果。k-means把初始中心的距离拉开方差能降一个数量级多写四行很值。GMM 用的是 EM 算法和 k-means 是同一个套路的软版本k-means 里每个样本硬性归到一个簇GMM 里每个样本对每个簇都有一个归属概率。def gmm_em(X, k, max_iter100, seed42): n, d X.shape rng np.random.RandomState(seed) mu X[rng.choice(n, k, replaceFalse)] # 初始均值随机挑 k 个样本 sigma np.array([np.eye(d)] * k) # 初始协方差单位阵 pi np.full(k, 1.0 / k) # 初始混合系数均匀 for _ in range(max_iter): # E 步算每个样本属于每个高斯分量的后验概率责任度 resp np.zeros((n, k)) for j in range(k): diff X - mu[j] inv np.linalg.inv(sigma[j] 1e-6 * np.eye(d)) # 加正则防奇异 coef 1.0 / np.sqrt((2 * np.pi) ** d * np.linalg.det(sigma[j])) resp[:, j] pi[j] * coef * np.exp(-0.5 * np.sum(diff inv * diff, axis1)) resp / resp.sum(axis1, keepdimsTrue) 1e-12 # 归一化成概率 # M 步用责任度加权更新三个参数 Nk resp.sum(axis0) mu (resp.T X) / Nk[:, None] for j in range(k): diff X - mu[j] sigma[j] (resp[:, j][:, None] * diff).T diff / Nk[j] 1e-6 * np.eye(d) pi Nk / n return mu, sigma, piE 步里np.linalg.det和np.linalg.inv同时出现是这个算法最脆弱的地方。当某个分量的协方差矩阵退化比如只有一两个样本归到它行列式趋近 0coef直接冲到无穷大接着resp全是nan。所以我在inv和sigma里都加了1e-6 * I的正则项——这在教材代码里通常没有但真实数据上不加几乎必炸。resp.sum(axis1, keepdimsTrue)的归一化分母后面跟了1e-12。理论上分母是各分量概率之和不会为 0但浮点下溢之后可能真的算出 0一整行除下来全是nan。这个 epsilon 就是后悔药。想直观感受 k-means 和 GMM 的区别构造一份两个长条形且互相斜交的二维数据分别跑一遍画出来。k-means 会得到一条直线切开的两个半圆GMM 会沿着椭圆的长轴切开——因为 GMM 的协方差矩阵能表达方向k-means 的欧氏距离不能。这个对比图比任何公式解释都清楚。4.4 降维PCA 走 SVD 比走协方差矩阵更稳PCA 的教材推导是从协方差矩阵的特征分解出发的但代码实现里更推荐直接对中心化后的数据做 SVDdef pca(X, n_components): X: (n, d)返回投影后的数据和解释方差比 Xc X - X.mean(axis0) # 中心化这一步漏了结果全错 # 对 Xc 做 SVDXc U S V^TV 的前 k 列就是主成分方向 U, S, Vt np.linalg.svd(Xc, full_matricesFalse) components Vt[:n_components] X_proj Xc components.T # 解释方差比 每个主成分对应的奇异值平方 / 全部奇异值平方和 explained (S ** 2) / np.sum(S ** 2) return X_proj, explained[:n_components], components为什么绕开协方差矩阵因为X^T X的构造过程本身就会放大数值误差——特征量纲差几个数量级时X^T X的条件数直接变成量级的平方特征分解出来的小特征值对应的方向基本是噪声。SVD 直接在原矩阵上算数值稳定性高得多而且full_matricesFalse时复杂度是O(nd·min(n,d))比构造d×d协方差再分解更省。中心化那一行是硬性要求。忘了减均值第一主成分会指向数据中心的方向而不是方差最大的方向画出来的图看着降维了其实降的是均值。这个错误在练习代码里极其常见因为教材数据集的两个特征均值本来就接近 0漏了也看不出来。explained数组是判断取几个主成分的依据。一般取累计解释方差达到 85%~95% 的维度数。教材代码里通常硬编码n_components1因为西瓜数据只有两个特征降到一维好画图用到真实数据上得改成按阈值自动选。5. 避坑与排查跑不通、收敛慢、结果不对的常见原因这份练习包的代码量不大但失败方式相当丰富。下面几条是我自己在跑、在帮别人看代码时反复遇到的按现象 → 原因 → 解决记下来。5.1 现象脚本一跑就SyntaxError或者除法结果全是 0原因基本只有一个——代码是从很早期的版本整理来的写的时候是 Python 2。print acc:, acc少了括号a / b在整数之间做的是整除range()返回的是列表而不是迭代器。这三处会让代码要么直接报语法错要么更隐蔽地算出全 0 的准确率。解决方式是一次性扫一遍# 找出所有没有括号的 printPython 2 风格 grep -rn print [^(] --include*.py . # 找出可疑的除法重点看整数相除的评估指标 grep -rn / --include*.py . | grep -v from __future__最省事的做法是在每个文件顶部加一行from __future__ import division, print_function能同时修好除法和 print。不过这行只在 Python 2 下有意义Python 3 里虽然不报错但也没作用长期看还是老老实实改代码。评估指标里的correct / total是最容易踩的地方——correct和total都是int整除之后准确率永远是 0 或 1你会以为模型彻底没学到东西。5.2 现象图能出来但标题和图例全是一堆方框原因很明确matplotlib 的默认字体不含中文遇到中文字符就渲染成tofu豆腐块。这个不影响计算结果但如果你要拿图给别人看方块图基本等于废图。解决方式是显式指定一个系统里存在的中文字体并且顺手把负号显示修好import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt # 按操作系统选一个存在的字体Windows 用 SimHeimacOS 用 PingFang SCLinux 用 Noto Sans CJK plt.rcParams[font.sans-serif] [SimHei, PingFang SC, Noto Sans CJK SC, DejaVu Sans] plt.rcParams[axes.unicode_minus] False # 否则负号会显示成方块字体列表按优先级排matplotlib 会挑第一个能找到的。最后的DejaVu Sans是保底至少在都没装的情况下不会报警告。axes.unicode_minus那行很多人会漏结果坐标轴上-1、-0.5这些刻度全变成方块看着像字体没生效其实是负号用的另一个字符集。5.3 现象决策树递归到一半报RecursionError或者建出来的树深得离谱原因是划分条件设计得不严。教材里的伪代码写的是对每个特征的每个取值划分选信息增益最大的但没写清楚两件事一是划分后某个子集为空怎么办二是所有特征的信息增益都为 0或者特征已经用完怎么办。少了这两个判断递归会在同一份数据上反复切分永远到不了叶节点。解决方式是在build_tree里补两道闸门划分前检查左右子集是否都非空划分后检查增益是否大于一个极小阈值比如1e-9任一条件不满足就地转成叶节点标签取当前集合的众数。除此之外再加一个max_depth兜底即使逻辑写错了也不至于把栈打爆。这三道闸门我都建议同时保留——它们防的是不同的病空子集防死循环零增益防无意义分裂深度上限防意外。5.4 现象训练损失一开始是nan或者几轮之后突然变成nan这是数值问题里最常见的一类来源有三个。第一个是指数溢出sigmoid里np.exp(-z)在z很负时上溢softmax 里减最大值那一步漏了。第二个是除零高斯混合模型里协方差矩阵行列式为 0或者 k-means 里出现空簇。第三个是学习率过大梯度下降一步跨过谷底损失指数增长然后溢出。排查顺序建议从数据入手而不是从代码入手。先在训练前打印X.mean()、X.std()、X.max()如果标准差比均值大两三个数量级先做标准化大部分nan就此消失。数据没问题再怀疑学习率把它调小十倍跑十轮看损失是不是稳定下降。最后才去加1e-8之类的 epsilon——那是兜底不是治疗。一个具体的习惯在fit循环里每轮检查np.isfinite(loss)一旦不是有限值就打印当前轮次和参数直接raise。让错误在你指定的地方爆出来比让它漂移到评估阶段变成一串莫名其妙的数字要好排查得多。5.5 现象SVM 跑完预测全是一类或者alpha大部分顶到C先看sum(alpha * y)是不是接近 0。如果不接近说明等式约束被破坏了问题一定出在ai_new的推导上——最常见的是把ai_old y_i * y_j * (aj_old - aj_new)里的符号写反了或者把y[i] ! y[j]和y[i] y[j]两个分支的L、H公式对调了。如果等式约束没问题但alpha大量顶到C检查tol和C的取值。C太小比如 0.01时所有乘子都被上界卡住模型退化成硬性限制的线性分类器稍微有点噪声的数据就会全预测成一类。先把C调到 1.0 或 10.0用线性可分的数据验证一遍逻辑再回到真实数据上调。b的更新也值得单独验一次。当两个乘子都落在(0, C)内部时用b1或b2都行都落在边界上时取两者的平均。如果不论什么情况都只取b1在某些数据分布下偏置会偏掉表现出来就是决策边界整体平移准确率莫名低几个点——这种偏差不打印中间量根本看不出来。6. 二次开发把这个练习包改成能跑实验的骨架这份代码最大的价值不是能跑通而是它给了你每个算法最朴素的实现。想再往前走一步我建议做三件事成本都不高但能让它从练手脚本集合变成你自己的实验底座。第一件事是统一接口。现在每个文件的函数签名五花八门有的返回(w, losses)有的返回(labels, centers)。抽一个极简基类只约定fit和predict两个方法之后所有算法都往外层套一层import numpy as np class BaseModel: 所有手写模型的最小接口约定不依赖 sklearn 也能对拍 def fit(self, X, yNone): raise NotImplementedError def predict(self, X): raise NotImplementedError def score(self, X, y): return float(np.mean(self.predict(X) y)) class LogisticRegression(BaseModel): def __init__(self, lr0.1, epochs500): self.lr, self.epochs lr, epochs def fit(self, X, y): self.w_ np.zeros(X.shape[1] 1) Xb np.hstack([np.ones((len(X), 1)), X]) # 补一列 1 对应偏置 for _ in range(self.epochs): p 1 / (1 np.exp(-np.clip(Xb self.w_, -30, 30))) self.w_ - self.lr * Xb.T (p - y) / len(y) return self def predict(self, X): Xb np.hstack([np.ones((len(X), 1)), X]) return (Xb self.w_ 0).astype(int)这个改动的收益在验证环节立刻体现。有了统一接口你可以写一个run_experiment(models, X, y)把 k 折交叉验证、准确率、耗时一次性跑完所有算法共用同一套评估逻辑不会再出现这个文件用留出法、那个文件用全量评估导致的横向对比失真。第二件事是加对拍测试。这是我认为性价比最高的一步给自己写的每个模型配一个断言和 sklearn 的同名实现比准确率差太多就 fail。import numpy as np from sklearn.tree import DecisionTreeClassifier from sklearn.datasets import load_iris from sklearn.model_selection import train_test_split def test_tree_matches_sklearn(): X, y load_iris(return_X_yTrue) Xtr, Xte, ytr, yte train_test_split(X, y, test_size0.3, random_state0) mine build_tree(Xtr, ytr, max_depth4) # 这里需要写一个从嵌套字典做预测的辅助函数 tree_predict acc_mine np.mean(tree_predict(mine, Xte) yte) acc_ref DecisionTreeClassifier(max_depth4, random_state0).fit(Xtr, ytr).score(Xte, yte) # 手写版用的是信息增益ID3sklearn 默认是 gini允许 15 个点的差距 assert acc_mine acc_ref - 0.15, f手写 {acc_mine:.3f} vs sklearn {acc_ref:.3f}那个 15 个点的容差不是随便定的。ID3 用信息增益、CART 用基尼指数划分准则不同长出来的树本来就不一样要求数值完全相等是不现实的。但只要差距在这个量级以内就说明你的递归、终止条件、预测路径三处逻辑没有结构性错误。真正的 bug比如划分时用错了特征索引会让准确率直接掉到 0.33三分类瞎猜的水平一眼就能看出来。第三件事是把参数从硬编码里解放出来。现在的代码里C1.0、k3、lr0.01全都是写死的想换个值得改源码。抽成argparse或者一个 YAML 配置一次实验对应一行命令python run.py --algo svm --C 10 --kernel linear --data data/watermelon.csv python run.py --algo kmeans --k 4 --init kmeans --seed 0改完这三步你会得到一个有点粗糙但完全自主可控的实验框架。它的算法实现肯定不如成熟库快但每一步你都能打断点、能改公式、能打印中间矩阵——这种掌控感是调包调不出来的。说个我自己的教训。刚开始跑这些练习代码时我习惯直接改原文件做实验改完也不留备份结果有一次把一个for循环的范围从range(n)改成range(d)决策树准确率从 0.8 掉到 0.3我花了整个晚上怀疑是信息增益公式推错了最后发现是那个范围。从那以后我每次动这类学习型代码之前都强制先git init提交一次原始版本再开一个分支改跑出来的结果和原版对不上就直接git diff看差异。这个习惯救过我很多次也希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/11/6 10:12:58 · 阅读 1.2k NumPy手写经典机器学习算法练习包：原理拆解与避坑实战
简介这份代码练习包面向正在学习周志华《机器学习》的读者适合希望把教材公式落到代码层面的初学者与课程实验者用于巩固经典算法的实现思路、训练流程与评估方式。压缩包共16个文件以15个Python脚本为主另附1份README说明整体约20KB体量轻便便于按章节逐个运行和修改。内容围绕线性模型、决策树、集成学习、支持向量机以及偏差—方差分析等主题展开包含逻辑回归、留一法与交叉验证评估、LDA、Iterated Bagging、MultiBoosting、简化SVM与SVR对比等练习可帮助读者对照教材理解算法细节、观察参数与泛化表现。已有510人学习适合作为课程配套练习和自学复现的参考。1. 从一行for循环开始的机器学习练习第一次把自己手写的感知机训练过程打印成日志看着损失从 0.69 一路抖到 0.02那种感觉和调通一个线上接口完全不是一回事——调接口你知道成功就在眼前训练模型你只能靠数字判断它有没有在学。这个练习包做的就是这件事它把一本经典机器学习教材里每一章的算法用 NumPy 从零写一遍线性回归、对数几率回归、决策树、支持向量机、AdaBoost、k-means、PCA 一个不落每个文件都能单独跑出数值和图表。它适合两类人。一类是刚啃完理论、怀疑自己对公式的理解有偏差的跑一遍就知道梯度推得对不对另一类是工作几年、调包调出手感但心里没底的想看看那些fit()背后到底在算什么。这个包不负责教你调参它负责让你把拉格朗日乘子、信息增益、EM 的 E 步真正写成能跑的循环。2. 练习包的目录结构与技术选型动手之前先摸清骨架拿到一个算法练习包最忌讳的操作是直接python xxx.py跑一遍报错了再从头翻。正确顺序是先看目录——目录结构决定了这份代码是能复用的工程还是一次性脚本两者的阅读方式完全不同。2.1 从章节目录看代码组织一份可预测的文件地图这类配套代码的组织方式高度趋同基本是一章一个目录、一个算法一个文件。整理出来大致长这样ml-practice/ ├── ch03_linear_model/ │ ├── linear_regression.py # 最小二乘 梯度下降 │ ├── logistic_regression.py # 对数几率回归牛顿法 │ └── lda.py # 线性判别分析 ├── ch04_decision_tree/ │ └── tree.py # ID3/C4.5 的划分与递归建树 ├── ch05_neural_network/ │ └── bp.py # 单隐层 BP误差逆传播 ├── ch06_svm/ │ └── smo.py # 简化版 SMO 求解对偶问题 ├── ch07_bayes/ │ └── naive_bayes.py # 朴素贝叶斯 拉普拉斯平滑 ├── ch08_ensemble/ │ ├── adaboost.py # 加权错误率与样本权重更新 │ └── bagging.py # 自助采样 基学习器投票 ├── ch09_cluster/ │ ├── kmeans.py │ ├── gmm.py # EM 求解高斯混合 │ └── dbscan.py ├── ch10_dim_reduce/ │ ├── pca.py # 协方差矩阵特征分解 / SVD │ └── lda.py # 有监督的线性判别降维 ├── data/ │ └── watermelon.csv # 教材自带的小样本数据 └── utils/ ├── metrics.py # 准确率、混淆矩阵、轮廓系数 └── plot.py # 决策边界、聚类散点、收敛曲线看清这个骨架之后你会发现一个规律每个算法文件内部基本都是三段式——顶部构造或加载数据、中间是算法本体一个类或一组函数、底部if __name__ __main__:里调用并画图。这个结构看着土但对学习型代码是最优解单文件可独立运行改坏了不影响别人也方便你只关心某一个算法时把它单独拎出来。utils/的存在是个信号。如果这个目录里的metrics.py被多个章节 import说明作者已经做过一轮解耦你可以直接在这些工具函数上加自己的评估指标如果每个算法文件里都各自写了一份准确率计算那说明代码是各写各的二次开发时第一件事就是把重复逻辑抽出来。我一般拿到这类包会先做一次grep -rn def --include*.py把所有函数名列出来扫一眼。命名风格能告诉你作者的抽象粒度如果全是train_xxx/predict_xxx这种成对命名说明已经有一层朴素接口约定如果全是main、run、test1那就只能一个文件一个文件读。2.2 技术选型NumPy 手写为主sklearn 只做对照这份代码的依赖清单非常短短到有点反直觉依赖用途是否必需numpy矩阵运算、linalg求逆与特征分解必需matplotlib决策边界、收敛曲线、聚类散点必需不需要绘图时可换 Agg 后端scipyspecial.expit之外的优化工具、距离矩阵可选pandas读取 csv、看数据结构可选np.loadtxt能完全替代scikit-learn只用来对拍结果不参与实现可选但强烈建议装核心算法一律用 NumPy 手写这是这份练习包最大的价值点。sklearn 在这里的角色是裁判而不是选手你写完一个决策树用同一个数据集跑一遍sklearn.tree.DecisionTreeClassifier如果准确率差出十几个点那一定是你的划分准则或者递归终止条件写错了。这种对拍是排查算法实现错误最快的手段比盯着公式看半小时有效得多。有一点要提前说明不同年份整理出来的练习代码NumPy 版本差异会带来 API 变化。老代码里常见的np.float、np.int、np.bool在 NumPy 1.24 之后已经彻底移除np.matrix也早就不推荐用了。如果你拿到的是几年前整理的版本第一件事就是全局搜一下这几个词替换成内置的float、int、bool以及用np.ndarray配合运算符代替np.matrix的*乘法。2.3 环境准备把 Python 版本和依赖钉死不要直接在系统 Python 里装依赖。这类练习包往往对 NumPy 版本有隐性要求比如某些写法依赖 1.x 的行为装到全局环境里后面做别的项目时容易翻车。用一个干净的虚拟环境# 创建并激活虚拟环境Python 3.9~3.11 兼容性最好 python -m venv .venv source .venv/bin/activate # Windows 下是 .venv\Scripts\activate # 装最小依赖集先不锁版本跑通再冻结 pip install numpy matplotlib scipy pandas scikit-learn # 跑通之后把当前版本固化下来方便复现 pip freeze requirements.txt装完之后先做一次冒烟测试确认这几个库能正常 import 且版本匹配import numpy as np import matplotlib matplotlib.use(Agg) # 无显示器环境服务器/容器必须设置否则 import pyplot 会卡住或报错 import matplotlib.pyplot as plt print(np.__version__, matplotlib.__version__) print(np.linalg.inv(np.array([[2.0, 1.0], [1.0, 3.0]]))) # 验证 linalg 可用matplotlib.use(Agg)这一行值得单独拎出来说。很多人在服务器上跑练习代码卡在import matplotlib.pyplot as plt十几秒甚至直接报no display name其实不是代码问题是默认后端在找图形界面。设置成Agg之后只出文件不弹窗plt.savefig()照常工作。如果要在 Jupyter 里看图再单独用%matplotlib inline两者不冲突。环境准备好之后建议按从线性模型到集成学习的顺序跑而不是按目录顺序。线性回归和对数几率回归是最短的完整闭环构造数据 → 定义损失 → 优化 → 评估 → 画图跑通这两个后面所有算法的阅读模式都是一样的。3. 核心算法模块拆解线性模型、对数几率回归与决策树这三块是整份练习代码的地基。后面的 SVM 对偶问题、AdaBoost 的指数损失、神经网络的链式求导全都能在这三章里找到原型。3.1 线性回归闭式解和梯度下降两条路都得会线性回归的代码一般会写两遍一遍用正规方程闭式解一遍用梯度下降。不是为了炫技是因为这两条路对应的数学直觉完全不同闭式解告诉你最优解长什么样梯度下降告诉你怎么一步步走到那。import numpy as np def linear_regression_closed_form(X, y): 最小二乘闭式解w (X^T X)^{-1} X^T y X: (n_samples, n_features1)第一列已补 1 y: (n_samples,) # 用 pinv 而不是 inv特征共线时 inv 会炸pinv 给最小范数解 w np.linalg.pinv(X.T X) X.T y return w def linear_regression_gd(X, y, lr0.01, epochs1000, tol1e-8): 批量梯度下降版返回权重和每轮损失方便画收敛曲线 n, d X.shape w np.zeros(d) losses [] for _ in range(epochs): y_pred X w err y_pred - y loss (err err) / (2 * n) # 均方误差除以 2 是为了求导后系数变干净 losses.append(loss) grad X.T err / n # dL/dw X^T(Xw - y)/n w - lr * grad if np.linalg.norm(grad) tol: # 梯度足够小就提前退出 break return w, losses参数上最容易出问题的是lr。对未归一化的特征0.01 这个学习率可能直接发散成nan因为损失对某个量纲大的特征求导后数值巨大一步就跨过了谷底。我一般会在gd函数内部先做一次特征标准化减均值除标准差训练完再把权重按同样比例还原回去这样学习率的取值就稳定在 0.01~0.1 这个区间不用一个数据集试一次。np.linalg.pinv代替inv是个血泪经验。教材里的数据往往只有两个特征、十几条样本X^T X很容易接近奇异比如两个特征线性相关或者样本数少于特征数用inv会抛LinAlgError: Singular matrix用pinv走 SVD 求伪逆得到的是最小范数解不会崩。两种写法在良态数据上结果完全一致你没有理由不用更稳的那个。收敛曲线一定要画。只看最终 loss 的数字你没法区分已经收敛和在震荡画出来一目了然——正常情况下应该是单调下降然后变平如果出现锯齿说明学习率偏大如果一千轮还在缓慢下降说明轮数不够或者学习率偏小。3.2 对数几率回归为什么这里用牛顿法而不是梯度下降对数几率回归名字叫回归干的是分类的活的代码通常会给出两种优化方式教材里重点讲的是牛顿法。原因很实际这个问题的目标函数是二阶可导的凸函数牛顿法利用二阶信息收敛速度是二次的几轮就能收敛梯度下降只有一阶信息在数据条件数不好的时候慢得让人怀疑人生。import numpy as np def sigmoid(z): # 数值稳定的写法避免 exp 溢出z 很大时 exp(-z) 上溢 return np.where(z 0, 1.0 / (1.0 np.exp(-z)), np.exp(z) / (1.0 np.exp(z))) def logistic_newton(X, y, epochs20, tol1e-8): 牛顿法求解对数几率回归 X: (n, d) 已补 1 列 y: (n,) 取值 0/1 n, d X.shape beta np.zeros(d) # 待求参数 (b; w) for i in range(epochs): p sigmoid(X beta) # 预测为正类的概率 grad X.T (p - y) / n # 一阶导 # 二阶导海森矩阵X^T diag(p(1-p)) X W np.diag(p * (1 - p)) H X.T W X / n # 加一个极小的对角项防止海森矩阵病态时求逆失败 step np.linalg.solve(H 1e-8 * np.eye(d), grad) beta - step if np.linalg.norm(step) tol: # 参数几乎不再变化就停 break return beta两个细节值得掰开说。第一是sigmoid里的分支写法直接写1 / (1 np.exp(-z))当z -800时np.exp(800)会溢出成inf接着整条损失链变成nan。分成正负两支处理后exp的输入永远是负数不会上溢。这个坑在小数据集上几乎遇不到一旦换成真实数据、特征没归一化第一次迭代就能炸。第二是海森矩阵那一行。教材上的形式是X^T R X其中R是对角阵对角元是p_i(1-p_i)。用np.diag显式构造一个n×n的矩阵再相乘在样本量大时会吃掉大量内存。规范的写法是直接用广播计算(X * (p * (1 - p))[:, None]).T X避免构造稠密对角阵。练习代码为了可读性通常选前者你自己扩展的时候记得换掉。牛顿法的停止条件用的是参数变化量||step||而不是梯度范数。这在接近最优解时更灵敏因为梯度可能已经很小但参数还在缓慢漂移。3.3 决策树递归划分与剪枝的落点决策树的代码是整个练习包里最容易写错、也最容易写对之后最有成就感的一块。核心就三件事选哪个特征划分、怎么递归、什么时候停。import numpy as np from collections import Counter def entropy(y): 信息熵用 log2 保持和教材一致单位 bit probs np.array(list(Counter(y).values())) / len(y) return -np.sum(probs * np.log2(probs 1e-12)) def info_gain(X, y, feat_idx): 按第 feat_idx 个特征做二分划分的信息增益 base entropy(y) vals np.unique(X[:, feat_idx]) if len(vals) 2: return 0.0, None # 连续特征取相邻取值中点作为候选切分点遍历找最优 best_gain, best_thr -1, None for i in range(len(vals) - 1): thr (vals[i] vals[i 1]) / 2 left, right y[X[:, feat_idx] thr], y[X[:, feat_idx] thr] if len(left) 0 or len(right) 0: continue cond len(left) / len(y) * entropy(left) len(right) / len(y) * entropy(right) gain base - cond if gain best_gain: best_gain, best_thr gain, thr return best_gain, best_thr def build_tree(X, y, depth0, max_depth5, min_samples2): 递归建树返回嵌套字典表示的树 # 终止条件一全部同类 if len(set(y)) 1: return {leaf: True, label: y[0]} # 终止条件二达到深度上限或样本太少预剪枝 if depth max_depth or len(y) min_samples: return {leaf: True, label: Counter(y).most_common(1)[0][0]} best (-1, None, None) for f in range(X.shape[1]): gain, thr info_gain(X, y, f) if gain best[0]: best (gain, f, thr) # 找不到有效划分增益为 0 或所有特征取值唯一退化成叶节点 if best[1] is None or best[0] 1e-9: return {leaf: True, label: Counter(y).most_common(1)[0][0]} _, f, thr best mask X[:, f] thr return { leaf: False, feat: f, thr: thr, left: build_tree(X[mask], y[mask], depth 1, max_depth, min_samples), right: build_tree(X[~mask], y[~mask], depth 1, max_depth, min_samples), }这段代码里有三个点是你在别人的实现里也一定会撞见的。entropy里的 1e-12是为了防止log2(0)。理论上当某一类占比为 0 时该类别根本不进Counter但浮点误差会让概率算出极小的正数log之后变成大负数最后熵可能出现-0.0000。加上这个 epsilon结果稳定在 0 附近。连续特征的处理用的是排序后取相邻中点这是 C4.5 的做法代价是每个特征要扫O(n)个候选切分点每算一次增益都是O(n)整体O(n^2)。教材数据集只有十几条样本随便跑换成几百条以上最好先对特征做一次分箱把候选点降到十几二十个。预剪枝用max_depth和min_samples两个参数控制写在递归入口。后剪枝先长满再自底向上砍在练习代码里通常不实现因为它需要一个独立的验证集来评估剪之前和剪之后哪个泛化更好工程量比建树本身还大。真想练手的话我的做法是留出 20% 数据当验证集写一个prune(node, X_val, y_val)递归函数用剪枝后的验证准确率是否下降作为判据——这个改动只涉及一个函数但能让你彻底理解为什么教材里说预剪枝有欠拟合风险。4. 核心算法模块拆解SVM、集成学习、聚类与降维过了线性模型这道坎后面的算法在代码结构上是同一套模板一个fit循环 一个predict分支。差别在于循环里更新什么和用什么指标判断收敛。4.1 支持向量机简化版 SMO 的两个循环怎么读SVM 的代码是整份练习包里最吓人的翻开就是两个嵌套循环加一堆if看着像天书。但它的本质很朴素每次挑两个拉格朗日乘子固定其他所有乘子把这两个调到一个让目标函数最小的位置反复直到满足 KKT 条件。import numpy as np def smo_simple(X, y, C1.0, tol1e-3, max_iter40): 简化版 SMO X: (n, d)y: (n,) 取值 -1/1 返回 alpha, b以及支持向量的索引 n, d X.shape alpha np.zeros(n) b 0.0 K X X.T # 线性核先算好避免循环里重复计算 it 0 while it max_iter: alpha_pairs_changed 0 for i in range(n): # 第 i 个样本的预测值sum_j alpha_j y_j K[j,i] b Ei (alpha * y) K[:, i] b - y[i] # KKT 违反判定乘子在边界内却偏离间隔或在边界上却没贴合 if (y[i] * Ei -tol and alpha[i] C) or (y[i] * Ei tol and alpha[i] 0): j i while j i: # 随便挑一个不等于 i 的 j j np.random.randint(n) Ej (alpha * y) K[:, j] b - y[j] ai_old, aj_old alpha[i], alpha[j] # 根据 y_i 与 y_j 是否同号计算 alpha_j 的上下界 L 和 H if y[i] ! y[j]: L, H max(0, aj_old - ai_old), min(C, C aj_old - ai_old) else: L, H max(0, ai_old aj_old - C), min(C, ai_old aj_old) if L H: continue # 目标函数对 alpha_j 的二阶导二阶差分形式 eta 2 * K[i, j] - K[i, i] - K[j, j] if eta 0: # 核函数不满足 Mercer 条件时的兜底 continue aj_new aj_old - y[j] * (Ei - Ej) / eta aj_new np.clip(aj_new, L, H) # 裁剪回合法区间 if abs(aj_new - aj_old) 1e-5: continue ai_new ai_old y[i] * y[j] * (aj_old - aj_new) # 由约束 sum alpha*y0 推出 alpha[i], alpha[j] ai_new, aj_new # 重新计算偏置 b两个乘子都在界内时取平均否则取其中一个 b1 b - Ei - y[i] * (ai_new - ai_old) * K[i, i] - y[j] * (aj_new - aj_old) * K[j, i] b2 b - Ej - y[i] * (ai_new - ai_old) * K[i, j] - y[j] * (aj_new - aj_old) * K[j, j] if 0 ai_new C: b b1 elif 0 aj_new C: b b2 else: b (b1 b2) / 2 alpha_pairs_changed 1 # 一整轮没有任何乘子被更新说明已经收敛 it it 1 if alpha_pairs_changed 0 else 0 sv np.where(alpha 1e-6)[0] return alpha, b, sv读这段代码的关键是抓住不变量sum(alpha_i * y_i) 0和0 alpha_i C必须始终成立。所以aj_new算出来之后必须立刻裁剪到[L, H]区间再反推ai_new。如果顺序反了或者裁剪时忘了考虑y_i y_j的符号跑出来的alpha会越界最后预测全是一类——这个现象很典型后面排查章节还会再提。eta 0这个判断容易被忽略。线性核下eta恒为负永远不会触发一旦你换成rbf核但核矩阵写得不对比如忘了距离平方eta就可能非负此时除法直接让alpha变成inf。加这一行兜底的成本极低。tol控制 KKT 违反的容忍度C控制软间隔的惩罚力度。数据集线性可分时把C调大接近硬间隔数据有噪声时把C调小允许更多样本落在间隔内。练习时我建议先把C设成 1.0 跑一遍画出决策边界然后改成 0.1 和 10 各跑一遍看支持向量数量怎么变——支持向量越多说明模型对噪声越宽容这个直观感受比看公式管用。4.2 集成学习AdaBoost 的权重更新与 Bagging 的自助采样集成学习这两段代码放在一个目录里特别合适因为它们刚好代表了两种相反的思路AdaBoost 是串行的每个基学习器盯着前一个犯的错改Bagging 是并行的每个基学习器看一份随机抽出来的数据。import numpy as np def adaboost_stump(X, y, n_estimators10): AdaBoost 决策树桩单层决策树 y 取值 -1/1 n len(y) w np.full(n, 1.0 / n) # 初始样本权重均匀 models [] # 保存 (特征, 阈值, 方向, alpha) for _ in range(n_estimators): best {err: np.inf} # 树桩就是找一个特征一个阈值把样本切成两半 for f in range(X.shape[1]): for thr in np.unique(X[:, f]): for direction in (-1, 1): pred np.ones(n) pred[direction * X[:, f] direction * thr] -1 err w[pred ! y].sum() # 加权错误率 if err best[err]: best {err: err, feat: f, thr: thr, dir: direction} # 错误率 0.5 说明这个树桩还不如瞎猜直接终止 if best[err] 0.5 - 1e-9: break eps max(best[err], 1e-10) # 防止 log(0) alpha 0.5 * np.log((1 - eps) / eps) # 基学习器权重 pred np.ones(n) pred[best[dir] * X[:, best[feat]] best[dir] * best[thr]] -1 w w * np.exp(-alpha * y * pred) # 分错的样本权重上升 w w / w.sum() # 必须重新归一化成概率分布 models.append({**best, alpha: alpha}) return models def predict_adaboost(models, X): 加权投票sign(sum alpha_m * h_m(x)) agg np.zeros(len(X)) for m in models: pred np.ones(len(X)) pred[m[dir] * X[:, m[feat]] m[dir] * m[thr]] -1 agg m[alpha] * pred return np.sign(agg)权重更新那一行w w * exp(-alpha * y * pred)是 AdaBoost 的核心。当预测正确时y*pred 1指数项小于 1权重被压低预测错误时y*pred -1指数项大于 1权重被抬高。归一化那一行绝对不能省——省了之后权重总和会指数级衰减几轮之后就全部下溢成 0err计算出来永远是 0模型直接在第一个基学习器上完美收敛你还会以为是代码写对了。err 0.5的提前终止也值得留意。在二分类问题里基学习器的错误率超过一半意味着它比随机猜还差此时alpha会变成负数加权投票的逻辑就反过来了。教材的写法是直接跳出循环更严谨的做法是把这类基学习器整个丢弃然后继续找下一个。Bagging 那边的代码短得多核心只有一句自助采样def bagging_fit(X, y, base_learner, n_estimators10, seed42): 自助采样有放回地抽 n 条约 63.2% 的样本会被抽到至少一次 rng np.random.RandomState(seed) models [] for _ in range(n_estimators): idx rng.choice(len(X), sizelen(X), replaceTrue) # 有放回 models.append(base_learner().fit(X[idx], y[idx])) return modelsreplaceTrue是 Bagging 和普通交叉验证的分水岭。每轮抽样大约有1 - (1-1/n)^n ≈ 63.2%的样本被抽中剩下 36.8% 的袋外样本天然构成一个验证集可以直接拿来估泛化误差不需要额外划分。这个性质在练习代码里往往没写出来但它正是 Bagging 比单模型稳的数学来源。想验证的话在bagging_fit里顺便把没被抽中的索引存下来用这些样本算一遍准确率你会发现它和留出法估出来的数字非常接近。4.3 聚类k-means 的初始化陷阱与 GMM 的 EM 两步k-means 的代码是全篇最短的二十行能写完但它的失败模式也是最多的。import numpy as np def kmeans(X, k, max_iter100, seed42): rng np.random.RandomState(seed) # k-means 初始化第一个中心随机之后按距离平方成正比的概率选 centers [X[rng.randint(len(X))]] for _ in range(1, k): d2 np.min([np.sum((X - c) ** 2, axis1) for c in centers], axis0) probs d2 / d2.sum() centers.append(X[rng.choice(len(X), pprobs)]) centers np.array(centers) for _ in range(max_iter): # E 步每个样本归到最近的中心 dist np.linalg.norm(X[:, None, :] - centers[None, :, :], axis2) labels dist.argmin(axis1) # M 步中心更新为簇内均值 new_centers np.array([ X[labels j].mean(axis0) if np.any(labels j) else centers[j] for j in range(k) ]) if np.allclose(new_centers, centers, atol1e-6): break centers new_centers return labels, centers两个地方必须防。一是空簇如果某个中心初始化得离所有样本都很远E 步之后没有任何样本归它M 步算mean会得到nan然后整个矩阵污染。上面代码里用if np.any(labels j)兜底保留原中心至少不会崩。二是初始化敏感纯随机挑 k 个样本点当初始中心跑十次能得到十种不同的聚类结果。k-means把初始中心的距离拉开方差能降一个数量级多写四行很值。GMM 用的是 EM 算法和 k-means 是同一个套路的软版本k-means 里每个样本硬性归到一个簇GMM 里每个样本对每个簇都有一个归属概率。def gmm_em(X, k, max_iter100, seed42): n, d X.shape rng np.random.RandomState(seed) mu X[rng.choice(n, k, replaceFalse)] # 初始均值随机挑 k 个样本 sigma np.array([np.eye(d)] * k) # 初始协方差单位阵 pi np.full(k, 1.0 / k) # 初始混合系数均匀 for _ in range(max_iter): # E 步算每个样本属于每个高斯分量的后验概率责任度 resp np.zeros((n, k)) for j in range(k): diff X - mu[j] inv np.linalg.inv(sigma[j] 1e-6 * np.eye(d)) # 加正则防奇异 coef 1.0 / np.sqrt((2 * np.pi) ** d * np.linalg.det(sigma[j])) resp[:, j] pi[j] * coef * np.exp(-0.5 * np.sum(diff inv * diff, axis1)) resp / resp.sum(axis1, keepdimsTrue) 1e-12 # 归一化成概率 # M 步用责任度加权更新三个参数 Nk resp.sum(axis0) mu (resp.T X) / Nk[:, None] for j in range(k): diff X - mu[j] sigma[j] (resp[:, j][:, None] * diff).T diff / Nk[j] 1e-6 * np.eye(d) pi Nk / n return mu, sigma, piE 步里np.linalg.det和np.linalg.inv同时出现是这个算法最脆弱的地方。当某个分量的协方差矩阵退化比如只有一两个样本归到它行列式趋近 0coef直接冲到无穷大接着resp全是nan。所以我在inv和sigma里都加了1e-6 * I的正则项——这在教材代码里通常没有但真实数据上不加几乎必炸。resp.sum(axis1, keepdimsTrue)的归一化分母后面跟了1e-12。理论上分母是各分量概率之和不会为 0但浮点下溢之后可能真的算出 0一整行除下来全是nan。这个 epsilon 就是后悔药。想直观感受 k-means 和 GMM 的区别构造一份两个长条形且互相斜交的二维数据分别跑一遍画出来。k-means 会得到一条直线切开的两个半圆GMM 会沿着椭圆的长轴切开——因为 GMM 的协方差矩阵能表达方向k-means 的欧氏距离不能。这个对比图比任何公式解释都清楚。4.4 降维PCA 走 SVD 比走协方差矩阵更稳PCA 的教材推导是从协方差矩阵的特征分解出发的但代码实现里更推荐直接对中心化后的数据做 SVDdef pca(X, n_components): X: (n, d)返回投影后的数据和解释方差比 Xc X - X.mean(axis0) # 中心化这一步漏了结果全错 # 对 Xc 做 SVDXc U S V^TV 的前 k 列就是主成分方向 U, S, Vt np.linalg.svd(Xc, full_matricesFalse) components Vt[:n_components] X_proj Xc components.T # 解释方差比 每个主成分对应的奇异值平方 / 全部奇异值平方和 explained (S ** 2) / np.sum(S ** 2) return X_proj, explained[:n_components], components为什么绕开协方差矩阵因为X^T X的构造过程本身就会放大数值误差——特征量纲差几个数量级时X^T X的条件数直接变成量级的平方特征分解出来的小特征值对应的方向基本是噪声。SVD 直接在原矩阵上算数值稳定性高得多而且full_matricesFalse时复杂度是O(nd·min(n,d))比构造d×d协方差再分解更省。中心化那一行是硬性要求。忘了减均值第一主成分会指向数据中心的方向而不是方差最大的方向画出来的图看着降维了其实降的是均值。这个错误在练习代码里极其常见因为教材数据集的两个特征均值本来就接近 0漏了也看不出来。explained数组是判断取几个主成分的依据。一般取累计解释方差达到 85%~95% 的维度数。教材代码里通常硬编码n_components1因为西瓜数据只有两个特征降到一维好画图用到真实数据上得改成按阈值自动选。5. 避坑与排查跑不通、收敛慢、结果不对的常见原因这份练习包的代码量不大但失败方式相当丰富。下面几条是我自己在跑、在帮别人看代码时反复遇到的按现象 → 原因 → 解决记下来。5.1 现象脚本一跑就SyntaxError或者除法结果全是 0原因基本只有一个——代码是从很早期的版本整理来的写的时候是 Python 2。print acc:, acc少了括号a / b在整数之间做的是整除range()返回的是列表而不是迭代器。这三处会让代码要么直接报语法错要么更隐蔽地算出全 0 的准确率。解决方式是一次性扫一遍# 找出所有没有括号的 printPython 2 风格 grep -rn print [^(] --include*.py . # 找出可疑的除法重点看整数相除的评估指标 grep -rn / --include*.py . | grep -v from __future__最省事的做法是在每个文件顶部加一行from __future__ import division, print_function能同时修好除法和 print。不过这行只在 Python 2 下有意义Python 3 里虽然不报错但也没作用长期看还是老老实实改代码。评估指标里的correct / total是最容易踩的地方——correct和total都是int整除之后准确率永远是 0 或 1你会以为模型彻底没学到东西。5.2 现象图能出来但标题和图例全是一堆方框原因很明确matplotlib 的默认字体不含中文遇到中文字符就渲染成tofu豆腐块。这个不影响计算结果但如果你要拿图给别人看方块图基本等于废图。解决方式是显式指定一个系统里存在的中文字体并且顺手把负号显示修好import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt # 按操作系统选一个存在的字体Windows 用 SimHeimacOS 用 PingFang SCLinux 用 Noto Sans CJK plt.rcParams[font.sans-serif] [SimHei, PingFang SC, Noto Sans CJK SC, DejaVu Sans] plt.rcParams[axes.unicode_minus] False # 否则负号会显示成方块字体列表按优先级排matplotlib 会挑第一个能找到的。最后的DejaVu Sans是保底至少在都没装的情况下不会报警告。axes.unicode_minus那行很多人会漏结果坐标轴上-1、-0.5这些刻度全变成方块看着像字体没生效其实是负号用的另一个字符集。5.3 现象决策树递归到一半报RecursionError或者建出来的树深得离谱原因是划分条件设计得不严。教材里的伪代码写的是对每个特征的每个取值划分选信息增益最大的但没写清楚两件事一是划分后某个子集为空怎么办二是所有特征的信息增益都为 0或者特征已经用完怎么办。少了这两个判断递归会在同一份数据上反复切分永远到不了叶节点。解决方式是在build_tree里补两道闸门划分前检查左右子集是否都非空划分后检查增益是否大于一个极小阈值比如1e-9任一条件不满足就地转成叶节点标签取当前集合的众数。除此之外再加一个max_depth兜底即使逻辑写错了也不至于把栈打爆。这三道闸门我都建议同时保留——它们防的是不同的病空子集防死循环零增益防无意义分裂深度上限防意外。5.4 现象训练损失一开始是nan或者几轮之后突然变成nan这是数值问题里最常见的一类来源有三个。第一个是指数溢出sigmoid里np.exp(-z)在z很负时上溢softmax 里减最大值那一步漏了。第二个是除零高斯混合模型里协方差矩阵行列式为 0或者 k-means 里出现空簇。第三个是学习率过大梯度下降一步跨过谷底损失指数增长然后溢出。排查顺序建议从数据入手而不是从代码入手。先在训练前打印X.mean()、X.std()、X.max()如果标准差比均值大两三个数量级先做标准化大部分nan就此消失。数据没问题再怀疑学习率把它调小十倍跑十轮看损失是不是稳定下降。最后才去加1e-8之类的 epsilon——那是兜底不是治疗。一个具体的习惯在fit循环里每轮检查np.isfinite(loss)一旦不是有限值就打印当前轮次和参数直接raise。让错误在你指定的地方爆出来比让它漂移到评估阶段变成一串莫名其妙的数字要好排查得多。5.5 现象SVM 跑完预测全是一类或者alpha大部分顶到C先看sum(alpha * y)是不是接近 0。如果不接近说明等式约束被破坏了问题一定出在ai_new的推导上——最常见的是把ai_old y_i * y_j * (aj_old - aj_new)里的符号写反了或者把y[i] ! y[j]和y[i] y[j]两个分支的L、H公式对调了。如果等式约束没问题但alpha大量顶到C检查tol和C的取值。C太小比如 0.01时所有乘子都被上界卡住模型退化成硬性限制的线性分类器稍微有点噪声的数据就会全预测成一类。先把C调到 1.0 或 10.0用线性可分的数据验证一遍逻辑再回到真实数据上调。b的更新也值得单独验一次。当两个乘子都落在(0, C)内部时用b1或b2都行都落在边界上时取两者的平均。如果不论什么情况都只取b1在某些数据分布下偏置会偏掉表现出来就是决策边界整体平移准确率莫名低几个点——这种偏差不打印中间量根本看不出来。6. 二次开发把这个练习包改成能跑实验的骨架这份代码最大的价值不是能跑通而是它给了你每个算法最朴素的实现。想再往前走一步我建议做三件事成本都不高但能让它从练手脚本集合变成你自己的实验底座。第一件事是统一接口。现在每个文件的函数签名五花八门有的返回(w, losses)有的返回(labels, centers)。抽一个极简基类只约定fit和predict两个方法之后所有算法都往外层套一层import numpy as np class BaseModel: 所有手写模型的最小接口约定不依赖 sklearn 也能对拍 def fit(self, X, yNone): raise NotImplementedError def predict(self, X): raise NotImplementedError def score(self, X, y): return float(np.mean(self.predict(X) y)) class LogisticRegression(BaseModel): def __init__(self, lr0.1, epochs500): self.lr, self.epochs lr, epochs def fit(self, X, y): self.w_ np.zeros(X.shape[1] 1) Xb np.hstack([np.ones((len(X), 1)), X]) # 补一列 1 对应偏置 for _ in range(self.epochs): p 1 / (1 np.exp(-np.clip(Xb self.w_, -30, 30))) self.w_ - self.lr * Xb.T (p - y) / len(y) return self def predict(self, X): Xb np.hstack([np.ones((len(X), 1)), X]) return (Xb self.w_ 0).astype(int)这个改动的收益在验证环节立刻体现。有了统一接口你可以写一个run_experiment(models, X, y)把 k 折交叉验证、准确率、耗时一次性跑完所有算法共用同一套评估逻辑不会再出现这个文件用留出法、那个文件用全量评估导致的横向对比失真。第二件事是加对拍测试。这是我认为性价比最高的一步给自己写的每个模型配一个断言和 sklearn 的同名实现比准确率差太多就 fail。import numpy as np from sklearn.tree import DecisionTreeClassifier from sklearn.datasets import load_iris from sklearn.model_selection import train_test_split def test_tree_matches_sklearn(): X, y load_iris(return_X_yTrue) Xtr, Xte, ytr, yte train_test_split(X, y, test_size0.3, random_state0) mine build_tree(Xtr, ytr, max_depth4) # 这里需要写一个从嵌套字典做预测的辅助函数 tree_predict acc_mine np.mean(tree_predict(mine, Xte) yte) acc_ref DecisionTreeClassifier(max_depth4, random_state0).fit(Xtr, ytr).score(Xte, yte) # 手写版用的是信息增益ID3sklearn 默认是 gini允许 15 个点的差距 assert acc_mine acc_ref - 0.15, f手写 {acc_mine:.3f} vs sklearn {acc_ref:.3f}那个 15 个点的容差不是随便定的。ID3 用信息增益、CART 用基尼指数划分准则不同长出来的树本来就不一样要求数值完全相等是不现实的。但只要差距在这个量级以内就说明你的递归、终止条件、预测路径三处逻辑没有结构性错误。真正的 bug比如划分时用错了特征索引会让准确率直接掉到 0.33三分类瞎猜的水平一眼就能看出来。第三件事是把参数从硬编码里解放出来。现在的代码里C1.0、k3、lr0.01全都是写死的想换个值得改源码。抽成argparse或者一个 YAML 配置一次实验对应一行命令python run.py --algo svm --C 10 --kernel linear --data data/watermelon.csv python run.py --algo kmeans --k 4 --init kmeans --seed 0改完这三步你会得到一个有点粗糙但完全自主可控的实验框架。它的算法实现肯定不如成熟库快但每一步你都能打断点、能改公式、能打印中间矩阵——这种掌控感是调包调不出来的。说个我自己的教训。刚开始跑这些练习代码时我习惯直接改原文件做实验改完也不留备份结果有一次把一个for循环的范围从range(n)改成range(d)决策树准确率从 0.8 掉到 0.3我花了整个晚上怀疑是信息增益公式推错了最后发现是那个范围。从那以后我每次动这类学习型代码之前都强制先git init提交一次原始版本再开一个分支改跑出来的结果和原版对不上就直接git diff看差异。这个习惯救过我很多次也希望帮到你。本文还有配套的精品资源点击获取