roots-multi.lisp 20 KB
Newer Older
1 2
;;; Multivariate roots.                
;;; Liam Healy 2008-01-12 12:49:08
3
;;; Time-stamp: <2012-01-13 12:01:17EST roots-multi.lisp>
4
;;
5
;; Copyright 2008, 2009, 2011 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; 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 3 of the License, or
;; (at your option) any later version.
;;
;; This program 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.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20 21 22

(in-package :gsl)

23 24
;;; /usr/include/gsl/gsl_multiroots.h

25 26
;;; Currently, functions defined for root solving will be passed
;;; scalars and should return scalars as multiple values.  A possible
27 28 29
;;; future enhancement is to optionally pass grid:foreign-arrays and return
;;; grid:foreign-arrays instead.  This would allow directly manipulation of
;;; grid:foreign-arrays by the user function.  Notes Mon Jan 19 2009.
30

31 32 33 34
;;;;****************************************************************************
;;;; Initialization
;;;;****************************************************************************

35
(defmobject multi-dimensional-root-solver-f "gsl_multiroot_fsolver"
36
  ((type :pointer) ((first dimensions) :sizet))
37 38
  "multi-dimensional root solver with function only"
  :documentation			; FDL
39 40 41
  "Make an instance of a solver of the type specified for a system of
   the specified number of dimensions.  Optionally
   set or reset an existing solver to use the function and the
42
   initial guess gsl-vector.  If scalarsp is T, the functions will
43
   be supplied scalars, and should return scalars."
44
  :initialize-suffix "set"
45
  :initialize-args ((callback :pointer) ((mpointer initial) :pointer))
46
  :callbacks
47
  (callback (:struct fnstruct-dimension) (dimension)
48 49
	    (function
	     :success-failure
50 51
	     (:input :double :foreign-array dim0) :slug
	     (:output :double :foreign-array dim0)))
52
  :arglists-function
53
  (lambda (set)
54
    `((type &optional function-or-dimension (initial nil ,set) (scalarsp t))
55
      (:type type
56
	     :dimensions
57
	     (if ,set (grid:dimensions initial) function-or-dimension))
58 59
      (:functions
       (list function-or-dimension) :initial initial :scalarsp scalarsp)))
60
  :inputs (initial))
61 62

(defmobject multi-dimensional-root-solver-fdf "gsl_multiroot_fdfsolver"
63
  ((type :pointer) ((first dimensions) :sizet))
64 65
  "multi-dimensional root solver with function and derivative"
  :documentation			; FDL
66 67 68
  "Make an instance of a derivative solver of the type specified for
   a system of the specified number of dimensions.  Optionally
   set or reset an existing solver to use the function and derivative
69 70
   (fdf) and the initial guess.  If scalarsp is T, the functions will
   be supplied, and should return scalars."
71
  :initialize-suffix "set"
72
  :initialize-args ((callback :pointer) ((mpointer initial) :pointer))
73
  :callbacks
74
  (callback (:struct fnstruct-dimension-fdf) (dimension)
75
	    (function :success-failure
76
		      (:input :double :foreign-array dim0)
77
		      :slug
78
		      (:output :double :foreign-array dim0))
79
	    (df :success-failure
80
		(:input :double :foreign-array dim0)
81
		:slug
82
		(:output :double :foreign-array dim0 dim0))
83
	    (fdf :success-failure
84
		 (:input :double :foreign-array dim0)
85
		 :slug
86 87
		 (:output :double :foreign-array dim0)
		 (:output :double :foreign-array dim0 dim0)))
88
  :arglists-function
89
  (lambda (set)
90 91
    `((type &optional function-or-dimension (initial nil ,set) 
	    (scalarsp t))
92 93
      (:type type
	     :dimensions
94
	     (if ,set (grid:dimensions initial) function-or-dimension))
95
      (:functions function-or-dimension :initial initial :scalarsp scalarsp)))
96
  :inputs (initial))
97

98
(defmfun name ((solver multi-dimensional-root-solver-f))
99
  "gsl_multiroot_fsolver_name"
100 101
  (((mpointer solver) :pointer))
  :definition :method
102
  :c-return :string
103
  :documentation			; FDL
104 105
  "The name of the solver.")

