亲手实现反向传播

1. 开场：把九课的工具拼成一台机器

上一课我们把矩阵求导 matrix calculus讲透了：知道了 \(\partial L/\partial W\) 长什么形状、怎么用链式法则一层层往回传。这一课是模块 1 的压轴实战——我们不再做孤立的公式练习，而是把前九课的全部工具（向量、矩阵变换、内积、ReLU、softmax、链式法则、矩阵求导）拼装成一台真正会"学习"的机器：一个两层全连接神经网络 two-layer fully-connected network，并亲手实现它的反向传播 backpropagation。

读完本课，你会从头到尾推导出每一个梯度公式、用 numpy 把它们写成几行代码、用梯度检查 gradient check证明它们一字不差，然后看着 loss 从高到低真实地降下来。这就是 PyTorch 里那句 loss.backward() 背后的全部真相——理解了它，以后遇到 NaN、梯度爆炸、梯度消失，你才知道去哪儿找毛病。

2. 网络长什么样：前向传播

我们要搭的网络结构非常标准：输入向量先过一个线性层，激活，再过一个线性层，最后用 softmax 变成概率，用交叉熵衡量它和真实标签差多少。画成计算图就是一条链：

先看单个样本（列向量 \(x\in\mathbb{R}^{D}\)）的前向公式。设隐藏层有 \(H\) 个单元、输出 \(C\) 个类别：

这里 \(y\) 是独热向量 one-hot vector：真实类别那一位是 1，其余是 0。所以交叉熵其实只剩一项 \(L=-\log\hat{y}_{\text{正确类}}\)——模型给正确类的概率越高，loss 越小。\(\mathrm{softmax}\) 的第 \(c\) 个分量是

它把任意一组实数"挤"成一个和为 1 的概率分布。

3. 数值稳定：softmax 的 log-sum-exp 技巧

在写代码之前必须先排一个雷。softmax 里有 \(e^{z}\)，如果某个 \(z\) 比较大（比如 \(z=1000\)），\(e^{1000}\) 会直接上溢 overflow 成 inf，整个计算变成 nan。解决办法是一个简单又精妙的恒等式：分子分母同乘 \(e^{-m}\) 不改变结果，

我们取 \(m=\max_k z_k\)（每一行各自的最大值）。这样减完之后，指数的输入最大是 0，\(e^{0}=1\) 不会溢出；最小的那些变成很负的数，\(e^{-\text{大}}\to 0\) 安全地下溢成 0。这就是log-sum-exp 稳定化。代码里永远写 Z - Z.max(axis=1, keepdims=True)，一行都不能省。

4. 反向传播：从损失一路倒推回去

现在是本课的心脏。我们要的是 \(\partial L/\partial W_1,\ \partial L/\partial b_1,\ \partial L/\partial W_2,\ \partial L/\partial b_2\)。策略是沿计算图反向走，每经过一个节点，就用上一课的链式法则把梯度往前一站传递。绿色是前向数据流，橙色是反向梯度流：

第一步（最关键）：softmax + 交叉熵合并求导

单独对 softmax 求导很麻烦（要算雅可比矩阵），但 softmax 紧接着交叉熵时，会发生一次"奇迹般"的化简。结论极其干净：

也就是预测概率减去真实独热标签。直觉上太美了：如果模型对正确类已经预测为 1（\(\hat y=y\)），梯度就是 0，不用再改；预测得越离谱，梯度越大，改得越狠。这个结果第 9 课已经推过，这里直接拿来当反向传播的"起手式"。我们记 \(\delta_2 := \partial L/\partial z_2\)。

第二步：穿过第二层线性 \(z_2=W_2 a_1+b_2\)

这是上一课的标准结论。对一个线性层 \(z=Wa+b\)，已知上游梯度 \(\delta=\partial L/\partial z\)，则（分母布局下）：

记住形状自检法：\(\delta_2\) 是 \(C\times1\)，\(a_1^\top\) 是 \(1\times H\)，外积得到 \(C\times H\)，正好和 \(W_2\) 同形——分母布局的约定让一切自动对齐。最后一项 \(W_2^\top\delta_2\)（\(H\times C\) 乘 \(C\times 1\) 得 \(H\times1\)）就是把梯度回传给 \(a_1\)，准备穿过 ReLU。

