1. 为什么非线性最小二乘值得单独整理如果你做过视觉SLAM、点云配准、曲线拟合、相机标定或者任何需要让模型输出逼近真实观测的工程任务那你大概率已经和非线性最小二乘打过照面了。它要解决的问题非常朴素给定一组观测数据和一个带参数的模型找一组参数让模型预测和真实观测之间的残差平方和最小。线性情形有闭式解一个正规方程就搞定了但一旦模型对参数是非线性的事情就变得麻烦起来——没有闭式解只能迭代逼近。我整理这份东西的初衷是因为在实际项目里反复踩坑同一个问题用高斯-牛顿能收敛换一组初值就发散LM调了半天阻尼因子结果发现是雅可比矩阵的数值精度在捣鬼好不容易收敛了结果掉进局部极小参数物理意义完全不对。这些坑单看某一篇论文或某一个库的文档是串不起来的必须把方法体系从头捋一遍。这份整理适合三类人一是刚接触优化、需要把理论落到代码上的工程师二是做视觉/机器人/测绘天天和BA、配准打交道但没系统梳理过方法的人三是想自己手写一个轻量优化器、不想被庞大框架绑架的开发者。我会从问题定义讲到主流方法高斯-牛顿、LM、Dog-Leg再到雅可比计算、收敛判据、工程实现细节最后给一份常见问题速查表。全程用大白话加公式加代码尽量让你看完就能动手。2. 问题定义与整体思路拆解2.1 非线性最小二乘到底在优化什么先把数学形式摆清楚。我们有一组观测通常写成残差向量的形式r(x) [r1(x), r2(x), ..., rm(x)]^T其中 x 是 n 维待优化参数每个 ri(x) 是第 i 个观测对应的残差预测值减观测值或者反过来符号约定统一即可。目标函数是F(x) 0.5 * ||r(x)||^2 0.5 * sum(ri(x)^2)前面那个 0.5 是习惯求导时能约掉平方的 2让梯度形式更干净。非线性最小二乘的核心任务就是找 x使 F(x) 最小*。为什么强调非线性因为如果 r(x) Ax - b 是线性的那 F(x) 是二次函数直接令梯度为零得到正规方程 A^T A x A^T b一步求解。非线性情况下F(x) 可能有很多局部极小、鞍点梯度为零只是必要条件而且方程本身也解不出来只能迭代。这里有个关键认知非线性最小二乘本质上是无约束优化问题的一个特例只不过目标函数有平方和这个特殊结构。这个结构非常重要它让我们可以利用雅可比矩阵而不是完整的海森矩阵省下大量计算。后面讲高斯-牛顿时会重点说这一点。2.2 迭代法的通用框架所有非线性最小二乘方法都遵循同一个套路我把它抽象成五步给定初值 x0这一步的选取极其关键后面会专门讲。在当前点计算残差 r(xk) 和雅可比 J(xk)J 是 m×n 矩阵第 i 行第 j 列是 ∂ri/∂xj。求解一个子问题得到增量 Δx不同方法的区别几乎全在这一步。更新 x_{k1} xk Δx或 xk α·Δx带线搜索。判断收敛满足条件就停否则回到第 2 步。这个框架看起来简单但每一步都有讲究。初值选不好直接发散雅可比算错收敛慢甚至不收敛子问题求解方式决定了收敛速度和鲁棒性收敛判据设不好要么提前停要么死循环。我见过太多项目在这几个环节出问题所以后面会逐个拆。2.3 方法选型的核心权衡速度 vs 鲁棒性非线性最小二乘的方法谱系本质上是在两个极端之间找平衡高斯-牛顿GN在残差接近零或者问题接近线性时收敛速度是二次的非常快。但它对初值敏感雅可比病态时增量会爆炸残差大时可能根本不下降。梯度下降最鲁棒保证下降方向但收敛是线性的慢得让人抓狂尤其在最优点附近几乎不动。列文伯格-马夸尔特LM在高斯-牛顿和梯度下降之间插一个阻尼因子残差大时偏向梯度下降稳残差小时偏向高斯-牛顿快。这是工程上最常用的默认选择。Dog-Leg另一种信赖域方法结合了高斯-牛顿步和梯度下降步用一条折线找信赖域内的最优步。理论上比 LM 更优雅实际中也有不少库在用。选型的经验法则如果你不确定用哪个先用 LM。它对初值的要求比 GN 宽松得多收敛也够快。如果你的问题残差能降到很小、雅可比条件数也好那 GN 会更快。Dog-Leg 适合那些 LM 调参调到崩溃的场景但实现复杂度更高。3. 核心方法深度解析与实操要点3.1 高斯-牛顿法快但别乱用高斯-牛顿的推导很漂亮。对残差做一阶泰勒展开r(x Δx) ≈ r(x) J(x) Δx代进目标函数F(x Δx) ≈ 0.5 * ||r(x) J Δx||^2对 Δx 求导令其为零得到高斯-牛顿方程(J^T J) Δx -J^T r这个方程的形式和线性最小二乘的正规方程一模一样左边 J^T J 是 n×n 的近似海森矩阵右边 -J^T r 是梯度。解出来就是 GN 步。为什么 GN 快因为它用 J^T J 近似了真实海森矩阵丢掉了残差的二阶导数项。当残差接近零时这个近似非常准所以收敛接近二次。但反过来当残差很大时丢掉的项不能忽略GN 步可能不是下降方向甚至让目标函数变大。实操中 GN 有几个致命陷阱J^T J 奇异或接近奇异比如参数冗余、观测不足、或者某些参数在当前点不可辨识。这时候解方程会得到巨大的 Δx直接把参数甩飞。解决办法是加一个小的对角正则项本质上是往 LM 靠或者做参数化精简。初值离最优解太远GN 没有全局收敛保证初值不好直接发散。我踩过的坑是相机标定里焦距初值给错GN 迭代三次就 NaN 了。残差量纲差异大比如一个残差是像素级几十另一个是角度零点几J^T J 的条件数会很难看。这时候要么归一化残差要么上 LM。提示GN 适合作为最后精修的手段。先用 LM 或带线搜索的方法把参数拉到最优点附近再切 GN 快速收敛这个组合在实际项目里非常好用。3.2 列文伯格-马夸尔特工程默认选择LM 的核心思想是给 GN 方程加一个阻尼项(J^T J λ I) Δx -J^T r或者更一般地用对角矩阵 D 代替 I(J^T J λ D) Δx -J^T rλ 是阻尼因子D 通常取 J^T J 的对角元素这叫 Marquardt 缩放比单位阵效果好因为它考虑了各参数的量纲差异。λ 的调节逻辑是 LM 的灵魂λ → 0退化成 GN快但激进。λ → ∞Δx ≈ -(1/λ) J^T r方向变成梯度下降方向步长很小稳但慢。实际算法里λ 是自适应调整的。一个典型的策略是用当前 λ 算 Δx计算实际的目标函数下降量。如果下降量大于某个阈值说明这一步有效接受更新并减小 λ比如除以 3 或 10让下一步更接近 GN。如果下降量不够甚至上升拒绝更新增大 λ比如乘以 3 或 10重新算 Δx。这个试探-接受/拒绝-调 λ的循环就是 LM 鲁棒性的来源。它保证了每次接受的步都让目标函数下降在数值精度允许范围内所以不会像 GN 那样乱飞。LM 的实操要点初始 λ 的选取常见做法是取 J^T J 对角元素的某个比例比如 λ0 τ · max(diag(J^T J))τ 取 1e-3 到 1e-6。取太小前期可能不稳取太大前期收敛慢。λ 的更新因子增大用 10减小用 0.1 是经典配置。也有用 2 和 0.5 的收敛更平滑但迭代次数多。我一般用 10/0.1配合最大迭代次数兜底。收敛判据不能只看 Δx 小还要看梯度小、目标函数变化小。三个判据同时用比较稳。数值稳定性J^T J λD 的求解用 Cholesky 分解比直接求逆快且稳。如果 Cholesky 失败矩阵非正定说明 λ 太小直接增大 λ 重试。3.3 Dog-Leg 与信赖域思路Dog-Leg 属于信赖域方法思路和高斯-牛顿、LM 都不太一样。它不直接调阻尼因子而是先定义一个信赖域半径 Δ认为在这个半径内二阶泰勒近似是可信的然后在这个球内找最优步。Dog-Leg 的名字来自它的求解路径它构造一条从当前点出发的折线先沿梯度下降方向走再折向高斯-牛顿步像狗腿一样。具体来说计算高斯-牛顿步 Δx_gn。如果 ||Δx_gn|| ≤ Δ直接用 GN 步说明 GN 步在信赖域内可信。否则计算梯度下降方向的最速下降步 Δx_sd通常取到信赖域边界然后在这两个方向张成的二维子空间里找最优步。Dog-Leg 的理论保证比 LM 更清晰它明确控制每步的可信程度收敛性分析更完整。但实现上比 LM 复杂需要处理信赖域半径的更新、折线求交等细节。实际工程里LM 因为实现简单、效果够好用得更多Dog-Leg 在一些对收敛性要求极高的库比如某些标定工具里能见到。信赖域半径 Δ 的更新逻辑和 LM 的 λ 类似但方向相反如果实际下降和预测下降的比值 ρ 接近 1说明近似很准扩大 Δ如果 ρ 很小甚至为负说明近似不准缩小 Δ。这个 ρ 叫增益比是信赖域方法的核心指标。3.4 三种方法的对比与选型表方法收敛速度鲁棒性实现复杂度适用场景高斯-牛顿接近二次残差小时差低初值好、残差小、精修阶段列文伯格-马夸尔特介于线性与二次之间好中通用默认初值一般Dog-Leg接近二次信赖域内好高对收敛性要求高、LM 调参困难选型建议新项目直接上 LM把 λ 的自适应逻辑写对90% 的场景够用。如果发现 LM 收敛慢检查是不是 λ 初始值太大或者更新因子太保守。如果 LM 都不收敛那大概率是雅可比算错了或者模型本身有问题换 Dog-Leg 也救不了。4. 雅可比矩阵算不对全白搭4.1 解析雅可比 vs 数值雅可比雅可比矩阵 J 是整套方法的命脉。它的计算方式直接决定收敛速度和稳定性。两条路解析雅可比手推导数公式代码里直接写。快、准但推导容易出错模型一改就得重推。数值雅可比用有限差分近似∂ri/∂xj ≈ (ri(x h·ej) - ri(x)) / h。实现简单但精度受 h 影响计算量是解析的 n 倍每个参数都要扰动一次。我的经验是原型阶段用数值雅可比快速验证性能敏感或精度要求高的场景换解析雅可比。数值雅可比的 h 选取很讲究太大截断误差大太小舍入误差大。理论最优 h 大约是 sqrt(ε) · |xj|ε 是机器精度double 约 2.2e-16所以 h 大概在 1e-8 量级。但实际中参数尺度差异大我一般用 h 1e-6 · max(|xj|, 1)实测比较稳。中心差分比前向差分精度高一个量级∂ri/∂xj ≈ (ri(x h·ej) - ri(x - h·ej)) / (2h)代价是计算量翻倍。如果残差函数计算不贵中心差分很划算。4.2 雅可比的计算成本与稀疏性利用在 BA光束法平差这类问题里参数可能上万观测几十万雅可比是巨大的稀疏矩阵。如果当成稠密矩阵算和存内存直接爆。这时候必须利用稀疏性J^T J 的稀疏结构BA 里 J^T J 通常是块对角加边箭头形可以用 Schur 补先消元把大问题拆成小问题。雅可比的块结构每个观测只依赖少数几个参数J 的每一行只有几个非零块。存储和计算都按块来。这些优化能把 BA 的复杂度从 O(n^3) 降到接近 O(n)是工程上必须掌握的。不过这是个大话题这里先点一下核心思想就是别把稀疏问题当稠密算。4.3 数值雅可比的步长选择实操我写过一个数值雅可比的调试工具专门用来对比解析和数值的差异。做法是对每个参数分别用解析和数值算导数打印相对误差。如果某个参数的误差超过 1e-4基本可以断定解析推导有错。步长选择的实操建议先用 h 1e-6 试看收敛情况。如果收敛慢或震荡试试 1e-8 和 1e-4对比迭代次数。对量纲差异大的参数用相对步长 h · |xj| 而不是绝对步长。中心差分优先除非计算成本实在扛不住。注意数值雅可比在参数接近零时容易出问题因为相对步长退化成绝对步长。这时候要么给参数加一个偏移要么用绝对步长兜底。5. 完整实操流程与关键环节实现5.1 从零实现一个 LM 优化器我用 Python 写一个最小可用的 LM把前面讲的逻辑串起来。假设残差函数和雅可比函数由用户提供。import numpy as np def lm_optimize(residual, jacobian, x0, max_iter100, lambda01e-3, nu10.0, tol_grad1e-8, tol_step1e-8, tol_cost1e-10): x np.array(x0, dtypefloat) n len(x) lam lambda0 r residual(x) cost 0.5 * np.dot(r, r) for it in range(max_iter): J jacobian(x) g J.T r # 梯度 H J.T J # 近似海森 D np.diag(np.diag(H)) # Marquardt 缩放 D np.where(D 0, 1.0, D) # 防止除零 # 收敛判据梯度足够小 if np.linalg.norm(g, np.inf) tol_grad: break improved False for _ in range(20): # 内层调 lambda A H lam * D try: dx np.linalg.solve(A, -g) except np.linalg.LinAlgError: lam * nu continue x_new x dx r_new residual(x_new) cost_new 0.5 * np.dot(r_new, r_new) # 实际下降 vs 预测下降 predicted -g dx - 0.5 * dx H dx actual cost - cost_new rho actual / predicted if predicted ! 0 else -1 if rho 0 and cost_new cost: x x_new r r_new cost cost_new lam max(lam / nu, 1e-12) improved True break else: lam * nu if not improved: break # 无法下降退出 # 收敛判据步长和目标变化 if np.linalg.norm(dx) tol_step * (np.linalg.norm(x) tol_step): break if abs(cost - cost_new) tol_cost: break return x, cost, it这段代码有几个关键点值得说增益比 rho 的计算预测下降用二阶泰勒近似 -g·dx - 0.5·dx^T H dx实际下降是 cost - cost_new。rho 接近 1 说明近似准接近 0 或负说明不准。这个指标比单纯看 cost 是否下降更有信息量。内层循环调 λ最多试 20 次每次失败就乘 nu。如果 20 次都不行说明当前点附近找不到下降方向直接退出。D 矩阵的处理用 J^T J 的对角元素做缩放但要对零元素兜底否则除零。收敛判据三重梯度、步长、目标变化任一满足就停。实际中梯度判据最可靠步长判据容易在病态问题里误判。5.2 一个曲线拟合的完整案例拿一个经典的指数拟合练手。模型是 y a · exp(b · t) c观测数据带噪声。参数 x [a, b, c]。残差函数def residual(x): a, b, c x return a * np.exp(b * t) c - y_obs def jacobian(x): a, b, c x e np.exp(b * t) J np.zeros((len(t), 3)) J[:, 0] e # d/da J[:, 1] a * t * e # d/db J[:, 2] 1.0 # d/dc return J这个问题的坑在于b 的初值极其敏感。如果 b 初值符号错了指数项要么爆炸要么衰减到零LM 很难拉回来。我的做法是先用对数变换或者网格搜索给一个粗略初值再上 LM 精修。实测下来初值 [a1, b0.1, c0] 在大多数数据上能收敛迭代 15 到 30 次最终 cost 在 1e-6 量级。如果初值给 [1, -0.1, 0]大概率发散或者收敛到错误的局部极小。5.3 收敛判据的设置与调试收敛判据设不好是很多优化不收敛问题的真凶。我总结了几条梯度判据||g||_inf 1e-8 是常用值。但注意如果残差量纲很大比如像素级梯度天然就大这个阈值要相应放大。更好的做法是用相对判据比如 ||g|| / (1 |cost|) 1e-8。步长判据||dx|| 1e-8 · (||x|| 1e-8)。加 1e-8 是防止 x 接近零时判据失效。目标函数判据|cost_new - cost| 1e-10 · (1 cost)。同样用相对形式。最大迭代次数一定要设100 到 200 是常见值。防止死循环。调试技巧把每次迭代的 cost、||g||、||dx||、λ 都打印出来。如果 cost 震荡不降看 λ 是不是太小如果 cost 几乎不变但 ||g|| 还很大看雅可比是不是算错了如果 ||dx|| 突然变得巨大看 J^T J 是不是奇异了。6. 常见问题与排查技巧实录6.1 不收敛、发散、震荡的排查路径我把实际项目中遇到的典型问题整理成一张速查表现象可能原因排查方法解决手段迭代几步后 NaN步长爆炸J^T J 奇异打印cost 震荡不降λ 太小步太激进看 rho 是否频繁为负增大 λ 更新因子加线搜索cost 几乎不变雅可比错误或参数不可辨识对比解析/数值雅可比修正雅可比精简参数收敛到错误解局部极小或初值差多组初值对比换初值加先验约束收敛慢λ 太大或缩放不当看 λ 变化曲线调小 λ0用 Marquardt 缩放参数物理意义不对参数冗余或过拟合看 J^T J 条件数加正则减参数这张表是我踩坑踩出来的基本覆盖了 80% 的现场问题。排查的第一原则是打印中间量别靠猜。cost、梯度范数、步长范数、λ、rho这五个量打出来问题基本就定位了。6.2 初值选取的实战经验初值选取没有万能公式但有几条通用策略物理意义约束参数有物理含义时初值必须落在合理范围。比如焦距不能为负深度不能为零。网格搜索参数少2 到 3 个时粗网格搜一遍取 cost 最小的点做初值。成本可接受效果很好。分步求解先固定一部分参数解另一部分再联合优化。比如标定里先解内参再解外参。多初值重启随机撒几组初值各自跑 LM取最优。适合参数少、单次优化快的场景。利用先验如果有历史数据或相似问题的解直接拿来当初值往往比随机猜好得多。我个人的习惯是任何非线性优化任务先花 20% 的时间在初值上。初值选对了后面 80% 的调参时间都能省下来。6.3 数值精度与尺度问题的避坑数值精度问题很隐蔽但杀伤力大。几个典型场景残差量纲差异大比如同时优化像素误差和角度误差。解决方法是给每个残差乘一个权重让它们的量级接近。权重怎么定看你对哪个更在意或者用观测噪声的倒数。参数尺度差异大比如一个参数是 1e6另一个是 1e-3。J^T J 的条件数会爆炸。解决方法是参数归一化或者用 Marquardt 缩放D 取对角元素。数值雅可比步长不当前面讲过h 取 1e-6 到 1e-8中心差分优先。累加误差迭代几百次后浮点误差累积可能让 cost 不再下降。这时候用双精度或者定期重新计算残差。提示如果 J^T J 的条件数超过 1e12基本可以断定有尺度问题或者参数冗余。先做参数归一化再考虑优化算法。6.4 局部极小与全局优化的取舍非线性最小二乘只能保证收敛到局部极小。全局最优需要额外手段多初值重启最简单粗暴适合低维问题。模拟退火/遗传算法全局搜索能力强但慢通常只用来找初值再用 LM 精修。凸松弛某些特定问题如旋转平均有凸松弛方法能保证全局最优但适用范围窄。问题重构换个参数化方式让目标函数更接近凸。比如用四元数代替欧拉角避免万向锁和周期跳变。实际工程里多初值重启 LM 精修是性价比最高的组合。全局优化算法除非问题特别难否则不划算。7. 工程落地的一些个人体会写到这里方法体系基本讲完了。最后分享几个我在实际项目里反复验证过的经验。第一别迷信算法先检查模型和雅可比。我遇到过好几次LM 不收敛最后发现是雅可比推导错了一个符号。花十分钟对比解析和数值雅可比比调一天 λ 有用得多。第二λ 的自适应逻辑比初始值重要。初始 λ 大一点没关系只要更新逻辑对几步就能降下来。但如果更新因子太保守比如用 2 和 0.5收敛会慢很多。我一般用 10 和 0.1配合增益比判断实测很稳。第三收敛判据要松紧适度。太松提前停结果精度不够太紧死循环浪费时间。我的默认配置是梯度 1e-8、步长 1e-8、目标 1e-10配合最大 100 次迭代覆盖大多数场景。第四把优化过程可视化。cost 随迭代的曲线、参数轨迹、残差分布画出来一眼就能看出问题。纯看数字容易漏掉模式。第五留好退路。优化失败时返回最后一次有效的参数和 cost别返回 NaN。工程代码里鲁棒性比理论最优更重要。这个内容后续还可以往几个方向扩展一是稀疏雅可比和 Schur 补在 BA 里的具体实现二是鲁棒核函数Huber、Cauchy怎么处理外点三是自动微分替代手推雅可比。这几个都是实际项目里绕不开的话题有机会再单独整理。