神经网络矩阵微积分：从 y = Wx + b 推导 MSE 梯度

线性层的梯度看起来只有两行公式，实际调试时却容易出错：偏置少一个维度，程序可能仍然运行；损失多除一次输出数，梯度方向不变但幅度变了；把有限差分步长不断缩小，误差反而可能增大。本文用一组固定数字，把这些问题逐个算出来。

2026-09-08 修订：正文代码已与手算使用同一组输入，增加偏置和输入梯度检查，并区分精确分数计算与浮点近似。旧版未设种子的随机示例和另一组二维数字的动画，不再作为本次结果的证据。原下载文件保持不变。

顶部运行说明和下方程序均使用本次线性层梯度实验包。旧 deep-learning-math-lab 作为历史参考保留，不要混用两个目录的入口脚本。

一、先固定输入和损失定义

先只看一个样本：输入 x 为 3×1，权重 W 为 2×3，偏置 b 和目标 y 都为 2×1。这里的损失是平方误差和的一半，不是不加说明地调用一个名为 MSE 的函数。

W = [[ 0.2, -0.4,  0.1],
     [ 0.7,  0.3, -0.2]]
b = [[ 0.05], [-0.10]]
x = [[ 1.50], [-2.00], [0.50]]
y = [[ 1.00], [-1.00]]

y_hat = W x + b = [[1.20], [0.25]]
e = y_hat - y   = [[0.20], [1.25]]
L = 0.5 * sum(e**2) = 641/800 = 0.80125

例如第一个输出为 0.2×1.5 + (-0.4)×(-2) + 0.1×0.5 + 0.05 = 1.2。每个小数都来自上面的输入，而不是先假设一个误差，再让代码随机生成另一组参数。

二、从标量微分读出三个梯度

不必先构造一个很大的 Jacobian。把误差写成分量，再观察每个微小变化前的系数：

e_i = sum_j W_ij x_j + b_i - y_i
dL  = sum_i e_i de_i
    = sum_i e_i (sum_j x_j dW_ij + sum_j W_ij dx_j + db_i)

dL/dW_ij = e_i x_j
dL/db_i  = e_i
dL/dx_j  = sum_i W_ij e_i

整理成矩阵形式就是 dW = e x^T、db = e、dx = W^T e。前两个用于更新本层参数，最后一个把梯度传给前一层。梯度与被求导变量同形状是必要检查，但形状正确不代表数值正确。

dW = [[ 0.300, -0.400,  0.100],
      [ 1.875, -2.500,  0.625]]
db = [[0.200], [1.250]]
dx = [[0.915], [0.295], [-0.230]]

以 W[1,0] 为例，增加它会让第二个输出增加 1.5×dW[1,0]。当前第二个输出误差为 1.25，所以损失的一阶变化系数为 1.25×1.5=1.875。这同时说明了外积中行和列各自对应什么。

三、批量计算时到底除以什么

下面的程序约定样本按列排列：X=(D,B)，Y=(M,B)，偏置仍为 (M,1)。其中 B 是样本数，M 是输出数。目标函数平均样本，不平均输出：

E  = W X + b - Y
L  = sum(E**2) / (2 B)
G  = E / B
dW = G X^T
db = sum(G, axis=1, keepdims=True)
dX = W^T G

偏置对每一个样本都产生贡献，因此必须沿样本轴求和；keepdims=True 保留列向量形状，但不能替你选对求和轴。

为暴露“两输出时碰巧相同”的问题，复核包另有下面这个三输出、两样本的固定对照：

W3 = [[0.2,-0.4,0.1], [0.7,0.3,-0.2], [-0.3,0.5,0.2]]
b3 = [[0.05], [-0.1], [0.1]]
X3 = [[1.5,0.25], [-2.0,1.0], [0.5,-1.5]]
Y3 = [[1.0,-0.5], [-1.0,0.25], [0.2,-0.4]]

