最小二乘问题详解8：Levenberg-Marquardt方法

1 引言

对于非线性最小二乘问题的求解来说，除了Gauss-Newton方法（以下简称GN方法）和梯度下降法，另外一种更加实用的求解算法就是Levenberg-Marquardt方法（以下简称LM方法）了。LM方法综合了GN方法和梯度下降法的特性，在实践中表现出极强的鲁棒性和收敛性。在阅读本文之前，至少需要阅读以下三篇前置文章：

1. 《最小二乘问题详解4：非线性最小二乘》
2. 《最小二乘问题详解6：梯度下降法》
3. 《最小二乘问题详解7：正则化最小二乘》

2 求解

2.1 基本原理

先复习一下GN方法的关键点，也就是求解线性最小二乘子问题：

$$\underset{\Delta \theta }{\min }\Vert {\mathbf{r}}_k-J_k\Delta \theta {\Vert }^2$$

正则方程是：

$$J_k^T{J}_k\Delta \theta =J_k^T{\mathbf{r}}_k$$

其解为：

$$\Delta \theta =(J_k^T{J}_k{)}^{-1}{J}_k^T{\mathbf{r}}_k$$

这里可以看到，GN方法有几个致命弱点：

问题	原因
$J_k^T{J}_k$ 可能奇异或病态	当雅可比矩阵 $J_k$ 列相关或接近秩亏时，$(J_k^T{J}_k{)}^{-1}$ 不存在或数值不稳定。
步长过大导致发散	GN 假设局部线性近似足够好，但如果 $\Delta \theta$ 太大，真实残差可能严重偏离线性模型，导致目标函数反而上升。
不适用于远离最优解的情况	在远离极小点时，高阶项不可忽略，一阶近似失效，GN 可能震荡甚至发散。

而LM方法正是为了克服这些问题而设计的，解决思路正是采用正则化最小二乘的思想，在GN方法的正规方程中加入一个正则化项（也被称为阻尼项），使得在 $J^{T}J$ 病态或线性化不良时步长更保守，从而提高稳定性：

$$\boxed{(J_k^T{J}_k+{\lambda }_k{D}_k)\Delta \theta =J_k^T{r}_k}$$

其中：

• $J_k=J({\theta }_k)\in {\mathbb{R}}^{N\times p}$（每个块 $J_i$ 是 $\frac{\partial f(x_i;\theta )}{\partial {\theta }^T})$。
• $r_k=r({\theta }_k)=y-f(x;{\theta }_k)\in {\mathbb{R}}^N$。
• ${\lambda }_k>0$：阻尼参数（damping parameter），控制正则化强度。
• $D_k$ 是对角正定矩阵，常见选择为单位矩阵 $I$ 或 $\mathrm{diag}(J_k^T{J}_k)$。
• $\Delta \theta$：待求的参数增量，LM方法的候选步长。

这个修改看起来简单，但具有很深刻的意义。我们可以观察到随着阻尼项 ${\lambda }_k{D}_k$ 的变化，会自动调节搜索方向：

• 当 ${\lambda }_k\to 0$ 时，退化为标准 Gauss–Newton：$(J^{T}J)\Delta \theta =J^{T}r$，适合接近收敛时使用，增加收敛速度。
• 当 ${\lambda }_k$ 很大且用 $D_k=I$ 时，方程近似为 ${\lambda }_{k}I\Delta \theta \approx J^{T}r$，即 $\Delta \theta \approx \frac{1}{{\lambda }_k}{J}^{T}r$，类似梯度下降小步长方向，适合初始阶段或不稳定情况，稳定但慢。
• 用 $D_k=\mathrm{diag}(J^{T}J)$可实现对不同参数尺度的自适应阻尼。

文章《最小二乘问题详解6：梯度下降法》里的雅可比矩阵是对残差向量$r$求偏导，而这里是雅可比矩阵式是对模型函数 $f(x;\theta )$ 求偏导，两者求偏导的结果方向相反。

所以，LM 实现了在迭代过程中智能平衡：在平坦、可信的区域大胆走 GN 的大步；在崎岖、不可信的区域小心走梯度下降的小步。

2.2 可信度比

既然LM方法可以自动调节搜索方向，那么关键就在于控制调节搜索方向的参数——也就是模型可信度比。在迭代逼近过程中，使用真实非线性模型计算目标函数的 实际减少（actual） 量是：

