反向传播计算图：两层 MLP 的前向、局部梯度和反向传播

反向传播（Backpropagation）常被渲染成深度学习的神秘引擎，但从本质上讲，它只是反向模式自动微分（Reverse-mode Automatic Differentiation）在计算图上的应用。它并不是神经网络专属的魔法，而是一种高度优化、基于微积分链式法则的程序化求导方法。

本文使用一个隐藏层、两个带参数的仿射层：x @ W1+b1 → ReLU → h @ W2+b2 → 平均交叉熵。输入按行组成 batch，先推导梯度，再执行完整的 NumPy 示例，最后逐个参数用有限差分核对。这里的随机数据是数值测试夹具，不用于证明模型具备分类能力。

实验版本说明：页首“代码运行说明”中的 src/mlp_backprop.py 属于旧系列包，使用单样本列向量。要复现本文的 32 样本 batch 与 1029 参数检查，请使用下方 Neural Gradient Lab 的结果、下载和运行命令；两套输入不同，数值不能直接对照。

一、揭开计算图的面纱

为了系统地计算导数，我们把复杂的神经网络拆解成由基础算子构成的有向无环图（DAG）。每个节点代表一个简单的数学操作（如矩阵乘法或 ReLU），边代表张量（Tensor）的流动。至关重要的是，每个节点必须具备两种能力：在前向时计算输出，在反向时计算局部的向量雅可比乘积（Vector-Jacobian Product, VJP）。

graph TD
    x["输入 x"] --> z1["z1 = x @ W1 + b1"]
    W1["权重 W1"] --> z1
    b1["偏置 b1"] --> z1
    z1 --> h["h = ReLU(z1)"]
    h --> logits["logits = h @ W2 + b2"]
    W2["权重 W2"] --> logits
    b2["偏置 b2"] --> logits
    logits --> p["p = Softmax(logits)"]
    p --> L["Loss = CrossEntropy(p, target)"]
    target["目标 y"] --> L

    classDef fwd fill:#e1f5fe,stroke:#039be5,stroke-width:2px;
    classDef param fill:#fce4ec,stroke:#d81b60,stroke-width:2px;
    class z1,h,logits,p,L fwd;
    class W1,b1,W2,b2 param;

反向传播需要访问前向中间值，例如 x、z1 和 h。可以保存它们，也可以在反向阶段重新计算部分激活。训练通常还需要梯度和优化器状态，但是否受显存容量、带宽或计算量限制，必须通过具体模型的性能分析判断。

二、Softmax Cross-Entropy 的关键简化

理论上，你可以先求交叉熵损失对 Softmax 概率的雅可比矩阵，再乘上 Softmax 对 logits 的雅可比矩阵。但在工程实践中，显式地这样做既会引发数值灾难，又会浪费大量算力。

对一个样本、整数类别标签、无类别权重、无 label smoothing 的交叉熵，logits 梯度如下。若对 N 个样本取平均损失，每行还要除以 N；若求和则不除。加权或忽略标签时，应重新核对归一化规则，不能机械套用。

dL/dlogits = p - one_hot(target)

类别从 0 开始编号时，目标为类 1、概率为 [0.1,0.7,0.2]，单样本梯度是 [0.1,-0.3,0.2]。注意 PyTorch CrossEntropyLoss 接收未归一化 logits，不应先做 softmax 再传入；概率只是本文推导的中间量。

三、两层 MLP 的反向公式推导

一旦我们拿到了损失对 logits 的梯度（dlogits），就可以利用链式法则将其向后传播。注意，我们在实际计算中从不实例化完整的雅可比矩阵，而是利用矩阵转置来高效计算 Vector-Jacobian Products。

# 第二层梯度
dW2 = h^T dlogits
db2 = sum(dlogits, axis=0)
dh  = dlogits W2^T

# 第一层梯度
dz1 = dh * ReLU'(z1)  # 逐元素相乘
dW1 = x^T dz1
db1 = sum(dz1, axis=0)

梯度的 shape 必须与对应参数完全一致，偏置梯度要沿 batch 维归约。范数可以提示异常尺度，但单次范数既不能证明梯度正确，也不能排除整个训练过程的消失或爆炸问题。下面给出固定输入下的实测范数和独立数值对照。

四、可直接运行的 NumPy 实现

