有过一段日子，我卡在一道经典机器学习公开课的编程练习上：公式不是看不懂，可一旦落到 Python 代码里，输出总是不对。这个练习要求徒手实现一个三层神经网络，核心就是反向传播算法。前向传播还能照着公式写，真正折磨人的是误差从输出层一层层传回去的过程，再加上偏置、正则、向量化维度这些问题，很多人在这一步直接心态爆炸。如果你正在学机器学习，或者刚好做到手写数字识别、多层感知机相关的练习，那这篇应该能帮上忙。我会从题目拆解讲起，说清楚反向传播每一步到底在算什么，然后给出完整可运行的 Python 代码，最后把常见坑位和调试技巧一并列出来。

1. 任务拆解：课后习题到底在问什么

1.1 核心需求：前向传播、反向传播、代价函数一个都不能少

经典课后习题4的题面，一般会给你一批手写数字图片。每张图片是 20 乘 20 的灰度图，展开成 400 个像素。目标是用神经网络把这些图片分成 10 类，也就是数字 0 到 9。

网络结构是固定的三层：输入层 400 个神经元，隐藏层 25 个神经元，输出层 10 个神经元。题面明确要求实现三件事：第一，带正则化的代价函数；第二，前向传播，计算每一层的激活值；第三，反向传播，求出代价函数对两个参数矩阵 Theta1、Theta2 的偏导数。最后还需要用梯度下降把参数迭代出来。

这里有个很常见的误解：题目虽然把“梯度下降”也放进练习里，但检查重点其实是“算出来的梯度到底正不正确”。这个梯度不是线性回归那种能直接写解析式的，它必须靠反向传播一点一点算出来。所以我在写代码时，强烈建议把“算梯度”和“用梯度更新参数”拆成两个独立环节。这样万一准确率上不去，你可以先用梯度校验定位问题，到底出在梯度计算还是出在优化策略。

1.2 先画数据流图，把维度写在每条边上

反向传播最烦人的不是数学，而是维度。有一个方法能极大降低出错率：在写任何代码之前，先在纸上画出数据流，每一步都标好张量的形状。不要嫌麻烦，这比报错后瞎猜高效得多。

具体到这个题：输入矩阵 X 是 (5000, 400)，5000 是样本数，400 是特征数。Theta1 设计成 (25, 401)，多出来的一列对应隐藏层偏置。前向传播的第一步是在 X 左边拼一列 1，变成 (5000, 401)，再与 Theta1 的转置做矩阵乘法，得到 z2，维度是 (5000, 25)。经过 sigmoid 激活得到 a2，维度不变。接着再给 a2 拼一列 1，变成 (5000, 26)，再与 Theta2 的转置乘，Theta2 是 (10, 26)，得到 z3 和 a3，维度是 (5000, 10)。

画完图就能发现，所有矩阵乘法的维度都严格对得上。后面写反向传播时，再判断谁该转置，直接看目标结果的维度就行。比如你想算 Delta2，它的形状必须和 Theta2 一样是 (10, 26)，已知 delta3 是 (5000, 10)，a2 是 (5000, 26)，唯一对得上的乘法就是 delta3.T.dot(a2)。

另一个容易被忽略的点是标签编码。输出层有 10 个神经元，所以不能用 y=3 这种整数直接和 a3 比较，得先把标签转成 one-hot 形式，也就是一个 10 维向量，对应位置是 1，其它是 0。我第一次做的时候就因为没转，导致梯度完全算错，损失看起来在下降，准确率却一直卡在 10% 左右。

2. 反向传播原理与公式推导

2.1 符号约定：从向量化角度重新看网络

为了不让公式和代码各说各话，先把符号统一起来。用 m 表示样本数，n_in 表示输入维度，n_hidden 表示隐藏层神经元数，n_out 表示输出类别数。Xi 是第 i 个样本，但在向量化实现里，我们一次性输入整个 X 矩阵。

Theta1 是隐藏层权重，维度是 (n_hidden, n_in + 1)；Theta2 是输出层权重，维度是 (n_out, n_hidden + 1)。注意两个矩阵的行数都等于当前层神经元数，列数等于上一层输入数加一。这个约定和很多教材里“Theta 的每一行是某个神经元所有输入边的权重”是一致的。

