看到标题里的 Jacobin我第一反应是拼错了——大概率是 Jacobian。这个笔误几乎每周都能在技术群里见到但问题不大真正要搞懂的是 Jacobian 矩阵和 Hessian 矩阵这一对组合到底在干什么。搞数值优化、机器人运动学、参数辨识、三维重建和机器学习的人基本天天和它们打交道但很多新手要么停留在“雅可比是个矩阵海森是二阶导”的背诵层面要么写代码时把维度搞反让优化器原地爆炸。这篇文章我想换个讲法不按教科书顺序罗列定义而是从工程视角把两个矩阵讲透。1. 雅可比矩阵的核心它把局部变化翻译成了矩阵1.1 先来一个可以直接算的例子很多人一看到 Jacobian 的定义就头疼因为它总是用一大串符号J[i][j] ∂f_i / ∂x_j翻译成人话就是我的第i个输出对第j个输入求偏导然后把所有偏导排成一张表。举个最简单的二维映射f1(x, y) x^2 sin(y) f2(x, y) x*y e^y这个f有两个输入x, y也有两个输出f1, f2所以 Jacobian 是 2×2 矩阵J [[2x, cos(y)], [y, x e^y]]拿点(x1, y0)代入具体得到J [[2, 1], [0, 2]]现在做一个小实验让x增加 0.01y不变。真实输出的变化约等于J第一列乘以 0.01也就是[0.02, 0.00]。让y增加 0.01真实输出的变化约等于J第二列乘以 0.01也就是[0.01, 0.02]。这个理解非常重要Jacobian 的每一列就是“某个输入分量变化时输出向量会朝哪个方向移动”。输入变一点输出就沿着对应的列方向滑出去输入几个分量同时变输出变化就是这些列的线性组合。哪怕函数再非线性只要扰动足够小这种线性翻译都是成立的。1.2 工程里真正关心的往往是 Jacobian 的列空间在参数辨识、机器人运动学、优化这些场景里我们很少只关心 Jacobian 的某一项而是关心整张矩阵里“列和列之间的关系”。简单说如果 Jacobian 的两列方向差不多意味着两个不同参数对输出的影响方式非常相似这时候数据很难把这两个参数区分开。比如某个模型里两个系数都以几乎相同的方式出现在输出里那么从观测数据反推参数时你会得到一个非常病态的问题参数组合有无穷多个拟合曲线却几乎一样。我在实际项目里碰到过类似情况。一开始用梯度下降怎么调都不收敛后来把 Jacobian 的列拿出来一看两列相关性高得离谱本质上是模型本身的可辨识性不足不是优化器不行。这时候要做的不是硬调迭代步长而是改模型结构、加约束或者换一组信息量更大的输入数据。如果处理的是方阵det(J) 0的点叫奇异位形。机器人的例子最直观末端速度v和关节速度q_dot的关系是v J(q) * q_dot。当机器人某个关节角度组合让 Jacobian 行列式变成零时某些方向的末端速度根本无法用关节运动产生这就是机械臂“卡死”或失去可控性的位置。很多控制器在这些位形附近会突然输出特别大的关节力矩原因就是 Jacobian 的逆接近无穷大微小误差被放大。1.3 从雅可比到海森一阶到二阶只差一次求导Jacobian 和 Hessian 不是两个孤立的概念。如果函数输出只有一个分量也就是说f是一个标量函数S(x)那它的 Jacobian 退化成一行偏导这就是我们常说的梯度。现在把梯度看作一个新的向量函数再对它求一次 Jacobian得到的就是 Hessian 矩阵H[i][j] ∂²S / ∂x_i ∂x_j我习惯这样记忆Hessian 是梯度的 Jacobian。例如标量函数S(x, y) 3x²y e^y先算梯度∇S [6xy, 3x² e^y]再对梯度求偏导H [[6y, 6x], [6x, e^y]]这个矩阵一定对称。原因是混合偏导和求导顺序无关∂²S/∂x∂y ∂²S/∂y∂x。工程上看到不对称的 Hessian第一反应应该是公式或代码写错了而不是“数学本来就如此”。2. Hessian 矩阵曲率不是玄学它直接决定优化方向2.1 光看梯度为零根本判断不了极值很多刚入门的朋友以为“梯度为零就是极小值”这在多元函数里完全不够。梯度为零只说明这是一个驻点它可能是极小值、极大值也可能是鞍点。最经典的例子是S(x, y) x² - y²在(0,0)处梯度是[0,0]。但沿x方向函数是下凸的像是谷底沿y方向函数是上凸的像是山顶。这个点既不是极小也不是极大它是一个鞍点。Hessian 能帮你区分在驻点处如果 Hessian 正定也就是所有特征值都大于零那这是局部极小值点如果 Hessian 负定所有特征值都小于零就是局部极大值点如果特征值有正有负那就是鞍点。对应上面例子Hessian 是diag(2, -2)特征值一正一负所以一眼就看出是鞍点。这也是为什么只靠梯度下降容易卡在鞍点附近尤其在高维神经网络损失函数里鞍点远比极小值点常见。2.2 牛顿步为什么比梯度下降更“聪明”把标量函数S在某个点x附近做二阶泰勒展开S(x δ) ≈ S(x) ∇S^T δ 0.5 * δ^T H δ如果 Hessian 正定极小点对应的δ满足H δ -∇S所以δ -H^{-1} ∇S这就是牛顿法。梯度下降的更新是δ -η ∇S区别在哪梯度下降用一个固定的标量学习率η给所有方向统一缩放。牛顿法则用H^{-1}对梯度做了一次线性变换相当于给每个方向配了不同的步长曲率大的方向走小步曲率小的地方走大步。我常用一个生活类比梯度告诉你在房间里往哪边走能下坡Hessian 告诉你这个坡是越来越陡还是越来越缓甚至前面是不是有个拐弯。只看梯度你可能会在一个长而平的“峡谷”里来回震荡带上曲率信息就能沿着谷底更快地滑下去。2.3 为什么现实里几乎没人显式存一个完整 Hessian理论很漂亮但工程里 Hessian 常常存不下来。原因很简单参数数量n一旦上来n × n的双精度浮点矩阵就吃不消了。参数数量 nHessian 元素个数约占用内存1001e480 KB1 万1e8800 MB100 万1e128 TB一百万个参数在机器学习里很常见8 TB 的矩阵基本告别常规服务器。所以现实工程中很少直接构造和存储完整 Hessian而会改用以下三类替代思路拟牛顿法比如 L-BFGS用最近几步的梯度差分信息隐式逼近H^{-1}只需要存一小段历史向量。Gauss-Newton 近似在最小二乘问题里用J^T J近似 Hessian后面会详细展开。Hessian-vector product利用自动微分算H v不需要把H形成矩阵。这些手段本质上都是在和“Hessian 太贵”这件事博弈。你能理解 Hessian 是什么再看这些方法就会知道它们到底在近似什么、牺牲了什么。3. 一个完整的拟合案例Jacobian 和 Hessian 如何联手干活3.1 建模先写出残差和解析 Jacobian理论知识说再多不如把一个具体问题跑通。假设我有若干组观测数据(x_i, y_i)希望用衰减指数模型去拟合y a * e^{-b x} c要估计的参数是p [a, b, c]。定义每个数据点的残差r_i(p) y_i - (a * e^{-b x_i} c)目标函数是S(p) 0.5 * Σ r_i^2这里0.5纯粹是为了让后面的梯度表达式少个系数完全不影响最优解。如果直接对r_i求偏导可以得到每一行残差的解析 Jacobian∂r_i / ∂a -e^{-b x_i} ∂r_i / ∂b a * x_i * e^{-b x_i} ∂r_i / ∂c -1如果有一千个观测点J就是一千行乘三列的矩阵。每一行对应一个残差每一列对应一个参数。很多初学者会习惯用有限差分去算 Jacobian省事但更好的是写出解析形式。解析求导速度快精度也高尤其是接近最优值时数值差分的误差会明显干扰收敛。3.2 高斯牛顿和 LMJ 怎么变出近似 Hessian最小二乘目标函数的海森矩阵其实可以写成两部分H Σ (J_i^T J_i r_i * ∇²r_i)第一项是J^T J第二项是残差对参数的二阶导乘以残差值本身。高斯牛顿法的核心近似就是把第二项丢掉直接令H ≈ J^T J为什么可以丢接近最优解时残差r_i本身已经很小乘以二阶导后就更小所以第二项相对于第一项不关键。但在迭代刚开始、残差还很大时这个近似就不太靠谱步子容易迈得太大导致发散。于是工程上几乎都用列文伯格-马夸尔特方法也就是加一个阻尼项A J^T J λ Iλ很大时A ≈ λ I更新步长约等于梯度下降的一小步稳定但慢。λ很小时A ≈ J^T J退化为高斯牛顿接近最优解时收敛很快。你只需要提供一次 Jacobian就能在迭代中通过调节λ自动在“稳妥的梯度下降”和“快速的高斯牛顿”之间切换。这就是为什么很多非线性最小二乘求解器让你返回 Jacobian而不要求你返回 Hessian。3.3 一个 10 行左右的 LM 迭代骨架下面是一个极简实现矩阵运算用常规线性代数库就能完成重点看迭代结构# Python 风格伪代码矩阵库负责 transpose、solve、norm 等操作 def residual(p, x, y): a, b, c p return y - (a * exp(-b * x) c) def jac_res(p, x): a, b, c p J empty((len(x), 3)) J[:, 0] -exp(-b * x) J[:, 1] a * x * exp(-b * x) J[:, 2] -1.0 return J def lm_fit(x, y, p0, lam1e-3, max_iter100): p p0[:] for _ in range(max_iter): r residual(p, x, y) J jac_res(p, x) A transpose(J) J lam * identity(3) g transpose(J) r step solve_linear(A, -g) new_p p step new_r residual(new_p, x, y) if norm(new_r) norm(r): p new_p lam * 0.7 else: lam * 2.0 if norm(step) 1e-12: break return p这里A就是近似的 Hessiang是目标函数的梯度。每次迭代如果残差下降就接受步长并减小λ否则增大λ让更新更保守。实际经验里有几个点值得留意初始λ可以从1e-3或1e-2开始太小容易一开始就发散太大则前几步很像梯度下降收敛慢。如果数据量级很大建议先做归一化否则J^T J的条件数会很难看。步长范数的停机阈值不要拍脑袋定最好看你参数的物理量级参数本身在1e-6量级时阈值设在1e-12会永远不触发。4. 工程里的四道坎从拼写到稀疏结构都有人栽过4.1 行和列方向搞反是最常见的“公式致命伤”把 Jacobian 拼写成 Jacobin 只是小尴尬真正致命的是维度方向写反。绝大多数优化器期望的 Jacobian 是“观测数 × 参数数”也就是每一行对应一个残差每一列对应一个参数。我见过很多人的代码里把J返回到参数数乘观测数之后J^T J的尺寸完全不对。最诡异的是一部分矩阵库并不会立刻报维度错误而是算出一个“莫名其妙”的结果让你绕很久才发现问题。这里给一个快速检查办法看J^T J的维度。如果参数数量是m那么J^T J必须是m × m。算出来的矩阵不是这个尺寸就说明 Jacobian 的维度反了。另一个更直接的检查是只让某一个参数变化一点点看残差向量的变化方向是不是和 Jacobian 里对应的那一列一致。这个检查不需要任何高级工具几行代码就能做。4.2 用有限差分验证解析导数别盲信手推公式手推 Jacobian 和 Hessian 很容易出错尤其是链式法则套来套去的时候。我的习惯是凡是写进代码里的导数公式先拿有限差分做一次交叉验证再放进优化器。中心差分的公式很简单df/dx ≈ (f(x h) - f(x - h)) / (2h)更粗糙的单侧差分(f(xh)-f(x))/h精度差一截能不用就不用。h并不是越小越好。受浮点舍入误差限制取h为sqrt(eps) * max(1, |x|)通常量级在1e-8左右太大会有截断误差太小会全是舍入噪声。验证 Hessian 的时候我一般先验证梯度再对梯度做中心差分得到近似 Hessian 和解析值对比。误差在1e-5这个量级很正常毕竟是差分近似但如果误差大到1e-2基本可以确定公式写错了。还有一个隐藏问题有些优化算法内部会用分解算法要求对称矩阵。哪怕你解析 Hessian 在数学上对称浮点舍入也可能让它不对称。拿到 Hessian 后顺手做一次(H H^T) / 2强制对称能省掉许多莫名其妙的线性代数报错。4.3 大规模稀疏结构J 千万别硬生生转成稠密矩阵大规模优化里最典型的例子是三维重建中的集束调整问题。每个重投影残差只与“当前相机的位姿”和“被观测的那个三维点”有关和场景里其他相机、其他点没有关系。于是整个 Jacobian 非常稀疏每一行只有少数几个非零块。但如果用最无脑的方式把完整J存成稠密矩阵再算J^T J内存和计算量都会爆炸。正确做法是保留稀疏结构用稀疏线性代数求解集束调整还会用分块消元利用J^T J的箭头形稀疏结构高效求解。很多新人第一次接触这种感觉是在跑大型非线性优化时节点一多机器突然内存爆掉。卡了很久才发现是稠密矩阵的问题。遇到这类问题第一反应不要是“优化器不够好”而是检查有没有过度稠密化。4.4 参数多到 Hessian 放不进内存时怎么办当模型参数到了百万级完整 Hessian 根本没法显式构造。这时候常见的选择有三类L-BFGS只保存最近几步的梯度差和参数差隐式构造逼近H^{-1}的矩阵内存开销是O(n)不是O(n²)。Hessian-vector product借助自动微分直接算H v不需要形成完整H。很多二阶优化算法真正需要的往往只是H v不是每个元素。分块近似把 Hessian 按网络结构分块每块用低秩或 Kronecker 近似这类方法是很多神经网络二阶优化器的理论基础。在这些方法里你依然需要通过 Jacobian 或梯度去描述一阶信息再通过一些操作去逼近二阶信息。所以把 Jacobian 和 Hessian 的基本概念吃透再去看这些高阶工具思路会清晰很多。5. 最后聊一点学习路径和实操心态每次讲完矩阵总有人问我“该从哪里开始刷题”。我的建议一直很直接不要先去背性质先找一个最小的二维例子比如本文里那个f1x²sin(y), f2xye^y自己动手写出 Jacobian再挑一个标量函数写出 Hessian。能用手算清楚两个小例子后面再大的矩阵也只是维度变大而已。等你能手算了就去做一个最简单的指数拟合。这个项目几乎是最合适的练功房模型不复杂解析 Jacobian 好写迭代过程又完整覆盖了“一阶信息→近似二阶信息→步长控制”这条主线。跑通之后再去挑战机器人运动学或三维重建里的稀疏 Jacobian会顺畅许多。最后补一个我自己的经验写完任何导数代码先跑五分钟有限差分验证再丢进优化器。这个习惯帮我省下的调试时间远超验证本身花掉的时间。遇到优化不收敛也不要急着怀疑求解器先看J^T J的条件数、再看 Hessian 的特征值符号多数问题在这一步就暴露了。想清楚每一步在干什么比记住一堆矩阵公式重要得多。