大家晚上好呀，这是我《科学计算与仿真》的大作业，上传到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)')
break
end
X=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》，需要原稿的可以私信我。
7617
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
