非线性最小二乘法之Gauss Newton、L-M、Dog-Leg9
最优化理论·非线性最小二乘
最优化理论·非线性最小二乘_努力努力努力-CSDN博客
最快下降法
理解:
a*b = |a|*|b|*cos()
最小二乘问题
通常的最小二乘问题都可以表示为:
GaussNewton
close all;
clc;a=1;b=2;c=1; %待求解的系数x=(0:0.01:1)';
w=rand(length(x),1)*2-1; %生成噪声
y_noise=exp(a*x.^2+b*x+c)+w; %带噪声的模型 pre=rand(3,1); %步骤1,待拟合的参数,初始值是随机数;for i=1:1000 f = exp(pre(1)*x.^2+pre(2)*x+pre(3)); % g = y_noise-f; % 步骤2,构建拟合的函数和真实函数的误差;视觉中,也是根据观测的值减去预测的值,然后优化其中的参数,使两者最小; p1 = exp(pre(1)*x.^2+pre(2)*x+pre(3)).*x.^2; % 残差函数g对a求偏导,在ceres中也是有对应的偏导项的p2 = exp(pre(1)*x.^2+pre(2)*x+pre(3)).*x; % 残差函数g对b求偏导;p3 = exp(pre(1)*x.^2+pre(2)*x+pre(3)); % 残差函数g对c求偏导;J = [p1 p2 p3]; % 步骤2中的雅克比矩阵delta = inv(J'*J)*J'* g; %步骤3,inv(J'*J)*J'为H的逆pcur = pre+delta; %步骤4if norm(delta) <1e-16break;endpre = pcur;
endyreal = exp(a*x.^2+b*x+c);
ypre=exp(pre(1)*x.^2+pre(2)*x+pre(3)); % 拟合的曲线figure(1)
plot(x,y_noise,'b')
hold on
plot(x,yreal,'r')
hold on
plot(x,ypre,'k');
grid on;
set(gca, 'GridLineStyle', '--'); % 设置为虚线
set(gca, 'GridAlpha', 1.0); % 设置透明度
xlabel('X');
ylabel('Y');
legend('带有噪声的真实曲线','不带噪声的真实曲线','拟合曲线')figure(2)
plot(x,yreal - ypre,'b')
grid on;
set(gca, 'GridLineStyle', '--'); % 设置为虚线
set(gca, 'GridAlpha', 1.0); % 设置透明度
xlabel('X');
ylabel('Y');
legend('真值和优化值之间的误差曲线')%比较一下
[a b c]
pre'
LM阻尼最小二乘法
Dog-Leg最小二乘法
举例
Gauss Newton代码:
double func(const VectorXd& input, const VectorXd& output, const VectorXd& params, double objIndex)
{// obj = A * sin(Bx) + C * cos(D*x) - Fdouble x1 = params(0);double x2 = params(1);double x3 = params(2);double x4 = params(3);double t = input(objIndex);double f = output(objIndex);return x1 * sin(x2 * t) + x3 * cos( x4 * t) - f;
}//return vector make up of func() element.
VectorXd objF(const VectorXd& input, const VectorXd& output, const VectorXd& params)
{VectorXd obj(input.rows());for(int i = 0; i < input.rows(); i++)obj(i) = func(input, output, params, i);return obj;
}//F = (f ^t * f)/2
double Func(const VectorXd& obj)
{return obj.squaredNorm()/2;
}double Deriv(const VectorXd& input, const VectorXd& output, int objIndex, const VectorXd& params,int paraIndex)
{VectorXd para1 = params;VectorXd para2 = params;para1(paraIndex) -= DERIV_STEP;para2(paraIndex) += DERIV_STEP;double obj1 = func(input, output, para1, objIndex);double obj2 = func(input, output, para2, objIndex);return (obj2 - obj1) / (2 * DERIV_STEP);
}MatrixXd Jacobin(const VectorXd& input, const VectorXd& output, const VectorXd& params)
{int rowNum = input.rows();int colNum = params.rows();MatrixXd Jac(rowNum, colNum);for (int i = 0; i < rowNum; i++){for (int j = 0; j < colNum; j++){Jac(i,j) = Deriv(input, output, i, params, j);}}return Jac;
}void gaussNewton(const VectorXd& input, const VectorXd& output, VectorXd& params)
{int errNum = input.rows(); //error numint paraNum = params.rows(); //parameter numVectorXd obj(errNum);double last_sum = 0;int iterCnt = 0;while (iterCnt < MAX_ITER){obj = objF(input, output, params);double sum = 0;sum = Func(obj);cout << "Iterator index: " << iterCnt << endl;cout << "parameter: " << endl << params << endl;cout << "error sum: " << endl << sum << endl << endl;if (fabs(sum - last_sum) <= 1e-12)break;last_sum = sum;MatrixXd Jac = Jacobin(input, output, params);VectorXd delta(paraNum);delta = (Jac.transpose() * Jac).inverse() * Jac.transpose() * obj;params -= delta;iterCnt++;}
}
LM代码:
double maxMatrixDiagonale(const MatrixXd& Hessian)
{int max = 0;for(int i = 0; i < Hessian.rows(); i++){if(Hessian(i,i) > max)max = Hessian(i,i);}return max;
}//L(h) = F(x) + h^t*J^t*f + h^t*J^t*J*h/2
//deltaL = h^t * (u * h - g)/2
double linerDeltaL(const VectorXd& step, const VectorXd& gradient, const double u)
{double L = step.transpose() * (u * step - gradient);return L/2;
}void levenMar(const VectorXd& input, const VectorXd& output, VectorXd& params)
{int errNum = input.rows(); //error numint paraNum = params.rows(); //parameter num//initial parameter VectorXd obj = objF(input,output,params);MatrixXd Jac = Jacobin(input, output, params); //jacobinMatrixXd A = Jac.transpose() * Jac; //HessianVectorXd gradient = Jac.transpose() * obj; //gradient//initial parameter tao v epsilon1 epsilon2double tao = 1e-3;long long v = 2;double eps1 = 1e-12, eps2 = 1e-12;double u = tao * maxMatrixDiagonale(A);bool found = gradient.norm() <= eps1;if(found) return;double last_sum = 0;int iterCnt = 0;while (iterCnt < MAX_ITER){VectorXd obj = objF(input,output,params);MatrixXd Jac = Jacobin(input, output, params); //jacobinMatrixXd A = Jac.transpose() * Jac; //HessianVectorXd gradient = Jac.transpose() * obj; //gradientif( gradient.norm() <= eps1 ){cout << "stop g(x) = 0 for a local minimizer optimizer." << endl;break;}cout << "A: " << endl << A << endl; VectorXd step = (A + u * MatrixXd::Identity(paraNum, paraNum)).inverse() * gradient; //negtive Hlm.cout << "step: " << endl << step << endl;if( step.norm() <= eps2*(params.norm() + eps2) ){cout << "stop because change in x is small" << endl;break;} VectorXd paramsNew(params.rows());paramsNew = params - step; //h_lm = -step;//compute f(x)obj = objF(input,output,params);//compute f(x_new)VectorXd obj_new = objF(input,output,paramsNew);double deltaF = Func(obj) - Func(obj_new);double deltaL = linerDeltaL(-1 * step, gradient, u);double roi = deltaF / deltaL;cout << "roi is : " << roi << endl;if(roi > 0){params = paramsNew;u *= max(1.0/3.0, 1-pow(2*roi-1, 3));v = 2;}else{u = u * v;v = v * 2;}cout << "u = " << u << " v = " << v << endl;iterCnt++;cout << "Iterator " << iterCnt << " times, result is :" << endl << endl;}
}
Dog-Leg代码:
void dogLeg(const VectorXd& input, const VectorXd& output, VectorXd& params)
{int errNum = input.rows(); //error numint paraNum = params.rows(); //parameter numVectorXd obj = objF(input, output, params);MatrixXd Jac = Jacobin(input, output, params); //jacobinVectorXd gradient = Jac.transpose() * obj; //gradient//initial parameter tao v epsilon1 epsilon2double eps1 = 1e-12, eps2 = 1e-12, eps3 = 1e-12;double radius = 1.0;bool found = obj.norm() <= eps3 || gradient.norm() <= eps1;if(found) return;double last_sum = 0;int iterCnt = 0;while(iterCnt < MAX_ITER){VectorXd obj = objF(input, output, params);MatrixXd Jac = Jacobin(input, output, params); //jacobinVectorXd gradient = Jac.transpose() * obj; //gradientif( gradient.norm() <= eps1 ){cout << "stop F'(x) = g(x) = 0 for a global minimizer optimizer." << endl;break;}if(obj.norm() <= eps3){cout << "stop f(x) = 0 for f(x) is so small";break;}//compute how far go along stepest descent direction.double alpha = gradient.squaredNorm() / (Jac * gradient).squaredNorm();//compute gauss newton step and stepest descent step.VectorXd stepest_descent = -alpha * gradient;VectorXd gauss_newton = (Jac.transpose() * Jac).inverse() * Jac.transpose() * obj * (-1);double beta = 0;//compute dog-leg step.VectorXd dog_leg(params.rows());if(gauss_newton.norm() <= radius)dog_leg = gauss_newton;else if(alpha * stepest_descent.norm() >= radius)dog_leg = (radius / stepest_descent.norm()) * stepest_descent;else{VectorXd a = alpha * stepest_descent;VectorXd b = gauss_newton;double c = a.transpose() * (b - a);beta = (sqrt(c*c + (b-a).squaredNorm()*(radius*radius-a.squaredNorm()))-c)/(b-a).squaredNorm();dog_leg = alpha * stepest_descent + beta * (gauss_newton - alpha * stepest_descent);}cout << "dog-leg: " << endl << dog_leg << endl;if(dog_leg.norm() <= eps2 *(params.norm() + eps2)){cout << "stop because change in x is small" << endl;break;}VectorXd new_params(params.rows());new_params = params + dog_leg;cout << "new parameter is: " << endl << new_params << endl;//compute f(x)obj = objF(input,output,params);//compute f(x_new)VectorXd obj_new = objF(input,output,new_params);//compute delta F = F(x) - F(x_new)double deltaF = Func(obj) - Func(obj_new);//compute delat L =L(0)-L(dog_leg)double deltaL = 0;if(gauss_newton.norm() <= radius)deltaL = Func(obj);else if(alpha * stepest_descent.norm() >= radius)deltaL = radius*(2*alpha*gradient.norm() - radius)/(2.0*alpha);else{VectorXd a = alpha * stepest_descent;VectorXd b = gauss_newton;double c = a.transpose() * (b - a);beta = (sqrt(c*c + (b-a).squaredNorm()*(radius*radius-a.squaredNorm()))-c)/(b-a).squaredNorm();deltaL = alpha*(1-beta)*(1-beta)*gradient.squaredNorm()/2.0 + beta*(2.0-beta)*Func(obj);}double roi = deltaF / deltaL;if(roi > 0){params = new_params;}if(roi > 0.75){radius = max(radius, 3.0 * dog_leg.norm());}else if(roi < 0.25){radius = radius / 2.0;if(radius <= eps2*(params.norm()+eps2)){cout << "trust region radius is too small." << endl;break;}}cout << "roi: " << roi << " dog-leg norm: " << dog_leg.norm() << endl;cout << "radius: " << radius << endl;iterCnt++;cout << "Iterator " << iterCnt << " times" << endl << endl;}
}
main()
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Sparse>
#include <iostream>
#include <iomanip>
#include <math.h>using namespace std;
using namespace Eigen;const double DERIV_STEP = 1e-5;
const int MAX_ITER = 100;#define max(a,b) (((a)>(b))?(a):(b))int main(int argc, char* argv[])
{// obj = A * sin(Bx) + C * cos(D*x) - F//there are 4 parameter: A, B, C, D.int num_params = 4;//generate random data using these parameterint total_data = 100;VectorXd input(total_data);VectorXd output(total_data);double A = 5, B= 1, C = 10, D = 2;//load observation datafor (int i = 0; i < total_data; i++){//generate a random variable [-10 10]double x = 20.0 * ((random() % 1000) / 1000.0) - 10.0;double deltaY = 2.0 * (random() % 1000) /1000.0;double y = A*sin(B*x)+C*cos(D*x) + deltaY;input(i) = x;output(i) = y;}//gauss the parametersVectorXd params_gaussNewton(num_params);//init gaussparams_gaussNewton << 1.6, 1.4, 6.2, 1.7;VectorXd params_levenMar = params_gaussNewton;VectorXd params_dogLeg = params_gaussNewton;gaussNewton(input, output, params_gaussNewton);levenMar(input, output, params_levenMar);dogLeg(input, output, params_dogLeg);cout << "gauss newton parameter: " << endl << params_gaussNewton << endl << endl << endl;cout << "Levenberg-Marquardt parameter: " << endl << params_levenMar << endl << endl << endl;cout << "dog-leg parameter: " << endl << params_dogLeg << endl << endl << endl;
}
相关资料
最小二乘问题的四种解法——牛顿法,梯度下降法,高斯牛顿法和列文伯格-马夸特法的区别和联系 - 知乎
非线性最小二乘法之Gauss Newton、L-M、Dog-Leg9相关推荐
- 非线性最小二乘法之Gauss Newton、L-M、Dog-Leg
非线性最小二乘法之Gauss Newton.L-M.Dog-Leg 最快下降法 假设hTF′(x)<0h^TF'(x) ,则h是F(x)F(x)下降方向,即对于任意足够小的α>0\alph ...
- 【scipy】Python调用非线性最小二乘法
文章目录 简介与构造函数 迭代策略 雅可比矩阵 测试 简介与构造函数 在scipy中,非线性最小二乘法的目的是找到一组函数,使得误差函数的平方和最小,可以表示为如下公式 arg minfiF(x)= ...
- 高斯牛顿(Gauss Newton)、列文伯格-马夸尔特(Levenberg-Marquardt)最优化算法与VSLAM
转载请说明出处:http://blog.csdn.net/zhubaohua_bupt/article/details/74973347 在VSLAM优化部分,我们多次谈到,构建一个关于待优化位姿的误 ...
- 非线性最小二乘法 python_[数值计算] 数据拟合——非线性最小二乘法
前面分析了 线性最小二乘 的解法YcoFlegs:[数值计算] 数据拟合--线性最小二乘法zhuanlan.zhihu.com 现在来看另一个问题:非线性最小二乘法 1. 定义 首先是如何定义这里这 ...
- Ceres Solver 官方教程学习笔记(十二)——非线性最小二乘法建模Modeling Non-linear Least Squares (下)
这一部分主要是最后的Problem比较重要. 带条件的代价函数ConditionedCostFunction 这个类使用户可以在使用封装好的代价函数的同时,对残差值加入一定的条件限制.举个例子,现在已 ...
- 非线性最小二乘法原理
一.一些基础概念[1] 1.最小二乘法定义 找到一个x∗{x^*}x∗,使得下面的式子到达局部最小值, F(x)=12∑i=1m(fi(x))2(1-1)F(x) = \frac{1}{2}\sum\ ...
- PCL 非线性最小二乘法拟合圆柱
文章目录 一.简介 二.实现代码 三.实现效果 参考资料 一.简介 这里通过非线性最小二乘的方法来实现圆柱体的拟合,具体的计算过程如下所述: 图中, p p p为输入数据的点位置,求解的参数为柱体的轴 ...
- 【矩阵基础与维度分析】【公式细节推导】矩阵非线性最小二乘法泰勒展开
最小二乘法的一般形式 fif_ifi为某一时刻的误差 比如预测值与预测值的差值 将残差写成向量的形式 和上面的FFF相比 每个fff计算的还是某一时刻的误差 将全部的误差求平方再求和 就相当于它的转 ...
- 求解非线性最小二乘法 Eigen
// 利用Eigen 求解非线性最小二乘; // 示例:<span style="font-family: Arial, Helvetica, sans-serif;"> ...
最新文章
- java mvc 绑定_关于Java:Spring MVC:将请求属性绑定到控制器方法参数
- 解决mysql“Access denied for user root@localhost”
- mvp 在 flutter 中的应用
- myeclipse使用太卡问题
- 计算机专业实训计划,计算机专业实习计划范文
- 2018年的AI/ML惊喜及预测19年的走势(二)
- Java 获取向量中的最大值
- 小强的HTML5移动开发之路(27)—— JavaScript回顾2
- Kubernetes中StatefulSet介绍
- C#获取%AppData%路径的方法
- array(numpy) 的size 和Image(PIL) 的shape比较
- python3计算运行时间_性能分析之代码运行时间计算——Python timeit 模块 介绍
- LDA-math-认识Beta/Dirichlet分布
- 【数学建模】相关性分析 - 皮尔逊相关系数 斯皮尔曼相关系数
- 毕业设计-基于微信小程序的医生随访系统
- 微信小程序连接mysql
- 基于CPT构建网络,熟悉各层协议
- PARAMETERS、SELECTION OPTIONS和SELECTION SCREEN的全部用法
- gcc: error: : No such file or directory
- Linux | 将文件移动到指定文件夹中
热门文章
- Maekdown使用
- linux网卡通信延迟高,Linux性能优化-网络请求延迟变大
- 运维管理系统,效果稳定、可靠运行
- 我的世界java版怎么加整合包_我的世界完美世界整合包
- Cesium出现Unknown crs name: urn:ogc:def:crs:EPSG:xxxx
- 《820事件:钻哥对AI人工智能下重手!!!》(Yanlz+AI+人工智能+元宇宙+Python+神经网络+深度学习+语言识别+NFT+VR+AR+MR+XR+虚拟现实+AI2026立钻哥哥++==)
- 【猿哥学二建】第一周 考试准备
- 打造爆款产品需要注意什么
- 兔子机器人Blossom成为萌宠,软体机器人将会是设计新方向?
- IBM安全报告:合法网站成企业IT安全的罩门