LM、GN算法原理及实现
2024-05-15 02:46:34
1、LM算法原理
2、 高斯-牛顿算法
将进行一阶展开为：
。
目标函数可转化为
上式对求导得：
简化为,即为高斯-牛顿算法的更新方程
3、GN与LM的区别
LM会拒绝可能导致更大的残差平方和的更新，速度也可能会慢些，而使用GN必须保证H的可逆性。
4、LM的c++简单实现
#include <bits/stdc++.h>
#include <eigen3/Eigen/Dense>using namespace std;
using namespace Eigen;vector<Vector2d> obs_;
int max_iterations_ = 100;
int max_inner_iterations_ = 10;
int lambda_;
VectorXd errors_;
Matrix3d hessian_;
Vector3d hb_;
MatrixXd jacobian_;
double tau_ = 1e-5;
double epsilon_ = 1e-5;
double current_squared_error;
double a_,b_,c_;void computeJacobian();
void computeHb();
void initLambda();int main(){double a = 1, b = 2 , c= 3;a_ = 0, b_= 0, c_=0;std::default_random_engine generator;std::normal_distribution<double> noise(0,0.1);for(int i=0;i<100;i++){Vector2d ob;ob(0) = i;ob(1) = a*i*i+b*i+c+noise(generator);obs_.push_back(ob);}computeJacobian();computeHb();initLambda();int iter = 0;double delta_norm; //每次优化增量的大小double v = 2.0;while (iter++<max_iterations_){current_squared_error = errors_.squaredNorm();int inner_iter = 0;while (inner_iter++ < max_inner_iterations_){Matrix3d damper = lambda_ * Matrix3d::Identity();Vector3d delta = (hessian_ + damper).inverse()*hb_;double aa = a_ + delta(0);double bb = b_ + delta(1);double cc = c_ + delta(2);delta_norm = delta.norm();double new_squared_error = 0.0;for(int i = 0;i<obs_.size();i++){double xi = obs_[i](0);double yi = obs_[i](1);double yy = aa*xi*xi+bb*xi+cc;new_squared_error += pow(yy-yi,2);}double rou = (current_squared_error-new_squared_error)/ (0.5*delta.transpose()*(lambda_*delta+hb_)+1e-3);if(rou>0){a_= aa,b_ = bb,c_=cc;lambda_ = lambda_*max(1.0/3.0, 1-pow((2*rou-1),3));v = 2;break;}else{lambda_ = lambda_ *v;v *= 2;}}if(delta_norm< epsilon_){break;}computeJacobian();computeHb();}cout<< a_<<" "<<b_<<" "<<c_<<endl;return 0;
}void initLambda(){double max_diagnal = 0.;for(int i=0;i<3;i++){max_diagnal = max(fabs(hessian_(i,i)), max_diagnal);}lambda_ = tau_ * max_diagnal;
}void computeJacobian(){jacobian_.resize(obs_.size(),3);errors_.resize(obs_.size());for(int i=0;i<obs_.size();i++){double xi = obs_[i][0];double yi = obs_[i][1];Matrix<double,1,3> jac;double yy = a_*xi*xi+b_*xi+c_;jac(0,0) = xi*xi;jac(0,1) = xi;jac(0,2) = 1;jacobian_.row(i) = jac;errors_(i) = yy-yi;}
}void computeHb(){hessian_ = jacobian_.transpose()*jacobian_;hb_ = -jacobian_.transpose()*errors_;
}5、Ceres示例
#include<bits/stdc++.h>
#include<eigen3/Eigen/Dense>
#include<ceres/ceres.h>
using namespace std;
using namespace Eigen;struct costFunc
{costFunc(double x, double y):x_(x),y_(y){}template<typename T>bool operator()(const T* const a,const T* const b, const T* const c, T* residual) const{residual[0] = T(y_) - (a[0]*x_*x_+b[0]*x_+c[0]);return true;}const double x_;const double y_;
};int main(){vector<Vector2d> obs;std::default_random_engine generator;std::normal_distribution<double> noise(0,0.1);double a = 1.,b = 2., c= 3.;double a0 = 0,b0= 0,c0 = 0;cout<< "original: " << a0<<" "<<b0<<" "<<c0<<endl;for(auto x=1.0;x<100;x+=1.){double y = a*x*x+b*x+c + noise(generator);obs.push_back(Vector2d(x,y));}ceres::Problem problem;for(int i=0;i<obs.size();i++){ceres::CostFunction* pCostFunction = new ceres::AutoDiffCostFunction<costFunc,1,1,1,1>(new costFunc(obs[i][0],obs[i][1]));problem.AddResidualBlock(pCostFunction, nullptr, &a0,&b0,&c0);}ceres::Solver::Options options;options.linear_solver_type = ceres::DENSE_QR;options.minimizer_progress_to_stdout = true;ceres::Solver::Summary summary;ceres::Solve(options,&problem,&summary);cout<< summary.BriefReport()<<endl;cout<< "optimized: " << a0<<" "<<b0<<" "<<c0<<endl;return 0;
}
解析求导形式：
#include<bits/stdc++.h>
#include<eigen3/Eigen/Dense>
#include<ceres/ceres.h>
using namespace std;
using namespace Eigen;class CostFunc: public ceres::SizedCostFunction<1,3>{
public:CostFunc(double x , double y):x_(x),y_(y){}virtual ~CostFunc(){}virtual bool Evaluate(double const* const* parameters, double* residuals, double** jacobians) const{const double a = parameters[0][0];const double b = parameters[0][1];const double c = parameters[0][2];residuals[0] = a*x_*x_+b*x_+c -y_;if(!jacobians) return true;double* jac = jacobians[0];if(!jac) return true;jac[0] = x_*x_;jac[1] = x_;jac[2] = 1;return true;}private:const double x_;const double y_;
};int main(){vector<Vector2d> obs;std::default_random_engine generator;std::normal_distribution<double> noise(0,0.1);double a = 1.,b = 2., c= 3.;double a0 = 0,b0= 0,c0 = 0;cout<< "original: " << a0<<" "<<b0<<" "<<c0<<endl;for(auto x=1.0;x<100;x+=1.){double y = a*x*x+b*x+c + noise(generator);obs.push_back(Vector2d(x,y));}double param[3] = {a0,b0,c0};ceres::Problem problem;for(int i=0;i<obs.size();i++){ceres::CostFunction* pCostFunction = new CostFunc(obs[i][0],obs[i][1]);problem.AddResidualBlock(pCostFunction, nullptr, param);}ceres::Solver::Options options;options.linear_solver_type = ceres::DENSE_QR;options.minimizer_progress_to_stdout = true;ceres::Solver::Summary summary;ceres::Solve(options,&problem,&summary);cout<< summary.BriefReport()<<endl;cout<< "optimized: " << param[0]<<" "<<param[1]<<" "<<param[2]<<endl;return 0;
}
LM、GN算法原理及实现相关推荐
- LM(Levenberg–Marquardt)算法原理及其python自定义实现
LM算法原理及其python自定义实现 LM(Levenberg–Marquardt)算法原理 LM算法python实现 实现步骤: 代码: 运行结果: LM(Levenberg–Marquardt) ...
- mlp神经网络和bp神经网络,bp神经网络lm算法原理
MATLAB中训练LM算法的BP神经网络 1.初始权值不一样,如果一样,每次训练结果是相同的 2.是 3.在train之前修改权值,IW,LW,b,使之相同 4.取多次实验的均值 一点浅见,仅供参考 ...
- 手把手实战：利用LM神经网络算法自动识别窃电用户（附代码）
来源:数据路 本文约3000字,建议阅读7分钟. 通过本文给大家介绍利用LM神经网络算法进行电力窃漏电用户的自动识别. 背景与挖掘目标 背景 传统的防窃漏电方法主要通过定期巡检.定期校验电表.用户举报 ...
- 加权GN算法的Java实现
本文为加权GN算法的Java实现,具体算法原理请参考前一篇文章GN算法的简介,整个代码可从http://download.csdn.net/detail/u012483487/9447115下载,如有 ...
- 利用GN算法进行社区发现
利用GN算法进行社区发现 算法原理 评价方法 实验结果 一种高效实现 算法原理 这个算法是这样做的.根据边的重要性将边一条一条的删去,那么随着边的删除,那些节点也会慢慢变成独立的子连通图.如下图: 所 ...
- MUSIC算法原理及MATLAB代码 阵列信号处理
MUSIC算法原理及MATLAB代码 阵列信号处理 MUSIC(multiple signal classification algorithm)算法是一种基于矩阵特征空间分解的方法.从几何角度讲,信 ...
- 神经网络算法实例说明，简单神经网络算法原理
神经网络算法实例说明有哪些? 在网络模型与算法研究的基础上,利用人工神经网络组成实际的应用系统,例如,完成某种信号处理或模式识别的功能.构作专家系统.制成机器人.复杂系统控制等等. 纵观当代新兴科学技 ...
- 社区发现算法——GN算法与FN算法
GN算法 本算法的具体内容请参考Finding and evaluating community structure in networks(Newman and Girvan). 重要概念 边介数( ...
- 极限学习机（ELM）算法原理及C++代码实现
从事机器学习研究两年多了,第一个用的算法就是极限学习机,CSDN算是我的领路人,在此感谢一下CSDN上分享知识的大神们,即将毕业,所学颇杂,想在此开始总结一下,顺便也为即将入坑的新手们做一些贡献.第一 ...
最新文章
- Javascript实现网页水印(非图片水印)
- 在RedHat4 64位操作系统下，安装Oracle 10g
- 以下关于java中布局管理说法错误的是_对于 Java 中的布局管理器，以下说法中错误的是（ ）。_2019复习答案_学小易找答案...
- REST Framework 的用户认证组件
- 分布式实物实现方式_这是您完成实物产品设计任务的方式
- golang常用手册：数据类型、变量和常量
- brctl 设置ip_Linux 网桥配置命令：brctl
- Python生成词云
- 看似简单的hashCode和equals面试题，竟然有这么多坑！
- Ensemble查看基因的外显子信息，并根据染色体位点判断是第几号外显子
- 【STL】C++ STL超全总结
- 倡导低碳低成本出行，神州租车用实力说话
- minisforum HX90G/HX99G miniPC-Hackintosh-Opencore 黑苹果efi引导文件
- 基于SaaS的教务系统平台设计构想
- 中国电子学会2022年12月份青少年软件编程Python等级考试试卷四级真题(含答案)
- 新任务管理系统YYSchedule-介绍-引擎执行机制及结果回收机制
- EndNote X9 闪退解决办法（最简版）
- android 8.1 9.0 10.0 默认允许安装第三方app去掉未知来源弹窗直接安装apk
- 关于在线使用remix-ide
- 项目管理（PJM）的定义