levenberg marquardt method
I made levenberg marquardt method script on maxima.
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);
output
[2.566454115081433, 2.566454115081433] - .9128770877033535 1.e-5 [2.342857486555295, 2.342857486555295] - .9996442697606387 1.0000000000000002e-6 [2.356197647528097, 2.356197647528097] - .9999999999800625 1.0000000000000002e-7 (%o16) [2.356197647528097, 2.356197647528097]
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 );