01 LM 算法及 Cpp 实现

已于 2023-11-24 10:01:29 修改
阅读量1.1k 收藏 17
点赞数 14
分类专栏： 视觉SLAM十四讲 文章标签： 算法 SLAM
于 2023-11-24 10:01:05 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_44175983/article/details/134591886
版权

文章目录

01 LM 算法及 Cpp 实现

1.1 应用

LM 算法用于解决非线性最小二乘问题

min ⁡ x F ( x ) = 1 2 ∥ f ( x ) ∥ 2 2 (1) \min _x F(x)=\frac{1}{2}\|f(\boldsymbol{x})\|_2^2 \tag{1} xmin​F(x)=21​∥f(x)∥22​(1)

其中， f ( x ) f(\boldsymbol{x}) f(x) 是指各项的残差。

LM 算法有两种推导方式，即 阻尼法 和 置信域 法（见《十四讲》），这里用阻尼法进行推导。

1.2 阻尼法推导

（1）一阶泰勒展开近似

将 f ( x ) f(\boldsymbol{x}) f(x) 在 x n \boldsymbol{x_n} xn​ 处一阶泰勒展开（把 Δ x \Delta \boldsymbol{x} Δx 看做未知数），

f ( x n + Δ x ) ≈ f ( x n ) + J ( x n ) Δ x (2) f(\boldsymbol{x_n}+\Delta \boldsymbol{x}) \approx f(\boldsymbol{x_n})+\boldsymbol{J}\left(\boldsymbol{x}_n\right) \Delta \boldsymbol{x} \tag{2} f(xn​+Δx)≈f(xn​)+J(xn​)Δx(2)

那么问题转化为，对于每次迭代，求解最优的 Δ x \Delta \boldsymbol{x} Δx。

Δ x ∗ = arg ⁡ min ⁡ Δ x 1 2 ∥ f ( x n ) + J n ( x n ) Δ x ∥ 2 (3) \Delta \boldsymbol{x}^*=\arg \min _{\Delta \boldsymbol{x}} \frac{1}{2}\left\|f\left(\boldsymbol{x}_{\boldsymbol{n}}\right)+\boldsymbol{J}_{\boldsymbol{n}}\left(\boldsymbol{x}_{\boldsymbol{n}}\right) \Delta \boldsymbol{x}\right\|^2 \tag{3} Δx∗=argΔxmin​21​∥f(xn​)+Jn​(xn​)Δx∥2(3)

（2）加入阻尼项

Δ x ∗ = arg ⁡ min ⁡ Δ x M ( Δ x ) = 1 2 ∥ f ( x n ) + J n ( x n ) Δ x ∥ 2 + 1 2 μ Δ x T Δ x (4) \Delta \boldsymbol{x}^*=\arg \min _{\Delta \boldsymbol{x}} M(\Delta \boldsymbol{x})=\frac{1}{2}\left\|f\left(\boldsymbol{x}_{\boldsymbol{n}}\right)+\boldsymbol{J}_{\boldsymbol{n}}\left(\boldsymbol{x}_{\boldsymbol{n}}\right) \Delta \boldsymbol{x}\right\|^2+\frac{1}{2} \mu \Delta \boldsymbol{x}^T \Delta \boldsymbol{x} \tag{4} Δx∗=argΔxmin​M(Δx)=21​∥f(xn​)+Jn​(xn​)Δx∥2+21​μΔxTΔx(4)

其中， μ \mu μ 为阻尼系数，阻尼项 1 2 μ Δ x T Δ x \frac{1}{2} \mu \Delta \boldsymbol{x}^T \Delta \boldsymbol{x} 21​μΔxTΔx 可以看做是对于过大的 Δ x \Delta \boldsymbol{x} Δx 的惩罚。

（3）求极值

对 Δ x \Delta \boldsymbol{x} Δx 求导，并令其等于零，

