```Update of /cvsroot/octave/octave-forge/main/optim/Optimize In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv21597/main/optim/Optimize Modified Files: BFGSMin.m Battery.m BisectionStep.m ExampleBFGSMin.m ExampleNumGradient.m NewtonStep.m NumGradient.m NumHessian.m QuadInterpStep.m Log Message: new (stable as of now) interface for BFGSMin Index: BFGSMin.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/BFGSMin.m,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- BFGSMin.m 17 Feb 2004 09:42:35 -0000 1.4 +++ BFGSMin.m 26 Feb 2004 10:18:42 -0000 1.5 @@ -19,15 +19,22 @@ #============================ BFGSMin.m ===================================== # Minimize a function of the form # -# value = f(args) +# Value = f(args) # # or # -# [value, gradient] = f(args) +# ValueGradient = f(args) # -# where args is a cell array. +# where # -# Minimization is with respect to the FIRST element of args. +# "args" is a cell array, minimization is w.r.t. FIRST element +# "Value" is the scalar value of f +# "ValueGradient" = {value, gradient} is a cell array +# containing the value and the gradient of f. +# +# That is, your function can either return a scalar value, +# or it can return the value and the gradient vector together +# in a cell array. # # Brief summary of syntax follows, for a more extensive example, # that shows how to write objective functions, see the file ExampleBFGSMin.m @@ -35,15 +42,15 @@ # # Calling syntax: # -# [theta, obj_value, iters, convergence] = BFGSMin(func, args, control, gradient) +# [theta, obj_value, iters, convergence] = BFGSMin(func, args, control) # -# NOTE: the last two arguments are optional +# NOTE: the last argument is optional # # # Input arguments: # func - (string) the function to minimize # args - arguments, in a cell array. -# control - OPTIONAL 3x1 vector OR a string, e.g., "default" +# control - OPTIONAL 3x1 vector # * 1st elem. controls maximum iterations # scalar > 0 - max. iters # or -1 for infinity (default) @@ -56,23 +63,6 @@ # are all tested (default) # 0 = only function conv tested # -# OR -# -# If you want defaults for control but want to set gradient -# (the following argument), set control to a string, -# e.g., "defaults". This is just for convenience, so that -# you don't have to remember the defaults. -# -# gradient - OPTIONAL: controls for analytic or numeric gradient -# * not provided, or not a string: use numeric gradient (default), -# function can be of form -# value = f(args) -# OR -# [value, gradient] = f(args) -# but numeric gradient will be used -# * any string - e.g., "yes" - function provides gradient, e.g., -# [value, gradient] = f(args) -# # Outputs: # theta - the minimizer # obj_value - the minimized funtion value @@ -82,7 +72,7 @@ function [theta, obj_value, iters, convergence] = BFGSMin(func, args, control, gradient_type) # Check number of args - if (nargin > 4) error("\nBFGSMin: you have provided more than 4 arguments\n"); endif + if (nargin > 3) error("\nBFGSMin: you have provided more than 3 arguments\n"); endif # Check types of required arguments if !ischar(func) error("\nBFGSMin: first argument must be a string that gives name of objective function\n"); endif @@ -99,40 +89,37 @@ max_iters = inf; convergence = -1; # if this doesn't change, it means that maxiters were exceeded verbosity = 0; - numeric_gradient = 1; # Now use options that are user-provided # Control verbosity, max iterations, and arg wrt which we minimize (optional) if nargin > 2 - if !ischar(control); - if (rows(control) != 3) - error("\nBFGSMin: 3rd argument (control) must be a 3x1 vector\n"); - else; - max_iters = control(1,:); - if max_iters == -1, max_iters = inf; endif - verbosity = control(2,:); - criterion = control(3,:); - endif - endif + if !isvector(control) error("\nBFGSMin: 3rd argument must be a vector\n"); endif + control = control(:); + if (rows(control) != 3) + error("\nBFGSMin: 3rd argument must be a 3x1 vector\n"); + else + max_iters = control(1,:); + if max_iters == -1, max_iters = inf; endif + verbosity = control(2,:); + criterion = control(3,:); + endif endif - # Check if analytic gradient is provided (optional - default is to use numeric) - if nargin > 3 - numeric_gradient = !ischar(gradient_type); - numeric_gradient - endif - # initialize things theta = args{1}; thetain = theta; H = eye (prod (sz = size (theta))); # Initial gradient and obj_value - if numeric_gradient - last_obj_value = feval(func, args); + temp = feval(func, args); + if isscalar(temp) + numeric_gradient = 1; # remember for future iterations + last_obj_value = temp; g = NumGradient(func, args); else - [last_obj_value, g] = feval(func, args); + numeric_gradient = 0; + last_obj_value = temp{1}; + g = temp{2}; endif g = g(:); # make sure it's a column vector @@ -208,7 +195,8 @@ if numeric_gradient g_new = NumGradient(func, args); else - [dummy, g_new] = feval(func, args); + temp = feval(func, args); + g_new = temp{2}; endif g_new = g_new(:); Index: Battery.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/Battery.m,v retrieving revision 1.2 retrieving revision 1.3 diff -u -d -r1.2 -r1.3 --- Battery.m 17 Feb 2004 09:43:00 -0000 1.2 +++ Battery.m 26 Feb 2004 10:18:42 -0000 1.3 @@ -28,7 +28,7 @@ [k,trials] = size(startvals); bestobj = inf; besttheta = zeros(k,1); -bfgscontrol = [maxiters;0;1]; +bfgscontrol = [maxiters;0;0]; # now try the supplied start values, and optionally the random start values for i = 1:trials args{1} = startvals(:,i); Index: BisectionStep.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/BisectionStep.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- BisectionStep.m 26 Jan 2004 10:29:04 -0000 1.1 +++ BisectionStep.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -18,7 +18,7 @@ #============================ BisectionStep ===================================== # -# this is for use by BFGSMin +# this is used by BFGSMin # # Uses bisection to find a stepsize that leads to a decrease, then # continues until no further improvement @@ -26,15 +26,16 @@ # 13/01/2004 # # usage: -# [a, obj_value] = BisectionStep(f, dx, args) +# +# [a, obj_value] = BisectionStep(f, dx, args) +# # inputs: -# f: the objective function -# dx: the direction -# args: the arguments of the function, in a cell array. -# The first argument is the one w.r.t. -# which we are minimizing. - - +# +# f: the objective function +# dx: the direction +# args: the arguments of the function, in a cell array. +# The first argument is the one w.r.t. +# which we are minimizing. function [a, obj] = BisectionStep(f, dx, args) @@ -43,12 +44,16 @@ args_in = args; obj_0 = feval(f, args); + if iscell(obj_0) obj_0 = obj_0{1}; endif + a = 1; # this first loop goes until an improvement is found while a > 2*eps # limit iterations args{1} = x + a*dx; obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif + if (obj > obj_0) || isnan(obj) # reduce stepsize if worse, or if function can't be evaluated a = 0.5 * a; else @@ -63,6 +68,7 @@ a = 0.5*a; args{1} = x + a*dx; obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif # if improved, record new best and try another step if ((obj < obj_0) & !isnan(obj)) Index: ExampleBFGSMin.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/ExampleBFGSMin.m,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- ExampleBFGSMin.m 16 Feb 2004 15:17:01 -0000 1.4 +++ ExampleBFGSMin.m 26 Feb 2004 10:18:42 -0000 1.5 @@ -25,11 +25,11 @@ # The argument w.r.t. which we minimize must be the FIRST argument # in the cell array that is the argument to BFGSMin. -1; - # Function value and gradient vector of the rosenbrock function # The minimizer is at the vector (1,1,..,1), # and the minimized value is 0. +1; + function [obj_value, gradient] = rosenbrock(x); dimension = length(x); obj_value = sum(100*(x(2:dimension)-x(1:dimension-1).^2).^2 + (1-x(1:dimension-1)).^2); @@ -43,11 +43,12 @@ # The example objective functions # example obj. fn. - this shows how to use analytic gradient -function [obj_value, grad] = objective(args) +function output = objective(args) theta = args{1}; location = args{2}; x = theta - location + ones(rows(theta),1); # move minimizer to "location" [obj_value, grad] = rosenbrock(x); + output = {obj_value, grad}; endfunction # example obj. fn. - this shows how to use numerical grafient @@ -76,7 +77,7 @@ printf("a vector of even steps from 0 to 5\n\n"); printf("ANALYTIC GRADIENT\n"); -[theta, obj_value, iterations, convergence] = BFGSMin("objective", args, control, "analytic"); +[theta, obj_value, iterations, convergence] = BFGSMin("objective", args, control); t = cputime() - t; printf("Elapsed time = %f\n",t); Index: ExampleNumGradient.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/ExampleNumGradient.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- ExampleNumGradient.m 26 Jan 2004 10:29:04 -0000 1.1 +++ ExampleNumGradient.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -21,6 +21,7 @@ # The Rosenbrock function - used to define objective functions # Function value and gradient vector of the rosenbrock function # The minimizer is x = [1; 1;...;1], and the minimized value is 0. +1; function [obj_value, gradient] = rosenbrock(x); dimension = length(x); obj_value = sum(100*(x(2:dimension)-x(1:dimension-1).^2).^2 + (1-x(1:dimension-1)).^2); @@ -34,16 +35,16 @@ # The example objective functions # example obj. fn. - this shows how to use analytic gradient -function [obj_value, grad] = objective(args) +function output = objective(args) theta = args{1}; [obj_value, grad] = rosenbrock(theta); + output = {obj_value, grad}; endfunction # example obj. fn. - this shows how to use numerical grafient -function [obj_value, grad] = objective2(args) +function obj_value = objective2(args) theta = args{1}; obj_value = rosenbrock(theta); - grad = "na"; endfunction # Here's how to get the numeric gradient @@ -51,7 +52,8 @@ ngradient = NumGradient("objective", args); # Analytic gradient, to verify -[o, agradient] = objective(args); +output = objective(args); +agradient = output{2}; ngradient = ngradient'; printf("\nCheck NumGradient\n"); printf("\nDerivatives of the Rosenbrock function, random start values\n\n"); Index: NewtonStep.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/NewtonStep.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- NewtonStep.m 26 Jan 2004 10:29:04 -0000 1.1 +++ NewtonStep.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -42,8 +42,10 @@ gradient = 1; evaluations = 0; - obj_0 = obj = feval(f, args); - + obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif + obj_0 = obj; + delta = 0.001; # experimentation show that this is a good choice x_right = x + delta*dx; @@ -51,8 +53,10 @@ args{1} = x_right; obj_right = feval(f, args); + if iscell(obj_right) obj_right = obj_right{1}; endif args{1} = x_left; obj_left = feval(f, args); + if iscell(obj_left) obj_left = obj_left{1}; endif gradient = (obj_right - obj_left) / (2*delta); # take central difference hessian = (obj_right - 2*obj + obj_left) / (delta^2); @@ -60,10 +64,13 @@ if hessian <= eps, hessian = 1; endif # avoid div by zero a = - gradient/hessian; # hessian inverse gradient: the Newton step - + test = a < 5; + a = a*test + 5*(1-test); # Let's avoid extreme steps that might cause crashes # check that this is improvement args{1} = x+a*dx; obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif + # if not, fall back to bisection if ((obj > obj_0) || isnan(obj)) [a, obj] = BisectionStep(f, dx, args_in); Index: NumGradient.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/NumGradient.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- NumGradient.m 26 Jan 2004 10:29:04 -0000 1.1 +++ NumGradient.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -24,6 +24,8 @@ # # * Allows diff. of vector-valued function # * Uses systematic finite difference +# +# See ExampleNumGradient.m function derivative = NumGradient(f, args) @@ -33,7 +35,8 @@ k = rows(parameter); obj_value = feval(f, args); - + if iscell(obj_value) obj_value = obj_value{1}; endif + n = rows(obj_value); derivative = zeros(n, k); @@ -46,12 +49,14 @@ delta_right = d - p; args{1} = parameter; obj_right = feval(f, args); - - # left size + if iscell(obj_right) obj_right = obj_right{1}; endif + + # left size parameter(i) = d = p - delta; delta_left = p - d; args{1} = parameter; obj_left = feval(f, args); + if iscell(obj_left) obj_left = obj_left{1}; endif parameter(i) = p; # restore original parameter Index: NumHessian.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/NumHessian.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- NumHessian.m 26 Jan 2004 10:29:04 -0000 1.1 +++ NumHessian.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -25,7 +25,7 @@ k = rows(parameter); obj_value = feval(f, args); - + if iscell(obj_value) obj_value = obj_value{1}; endif derivative = zeros(k, k); for i = 1:k # approximate 2nd deriv. by central difference @@ -39,37 +39,40 @@ hia = di - pi; hja = dj - pj; args{1} = parameter; fpp = feval(f, args); - + if iscell(fpp) fpp = fpp{1}; endif parameter(i) = di = pi - hi; parameter(j) = dj = pj - hj; # -1 -1 hia = hia + pi - di; hja = hja + pj - dj; args{1} = parameter; fmm = feval(f, args); - + if iscell(fmm) fmm = fmm{1}; endif + parameter(i) = pi + hi; parameter(j) = pj - hj; # +1 -1 args{1} = parameter; fpm = feval(f, args); - + if iscell(fpm) fpm = fpm{1}; endif parameter(i) = pi - hi; parameter(j) = pj + hj; # -1 +1 args{1} = parameter; fmp = feval(f, args); - + if iscell(fmp) fmp = fmp{1}; endif derivative(j,i) = ((fpp - fpm) + (fmm - fmp)) / (hia * hja); derivative(i,j) = derivative(j,i); parameter(j) = pj; endfor - # diagonal elements + + # diagonal elements parameter(i) = di = pi + 2 * hi; # +1 +1 args{1} = parameter; fpp = feval(f, args); + if iscell(fpp) fpp = fpp{1}; endif hia = (di - pi) / 2; parameter(i) = di = pi - 2 * hi; # -1 -1 args{1} = parameter; fmm = feval(f, args); - + if iscell(fmm) fmm = fmm{1}; endif hia = hia + (pi - di) / 2; derivative(i,i) = ((fpp - obj_value) + (fmm - obj_value)) / (hia * hia); Index: QuadInterpStep.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/QuadInterpStep.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- QuadInterpStep.m 26 Jan 2004 10:29:04 -0000 1.1 +++ QuadInterpStep.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -42,17 +42,20 @@ evaluations = 0; obj_0 = feval(f, args); - + if iscell(obj_0) obj_0 = obj_0{1}; endif left = 0.1; right = 3; center = 1; args{1} = x + left*dx; obj_left = feval(f, args); + if iscell(obj_left) obj_left = obj_left{1}; endif args{1} = x + center*dx; obj_center = feval(f, args); + if iscell(obj_center) obj_center = obj_center{1}; endif args{1} = x + right*dx; obj_right = feval(f, args); + if iscell(obj_right) obj_right = obj_right{1}; endif # Best is on L extreme if (obj_left < obj_center) & (obj_left < obj_right) @@ -69,6 +72,8 @@ a = 0.5 * a / (obj_left*(center - right) + obj_center*(right - left) + obj_right*(left - center)); args{1} = x + a*dx; obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif + endif # fall back to bisection if this is not improvement or there is some sort of crash ```
 ```Update of /cvsroot/octave/octave-forge/main/optim/Optimize In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv21597/main/optim/Optimize Modified Files: BFGSMin.m Battery.m BisectionStep.m ExampleBFGSMin.m ExampleNumGradient.m NewtonStep.m NumGradient.m NumHessian.m QuadInterpStep.m Log Message: new (stable as of now) interface for BFGSMin Index: BFGSMin.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/BFGSMin.m,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- BFGSMin.m 17 Feb 2004 09:42:35 -0000 1.4 +++ BFGSMin.m 26 Feb 2004 10:18:42 -0000 1.5 @@ -19,15 +19,22 @@ #============================ BFGSMin.m ===================================== # Minimize a function of the form # -# value = f(args) +# Value = f(args) # # or # -# [value, gradient] = f(args) +# ValueGradient = f(args) # -# where args is a cell array. +# where # -# Minimization is with respect to the FIRST element of args. +# "args" is a cell array, minimization is w.r.t. FIRST element +# "Value" is the scalar value of f +# "ValueGradient" = {value, gradient} is a cell array +# containing the value and the gradient of f. +# +# That is, your function can either return a scalar value, +# or it can return the value and the gradient vector together +# in a cell array. # # Brief summary of syntax follows, for a more extensive example, # that shows how to write objective functions, see the file ExampleBFGSMin.m @@ -35,15 +42,15 @@ # # Calling syntax: # -# [theta, obj_value, iters, convergence] = BFGSMin(func, args, control, gradient) +# [theta, obj_value, iters, convergence] = BFGSMin(func, args, control) # -# NOTE: the last two arguments are optional +# NOTE: the last argument is optional # # # Input arguments: # func - (string) the function to minimize # args - arguments, in a cell array. -# control - OPTIONAL 3x1 vector OR a string, e.g., "default" +# control - OPTIONAL 3x1 vector # * 1st elem. controls maximum iterations # scalar > 0 - max. iters # or -1 for infinity (default) @@ -56,23 +63,6 @@ # are all tested (default) # 0 = only function conv tested # -# OR -# -# If you want defaults for control but want to set gradient -# (the following argument), set control to a string, -# e.g., "defaults". This is just for convenience, so that -# you don't have to remember the defaults. -# -# gradient - OPTIONAL: controls for analytic or numeric gradient -# * not provided, or not a string: use numeric gradient (default), -# function can be of form -# value = f(args) -# OR -# [value, gradient] = f(args) -# but numeric gradient will be used -# * any string - e.g., "yes" - function provides gradient, e.g., -# [value, gradient] = f(args) -# # Outputs: # theta - the minimizer # obj_value - the minimized funtion value @@ -82,7 +72,7 @@ function [theta, obj_value, iters, convergence] = BFGSMin(func, args, control, gradient_type) # Check number of args - if (nargin > 4) error("\nBFGSMin: you have provided more than 4 arguments\n"); endif + if (nargin > 3) error("\nBFGSMin: you have provided more than 3 arguments\n"); endif # Check types of required arguments if !ischar(func) error("\nBFGSMin: first argument must be a string that gives name of objective function\n"); endif @@ -99,40 +89,37 @@ max_iters = inf; convergence = -1; # if this doesn't change, it means that maxiters were exceeded verbosity = 0; - numeric_gradient = 1; # Now use options that are user-provided # Control verbosity, max iterations, and arg wrt which we minimize (optional) if nargin > 2 - if !ischar(control); - if (rows(control) != 3) - error("\nBFGSMin: 3rd argument (control) must be a 3x1 vector\n"); - else; - max_iters = control(1,:); - if max_iters == -1, max_iters = inf; endif - verbosity = control(2,:); - criterion = control(3,:); - endif - endif + if !isvector(control) error("\nBFGSMin: 3rd argument must be a vector\n"); endif + control = control(:); + if (rows(control) != 3) + error("\nBFGSMin: 3rd argument must be a 3x1 vector\n"); + else + max_iters = control(1,:); + if max_iters == -1, max_iters = inf; endif + verbosity = control(2,:); + criterion = control(3,:); + endif endif - # Check if analytic gradient is provided (optional - default is to use numeric) - if nargin > 3 - numeric_gradient = !ischar(gradient_type); - numeric_gradient - endif - # initialize things theta = args{1}; thetain = theta; H = eye (prod (sz = size (theta))); # Initial gradient and obj_value - if numeric_gradient - last_obj_value = feval(func, args); + temp = feval(func, args); + if isscalar(temp) + numeric_gradient = 1; # remember for future iterations + last_obj_value = temp; g = NumGradient(func, args); else - [last_obj_value, g] = feval(func, args); + numeric_gradient = 0; + last_obj_value = temp{1}; + g = temp{2}; endif g = g(:); # make sure it's a column vector @@ -208,7 +195,8 @@ if numeric_gradient g_new = NumGradient(func, args); else - [dummy, g_new] = feval(func, args); + temp = feval(func, args); + g_new = temp{2}; endif g_new = g_new(:); Index: Battery.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/Battery.m,v retrieving revision 1.2 retrieving revision 1.3 diff -u -d -r1.2 -r1.3 --- Battery.m 17 Feb 2004 09:43:00 -0000 1.2 +++ Battery.m 26 Feb 2004 10:18:42 -0000 1.3 @@ -28,7 +28,7 @@ [k,trials] = size(startvals); bestobj = inf; besttheta = zeros(k,1); -bfgscontrol = [maxiters;0;1]; +bfgscontrol = [maxiters;0;0]; # now try the supplied start values, and optionally the random start values for i = 1:trials args{1} = startvals(:,i); Index: BisectionStep.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/BisectionStep.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- BisectionStep.m 26 Jan 2004 10:29:04 -0000 1.1 +++ BisectionStep.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -18,7 +18,7 @@ #============================ BisectionStep ===================================== # -# this is for use by BFGSMin +# this is used by BFGSMin # # Uses bisection to find a stepsize that leads to a decrease, then # continues until no further improvement @@ -26,15 +26,16 @@ # 13/01/2004 # # usage: -# [a, obj_value] = BisectionStep(f, dx, args) +# +# [a, obj_value] = BisectionStep(f, dx, args) +# # inputs: -# f: the objective function -# dx: the direction -# args: the arguments of the function, in a cell array. -# The first argument is the one w.r.t. -# which we are minimizing. - - +# +# f: the objective function +# dx: the direction +# args: the arguments of the function, in a cell array. +# The first argument is the one w.r.t. +# which we are minimizing. function [a, obj] = BisectionStep(f, dx, args) @@ -43,12 +44,16 @@ args_in = args; obj_0 = feval(f, args); + if iscell(obj_0) obj_0 = obj_0{1}; endif + a = 1; # this first loop goes until an improvement is found while a > 2*eps # limit iterations args{1} = x + a*dx; obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif + if (obj > obj_0) || isnan(obj) # reduce stepsize if worse, or if function can't be evaluated a = 0.5 * a; else @@ -63,6 +68,7 @@ a = 0.5*a; args{1} = x + a*dx; obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif # if improved, record new best and try another step if ((obj < obj_0) & !isnan(obj)) Index: ExampleBFGSMin.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/ExampleBFGSMin.m,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- ExampleBFGSMin.m 16 Feb 2004 15:17:01 -0000 1.4 +++ ExampleBFGSMin.m 26 Feb 2004 10:18:42 -0000 1.5 @@ -25,11 +25,11 @@ # The argument w.r.t. which we minimize must be the FIRST argument # in the cell array that is the argument to BFGSMin. -1; - # Function value and gradient vector of the rosenbrock function # The minimizer is at the vector (1,1,..,1), # and the minimized value is 0. +1; + function [obj_value, gradient] = rosenbrock(x); dimension = length(x); obj_value = sum(100*(x(2:dimension)-x(1:dimension-1).^2).^2 + (1-x(1:dimension-1)).^2); @@ -43,11 +43,12 @@ # The example objective functions # example obj. fn. - this shows how to use analytic gradient -function [obj_value, grad] = objective(args) +function output = objective(args) theta = args{1}; location = args{2}; x = theta - location + ones(rows(theta),1); # move minimizer to "location" [obj_value, grad] = rosenbrock(x); + output = {obj_value, grad}; endfunction # example obj. fn. - this shows how to use numerical grafient @@ -76,7 +77,7 @@ printf("a vector of even steps from 0 to 5\n\n"); printf("ANALYTIC GRADIENT\n"); -[theta, obj_value, iterations, convergence] = BFGSMin("objective", args, control, "analytic"); +[theta, obj_value, iterations, convergence] = BFGSMin("objective", args, control); t = cputime() - t; printf("Elapsed time = %f\n",t); Index: ExampleNumGradient.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/ExampleNumGradient.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- ExampleNumGradient.m 26 Jan 2004 10:29:04 -0000 1.1 +++ ExampleNumGradient.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -21,6 +21,7 @@ # The Rosenbrock function - used to define objective functions # Function value and gradient vector of the rosenbrock function # The minimizer is x = [1; 1;...;1], and the minimized value is 0. +1; function [obj_value, gradient] = rosenbrock(x); dimension = length(x); obj_value = sum(100*(x(2:dimension)-x(1:dimension-1).^2).^2 + (1-x(1:dimension-1)).^2); @@ -34,16 +35,16 @@ # The example objective functions # example obj. fn. - this shows how to use analytic gradient -function [obj_value, grad] = objective(args) +function output = objective(args) theta = args{1}; [obj_value, grad] = rosenbrock(theta); + output = {obj_value, grad}; endfunction # example obj. fn. - this shows how to use numerical grafient -function [obj_value, grad] = objective2(args) +function obj_value = objective2(args) theta = args{1}; obj_value = rosenbrock(theta); - grad = "na"; endfunction # Here's how to get the numeric gradient @@ -51,7 +52,8 @@ ngradient = NumGradient("objective", args); # Analytic gradient, to verify -[o, agradient] = objective(args); +output = objective(args); +agradient = output{2}; ngradient = ngradient'; printf("\nCheck NumGradient\n"); printf("\nDerivatives of the Rosenbrock function, random start values\n\n"); Index: NewtonStep.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/NewtonStep.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- NewtonStep.m 26 Jan 2004 10:29:04 -0000 1.1 +++ NewtonStep.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -42,8 +42,10 @@ gradient = 1; evaluations = 0; - obj_0 = obj = feval(f, args); - + obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif + obj_0 = obj; + delta = 0.001; # experimentation show that this is a good choice x_right = x + delta*dx; @@ -51,8 +53,10 @@ args{1} = x_right; obj_right = feval(f, args); + if iscell(obj_right) obj_right = obj_right{1}; endif args{1} = x_left; obj_left = feval(f, args); + if iscell(obj_left) obj_left = obj_left{1}; endif gradient = (obj_right - obj_left) / (2*delta); # take central difference hessian = (obj_right - 2*obj + obj_left) / (delta^2); @@ -60,10 +64,13 @@ if hessian <= eps, hessian = 1; endif # avoid div by zero a = - gradient/hessian; # hessian inverse gradient: the Newton step - + test = a < 5; + a = a*test + 5*(1-test); # Let's avoid extreme steps that might cause crashes # check that this is improvement args{1} = x+a*dx; obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif + # if not, fall back to bisection if ((obj > obj_0) || isnan(obj)) [a, obj] = BisectionStep(f, dx, args_in); Index: NumGradient.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/NumGradient.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- NumGradient.m 26 Jan 2004 10:29:04 -0000 1.1 +++ NumGradient.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -24,6 +24,8 @@ # # * Allows diff. of vector-valued function # * Uses systematic finite difference +# +# See ExampleNumGradient.m function derivative = NumGradient(f, args) @@ -33,7 +35,8 @@ k = rows(parameter); obj_value = feval(f, args); - + if iscell(obj_value) obj_value = obj_value{1}; endif + n = rows(obj_value); derivative = zeros(n, k); @@ -46,12 +49,14 @@ delta_right = d - p; args{1} = parameter; obj_right = feval(f, args); - - # left size + if iscell(obj_right) obj_right = obj_right{1}; endif + + # left size parameter(i) = d = p - delta; delta_left = p - d; args{1} = parameter; obj_left = feval(f, args); + if iscell(obj_left) obj_left = obj_left{1}; endif parameter(i) = p; # restore original parameter Index: NumHessian.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/NumHessian.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- NumHessian.m 26 Jan 2004 10:29:04 -0000 1.1 +++ NumHessian.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -25,7 +25,7 @@ k = rows(parameter); obj_value = feval(f, args); - + if iscell(obj_value) obj_value = obj_value{1}; endif derivative = zeros(k, k); for i = 1:k # approximate 2nd deriv. by central difference @@ -39,37 +39,40 @@ hia = di - pi; hja = dj - pj; args{1} = parameter; fpp = feval(f, args); - + if iscell(fpp) fpp = fpp{1}; endif parameter(i) = di = pi - hi; parameter(j) = dj = pj - hj; # -1 -1 hia = hia + pi - di; hja = hja + pj - dj; args{1} = parameter; fmm = feval(f, args); - + if iscell(fmm) fmm = fmm{1}; endif + parameter(i) = pi + hi; parameter(j) = pj - hj; # +1 -1 args{1} = parameter; fpm = feval(f, args); - + if iscell(fpm) fpm = fpm{1}; endif parameter(i) = pi - hi; parameter(j) = pj + hj; # -1 +1 args{1} = parameter; fmp = feval(f, args); - + if iscell(fmp) fmp = fmp{1}; endif derivative(j,i) = ((fpp - fpm) + (fmm - fmp)) / (hia * hja); derivative(i,j) = derivative(j,i); parameter(j) = pj; endfor - # diagonal elements + + # diagonal elements parameter(i) = di = pi + 2 * hi; # +1 +1 args{1} = parameter; fpp = feval(f, args); + if iscell(fpp) fpp = fpp{1}; endif hia = (di - pi) / 2; parameter(i) = di = pi - 2 * hi; # -1 -1 args{1} = parameter; fmm = feval(f, args); - + if iscell(fmm) fmm = fmm{1}; endif hia = hia + (pi - di) / 2; derivative(i,i) = ((fpp - obj_value) + (fmm - obj_value)) / (hia * hia); Index: QuadInterpStep.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/optim/Optimize/QuadInterpStep.m,v retrieving revision 1.1 retrieving revision 1.2 diff -u -d -r1.1 -r1.2 --- QuadInterpStep.m 26 Jan 2004 10:29:04 -0000 1.1 +++ QuadInterpStep.m 26 Feb 2004 10:18:42 -0000 1.2 @@ -42,17 +42,20 @@ evaluations = 0; obj_0 = feval(f, args); - + if iscell(obj_0) obj_0 = obj_0{1}; endif left = 0.1; right = 3; center = 1; args{1} = x + left*dx; obj_left = feval(f, args); + if iscell(obj_left) obj_left = obj_left{1}; endif args{1} = x + center*dx; obj_center = feval(f, args); + if iscell(obj_center) obj_center = obj_center{1}; endif args{1} = x + right*dx; obj_right = feval(f, args); + if iscell(obj_right) obj_right = obj_right{1}; endif # Best is on L extreme if (obj_left < obj_center) & (obj_left < obj_right) @@ -69,6 +72,8 @@ a = 0.5 * a / (obj_left*(center - right) + obj_center*(right - left) + obj_right*(left - center)); args{1} = x + a*dx; obj = feval(f, args); + if iscell(obj) obj = obj{1}; endif + endif # fall back to bisection if this is not improvement or there is some sort of crash ```