conjgrad_min

(tmp_for_tar/deriv_min.doc/conjgrad_min.m)

Function Synopsis

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

Help text

 [x,v,n] = conjgrad_min(f,df,x,opt) - Conjugate Gradient minimization

 Conjugate Gradient method for minimizing the function 'f', whose
 derivative is 'df', starting from 'xi'.

 ARGUMENTS : 
 f  : string : Name of function, whose synopsis is 'val = f (x)' 
               where x is Nx1, val is 1x1.
 df : string : Name of derivative of f Synopsis is 'dv = df (x)'
               where dv is 1xN.
 x  : N      : Starting point

 RETURNED VALUES :
 x  : Nx1    : Found minimum
 v  : 1      : Value at minimum
 n  : 1      : Number of iterations

 TODO : Ease the following restrictions

 - x  is a N-by-1 column vector. 
 - f  takes a single (column vector) argument. 
 - df 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/conjgrad_min.m

## [x,v,n] = conjgrad_min(f,df,x,opt) - Conjugate Gradient minimization
##
## Conjugate Gradient method for minimizing the function 'f', whose
## derivative is 'df', starting from 'xi'.
##
## ARGUMENTS : 
## f  : string : Name of function, whose synopsis is 'val = f (x)' 
##               where x is Nx1, val is 1x1.
## df : string : Name of derivative of f Synopsis is 'dv = df (x)'
##               where dv is 1xN.
## x  : N      : Starting point
##
## RETURNED VALUES :
## x  : Nx1    : Found minimum
## v  : 1      : Value at minimum
## n  : 1      : Number of iterations
##
## TODO : Ease the following restrictions
##
## - x  is a N-by-1 column vector. 
## - f  takes a single (column vector) argument. 
## - df 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] = conjgrad_min(func,dfunc,x,...)

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

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

niter = 0 ;
maxstep = 100 ;

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

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

va_start() ;
nargin = nargin - 3 ;
while nargin>0 ,
  tmp = va_arg() ; nargin-- ;
  if ! isstr(tmp) , break end	# Last option has been read. 'tmp' is
				# first arg.

  if index(opt2,[" ",tmp," "]) ,
    
    tmp2 = va_arg() ; nargin-- ;
    nargin-- ;
    eval([tmp,"=tmp2;"]) ;
    if verbose ,
      printf("conjgrad_min : Read option : %s. Value is :\n",tmp) ;
      tmp2 
    end

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


N = prod(size(x)) ;
x = reshape(x,N,1) ;
x0 = x ;
if isnan(vinit) , v = feval( func , x ) ; 
else              v = vinit ; 
end

xi = h = df = -feval( dfunc, x ) ;

while niter++ <= maxiter ,
  xprev = x ;
  [w,vnew,ok] = brent0( xi', func, x, 1 ) ;
  ## [w,vnew,ok] = brent( "auto", xi', func, x ) ;
  x = x + w*xi' ; 
  if 2*abs(vnew-v) <= gtol*( abs(vnew)+abs(v)+100*eps ),
    return
  end
  if vnew>v ,
    printf("conjgrad_min: step increased cost function\n");
    keyboard
  end
  
  ## if abs(1-(x-xprev)'*xi'/norm(xi)/norm(x-xprev))>1000*eps,
  ##  printf("conjgrad_min: step is not in the right direction\n");
  ##  keyboard
  ## end
  v = feval(func,x) ;
  if verbose, printf("conjgrad_min : niter=%4i, v=%8.3g\n",niter,v) ; end
  ## if v2>vnew ,			# Should hold, imho
  ##  printf("conjgrad_min: Whoa!! \n");
  ##  keyboard
  ## end
  xi = feval(dfunc,x) ;
  gg = df*df' ;
  if gg == 0,
    return
  end
  dgg = (df+xi)*xi' ;
  gam = dgg/gg ;
  df = -xi ;
  xi = h = df+gam*h ;
end
printf("conjgrad_min: Iterates forever!\n");