非线性最小二乘问题的分析与理解(附高斯牛顿法matlab代码)
大家晚上好呀,这是我《科学计算与仿真》的大作业,上传到c站上供大家参考。我在做这份作业的时候也查找了不少资料,希望以我通俗的语言大家能够大概的理解非线性最小二乘法的高斯牛顿解法以及实际应用
- 问题的提出
经调查资料、查阅文献,我认为在很多需要求解非线性方程组时,会出现多变量但是条件有限的情况。就比如说解决六元二次方程时,每一个元都代表着传感器测量到的当前值,如果我们能够测出六组数据,利用高斯消元法就能解决,但是关键难点在于,即使测量了六组数据,每一组数据都有噪声误差,测量的组数越多,最后结果的误差越大。所以这时候就需要利用最小二乘法来得到最优解。它通过最小化误差的平方和寻找数据的最佳函数匹配,可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。
此处参考:姿态篇:四.非线性最小二乘与飞控传感器校准 - 简书 (jianshu.com)
在实际情况中,我们调参数的时候就经常会遇到这样的情况。比如在调飞行器的参数时,我们需要合适的参数以减小误差。有如,飞控中存在误差有零偏误差与刻度误差。
我们设这六个误差分别为Ox、Oy、Oz、Sx、Sy、Sz ,设当前时刻真实加速度向量为α=(αx , αy , αz),由于加速度向量模值始终为1,则得到方程:
我们需要求解Ox、Oy、Oz、Sx、Sy、Sz。
- 问题的分析
由于需要求解的六个参数均为非线性函数,则无法使用常规的高斯消元法进行计算,所以我们使用非线性最小二乘以求得最优解。
经分析,我们无法求得唯一解的原因就在于获取的数据不是呈线性关系,那么便引入残差向量r,r为每一个参数的估计值与实际值的差的向量:
我们的目标是得到最接近真实值的结果,所以我们需要让r(x)尽可能地小,这时我们就将非线性最小二乘问题转化为多元函数的极值问题。
记
首先我们设定一个初始值x0,将S(x)在初始值处进行一阶泰勒展开,得到S(x)的线性近似:
两边同时求导得到:
令上式为0,可求得:
我们得到了估计误差Δx。接着由xnext=x0+Δx,带入求得对应的Δx,不断循环,直至误差小于能接受的误差阈值。
总而言之,在我的实操过程中,我发现最重要的就是求解残差函数与雅可比矩阵,其他参数,如初始值由系统产生随机数给出,迭代次数与最大可接受阈值依据题目来决定。
3高斯-牛顿算法
3.1算法的原理
高斯牛顿算法分为高斯与牛顿算法两部分,高斯牛顿法是对牛顿法的改进。
我们首先介绍牛顿法,牛顿法将求极值点转化为求f ‘ (x)=0。我们首先取初始点,接着在初始点处做切线交x轴于xnext,接着从xnext重复以上过程,使切线的斜率不断逼近于0,则最后得到最优解。但是牛顿法的缺点在于对初始点的选择要求较高,并且目标函数的Hessian矩阵的逆矩阵计算较为复杂,所以我们采用高斯牛顿法简化Hessian矩阵的计算。
其中Jr为残差函数的雅可比矩阵。
则
则
我们即能求出Δx接着由xnext=x0+Δx,带入求得对应的Δx,不断循环,直至误差小于能接受的误差阈值。
3.2算法的实现
我首先使用熟悉的c++来编写程序,发现自己对c++的eigen库和opencv库的了解还是有欠缺,所以我打算先用matlab来编写程序,因为matlab对于矩阵的操作较为简单,比如Jacobian矩阵只需要输入J=jacobian([f1;f2],[x1;x2])即可得到结果,避免了繁琐的公式。
如果你想找到c++的参考可传送至:
(101条消息) 非线性最小二乘求解方法详解_陈建驱的博客-CSDN博客_非线性最小二乘法估计参数的步骤
对于实际应用问题,我的思路是:
1.根据已知数据的点图来判断该模型近似于哪一种分布。
2.假设该模型为此分布,并列出方程。。
3.使用矩阵函数求出R与J。
4.使用高斯牛顿的公式求出Δx并继续进行迭代,直至误差小于能接受的误差阈值。
但是对于[2]中1.1节的例子里,θ与x的方程容易混淆,我花了一些时间去翻译分清两者的关系。结果由于该例难度(四个未知参数)超乎了我的水平,我选择换一个例子实现代码。但是在摸索该例的过程中我收集了几十个报错,通过查阅资料也都基本解决,但是最终让我换方向的原因是:在迭代的过程中,由于例中a b c d四个待求参数不是一个量级的,并且预估函数是e上加e,难上加难,我认为是这些原因导致得到的雅可比矩阵在第二次迭代就变成了奇异矩阵,无法再运行下去。
我将上述例子粘贴在文章的最下方,感兴趣的读者可以与我交流。
我也了解了两个求逆函数inv与pinv,虽然pinv能求奇异矩阵的伪逆矩阵,但是在之后的迭代过程中,J与x就不会再发生改变了,依旧不能得到结果。
在求J的过程中,由于需要代入的公式的值有很多,所以我设计了自动带入的循环,如下图:
虽然可以这么求,但是依据图上的结果,您可以清楚看出,矩阵保存的数据位数是有限的,所以最终还是失败了。
在边做边学的过程中,我还由求雅可比矩阵的公式中的求偏导学会了使用matlab的求多元函数的偏导函数公式diff,求出了e上加e对于各参数的偏导:
与此同时,在整理了所有的报错之后,我学会了对矩阵的定义,对参数的定义等等。
虽然对于解题是失败的,但是我在失败的项目中学到了不少东西,以至于最后我决定举一个稍微简单些的例子来完成这次作业。
4 应用实例
例子以及解法我参考了(小学期太短了,手撸上述未解开的题目花了太长时间,理解了方法但是没时间再手撸新题目了,所以我就多写点注释以表歉意):
(101条消息) 非线性最小二乘问题的高斯-牛顿算法_work_harderer的博客-CSDN博客
没想到居然是我的直系学长哈哈哈。
假设某个初创科技企业,在一段时间内年营业额增长曲线类似于指数分布,收集该企业近十年的数据,绘制成以下表格
年份 |
2014 |
2015 |
2016 |
2017 |
2018 |
2019 |
2020 |
2021 |
年营业额(单位:万元) |
80 |
112 |
180 |
300 |
480 |
700 |
1200 |
2000 |
假定该企业营业额函数为:
运行代码作图得到最佳参数x1与x2为:
曲线为:
5参考文献
[1]. Methods for Non-linear Least Squares Problems
[2]. 视觉SLAM十四讲
[3]. Multiple View Geometry in Computer Vision, Appendix 6, Iterative Estimation Methods
[4]. 机器学习.Lecture 2(吴恩达)
6附录(仅适用于特定题)
format long; %默认有效数字16位
k=0; %初始迭代次数
s=1e-3; %最小可接受的阈值
X=[1,1]; %参数矩阵
K=100; %最大迭代次数
x=1:1:8; %此x为未知数x
A=X(1); %A等于参数x1,B等于参数x2
B=X(2);
p=[80 112 180 300 480 700 1200 2000]; %实际值,参见表格数据while k<K %当k未到达最大迭代次数并未满足小于最小可接受阈值时一直迭代R=[A.*exp(B)-80;A.*exp(2.*B)-112;A.*exp(3.*B)-180;A.*exp(4.*B)-300;A.*exp(5.*B)-480;A.*exp(6.*B)-700;A.*exp(7.*B)-1200;A.*exp(8.*B)-2000];J=[exp(B),A.*exp(B);exp(2.*B),2.*A.*exp(2.*B);exp(3.*B),3.*A.*exp(3.*B);exp(4.*B),4.*A.*exp(4.*B);exp(5.*B),5.*A.*exp(5.*B);exp(6.*B),6.*A.*exp(6.*B);exp(7.*B),7.*A.*exp(7.*B);exp(8.*B),8.*A.*exp(8.*B)];delta=inv(J'*J)*J'*R if norm(delta,2)<s %返回A的最大奇异值 fprintf('经过%d次迭代:最优参数X=[%f,%f]\n',k,X(1),X(2)')breakendX=X-delta; %更新k=k+1;fprintf('第%d迭代:X=[%f,%f]\n',k,X(1),X(2));A=X(1);B=X(2);
y=A.*exp(B.*x);
subplot(3,4,k);
plot(x,p,'*');
hold on
plot(x,y);
xlabel('time/year');
ylabel('annual turnover/Ten thousand yuan');
title(['第',num2str(k),'次迭代拟合曲线图']);
end
7解不开的例子
该例截自书《2004 Statistical Tools for Nonlinear Regression》,需要原稿的可以私信我。
非线性最小二乘问题的分析与理解(附高斯牛顿法matlab代码)相关推荐
- 【综合评价分析】topsis评价 原理+完整MATLAB代码+详细注释+操作实列
[综合评价分析]topsis评价 原理+完整MATLAB代码+详细注释+操作实列 文章目录 1.TOPSIS法的原理 2.TOPSIS法案例分析 3.建立模型并求解 3.1数据预处理 3.2代码实现数 ...
- MATLAB算法实战应用案例精讲-【数模应用】小批量梯度下降(MBGD)(附Python和MATLAB代码)
目录 前言 几个相关概念 1.Batch 2. Iteration 3. Epoch 知识储备 损失函数
- 【ESN-PSO】基于PSO的回波状态网络参数分析,用于时间序列预测(Matlab代码实现)
- 现代信号处理之手机加速度传感器步态数据采集、分析(采集的数据和MATLAB代码见CSDN同名资源)
一.实验目的 通过实际数据采集.处理加深对理论知识的理解和掌握,提高学生动手能力. 二.实验原理 零漂处理.降噪 谱分析 滤波 三.实验内容与结果 3.1 数据采集 下载MATLAB APP或其它手机 ...
- 【故障诊断分析】基于 FFT轴承故障诊断matlab代码
1 简介 快速傅里叶变换( FFT) 是 1965 年 J W Cooley 和 J W Tukey 巧妙地利用 WN 因子的周期性和对称性,构造了离散傅里叶变换( DFT) 的快速算法,即快速离散傅 ...
- MATLAB算法实战应用案例精讲-【人工智能】机器视觉(概念篇)(补充篇)(附python和matlab代码实现)
目录 前言 知识储备 机器视觉成像 2. 镜头 3. 为机器视觉项目选择镜头和相机的简化流程
- MATLAB算法实战应用案例精讲-【智能优化算法】斑点鬣狗优化-SHO(附python和matlab代码)
目录 前言 算法原理 算法思想 1.1包围机制 1.2 狩猎机制 1
- 非线性最小二乘问题的高斯-牛顿算法
@非线性最小二乘问题的高斯-牛顿算法 非线性最小二乘与高斯-牛顿算法 开始做这个东西还是因为学校里的一次课程设计任务,找遍了全网好像也没有特别好用的,于是就自己写了一个.仅供参考. 首先,介绍下非线性 ...
- 层次分析法及matlab代码
数学建模算法(一) 层次分析法 The analytic hierarchy process(AHP) [清风数学建模课程笔记] 文章目录 数学建模算法(一) 层次分析法 The analytic h ...
- 谈谈自己对线性最小二乘和非线性最小二乘之间关系的理解~
最小二乘问题,是最基本.最实用.最应用广泛的的数学模型,在机器学习中更是得到了极大的应用(公式1),比如说我们的PCA,最经典的解释就是:最小化样本与其重构样本之间的误差和.采用的公式我不用写出来大家 ...
最新文章
- 502 Server dropped connection
- Java中常见的几种类型转换
- python开源项目及示例代码
- 对于指针传入函数,却最终没有改变指针的值的问题
- 使用ffmpeg裁剪和合并视频
- SAP FI模块与SD、MM的接口配置
- Kafka 安装和搭建 (一)
- 分治法在二叉树遍历中的应用(JAVA)--二叉查找树高度、前序遍历、中序遍历、后序遍
- centos 6.3 64bit 安装VMware workstation 9.1 64bit
- PyTorch 入坑七:模块与nn.Module学习
- 文献管理三剑客之Noteexpress:更新一次文献后把文献自动插一次
- 工科数学分析序言及索引(不断更新中)
- 什么是BLOB URL,为什么要使用它?
- Python游戏编程(五)Tic Tac Toe
- Android 仿微信实现语音聊天功能
- MP4文件格式简要解析---代码解析
- 桥梁通服务器物理连接成功,ZStack 实践汇|OSPF搭建与物理网络通信的“桥梁”
- 算法笔记——基数排序
- c#操作word文档之简历导出
- 【解决方案】“/usr/bin/nvcc“ is not able to compile a simple test program解决方案
热门文章
- python判断身份证号码是否合法_怎样使用 Python 判断身份证号码是否正确-阿里云开发者社区...
- wps目录怎么加一条_wps目录怎么自动生成目录?目录自动生出方法介绍
- Linux驱动(并发):02---编译乱序、执行乱序(屏障API(bm、rmb、wmb、__iormb、__iowmb))
- dtu转发虚拟服务器,DTU转发云服务器
- Pillow教程(一)
- 使用Python-OpenCV将图片批量转换为jpg格式
- 微信充错服务器,微信话费充错了怎么追回来(只需一招即可追回)
- web前端期末大作业 HTML+CSS+JavaScript仿安踏
- 台式计算机时间不能同步,台式电脑时间同步不了?一分钟就能快速解决
- 如何用Python制作词云,对1000首古诗做词云分析!