从零推导反向传播：计算图、自动微分与梯度排查实战

第一次把反向传播的公式从头推完，是在一个加班到凌晨的夜里。当时手上有个两层的小网络，训练损失死死卡在0.69不动——玩过二分类的人都知道，这是交叉熵在瞎猜时的理论值。我一边打印每层梯度一边排查，发现输出层的梯度量级正常，第一层的权重梯度却小到1e-8。问题既不在优化器也不在学习率，而是我压根没弄明白梯度是怎么一层一层传回去的。后来把链式法则老老实实手推了三遍，又用几百行NumPy重写了一遍前向与反向，那种“原来如此”的感觉才真正落地——这篇文章就是把那次踩坑和后续积累的整套经验完整摊开讲。

先把定位讲清楚：**反向传播（Backpropagation）**不是某种网络结构，也不是优化算法本身，它是深度学习里计算损失函数对所有参数梯度的核心方法。有了它，一个百万级参数的网络才能在一次前向加一次反向里拿到全部梯度；没有它，你只能用数值差分一个参数一个参数地试，代价高到完全不可接受。这篇内容适合三类人：刚学完吴恩达课程但一写代码就懵的入门者、能跑通框架却说不清梯度流向的中间层，以及面试前想系统梳理一遍细节的同学。我会从计算图的视角讲透“为什么是反向”，再手推两层网络的完整梯度，然后从零实现一个能真正收敛的自动微分内核，最后集中处理梯度消失、梯度校验和局部最小值这几个高频疑问。

1. 反向传播到底在解决什么问题

1.1 从“调参靠猜”到梯度可算

要理解反向传播的价值，先得想象没有它的世界长什么样。假设你有一个一百万个参数的网络，想通过试错找到让损失下降的方向，最朴素的做法是数值差分：把第 i 个参数加上一个极小量，重新跑一次前向，看损失变了多少，这个变化率就是该参数的近似梯度。问题在于，算一个参数的梯度就要一次完整前向，算一百万个参数就要一百万次前向。一次前向在真实数据上可能要几百毫秒，一百万次就是几十个小时，而且是每更新一步都这样。这条路在工程上直接封死。

梯度下降之所以能跑起来，靠的是一个关键事实：我们真正需要的不是损失值本身，而是损失关于每个参数的偏导数。有了这组偏导，参数更新就变成了简单的一维问题——每个参数沿着自己那个方向按比例走一小步。反向传播把“求全部偏导”这件事的代价从 O(N) 次前向压到了 O(1) 次前向加 O(1) 次反向，一次遍历拿到全部梯度，这才是深度学习能规模化的根本原因。

打个比方，你在浓雾里下山，看不见谷底在哪。你能做的只有弯腰摸一摸脚边地面的倾斜方向，然后朝着最陡的下坡方向迈一小步。梯度就是那个“倾斜方向”，反向传播就是让你不用把整座山的每一寸地面都摸一遍，只沿着你走过的路径回推一次，就能知道每个落脚点该往哪偏。

1.2 计算图视角：为什么一定是“反向”

很多人背下了链式法则，却没想过为什么这套流程叫“反向”传播。把它放到计算图里看会很清楚。前向传播时，每一步运算都是一个节点：输入乘以权重得到中间变量，经过激活函数得到输出，再到损失。这些节点连成一张有向无环图，方向从输入指向损失。所谓求导，就是在这张图上把“损失对某个参数的敏感度”一路乘回去。

链式法则本身：若 z 依赖 y，y 依赖 x，那么 dz/dx = dz/dy · dy/dx。要把损失对最底层参数的导数算出来，就得沿着图从损失往回走，每经过一个节点就乘上该节点的局部导数。这个“往回走”的顺序，正是“反向”两个字的由来。

那为什么不用前向模式？前向模式是从输入出发，同时传播“这个输入的变化会引起输出多大变化”。它的问题在于：每换一个输入变量作为起点，就得重新完整遍历一次计算图。对于神经网络这种“输入千万个参数、输出只有一个标量损失”的结构，前向模式要遍历一千万次，反向模式一次就够。这背后其实是雅可比矩阵和向量乘积的结合——反向传播每一步算的都是“某个局部雅可比转置乘以回传梯度”，把矩阵和向量相乘的顺序安排得恰到好处，避免了大矩阵显式构造。