关于激活函数，这个题用的不是 ReLU，而是 sigmoid。隐藏层、输出层都用 sigmoid。输出层不用 softmax，这一点让很多人困惑。我理解课程选择 sigmoid，是为了让输出层的每个神经元可以看作一个独立的逻辑回归模型，这样在推导交叉熵代价对 sigmoid 的梯度时，二项交叉熵的形式可以直接套用。虽然 softmax 更常用于多分类，但在这里两者的思想是一致的，只是误差项表达稍有不同。

如果换成 softmax，反向传播的 delta3 会变成 a3 - Y，和 sigmoid 加交叉熵的结果一模一样。所以用这个练习理解多分类神经网络，完全不会走偏。

2.2 三条核心公式，从 delta3 到梯度更新

反向传播的基础公式其实只有几个。第一，输出层误差 delta3。在代价函数为交叉熵、激活函数为 sigmoid 的设定下，delta3 = a3 - Y，Y 是 one-hot 标签矩阵。为什么不是 a3 乘以 (1-a3) 之类的东西？因为交叉熵的导数恰好把 sigmoid 的导数抵消了。也就是说，这个简洁形式是有条件的，换一个代价函数就不成立。

第二，隐藏层误差 delta2。delta2 = (delta3 乘以 Theta2 去掉偏置列) 然后逐元素乘以 sigmoid 的导数。具体写出来就是 delta3.dot(theta2[:, 1:]) * a2 * (1 - a2)。这里逐元素相乘不能写成矩阵乘法，因为每个样本在每个神经元上的误差是独立的。

第三，梯度累加。Delta1 = delta2.T.dot(a1) / m，Delta2 = delta3.T.dot(a2) / m。如果要带正则化，还要在对应矩阵除第一列外加上 lambda/m 倍的权重。我的做法是先把不包含正则的梯度算出来，再单独加正则项，这样梯度校验时能分开排查。

参数更新就用最朴素的批量梯度下降：Theta1 -= alpha * Delta1，Theta2 -= alpha * Delta2。这个题没有要求用 Adam，纯梯度下降就能收敛到不错的效果。

2.3 理解“反向”到底反的是什么

很多人把反向传播当复杂黑盒，其实它的核心就是链式法则。前向传播像河流从上游往下游，代价函数是下游的湖泊；反向传播则从湖泊出发，逆流而上，把每个闸门受到的影响算清楚。

具体来说，代价函数对某一层权重 Theta 的偏导数，可以拆成“误差项”和“该层输入激活值”的乘积。误差项是从输出层一层层往前传的，它记录的是“如果这个神经元输出变化一点点，最终代价会变化多少”。有了这些误差项，计算梯度的效率就很高，不需要对每个权重单独扰动再跑一次前向传播，一次前向加一次反向就能拿到所有梯度。

这也是为什么讲神经网络时总说“反向传播是计算梯度的引擎”。它本身不负责优化，优化是梯度下降的事情。

3. Python 实现的关键步骤

3.1 数据准备：归一化、one-hot、维度检查

拿到数据第一步是检查范围。如果像素值在 0 到 255 之间，需要除以 255，让特征落入 [0, 1]。不做这步，输入到 sigmoid 的加权和很容易变得很大，激活值饱和，梯度接近 0。如果数据已经是 [0, 1]，就跳过。

第二步是标签转换。我一般用 np.eye(10)[y]，y 是形状为 (m,) 的整数数组。这样得到的 Y 是 (m, 10) 的 one-hot 矩阵。如果标签不是从 0 开始，需要先减掉最小值，确保类别索引正确。这里有个常见坑：数据里标签可能是 1 到 10，而 np.eye(10) 只支持 0 到 9，直接用会越界。所以转换前务必确认标签范围。

第三步是维度检查。X.shape[1] 应该等于 400，Y.shape[1] 应该等于 10。任何不一致都说明数据读取有误。我习惯在一开始就加 assert，让错误尽早暴露。