$$\text{ared}=S({\theta }_k)-S({\theta }_k+\Delta \theta )=\Vert r_k{\Vert }^2-\Vert r({\theta }_k+\Delta \theta ){\Vert }^2.$$

在GN方法中，将迭代过程的线性最小二乘子问题模型展开：

$$m_k(\Delta \theta )\equiv \Vert r_k-J_k\Delta \theta {\Vert }^2=r_k^T{r}_k-2r_k^T{J}_k\Delta \theta +\Delta {\theta }^T(J_k^T{J}_k)\Delta \theta .$$

那么基于线性模型近似的 预期减少（predicted） 量是：

$$\text{pred}=m_k(0)-m_k(\Delta \theta )=2r_k^T{J}_k\Delta \theta -\Delta {\theta }^T(J_k^T{J}_k)\Delta \theta .$$

可定义模型可信度比：

$$\rho =\frac{\text{ared}}{\text{pred}}$$

用来判断步长是否有效。具体来说：

• 若 $\rho$ 大（例如 $\rho >0$ 且远离 0），说明真实下降与模型预测一致或更好，应接受步并减少 $\lambda$；
• 若 $\rho \le 0$ 或很小，说明模型预测不可靠，应拒绝步进并增大 $\lambda$。

其实模型可信度比 $\rho$ 是来源于 信任域(Trust Region) 中的概念，线性模型 $\mathbf{r}(\theta )\approx {\mathbf{r}}_k-J_k\Delta \theta$ 只在某个“信任区域”内有效，阻尼参数 ${\lambda }_k$ 则控制了这个区域的大小。不过这里也就是提一提，笔者理解的也不是很深入，以后有机会再深入探讨。

2.3 算法流程

初始化如下参数：

• 初始参数猜测 ${\theta }_0$
• 初始阻尼参数 ${\lambda }_0$（例如 $10^{-3}$ 或基于 $\text{diag}(J_0^T{J}_0)$）
• 阻尼更新因子 $\nu >1$（常用 $\nu =10$）
• 阻尼项 $D_k$ 选择（$I$ 或 $\text{diag}(J_k^T{J}_k)$）
• 收敛阈值 $ϵ$

进行迭代逼近，对 $k=0,1,2,\dots$

1. 计算残差 $r_k=y-f(x;{\theta }_k)$,维度为 $N$
2. 计算雅可比 $J_k=J({\theta }_k)$ （维度 $N\times p$，$J_k={\frac{\partial f}{\partial {\theta }^T}|}_{\theta ={\theta }_k}$）
3. 构造 $A=J_k^T{J}_k$, $g=J_k^T{r}_k$ （$A$ 是 $p\times p$，$g$ 是 $p\times 1$）
4. 构造阻尼矩阵 $B=A+{\lambda }_k{D}_k$
5. 求解线性系统 $B\Delta \theta =g$，得到候选步进值 $\Delta \theta$，可使用Cholesky分解/LDLT分解
6. 计算实际减少：
   • $S_{old}=||r_k|{|}^2$
   • $S_{new}=||r({\theta }_k+\Delta \theta )|{|}^2$
   • $ared=S_{old}-S_{new}$
7. 计算预测减少: $pred=2g^T\Delta \theta -\Delta {\theta }^{T}A\Delta \theta$
8. 计算模型可信度比 $\rho =ared/pred$
9. 如果 ${\rho }_k>0$，表示更新有效：
   • 接受更新：${\theta }_{k+1}={\theta }_k+\Delta \theta$
   • 如果满足如下收敛条件之一，则停止更新并返回 ${\theta }_k$
     • $\Vert \Delta {\theta }_k\Vert <{ϵ}_1$
     • $\Vert \nabla S({\theta }_k)\Vert =\Vert g\Vert <{ϵ}_2$
     • $|S_{new}-S_{old}|<{ϵ}_3\ast S_{old}$，使用相对变化判据，避免不同尺度下的误判
   • 减小 ${\lambda }_k$（例如 ${\lambda }_{k+1}={\lambda }_k/\nu$），更接近 GN，加快收敛
10. 否则${\rho }_k\le 0$，模型预测失败：
    • 拒绝更新：${\theta }_{k+1}={\theta }_k$
    • 增大 ${\lambda }_k$（例如 ${\lambda }_{k+1}={\lambda }_k\cdot \nu$），更接近梯度下降，更保守

