minimization-multi.lisp 15.1 KB
Newer Older
1 2
;; Multivariate minimization.
;; Liam Healy  <Tue Jan  8 2008 - 21:28>
3
;; Time-stamp: <2009-01-22 21:35:26EST minimization-multi.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

8 9
;;; /usr/include/gsl/gsl_multimin.h

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.
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)
31
(defmacro def-minimization-functions (function dimensions &optional df fdf array)
32 33 34
  "Setup functions for multivariate minimization.
   The CL functions name and derivative should be defined previously
   with defuns."
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)))))))
49 50 51 52 53

;;;;****************************************************************************
;;;; Initialization
;;;;****************************************************************************

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"
64
  ;; Could have one fewer argument: dimension=(dim0 initial)
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"
83
  ;; Could have one fewer argument: dimension=(dim0 initial)
84 85
  ((function-derivative :pointer) ((mpointer initial) :pointer)
   (step-size :double) (tolerance :double)))
86

87
(defmfun name ((minimizer multi-dimensional-minimizer-f))
88
  "gsl_multimin_fminimizer_name"
89 90
  (((mpointer minimizer) :pointer))
  :definition :method
91
  :c-return :string
liam's avatar
liam committed
92
  :documentation			; FDL
93 94
  "The name of the minimizer.")

95
(defmfun name ((minimizer multi-dimensional-minimizer-fdf))
96
  "gsl_multimin_fdfminimizer_name"
97 98
  (((mpointer minimizer) :pointer))
  :definition :method
99
  :c-return :string
liam's avatar
liam committed
100
  :documentation			; FDL
101 102 103 104 105 106
  "The name of the minimizer.")

;;;;****************************************************************************
;;;; Iteration
;;;;****************************************************************************

