1 引言

對于非線性最小二乘問題的求解來說，除了Gauss-Newton方法（以下簡稱GN方法）和梯度下降法，另外一種更加實用的求解算法就是Levenberg-Marquardt方法（以下簡稱LM方法）了。LM方法綜合了GN方法和梯度下降法的特性，在實踐中表現出極強的魯棒性和收斂性。在閱讀本文之前，至少需要閱讀以下三篇前置文章：

1. 《最小二乘問題詳解4：非線性最小二乘》
2. 《最小二乘問題詳解6：梯度下降法》
3. 《最小二乘問題詳解7：正則化最小二乘》

2 求解

2.1 基本原理

先復習一下GN方法的關鍵點，也就是求解線性最小二乘子問題：

正則方程是：

其解為：

這里可以看到，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||2\)
   • \(S_{new} = ||r(θ_k + Δθ)||2\)
   • \(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\)，使其調整幅度更柔和。常見做法如下：

這種分級調整在實踐中能顯著提高收斂穩定性，尤其是在存在噪聲或殘差面高度非線性的情況下。

5 總結

Levenberg–Marquardt 方法的最大魅力在于：它不是在梯度下降和高斯牛頓之間取折中，而是根據模型的“可信度”在兩者之間智能切換。這使得 LM 成為現代非線性最小二乘的工業標準算法，也讓它成為理解信任域思想的入門之路。