第三步：穿过 ReLU

ReLU 是逐元素的 \(a_1=\max(0,z_1)\)，它的导数是一个"开关"：输入为正时导数为 1（信号原样通过），输入为负时导数为 0（信号被掐断）。所以

其中 \(\odot\) 是逐元素乘 element-wise product，\(\mathbf{1}[z_1>0]\) 是一个由 0/1 组成的掩码（mask）。记 \(\delta_1:=\partial L/\partial z_1\)。

第四步：穿过第一层线性 \(z_1=W_1 x+b_1\)

和第二步一模一样的模式，只是把 \(a_1\) 换成输入 \(x\)：

到这里四个梯度全部到手。注意整个过程只用了前向时存下来的中间量（\(x, z_1, a_1, \hat y\)）和一路传下来的 \(\delta\)，没有任何重复计算——这正是反向传播比"逐个参数算数值导数"快上千万倍的原因。

5. 例题：手推一遍全部梯度（带真实数字）

6. batch 维度：N 个样本一起算

真实训练不会一次只喂一个样本，而是一批 batch \(N\) 个一起。沿用"行=样本"的约定，把数据堆成矩阵 \(X\in\mathbb{R}^{N\times D}\)（每行一个 \(x^\top\)）。前向时所有公式从"矩阵×列向量"变成"矩阵×矩阵"，且权重摆在右边：

（这里 \(W_1\) 取 \(D\times H\)、\(W_2\) 取 \(H\times C\)，正好让"行样本"流畅地左乘；偏置 \(b\) 靠广播加到每一行。）损失是整批的平均：\(L=\frac1N\sum_{i}\mathrm{CE}(\hat y_i,y_i)\)。反向时的关键约定是：

• 权重梯度对 batch 求和、再除以 \(N\)（即平均）。\(\partial L/\partial Z_2=(\hat Y-Y)/N\)，于是 \(\partial L/\partial W_2=A_1^\top\,(\hat Y-Y)/N\)。这个矩阵乘法里的求和 \(\sum_i\) 自动把所有样本的贡献加起来，除以 \(N\) 就是平均。
• 偏置梯度沿样本维 axis=0 求和，因为同一个 \(b\) 被加到了每个样本上，它的总梯度就是各样本贡献之和。注意求和与平均在这里只差一个"已经提前除掉的 \(N\)"：我们传进来的本就是带 \(1/N\) 因子的 \((\hat Y-Y)/N\)，所以 db = ((Ŷ−Y)/N).sum(axis=0) 沿样本维加完，得到的恰好是 \(N\) 个样本梯度的平均——和权重梯度口径一致，无需再额外除 \(N\)。

形状上：\(A_1^\top\) 是 \(H\times N\)，\((\hat Y-Y)/N\) 是 \(N\times C\)，乘出来 \(H\times C\) 和 \(W_2\) 同形。漂亮地对齐。

7. 梯度检查：怎么证明你没推错

推导和写代码都极易出符号错、转置错、漏除 \(N\)。怎么自证清白？用有限差分 finite difference 去逼近"真"梯度，和你的解析梯度比。对任意一个参数分量 \(\theta\)，用中心差分（误差是 \(O(\epsilon^2)\)，比单边差分准得多）：

然后比较解析梯度 \(g_a\) 与数值梯度 \(g_n\) 的相对误差 relative error：

如果你的反向传播写对了，相对误差应当小于 \(10^{-7}\)（用 float64 时常能到 \(10^{-9}\)）。如果某个参数误差是 \(10^{-2}\) 量级，那里几乎肯定有 bug。下面代码里 grad_check() 抽查若干分量，实测最大相对误差约 \(5\times10^{-9}\)——过关。

8. 完整代码：合成数据 → 训练 → 梯度检查全过关

下面是一段纯 numpy、可直接在本页运行的完整实现：现场合成三团高斯点做三分类，跑梯度检查，再训练 800 步，打印 loss 下降与准确率。几秒内跑完。

import numpy as np
rng = np.random.default_rng(42)

