dfp_min

(tmp_for_tar/deriv_min.doc/dfp_min.m)

Function Synopsis

[x,v,niter,h] = dfp_min(func,dfunc,x,...)

Help text

       [x,v,niter,h] = dfp_min(func,dfunc,xinit,options)

 Broyden-Fletcher-Goldfarb-Shanno Variable metric method for
 minimizing the function 'func', whose derivative is 'dfunc', starting
 from 'xinit'. 

 TODO : Ease the following restrictions

 - xinit is a N-by-1 column vector. 
 - func  takes a single (column vector) argument.
 - dfunc takes a single (column vector) argument and returns a row
   vector. 

 OPTIONS

 'vinit',v     : Value of the function at the starting xinit.
 'maxiter', m  : At most m iterations are done (default : 200) 
 'gtol',gtol   : Set the threshold for the stopping criterion  

              gtol > max { df(i)*max(|x(i)|,1)/max(v,1) | i in 1..N }

                 where x is the current minimum, v is func(x) and df
                 is dfunc(x).
                 Default value is 10*eps

  'tol', tol    : Set the threshold for the stopping criterion  

              tol > max { dx(i)/max(|x(i)|,1) | i in 1..N }

                 where  dx is the change in the x that occured in the
                 last iteration.
                 Default value is 10*eps

 Last modified: October 2000

Listing of function file tmp_for_tar/deriv_min.doc/dfp_min.m

##       [x,v,niter,h] = dfp_min(func,dfunc,xinit,options)
##
## Broyden-Fletcher-Goldfarb-Shanno Variable metric method for
## minimizing the function 'func', whose derivative is 'dfunc', starting
## from 'xinit'. 
##
## TODO : Ease the following restrictions
## 
## - xinit is a N-by-1 column vector. 
## - func  takes a single (column vector) argument.
## - dfunc takes a single (column vector) argument and returns a row
##   vector. 
## 
## OPTIONS
## 
## 'vinit',v     : Value of the function at the starting xinit.
## 'maxiter', m  : At most m iterations are done (default : 200) 
## 'gtol',gtol   : Set the threshold for the stopping criterion  
##
##              gtol > max { df(i)*max(|x(i)|,1)/max(v,1) | i in 1..N }
##
##                 where x is the current minimum, v is func(x) and df
##                 is dfunc(x).
##                 Default value is 10*eps
##
##  'tol', tol    : Set the threshold for the stopping criterion  
##
##              tol > max { dx(i)/max(|x(i)|,1) | i in 1..N }
##
##                 where  dx is the change in the x that occured in the
##                 last iteration.
##                 Default value is 10*eps
##
##

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

function [x,v,niter,h] = dfp_min(func,dfunc,x,...)

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

report = 0 ;			# Never report
verbose = 0 ;			# Not even at end
maxiter = 500 ;			

niter = 0 ;

tol = 3*eps ;			# too small a step
alf = 1e-1 ;

gtol = 3*eps ; ## 1e-12 ;

maxstep = 100 ;
ok = 1 ;

# ####################################################################
# Read the options ###################################################
# ####################################################################
# Options with a value
opt1 = " vinit maxiter gtol tol " ;
# Boolean options 
opt0 = " verbose " ;

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;"]) ;
    if verbose ,
      printf("dfp_min : Read option : %s. Value is :\n",tmp) ;
      tmp2 
    end

  elseif index(opt0,[" ",tmp," "]) ,
    
    eval([tmp,"=1;"]) ;
    if verbose ,
      printf("dfp_min : Read boolean option : %s\n",tmp) ;
    end
  else
    printf("dfp_min : Unknown option : %s\n",tmp) ;
    keyboard
  end
endwhile


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

if isnan(vinit) , v = feval( func , x ) ; 
else              v = vinit ; 
end
## tol
## gtol
df = feval( dfunc, x ) ;
h = eye(N) ;
xi = -df' ;
maxstep = maxstep*max(sum(df(:).^2),N) ;

while niter++ <= maxiter ,
  xold = x ;
  xiold = xi ;
  ## keyboard
  [xnew,v,ok] = linesearch(df',xi,func,x,"maxstep",maxstep,"vinit",v ) ;
  xi = xnew - x ;
  x = xnew ;
  ## xi' xnew'
  test1 = max( abs(xi(:)) ./ max( 1,abs(x(:)) ) ) ;
  if test1 < tol,
    if verbose, printf("test1 = %8.3g, tol=%8.3g\n",test1,tol) ; end
    ## keyboard
    return ;
  end
  dfold = df2 = df ;
  df = feval(dfunc,x) ;
  if verbose, printf("Norm df=%8.3g\n",norm(df)); end
  den = max( v, 1 );
  test2 = max( abs(df(:)) .* max( abs(x(:)),1 ) ) / den ; ## max(max(v(:)),1) ;
  if test2 < gtol, 
    if verbose, printf("test2 = %8.3g, gtol=%8.3g\n",test1,gtol) ; end
    return ;
  end
  if verbose, printf("dfp_min : niter=%4i, v=%8.3g\n",niter,v) ; end
  df2 = df - df2 ;
  if any( abs(h-h') ),
    printf("dfp_min : h ain't symmetric, goddamit! %8.3g (1)\n",
	   max(max(abs(h-h'))));
    keyboard
  end
  hdf2 = df2*h ;
  fac = df2*xi ;
  fae = df2*hdf2' ;
  sdf2 = sum(df2.^2) ;
  sxi = sum(xi.^2) ;
  if fac^2 > eps*sdf2*sxi ,
    ## fac = 1/fac ;
    ## fad = 1/fae ;
    df2 = xi'/fac - hdf2/fae ;
    if any( abs(h-h') ),
      printf("dfp_min : h ain't symmetric, goddamit! %8.3g (3)\n",
	     max(max(abs(h-h'))));
      keyboard
    end
    hh = (xi*xi')/fac - (hdf2'*hdf2)/fae + fae*(df2'*df2) ;
    if any( abs(hh-hh') ),
      printf("dfp_min : hh ain't symmetric, goddamit! %8.3g (5)\n",
	     max(max(abs(hh-hh'))));
      keyboard
    end
    h = h+hh ;
    if any( abs(h-h') ),
      printf("dfp_min : h ain't symmetric, goddamit! %8.3g (4)\n",
	     max(max(abs(h-h'))));
      keyboard
    end
  end
  xi = - h*df' ;
  if any( abs(h-h') ),
    printf("dfp_min : h ain't symmetric, goddamit! %8.3g (2)\n",
	   max(max(abs(h-h'))));
    keyboard
  end
  ## keyboard
  
end