half_sample_mean = 1.0696875
element_mean_MSE = 0.713125
element_mean_MSE / half_sample_mean = 2/3

如果迁移到框架后梯度整体差一个固定倍数，先核对 loss reduction；不要立即靠调整学习率掩盖目标函数已经变化的事实。

四、与手算一致的完整 NumPy 程序

下面就是下载包中的 linear_example.py，不需要随机数或外部数据。默认用 float64；它检查 W、b、X 的全部 11 个梯度坐标，并在独立副本上做正负扰动，避免检查过程中改坏原输入。Y 在这个例子中视为固定目标。

"""A fixed linear-layer example; samples occupy columns, loss averages samples."""
import json
import numpy as np


def fixture(dtype=np.float64):
    return tuple(np.array(a, dtype=dtype) for a in (
        [[0.2, -0.4, 0.1], [0.7, 0.3, -0.2]],
        [[0.05], [-0.1]], [[1.5], [-2.0], [0.5]], [[1.0], [-1.0]]))


def evaluate(W, b, X, Y):
    arrays = (W, b, X, Y)
    if any(not isinstance(a, np.ndarray) or a.ndim != 2 for a in arrays):
        raise ValueError("Use two-dimensional NumPy arrays")
    if W.dtype not in (np.dtype('float32'), np.dtype('float64')) or any(a.dtype != W.dtype for a in arrays):
        raise ValueError("Use one common float32 or float64 dtype")
    outputs, inputs = W.shape
    batch = X.shape[1]
    if min(outputs, inputs, batch) == 0 or X.shape[0] != inputs or b.shape != (outputs, 1) or Y.shape != (outputs, batch):
        raise ValueError("Expected W=(M,D), b=(M,1), X=(D,B), Y=(M,B)")
    if any(not np.isfinite(a).all() for a in arrays):
        raise ValueError("Inputs must be finite")
    with np.errstate(over='raise', invalid='raise', divide='raise'):
        error = W @ X + b - Y
        loss = float(np.sum(error * error) * W.dtype.type(0.5 / batch))
        G = error / batch
        gradients = (G @ X.T, G.sum(axis=1, keepdims=True), W.T @ G)
    return loss, gradients


def finite_difference(arrays, slot, h=1e-5):
    if slot not in (0, 1, 2) or not np.isscalar(h) or not np.isfinite(h) or h <= 0:
        raise ValueError("Choose W/b/X (0/1/2) and a finite positive step")
    evaluate(*arrays)
    numerical = np.empty(arrays[slot].shape, dtype=np.float64)
    for index in np.ndindex(numerical.shape):
        plus = [a.copy() for a in arrays]
        minus = [a.copy() for a in arrays]
        plus[slot][index] += h
        minus[slot][index] -= h
        numerical[index] = (evaluate(*plus)[0] - evaluate(*minus)[0]) / (2 * h)
    return numerical


def main():
    W, b, X, Y = arrays = fixture()
    loss, gradients = evaluate(*arrays)
    checks = {}
    for slot, name in enumerate(('W', 'b', 'X')):
        numerical = finite_difference(arrays, slot)
        np.testing.assert_allclose(gradients[slot], numerical, atol=1e-8, rtol=1e-8)
        checks[name] = float(np.max(np.abs(gradients[slot] - numerical)))
    next_loss, _ = evaluate(W - 0.1 * gradients[0], b - 0.1 * gradients[1], X, Y)
    print(json.dumps({'loss': loss, 'error': (W @ X + b - Y).tolist(),
        'dW': gradients[0].tolist(), 'db': gradients[1].tolist(), 'dX': gradients[2].tolist(),
        'max_abs_difference': checks, 'loss_after_one_update': next_loss}, indent=2))
    print('LINEAR_EXAMPLE_OK')


if __name__ == '__main__':
    main()

输出中的 0.8012499999999998 与手算的 0.80125 之间是二进制浮点舍入差异。程序使用 atol=1e-8, rtol=1e-8 比较 float64 梯度，而不是要求打印出来的每一位都等于十进制手算。

