summaryrefslogtreecommitdiffstats
path: root/minimize.scm
blob: e28568afb3262142b23b477df5696a9fb2e59f08 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
;;; "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 (= 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))))))))))))