矩阵求导

上一课我们用计算图，把损失对每个标量的偏导一个一个地反向传播出来。那很扎实，但当一层网络有成千上万个权重时，逐元素写偏导既慢又容易错。这一课我们把同样的事情向量化：直接对整个向量、整个矩阵求导，让一行公式 \(\partial L/\partial W=\delta\,x^\top\) 替代成千上万行链式法则。

1. 先把"布局"这件事一次说清

"矩阵求导"之所以让初学者头疼，90% 的痛苦来自一个根本没说清的约定：导数排成什么形状。同一个 \(\partial(a^\top x)/\partial x\)，有的书得到行向量，有的书得到列向量——它们都没错，只是布局约定（layout convention）不同。

有两种约定：

• 分子布局（numerator layout）：结果的形状跟随分子。\(\partial y/\partial x\) 中，\(y\) 决定行、\(x\) 决定列。
• 分母布局（denominator layout）：结果的形状跟随分母。这正是深度学习的惯例，因为我们最关心 \(\partial L/\partial W\)（标量损失对参数），希望它和参数同形，这样才能直接做梯度下降 \(W \leftarrow W-\eta\,\partial L/\partial W\)。

1.1 四类导数的形状总览

把"谁对谁求导"分成四类。设标量 \(a\)、列向量 \(x\in\mathbb{R}^n\)、\(y\in\mathbb{R}^m\)、矩阵 \(W\in\mathbb{R}^{m\times n}\)。在分母布局下：

注意第二行第二列：\(\partial y/\partial x\)。上一课我们把雅可比（Jacobian）\(J\) 定义成 \(m\times n\)（\(J_{ij}=\partial y_i/\partial x_j\)，分子布局）。在分母布局下，\(\partial y/\partial x\) 是它的转置 \(J^\top\in\mathbb{R}^{n\times m}\)。这正是上面"总原则"里说的唯一例外——也是两套约定唯一容易"打架"的地方，下面会专门点出。

2. 从第一性原理推导四个必备公式

我们不背表，而是从定义推。所有推导只用一件事：标量对向量的导数，第 \(k\) 个分量就是对 \(x_k\) 的偏导。

2.1 线性型 \(\partial(a^\top x)/\partial x = a\)

设 \(f(x)=a^\top x=\sum_{k=1}^n a_k x_k\)，这里 \(a\) 是常向量。对第 \(i\) 个分量求偏导：

把所有 \(i\) 堆成一个列向量（分母布局），得到

