levenberg marquardt method

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
);