下面代码固定 seed=42，输入为 32×10，隐藏层 64 单元，输出 5 类；NumPy 默认生成 float64 数据。它计算平均交叉熵、四组参数梯度和一次更新后的损失，代码与下载包的 mlp_example.py 一致。减最大值用于稳定指数运算，损失从 log-sum-exp 形式计算，避免先得到零概率再求 log。这是有明确输入约束的教学实现，不是通用训练库。

"""One-hidden-layer ReLU MLP, row-batch convention, mean unweighted CE."""

import numpy as np


def fixture(seed=42):
    rng = np.random.default_rng(seed)
    x = rng.standard_normal((32, 10))
    params = {
        "W1": rng.standard_normal((10, 64)) * 0.1,
        "b1": np.zeros((1, 64)),
        "W2": rng.standard_normal((64, 5)) * 0.1,
        "b2": np.zeros((1, 5)),
    }
    targets = rng.integers(0, 5, size=32)
    return x, targets, params


def cross_entropy(logits, targets):
    if logits.ndim != 2 or min(logits.shape) < 1:
        raise ValueError("Expected nonempty (batch, classes) logits")
    if targets.shape != (logits.shape[0],) or targets.dtype.kind not in "iu":
        raise ValueError("Expected one integer class index per row")
    if not np.isfinite(logits).all() or np.any(targets < 0) or np.any(targets >= logits.shape[1]):
        raise ValueError("Nonfinite logits or out-of-range target")
    shifted = logits - logits.max(axis=1, keepdims=True)
    exps = np.exp(shifted)
    normalizer = exps.sum(axis=1, keepdims=True)
    probs = exps / normalizer
    # Compute the loss from logits, not log(probs), which can underflow to log(0).
    loss = np.mean(np.log(normalizer[:, 0]) - shifted[np.arange(len(targets)), targets])
    dlogits = probs.copy()
    dlogits[np.arange(len(targets)), targets] -= 1
    dlogits /= len(targets)
    return float(loss), probs, dlogits


def loss_and_grads(x, targets, params):
    W1, b1, W2, b2 = (params[k] for k in ("W1", "b1", "W2", "b2"))
    if x.ndim != 2 or W1.ndim != 2 or W2.ndim != 2:
        raise ValueError("Expected matrices for x, W1, W2")
    if x.shape[1] != W1.shape[0] or W1.shape[1] != W2.shape[0]:
        raise ValueError("Incompatible layer shapes")
    if b1.shape != (1, W1.shape[1]) or b2.shape != (1, W2.shape[1]):
        raise ValueError("Biases must have shape (1, units)")
    if not all(np.isfinite(a).all() for a in (x, W1, b1, W2, b2)):
        raise ValueError("Nonfinite input or parameter")
    z1 = x @ W1 + b1
    h = np.maximum(z1, 0)
    logits = h @ W2 + b2
    loss, probs, dlogits = cross_entropy(logits, targets)
    dh = dlogits @ W2.T
    dz1 = dh * (z1 > 0)  # Choose derivative 0 at the ReLU kink.
    grads = {
        "W1": x.T @ dz1,
        "b1": dz1.sum(axis=0, keepdims=True),
        "W2": h.T @ dlogits,
        "b2": dlogits.sum(axis=0, keepdims=True),
    }
    for key in params:
        if grads[key].shape != params[key].shape:
            raise AssertionError("Gradient shape mismatch: " + key)
    return loss, grads, {"z1": z1, "h": h, "logits": logits, "probs": probs}


def step(params, grads, learning_rate):
    return {key: value - learning_rate * grads[key] for key, value in params.items()}


if __name__ == "__main__":
    x, targets, params = fixture()
    loss, grads, cache = loss_and_grads(x, targets, params)
    next_loss = loss_and_grads(x, targets, step(params, grads, 0.1))[0]
    print(f"seed=42 dtype={x.dtype} batch=32, features=10, hidden=64, classes=5")
    print(f"mean CE: {loss:.12f} -> {next_loss:.12f}")
    for key, grad in grads.items():
        print(f"{key}: shape={grad.shape}, norm={np.linalg.norm(grad):.12f}")

五、从数值失败定位问题

以下结论来自本页可运行的测试与公开文档，不依赖个人履历或未经展示的 CUDA 性能经验。

