Gauss-Newton法matlab求解

sanliujiu222

已于 2023-10-07 20:41:50 修改

阅读量2.2k

点赞数

文章标签： matlab 最小二乘法

于 2022-07-14 21:18:03 首次发布

本文链接：https://blog.csdn.net/sanliujiu222/article/details/125786692

版权

目标函数 $y=Ae^{Bx}$

输入变量 $x=\begin{bmatrix} 1 &2 &3 &4 &5 &6 &7 &8 \end{bmatrix}$

目标值 $y=\begin{bmatrix} 8.3 &11.0 &14.7 &19.7 &26.7 &35.2 &44.4 &55.9 \end{bmatrix}$

简单介绍一下高斯-牛顿迭代求解非线性方程组的过程。

高斯-牛顿迭代可以说是牛顿迭代的升级版->就是将误差的二范数进行最小值的优化求解。

这个方法的原理是，利用二阶泰勒展开去近似线性回归模型，通过迭代更新系数，令残差最小。

模型函数： $f(x,\beta )$ ,目标函数: $min \, r=(f-y)^{2}$ ,此为单变量目标函数（单变量的意思就是 $f$ 是一个单一变量），但是实际问题中肯定不止单变量，会出现多变量，因此会出现： $min\, r=(f-y)^{T}(f-y)$ 。

第i个输入变量的残差可以表示为： $r_{i}=f(x_{i},\beta )-y_{i}$

又因为在目标函数 $y=Ae^{Bx}$ 中，我们需要求A,B这两个需要更新的数值，因此对函数f中的A,B分别进行偏微分运算，得到对应的雅可比矩阵 $J^{T}=\begin{bmatrix} \frac{\delta J}{\delta A}\\ \frac{\delta J}{\delta B} \end{bmatrix}$

在牛顿迭代中，牛顿迭代的公式为： $x^{k+1}=x^{k}-H^{-1}\bigtriangledown f$

其中 $H^{-1}$ 是f的二阶梯度，在高斯-牛顿法中，省略二阶项，只留下一阶项的二次方，刚好与梯度 $\bigtriangledown f$ 吻合

因此，得到最终的迭代公式： $x^{k+1}=x^{k}+(J^{T}J)^{-1}J^{T}r_{i}$

代码如下：

%object function y=A*e^(b*t)

clear
clc

deriv_step=1e-5;
iternum=100;

input=[1;2;3;4;5;6;7;8];
output=[8.3;11.0;14.7;19.7;26.7;35.2;44.4;55.9];

%gauss_newton
A=6;
B=0.3;
loop_num=100;
tt=1e-8;
for i=1:loop_num
    y=A*exp(B*input);
    r=output-y;
    JfA=exp(B*input);
    JfB=A*exp(B*input).*input;
    Jf=[JfA,JfB];
    betha=inv(Jf'*Jf)*Jf'*r;
    Anew=A+betha(1);
    Bnew=B+betha(2);
    if norm(r)<tt
        break;
    end
    A=Anew;
    B=Bnew;
end
subplot(1,2,1)
plot(input,output,'--bo',input,y,'rx');
subplot(1,2,2)
plot(input,r,'--b');
% xlim([-10 10]);
% ylim([0 100]);

得到结果如图：