106
(defmfun name ((solver multi-dimensional-root-solver-fdf))
107
  "gsl_multiroot_fdfsolver_name"
108 109
  (((mpointer solver) :pointer))
  :definition :method
110
  :c-return :string
111
  :documentation			; FDL
112 113 114 115 116 117
  "The name of the solver.")

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

118
(defmfun iterate ((solver multi-dimensional-root-solver-f))
119
  "gsl_multiroot_fsolver_iterate"
120
  (((mpointer solver) :pointer))
121
  :definition :method
122
  :callback-object solver
123
  :documentation			; FDL
124 125 126 127 128 129
  "Perform a single iteration of the solver.  The following errors may
   be signalled: 'bad-function-supplied, the iteration encountered a
   singular point where the function or its derivative evaluated to
   infinity or NaN, or 'gsl-division-by-zero, the derivative of the
   function vanished at the iteration point, preventing the algorithm
   from continuing without a division by zero.")
130

131
(defmfun iterate ((solver multi-dimensional-root-solver-fdf))
132
  "gsl_multiroot_fdfsolver_iterate"
133
  (((mpointer solver) :pointer))
134
  :definition :method
135
  :callback-object solver
136
  :documentation			; FDL
137 138 139 140 141 142
  "Perform a single iteration of the solver.  The following errors may
   be signalled: 'bad-function-supplied, the iteration encountered a
   singular point where the function or its derivative evaluated to
   infinity or NaN, or 'gsl-division-by-zero, the derivative of the
   function vanished at the iteration point, preventing the algorithm
   from continuing without a division by zero.")
143

144
(defmfun solution ((solver multi-dimensional-root-solver-f))
145
  "gsl_multiroot_fsolver_root"
146
  (((mpointer solver) :pointer))
147
  :definition :method
148
  :callback-object solver
149
  :c-return (crtn :pointer)
150
  :return ((make-foreign-array-from-mpointer crtn))
151
  :documentation			; FDL
152 153
  "The current estimate of the root for the solver.")

154
(defmfun solution ((solver multi-dimensional-root-solver-fdf))
155
  "gsl_multiroot_fdfsolver_root"
156
  (((mpointer solver) :pointer))
157
  :definition :method
158
  :callback-object solver
159
  :c-return (crtn :pointer)
160
  :return ((make-foreign-array-from-mpointer crtn))
161 162 163
  :documentation
  "The current estimate of the root for the solver.")

164
(defmfun function-value ((solver multi-dimensional-root-solver-f))
165
  "gsl_multiroot_fsolver_f"
166
  (((mpointer solver) :pointer))
167
  :definition :method
168
  :c-return (crtn :pointer)
169
  :return ((make-foreign-array-from-mpointer crtn))
170
  :documentation			; FDL
171 172
  "The function value f(x) at the current estimate x of the root for the solver.")

173
(defmfun function-value ((solver multi-dimensional-root-solver-fdf))
174
  "gsl_multiroot_fdfsolver_f"
175
  (((mpointer solver) :pointer))
176
  :definition :method
177
  :c-return (crtn :pointer)
178
  :return ((make-foreign-array-from-mpointer crtn))
179
  :documentation			; FDL
180 181
  "The function value f(x) at the current estimate x of the root for the solver.")

182
(defmfun last-step ((solver multi-dimensional-root-solver-f))
183
  "gsl_multiroot_fsolver_dx"
184
  (((mpointer solver) :pointer))
185
  :definition :method
186
  :c-return (crtn :pointer)
187
  :return ((make-foreign-array-from-mpointer crtn))
188
  :documentation			; FDL
189 190
  "The last step dx taken by the solver.")

191
(defmfun last-step ((solver multi-dimensional-root-solver-fdf))
192
  "gsl_multiroot_fsolver_dx"
193
  (((mpointer solver) :pointer))
194
  :definition :method
195
  :c-return (crtn :pointer)
196
  :return ((make-foreign-array-from-mpointer crtn))
197
  :documentation			; FDL
198 199 200 201 202 203 204
  "The last step dx taken by the solver.")

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

;;; The only place we need to pick apart the gsl_multiroot_fsolver
205
;;; struct is here.  We could use last-step etc., but then we'd have
206 207 208
;;; to discriminate on mfsolver vs. mfdfsolver.

(defun multiroot-slot (solver slot)
209
  (cffi:foreign-slot-value (mpointer solver) '(:struct gsl-multiroot-fsolver) slot))
210

liam's avatar
liam committed
211
(defmfun multiroot-test-delta (solver absolute-error relative-error)
212 213 214 215 216
  "gsl_multiroot_test_delta"
  (((multiroot-slot solver 'dx) :pointer)
   ((multiroot-slot solver 'x) :pointer)
   (absolute-error :double) (relative-error :double))
  :c-return :success-continue
217
  :documentation			; FDL
218 219 220 221 222 223 224
  "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
225
(defmfun multiroot-test-residual (solver absolute-error)
226 227
  "gsl_multiroot_test_residual"
  (((multiroot-slot solver 'f) :pointer) (absolute-error :double))
228
  :c-return :success-continue
229
  :documentation			; FDL
230 231 232 233 234 235 236 237 238 239 240 241
  "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 Healy's avatar
Liam Healy committed
242
(defmpar +powells-hybrid+ "gsl_multiroot_fdfsolver_hybridsj"
243
  ;; FDL
244
  "This is a modified version of Powell's Hybrid method as implemented in
liam's avatar
liam committed
245
   the hybrj algorithm in minpack.  Minpack was written by Jorge
246 247 248 249 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
   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,
281
   'no-progress
282 283
   the iteration is not making any progress, preventing the algorithm from
   continuing.
284
   'jacobian-not-improving
285 286 287
   re-evaluations of the Jacobian indicate that the iteration is not
   making any progress, preventing the algorithm from continuing.")

Liam Healy's avatar
Liam Healy committed
288
(defmpar +powells-hybrid-unscaled+ "gsl_multiroot_fdfsolver_hybridj"
289
  ;; FDL
290 291 292 293 294
  "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 Healy's avatar
Liam Healy committed
295
(defmpar +newton-mfdfsolver+ "gsl_multiroot_fdfsolver_newton"
296
  ;; FDL
297 298 299 300 301 302 303 304 305 306 307 308
  "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 Healy's avatar
Liam Healy committed
309
(defmpar +gnewton-mfdfsolver+ "gsl_multiroot_fdfsolver_gnewton"
310
  ;; FDL
311 312 313 314 315 316 317 318 319 320 321 322 323
  "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 Healy's avatar
Liam Healy committed
324
(defmpar +hybrid-scaled+ "gsl_multiroot_fsolver_hybrids"
325 326
  ;; FDL
  "This is a version of the Hybrid algorithm which replaces calls to the
327 328 329 330 331 332
     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 Healy's avatar
Liam Healy committed
333
(defmpar +hybrid-unscaled+ "gsl_multiroot_fsolver_hybrid"
334
  ;; FDL
335 336 337
  "A finite difference version of the Hybrid algorithm without
   internal scaling.")

Liam Healy's avatar
Liam Healy committed
338
(defmpar +discrete-newton+ "gsl_multiroot_fsolver_dnewton"
339
  ;; FDL
340 341 342 343 344 345 346 347 348 349 350 351 352 353 354
  "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 Healy's avatar
Liam Healy committed
355
(defmpar +broyden+ "gsl_multiroot_fsolver_broyden"
356
  ;; FDL
357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378
  "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
;;;;****************************************************************************

379 380 381 382 383 384 385
;;; This is the example given in GSL manual, Sec. 35.8.
;;; http://www.gnu.org/software/gsl/manual/html_node/Example-programs-for-Multidimensional-Root-finding.html

;;; These examples use use scalarsp=T in the
;;; multi-dimensional-root-solver argument.  To see how vectors would
;;; be used, see minimization-multi.

386
(defparameter *powell-A* 1.0d4)
387
(defun powell (arg0 arg1)
388
  "Powell's test function."
389 390 391
  (values
   (- (* *powell-A* arg0 arg1) 1)
   (+ (exp (- arg0)) (exp (- arg1)) (- (1+ (/ *powell-A*))))))
392
;; not used?
393 394 395 396

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

397
(defun rosenbrock (arg0 arg1)
398
  "Rosenbrock test function."
399 400 401
  (values
   (* *rosenbrock-a* (- 1 arg0))
   (* *rosenbrock-b* (- arg1 (expt arg0 2)))))
402

403
(defun roots-multi-example-no-derivative
Liam Healy's avatar
Liam Healy committed
404
    (&optional (method +hybrid-scaled+) (print-steps t))
405 406 407
  "Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual."
  (let ((max-iter 1000)
	(solver (make-multi-dimensional-root-solver-f 
408 409 410
		 method 'rosenbrock
		 (grid:make-foreign-array
		  'double-float :initial-contents '(-10.0d0 -5.0d0)))))
411
    (loop for iter from 0
412 413 414 415 416
	  with fnval and argval
	  while (and (< iter max-iter)
		     (or (zerop iter)
			 (not (multiroot-test-residual solver 1.0d-7))))
	  do
417 418 419 420 421 422
       (iterate solver)
       (setf fnval (function-value solver)
	     argval (solution solver))
       (when print-steps
	 (format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&"
		 iter
423 424 425 426
		 (grid:aref argval 0)
		 (grid:aref argval 1)
		 (grid:aref fnval 0)
		 (grid:aref fnval 1)))
427 428 429 430 431
	  finally (return
		    (values (grid:aref argval 0)
			    (grid:aref argval 1)
			    (grid:aref fnval 0)
			    (grid:aref fnval 1))))))
432

433
(defun rosenbrock-df (arg0 arg1)
434
  "The partial derivatives of the Rosenbrock functions."
435 436 437 438 439 440 441
  (declare (ignore arg1))
  (values (- *rosenbrock-a*)
	  0.0d0
	  (* -2 *rosenbrock-b* arg0)
	  *rosenbrock-b*))

;;; Why is it necessary to define a function that calls the two other functions?
442
(defun rosenbrock-fdf (arg0 arg1)
443 444 445 446 447
  (multiple-value-bind (v0 v1)
      (rosenbrock arg0 arg1)
    (multiple-value-bind (j0 j1 j2 j3)
	(rosenbrock-df arg0 arg1)
      (values v0 v1 j0 j1 j2 j3))))
448

449
(defun roots-multi-example-derivative
Liam Healy's avatar
Liam Healy committed
450
    (&optional (method +gnewton-mfdfsolver+) (print-steps t))
451 452 453
  "Solving Rosenbrock with derivatives, the example given in Sec. 34.8
   of the GSL manual."
  (flet ((print-state (iter argval fnval)
454 455 456
	   (when print-steps
	     (format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&"
		     iter
457 458 459 460
		     (grid:aref argval 0)
		     (grid:aref argval 1)
		     (grid:aref fnval 0)
		     (grid:aref fnval 1)))))
461 462
    (let ((max-iter 1000)
	  (solver (make-multi-dimensional-root-solver-fdf
463
		   method
464
		   '(rosenbrock rosenbrock-df rosenbrock-fdf)
465 466
		   (grid:make-foreign-array
		    'double-float :initial-contents '(-10.0d0 -5.0d0)))))
467
      (loop for iter from 0
468 469 470 471 472 473
	    with fnval = (function-value solver)
	    and argval = (solution solver)
	    while (and (< iter max-iter)
		       (not (multiroot-test-residual solver 1.0d-7)))
	      initially (print-state iter argval fnval)
	    do
474 475 476 477
	 (iterate solver)
	 (setf fnval (function-value solver)
	       argval (solution solver))
	 (print-state iter argval fnval)
478 479 480 481 482
	    finally (return
		      (values (grid:aref argval 0)
			      (grid:aref argval 1)
			      (grid:aref fnval 0)
			      (grid:aref fnval 1)))))))
483

484 485
;; To see step-by-step information as the solution progresses, make
;; the last argument T.
486
(save-test roots-multi
Liam Healy's avatar
Liam Healy committed
487 488 489 490 491 492 493 494
 (roots-multi-example-no-derivative +hybrid-unscaled+ nil)
 (roots-multi-example-no-derivative +hybrid-scaled+ nil)
 (roots-multi-example-no-derivative +discrete-newton+ nil)
 (roots-multi-example-no-derivative +broyden+ nil)
 (roots-multi-example-derivative +newton-mfdfsolver+ nil)
 (roots-multi-example-derivative +gnewton-mfdfsolver+ nil)
 (roots-multi-example-derivative +powells-hybrid+ nil)
 (roots-multi-example-derivative +powells-hybrid-unscaled+ nil))