3.2 参数初始化：随机但不离谱

初始化对神经网络影响很大。如果 Theta1 和 Theta2 全零，同一层的每个神经元会收到完全相同的梯度，无论训练多少轮，它们都像一个人，网络实际上退化成一个很窄的模型。如果初始值太大，激活值一开始就饱和。所以要用小随机数。

这个练习给出的 epsilon 范围大约是 0.12，我会直接用 np.random.uniform(-0.12, 0.12, size=(25, 401))。在实际项目中，可以用 Xavier 规则：epsilon = sqrt(6) / sqrt(n_in + n_hidden)。算出来大概也是 0.11，和 0.12 很接近。

随机种子建议固定。不是所有随机初始化效果都一样好，固定种子才能复现实验结果，排查问题时也才能保证每次对比都在同一条件下。

3.3 前向传播：缓存每一层的 z 和 a

前向传播不能只返回最终预测，还必须把中间的 z 和 a 缓存下来。因为反向传播要用到它们。我习惯在函数里计算完每一层后，存到列表或字典里，最后一次性返回。

具体流程是：a1 = np.c_[np.ones((m, 1)), X]；z2 = a1.dot(theta1.T)；a2 = sigmoid(z2)；然后 a2 = np.c_[np.ones((m, 1)), a2]；z3 = a2.dot(theta2.T)；a3 = sigmoid(z3)。注意，拼接偏置这个操作要重复两次，第二次很容易漏掉。我的记忆方法是：每进入一个新的线性层之前，都要给上一层的激活值加上偏置列。输入层要加，隐藏层激活后也要加，输出层之后就不用加了。

这里选择用 np.c_ 来拼列，避免 reshape 带来的麻烦。如果你习惯 np.hstack，效果一样。

3.4 反向传播：四个关键操作背后的直觉

反向传播函数接受前向传播缓存下来的 a 和 z，以及 one-hot 标签，返回 D1、D2。第一步计算 delta3 = a3 - Y。这是输出层误差。

第二步计算 delta2 = delta3.dot(theta2[:, 1:]) * sigmoid_gradient(z2)。这里的 theta2[:, 1:] 表示去掉每一行中的偏置列，因为偏置对应的权重不参与上一层误差传递。很多初学的人会问，为什么这里不去掉 delta3 中的偏置列？因为偏置是加在输入侧的，输出侧没有偏置。去掉 theta2 的偏置列后，得到的线性组合乘以 sigmoid 导数，就是隐藏层误差。

第三步计算 D1 = delta2.T.dot(a1) / m。这里 a1 是包含偏置的。D2 = delta3.T.dot(a2) / m，这里 a2 也是包含偏置的。这一步的原则是：反传时去掉权重偏置列，算梯度时保留激活偏置列。

第四步，如果正则化系数 lambda 大于 0，就在 D1 和 D2 的除第 0 列以外的列上加上 lambda/m 乘以对应的权重。这相当于把权重衰减梯度加到普通梯度里。注意不能加在偏置对应的第 0 列上。

3.5 梯度校验：用数值梯度当照妖镜

反向传播代码写完，不应该直接拿去训练。第一件事是梯度校验。数值梯度的思想是用导数的定义近似：J(theta + epsilon) - J(theta - epsilon) 除以 2epsilon。如果反向传播算出来的梯度和数值梯度接近，说明公式和实现没问题。

我建议随机抽取几个参数位置做校验，用同一份小样本数据，lambda 可以先设为 0 或 0.1。对比时看相对误差：|numerical - analytic| / max(1, |numerical|, |analytic|)。小于 1e-4 基本合格，大于 1e-2 必然有 bug。

实现细节是这样的：要复制 theta 再加或减 epsilon，不要在原数组上直接改。而且要保证 epsilon 比数值精度足够大，我常用 1e-5。梯度校验跑起来很慢，所以只针对少量位置，不要全参数扫描。

3.6 训练循环与收敛判断

训练循环其实很简单：用当前 theta 算代价 J，用 backprop 算梯度，更新参数，记录历史损失。我在每 10 轮打印一次 loss，方便观察。刚开始跑的时候，loss 会有一个明显下降，之后逐渐平缓。