提示：判断一个自动微分实现是前向还是反向，看它从哪端开始传播即可。深度学习框架全部采用反向模式，原因就是上面这个维度不对称。

1.3 三个必须先统一的约定

动手推公式之前，有几个约定不统一，后面必然乱。第一，梯度的形状永远和被求导张量的形状一致。权重矩阵 W 是 (m, n)，那它的梯度也必须是 (m, n)，这样更新时才能逐元素相减。第二，反向传播里所有加法都是累加而不是覆盖。同一个张量如果在前向中被用了两次，它的梯度必须来自两条路径的和，这也是为什么实现里写的是grad += ...。第三，标量损失才有意义。整个反向传播的起点是损失这一个标量，它的初始回传梯度恒为 1，相当于把“损失对自身的导数”设定为基准。

这三条看着琐碎，但几乎每个初学者的第一次实现翻车，都出在其中某一条上。尤其是第二条，用覆盖而不是累加，网络在小数据集上还能凑合收敛，一旦遇到有分支的结构就会静默出错，非常难查。

2. 手推一遍才算真懂

2.1 单个神经元的梯度长什么样

先从最小的单元开始。一个二分类神经元，前向是三步：线性加权 z = wᵀx + b，激活 a = σ(z)，损失用二分类交叉熵 L = -[y·ln a + (1-y)·ln(1-a)]。反向要回答的是：w 和 b 各自变化一点点，L 会怎么变。

先用链式法则拆 L 对 a 的导数。对 a 求导，ln a 那项给出 -y/a，ln(1-a) 那项给出 (1-y)/(1-a)，整体是 (a - y) / [a(1-a)]。再看 a 对 z 的导数，sigmoid 的导数有个漂亮的性质：σ'(z) = σ(z)(1 - σ(z)) = a(1-a)。把这两项相乘，分母上的 a(1-a) 正好被消掉，剩下 dz 上的回传梯度就是简洁的(a - y)。

再往下就简单了：z 对 w 的导数是 x，z 对 b 的导数是 1。所以 dL/dw = (a - y)·x，dL/db = (a - y)。这三个字母看起来平淡，但它揭示了一个非常实用的规律——sigmoid 配交叉熵，输出层的梯度就是“预测减真值”。这个结论对后面的 softmax 也成立，记住它，很多推导能省一大半力气。

注意：这个简洁结果依赖“交叉熵 + sigmoid”这对组合，如果你换成均方误差，梯度里会多出 a(1-a) 这个因子，量级立刻变小一截。输出层用什么损失，直接决定了回传梯度的尺度。

2.2 两层网络的完整手推过程

把规模抬到两层，用带批量的矩阵形式，这样推导结果能直接落到代码里。设输入批次 X 的形状为 (n_in, B)，第一层权重 W1 为 (n_h, n_in)，偏置 b1 为 (n_h, 1)，激活用 ReLU。第一层前向：Z1 = W1·X + b1，A1 = ReLU(Z1)。第二层 W2 为 (n_out, n_h)，b2 为 (n_out, 1)，输出 logits：Z2 = W2·A1 + b2。经过 softmax 得到概率 P，再配交叉熵算平均损失 L。

反向从 L 开始。前面提过 softmax 配交叉熵的组合梯度极其干净：dL/dZ2 = (P - Y) / B，其中 Y 是标签的独热形式，B 是批量大小，除以 B 是因为用了平均损失。这一步是整个推导里收益最高的结论，值得单独记住。

拿到 dZ2 之后逐层回推。Z2 = W2·A1 + b2，所以 dL/dW2 = dZ2 · A1ᵀ，形状是 (n_out, n_h) 对上了；dL/db2 是把 dZ2 沿批量维度求和，得到 (n_out, 1)；dL/dA1 = W2ᵀ · dZ2，形状回到 (n_h, B)。接下来穿过 ReLU：ReLU 的导数在输入大于零处为 1，否则为 0，所以dL/dZ1 = dL/dA1 ⊙ (Z1 > 0)，这里的 ⊙ 表示逐元素相乘。最后一步形式和第一层对称：dL/dW1 = dL/dZ1 · Xᵀ，dL/db1 是 dL/dZ1 沿批量方向求和。