107
(defmfun iterate ((minimizer multi-dimensional-minimizer-f))
108
  "gsl_multimin_fminimizer_iterate"
109
  (((mpointer minimizer) :pointer))
110
  :definition :method
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.")

116
(defmfun iterate ((minimizer multi-dimensional-minimizer-fdf))
117
  "gsl_multimin_fdfminimizer_iterate"
118
  (((mpointer minimizer) :pointer))
119
  :definition :method
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.")

125
(defmfun solution ((minimizer multi-dimensional-minimizer-f))
126
  "gsl_multimin_fminimizer_x"
127
  (((mpointer minimizer) :pointer))
128
  :definition :method
129 130
  :c-return (crtn :pointer)
  :return ((make-marray 'double-float :from-pointer crtn))
131 132 133
  :documentation			; FDL
  "The current best estimate of the location of the minimum.")

134
(defmfun solution ((minimizer multi-dimensional-minimizer-fdf))
135
  "gsl_multimin_fdfminimizer_x"
136
  (((mpointer minimizer) :pointer))
137
  :definition :method
138 139
  :c-return (crtn :pointer)
  :return ((make-marray 'double-float :from-pointer crtn))
140 141 142
  :documentation			; FDL
  "The current best estimate of the location of the minimum.")

143
(defmfun function-value ((minimizer multi-dimensional-minimizer-f))
144
  "gsl_multimin_fminimizer_minimum"
145
  (((mpointer minimizer) :pointer))
146
  :definition :method
147 148 149 150
  :c-return :double
  :documentation			; FDL
  "The current best estimate of the value of the minimum.")

151
(defmfun function-value ((minimizer multi-dimensional-minimizer-fdf))
152
  "gsl_multimin_fdfminimizer_minimum"
153
  (((mpointer minimizer) :pointer))
154
  :definition :method
155 156 157 158
  :c-return :double
  :documentation			; FDL
  "The current best estimate of the value of the minimum.")

liam's avatar
liam committed
159
(defmfun mfminimizer-size (minimizer)
160
  "gsl_multimin_fminimizer_size"
161
  (((mpointer minimizer) :pointer))
162 163 164 165
  :c-return :double
  :documentation			; FDL
  "A minimizer-specific characteristic size for the minimizer.")

liam's avatar
liam committed
166
(defmfun mfdfminimizer-gradient (minimizer)
167
  "gsl_multimin_fdfminimizer_gradient"
168
  (((mpointer minimizer) :pointer))
169 170
  :c-return (crtn :pointer)
  :return ((make-marray 'double-float :from-pointer crtn))
171 172 173
  :documentation			; FDL
  "The current best estimate of the gradient for the minimizer.")

liam's avatar
liam committed
174
(defmfun mfdfminimizer-restart (minimizer)
175
  "gsl_multimin_fdfminimizer_restart"
176
  (((mpointer minimizer) :pointer))
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's avatar
liam committed
185
(defmfun min-test-gradient (gradient absolute-error)
186
  "gsl_multimin_test_gradient"
187
  (((mpointer gradient) :pointer) (absolute-error :double))
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's avatar
liam committed
198
(defmfun min-test-size (size absolute-error)
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's avatar
liam committed
204
   the used minimizer) against absolute tolerance absolute-error.  The test
205 206 207 208 209 210
   returns T if the size is smaller than tolerance, and NIL otherwise.")

;;;;****************************************************************************
;;;; Algorithms
;;;;****************************************************************************

liam's avatar
liam committed
211
(defmpar *conjugate-fletcher-reeves*
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's avatar
liam committed
229
(defmpar *conjugate-polak-ribiere*
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's avatar
liam committed
238
(defmpar *vector-bfgs*
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's avatar
liam committed
249
(defmpar *simplex-nelder-mead*
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.

284 285
;;; Example using derivatives, taking a vector argument.

286 287 288
(defparameter *parabaloid-center* #(1.0d0 2.0d0))

(defun parabaloid (gsl-vector-pointer)
289 290
  "A parabaloid function of two arguments, given in GSL manual Sec. 35.4.
   This version takes a vector-double-float argument."
291 292
  (let ((x (maref gsl-vector-pointer 0))
	(y (maref gsl-vector-pointer 1))
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)
301 302
  (let ((x (maref arguments-gv-pointer 0))
	(y (maref arguments-gv-pointer 1))
303 304
	(dp0 (aref *parabaloid-center* 0))
	(dp1 (aref *parabaloid-center* 1)))
305
    (setf (maref derivative-gv-pointer 0)
306
	  (* 20 (- x dp0))
307
	  (maref derivative-gv-pointer 1)
308 309
	  (* 40 (- y dp1)))))

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)))
315 316

(def-minimization-functions
317
    parabaloid 2 parabaloid-derivative parabaloid-and-derivative t)
318

319 320
(defun multimin-example-derivative
    (&optional (method *conjugate-fletcher-reeves*) (print-steps t))
321
  (let* ((initial #m(5.0d0 7.0d0))
322
	 (minimizer
323
	  (make-multi-dimensional-minimizer-fdf
324
	   method 2 parabaloid
325
	   initial 0.01d0 1.0d-4)))
326
    (loop with status = T
327 328 329
       for iter from 0 below 100
       while status
       do
330
       (iterate minimizer)
331 332 333 334
       (setf status
	     (not (min-test-gradient
		   (mfdfminimizer-gradient minimizer)
		   1.0d-3)))
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)))
355 356 357

(def-minimization-functions parabaloid-f 2)

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)))
361
    (set-all step-size 1.0d0)
362 363
    (let ((minimizer
	   (make-multi-dimensional-minimizer-f
364
	    method 2 parabaloid-f #m(5.0d0 7.0d0) step-size)))
365
      (loop with status = T and size
366 367
	 for iter from 0 below 100
	 while status
368
	 do (iterate minimizer)
369 370 371 372
	 (setf size
	       (mfminimizer-size minimizer)
	       status
	       (not (min-test-size size 1.0d-2)))
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))