雅可比与海森

上一课我们为标量函数 \(f:\mathbb{R}^n\to\mathbb{R}\) 装上了"梯度 \(\nabla f\)"这台指南针，它在每个点指向上升最快的方向。但神经网络里到处是向量进、向量出的映射：一层把激活向量变成另一个激活向量。这一课我们把导数推广到向量值函数，得到雅可比矩阵 Jacobian；再回头给标量函数装上"二阶导数"——海森矩阵 Hessian，它描述函数的弯曲程度（曲率 curvature），正是优化难易、鞍点、二阶方法的几何根源。

一、雅可比：向量值函数的导数

考虑一个把 \(n\) 维向量映成 \(m\) 维向量的函数

它有 \(m\) 个输出分量，每个分量又依赖 \(n\) 个输入变量。于是"导数"不再是一个数，也不再是一个向量，而是把所有偏导数 \(\partial f_i/\partial x_j\) 排成的一张表——这就是雅可比矩阵。

两个熟悉的特例帮你定位它在知识地图上的位置：

• 当 \(m=1\)（标量输出）时，\(J\) 只有一行，就是梯度的转置 \(\nabla f^\top\)（一个行向量）。所以梯度是雅可比的特例。（沿用上一课的约定：梯度 \(\nabla f\) 默认是列向量，它转置后才是雅可比那一行。）
• 当 \(n=1,m=1\) 时，\(J\) 退化成一个数——就是高中的导数 \(f'(x)\)。

雅可比 = 局部线性映射

把这句话写成公式，就是一阶线性化（first-order linearization）。设在点 \(x\) 处给一个微小扰动 \(\Delta\in\mathbb{R}^n\)，则

右边 \(J\Delta\) 是矩阵乘向量，结果是 \(m\) 维向量，与 \(f(x)\) 同形——量纲对上了。这就是"用一块切平面去贴住弯曲曲面"的多维推广。误差是 \(O(\|\Delta\|^2)\)，扰动越小越准。

方阵时：\(\det J\) 是局部体积缩放因子

当 \(m=n\) 时 \(J\) 是方阵，它的行列式有一个漂亮的几何意义。回忆线性代数：一个矩阵 \(A\) 把单位立方体映成平行多面体，体积变成 \(|\det A|\) 倍。既然 \(f\) 在局部就是线性映射 \(J\)，于是：

"翻面"是什么直觉？ 就像照镜子——左手照出来变成右手；在二维里，它相当于把图形沿某条线做一次镜像翻转，原本逆时针绕的小圈会被映成顺时针。所以 \(\det J\) 同时编码了两件事：绝对值 \(|\det J|\) 是面积/体积缩放倍数，符号则告诉你有没有镜像翻面。

这正是换元积分里那个雅可比行列式因子的来历（你在概率论里会反复见到）。设 \(y=f(x)\) 为正向变换，\(J=J_f(x)\) 是 \(f\) 在 \(x\) 处的雅可比。若随机向量 \(y=f(x)\) 是 \(x\) 的可逆变换，则概率密度满足变量替换公式：

体积被 \(f\) 放大 \(|\det J_f|\) 倍，概率要守恒（总和为 1），所以密度就得除以这个因子。注意指数 \(-1\) 是因为这里 \(J_f\) 是正向映射 \(x\to y\) 的雅可比；若改用逆向雅可比 \(J_{f^{-1}}(y)\)，由于 \(\det J_{f^{-1}}=1/\det J_f\)，指数就变成 \(+1\)。两种写法等价，关键是认清你算的是哪个方向的雅可比。这条公式是生成模型里"标准化流（normalizing flows）"的数学心脏。

二、海森：标量函数的二阶导数（曲率）

梯度告诉我们函数往哪儿走、走多陡；但它说不清函数弯得多厉害——是平缓的碗，还是狭长的山谷，还是马鞍。要回答这个，得求二阶导数。对标量函数 \(f:\mathbb{R}^n\to\mathbb{R}\)，把所有二阶偏导数排成方阵，就是海森矩阵。

海森是对称的：只要二阶偏导连续，混合偏导可交换次序（克莱罗定理 Clairaut's theorem），即 \(\dfrac{\partial^2 f}{\partial x_i\partial x_j}=\dfrac{\partial^2 f}{\partial x_j\partial x_i}\)，于是 \(H=H^\top\)。对称矩阵有实特征值、可正交对角化（这是对称矩阵谱定理的结论，可回顾线性代数那一课）——这一点后面判断曲率时至关重要。

二阶泰勒展开

一维泰勒展开 \(f(x+\Delta)\approx f+f'\Delta+\tfrac12 f''\Delta^2\) 的多维版本，把 \(f'\) 换成梯度、\(f''\) 换成海森：

曲率与海森的特征值

二次型 \(\Delta^\top H\Delta\) 的"脾气"由 \(H\) 的特征值（eigenvalue）决定。把 \(H\) 正交对角化 \(H=Q\Lambda Q^\top\)（即上面提到的对称矩阵谱定理，此处即引即用），沿特征向量（eigenvector）方向看，函数就是一族独立的抛物线，第 \(k\) 条的"曲率"正比于特征值 \(\lambda_k\)：

• 沿 \(\lambda_k>0\) 的方向，函数向上弯（像碗底，往两边都升高）；
• 沿 \(\lambda_k<0\) 的方向，函数向下弯（像山顶，往两边都降低）；
• \(\lambda_k=0\) 表示该方向是平的（局部退化）。

由此得到临界点（\(\nabla f=0\) 的点）分类，只看 \(H\) 的特征值符号：

凸性 ⇔ 海森半正定

这就是为什么线性回归、逻辑回归这类凸问题"好优化"，而深层神经网络的损失非凸、海森到处变号、布满鞍点，优化才那么棘手。

三、牛顿法：用海森修正梯度

梯度下降 \(x\leftarrow x-\eta\nabla f\) 只用了一阶信息，它不知道每个方向弯得多厉害，所以只能用一个全局学习率 \(\eta\) 凑合。牛顿法（Newton's method）更聪明：它在当前点用二阶泰勒展开搭一个抛物面，然后一步跳到这个抛物面的最低点。

对二阶近似 \(q(\Delta)=f+\nabla f^\top\Delta+\tfrac12\Delta^\top H\Delta\) 求最小，令其梯度为零：\(\nabla f+H\Delta=0\)，解出 \(\Delta=-H^{-1}\nabla f\)。于是牛顿更新为

四、向量-雅可比积（VJP）——下一课的钥匙

牛顿法展示了二阶信息多有用，但代价是要算 \(H^{-1}\)；回到现实，深度学习几乎全靠一阶方法——而支撑一阶反向传播的最小积木，就是下面要讲的 VJP。它把我们从"二阶（理想但昂贵）"自然带回"一阶（便宜且实用，下一课的主角）"。

反向传播的全部秘密，可以浓缩成一句话：它从不显式构造庞大的雅可比矩阵，只反复计算"向量乘雅可比"这个组合。给定某一层的雅可比 \(J\in\mathbb{R}^{m\times n}\)，以及一个暂时看作任意给定的 \(m\) 维向量 \(v\in\mathbb{R}^m\)（沿用全篇约定：向量默认是列向量），所谓向量-雅可比积（vector–Jacobian product, VJP）就是把 \(v\) 转置后左乘 \(J\)：

这两种写法只差一个转置：\(v^\top J\) 是行向量，\(J^\top v=(v^\top J)^\top\) 是列向量。机器学习库里习惯让"梯度与参数同形（列向量）"，所以最终拿到手的是 \(J^\top v\in\mathbb{R}^n\)。\(v\) 是哪来的？ 这里你只需把它看成一个给定的 \(m\) 维向量；在下一课它会是损失对本层输出的梯度，这里无需剧透，只要知道 VJP 是个比存下整个 \(J\) 便宜得多的运算即可。

为什么不直接存 \(J\)？因为在神经网络里 \(m,n\) 可能是几千几万，\(J\) 大到存不下；但 \(J^\top v\) 只是一个 \(n\) 维向量，而且对许多层（如逐元素激活、矩阵乘）有不依赖于显式 \(J\) 的快速算法。把每一层的 VJP 从输出端一路串到输入端，就是反向传播。

五、用 numpy 数值验证手算

有限差分（finite difference）能在不会求导的情况下逼近雅可比与海森，是检查手推公式的利器。下面这段只用 numpy，几秒跑完，把上面两道例题都验一遍。

import numpy as np

# ---------- 1) 数值雅可比：中心差分 ----------
def num_jacobian(f, x, h=1e-5):
    x = np.asarray(x, float)
    fx = np.asarray(f(x), float)
    m, n = fx.size, x.size
    J = np.zeros((m, n))
    for j in range(n):                 # 逐列：扰动第 j 个输入
        e = np.zeros(n); e[j] = h
        J[:, j] = (np.asarray(f(x + e)) - np.asarray(f(x - e))) / (2 * h)
    return J

# f: R^2 -> R^2 , f(x,y) = [x^2 y , x + sin y]
def f(v):
    x, y = v
    return np.array([x**2 * y, x + np.sin(y)])

p = np.array([1.0, np.pi / 2])
J_num = num_jacobian(f, p)
J_hand = np.array([[2*p[0]*p[1], p[0]**2],
                   [1.0,          np.cos(p[1])]])   # = [[pi,1],[1,0]]
print("数值雅可比 J =\n", np.round(J_num, 6))
print("手算雅可比 J =\n", np.round(J_hand, 6))
print("det(J) 数值 =", round(np.linalg.det(J_num), 6), " (手算 = -1)")

# ---------- 2) 数值海森：二阶中心差分 ----------
def num_hessian(g, x, h=1e-4):
    x = np.asarray(x, float); n = x.size
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            ei = np.zeros(n); ei[i] = h
            ej = np.zeros(n); ej[j] = h
            H[i, j] = (g(x+ei+ej) - g(x+ei-ej)
                       - g(x-ei+ej) + g(x-ei-ej)) / (4*h*h)
    return 0.5 * (H + H.T)             # 对称化，抹掉差分噪声

# g: R^2 -> R , g(x,y) = x^2 - y^2 ，临界点 (0,0)
def g(v):
    x, y = v
    return x**2 - y**2

H = num_hessian(g, np.array([0.3, -0.4]))   # 海森处处为常数阵
eig = np.linalg.eigvalsh(H)                 # 对称阵用 eigvalsh
print("\n数值海森 H =\n", np.round(H, 4), " (手算 = [[2,0],[0,-2]])")
print("特征值 =", np.round(eig, 4))
kind = ("极小" if np.all(eig > 0) else
        "极大" if np.all(eig < 0) else "鞍点")
print("临界点 (0,0) 类型 ->", kind)        # -> 鞍点

运行后你会看到数值结果与手算几乎逐位吻合：雅可比为 \([[\pi,1],[1,0]]\)、行列式 \(-1\)，海森为 \(\mathrm{diag}(2,-2)\)、特征值 \(\{2,-2\}\)，判定为鞍点。

调一调，观察现象

下面三个微任务都基于本课已验证的公式，改一个数、看一个数怎么变，把"曲率、行列式、条件数"从符号变成你亲眼跑出来的现象。每段都只用 numpy 和 print，几秒跑完。

import numpy as np
def f(v):
    x, y = v
    return np.array([x**2 * y, x + np.sin(y)])
p = np.array([1.0, np.pi/2])
J = np.array([[2*p[0]*p[1], p[0]**2],
              [1.0,          np.cos(p[1])]])   # = [[pi,1],[1,0]]
for d in [np.array([0.01, -0.02]), np.array([0.001, -0.002])]:
    err = np.linalg.norm(f(p+d) - (f(p) + J @ d))
    print("delta =", d, " 线性化误差 =", err)
# 第二行误差约为第一行的 1/100

import numpy as np
def classify(H):
    eig = np.linalg.eigvalsh(H)        # 对称阵用 eigvalsh
    if np.all(eig > 0): k = "极小"
    elif np.all(eig < 0): k = "极大"
    else: k = "鞍点"
    return np.round(eig, 2), k
for a, b in [(2, 2), (2, -2), (-2, -2)]:
    H = np.diag([a, b]).astype(float)
    print("对角元", (a, b), " 特征值/类型 ->", classify(H))

import numpy as np
def steps_to_converge(H, tol=1e-3, max_iter=100000):
    lam = np.diag(H)
    eta = 2.0 / (lam.max() + lam.min())   # 各自的最优学习率
    x = np.ones(H.shape[0])
    for k in range(max_iter):
        if np.linalg.norm(x) < tol:
            return k
        x = x - eta * (H @ x)             # 梯度 = H x
    return max_iter
for H in [np.diag([1., 1.]), np.diag([1., 10.]), np.diag([1., 100.])]:
    print(f"kappa={np.linalg.cond(H):6.1f}  收敛步数={steps_to_converge(H)}")
# 1 -> 37 -> 363：条件数越大，一阶法越慢

动手练习

1. 写雅可比。 对极坐标映射 \(f(r,\theta)=[\,r\cos\theta,\ r\sin\theta\,]\) 手推 \(J\)，并证明 \(\det J=r\)（这正是二重积分换元里的 \(r\,dr\,d\theta\) 因子）。用下面骨架在 \((r,\theta)=(2,\pi/3)\) 处验证：

   import numpy as np
   def f(v):
       r, t = v
       return np.array([r*np.cos(t), r*np.sin(t)])
   # TODO: 复用上文 num_jacobian，打印数值 J 与 det，并与手算 det=r=2 对比
   

2. 判类型。 对 \(f(x,y)=x^2+xy+y^2\)，手算梯度找临界点、写出海森、求特征值，判断类型；再用 num_hessian 验证特征值是否为 \(\{1,3\}\)。
3. 线性化的精度。 取第一道例题的 \(f\) 与点 \(p=(1,\pi/2)\)，令 \(\Delta=[0.01,-0.02]^\top\)，分别算真值 \(f(p+\Delta)\) 与线性近似 \(f(p)+J\Delta\)，打印两者之差的范数 np.linalg.norm(...)。再把 \(\Delta\) 缩小 10 倍，观察误差大约缩小多少倍（提示：应约为 \(100\) 倍，因为误差是 \(O(\|\Delta\|^2)\)）。
4. 条件数与优化。 我们要让"病态曲率拖累一阶法"这个现象真正显现，关键是给每个 \(H\) 都配上它各自能用的最优学习率，再比较收敛到 \(\|x\|<10^{-3}\) 各需多少步。对二次函数 \(f(x)=\tfrac12 x^\top H x\)（梯度 \(\nabla f=Hx\)），对角海森的最优学习率约为 \(\eta^\star=\dfrac{2}{\lambda_{\max}+\lambda_{\min}}\)。比较 \(H=\mathrm{diag}(1,100)\)（病态，\(\kappa=100\)）与 \(H=\mathrm{diag}(1,1)\)（理想，\(\kappa=1\)）：前者约需数百步、后者一步即到，对比鲜明。
   原理点睛：一阶法的步长被 \(\lambda_{\max}\) 卡死（步子再大就发散），收敛却被 \(\lambda_{\min}\) 拖死（平缓方向爬得慢），两者之比 \(\kappa=\lambda_{\max}/\lambda_{\min}\) 才是元凶——这与正文 ML 小节的条件数论述闭环。

   import numpy as np
   
   def steps_to_converge(H, tol=1e-3, max_iter=100000):
       lam = np.diag(H)
       eta = 2.0 / (lam.max() + lam.min())   # 每个 H 各自的最优学习率
       x = np.ones(H.shape[0])
       for k in range(max_iter):
           if np.linalg.norm(x) < tol:
               return k, eta
           x = x - eta * (H @ x)             # f = 0.5 x^T H x 的梯度 = H x
       return max_iter, eta
   
   for H in [np.diag([1.0, 1.0]), np.diag([1.0, 100.0])]:
       k, eta = steps_to_converge(H)
       print(f"kappa={np.linalg.cond(H):6.1f}  eta*={eta:.4f}  收敛步数={k}")
   # 你会看到 kappa=1 只要 1 步，kappa=100 要数百步：条件数越大，一阶法越慢
   

5. 牛顿一步到位。 同上二次函数（取正定的 \(H=\mathrm{diag}(1,100)\)），用牛顿更新 \(x\leftarrow x-H^{-1}\nabla f\) 走一步，验证无论 \(H\) 多病态都一步到达原点（打印更新后的 \(\|x\|\)，应≈0）。再思考：若把 \(H\) 换成不定的 \(\mathrm{diag}(1,-100)\)，牛顿一步会落到哪里？（提示：仍落到 \(\nabla f=0\) 的临界点，但那是鞍点，不是极小。）

掌握自检

• 给定 \(f:\mathbb{R}^n\to\mathbb{R}^m\)，我能立刻说出雅可比的形状（\(m\times n\)）、哪个下标是行哪个是列，并解释 \(f(x+\Delta)\approx f(x)+J\Delta\) 为何维度自洽。
• 我能解释 \(m=n\) 时 \(|\det J|\) 是体积缩放因子、\(\det J<0\) 代表镜像翻面，并说出它在概率密度变量替换中的作用（以及指数 \(\pm1\) 取决于用正向还是逆向雅可比）。
• 我能说清"海森 = 梯度的雅可比"，会写二阶泰勒展开，并知道它为何对称（克莱罗定理 + 谱定理）。
• 给一个 \(2\times2\) 海森，我能算特征值并判定极小/极大/鞍点；我知道"凸 ⇔ 海森半正定"。
• 我能用一句话说出牛顿法 \(x\leftarrow x-H^{-1}\nabla f\) 的直觉（仅在海森正定时一步到极小），以及条件数 \(\kappa=\lambda_{\max}/\lambda_{\min}\) 为何决定一阶方法的快慢。
• 我能复述"反向传播 = 一串 VJP（\(v^\top J\) 或等价的 \(J^\top v\)）"，说明为何不显式存雅可比，并知道梯度默认按列向量约定。

下一课，我们把这些一阶导数（每层的雅可比/VJP）用链式法则沿计算图从输出端一路乘回输入端——这串 VJP 的连乘，就是驱动整个深度学习的反向传播。