跳转到内容
← 返回算法
数值算法计算机科学 · 数值计算 · 优化约 17 分钟阅读

牛顿法与数值求解

Newton's Method and Numerical Analysis

1669 年，艾萨克·牛顿（Isaac Newton）在一份未发表的手稿中描述了一种求多项式根的迭代方法。三个多世纪后，这个方法仍然是科学计算的基石：每一次 GPS 定位、每一次卫星轨道修正、每一次机器学习模型的二阶优化，背后都有它的影子。 牛顿法在教科书里通常以"求方程 $f(x) = 0$ 的根"的形式出现，但这只…

牛顿法数值方法求根算法优化迭代

1669 年，艾萨克·牛顿（Isaac Newton）在一份未发表的手稿中描述了一种求多项式根的迭代方法。三个多世纪后，这个方法仍然是科学计算的基石：每一次 GPS 定位、每一次卫星轨道修正、每一次机器学习模型的二阶优化，背后都有它的影子。

破除误解：牛顿法不仅仅是"求根"

牛顿法在教科书里通常以"求方程 $f(x) = 0$ 的根"的形式出现，但这只是它的一个应用。更广义地看，牛顿法是利用函数的局部线性近似（切线）反复逼近目标点的框架。

它的三个主要用途： 1. 求根：找 $f(x) = 0$ 的解 2. 求极值：对 $f'(x) = 0$ 应用牛顿法，即求导函数的根 3. 求逆：如求平方根 a\sqrt{a}，等价于解 x2−a=0x^2 - a = 0

核心思想：用切线代替曲线

设我们要求 $f(x) = 0$，当前猜测值为 xnx_n。在 xnx_n 处作 $f$ 的切线：

y=f(xn)+f′(xn)(x−xn)y = f(x_n) + f'(x_n)(x - x_n)

令切线等于零，解得切线与 $x$ 轴的交点：

xn+1=xn−f(xn)f′(xn)x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}

这就是牛顿迭代公式。每次迭代用当前点的切线（一阶近似）替代原函数，找到新的猜测值，直到 ∣f(xn)∣|f(x_n)| 足够小。

收敛速度：二次收敛的奇迹

牛顿法最令人惊叹的性质是它的二次收敛（Quadratic Convergence）：在根的附近，每次迭代误差大约被平方。若当前误差 en=xn−x∗e_n = x_n - x^*，则：

