Re: Marquadt-Levenberg method. Message #3 Posted by Les Wright on 26 June 2006, 10:55 p.m., in response to message #1 by Ron Allen
Here is a routine I put together some years ago in Maple (an earlier version), with some inspiration from the NR routine and some other assistance and inspiration from a similar routine of a fellow named Dave Holmgren.
If you can "read" Maple you may find it informative. Seems to work, though it is slow. There is faster curve fitting software out there!
Les
LevMarq:= proc(f::name=algebraic,indvar::name,firstguess::list(name=numeric),xlist::list(numeric),ylist::list(numeric),tol)
> local x,params,reso,npar,npts,drvs,eqs,i,j, w, A, alpha, dy, beta, lambda,p,dynew,pnew,ssold,ssnew,niter,pnlist:
> w:=interface(warnlevel):interface(warnlevel=0):with(linalg):interface(warnlevel=w):
> reso:=10^10;
> x:=indvar;
> params:=map(lhs,firstguess);
> npar:=nops(params);
> npts:=nops(xlist);
> drvs:=grad(rhs(f),vector(params));
> A:=matrix(npts,npar,0);
> alpha:=matrix(npar,npar,0);
> dy:=vector(npts,0);
> dynew:=vector(npts,0);
> beta:=vector(npts,0);
> pnew:=vector(npar,0);
> lambda:=.001;
> p:=map(rhs,firstguess);
> eqs:=zip((s,t)->s=t,params,p);
>
> for niter from 1 to 200 do
> for i from 1 to npts do
> for j from 1 to npar do
> A[i,j]:=evalf(subs(eqs,subs(x=xlist[i],drvs[j])));
> od
> od;
> alpha:=evalm(transpose(A)&*A);
> for j from 1 to npar do alpha[j,j]:=(1+lambda)*alpha[j,j] od;
> for i from 1 to npts do
> dy[i]:=evalf(ylist[i]-subs(eqs,subs(x=xlist[i],rhs(f))))
> od:
> beta:=evalm(transpose(A)&*dy);
> pnew:=evalm(vector(p)+linsolve(alpha,beta)):
> pnlist:=[seq(pnew[i],i=1..npar)];
> eqs:=zip((s,t)->s=t,params,pnlist);
> for i from 1 to npts do
> dynew[i]:=evalf(ylist[i]-subs(eqs,subs(x=xlist[i],rhs(f))))
> od;
> ssold:=dotprod(dy,dy);
> ssnew:=dotprod(dynew,dynew);
> if ssold<=ssnew then
> lambda:=10*lambda;
> else
> lambda:=.1*lambda;
> p:=pnlist;
> fi;
> print(sprintf(`Iteration no. %a`,niter));
> print(`Old SSR`=ssold):
> print(`New SSR`=ssnew):
> eqs:=zip((s,t)->s=t,params,p);
> print(eqs);
> print(``);
> if abs(ssnew-reso) <= tol then break fi;
> reso:=ssnew;
> od:
> lhs(f)=subs(eqs,rhs(f)):
> end:
|