优化算法学习(LM算法)

优化算法学习(LM算法)LM算法可以理解为**Gauss-Newton算法与最速下降法的结合**

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE稳定放心使用

推荐书籍

建议学习,METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS:
http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
篇幅不长,容易理解
学习的时候可以参考另一篇,UNCONSTRAINED OPTIMIZATION:http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3217/pdf/imm3217.pdf

Numerical Optimization 2nd –Jorge Nocedal Stephen J. Wright:
http://www.bioinfo.org.cn/~wangchao/maa/Numerical_Optimization.pdf
《视觉SLAM十四讲》第六讲 https://github.com/gaoxiang12/slambook

理论理解

知乎上看到一个回答非常好:

LM算法可以理解为Gauss-Newton算法与最速下降法的结合,如果理解了如何用上述算法求解目标函数最小值的问题,自然也能理解LM。

其实算法的本质就是 a. 站在当前位置( x k x_k xk ),我们需要一个预言(oracle)告诉我们往哪走能找到目的地(最优解可能的方向,比如梯度方向);b. 我们沿着该方向走了一段距离之后(stepsize),更新当前位置信息( x k + 1 x_{k+1} xk+1 ),再问预言家我们下一步往哪走,以此反复。

所以,梯度下降法,给的 oracle 就是当前位置的梯度信息(损失方程关于变量的一阶导数):

x k + 1 = x k − α g k x_{k+1}=x_k-\alpha g_k xk+1=xkαgk

如果是牛顿法,给的 oracle 就是Hessian matrix(损失方程关于变量的二阶导数):

x k + 1 = x k − H k − 1 g k x_{k+1}=x_k-H_k^{-1}g_k xk+1=xkHk1gk(1)

为什么是一阶导数和二阶导数?因为我们知道,对于任意(处处可导的)方程,在其任意一点,我们都可以用泰勒展开式对其拟合,阶数越高,精度越高。但是,考虑到高阶导数的计算复杂度,以及三阶以上函数的非凸性,也不会使用高阶导数。

好了,那么LM算法的优势是什么?牛顿法虽然收敛速度快,但是需要计算 Hessian matrix,对于高维的问题,计算二阶导数会很复杂。因此我们有了Gauss-Newton算法。Gauss-Newton算法不直接计算Hessian matrix,而是通过 Jacobian matrix 对 Hessian matrix 进行拟合:

H ≈ J T J H\approx J^TJ HJTJ

但是,用 Jacobian matrix 拟合Hessian matrix,所计算出来的结果不一定可逆。所以在此基础上,我们引入了一个identity matrix:

H ≈ J T J + μ I H\approx J^TJ+\mu I HJTJ+μI

这也就得到了LM算法。如果我们把上述式子带入之前的公式(1),可以得到

x k + 1 = x k − ( J k T J k + μ I ) − 1 g k x_{k+1}=x_k-(J_k^TJ_k+\mu I)^{-1}g_k xk+1=xk(JkTJk+μI)1gk

所以我们发现,当 μ \mu μ接近于0时,这个算法近似于Gauss-Newton算法;当 μ \mu μ很大时,这个算法近似于最速下降法。因此,这也是为什么LM算法称为Gauss-Newton算法与最速下降法的结合。最后,上一张图表示几种算法之间的关系:
在这里插入图片描述
参考文献:Wilamowski, B. M., & Yu, H. (2010). Improved computation for Levenberg–Marquardt training. IEEE transactions on neural networks, 21(6), 930-937.

一个回答:Matlab 的话现成的代码也是很多的;比如,Solve nonlinear least-squares (nonlinear data-fitting) problems,或者 Levenberg-Marquardt-Fletcher algorithm for nonlinear least squares problems。你可以在网站里面搜搜有没有适合你的。

作者:Sixiang
链接:https://www.zhihu.com/question/269579938/answer/349205519
来源:知乎

这是最上面推荐的书,英文不难:
在这里插入图片描述

gradient matrix, hessian matrix, jacobian matrix:gradient hessian jacobian matrix
https://www.value-at-risk.net/functions/

csdn 有个不错的博客:数值优化(Numerical Optimization)学习系列-目录:https://blog.csdn.net/fangqingan_java/article/details/48951191

程序实现

视觉SLAM十四讲里推荐了**Ceres库**,Ceres solver 是谷歌开发的一款用于非线性优化的库,在谷歌的开源激光雷达slam项目cartographer中被大量使用。
安装和使用参考:
https://zhaoxuhui.top/blog/2018/04/04/ceres&ls.html
下面把关键操作贴出来:

