1 ## Copyright (C) 2000 Ben Sapp <bsapp@lanl.gov>
2 ## Copyright (C) 2002 Etienne Grossmann <etienne@egdn.net>
3 ## Copyright (C) 2011 Nir Krakauer nkrakauer@ccny.cuny.edu
5 ## This program is free software; you can redistribute it and/or modify it under
6 ## the terms of the GNU General Public License as published by the Free Software
7 ## Foundation; either version 3 of the License, or (at your option) any later
10 ## This program is distributed in the hope that it will be useful, but WITHOUT
11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
15 ## You should have received a copy of the GNU General Public License along with
16 ## this program; if not, see <http://www.gnu.org/licenses/>.
18 ## [a,fx,nev] = line_min (f, dx, args, narg, h, nev_max) - Minimize f() along dx
21 ## f : string : Name of minimized function
22 ## dx : matrix : Direction along which f() is minimized
23 ## args : cell : Arguments of f
24 ## narg : integer : Position of minimized variable in args. Default=1
25 ## h : scalar : Step size to use for centered finite difference
26 ## approximation of first and second derivatives. Default=1E-3.
27 ## nev_max : integer : Maximum number of function evaluations. Default=30
30 ## a : scalar : Value for which f(x+a*dx) is a minimum (*)
31 ## fx : scalar : Value of f(x+a*dx) at minimum (*)
32 ## nev : integer : Number of function evaluations
34 ## (*) The notation f(x+a*dx) assumes that args == {x}.
36 ## Reference: David G Luenberger's Linear and Nonlinear Programming
38 function [a,fx,nev] = line_min (f, dx, args, narg, h, nev_max)
42 if (nargin < 4) narg = 1; endif
43 if (nargin < 5) h = 0.001; endif
44 if (nargin < 6) nev_max = 30; endif
50 min_velocity_change = 0.000001;
52 while (abs (velocity) > min_velocity_change && nev < nev_max)
53 fx = feval (f,args{1:narg-1}, x+a*dx, args{narg+1:end});
54 fxph = feval (f,args{1:narg-1}, x+(a+h)*dx, args{narg+1:end});
55 fxmh = feval (f,args{1:narg-1}, x+(a-h)*dx, args{narg+1:end});
60 velocity = (fxph - fxmh)/(2*h);
61 acceleration = (fxph - 2*fx + fxmh)/(h^2);
62 if abs(acceleration) <= eps, acceleration = 1; end # Don't do div by zero
63 # Use abs(accel) to avoid problems due to
65 a = a - velocity/abs(acceleration);
69 fx = feval (f, args{1:narg-1}, x+a*dx, args{narg+1:end});
71 if fx >= fx0 # if no improvement, return the starting value
77 disp ("line_min: maximum number of function evaluations reached")
82 ## Rem : Although not clear from the code, the returned a always seems to
83 ## correspond to (nearly) optimal fx.