lm opencv 算法_Levenberg–Marquardt算法学习(和matlab的LM算法对比)
回顾高斯牛顿算法,引入LM算法
惩罚因子的计算(迭代步子的计算)
完整的算法流程及代码样例
1. 回顾高斯牛顿,引入LM算法 根据之前的博文:Gauss-Newton算法学习
假设我们研究如下形式的非线性最小二乘问题:
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(xk+1)-f(xk)小于阈值ε时,算法退出,否则(3)
(3)xk+1=xk+hlm,代入f,返回(1)
两个例子还是沿用之前的。
例子1,根据美国1815年至1885年数据,估计人口模型中的参数A和B。如下表所示,已知年份和人口总量,及人口模型方程,求方程中的参数。
// A simple demo of Gauss-Newton algorithm on a user defined function
#include
#include
#include
usingnamespacestd;
usingnamespacecv;
constdoubleDERIV_STEP = 1e-5;
constintMAX_ITER = 100;
voidLM(double(*Func)(constMat &input,constMat ¶ms),// function pointer
constMat &inputs,constMat &outputs, Mat ¶ms);
doubleDeriv(double(*Func)(constMat &input,constMat ¶ms),// function pointer
constMat &input,constMat ¶ms,intn);
// The user defines their function here
doubleFunc(constMat &input,constMat ¶ms);
intmain()
{
// For this demo we're going to try and fit to the function
// F = A*exp(t*B)
// There are 2 parameters: A B
intnum_params = 2;
// Generate random data using these parameters
inttotal_data = 8;
Mat inputs(total_data, 1, CV_64F);
Mat outputs(total_data, 1, CV_64F);
//load observation data
for(inti=0; i
inputs.at(i,0) = i+1;//load year
}
//load America population
outputs.at(0,0)= 8.3;
outputs.at(1,0)= 11.0;
outputs.at(2,0)= 14.7;
outputs.at(3,0)= 19.7;
outputs.at(4,0)= 26.7;
outputs.at(5,0)= 35.2;
outputs.at(6,0)= 44.4;
outputs.at(7,0)= 55.9;
// 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(0,0) = 6;
params.at(1,0) = 0.3;
LM(Func, inputs, outputs, params);
printf("Parameters from GaussNewton: %lf %lf\n", params.at(0,0), params.at(1,0));
return0;
}
doubleFunc(constMat &input,constMat ¶ms)
{
// Assumes input is a single row matrix
// Assumes params is a column matrix
doubleA = params.at(0,0);
doubleB = params.at(1,0);
doublex = input.at(0,0);
returnA*exp(x*B);
}
//calc the n-th params' partial derivation , the params are our final target
doubleDeriv(double(*Func)(constMat &input,constMat ¶ms),constMat &input,constMat ¶ms,intn)
{
// 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(n,0) -= DERIV_STEP;
params2.at(n,0) += DERIV_STEP;
doublep1 = Func(input, params1);
doublep2 = Func(input, params2);
doubled = (p2 - p1) / (2*DERIV_STEP);
returnd;
}
voidLM(double(*Func)(constMat &input,constMat ¶ms),
constMat &inputs,constMat &outputs, Mat ¶ms)
{
intm = inputs.rows;
intn = inputs.cols;
intnum_params = params.rows;
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();
doublelast_mse = 0;
floatu = 1, v = 2;
Mat I = Mat::ones(num_params, num_params, CV_64F);//construct identity matrix
inti =0;
for(i=0; i
doublemse = 0;
doublemse_temp = 0;
for(intj=0; j
for(intk=0; k
input.at(0,k) = inputs.at(j,k);
}
r.at(j,0) = outputs.at(j,0) - Func(input, params);//diff between previous estimate and observation population
mse += r.at(j,0)*r.at(j,0);
for(intk=0; k
Jf.at(j,k) = Deriv(Func, input, params, k);//construct jacobian matrix
}
}
mse /= m;
params_tmp = params.clone();
Mat hlm = (Jf.t()*Jf + u*I).inv()*Jf.t()*r;
params_tmp += hlm;
for(intj=0; j
r_tmp.at(j,0) = outputs.at(j,0) - Func(input, params_tmp);//diff between current estimate and observation population
mse_temp += r_tmp.at(j,0)*r_tmp.at(j,0);
}
mse_temp /= m;
Mat q(1,1,CV_64F);
q = (mse - mse_temp)/(0.5*hlm.t()*(u*hlm-Jf.t()*r));
doubleq_value = q.at(0,0);
if(q_value>0)
{
doubles = 1.0/3.0;
v = 2;
mse = mse_temp;
params = params_tmp;
doubletemp = 1 - pow(2*q_value-1,3);
if(s>temp)
{
u = u * s;
}else
{
u = u * temp;
}
}else
{
u = u*v;
v = 2*v;
params = params_tmp;
}
// The difference in mse is very small, so quit
if(fabs(mse - last_mse)
break;
}
//printf("%d: mse=%f\n", i, mse);
printf("%d %lf\n", i, mse);
last_mse = mse;
}
}
A=7.0,B=0.26(初始值,A=6,B=0.3),100次迭代到第7次收敛了。和之前差不多,但是LM对于初始的选择是非常敏感的,如果A=6,B=6,则拟合失败!
我调用了matlab的接口跑LM,结果也是一样错误的,图片上可以看到拟合失败
clc;
clear;
a0=[6,6];
options=optimset('Algorithm','Levenberg-Marquardt','Display','iter');
data_1=[1 2 3 4 5 6 7 8];
obs_1=[8.3 11.0 14.7 19.7 26.7 35.2 44.4 55.9];
a=lsqnonlin(@myfun,a0,[],[],options,data_1,obs_1);
plot(data_1,obs_1,'o');
hold on
plot(data_1,a(1)*exp(a(2)*data_1),'b');
plot(data_1,7*exp(a(2)*data_1),'b');
%hold off
a(1)
a(2)
function E = myfun(a, x,y)
%这是一个测试文件用于测试 lsqnonlin
% Detailed explanation goes here
x=x(:);
y=y(:);
Y=a(1)*exp(a(2)*x);
E=y-Y;
end
最后一个点拟合失败的,所以函数不对的
因此虽然莱文博格-马夸特迭代法能够自适应的在高斯牛顿和最速下降法之间调整,既可保证在收敛较慢时迭代过程总是下降的,又可保证迭代过程在解的邻域内迅速收敛。但是,LM对于初始点选择还是比较敏感的!
例子2:我想要拟合如下模型,
由于缺乏观测量,就自导自演,假设4个参数已知A=5,B=1,C=10,D=2,构造100个随机数作为x的观测值,计算相应的函数观测值。然后,利用这些观测值,反推4个参数。
// A simple demo of Gauss-Newton algorithm on a user defined function
#include
#include
#include
usingnamespacestd;
usingnamespacecv;
constdoubleDERIV_STEP = 1e-5;
constintMAX_ITER = 100;
voidLM(double(*Func)(constMat &input,constMat ¶ms),// function pointer
constMat &inputs,constMat &outputs, Mat ¶ms);
doubleDeriv(double(*Func)(constMat &input,constMat ¶ms),// function pointer
constMat &input,constMat ¶ms,intn);
// The user defines their function here
doubleFunc(constMat &input,constMat ¶ms);
intmain()
{
// For this demo we're going to try and fit to the function
// F = A*sin(Bx) + C*cos(Dx)
// There are 4 parameters: A, B, C, D
intnum_params = 4;
// Generate random data using these parameters
inttotal_data = 100;
doubleA = 5;
doubleB = 1;
doubleC = 10;
doubleD = 2;
Mat inputs(total_data, 1, CV_64F);
Mat outputs(total_data, 1, CV_64F);
for(inti=0; i
doublex = -10.0 + 20.0* rand() / (1.0 + RAND_MAX);// random between [-10 and 10]
doubley = A*sin(B*x) + C*cos(D*x);
// Add some noise
// y += -1.0 + 2.0*rand() / (1.0 + RAND_MAX);
inputs.at(i,0) = x;
outputs.at(i,0) = y;
}
// 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);
params.at(0,0) = 1;
params.at(1,0) = 1;
params.at(2,0) = 8;// changing to 1 will cause it not to find the solution, too far away
params.at(3,0) = 1;
LM(Func, inputs, outputs, params);
printf("True parameters: %f %f %f %f\n", A, B, C, D);
printf("Parameters from GaussNewton: %f %f %f %f\n", params.at(0,0), params.at(1,0),
params.at(2,0), params.at(3,0));
return0;
}
doubleFunc(constMat &input,constMat ¶ms)
{
// Assumes input is a single row matrix
// Assumes params is a column matrix
doubleA = params.at(0,0);
doubleB = params.at(1,0);
doubleC = params.at(2,0);
doubleD = params.at(3,0);
doublex = input.at(0,0);
returnA*sin(B*x) + C*cos(D*x);
}
//calc the n-th params' partial derivation , the params are our final target
doubleDeriv(double(*Func)(constMat &input,constMat ¶ms),constMat &input,constMat ¶ms,intn)
{
// 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(n,0) -= DERIV_STEP;
params2.at(n,0) += DERIV_STEP;
doublep1 = Func(input, params1);
doublep2 = Func(input, params2);
doubled = (p2 - p1) / (2*DERIV_STEP);
returnd;
}
voidLM(double(*Func)(constMat &input,constMat ¶ms),
constMat &inputs,constMat &outputs, Mat ¶ms)
{
intm = inputs.rows;
intn = inputs.cols;
intnum_params = params.rows;
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();
doublelast_mse = 0;
floatu = 1, v = 2;
Mat I = Mat::ones(num_params, num_params, CV_64F);//construct identity matrix
inti =0;
for(i=0; i
doublemse = 0;
doublemse_temp = 0;
for(intj=0; j
for(intk=0; k
input.at(0,k) = inputs.at(j,k);
}
r.at(j,0) = outputs.at(j,0) - Func(input, params);//diff between estimate and observation population
mse += r.at(j,0)*r.at(j,0);
for(intk=0; k
Jf.at(j,k) = Deriv(Func, input, params, k);//construct jacobian matrix
}
}
mse /= m;
params_tmp = params.clone();
Mat hlm = (Jf.t()*Jf + u*I).inv()*Jf.t()*r;
params_tmp += hlm;
for(intj=0; j
r_tmp.at(j,0) = outputs.at(j,0) - Func(input, params_tmp);//diff between estimate and observation population
mse_temp += r_tmp.at(j,0)*r_tmp.at(j,0);
}
mse_temp /= m;
Mat q(1,1,CV_64F);
q = (mse - mse_temp)/(0.5*hlm.t()*(u*hlm-Jf.t()*r));
doubleq_value = q.at(0,0);
if(q_value>0)
{
doubles = 1.0/3.0;
v = 2;
mse = mse_temp;
params = params_tmp;
doubletemp = 1 - pow(2*q_value-1,3);
if(s>temp)
{
u = u * s;
}else
{
u = u * temp;
}
}else
{
u = u*v;
v = 2*v;
params = params_tmp;
}
// The difference in mse is very small, so quit
if(fabs(mse - last_mse)
break;
}
//printf("%d: mse=%f\n", i, mse);
printf("%d %lf\n", i, mse);
last_mse = mse;
}
}
我们看到迭代了100次,结果几何和高斯牛顿算出来是一样的。我们绘制LM和高斯牛顿的残差函数收敛过程,发现LM一直是总体下降的,没有太多反复。
高斯牛顿:
LM:
lm opencv 算法_Levenberg–Marquardt算法学习(和matlab的LM算法对比)相关推荐
- matlab切割肿瘤算法,ML之RF:基于Matlab利用RF算法实现根据乳腺肿瘤特征向量高精度(better)预测肿瘤的是恶性还是良性...
ML之RF:基于Matlab利用RF算法实现根据乳腺肿瘤特征向量高精度(better)预测肿瘤的是恶性还是良性 目录 输出结果 实现代码 输出结果 更新-- 实现代码 %RF:RF实现根据乳腺肿瘤特征 ...
- matlab粒子群算法求解无约束最小值,pso matlab粒子群算法和遗传 是解决约束优化问题,无 和多目标 的优 259万源代码下载- www.pudn.com...
文件名称: pso下载 收藏√ [ 5 4 3 2 1 ] 开发工具: matlab 文件大小: 51 KB 上传时间: 2016-06-01 下载次数: 0 提 供 者: 孙志勇 详细说 ...
- 算法工程师的必备学习资料,《AI算法工程师手册》正式开源了
2019-05-14 23:41:00 前言 最近前阿里的一位工程师开源了一份网页版的算法工程师学习手册,没有纸质版的图书,直接在线开源,小编去看看了一下,总结的非常到位,几乎涵盖的机器学习.深度学习 ...
- 计算机视觉与深度学习 | 基于MATLAB的Vibe算法消除鬼影(代码版)
/********************************************************** github:https://github.com/MichaelBeechan ...
- python数据结构和算法讲解_【学习】python数据结构和算法
二.算法分析 2.2 什么是算法分析 大O表示法 image.png 2.3 python数据结构的性能 列表 image.png 字典 image.png 说一下list[index]的o(1)原理 ...
- 【CT算法,radon变换】基于MATLAB的CT算法,radon变换的三维建模仿真
1.软件版本 MATLAB2021a 2.本算法理论知识 1.输入:T(x,y,z) 使用stl读取函数完成T的导入工作 2.做Radon变换,得投影图:P 正常Radon变换即可. 3.对P:应用斜 ...
- 基于matlab的捕食算法,【优化求解】基于matlab细菌觅食算法的函数优化分析【含Matlab源码 217期】...
一.简介 实际生活需求促进了最优化方法的发展.近半个多世纪以来,由于传统优化方法的不足,一些具有全局优化性能且通用性强的进化算法,因其高效的优化性能.无需问题精确描述信息等优点,受到各领域广泛的关注和 ...
- matlab的多变量dmc源程序,基于MATLAB多变量DMC算法的仿真技术研究
基于MATLAB多变量DMC算法的仿真技术研究 基于MATLAB多变量DMC算法的仿真技术研究 作者:李凤霞 于佐军 来源:<科技创新导报>2011年第17期 摘 要:利用MATLAB开发 ...
- 用matlab实现理查森外推算法,Matlab数值积分(2)
实验目的: 掌握理查森外推法 实验要求: 1. 给出理查森外推算法 2. 用Matlab实现理查森外推算法 3. 用Matlab实现自适应积分算法 实验内容: 1. 理查森外推算法,数学知识:利用Ri ...
最新文章
- python使用fpdf创建页眉、页脚并嵌入图片
- VMVMware-workstation以及CentOS-7安装
- GTX1060 6G是低端电脑显卡吗?
- Luogu1613 跑路
- 扫雷游戏网页版_借“买量”造爆款,梦幻西游网页版击穿H5游戏天花板
- 【系统分析师之路】2010年系统分析师上午综合知识真题
- ios 手势返回监听方法
- 【小技巧】2345——劫持IE浏览器主页
- Zigbee 协议栈网络管理
- WPS Office宏病毒实现shell反弹
- NAS网络配置、资源管理和用户访问权限
- Python之数据容器
- java 发送notes_JAVA使用B/S模式(网页)发送Notes邮件
- 腾讯 Techo Hub 2022 年首站落地福州|723,与开发者们探讨工业数字化!
- 腾讯云内容生态助力猿辅导,线上线下全方位推动教育云进入快车道
- Excel批注教学:一键给多个单元格添加相同批注
- WebView部分源码概览
- python数字转大写字母_python变量名称如何转化为大写字母?
- 八、CSS3的美化背景与边框
- agent开发之oneAgent
热门文章
- POJ1917 UVA10361 Automatic Poetry【文本】
- POJ3070 Fibonacci【矩阵快速幂】
- 基础训练(六~十)题解
- CCF NOI1011 正方形
- HDU2021 发工资咯:)【入门】
- TensorFlow 学习(二)—— tf.Graph、tf.Session() 与 tf.Session().run()
- fatal error: caffe/proto/caffe.pb.h: No such file or directory
- 仿射变换(Affine transformation)与python实践
- incrby redis 最大值_Redis 的 8 大数据类型,写得非常好!
- python零基础入门教程-零基础入门Python爬虫不知道怎么学?这是入门的完整教程...