J ( x n ) T f ( x n ) + J ( x n ) T J ( x n ) Δ x + μ Δ x = 0 (5) \boldsymbol{J}\left(\boldsymbol{x}_{\boldsymbol{n}}\right)^T f\left(\boldsymbol{x}_{\boldsymbol{n}}\right)+\boldsymbol{J}\left(\boldsymbol{x}_{\boldsymbol{n}}\right)^T \boldsymbol{J}\left(\boldsymbol{x}_{\boldsymbol{n}}\right) \Delta \boldsymbol{x}+\mu \Delta \boldsymbol{x}=\mathbf{0} \tag{5} J(xn​)Tf(xn​)+J(xn​)TJ(xn​)Δx+μΔx=0(5)

得

Δ x ∗ = − ( J ( x n ) T J ( x n ) + μ I ) − 1 J ( x n ) T f ( x n ) (6) \Delta \boldsymbol{x}^*=-\left(\boldsymbol{J}\left(\boldsymbol{x}_{\boldsymbol{n}}\right)^T \boldsymbol{J}\left(\boldsymbol{x}_{\boldsymbol{n}}\right)+\mu \boldsymbol{I}\right)^{-1} \boldsymbol{J}\left(\boldsymbol{x}_{\boldsymbol{n}}\right)^T f\left(\boldsymbol{x}_{\boldsymbol{n}}\right) \tag{6} Δx∗=−(J(xn​)TJ(xn​)+μI)−1J(xn​)Tf(xn​)(6)

令

J ( x n ) T J ( x n ) = H J ( x n ) T f ( x n ) = g \begin{aligned} \boldsymbol{J}\left(\boldsymbol{x}_{\boldsymbol{n}}\right)^T \boldsymbol{J}\left(\boldsymbol{x}_{\boldsymbol{n}}\right) &=\boldsymbol{H} \\ \boldsymbol{J}\left(\boldsymbol{x}_{\boldsymbol{n}}\right)^T f\left(\boldsymbol{x}_{\boldsymbol{n}}\right) &=\boldsymbol{g} \end{aligned} J(xn​)TJ(xn​)J(xn​)Tf(xn​)​=H=g​

则式（6）可写为

Δ x ∗ = − ( H + μ I ) − 1 g (7) \Delta \boldsymbol{x}^*=-(\boldsymbol{H}+\mu \boldsymbol{I})^{-1} \boldsymbol{g} \tag{7} Δx∗=−(H+μI)−1g(7)

（4）阻尼系数 μ \mu μ 的调节

定义增益系数 ρ \rho ρ

ρ = f ( x + Δ x ) − f ( x ) J ( x ) T Δ x (8) \rho=\frac{f(\boldsymbol{x}+\boldsymbol{\Delta x})-f(\boldsymbol{x})}{\boldsymbol{J}(\boldsymbol{x})^T\Delta \boldsymbol{x}} \tag{8} ρ=J(x)TΔxf(x+Δx)−f(x)​(8)

其中，分子是实际下降的值，分母是近似下降的值。若 ρ \rho ρ 接近于 1 ，则近似效果好；若 ρ \rho ρ 太小，则说明实际减小的值远小于近似减小的值，即近似效果较差，需缩小近似范围，即增大阻尼系数 μ \mu μ； 若 ρ \rho ρ 太大，则说明实际减小的值大于近似减小的值，则需放大近似范围，即减小 μ \mu μ。

两种调节方法：

 if  ρ < 0.25 μ : = μ ∗ 2  else if  ρ > 0.75 μ : = μ / 3 \begin{aligned} \text { if } & \rho<0.25 \\ & \mu:=\mu * 2 \\ \text { else if } & \rho>0.75 \\ & \mu:=\mu / 3 \end{aligned}  if  else if ​ρ<0.25μ:=μ∗2ρ>0.75μ:=μ/3​

也就是， ρ < 0.25 \rho<0.25 ρ<0.25 时增大阻尼系数； ρ > 0.75 \rho>0.75 ρ>0.75 时，减小阻尼系数。

或者：

i f ρ > 0 μ : = μ ∗ max ⁡ { 1 3 , 1 − ( 2 ρ − 1 ) 3 } ; ν : = 2 e l s e μ : = μ ∗ ν ; ν : = 2 ∗ ν \begin{aligned} if \rho>0 \\ &\mu:=\mu * \max \left\{\frac{1}{3}, 1-(2 \rho-1)^3\right\} ; \quad \nu:=2\\ else\\ &\mu:=\mu * \nu ; \quad \nu:=2 * \nu \end{aligned} ifρ>0else​μ:=μ∗max{31​,1−(2ρ−1)3};ν:=2μ:=μ∗ν;ν:=2∗ν​