# ---------- 1. 合成数据：3 团高斯点（故意让它们有重叠，问题才不平凡）----------
def make_blobs(n_per=200):
    centers = np.array([[ 1.3,  0.0],
                        [-1.0,  1.1],
                        [-0.7, -1.1]])
    X, y = [], []
    for k, c in enumerate(centers):
        pts = rng.normal(c, 0.95, size=(n_per, 2))   # 标准差大 -> 三团互相重叠
        X.append(pts); y.append(np.full(n_per, k))
    X = np.vstack(X); y = np.concatenate(y)
    idx = rng.permutation(len(y))
    return X[idx], y[idx]

X, y = make_blobs()
N, D = X.shape          # N 样本, D=2 特征（行=样本）
C = 3                   # 类别数
# np.eye(C) 是 C×C 单位矩阵；用标签数组 y 索引它的行，
# 把每个整数标签 k 取成第 k 行（第 k 位为 1 的独热向量），
# 一次性得到形状 (N, C) 的独热标签矩阵 Y。
Y = np.eye(C)[y]
H = 16                  # 隐藏单元数

# ---------- 2. 参数初始化（He 初始化，配 ReLU）----------
# He 初始化：把权重的随机尺度设成 sqrt(2/扇入)，目的是让前向信号
# 不逐层放大或缩小、训练能顺利起步。为什么恰好是 2/n 留到优化模块讲，
# 现在只需照用即可。
def init():
    return (rng.normal(0, np.sqrt(2/D), (D, H)),  # W1 (D,H)
            np.zeros(H),                            # b1
            rng.normal(0, np.sqrt(2/H), (H, C)),   # W2 (H,C)
            np.zeros(C))                            # b2

# ---------- 3. 前向 + 损失 ----------
def softmax(Z):
    Z = Z - Z.max(axis=1, keepdims=True)   # log-sum-exp 稳定化（按行减最大值）
    e = np.exp(Z)
    return e / e.sum(axis=1, keepdims=True)

def forward(X, W1, b1, W2, b2):
    z1 = X @ W1 + b1          # (N,H)
    a1 = np.maximum(z1, 0)    # ReLU
    z2 = a1 @ W2 + b2         # (N,C)
    yhat = softmax(z2)
    return yhat, (X, z1, a1, z2, yhat)

def loss_fn(yhat, Y):
    N = Y.shape[0]
    return -np.sum(Y * np.log(yhat + 1e-12)) / N   # 加 1e-12 防 log(0)

# ---------- 4. 反向（本课的核心，对照第 4、6 节公式）----------
def backward(cache, Y, W2):
    X, z1, a1, z2, yhat = cache
    N = X.shape[0]
    dz2 = (yhat - Y) / N        # 起手式 ŷ-y，并对 batch 平均
    dW2 = a1.T @ dz2            # (H,C)  = A1ᵀ δ2
    db2 = dz2.sum(axis=0)       # (C,)   偏置沿样本维求和
    da1 = dz2 @ W2.T            # (N,H)  回传到 a1
    dz1 = da1 * (z1 > 0)        # 过 ReLU：逐元素乘 0/1 掩码
    dW1 = X.T @ dz1             # (D,H)  = Xᵀ δ1
    db1 = dz1.sum(axis=0)       # (H,)
    return dW1, db1, dW2, db2

# ---------- 5. 梯度检查：有限差分 vs 解析梯度 ----------
def grad_check():
    W1, b1, W2, b2 = init()
    Xs, Ys = X[:8], Y[:8]                    # 小批量足够
    _, cache = forward(Xs, W1, b1, W2, b2)
    grads  = dict(zip(['W1','b1','W2','b2'], backward(cache, Ys, W2)))
    params = {'W1':W1, 'b1':b1, 'W2':W2, 'b2':b2}
    eps, worst = 1e-5, 0.0
    for name, P in params.items():
        g_ana = grads[name]
        for _ in range(20):                  # 每个参数抽查 20 个分量
            i = tuple(rng.integers(0, s) for s in P.shape)
            old = P[i]
            P[i] = old + eps; Lp = loss_fn(forward(Xs,W1,b1,W2,b2)[0], Ys)
            P[i] = old - eps; Lm = loss_fn(forward(Xs,W1,b1,W2,b2)[0], Ys)
            P[i] = old
            g_num = (Lp - Lm) / (2*eps)       # 中心差分
            rel = abs(g_num - g_ana[i]) / max(1e-12, abs(g_num) + abs(g_ana[i]))
            worst = max(worst, rel)
    print(f"梯度检查最大相对误差 = {worst:.2e}   (应 < 1e-7)")