整个链条的手感是：每一个线性层都有一个“权重梯度 = 回传梯度 × 输入转置”和“偏置梯度 = 回传梯度按批量求和”的固定套路，激活函数只负责在回传路径上乘一个逐元素的导数。把这四行记住，任何前馈网络的梯度你都能在纸上推出来。

2.3 维度对齐是最有效的自查手段

推完公式后，强烈建议把每一步的形状写在旁边。以上面为例：dZ2 是 (n_out, B)，A1ᵀ 是 (B, n_h)，两者相乘得到 (n_out, n_h)，和 W2 一致，通过。如果某一步算出来的形状和参数对不上，那一定是转置或者乘法顺序错了，不用怀疑数学，直接查形状。

这套方法在排查真实代码时特别管用。框架报错往往只告诉你“维度不匹配”，但不会告诉你逻辑哪里错了；而如果你在脑子里已经跑过一遍形状推导，定位通常不超过两分钟。我自己的习惯是，写任何自定义层之前先在注释里把输入输出形状标好，正反向都标，这样连调试的时间都能省掉一大块。

3. 从零写一个能跑通的反向传播

3.1 自动微分内核怎么设计

搞清楚公式之后，我想真正理解自动微分的工作方式，于是写了这版最小内核。核心思想是把每个张量包装成一个节点，节点内部记录三样东西：数据、梯度、以及生成自己的子节点和反向函数。前向运算时，每做一次加法或乘法，就顺手把“如何反向”这个闭包挂到输出节点上。反向时先对整张图做一次拓扑排序，再从损失出发逆序遍历，逐个调用反向函数。

拓扑排序这一步很关键。计算图本质上是有向无环图，直接递归反向可能在某些节点还没收集完所有回传梯度时就提前计算，导致梯度少算。用后序遍历保证“所有依赖它的节点都处理完之后才处理它”，逆序执行就天然满足反向传播的顺序要求。

import numpy as np def _unbroadcast(grad, shape): """把广播后的梯度还原成原始形状""" while grad.ndim > len(shape): grad = grad.sum(axis=0) for i, s in enumerate(shape): if s == 1 and grad.shape[i] != 1: grad = grad.sum(axis=i, keepdims=True) return grad class Tensor: def __init__(self, data, requires_grad=False, _children=(), _op=""): self.data = np.asarray(data, dtype=np.float64) self.requires_grad = requires_grad self.grad = None self._prev = tuple(_children) self._op = _op self._backward = lambda: None def __add__(self, other): other = other if isinstance(other, Tensor) else Tensor(other) out = Tensor(self.data + other.data, self.requires_grad or other.requires_grad, (self, other), "+") def _backward(): g = _unbroadcast(out.grad, self.data.shape) self.grad = g if self.grad is None else self.grad + g g2 = _unbroadcast(out.grad, other.data.shape) other.grad = g2 if other.grad is None else other.grad + g2 out._backward = _backward return out def __mul__(self, other): other = other if isinstance(other, Tensor) else Tensor(other) out = Tensor(self.data * other.data, self.requires_grad or other.requires_grad, (self, other), "*") def _backward(): g = _unbroadcast(out.grad * other.data, self.data.shape) self.grad = g if self.grad is None else self.grad + g g2 = _unbroadcast(out.grad * self.data, other.data.shape) other.grad = g2 if other.grad is None else other.grad + g2 out._backward = _backward return out def __matmul__(self, other): other = other if isinstance(other, Tensor) else Tensor(other) out = Tensor(self.data @ other.data, self.requires_grad or other.requires_grad, (self, other), "@") def _backward(): g = out.grad @ other.data.T self.grad = g if self.grad is None else self.grad + g g2 = self.data.T @ out.grad other.grad = g2 if other.grad is None else other.grad + g2 out._backward = _backward return out def relu(self): out = Tensor(np.maximum(self.data, 0.0), self.requires_grad, (self,), "relu") def _backward(): g = out.grad * (self.data > 0) self.grad = g if self.grad is None else self.grad + g out._backward = _backward return out def softmax_cross_entropy(self, labels): """logits: (C, B)，labels: (B,)，返回标量损失""" z = self.data - self.data.max(axis=0, keepdims=True) e = np.exp(z) p = e / e.sum(axis=0, keepdims=True) B = labels.shape[0] loss = float(-np.log(p[labels, np.arange(B)] + 1e-12).mean()) out = Tensor(loss, self.requires_grad, (self,), "sce") def _backward(): dp = p.copy() dp[labels, np.arange(B)] -= 1.0 g = dp / B self.grad = g if self.grad is None else self.grad + g out._backward = _backward return out def backward(self): topo, visited = [], set() def build(v): if id(v) not in visited: visited.add(id(v)) for c in v._prev: build(c) topo.append(v) build(self) for v in topo: v.grad = None self.grad = np.ones_like(self.data) for v in reversed(topo): v._backward()

