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

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)\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\) 或 \(\operatorname{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=\operatorname{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) = \|r_k\|^2 - \|r(\theta_k+\Delta\theta)\|^2. \]

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

\[m_k(\Delta\theta) \equiv \| r_k - J_k \Delta\theta \|^2 = r_k^T r_k - 2 r_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) = 2 r_k^T J_k \Delta\theta - \Delta\theta^T (J_k^T J_k)\Delta\theta. \]

可定义模型可信度比：

\[\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)\Delta\theta = J^T r \]

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

4.1 阻尼项“自适应缩放”

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

\[(J^T J + \lambda \operatorname{diag}(J^T J)) \Delta\theta = 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 \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 成为现代非线性最小二乘的工业标准算法，也让它成为理解信任域思想的入门之路。

posted @ 2025-11-04 20:55 charlee44 阅读(24) 评论(1) 收藏举报

刷新页面返回顶部

人生海海山山而川

分说，不分说，不由分说！

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

1 引言

2 求解

2.1 基本原理

2.2 可信度比

2.3 算法流程

3 实例

4 改进

4.1 阻尼项“自适应缩放”

4.2 阻尼因子动态调整

4.3 多级分级的ρ判定策略

5 总结

公告

人生海海 山山而川

分说，不分说，不由分说！

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

1 引言

2 求解

2.1 基本原理

2.2 可信度比

2.3 算法流程

3 实例

4 改进

4.1 阻尼项“自适应缩放”

4.2 阻尼因子动态调整

4.3 多级分级的ρ判定策略

5 总结

公告

人生海海山山而川