如果你在数学建模竞赛里解过带非线性目标函数的优化问题大概率经历过两种痛苦用梯度下降吧步长调不好收敛慢得让人怀疑人生用现成的优化库吧又不知道它内部到底在做什么论文里写算法部分心里发虚。我第一次在国赛里硬啃一个参数估计问题目标函数是非线性最小二乘的形式梯度下降跑了三百多步还在原地打转后来老老实实把多变量最优化计算里的牛顿法补上才觉得这块真正通了一点。这篇笔记不打算写成教科书式的推导罗列就按我自己从单变量到多变量、从理论到代码、从踩坑到改进的真实过程来写涉及的核心关键词是所有优化方法里绕不开的牛顿法希望对准备数模或者做工程优化的朋友有直接的参考价值。1. 我为什么在数模里最终选了牛顿法——它到底在算什么1.1 一个容易被忽略的前提牛顿法解的其实是个方程先说个特别容易混淆的点。很多人一看到单变量牛顿法第一反应是高中课本里那个“用切线逐步逼近零点”的公式$$x_{k1} x_k - \frac{f(x_k)}{f(x_k)}$$这个公式是拿来解方程 f(x)0 的不是直接求极值。那多变量最优化里怎么用它关键在于求一个光滑函数的极小值点本质上是求梯度为零的点。也就是说我们要解的是 ∇f(x)0 这个方程组而不是直接在 f(x) 上做牛顿迭代。这个转换看起来只是换了个函数实际上影响巨大。你要处理的方程个数和变量个数一样多迭代公式里的分母也变成了一整个矩阵——二阶偏导组成的海森矩阵 H。之前很多教程把这步跳过了直接丢出多维迭代公式读者看得一头雾水其实原因就在这里牛顿法求极值逻辑内核是“求梯度的零点”。1.2 从切线到曲面单变量的直觉怎么推广到多维单变量的几何直觉是在当前位置放一条切线让切线延伸到和 x 轴相交把交点作为下一步。多变量情况下目标函数是一个 n 维曲面你不可能放一条线你应该放一个弯曲的曲面。什么曲面最合适二次曲面因为任何一个光滑函数在足够小的邻域里都能用二次函数去逼近。所以多变量牛顿法干的事情是每一步都用原函数在这个点的二次近似模型来代替原函数然后直接跳到这个二次模型的极值点。如果原函数本身就是二次函数一步就到如果不太像二次函数那就多迭代几次每一次都在新的点上重新拟合一个二次曲面。这比梯度下降聪明的地方在于梯度下降只用了当前位置的坡度信息一阶导数它知道哪个方向是下坡但不知道坡面有多弯曲所以只能拿一个固定的学习率小心翼翼地走。牛顿法用了坡度加曲率信息二阶导数它不仅能判断方向还能根据曲率预估出“再往前走多远就应该到谷底”于是它可以大胆地一次性跨出一大步。1.3 二次收敛是什么概念我当年第一次看到“牛顿法具有二次收敛性”这个表述时毫无感觉直到自己在 R 语言里跑了一组误差数据才真正被震撼。把误差记为 e_k二次收敛意味着误差按这样的规律缩小$$e_{k1} \approx C e_k^2$$假设 C 不是特别离谱误差序列可能是 0.1 → 0.001 → 10^(-8) → 10^(-16)。每一步的误差大约是上一步的平方。也就是说一旦进入牛顿法的收敛邻域精度提升是指数级的三到五步就能把误差压到双精度浮点的极限。对比梯度下降的线性收敛e_{k1}≈C e_k每一步误差只缩一个固定比例要磨很久才能到高精度。所以如果你在数模论文里需要展示算法效率把牛顿法和梯度下降的迭代次数放在同一张表里效果非常直观评审老师看一眼就明白你选的为什么是牛顿法。2. 多变量牛顿法的数学推导每一步都不是魔法2.1 从泰勒展开到驻点方程假设目标函数 f: R^n → R在当前迭代点 x_k 附近足够光滑。我们把它做二阶泰勒展开$$f(x) \approx f(x_k) g_k^{\top}(x - x_k) \frac{1}{2}(x - x_k)^{\top} H_k (x - x_k)$$其中g_k ∇f(x_k) 是梯度向量n 维H_k ∇²f(x_k) 是海森矩阵n×n 对称矩阵第 i 行第 j 列的元素就是 ∂²f/∂x_i ∂x_j 在 x_k 处的值。现在对这个二次模型求梯度并让它等于零。对向量 x 求梯度得到的是一个线性方程组$$g_k H_k (x - x_k) 0$$如果 H_k 是非奇异的就可以解出下一步的迭代点$$x_{k1} x_k - H_k^{-1} g_k$$整个多变量最优化计算里最重要的公式其实就是这一条。你看它没有任何魔法就是从“梯度为零”这个极值必要条件出发加上二阶近似推导出来的结果。2.2 为什么迭代方向是 H^{-1}g 而不是 g这里有个常见的疑问为什么牛顿方向是 H_k^{-1} g_k而不是直接用梯度方向直观理解是梯度方向只告诉你“往下坡走”但没有考虑每个方向上的曲率。海森矩阵编码的是局部曲率信息用它的逆矩阵去左乘梯度相当于对坐标做了一次与曲率成反比的伸缩变换。我经常用一个类比你在一个狭长的山谷里谷底的走向是斜的。梯度方向可能指向谷壁如果在梯度方向上走大步会撞到山坡而海森矩阵知道了“垂直谷底方向曲率很大”它会把这一步幅度压小把主要步长放到“沿着谷底那个曲率小的方向”上去。这就是为什么二次函数上一牛顿法一步到位——因为它把曲率信息全部利用了相当于直接把谷底抛物线方程解出来了。2.3 手算一个二维例子的海森矩阵以最经典的优化测试函数 Rosenbrock 函数为例$$f(x_1, x_2) (1-x_1)^2 100(x_2 - x_1^2)^2$$它的梯度是$$\frac{\partial f}{\partial x_1} -2(1-x_1) - 400x_1(x_2 - x_1^2)$$$$\frac{\partial f}{\partial x_2} 200(x_2 - x_1^2)$$海森矩阵的四个分量$$\frac{\partial^2 f}{\partial x_1^2} 2 - 400(x_2 - x_1^2) 800x_1^2$$$$\frac{\partial^2 f}{\partial x_1 \partial x_2} -400x_1$$$$\frac{\partial^2 f}{\partial x_2 \partial x_1} -400x_1$$$$\frac{\partial^2 f}{\partial x_2^2} 200$$注意海森矩阵是对称的两个混合偏导相等。这是检验你手算是否出错的一个天然手段。我建议在纸面推导时每次都检查 H_12 和 H_21 是否一致如果不等八成是某个链式法则求错了。3. 手写一个能用的多变量牛顿法代码与循环条件3.1 迭代里用“解方程组”而不是“求逆”理论上迭代公式是 x_{k1} x_k - H_k^{-1} g_k但实际代码里千万不要直接去计算矩阵的逆。原因有两个第一求逆的计算量是 O(n^3)而且常数很大第二数值上求逆再乘向量误差比直接解线性方程组大得多。正确的做法是引入方向向量 d_k解这个线性方程组$$H_k d_k -g_k$$然后令 x_{k1} x_k d_k。在 Python 里调用 numpy 的np.linalg.solve一行就能搞定它内部用 LU 分解速度快而且数值稳定性更好。这个细节在数模论文里的算法描述环节值得写进去显得你真正实现过而不是抄了一遍公式。3.2 一个完整的 Python 实现我用 Rosenbrock 函数做测试手写一个纯正的多变量牛顿法。这个函数的最优解是 (1, 1)最优值是 0公认的优化测试函数用来验证算法很合适。import numpy as np def rosen(x): x1, x2 x return (1 - x1)**2 100.0 * (x2 - x1**2)**2 def rosen_grad(x): x1, x2 x g np.zeros(2) g[0] -2.0 * (1 - x1) - 400.0 * x1 * (x2 - x1**2) g[1] 200.0 * (x2 - x1**2) return g def rosen_hess(x): x1, x2 x h11 2.0 - 400.0 * (x2 - x1**2) 800.0 * x1**2 h12 -400.0 * x1 h22 200.0 return np.array([[h11, h12], [h12, h22]]) def newton_minimize(grad, hess, x0, tol1e-8, max_iter100): x np.array(x0, dtypefloat) for i in range(max_iter): g grad(x) H hess(x) grad_norm np.linalg.norm(g, ordnp.inf) print(fiter {i:3d}: x ({x[0]:.10f}, {x[1]:.10f}), ||g||inf {grad_norm:.3e}) if grad_norm tol: return x, i, grad_norm # 关键步骤解 H d -g而不是求逆 d np.linalg.solve(H, -g) x x d return x, max_iter, np.linalg.norm(grad(x), ordnp.inf) x_opt, it, gn newton_minimize(rosen_grad, rosen_hess, [-1.2, 1.0]) print(f\n最优解: {x_opt}) print(f迭代次数: {it}) print(f最终梯度无穷范数: {gn:.3e})我第一次跑这段代码的时候从初始点 (-1.2, 1.0) 出发大概只迭代了五六步就到了 (1, 1) 附近。这就是二次收敛的威力。如果用梯度下降同样的初始点、同样的精度通常需要上千步。3.3 停机条件怎么设才不翻车写优化算法停机条件是特别容易被忽略但特别重要的部分。我自己常用的有三条优先级从高到低梯度范数足够小||g||_inf tol。这是最优性的一阶必要条件无约束光滑问题里它趋近于零说明已经接近驻点。tol 通常取 1e-6 到 1e-8具体看你的目标函数数值尺度。如果函数值本身很大梯度范数阈值也相应放大一点。迭代点变化量足够小||x_{k1} - x_k|| tol。这个条件能防止梯度范数卡在一个不小的值上但迭代点已经动不了的情况——虽然这不常见但碰到海森矩阵接近奇异时会遇到。最大迭代次数兜底max_iter。纯牛顿法有时候会震荡甚至发散必须设一个上限防止死循环。三条都用不要省。尤其第二条我见过不少人在梯度过零点附近来回跳动但梯度范数始终降不到 1e-8 以下的情况没有最大迭代数限制程序就卡死了。4. 牛顿法翻车实录海森矩阵奇异、不下降、初值敏感4.1 最直接的崩溃海森矩阵奇异解不出来多变量牛顿法看起来很美但现实里第一个坑就是海森矩阵可能奇异。什么叫奇异就是行列式为零方程组 H d -g 没有唯一解。几何上对应着当前点所在的局部二次模型是个“平的”或“有个方向没有曲率”的曲面。这种情况在球面型函数或者一些对称性较强的函数里很容易发生。比如 f(x)x1^4 x2^4在某些初始点附近海森矩阵可能是奇异的。代码里表现为np.linalg.solve抛出LinAlgError。如果你只是想要一个能跑的成绩最简单的方式是给海森矩阵加上一个小的单位矩阵的倍数把迭代公式改成$$(H_k \lambda I) d_k -g_k$$\lambda 取一个很小的正数比如 1e-6。这个操作本质上是保证矩阵正定也是后面要讲的阻尼牛顿法和 LM 算法的核心思想。4.2 牛顿方向不一定下降这个坑藏得更深。牛顿法推导时默认了海森矩阵是正定的但实际函数在海森正定以外的区域牛顿方向 d -H^{-1} g 可能与真正的梯度下降方向夹角超过 90 度也就是说沿着这个方向走目标函数不仅不下降还会上升。原因在于海森矩阵非正定时局部二次模型不是“碗型”而是“马鞍形”二次模型的极值点可能是一个鞍点甚至极大值点你跳过去自然不下降。判断方法很简单每次算出 d 之后检查一下 g^T d 是否小于零。如果 g^T d 0说明 d 不是下降方向这时候就不要固执地用牛顿方向了否则数值会瞬间爆炸。4.3 初值敏感一步跳飞牛顿法在局部收敛邻域内是神出了邻域直接变鬼。Rosenbrock 函数还算友好从 (-1.2, 1.0) 出发没问题。但如果你从更远的地方出发比如 (10, 10)二次模型完全看不出真实谷底在哪一步迭代可能直接跳到函数值大到不可描述的区域。我在自己的试验里见过最典型的情况初始点选得不合适第一步 d 的长度超过 1e4第二步目标函数值变成 10^200程序当场溢出。解决思路有两个一是用多点随机初始化比如均匀撒 50 个初始点每个点跑一遍牛顿法收敛后再看极值分布二是先用梯度下降或单纯形法粗跑几十步得到一个较好的初值再切换到牛顿法收尾。这种“粗糙全局精细局部”的组合拳在数模实战里非常有效。4.4 计算成本与数值误差多变量牛顿法每次迭代都要计算海森矩阵并求解一个 n 阶线性方程组复杂度 O(n^3)。当变量个数 n 到几千上万时这个代价大到基本不可用。所以牛顿法在小规模问题上表现优秀大规模问题上基本没有优势这个一定要在选型时想清楚。数值误差方面海森矩阵如果是手动推导的非常容易出错。我的经验是无论用链式法则手推时觉得多简单都要利用“海森矩阵必须对称”这条性质检查一遍。另外可以用数值差分法验证梯度和海森的正确性。所谓数值差分就是让某个变量动一个小量 ε用 (f(xεe_i) - f(x-εe_i)) / (2ε) 近似偏导数然后和解析表达式对比。这个方法虽然跑起来慢但只用来测试正确性一锤定音。5. 从牛顿法出发的改进路线阻尼、高斯-牛顿、拟牛顿怎么选5.1 阻尼牛顿法一行代码救回大半场景理解了“牛顿方向不下降”的问题之后最自然的修复方案就是给牛顿方向加一个步长因子变成$$x_{k1} x_k \alpha_k d_k$$其中 \alpha_k 通过一维线搜索确定比如 Armijo 准则从 \alpha1 开始不断减半直到目标函数充分下降。这个方法叫阻尼牛顿法本质是保留海森矩阵提供的高质量方向信息同时利用线搜索保证每一步都实际下降。阻尼牛顿法在代码上的改动很小主框架不变只是每次迭代里多一个回溯循环。但它能解决牛顿法一大半的不稳定问题。我建议在数模论文里如果涉及无约束优化算法实现直接写阻尼牛顿法比纯牛顿法体面得多。5.2 高斯-牛顿和 LM最小二乘问题的专用加速数模里大量问题最终都归结为最小二乘比如曲线拟合、参数估计。假设目标函数是$$f(x) \frac{1}{2} \sum_{i1}^{m} r_i(x)^2$$其中 r_i 是第 i 个残差。直接算海森矩阵太麻烦但有个漂亮的近似海森矩阵可以写成 J^T J 残差的相关项其中 J 是残差向量对 x 的雅可比矩阵。当残差比较小、模型接近线性时可以用 J^T J 近似海森矩阵于是迭代公式变成$$(J_k^{\top} J_k) d_k -J_k^{\top} r_k$$这就是高斯-牛顿法。它在工程上非常常用因为只需要一阶导数信息却能有接近牛顿法的收敛速度。如果再加一个阻尼因子 \lambda让迭代公式变成 (J^T J \lambda I) d -J^T r就是 LM 算法。\lambda 大时算法接近梯度下降variance 小步幅可控\lambda 小时算法接近高斯-牛顿收敛快。我遇到拟合问题首选就是 LM。5.3 拟牛顿法不想算海森就用逐步逼近还有一种情况目标函数的解析表达式特别复杂手动推海森矩阵不现实数值差分计算海森代价也高。这时候可以用拟牛顿法。它的思路是不直接计算海森矩阵而是维护一个近似矩阵 B_k每次迭代后利用梯度变化信息做低秩更新让 B_k 逐渐逼近真实的海森矩阵。最著名的家族是 DFP 和 BFGS其中 BFGS 在数值稳定性上通常优于 DFP。拟牛顿法的收敛速度介于梯度下降和牛顿法之间是超线性收敛工程上已经足够快。Python 里 scipy.optimize.minimize 的methodBFGS就是这个算法你只需要提供目标函数和梯度连梯度都不用提供可以用数值微分。我在数模里如果不想手推海森就会把 BFGS 当备选项。5.4 一张表帮你选型很多朋友看到多种优化方法后会纠结到底用哪个。我把常用方法的特性整理成一张表你在数模赛场上直接按表取用方法单次迭代成本收敛速度对初值要求适用场景梯度下降低线性宽松大规模问题粗略优化纯牛顿法高需要海森矩阵二次苛刻低维、光滑强凸、海森易求阻尼牛顿法高二次局部中等数模中最推荐手写的牛顿法变体高斯-牛顿法中接近二次中等最小二乘拟合与参数估计LM 算法中接近二次中等带非线性的最小二乘、曲线拟合BFGS低超线性中等海森解析难求时注意一点表格说“对初值要求”是相对而言的。真实比赛里条件不好的函数再好的方法也会因为初值不好而失败所以多起点启动是通用的增强手段。6. 数模实战里的几条个人经验最后聊几句真实的比赛经验不写虚的。第一动笔写代码之前先明确你的目标函数是什么形式。如果是最小二乘结构别死磕手写纯牛顿直接用高斯-牛顿或者 LM省掉海森矩阵推导一整个环节而且更稳。如果是小型低维的普通无约束函数手写阻尼牛顿法也不过四五十行代码完全可控。第二写牛顿类算法的时候一定先用小的二维测试函数验证梯度和海森的解析表达式。我吃过一次亏在三维问题上某个混合偏导少乘了一项结果算法表现比梯度下降还差查了半天才发现是符号推错。后来我无论多简单都先跑一次数值差分校验把校验函数放到工具模块里一劳永逸。第三警惕多峰函数的全局最优陷阱。牛顿法是局部优化算法遇到多峰函数初值落在哪个峰附近它就会收敛到哪个峰。稳妥的做法是先画一个粗略的网格图或者用随机采样看目标函数的形态判断它是不是多峰。如果多峰就采样几十个初始点分别跑最后对结果做聚类而不是直接相信某一次的收敛值。第四论文写作上有一点小技巧算法推导部分写清楚泰勒展开到二次项、令梯度为零、得到迭代公式这个逻辑链条就够了不用贴完整代码。然后放一张迭代收敛表列“迭代次数、目标函数值、梯度范数”三列最好再附上 MATLAB 或 Python 画的收敛曲线图。评审老师看多的都是梯度下降你拿出一张几步收敛到 1e-8 的表说服力一下就出来了。我个人在实际操作中的体会是多变量最优化计算里牛顿法的核心不是那个公式而是它背后“用二次模型代替真实函数”的思维方式。当你理解了这一点再去看阻尼牛顿、高斯-牛顿、LM、拟牛顿这些变体会发现它们全是围绕“海森矩阵昂贵且不稳定”这一痛点做的妥协。把这个主线抓住以后遇到任何优化问题你都不会慌——先判断结构再选方法最后用测试函数验证。这条路我走过确实可靠。 
长河君 · 2026/9/16 2:40:12 · 阅读 1.2k 多变量牛顿法全解析：从公式推导到Python实现与改进
如果你在数学建模竞赛里解过带非线性目标函数的优化问题大概率经历过两种痛苦用梯度下降吧步长调不好收敛慢得让人怀疑人生用现成的优化库吧又不知道它内部到底在做什么论文里写算法部分心里发虚。我第一次在国赛里硬啃一个参数估计问题目标函数是非线性最小二乘的形式梯度下降跑了三百多步还在原地打转后来老老实实把多变量最优化计算里的牛顿法补上才觉得这块真正通了一点。这篇笔记不打算写成教科书式的推导罗列就按我自己从单变量到多变量、从理论到代码、从踩坑到改进的真实过程来写涉及的核心关键词是所有优化方法里绕不开的牛顿法希望对准备数模或者做工程优化的朋友有直接的参考价值。1. 我为什么在数模里最终选了牛顿法——它到底在算什么1.1 一个容易被忽略的前提牛顿法解的其实是个方程先说个特别容易混淆的点。很多人一看到单变量牛顿法第一反应是高中课本里那个“用切线逐步逼近零点”的公式$$x_{k1} x_k - \frac{f(x_k)}{f(x_k)}$$这个公式是拿来解方程 f(x)0 的不是直接求极值。那多变量最优化里怎么用它关键在于求一个光滑函数的极小值点本质上是求梯度为零的点。也就是说我们要解的是 ∇f(x)0 这个方程组而不是直接在 f(x) 上做牛顿迭代。这个转换看起来只是换了个函数实际上影响巨大。你要处理的方程个数和变量个数一样多迭代公式里的分母也变成了一整个矩阵——二阶偏导组成的海森矩阵 H。之前很多教程把这步跳过了直接丢出多维迭代公式读者看得一头雾水其实原因就在这里牛顿法求极值逻辑内核是“求梯度的零点”。1.2 从切线到曲面单变量的直觉怎么推广到多维单变量的几何直觉是在当前位置放一条切线让切线延伸到和 x 轴相交把交点作为下一步。多变量情况下目标函数是一个 n 维曲面你不可能放一条线你应该放一个弯曲的曲面。什么曲面最合适二次曲面因为任何一个光滑函数在足够小的邻域里都能用二次函数去逼近。所以多变量牛顿法干的事情是每一步都用原函数在这个点的二次近似模型来代替原函数然后直接跳到这个二次模型的极值点。如果原函数本身就是二次函数一步就到如果不太像二次函数那就多迭代几次每一次都在新的点上重新拟合一个二次曲面。这比梯度下降聪明的地方在于梯度下降只用了当前位置的坡度信息一阶导数它知道哪个方向是下坡但不知道坡面有多弯曲所以只能拿一个固定的学习率小心翼翼地走。牛顿法用了坡度加曲率信息二阶导数它不仅能判断方向还能根据曲率预估出“再往前走多远就应该到谷底”于是它可以大胆地一次性跨出一大步。1.3 二次收敛是什么概念我当年第一次看到“牛顿法具有二次收敛性”这个表述时毫无感觉直到自己在 R 语言里跑了一组误差数据才真正被震撼。把误差记为 e_k二次收敛意味着误差按这样的规律缩小$$e_{k1} \approx C e_k^2$$假设 C 不是特别离谱误差序列可能是 0.1 → 0.001 → 10^(-8) → 10^(-16)。每一步的误差大约是上一步的平方。也就是说一旦进入牛顿法的收敛邻域精度提升是指数级的三到五步就能把误差压到双精度浮点的极限。对比梯度下降的线性收敛e_{k1}≈C e_k每一步误差只缩一个固定比例要磨很久才能到高精度。所以如果你在数模论文里需要展示算法效率把牛顿法和梯度下降的迭代次数放在同一张表里效果非常直观评审老师看一眼就明白你选的为什么是牛顿法。2. 多变量牛顿法的数学推导每一步都不是魔法2.1 从泰勒展开到驻点方程假设目标函数 f: R^n → R在当前迭代点 x_k 附近足够光滑。我们把它做二阶泰勒展开$$f(x) \approx f(x_k) g_k^{\top}(x - x_k) \frac{1}{2}(x - x_k)^{\top} H_k (x - x_k)$$其中g_k ∇f(x_k) 是梯度向量n 维H_k ∇²f(x_k) 是海森矩阵n×n 对称矩阵第 i 行第 j 列的元素就是 ∂²f/∂x_i ∂x_j 在 x_k 处的值。现在对这个二次模型求梯度并让它等于零。对向量 x 求梯度得到的是一个线性方程组$$g_k H_k (x - x_k) 0$$如果 H_k 是非奇异的就可以解出下一步的迭代点$$x_{k1} x_k - H_k^{-1} g_k$$整个多变量最优化计算里最重要的公式其实就是这一条。你看它没有任何魔法就是从“梯度为零”这个极值必要条件出发加上二阶近似推导出来的结果。2.2 为什么迭代方向是 H^{-1}g 而不是 g这里有个常见的疑问为什么牛顿方向是 H_k^{-1} g_k而不是直接用梯度方向直观理解是梯度方向只告诉你“往下坡走”但没有考虑每个方向上的曲率。海森矩阵编码的是局部曲率信息用它的逆矩阵去左乘梯度相当于对坐标做了一次与曲率成反比的伸缩变换。我经常用一个类比你在一个狭长的山谷里谷底的走向是斜的。梯度方向可能指向谷壁如果在梯度方向上走大步会撞到山坡而海森矩阵知道了“垂直谷底方向曲率很大”它会把这一步幅度压小把主要步长放到“沿着谷底那个曲率小的方向”上去。这就是为什么二次函数上一牛顿法一步到位——因为它把曲率信息全部利用了相当于直接把谷底抛物线方程解出来了。2.3 手算一个二维例子的海森矩阵以最经典的优化测试函数 Rosenbrock 函数为例$$f(x_1, x_2) (1-x_1)^2 100(x_2 - x_1^2)^2$$它的梯度是$$\frac{\partial f}{\partial x_1} -2(1-x_1) - 400x_1(x_2 - x_1^2)$$$$\frac{\partial f}{\partial x_2} 200(x_2 - x_1^2)$$海森矩阵的四个分量$$\frac{\partial^2 f}{\partial x_1^2} 2 - 400(x_2 - x_1^2) 800x_1^2$$$$\frac{\partial^2 f}{\partial x_1 \partial x_2} -400x_1$$$$\frac{\partial^2 f}{\partial x_2 \partial x_1} -400x_1$$$$\frac{\partial^2 f}{\partial x_2^2} 200$$注意海森矩阵是对称的两个混合偏导相等。这是检验你手算是否出错的一个天然手段。我建议在纸面推导时每次都检查 H_12 和 H_21 是否一致如果不等八成是某个链式法则求错了。3. 手写一个能用的多变量牛顿法代码与循环条件3.1 迭代里用“解方程组”而不是“求逆”理论上迭代公式是 x_{k1} x_k - H_k^{-1} g_k但实际代码里千万不要直接去计算矩阵的逆。原因有两个第一求逆的计算量是 O(n^3)而且常数很大第二数值上求逆再乘向量误差比直接解线性方程组大得多。正确的做法是引入方向向量 d_k解这个线性方程组$$H_k d_k -g_k$$然后令 x_{k1} x_k d_k。在 Python 里调用 numpy 的np.linalg.solve一行就能搞定它内部用 LU 分解速度快而且数值稳定性更好。这个细节在数模论文里的算法描述环节值得写进去显得你真正实现过而不是抄了一遍公式。3.2 一个完整的 Python 实现我用 Rosenbrock 函数做测试手写一个纯正的多变量牛顿法。这个函数的最优解是 (1, 1)最优值是 0公认的优化测试函数用来验证算法很合适。import numpy as np def rosen(x): x1, x2 x return (1 - x1)**2 100.0 * (x2 - x1**2)**2 def rosen_grad(x): x1, x2 x g np.zeros(2) g[0] -2.0 * (1 - x1) - 400.0 * x1 * (x2 - x1**2) g[1] 200.0 * (x2 - x1**2) return g def rosen_hess(x): x1, x2 x h11 2.0 - 400.0 * (x2 - x1**2) 800.0 * x1**2 h12 -400.0 * x1 h22 200.0 return np.array([[h11, h12], [h12, h22]]) def newton_minimize(grad, hess, x0, tol1e-8, max_iter100): x np.array(x0, dtypefloat) for i in range(max_iter): g grad(x) H hess(x) grad_norm np.linalg.norm(g, ordnp.inf) print(fiter {i:3d}: x ({x[0]:.10f}, {x[1]:.10f}), ||g||inf {grad_norm:.3e}) if grad_norm tol: return x, i, grad_norm # 关键步骤解 H d -g而不是求逆 d np.linalg.solve(H, -g) x x d return x, max_iter, np.linalg.norm(grad(x), ordnp.inf) x_opt, it, gn newton_minimize(rosen_grad, rosen_hess, [-1.2, 1.0]) print(f\n最优解: {x_opt}) print(f迭代次数: {it}) print(f最终梯度无穷范数: {gn:.3e})我第一次跑这段代码的时候从初始点 (-1.2, 1.0) 出发大概只迭代了五六步就到了 (1, 1) 附近。这就是二次收敛的威力。如果用梯度下降同样的初始点、同样的精度通常需要上千步。3.3 停机条件怎么设才不翻车写优化算法停机条件是特别容易被忽略但特别重要的部分。我自己常用的有三条优先级从高到低梯度范数足够小||g||_inf tol。这是最优性的一阶必要条件无约束光滑问题里它趋近于零说明已经接近驻点。tol 通常取 1e-6 到 1e-8具体看你的目标函数数值尺度。如果函数值本身很大梯度范数阈值也相应放大一点。迭代点变化量足够小||x_{k1} - x_k|| tol。这个条件能防止梯度范数卡在一个不小的值上但迭代点已经动不了的情况——虽然这不常见但碰到海森矩阵接近奇异时会遇到。最大迭代次数兜底max_iter。纯牛顿法有时候会震荡甚至发散必须设一个上限防止死循环。三条都用不要省。尤其第二条我见过不少人在梯度过零点附近来回跳动但梯度范数始终降不到 1e-8 以下的情况没有最大迭代数限制程序就卡死了。4. 牛顿法翻车实录海森矩阵奇异、不下降、初值敏感4.1 最直接的崩溃海森矩阵奇异解不出来多变量牛顿法看起来很美但现实里第一个坑就是海森矩阵可能奇异。什么叫奇异就是行列式为零方程组 H d -g 没有唯一解。几何上对应着当前点所在的局部二次模型是个“平的”或“有个方向没有曲率”的曲面。这种情况在球面型函数或者一些对称性较强的函数里很容易发生。比如 f(x)x1^4 x2^4在某些初始点附近海森矩阵可能是奇异的。代码里表现为np.linalg.solve抛出LinAlgError。如果你只是想要一个能跑的成绩最简单的方式是给海森矩阵加上一个小的单位矩阵的倍数把迭代公式改成$$(H_k \lambda I) d_k -g_k$$\lambda 取一个很小的正数比如 1e-6。这个操作本质上是保证矩阵正定也是后面要讲的阻尼牛顿法和 LM 算法的核心思想。4.2 牛顿方向不一定下降这个坑藏得更深。牛顿法推导时默认了海森矩阵是正定的但实际函数在海森正定以外的区域牛顿方向 d -H^{-1} g 可能与真正的梯度下降方向夹角超过 90 度也就是说沿着这个方向走目标函数不仅不下降还会上升。原因在于海森矩阵非正定时局部二次模型不是“碗型”而是“马鞍形”二次模型的极值点可能是一个鞍点甚至极大值点你跳过去自然不下降。判断方法很简单每次算出 d 之后检查一下 g^T d 是否小于零。如果 g^T d 0说明 d 不是下降方向这时候就不要固执地用牛顿方向了否则数值会瞬间爆炸。4.3 初值敏感一步跳飞牛顿法在局部收敛邻域内是神出了邻域直接变鬼。Rosenbrock 函数还算友好从 (-1.2, 1.0) 出发没问题。但如果你从更远的地方出发比如 (10, 10)二次模型完全看不出真实谷底在哪一步迭代可能直接跳到函数值大到不可描述的区域。我在自己的试验里见过最典型的情况初始点选得不合适第一步 d 的长度超过 1e4第二步目标函数值变成 10^200程序当场溢出。解决思路有两个一是用多点随机初始化比如均匀撒 50 个初始点每个点跑一遍牛顿法收敛后再看极值分布二是先用梯度下降或单纯形法粗跑几十步得到一个较好的初值再切换到牛顿法收尾。这种“粗糙全局精细局部”的组合拳在数模实战里非常有效。4.4 计算成本与数值误差多变量牛顿法每次迭代都要计算海森矩阵并求解一个 n 阶线性方程组复杂度 O(n^3)。当变量个数 n 到几千上万时这个代价大到基本不可用。所以牛顿法在小规模问题上表现优秀大规模问题上基本没有优势这个一定要在选型时想清楚。数值误差方面海森矩阵如果是手动推导的非常容易出错。我的经验是无论用链式法则手推时觉得多简单都要利用“海森矩阵必须对称”这条性质检查一遍。另外可以用数值差分法验证梯度和海森的正确性。所谓数值差分就是让某个变量动一个小量 ε用 (f(xεe_i) - f(x-εe_i)) / (2ε) 近似偏导数然后和解析表达式对比。这个方法虽然跑起来慢但只用来测试正确性一锤定音。5. 从牛顿法出发的改进路线阻尼、高斯-牛顿、拟牛顿怎么选5.1 阻尼牛顿法一行代码救回大半场景理解了“牛顿方向不下降”的问题之后最自然的修复方案就是给牛顿方向加一个步长因子变成$$x_{k1} x_k \alpha_k d_k$$其中 \alpha_k 通过一维线搜索确定比如 Armijo 准则从 \alpha1 开始不断减半直到目标函数充分下降。这个方法叫阻尼牛顿法本质是保留海森矩阵提供的高质量方向信息同时利用线搜索保证每一步都实际下降。阻尼牛顿法在代码上的改动很小主框架不变只是每次迭代里多一个回溯循环。但它能解决牛顿法一大半的不稳定问题。我建议在数模论文里如果涉及无约束优化算法实现直接写阻尼牛顿法比纯牛顿法体面得多。5.2 高斯-牛顿和 LM最小二乘问题的专用加速数模里大量问题最终都归结为最小二乘比如曲线拟合、参数估计。假设目标函数是$$f(x) \frac{1}{2} \sum_{i1}^{m} r_i(x)^2$$其中 r_i 是第 i 个残差。直接算海森矩阵太麻烦但有个漂亮的近似海森矩阵可以写成 J^T J 残差的相关项其中 J 是残差向量对 x 的雅可比矩阵。当残差比较小、模型接近线性时可以用 J^T J 近似海森矩阵于是迭代公式变成$$(J_k^{\top} J_k) d_k -J_k^{\top} r_k$$这就是高斯-牛顿法。它在工程上非常常用因为只需要一阶导数信息却能有接近牛顿法的收敛速度。如果再加一个阻尼因子 \lambda让迭代公式变成 (J^T J \lambda I) d -J^T r就是 LM 算法。\lambda 大时算法接近梯度下降variance 小步幅可控\lambda 小时算法接近高斯-牛顿收敛快。我遇到拟合问题首选就是 LM。5.3 拟牛顿法不想算海森就用逐步逼近还有一种情况目标函数的解析表达式特别复杂手动推海森矩阵不现实数值差分计算海森代价也高。这时候可以用拟牛顿法。它的思路是不直接计算海森矩阵而是维护一个近似矩阵 B_k每次迭代后利用梯度变化信息做低秩更新让 B_k 逐渐逼近真实的海森矩阵。最著名的家族是 DFP 和 BFGS其中 BFGS 在数值稳定性上通常优于 DFP。拟牛顿法的收敛速度介于梯度下降和牛顿法之间是超线性收敛工程上已经足够快。Python 里 scipy.optimize.minimize 的methodBFGS就是这个算法你只需要提供目标函数和梯度连梯度都不用提供可以用数值微分。我在数模里如果不想手推海森就会把 BFGS 当备选项。5.4 一张表帮你选型很多朋友看到多种优化方法后会纠结到底用哪个。我把常用方法的特性整理成一张表你在数模赛场上直接按表取用方法单次迭代成本收敛速度对初值要求适用场景梯度下降低线性宽松大规模问题粗略优化纯牛顿法高需要海森矩阵二次苛刻低维、光滑强凸、海森易求阻尼牛顿法高二次局部中等数模中最推荐手写的牛顿法变体高斯-牛顿法中接近二次中等最小二乘拟合与参数估计LM 算法中接近二次中等带非线性的最小二乘、曲线拟合BFGS低超线性中等海森解析难求时注意一点表格说“对初值要求”是相对而言的。真实比赛里条件不好的函数再好的方法也会因为初值不好而失败所以多起点启动是通用的增强手段。6. 数模实战里的几条个人经验最后聊几句真实的比赛经验不写虚的。第一动笔写代码之前先明确你的目标函数是什么形式。如果是最小二乘结构别死磕手写纯牛顿直接用高斯-牛顿或者 LM省掉海森矩阵推导一整个环节而且更稳。如果是小型低维的普通无约束函数手写阻尼牛顿法也不过四五十行代码完全可控。第二写牛顿类算法的时候一定先用小的二维测试函数验证梯度和海森的解析表达式。我吃过一次亏在三维问题上某个混合偏导少乘了一项结果算法表现比梯度下降还差查了半天才发现是符号推错。后来我无论多简单都先跑一次数值差分校验把校验函数放到工具模块里一劳永逸。第三警惕多峰函数的全局最优陷阱。牛顿法是局部优化算法遇到多峰函数初值落在哪个峰附近它就会收敛到哪个峰。稳妥的做法是先画一个粗略的网格图或者用随机采样看目标函数的形态判断它是不是多峰。如果多峰就采样几十个初始点分别跑最后对结果做聚类而不是直接相信某一次的收敛值。第四论文写作上有一点小技巧算法推导部分写清楚泰勒展开到二次项、令梯度为零、得到迭代公式这个逻辑链条就够了不用贴完整代码。然后放一张迭代收敛表列“迭代次数、目标函数值、梯度范数”三列最好再附上 MATLAB 或 Python 画的收敛曲线图。评审老师看多的都是梯度下降你拿出一张几步收敛到 1e-8 的表说服力一下就出来了。我个人在实际操作中的体会是多变量最优化计算里牛顿法的核心不是那个公式而是它背后“用二次模型代替真实函数”的思维方式。当你理解了这一点再去看阻尼牛顿、高斯-牛顿、LM、拟牛顿这些变体会发现它们全是围绕“海森矩阵昂贵且不稳定”这一痛点做的妥协。把这个主线抓住以后遇到任何优化问题你都不会慌——先判断结构再选方法最后用测试函数验证。这条路我走过确实可靠。