一、思想

高斯牛顿法的对象是最小二乘法。

采用一定的方法对Hession 矩阵进行近似，这样的话可以减少计算量，只需要计算一阶偏导数得到雅可比矩阵即可。

               minF(x)=|| f(x)||^2  

那么x在xk处的增量Δxk出的最小二乘法为 

    minF(xk+Δxk)=∣∣f(xk​+Δxk​)∣∣^2

同理：

     F(xk+Δxk)=≈∣∣f(xk​)+J(xk​)TΔxk​∣∣^2

J(xk​) 是一阶导数，也就是雅可比矩阵，xk是一个已知数据，只有Δxk​是未知数，那么我们要计算的就是这个Δxk​。

二、算法流程

1、确定函数的模型,设定迭代精度m

2、设定一个初始值  y=A*exp(Bx),  已知x y 数据，所以拟合的就是参数  A 和B，那么就要给A  和B 设定一个初始值。

3、对与k次迭代，计算J(xk)雅可比矩阵 和 误差r 

 y'=A*exp(Bx),  r=y-y'

4、计算 Δxk=inv(J(xk)'J(xk)T)*(-J(xk)*r)

5、计算 norm(Δxk) <m？Ok：x(k+1)=x(x)+Δxk    ---》返回3

eg1 :

clear ;
clc;
iterator_num=500;
segma=1e-8;% 迭代精度
% y=a*exp(bx)+c 函数模型
% 设置其实参数值
a=6;
b=0.3;
c=0.1;
act_x=[1;2;3;4;5;6;7;8];
act_y=[8.3;11.0;14.7;19.7;26.7;35.2;44.4;55.9];
% 开始迭代
for i=1:iterator_num
    y=a*exp(b*act_x);
    %计算误差 r=act_y-y
    r=act_y-y;
    % 计算偏导数a 的雅可比矩阵
    Jacobian_a=exp(b*act_x);
    % 计算偏导数b 的雅可比矩阵
    Jacobian_b=a*exp(b*act_x).*act_x;
  
    % 联立雅可比方阵
    Jf=[Jacobian_a,Jacobian_b];
    % =====
    delta_abc=inv(Jf'*Jf)*Jf'*r;
    if norm(delta_abc)<segma
        break;
    end
    a_iterator=a+delta_abc(1);
    b_iterator=b+delta_abc(2);
    % 判断精度
    if norm(r)<segma
        break;
    end
    a=a_iterator;
    b=b_iterator;
end
it_y=a*exp(b*act_x);
plot(act_x,act_y,'.',act_x,it_y,'-');
xlim([0 10]);
ylim([0 70]);
legend('act','fit','Location','southoutside','Orientation','horizontal')

eg2；

Matlab Gauss 拟合_Σίσυφος1900的博客-CSDN博客

从上面的链接可以得到一些数据x y

 函数的模型是：y(i)=a*exp(-(x(i)-b).^2/c.^2)+0.1*rand(1);

补充：求导数

>> syms  act_x
>> y=a*exp(-(act_x-b).^2/c.^2)
 
y =
 
a*exp(-(act_x - b)^2/c^2)
 
>> diff_a=diff(y,a)
 
diff_a =
 
exp(-(act_x - b)^2/c^2)
 
>> diff_b=diff(y,b)
 
diff_b =
 
(a*exp(-(act_x - b)^2/c^2)*(2*act_x - 2*b))/c^2
 
>> diff_c=diff(y,c)
 
diff_c =
 
(2*a*exp(-(act_x - b)^2/c^2)*(act_x - b)^2)/c^3
 
>>