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

1 引言

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

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

2 求解

2.1 基本原理

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

正则方程是：

其解为：

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

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

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

其中：

• $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） 量是：

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

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

可定义模型可信度比：

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

• 若 $\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$

1. 计算残差 $r_k = y - f(x; θ_k)$,维度为 $N$
2. 计算雅可比 $J_k = J(θ_k)$ （维度 $N×p$，$J_k = \left.\frac{\partial f}{\partial \theta^T}\right|_{\theta = \theta_k}$）
3. 构造 $A = J_k^T J_k$, $g = J_k^T r_k$ （$A$ 是 $p×p$，$g$ 是 $p×1$）
4. 构造阻尼矩阵 $B = A + λ_k D_k$
5. 求解线性系统 $B Δθ = g$，得到候选步进值 $Δθ$，可使用Cholesky分解/LDLT分解
6. 计算实际减少：
   • $S_{old} = ||r_k||²$
   • $S_{new} = ||r(θ_k + Δθ)||²$
   • $ared = S_{old} - S_{new}$
7. 计算预测减少: $pred = 2 g^T Δθ - Δθ^T A Δθ$
8. 计算模型可信度比 $ρ = ared / pred$
9. 如果 $\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，加快收敛
10. 否则$\rho_k \leq 0$，模型预测失败：
    • 拒绝更新：$θ_{k+1} = θ_k$
    • 增大 $\lambda_k$（例如 $\lambda_{k+1} = \lambda_k \cdot \nu$），更接近梯度下降，更保守

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

1. 在更新 $\lambda$ 时加边界保护：
   • $λ_{k+1} = \text{max}(λ_k / ν, λ_{min})$ # 防止 λ → 0 导致 GN 不稳定
   • $λ_{k+1} = \text{min}(λ_k * ν, λ_{max})$ # 防止 λ → ∞ 导致步长太小
2. 公式 $ρ = ared / pred$ 在 $pred <= 0$ 时会导致误判或除零，此时需要将这次逼近视为不可信，对应同 $ρ <= 0$ 的处理（拒绝步、增大 $λ$）。
3. 推荐初始 $λ$ 取法：(\lambda_0 = \tau \cdot \max(\operatorname{diag}(A_0)))，$τ$ 可取 1e-3 至 1e-1。这样能自动按参数尺度调整初始阻尼。
4. $D_k = \text{diag}(A)$ 通常比 $I$ 更稳健（参数尺度不同会导致不同的步长）。
5. 典型的$\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方法中加入阻尼项：

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

4.1 阻尼项“自适应缩放”

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

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

4.2 阻尼因子动态调整

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

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

同时还引入因子 $\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 成为现代非线性最小二乘的工业标准算法，也让它成为理解信任域思想的入门之路。