关于收敛判断，应该先看训练集准确率，而不是死盯 loss。如果 loss 降到 0.2 左右，但准确率只有 80%，很可能是隐藏层规模不够，或者数据本身有噪声。在这个练习中，固定结构下 50 到 100 轮基本能到 95% 以上。如果 30 轮准确率还低于 70%，通常不是轮数问题，而是前面哪个环节错了。不要盲目加 epochs。

4. 踩坑与排查实录

4.1 常见问题速查表

下面是我实际调试中遇到比较多的情况，整理成表方便你对照。

4.2 梯度消失与学习率的选择

这个练习用的是 sigmoid，sigmoid 在两端饱和。如果 z 的绝对值大于 4，导数已经非常接近 0。你可以打印 z2 和 z3 的最大最小值，如果都很大，说明初始化范围太大或特征没缩放。这时梯度传到前面层会变得很小，前面权重几乎不更新。

学习率方面，我通常从 alpha=1.0 开始试，跑 10 轮如果 loss 单调下降，就继续用；如果震荡，就降低到 0.3；如果下降太慢，可以考虑加大到 3。梯度下降对学习率很敏感，但在这个小规模任务里，0.3 到 1 之间通常能找到稳定解。

4.3 初始化范围与随机种子敏感度

即使代码完全正确，随机初始化不同也会导致最终准确率在小范围内波动。这很正常。但如果某一次跑到 90%，另一次只有 80%，说明初始化或者学习率对网络影响太大。一种做法是多次运行取平均，另一种是用固定种子在调参阶段保证对比公平。

我见过有人为了追求稳定，把所有权重视为相同小值。这样虽然梯度不为零，但每个隐藏神经元学到的特征仍然一样，网络容量被浪费了。所以随机初始化必须保留，固定种子就够用。

4.4 独门调试技巧：四个立竿见影的办法

第一个办法：小数据过拟合测试。选 100 条样本，把 lambda 设为 0，训练 20 轮。如果梯度计算正确，训练集准确率应该接近 100%。如果达不到，问题基本在实现，不在数据。

第二个办法：线性网络降级测试。把 sigmoid 替换成线性激活，网络就退化成多层线性回归。如果梯度实现里有概念性错误，这种退化结构会很容易暴露。

第三个办法：梯度校验分阶段。第一次只校验不带正则的 D1，第二次校验带正则的 D1，第三次再校验 D2。不要一次性都打开，否则很难定位。

第四个办法：冻结前层只调后层。把 Theta1 固定，只用反向传播更新 Theta2。如果这个子模型能收敛，说明后层实现没问题，问题多半出在前面。

5. 完整代码：这份实现可以直接运行

5.1 代码整体结构说明

下面这份代码不依赖任何深度学习框架，只用 numpy。我把整个流程拆成多个小函数，方便你修改和调试。每个函数只负责一件事：数据加载、前向传播、代价函数、反向传播、梯度校验、训练、预测。

load_data 目前用随机数占位，目的是让代码可以立刻跑通。你把它替换成自己的数据读取逻辑后，才能看到真实分类效果。替换方法我会写在后面。

5.2 完整代码块