grad_check()

# ---------- 6. 训练循环 ----------
W1, b1, W2, b2 = init()
eta = 0.4                                     # 学习率
print("开始训练：")
for step in range(1, 801):
    yhat, cache = forward(X, W1, b1, W2, b2)
    L = loss_fn(yhat, Y)
    dW1, db1, dW2, db2 = backward(cache, Y, W2)
    W1 -= eta*dW1; b1 -= eta*db1             # 梯度下降一步
    W2 -= eta*dW2; b2 -= eta*db2
    if step in (1, 50, 100, 200, 400, 800):
        acc = (yhat.argmax(1) == y).mean()
        print(f"  step {step:4d}   loss {L:.4f}   acc {acc:.3f}")

yhat, _ = forward(X, W1, b1, W2, b2)
print(f"最终训练准确率 = {(yhat.argmax(1)==y).mean():.3f}")

运行后你会看到类似输出（数字因随机种子固定而可复现）：梯度检查约 5e-09 过关；loss 从 step 1 的 1.8 左右一路降下来，前 50 步降得最猛（50 步时已到 0.42 上下、基本收敛），之后缓慢微降，到 800 步稳定在 0.40 附近，准确率到 0.84 左右。三团点故意做成重叠的，所以达不到 100%——这恰恰说明它在真学，而不是死记。下图是 loss 随步数下降的示意（前陡后缓）：

试着把 eta 调到 5.0（loss 可能震荡或变 nan），或调到 0.01（降得极慢），亲手感受学习率的作用；再把 H 改成 2，看看模型容量不够时准确率怎么掉下来。

9. micrograd：把"手推链式"升级成"自动反向"

上面我们是手工为这个特定网络写 backward。但框架（PyTorch/TensorFlow）不可能为每种网络手写——它们用自动微分 automatic differentiation。核心思想由 Karpathy 的玩具库 micrograd 讲得最透彻：把每个标量包成一个带 grad 字段的节点，记住它是由哪些节点、经什么运算得来的；每个运算自己知道怎么把上游梯度传给输入。前向时顺手把计算图连起来，调用 .backward() 时按拓扑序倒着走一遍，每个节点执行自己的局部链式法则——这正是我们第 4 节做的事，只是自动化了。

import numpy as np

class Value:
    """一个标量 + 它的梯度 + 它怎么算出来的（自动微分最小内核）"""
    def __init__(self, data, _children=()):
        self.data = float(data)
        self.grad = 0.0
        self._backward = lambda: None     # 该节点把梯度回传给输入的局部规则
        self._prev = set(_children)

    def __add__(self, o):
        o = o if isinstance(o, Value) else Value(o)
        out = Value(self.data + o.data, (self, o))
        def _bw():                        # d(a+b)/da = 1, d(a+b)/db = 1
            self.grad += out.grad
            o.grad    += out.grad
        out._backward = _bw
        return out

    def __mul__(self, o):
        o = o if isinstance(o, Value) else Value(o)
        out = Value(self.data * o.data, (self, o))
        def _bw():                        # d(a*b)/da = b, d(a*b)/db = a
            self.grad += o.data * out.grad
            o.grad    += self.data * out.grad
        out._backward = _bw
        return out

    def relu(self):
        out = Value(max(0.0, self.data), (self,))
        def _bw():                        # ReLU 导数：输入>0 时 1，否则 0
            self.grad += (out.data > 0) * out.grad
        out._backward = _bw
        return out

    def backward(self):
        # 1) 拓扑排序：保证回传一个节点前，它的所有"下游"已处理完。
        #    直觉像倒水：必须先处理离输出最近的节点，再处理喂给它的上游节点；
        #    这样轮到某节点时，所有用到它的下游节点都已把各自梯度累加进来。
        topo, seen = [], set()
        def build(v):
            if v not in seen:
                seen.add(v)
                for ch in v._prev:
                    build(ch)
                topo.append(v)
        build(self)
        # 2) 从输出开始，梯度设 1，倒序执行每个节点的局部规则
        self.grad = 1.0
        for v in reversed(topo):
            v._backward()

