minimization-multi.lisp 15.1 KB
 liam committed Jan 21, 2008 1 2 ;; Multivariate minimization. ;; Liam Healy  Liam Healy committed Jan 23, 2009 3 ;; Time-stamp: <2009-01-22 21:35:26EST minimization-multi.lisp>  lhealy committed Jul 25, 2008 4 ;; $Id$  liam committed Jan 21, 2008 5 6 7  (in-package :gsl)  Liam Healy committed Dec 27, 2008 8 9 ;;; /usr/include/gsl/gsl_multimin.h  Liam Healy committed Jan 23, 2009 10 11 12 13 14 15 ;;; In the parabaloid example, I notice that the consruct ;;; (min-test-gradient (mfdfminimizer-gradient minimizer) 1.0d-3) ;;; is constructing a CL vector-double-float (in mfdfminimizer-gradient) and ;;; then immediately pulling out the pointer (in min-test-gradient). It ;;; is easy enough to eliminate this, but then mfdfminimizer-gradient ;;; would not be useful to a CL user.  liam committed Jan 21, 2008 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30  ;;;;**************************************************************************** ;;;; Function definition ;;;;**************************************************************************** ;;; The structures gsl-mfunction and gsl-mfunction-fdf are, from the ;;; CFFI point of view, equally valid for gsl_multimin_function and ;;; gsl_multimin_function_fdf defined in ;;; /usr/include/gsl/gsl_multimin.h as they are for ;;; gsl_multiroot_function and gsl_multiroot_function_fdf defined in ;;; /usr/include/gsl/gsl_multiroots.h. As far as CFFI is concerned, a ;;; pointer is a pointer, even though the C definition the functions ;;; they point to have different signatures. (export 'def-minimization-functions)  Liam Healy committed Jan 23, 2009 31 (defmacro def-minimization-functions (function dimensions &optional df fdf array)  liam committed Jan 21, 2008 32 33 34  "Setup functions for multivariate minimization. The CL functions name and derivative should be defined previously with defuns."  Liam Healy committed Jan 23, 2009 35 36 37 38 39 40 41 42 43 44 45 46 47 48  (let ((vdec (if array :pointer ((:double ,dimensions)))) (svdec (if array :pointer (:set :double ,dimensions)))) (progn (defmcallback ,function :double ,vdec nil t) ,@(when df ((defmcallback ,df :void ,vdec (,svdec) t) (defmcallback ,fdf :void ,vdec ((:set :double -1) ,svdec) t))) ,(if df (defcbstruct (,function function ,df df ,fdf fdf) gsl-mfunction-fdf ((dimensions ,dimensions))) (defcbstruct (,function function) gsl-mfunction ((dimensions ,dimensions)))))))  liam committed Jan 21, 2008 49 50 51 52 53  ;;;;**************************************************************************** ;;;; Initialization ;;;;****************************************************************************  Liam Healy committed Dec 25, 2008 54 55 56 57 58 59 60 61 62 63 (defmobject multi-dimensional-minimizer-f "gsl_multimin_fminimizer" ((type :pointer) (dimension sizet)) "multi-dimensional minimizer with function only" ; FDL "Make an instance of a minimizer of the given for an function of the given dimensions. Optionally initialize the minimizer to minimize the function starting from the initial point. The size of the initial trial steps is given in vector step-size. The precise meaning of this parameter depends on the method used." "set"  Liam Healy committed Dec 26, 2008 64  ;; Could have one fewer argument: dimension=(dim0 initial)  Liam Healy committed Dec 25, 2008 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82  ((function :pointer) ((mpointer initial) :pointer) ((mpointer step-size) :pointer))) (defmobject multi-dimensional-minimizer-fdf "gsl_multimin_fdfminimizer" ((type :pointer) (dimension sizet)) "multi-dimensional minimizer with function and derivative" ; FDL "Make an instance of a derivative-based minimizer of the given for an function of the given dimensions. Optionally initialize the minimizer to minimize the function starting from the initial point. The size of the first trial step is given by step-size. The accuracy of the line minimization is specified by tolernace. The precise meaning of this parameter depends on the method used. Typically the line minimization is considered successful if the gradient of the function g is orthogonal to the current search direction p to a relative accuracy of tolerance, where dot(p,g) < tol |p| |g|." "set"  Liam Healy committed Dec 26, 2008 83  ;; Could have one fewer argument: dimension=(dim0 initial)  Liam Healy committed Dec 25, 2008 84 85  ((function-derivative :pointer) ((mpointer initial) :pointer) (step-size :double) (tolerance :double)))  liam committed Jan 21, 2008 86   Liam Healy committed Dec 26, 2008 87 (defmfun name ((minimizer multi-dimensional-minimizer-f))  liam committed Jan 21, 2008 88  "gsl_multimin_fminimizer_name"  Liam Healy committed Dec 26, 2008 89 90  (((mpointer minimizer) :pointer)) :definition :method  liam committed Jan 21, 2008 91  :c-return :string  liam committed Feb 18, 2008 92  :documentation ; FDL  liam committed Jan 21, 2008 93 94  "The name of the minimizer.")  Liam Healy committed Dec 26, 2008 95 (defmfun name ((minimizer multi-dimensional-minimizer-fdf))  liam committed Jan 21, 2008 96  "gsl_multimin_fdfminimizer_name"  Liam Healy committed Dec 26, 2008 97 98  (((mpointer minimizer) :pointer)) :definition :method  liam committed Jan 21, 2008 99  :c-return :string  liam committed Feb 18, 2008 100  :documentation ; FDL  liam committed Jan 21, 2008 101 102 103 104 105 106  "The name of the minimizer.") ;;;;**************************************************************************** ;;;; Iteration ;;;;****************************************************************************  Liam Healy committed Jan 03, 2009 107 (defmfun iterate ((minimizer multi-dimensional-minimizer-f))  liam committed Jan 21, 2008 108  "gsl_multimin_fminimizer_iterate"  Liam Healy committed Dec 26, 2008 109  (((mpointer minimizer) :pointer))  Liam Healy committed Jan 03, 2009 110  :definition :method  liam committed Jan 21, 2008 111 112 113 114 115  :documentation ; FDL "Perform a single iteration of the minimizer. If the iteration encounters an unexpected problem then an error code will be returned.")  Liam Healy committed Jan 03, 2009 116 (defmfun iterate ((minimizer multi-dimensional-minimizer-fdf))  liam committed Jan 21, 2008 117  "gsl_multimin_fdfminimizer_iterate"  Liam Healy committed Dec 26, 2008 118  (((mpointer minimizer) :pointer))  Liam Healy committed Jan 03, 2009 119  :definition :method  liam committed Jan 21, 2008 120 121 122 123 124  :documentation ; FDL "Perform a single iteration of the minimizer. If the iteration encounters an unexpected problem then an error code will be returned.")  Liam Healy committed Jan 03, 2009 125 (defmfun solution ((minimizer multi-dimensional-minimizer-f))  liam committed Jan 21, 2008 126  "gsl_multimin_fminimizer_x"  Liam Healy committed Dec 26, 2008 127  (((mpointer minimizer) :pointer))  Liam Healy committed Jan 03, 2009 128  :definition :method  Liam Healy committed Dec 29, 2008 129 130  :c-return (crtn :pointer) :return ((make-marray 'double-float :from-pointer crtn))  liam committed Jan 21, 2008 131 132 133  :documentation ; FDL "The current best estimate of the location of the minimum.")  Liam Healy committed Jan 03, 2009 134 (defmfun solution ((minimizer multi-dimensional-minimizer-fdf))  liam committed Jan 21, 2008 135  "gsl_multimin_fdfminimizer_x"  Liam Healy committed Dec 26, 2008 136  (((mpointer minimizer) :pointer))  Liam Healy committed Jan 03, 2009 137  :definition :method  Liam Healy committed Dec 29, 2008 138 139  :c-return (crtn :pointer) :return ((make-marray 'double-float :from-pointer crtn))  liam committed Jan 21, 2008 140 141 142  :documentation ; FDL "The current best estimate of the location of the minimum.")  Liam Healy committed Jan 03, 2009 143 (defmfun function-value ((minimizer multi-dimensional-minimizer-f))  liam committed Jan 21, 2008 144  "gsl_multimin_fminimizer_minimum"  Liam Healy committed Dec 26, 2008 145  (((mpointer minimizer) :pointer))  Liam Healy committed Jan 03, 2009 146  :definition :method  liam committed Jan 21, 2008 147 148 149 150  :c-return :double :documentation ; FDL "The current best estimate of the value of the minimum.")  Liam Healy committed Jan 03, 2009 151 (defmfun function-value ((minimizer multi-dimensional-minimizer-fdf))  liam committed Jan 21, 2008 152  "gsl_multimin_fdfminimizer_minimum"  Liam Healy committed Dec 26, 2008 153  (((mpointer minimizer) :pointer))  Liam Healy committed Jan 03, 2009 154  :definition :method  liam committed Jan 21, 2008 155 156 157 158  :c-return :double :documentation ; FDL "The current best estimate of the value of the minimum.")  liam committed Feb 18, 2008 159 (defmfun mfminimizer-size (minimizer)  liam committed Jan 21, 2008 160  "gsl_multimin_fminimizer_size"  Liam Healy committed Dec 26, 2008 161  (((mpointer minimizer) :pointer))  liam committed Jan 21, 2008 162 163 164 165  :c-return :double :documentation ; FDL "A minimizer-specific characteristic size for the minimizer.")  liam committed Feb 18, 2008 166 (defmfun mfdfminimizer-gradient (minimizer)  liam committed Jan 21, 2008 167  "gsl_multimin_fdfminimizer_gradient"  Liam Healy committed Dec 26, 2008 168  (((mpointer minimizer) :pointer))  Liam Healy committed Dec 29, 2008 169 170  :c-return (crtn :pointer) :return ((make-marray 'double-float :from-pointer crtn))  liam committed Jan 21, 2008 171 172 173  :documentation ; FDL "The current best estimate of the gradient for the minimizer.")  liam committed Feb 18, 2008 174 (defmfun mfdfminimizer-restart (minimizer)  liam committed Jan 21, 2008 175  "gsl_multimin_fdfminimizer_restart"  Liam Healy committed Dec 26, 2008 176  (((mpointer minimizer) :pointer))  liam committed Jan 21, 2008 177 178 179 180 181 182 183 184  :documentation ; FDL "Reset the minimizer to use the current point as a new starting point.") ;;;;**************************************************************************** ;;;; Stopping criteria ;;;;****************************************************************************  liam committed Feb 18, 2008 185 (defmfun min-test-gradient (gradient absolute-error)  liam committed Jan 21, 2008 186  "gsl_multimin_test_gradient"  Liam Healy committed Dec 29, 2008 187  (((mpointer gradient) :pointer) (absolute-error :double))  liam committed Jan 21, 2008 188 189 190 191 192 193 194 195 196 197  :c-return :success-continue :documentation ; FDL "Test the norm of the gradient against the absolute tolerance absolute-error. The gradient of a multidimensional function goes to zero at a minimum. The test returns T if |g| < epsabs is achieved, and NIL otherwise. A suitable choice of absolute-error can be made from the desired accuracy in the function for small variations in x. The relationship between these quantities \delta f = g \delta x.")  liam committed Feb 18, 2008 198 (defmfun min-test-size (size absolute-error)  liam committed Jan 21, 2008 199 200 201 202 203  "gsl_multimin_test_size" ((size :double) (absolute-error :double)) :c-return :success-continue :documentation ; FDL "Test the minimizer specific characteristic size (if applicable to  liam committed Feb 18, 2008 204  the used minimizer) against absolute tolerance absolute-error. The test  liam committed Jan 21, 2008 205 206 207 208 209 210  returns T if the size is smaller than tolerance, and NIL otherwise.") ;;;;**************************************************************************** ;;;; Algorithms ;;;;****************************************************************************  liam committed Feb 18, 2008 211 (defmpar *conjugate-fletcher-reeves*  liam committed Jan 21, 2008 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228  "gsl_multimin_fdfminimizer_conjugate_fr" ;; FDL "The Fletcher-Reeves conjugate gradient algorithm. The conjugate gradient algorithm proceeds as a succession of line minimizations. The sequence of search directions is used to build up an approximation to the curvature of the function in the neighborhood of the minimum. An initial search direction p is chosen using the gradient, and line minimization is carried out in that direction. The accuracy of the line minimization is specified by the parameter tol. The minimum along this line occurs when the function gradient g and the search direction p are orthogonal. The line minimization terminates when dot(p,g) < tol |p| |g|. The search direction is updated using the Fletcher-Reeves formula p' = g' - \beta g where \beta=-|g'|^2/|g|^2, and the line minimization is then repeated for the new search direction.")  liam committed Feb 18, 2008 229 (defmpar *conjugate-polak-ribiere*  liam committed Jan 21, 2008 230 231 232 233 234 235 236 237  "gsl_multimin_fdfminimizer_conjugate_pr" ;; FDL "The Polak-Ribiere conjugate gradient algorithm. It is similar to the Fletcher-Reeves method, differing only in the choice of the coefficient \beta. Both methods work well when the evaluation point is close enough to the minimum of the objective function that it is well approximated by a quadratic hypersurface.")  liam committed Feb 18, 2008 238 (defmpar *vector-bfgs*  liam committed Jan 21, 2008 239 240 241 242 243 244 245 246 247 248  "gsl_multimin_fdfminimizer_vector_bfgs" ;; FDL "The vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) conjugate gradient algorithm. It is a quasi-Newton method which builds up an approximation to the second derivatives of the function using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region.")  liam committed Feb 18, 2008 249 (defmpar *simplex-nelder-mead*  liam committed Jan 21, 2008 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283  "gsl_multimin_fminimizer_nmsimplex" ;; FDL "The Simplex algorithm of Nelder and Mead. It constructs n vectors p_i from the starting vector initial and the vector step-size as follows: p_0 = (x_0, x_1, ... , x_n) p_1 = (x_0 + step_size_0, x_1, ... , x_n) p_2 = (x_0, x_1 + step_size_1, ... , x_n) ... = ... p_n = (x_0, x_1, ... , x_n+step_size_n) These vectors form the n+1 vertices of a simplex in n dimensions. On each iteration the algorithm tries to improve the parameter vector p_i corresponding to the highest function value by simple geometrical transformations. These are reflection, reflection followed by expansion, contraction and multiple contraction. Using these transformations the simplex moves through the parameter space towards the minimum, where it contracts itself. After each iteration, the best vertex is returned. Note, that due to the nature of the algorithm not every step improves the current best parameter vector. Usually several iterations are required. The routine calculates the minimizer specific characteristic size as the average distance from the geometrical center of the simplex to all its vertices. This size can be used as a stopping criteria, as the simplex contracts itself near the minimum. The size is returned by the function #'mfminimizer-size.") ;;;;**************************************************************************** ;;;; Examples ;;;;**************************************************************************** ;;; Examples from Sec. 35.8.  Liam Healy committed Jan 23, 2009 284 285 ;;; Example using derivatives, taking a vector argument.  liam committed Jan 21, 2008 286 287 288 (defparameter *parabaloid-center* #(1.0d0 2.0d0)) (defun parabaloid (gsl-vector-pointer)  Liam Healy committed Jan 23, 2009 289 290  "A parabaloid function of two arguments, given in GSL manual Sec. 35.4. This version takes a vector-double-float argument."  Liam Healy committed Aug 31, 2008 291 292  (let ((x (maref gsl-vector-pointer 0)) (y (maref gsl-vector-pointer 1))  liam committed Jan 21, 2008 293 294 295 296 297 298 299 300  (dp0 (aref *parabaloid-center* 0)) (dp1 (aref *parabaloid-center* 1))) (+ (* 10 (expt (- x dp0) 2)) (* 20 (expt (- y dp1) 2)) 30))) (defun parabaloid-derivative (arguments-gv-pointer derivative-gv-pointer)  Liam Healy committed Aug 31, 2008 301 302  (let ((x (maref arguments-gv-pointer 0)) (y (maref arguments-gv-pointer 1))  liam committed Jan 21, 2008 303 304  (dp0 (aref *parabaloid-center* 0)) (dp1 (aref *parabaloid-center* 1)))  Liam Healy committed Aug 31, 2008 305  (setf (maref derivative-gv-pointer 0)  liam committed Jan 21, 2008 306  (* 20 (- x dp0))  Liam Healy committed Aug 31, 2008 307  (maref derivative-gv-pointer 1)  liam committed Jan 21, 2008 308 309  (* 40 (- y dp1)))))  Liam Healy committed Jan 23, 2009 310 311 312 313 314 (defun parabaloid-and-derivative (arguments-gv-pointer derivative-gv-pointer) (prog1 (parabaloid arguments-gv-pointer) (parabaloid-derivative arguments-gv-pointer derivative-gv-pointer)))  liam committed Jan 21, 2008 315 316  (def-minimization-functions  Liam Healy committed Jan 23, 2009 317  parabaloid 2 parabaloid-derivative parabaloid-and-derivative t)  liam committed Jan 21, 2008 318   Liam Healy committed Jan 23, 2009 319 320 (defun multimin-example-derivative (&optional (method *conjugate-fletcher-reeves*) (print-steps t))  Liam Healy committed Dec 26, 2008 321  (let* ((initial #m(5.0d0 7.0d0))  lhealy committed Jul 25, 2008 322  (minimizer  Liam Healy committed Dec 26, 2008 323  (make-multi-dimensional-minimizer-fdf  Liam Healy committed Jan 23, 2009 324  method 2 parabaloid  Liam Healy committed Dec 26, 2008 325  initial 0.01d0 1.0d-4)))  liam committed Feb 04, 2008 326  (loop with status = T  Liam Healy committed Dec 26, 2008 327 328 329  for iter from 0 below 100 while status do  Liam Healy committed Jan 03, 2009 330  (iterate minimizer)  Liam Healy committed Dec 26, 2008 331 332 333 334  (setf status (not (min-test-gradient (mfdfminimizer-gradient minimizer) 1.0d-3)))  Liam Healy committed Jan 23, 2009 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354  (when print-steps (let ((x (solution minimizer))) (format t "~d~6t~10,6f~18t~10,6f~28t~12,9f~&" iter (maref x 0) (maref x 1) (function-value minimizer)))) finally (return (let ((x (solution minimizer))) (values (maref x 0) (maref x 1) (function-value minimizer))))))) ;;; Example without derivatives, same function but now defined using (defun parabaloid-f (x y) "A parabaloid function of two arguments, given in GSL manual Sec. 35.4. This version takes scalar arguments." (let ((dp0 (aref *parabaloid-center* 0)) (dp1 (aref *parabaloid-center* 1))) (+ (* 10 (expt (- x dp0) 2)) (* 20 (expt (- y dp1) 2)) 30)))  liam committed Jan 21, 2008 355 356 357  (def-minimization-functions parabaloid-f 2)  Liam Healy committed Jan 23, 2009 358 359 360 (defun multimin-example-no-derivative (&optional (method *simplex-nelder-mead*) (print-steps t)) (let ((step-size (make-marray 'double-float :dimensions 2)))  liam committed Feb 04, 2008 361  (set-all step-size 1.0d0)  Liam Healy committed Dec 26, 2008 362 363  (let ((minimizer (make-multi-dimensional-minimizer-f  Liam Healy committed Jan 23, 2009 364  method 2 parabaloid-f #m(5.0d0 7.0d0) step-size)))  liam committed Feb 04, 2008 365  (loop with status = T and size  Liam Healy committed Dec 26, 2008 366 367  for iter from 0 below 100 while status  Liam Healy committed Jan 03, 2009 368  do (iterate minimizer)  Liam Healy committed Dec 26, 2008 369 370 371 372  (setf size (mfminimizer-size minimizer) status (not (min-test-size size 1.0d-2)))  Liam Healy committed Jan 23, 2009 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388  (when print-steps (let ((x (solution minimizer))) (format t "~d~6t~10,6f~18t~10,6f~28t~12,9f~40t~8,3f~&" iter (maref x 0) (maref x 1) (function-value minimizer) size))) finally (return (let ((x (solution minimizer))) (values (maref x 0) (maref x 1) (function-value minimizer)))))))) (save-test minimization-multi (multimin-example-no-derivative *simplex-nelder-mead* nil) (multimin-example-derivative *conjugate-fletcher-reeves* nil) (multimin-example-derivative *conjugate-polak-ribiere* nil) (multimin-example-derivative *vector-bfgs* nil))