在实践时，有如下问题需要注意：

1. 在更新 $\lambda$ 时加边界保护：
   • ${\lambda }_{k+1}=\text{max}({\lambda }_k/\nu ,{\lambda }_{min})$ # 防止 λ → 0 导致 GN 不稳定
   • ${\lambda }_{k+1}=\text{min}({\lambda }_k\ast \nu ,{\lambda }_{max})$ # 防止 λ → ∞ 导致步长太小
2. 公式 $\rho =ared/pred$ 在 $pred<=0$ 时会导致误判或除零，此时需要将这次逼近视为不可信，对应同 $\rho <=0$ 的处理（拒绝步、增大 $\lambda$）。
3. 推荐初始 $\lambda$ 取法：(\lambda_0 = \tau \cdot \max(\operatorname{diag}(A_0)))，$\tau$ 可取 1e-3 至 1e-1。这样能自动按参数尺度调整初始阻尼。
4. $D_k=\text{diag}(A)$ 通常比 $I$ 更稳健（参数尺度不同会导致不同的步长）。
5. 典型的$ϵ$默认值是：${ϵ}_1=1e-6,{ϵ}_2=1e-8,{ϵ}_3=1e-12$，也可以按照问题尺度进行调整。

3 实例

改进《最小二乘问题详解5：非线性最小二乘求解实例》中的实例，将原来的GN方法改进成LM方法：