# 验证：算一个神经元 n = relu(x*w + b) 对各输入的梯度
x, w, b = Value(2.0), Value(3.0), Value(1.0)
n = (x*w + b).relu()        # relu(2*3+1) = relu(7) = 7
n.backward()
print("n =", n.data)                  # 7.0
print("dx =", x.grad)                 # = w = 3.0
print("dw =", w.grad)                 # = x = 2.0
print("db =", b.grad)                 # = 1.0

# 换成会被 ReLU 掐断的情形：relu(-5)=0，梯度全为 0
x2, w2, b2 = Value(-2.0), Value(3.0), Value(1.0)
n2 = (x2*w2 + b2).relu()    # relu(-6+1) = relu(-5) = 0
n2.backward()
print("n2 =", n2.data, " dx =", x2.grad, " dw =", w2.grad, " db =", b2.grad)

运行会打印 n=7.0, dx=3.0, dw=2.0, db=1.0，第二组因 ReLU 截断全是 0——和我们手推的规则完全一致。把 Value 扩展到向量/矩阵、再加上 exp、log 等运算，就是一个迷你 PyTorch。第 8 节那个网络的每一步 .backward()，本质就是这套机制在矩阵上的批量版本。

10. 小结：你刚刚打通了整个模块 1

回头看，这一课把前九课串成了一条完整链路：向量与矩阵是数据和参数的载体，内积/线性变换是每一层的运算，导数与梯度是"往哪改"的方向，链式法则与矩阵求导是反向传播的引擎，而 softmax/交叉熵把它接到分类任务上。你现在能从损失一路手推到第一层权重的梯度、用代码实现它、用梯度检查证明它对、看着它把 loss 真的降下来——这是理解一切深度模型训练的地基。

调一调，观察现象

下面三个微任务都基于本课已验证的代码，改一个地方、跑几秒，亲手看反向传播的几个关键约定"坏掉"或"生效"。建议每个都先猜结果再运行。

import numpy as np
rng = np.random.default_rng(42)
def make_blobs(n_per=200):
    centers = np.array([[1.3,0.0],[-1.0,1.1],[-0.7,-1.1]])
    X,y=[],[]
    for k,c in enumerate(centers):
        X.append(rng.normal(c,0.95,size=(n_per,2))); y.append(np.full(n_per,k))
    X=np.vstack(X); y=np.concatenate(y); idx=rng.permutation(len(y))
    return X[idx],y[idx]
X,y=make_blobs(); N,D=X.shape; C=3; Y=np.eye(C)[y]; H=16
def softmax(Z):
    Z=Z-Z.max(axis=1,keepdims=True); e=np.exp(Z); return e/e.sum(axis=1,keepdims=True)
def forward(X,W1,b1,W2,b2):
    z1=X@W1+b1; a1=np.maximum(z1,0); z2=a1@W2+b2; return softmax(z2),(X,z1,a1,z2)
def loss_fn(yh,Y): return -np.sum(Y*np.log(yh+1e-12))/Y.shape[0]
def backward(c,Y,W2):
    X,z1,a1,z2=c; N=X.shape[0]
    dz2=(softmax(z2)-Y)/N; dW2=a1.T@dz2; db2=dz2.sum(0)
    dz1=(dz2@W2.T)*(z1>0); return X.T@dz1, dz1.sum(0), dW2, db2
for eta in [0.01,0.4,5.0]:
    r=np.random.default_rng(0)
    W1=r.normal(0,np.sqrt(2/D),(D,H)); b1=np.zeros(H)
    W2=r.normal(0,np.sqrt(2/H),(H,C)); b2=np.zeros(C)
    for _ in range(300):
        yh,c=forward(X,W1,b1,W2,b2); dW1,db1,dW2,db2=backward(c,Y,W2)
        W1-=eta*dW1; b1-=eta*db1; W2-=eta*dW2; b2-=eta*db2
    yh,_=forward(X,W1,b1,W2,b2)
    print(f"eta={eta:<5} loss={loss_fn(yh,Y):.4f} acc={(yh.argmax(1)==y).mean():.3f}")

