晚上好，这是我的科学计算和模拟作业，上传到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，连续循环，直到误差小于可接受的误差阈值。
总之，在我的实际操作过程中，我发现最重要的是解决残差函数和雅可比矩阵。其他参数，如初始值，由系统产生的随机数给出，迭代次数和最大可接受阈值根据主题确定。
高斯牛顿算法分为两部分：高斯牛顿算法和牛顿算法。
我们先介绍牛顿法，牛顿法把求极值点转化为求极值点f ‘ (x)=0。我们先取初始点，然后在初始点做切线交x轴xnext，接着从xnext重复上述过程，使切线斜率接近0，最终得到最优解。但牛顿法的缺点是对初始点的选择要求高，目标函数Hessian矩阵的逆矩阵计算比较复杂，所以采用高斯牛顿法简化Hessian计算矩阵。
其中Jr雅可比矩阵是残差函数。
则
则
我们可以求出Δx接着由xnext=x0 Δx，带入对应Δx，连续循环，直到误差小于可接受的误差阈值。
我先用熟悉的c 写程序，发现自己对c 的eigen库和opencv库的了解还是有欠缺，所以我打算先用matlab因为matlab矩阵操作相对简单，比如Jacobian只需输入矩阵J=jacobian（[f1;f2],[x1;x2])可以得到结果，避免繁琐的公式。
如果你想找到c 参考可传输至：
(101条消息) 非线性最小二乘的解释方法_陈建驱的博客-CSDN博客_非线性最小二乘法估计参数的步骤
我的想法是：
1.根据已知数据的点图来判断该模型近似于哪一种分布。
假设模型为此分布，并列出方程。
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对于各参数的偏导：
与此同时，在整理了所有的报错之后，我学会了对矩阵的定义，对参数的定义等等。
虽然对于解题是失败的，但是我在失败的项目中学到了不少东西，以至于最后我决定举一个稍微简单些的例子来完成这次作业。
例子以及解法我参考了（小学期太短了，手撸上述未解开的题目花了太长时间，理解了方法但是没时间再手撸新题目了，所以我就多写点注释以表歉意）：
(101条消息) 非线性最小二乘问题的高斯-牛顿算法_work_harderer的博客-CSDN博客
没想到居然是我的直系学长哈哈哈。
假设某个初创科技企业，在一段时间内年营业额增长曲线类似于指数分布，收集该企业近十年的数据，绘制成以下表格
| 年份 |
2014 |
2015 |
2016 |
2017 |
2018 |
2019 |
2020 |
2021 |
| 年营业额（单位：万元） |
80 |
112 |
180 |
300 |
480 |
700 |
1200 |
2000 |
假定该企业营业额函数为：
运行代码作图得到最佳参数x1与x2为：
曲线为：
[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（吴恩达）
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
该例截自书《2004 Statistical Tools for Nonlinear Regression》，需要原稿的可以私信我。