1.3 Cpp 算法实现

两种方法差别在于 ρ \rho ρ 的分母的计算（我的更快？）。

优化目标是待定系数，把他们看做 未知数，需要计算的 step 就是这些待定系数的 step。

（1）我的方法

/***********************************************************                                          *
* Time: 2022/11/3
* Author: xiaocong
* Function: LM 算法
 *
 * @ 解决的是最小二乘问题，也就是找到最优的系数使得残差最小
 * @ obj = A * sin(Bx) + C * cos(D*x) - F
 *      A * sin(Bx) + C * cos(D*x) 是目标函数
 *      F 是实际值
 *      目标是找到使 obj 最小的系数 A B C D
***********************************************************/

#include <Eigen/Dense>      // 稠密矩阵
#include <Eigen/Sparse>      // 稀疏矩阵
#include <iostream>
#include <cmath>

using namespace std;
using namespace Eigen;

const double DERIV_STEP = 1e-5;
const int MAX_INTER = 100;

#define max(a, b) (((a)>(b))?(a):(b))

// 定义目标函数
double func(const VectorXd input, const VectorXd &output, const VectorXd &params, int objIndex)
{
    // obj = A * sin(Bx) + C * cos(D*x) - F
    double x1 = params(0);   // params 中存储的是系数
    double x2 = params(1);
    double x3 = params(2);
    double x4 = params(3);

    double t = input(objIndex);    // input 中存储的是 x
    double f = output(objIndex);   // output 中存储的是实际值

    return x1 * sin(x2 * t) + x3 * cos(x4 * t) - f;    // 返回 objIndex 项的残差
}

// 计算残差矩阵
VectorXd objF(const VectorXd &input, const VectorXd &output, const VectorXd &params)
{
    VectorXd obj(input.rows());    // 存储所有残差，也就是残差矩阵
    for (int i = 0; i < input.rows(); i++)
        obj(i) = func(input, output, params, i);
    return obj;   // 返回残差矩阵
}


// 残差平方和
double Func(const VectorXd &obj)
{
    return obj.squaredNorm() / 2.0;
}

// 求某个系数在某点的导数
double Deriv(const VectorXd &input, const VectorXd &output, int objIndex, const VectorXd &params, int paraIndex)
{
    VectorXd para1 = params;
    VectorXd para2 = params;

    para1(paraIndex) -= DERIV_STEP;
    para2(paraIndex) += DERIV_STEP;

    double obj1 = func(input, output, para1, objIndex);
    double obj2 = func(input, output, para2, objIndex);

    return (obj2 - obj1) / (2 * DERIV_STEP);     // 该点处的导数，为求雅克比矩阵做准备
}


// 计算雅克比矩阵
/****************************
 * 我们优化的是系数 params，把他们看做未知数，分别求导，得到雅克比矩阵
 * 维度：（input.rows() x params.rows()）
 * [[df/dA df/dB df/dC df/dD]    <--- x1
 *  [df/dA df/dB df/dC df/dD]    <--- x2
 *  [.......................]
 *  [df/dA df/dB df/dC df/dD]]    <--- xn
 ****************************/
MatrixXd Jacobian(const VectorXd &input, const VectorXd &output, const VectorXd &params)
{
    int rowNum = input.rows();
    int colNum = params.rows();

    MatrixXd Jac(rowNum, colNum);

    for (int i = 0; i < rowNum; i++)
        for (int j = 0; j < colNum; j++)
            Jac(i, j) = Deriv(input, output, i, params, j);

    return Jac;
}


//求 Hessian 矩阵对角线最大值
// Hessian 矩阵：二阶导数
double maxMatrixDiagonale(const MatrixXd &Hessian)
{
    int max = 0;
    for (int i = 0; i < Hessian.rows(); i++)
    {
        if(Hessian(i, i) > max)
            max = Hessian(i, i);
    }
    return max;
}