这段代码不长，但把自动微分的骨架全交代了。几个细节值得展开：_unbroadcast处理的是偏置加法时的广播问题——(n_h, 1)的偏置加到(n_h, B)上，反向时梯度必须沿批量维度求和还原成(n_h, 1)，否则更新时形状对不上；softmax_cross_entropy里先减去每列最大值，是防止 exp 溢出，这个减常数操作不影响最终梯度，因为 softmax 本身对平移不变；反向函数里一律是+=而不是=，保证同一个张量被多次引用时梯度能正确累加。

3.2 训练一个两层网络验证正确性

为了验证这套内核真能学习，我用双月牙数据搭一个两层网络。这个数据集是非线性可分的，线性模型最多只能到八成左右，能跑上去说明隐藏层和反向传播确实在起作用。

rng = np.random.default_rng(42) def make_moons(n=400, noise=0.12, seed=0): r = np.random.default_rng(seed) n1 = n // 2 t1 = np.linspace(0, np.pi, n1) x1 = np.stack([np.cos(t1), np.sin(t1)], axis=1) t2 = np.linspace(0, np.pi, n - n1) x2 = np.stack([1 - np.cos(t2), 1 - np.sin(t2) - 0.5], axis=1) X = np.concatenate([x1, x2], 0) + r.normal(0, noise, (n, 2)) y = np.concatenate([np.zeros(n1), np.ones(n - n1)]).astype(int) return X, y X, y = make_moons(400, 0.12, 0) X = (X - X.mean(0)) / X.std(0) W1 = Tensor(rng.normal(0, 0.5, (32, 2)), requires_grad=True) b1 = Tensor(np.zeros((32, 1)), requires_grad=True) W2 = Tensor(rng.normal(0, 0.5, (2, 32)), requires_grad=True) b2 = Tensor(np.zeros((2, 1)), requires_grad=True) params = [W1, b1, W2, b2] lr = 0.1 for epoch in range(2000): xb = Tensor(X.T) # (2, B) h = (W1 @ xb + b1).relu() # (32, B) logits = W2 @ h + b2 # (2, B) loss = logits.softmax_cross_entropy(y) loss.backward() for p in params: p.data -= lr * p.grad if epoch % 200 == 0: acc = (logits.data.argmax(0) == y).mean() print(f"epoch {epoch:4d} loss {loss.data:.4f} acc {acc:.3f}")

实测下来，初始损失大约在0.693（也就是随机猜的水平），几十轮之后快速下降，两千轮时准确率通常能到九成五以上。这个结果和官方框架跑出来的几乎一致，说明手写的反向传播逻辑是站得住的。学习率这里定的是0.1，配合全部样本做批梯度下降比较稳；如果你改成小批量，学习率可以适当放大，但要注意把批量大小体现在损失的平均里，否则梯度尺度会跟着批量大小变化。

注意：每次迭代之后并不需要手动清零梯度，因为backward()一开始会把拓扑图中所有节点的梯度重置为 None。这个设计省去了一处常见的手写错误，但代价是每轮都要重新分配一批数组，训练大规模网络时会有额外开销，真实框架里是用 zero_grad 单独做的。

3.3 参数选择的来龙去脉