import numpy as np def load_data(): # 这里放你的数据读取逻辑 # 返回 X: (m, 400), y: (m,) 整数0-9 # 先用随机数据模拟，保证流程能跑通 m = 5000 X = np.random.rand(m, 400) y = np.random.randint(0, 10, m) return X, y def sigmoid(z): return 1.0 / (1.0 + np.exp(-z)) def sigmoid_gradient(z): a = sigmoid(z) return a * (1 - a) def forward(X, theta1, theta2): m = X.shape[0] a1 = np.c_[np.ones((m, 1)), X] # (m, 401) z2 = a1.dot(theta1.T) # (m, 25) a2 = sigmoid(z2) a2 = np.c_[np.ones((m, 1)), a2] # (m, 26) z3 = a2.dot(theta2.T) # (m, 10) a3 = sigmoid(z3) return a1, z2, a2, z3, a3 def cost(theta1, theta2, X, y, lam): m = X.shape[0] _, _, _, _, h = forward(X, theta1, theta2) Y = np.eye(10)[y] J = -np.mean(Y * np.log(h + 1e-10) + (1 - Y) * np.log(1 - h + 1e-10)) reg = (lam / (2 * m)) * (np.sum(theta1[:, 1:] ** 2) + np.sum(theta2[:, 1:] ** 2)) return J + reg def backprop(X, y, theta1, theta2, lam): m = X.shape[0] a1, z2, a2, z3, a3 = forward(X, theta1, theta2) Y = np.eye(10)[y] delta3 = a3 - Y # (m, 10) delta2 = (delta3.dot(theta2[:, 1:])) * sigmoid_gradient(z2) D1 = delta2.T.dot(a1) / m D2 = delta3.T.dot(a2) / m D1[:, 1:] += (lam / m) * theta1[:, 1:] D2[:, 1:] += (lam / m) * theta2[:, 1:] return D1, D2 def numerical_grad_check(X, y, theta1, theta2, lam, idx): eps = 1e-5 theta_plus = theta1.copy() theta_minus = theta1.copy() i, j = idx theta_plus[i, j] += eps theta_minus[i, j] -= eps num = (cost(theta_plus, theta2, X, y, lam) - cost(theta_minus, theta2, X, y, lam)) / (2 * eps) return num def train(X, y, n_hidden=25, alpha=0.5, lam=0.1, epochs=60, seed=42): rng = np.random.default_rng(seed) n_in = X.shape[1] n_out = 10 eps_init = 0.12 theta1 = rng.uniform(-eps_init, eps_init, (n_hidden, n_in + 1)) theta2 = rng.uniform(-eps_init, eps_init, (n_out, n_hidden + 1)) history = [] for epoch in range(epochs): J = cost(theta1, theta2, X, y, lam) D1, D2 = backprop(X, y, theta1, theta2, lam) theta1 -= alpha * D1 theta2 -= alpha * D2 history.append(J) if (epoch + 1) % 10 == 0: print(f"Epoch {epoch+1:3d}, cost = {J:.6f}") return theta1, theta2, history def predict(X, theta1, theta2): _, _, _, _, h = forward(X, theta1, theta2) return np.argmax(h, axis=1) def accuracy(X, y, theta1, theta2): p = predict(X, theta1, theta2) return np.mean(p == y) if __name__ == "__main__": X, y = load_data() theta1, theta2, history = train(X, y) print("Accuracy:", accuracy(X, y, theta1, theta2))

5.3 运行结果与数据替换说明

用随机数据训练不会得到有意义准确率，但代码流程是完整的。load_data 设计成可替换接口，你把自己的数据读进来后，只需要把 X、y 赋值给它。常见读取方式有 np.loadtxt 或 pandas.read_csv。读进去之后确认 X 是二维数组 (m, n)，y 是一维整数数组。

如果你使用课程附带的 5000 张手写数字子集，这个网络结构在几十轮迭代后，训练集准确率一般能到 95% 以上，部分初始化下能到 98%。100 轮之后基本稳定。建议你把准确率作为验收标准：低于 95%，就要回去查梯度。

如果想换隐藏层大小，把 train 函数里的 n_hidden 改成其他值就行。但注意，隐藏层大小改变后，初始化范围、迭代次数都可能要重新调，不是越大越好。过大的隐藏层在这个任务里会让训练变慢，边际收益递减。

5.4 扩展建议与最终提醒

我最大的体会是：这个习题的关键不是背公式，而是建立一套“工程调试”直觉。代码写完后，如果能在一分钟内想到“准确率 10% 先查标签编码，损失 nan 先查学习率，梯度对不上先查偏置列”，那这道题才算真正消化了。代码可以抄，但调试的经验一定得自己踩出来。

最后再建议你做一个很小的改动：把代码里的隐藏层从 25 改成 50，然后观察准确率和训练时间的变化。这种动手测试比单纯看书更能帮助你理解网络容量对模型的影响。踩过几次坑之后，再回头看反向传播，你会觉得它只是链式法则加一点矩阵运算，没有想象中那么难。