我从不建议第一个模型就直接上框架。框架封装得太干净反向传播在底层把梯度算好用户只看到一个.backward()调用真正的数学过程被藏得死死的。一旦遇到梯度爆炸、训练不收敛这类问题没有底层逻辑的人只能瞎猜学习率。这篇就把我用 NumPy 在 MNIST 上从零实现两层神经网络、并把反向传播的每一步数学推导和代码逐行对应的过程完整记录下来包括推导细节、实现代码和调试踩坑记录适合已经会基础 Python、想真正吃透反向传播原理的读者。1. 整体设计为什么必须手写一遍反向传播1.1 框架黑盒带来的认知断层用框架训练模型太容易了。定义网络结构、选损失函数、调用反向传播三行代码模型就跑起来了。但“能跑”和“理解”是两回事。我在带新人时发现一个典型现象问激活函数的导数是什么能答上来问这一层传到下一层的梯度矩阵到底是什么形状很多人就开始犹豫。这说明对反向传播缺少真正可操作的理解。手写 NumPy 实现的本质是让自己被迫面对张量的每一个维度变化。(64, 784)的权重矩阵和(784,)的输入向量怎么相乘结果传给谁反向的时候梯度怎么从 10 维的输出层流回 784 维的输入层——这些在框架里都被自动处理了手写的时候每一步都要自己算清楚。这种经验一旦建立后面学卷积、循环网络里的反向传播都会顺畅很多。1.2 网络结构选型MNIST 是 28×28 的灰度图每个像素值 0 到 255展平后就是 784 维的输入向量。手写数字有 10 类所以输出层是 10 维。中间的隐藏层我选了 64 个神经元激活函数用 ReLU。这个网络结构很小单张显卡都用不上普通 CPU 几十秒就能跑完一个 epoch特别适合反复实验验证梯度正确性。整个前向过程可以概括成两条线性变换、两次激活。第一层把 784 维的输入映射到 64 维的隐藏空间第二层把 64 维的隐藏表示映射到 10 维的类别得分。注意输出层用的是 softmax 激活把得分转换成概率分布。隐藏层用 ReLU 而不是 sigmoid主要原因是 ReLU 的导数在正区间恒为 1梯度可以顺利传导到第一层权重不会像 sigmoid 那样在两端饱和、梯度趋近于零。2. 反向传播的完整数学推导2.1 符号约定与前向过程后面所有公式都遵循这套记号先统一再说推导否则很容易看乱X输入矩阵形状(m, 784)m 是批大小W1第一层权重形状(784, 64)b1形状(64,)z1第一层线性输出z1 X W1 b1形状(m, 64)a1第一层激活输出a1 ReLU(z1)形状(m, 64)W2第二层权重形状(64, 10)b2形状(10,)z2第二层线性输出z2 a1 W2 b2形状(m, 10)a2softmax 输出形状(m, 10)每一行对应一个样本在 10 个类别上的概率Y真实标签形状(m,)取值 0 到 9Y_onehotone-hot 编码后的标签形状(m, 10)尖括号里的矩阵乘法我用np.dot实现统一处理。需要注意的是z1 X W1 b1这一步b1是向量加的时候利用 NumPy 广播机制加到每一行上。很多新手在这里困惑以为维度不匹配其实广播规则会帮我们把偏置自动扩展。2.2 损失函数的选取与化简分类任务最常用的是交叉熵损失。对单个样本损失定义为L_i -Σ_k y_k * log(a2_k)因为y是 one-hot 向量只有真实类别那一项是 1其余都是 0所以这个式子直接化简成L_i -log(a2_y)也就是说每个样本的损失只取决于模型给真实类别分配的概率。预测概率越接近 1损失越小越接近 0损失越大。从信息论角度看这就是衡量预测分布和真实分布之间的差异但实际计算时用化简后的形式就够了。整个批次的损失是所有样本损失的平均L (1/m) * Σ_i L_i后面梯度计算时记得除以 m否则梯度会随批大小变化不同批大小下的学习率就不具可比性。2.3 softmax 交叉熵的梯度推导这一节是整个反向传播里最核心的一步也是唯一需要一点数学技巧的地方。先看损失对z2的梯度。交叉熵损失对a2_y求导是-1 / a2_y。但a2是 softmax 的输出它依赖所有的z2分量不能直接对某一个z2_k求导。直接暴力推导会涉及对i j和i ! j两种情况的分类讨论比较麻烦。这里用一个小技巧。既然已经化简成L_i -z2_y log(Σ_k exp(z2_k))这是对数运算法则直接展开的结果可以逐项求导∂L_i / ∂z2_j -1_jy exp(z2_j) / (Σ_k exp(z2_k))右边第二项恰好就是 softmax 的第 j 个输出。所以得到非常简洁的形式∂L_i / ∂z2 softmax(z2) - y_onehot向量化后写成dz2 a2 - Y_onehot。对一批数据就是逐行相减。这也是为什么网上很多代码里反向传播的第一步直接就是dz2 a2 - y_onehot看似凭空出现其实背后是交叉熵和 softmax 结合后的化简结果。这一步单独记忆即可推导一次以后就不用重复了。2.4 梯度向隐藏层传播拿到输出层的梯度之后就可以沿着计算图往回传。第二层权重W2的梯度用链式法则拆成两部分损失对z2的梯度和z2对W2的梯度。∂L / ∂W2 (a1)^T · (a2 - Y_onehot) / m ∂L / ∂b2 mean(a2 - Y_onehot, axis0)这里注意维度的匹配(a1)^T形状是(64, m)乘上dz2形状(m, 10)结果就是(64, 10)正好和W2的形状一致。偏置的梯度更简单每个 batch 上对偏置求导就是对所有样本的梯度取平均所以用np.mean(dz2, axis0)。继续往前梯度要穿过 ReLU 激活层回到z1。ReLU 的导数是分段函数ReLU(z) 1 (当 z 0)否则 0所以∂L / ∂z1 ((a2 - Y_onehot) · (W2)^T) ⊙ ReLU(z1)∂L / ∂W1 (X)^T · ∂L/∂z1 / m ∂L / ∂b1 mean(∂L/∂z1, axis0)其中 ⊙ 表示逐元素相乘。为什么这里需要ReLU(z1)因为梯度穿过激活函数时只有对输入大于 0 的神经元梯度才能继续流动小于等于 0 的神经元相当于被“关闭”了。如果某一层的 ReLU 输出全部为零那这一层往后就没有梯度了这就是所谓“死亡 ReLU”。2.5 权重初始化的作用权重初始化看起来很不起眼但对训练影响非常大。我见过有人直接用np.random.randn(784, 64)初始化结果训了几个 epoch 损失都不降。原因在于如果初始化权重过大线性输出z1的方差会被放大很多倍ReLU 一下把所有负值置零梯度传不回去。使用 He 初始化可以缓解这个问题。它的原理是从输入维度出发让每层输出的方差保持稳定。ReLU 激活因为会把一半神经元置零输出方差减半所以 He 初始化用np.sqrt(2 / fan_in)作为标准差self.W1 np.random.randn(784, hidden_size) * np.sqrt(2.0 / 784) self.b1 np.zeros(hidden_size)输出层可以用np.sqrt(2.0 / hidden_size)类似的思路初始化。偏置直接用零向量即可因为后面梯度更新会自己调整。3. 代码实现细节一行一行对应公式3.1 数据加载与预处理MNIST 数据不一定要用框架自带接口。原始文件格式是 IDX 格式能帮我们熟悉 NumPy 的二进制读取能力。这里用一个极简的加载函数直接读训练集和测试集。import numpy as np def load_mnist_images(filename): with open(filename, rb) as f: data np.frombuffer(f.read(), dtypenp.uint8, offset16) return data.reshape(-1, 784) def load_mnist_labels(filename): with open(filename, rb) as f: data np.frombuffer(f.read(), dtypenp.uint8, offset8) return data这个代码有个细节很值得注意图像文件的头部占 16 字节包含魔数、图像数量、行数、列数各 4 字节标签文件头部占 8 字节包含魔数和标签数量。offset参数就是为了跳过这些头部信息。如果忘记 offset 直接读会把文件头当成数据导致维度完全错乱。拿到原始数据后有两个预处理步骤必不可少。第一是把像素值从 uint8 转成 float 并归一化到 [0, 1] 区间直接除以 255.0。这一步看起来简单但如果不做输入数据的量纲太大第一层加权求和结果会非常大没有归一化的数据会让梯度更新很不稳定。第二是把标签转成 one-hot 编码方便和 softmax 输出直接相减def one_hot(y, num_classes10): m y.shape[0] Y np.zeros((m, num_classes)) Y[np.arange(m), y] 1 return Y最常用的 trick 是用np.arange(m)生成行索引用原始标签生成列索引一行代码完成填充比循环快得多。3.2 两层网络的类定义与前向传播我把整个网络封装成一个类成员变量保存各层的参数和中间变量。之所以把中间变量z1、a1、z2都存下来是因为反向传播需要用到它们少了哪一个都没法求梯度。这是手写实现的惯用模式。class TwoLayerNet: def __init__(self, input_size, hidden_size, output_size, lr0.1): self.lr lr self.W1 np.random.randn(input_size, hidden_size) * np.sqrt(2.0 / input_size) self.b1 np.zeros(hidden_size) self.W2 np.random.randn(hidden_size, output_size) * np.sqrt(2.0 / hidden_size) self.b2 np.zeros(output_size) def forward(self, x): self.z1 np.dot(x, self.W1) self.b1 self.a1 np.maximum(0, self.z1) self.z2 np.dot(self.a1, self.W2) self.b2 exp_z2 np.exp(self.z2 - np.max(self.z2, axis1, keepdimsTrue)) self.a2 exp_z2 / np.sum(exp_z2, axis1, keepdimsTrue) return self.a2视角在原 calc_forward 里softmax 实现那部分有个数值稳定的关键操作np.exp(self.z2 - np.max(self.z2, axis1, keepdimsTrue))。如果不减最大值当 z2 中某些分量为较大正数时np.exp(z2)会出现溢出结果变成 inf除以 inf 就得到 NaN。减去最大值之后输入到 exp 的函数值最大为 0再大的 z2 也不会溢出。注意这里的np.max的轴方向。keepdimsTrue保留维度使得(m, 10)的矩阵减去(m, 1)的结果时正确广播逐行减去每行的最大值而不是减整个矩阵的最大值。3.3 反向传播与参数更新反向传播按推导向后封装即可核心代码很短但每一行的含义都要清楚。def backward(self, x, y_onehot): m x.shape[0] dz2 self.a2 - y_onehot dW2 np.dot(self.a1.T, dz2) / m db2 np.sum(dz2, axis0) / m drelu (self.z1 0).astype(np.float64) dz1 np.dot(dz2, self.W2.T) * drelu dW1 np.dot(x.T, dz1) / m db1 np.sum(dz1, axis0) / m return dW1, db1, dW2, db2 def update(self, dW1, db1, dW2, db2): self.W1 - self.lr * dW1 self.b1 - self.lr * db1 self.W2 - self.lr * dW2 self.b2 - self.lr * db2反向传播这块有几个容易写错的地方。首先是drelu用(self.z1 0)的布尔矩阵转成浮点型本质上就是一个掩码矩阵在反向传播中负责把梯度乘上 ReLU 的导数值。np.dot(dz2, self.W2.T)这一步是把输出层的梯度“从后往前投影”到隐藏层空间再逐元素乘上 ReLU 掩码得到dz1。最后dW1用x.T乘dz1再除以 mm 是批大小。除以 m 是和前向损失取平均值相对应的。如果这一步忘了除以 m在不同 batch size 下梯度量级完全不同哪怕模型结构相同也需要重新调学习率。参数更新用的是最简单的小批量随机梯度下降。为了避免在一条线上来回走位后面可以考虑加动量我会在第五节里展开讨论。3.4 训练循环与评估训练循环分成几个常规步骤前向算出预测、算损失、反向求梯度、更新参数。我额外增加了每批记录一次损失和累计正确数这样训练结束后可以直接观察到损失下降曲线和准确率变化。def train(self, X, Y, epochs, batch_size64): m X.shape[0] indices np.arange(m) for epoch in range(epochs): np.random.shuffle(indices) total_loss 0.0 correct 0 for i in range(0, m, batch_size): batch_idx indices[i:ibatch_size] x_batch X[batch_idx] y_batch Y[batch_idx] y_onehot one_hot(y_batch, 10) pred self.forward(x_batch) loss self.loss(pred, y_batch) total_loss loss * x_batch.shape[0] dW1, db1, dW2, db2 self.backward(x_batch, y_onehot) self.update(dW1, db1, dW2, db2) correct np.sum(np.argmax(pred, axis1) y_batch) train_acc correct / m test_acc self.evaluate(X_test, Y_test) print(fEpoch {epoch1}/{epochs}, loss{total_loss/m:.4f}, train_acc{train_acc:.4f}, test_acc{test_acc:.4f})这里total_loss的累加考虑了批大小用x_batch.shape[0]乘以每个样本的平均损失把不同批次的损失加总最后再除以总样本数得到整个 epoch 的平均损失。每次都去求 loss 其实有点浪费但为方便观察还是保留。测试评估函数类似前向不需要反向、不需要存中间变量def evaluate(self, X, Y): pred self.forward(X) return np.mean(np.argmax(pred, axis1) Y)后续可以加一个简单的 mini-batch 构造器方便在同一数据集上反复打乱迭代。上面这个循环已经把 shuffle 写进了 epoch 里用数组索引取子集不复制大数组内存占用很低。4. 常见问题与排查技巧实录4.1 先用数值梯度验证反向传播正确性手写反向传播最常见的 bug 是维度对不上、或者某个矩阵转置写错导致训练时损失看似下降但准确率奇低。一个非常有效的方法是做梯度检查用数值方式近似计算损失对某个参数的梯度和解析梯度对比。数值梯度的原理很简单参数 $\theta$ 的梯度可以用中心差分近似dL/dθ ≈ (L(θ h) - L(θ - h)) / (2h)取h 1e-5即可。实现时可以这样写def numerical_gradient(net, X, Y, param_name, epsilon1e-5): param getattr(net, param_name) grad np.zeros_like(param) it np.nditer(param, flags[multi_index]) while not it.finished: idx it.multi_index old param[idx] param[idx] old epsilon loss_plus net.compute_loss(X, Y) param[idx] old - epsilon loss_minus net.compute_loss(X, Y) param[idx] old grad[idx] (loss_plus - loss_minus) / (2 * epsilon) it.iternext() return grad验证的时候取一个小数据集比如 32 个样本分别计算解析梯度和数值梯度计算相对误差relative_error |grad_analytic - grad_numeric| / (|grad_analytic| |grad_numeric|)如果误差在1e-7到1e-5量级基本说明反向传播是对的。我每次搭新网络都会跑这个检查因为它可以一次性把所有权重矩阵的形状、转置、除法写错的问题抓出来。要注意的是数值梯度每更新一个参数就要重新走一遍前向矩阵大的时候非常慢所以只在小模型上做一次验证就够了不用每次训练都带。4.2 学习率与收敛速度的关系学习率是整个训练过程里最需要调的参数。拿我初始用的lr0.1在二层网络加 ReLU 的条件下loss 下降是比较平稳的。试过把学习率调到1.0结果 loss 在几百次迭代内变成inf因为更新步长太大一步跨过了最优区域在极值点附近来回震荡甚至发散。将学习率调小到0.01之后能收敛但速度明显变慢第一个 epoch 结束训练准确率只有 60% 上下用0.1能到 90% 左右。这说明学习率不能拍脑袋定可以从0.1开始如果 loss 振荡就降低一个量级再试如果 loss 下降太慢就适当增大。实操时还可以做学习率衰减每几个 epoch 把学习率乘一个小于 1 的系数让训练初期快速收敛、后期精细微调。4.3 softmax 的数值稳定性如果直接写np.exp(self.z2) / np.sum(np.exp(self.z2), axis1, keepdimsTrue)训练中一旦z2里出现较大的正值比如 100exp(100)会溢出为inf。一旦出现 inf交叉熵算出来是 NaN后面梯度也跟着变成 NaN整个训练就废了。解决办法就是减掉每行的最大值前面代码已经写了。这一步不会改变 softmax 的数学结果因为分子分母同时乘以exp(-max)概率不变但数值上稳定得多。这是一个非常简单但能救命的处理框架里已经内置了手写时一定要记得。4.4 训练准确率高但测试准确率低隐藏层设成 128 的时候训练集准确率能到 99% 以上但测试集只有 96% 左右。这是典型的过拟合表现。MNIST 本身比较简单过拟合空间不大但这个现象值得注意因为原理和更复杂的任务是一致的。缓解过拟合有几个立即可做的操作第一是减小模型容量比如隐藏层从 128 降到 32观察测试准确率是否反弹第二是加正则化最简单的 L2 正则化就是给每个权重矩阵的梯度加上一个lambda * W的惩罚项第三是加 dropout在前向和反向传播时随机屏蔽掉一部分隐藏层神经元让模型不敢过度依赖某一个特征。用手写实现 dropout 也不难但要注意训练时对激活输出除以保留概率做缩放推理时不缩放。4.5 输出全为零或者梯度消失的排查ReLU 网络一个常见副作用是“死亡 ReLU”就是某层大量神经元输出恒为 0反向传播时梯度被掩码屏蔽权重永远得不到更新。判断方法很简单训练几轮后打印隐藏层输出的均值如果一组数据全部是 0就要检查初始化方法是不是有问题。用 He 初始化能大幅减少这种情况。如果用np.random.randn裸初始化标准差是 1线性输出 z 的方差大约等于输入维度ReLU 后负值全部变零正值的方差也被放大非常容易造成神经元饱和或死亡。所以初始化权重不是“运气问题”而是有明确数学依据的。4.6 常见问题速查表现象可能原因排查方法loss 为 NaN 或 infsoftmax 数值溢出检查是否做了 max 减法稳定处理loss 震荡严重学习率过大降低学习率一个量级训练不下降权重初始化不当换用 He 初始化训练/测试准确率差距大过拟合加 dropout 或 L2 正则化隐藏层输出全是 0ReLU 死亡调初始化、适当调低学习率梯度检查误差大反向传播实现错误核对每个矩阵转置和除法5. 扩展优化从能跑到跑得好5.1 动量与小批量梯度下降纯 SGD 的问题在于每一步只看了当前批次更新方向噪声很大loss 曲线容易抖动。最简单的改进是加动量项维护一个历史梯度移动平均每次更新时不直接用当前梯度而是把历史方向和当前方向叠加。def update_with_momentum(self, dW1, db1, dW2, db2, gamma0.9): self.vW1 gamma * self.vW1 self.lr * dW1 self.W1 - self.vW1 self.vW2 gamma * self.vW2 self.lr * dW2 self.W2 - self.vW2初次使用vW1时要在__init__里初始化为零矩阵。动量项的直观意思像小球滚下坡当前梯度方向一致时加速方向相反时减速所以能较快通过平坦区域同时抑制震荡。实测加了这个之后同样 10 个 epoch准确率能稳定提升零点几个百分点而且 loss 曲线更平滑。5.2 更复杂的优化器如果还想再进一步可以在 NumPy 里手写 Adam。Adam 融合了动量项和自适应学习率对每个参数单独调整学习率分量数学一点都不复杂m_t β1 * m_{t-1} (1 - β1) * g_t v_t β2 * v_{t-1} (1 - β2) * g_t^2 m_hat m_t / (1 - β1^t) v_hat v_t / (1 - β2^t) param - lr * m_hat / (sqrt(v_hat) eps)其中 β1 常用 0.9β2 常用 0.999eps 取1e-8。手写 Adam 也不难但第一次做网络时我建议先把 SGD 玩透再上 Adam。很多初学者直接 Adam 一把梭反而学不到调参的手感。5.3 训练过程可视化手写实现的一大价值就是能看到每一步的内部变化。我一般每层都存下激活值的统计均值、方差、极端值打印出来观察。比如发现激活均值偏向 0说明 ReLU 很多输出是 0隐藏层表达能力变弱就要注意初始化或者学习率是否过大。下面这段代码可以打印每一层激活的统计信息我一般每隔几个 batch 执行一次用来感知训练过程是否正常def inspect_activation(self, x): self.forward(x) print(fz1 mean{self.z1.mean():.4f}, std{self.z1.std():.4f}) print(fa1 mean{self.a1.mean():.4f}, zero_ratio{(self.z1 0).mean():.4f}) print(fa2 max{self.a2.max():.4f}, min{self.a2.min():.4f})zero_ratio这个指标非常有用。如果第一轮训练中隐藏层神经元死亡比例超过 0.5说明初始化可能有问题梯度流动不畅模型很难学到特征。5.4 从两层到多层的自然演进这篇基于两层网络讲透了反向传播扩展到三层或更深网络时其实只需要在反向传播代码里多写一组相同模式的层。具体做法是前向传播时依次计算z - a - z - a把每一层的中间结果存进一个列表反向传播时从最后一层的dz a2 - y开始一层一层往前乘以W.T和激活函数的导数。你可以把核心反向代码抽象成一个循环dZ self.a_layers[-1] - y_onehot for i in reversed(range(num_layers)): dW np.dot(self.activations[i].T, dZ) / m db np.sum(dZ, axis0) / m dZ np.dot(dZ, self.W_layers[i].T) * (self.z_layers[i-1] 0)不管网络多深这两件事永远不变先求损失对当前层线性输出的梯度再把它传回上一层的激活输入。这个规律弄明白了学卷积网络、循环网络都是同一个思路只是每个模块的局部梯度不同。6. 最后的几点体会把整个项目跑通的完整过程我觉得比跑通框架里一条命令的收获大得多。数值梯度检查一定要做它能在五分钟内定位反向传播所有形状错误省下后面无数排查时间批次归一化虽然是现代网络的标配但第一步不要加先把裸 BP 的数学机制彻彻底底跑一遍再逐步加技术手段每一步改进才有的放矢。最后分享一个我训练时的小技巧不要只看最终准确率要在中间打印每一轮训练集和测试集准确率。如果你发现训练集涨、测试集停滞过拟合的信号就出来了如果两边都不涨更可能是学习率或初始化问题。反向传播就是这么个东西数学上一次性推导清楚代码上反复验证它就从“高深名词”变成一个可靠的工程工具。之后再去学任何现代框架你会知道每一个.backward()背后到底在算什么。