ceres安装

  • 下载源码
    git clone https://github.com/ceres-solver/ceres-solver.git
  • 安装依赖:
sudo apt-get install cmake
# google-glog + gflags
sudo apt-get install libgoogle-glog-dev libgtest-dev libgflags-dev
# BLAS & LAPACK
sudo apt-get install libatlas-base-dev liblapack-dev
# Eigen3
sudo apt-get install libeigen3-dev
# SuiteSparse and CXSparse (optional)
sudo apt-get install libsuitesparse-dev libcxsparse3.1.4

这里libcxsparse可能存在版本问题(出现找不到对应版本),解决办法:

sudo apt-get install bash-completion
sudo gedit /etc/bash.bashrc

将这一部分取消注释,并保存,即可自动补全:
在这里插入图片描述
sudo apt-get install libsuitesparse-dev libcxsparse(按tab)

cd ceres-solver/
mkdir build
cd build
cmake ..
make -j3
sudo make install

代码:

利用Ceres简单实现最小二乘曲线拟合。首先需要生成数据,这里采用OpenCV的随机数生成器生成误差。

#include <iostream>
#include <opencv2/core/core.hpp>
#include <ceres/ceres.h>

using namespace std;
using namespace cv;
using namespace ceres;

//vector,用于存放x、y的观测数据
//待估计函数为y=3.5x^3+1.6x^2+0.3x+7.8
vector<double> xs;
vector<double> ys;

//定义CostFunctor结构体用于描述代价函数
struct CostFunctor{ 
   
  
  double x_guan,y_guan;
  
  //构造函数,用已知的x、y数据对其赋值
  CostFunctor(double x,double y)
  { 
   
    x_guan = x;
    y_guan = y;
  }
  
  //重载括号运算符,两个参数分别是估计的参数和由该参数计算得到的残差
  //注意这里的const,一个都不能省略,否则就会报错
  template <typename T>
  bool operator()(const T* const params,T* residual)const
  { 
   
    residual[0]=y_guan-(params[0]*x_guan*x_guan*x_guan+params[1]*x_guan*x_guan+params[2]*x_guan+params[3]);
    return true;  
  }
};

//生成实验数据
void generateData()
{ 
   
  RNG rng;
  //RNG::gaussian( σ) 返回一个均值为0,标准差为σ的随机数。
  double w_sigma = 1.0;
  for(int i=0;i<100;i++)
  { 
   
    double x = i;
    double y = 3.5*x*x*x+1.6*x*x+0.3*x+7.8;
    xs.push_back(x);
    ys.push_back(y+rng.gaussian(w_sigma));
  }
  for(int i=0;i<xs.size();i++)
  { 
   
    cout<<"x:"<<xs[i]<<" y:"<<ys[i]<<endl;
  }
}

//简单描述我们优化的目的就是为了使我们估计参数算出的y'和实际观测的y的差值之和最小
//所以代价函数(CostFunction)就是y'-y,其对应每一组观测值与估计值的残差。
//由于我们优化的是残差之和,因此需要把代价函数全部加起来,使这个函数最小,而不是单独的使某一个残差最小
//默认情况下,我们认为各组的残差是等权的,也就是核函数系数为1。
//但有时可能会出现粗差等情况,有可能不等权,但这里不考虑。
//这个求和以后的函数便是我们优化的目标函数
//通过不断调整我们的参数值,使这个目标函数最终达到最小,即认为优化完成
int main(int argc, char **argv) { 
   
  
  generateData();
  
  //创建一个长度为4的double数组用于存放参数
  double params[4]={ 
   1.0};

  //第一步,创建Problem对象,并对每一组观测数据添加ResidualBlock
  //由于每一组观测点都会得到一个残差,而我们的目的是最小化所有残差的和
  //所以采用for循环依次把每个残差都添加进来
  Problem problem;
  for(int i=0;i<xs.size();i++)
  { 
   
    //利用我们之前写的结构体、仿函数,创建代价函数对象,注意初始化的方式
    //尖括号中的参数分别为误差类型,输出维度(因变量个数),输入维度(待估计参数的个数)
    CostFunction* cost_function = new AutoDiffCostFunction<CostFunctor,1,4>(new CostFunctor(xs[i],ys[i]));
    //三个参数分别为代价函数、核函数和待估参数
    problem.AddResidualBlock(cost_function,NULL,params);
  }
  
  //第二步,配置Solver
  Solver::Options options;
  //配置增量方程的解法
  options.linear_solver_type=ceres::DENSE_QR;
  //是否输出到cout
  options.minimizer_progress_to_stdout=true;
  
  //第三步,创建Summary对象用于输出迭代结果
  Solver::Summary summary;
  
  //第四步,执行求解
  Solve(options,&problem,&summary);
  
  //第五步,输出求解结果
  cout<<summary.BriefReport()<<endl;
  
  cout<<"p0:"<<params[0]<<endl;
  cout<<"p1:"<<params[1]<<endl;
  cout<<"p2:"<<params[2]<<endl;
  cout<<"p3:"<<params[3]<<endl;
  return 0;
}

