roots-multi.lisp 18.2 KB
Newer Older
1 2
;;; Multivariate roots.                
;;; Liam Healy 2008-01-12 12:49:08
3
;;; Time-stamp: <2008-08-23 23:04:50EDT roots-multi.lisp>
lhealy's avatar
lhealy committed
4
;;; $Id$
5 6 7 8 9 10 11 12 13 14 15

(in-package :gsl)

;;;;****************************************************************************
;;;; Function definition
;;;;****************************************************************************

(cffi:defcstruct gsl-mfunction
  ;; See /usr/include/gsl/gsl_multiroots.h
  "The definition of a function for multiroot finding in GSL."
  (function :pointer)
16
  (dimensions sizet)
17 18 19 20 21
  (parameters :pointer))

(export 'def-mfunction)
(defmacro def-mfunction (name dimensions)
  "Define a function for multivariate root solving."
22
  `(def-single-function ,name :success-failure :pointer gsl-mfunction
23
    ((dimensions ,dimensions))
24
    (gsl-vector-c)))
25

26 27 28 29 30 31 32
(cffi:defcstruct gsl-mfunction-fdf
  ;; See /usr/include/gsl/gsl_multiroots.h
  "The definition of a function and its derivatives for multiroot
   finding in GSL."
  (function :pointer)
  (df :pointer)
  (fdf :pointer)
33
  (dimensions sizet)
34 35
  (parameters :pointer))

36 37 38 39
;;;;****************************************************************************
;;;; Initialization
;;;;****************************************************************************

40
(defgo mfsolver (type function-derivative initial)
41 42 43 44 45 46
  (list
   `(allocate-mfsolver ,type (dim0 ,initial))
   'free-mfsolver
   (lambda (symb)
     `(set-mfsolver ,symb ,function-derivative ,initial))))

47
(defgo mfdfsolver (type function-derivative initial)
48 49 50 51 52
  (list
   `(allocate-mfdfsolver ,type (dim0 ,initial))
   'free-mfdfsolver
   (lambda (symb)
     `(set-mfdfsolver ,symb ,function-derivative ,initial))))
53

liam's avatar
liam committed
54
(defmfun allocate-mfsolver (type dimension)
55
  "gsl_multiroot_fsolver_alloc"
56
  ((type :pointer) (dimension sizet))
57
  :c-return :pointer
58
  :export nil
59
  :index (letm mfsolver)
60
  :documentation			; FDL
61 62 63
  "Allocate an instance of a solver of the type specified for a system of
   the specified number of dimensions.")

liam's avatar
liam committed
64
(defmfun allocate-mfdfsolver (type dimension)
65
  "gsl_multiroot_fdfsolver_alloc"
66
  ((type :pointer) (dimension sizet))
67
  :c-return :pointer
68
  :export nil
69
  :index (letm mfdfsolver)
70
  :documentation			; FDL
71 72 73
  "Allocate an instance of a derivative solver of the type specified for
   a system of the specified number of dimensions.")

liam's avatar
liam committed
74
(defmfun set-mfsolver (solver function initial)
75
  "gsl_multiroot_fsolver_set"
76
  ((solver :pointer) (function :pointer) ((mpointer initial) :pointer))
77
  :export nil
78
  :index (letm mfsolver)
79
  :documentation			; FDL
80 81 82
  "Set or reset an existing solver to use the function and the
   initial guess gsl-vector.")

liam's avatar
liam committed
83
(defmfun set-mfdfsolver (solver function-derivative initial)
84 85
  "gsl_multiroot_fdfsolver_set"
  ((solver :pointer) (function-derivative :pointer)
86
   ((mpointer initial) :pointer))
87
  :export nil
88
  :index (letm mfdfsolver)
89
  :documentation			; FDL
90 91 92
  "Set or reset an existing solver to use the function and derivative
   (fdf) and the initial guess.")

liam's avatar
liam committed
93
(defmfun free-mfsolver (solver)
94 95 96
  "gsl_multiroot_fsolver_free"
  ((solver :pointer))
  :c-return :void
97
  :export nil
98
  :index (letm mfsolver)
99
  :documentation			; FDL
100 101
  "Free all the memory associated with the solver.")

liam's avatar
liam committed
102
(defmfun free-mfdfsolver (solver)
103 104 105
  "gsl_multiroot_fdfsolver_free"
  ((solver :pointer))
  :c-return :void
106
  :export nil
107
  :index (letm mfdfsolver)
108
  :documentation			; FDL
109 110
  "Free all the memory associated with the solver.")

liam's avatar
liam committed
111
(defmfun mfsolver-name (solver)
112 113 114
  "gsl_multiroot_fsolver_name"
  ((solver :pointer))
  :c-return :string
115
  :documentation			; FDL
116 117
  "The name of the solver.")

liam's avatar
liam committed
118
(defmfun mfdfsolver-name (solver)
119 120 121
  "gsl_multiroot_fdfsolver_name"
  ((solver :pointer))
  :c-return :string
122
  :documentation			; FDL
123 124 125 126 127 128
  "The name of the solver.")

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

liam's avatar
liam committed
129
(defmfun iterate-mfsolver (solver)
130 131
  "gsl_multiroot_fsolver_iterate"
  ((solver :pointer))
132
  :documentation			; FDL
133 134 135 136 137 138 139
  "Perform a single iteration of the solver.  The following
   errors may be signalled: :EBADFUNC,
   the iteration encountered a singular point where the function or its
   derivative evaluated to infinity or NaN, or
   :EZERODIV, the derivative of the function vanished at the iteration point,
   preventing the algorithm from continuing without a division by zero.")

liam's avatar
liam committed
140
(defmfun iterate-mfdfsolver (solver)
141 142
  "gsl_multiroot_fdfsolver_iterate"
  ((solver :pointer))
143
  :documentation			; FDL
144 145 146 147 148 149 150
  "Perform a single iteration of the solver.  The following
   errors may be signalled: :EBADFUNC,
   the iteration encountered a singular point where the function or its
   derivative evaluated to infinity or NaN, or
   :EZERODIV, the derivative of the function vanished at the iteration point,
   preventing the algorithm from continuing without a division by zero.")

liam's avatar
liam committed
151
(defmfun mfsolver-root (solver)
152 153
  "gsl_multiroot_fsolver_root"
  ((solver :pointer))
154 155
  :c-return :pointer
  :return (:c-return)
156
  :documentation			; FDL
157 158
  "The current estimate of the root for the solver.")

liam's avatar
liam committed
159
(defmfun mfdfsolver-root (solver)
160 161
  "gsl_multiroot_fdfsolver_root"
  ((solver :pointer))
162 163
  :c-return :pointer
  :return (:c-return)
164 165 166
  :documentation
  "The current estimate of the root for the solver.")

liam's avatar
liam committed
167
(defmfun mfsolver-f (solver)
168 169
  "gsl_multiroot_fsolver_f"
  ((solver :pointer))
170 171
  :c-return :pointer
  :return (:c-return)
172
  :documentation			; FDL
173 174
  "The function value f(x) at the current estimate x of the root for the solver.")

liam's avatar
liam committed
175
(defmfun mfdfsolver-f (solver)
176 177
  "gsl_multiroot_fdfsolver_f"
  ((solver :pointer))
178 179
  :c-return :pointer
  :return (:c-return)
180
  :documentation			; FDL
181 182
  "The function value f(x) at the current estimate x of the root for the solver.")

liam's avatar
liam committed
183
(defmfun mfsolver-dx (solver)
184 185
  "gsl_multiroot_fsolver_dx"
  ((solver :pointer))
186 187
  :c-return :pointer
  :return (:c-return)
188
  :documentation			; FDL
189 190
  "The last step dx taken by the solver.")

liam's avatar
liam committed
191
(defmfun mfdfsolver-dx (solver)
192 193
  "gsl_multiroot_fsolver_dx"
  ((solver :pointer))
194 195
  :c-return :pointer
  :return (:c-return)
196
  :documentation			; FDL
197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217
  "The last step dx taken by the solver.")

;;;;****************************************************************************
;;;; Search stopping conditions
;;;;****************************************************************************

;;; The only place we need to pick apart the gsl_multiroot_fsolver
;;; struct is here.  We could use mfsolver-dx etc., but then we'd have
;;; to discriminate on mfsolver vs. mfdfsolver.
(cffi:defcstruct gsl-multiroot-fsolver
  ;; See /usr/include/gsl/gsl_multiroots.h
  (type :pointer)
  (function :pointer)
  (x :pointer)
  (f :pointer)
  (dx :pointer)
  (state :pointer))

(defun multiroot-slot (solver slot)
  (cffi:foreign-slot-value solver 'gsl-multiroot-fsolver slot))

liam's avatar
liam committed
218
(defmfun multiroot-test-delta (solver absolute-error relative-error)
219 220 221 222 223
  "gsl_multiroot_test_delta"
  (((multiroot-slot solver 'dx) :pointer)
   ((multiroot-slot solver 'x) :pointer)
   (absolute-error :double) (relative-error :double))
  :c-return :success-continue
224
  :documentation			; FDL
225 226 227 228 229 230 231
  "Test for the convergence of the sequence by comparing the
   last step dx with the absolute error and relative
   errors given to the current position x.  The test returns
   T if the following condition is achieved:
   |dx_i| < epsabs + epsrel |x_i|
   for each component of x and returns NIL otherwise.")

liam's avatar
liam committed
232
(defmfun multiroot-test-residual (solver absolute-error)
233 234 235
  "gsl_multiroot_test_residual"
  (((multiroot-slot solver 'f) :pointer) (absolute-error :double))
  :c-return :success-failure
236
  :documentation			; FDL
237 238 239 240 241 242 243 244 245 246 247 248
  "Test the residual value f against the absolute error,
   returning T if the following condition is achieved:
   \sum_i |f_i| < absolute_error
   and returns NIL otherwise.  This criterion is suitable
   for situations where the precise location of the root x is
   unimportant provided a value can be found where the
   residual is small enough.")

;;;;****************************************************************************
;;;; Algorithms using derivatives
;;;;****************************************************************************

liam's avatar
liam committed
249
(defmpar *powells-hybrid* "gsl_multiroot_fdfsolver_hybridsj"
250
  ;; FDL
251
  "This is a modified version of Powell's Hybrid method as implemented in
liam's avatar
liam committed
252
   the hybrj algorithm in minpack.  Minpack was written by Jorge
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 284 285 286 287 288 289 290 291 292 293 294
   J. More, Burton S. Garbow and Kenneth E. Hillstrom.  The Hybrid
   algorithm retains the fast convergence of Newton's method but will also
   reduce the residual when Newton's method is unreliable. 

   The algorithm uses a generalized trust region to keep each step under
   control.  In order to be accepted a proposed new position x' must
   satisfy the condition |D (x' - x)| < \delta, where D is a
   diagonal scaling matrix and \delta is the size of the trust
   region.  The components of D are computed internally, using the
   column norms of the Jacobian to estimate the sensitivity of the residual
   to each component of x.  This improves the behavior of the
   algorithm for badly scaled functions.

   On each iteration the algorithm first determines the standard Newton
   step by solving the system J dx = - f.  If this step falls inside
   the trust region it is used as a trial step in the next stage.  If not,
   the algorithm uses the linear combination of the Newton and gradient
   directions which is predicted to minimize the norm of the function while
   staying inside the trust region,
   dx = - \alpha J^{-1} f(x) - \beta \nabla |f(x)|^2.
   This combination of Newton and gradient directions is referred to as a
   dogleg step.

   The proposed step is now tested by evaluating the function at the
   resulting point, x'.  If the step reduces the norm of the function
   sufficiently then it is accepted and size of the trust region is
   increased.  If the proposed step fails to improve the solution then the
   size of the trust region is decreased and another trial step is
   computed.

   The speed of the algorithm is increased by computing the changes to the
   Jacobian approximately, using a rank-1 update.  If two successive
   attempts fail to reduce the residual then the full Jacobian is
   recomputed.  The algorithm also monitors the progress of the solution
   and returns an error if several steps fail to make any improvement,
   :ENOPROG
   the iteration is not making any progress, preventing the algorithm from
   continuing.
   :ENOPROGJ
   re-evaluations of the Jacobian indicate that the iteration is not
   making any progress, preventing the algorithm from continuing.")

liam's avatar
liam committed
295
(defmpar *powells-hybrid-unscaled* "gsl_multiroot_fdfsolver_hybridj"
296
  ;; FDL
297 298 299 300 301
  "This algorithm is an unscaled version of *powells-hybrid*.  The steps are
   controlled by a spherical trust region |x' - x| < \delta, instead
   of a generalized region.  This can be useful if the generalized region
   estimated by *powells-hybrid* is inappropriate.")

liam's avatar
liam committed
302
(defmpar *newton-mfdfsolver* "gsl_multiroot_fdfsolver_newton"
303
  ;; FDL
304 305 306 307 308 309 310 311 312 313 314 315
  "Newton's Method is the standard root-polishing algorithm.  The algorithm
   begins with an initial guess for the location of the solution.  On each
   iteration a linear approximation to the function F is used to
   estimate the step which will zero all the components of the residual.
   The iteration is defined by the following sequence,
   x -> x' = x - J{-1} f(x)
   where the Jacobian matrix J is computed from the derivative
   functions provided by f.  The step dx is obtained by solving
   the linear system,
   J dx = - f(x)
   using LU decomposition.")

liam's avatar
liam committed
316
(defmpar *gnewton-mfdfsolver* "gsl_multiroot_fdfsolver_gnewton"
317
  ;; FDL
318 319 320 321 322 323 324 325 326 327 328 329 330
  "A modified version of Newton's method which attempts to improve
   global convergence by requiring every step to reduce the Euclidean norm
   of the residual, |f(x)|.  If the Newton step leads to an increase
   in the norm then a reduced step of relative size,
   t = (\sqrt(1 + 6 r) - 1) / (3 r)
   is proposed, with r being the ratio of norms
   |f(x')|^2/|f(x)|^2.  This procedure is repeated until a suitable step
   size is found.")

;;;;****************************************************************************
;;;; Algorithms without derivatives
;;;;****************************************************************************

liam's avatar
liam committed
331
(defmpar *hybrid-scaled* "gsl_multiroot_fsolver_hybrids"
332 333
  ;; FDL
  "This is a version of the Hybrid algorithm which replaces calls to the
334 335 336 337 338 339
     Jacobian function by its finite difference approximation.  The finite
     difference approximation is computed using gsl_multiroots_fdjac
     with a relative step size of GSL_SQRT_DBL_EPSILON.")
;; Where is this function and parameter?  Only thing that shows in the
;; library is gsl_multiroot_fdjacobian.
 
liam's avatar
liam committed
340
(defmpar *hybrid-unscaled* "gsl_multiroot_fsolver_hybrid"
341
  ;; FDL
342 343 344
  "A finite difference version of the Hybrid algorithm without
   internal scaling.")

liam's avatar
liam committed
345
(defmpar *discrete-newton* "gsl_multiroot_fsolver_dnewton"
346
  ;; FDL
347 348 349 350 351 352 353 354 355 356 357 358 359 360 361
  "The discrete Newton algorithm is the simplest method of solving a
   multidimensional system.  It uses the Newton iteration
   x -> x - J^{-1} f(x)
   where the Jacobian matrix J is approximated by taking finite
   differences of the function f.  The approximation scheme used by
   this implementation is
   J_{ij} = (f_i(x + \delta_j) - f_i(x)) /  \delta_j
   where \delta_j is a step of size \sqrt\epsilon |x_j| with
   \epsilon being the machine precision 
   (\epsilon \approx 2.22 \times 10^-16}).
   The order of convergence of Newton's algorithm is quadratic, but the
   finite differences require n^2 function evaluations on each
   iteration.  The algorithm may become unstable if the finite differences
   are not a good approximation to the true derivatives.")

liam's avatar
liam committed
362
(defmpar *broyden* "gsl_multiroot_fsolver_broyden"
363
  ;; FDL
364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388
  "The Broyden algorithm is a version of the discrete Newton
   algorithm which attempts to avoids the expensive update of the Jacobian
   matrix on each iteration.  The changes to the Jacobian are also
   approximated, using a rank-1 update,
   J^{-1} \to J^{-1} - (J^{-1} df - dx) dx^T J^{-1} / dx^T J^{-1} df
   where the vectors dx and df are the changes in x
   and f.  On the first iteration the inverse Jacobian is estimated
   using finite differences, as in the discrete Newton algorithm.
    
   This approximation gives a fast update but is unreliable if the changes
   are not small, and the estimate of the inverse Jacobian becomes worse as
   time passes.  The algorithm has a tendency to become unstable unless it
   starts close to the root.  The Jacobian is refreshed if this instability
   is detected (consult the source for details).

   This algorithm is included only for demonstration purposes, and is not
   recommended for serious use.")

;;;;****************************************************************************
;;;; Examples
;;;;****************************************************************************

(defparameter *powell-A* 1.0d4)
(defun powell (argument return)
  "Powell's test function."
389 390
  (setf (maref return 0)
	(- (* *powell-A* (maref argument 0) (maref argument 1))
391
	   1)
392 393
	(maref return 1)
	(+ (exp (- (maref argument 0))) (exp (- (maref argument 1)))
394 395 396 397 398 399 400 401 402
	   (- (1+ (/ *powell-A*))))))

;;; (def-mfunction powell 2)

;;; This is the example given in Sec. 34.8.

(defparameter *rosenbrock-a* 1.0d0)
(defparameter *rosenbrock-b* 10.0d0)

403
#|
404
;;; One alternative way of writing the function, not recommended.
405 406
(defun rosenbrock (argument return)
  "Rosenbrock test function."
407 408
  (with-c-doubles (((vector-data argument) x0 x1)
		   ((vector-data return) f0 f1))
409 410
    (setf f0 (* *rosenbrock-a* (- 1 x0))
	  f1 (* *rosenbrock-b* (- x1 (expt x0 2))))))
411
|#
412

413 414 415
;;; The recommended alternative
(defun rosenbrock (argument return)
  "Rosenbrock test function."
416 417 418 419
  (setf (maref return 0)
	(* *rosenbrock-a* (- 1 (maref argument 0)))
	(maref return 1)
	(* *rosenbrock-b* (- (maref argument 1) (expt (maref argument 0) 2)))))
420

421 422 423 424 425
(def-mfunction rosenbrock 2)

(defun roots-multi-example ()
  "Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual."
  (let ((max-iter 1000))
426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447
    (letm ((vect (vector-double-float (a -10.0d0 -5.0d0)))
	   (solver (mfsolver *hybrid-scaled* rosenbrock vect)))
      (loop for iter from 0
	 with fnval and argval
	 while (and (< iter max-iter)
		    (not (multiroot-test-residual solver 1.0d-7)))
	   
	 do
	 (iterate-mfsolver solver)
	 (setf fnval (cl-array (mfsolver-f solver))
	       argval (cl-array (mfsolver-root solver)))
	 (format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g"
		 iter
		 (aref argval 0)
		 (aref argval 1)
		 (aref fnval 0)
		 (aref fnval 1))
	 finally (return
		   (values (aref argval 0)
			   (aref argval 1)
			   (aref fnval 0)
			   (aref fnval 1)))))))
448

449 450
(defun rosenbrock-df (argument jacobian)
  "The partial derivatives of the Rosenbrock functions."
451 452 453 454
  (setf (maref jacobian 0 0) (- *rosenbrock-a*)
	(maref jacobian 0 1) 0.0d0
	(maref jacobian 1 0) (* -2 *rosenbrock-b* (maref argument 0))
	(maref jacobian 1 1) *rosenbrock-b*))
455 456 457 458 459

(defun rosenbrock-fdf (argument value jacobian)
  (rosenbrock argument value)
  (rosenbrock-df argument jacobian))

460
;;; Because def-solver-functions and def-single-function bind a symbol
461 462 463 464 465 466 467 468 469 470 471 472 473
;;; of the same name as the first function, and we want both to run,
;;; we'll make an alias function so we can use both.  
(eval-when (:load-toplevel :execute)
  (setf (fdefinition 'rosenbrock-f) #'rosenbrock))

(def-solver-functions rosenbrock-f rosenbrock-df rosenbrock-fdf 2)

(defun roots-multi-example-df ()
  "Solving Rosenbrock with derivatives, the example given in Sec. 34.8
   of the GSL manual."
  (flet ((print-state (iter argval fnval)
	   (format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g"
		   iter
474 475 476 477
		   (maref argval 0)
		   (maref argval 1)
		   (maref fnval 0)
		   (maref fnval 1))))
478
    (let ((max-iter 1000))
479
      (letm ((vect (vector-double-float (a -10.0d0 -5.0d0))))
480 481
	(letm
	    ((solver (mfdfsolver *gnewton-mfdfsolver* rosenbrock-f vect)))
482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497
	  (loop for iter from 0
	     with fnval = (mfdfsolver-f solver)
	     and argval = (mfdfsolver-root solver)
	     while (and (< iter max-iter)
			(not (multiroot-test-residual solver 1.0d-7)))
	     initially (print-state iter argval fnval)
	     do
	     (iterate-mfdfsolver solver)
	     (setf fnval (mfdfsolver-f solver)
		   argval (mfdfsolver-root solver))
	     (print-state iter argval fnval)
	     finally (return
		       (values (maref argval 0)
			       (maref argval 1)
			       (maref fnval 0)
			       (maref fnval 1)))))))))