初始化的标准差选0.5、隐藏层宽度选32、学习率0.1，这些不是随手写的。隐藏层宽度32对于双月牙这种二维问题是足够的，太窄（比如4）会导致决策边界不够弯曲，太宽则容易在小样本上过拟合。初始化用正态分布而不是全零，是因为全零会让同一层的所有神经元在任何输入下都输出相同结果，反向梯度也完全相同，网络永远学不出差异；标准差过大又会让激活值分布发散，梯度爆炸。0.5这个量级在实践中对浅层网络很少出问题，深层网络则需要更精细的方案，后面会讲到。

学习率0.1对应的是全批量梯度下降。这个配置下损失曲线平滑、单调下降，适合验证反向传播实现的正确性。如果想要更快的收敛速度，可以换成带动量的更新，但那属于优化器层面的改进，不改变反向传播本身的逻辑。

4. 梯度问题排查与校验实录

4.1 数值梯度校验的标准做法

写完了反向传播，怎么知道自己推的公式没错？最可靠的办法是数值梯度校验：对每个参数施加一个极小的扰动，用中心差分近似导数，再跟解析梯度比较相对误差。中心差分比单边差分精度高一个量级，用的是 (f(x+ε) - f(x-ε)) / (2ε)。

def numerical_grad(f, x, eps=1e-6): """f: 接受 numpy 数组返回标量的函数, x: numpy 数组""" grad = np.zeros_like(x) it = np.nditer(x, flags=['multi_index']) while not it.finished: idx = it.multi_index old = x[idx] x[idx] = old + eps f1 = f(x) x[idx] = old - eps f2 = f(x) x[idx] = old grad[idx] = (f1 - f2) / (2 * eps) it.iternext() return grad

实践中的判据是相对误差|数值梯度 - 解析梯度| / max(1e-8, |数值梯度| + |解析梯度|)：小于1e-7说明实现完全正确，1e-4到1e-7之间可接受，超过1e-4就要警惕。用的时候有几个坑要避开：epsilon取1e-4到1e-6之间比较合适，太大会被截断误差污染，太小会被浮点精度吃掉；校验时应该切换到双精度，单精度下差分本身就不可靠；如果损失里含 ReLU，扰动点恰好落在零点附近会让数值梯度变得很怪，这种位置直接跳过。

还有一个容易被忽略的细节：梯度校验要在随机初始化之后、训练之前做。训练几轮之后参数已经跑到某些激活的死亡区域，很多神经元输出恒为零，梯度信息本来就少，校验反而看不出问题。我自己习惯是每次改完一个自定义层的反向实现，先跑一遍小批量数值校验，通过了再上真实数据。

4.2 梯度消失和梯度爆炸的成因

前面手推时已经看到，反向传播每穿过一层，梯度就要乘上那一层激活函数的导数。sigmoid 的导数最大值只有0.25，如果网络有十层，理想情况下梯度最大也要乘0.25的十次方，大约是1e-6；实际情况更糟，因为输入很少恰好落在导数峰值处。这就是梯度消失：靠近输入层的参数几乎收不到有效的更新信号，训练表现为损失卡住、前几层权重几乎不变。

反过来，如果每一层的回传因子都大于1，连乘之后梯度指数增长，一次更新就把参数推到很远的地方，损失直接变成 NaN，这是梯度爆炸。它在循环网络里尤其常见，因为同一组权重被反复使用，连乘的层数等于序列长度。

He 初始化针对 ReLU 设计，方差取 2/fan_in，正好抵消 ReLU 把一半激活置零带来的方差减半；Xavier 初始化面向 tanh，方差取 1/fan_in 和 1/fan_out 的调和平均。这两个推导都基于“保持每一层前向输出方差不变”的目标，理解了目标，选择哪种初始化就不难判断。梯度裁剪则是硬性兜底：算出梯度后先看它的范数，超过阈值就按比例缩回去，代价是可能改变梯度方向，所以阈值要设得比正常梯度略高一点。

4.3 常见报错速查