// 近似下降值
double linerDeltaL(const VectorXd &step, const MatrixXd &Jac)
{
    VectorXd L = Jac * step;
    return L.norm();
}


void levenMar(const VectorXd &input, const VectorXd &output, VectorXd &params)
{
    double epsilon = 1e-12;
    int iterCnt = 0;         // 迭代计数
    double tao = 1e-3;
    long long v = 2;

    // 求初始的 u
    MatrixXd Jac = Jacobian(input, output, params);    // 得到雅克比矩阵
    MatrixXd H = Jac.transpose() * Jac;            // 得到 Hessian 矩阵，4x4
    double u = tao * maxMatrixDiagonale(H);     // Hessian 矩阵对角线最大值乘 tao


    while (iterCnt < MAX_INTER)     // double 类型不能直接作比较
    {
        VectorXd obj = objF(input, output, params);    // 误差矩阵
        MatrixXd Jac = Jacobian(input, output, params);    // 得到雅克比矩阵
        MatrixXd H = Jac.transpose() * Jac;            // 得到 Hessian 矩阵，4x4
        VectorXd g = Jac.transpose() * obj;     // 也就是 g，4x1

        VectorXd step = (H + u * MatrixXd::Identity(H.rows(), H.cols())).inverse() * g;
        if(step.norm() < epsilon)
            break;

        cout << "H：" << endl << H << endl;
        cout << "step: " << endl << step << endl;

        VectorXd paramsNew(params.rows());
        paramsNew = params - step;            // 更新 params

        // 计算 params 误差
        obj = objF(input, output, params);

        // 计算 paramsNew 误差
        VectorXd  obj_new = objF(input, output, paramsNew);

        double deltaF = Func(obj) - Func(obj_new);     // 求差
        double deltaL = linerDeltaL(step, Jac);

        // 计算增益系数 rho
        double rho = deltaF / deltaL;          // 实际下降值 / 近似下降值
        cout << "rho is； " << rho <<endl;

        if(rho > 0)
        {
            params = paramsNew;
            u *= max(1.0 / 3.0, 1 - pow(2 * rho - 1, 3));
            v = 2;
        } else
        {
            u = u * v;
            v = v * 2;
        }
        cout << "u= " << u << "\tv= " << v << endl;

        iterCnt ++;
        cout << "Iteration " << iterCnt << " times, result is :" << endl<< params << endl;

    }
}


int main()
{
    int params_num = 4;
    int total_data = 100;

    VectorXd input(total_data);
    VectorXd output(total_data);

    double A = 5, B = 1, C = 10, D = 2;    // 初始化

    // 生成数据
    for (int i = 0; i < total_data; i++)
    {
        double x = 20.0 * ((rand() % 1000) / 1000.0) - 10.0;    // [-10, 10]
        double deltaY = 2.0 * (rand() % 1000) / 1000.0;         // 随机噪声，[0, 2]
        double y = A * sin(B*x) + C * cos(D*x) + deltaY;
        input(i) = x;
        output(i) = y;
    }


    VectorXd params_levenMar(params_num);
    params_levenMar << 3.6, 1.3, 7.2, 1.7;

    levenMar(input, output, params_levenMar);
    cout << "**********************************************" << endl;
    cout << "Levenberg-Marquardt parameter: " << endl << params_levenMar << endl;

}

输出

Levenberg-Marquardt parameter:
 4.85628
0.997904
 10.0524
   2.003

（2）网络方法

/***********************************************************                                          *
* Time: 2022/11/2
* Author: xiaocong
* Function: LM 算法
 *
 * @ 解决的是最小二乘问题，也就是找到最优的系数使得残差最小
 * @ obj = A * sin(Bx) + C * cos(D*x) - F
 *      A * sin(Bx) + C * cos(D*x) 是目标函数
 *      F 是实际值
 *      目标是找到使 obj 最小的系数 A B C D
***********************************************************/

#include <Eigen/Dense>      // 稠密矩阵
#include <Eigen/Sparse>      // 稀疏矩阵
#include <iostream>
#include <iomanip>    // 控制输入输出格式等
#include <cmath>

