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

charlee44

发布于 2025-11-05 08:36:31

1370

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 可能震荡甚至发散。

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

是对角正定矩阵，常见选择为单位矩阵

或

\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：梯度下降法》里的雅可比矩阵是对残差向量

求偏导，而这里是雅可比矩阵式是对模型函数

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

选择（

或

\text{diag}(J_k^T J_k)

）

收敛阈值

\epsilon

进行迭代逼近，对

k = 0, 1, 2, \dots

计算残差

r_k = y - f(x; θ_k)

,维度为

计算雅可比

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

（

是

p×p

，

是

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

减小

\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)

通常比

更稳健（参数尺度不同会导致不同的步长）。

典型的

\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 观察到，单纯使用

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

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

这样能让每个参数方向按自身曲率大小进行调节，避免大尺度参数步长过大、小尺度参数步长过小的问题。也就是前面算法中“

D_k = \text{diag}(A)

比

更稳健”的由来。

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\)