linesearch

(tmp_for_tar/deriv_min.doc/linesearch.m)

Function Synopsis

[x,v,ok] = linesearch(df,dx,func,xinit,...)

Help text

       [x,v] = linesearch(df,dx,func,xinit,options)

 Last modified: October 2000

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

##       [x,v] = linesearch(df,dx,func,xinit,options)
##


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

function [x,v,ok] = linesearch(df,dx,func,xinit,...)

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

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

debug.dummy = 0 ;		# Dummy structure element
niter = 0 ;

tol = 1e-7 ;			# too small a step
alf = 1e-4 ;

maxstep = nan ;
ok = 1 ;

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

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(opt2,[" ",tmp," "]) ,
    
    tmp2 = va_arg() ; nargin-- ;
    nargin-- ;
    eval([tmp,"=tmp2;"]) ;
    if verbose ,
      printf("linesearch : Read option : %s. Value is :\n",tmp) ;
      ## tmp2 
    end

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

if isnan(vinit) , vinit = feval( func , xinit ) ; end

N = prod(size(xinit)) ;

sm = sqrt(sum(dx.^2)) ;
if !isnan(maxstep) && sm>maxstep, 
  dx = dx*(maxstep/sm) ;
end
slope = sum(df(:).*dx(:)) ;	# Avoid matrix multiplication
if slope > 0 ,
  printf("linesearch : slope is > 0\n") ; 
  ##   dx = -dx ;
  ## slope = -slope ;
  keyboard
end

tmp = abs(dx)./max(abs(xinit),1) ;
test = max(tmp(:)) ;

lamin = tol/test ;
lam = 1;
lam2 = lam ;

while niter++>=0,		# Never stop
  ## dx2 = (dx*lam +(lam2-lam)*df)/lam2 ;
  ## df'*dx2/sqrt(dx2'*dx2+df'*df)
  ## sm = sqrt(sum(dx2.^2)) ;
  ## if !isnan(maxstep) && sm>maxstep, 
  ##   dx2 = dx2*(maxstep/sm) ;
  ## end
  ## slope = sum(df(:).*dx2(:)) ;
  x = xinit+lam*dx ;
  v = feval(func,x) ;
  if report && rem(niter-1,report) == 0 ,
    printf("linesearch : niter=%d, v=%8.3g,vinit=%8.3g, lam=%8.3g\n",...
	   niter,v,vinit,lam); 
  end
  if lam < lamin ,		# Failure
    x = xinit ;
    ok = 0 ;
    return 
				# Success
  elseif v < vinit+alf*lam*slope ,
    return 
				# Continue
  elseif lam==1			# 1st time around
    tmp = -slope/(2*(v-vinit-slope)) ;
  else				# Not 1st time around
    rh1 = v-vinit-lam*slope ;
    rh2 = v2-vinit2-lam2*slope ;
    a = ( rh1/lam^2 - rh2/lam2^2 )/(lam-lam2) ;
    b = (-lam2*rh1/lam^2+lam*rh2/lam2^2)/(lam-lam2) ;
    if a==0, 
      tmp = -slope/(2*b) ;
    else   
      d = b^2-3*a*slope ;
      if d < 0,
	printf("linesearch : Roundoff problem\n");
      else
	tmp = (-b+sqrt(d))/(3*a) ;
      end
    end
  end
  if tmp>0.5*lam, tmp = 0.5*lam ; end
  lam2 = lam ;
  v2 = v ;
  vinit2 = vinit ;
  lam = max(tmp,0.1*lam);
end