调试反向传播时，报错往往集中在几类。第一类是形状不匹配，八成是转置漏了或者乘法顺序反了，按第2.3节的形状表逐项核对，通常几分钟能定位。第二类是梯度全是零，常见原因是激活函数的掩码写错，比如把(z > 0)写成(z >= 0)影响不大，但如果写成(self.data > 0)而 self 是激活后的输出，就会全错。第三类是损失不收敛但也不报错，这时候先怀疑学习率，再怀疑标签编码是否和损失函数匹配——用 softmax 交叉熵却传了独热标签而不是索引，是高频错误。

还有一种隐蔽情况：某些批次损失正常，某些批次突然爆掉。这往往和输入没有归一化有关，个别样本的数值范围远大于其他样本，反传时梯度被放大。解决方式是统一做标准化，或者把异常样本剔除。我自己的经验是，任何梯度异常先查数据，再看公式，最后才怀疑框架，顺序反过来会浪费大量时间。

5. 反向传播能解决局部最小值问题吗

5.1 局部最小值、鞍点和高原的区别

这个问题在面试和考试里出现频率极高，也是理解深浅的分水岭。直接给结论：反向传播不能解决局部最小值问题，因为它根本不是用来解决这个问题的。反向传播只负责一件事——给定当前参数，算出损失对每个参数的梯度。至于这个梯度把你带到哪里，是优化算法的事。梯度下降遇到局部最小值会停，遇到鞍点会减速，遇到高原会缓慢挪动，这些都发生在“拿到梯度之后”，跟反向传播无关。

更值得说的是，在高维参数空间里，真正的局部最小值其实很少见。一维函数里局部最小值是常见的，但参数维度上升到百万级，损失曲面上绝大多数梯度为零的点都是鞍点——某些方向往上、某些方向往下，只是整体梯度恰好为零。维数越高，每个方向都恰好朝上的概率越低，鞍点的比例就越大。而随机初始化基本不会恰好落在某个局部最小值上，加上训练时用的是小批量，每一步的损失曲面都略有不同，噪声本身就会把参数从浅坑里推出来。

真正让人头疼的其实是高原区域，也就是梯度长时间接近零但损失仍然很高的平坦地带。它不像局部最小值那样“看起来到底了”，但训练曲线会表现为长时间纹丝不动，容易让人误判成模型容量不够。

5.2 工程上真正有效的几招

既然局部最小值大多不是瓶颈，那实际训练中该关注什么？第一是动量。它把历史梯度按指数加权累积起来，在鞍点附近即使当前梯度很小，累积的动量也能把参数推过去。动量还有一个副作用是能抑制震荡，在峡谷型曲面上收敛更快。第二是自适应学习率方法，AdaGrad、RMSProp、Adam 这一族会为每个参数单独调整步长，梯度一直很小的维度会获得更大的有效步长，这在特征尺度差异大的场景里帮助明显。

第三是随机性本身。小批量梯度下降每个批次用的是不同样本，损失曲面不断变化，相当于给优化过程加了噪声，参数不容易卡死在某个点上。批量越小噪声越大，跳出浅坑的能力越强，但收敛会更抖；这是一个需要根据任务调的平衡。第四是残差连接和归一化，它们从结构上保证梯度能顺畅回传，让深层网络的优化曲面更平滑，这在实践中带来的收益往往比调优化器更大。

第五个是学习率调度。训练初期用较大的学习率快速下降到低损失区域，后期逐渐减小让参数稳定下来。常见的做法是余弦退火或者阶梯式衰减，配合预热处理。这些手段组合起来，实际训练中几乎不会因为“卡在局部最小值”而失败，反倒是梯度消失、数据质量、超参数设置这几项更容易出问题。

所以回到标题里的那个疑问：与其把精力放在“怎么逃出局部最小值”，不如先把反向传播实现正确、把梯度量级监控起来、把数据预处理做好。这三点做到位，绝大多数训练问题都会迎刃而解。

最后分享一个我自己一直在用的小习惯：训练脚本里加一段代码，每轮打印所有参数梯度的全局范数和每层范数。梯度范数突然增大说明学习率偏高或者数据有问题，持续减小说明梯度在衰减，长期停在极小值说明前几层根本没参与学习。这个信息比损失曲线更早暴露问题，在很多次排查里帮我省下了大把时间。