The Museum of HP Calculators

HP Forum Archive 16

 Marquadt-Levenberg method.Message #1 Posted by Ron Allen on 26 June 2006, 2:16 a.m. Would someone be interested in sharing some knowledge about the Marquate-Levenberg method of optimization? Source for reading or source code for the 48gx or successor. I am probably among the first to write RPN for Newton-Raphson for the 48 and earler TI 59, etc., so I umderstand the subject reasonably well, just haven't studied the techniques involved. I use mathcad version 13 in my desktop which uses the method imbedded in machine language. so does excel for that matter. Anyway, I will appreciate any information you want to share. As far as I can tell from using it, the method is suited for simultaneous solutions of complex non-linear models. Ron

 Re: Marquadt-Levenberg method.Message #2 Posted by Thomas Radtke on 26 June 2006, 5:54 a.m.,in response to message #1 by Ron Allen Here is everything you need: LM should be covered in 15.5 don't remember exactly. Good luck, Thomas

 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: ```

 Re: Marquadt-Levenberg method.Message #4 Posted by Ron on 27 June 2006, 12:56 a.m.,in response to message #3 by Les Wright Thanks, Les. I went with MathCad which borrows much from Maple almost a subset, but it was the difficulty of using Maple at the time for me (DOS). I think it was about Maple 5 that I found it too cumbersome for my needs. Thank you again for the code. I can learn a lot just attempting to decipher the code in a structured program. This one looks almost more complex than the problems to be solved, but I know that a lot of the extra code is for general application. Ron

 Re: Marquadt-Levenberg method.Message #5 Posted by Les Wright on 27 June 2006, 12:25 p.m.,in response to message #4 by Ron I agree it looks complex but fortunately a CAS like Maple is well suited to coding the algorithm since one can use the built in calculus functions like grad(). Linear algebra and matrix things are built in and come naturally--things like matrix inversion, solution of a linear system, vector products. I am not much of a programmer, so I was easily bewildered at first by the NR C and Pascal code. Maple removes the burden of creating matrix and vector handling routines from scratch. I look forward to getting a look a Namir's stuff once I can get the link to work :) Les

 Re: Marquadt-Levenberg method.Message #6 Posted by Namir on 27 June 2006, 8:36 a.m.,in response to message #1 by Ron Allen Ron, If you got to my web site you can look at various implementations of the Marquardt algorithm (MATLAB, True BASIC, VB .Net, and Visual C# .Net). Namir

 Re: Marquadt-Levenberg method.Message #7 Posted by bennyboy6 on 27 June 2006, 10:15 p.m.,in response to message #1 by Ron Allen I've written a Mathematica program for it too; easy with that. Numerical recipes is good, I see it was linked. A similar method is the Newton method. I don't know of any published Levenberg-Marquart or Newton calculator methods, but one gentleman published HP-41C and TI-58/59 versions of a nonlinear regression using Gauss-Newton methods. See the following: "Nonlinear regression on a pocket calculator", Brian W. Clare, Chemical Engineering, August 23, 1982 It may be of some help? It's a neat article in any case.

 Re: Marquadt-Levenberg method.Message #8 Posted by bennyboy6 on 27 June 2006, 10:37 p.m.,in response to message #7 by bennyboy6 Incidentally this article is found in Volume II of "Calculator Programs for Chemical Engineers" I may have found it at amazon.com But I don't know if those sellers are selling volume I or II. You should ask first. Volume II is the yellow one. Lots of HP calculator programs in these books! Edited: 27 June 2006, 10:45 p.m.

 Re: Marquadt-Levenberg method.Message #9 Posted by Ron Allen on 28 June 2006, 12:15 p.m.,in response to message #7 by bennyboy6 Thanks, Brian. Yes, I am familiar with Newton-Raphson. So neat and compact, especially if it applies and you have a closed form for the first derivative. Not bad otherwise even. I use my laptop and MathCad when I can. One optimizing routine that MathCad uses is labeled "GenFit." This curve-fitting routine gives the user a choice of different methods, like a least squares, M.L.. etc. I have used a mini/max routine routine in Excel to handle an entire spreadsheet making the problem pretty complex. Being naturally curios, I wanted to read up on the logic behind the method and answer questions like, "When it requires derivatives, etc., how does it use them. I get some help in the field with the 48/49 calculators performing what would have seemed like miracles in the 60's, just curios I guess. There are so many functions and subroutines built in that we no longer need the fundamentals to produce results. Kind of sad in a way, when there are still a few of us able to calculate an accurate mortgage payment on a slide rule with adequate log scales. Thanks again, Ron

 Re: Marquadt-Levenberg method.Message #10 Posted by bennyboy6 on 28 June 2006, 1:09 p.m.,in response to message #9 by Ron Allen This is not Newton-Raphson. This is Gauss-Newton. Newton-Raphson is a root-finding method. Gauss-Newton is, to quote wikipedia, a nonlinear least squares method. Marquardt method is a combination of steepest descent and Gauss-Newon methods.

 Re: Marquadt-Levenberg method.Message #11 Posted by Ron Allen on 28 June 2006, 11:51 p.m.,in response to message #10 by bennyboy6 I do appreciate your interest, Bennyboy6, but I don't believe that I showed any confusion about Newton. One of my responses was directed to the person who made a comment about it as an alternative. I tried to set the record straight for all those who came behind the first to respond. Guess I missed something, for I am very familiar with the very compact Newton-Raphson, usually applicable to a single equation as one of (as you pointed out) many root finders available to devices with few resources. I wanted a description of M&L for personal reasons and more than one kind person has answered, some with coded subroutines. I appreciate the willingness for many to help without the self-serving side comments we tend to suffer from less courteous forums. As a very qualified professor and mentor reminded me more than once, "There are few foolish questions, but many foolish answers." Ron

 Re: Marquadt-Levenberg method.Message #12 Posted by bennyboy6 on 29 June 2006, 11:19 a.m.,in response to message #11 by Ron Allen No problem Ron. If you live near a university library, check out the Optimization of Chemical Processes book I mentioned. It really does have a good explanation and algorithm, though the above mentioned Numerical Recipes is also a great reference.

 Re: Marquadt-Levenberg method.Message #13 Posted by Ron Allen on 15 July 2006, 10:17 p.m.,in response to message #12 by bennyboy6 Thanks again for all your assistance, Bennyboy. I use Mathcad version 13 which has internal solver routines internally planted for various functions including some curve fitting routines and differential equations. They seem to work well and I was wondering about the logic. You mentioned earlier that higher derivatives are used to employ steeper descent. The "GENFIT" function, for example offers choices to the user for M&L, or an alternative "Least Squares," raising my curiosity as to what advantages might I be passing up accepting one over the other. From your responses it sounds like the more unknowns and or parameters involved and or the expectation of asymptopes in the vicinity of solutions, the more I should lean towards M&L as the choice. On the other hand, I should lean towards other alternatives for less complex solutions to save time. Ron

 Re: Marquadt-Levenberg method.Message #14 Posted by Les Wright on 28 June 2006, 1:12 p.m.,in response to message #9 by Ron Allen I would be very interested if you could share the Excel spreadsheet that does minimax fitting. The email given is a real one, since I quite foolishly like tempting the spambots. Besides, everyone in this forum is pretty harmless, it seems, and I like to make it easy to be reached :) Les