levmar_min

(tmp_for_tar/levmar/levmar_min.m)

Function Synopsis

[xbest,vbest,niter,hbest,y] = \
      levmar_min(u,func,d2func,xinit,...)

Help text

       [x,v,niter,h,y] = levmar_min(u,func,d2func,xinit,options)
        ^ ^  ^    ^ ^
        | |  |    | \ 
        | |  |    \  See below
        | |   \    Pseudo-inverse of Hessian
        |  \   Number of iterations
        |   Value at x
 Best "x" found

 Not-quite Levenberg-Marquardt method for minimizing the norm of
 func(u,x) starting from "xinit".

 "Not quite" because it minimizes a general function, given an
 implementation of the function, and of its 1st (gradient vector) and
 2nd (Hessian matrix) derivatives, while Levenberg-Marquart is for
 minimizing a sum-of-squares, given a (vector-valued) function and its
 1st derivatives.

 Somewhat like Levenberg-Marquardt in the way of avoiding the pitfall
 of Newton methods that are not garanteed to be non-increasing.

 So it's some sort of bastard optimization function ...

 Arguments :

 - xinit is a N-by-1 column vector. 

 - "func" is the cost function : v = func(u,x), or v = func(u,x,y) if
   an extra argument "y" is specified amongst the options.

 - "d2func" returns the cost (1x1), its differential (1xN) and the
   *pseudo-inverse* of its second differential (NxN) :

    [v,dv,d2v]=d2func(u,x[,y]). "y" is again optional.

 Options :

 "vinit",v     : Value of the function at xinit. Saves one evaluation.

 "maxouter", m : Max number of outer loop iterations (default : 500)

 "maxinner", m : Max number of inner loop iterations (default :  30)

 "mimprov", w  : Stop if best_value over new_best_value is above w,
                 which should be 0<w<1.              (default:0.999)

 "mingrad", w  : Stop if gradient has 2-norm smaller than w. 
                 (default : eps) 

 "y",       y  : Extra argument to func,d2func.

 "code",code   : The passed string will be evaluated after each outer
                 loop that produced some (any) improvement. "x", the
                 best parameter found, "v" the best value and
                 eventually "y" are visible from the code.

 Last modified: March 2000

Cross-Reference Information

sayif     tmp_for_tar/levmar/sayif.m

Listing of function file tmp_for_tar/levmar/levmar_min.m

##       [x,v,niter,h,y] = levmar_min(u,func,d2func,xinit,options)
##        ^ ^  ^    ^ ^
##        | |  |    | \ 
##        | |  |    \  See below
##        | |   \    Pseudo-inverse of Hessian
##        |  \   Number of iterations
##        |   Value at x
## Best "x" found
##
## Not-quite Levenberg-Marquardt method for minimizing the norm of
## func(u,x) starting from "xinit".
##
## "Not quite" because it minimizes a general function, given an
## implementation of the function, and of its 1st (gradient vector) and
## 2nd (Hessian matrix) derivatives, while Levenberg-Marquart is for
## minimizing a sum-of-squares, given a (vector-valued) function and its
## 1st derivatives.
##
## Somewhat like Levenberg-Marquardt in the way of avoiding the pitfall
## of Newton methods that are not garanteed to be non-increasing.
##
## So it's some sort of bastard optimization function ...
## 
## Arguments :
##
## - xinit is a N-by-1 column vector. 
##
## - "func" is the cost function : v = func(u,x), or v = func(u,x,y) if
##   an extra argument "y" is specified amongst the options.
##
## - "d2func" returns the cost (1x1), its differential (1xN) and the
##   *pseudo-inverse* of its second differential (NxN) :
##
##    [v,dv,d2v]=d2func(u,x[,y]). "y" is again optional.
##
## Options :
##
## "vinit",v     : Value of the function at xinit. Saves one evaluation.
##
## "maxouter", m : Max number of outer loop iterations (default : 500)
##
## "maxinner", m : Max number of inner loop iterations (default :  30)
##
## "mimprov", w  : Stop if best_value over new_best_value is above w,
##                 which should be 0<w<1.              (default:0.999)
##
## "mingrad", w  : Stop if gradient has 2-norm smaller than w. 
##                 (default : eps) 
##
## "y",       y  : Extra argument to func,d2func.
##
## "code",code   : The passed string will be evaluated after each outer
##                 loop that produced some (any) improvement. "x", the
##                 best parameter found, "v" the best value and
##                 eventually "y" are visible from the code.

