From octave-sources-request at bevo dot che dot wisc dot edu Fri Oct 31 09:52:17 2003 Subject: Line minimization routines - proposed modification to octave forge line_min From: Michael Creel To: octave-sources at bevo dot che dot wisc dot edu Date: Fri, 31 Oct 2003 17:51:14 +0100 --Boundary_(ID_ZFiIUjqEC5tBth+oWVtnOQ) Content-type: text/plain; charset=us-ascii Content-disposition: inline Hello all, Getting ready to do some large scale minimizations, I've been looking into getting an efficient and robust line minimization routine working. Up until now I've been using my own routine, line_search.m: line_search.m: * This is a robust routine that checks to see that it really finds a stepsize that leads to an improvement. On the other hand, there is line_min in octave forge. * this is "efficient" in that it tries to find the best stepsize. Seems like golden step, if my neurons aren't failing A couple of problems with this code: * It does not limit the number of iterations used to find the best stepsize * it does not check that the final stepsize is in fact better than a stepsize of zero. This is a problem that occurs in practice with this code, not just a possibility. Now, there is no reason to try to get the exact best stepsize, if it takes a zillion iterations. We're trying to minimize the time needed to minimize the objective function. The many function evaluations might better be used to re-evaluate the direction, rather than find a really optimal stepsize. So, there is a tradeoff. The routine mc_line_min.m combines the two approaches. It is basically line_min, but with a limit on the maximum number of function evaluations used to get the best stepsize that is proportional to the number of parameters. Also, checking for actual improvement is enforced, and it falls back to the line_search routine if none is found. In some benchmarking, I have found that mc_line_min uses about 0.6 the time that line_search does. Furthermore, the line_min that comes with octave forge is prone to getting hung up doing a zillion iterations, and in my work, it is not usable. So, for your minimizing pleasure, I attach the following so you can try it out yourself. Regards, Michael --Boundary_(ID_ZFiIUjqEC5tBth+oWVtnOQ) Content-type: text/octave; name=line_search.m; charset=us-ascii Content-disposition: attachment; filename=line_search.m Content-transfer-encoding: 7BIT #============================ line_search ===================================== # own inspiration, no credits necessary # line_search : Subject to using few iterations, # find value of a that leads to greatest reduction # in f, moving in direction dx function [a, fmin] = line_search (f, dx, args, narg) x = nth (args, narg); fin = leval(f, splice(args, narg, 1, list(x))); a = 1; f0 = fin; # this first loop goes until an improvement is found while a > 2*eps # limit iterations fx = leval (f,splice (args, narg, 1, list (x+a*dx))); if (fx > f0) || isnan(fx) # reduce stepsize if worse, or if function can't be evaluated a = 0.5 * a; else f0 = fx; break; endif endwhile # now keep going until we no longer improve, or reach max trials while a > 2*eps a = 0.5*a; fx = leval (f,splice (args, narg, 1, list (x+a*dx))); # if improved, record new best and try another step if (fx < f0) & !isnan(fx) f0 = fx; else a = a / 0.5; # put it back to best found break; endif; endwhile fmin = f0; endfunction #============================ end line_search ================================== === --Boundary_(ID_ZFiIUjqEC5tBth+oWVtnOQ) Content-type: text/octave; name=mc_line_min.m; charset=us-ascii Content-disposition: attachment; filename=mc_line_min.m Content-transfer-encoding: 7BIT ## Copyright (C) 2000 Ben Sapp. All rights reserved. ## ## This program is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by the ## Free Software Foundation; either version 2, or (at your option) any ## later version. ## ## This is distributed in the hope that it will be useful, but WITHOUT ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License ## for more details. ## [a,fx,nev] = line_min (f, dx, args, narg) - Minimize f() along dx ## ## INPUT ---------- ## f : string : Name of minimized function ## dx : matrix : Direction along which f() is minimized ## args : list : List of argument of f ## narg : integer : Position of minimized variable in args. Default=1 ## ## OUTPUT --------- ## a : scalar : Value for which f(x+a*dx) is a minimum (*) ## fx : scalar : Value of f(x+a*dx) at minimum (*) ## nev : integer : Number of function evaluations ## ## (*) The notation f(x+a*dx) assumes that args == list (x). ## Author: Ben Sapp ## Reference: David G Luenberger's Linear and Nonlinear Programming ## ## Changelog : ----------- ## 2002-01-28 Paul Kienzle ## * save two function evaluations by inlining the derivatives ## * pass through varargin{:} to the function ## 2002-03-13 Paul Kienzle ## * simplify update expression ## 2002-04-17 Etienne Grossmann ## * Rename nrm.m to line_min.m (in order not to break dfp, which uses nrm) ## * Use list of args, suppress call to __pseudo_func__ ## * Add nargs argument, assume args is a list ## * Change help text ## 2003-10-31 -- boo! Michael Creel ## * limit number of function evaluations used, and fall back to ## * a robust brute force method when no improvement found. function [a,fx,nev] = mc_line_min(f, dx, args, narg) velocity = 1; acceleration = 1; if nargin < 4, narg = 1; end nev = 0; h = 0.001; # Was 0.01 here x = nth (args,narg); a = 0; f0 = leval (f,splice (args, narg, 1, list (x+a*dx))); while ((abs (velocity) > 0.000001) & nev < 3*rows(dx)) fx = leval (f,splice (args, narg, 1, list (x+a*dx))); if (fx > f0) break; endif fxph = leval (f,splice (args, narg,1,list (x+(a+h)*dx))); fxmh = leval (f,splice (args, narg,1,list (x+(a-h)*dx))); velocity = (fxph - fxmh)/(2*h); acceleration = (fxph - 2*fx + fxmh)/(h^2); if abs(acceleration) <= eps, acceleration = 1; end # Don't do div by zero # Use abs(accel) to avoid problems due to # concave function a = a - velocity/abs(acceleration); nev += 3; endwhile if (fx > f0) printf("caught you!\n"); a = line_search(f, dx, args,narg); endif endfunction ## Rem : Although not clear from the code, the returned a always seems to ## correspond to (nearly) optimal fx. --Boundary_(ID_ZFiIUjqEC5tBth+oWVtnOQ)--