LM拟合算法

最新推荐文章于 2024-09-03 23:15:24 发布

weixin_30383279

最新推荐文章于 2024-09-03 23:15:24 发布

阅读量1.9k

点赞数

文章标签： c#

原文链接：http://www.cnblogs.com/ruo-li-suo-yi/p/8619150.html

版权

本文详细介绍了Levenberg-Marquardt(LM)算法，包括其在函数y=a*e.^(-b*x)和y=a*x.^b+c形式的拟合应用，强调了参数初始化的重要性。LM算法通过线性近似转化为线性最小二乘问题，适用于快速收敛。文章还提及了nlinfit拟合、最小二乘法和信赖域法的概念，并指出在复杂函数和大量参数情况下，可能需要考虑其他算法。此外，文章讨论了数值求导、偏导数计算以及LM算法的实现步骤，推荐在条件允许时手动计算偏导数以获取更高精度。最后，提到了程序实现和相关资源链接。

摘要由CSDN通过智能技术生成

一、 Levenberg-Marquardt算法

（1）y=a*e.^(-b*x)形式拟合

clear all
% 计算函数f的雅克比矩阵，是解析式
syms a b y x real;
f=a*exp(-b*x);
Jsym=jacobian(f,[a b]);
% 拟合用数据。参见《数学试验》，p190，例2
% data_1=[0.25 0.5 1 1.5 2 3 4 6 8];
% obs_1=[19.21 18.15 15.36 18.10 12.89 9.32 7.45 5.24 3.01];
load fla
% data_1=[0.25 0.5 1 1.5 2 3 4 6 8];
% obs_1=[19.21 18.15 15.36 18.10 12.89 9.32 7.45 5.24 3.01];
data_1=1:26;
obs_1=fla3_1;%y轴的值

% 2. LM算法
% 初始猜测s
a0=1e+09; b0=0;
y_init = a0*exp(-b0*data_1);
% 数据个数
Ndata=length(obs_1);
% 参数维数
Nparams=2;
% 迭代最大次数
n_iters=50;
% LM算法的阻尼系数初值
lamda=0.01;
% step1: 变量赋值
updateJ=1;
a_est=a0;
b_est=b0;
%% 
% step2: 迭代
for it=1:n_iters %迭代次数
   if updateJ==1
       % 根据当前估计值，计算雅克比矩阵
       J=zeros(Ndata,Nparams);%数据的维度，系数的维度
       for i=1:length(data_1)
           J(i,:)=[exp(-b_est*data_1(i)) -a_est*data_1(i)*exp(-b_est*data_1(i))];%雅克比矩阵的生成
                                                                                 % 对系数a和b求导
       end
       % 根据当前参数，得到函数值
       y_est = a_est*exp(-b_est*data_1);
       % 计算误差
       d=obs_1-y_est;
       % 计算（拟）海塞矩阵
       H=J'*J;
       % 若是第一次迭代，计算误差
       if it==1
           e=dot(d,d);%向量点乘
       end
   end
   % 根据阻尼系数lamda混合得到H矩阵
   H_lm=H+(lamda*eye(Nparams,Nparams)); % J'*J+lamda*I
   % 计算步长dp，并根据步长计算新的可能的参数估计值
   dp=inv(H_lm)*(J'*d(:)); %deltp=[J'*J+lamda*I]^(-1)*(J'*εk)矩阵求逆
   g = J'*d(:);
   a_lm=a_est+dp(1);
   b_lm=b_est+dp(2);
   % 计算新的可能估计值对应的y和计算残差e
   y_est_lm = a_lm*exp(-b_lm*data_1);
   d_lm=obs_1-y_est_lm;
   e_lm=dot(d_lm,d_lm);
   % 根据误差，决定如何更新参数和阻尼系数
   if e_lm<e
       lamda=lamda/10;
       a_est=a_lm;
       b_est=b_lm;
       e=e_lm;
       disp(e);
       updateJ=1;
   else
       updateJ=0;
       lamda=lamda*10;
   end
end
%% 
%显示优化的结果
a_est
b_est
%从图形上观察拟合程度
plot(data_1,obs_1,'*')
hold on
x=1:54;
y=a_est*exp(-b_est*x);
plot(x,y,'o-')

（2）y=a*x.^b+c形式

注意公式变化后，求导的变化

clear all
% 计算函数f的雅克比矩阵，是解析式
syms a b c x ;
f=a*x.^b+c;
Jsym=jacobian(f,[a b c]);
% 拟合用数据。参见《数学试验》，p190，例2
% data_1=[0.25 0.5 1 1.5 2 3 4 6 8];
% obs_1=[19.21 18.15 15.36 18.10 12.89 9.32 7.45 5.24 3.01];
load fla
% data_1=[0.25 0.5 1 1.5 2 3 4 6 8];
% obs_1=[19.21 18.15 15.36 18.10 12.89 9.32 7.45 5.24 3.01];
data_1=1:length(fla3_1);
obs_1=fla3_1;%y轴的值

% 2. LM算法
% 初始猜测s
a0=1e+09; b0=0;c0=1.4e+09;
y_init = a0*data_1.^b0+c0;
% 数据个数
Ndata=length(obs_1);
% 参数维数
Nparams=3;
% 迭代最大次数
n_iters=50;
% LM算法的阻尼系数初值
lamda=0.01;
% step1: 变量赋值
updateJ=1;
a_est=a0;