Levenberg-Marquardt算法（LMA）是一种优化算法，常用于最小化非线性最小二乘问题。该算法结合了高斯-牛顿法和梯度下降法的优点，因此在处理非线性最小二乘问题时，既能够快速收敛，也能避免高斯-牛顿法在某些情况下不稳定的情况。LMA特别适用于拟合问题，例如曲线拟合和非线性回归。

Levenberg-Marquardt算法概述

给定一个优化问题，我们的目标是最小化以下目标函数（残差平方和）：

其中：

• xi​ 是输入数据，yi​ 是目标值。
• f(xi​,p) 是拟合模型，通常是非线性的。
• p 是模型的参数，目标是通过最小化 S(p) 来求解。

Levenberg-Marquardt算法通过迭代的方式更新参数 p，每次更新使用以下公式：

其中：

• J 是雅可比矩阵。
• r 是残差向量，即 r=y−f(x,p)。
• λ 是一个正则化参数，控制算法是更倾向于高斯-牛顿法还是梯度下降法。
• I 是单位矩阵。

算法步骤

1. 初始化：设定初始参数值 p0，选择初始正则化参数 λ。
2. 计算雅可比矩阵：计算残差向量和雅可比矩阵。
3. 求解更新步骤：根据Levenberg-Marquardt公式更新参数。
4. 调整正则化参数 λ\lambdaλ：如果更新步骤使得目标函数减小，则减小 λ；否则，增大 λ。
5. 收敛判断：如果参数更新的幅度足够小或目标函数变化足够小，则认为算法收敛，停止迭代。

C++实现Levenberg-Marquardt算法

以下是一个简单的C++实现，其中使用Levenberg-Marquardt算法拟合一个非线性模型 y = a e^{bx} + c。

#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

// 定义拟合函数：y = a * exp(bx) + c
double model(double x, const VectorXd& p) {
    return p(0) * exp(p(1) * x) + p(2);
}

// 计算残差：r = y - f(x, p)
VectorXd calculateResiduals(const vector<double>& x, const vector<double>& y, const VectorXd& p) {
    int m = x.size();
    VectorXd residuals(m);
    for (int i = 0; i < m; ++i) {
        residuals(i) = y[i] - model(x[i], p);
    }
    return residuals;
}

// 计算雅可比矩阵：J = [df/dp_0, df/dp_1, ..., df/dp_n]
MatrixXd calculateJacobian(const vector<double>& x, const VectorXd& p) {
    int m = x.size();
    MatrixXd J(m, p.size());
    for (int i = 0; i < m; ++i) {
        double xi = x[i];
        double exp_term = exp(p(1) * xi);
        J(i, 0) = exp_term;       // df/dp_0
        J(i, 1) = p(0) * xi * exp_term;  // df/dp_1
        J(i, 2) = -1;             // df/dp_2
    }
    return J;
}

// Levenberg-Marquardt算法
VectorXd levenbergMarquardt(const vector<double>& x, const vector<double>& y, VectorXd p_init, double lambda_init = 0.01, int max_iter = 100, double tol = 1e-6) {
    VectorXd p = p_init;
    double lambda = lambda_init;
    
    // 迭代过程
    for (int iter = 0; iter < max_iter; ++iter) {
        VectorXd residuals = calculateResiduals(x, y, p);  // 计算残差
        MatrixXd J = calculateJacobian(x, p);               // 计算雅可比矩阵

        // 计算正规方程 (J^T * J + lambda * I) * delta_p = J^T * residuals
        MatrixXd JTJ = J.transpose() * J;
        MatrixXd lambdaI = lambda * MatrixXd::Identity(p.size(), p.size());
        MatrixXd A = JTJ + lambdaI;
        VectorXd delta_p = A.ldlt().solve(J.transpose() * residuals);

        // 更新参数
        VectorXd p_new = p - delta_p;

        // 判断是否收敛
        if ((p_new - p).norm() < tol) {
            cout << "收敛, 迭代次数: " << iter << endl;
            break;
        }

        // 判断目标函数是否减小，如果减小，减小lambda，否则增大lambda
        double cost_old = residuals.squaredNorm();
        residuals = calculateResiduals(x, y, p_new);
        double cost_new = residuals.squaredNorm();

        if (cost_new < cost_old) {
            p = p_new;
            lambda /= 10;  // 减小lambda
        } else {
            lambda *= 10;  // 增大lambda
        }
    }

    return p;
}

int main() {
    // 数据点（假设这是实验数据）
    vector<double> x = {0, 1, 2, 3, 4, 5};
    vector<double> y = {2.5, 3.1, 4.8, 7.3, 11.0, 15.0};

    // 初始猜测参数 (a, b, c)
    VectorXd p_init(3);
    p_init << 1.0, 0.5, 1.0;  // 初始猜测

    // 使用Levenberg-Marquardt算法进行拟合
    VectorXd p_opt = levenbergMarquardt(x, y, p_init);

    // 输出结果
    cout << "优化后的参数: a = " << p_opt(0) << ", b = " << p_opt(1) << ", c = " << p_opt(2) << endl;

    // 计算拟合值并输出
    for (size_t i = 0; i < x.size(); ++i) {
        cout << "x = " << x[i] << ", y = " << y[i] << ", 拟合值 = " << model(x[i], p_opt) << endl;
    }

    return 0;
}

代码解析

1. 拟合模型：假设我们拟合的函数是 y = a e^{bx} + c。model 函数用于计算给定参数 p=[a,b,c] 和输入数据 x 对应的输出 y。

2. 计算残差：calculateResiduals 函数计算每个数据点的残差，即目标值 y_i 与模型输出 f(x_i, p) 之间的差异。

3. 计算雅可比矩阵：calculateJacobian 函数计算雅可比矩阵。对于每个参数，雅可比矩阵的元素是模型输出对该参数的偏导数。

4. Levenberg-Marquardt优化：levenbergMarquardt 函数实现了Levenberg-Marquardt算法。在每次迭代中，它计算残差、雅可比矩阵，并利用正规方程更新参数 p。

5. 收敛判断：如果参数更新幅度足够小（即 小于阈值），则认为算法收敛。

6. 调整正则化参数：如果当前的参数更新能减小目标函数，则减小 λ；否则，增大 λ。

运行示例

假设数据点为：

x = {0, 1, 2, 3, 4, 5}
y = {2.5, 3.1, 4.8, 7.3, 11.0, 15.0}

输出结果可能为：

优化后的参数: a = 1.00488, b = 0.568993, c = 1.03269
x = 0, y = 2.5, 拟合值 = 2.54419
x = 1, y = 3.1, 拟合值 = 3.14998
x = 2, y = 4.8, 拟合值 = 4.80998
x = 3, y = 7.3, 拟合值 = 7.39988
x = 4, y = 11, 拟合值 = 11.0201
x = 5, y = 15, 拟合值 = 14.9962

总结

Levenberg-Marquardt算法是一种高效的优化方法，适用于非线性最小二乘问题。它结合了高斯-牛顿法和梯度下降法的优点，在实际应用中能够快速收敛并避免不稳定的情况。通过合理调整正则化参数 λ\lambdaλ，算法能够在不同问题中保持稳定性。