using namespace std;
using namespace Eigen;

const double DERIV_STEP = 1e-5;
const int MAX_INTER = 100;

#define max(a, b) (((a)>(b))?(a):(b))

// 定义目标函数
double func(const VectorXd input, const VectorXd &output, const VectorXd &params, int objIndex)
{
    // obj = A * sin(Bx) + C * cos(D*x) - F
    double x1 = params(0);   // params 中存储的是系数
    double x2 = params(1);
    double x3 = params(2);
    double x4 = params(3);

    double t = input(objIndex);    // input 中存储的是 x
    double f = output(objIndex);   // output 中存储的是实际值

    return x1 * sin(x2 * t) + x3 * cos(x4 * t) - f;    // 返回 objIndex 项的残差
}

// 计算残差矩阵
VectorXd objF(const VectorXd &input, const VectorXd &output, const VectorXd &params)
{
    VectorXd obj(input.rows());    // 存储所有残差，也就是残差矩阵
    for (int i = 0; i < input.rows(); i++)
        obj(i) = func(input, output, params, i);
    return obj;   // 返回残差矩阵
}


// 残差平方和
double Func(const VectorXd &obj)
{
    return obj.squaredNorm() / 2.0;
}

// 求某个系数在某点的导数
double Deriv(const VectorXd &input, const VectorXd &output, int objIndex, const VectorXd &params, int paraIndex)
{
    VectorXd para1 = params;
    VectorXd para2 = params;

    para1(paraIndex) -= DERIV_STEP;
    para2(paraIndex) += DERIV_STEP;

    double obj1 = func(input, output, para1, objIndex);
    double obj2 = func(input, output, para2, objIndex);

    return (obj2 - obj1) / (2 * DERIV_STEP);     // 该点处的导数，为求雅克比矩阵做准备
}


// 计算雅克比矩阵
/****************************
 * 我们优化的是系数 params，把他们看做未知数，分别求导，得到雅克比矩阵
 * 维度：（input.rows() x output.rows()）
 * [[df/dA df/dB df/dC df/dD]    <--- x1
 *  [df/dA df/dB df/dC df/dD]    <--- x2
 *  [.......................]
 *  [df/dA df/dB df/dC df/dD]]    <--- xn
 ****************************/
MatrixXd Jacobian(const VectorXd &input, const VectorXd &output, const VectorXd &params)
{
    int rowNum = input.rows();
    int colNum = params.rows();

    MatrixXd Jac(rowNum, colNum);

    for (int i = 0; i < rowNum; i++)
        for (int j = 0; j < colNum; j++)
            Jac(i, j) = Deriv(input, output, i, params, j);

    return Jac;
}


//求 Hessian 矩阵对角线最大值
// Hessian 矩阵：二阶导数
double maxMatrixDiagonale(const MatrixXd &Hessian)
{
    int max = 0;
    for (int i = 0; i < Hessian.rows(); i++)
    {
        if(Hessian(i, i) > max)
            max = Hessian(i, i);
    }
    return max;
}


double linerDeltaL(const VectorXd &step, const VectorXd &gradient, const double u)
{
    double L = step.transpose() * (u * step - gradient);
    return L;
}