## Author:        Etienne Grossmann  <etienne@isr.ist.utl.pt>
## Last modified: March 2000

function [xbest,vbest,niter,hbest,y] = \
      levmar_min(u,func,d2func,xinit,...)

## Sometimes I suspect modified sources are not recompiled. This
## static variable is displayed each time this func is called, and
## should change at (almost) each recompilation.

## static dummy_levmar = floor(100*rand());
## dummy_levmar 

# Set the defaults (not all) #########################################
vinit = nan ;			# maximum distance

maxouter = 500 ;			
maxinner = 30 ;
mimprov = 0.999 ;
mingrad = eps ;

tcoeff = 0.5 ;			# Discount on total weight
ncoeff = 0.5 ;			# Discount on weight of newton
ocoeff = 1.5 ;			# Factor for outwards searching

report = 0 ;			# Never report
verbose = 0 ;			# Be quiet
prudent = 0 ;			# Check coherence of d2func and func?

niter = 0 ;
neval = 0 ;


## d2func = dfunc = 0 ;

# ####################################################################
# Read the options ###################################################
# ####################################################################
# Options with a value
opt1 = " vinit maxouter maxinner report mimprov mingrad code y " ;
# Boolean options 
opt0 = " verbose prudent " ;

va_start() ;
nargin = nargin - 4 ;
while nargin>0 ,
  tmp = va_arg() ; nargin-- ;
  if ! isstr(tmp) ,
    break			# Last option has been read. "tmp" is
				# first arg.
  end
  if index(opt1,[" ",tmp," "]) ,
    
    tmp2 = va_arg() ; nargin-- ;
    nargin-- ;
    eval([tmp,"=tmp2;"]) ;

    sayif(verbose ,"levmar_min : Read option : %s.\n",tmp) ;

  elseif index(opt0,[" ",tmp," "]) ,
    
    eval([tmp,"=1;"]) ;
    sayif(verbose ,"levmar_min : Read boolean option : %s\n",tmp) ;

  else
    printf("levmar_min : Unknown option : %s\n",tmp) ;
    keyboard
  end
endwhile

############################## Checking ##############################
if mimprov<=0 || mimprov>=1,
  printf("levmar_min : mimprov=%8.3g is not in ]0,1[ \n",mimprov) ;
  keyboard
end

if !isstr(d2func) && !isstr(dfunc),
  printf("levmar_min : I need dfunc and d2func to be strings!\n");
  keyboard
end

x = xinit ;
N = prod(size(x)) ;
x = reshape(x,N,1) ;

if isnan(vinit) , 
  if exist("y")==1,  v = feval( func , u, x, y ) ; 
  else               v = feval( func , u, x ) ; 
  end
else                 v = vinit ; 
end
## keyboard

xbest = x ;
vold = vbest = nan ;
hbest = nan ;

sayif( verbose,"levmar_min : Initially, v=%8.3g\n",v);

