资讯动态

手写机器学习算法:从NumPy推导到梯度更新的完整实现

发布时间:2026/10/14 11:59:01 来源:尧图企业网站定制
简介本资源是一套面向Python初学者与机器学习入门者的算法实践代码库聚焦于核心算法原理理解与动手实现适用于课程学习、课设开发及算法复现训练。压缩包共54个文件含31个Python源码如logistic_regression_train.py、svm.py、kmeans.py等、19个文本类数据与说明文件含heart_scale、svm_test_data等典型数据集以及test_data、data等结构化测试资源整体仅323KB轻量易解压运行。已有361人下载学习反映出其在算法教学与自学场景中的实用价值。读者可直接运行各章节独立模块如Chapter_1逻辑回归、Chapter_4 SVM、Chapter_10 KMeans等结合配套数据与简明文档快速验证算法效果目录按算法分章组织命名规范、模块解耦清晰便于逐个研读源码逻辑、调试参数或拓展实验是夯实机器学习工程实现能力的优质实践素材。1. 这不是“抄个sklearn就完事”的源码它要你亲手推导损失函数、手写梯度更新、在NumPy里把反向传播跑通“基于Python的机器学习算法设计源码”——这个标题在招聘JD、课程大纲和GitHub仓库名里高频出现但它常被误读为“用Python调包实现几个模型”。真相是真正值得称作‘算法设计源码’的必须脱离scikit-learn、PyTorch等高层封装在纯NumPy或少量SciPy中从数学定义出发逐行写出前向计算、损失构建、梯度推导与参数更新全过程。它解决的不是“怎么快速出结果”而是“为什么这个公式能收敛”“为什么学习率设0.01不设0.1”“为什么sigmoid在深层网络里会消失”。适合三类人刚学完《统计学习方法》想验证公式的研究生、面试前需手撕逻辑回归推导的算法岗求职者、以及带学生做课程设计的某高校导师——他们需要的不是黑匣子API而是一个可打断点、可改符号、可替换激活函数、可观察每轮权重变化的透明沙盒。这类源码的价值不在运行速度而在可控性、可解释性与教学穿透力你能看着w w - lr * (X.T (y_pred - y))这一行讲清它如何从最小二乘目标函数对w求偏导而来也能把sigmoid换成tanh后立刻在loss曲线里看到震荡加剧——这才是“设计”的起点。2. 从零构建线性回归不调用LinearRegression只靠NumPy完成推导闭环2.1 为什么必须手写——看清“最小二乘”不是魔法而是可微分的数学操作很多初学者认为线性回归就是model.fit(X, y)但真正的设计起点是目标函数$$J(w,b) \frac{1}{2m}\sum_{i1}^{m}(w^Tx^{(i)} b - y^{(i)})^2$$注意这里用了$\frac{1}{2m}$而非$\frac{1}{m}$是为了求导后消去平方项的系数2让梯度更简洁。这个细节在sklearn源码里被封装掉了但手写时你必须主动选择——它直接影响后续梯度表达式的干净程度。常见误区是直接套用矩阵形式$J(w)\frac{1}{2m}|Xw-y|^2$却忽略b的处理。正确做法是将偏置b并入权重向量令$\tilde{w} [w; b]$$\tilde{X} [X ; \mathbf{1}]$在X右侧拼一列全1向量则目标函数简化为$J(\tilde{w}) \frac{1}{2m}|\tilde{X}\tilde{w} - y|^2$。这一步看似简单却是后续所有梯度推导的基石。若跳过此步后续手动求$\frac{\partial J}{\partial b}$时极易出错。2.2 前向传播与损失计算用NumPy写得像数学公式一样直白import numpy as np def linear_forward(X, w): 前向传播X为(m,n)矩阵w为(n1,)向量含偏置 # X_aug: (m, n1)最后一列为1用于承载偏置项 X_aug np.hstack([X, np.ones((X.shape[0], 1))]) return X_aug w # 返回(m,)向量 def mse_loss(y_pred, y_true): 均方误差损失标量输出 m len(y_true) return np.sum((y_pred - y_true) ** 2) / (2 * m) # 示例数据生成模拟真实场景100个样本3个特征 np.random.seed(42) X np.random.randn(100, 3) y_true 2.5 * X[:, 0] - 1.3 * X[:, 1] 0.8 * X[:, 2] 4.2 np.random.randn(100) * 0.5 # 初始化权重含偏置共4维 w_init np.random.randn(4) * 0.01 y_pred linear_forward(X, w_init) loss mse_loss(y_pred, y_true) print(f初始损失: {loss:.4f}) # 输出类似初始损失: 28.6732关键说明linear_forward中np.hstack([X, np.ones(...)])是核心技巧——它把偏置b的线性组合b自然融入矩阵乘法避免单独加法操作。mse_loss分母用2*m而非m正是为了匹配目标函数定义确保后续梯度为$\frac{1}{m}\tilde{X}^T(\tilde{X}w - y)$无额外系数干扰。2.3 梯度推导与更新从公式到代码的逐项映射损失函数对$\tilde{w}$的梯度为$$\nabla_{\tilde{w}} J(\tilde{w}) \frac{1}{m}\tilde{X}^T(\tilde{X}\tilde{w} - y)$$这个公式必须手算一遍对$\tilde{w}_j$求偏导不能只背结果。代码实现需严格对应def compute_gradient(X, w, y_true): 计算梯度返回(n1,)形状的梯度向量 X_aug np.hstack([X, np.ones((X.shape[0], 1))]) y_pred X_aug w error y_pred - y_true # (m,)向量 grad (X_aug.T error) / len(y_true) # (n1, m) (m,) - (n1,) return grad def gradient_descent(X, y_true, w_init, lr0.01, epochs1000): 完整梯度下降训练循环 w w_init.copy() losses [] for epoch in range(epochs): grad compute_gradient(X, w, y_true) w w - lr * grad # 标准更新w : w - lr * ∇J(w) y_pred linear_forward(X, w) loss mse_loss(y_pred, y_true) losses.append(loss) if epoch % 200 0: print(fEpoch {epoch}, Loss: {loss:.6f}) return w, losses # 执行训练 w_trained, loss_history gradient_descent(X, y_true, w_init, lr0.01, epochs1000) print(f训练后权重: {w_trained}) # 输出类似训练后权重: [ 2.498 -1.299 0.799 4.198]参数说明lr0.01是经验起点过大导致loss震荡甚至发散见后文避坑epochs1000需配合loss曲线判断是否收敛w_init用小随机数如*0.01而非全零避免对称性导致梯度为零。注意compute_gradient中X_aug.T error是矩阵乘法不是点积——这是NumPy高效实现的关键也是公式到代码最易出错的映射点。3. 逻辑回归的完整设计Sigmoid、交叉熵与数值稳定性的硬核落地3.1 为什么不用MSE——交叉熵才是分类问题的理论最优解对二分类若仍用线性回归的MSE损失会出现两个致命问题一是预测值y_pred无界可为负或远大于1而真实标签y_true只有0/1导致loss在边界处梯度极小“梯度消失”二是MSE不是概率框架下的最大似然估计。正确目标是二元交叉熵Binary Cross-Entropy, BCE$$J(w) -\frac{1}{m}\sum_{i1}^{m}\left[y^{(i)}\log(\sigma(z^{(i)})) (1-y^{(i)})\log(1-\sigma(z^{(i)}))\right]$$其中$z^{(i)} w^Tx^{(i)}$$\sigma(z) \frac{1}{1e^{-z}}$。这个公式必须手写推导其梯度$$\nabla_w J(w) \frac{1}{m}X^T(\sigma(Xw) - y)$$惊人地简洁且梯度大小正比于预测误差$(\sigma(z)-y)$天然具备自适应学习能力。这正是逻辑回归优于线性回归分类的根本原因。3.2 Sigmoid的数值稳定性避免exp(-z)下溢与exp(z)上溢直接写1/(1np.exp(-z))在z很大如700时np.exp(-z)下溢为0结果为1z很小如-700时np.exp(-z)上溢为inf结果为nan。工业级实现必须用scipy.special.expit或手动裁剪def sigmoid_stable(z): 数值稳定的sigmoid处理z极大/极小情况 # 对z0计算 1/(1exp(-z)) # 对z0计算 exp(z)/(1exp(z))避免exp(-z)上溢 pos_mask (z 0) neg_mask ~pos_mask result np.empty_like(z, dtypefloat) result[pos_mask] 1 / (1 np.exp(-z[pos_mask])) result[neg_mask] np.exp(z[neg_mask]) / (1 np.exp(z[neg_mask])) return result # 验证稳定性 z_test np.array([700, -700, 0, 10]) print(sigmoid_stable(z_test)) # [1. 0. 0.5 0.9999546]提示sigmoid_stable中pos_mask/neg_mask的分离是关键。若统一用1/(1np.exp(-z))z-700时np.exp(700)直接溢出。此技巧在所有涉及sigmoid/tanh的源码中必须强制应用否则训练过程会在某次迭代突然崩坏。3.3 交叉熵损失与梯度的NumPy实现一行代码背后的数学def sigmoid_forward(X, w): 前向计算logits z Xw再经sigmoid X_aug np.hstack([X, np.ones((X.shape[0], 1))]) z X_aug w return sigmoid_stable(z) def bce_loss(y_pred_proba, y_true): 二元交叉熵损失y_pred_proba为sigmoid输出shape(m,) # 加极小值防止log(0) eps 1e-15 y_pred_clipped np.clip(y_pred_proba, eps, 1 - eps) return -np.mean( y_true * np.log(y_pred_clipped) (1 - y_true) * np.log(1 - y_pred_clipped) ) def bce_gradient(X, w, y_true): BCE梯度∇J (1/m) * X_aug.T (sigmoid(X_augw) - y_true) X_aug np.hstack([X, np.ones((X.shape[0], 1))]) y_pred sigmoid_forward(X, w) # 注意此处调用的是带stable的sigmoid error y_pred - y_true # (m,)即∂J/∂z grad (X_aug.T error) / len(y_true) # (n1, m) (m,) - (n1,) return grad # 生成二分类数据用线性可分模拟 np.random.seed(42) X_cls np.random.randn(200, 2) y_cls (X_cls[:, 0] X_cls[:, 1] 0).astype(int) # 标签0/1 # 初始化并训练 w_cls_init np.random.randn(3) * 0.01 w_cls_trained, _ gradient_descent( X_cls, y_cls, w_cls_init, lr0.1, # 分类问题常用更大lr epochs2000, loss_funcbce_loss, grad_funcbce_gradient )注意bce_loss中np.clip(y_pred_proba, eps, 1-eps)是防log(0)的必需操作gradient_descent函数需支持传入自定义loss_func和grad_func此处未展开但实际源码中必须模块化。lr0.1比回归任务大10倍因BCE梯度天然更平滑——这是通过亲手调试才能建立的直觉。4. 决策树ID3算法信息增益、递归分裂与离散化陷阱4.1 为什么ID3是理解树模型的“第一块砖”sklearn的DecisionTreeClassifier背后是CART分类回归树它使用基尼不纯度或MSE而ID3仅用信息增益Information Gain且要求输入特征为离散型。这看似局限实则是绝佳的教学入口它强迫你直面“如何量化一个特征对分类的帮助”——答案就是香农熵。ID3的每一步分裂都是在回答“按特征A的每个取值分组后组内标签的混乱程度比不分组时降低了多少”这个思想是所有树模型包括XGBoost、LightGBM的源头。4.2 熵与信息增益的手动计算用NumPy向量化实现给定标签数组yshape(m,)其熵定义为$$H(y) -\sum_{k1}^{K}p_k \log_2 p_k$$其中$p_k$是第k类占比。NumPy实现需避开循环def entropy(y): 计算标签y的香农熵 if len(y) 0: return 0 # 统计各类别频次 _, counts np.unique(y, return_countsTrue) probs counts / len(y) # 计算 -sum(p * log2(p))log2(0)定义为0 return -np.sum([p * np.log2(p) for p in probs if p 0]) def information_gain(X, y, feature_idx, thresholdNone): 计算按feature_idx特征分裂的信息增益 feature_vals X[:, feature_idx] if threshold is None: # 离散特征按每个唯一值分组 unique_vals np.unique(feature_vals) weighted_entropy 0 for val in unique_vals: mask (feature_vals val) subset_y y[mask] weight len(subset_y) / len(y) weighted_entropy weight * entropy(subset_y) else: # 连续特征按threshold二分ID3本不支持但实践中常需 mask_left (feature_vals threshold) subset_y_left y[mask_left] subset_y_right y[~mask_left] weight_left len(subset_y_left) / len(y) weight_right len(subset_y_right) / len(y) weighted_entropy ( weight_left * entropy(subset_y_left) weight_right * entropy(subset_y_right) ) return entropy(y) - weighted_entropy # 示例计算某特征的信息增益 X_demo np.array([[1, 2], [1, 3], [2, 2], [2, 3]]) y_demo np.array([0, 0, 1, 1]) print(f整体熵: {entropy(y_demo):.4f}) # 1.0000 print(f按特征0分裂增益: {information_gain(X_demo, y_demo, 0):.4f}) # 1.0000关键说明entropy中if p 0过滤掉零概率项避免log2(0)报错information_gain支持离散thresholdNone和连续指定threshold两种模式后者虽非ID3原生但工程中必备。注意np.unique返回的counts顺序与unique_vals严格对应这是向量化统计的基础。4.3 ID3递归建树终止条件与字典结构的巧妙表达ID3树本质是嵌套字典{feature: i, threshold: t, children: {val1: subtree1, val2: subtree2}}。递归终止条件有三1所有标签相同 → 返回该标签2无特征可分 → 返回众数标签3达到最大深度。手动实现强调可读性def build_id3_tree(X, y, features, max_depth5, depth0): 递归构建ID3树 # 终止条件1所有y相同 if len(np.unique(y)) 1: return {type: leaf, class: y[0]} # 终止条件2无特征或达最大深度 if len(features) 0 or depth max_depth: # 返回众数类 classes, counts np.unique(y, return_countsTrue) majority_class classes[np.argmax(counts)] return {type: leaf, class: majority_class} # 寻找最佳分裂特征遍历所有features best_gain -1 best_feature None best_threshold None for f_idx in features: # 离散特征遍历所有唯一值作为分割点ID3标准 unique_vals np.unique(X[:, f_idx]) for val in unique_vals: gain information_gain(X, y, f_idx, thresholdval) if gain best_gain: best_gain gain best_feature f_idx best_threshold val # 若无正增益退化为叶节点 if best_gain 0: classes, counts np.unique(y, return_countsTrue) majority_class classes[np.argmax(counts)] return {type: leaf, class: majority_class} # 递归构建子树 feature_vals X[:, best_feature] tree { type: split, feature: best_feature, threshold: best_threshold, children: {} } # 按best_threshold分割数据 mask (feature_vals best_threshold) X_left, y_left X[mask], y[mask] X_right, y_right X[~mask], y[~mask] # 递归构建左右子树 remaining_features [f for f in features if f ! best_feature] tree[children][left] build_id3_tree( X_left, y_left, remaining_features, max_depth, depth1 ) tree[children][right] build_id3_tree( X_right, y_right, remaining_features, max_depth, depth1 ) return tree # 构建示例树 tree build_id3_tree(X_demo, y_demo, features[0,1], max_depth3) print(tree) # 输出为嵌套字典清晰展示分裂逻辑提示build_id3_tree中remaining_features的剔除确保每个特征只用一次符合ID3特性mask (feature_vals best_threshold)是离散分裂的核心若处理连续特征此处应改为和。树结构用字典而非类便于print调试和json.dump保存——这是教学源码的实用主义选择。5. 避坑指南那些让算法在第37次迭代突然崩溃的隐藏雷区5.1 学习率设置不当不是越小越好也不是越大越快现象训练初期loss下降极慢1000轮后仍高于初始值或loss在几轮内骤降后剧烈震荡甚至NaN。原因学习率lr过小导致收敛缓慢过大则梯度更新跨过最优解在损失曲面“弹跳”。更隐蔽的是不同算法对lr敏感度差异巨大线性回归MSE梯度较平缓lr0.01安全逻辑回归BCE梯度在预测错误时较大lr0.1更合适而神经网络隐藏层梯度可能爆炸需lr0.001。解决初始化lr0.01观察前100轮loss曲线若单调下降可尝试×10若震荡×0.1实现学习率衰减lr_t lr_0 / (1 decay_rate * t)其中decay_rate0.01是常用起点绝对不要在训练中固定lr后不做监控——每次修改算法必须重调lr。5.2 特征未标准化梯度下降变成“瘸腿走路”现象某特征如房屋面积单位平方米数值在10^3量级另一特征如房间数在10^0量级训练loss下降异常缓慢或权重w中面积对应项极小如1e-5房间数对应项极大如1e2。原因梯度$\nabla_w J$中大尺度特征主导梯度方向小尺度特征更新几乎停滞。几何上损失曲面呈狭长椭圆GD需沿“之”字形路径爬行。解决from sklearn.preprocessing import StandardScaler scaler StandardScaler() X_scaled scaler.fit_transform(X) # 训练集上fit测试集上transform注意StandardScaler必须仅在训练集上fit测试集用同一scaler的transform否则数据泄露。手写源码中可手动计算X_mean X_train.mean(axis0),X_std X_train.std(axis0), 然后X_train_norm (X_train - X_mean) / X_std。5.3 Sigmoid/Tanh的梯度消失模型“失忆”的根源现象深层网络如3层以上全连接训练时前层权重几乎不更新loss下降停滞或训练后期准确率卡在50%附近不上升。原因Sigmoid导数$\sigma(z) \sigma(z)(1-\sigma(z))$最大值仅0.25当z远离0时如|z|3导数趋近于0。反向传播中梯度经多层连乘迅速衰减至1e-10以下无法驱动更新。解决换激活函数用ReLUmax(0,z)其导数在z0时为1彻底避免消失权重初始化用He初始化w ~ N(0, 2/n_in)替代随机小数使z初始分布集中在0附近批归一化BatchNorm在每层后加入BN层强制输入分布稳定——但这已超出纯NumPy范畴属进阶技巧。5.4 标签编码错误把多分类当二分类训现象训练时loss正常下降但预测准确率始终低于随机如3分类得33%或predict_proba输出所有类别概率接近相等。原因多分类任务中y_true应为整数标签如[0,1,2]或one-hot编码如[[1,0,0],[0,1,0]]但误用二分类的[0,1]标签训练softmax。此时模型被强制学习“非此即彼”丢失第三类信息。解决确认任务类型二分类用BCE损失sigmoid输出多分类用交叉熵损失softmax输出检查y维度y.shape应为(m,)整数标签或(m, k)one-hot手写softmax时务必用np.exp(z - np.max(z, axis1, keepdimsTrue))做数值稳定否则exp(1000)直接溢出。6. 验证你的源码是否真“可设计”三个不可妥协的检验动作6.1 梯度检验Gradient Checking用有限差分法给你的解析梯度“照X光”无论你推导的梯度公式多么优美都必须用数值方法验证。原理很简单对权重w_i施加微小扰动h计算损失变化率$(J(wh)-J(w-h))/(2h)$应与解析梯度grad[i]高度一致相对误差1e-7。这是所有可信算法源码的出厂必检项def gradient_check(X, y, w, loss_func, grad_func, h1e-5): 梯度检验对比解析梯度与数值梯度 grad_analytic grad_func(X, w, y) grad_numeric np.zeros_like(w) for i in range(len(w)): # 扰动w[i] w_plus w.copy() w_minus w.copy() w_plus[i] h w_minus[i] - h loss_plus loss_func(linear_forward(X, w_plus), y) loss_minus loss_func(linear_forward(X, w_minus), y) grad_numeric[i] (loss_plus - loss_minus) / (2 * h) # 计算相对误差 diff np.linalg.norm(grad_analytic - grad_numeric) norm np.linalg.norm(grad_analytic) np.linalg.norm(grad_numeric) rel_error diff / norm if norm 0 else 0 print(f梯度检验相对误差: {rel_error:.2e}) return rel_error 1e-7 # 对线性回归执行检验 is_correct gradient_check(X, y_true, w_init, mse_loss, compute_gradient) print(f梯度检验通过: {is_correct}) # 应输出True血泪经验我曾在一个自定义损失函数中漏掉一个负号梯度检验误差高达1e-1而模型在训练集上loss下降完美——这证明loss下降≠梯度正确。没有梯度检验的源码就像没装刹车的车表面能跑实则危险。6.2 参数敏感性分析画出lr、batch_size对收敛的影响热力图一个“可设计”的源码必须能快速回答“如果我把学习率从0.01改成0.05会怎样” 而不是重跑10次看运气。建议用itertools.product遍历参数组合记录每组的最终loss和收敛轮次import itertools import pandas as pd lrs [0.001, 0.01, 0.1, 1.0] batch_sizes [16, 32, 64, 128] results [] for lr, bs in itertools.product(lrs, batch_sizes): # 此处插入你的训练函数支持batch_size参数 final_loss, epochs_to_converge train_with_batch(X, y_true, lrlr, batch_sizebs) results.append({ lr: lr, batch_size: bs, final_loss: final_loss, epochs: epochs_to_converge }) df_results pd.DataFrame(results) # 用seaborn画热力图此处略去绘图代码 # 关键结论通常lr0.01 bs32是稳健起点lr0.1时bs需增大以稳定梯度提示热力图中若发现lr0.1在bs16时loss发散但在bs128时收敛良好说明小批量梯度噪声太大——这直接指导你后续实验应优先调bs而非盲目降lr。6.3 替换组件实验把Sigmoid换成Tanh观察决策边界如何变形“可设计”的终极体现是能像搭乐高一样替换核心组件。例如将逻辑回归的sigmoid_stable换成tanh并可视化决策边界变化def tanh_stable(z): 数值稳定的tanh return np.tanh(z) # numpy的tanh已内置稳定性处理 def tanh_forward(X, w): X_aug np.hstack([X, np.ones((X.shape[0], 1))]) z X_aug w return tanh_stable(z) # 训练tanh版逻辑回归损失函数需相应调整此处略 # 然后用meshgrid生成决策边界 xx, yy np.meshgrid(np.linspace(-3, 3, 100), np.linspace(-3, 3, 100)) grid_points np.c_[xx.ravel(), yy.ravel()] z_grid tanh_forward(grid_points, w_tanh_trained) pred_grid (z_grid 0).reshape(xx.shape) # tanh输出[-1,1]0判为正类 plt.contourf(xx, yy, pred_grid, alpha0.3, cmapRdYlBu) plt.scatter(X_cls[:, 0], X_cls[:, 1], cy_cls, cmapRdYlBu, edgecolorsk) plt.title(Tanh决策边界对比Sigmoid)玄学时刻你会发现tanh的边界比sigmoid更“锐利”在数据稀疏区延伸更远——这是因为tanh在z0处导数更大1 vs 0.25对边界附近的样本更敏感。这种直观对比是调包无法提供的认知深度。我坚持在每份新算法源码中执行这三项检验梯度检验是底线参数热力图是导航仪组件替换是显微镜。它们不增加最终模型性能但让每一次调试都成为对数学本质的确认。希望帮到你。本文还有配套的精品资源点击获取

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