test_LM.zip_LEVMAR_LM_levenberg marquardt_lm.ref_marquardt

function [p,res_p]=test_LM(p0,xdata,ydata) % this m-file is used to test the simplest Levenberg-Marquardt method % Reference: % A Brief Description of the Levenberg-Marquardt Algorithm Implemented by % levmar % Ref: http://www.ics.forth.gr/~lourakis/levmar % Edited by LUXP % Date: 2011.1.17 % Initial the sample: f(x)=acos(x)+bsin(x) for given x=1,2,3,fx=fx1,fx2,fx3 % and initial test p0=[0,0]' % Minimizing ||fx-f(a,b)||^2 to compute the optimizal a,b by LM algorithm %fx=fun_LM(x,p); %given x and fx values % define the threshold parameters: kmax=100; % the largest iteration steps eps1=10^(-15); %check whether the ||g||_inf is close to 0 eps2=10^(-15); % check eps3=10^(-15); % check tao=10^(-3); factor=2; %the diagnal values of J'J is enlargered by factor %Initial values p=p0; J=Jacobi_LM(p,xdata); %Jacobi matrix of fx on parameters p A=J'*J; res_p=ydata-fun_LM(p,xdata); %residual of p: ydata-f(x,p) g=J'*res_p; %solve the linear equation: Ndp=g to find dp:delta P:change of P; % main iteration steps if norm(g,inf)<=eps1 disp('Jacobi*residual=0, the initial p is wanted!'); else eye_mu=eye(size(A)); mu=tao*max(diag(A)); for k=1:kmax A=A+mu.*eye_mu; dp=A\g; if norm(dp)<=eps2*norm(p) disp('delta P is large less than P, the p is wanted!' ); break; end p_new=p+dp; rou=(norm(res_p)^2-norm(ydata-fun_LM(p_new,xdata))^2)/(dp'*(mu*dp+g)); if rou>0 p=p_new; J=Jacobi_LM(p,xdata); A=J'*J; res_p=ydata-fun_LM(p,xdata); g=J'*res_p; if (norm(g,inf)<=eps1) || (norm(res_p)^2<=eps3) disp('norm g is less then eps1 or norm of res_p less than eps3'); break; end; mu=mu*max(1/3,1-(2*rou-1)^3); factor=2; else mu=mu*factor; factor=factor*2; end end end p=p; res_p=res_p; end function y=fun_LM(p,x) y=p(1).*cos(x)+p(2).*sin(x); end function J=Jacobi_LM(p,x) J=cos(x); J=[J,sin(x)]; end