Levenberg–Marquardt算法学习
本次是对Levenberg–Marquardt的学习总结，是为之后看懂sparse bundle ajdustment打基础。这篇笔记包含如下内容：
回顾⾼斯⽜顿算法，引⼊LM算法
惩罚因⼦的计算(迭代步⼦的计算)
完整的算法流程及代码样例
1.      回顾⾼斯⽜顿，引⼊LM算法
根据之前的博⽂：
假设我们研究如下形式的⾮线性最⼩⼆乘问题：
r(x)为某个问题的残差residual，是关于x的⾮线性函数。我们知道⾼斯⽜顿法的迭代公式：
Levenberg–Marquardt算法是对⾼斯⽜顿的改进，在迭代步长上略有不同：
最速下降法对初始点没有特别要求具有整体收敛性，但是相邻两次的搜索⽅向是相互垂直的，所以收敛并不⼀定快。总⽽⾔之就是：当⽬标函数的等值线接近于圆(球)时，下降较快；等值线类似于扁长的椭球时，⼀开始快，后来很慢。This is good if the current iterate is far from the solution.
c.  如果μ的值很⼩，那么hlm成了⾼斯⽜顿法的⽅向(适合迭代的最后阶段，⾮常接近最优解，避免了最速下降的震荡)
由此可见，惩罚因⼦既下降的⽅向⼜影响下降步⼦的⼤⼩。
2.    惩罚因⼦的计算[迭代步长计算]
我们的⽬标是求f的最⼩值，我们希望迭代开始时，惩罚因⼦μ被设定为较⼩的值，若
信赖域⽅法与线搜索技术⼀样，也是优化算法中的⼀种保证全局收敛的重要技术. 它们的功能都是在优化算法中求出每次迭代的位移, 从⽽确定新的迭代点.所不同的是: 线搜索技术是先产⽣位移⽅向(亦称为搜索⽅向), 然后确定位移的长度(亦称为搜索步长)。⽽信赖域技术则是直接确定位移, 产⽣新的迭代点。
信赖域⽅法的基本思想是:⾸先给定⼀个所谓的“信赖域半径”作为位移长度的上界，并以当前迭代点为中⼼以此“上界”为半径确定⼀个称之为“信赖域”的闭球区域。然后，通过求解这个区域内的“信赖域⼦问题”(⽬标函数的⼆次近似模型) 的最优点来确定“候选位移”。若候选位移能使⽬标函数值有充分的下降量, 则接受该候选位移作为新的位移，并保持或扩⼤信赖域半径, 继续新的迭代。否则, 说明⼆次模型与⽬标函数的近似度不够理想，需要缩⼩信赖域半径，再通过求解新的信赖域内的⼦问题得到新的候选位移。如此重复下去，直到满⾜迭代终⽌条件。
现在⽤信赖域⽅法解决之前的⽆约束线性规划：
如果q很⼤，说明L(h)⾮常接近F(x+h)，我们可以减少惩罚因⼦μ，以便于下次迭代此时算法更接近⾼斯⽜顿算法。如果q很⼩或者是负的，说明是poor approximation，我们需要增⼤惩罚因⼦，减少步长，此时算法更接近最速下降法。具体来说，
a.当q⼤于0时，此次迭代有效：
b．当q⼩于等于0时，此次迭代⽆效：
3.完整的算法流程及代码距离
LM的算法流程和⾼斯⽜顿⼏乎⼀样，只是迭代步长求法利⽤信赖域法
(1)给定初始点x(0)，允许误差ε>0，置k=0
(2)当f(x k+1)-f(x k)⼩于阈值ε时，算法退出，否则(3)
(3)x k+1=x k+hlm，代⼊f，返回(1)
两个例⼦还是沿⽤之前的。
例⼦1，根据美国1815年⾄1885年数据，估计⼈⼝模型中的参数A和B。如下表所⽰，已知年份和⼈⼝总量，及⼈⼝模型⽅程，求⽅程中的参数。
// A simple demo of Gauss-Newton algorithm on a user defined function
#include <cstdio>
#include <vector>
#include <opencv2/core/core.hpp>
using namespace std;
using namespace cv;
const double DERIV_STEP = 1e-5;
const int MAX_ITER = 100;
void LM(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer
const Mat &inputs, const Mat &outputs, Mat ¶ms);
double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer    const Mat &input, const Mat ¶ms, int n);
// The user defines their function here
double Func(const Mat &input, const Mat ¶ms);
int main()
{
// For this demo we're going to try and fit to the function
// F = A*exp(t*B)
// There are 2 parameters: A B
int num_params = 2;
// Generate random data using these parameters
int total_data = 8;
Mat inputs(total_data, 1, CV_64F);
Mat outputs(total_data, 1, CV_64F);
//load observation data
for(int i=0; i < total_data; i++) {
inputs.at<double>(i,0) = i+1;  //load year
}
//load America population
outputs.at<double>(0,0)= 8.3;
outputs.at<double>(1,0)= 11.0;
outputs.at<double>(2,0)= 14.7;
outputs.at<double>(3,0)= 19.7;
outputs.at<double>(4,0)= 26.7;
outputs.at<double>(4,0)= 26.7;
outputs.at<double>(5,0)= 35.2;
outputs.at<double>(6,0)= 44.4;
outputs.at<double>(7,0)= 55.9;
identity matrix是什么意思
// Guess the parameters, it should be close to the true value, else it can fail for very sensitive functions!    Mat params(num_params, 1, CV_64F);
//init guess
params.at<double>(0,0) = 6;
params.at<double>(1,0) = 0.3;
LM(Func, inputs, outputs, params);
printf("Parameters from GaussNewton: %lf %lf\n", params.at<double>(0,0), params.at<double>(1,0));
return 0;
}
double Func(const Mat &input, const Mat ¶ms)
{
// Assumes input is a single row matrix
// Assumes params is a column matrix
double A = params.at<double>(0,0);
double B = params.at<double>(1,0);
double x = input.at<double>(0,0);
return A*exp(x*B);
}
//calc the n-th params' partial derivation ， the params are our  final target
double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), const Mat &input, const Mat ¶ms, int n) {
// Assumes input is a single row matrix
// Returns the derivative of the nth parameter
Mat params1 = params.clone();
Mat params2 = params.clone();
// Use central difference  to get derivative
params1.at<double>(n,0) -= DERIV_STEP;
params2.at<double>(n,0) += DERIV_STEP;
double p1 = Func(input, params1);
double p2 = Func(input, params2);
double d = (p2 - p1) / (2*DERIV_STEP);
return d;
}
void LM(double(*Func)(const Mat &input, const Mat ¶ms),
const Mat &inputs, const Mat &outputs, Mat ¶ms)
{
int m = ws;
int n = ls;
int num_params = ws;
Mat r(m, 1, CV_64F); // residual matrix
Mat r_tmp(m, 1, CV_64F);
Mat Jf(m, num_params, CV_64F); // Jacobian of Func()
Mat input(1, n, CV_64F); // single row input
Mat params_tmp = params.clone();

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。