void levenMar(const VectorXd &input, const VectorXd &output, VectorXd &params)
{
    int errNum = input.rows();
    int paraNum = params.rows();

    // initial parameters
    VectorXd obj = objF(input, output, params);    // 得到误差矩阵
    MatrixXd Jac = Jacobian(input, output, params);    // 得到雅克比矩阵
    MatrixXd H = Jac.transpose() * Jac;            // 得到 Hessian 矩阵，4x4
    VectorXd gradient = Jac.transpose() * obj;     // 也就是 g，4x1

    double tao = 1e-3;
    long  long v = 2;
    double epsilon1 = 1e-12, epsilon2 = 1e-12;
    double u = tao * maxMatrixDiagonale(H);     // Hessian 矩阵对角线最大值乘 tao
    bool found = gradient.norm() <= epsilon1;
    if (found) return;        // 直接退出程序，不再执行后面的程序

    double last_sum = 0;
    int iterCnt = 0;         // 迭代计数

    while (iterCnt < MAX_INTER)
    {
        VectorXd obj = objF(input, output, params);    // 误差矩阵

        MatrixXd Jac = Jacobian(input, output, params);    // 得到雅克比矩阵
        MatrixXd H = Jac.transpose() * Jac;            // 得到 Hessian 矩阵，4x4
        VectorXd gradient = Jac.transpose() * obj;     // 也就是 g，4x1

        if(gradient.norm() < epsilon1)
        {
            cout << "stop g(x) = 0 for a local minimizer optimizer." << endl;
            break;
        }

        cout << "H：" << endl << H << endl;

        VectorXd step = (H + u * MatrixXd::Identity(paraNum, paraNum)).inverse() * gradient;
        // 求 Delta x = (H + uI)^{-1}g      注意：step 维度(4x1)

        cout << "step: " << endl << step << endl;

        if(step.norm() <= epsilon2 * (params.norm()) + epsilon2)
        {
            cout << "stop because change in x is small" << endl;
            break;
        }

        VectorXd paramsNew(params.rows());
        paramsNew = params - step;            // 更新 params

        // 计算 params 误差
        obj = objF(input, output, params);

        // 计算 paramsNew 误差
        VectorXd  obj_new = objF(input, output, paramsNew);

        double deltaF = Func(obj) - Func(obj_new);     // 求差
        double deltaL = linerDeltaL(-1 * step, gradient, u);

        // 计算增益系数 rho
        double rho = deltaF / deltaL;          // 实际下降值 / 近似下降值
        cout << "rho is； " << rho <<endl;

        if(rho > 0)
        {
            params = paramsNew;
            u *= max(1.0 / 3.0, 1 - pow(2 * rho - 1, 3));
            v = 2;
        } else
        {
            u = u * v;
            v = v * 2;
        }
        cout << "u= " << u << "\tv= " << v << endl;

        iterCnt ++;
        cout << "Iteration " << iterCnt << " times, result is :" << endl<< params << endl;


    }
}