import numpy as np
z = np.array([[1.0, 2.0, 1000.0]])
e_bad = np.exp(z)
print("unstable:", e_bad / e_bad.sum(axis=1, keepdims=True))
zs = z - z.max(axis=1, keepdims=True)
e = np.exp(zs)
print("stable:  ", e / e.sum(axis=1, keepdims=True))

import numpy as np
rng = np.random.default_rng(1)
N, C = 8, 3
dz2 = rng.normal(size=(N, C)) / N
print("correct db (axis=0):", dz2.sum(axis=0).round(4), "shape", dz2.sum(axis=0).shape)
print("wrong   db (scalar):", round(float(dz2.sum()), 4))

动手练习

1. 手推 + 复核。 用第 5 节的同一组数字，但把真实标签改成第 1 类（\(y=[0,1]^\top\)），手推出 \(\delta_2,\ \partial L/\partial W_2,\ \partial L/\partial W_1\)，再用下面骨架核对。进阶：再换一组让 \(z_1\) 全为正的输入（例如 \(x=[2,1]^\top\)）重推一遍，体会当所有隐藏单元都被激活时，各梯度项都非零、外积是"满"的——不像例题里那样退化得稀疏。

   import numpy as np
   x=np.array([1.,-2.]); 
   W1=np.array([[.5,-.3],[.1,.2],[-.4,.6]]); b1=np.array([0.,.1,-.2])
   W2=np.array([[.2,-.1,.4],[.3,.5,-.2]]); b2=np.array([0.,.1])
   y=np.array([0.,1.])
   z1=W1@x+b1; a1=np.maximum(z1,0); z2=W2@a1+b2
   e=np.exp(z2-z2.max()); yh=e/e.sum()
   dz2=yh-y
   print("dz2",dz2)
   print("dW2",np.outer(dz2,a1))
   da1=W2.T@dz2; dz1=da1*(z1>0)
   print("dW1",np.outer(dz1,x))
   

2. 去掉稳定化看溢出。 把第 8 节 softmax 里的 Z - Z.max(...) 删掉，并把某个输入特征乘以 1000（X = X*1000），运行看是否出现 nan；再加回稳定化，确认恢复正常。
3. 偏置约定实验。 把 backward 里的 db2 = dz2.sum(axis=0) 误写成 db2 = dz2.sum()（求成标量），重新跑 grad_check()，观察相对误差从约 5e-9 跳到多大——亲手见证"偏置要沿样本维求和"。
4. 学习率扫描。 写个循环对 eta in [0.01, 0.1, 0.4, 1.5, 5.0] 各训练 300 步，打印最终 loss 与准确率，找出最好的 eta，并解释太大/太小各发生了什么。
5. 给 micrograd 加新运算。 给第 9 节的 Value 类加一个 tanh 方法（前向 np.tanh(x)，反向局部导数 1 - tanh(x)**2），用 backward() 验证它对一个标量输入的梯度，与中心差分对比。

掌握自检

• 我能不看讲义，从 \(\delta_2=\hat y-y\) 起手，独立推出四个梯度，并说出每个矩阵乘法里转置和左右摆放的理由（靠形状自检）。
• 我能解释 log-sum-exp 稳定化为什么不改变 softmax 结果，以及为什么要按行减最大值。
• 我能说清 batch 维下"权重梯度对 \(N\) 求平均、偏置沿样本维求和"的来历，知道偏置的求和与平均只差那个已被 \((\hat Y-Y)/N\) 提前吸收的 \(N\)，并知道漏掉除 \(N\) 会怎样。
• 我能写出中心差分公式，并知道梯度检查的相对误差应当 \(<10^{-7}\)，否则代码有 bug。
• 我能用一句话说明 loss.backward() 的本质（计算图 + 拓扑逆序 + 局部链式法则 + 梯度累加），并解释 ReLU 掩码如何引起梯度消失、\(W^\top\) 连乘如何引起梯度消失或爆炸。

到这里，模块 1（线性代数与微积分）全部完成——你已经备齐了驱动一切模型的数学引擎。下一模块我们换一套语言：进入概率、信息论与最优化（含强化学习速成），去回答"模型该有多自信""信息怎么量化""梯度下降之外还有哪些更聪明的优化器"。