LM算法程序

最近科研中用到LM算法，在研究，其基本原理可以再如下网页查看到：讲解比较深刻。基本能看懂。用matlab也能跑通。下面的代码是这个博客上的，也是别的很多地方的LM范例。Levenberg-Marquardt快速入门教程（荐）例子程序（MATLAB源程序）本程序不到100行，实现了求雅克比矩阵的解析解，Levenberg-Marquardt最优化迭代，演示了如何求解拟合问题。采用萧树铁主编的《数学试验》（第二版）（高等教育出版社）中p190例2（血药浓度）来演示。在MATLAB中可直接运行得到最优解。*************************************************************************%计算函数f的雅克比矩阵，是解析式symsabyxreal;f=a*exp(-b*x);Jsym=jacobian(f,[ab])%拟合用数据。参见《数学试验》，p190，例2data_1=[0.250.511.523468];obs_1=[19.2118.1515.3614.1012.899.327.455.243.01];%2.LM算法%初始猜测sa0=10;b0=0.5;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:迭代forit=1:n_itersifupdateJ==1%根据当前估计值，计算雅克比矩阵J=zeros(Ndata,Nparams);fori=1:length(data_1)J(i,:)=[exp(-b_est*data_1(i))-a_est*data_1(i)*exp(-b_est*data_1(i))];end%根据当前参数，得到函数值y_est=a_est*exp(-b_est*data_1);%计算误差d=obs_1-y_est;%计算（拟）海塞矩阵H=J'*J;%若是第一次迭代，计算误差ifit==1e=dot(d,d);endend%根据阻尼系数lamda混合得到H矩阵H_lm=H+(lamda*eye(Nparams,Nparams));%计算步长dp，并根据步长计算新的可能的\参数估计值dp=inv(H_lm)*(J'*d(:));g=J'*d(:);a_lm=a_est+dp(1);b_lm=b_est+dp(2);%计算新的可能估计值对应的y和计算残差ey_est_lm=a_lm*exp(-b_lm*data_1);d_lm=obs_1-y_est_lm;e_lm=dot(d_lm,d_lm);%根据误差，决定如何更新参数和阻尼系数ife_lmlamda=lamda/10;a_est=a_lm;b_est=b_lm;e=e_lm;disp(e);updateJ=1;elseupdateJ=0;lamda=lamda*10;endend%显示优化的结果a_estb_est************************************************************然后自己用matlab工具中的优化函数：lsqnonlin试了一下，也能跑通，效果不错：代码如下：（下面的代码直接粘贴在上面的后面就可以了。）a0=[1,2];options=optimset('Algorithm','Levenberg-Marquardt',...'Display','iter');a=lsqnonlin(@EFunctionEst,a0,[],[],options,data_1,obs_1)plot(data_1,obs_1,'o');holdonplot(data_1,a(1)*exp(-a(2)*data_1));holdoff其中EFunctionEst的代码如下：functionE=EFunctionEst(a,x,y)%这是一个测试文件用于测试lsqnonlin%Detailedexplanationgoesherex=x(:);y=y(:);Y=a(1)*exp(-a(2)*x);E=y-Y;end迭代13次就可以完成了。结果如下：First-OrderNormofIterationFunc-countResidualoptimalityLambdastep03143620.80.01191277.5764.1102.17442131125.581661000.6872933171058.0722710000.167484420932.6713401000.216115524827.55732710000.225789627810.12463.51000.0731397730737.76860.2100.598774833402.1934113.8916793651.96412650.18.2791810391.155236.360.013.7611911421.06590.05970.0010.20884812451.065890.003050.00010.0027802813481.065896.09e-0051e-0053.51923e-005