summaryrefslogtreecommitdiffstats
path: root/minimize.scm
diff options
context:
space:
mode:
Diffstat (limited to 'minimize.scm')
-rw-r--r--minimize.scm114
1 files changed, 114 insertions, 0 deletions
diff --git a/minimize.scm b/minimize.scm
new file mode 100644
index 0000000..50a7e65
--- /dev/null
+++ b/minimize.scm
@@ -0,0 +1,114 @@
+;;; "minimize.scm" finds minimum f(x) for x0 <= x <= x1.
+;;; Author: Lars Arvestad
+;;;
+;;; This code is in the public domain.
+
+;;@noindent
+;;
+;;The Golden Section Search
+;;@footnote{David Kahaner, Cleve Moler, and Stephen Nash
+;;@cite{Numerical Methods and Software}
+;;Prentice-Hall, 1989, ISBN 0-13-627258-4}
+;;algorithm finds minima of functions which
+;;are expensive to compute or for which derivatives are not available.
+;;Although optimum for the general case, convergence is slow,
+;;requiring nearly 100 iterations for the example (x^3-2x-5).
+;;
+;;@noindent
+;;
+;;If the derivative is available, Newton-Raphson is probably a better
+;;choice. If the function is inexpensive to compute, consider
+;;approximating the derivative.
+
+;;@body
+;;
+;;@var{x_0} are @var{x_1} real numbers. The (single argument)
+;;procedure @var{f} is unimodal over the open interval (@var{x_0},
+;;@var{x_1}). That is, there is exactly one point in the interval for
+;;which the derivative of @var{f} is zero.
+;;
+;;@0 returns a pair (@var{x} . @var{f}(@var{x})) where @var{f}(@var{x})
+;;is the minimum. The @var{prec} parameter is the stop criterion. If
+;;@var{prec} is a positive number, then the iteration continues until
+;;@var{x} is within @var{prec} from the true value. If @var{prec} is
+;;a negative integer, then the procedure will iterate @var{-prec}
+;;times or until convergence. If @var{prec} is a procedure of seven
+;;arguments, @var{x0}, @var{x1}, @var{a}, @var{b}, @var{fa}, @var{fb},
+;;and @var{count}, then the iterations will stop when the procedure
+;;returns @code{#t}.
+;;
+;;Analytically, the minimum of x^3-2x-5 is 0.816497.
+;;@example
+;;(define func (lambda (x) (+ (* x (+ (* x x) -2)) -5)))
+;;(golden-section-search func 0 1 (/ 10000))
+;; ==> (816.4883855245578e-3 . -6.0886621077391165)
+;;(golden-section-search func 0 1 -5)
+;; ==> (819.6601125010515e-3 . -6.088637561916407)
+;;(golden-section-search func 0 1
+;; (lambda (a b c d e f g ) (= g 500)))
+;; ==> (816.4965933140557e-3 . -6.088662107903635)
+;;@end example
+
+(define golden-section-search
+ (let ((gss 'golden-section-search:)
+ (r (/ (- (sqrt 5) 1) 2))) ; 1 / golden-section
+ (lambda (f x0 x1 prec)
+ (cond ((not (procedure? f)) (slib:error gss 'procedure? f))
+ ((not (number? x0)) (slib:error gss 'number? x0))
+ ((not (number? x1)) (slib:error gss 'number? x1))
+ ((>= x0 x1) (slib:error gss x0 'not '< x1)))
+ (let ((stop?
+ (cond
+ ((procedure? prec) prec)
+ ((number? prec)
+ (if (>= prec 0)
+ (lambda (x0 x1 a b fa fb count) (<= (abs (- x1 x0)) prec))
+ (if (integer? prec)
+ (lambda (x0 x1 a b fa fb count) (>= count (- prec)))
+ (slib:error gss 'integer? prec))))
+ (else (slib:error gss 'procedure? prec))))
+ (a0 (+ x0 (* (- x1 x0) (- 1 r))))
+ (b0 (+ x0 (* (- x1 x0) r)))
+ (delta #f)
+ (fmax #f)
+ (fmin #f))
+ (let loop ((left x0)
+ (right x1)
+ (a a0)
+ (b b0)
+ (fa (f a0))
+ (fb (f b0))
+ (count 1))
+ (define finish
+ (lambda (x fx)
+ (if (> fx fmin) (slib:warn gss fx 'not 'min (list '> fmin)))
+ (if (and (> count 9) (or (eqv? x0 left) (eqv? x1 right)))
+ (slib:warn gss 'min 'not 'found))
+ (cons x fx)))
+ (case count
+ ((1)
+ (set! fmax (max fa fb))
+ (set! fmin (min fa fb)))
+ ((2)
+ (set! fmin (min fmin fa fb))
+ (if (eqv? fmax fa fb) (slib:error gss 'flat? fmax)))
+ (else
+ (set! fmin (min fmin fa fb))))
+ (cond ((stop? left right a b fa fb count)
+ (if (< fa fb)
+ (finish a fa)
+ (finish b fb)))
+ ((< fa fb)
+ (let ((a-next (+ left (* (- b left) (- 1 r)))))
+ (cond ((and delta (< delta (- b a)))
+ (finish a fa))
+ (else (set! delta (- b a))
+ (loop left b a-next a (f a-next) fa
+ (+ 1 count))))))
+ (else
+ (let ((b-next (+ a (* (- right a) r))))
+ (cond ((and delta (< delta (- b a)))
+ (finish b fb))
+ (else (set! delta (- b a))
+ (loop a right b b-next fb (f b-next)
+ (+ 1 count))))))))))))