int main()
{
    int params_num = 4;
    int total_data = 100;

    VectorXd input(total_data);
    VectorXd output(total_data);

    double A = 5, B = 1, C = 10, D = 2;    // 初始化

    // 生成数据
    for (int i = 0; i < total_data; i++)
    {
        double x = 20.0 * ((rand() % 1000) / 1000.0) - 10.0;    // [-10, 10]
        double deltaY = 2.0 * (rand() % 1000) / 1000.0;         // 随机噪声，[0, 2]
        double y = A * sin(B*x) + C * cos(D*x) + deltaY;
        input(i) = x;
        output(i) = y;
    }


    VectorXd params_levenMar(params_num);
    params_levenMar << 3.6, 1.3, 7.2, 1.7;

    levenMar(input, output, params_levenMar);
    cout << "Levenberg-Marquardt parameter: " << endl << params_levenMar << endl << endl << endl;
    cout << "**********************************************" << endl;

}
算法导航
关注
  • 14
    点赞
  • 踩
  • 17
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
LM算法推导
qq_43526137的博客
09-08 286
最小二乘
LM算法的C语言源代码
09-15
非线性优化中的经典算法，用C语言实现。内附Matlab源代码和相应的参考文献。
matlablm算法代码-lsnn:最小二乘神经网络算法
05-26
matlab lm算法代码lsnn（最小二乘神经网络） 语言：MATLAB 该项目用于非迭代算法开发，该算法主要基于最小二乘法。 以下所有这些算法均基于单隐藏层前馈神经网络（SLFN）结构。 算法列表： ELM（极限学习机），PCA-ELM（主要成分分析），I-ELM（增量），EI-ELM（增强增量），DP-ELM（破坏性简约）和CP-ELM（建设性简约）。 AIL（肛门增量学习） LSM（本地Sigmoid方法） BP（反向传播，LM）matlab包装器。 如何进行所有算法的比较 回归案例：运行runcvreg.m 。 分类案例：运行runcvclass.m 。 要查找超参数：运行hpsearching.m ，您可以将大小写更改为回归或分类。 要总结指标（准确性，＃nodes，时间）：运行`perfsummay.m，您需要手动更改指标-请参阅/阅读代码。
[原创]LM算法的实现
热门推荐
编码无悔
08-09 2万+
完整文章请查看这里。转载请注明出处：本文来自learnhard的博客：http://www.codelast.com/ & http://blog.csdn.net/learnhard/，并保持文章的完整性。 LM算法可用于解决非线性最小二乘问题。多用于曲线拟合等场合。LM算法的实现并不难，这里不讨论使用MATLAB等工具直接得到结果的过程，使用那些工具对于算法编程能力的提高无任何益处。 LM算法
手写系列之手写LM(Levenberg–Marquardt)算法(基于eigen)
卫少东个人博客
07-08 2391
紧接上次的手写高斯牛顿,今天顺便将LM算法也进行了手写实现,并且自己基于eigen的高斯牛顿进行了对比,可以很明显的看到,LM的算法收敛更快,精度也略胜一筹,这次高博的书不够用了,参考网上伪代码进行实现.同学们有什么问题都可以讨论,理论部分的细节后面会补上,最近确实有点忙. 伪代码 程序源码 #include <eigen3/Eigen/Core> #include <eigen3/Eigen/Dense> #include"opencv2/opencv.hpp" #include
都是被逼的... ，LM算法的具体实现python和C++代码
分享型博主
06-14 2370
L-M方法全称Levenberg-Marquardt方法，是一种非线性最小二乘优化算法，它通过同时利用高斯牛顿方法和梯度下降方法来解决非线性最小二乘问题。其核心思想是在每次迭代中，根据当前参数估计计算目标函数的梯度和海森矩阵，并使用这些信息来更新模型参数。LM算法适用于解决各种非线性最小二乘问题，例如数据拟合、无约束非线性优化等。LM算法相对于其他算法的优势在于其自适应调整步长的能力，使得模型更容易收敛到最优解，并且可以避免梯度爆炸或消失的问题。
Levenberg-Marquardt (LM) 算法进行非线性拟合
最新发布
Jurio的博客
04-24 1749
LM算法是一种非线性最小二乘优化算法，用于求解非线性最小化问题。
LM算法详解
penuel
12-11 1万+
1. 高斯牛顿法 残差函数f(x)为非线性函数，对其一阶泰勒近似有： 这里的J是残差函数f的雅可比矩阵，带入损失函数的： 令其一阶导等于0，得： 这就是论文里常看到的normal equation。 2.LM LM是对高斯牛顿法进行了改进，在求解过程中引入了阻尼因子： 2.1 阻尼因子的作用： 2.2 阻尼因子的初始值选取： 一个简单的策略就是： 2.3 阻尼因子的更新策略 3....
从Gauss-Newton算法到 LM算法 （详细推导及MATLAB实现、多自变量问题）
强强的博客
12-18 1万+
Gauss-Newton算法 MATLAB实现结果回顾算法实现总结 结果回顾 Gauss-Newton算法对Gauss-newton算法做了详细的解释，并且使用C++做了实例程序。 但是程序其实有微小错误，实际的坐标并不是年代1815—1885,而是1—8，否则 p=A∗exp(B∗t)p = A*exp(B*t)p=A∗exp(B∗t)拟合时将会迅速增大，也得不到 A=0.7A=0.7A=0.7...
牛顿法，高斯牛顿法，列文伯格-马夸尔特(LM)法 (含代码)【最小二乘非线性求解】
weixin_43763292的博客
11-29 5943
牛顿法是一种函数逼近法,它的基本思想是：在极小点附近用x(k)x^{(k)}x(k)点的二阶泰勒多项式来近似目标函数f(x)f(x)f(x),并用选代点x(k)x^{(k)}x(k)处指向近似二次函数的极小点方向作为搜索方向p(k)p^{(k)}p(k)。minf(x),x∈Rnminf(x),x∈Rn其中f(x)f(x)f(x)在点x(k)x^{(k)}x(k。
LM算法的实现
05-27
Levenberg_Marquardt算法，可进一步优化，得到更优的结果，仅作参考
matlab编写的LM优化算法
05-23
使用matlab编写的LM优化算法，基于imm3215，不是matlab内部函数。FK.m是原函数，JFK是一阶Jaccobi导数
lm算法的matlab实现
12-14
该算法是自己编的一个lm算法，是基于MATLAB的程序
LM算法+PPT讲解
06-04
MATLAB中LM算法的实现+PPT讲解，希望对大家有帮助！LM算法非线性最小二乘法
LM.zip_LM_lm4798_lm算法_opencv lm算法_opencv实现lm
09-23
1. `LM.cpp`或`LM.hpp`: OpenCV的C++实现，包含LM算法的函数声明和实现。 2. `LM.m`: Matlab版本的LM算法实现，可能包含非线性模型的定义、残差计算和雅可比矩阵的代码。 3. 示例数据和测试用例：用于验证算法正确性...
CPP-LM-BP.rar_BP LM_LM BP_LM神经网络_bp神经网络 lm_优化算法BP
07-15
利用LM训练算法来优化BP神经网络算法。
c++实现LM和NTLM哈希
10-31
通过以上步骤，你可以在C++中实现LM和NTLM哈希算法，并在VS2008环境中编译运行。理解这些概念和技术对于网络安全和身份验证领域的开发者来说是非常重要的，同时也能帮助你更好地理解密码学的基本原理。
关于基于VS的Levenberg-Marquardt算法参考
10-08
- `LMAlgorithm.cpp/h`：LM算法的核心实现。 - `NonlinearFunction.cpp/h`：定义非线性函数的类，包括函数值计算和雅可比矩阵计算。 - `main.cpp`：主程序，调用LM算法并输入初始参数。 - `Data.h`：可能包含数据...
NTLM口令加密算法C++实现
05-23
NT（Windows NT LAN Manager）是一种口令加密算法，它用于计算Windows操作系统中的用户口令哈希值。以下是NTLM口令加密算法的C++实现： ```cpp #include <iostream> #include <cstring> #include <Windows.h> using namespace std; void LMHash(const char* pwd, unsigned char* lmhash) { unsigned char* p = lmhash; char password[14] = {0}; strncpy(password, pwd, 14); CharUpperA(password); memset(p, 0, 16); for (int i = 0; i < 14; i += 7) { for (int j = 0; j < 7; j++) { if (password[i + j]) { unsigned char c = password[i + j]; for (int k = 0; k < 8; k++) { unsigned char bit = (c >> k) & 1; *p++ = bit ? 0xff : 0; } } else { *p++ = 0; } } } } void NTLMHash(const char* pwd, unsigned char* ntlmhash) { char password[256]; strncpy(password, pwd, 256); int len = strlen(password); memset(&password[len], 0, 256 - len); HCRYPTPROV hProv; HCRYPTHASH hHash; DWORD dwLen = 16; BYTE rgbHash[16]; if (CryptAcquireContext(&hProv, NULL, NULL, PROV_RSA_FULL, CRYPT_VERIFYCONTEXT)) { if (CryptCreateHash(hProv, CALG_MD4, 0, 0, &hHash)) { if (CryptHashData(hHash, (BYTE*)password, strlen(password), 0)) { if (CryptGetHashParam(hHash, HP_HASHVAL, rgbHash, &dwLen, 0)) { memcpy(ntlmhash, rgbHash, 16); } } CryptDestroyHash(hHash); } CryptReleaseContext(hProv, 0); } } int main() { const char* pwd = "password"; unsigned char lmhash[16]; unsigned char ntlmhash[16]; LMHash(pwd, lmhash); NTLMHash(pwd, ntlmhash); cout << "LM hash: "; for (int i = 0; i < 16; i++) { printf("%02X", lmhash[i]); } cout << endl; cout << "NTLM hash: "; for (int i = 0; i < 16; i++) { printf("%02X", ntlmhash[i]); } cout << endl; return 0; } ``` 该代码使用Windows API中的加密函数计算口令的LM哈希值和NTLM哈希值。在计算LM哈希值时，密码将被转换为大写，并且长度必须小于14个字符。在计算NTLM哈希值时，使用MD4哈希算法计算密码的哈希值，并将结果拷贝到NTLM哈希值输出参数中。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值