• 偏置广播。本代码要求 b1 的 shape 为 (1,64)，所以 db1 使用 dz1.sum(axis=0, keepdims=True)。如果参数用一维 (64,)，相应梯度也可用一维；重点是匹配参数，而非无条件要求 keepdims=True。测试会拒绝错误的 (32,64) batch 形状偏置。
• 差分检查。使用中心差分 (L(p+eps)-L(p-eps))/(2eps)，每次只改变一个坐标并恢复。误差还取决于步长、精度和 ReLU 是否跨越 0，不能看到大于 1e-4 就直接断言代码错误。PyTorch gradcheck 文档也明确提醒双精度和不可微点问题。
• 稳定损失。对 logits=[1000,-1000,0]、目标类 1，减最大值的 softmax 仍会把目标概率下溢为 0。直接 -log(probability) 得到 Inf；本代码的损失是有限的 2000，梯度为 [1,-1,0]。这不代表任意极大有限输入都不会造成中间运算溢出。
• 激活重算。Gradient checkpointing 用重算换取较少保存的激活，并不证明所有训练瓶颈都在带宽；是否划算要测实际时间和峰值内存。

六、如何观看演示动画

在观看动画时，不要只盯着箭头的方向。仔细观察每个节点在前向时保存了什么值，以及为什么反向传播时必须复用这些值。在下一篇文章中，我们将研究这些梯度是如何驱动参数更新的，以及不同的优化器为什么会走出完全不同的参数轨迹。

七、反向传播验证矩阵

复现这篇文章时，建议把一次前向和一次反向拆成可审计步骤。下面的表格把“公式是否正确”转化为可观测证据，避免只看到 loss 下降就误以为反向传播一定正确。

八、逐个参数检查的实际结果

运行环境为 Python 3.13.9、NumPy 2.3.5，CPU、float64、seed=42。初始平均 CE 为 1.647801573097；所有梯度取自同一参数快照，学习率 0.1 更新后为 1.628354091165。随机标签没有训练/测试划分，这不是模型准确率或泛化实验。

总共检查 1029 个参数坐标，而不是只抽查权重、漏掉偏置。每个坐标的通过条件为 abs(a-n) ≤ 1e-7 + 1e-5 × max(abs(a),abs(n))，a 为解析梯度，n 为数值梯度。容差属于这个双精度测试，不是所有模型通用的合格线。

前两行存在跨 ReLU 拐点的探针，不能当作光滑点局部导数的验证；后两行没有跨点，均为 1029/1029 坐标通过。表格对显示值做了四舍五入，完整精度及每个坐标的结果保留在下载包的 CSV 和 JSON 中。

最小 |z1| 约为 1.34039e-5，但扰动 W1 的一个坐标会让 z1 改变 x_j × eps，不一定只改变 eps。因此 1e-5 也可能越过拐点。更小步长又会放大减法舍入误差，表中 1e-7 的误差反而比 1e-6 大。在恰好 z=0 时，ReLU 的中心差分为 0.5，本实现选择的次梯度为 0；这里不能用一个固定阈值判谁“算错”。

为检查测试真的能发现错误，实验分别把反向梯度整体反号、故意漏除 batch 大小，两种错误都被拒绝；把同一 batch 重复一次，正确的平均损失和梯度保持一致。额外检查拒绝了 6 类非法标签、数值或形状输入。这些都是具体夹具的通过证据，不是对所有网络的形式化证明。

下载与复现

下载反向传播梯度检查实验包。其中 reference/mlp-arrays.npz 保存本次输入、标签、参数、前向中间量和解析梯度；四份 CSV 列出每个差分坐标，reference/audit.json 保存摘要与源文件哈希。旧系列包的单样本列向量示意结果仍保留，但不能与本文 32 样本 batch 的数值混用。

python3 -m venv .venv
source .venv/bin/activate
python -m pip install -r requirements.txt
python mlp_example.py
python audit_gradients.py --output run-results

成功结束会打印 NEURAL_GRADIENT_AUDIT_PASSED。审计本身无需网络、训练数据或 GPU；安装依赖可能访问软件包源。不同 BLAS、CPU 或 NumPy 版本的末位浮点数可能不同，应按误差容限比较，不要求跨机器 JSON 字节一致。本轮没有运行框架 autograd 或 CUDA 对照。

需要更小的手算例子，可先读 神经网络基础中的五参数更新与初始化反例。