#include <Eigen/Dense>
#include <cmath>
#include <iostream>
#include <random>
#include <vector>using namespace std;
using namespace Eigen;// 模型函数: y = exp(a*x^2 + b*x + c)
double model(double x, const Vector3d& theta) {double a = theta(0);double b = theta(1);double c = theta(2);double exponent = a * x * x + b * x + c;// 防止溢出if (exponent > 300) return exp(300);if (exponent < -300) return exp(-300);return exp(exponent);
}// 计算 Jacobian 矩阵（数值导数或解析导数）
MatrixXd computeJacobian(const vector<double>& x_data,const vector<double>& y_data, const Vector3d& theta) {int N = x_data.size();MatrixXd J(N, 3);  // 每行对应一个点，三列：∂f/∂a, ∂f/∂b, ∂f/∂cfor (int i = 0; i < N; ++i) {double x = x_data[i];double exponent = theta(0) * x * x + theta(1) * x + theta(2);double f = model(x, theta);  // 当前预测值// 解析导数（链式法则）double df_de = f;  // d(exp(u))/du = exp(u)double de_da = x * x;double de_db = x;double de_dc = 1.0;J(i, 0) = df_de * de_da;  // ∂f/∂aJ(i, 1) = df_de * de_db;  // ∂f/∂bJ(i, 2) = df_de * de_dc;  // ∂f/∂c}return J;
}// 计算残差向量 r_i = y_i - f(x_i; theta)
VectorXd computeResiduals(const vector<double>& x_data,const vector<double>& y_data, const Vector3d& theta) {int N = x_data.size();VectorXd residuals(N);for (int i = 0; i < N; ++i) {double pred = model(x_data[i], theta);residuals(i) = y_data[i] - pred;}return residuals;
}int main() {// ========================// 1. 设置真实参数// ========================Vector3d true_params;true_params << 0.05, -0.4, 1.0;  // a, b, cdouble a_true = true_params(0), b_true = true_params(1),c_true = true_params(2);cout << "真实参数: a=" << a_true << ", b=" << b_true << ", c=" << c_true<< endl;// ========================// 2. 生成观测数据（带高斯噪声）// ========================int N = 50;vector<double> x_data(N), y_data(N);random_device rd;mt19937 gen(rd());uniform_real_distribution<double> x_dis(-5.0, 5.0);  // x 在 [-5, 5]normal_distribution<double> noise_gen(0.0, 0.1);     // 噪声 ~ N(0, 0.1)for (int i = 0; i < N; ++i) {x_data[i] = x_dis(gen);double y_true = model(x_data[i], true_params);y_data[i] = y_true + noise_gen(gen);  // 添加噪声}// ========================// 3. 初始化参数（随便猜）// ========================Vector3d theta = Vector3d::Zero();  // 初始猜测: a=0, b=0, c=0cout << "初始猜测: a=" << theta(0) << ", b=" << theta(1) << ", c=" << theta(2)<< endl;// ========================// 4. Levenberg-Marquardt 迭代// ========================int max_iter = 50;double tau = 1e-3;double lambda = 1e-3;double nu = 10;double epsilon_1 = 1e-6;double epsilon_2 = 1e-8;double epsilon_3 = 1e-12;cout << "\n开始 Levenberg-Marquardt 迭代...\n";for (int iter = 0; iter < max_iter; ++iter) {// 计算残差 rVectorXd r = computeResiduals(x_data, y_data, theta);// 计算代价函数 ||r||^2double S_old = r.squaredNorm();cout << "迭代 " << iter << ": 残差平方和 = " << S_old << endl;// 计算 Jacobian 矩阵MatrixXd J = computeJacobian(x_data, y_data, theta);// A = J_k^T J_kMatrixXd A = J.transpose() * J;// g = J_k ^ T r_kVectorXd g = J.transpose() * r;// D_kMatrixXd D = A.diagonal().asDiagonal();// 自适应初始阻尼if (iter == 0) {lambda = tau * A.diagonal().maxCoeff();}// B = A + λ_k D_kMatrixXd B = A + lambda * D;// 求解线性系统 BΔθ = gVectorXd delta = B.colPivHouseholderQr().solve(g);// 计算实际减少VectorXd r_new = computeResiduals(x_data, y_data, theta + delta);double S_new = r_new.squaredNorm();double ared = S_old - S_new;// 计算预测减少pred = 2 g^T Δθ - Δθ^T A Δθdouble pred = 2.0 * g.dot(delta) - delta.dot(A * delta);if (pred <= 0) {  // 模型预测无效（可能是数值误差或矩阵病态）cout << "预测减少量 <= 0，拒绝更新" << endl;lambda *= nu;lambda = std::min(lambda, 1e12);  // 防止 lambda 太大continue;}// 模型可信度比double rho = ared / pred;if (rho > 0) {cout << "接受更新" << endl;theta += delta;bool stop = (delta.norm() < epsilon_1) || (g.norm() < epsilon_2) ||(fabs(ared) < epsilon_3 * S_old);if (stop) {break;}lambda /= nu;lambda = std::max(lambda, 1e-10);  // 防止 lambda 太小} else {cout << "拒绝更新" << endl;lambda *= nu;lambda = std::min(lambda, 1e12);  // 防止 lambda 太大}}// ========================// 5. 输出结果// ========================cout << "\n--- 拟合完成 ---" << endl;cout << "估计参数: a=" << theta(0) << ", b=" << theta(1) << ", c=" << theta(2)<< endl;cout << "真实参数: a=" << a_true << ", b=" << b_true << ", c=" << c_true<< endl;// 最终残差VectorXd final_r = computeResiduals(x_data, y_data, theta);cout << "最终残差平方和: " << final_r.squaredNorm() << endl;// ========================// 6. （可选）计算参数协方差与标准差// ========================MatrixXd J_final = computeJacobian(x_data, y_data, theta);int n = N, p = 3;double sigma_squared = final_r.squaredNorm() / (n - p);  // 估计噪声方差MatrixXd cov_theta =sigma_squared * (J_final.transpose() * J_final).inverse();Vector3d std_error;std_error << sqrt(cov_theta(0, 0)), sqrt(cov_theta(1, 1)),sqrt(cov_theta(2, 2));cout << "\n参数标准差 (近似):" << endl;cout << "a: ±" << std_error(0) << endl;cout << "b: ±" << std_error(1) << endl;cout << "c: ±" << std_error(2) << endl;return 0;
}

应该来说，LM算法的关键点全部都已经在第2节中已经说明，也没什么值得额外注意的。运行结果如下：