最后一行更新只作用于 W 和 b，X、Y 固定。学习率为 0.1 时，新损失是 0.050078125。这组单样本中，更新后的误差为原来的 0.25 倍，因而损失为原来的 1/16；它不是任意网络或任意学习率都会下降的承诺。

五、哪些结果已经实际核对

2026-09-08 在 arm64 macOS、Python 3.13.9、NumPy 2.3.5 上运行了三个固定样例。除中心差分外，还用 Python Fraction 和标量循环独立计算损失与导数，不复用 NumPy 矩阵乘法作为参考。

共检查 46 个梯度坐标。分数计算还对每个坐标使用 1/10 和 1/100000 两种精确步长做中心差分，得到 92 次精确相等检查；NumPy 解析梯度与分数参考的最大绝对差约为 4.44e-16。这些检查覆盖指定样例，不是对所有输入的正确性证明。

完整梯度检查 CSV保留浮点差值，精确导数 CSV保留分数。旧 CSV 只记录六个权重坐标，差值按六位小数显示为 0.000000，偏置的数值检查栏为空；新包没有把这些空栏误称为已经验证。

六、步长越小，为什么可能越不准

中心差分使用 [L(θ+h)-L(θ-h)]/(2h)。一般光滑非线性函数存在截断误差和浮点舍入误差的权衡。但本例固定其他变量、只扰动一个坐标时，平方损失是该坐标的二次多项式；在精确算术里，中心差分的截断误差恰好为零。不能把别处常见的 U 形误差曲线硬套到它上面。

查看原尺寸误差图；下载全部 24 行扫描数据。

有限差分只是数值检查工具，不是绝对 ground truth。非光滑点、随机算子、舍入误差和共享存储都可能干扰解释。PyTorch gradcheck 文档也提醒其默认容差面向双精度，并说明非可微点可能导致检查失败；本轮没有安装或运行 PyTorch，不能把 NumPy 检查称为框架交叉验证。

七、广播错误要用失败输入检查

把本例的 b 从 (2,1) 压成 (2,)，表达式 W @ X + b.ravel() - Y 仍能执行，却会产生 (2,2) 的误差数组。损失从正确的 0.80125 变成 1.7825。这是一个可以直接复现的语义错误，不需要夸张成“必然显存溢出”。

NumPy 广播规则从右侧维度开始比较，维度相同或其中一个为 1 才能匹配。因此，(64,1)+(64,) 会形成 64×64 结果；一个 float32 结果载荷为 16,384 字节，即 16 KiB。这个小例子本身不足以证明 GPU OOM，大规模内存后果需要另测。

失败输入用于确认检查器能发现具体错误，不是把“通过测试的数量”当作正确性证明。微分推导、形状检查、独立参考与数值误差分析各自解决不同的问题。

八、下载与继续验证

下载线性层梯度实验包，在解压后的目录运行：

python3 -m venv .venv
.venv/bin/python -m pip install -r requirements.txt
.venv/bin/python linear_example.py
.venv/bin/python audit_matrix.py --output reproduced

两个成功标记分别是 LINEAR_EXAMPLE_OK 和 MATRIX_GRADIENT_AUDIT_OK。结果 JSON记录环境、源码哈希、全部样例、失败对照和边界。绘图是可选步骤，命令和依赖在包内 README 中。

本轮没有进行 GPU 吞吐、显存带宽、分布式训练或混合精度性能测试。有限差分中的 h 与混合精度训练的 loss scaling 不是同一机制，不能用调大 loss scale 来替代检查扰动是否被浮点表示保留。

下一步可阅读两层 MLP 的反向传播与梯度检查，观察这里的 dX 如何变成前一层的上游梯度；再结合梯度下降、Momentum 与 Adam 的固定实验，区分“梯度计算正确”和“更新过程按预期收敛”。