CMakeLists.txt:

cmake_minimum_required(VERSION 2.6)
project(ceres_test)

set( CMAKE_CXX_FLAGS "-std=c++11 -O3" )

# 添加cmake模块以使用ceres库
list( APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake_modules )

# 寻找Ceres库并添加它的头文件
find_package( Ceres REQUIRED )
include_directories( ${CERES_INCLUDE_DIRS} )

# OpenCV
find_package( OpenCV REQUIRED )
include_directories( ${OpenCV_DIRS} )

add_executable(ceres_test main.cpp)

# 与Ceres和OpenCV链接
target_link_libraries( ceres_test ${CERES_LIBRARIES} ${OpenCV_LIBS} )

install(TARGETS ceres_test RUNTIME DESTINATION bin)

另外,有个levmar的C/C++的库:(这个还不会用)
levmar : Levenberg-Marquardt nonlinear least squares algorithms in C/C++
http://users.ics.forth.gr/~lourakis/levmar/index.html#download
http://users.ics.forth.gr/~lourakis/sparseLM/

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/187088.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)


相关推荐

  • 邮件群发怎么设置_qq邮箱邮件怎么群发

    邮件群发怎么设置_qq邮箱邮件怎么群发大家都说30岁的女人一枝花,但是我就狠狠的被家里人催婚了。老妈让我去参加一个相亲,虽然心里不想去,但是为了让老妈开心,还是去参加了这场相亲局。当相亲那天来临时,我进入了跟人约好的咖啡馆,见面时寒暄了几句,就进入了无声的沉默,后来我们聊天时提起了我们的职业,我说我是外企HR,他跟我说他是会展公司的市场部部员,然后我问他工作具体是干什么的,然后他跟我说他是具体用邮件来开发客户,介绍会展公司承包的展览这种工作。我说好厉害的样子,你们是不是也需要邮件群发啊?最近我的邮箱有限,不是特别好用,刚好想换一个邮箱,你平常使

    2022年10月30日
  • rpm卸载多个有依赖的rpm包[通俗易懂]

    rpm卸载多个有依赖的rpm包[通俗易懂][root@devOOo_3.1.0_src]#rpm-qlibxml2[root@devOOo_3.1.0_src]#rpm-qalibxml2*[root@dev~]#rpm-qa|greplibxml2libxml2-python-2.6.26-2.1.12libxml2-devel-2.6.26-2.1.12libxml2-2.6.26libxm…

  • c++实现远程开关机「建议收藏」

    c++实现远程开关机「建议收藏」把远程开、关机写成了一个类的两个静态函数。这两个功能的实现都需要事先对目标主机进行一些设置。其中远程开机需要目标主机主板支持,并且插上网线。部分主机的设置已经写明。另可实现方法:https://www.cnblogs.com/findumars/p/6009474.html原理参考:https://blog.csdn.net/smstong/article/details/16879347…

  • generic host process for win32_hostunreachable怎么解决

    generic host process for win32_hostunreachable怎么解决本人使用的windows2003sp1英文版。昨天开始总是莫明其妙出现GenericHostProcess进程出错提示框,紧跟着svchost内存出错提示框, 之后一些service就停止工作,比如WindowsAudio,必须手动重启才能听音乐;网络连接假死,tcp连接需要重连,不胜其扰。上网搜了些文章,基本上分为三个原因:1,木马病毒。2,系统漏洞。3,硬件驱动问题

    2022年10月10日
  • Invalidate介绍[通俗易懂]

    Invalidate介绍[通俗易懂]1、Invalidate介绍  voidInvalidate(BOOLbErase=TRUE);  该函数的作用是使整个窗口客户区无效。窗口的客户区无效意味着需要重绘,例如,如果一个被其它窗口遮住的窗口变成了前台窗口,那么原来被遮住的部分就是无效的,需要重绘。这时Windows会在应用程序的消息队列中放置WM_PAINT消息。MFC为窗口类提供了WM_PAINT的消息处理函数OnPaint,OnPaint负责重绘窗口。视图类有一些例外,在视图类的OnPaint函数中调用了OnDraw函数,实际

  • 机器学习影响现代云计算的五种方式

    机器学习影响现代云计算的五种方式

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号