Gauss-Newton算法学习
Gauss-Newton算法是解决非线性最优问题的常见算法之一,最近研读开源项目代码,又碰到了,索性深入看下。本次讲解内容如下:
- 基本数学名词识记
- 牛顿法推导、算法步骤、计算实例
- 高斯牛顿法推导(如何从牛顿法派生)、算法步骤、编程实例
- 高斯牛顿法优劣总结
一、基本概念定义
1.非线性方程定义及最优化方法简述
指因变量与自变量之间的关系不是线性的关系,比如平方关系、对数关系、指数关系、三角函数关系等等。对于此类方程,求解n元实函数f在整个n维向量空间Rn上的最优值点往往很难得到精确解,经常需要求近似解问题。
求解该最优化问题的方法大多是逐次一维搜索的迭代算法,基本思想是在一个近似点处选定一个有利于搜索方向,沿这个方向进行一维搜索,得到新的近似点。如此反复迭代,知道满足预定的精度要求为止。根据搜索方向的取法不同,这类迭代算法可分为两类:
解析法:需要用目标函数的到函数,
梯度法:又称最速下降法,是早期的解析法,收敛速度较慢
牛顿法:收敛速度快,但不稳定,计算也较困难。高斯牛顿法基于其改进,但目标作用不同
共轭梯度法:收敛较快,效果好
变尺度法:效率较高,常用DFP法(Davidon Fletcher Powell)
直接法:不涉及导数,只用到函数值。有交替方向法(又称坐标轮换法)、模式搜索法、旋转方向法、鲍威尔共轭方向法和单纯形加速法等。
2.非线性最小二乘问题
非线性最小二乘问题来自于非线性回归,即通过观察自变量和因变量数据,求非线性目标函数的系数参数,使得函数模型与观测量尽量相似。
高斯牛顿法解决非线性最小二乘问题的最基本方法,并且它只能处理二次函数。(使用时必须将目标函数转化为二次的)
Unlike Newton'smethod, the Gauss–Newton algorithm can only be used to minimize a sum ofsquared function values
3.基本数学表达
a.梯度gradient,由多元函数的各个偏导数组成的向量
以二元函数为例,其梯度为:
b.黑森矩阵Hessian matrix,由多元函数的二阶偏导数组成的方阵,描述函数的局部曲率,以二元函数为例,
c.雅可比矩阵 Jacobian matrix,是多元函数一阶偏导数以一定方式排列成的矩阵,体现了一个可微方程与给出点的最优线性逼近。以二元函数为例,
如果扩展多维的话F: Rn-> Rm,则雅可比矩阵是一个m行n列的矩阵:
雅可比矩阵作用,如果P是Rn中的一点,F在P点可微分,那么在这一点的导数由JF(P)给出,在此情况下,由F(P)描述的线性算子即接近点P的F的最优线性逼近:
d.残差 residual,表示实际观测值与估计值(拟合值)之间的差
二、牛顿法
牛顿法的基本思想是采用多项式函数来逼近给定的函数值,然后求出极小点的估计值,重复操作,直到达到一定精度为止。
1.考虑如下一维无约束的极小化问题:
因此,一维牛顿法的计算步骤如下:
需要注意的是,牛顿法在求极值的时候,如果初始点选取不好,则可能不收敛于极小点
2.下面给出多维无约束极值的情形:
若非线性目标函数f(x)具有二阶连续偏导,在x(k)为其极小点的某一近似,在这一点取f(x)的二阶泰勒展开,即:
如果f(x)是二次函数,则其黑森矩阵H为常数,式(1)是精确的(等于号),在这种情况下,从任意一点处罚,用式(2)只要一步可求出f(x)的极小点(假设黑森矩阵正定,所有特征值大于0)
如果f(x)不是二次函数,式(1)仅是一个近似表达式,此时,按式(2)求得的极小点,只是f(x)的近似极小点。在这种情况下,常按照下面选取搜索方向:
牛顿法收敛的速度很快,当f(x)的二阶导数及其黑森矩阵的逆矩阵便于计算时,这一方法非常有效。【但通常黑森矩阵很不好求】
3.下面给出一个实际计算例子。
例:试用牛顿法求的极小值
解:
【f(x)是二次函数,H矩阵为常数,只要任意点出发,只要一步即可求出极小点】
三、牛顿高斯法
1. gauss-newton是如何由上述派生的
有时候为了拟合数据,比如根据重投影误差求相机位姿(R,T为方程系数),常常将求解模型转化为非线性最小二乘问题。高斯牛顿法正是用于解决非线性最小二乘问题,达到数据拟合、参数估计和函数估计的目的。
假设我们研究如下形式的非线性最小二乘问题:
这两个位置间残差(重投影误差):
如果有大量观测点(多维),我们可以通过选择合理的T使得残差的平方和最小求得两个相机之间的位姿。机器视觉这块暂时不扩展,接着说怎么求非线性最小二乘问题。
若用牛顿法求式3,则牛顿迭代公式为:
看到这里大家都明白高斯牛顿和牛顿法的差异了吧,就在这迭代项上。经典高斯牛顿算法迭代步长λ为1.
那回过头来,高斯牛顿法里为啥要舍弃黑森矩阵的二阶偏导数呢?主要问题是因为牛顿法中Hessian矩阵中的二阶信息项通常难以计算或者花费的工作量很大,而利用整个H的割线近似也不可取,因为在计算梯度时已经得到J(x),这样H中的一阶信息项JTJ几乎是现成的。鉴于此,为了简化计算,获得有效算法,我们可用一阶导数信息逼近二阶信息项。注意这么干的前提是,残差r接近于零或者接近线性函数从而接近与零时,二阶信息项才可以忽略。通常称为“小残量问题”,否则高斯牛顿法不收敛。
3. 举例
接下来的代码里并没有保证算法收敛的机制,在例子2的自嗨中可以看到劣势。关于自变量维数,代码可以支持多元,但两个例子都是一维的,比如例子1中只有年份t,其实可以增加其他因素的,不必在意。
例子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 GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms), // function pointerconst Mat &inputs, const Mat &outputs, Mat ¶ms);double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), // function pointerconst 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 Bint num_params = 2;// Generate random data using these parametersint total_data = 8;Mat inputs(total_data, 1, CV_64F);Mat outputs(total_data, 1, CV_64F);//load observation datafor(int i=0; i < total_data; i++) {inputs.at<double>(i,0) = i+1; //load year}//load America populationoutputs.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>(5,0)= 35.2;outputs.at<double>(6,0)= 44.4;outputs.at<double>(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 guessparams.at<double>(0,0) = 6;params.at<double>(1,0) = 0.3;GaussNewton(Func, inputs, outputs, params);printf("Parameters from GaussNewton: %f %f\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 matrixdouble 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 parameterMat params1 = params.clone();Mat params2 = params.clone();// Use central difference to get derivativeparams1.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 GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms),const Mat &inputs, const Mat &outputs, Mat ¶ms)
{int m = inputs.rows;int n = inputs.cols;int num_params = params.rows;Mat r(m, 1, CV_64F); // residual matrixMat Jf(m, num_params, CV_64F); // Jacobian of Func()Mat input(1, n, CV_64F); // single row inputdouble last_mse = 0;for(int i=0; i < MAX_ITER; i++) {double mse = 0;for(int j=0; j < m; j++) {for(int k=0; k < n; k++) {//copy Independent variable vector, the yearinput.at<double>(0,k) = inputs.at<double>(j,k);}r.at<double>(j,0) = outputs.at<double>(j,0) - Func(input, params);//diff between estimate and observation populationmse += r.at<double>(j,0)*r.at<double>(j,0);for(int k=0; k < num_params; k++) {Jf.at<double>(j,k) = Deriv(Func, input, params, k);}}mse /= m;// The difference in mse is very small, so quitif(fabs(mse - last_mse) < 1e-8) {break;}Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;params += delta;//printf("%d: mse=%f\n", i, mse);printf("%d %f\n", i, mse);last_mse = mse;}
}
运行结果:
A=7.0,B=0.26 (初始值,A=6,B=0.3),100次迭代到第4次就收敛了。
若初始值A=1,B=1,则要迭代14次收敛。
下图为根据上面得到的A、B系数,利用matlab拟合的人口模型曲线
例子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 <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 GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms), // function pointerconst Mat &inputs, const Mat &outputs, Mat ¶ms);double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), // function pointerconst 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*sin(Bx) + C*cos(Dx)// There are 4 parameters: A, B, C, Dint num_params = 4;// Generate random data using these parametersint total_data = 100;double A = 5;double B = 1;double C = 10;double D = 2;Mat inputs(total_data, 1, CV_64F);Mat outputs(total_data, 1, CV_64F);for(int i=0; i < total_data; i++) {double x = -10.0 + 20.0* rand() / (1.0 + RAND_MAX); // random between [-10 and 10]double y = A*sin(B*x) + C*cos(D*x);// Add some noise// y += -1.0 + 2.0*rand() / (1.0 + RAND_MAX);inputs.at<double>(i,0) = x;outputs.at<double>(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<double>(0,0) = 1;params.at<double>(1,0) = 1;params.at<double>(2,0) = 8; // changing to 1 will cause it not to find the solution, too far awayparams.at<double>(3,0) = 1;GaussNewton(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<double>(0,0), params.at<double>(1,0),params.at<double>(2,0), params.at<double>(3,0));return 0;
}double Func(const Mat &input, const Mat ¶ms)
{// Assumes input is a single row matrix// Assumes params is a column matrixdouble A = params.at<double>(0,0);double B = params.at<double>(1,0);double C = params.at<double>(2,0);double D = params.at<double>(3,0);double x = input.at<double>(0,0);return A*sin(B*x) + C*cos(D*x);
}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 parameterMat params1 = params.clone();Mat params2 = params.clone();// Use central difference to get derivativeparams1.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 GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms),const Mat &inputs, const Mat &outputs, Mat ¶ms)
{int m = inputs.rows;int n = inputs.cols;int num_params = params.rows;Mat r(m, 1, CV_64F); // residual matrixMat Jf(m, num_params, CV_64F); // Jacobian of Func()Mat input(1, n, CV_64F); // single row inputdouble last_mse = 0;for(int i=0; i < MAX_ITER; i++) {double mse = 0;for(int j=0; j < m; j++) {for(int k=0; k < n; k++) {input.at<double>(0,k) = inputs.at<double>(j,k);}r.at<double>(j,0) = outputs.at<double>(j,0) - Func(input, params);mse += r.at<double>(j,0)*r.at<double>(j,0);for(int k=0; k < num_params; k++) {Jf.at<double>(j,k) = Deriv(Func, input, params, k);}}mse /= m;// The difference in mse is very small, so quitif(fabs(mse - last_mse) < 1e-8) {break;}Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;params += delta;//printf("%d: mse=%f\n", i, mse);printf("%f\n",mse);last_mse = mse;}
}
运行结果,得到的参数并不够理想,50次后收敛了
下图中,每次迭代残差并没有持续减少,有反复
4.优缺点分析
优点:
对于零残量问题,即r=0,有局部二阶收敛速度
对于小残量问题,即r较小或接近线性,有快的局部收敛速度
对于线性最小二乘问题,一步达到极小点
缺点:
对于不是很严重的大残量问题,有较慢的局部收敛速度
对于残量很大的问题或r的非线性程度很大的问题,不收敛
不一定总体收敛
如果J不满秩,则方法无定义
对于它的缺点,我们通过增加线性搜索策略,保证目标函数每一步下降,对于几乎所有非线性最小二乘问题,它都具有局部收敛性及总体收敛,即所谓的阻尼高斯牛顿法。
其中,ak为一维搜索因子
Gauss-Newton算法学习相关推荐
- 高斯牛顿(Gauss Newton)、列文伯格-马夸尔特(Levenberg-Marquardt)最优化算法与VSLAM
转载请说明出处:http://blog.csdn.net/zhubaohua_bupt/article/details/74973347 在VSLAM优化部分,我们多次谈到,构建一个关于待优化位姿的误 ...
- 非线性最小二乘法之Gauss Newton、L-M、Dog-Leg
非线性最小二乘法之Gauss Newton.L-M.Dog-Leg 最快下降法 假设hTF′(x)<0h^TF'(x) ,则h是F(x)F(x)下降方向,即对于任意足够小的α>0\alph ...
- 拿下斯坦福和剑桥双offer,00后的算法学习之路
董文馨,00后,精通英语,西班牙语.斯坦福大学计算机系和剑桥大学双Offer,秋季将进入斯坦福大学学习. 10岁开始在国外上学:12岁学Scratch: 13岁学HTML & CSS: 14岁 ...
- 好久没有看到这么有建设性德文章,由衷地赞叹《知其所以然地学习(以算法学习为例)》-By 刘未鹏(pongba)
知其所以然地学习(以算法学习为例) By 刘未鹏(pongba) C++的罗浮宫(http://blog.csdn.net/pongba) Updated(2008-7-24):更新见正文部分,有标注 ...
- 原创 | 初学者友好!最全算法学习资源汇总(附链接)
在计算机发展飞速的今天,也许有人会问,"今天计算机这么快,算法还重要吗?"其实永远不会有太快的计算机,因为我们总会想出新的应用.虽然在摩尔定律的作用下,计算机的计算能力每年都在飞快 ...
- 基本算法学习(一)之希尔排序(JS)
参考书: 严蔚敏-数据结构 希尔排序(Shell's Sort) 希尔排序又称"缩小增量排序",归属于插入排序一类,简单来说,和我们的插入排序比,它更快. 奇妙的记忆点: 内排序( ...
- 大顶堆删除最大值_算法学习笔记(47): 二叉堆
堆(Heap)是一类数据结构,它们拥有树状结构,且能够保证父节点比子节点大(或小).当根节点保存堆中最大值时,称为大根堆:反之,则称为小根堆. 二叉堆(Binary Heap)是最简单.常用的堆,是一 ...
- Surf算法学习心得(一)——算法原理
Surf算法学习心得(一)--算法原理 写在前面的话: Surf算法是对Sift算法的一种改进,主要是在算法的执行效率上,比Sift算法来讲运行更快!由于我也是初学者,刚刚才开始研究这个算法,然而网上 ...
- 算法学习:后缀自动机
[前置知识] AC自动机(没有什么关联,但是看懂了会对后缀自动机有不同的理解) [解决问题] 各种子串的问题 [算法学习] 学习后缀自动机的过程中,看到了许多相关性质和证明,但是奈何才疏学浅(lan) ...
最新文章
- 数学仍然是人类的“火炬”
- 孙正义下重金的机械臂独角兽梦碎:估值最高40亿美元,做披萨太难吃,只好去做披萨盒...
- ESXI使用记录---安装vSphere(VCSA)
- XCode6 生成prefix.pch文件
- TensorFlow 使用例子-LSTM实现序列标注
- Hibernate中使用Criteria查询及注解——(Dept.hbm.xml)
- matlab imwrite将图像保存到其他目录
- aix oracle 10.2.0.1 升级 10.2.0.4,AIX Oracle RAC 升级到10.2.0.4.0要特别注意的问题 - 爱肯的专栏 ......
- GDPR到底是如何影响机器学习的?
- mysql通配符like,不吃透都对不起自己
- php 多维数组怎么表达,php 对多维数组的操作,该怎么解决
- 「一本通 4.1 例 3」校门外的树 (loj10115)
- (python)域名查询服务(whois)程序+检查5位以内域名到期时间邮件报警抢注域名
- 打造立体文案矩阵库之二:直复式营销文案
- 清理服务器 归档日志文件,服务器归档日志模式
- 现在可以把小程序交给第三方开发或管理了
- 软件开发中,做产品与做项目有什么区别?
- 男人就应该对自己狠一点
- 华为防火墙用户与认证
- 如何在官网下载tomcat
热门文章
- treasure what you have now
- what you should do in the morning?
- iPad+MacBook+安卓手机的图书馆工作方案!超高效率!堪比移动工作站!
- net.sf包JSONArray与JSONObject遍历
- 第四章:滚动堆栈(1)
- 开发:随笔记录之 Json字符串和对象的相互转换
- 五:二叉树中和为某一直的路径
- REMOTE_ADDR,HTTP_CLIENT_IP,HTTP_X_FORWARDED
- 如何在linux系统自动mount一个NTFS分区只读方式挂载
- symbian 视频播放解决方案