en+1≈−f′′(x∗)2f′(x∗)en2e_{n+1} \approx -\frac{f''(x^*)}{2f'(x^*)} e_n^2

这意味着：从 4 位有效数字开始，下一步得到约 8 位，再下一步约 16 位。实践中，牛顿法求解到机器精度通常只需要 5–10 次迭代。对比梯度下降的线性收敛（每步误差乘以固定常数 $< 1$），二次收敛是质的飞跃。

经典应用：平方根的快速计算

计算 a\sqrt{a}，令 f(x)=x2−af(x) = x^2 - a，则 $f'(x) = 2x$，迭代公式变为：

xn+1=12(xn+axn)x_{n+1} = \frac{1}{2}\left(x_n + \frac{a}{x_n}\right)

以 $a = 2$，x0=1x_0 = 1 为例：

迭代次数xnx_n误差
01.000000−0.414213…
11.5000000.085786…
21.4166670.002454…
31.4142160.0000021…
41.414213562…<10−12< 10^{-12}

4 步达到 12 位精度——这就是二次收敛的威力。这个迭代式（xn+1=12(xn+a/xn)x_{n+1} = \frac{1}{2}(x_n + a/x_n)）比牛顿本人早了一千六百多年：它就是古希腊的海伦迭代法，巴比伦泥板上也出现过类似算法。牛顿的贡献不是这个特例，而是指出任何可微函数都能这样干。

现场：一次迭代就够的时候

1999 年发布的《雷神之锤 III 竞技场》（Quake III Arena，源码于 2005 年由 id Software 以 GPL 开源）里有一段后来极为出名的代码，用来算 1/x1/\sqrt{x}——3D 渲染每帧要对成千上万个法向量做归一化，这个运算是热点中的热点：

c
float Q_rsqrt(float number) {
    long i; float x2, y;
    x2 = number * 0.5F;
    y  = number;
    i  = *(long*)&y;              // 把浮点数的二进制位当整数读
    i  = 0x5f3759df - (i >> 1);   // 神秘常数：直接算出一个初始猜测
    y  = *(float*)&i;
    y  = y * (1.5F - (x2 * y * y)); // ← 一次牛顿迭代，就一次
    return y;
}
```

两步：先用位运算+魔数 0x5f3759df 造一个粗糙初值（利用 IEEE 754 里指数位近似等于 log⁡2x\log_2 x 这个事实），再做一次牛顿迭代收尾。对 f(y)=1/y2−xf(y) = 1/y^2 - x 应用牛顿法，迭代式化简后正是 y * (1.5 - x/2 * y * y)——代码最后一行的来源。

一次迭代就把最大相对误差压到约 1.75×10−31.75 \times 10^{-3}，对游戏画面完全够用。这是二次收敛的另一种用法：初值给得好，一步就到位。（后来的分析指出最优魔数是 0x5F375A86，误差略低一点点；这个位操作技巧的源头可追到 Kahan 与 Ng 1986 年的一份未发表笔记，经 Greg Walsh 之手进入游戏代码。）

一个失败现场：它也会原地打转

牛顿法的教科书版本给人一种"必然收敛"的错觉。看这个例子：

f(x)=x3−2x+2,f′(x)=3x2−2f(x) = x^3 - 2x + 2, \qquad f'(x) = 3x^2 - 2

从 x0=0x_0 = 0 出发：

  • $f(0) = 2$，$f'(0) = -2$ → x1=0−2−2=1x_1 = 0 - \frac{2}{-2} = 1
  • $f(1) = 1$，$f'(1) = 1$ → x2=1−11=0x_2 = 1 - \frac{1}{1} = 0

0→1→0→1→⋯0 \to 1 \to 0 \to 1 \to \cdots，无限二循环，永远到不了那个真实的根 x∗≈−1.7693x^* \approx -1.7693。它不是收敛慢，是根本不动。

再看文献里常提的 f(x)=x1/3f(x) = x^{1/3}（根为 0）。代入迭代式：

xn+1=xn−xn1/313xn−2/3=xn−3xn=−2xnx_{n+1} = x_n - \frac{x_n^{1/3}}{\frac13 x_n^{-2/3}} = x_n - 3x_n = -2x_n

每一步都跳到反向、并且距离翻倍——从任何非零点出发都必然发散。原因是 $f'(0)$ 不存在，二次收敛的前提（根附近 $f'$ 有界且非零、$f''$ 有界）被打破了。

结论：二次收敛是一条局部定理。牛顿法的实际工程难点从来不是那个漂亮的迭代式，而是"怎么保证进入收敛域"。

多维推广：牛顿-拉弗森法

在高维空间中，对向量函数 F(x)=0\mathbf{F}(\mathbf{x}) = \mathbf{0}，牛顿法推广为：

xn+1=xn−J(xn)−1F(xn)\mathbf{x}_{n+1} = \mathbf{x}_n - J(\mathbf{x}_n)^{-1} \mathbf{F}(\mathbf{x}_n)

其中 $J$ 是雅可比矩阵（Jacobian Matrix），每个元素 Jij=∂Fi/∂xjJ_{ij} = \partial F_i / \partial x_j。实践中不直接求逆，而是解线性方程组 JΔx=−FJ \Delta x = -F，再令 xn+1=xn+Δxx_{n+1} = x_n + \Delta x。

优化中的牛顿法

在优化（求 $f$ 的极小值）中，目标是令梯度 ∇f=0\nabla f = 0，对梯度应用牛顿法：

xn+1=xn−H(xn)−1∇f(xn)\mathbf{x}_{n+1} = \mathbf{x}_n - H(\mathbf{x}_n)^{-1} \nabla f(\mathbf{x}_n)

其中 $H$ 是海森矩阵（Hessian Matrix）——$f$ 的二阶偏导矩阵。这比梯度下降收敛快得多，但代价是每步需要计算和求逆 n×nn \times n 的 Hessian 矩阵，对大模型（n∼109n \sim 10^9）完全不可行。

方法每步计算代价收敛速度
梯度下降$O(n)$线性
牛顿法O(n3)O(n^3)（Hessian 逆）二次
拟牛顿法（BFGS）O(n2)O(n^2)超线性

拟牛顿法（Quasi-Newton）：用梯度信息逐步逼近 Hessian 的逆，避免直接计算二阶导数。L-BFGS（Limited-memory BFGS）是科学计算和中小规模机器学习中最常用的优化算法之一。

实践中的牛顿法：没人用教科书那一版

上面那个"原地打转"的例子说明，工业级求解器不可能裸用 xn+1=xn−f/f′x_{n+1} = x_n - f/f'。真实代码里通常裹着三层护具：

1. 阻尼与线搜索（Damped Newton / Line Search）。牛顿步 Δx=−H−1∇f\Delta x = -H^{-1}\nabla f 只决定方向，步长交给一维搜索：先试 α=1\alpha = 1，若目标函数下降得不够多就把 α\alpha 减半重试。判据通常是 Armijo 条件

f(x+αΔx)≤f(x)+c1α∇f⊤Δx,c1≈10−4f(x + \alpha \Delta x) \le f(x) + c_1 \alpha \nabla f^\top \Delta x, \quad c_1 \approx 10^{-4}

即"实际下降量至少达到线性预测下降量的 c1c_1 倍"。这一条就足以排除掉"跳得太远反而变差"的情形，把全局收敛性救回来，代价是靠近解时退化成几次多余的函数求值。

2. 信赖域（Trust Region）。换个思路：不问"沿这个方向走多远"，而是先划一个半径 Δ\Delta 的球，只在球内相信二次模型，在球内求最优步。若这一步实际效果与模型预测吻合就扩大 Δ\Delta，否则缩小。它比线搜索更能处理 Hessian 不正定的情形——牛顿方向此时可能指向鞍点甚至上坡，而信赖域子问题依然有解。

3. 正则化：Levenberg–Marquardt。求解非线性最小二乘 min⁡∑ri(x)2\min \sum r_i(x)^2 时，用 J⊤J+λIJ^\top J + \lambda I 替代 Hessian：λ→0\lambda \to 0 时退化为高斯-牛顿（快），λ\lambda 很大时退化为小步长梯度下降（稳）。算法在迭代中自动调 λ\lambda，在"快"和"稳"之间连续滑动。这是曲线拟合、相机标定、SLAM 后端优化里的默认工具。

二阶方法在深度学习里的回归

前面那张表给了一个悲观结论：Hessian 求逆 O(n3)O(n^3)，参数量 10910^9 的模型没戏。所以十年来深度学习几乎全用一阶方法（SGD、Adam）。但故事在改写——关键在于不求完整 Hessian，只求它的结构化近似：

  • K-FAC（Martens & Grosse, 2015）：把 Fisher 信息矩阵近似成层内两个小矩阵的 Kronecker 积，于是"求逆一个百万维矩阵"变成"求逆两个几百维矩阵"。
  • Shampoo（Gupta, Koren & Singer, 2018）：对每个参数张量的每个维度各维护一个预条件矩阵，同样用 Kronecker 结构把代价压下来。注意它预条件的是梯度二阶矩而非严格 Hessian，属于"非对角预条件"而非纯粹的牛顿法。

2024 年 MLCommons 举办的首届 AlgoPerf 训练算法竞赛给出了一个有分量的实证结果：Meta 提交的 Distributed Shampoo 拿下外部调参赛道第一，在跨多个真实工作负载的墙钟时间上比强基线快约 28%。这一点很关键——以往二阶方法总能宣称"每步进展更大"，但算上每步的额外开销就不划算了；AlgoPerf 是在同一块墙钟表上比的。

牛顿的那条切线，绕了 355 年重新回到了最前沿的优化器里。

代价与争议

收敛只在局部保证：见上文的两个失败现场。这不是可以靠"写好代码"消除的缺陷，而是方法的内在性质，只能靠阻尼、信赖域或更好的初值来管理。

需要导数：需要 $f'(x)$（或 Hessian）的解析表达式，某些实际问题中难以得到。割线法（Secant Method）用差商近似导数，收敛阶降到约 1.618（黄金比例）——比线性快、比二次慢。自动微分的普及大幅削弱了这条限制：今天求导常常是编译器的事，而不是人的事。

重根退化：若 f′(x∗)=0f'(x^*) = 0，即根是重根，牛顿法退化为线性收敛，效率大幅下降。修正方法是使用 xn+1=xn−m⋅f(xn)/f′(xn)x_{n+1} = x_n - m \cdot f(x_n)/f'(x_n)（$m$ 为根的重数）——但重数通常事先不知道，这让修正在实践中难用。

混沌与分形：把牛顿法搬到复平面上，行为会变得极其古怪。取最简单的 f(z)=z3−1f(z) = z^3 - 1，它有三个根（$1$ 和两个复根）。给平面上每个初始点按"最终收敛到哪个根"染色，得到的不是三块规整区域，而是三色互相缠绕、边界处处自相似的牛顿分形（Newton Fractal）——边界上任意一点的任意小邻域内，三种颜色都同时出现。也就是说，初值哪怕挪动 10−1510^{-15}，结果也可能换一个根。

这件事的分量超出数值分析：一个每步都由确定公式给出、连"随机"都没有的算法，它的全局行为却对初值敏感依赖到无法预测。这正是确定性 ≠ 可预测性的教科书式例证，和 chaos-theory 讲的是同一件事。

跨域连接

  • 混沌理论：把迭代搬到复平面上，按"最终收敛到哪个根"给初始点染色，得到的是三色互相缠绕、边界处处自相似的分形——边界上任一点的任意小邻域内三种颜色同时出现。一个每步都由确定公式给出的算法，全局行为却对初值敏感依赖到无法预测。确定性不等于可预测性，这是最短的例证。
  • 梯度下降与反向传播：一阶方法只知道往哪下坡，二阶还知道坡的弯曲程度，因而能自动定步长、在病态方向上不来回打转。代价是曲率信息的规模：完整的二阶矩阵在十亿参数上不可能求逆，所以现代做法一律是结构化近似——按层分块、用克罗内克积分解，把"求逆一个百万维矩阵"换成求逆两个几百维的。
  • 大地测量与大地水准面：卫星定位要解一组关于三维位置与接收机钟差的非线性方程，标准解法就是这套迭代。雅可比矩阵的条件数直接决定定位精度：卫星在天空中挤成一团时矩阵接近奇异，同样的测距误差被放大好几倍，这就是几何精度因子。误差来自几何而非算法，换求解器救不了。
  • 化学平衡：多组分体系要同时满足元素守恒与多个平衡常数，这是一组强非线性方程。麻烦在于各物种浓度可以相差十几个数量级，直接对浓度迭代会让雅可比病态，还容易迭出负浓度。实践中改对浓度的对数迭代，既天然保证为正，又把量级差异压平——变量替换在这里比换求解器有用得多。
  • 期权、衍生品与布莱克-斯科尔斯：定价公式能从波动率算出价格，却无法反解出波动率。隐含波动率是把市场价当作方程右端、对波动率现算出来的，而它对波动率的导数恰好有闭式解，正好当分母。深度虚值期权是这条流程的失败现场：那里导数接近零，除法放大误差，迭代不稳定，必须换成有区间保证的方法兜底。

参考文献

  • Newton, I. De analysi per aequationes numero terminorum infinitas. 1669（手稿，1711 年出版）。
  • Raphson, J. Analysis Aequationum Universalis. 1690.
  • Nocedal, J. & Wright, S. Numerical Optimization. 2nd ed. Springer, 2006.（第 3 章线搜索与 Armijo 条件、第 4 章信赖域、第 6 章拟牛顿法、第 10 章 Levenberg–Marquardt）
  • Martens, J. & Grosse, R. Optimizing Neural Networks with Kronecker-factored Approximate Curvature. ICML 2015. arXiv:1503.05671.
  • Gupta, V., Koren, T. & Singer, Y. Shampoo: Preconditioned Stochastic Tensor Optimization. ICML 2018. arXiv:1802.09568.
  • Kasimbeg, P. et al. Accelerating Neural Network Training: An Analysis of the AlgoPerf Competition. ICLR 2025. arXiv:2502.15015.（首届 AlgoPerf 竞赛结果与墙钟时间口径）

延伸阅读

  • Press, W. H. et al. Numerical Recipes: The Art of Scientific Computing. 3rd ed. Cambridge University Press, 2007.（第 9 章求根、第 10 章极小化）
  • Trefethen, L. N. & Bau, D. Numerical Linear Algebra. SIAM, 1997.（理解牛顿步为何最终归结为"解一个线性系统"）

正文中引用本条 · 2 · 其中 2 条来自其他领域

数学代数基本定理刘徽