真实参数: a=0.05, b=-0.4, c=1
初始猜测: a=0, b=0, c=0开始 Levenberg-Marquardt 迭代...
迭代 0: 残差平方和 = 17583
拒绝更新
迭代 1: 残差平方和 = 17583
接受更新
迭代 2: 残差平方和 = 16798.6
拒绝更新
迭代 3: 残差平方和 = 16798.6
接受更新
迭代 4: 残差平方和 = 15697.2
接受更新
迭代 5: 残差平方和 = 4748.15
接受更新
迭代 6: 残差平方和 = 471.045
接受更新
迭代 7: 残差平方和 = 58.1985
接受更新
迭代 8: 残差平方和 = 10.2766
接受更新
迭代 9: 残差平方和 = 0.674626
接受更新
迭代 10: 残差平方和 = 0.356372
接受更新
迭代 11: 残差平方和 = 0.35541
接受更新
迭代 12: 残差平方和 = 0.35541
拒绝更新
迭代 13: 残差平方和 = 0.35541
拒绝更新
迭代 14: 残差平方和 = 0.35541
拒绝更新
迭代 15: 残差平方和 = 0.35541
拒绝更新
迭代 16: 残差平方和 = 0.35541
接受更新--- 拟合完成 ---
估计参数: a=0.0504625, b=-0.396944, c=1.00441
真实参数: a=0.05, b=-0.4, c=1
最终残差平方和: 0.35541参数标准差 (近似):
a: ±0.000332532
b: ±0.00203305
c: ±0.00425019

可以多运行几次看看不同随机数的结果，可以看到改进后的LM算法运行结果非常稳定，基本每次都能收敛；而原来的GN方法总是有一定概率不能收敛。以这个例子来说，LM方法解决了GN方法初值太差、局部线性近似不足导致发散的问题，表现除了极强的稳健性。

4 改进

Levenberg 最早提出在GN方法中加入阻尼项：

$$(J^{T}J+\lambda I)\Delta \theta =J^{T}r$$

其本质是让矩阵始终可逆，并在梯度下降与高斯牛顿之间做插值。但这种方法存在收敛速度慢、阻尼调整不灵活的问题。Marquardt 对其进行了关键改进，使算法在工程实践中更加高效与稳健：

4.1 阻尼项“自适应缩放”

Marquardt 观察到，单纯使用 $I$ 作为阻尼方向可能破坏不同参数的尺度关系。因此，他提出使用矩阵的对角项进行缩放：

$$(J^{T}J+\lambda \mathrm{diag}(J^{T}J))\Delta \theta =J^{T}r$$

这样能让每个参数方向按自身曲率大小进行调节，避免大尺度参数步长过大、小尺度参数步长过小的问题。也就是前面算法中“$D_k=\text{diag}(A)$ 比 $I$ 更稳健”的由来。

4.2 阻尼因子动态调整

Levenberg 原始算法只使用“接受则除以ν、拒绝则乘以ν”的二元调整策略：

$${\lambda }_{k+1}=\{\begin{array}{ll}{\lambda }_k/\nu ,& {\rho }_k>0\\ {\lambda }_k\cdot \nu ,& {\rho }_k\le 0\end{array}$$

而 Marquardt 提出更细腻的自适应方案，使 $\lambda$ 的变化与模型可信度 $\rho$ 连续相关：

$${\lambda }_{k+1}={\lambda }_k\cdot \max \left(\frac{1}{3},1-(2{\rho }_k-1{)}^3\right)$$

同时还引入因子 ${\nu }_{k+1}=2$，当 ${\rho }_k<0.25$ 时再增大 $\lambda$。
这一改进让算法能平滑地在“梯度下降模式”与“高斯牛顿模式”之间过渡，避免震荡与过调。

4.3 多级分级的ρ判定策略

现代实现（如 Ceres Solver、MPFIT、g2o）通常采用分级控制策略来调整 $\lambda$，使其调整幅度更柔和。常见做法如下：

${\rho }_k$ 区间	含义	$\lambda$ 调整策略
${\rho }_k>0.75$	模型拟合非常好	${\lambda }_{k+1}={\lambda }_k/3$
$0.25<{\rho }_k\le 0.75$	拟合良好	${\lambda }_{k+1}={\lambda }_k/2$
$0<{\rho }_k\le 0.25$	拟合一般	${\lambda }_{k+1}={\lambda }_k$（保持）
${\rho }_k\le 0$	拟合失败	${\lambda }_{k+1}={\lambda }_k\cdot \nu$

这种分级调整在实践中能显著提高收敛稳定性，尤其是在存在噪声或残差面高度非线性的情况下。

5 总结

Levenberg–Marquardt 方法的最大魅力在于：它不是在梯度下降和高斯牛顿之间取折中，而是根据模型的“可信度”在两者之间智能切换。这使得 LM 成为现代非线性最小二乘的工业标准算法，也让它成为理解信任域思想的入门之路。