quad_min

(tmp_for_tar/quad_min.doc/quad_min.m)

Function Synopsis

[xbest,vbest,niter,hbest,ybest]=quad_min(func,d2func,xinit,...)

Help text

 [x,v,niter,h,y] = quad_min(f,d2f,init,...) - Newton-like minimization

               Minimized function.          |Returns Value, 1st and 
                                  \         |2nd derivatives of func
                                   \       /      .Starting point
                                    |     |      /
                                    v     v     v
       [x,v,niter,h,y] = quad_min(func,d2func,xinit,options)
        ^ ^  ^    ^ ^
        | |  |    |  \ 
        | |  |     \  See below
        | |   \     Pseudo-inverse of Hessian
        |  \   Number of iterations
        |   Best value of "func"
 Best "x" found

 Not-quite like Newton iterations to find minimum evenberg-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.

Cross-Reference Information

sayif         tmp_for_tar/quad_min.doc/sayif.m

Listing of function file tmp_for_tar/quad_min.doc/quad_min.m

## [x,v,niter,h,y] = quad_min(f,d2f,init,...) - Newton-like minimization
## 
##               Minimized function.          |Returns Value, 1st and 
##                                  \         |2nd derivatives of func
##                                   \       /      .Starting point
##                                    |     |      /
##                                    v     v     v
##       [x,v,niter,h,y] = quad_min(func,d2func,xinit,options)
##        ^ ^  ^    ^ ^
##        | |  |    |  \ 
##        | |  |     \  See below
##        | |   \     Pseudo-inverse of Hessian
##        |  \   Number of iterations
##        |   Best value of "func"
## Best "x" found
##
## Not-quite like Newton iterations to find minimum evenberg-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>
##
function [xbest,vbest,niter,hbest,ybest]=quad_min(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 ;

code = "" ;
## 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 - 3 ;
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 ,"quad_min : Read option : %s.\n",tmp) ;

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

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

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

if !isstr (d2func) || !isstr (func),
  printf("quad_min : I need func and d2func to be strings!\n");
  keyboard
end

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


exist_y = exist ("y","var");
if isnan (vinit) , 
  if exist_y, v = feval (func, x, y); ybest = y;
  else        v = feval (func, x);  
  end
else          v = vinit ; 
end
## keyboard

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

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

while niter++ <= maxouter ,

  if exist_y, [v,d,h] = feval(d2func,x,y) ;
  else        [v,d,h] = feval(d2func,x) ;
  end

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

  ## printf("vold=%8.3g , v=%8.3g *****************************\n",vold,v);
  
  sayif( verbose || report && (rem(niter,max(report,1)) == 1) ,...
	"quad_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 ,...
	  "quad_min : exiting 'cause low gradient\n");
    ## keyboard
    break			# Exit outer loop
  end

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

  ninner = 0 ;
  while ninner++ < maxinner,	# Inner loop #########################
		       
    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("quad_min : Whoa!! any(isnan(xnew)) (1)\n"); 
      keyboard
    end
    
    if exist_y, vnew = feval (func, xnew, y); 
    else        vnew = feval (func, xnew); 
    end

    neval++ ;
    if vnew<vbest ,
      dbest = dx ; vbest = vnew; xbest = xnew; 
      done_inner = 1 ;		# Will go out at next increase
      sayif( verbose,"quad_min : going down\n" );
      
    elseif done_inner == 1,	# Time to go out
      sayif( verbose,"quad_min : quitting %d th inner loop\n",niter); 
      break			# out of inner loop
    end
    wt = wt*tcoeff ;
    wn = wn*ncoeff ;
  end				# End of inner loop ##################

  if ninner >= maxinner ,	# There was a problem
    sayif(verbose, "quad_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("quad_min : Whoa!! any(isnan(xnew)) (1)\n"); 
      keyboard
    end
    if exist_y, vnew = feval (func, xnew, y) ; 
    else        vnew = feval (func, xnew) ; 
    end

    while vnew < vbest,
      vbest = vnew; xbest = xnew; 
      wn = wn*ocoeff ;
      xnew = x+wn*dbest;
      if exist_y, vnew = feval (func, xnew, y); 
      else        vnew = feval (func, xnew); 
      end
    end
  end				# End of improving along dbest

  if verbose || rem(niter,max(report,1)) == 1,
    if vold,
      sayif(verbose, "quad_min : done inner vbest=%8.5g, vbest/vold=%8.5g\n",...
	     vbest,vbest/vold);
    else
      sayif(verbose, "quad_min : done inner vbest=%8.5g, vold=0\n",...
	     vbest);
    end
  end
  if vbest < vold ,
    ## "improvement found"
    if prudent,
      if exist_y, tmpv = feval (func, xbest, y) ; 
      else        tmpv = feval (func, xbest) ; 
      end
      if abs(tmpv-vbest)>eps,
	printf("quad_min : Whoa! best value is changing\n");
	keyboard
      end
    end
    v = vbest; x = xbest;
    if ! isempty (code),
      sayif (verbose,"quad_min : Gonna eval("%s"\n",code); 
      eval (code);
      xbest = x; 
      if exist_y, ybest = y; end
      if prudent,
	tmpv = feval (func, x, y);
	if abs (tmpv-vbest) > eps,
	  printf("quad_min : Whoa! best value changes after eval(code)\n");
	  keyboard
	end
      end
    end
  end

  if vbest > mimprov*vold , 
    if verbose || report ,
      printf("quad_min : quitting, niter=%-3d v=%8.3g, ",niter,v);
      if vold, printf("v/vold=%8.3g \n",v/vold);
      else     printf("vold  =0     \n",v);
      end
    end
    ## keyboard
    break ;			# out of outer loop
  end   
end				# End of outer loop ##################

if niter >= maxouter ,
  sayif(verbose, "quad_min : outer looping forever\n");
end

if exist_y, err = feval (func, xbest, ybest); 
else        err = feval (func, xbest); 
end

if abs (err-vbest) > eps,
  printf("quad_min : Whoa!! xbest does not eval to vbest\n");
  printf("         : %8.3e - %8.3e = %8.3e != 0\n",err,vbest,err-vbest);
  keyboard
end