""" ============================================================================= 第 3 讲 · 代码 1：矩阵微积分 —— 雅可比矩阵、形状约定、一层神经网络的梯度 ============================================================================= 对应讲义：lectures/lecture03_neural_networks_backprop.md 第 4 节 对应作业：CS224N 作业 2 的书面部分（Neural Network Foundations / Tensor Derivatives） ----------------------------------------------------------------------------- 【这个文件要解决什么问题？】 ----------------------------------------------------------------------------- 第 1、2 讲我们一个一个地推导了标量、向量的梯度。神经网络里到处是"矩阵乘向量、再过一个 逐元素的非线性函数"，逐个元素地推导会非常繁琐。矩阵微积分给出了一套"成块"推导的规则： 雅可比矩阵（Jacobian）：函数 f: R^n → R^m 的所有偏导数排成的 m × n 矩阵， 第 i 行第 j 列 = ∂f_i / ∂x_j。 链式法则 ：复合函数的雅可比矩阵 = 各层雅可比矩阵的乘积。 本文件用**数值方法**逐一验证讲义第 4 节的结论： 演示 1：几个基本函数的雅可比矩阵（线性层、逐元素函数、softmax） 演示 2：逐元素函数的雅可比矩阵是对角阵 —— 为什么实现时永远不要真的把它造出来 演示 3：softmax + 交叉熵的梯度 = ŷ - y（第 1 讲老朋友） 演示 4：一层神经网络 s = uᵀ tanh(Wx + b) 的全部梯度，以及"形状约定" 演示 5：窗口分类中，梯度怎样传回到每个词向量 ----------------------------------------------------------------------------- 【运行方式】 ----------------------------------------------------------------------------- python 01_matrix_calculus.py 需要 numpy。运行时间约 1 秒。 """ import time import numpy as np def print_header(title): """打印一个带分隔线的标题。""" print() print("=" * 70) print(f" {title}") print("=" * 70) # ============================================================================= # 第 1 部分：数值雅可比矩阵 # ============================================================================= def numerical_jacobian(f, x, eps=1e-6): """ 用中心差分数值地计算雅可比矩阵 J[i, j] ≈ ∂f_i / ∂x_j。 【思路】 每次只扰动 x 的一个分量 x_j（加减 eps），看输出向量的每个分量变了多少， 得到雅可比矩阵的第 j **列**。对所有 j 重复一遍。 【参数】 f : 函数，输入一维数组 (n,)，输出一维数组 (m,)（标量输出请包成长度 1 的数组） x : np.ndarray (n,) eps : 扰动大小 【返回】 np.ndarray (m, n) 【边界】 - 需要 2n 次调用 f，只适合小规模检查。 - x 会被复制，不会被修改。 """ x = np.array(x, dtype=float) m = len(np.atleast_1d(f(x))) J = np.zeros((m, len(x))) for j in range(len(x)): xp, xm = x.copy(), x.copy() xp[j] += eps xm[j] -= eps J[:, j] = (np.atleast_1d(f(xp)) - np.atleast_1d(f(xm))) / (2 * eps) return J def max_abs_diff(a, b): """两个数组对应元素之差的最大绝对值。""" return float(np.max(np.abs(np.asarray(a) - np.asarray(b)))) def softmax(z): """数值稳定的 softmax（见第 1 讲代码 04）。""" e = np.exp(z - z.max()) return e / e.sum() # ============================================================================= # 演示 # ============================================================================= def demo_basic_jacobians(): """演示 1：线性层、逐元素函数、softmax 的雅可比矩阵。""" print_header("演示 1：几个基本函数的雅可比矩阵（解析公式 vs 数值计算）") rng = np.random.default_rng(0) n, m = 4, 3 W, b, x = rng.normal(size=(m, n)), rng.normal(size=m), rng.normal(size=n) # (a) 线性层 z = Wx + b：∂z/∂x = W J_num = numerical_jacobian(lambda v: W @ v + b, x) print(f"\n (a) z = Wx + b（W 形状 {W.shape}）：∂z/∂x = W，数值结果与 W 的最大差 {max_abs_diff(J_num, W):.1e}") # (b) 逐元素 h = tanh(z)：∂h/∂z = diag(1 - tanh(z)²) z = rng.normal(size=m) J_num = numerical_jacobian(np.tanh, z) J_ana = np.diag(1 - np.tanh(z) ** 2) print(f" (b) h = tanh(z)（逐元素）：∂h/∂z = diag(1 - tanh²(z))，最大差 {max_abs_diff(J_num, J_ana):.1e}") print(" 数值算出的雅可比矩阵（非对角元素都是 0）：") for row in J_num: print(" " + " ".join(f"{v:+.3f}" for v in row)) # (c) softmax：∂p/∂z = diag(p) - p pᵀ p = softmax(z) J_num = numerical_jacobian(softmax, z) J_ana = np.diag(p) - np.outer(p, p) print(f" (c) p = softmax(z)：∂p/∂z = diag(p) - ppᵀ，最大差 {max_abs_diff(J_num, J_ana):.1e}") print(f" 每一列的和 = {np.round(J_num.sum(axis=0), 8)}（概率之和恒为 1，所以怎么改 z，变化量加起来都是 0）") # (d) 标量对向量：s = uᵀh，∂s/∂h = uᵀ（1 × m 的行向量） u = rng.normal(size=m) J_num = numerical_jacobian(lambda v: np.array([u @ v]), z) print(f" (d) s = uᵀh：∂s/∂h = uᵀ（形状 1 × {m}），最大差 {max_abs_diff(J_num, u[None, :]):.1e}") def demo_elementwise_cost(): """演示 2：逐元素函数的雅可比矩阵是对角阵，实现时用逐元素乘法代替矩阵乘法。""" print_header("演示 2：逐元素函数 → 对角雅可比矩阵 → 用逐元素乘法代替") rng = np.random.default_rng(1) for n in [100, 1000, 4000]: z, upstream = rng.normal(size=n), rng.normal(size=n) t0 = time.perf_counter() J = np.diag(1 - np.tanh(z) ** 2) # 真的造出 n × n 的对角阵 g1 = upstream @ J t_full = time.perf_counter() - t0 t0 = time.perf_counter() g2 = upstream * (1 - np.tanh(z) ** 2) # 逐元素乘法 t_elem = time.perf_counter() - t0 assert np.allclose(g1, g2) print(f" n = {n:>5d}：造对角阵再相乘 {t_full * 1e3:7.2f} 毫秒（{n * n * 8 / 1e6:6.1f} MB），" f"逐元素乘法 {t_elem * 1e3:6.3f} 毫秒，结果相同") print(""" 观察：数学上写 "上游梯度 × diag(f'(z))"，实现时永远写成 "上游梯度 ∘ f'(z)"（∘ 是逐元素乘法）。 对角阵的 n² 个元素里只有 n 个非零，造出来既浪费内存又浪费时间。 深度学习框架里所有逐元素运算（激活函数、加法、乘以常数……）的反向传播都是这样实现的。""") def demo_softmax_cross_entropy(): """演示 3：softmax + 交叉熵，对 logits 的梯度是 ŷ - y。""" print_header("演示 3：softmax + 交叉熵的梯度 = ŷ - y") rng = np.random.default_rng(2) z = rng.normal(size=5) y_idx = 3 loss = lambda v: np.array([-np.log(softmax(v)[y_idx])]) g_num = numerical_jacobian(loss, z)[0] g_ana = softmax(z) - np.eye(5)[y_idx] print(f"\n logits z = {np.round(z, 3)}，真实类别 = {y_idx}") print(f" 解析梯度 ŷ - y = {np.round(g_ana, 4)}") print(f" 数值梯度 = {np.round(g_num, 4)}，最大差 {max_abs_diff(g_num, g_ana):.1e}") print(""" 推导思路（讲义第 4.4 节）：用链式法则，∂J/∂z = ∂J/∂p · ∂p/∂z = (-1/p_y) e_yᵀ · (diag(p) - ppᵀ)， 化简后正好是 pᵀ - e_yᵀ。两个"复杂"的雅可比矩阵相乘，得到一个极其简单的结果 —— 这也是为什么所有框架都把 softmax 和交叉熵合并成一个函数（数值上也更稳定）。""") def demo_one_layer_network(): """演示 4：s = uᵀ tanh(Wx + b) 的全部梯度，以及形状约定。""" print_header("演示 4：一层神经网络 s = uᵀ tanh(Wx + b) 的梯度") rng = np.random.default_rng(3) n, m = 4, 3 W, b, x, u = rng.normal(size=(m, n)), rng.normal(size=m), rng.normal(size=n), rng.normal(size=m) def forward(W, b, x, u): return u @ np.tanh(W @ x + b) # 解析梯度（讲义第 4.5 节） z = W @ x + b h = np.tanh(z) delta = u * (1 - h ** 2) # δ = ∂s/∂z = uᵀ diag(f'(z))，写成逐元素乘法 grads = { "∂s/∂u（= h）": (h, lambda v: forward(W, b, x, v), u), "∂s/∂b（= δ）": (delta, lambda v: forward(W, v, x, u), b), "∂s/∂x（= Wᵀδ）": (W.T @ delta, lambda v: forward(W, b, v, u), x), } print(f"\n 形状：x ∈ R^{n}，W ∈ R^{m}×{n}，b、u ∈ R^{m}；δ = u ∘ tanh'(z) ∈ R^{m}\n") for name, (ana, f, at) in grads.items(): num = numerical_jacobian(lambda v: np.array([f(v)]), at)[0] print(f" {name:<16s}最大差 {max_abs_diff(ana, num):.1e}") # ∂s/∂W：形状约定 —— 梯度和参数同形状 gW_ana = np.outer(delta, x) # δ xᵀ，形状 (m, n)，和 W 一样 num_flat = numerical_jacobian(lambda v: np.array([forward(v.reshape(m, n), b, x, u)]), W.ravel()) print(f" {'∂s/∂W（= δxᵀ）':<16s}最大差 {max_abs_diff(gW_ana, num_flat.reshape(m, n)):.1e}") print(f""" 形状约定（shape convention）： 严格按雅可比矩阵的定义，标量 s 对 {m}×{n} 的矩阵 W 求导，得到的是 1 × {m * n} 的长行向量（数值方法算出的就是它）。 但为了做梯度下降 W ← W - α ∂s/∂W，我们希望梯度和 W **同形状**。 所以按惯例把它排成 {m}×{n} 的矩阵：第 i 行第 j 列 = ∂s/∂W_ij = δ_i · x_j，也就是外积 δ xᵀ。 直觉：W_ij 连接输入 x_j 和隐藏单元 i，它的梯度 = "上游传到单元 i 的误差 δ_i" × "这条边的输入 x_j"。""") # 复用：δ 被 ∂s/∂b、∂s/∂W、∂s/∂x 共同使用 print(" 注意 δ 在 ∂s/∂b、∂s/∂W、∂s/∂x 中被重复使用 —— 算一次、到处用，这就是反向传播高效的原因（代码 02）。") def demo_window_gradient(): """演示 5：窗口分类中，梯度怎样传回每个词向量。""" print_header("演示 5：窗口分类 —— 梯度传回到每个词向量") rng = np.random.default_rng(4) d, m, window = 3, 4, ["museums", "in", "paris", "are", "amazing"] E = {w: rng.normal(size=d) for w in window} # 每个词一个 d 维向量 x = np.concatenate([E[w] for w in window]) # 拼接成 5d 维的窗口向量 W, b, u = rng.normal(size=(m, len(x))), rng.normal(size=m), rng.normal(size=m) h = np.tanh(W @ x + b) delta = u * (1 - h ** 2) grad_x = W.T @ delta # (5d,) print(f"\n 窗口：{window}，每个词向量 {d} 维，拼接后 x 是 {len(x)} 维") print(f" ∂s/∂x 是 {len(grad_x)} 维，把它按 {d} 个一组切开，就是每个词向量的梯度：") for k, w in enumerate(window): print(f" ∂s/∂v_{w:<8s}= {np.round(grad_x[k * d:(k + 1) * d], 3)}") num = numerical_jacobian(lambda v: np.array([u @ np.tanh(W @ v + b)]), x)[0] print(f" 数值检查：最大差 {max_abs_diff(grad_x, num):.1e}") print(""" 观察：拼接（concatenate）在反向传播时就是"切开"（split）。 如果一个词在窗口里出现了两次，它的两段梯度要**加起来**（同一个参数被用了两次，梯度累加）。 是否要用这些梯度去更新词向量本身，是一个需要权衡的问题（讲义第 4.6 节、代码 04 演示 4）。""") if __name__ == "__main__": demo_basic_jacobians() demo_elementwise_cost() demo_softmax_cross_entropy() demo_one_layer_network() demo_window_gradient()