数学建模的题目做到后期十有八九会撞到优化问题。曲线拟合要调参数资源分配要定方案运筹调度要找最优这些本质都是“在某个目标函数上求极值”。真到了多变量场景很多教材翻来覆去只讲梯度下降讲到牛顿法也就给个一维公式等到Hessian矩阵一出现大家基本就绕道走了。其实多变量最优化计算里的牛顿法是理解数值优化很关键的一环理解了它再看拟牛顿法、信赖域法都会顺畅很多。这篇数模笔记主要整理牛顿法在多变量优化里的完整推导、手算过程、Python实现以及实操中Hessian矩阵“翻车”时的补救方案。建议正在准备数学建模竞赛、算法竞赛或者刚接触数值优化计算的同学重点收藏这套笔记是按“公式推导—手算验证—代码实现—避坑总结”的流程走的可以直接照着用。1. 项目概述多变量牛顿法是什么能解决哪些问题1.1 从“求根”到“求极值”牛顿法的原始身份很多人第一眼看到牛顿法是在数值分析里解非线性方程的部分所谓Newton-Raphson迭代。当时的场景很简单给定一个一元函数f(x)我们想找到某个x*使f(x*) 0。做法是在当前点xk处用切线近似f(x)然后求切线与x轴的交点x_{k1} x_k - f(x_k) / f(x_k)这个迭代式看起来平平无奇但它的思想极其朴素任何光滑函数都可以在某一点的邻域内用低阶多项式逼近而在局部上用一次函数切线去逼近就已经能给出不错的下一步位置了。到了最优化计算里我们的目标变了不是找f(x) 0而是找f(x)取得极小值的点x*。但两者在数学上可以打一个“等号”如果f(x)可微那么在极值点处梯度g(x) ∇f(x) 0。这等于说优化问题可以转化为一个“梯度方程组求根”的问题。把牛顿法从“求f(x)0”迁移到“求∇f(x)0”就引出了多变量牛顿法。这也解释了为什么多变量牛顿法在收敛速度上有先天优势。它天生是冲着方程的根去的而极值点恰好是梯度向量的零点所以它对极值点的逼近是“二阶”级别的。梯度下降法只用了目标函数的一阶信息牛顿法把二阶信息Hessian矩阵也用了收敛快是自然的事。1.2 数学建模里哪些场景会用到它结合我自己带比赛的经验下面几类数模问题会频繁出现需要多变量优化牛顿法的场景第一类是最小二乘曲线拟合。你有一堆观测数据点想确定模型里的若干参数让模型输出和实际观测的残差平方和最小。这个残差函数就是一个多变量函数参数个数可能两三个也可能十几个用牛顿法可以直接迭代参数向量。最小二乘问题的目标函数有天然的平方结构Hessian矩阵在很多情形下可以通过Jacobian矩阵近似算起来也不慢。第二类是最大似然估计。比如排队论、可靠性分析中要估计分布参数把似然函数取对数后变成最大化问题求导出来的往往是非线性方程组这正是多变量牛顿法发挥优势的地方。第三类是带约束优化的子问题求解。很多带约束的问题会转成拉格朗日函数求驻点驻点方程是一组多元方程组这时候牛顿法作为一个高效的局部迭代器被嵌套在更大框架里。比赛中你不会直接看到“请用牛顿法”这种直白描述但一旦进入最优策略搜索底层就是这些算法在跑。需要提醒的是牛顿法不能保证收敛到全局最优。它是个局部优化方法对初始点非常敏感。数学建模中更稳妥的做法是先用网格搜索、遗传算法这类全局方法找到一块比较好的区域再在这个区域内用牛顿法做精细局部收敛俗称“全局找凹地牛顿爬到底”。2. 核心原理解析多变量牛顿法的公式是怎么推出来的2.1 直觉起点用二次曲面逼近函数要理解多变量牛顿法我的经验是先别急着把矩阵代数抄一遍而是把它看成一个几何问题。在一维情况下一个函数在某个点附近可以用抛物线来近似不仅知道斜率还知道弯曲程度二阶导所以可以找到抛物线的底部。多变量情形完全类似只不过“斜率”变成了梯度向量“弯曲程度”变成了一块矩阵也就是Hessian矩阵。我们手里有一个n维函数f(x)x ∈ R^n假设它在我们关心的当前点x_k附近是足够光滑的。利用多元泰勒展开舍去三阶及以上的项可以得到f(x) ≈ f(x_k) g_k^T (x - x_k) (1/2)(x - x_k)^T H_k (x - x_k)其中g_k是梯度向量H_k是Hessian矩阵即二阶偏导数矩阵。这个展开式可以理解为我们用了一个“凹陷”或“隆起”的二次曲面去近似原曲面。牛顿法的灵魂就在这一步要知道一个二次曲面在哪个位置达到极小就不用再做复杂的搜索。二次函数的极值有一个闭式解对它求梯度并令其等于零可以直接解出极小点。对上面的近似式关于x求梯度得到∇f(x) ≈ g_k H_k (x - x_k)。极小点满足梯度为零所以令它等于零g_k H_k (x_{k1} - x_k) 0这就是多变量牛顿法的核心方程。整理一下H_k (x_{k1} - x_k) -g_k换成更常见的迭代形式x_{k1} x_k - H_k^{-1} g_k2.2 为什么用Hessian矩阵而不是直接套一维公式有人会问既然一维牛顿法的格式是x_{k1}x_k - f(x_k)/f(x_k)那多变量版是不是用每个方向上的二阶导分别除一下就行这个理解是错的。多变量问题的关键难点在于变量之间会互相耦合。某个方向上不仅有自己的弯曲程度还会被其他方向的变化牵连。Hessian矩阵的非对角元素∂²f/∂x_i ∂x_j恰好刻画了这种“交叉弯曲”。举个例子如果目标函数是f(x,y)x²5xyy²变量x和y并不是独立的。单独看x方向二阶导是2单独看y方向二阶导也是2但混合偏导是5中间那个小山包的方向实际上沿着一个斜轴。如果不考虑交叉项把一个方向上的单变量二阶导拿来直接除迭代路径就会歪掉。用矩阵的语言来说我们求解的是H_k d_k -g_k这个线性方程组的解d_k就是本次迭代的搜索方向。它跟最速下降方向-g_k不一样是在Hessian张成的度量空间里做了重新标定。这正是牛顿法的几何意义它不是简单地沿着当前梯度负方向走而是根据局部的曲率信息修正了这个方向。通俗地讲梯度下降像闭眼往坡下走牛顿法睁着眼看地形的弯曲提前把弯路掰直了。2.3 二次收敛速度为什么它可以几步走完牛顿法在数学上的经典结论是局部二次收敛。假设目标函数足够光滑、在最优点附近Hessian矩阵正定且初始点离最优点足够近那么迭代误差e_k x_k - x*满足‖e_{k1}‖ ≤ C‖e_k‖²这里的C是一个与函数三阶导有关的常数。二次收敛是非常恐怖的收敛速度。这直观上怎么理解如果当前误差是0.1量级平方一下变成0.01再下次0.0001误差会像双重指数一样崩塌下来。对比梯度下降的线性收敛每次误差乘以一个小于1的常数q比如q0.5从0.1误差到1e-9需要大概30次迭代而对牛顿法来说可能只需要两三次。但“局部”二字要画重点。这个收敛速度的前提是初始点已经落在最优点某个邻域内Hessian矩阵在这一区域里还得是正定的。如果初始点离最优点太远或者Hessian矩阵奇异性严重牛顿法照样可能震荡甚至发散。这也是为什么在真实项目里很少直接用“纯净版”牛顿法从任意初始点硬跑而是通常加规则化项或进行线搜索。3. 手算推演二元二次函数用牛顿法迭代的完整过程3.1 先造一个函数看牛顿法怎么做在数学建模问题里我们经常会构造一些便于检验的测试函数。这里用一个二元二次函数f(x, y) 3x² 2xy 4y² - x - 3y这个函数的梯度向量是g [∂f/∂x, ∂f/∂y] [6x 2y - 1, 2x 8y - 3]Hessian矩阵是H [[∂²f/∂x², ∂²f/∂x∂y], [∂²f/∂y∂x, ∂²f/∂y²]] [[6, 2], [2, 8]]这是一个二次函数并且Hessian矩阵是正定的意味着函数是一个向上开口的椭圆抛物面有且仅有一个全局最小值点。把梯度方程组g 0解出来可以得到其解析最优点6x 2y 1 2x 8y 3解得y 4/11x 1/22所以精确极小点是(1/22, 4/11) ≈ (0.0455, 0.3636)。这是后面判断迭代是否正确的标准值。3.2 从初始点(1, 1)出发手算一次迭代现在从初始点x_0 (1, 1)开始。在该点处的梯度为g_0 [6×1 2×1 - 1, 2×1 8×1 - 3] [7, 7]Hessian矩阵是不随位置变化的常数矩阵H [[6, 2], [2, 8]]牛顿法迭代要求解线性方程组 H·d -g也就是6d_x 2d_y -7 2d_x 8d_y -7由第一个式子得d_x (-7 - 2d_y)/6代入第二个式子2×(-7 - 2d_y)/6 8d_y -7 (-7 - 2d_y)/3 8d_y -7 -7/3 - (2/3)d_y 8d_y -7 (22/3)d_y -14/3 d_y -7/11带回d_x的表达式d_x (-7 14/11)/6 (-63/11)/6 -21/22因此搜索方向是d (-21/22, -7/11)。更新位置x_1 x_0 d (1 - 21/22, 1 - 7/11) (1/22, 4/11)一步就到达了理论最优点。这个结果非常漂亮也揭示了一个很重要的性质对于二次目标函数只要Hessian正定且非奇异牛顿法一步就能收敛到精确解不需要迭代第二次。因为二次函数的泰勒展开没有三阶以上的尾项用二次模型去近似它是精确的这一步相当于直接解了极小点的线性方程。这也解释了为什么在最小二乘问题里如果模型是线性的那么“高斯-牛顿法”一次迭代就可以给出闭式解。而实际问题之所以要多次迭代正是因为目标函数不是二次函数舍去的那些高阶项导致需要反复修正。4. 实操用Python实现多变量牛顿法4.1 最小实现纯解析Hessian版本在数模竞赛里代码讲究“短平快”。下面给一个最基础的牛顿法实现模板用解析梯度解析Hessian可以处理任意n维问题import numpy as np def newton_method(f, grad, hess, x0, tol1e-8, max_iter100): x np.array(x0, dtypefloat) path [x.copy()] for i in range(max_iter): g grad(x) if np.linalg.norm(g, ordnp.inf) tol: print(f收敛于第{i}次迭代) break H hess(x) try: d np.linalg.solve(H, -g) except np.linalg.LinAlgError: print(Hessian矩阵奇异开始使用阻尼处理) break # 简单线搜索步长从1开始不满足下降条件就减半 alpha 1.0 while f(x alpha * d) f(x) 1e-4 * alpha * np.dot(g, d): alpha * 0.5 if alpha 1e-10: break x x alpha * d path.append(x.copy()) return x, len(path) def f(v): x, y v return 3*x**2 2*x*y 4*y**2 - x - 3*y def grad(v): x, y v return np.array([6*x 2*y - 1, 2*x 8*y - 3]) def hess(v): return np.array([[6.0, 2.0], [2.0, 8.0]]) x_opt, iters newton_method(f, grad, hess, [1, 1]) print(x_opt, iters)这段代码里有个容易被新手忽略的关键点用np.linalg.solve解线性方程而不是用np.linalg.inv(H)先求逆再乘。原因很实际solve底层用的是LU分解计算量小一个量级数值稳定性也更好。在比赛中数据维度可能不高但求逆操作容易放大舍入误差属于“能不用就不用”的坏习惯。另外一个细节是线搜索。很多教材讲牛顿法都不提步长默认α1。但在非二次函数上大步长经常导致迭代发散。我在代码里加了一个小回退线搜索条件用的是Armijo准则目标函数必须比当前值下降一个正比于梯度方向导数的量。如果在数值上出现新点函数值不降反升就说明α1这一步迈大了果断减半。这个技巧代码量小但对收敛稳定性帮助非常大。4.2 遇到复杂函数用数值Hessian代替解析推导实际建模时目标函数经常是一大堆代码拼出来的梯度都未必能手推更别提二阶偏导。比如你写的目标函数里有循环、有查表、有分段逻辑这时解析求导简直要命。一个现实的替代方案是用有限差分法数值估计梯度和Hessian。梯度可以用中心差分∂f/∂x_i ≈ (f(x h e_i) - f(x - h e_i)) / (2h)Hessian矩阵的对角元和非对角元也可以用二阶差分∂²f/∂x_i² ≈ (f(x h e_i) - 2f(x) f(x - h e_i)) / h² ∂²f/∂x_i∂x_j ≈ (f(x h e_i h e_j) - f(x h e_i - h e_j) - f(x - h e_i h e_j) f(x - h e_i - h e_j)) / (4h²)这里h的选取很讲究。数值微分天生有截断误差和舍入误差的权衡h太大截断误差大h太小浮点数相减导致有效数字丢失舍入误差大。我的经验是取h ≈ eps^(1/3)量级其中eps是机器精度双精度下大约1e-6到1e-5。实际做项目时不会直接从零写数值差分用scipy.optimize里现成的approx_fprime或者直接借用numdifftools之类的库更方便。但拼比赛中如果只依赖第三方库容易被环境限制所以我会在代码包里常备一个numpy版本的数值Hessian函数简单几行心里踏实。4.3 终止条件和容差设计程序判断要不要停一般有三个指标梯度范数、相邻两步之间的位移、目标函数下降量。我推荐看梯度的无穷范数因为最优点处的梯度应该是零向量梯度范数越小离最优点的距离通常也越近。if np.linalg.norm(g, ordnp.inf) tol: print(gradient norm is small enough) break关于tol不同问题差别很大。如果目标是比赛里的拟合建模一般取1e-6就够了再小可能陷入数值噪声区如果目标是想求一个高精度解比如用来判断算法内部逻辑正确性可以取到1e-10但要确保Hessian条件数别太高。条件数很高的时候梯度下降和牛顿法的浮点误差都会被放大这时候单纯收紧容差没有意义。迭代次数上限也要设置我习惯设100到200次。因为牛顿法如果效果好根本走不了几步如果走到几十步还不收敛说明初始点选得不行或Hessian出了问题继续死磕大概率是浪费算力。5. 改进策略Hessian奇异、不正定时的应对方案5.1 牛顿法最怕的Hessian问题理想情况下Hessian矩阵正定牛顿方向是下降方向迭代能顺畅推进。但实际问题永远不按剧本走。可能出现的情况有三种。第一种是Hessian奇异线性方程组解不出来第二种是Hessian非正定有负特征值此时牛顿方向可能根本不是下山方向函数值不降反升第三种是Hessian条件数极大虽能解方程但数值误差大到不可信。这几种情况在数值上就像开车遇到结冰路面刹车失灵又打滑。常规解决思路是“正则化”给Hessian矩阵加一个对角占优的扰动项把矩阵“掰”回正定。最常用的招是H_modified H λI其中λ是一个非负参数I是单位矩阵。这个技巧最早在Levenberg-Marquardt算法里被系统使用实际上就是在牛顿方向和最速下降方向之间做插值。λ取0是纯牛顿法λ取无穷大时方向就趋于负梯度方向——也就是最速下降法。实际使用时如果Hessian检测到负特征值就把λ设成一个正数比如1e-4或者更大的值然后再解方程能显著提升稳定性。还有一点加的对角扰动量应该跟Hessian对角元的量级适配。如果Hessian元素本身在1e3量级却只加1e-6等于没加如果Hessian量级只有1e-6加个1反而把原本的信息都盖住了。比较通用的做法是λ从1e-3量级开始通过循环调整直到矩阵变为正定为止。5.2 阻尼牛顿法步长不是越大越好牛顿法更新式里的“步长”其实已经被Hessian隐式决定了大多数情况下方向是好的但步子通常迈得太猛。很多教材默认α1而实际问题告诉我们α1经常导致越界、震荡、甚至跳到一个函数值更大的区域。阻尼牛顿法的设计很简单就是不改变搜索方向只对步长α做一维搜索。每轮迭代先解牛顿方程拿到方向d_k然后在这一条射线上找使得f(x_k αd_k)最小的α。最简单好用的是回溯线搜索从α1开始如果函数值满足Armijo条件就接受否则α减半重复直到条件满足。我在代码里用的是Armijo条件f(x_k αd_k) ≤ f(x_k) c·α·g_k^T·d_k其中c是常数通常取1e-4。右边这个g_k^T·d_k是一个负数因为下降方向满足g_k^T·d_k 0所以右边表示“接受一个足够明显的下降”。如果这一步下降幅度连很小系数下的线性预估都达不到那就说明步长太大需要往回缩。这个看似简单的改动能让牛顿法在实际应用中脱胎换骨。我测试Rosenbrock函数时纯牛顿法在某些初始点直接飞出坐标范围加上回退线搜索后收敛稳定性和成功率都大幅提升。5.3 从牛顿法到拟牛顿法一劳永逸的扩展思路牛顿法每轮迭代都要组装Hessian矩阵并解线性方程计算量是O(n³)的量级。当变量数达到几百上千时这个成本就有点吃不消了。更重要的是高阶目标函数的Hessian矩阵并不总是容易获得。于是拟牛顿法登场不直接计算Hessian矩阵而是用每一步的梯度和位移信息去逼近Hessian的逆矩阵最著名的是BFGS公式。它的思路是用历史的梯度变化去迭代更新B_k ≈ H_k^{-1}避免每步求矩阵和求逆。数学上拟牛顿法保持了超线性收敛速度比梯度下降快得多又比牛顿法省计算。因为BFGS更新只涉及向量乘法和外积每轮复杂度降到O(n²)。在比赛实战中如果优化问题的变量数量不超过几十个我个人还是推荐直接用牛顿法配合数值Hessian因为逻辑直观、调参方便。一旦维度上了百级就果断切到BFGS或L-BFGS后者在内存受限时尤其有效。可以把牛顿法看成“精确但昂贵的探测器”拟牛顿法看成“高效但近似的替代品”两者在不同场景下各有所长。6. 常见问题与实操经验数模比赛里的牛顿法6.1 典型故障排查速查表迭代不收敛、报错矩阵奇异、函数值越界这些是新手用牛顿法时最常见的问题。下面这表是我在带比赛时整理的排查思路。现象可能原因优先检查项固定解法第1次迭代后函数值暴涨初始点离最优点太远或方向错误打印每次迭代的梯度范数、函数值加入回溯线搜索限制步长Hessian矩阵奇异或接近奇异目标函数在该点平坦或存在冗余参数检查条件数、特征值是否接近0加正则项HλI或改用拟牛顿法收敛到同一个错误点函数多峰初始点落错吸引域画出函数等高线或做多点初始化测试先用全局搜索找近似最优区域再局部收敛迭代很慢几十步不收敛Hessian条件数过大统计Hessian最大最小特征值比值考虑变量缩放或改用BFGS最终结果依赖初始点波动大目标函数非凸对比多组初始点结果用网格/遗传算法确定初值这里我要多说一句变量缩放的问题。实际建模中比如一个变量是温度量级300另一个变量是流量量级0.1这种情况下Hessian矩阵的对角元素会差好几个数量级矩阵条件数天然就大。数值计算里的浮点精度再高也扛不住这种尺度失衡。解决办法是建模时先做无量纲化或变量标准化把变量全部压缩到相近量级这在优化里几乎是必须做的预处理步骤却经常被学生忽略。6.2 比赛中的实战选择到底该不该用牛顿法数模竞赛和纯数值计算研究有个区别竞赛更看重“拿到一个合理结果并解释清楚”而不是追求极端精度。因此在比赛中牛顿法更像是一个“局部精修工具”不要一上来就指望它暴力收敛。我的惯用套路是分两步走。第一步做全局搜索把决策变量的可行域做一次初略网格扫描或跑一个种群算法比如遗传算法得到一个“大概差不多”的候选点。第二步从这个候选点出发用牛顿法去做局部快速收敛。这个组合拳比单纯用哪一个都稳。你换一个初始点局部结果可能会略有差异最后报告里写清楚算法流程、初始点选取方式、收敛容差评卷人就能看出你是在认真做数值实验而不是盲调黑盒。另一个经验是如果能手工写出梯度就别用数值差分梯度如果能写出Hessian那更好。手推公式虽然费时间但能让你对问题结构保持敏感。有一次我带学生做一个交通流参数估计问题目标函数涉及好几个求和项手推Hessian时发现其中一项在某个参数区间内符号会反转导致Hessian非正定进一步发现那个参数本身不可识别。这个结论如果完全依赖数值库很可能就被淹没了。6.3 一个替换思路牛顿法失效时用什么兜底如果做好了以上准备牛顿法还是在某个方向上不断触发异常我的兜底方案是按顺序降级牛顿法 → 阻尼牛顿法 → 高斯-牛顿法 → BFGS → 梯度下降。降级不是丢人而是理性判断。很多时候高斯-牛顿法已经够用了它在最小二乘问题里比标准牛顿法更便宜因为它只用Jacobian矩阵构造Hessian近似不用求二阶偏导。而梯度下降虽然收敛慢但胜在稳健只要步长取得够小几乎永远不会发散。这个降级顺序的核心思想是你永远要有一个“能出结果”的方案。比赛最怕的不是精度不够而是程序在截止时间前跑不出任何结果。我见过太多队伍在一个复杂算法上磨了几个小时最后交上去一份没握到底的代码。灵活切换算法比固执追求单一方法更符合实战需求。最后分享一个个人习惯。每次写牛顿法代码我都会顺手写一个测试函数——一个简单的二次函数和一个Rosenbrock函数跑通后再套真实数据。这个习惯帮我避开了很多因为索引错误、维度不匹配导致的低级bug。多变量牛顿法其实不难难的是在拟合、估计、优化这些实际场景中搞清楚它什么时候灵、什么时候不灵。如果你能把这篇笔记对应的推导流程亲手走一遍再跑通一个二维例子的代码后面再遇到所谓的“多变量最优化计算”你会发现自己已经不像之前那样没底了。 
长河君 · 2026/10/4 9:50:23 · 阅读 1.2k 多变量牛顿法全解析：从Hessian矩阵到Python实现与避坑指南
数学建模的题目做到后期十有八九会撞到优化问题。曲线拟合要调参数资源分配要定方案运筹调度要找最优这些本质都是“在某个目标函数上求极值”。真到了多变量场景很多教材翻来覆去只讲梯度下降讲到牛顿法也就给个一维公式等到Hessian矩阵一出现大家基本就绕道走了。其实多变量最优化计算里的牛顿法是理解数值优化很关键的一环理解了它再看拟牛顿法、信赖域法都会顺畅很多。这篇数模笔记主要整理牛顿法在多变量优化里的完整推导、手算过程、Python实现以及实操中Hessian矩阵“翻车”时的补救方案。建议正在准备数学建模竞赛、算法竞赛或者刚接触数值优化计算的同学重点收藏这套笔记是按“公式推导—手算验证—代码实现—避坑总结”的流程走的可以直接照着用。1. 项目概述多变量牛顿法是什么能解决哪些问题1.1 从“求根”到“求极值”牛顿法的原始身份很多人第一眼看到牛顿法是在数值分析里解非线性方程的部分所谓Newton-Raphson迭代。当时的场景很简单给定一个一元函数f(x)我们想找到某个x*使f(x*) 0。做法是在当前点xk处用切线近似f(x)然后求切线与x轴的交点x_{k1} x_k - f(x_k) / f(x_k)这个迭代式看起来平平无奇但它的思想极其朴素任何光滑函数都可以在某一点的邻域内用低阶多项式逼近而在局部上用一次函数切线去逼近就已经能给出不错的下一步位置了。到了最优化计算里我们的目标变了不是找f(x) 0而是找f(x)取得极小值的点x*。但两者在数学上可以打一个“等号”如果f(x)可微那么在极值点处梯度g(x) ∇f(x) 0。这等于说优化问题可以转化为一个“梯度方程组求根”的问题。把牛顿法从“求f(x)0”迁移到“求∇f(x)0”就引出了多变量牛顿法。这也解释了为什么多变量牛顿法在收敛速度上有先天优势。它天生是冲着方程的根去的而极值点恰好是梯度向量的零点所以它对极值点的逼近是“二阶”级别的。梯度下降法只用了目标函数的一阶信息牛顿法把二阶信息Hessian矩阵也用了收敛快是自然的事。1.2 数学建模里哪些场景会用到它结合我自己带比赛的经验下面几类数模问题会频繁出现需要多变量优化牛顿法的场景第一类是最小二乘曲线拟合。你有一堆观测数据点想确定模型里的若干参数让模型输出和实际观测的残差平方和最小。这个残差函数就是一个多变量函数参数个数可能两三个也可能十几个用牛顿法可以直接迭代参数向量。最小二乘问题的目标函数有天然的平方结构Hessian矩阵在很多情形下可以通过Jacobian矩阵近似算起来也不慢。第二类是最大似然估计。比如排队论、可靠性分析中要估计分布参数把似然函数取对数后变成最大化问题求导出来的往往是非线性方程组这正是多变量牛顿法发挥优势的地方。第三类是带约束优化的子问题求解。很多带约束的问题会转成拉格朗日函数求驻点驻点方程是一组多元方程组这时候牛顿法作为一个高效的局部迭代器被嵌套在更大框架里。比赛中你不会直接看到“请用牛顿法”这种直白描述但一旦进入最优策略搜索底层就是这些算法在跑。需要提醒的是牛顿法不能保证收敛到全局最优。它是个局部优化方法对初始点非常敏感。数学建模中更稳妥的做法是先用网格搜索、遗传算法这类全局方法找到一块比较好的区域再在这个区域内用牛顿法做精细局部收敛俗称“全局找凹地牛顿爬到底”。2. 核心原理解析多变量牛顿法的公式是怎么推出来的2.1 直觉起点用二次曲面逼近函数要理解多变量牛顿法我的经验是先别急着把矩阵代数抄一遍而是把它看成一个几何问题。在一维情况下一个函数在某个点附近可以用抛物线来近似不仅知道斜率还知道弯曲程度二阶导所以可以找到抛物线的底部。多变量情形完全类似只不过“斜率”变成了梯度向量“弯曲程度”变成了一块矩阵也就是Hessian矩阵。我们手里有一个n维函数f(x)x ∈ R^n假设它在我们关心的当前点x_k附近是足够光滑的。利用多元泰勒展开舍去三阶及以上的项可以得到f(x) ≈ f(x_k) g_k^T (x - x_k) (1/2)(x - x_k)^T H_k (x - x_k)其中g_k是梯度向量H_k是Hessian矩阵即二阶偏导数矩阵。这个展开式可以理解为我们用了一个“凹陷”或“隆起”的二次曲面去近似原曲面。牛顿法的灵魂就在这一步要知道一个二次曲面在哪个位置达到极小就不用再做复杂的搜索。二次函数的极值有一个闭式解对它求梯度并令其等于零可以直接解出极小点。对上面的近似式关于x求梯度得到∇f(x) ≈ g_k H_k (x - x_k)。极小点满足梯度为零所以令它等于零g_k H_k (x_{k1} - x_k) 0这就是多变量牛顿法的核心方程。整理一下H_k (x_{k1} - x_k) -g_k换成更常见的迭代形式x_{k1} x_k - H_k^{-1} g_k2.2 为什么用Hessian矩阵而不是直接套一维公式有人会问既然一维牛顿法的格式是x_{k1}x_k - f(x_k)/f(x_k)那多变量版是不是用每个方向上的二阶导分别除一下就行这个理解是错的。多变量问题的关键难点在于变量之间会互相耦合。某个方向上不仅有自己的弯曲程度还会被其他方向的变化牵连。Hessian矩阵的非对角元素∂²f/∂x_i ∂x_j恰好刻画了这种“交叉弯曲”。举个例子如果目标函数是f(x,y)x²5xyy²变量x和y并不是独立的。单独看x方向二阶导是2单独看y方向二阶导也是2但混合偏导是5中间那个小山包的方向实际上沿着一个斜轴。如果不考虑交叉项把一个方向上的单变量二阶导拿来直接除迭代路径就会歪掉。用矩阵的语言来说我们求解的是H_k d_k -g_k这个线性方程组的解d_k就是本次迭代的搜索方向。它跟最速下降方向-g_k不一样是在Hessian张成的度量空间里做了重新标定。这正是牛顿法的几何意义它不是简单地沿着当前梯度负方向走而是根据局部的曲率信息修正了这个方向。通俗地讲梯度下降像闭眼往坡下走牛顿法睁着眼看地形的弯曲提前把弯路掰直了。2.3 二次收敛速度为什么它可以几步走完牛顿法在数学上的经典结论是局部二次收敛。假设目标函数足够光滑、在最优点附近Hessian矩阵正定且初始点离最优点足够近那么迭代误差e_k x_k - x*满足‖e_{k1}‖ ≤ C‖e_k‖²这里的C是一个与函数三阶导有关的常数。二次收敛是非常恐怖的收敛速度。这直观上怎么理解如果当前误差是0.1量级平方一下变成0.01再下次0.0001误差会像双重指数一样崩塌下来。对比梯度下降的线性收敛每次误差乘以一个小于1的常数q比如q0.5从0.1误差到1e-9需要大概30次迭代而对牛顿法来说可能只需要两三次。但“局部”二字要画重点。这个收敛速度的前提是初始点已经落在最优点某个邻域内Hessian矩阵在这一区域里还得是正定的。如果初始点离最优点太远或者Hessian矩阵奇异性严重牛顿法照样可能震荡甚至发散。这也是为什么在真实项目里很少直接用“纯净版”牛顿法从任意初始点硬跑而是通常加规则化项或进行线搜索。3. 手算推演二元二次函数用牛顿法迭代的完整过程3.1 先造一个函数看牛顿法怎么做在数学建模问题里我们经常会构造一些便于检验的测试函数。这里用一个二元二次函数f(x, y) 3x² 2xy 4y² - x - 3y这个函数的梯度向量是g [∂f/∂x, ∂f/∂y] [6x 2y - 1, 2x 8y - 3]Hessian矩阵是H [[∂²f/∂x², ∂²f/∂x∂y], [∂²f/∂y∂x, ∂²f/∂y²]] [[6, 2], [2, 8]]这是一个二次函数并且Hessian矩阵是正定的意味着函数是一个向上开口的椭圆抛物面有且仅有一个全局最小值点。把梯度方程组g 0解出来可以得到其解析最优点6x 2y 1 2x 8y 3解得y 4/11x 1/22所以精确极小点是(1/22, 4/11) ≈ (0.0455, 0.3636)。这是后面判断迭代是否正确的标准值。3.2 从初始点(1, 1)出发手算一次迭代现在从初始点x_0 (1, 1)开始。在该点处的梯度为g_0 [6×1 2×1 - 1, 2×1 8×1 - 3] [7, 7]Hessian矩阵是不随位置变化的常数矩阵H [[6, 2], [2, 8]]牛顿法迭代要求解线性方程组 H·d -g也就是6d_x 2d_y -7 2d_x 8d_y -7由第一个式子得d_x (-7 - 2d_y)/6代入第二个式子2×(-7 - 2d_y)/6 8d_y -7 (-7 - 2d_y)/3 8d_y -7 -7/3 - (2/3)d_y 8d_y -7 (22/3)d_y -14/3 d_y -7/11带回d_x的表达式d_x (-7 14/11)/6 (-63/11)/6 -21/22因此搜索方向是d (-21/22, -7/11)。更新位置x_1 x_0 d (1 - 21/22, 1 - 7/11) (1/22, 4/11)一步就到达了理论最优点。这个结果非常漂亮也揭示了一个很重要的性质对于二次目标函数只要Hessian正定且非奇异牛顿法一步就能收敛到精确解不需要迭代第二次。因为二次函数的泰勒展开没有三阶以上的尾项用二次模型去近似它是精确的这一步相当于直接解了极小点的线性方程。这也解释了为什么在最小二乘问题里如果模型是线性的那么“高斯-牛顿法”一次迭代就可以给出闭式解。而实际问题之所以要多次迭代正是因为目标函数不是二次函数舍去的那些高阶项导致需要反复修正。4. 实操用Python实现多变量牛顿法4.1 最小实现纯解析Hessian版本在数模竞赛里代码讲究“短平快”。下面给一个最基础的牛顿法实现模板用解析梯度解析Hessian可以处理任意n维问题import numpy as np def newton_method(f, grad, hess, x0, tol1e-8, max_iter100): x np.array(x0, dtypefloat) path [x.copy()] for i in range(max_iter): g grad(x) if np.linalg.norm(g, ordnp.inf) tol: print(f收敛于第{i}次迭代) break H hess(x) try: d np.linalg.solve(H, -g) except np.linalg.LinAlgError: print(Hessian矩阵奇异开始使用阻尼处理) break # 简单线搜索步长从1开始不满足下降条件就减半 alpha 1.0 while f(x alpha * d) f(x) 1e-4 * alpha * np.dot(g, d): alpha * 0.5 if alpha 1e-10: break x x alpha * d path.append(x.copy()) return x, len(path) def f(v): x, y v return 3*x**2 2*x*y 4*y**2 - x - 3*y def grad(v): x, y v return np.array([6*x 2*y - 1, 2*x 8*y - 3]) def hess(v): return np.array([[6.0, 2.0], [2.0, 8.0]]) x_opt, iters newton_method(f, grad, hess, [1, 1]) print(x_opt, iters)这段代码里有个容易被新手忽略的关键点用np.linalg.solve解线性方程而不是用np.linalg.inv(H)先求逆再乘。原因很实际solve底层用的是LU分解计算量小一个量级数值稳定性也更好。在比赛中数据维度可能不高但求逆操作容易放大舍入误差属于“能不用就不用”的坏习惯。另外一个细节是线搜索。很多教材讲牛顿法都不提步长默认α1。但在非二次函数上大步长经常导致迭代发散。我在代码里加了一个小回退线搜索条件用的是Armijo准则目标函数必须比当前值下降一个正比于梯度方向导数的量。如果在数值上出现新点函数值不降反升就说明α1这一步迈大了果断减半。这个技巧代码量小但对收敛稳定性帮助非常大。4.2 遇到复杂函数用数值Hessian代替解析推导实际建模时目标函数经常是一大堆代码拼出来的梯度都未必能手推更别提二阶偏导。比如你写的目标函数里有循环、有查表、有分段逻辑这时解析求导简直要命。一个现实的替代方案是用有限差分法数值估计梯度和Hessian。梯度可以用中心差分∂f/∂x_i ≈ (f(x h e_i) - f(x - h e_i)) / (2h)Hessian矩阵的对角元和非对角元也可以用二阶差分∂²f/∂x_i² ≈ (f(x h e_i) - 2f(x) f(x - h e_i)) / h² ∂²f/∂x_i∂x_j ≈ (f(x h e_i h e_j) - f(x h e_i - h e_j) - f(x - h e_i h e_j) f(x - h e_i - h e_j)) / (4h²)这里h的选取很讲究。数值微分天生有截断误差和舍入误差的权衡h太大截断误差大h太小浮点数相减导致有效数字丢失舍入误差大。我的经验是取h ≈ eps^(1/3)量级其中eps是机器精度双精度下大约1e-6到1e-5。实际做项目时不会直接从零写数值差分用scipy.optimize里现成的approx_fprime或者直接借用numdifftools之类的库更方便。但拼比赛中如果只依赖第三方库容易被环境限制所以我会在代码包里常备一个numpy版本的数值Hessian函数简单几行心里踏实。4.3 终止条件和容差设计程序判断要不要停一般有三个指标梯度范数、相邻两步之间的位移、目标函数下降量。我推荐看梯度的无穷范数因为最优点处的梯度应该是零向量梯度范数越小离最优点的距离通常也越近。if np.linalg.norm(g, ordnp.inf) tol: print(gradient norm is small enough) break关于tol不同问题差别很大。如果目标是比赛里的拟合建模一般取1e-6就够了再小可能陷入数值噪声区如果目标是想求一个高精度解比如用来判断算法内部逻辑正确性可以取到1e-10但要确保Hessian条件数别太高。条件数很高的时候梯度下降和牛顿法的浮点误差都会被放大这时候单纯收紧容差没有意义。迭代次数上限也要设置我习惯设100到200次。因为牛顿法如果效果好根本走不了几步如果走到几十步还不收敛说明初始点选得不行或Hessian出了问题继续死磕大概率是浪费算力。5. 改进策略Hessian奇异、不正定时的应对方案5.1 牛顿法最怕的Hessian问题理想情况下Hessian矩阵正定牛顿方向是下降方向迭代能顺畅推进。但实际问题永远不按剧本走。可能出现的情况有三种。第一种是Hessian奇异线性方程组解不出来第二种是Hessian非正定有负特征值此时牛顿方向可能根本不是下山方向函数值不降反升第三种是Hessian条件数极大虽能解方程但数值误差大到不可信。这几种情况在数值上就像开车遇到结冰路面刹车失灵又打滑。常规解决思路是“正则化”给Hessian矩阵加一个对角占优的扰动项把矩阵“掰”回正定。最常用的招是H_modified H λI其中λ是一个非负参数I是单位矩阵。这个技巧最早在Levenberg-Marquardt算法里被系统使用实际上就是在牛顿方向和最速下降方向之间做插值。λ取0是纯牛顿法λ取无穷大时方向就趋于负梯度方向——也就是最速下降法。实际使用时如果Hessian检测到负特征值就把λ设成一个正数比如1e-4或者更大的值然后再解方程能显著提升稳定性。还有一点加的对角扰动量应该跟Hessian对角元的量级适配。如果Hessian元素本身在1e3量级却只加1e-6等于没加如果Hessian量级只有1e-6加个1反而把原本的信息都盖住了。比较通用的做法是λ从1e-3量级开始通过循环调整直到矩阵变为正定为止。5.2 阻尼牛顿法步长不是越大越好牛顿法更新式里的“步长”其实已经被Hessian隐式决定了大多数情况下方向是好的但步子通常迈得太猛。很多教材默认α1而实际问题告诉我们α1经常导致越界、震荡、甚至跳到一个函数值更大的区域。阻尼牛顿法的设计很简单就是不改变搜索方向只对步长α做一维搜索。每轮迭代先解牛顿方程拿到方向d_k然后在这一条射线上找使得f(x_k αd_k)最小的α。最简单好用的是回溯线搜索从α1开始如果函数值满足Armijo条件就接受否则α减半重复直到条件满足。我在代码里用的是Armijo条件f(x_k αd_k) ≤ f(x_k) c·α·g_k^T·d_k其中c是常数通常取1e-4。右边这个g_k^T·d_k是一个负数因为下降方向满足g_k^T·d_k 0所以右边表示“接受一个足够明显的下降”。如果这一步下降幅度连很小系数下的线性预估都达不到那就说明步长太大需要往回缩。这个看似简单的改动能让牛顿法在实际应用中脱胎换骨。我测试Rosenbrock函数时纯牛顿法在某些初始点直接飞出坐标范围加上回退线搜索后收敛稳定性和成功率都大幅提升。5.3 从牛顿法到拟牛顿法一劳永逸的扩展思路牛顿法每轮迭代都要组装Hessian矩阵并解线性方程计算量是O(n³)的量级。当变量数达到几百上千时这个成本就有点吃不消了。更重要的是高阶目标函数的Hessian矩阵并不总是容易获得。于是拟牛顿法登场不直接计算Hessian矩阵而是用每一步的梯度和位移信息去逼近Hessian的逆矩阵最著名的是BFGS公式。它的思路是用历史的梯度变化去迭代更新B_k ≈ H_k^{-1}避免每步求矩阵和求逆。数学上拟牛顿法保持了超线性收敛速度比梯度下降快得多又比牛顿法省计算。因为BFGS更新只涉及向量乘法和外积每轮复杂度降到O(n²)。在比赛实战中如果优化问题的变量数量不超过几十个我个人还是推荐直接用牛顿法配合数值Hessian因为逻辑直观、调参方便。一旦维度上了百级就果断切到BFGS或L-BFGS后者在内存受限时尤其有效。可以把牛顿法看成“精确但昂贵的探测器”拟牛顿法看成“高效但近似的替代品”两者在不同场景下各有所长。6. 常见问题与实操经验数模比赛里的牛顿法6.1 典型故障排查速查表迭代不收敛、报错矩阵奇异、函数值越界这些是新手用牛顿法时最常见的问题。下面这表是我在带比赛时整理的排查思路。现象可能原因优先检查项固定解法第1次迭代后函数值暴涨初始点离最优点太远或方向错误打印每次迭代的梯度范数、函数值加入回溯线搜索限制步长Hessian矩阵奇异或接近奇异目标函数在该点平坦或存在冗余参数检查条件数、特征值是否接近0加正则项HλI或改用拟牛顿法收敛到同一个错误点函数多峰初始点落错吸引域画出函数等高线或做多点初始化测试先用全局搜索找近似最优区域再局部收敛迭代很慢几十步不收敛Hessian条件数过大统计Hessian最大最小特征值比值考虑变量缩放或改用BFGS最终结果依赖初始点波动大目标函数非凸对比多组初始点结果用网格/遗传算法确定初值这里我要多说一句变量缩放的问题。实际建模中比如一个变量是温度量级300另一个变量是流量量级0.1这种情况下Hessian矩阵的对角元素会差好几个数量级矩阵条件数天然就大。数值计算里的浮点精度再高也扛不住这种尺度失衡。解决办法是建模时先做无量纲化或变量标准化把变量全部压缩到相近量级这在优化里几乎是必须做的预处理步骤却经常被学生忽略。6.2 比赛中的实战选择到底该不该用牛顿法数模竞赛和纯数值计算研究有个区别竞赛更看重“拿到一个合理结果并解释清楚”而不是追求极端精度。因此在比赛中牛顿法更像是一个“局部精修工具”不要一上来就指望它暴力收敛。我的惯用套路是分两步走。第一步做全局搜索把决策变量的可行域做一次初略网格扫描或跑一个种群算法比如遗传算法得到一个“大概差不多”的候选点。第二步从这个候选点出发用牛顿法去做局部快速收敛。这个组合拳比单纯用哪一个都稳。你换一个初始点局部结果可能会略有差异最后报告里写清楚算法流程、初始点选取方式、收敛容差评卷人就能看出你是在认真做数值实验而不是盲调黑盒。另一个经验是如果能手工写出梯度就别用数值差分梯度如果能写出Hessian那更好。手推公式虽然费时间但能让你对问题结构保持敏感。有一次我带学生做一个交通流参数估计问题目标函数涉及好几个求和项手推Hessian时发现其中一项在某个参数区间内符号会反转导致Hessian非正定进一步发现那个参数本身不可识别。这个结论如果完全依赖数值库很可能就被淹没了。6.3 一个替换思路牛顿法失效时用什么兜底如果做好了以上准备牛顿法还是在某个方向上不断触发异常我的兜底方案是按顺序降级牛顿法 → 阻尼牛顿法 → 高斯-牛顿法 → BFGS → 梯度下降。降级不是丢人而是理性判断。很多时候高斯-牛顿法已经够用了它在最小二乘问题里比标准牛顿法更便宜因为它只用Jacobian矩阵构造Hessian近似不用求二阶偏导。而梯度下降虽然收敛慢但胜在稳健只要步长取得够小几乎永远不会发散。这个降级顺序的核心思想是你永远要有一个“能出结果”的方案。比赛最怕的不是精度不够而是程序在截止时间前跑不出任何结果。我见过太多队伍在一个复杂算法上磨了几个小时最后交上去一份没握到底的代码。灵活切换算法比固执追求单一方法更符合实战需求。最后分享一个个人习惯。每次写牛顿法代码我都会顺手写一个测试函数——一个简单的二次函数和一个Rosenbrock函数跑通后再套真实数据。这个习惯帮我避开了很多因为索引错误、维度不匹配导致的低级bug。多变量牛顿法其实不难难的是在拟合、估计、优化这些实际场景中搞清楚它什么时候灵、什么时候不灵。如果你能把这篇笔记对应的推导流程亲手走一遍再跑通一个二维例子的代码后面再遇到所谓的“多变量最优化计算”你会发现自己已经不像之前那样没底了。