最近科研中用到LM算法,在研究,其基本原理可以再如下网页查看到:
http://blog.sina.com.cn/s/blog_4a8e595e01014tvb.html
讲解比较深刻。基本能看懂。用matlab也能跑通。下面的代码是这个博客上的,也是别的很多地方的LM范例。
Levenberg-Marquardt快速入门教程(荐)
例子程序(MATLAB源程序)
本程序不到100行,实现了求雅克比矩阵的解析解,Levenberg-Marquardt最优化迭代,演示了如何求解拟合问题。采用萧树铁主编的《数学试验》(第二版)(高等教育出版社)中p190例2(血药浓度)来演示。在MATLAB中可直接运行得到最优解。
*************************************************************************
% 计算函数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 14.10 12.89 9.32 7.45 5.24 3.01];
% 2. LM算法
% 初始猜测s
a0=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: 迭代
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))];
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));
% 计算步长dp,并根据步长计算新的可能的\参数估计值
dp=inv(H_lm)*(J'*d(:));
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 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
************************************************************
然后自己用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');
hold on
plot(data_1,a(1)*exp(-a(2)*data_1));
hold off
其中EFunctionEst的代码如下:
function E = EFunctionEst(a, x,y )
%这是一个测试文件用于测试 lsqnonlin
% Detailed explanation goes here
x=x(:);
y=y(:);
Y=a(1)*exp(-a(2)*x);
E=y-Y;
end
迭代13次就可以完成了。结果如下:
First-Order Norm of
Iteration Func-count Residual optimality Lambda step
0 3 1436 20.8 0.01
1 9 1277.57 64.1 10 2.1744
2 13 1125.58 166 100 0.687293
3 17 1058.07 227 1000 0.167484
4 20 932.671 340 100 0.216115
5 24 827.557 327 1000 0.225789
6 27 810.124 63.5 100 0.0731397
7 30 737.768 60.2 10 0.598774
8 33 402.19 341 1 3.89167
9 36 51.9641 265 0.1 8.27918
10 39 1.15523 6.36 0.01 3.76119
11 42 1.0659 0.0597 0.001 0.208848
12 45 1.06589 0.00305 0.0001 0.00278028
13 48 1.06589 6.09e-005 1e-005 3.51923e-005