levenberg marquardt method

How is this different from (or better than) minpack_lsquares ?


--- On Sun, 7/24/11, Niitsuma Hirotaka <hirotaka.niitsuma at gmail.com> wrote:

From: Niitsuma Hirotaka <hirotaka.niitsuma at gmail.com>
Subject: levenberg marquardt method
To: maxima at math.utexas.edu
Date: Sunday, July 24, 2011, 7:16 AM

I made levenberg marquardt method script on maxima.

http://d.hatena.ne.jp/niitsuma/20080325/1311505968


usage
----------
/*minimize f*/
f(xx,p):=sin(apply("+",map("*",xx,p)));
df(xx,p):=transpose(p)*cos(apply("+",map("*",xx,p)));
ddf(xx,p):= -( transpose(p) .p) * sin(apply("+",map("*",xx,p)));
xinit:[0.1,0.1];
param:[1,1];
eps:0.001;
levenberg_marquardt_minimize([f,df,ddf,levenberg_marquardt_dx2x_default],xinit,
eps,param);
----------



code
----------


load(diag);
matrix2list(m):=substpart("[",m,0);
rowmatrix2list(m):=matrix2list(m)[1];
colmatrix2list(m):=rowmatrix2list(transpose(m));


mat_en_diag(m):=diag(makelist(m[k][k],k,1,length(m)));

levenberg_marquardt_dx2x_default(x,dx):=x+colmatrix2list(dx);

levenberg_marquardt_minimize(fobj,x0,eps,param):=block(
? [n,m,f,df,ddf,dx2x,fv,fv_tmp,fv0,dfv,ddfv,dx,x_tmp
? ,
? x:x0
? ,
? c:0.0001
? ,
? c_factor:10
? /*c_factor:1.2*/
? ],
? [f,df,ddf,dx2x]:fobj
? ,
? dx:df(x0,param)
? ,
? fv:f(x,param)
? ,
? fv0:fv+100
? ,
? for n:1
? /*thru 1000000000*/
? ? while???abs(fv -???fv0) >? eps
? /*
? norm2mat(dx) > eps
? */
? do (
? ? /*disp(norm2_list_or_scalar(dx)) ,*/
? ? dfv:df(x,param)
? ? ,
? ? ddfv:ddf(x,param)
? ? ,
? ? fv_tmp:fv + 1000000000
? ? ,
? ? for? m:1 thru 1000000000???while fv_tmp > fv do (
? ? ? damped_ddfv: ddfv+c*mat_en_diag(ddfv)
? ? ? ,
? ? ? dx: - invert(damped_ddfv).dfv
? ? ? ,
? ? ? disp(x_tmp:dx2x(x,dx))
? ? ? ,
? ? ? fv_tmp:f(x_tmp,param)
? ? ? ,
? ? ? if fv_tmp > fv then c:c_factor*c
? ? ? )
? ? ,
? ? x:x_tmp,
? ? fv0:fv,
? ? disp(fv:fv_tmp)
? ? ,
? ? disp(c:c/c_factor)
? ? ,
? ? x
? ? ),
? x
);
_______________________________________________
Maxima mailing list
Maxima at math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima