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

技术分享 4个月前 (11-05) 0 999+

1 引言

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

2 求解

2.1 基本原理

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

[min_{Delta theta} | mathbf{r}_k - J_k Delta theta |^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)Deltatheta = J_k^T r_k} ]

其中：

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

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

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

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

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

2.2 可信度比

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

[text{ared} = S(theta_k) - S(theta_k+Deltatheta) = |r_k|^2 - |r(theta_k+Deltatheta)|^2. ]

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

[m_k(Deltatheta) equiv | r_k - J_k Deltatheta |^2 = r_k^T r_k - 2 r_k^T J_k Deltatheta + Deltatheta^T (J_k^T J_k)Deltatheta. ]

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

[text{pred} = m_k(0) - m_k(Deltatheta) = 2 r_k^T J_k Deltatheta - Deltatheta^T (J_k^T J_k)Deltatheta. ]

可定义模型可信度比：

[rho = dfrac{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))）
收敛阈值 (epsilon)

进行迭代逼近，对 (k = 0, 1, 2, dots)

计算残差 (r_k = y - f(x; θ_k)),维度为 (N)
计算雅可比 (J_k = J(θ_k)) （维度 (N×p)，(J_k = left.frac{partial f}{partial theta^T}right|_{theta = theta_k})）
构造 (A = J_k^T J_k), (g = J_k^T r_k) （(A) 是 (p×p)，(g) 是 (p×1)）
构造阻尼矩阵 (B = A + λ_k D_k)
求解线性系统 (B Δθ = g)，得到候选步进值 (Δθ)，可使用Cholesky分解/LDLT分解
计算实际减少：
- (S_{old} = ||r_k||²)
- (S_{new} = ||r(θ_k + Δθ)||²)
- (ared = S_{old} - S_{new})
计算预测减少: (pred = 2 g^T Δθ - Δθ^T A Δθ)
计算模型可信度比 (ρ = ared / pred)
如果 (rho_k > 0)，表示更新有效：
- 接受更新：(θ_{k+1} = θ_k + Δθ)
- 如果满足如下收敛条件之一，则停止更新并返回 (θ_k)
  - (|Delta theta_k| < epsilon_1)
  - (|nabla S(theta_k)| = |g| < epsilon_2)
  - (|S_{new} - S_{old}| < epsilon_3 * S_{old})，使用相对变化判据，避免不同尺度下的误判
- 减小 (lambda_k)（例如 (lambda_{k+1} = lambda_k / nu)），更接近 GN，加快收敛
否则(rho_k leq 0)，模型预测失败：
- 拒绝更新：(θ_{k+1} = θ_k)
- 增大 (lambda_k)（例如 (lambda_{k+1} = lambda_k cdot nu)），更接近梯度下降，更保守

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

在更新 (lambda) 时加边界保护：
- (λ_{k+1} = text{max}(λ_k / ν, λ_{min})) # 防止 λ → 0 导致 GN 不稳定
- (λ_{k+1} = text{min}(λ_k * ν, λ_{max})) # 防止 λ → ∞ 导致步长太小
公式 (ρ = ared / pred) 在 (pred <= 0) 时会导致误判或除零，此时需要将这次逼近视为不可信，对应同 (ρ <= 0) 的处理（拒绝步、增大 (λ)）。
推荐初始 (λ) 取法：(lambda_0 = tau cdot max(operatorname{diag}(A_0)))，(τ) 可取 1e-3 至 1e-1。这样能自动按参数尺度调整初始阻尼。
(D_k = text{diag}(A)) 通常比 (I) 更稳健（参数尺度不同会导致不同的步长）。
典型的(epsilon)默认值是：(epsilon_1=1e-6, epsilon_2=1e-8, epsilon_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/∂c    for (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/∂a     J(i, 1) = df_de * de_db;  // ∂f/∂b     J(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, c   double 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=0   cout << "初始猜测: 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) {     // 计算残差 r     VectorXd r = computeResiduals(x_data, y_data, theta);      // 计算代价函数 ||r||^2     double S_old = r.squaredNorm();     cout << "迭代 " << iter << ": 残差平方和 = " << S_old << endl;      // 计算 Jacobian 矩阵     MatrixXd J = computeJacobian(x_data, y_data, theta);      // A = J_k^T J_k     MatrixXd A = J.transpose() * J;      // g = J_k ^ T r_k     VectorXd g = J.transpose() * r;      // D_k     MatrixXd D = A.diagonal().asDiagonal();      // 自适应初始阻尼     if (iter == 0) {       lambda = tau * A.diagonal().maxCoeff();     }      // B = A + λ_k D_k     MatrixXd B = A + lambda * D;      // 求解线性系统 BΔθ = g     VectorXd 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)Deltatheta = J^T r ]

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

4.1 阻尼项“自适应缩放”

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

[(J^T J + lambda operatorname{diag}(J^T J)) Deltatheta = J^T r ]

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

4.2 阻尼因子动态调整

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

[lambda_{k+1} = begin{cases} lambda_k / nu, & rho_k > 0 \ lambda_k cdot nu, & rho_k le 0 \ end{cases} ]

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

[lambda_{k+1} = lambda_k cdot maxleft(frac{1}{3}, 1 - (2rho_k - 1)^3right) ]

同时还引入因子 (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)