1. 从"拟合不上"说起:非线性最小二乘到底在解决什么问题
做数据拟合、参数辨识、位姿估计或者曲线拟合的人,几乎都绕不开一个场景:手里有一堆观测数据,心里有一个带参数的模型,但参数怎么调都拟合不上。线性问题还能靠矩阵求逆一步到位,一旦模型对参数是非线性的,事情就变得棘手了。非线性最小二乘(Non-linear Least Squares,简称 NLLS)就是专门对付这类问题的工具族。
它的核心目标非常朴素:找到一组参数,让模型预测值和实际观测值之间的残差平方和最小。用数学语言写出来就是最小化 $S(\mathbf{x}) = \sum_{i=1}^{m} r_i(\mathbf{x})^2$,其中 $r_i$ 是第 $i$ 个残差,$\mathbf{x}$ 是待求参数向量。听起来简单,但难点在于 $r_i(\mathbf{x})$ 对 $\mathbf{x}$ 是非线性的,这意味着目标函数可能存在多个局部极小值,没有解析解,只能靠迭代逼近。
我最初接触这块内容是在做传感器标定的时候。当时用的是一个简单的多项式模型,本以为随便调调参数就能收敛,结果发现初值稍微偏一点,迭代就发散或者卡在明显不对的结果上。后来才明白,非线性最小二乘不是"调包就能用"的黑盒,方法选型、初值策略、收敛判据、数值稳定性,每一个环节都会直接影响最终结果。这也是为什么我把这份整理写下来——它适合正在做参数拟合、曲线拟合、SLAM后端优化、神经网络训练底层推导,或者任何需要"让模型贴合数据"的开发者。不管你是刚接触数值优化的学生,还是已经用过scipy.optimize.least_squares但说不清背后原理的工程师,这份整理都能帮你把方法脉络理清楚。
2. 高斯-牛顿法:最直觉的迭代思路,以及它为什么会崩
2.1 从线性最小二乘"借"来的迭代框架
高斯-牛顿法(Gauss-Newton,简称 GN)是理解所有非线性最小二乘方法的起点。它的思路非常直接:既然非线性难处理,那我就在当前参数点附近把残差函数线性化,用线性最小二乘的解作为下一次迭代的方向。
具体来说,在当前参数 $\mathbf{x}_k$ 处对残差做一阶泰勒展开:$r_i(\mathbf{x}_k + \Delta \mathbf{x}) \approx r_i(\mathbf{x}_k) + \mathbf{J}_i \Delta \mathbf{x}$,其中 $\mathbf{J}_i$ 是第 $i$ 个残差对参数的雅可比行向量。把所有残差堆起来得到雅可比矩阵 $\mathbf{J}$,那么最小化线性化后的残差平方和就变成一个标准线性最小二乘问题,其正规方程为:
$$\mathbf{J}^T \mathbf{J} \Delta \mathbf{x} = -\mathbf{J}^T \mathbf{r}$$
解出 $\Delta \mathbf{x}$ 后更新 $\mathbf{x}_{k+1} = \mathbf{x}_k + \Delta \mathbf{x}$,反复迭代直到收敛。这个式子是整个非线性最小二乘领域的"心脏",后面要讲的 LM 方法也是在这个基础上做修正。
我实测下来,GN 法在残差接近零、模型接近线性的区域收敛速度极快,通常是二次收敛——也就是说每次迭代有效数字大约翻倍。这也是它在 bundle adjustment 这类问题里被广泛使用的原因。
2.2 雅可比矩阵的推导与计算:手推还是自动微分
GN 法能不能用,很大程度上取决于你能不能拿到准确的雅可比矩阵。这里有三条路:
- 解析求导:手动推导每个残差对每个参数的偏导。精度最高、速度最快,但模型一复杂就容易推错,而且改模型就要重推。
- 数值差分:用有限差分近似偏导,比如 $\partial r_i / \partial x_j \approx (r_i(\mathbf{x} + h\mathbf{e}_j) - r_i(\mathbf{x})) / h$。实现简单,但步长 $h$ 的选取很讲究,太大截断误差大,太小舍入误差大。
- 自动微分:借助工具对表达式做符号级求导,兼顾精度和开发效率。
提示:数值差分的步长经验值一般在 $10^{-8}$ 到 $10^{-6}$ 之间(针对双精度),但更稳妥的做法是用中心差分 $(r_i(\mathbf{x}+h\mathbf{e}_j) - r_i(\mathbf{x}-h\mathbf{e}_j))/(2h)$,精度从一阶提升到二阶,代价是每次求导多算一次函数值。
我在做一个小型位姿估计 Demo 时,一开始图省事用前向差分,结果迭代后期残差下降变得很慢,换成中心差分后收敛明显顺畅。后来分析原因,是前向差分的截断误差在接近最优解时干扰了梯度方向。
2.3 GN 法的三个致命弱点
GN 法虽然优雅,但它的假设很脆弱,实际用起来经常出问题:
第一,雅可比必须满秩。如果 $\mathbf{J}^T\mathbf{J}$ 奇异或接近奇异,正规方程就解不出来或者解出来是个巨大的方向。这在参数冗余、观测不足、或者模型在某些点退化时非常常见。
第二,残差不能太大。GN 法本质上是忽略了泰勒展开的二阶项。当残差较大、模型非线性强时,这个忽略会带来严重误差,导致迭代方向不可靠,甚至发散。
第三,没有步长控制。GN 法直接接受 $\Delta \mathbf{x}$,如果这一步迈得太大,目标函数反而可能上升。它没有"往回退"的机制。
这三点合起来,就是为什么实际工程里很少直接用裸的 GN 法,而是用下面要讲的带阻尼的 LM 方法。
3. 列文伯格-马夸特法:给迭代装上"刹车"和"油门"
3.1 阻尼项的引入逻辑
列文伯格-马夸特法(Levenberg-Marquardt,简称 LM)是工程上最常用的非线性最小二乘算法,没有之一。它的核心改动是在正规方程里加了一个阻尼项:
$$(\mathbf{J}^T\mathbf{J} + \lambda \mathbf{I})\Delta \mathbf{x} = -\mathbf{J}^T\mathbf{r}$$
这里的 $\lambda$ 就是阻尼因子。它的作用可以这样理解:当 $\lambda$ 很小时,方程退化成 GN 法,收敛快;当 $\lambda$ 很大时,$\mathbf{J}^T\mathbf{J}$ 被 $\lambda\mathbf{I}$ 主导,$\Delta\mathbf{x} \approx -\frac{1}{\lambda}\mathbf{J}^T\mathbf{r}$,这实际上变成了最速下降法的一小步,虽然慢但保证目标函数下降。
所以 LM 法本质上是 GN 法和最速下降法之间的连续过渡。$\lambda$ 就像一个"信任度旋钮":模型可信时调小,让 GN 大步快跑;模型不可信时调大,退化成保守的小步下降。
3.2 阻尼因子的更新策略:为什么不能简单调参
$\lambda$ 的更新是 LM 法的灵魂。常见策略是这样的:
- 计算当前步 $\Delta\mathbf{x}$,尝试更新参数。
- 如果新参数下目标函数下降,说明这一步有效,接受更新,并减小 $\lambda$(比如除以 10),下次更激进。
- 如果目标函数上升,说明这一步迈错了,拒绝更新,并增大 $\lambda$(比如乘以 10),下次更保守。
这个"接受就减、拒绝就加"的反馈机制,让 LM 法在远离最优解时表现得像最速下降法一样稳健,接近最优解时又像 GN 法一样快速。
注意:$\lambda$ 的初始值选取很关键。太小容易第一步就发散,太大则前期收敛极慢。经验上可以从 $\lambda_0 = 10^{-3} \cdot \max(\text{diag}(\mathbf{J}^T\mathbf{J}))$ 开始,让阻尼和问题尺度匹配。
3.3 尺度不变性:LM 法的一个隐藏坑
标准 LM 法加的阻尼是 $\lambda\mathbf{I}$,这假设所有参数的尺度是一样的。但实际模型里,参数往往量纲不同——比如一个参数是角度(量级 1),另一个是距离(量级 1000)。这时候 $\lambda\mathbf{I}$ 对不同参数施加的阻尼效果完全不同,会导致某些参数几乎不动,另一些剧烈震荡。
解决办法是用 $\lambda \cdot \text{diag}(\mathbf{J}^T\mathbf{J})$ 代替 $\lambda\mathbf{I}$,也就是马夸特本人后来提出的改进版本。这样阻尼会按每个参数的曲率自动缩放,实现尺度不变性。我在做多参数标定时踩过这个坑:一个距离参数和一个角度参数混在一起,用 $\lambda\mathbf{I}$ 版本怎么都不收敛,换成 $\text{diag}$ 版本后立刻正常。
| 阻尼形式 | 表达式 | 适用场景 | 缺点 |
|---|---|---|---|
| 原始 LM | $\mathbf{J}^T\mathbf{J} + \lambda\mathbf{I}$ | 参数尺度相近 | 尺度敏感 |
| 马夸特改进 | $\mathbf{J}^T\mathbf{J} + \lambda,\text{diag}(\mathbf{J}^T\mathbf{J})$ | 参数量纲差异大 | 实现稍复杂 |
| 信赖域形式 | 约束 $|\mathbf{D}\Delta\mathbf{x}| \le \Delta$ | 理论更严谨 | 需要子问题求解器 |
4. 收敛判据与终止条件:什么时候该停下来
4.1 四类常用判据及其陷阱
迭代算法必须有个"停"的规则,否则要么浪费算力,要么永远不收敛。非线性最小二乘里常用的判据有四类:
- 梯度判据:$|\mathbf{J}^T\mathbf{r}|_\infty < \epsilon_1$。梯度接近零说明到了驻点,但驻点可能是局部极大或鞍点,不一定是极小。
- 步长判据:$|\Delta\mathbf{x}| < \epsilon_2(|\mathbf{x}| + \epsilon_2)$。这一步几乎没动,说明已经稳定。
- 目标函数判据:$|S(\mathbf{x}_{k+1}) - S(\mathbf{x}_k)| < \epsilon_3$。函数值不再明显下降。
- 相对误差判据:结合上面几种,用相对量而非绝对量。
单独用任何一个都容易误判。比如梯度判据在病态问题里可能很早触发,步长判据在参数尺度差异大时不可靠。实践中通常组合使用,任一满足就停,但要注意组合顺序。
4.2 为什么"收敛"不等于"找到正确答案"
这是新手最容易误解的地方。非线性最小二乘的收敛只意味着迭代稳定了,不代表你找到了全局最优。目标函数可能有多个局部极小,算法停在哪个取决于初值。
我做过一个实验:同一个模型,初值取两组不同的值,LM 法都"收敛"了,但残差平方和差了将近三倍。这说明两次都停在了不同的局部极小。所以判断结果好坏,不能只看"是否收敛",还要看残差是否合理、参数是否符合物理意义、多组初值是否收敛到同一结果。
提示:工程上稳妥的做法是跑多组随机初值,取残差最小的那组。如果多组初值收敛到不同结果,说明问题本身病态或模型有冗余,需要重新审视模型结构。
4.3 最大迭代次数:一个必须设置的保险
无论判据多完善,都要设一个最大迭代次数上限。原因很简单:有些问题会陷入"接受-拒绝"的循环,$\lambda$ 反复调整但始终找不到下降方向,迭代次数无限增长。设个上限(比如 100 到 200 次),超了就报警,让开发者介入检查模型或数据。
5. 初值、病态与数值稳定性:那些让算法"翻车"的细节
5.1 初值选取:非线性优化的第一道坎
非线性最小二乘对初值极其敏感。好的初值能让算法几步收敛,坏的初值直接发散。初值策略没有万能公式,但有几条实用原则:
- 利用物理先验:如果参数有明确物理意义,用测量值或理论估计值作为初值。
- 先粗后精:先用简单模型或网格搜索粗略定位,再用 NLLS 精修。
- 多起点:随机撒若干初值,取最优。
- 参数变换:把有约束的参数(如必须为正)通过变换变成无约束,比如用 $\log$ 变换处理正参数。
我在做一个衰减曲线拟合时,衰减系数初值给得太大,导致指数项几乎为零,梯度消失,算法直接卡住。后来改成先对数据做对数变换线性拟合得到初值,再喂给 NLLS,一次就收敛了。
5.2 病态问题与正则化
当 $\mathbf{J}^T\mathbf{J}$ 的条件数很大时,问题就是病态的。病态的表现是:参数估计对数据微小扰动极其敏感,解不稳定。LM 法的阻尼项本身就有正则化效果,能缓解病态。但如果病态严重,还需要额外手段:
- 参数精简:去掉冗余参数,降低问题维度。
- Tikhonov 正则化:在目标函数里加 $\mu|\mathbf{x}|^2$,等价于在正规方程里加 $\mu\mathbf{I}$。
- 变量投影:把线性参数和线性参数分离,只对非线性参数迭代。
5.3 数值实现中的几个坑
即使算法理论正确,实现细节也会翻车:
- 矩阵求逆不要显式做:解 $(\mathbf{J}^T\mathbf{J} + \lambda\mathbf{I})\Delta\mathbf{x} = -\mathbf{J}^T\mathbf{r}$ 应该用 Cholesky 分解或 QR 分解,而不是先求逆再乘。
- 残差和雅可比要同步更新:每次参数更新后,两者都要重算,不能复用旧值。
- 注意浮点溢出:残差平方和在大数据量下可能溢出,必要时用 Kahan 求和或缩放。
| 问题现象 | 可能原因 | 排查方向 |
|---|---|---|
| 迭代发散 | 初值太远 / $\lambda$ 太小 | 检查初值,增大初始 $\lambda$ |
| 收敛极慢 | $\lambda$ 太大 / 尺度不匹配 | 用 diag 阻尼,调小初始 $\lambda$ |
| 停在明显错误解 | 局部极小 | 多初值重启 |
| 参数剧烈震荡 | 病态 / 参数冗余 | 正则化,精简参数 |
| 残差下降后反弹 | 步长过大 | 引入线搜索或信赖域 |
6. 从理论到代码:一个可复现的 LM 实现骨架
6.1 核心循环的伪代码拆解
把前面所有内容串起来,一个完整的 LM 实现大致是这样的结构:
# 伪代码,展示 LM 核心逻辑 x = x0 # 初值 lam = 1e-3 # 初始阻尼 for k in range(max_iter): r = compute_residual(x) # 计算残差 J = compute_jacobian(x) # 计算雅可比 g = J.T @ r # 梯度 if norm(g, inf) < eps_g: # 梯度判据 break H = J.T @ J # 近似 Hessian # 用 diag 阻尼实现尺度不变 D = diag(diag(H)) delta = solve(H + lam * D, -g) # 解线性方程 x_new = x + delta if sum(compute_residual(x_new)**2) < sum(r**2): x = x_new # 接受 lam = max(lam / 10, 1e-12) # 减小阻尼 else: lam = min(lam * 10, 1e12) # 增大阻尼,拒绝这段代码虽然短,但每一行都对应前面讲的一个决策点:阻尼形式、接受准则、$\lambda$ 更新、终止判据。
6.2 用现成库还是自己写
实际项目里,除非有特殊需求,我一般直接用成熟库:
- Python:
scipy.optimize.least_squares,支持 LM、trf、dogbox 三种方法,自动处理边界约束。 - C++:Ceres Solver、g2o,适合大规模稀疏问题。
- MATLAB:
lsqnonlin,接口友好,适合快速验证。
自己写的价值在于理解原理和定制特殊需求,比如自定义阻尼策略、特殊稀疏结构利用。但如果只是常规拟合,用库能省下大量调试时间。
6.3 验证实现正确性的方法
写完或选好算法后,怎么确认它是对的?我通常做三件事:
- 合成数据测试:用已知参数生成带噪声数据,看算法能否恢复出接近真值的参数。
- 残差单调性检查:正常情况下残差平方和应该单调不增,如果出现上升,说明接受准则有问题。
- 与库对比:同一问题分别用自己的实现和成熟库跑,结果应该接近。
7. 方法选型与实战经验:什么时候用哪个
7.1 GN、LM、信赖域方法的适用边界
三种主流方法各有适用场景:
- GN 法:残差小、模型接近线性、雅可比满秩时首选,速度最快。典型场景是 bundle adjustment 的后期迭代。
- LM 法:通用性最强,绝大多数中小规模非线性最小二乘问题都能用。工程默认选择。
- 信赖域方法:理论最严谨,适合病态、强非线性问题,但实现复杂,通常用库。
选型的核心判断依据是:残差量级、模型非线性程度、参数尺度差异、问题规模。规模大到稀疏性可以利用时,还要考虑稀疏求解器。
7.2 我踩过的三个真实坑
坑一:忽略参数尺度。前面提过,多参数标定时用 $\lambda\mathbf{I}$ 导致不收敛,换 diag 阻尼解决。教训是:只要参数量纲不同,就必须用尺度不变的阻尼。
坑二:残差函数写错符号。残差定义成"观测减预测"还是"预测减观测",会影响雅可比符号,进而影响迭代方向。我因为符号搞反,算法朝反方向优化,越跑越远。后来养成习惯:先用一个简单线性案例验证残差和雅可比符号。
坑三:过早停止。默认的收敛阈值有时太松,算法在残差还很大时就停了。我现在的做法是把阈值调紧一档,同时设最大迭代次数兜底,观察残差曲线是否真的平了。
7.3 性能优化的一些思路
当问题规模变大,性能就成了瓶颈。几个实用方向:
- 利用稀疏性:雅可比通常是稀疏的,用稀疏矩阵存储和求解能大幅提速。
- 并行计算雅可比:每个残差的偏导可以并行算。
- 解析雅可比:比数值差分快一个数量级,值得手推。
- 减少函数调用:残差和雅可比尽量在一次遍历里同时算。
非线性最小二乘这块内容,理论不复杂,但工程细节极多。真正拉开差距的,往往不是知不知道 LM 公式,而是能不能在初值、尺度、收敛判据、数值稳定性这些地方做出正确判断。我个人的体会是:先把 GN 和 LM 的推导自己推一遍,再拿一个真实的小问题从头实现一遍,比看十篇教程都管用。后续如果要做大规模问题,再往稀疏求解和信赖域方向深入,路径会很顺。