while niter++ <= maxouter ,

  if exist("y")==1, [v,d,h] = feval(d2func,u,x,y) ;
  else              [v,d,h] = feval(d2func,u,x) ;
  end

  if prudent,
    if exist("y")==1, v2 = feval(func,u,x,y) ;
    else              v2 = feval(func,u,x) ;
    end
    if abs(v2-v)>0.001*sqrt(eps),
      printf("levmar_min : func and d2func disagree %8.3g\n",abs(v2-v));
      keyboard
    end
  end
  xbest = x ;
  vold = vbest = v ;
  hbest = h ;

  ## printf("vold=%8.3g , v=%8.3g *****************************\n",vold,v);
  
  sayif( verbose || report && (rem(niter,report) == 1) ,...
	"levmar_min : niter=%d, v=%8.3g\n",niter,v ); 
  wn = 1 ;			# Weight of Newton step
  wt = 1 ;			# Total weight

  if msq(d) < mingrad,
    sayif( verbose || report ,...
	  "levmar_min : exiting 'cause low gradient\n");
    ## keyboard
    break
  end

  dnewton = -h*d' ;
  done_inner = 0;		# 0=not found. 1=Ready to quit inner.

  ninner = 0 ;
  while ninner++ < maxinner,
		       
    dx = wt*(wn*dnewton - (1-wn)*d') ;
    xnew = x+dx ;

    sayif( verbose, ...
	  ["Weight : total=%8.3g, newtons's=%8.3g  vbest=%8.3g ",...
	   "Norm:Newton=%8.3g, deriv=%8.3g\n"],...
	  wt,wn,vbest,norm(wt*wn*dnewton),norm(wt*(1-wn)*d'));
    
    if any(isnan(xnew)),
      printf("levmar_min : Whoa!! any(isnan(xnew)) (1)\n"); 
      keyboard
    end
    
    if exist("y")==1, vnew = feval( func , u, xnew, y ) ; 
    else              vnew = feval( func , u, xnew ) ; 
    end
    ## vnew = feval(func,u,xnew) ;
    neval++ ;
    if vnew<vbest ,
      dbest = dx ; vbest = vnew; xbest = xnew; 
      done_inner = 1 ;		# Will go out at next increase
      sayif( verbose,"levmar_min : going down\n" );
      
    elseif done_inner == 1,	# Time to go out
      sayif( verbose,"levmar_min : quitting %d th inner loop\n",niter); 
      break			
    end
    wt = wt*tcoeff ;
    wn = wn*ncoeff ;
  end				# End of inner loop

  if ninner >= maxinner ,	# There was a problem
    printf("levmar_min : inner looping forever (vnew=%8.3g)\n",vnew);

  else				# Look for improvement along dbest
    wn = 1.5 ; 
    xnew = x+wn*dbest;
    if any(isnan(xnew)),
      printf("levmar_min : Whoa!! any(isnan(xnew)) (1)\n"); 
      keyboard
    end
    if exist("y")==1, vnew = feval( func , u, xnew, y ) ; 
    else              vnew = feval( func , u, xnew ) ; 
    end
    ## vnew = feval(func,u,xnew) ;
    while vnew < vbest,
      vbest = vnew; xbest = xnew; 
      wn = wn*ocoeff ;
      xnew = x+wn*dbest;
      if exist("y")==1, vnew = feval( func , u, xnew, y ) ; 
      else              vnew = feval( func , u, xnew ) ; 
      end
      ## vnew = feval(func,u,xnew) ;
    end
  end				# End of improving along dbest
  if verbose || rem(niter,report) == 1,...
    if vold,
      printf("levmar_min : done inner vbest=%8.5g, vbest/vold=%8.5g\n",...
	     vbest,vbest/vold);
    else
      printf("levmar_min : done inner vbest=%8.5g, vold=0\n",...
	     vbest);
    end
  end
  if vbest < vold ,
    if prudent,
      tmpv = feval(func,u,xbest,y);
      if abs(tmpv-vbest)>eps,
	printf("levmar_min : Whoa! best value is changing\n");
	keyboard
      end
    end
    v = vbest; x = xbest;
    if exist("code")==1, 
      sayif( verbose,"levmar_min : Gonna eval(\"%s\"\n",code); 
      eval(code); 
      if prudent,
	tmpv = feval(func,u,x,y);
	if abs(tmpv-vbest)>eps,
	  printf("levmar_min : Whoa! best value changes after eval(code)\n");
	  keyboard
	end
      end
    end
  end
  ## printf("vold=%8.3g , v=%8.3g, vbest=%8.3g\n",vold,v,vbest);
  ## [feval( func , u, x, y ), feval( d2func , u, x, y ) , v; ...
  ## feval( func , u, xbest, y ), feval( d2func , u, xbest, y ) , vbest  ]
  if vbest > mimprov*vold , 
    if verbose || report ,
      printf("levmar_min : quitting, niter=%-3d v=%8.3g, ",niter,v);
      if vold, printf("v/vold=%8.3g ",v/vold);
      else     printf("vold  =0     ",v);
      end
    end
    ## keyboard
    break ; 
  end   
end

if niter >= maxouter ,
  printf("levmar_min : outer looping forever\n");
end