Levenberg-Marquardt（LM）算法是求解非线性最小二乘问题的迭代优化方法，融合梯度下降与高斯 - 牛顿法优势。文章解析其数学原理，包括残差函数、雅可比矩阵、Hessian 近似及阻尼因子动态调节策略。重点讨论参数尺度归一化、收敛判断条件及工业级 C++ 实现细节，提供模块化类设计与核心求解流程代码，适用于参数估计、SLAM 重投影等工程场景。
在现代科学计算、机器人感知、计算机视觉和机器学习中，我们常常面临一个核心挑战：如何从带有噪声的观测数据中，反推出最能解释这些数据的模型参数？这不仅仅是一个数学问题，更是一场关于精度、效率与鲁棒性的系统性博弈。而在这条探索之路上，Levenberg-Marquardt（LM）算法就像一位经验老道的导航员，既能穿越崎岖不平的非线性地形，又能巧妙地在'大胆前进'与'谨慎试探'之间找到最佳平衡。
你可能已经见过它的标准公式：
$$ \mathbf{J}^T\mathbf{J} + \lambda \mathbf{I})\Delta\theta = -\mathbf{J}^T\mathbf{r} $$
但这个简洁的表达背后，隐藏着怎样的直觉？为什么它能在梯度下降和高斯 - 牛顿法之间游刃有余？更重要的是——当我们在 C++ 里真正实现它时，会遇到哪些教科书上不会告诉你的坑？
今天，我们就来一次彻底的拆解，从物理直觉到代码细节，从理论边界到实战技巧，把 LM 算法讲个明明白白。
让我们先放下那些复杂的公式，回到问题的本质。
假设你在做实验，测得了一组 $(x_i, y_i)$ 的数据点，你觉得它们应该符合某种规律，比如指数衰减、抛物线增长，或者某个复杂的动力学方程。你的目标是找到一组参数 $\theta$，让模型 $f(x; \theta)$ 尽可能贴近这些点。
听起来很简单？可一旦模型是非线性的，事情就变得棘手了。
在优化的世界里，'贴近'被翻译成一种叫做**残差（residual）**的东西：
$$ r_i(\theta) = y_i - f(x_i; \theta) $$
每一个数据点对应一个残差，表示当前模型预测值与真实值之间的差距。我们的目标不再是'让所有点都刚好落在曲线上'——那是理想主义；而是'让所有残差的平方和尽可能小'，也就是：
$$ \min_\theta F(\theta) = \frac{1}{2} \sum_{i=1}^m r_i(\theta)^2 = \frac{1}{2} | \mathbf{r}(\theta) |^2 $$
你可能会问：为啥要加个 $\frac{1}{2}$？其实是为了求导方便！当你对 $F(\theta)$ 求梯度时，$\frac{1}{2}$ 和平方项的 2 刚好抵消，得到一个干净的结果：
$$ \nabla F(\theta) = \mathbf{J}^T(\theta)\mathbf{r}(\theta) $$
其中 $\mathbf{J}$ 是雅可比矩阵，记录了每个残差对每个参数的变化率。这个结构太重要了——它意味着我们可以只关心一阶导数，而不必直接处理复杂的二阶 Hessian 矩阵。
📌 小贴士：如果你发现你的目标函数没有这种'半范数平方'的形式，那可能需要重新考虑建模方式。LM 算法最喜欢的就是这种结构清晰、梯度容易计算的问题。
别以为残差只能用来画曲线。在真实世界中，它的形态千变万化：
| 应用领域 | 残差定义 |
|---|---|
| 曲线拟合 | $y_i - (ae^{-bx_i} + c)$ |
| 相机标定 | $ |
| SLAM 重投影 | $(u_i^{\text{obs}} - u_i^{\text{proj}}, v_i^{\text{obs}} - v_i^{\text{proj}})$ |
| 化学反应动力学 | $C_i^{\text{meas}} - C(t_i; k)$ |
无论形式如何变化，关键在于：残差必须是连续可导的。如果你用了 abs()、max() 或者 if/else 条件判断，LM 算法很可能会罢工。
例如，Huber loss 虽然抗噪能力强，但它本身不可导。正确的做法是将其分解为加权残差项，或使用平滑近似如 $\sqrt{x^2 + \epsilon}$ 来代替绝对值。
想象一下你在山谷中徒步。你想尽快到达谷底，但周围雾气弥漫，看不清路。
这就是梯度下降法vs高斯 - 牛顿法的区别。
而 LM 算法呢？它像是一个聪明的徒步者，手里拿着一张模糊的地图和一个灵敏的姿态传感器。他既知道大致方向，也感知脚下地面的坡度变化。于是他决定：
'如果下一步走得稳，我就相信地图，大步向前； 如果踩空了，我就立刻收脚，改用小碎步试探。'
LM 的核心思想体现在那个著名的更新方程：
$$ (\mathbf{J}^T\mathbf{J} + \lambda \mathbf{I})\Delta\theta = -\mathbf{J}^T\mathbf{r} $$
这里的 $\lambda$ 就是那个'信任度调节器'。
所以你看，LM 不是简单地混合两种方法，而是通过一个参数实现了动态平衡。这才是它真正的智慧所在！
// 简化版 LM 步长计算示意（依赖 Eigen 库）
Eigen::VectorXd ComputeLMStep(const Eigen::MatrixXd& J, const Eigen::VectorXd& r, double lambda) {
Eigen::MatrixXd A = J.transpose() * J + lambda * Eigen::MatrixXd::Identity(J.cols(), J.cols());
Eigen::VectorXd b = -J.transpose() * r;
return A.ldlt().solve(b); // 使用 LDLT 分解求解对称正定系统
}
这段代码虽然短，却浓缩了 LM 的精髓。不过注意：这里用了 LDLT 分解，因为它适用于对称正定矩阵，且比 LU 更快更稳定。但在工业级实现中，还会考虑更多因素，比如稀疏性、内存布局等。
现在我们来看看 LM 背后的两个关键技术点：Hessian 近似和阻尼因子调节。
完整的目标函数 Hessian 其实是这样的：
$$ \mathbf{H} = \mathbf{J}^T\mathbf{J} + \sum_{i=1}^m r_i(\theta) \cdot \nabla^2 r_i(\theta) $$
第二项包含了残差函数的二阶导数，理论上更准确。但在实践中，我们通常忽略它，理由如下：
因此，采用 $\mathbf{J}^T\mathbf{J}$ 不仅合理，而且高效。这也是高斯 - 牛顿法和 LM 算法得以广泛应用的基础。
当然，这也带来了副作用：当初始残差很大或模型高度非线性时，这个近似可能失效，导致搜索方向错误。这时候，就得靠阻尼因子 $\lambda$ 来救场了。
固定 $\lambda$ 是行不通的。聪明的做法是根据每次迭代的效果动态调整它。
下面是我在实际项目中常用的策略：
double lambda = 1e-3;
double nu = 2.0;
for (int iter = 0; iter < max_iters; ++iter) {
// 构造带阻尼的矩阵
MatrixXd jtj = jacobian.transpose() * jacobian;
VectorXd diag_jtj = jtj.diagonal();
MatrixXd damping = lambda * diag_jtj.asDiagonal(); // 改进型：按主对角线缩放
MatrixXd A = jtj + damping;
VectorXd b = -jacobian.transpose() * residuals;
VectorXd delta = A.ldlt().solve(b);
// 尝试更新
VectorXd new_params = params + delta;
VectorXd new_residuals = ComputeResiduals(new_params);
double new_cost = new_residuals.squaredNorm();
double old_cost = residuals.squaredNorm();
// 计算预测下降量（来自二次模型）
double predicted_reduction = -(delta.dot(jtj * delta + b)); // 利用泰勒展开估算
// 计算实际下降与预测下降的比率
double rho = (old_cost - new_cost) / predicted_reduction;
if (rho > 1e-3) {
// 更新有效！接受新参数，并减小λ
params = new_params;
residuals = new_residuals;
// 动态缩减λ，加速收敛
lambda *= std::max(1.0/3.0, 1.0 - pow(2*rho - 1, 3));
nu = 2.0; // 重置增长倍数
} else {
// 更新失败！拒绝更新，增大λ以增强稳定性
lambda *= nu;
nu *= 2.0; // 指数增长，快速逃离危险区域
}
// 检查收敛
if (Converged(delta, residuals, old_cost)) break;
}
这套逻辑有几个精妙之处：
rho 是关键指标：它衡量了'我们的局部模型是否靠谱'。如果实际下降远小于预测，说明模型不准，必须保守一点。nu 指数增长：避免陷入无限循环，确保即使连续失败也能迅速跳出困境。💡 经验法则：初始 $\lambda$ 建议设为 $\tau \cdot \max(\mathrm{diag}(\mathbf{J}^T\mathbf{J}))$，其中 $\tau \in [10^{-6}, 10^{-3}]$。维度越高，$\tau$ 越小。
你有没有遇到过这种情况：明明算法跑起来了，但某些参数几乎不动，而另一些疯狂震荡？
这很可能是因为参数尺度差异太大！
举个例子，在 SLAM 中：
如果不做处理，雅可比矩阵的列就会严重不平衡，导致 $\mathbf{J}^T\mathbf{J}$ 病态，条件数巨大，求解极不稳定。
最直接的方法是对参数进行预处理：
$$ \theta'_j = \frac{\theta_j}{s_j} $$
其中 $s_j$ 是第 $j$ 个参数的典型尺度。例如：
这样所有参数都在相似的数量级上，优化过程更加平稳。
更优雅的方式是在 LM 方程中用 $\mathrm{diag}(\mathbf{J}^T\mathbf{J})$ 代替 $\mathbf{I}$：
$$ (\mathbf{J}^T\mathbf{J} + \lambda \cdot \mathrm{diag}(\mathbf{J}^T\mathbf{J})) \Delta\theta = -\mathbf{J}^T\mathbf{r} $$
这意味着每个参数的正则化强度与其自身的曲率成比例，实现了自动的'各向异性阻尼'。
这一技巧已被广泛应用于 Ceres Solver、g2o 等主流框架中，强烈推荐在跨尺度系统中使用。
优化不能无限循环下去。我们需要一套可靠的终止条件。
常见的几种判据包括：
$$ | \Delta\theta | < \epsilon_1 (|\theta| + \epsilon_1) $$
这是一个相对阈值，避免因参数尺度不同而导致误判。
bool CheckParamConvergence(
const VectorXd& delta,
const VectorXd& current,
double eps = 1e-8)
{
double norm_delta = delta.norm();
double norm_param = current.norm();
return norm_delta < eps * (norm_param + eps);
}
$$ \frac{|F_k - F_{k+1}|}{|F_k|} < \epsilon_2 $$
适用于残差本身较大的情况。
$$ | \nabla F | = | \mathbf{J}^T\mathbf{r} | < \epsilon_3 $$
表明已到达极小值附近。
bool HasConverged(
double cost,
double prev_cost,
const VectorXd& grad,
const VectorXd& delta)
{
double cost_diff = std::abs(cost - prev_cost);
double rel_cost_change = cost_diff / (std::abs(prev_cost) + 1e-8);
bool by_cost = rel_cost_change < 1e-6;
bool by_grad = grad.norm() < 1e-5;
bool by_step = delta.norm() < 1e-6 * (params.norm() + 1e-6);
return by_cost || by_grad || by_step;
}
同时设置最大迭代次数作为兜底机制：
int iter = 0;
while (!HasConverged() && iter < max_iter) {
// 主循环体
iter++;
}
if (iter >= max_iter) {
LOG(WARNING) << "LM terminated by max iterations";
}
一般默认设为 50~100 次，视问题复杂度调整。
说了这么多，终于到了动手环节！
下面是一个模块化、可复用的 LM 类设计：
class LevenbergMarquardt {
public:
using ResidualFunc = std::function<void(const Eigen::VectorXd&, Eigen::VectorXd&)>;
using JacobianFunc = std::function<void(const Eigen::VectorXd&, Eigen::MatrixXd&)>;
LevenbergMarquardt(ResidualFunc res_func, JacobianFunc jac_func)
: residual_func_(std::move(res_func)), jacobian_func_(std::move(jac_func)) {}
bool solve(Eigen::VectorXd& parameters);
void setMaxIterations(int max_iter) { max_iterations_ = max_iter; }
void setTolerance(double tol) { tolerance_ = tol; }
void setLambda(double lambda) { lambda_ = lambda; }
private:
ResidualFunc residual_func_;
JacobianFunc jacobian_func_;
int max_iterations_ = 100;
double tolerance_ = 1e-6;
double lambda_ = 0.01;
double lambda_factor_ = 10.0;
Eigen::VectorXd current_params_;
Eigen::VectorXd residuals_;
Eigen::MatrixXd jacobian_;
};
核心求解流程如下：
bool LevenbergMarquardt::solve(Eigen::VectorXd& parameters) {
current_params_ = parameters;
double last_cost = 0.0;
for (int iter = 0; iter < max_iterations_; ++iter) {
// 计算残差和雅可比
residual_func_(current_params_, residuals_);
jacobian_func_(current_params_, jacobian_);
// 构造矩阵
Eigen::MatrixXd jtj = jacobian_.transpose() * jacobian_;
Eigen::VectorXd b = -jacobian_.transpose() * residuals_;
// 添加阻尼
Eigen::MatrixXd A = jtj + lambda_ * Eigen::MatrixXd::Identity(jtj.rows(), jtj.cols());
// 求解增量
Eigen::VectorXd delta = A.ldlt().solve(b);
// 尝试更新
Eigen::VectorXd new_params = current_params_ + delta;
Eigen::VectorXd new_residuals;
residual_func_(new_params, new_residuals);
double new_cost = new_residuals.squaredNorm();
double old_cost = residuals_.squaredNorm();
// 判断是否接受
if (new_cost < old_cost) {
current_params_ = new_params;
residuals_ = new_residuals;
last_cost = new_cost;
lambda_ /= lambda_factor_; // 减小阻尼
} else {
lambda_ *= lambda_factor_; // 增加阻尼
}
// 检查收敛
if (delta.norm() < tolerance_) break;
}
parameters = current_params_;
return true;
}
通过上述模块化的设计与实现，你可以将 LM 算法灵活集成到自己的 C++ 项目中，解决各类非线性优化问题。
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
生成新的随机RSA私钥和公钥pem证书。 在线工具，RSA密钥对生成器在线工具，online
基于 Mermaid.js 实时预览流程图、时序图等图表，支持源码编辑与即时渲染。 在线工具，Mermaid 预览与可视化编辑在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online
将 Markdown（GFM）转为 HTML 片段，浏览器内 marked 解析；与 HTML转Markdown 互为补充。 在线工具，Markdown转HTML在线工具，online