Shape check：分母是 \(x\in\mathbb{R}^n\)，结果应是 \(\mathbb{R}^n\) 列向量——而 \(a\in\mathbb{R}^n\)，对上了。这就是标量对向量求导最基本的"砖块"，相当于一元微积分里的 \((ax)'=a\)。

2.2 矩阵乘向量 \(\partial(Wx)/\partial x\)

这里 \(y=Wx\) 是向量对向量——也就是 §1 里说的那个唯一例外。先别急着记结论：实践中我们几乎从不单独用这个导数，谁该转置最终由"形状必须对上"自动决定（见 §4），所以下面的方框只是为了和上一课的雅可比对照，不必背。

第 \(i\) 个输出 \(y_i=\sum_k W_{ik}x_k\)，所以 \(\partial y_i/\partial x_j=W_{ij}\)。雅可比（分子布局，\(m\times n\)）正是 \(W\) 本身：

而在分母布局下，\(\partial y/\partial x\) 是它的转置：

Shape check：分母 \(x\in\mathbb{R}^n\) 占行（\(n\)），分子 \(y\in\mathbb{R}^m\) 占列（\(m\)），所以是 \(n\times m\)——而 \(W^\top\) 正是 \(n\times m\)。

2.3 二次型 \(\partial(x^\top A x)/\partial x=(A+A^\top)x\)

设 \(f(x)=x^\top A x=\sum_{i}\sum_{j} x_i A_{ij} x_j\)，\(A\) 是常矩阵。要对 \(x_k\) 求偏导，注意 \(x_k\) 出现在两处：作为外层下标 \(i=k\)，以及作为内层下标 \(j=k\)。用乘积法则：

第一个和是 \((Ax)_k\)；第二个和是 \(\sum_i A_{ik}x_i=(A^\top x)_k\)。合起来：

当 \(A\) 对称（\(A=A^\top\)，机器学习里协方差矩阵、海森矩阵都对称）时简化为

Shape check：\((A+A^\top)\in\mathbb{R}^{n\times n}\) 乘 \(x\in\mathbb{R}^n\) 得 \(\mathbb{R}^n\)，与分母 \(x\) 同形。

2.4 平方范数 \(\partial\|x\|^2/\partial x=2x\)

这是 §2.3 的特例：\(\|x\|^2=x^\top x=x^\top I x\)，取 \(A=I\)（对称），立刻得

它是一元 \((x^2)'=2x\) 的向量版本，也是 L2 正则、最小二乘损失里最常见的梯度。

3. 微分法：一招通杀的 \(dL=\operatorname{tr}(G^\top dX)\Rightarrow \partial L/\partial X=G\)

逐元素求偏导对复杂表达式会越来越累。有一种更优雅、更不易错的工具——微分法（differential method）。它的核心要用到迹（trace），所以先花一句话把迹讲清楚。

有了迹，微分法的核心是一条"识别定理"：

为什么成立？把迹展开：\(\operatorname{tr}(G^\top dX)=\sum_{i,j}G_{ij}\,dX_{ij}\)。而由全微分定义 \(dL=\sum_{i,j}\dfrac{\partial L}{\partial X_{ij}}dX_{ij}\)。两式对每个独立的 \(dX_{ij}\) 都成立，逐项比较系数即得 \(G_{ij}=\partial L/\partial X_{ij}\)，也就是 \(G=\partial L/\partial X\)。向量是 \(n\times 1\) 矩阵的特例，公式照用：\(dL=g^\top dx\Rightarrow \partial L/\partial x=g\)。

三条够用的运算法则（\(d\) 像普通微分一样线性，且服从乘积法则）：

• \(d(AX)=A\,dX\)，\(d(XA)=(dX)A\)（常矩阵提到外面）；
• \(d(\operatorname{tr}(X))=\operatorname{tr}(dX)\)，且迹有循环性 \(\operatorname{tr}(ABC)=\operatorname{tr}(BCA)\)；
• \(\operatorname{tr}(A^\top)=\operatorname{tr}(A)\)，\(\operatorname{tr}(A^\top B)=\sum_{i,j}A_{ij}B_{ij}\)（这正是把"逐元素相乘再求和"写成迹）。

小试身手（用微分法重做平方范数）：\(L=\|x\|^2=x^\top x\)。两边微分，用乘积法则：

对照 \(dL=g^\top dx\)，立得 \(g=\partial L/\partial x=2x\)。和 §2.4 一致，但没碰任何下标。微分法的威力在矩阵情形（下一节）才真正显现。

4. 一层线性层的全部梯度

这是反向传播的"原子层"。设

后面接一个标量损失 \(L(y)\)。我们已经（在上一课或本层下游）拿到 上游梯度

目标：求 \(\partial L/\partial W,\ \partial L/\partial b,\ \partial L/\partial x\)。

4.1 用微分法一次推完

这里我们分别对每个变量求偏导：求 \(\partial L/\partial W\) 时，把 \(b\)、\(x\) 视作常数（即令 \(db=0,\ dx=0\)，注意 \(x\) 是该层的输入、\(b\) 是参数，它们并非恒为常量，只是在对 \(W\) 求偏导这一语境里暂时当常数）。于是 \(y=Wx+b\) 的微分只剩 \(dy=(dW)x\)。又 \(dL=\delta^\top dy\)（这是 \(\delta\) 的定义，标量对向量）。代入：

这是个标量，等于自己的迹（标量看成 \(1\times1\) 矩阵），再用循环性把 \(dW\) 凑到末尾：

对照识别公式 \(dL=\operatorname{tr}(G^\top dW)\)，读出 \(G=\delta x^\top\)：

同理，对 \(b\) 求偏导时把 \(W\)、\(x\) 当常数，\(dy=db\Rightarrow dL=\delta^\top db\)，得 \(\partial L/\partial b=\delta\)；对 \(x\) 求偏导时把 \(W\)、\(b\) 当常数，\(dy=W\,dx\Rightarrow dL=\delta^\top W\,dx=(W^\top\delta)^\top dx\)，得 \(\partial L/\partial x=W^\top\delta\)。汇总：

4.2 外积结构与 shape check

\(\partial L/\partial W=\delta x^\top\) 是一个外积（outer product）：列向量 \(\delta\in\mathbb{R}^{m\times 1}\) 乘行向量 \(x^\top\in\mathbb{R}^{1\times n}\)，结果 \(m\times n\)——恰好与 \(W\) 同形。这个"上游梯度 ⊗ 本层输入"的外积，是所有全连接层反向传播的统一模板。

再核对另外两个：\(\partial L/\partial b=\delta\in\mathbb{R}^m\) 与 \(b\) 同形 ✓；\(\partial L/\partial x=W^\top\delta\)，\(W^\top\) 是 \(n\times m\) 乘 \(\delta\) 的 \(m\times1\) 得 \(n\times1\)，与 \(x\) 同形 ✓。三个梯度的形状全部"自动对上"——这正是分母布局的好处。注意 §2.2 那个"该用 \(W\) 还是 \(W^\top\)"的纠结，在这里被 shape check 自动解决了：要让结果与 \(x\) 同形，就只能是 \(W^\top\delta\)。

5. softmax + 交叉熵：那个"好得不像话"的梯度

分类网络最后一层几乎总是 softmax 把打分（logits）\(z\in\mathbb{R}^K\) 变成概率，再用交叉熵 cross-entropy 与真标签比对。神奇的是，两者复合后梯度极其简洁：\(\partial L/\partial z=\hat y-y\)。我们来推它——而且照本课一贯的做法，softmax 的雅可比也从定义亲手推，不直接抛结论。

记 softmax 输出 \(\hat y_k=\dfrac{e^{z_k}}{\sum_{j}e^{z_j}}\)，真标签是 one-hot 向量 \(y\)（正确类 \(c\) 处为 1，其余 0）。交叉熵损失

第一步：softmax 的偏导。设 \(S=\sum_j e^{z_j}\)，则 \(\hat y_k=e^{z_k}/S\)。用商法则对 \(z_i\) 求导（注意 \(\partial S/\partial z_i=e^{z_i}\)，且分子 \(e^{z_k}\) 只在 \(k=i\) 时才依赖 \(z_i\)），分两种情况：

两种情况可以合写成一条（用 Kronecker 记号 \(\delta_{ki}\)：\(k=i\) 取 1，否则 0）：

第二步：链式法则。先用 \(L=-\sum_k y_k\log\hat y_k\)，\(\partial L/\partial\hat y_k=-y_k/\hat y_k\)，于是

拆成两个和：\(-\sum_k y_k\delta_{ki}=-y_i\)（只有 \(k=i\) 那项留下）；而第二项里 \(\hat y_i\) 不含求和指标 \(k\)，可提到求和号外面，得 \(\sum_k y_k\hat y_i=\hat y_i\sum_k y_k=\hat y_i\)（因为 one-hot 的 \(\sum_k y_k=1\)）。所以

堆成向量：

Shape check：\(\hat y,y\in\mathbb{R}^K\)，差还是 \(\mathbb{R}^K\)，与分母 \(z\) 同形 ✓。在实现里，这个 \(\hat y-y\) 就是喂给上一层（输出线性层）的 \(\delta\)，直接套 §4 的三大梯度。

6. 例题

7. 代码：用数值梯度验证我们的公式

口说无凭。验证矩阵求导公式最可靠的办法是有限差分数值梯度（finite-difference numerical gradient）：用 \(\dfrac{\partial L}{\partial \theta_i}\approx\dfrac{L(\theta+\epsilon e_i)-L(\theta-\epsilon e_i)}{2\epsilon}\)（中心差分，误差 \(O(\epsilon^2)\)）逐个分量逼近，再和解析公式比相对误差。下面三段一次性验证 §2.3、§4、§5 的全部公式，相对误差都应 \(<10^{-6}\)。

import numpy as np
rng = np.random.default_rng(0)
eps = 1e-6

def numgrad(f, theta):
    """对标量函数 f(theta) 在 theta 处做中心差分数值梯度（theta 任意形状）。"""
    g = np.zeros_like(theta, dtype=float)
    it = np.nditer(theta, flags=['multi_index'])
    for _ in it:
        idx = it.multi_index
        old = theta[idx]
        theta[idx] = old + eps; fp = f(theta)
        theta[idx] = old - eps; fm = f(theta)
        theta[idx] = old
        g[idx] = (fp - fm) / (2*eps)
    return g

def relerr(a, b):
    return np.linalg.norm(a-b) / (np.linalg.norm(a)+np.linalg.norm(b) + 1e-12)

# ---- 验证 1: d(x^T A x)/dx = (A + A^T) x （A 故意非对称）----
n = 4
A = rng.standard_normal((n, n))
x = rng.standard_normal(n)
f1 = lambda x: x @ A @ x
g_analytic = (A + A.T) @ x
g_numeric  = numgrad(f1, x.copy())
print("[1] x^T A x   相对误差 =", relerr(g_analytic, g_numeric))

# ---- 验证 2: 线性层 y=Wx+b, L=0.5||y-t||^2 ----
m, n = 3, 5
W = rng.standard_normal((m, n)); b = rng.standard_normal(m)
xv = rng.standard_normal(n);     t = rng.standard_normal(m)
def L2(W, b, xv):
    y = W @ xv + b
    return 0.5*np.sum((y - t)**2)
y = W @ xv + b
delta = (y - t)              # = dL/dy
dW = np.outer(delta, xv)     # = delta x^T  (m x n)
db = delta
dx = W.T @ delta             # (n,)
print("[2a] dL/dW    相对误差 =", relerr(dW, numgrad(lambda W: L2(W,b,xv), W.copy())))
print("[2b] dL/db    相对误差 =", relerr(db, numgrad(lambda b: L2(W,b,xv), b.copy())))
print("[2c] dL/dx    相对误差 =", relerr(dx, numgrad(lambda xv: L2(W,b,xv), xv.copy())))
print("     形状检查 dW.shape == W.shape :", dW.shape == W.shape)

# ---- 验证 3: softmax + 交叉熵, dL/dz = yhat - y ----
K = 5
z = rng.standard_normal(K)
y_oh = np.zeros(K); y_oh[2] = 1.0          # 正确类 = 第 2 类
def softmax(z):
    z = z - z.max()                         # 数值稳定
    e = np.exp(z); return e / e.sum()
def ce(z):
    p = softmax(z); return -np.sum(y_oh * np.log(p))
yhat = softmax(z)
g3 = yhat - y_oh
print("[3] softmax+CE 相对误差 =", relerr(g3, numgrad(ce, z.copy())))

运行后五个相对误差应当都在 \(10^{-9}\sim10^{-6}\) 量级，形状检查输出 True。这就是工程上验证任何反向传播实现是否正确的"金标准"——下一课写完两层网络，我们也会用同样的 numgrad 给它做体检。

调一调，观察现象

下面三个小实验都只用 numpy 和 print，几秒钟跑完。改一个数、看一个现象，把本课的公式从"看懂"变成"亲眼见过"。

import numpy as np
A = np.array([[1.,2.],[0.,3.]])   # 故意非对称
x = np.array([1.,2.])
print("正确 (A+A^T)x =", (A + A.T) @ x)   # [ 6. 14.]
print("偷懒 2Ax       =", 2 * A @ x)       # [10. 12.]  —— 不一样！
As = 0.5*(A + A.T)                          # 对称化
print("对称化后 2*As*x =", 2 * As @ x)      # [ 6. 14.]  —— 又对上了

import numpy as np
delta = np.array([0.5, -2.0])     # m=2 上游梯度
xv    = np.array([1., 0., -1.])   # n=3 本层输入
dW = np.outer(delta, xv)          # delta x^T
print("dW =\n", dW)
print("dW.shape =", dW.shape, " 应与 W 同形 (2,3)")
print("rank(dW) =", np.linalg.matrix_rank(dW))   # 永远是 1

import numpy as np
def softmax(z):
    z = z - z.max(); e = np.exp(z); return e/e.sum()
y = np.array([1.,0.,0.])               # 真标签=第0类
for z in ([2.,1.,0.], [5.,1.,0.]):
    g = softmax(np.array(z)) - y
    print("logits", z, "-> grad", np.round(g,3),
          " sum =", round(g.sum() + 0.0, 12))   # 两组都打印 0.0

动手练习

1. 对称化二次型。证明对任意 \(A\)，\(x^\top A x=x^\top A_{\text{sym}}x\)，其中 \(A_{\text{sym}}=\tfrac12(A+A^\top)\)。由此说明为何梯度公式里只出现 \(A+A^\top\)（"反对称部分对二次型无贡献"）。再用 numpy 取一个随机非对称 \(A\)，数值验证 \(x^\top A x=x^\top A_{\text{sym}}x\)。

   import numpy as np
   rng = np.random.default_rng(1)
   A = rng.standard_normal((4,4)); x = rng.standard_normal(4)
   As = 0.5*(A + A.T)
   print("两种二次型之差 =", x@A@x - x@As@x)   # 应≈0
   

2. L2 正则的梯度。设 \(L(W)=\tfrac{\lambda}{2}\|W\|_F^2=\tfrac{\lambda}{2}\sum_{ij}W_{ij}^2\)（Frobenius 范数平方）。手推 \(\partial L/\partial W\)，并解释它为什么导致"权重衰减 weight decay"。用 numgrad 验证你的解析式。（提示：答案是 \(\lambda W\)，与 \(W\) 同形。）
3. batch 版线性层。把 §4 推广到 mini-batch：\(Y=XW^\top+\mathbf{1}b^\top\)，\(X\in\mathbb{R}^{B\times n}\)（行=样本），其中 \(\mathbf{1}\in\mathbb{R}^B\) 是全 1 列向量、\(\mathbf{1}b^\top\in\mathbb{R}^{B\times m}\) 把偏置复制到每一行（即 numpy 的广播）。上游 \(\Delta=\partial L/\partial Y\in\mathbb{R}^{B\times m}\)。证明 \(\partial L/\partial W=\Delta^\top X\)、\(\partial L/\partial b=\Delta^\top\mathbf{1}\)（按样本求和）、\(\partial L/\partial X=\Delta W\)。用下面骨架做数值验证。

   import numpy as np
   rng = np.random.default_rng(2)
   B, n, m = 8, 5, 3
   X = rng.standard_normal((B, n)); W = rng.standard_normal((m, n)); b = rng.standard_normal(m)
   T = rng.standard_normal((B, m))
   ones = np.ones(B)               # 全 1 列向量 1
   def Lb(W, b, X):
       Y = X @ W.T + np.outer(ones, b)   # = X W^T + 1 b^T （广播写法等价于 X @ W.T + b）
       return 0.5*np.sum((Y - T)**2)
   Y = X @ W.T + np.outer(ones, b)
   Delta = (Y - T)                 # B x m
   dW = Delta.T @ X                # m x n  —— 你来核对维度
   db = Delta.T @ ones             # m      （= 按样本求和 Delta.sum(axis=0)）
   # TODO: 写一个 numgrad 验证 dW, db（可复用正文里的 numgrad）
   print("dW.shape =", dW.shape, " 应为", W.shape)
   

4. softmax 雅可比。不带交叉熵，直接写出 softmax 的雅可比 \(J_{ki}=\partial\hat y_k/\partial z_i=\hat y_k(\delta_{ki}-\hat y_i)\) 的矩阵形式 \(J=\operatorname{diag}(\hat y)-\hat y\hat y^\top\)。用 numpy 构造它，并验证 \(J\mathbf{1}=0\)（每行和为 0，因为概率归一）。
5. 从微分法推 \(\partial L/\partial x=W^\top\delta\)。仿照 §4.1，对 \(y=Wx\) 取 \(dy=W\,dx\)，代入 \(dL=\delta^\top dy\)，整理成 \(dL=g^\top dx\) 读出 \(g\)。写下每一步，确认 \(g=W^\top\delta\) 且形状为 \(\mathbb{R}^n\)。

掌握自检

• 我能一口说出本课的布局约定，并解释"为什么 \(\partial L/\partial W\) 必须与 \(W\) 同形"。
• 给我任意一个梯度表达式，我能立刻做 shape check，仅凭形状就判断它"不可能对"还是"可能对"。
• 我能从第一性原理（逐元素或微分法）独立推出 \(\partial(a^\top x)/\partial x\)、\(\partial(x^\top Ax)/\partial x\)、\(\partial\|x\|^2/\partial x\)，不靠背表。
• 我记得线性层三大梯度 \(\partial L/\partial W=\delta x^\top\)、\(\partial L/\partial b=\delta\)、\(\partial L/\partial x=W^\top\delta\)，并能说清那个外积为什么自动与 \(W\) 同形、为什么反向用 \(W^\top\)。
• 我能从 softmax 的商法则推导出雅可比 \(\hat y_k(\delta_{ki}-\hat y_i)\)，再得到 softmax+CE 的 \(\partial L/\partial z=\hat y-y\)，并解释它"预测越错梯度越大、分量和为 0"的直觉。
• 我会写中心差分 numgrad 去验证任意解析梯度，知道 \(\epsilon\) 不能取太小，也知道要避开 ReLU/绝对值/max 这类不可导点。

下一课，我们把第 1 模块的所有零件——向量、矩阵变换、链式法则、计算图、以及本课的矩阵求导——拼到一起，从零手写并训练一个两层神经网络，亲眼看着损失曲线降下来。