wolffd@0
|
1 function [br_min, br_mid, br_max, num_evals] = minbrack(f, a, b, fa, ...
|
wolffd@0
|
2 varargin)
|
wolffd@0
|
3 %MINBRACK Bracket a minimum of a function of one variable.
|
wolffd@0
|
4 %
|
wolffd@0
|
5 % Description
|
wolffd@0
|
6 % BRMIN, BRMID, BRMAX, NUMEVALS] = MINBRACK(F, A, B, FA) finds a
|
wolffd@0
|
7 % bracket of three points around a local minimum of F. The function F
|
wolffd@0
|
8 % must have a one dimensional domain. A < B is an initial guess at the
|
wolffd@0
|
9 % minimum and maximum points of a bracket, but MINBRACK will search
|
wolffd@0
|
10 % outside this interval if necessary. The bracket consists of three
|
wolffd@0
|
11 % points (in increasing order) such that F(BRMID) < F(BRMIN) and
|
wolffd@0
|
12 % F(BRMID) < F(BRMAX). FA is the value of the function at A: it is
|
wolffd@0
|
13 % included to avoid unnecessary function evaluations in the
|
wolffd@0
|
14 % optimization routines. The return value NUMEVALS is the number of
|
wolffd@0
|
15 % function evaluations in MINBRACK.
|
wolffd@0
|
16 %
|
wolffd@0
|
17 % MINBRACK(F, A, B, FA, P1, P2, ...) allows additional arguments to be
|
wolffd@0
|
18 % passed to F
|
wolffd@0
|
19 %
|
wolffd@0
|
20 % See also
|
wolffd@0
|
21 % LINEMIN, LINEF
|
wolffd@0
|
22 %
|
wolffd@0
|
23
|
wolffd@0
|
24 % Copyright (c) Ian T Nabney (1996-2001)
|
wolffd@0
|
25
|
wolffd@0
|
26 % Check function string
|
wolffd@0
|
27 f = fcnchk(f, length(varargin));
|
wolffd@0
|
28
|
wolffd@0
|
29 % Value of golden section (1 + sqrt(5))/2.0
|
wolffd@0
|
30 phi = 1.6180339887499;
|
wolffd@0
|
31
|
wolffd@0
|
32 % Initialise count of number of function evaluations
|
wolffd@0
|
33 num_evals = 0;
|
wolffd@0
|
34
|
wolffd@0
|
35 % A small non-zero number to avoid dividing by zero in quadratic interpolation
|
wolffd@0
|
36 TINY = 1.e-10;
|
wolffd@0
|
37
|
wolffd@0
|
38 % Maximal proportional step to take: don't want to make this too big
|
wolffd@0
|
39 % as then spend a lot of time finding the minimum inside the bracket
|
wolffd@0
|
40 max_step = 10.0;
|
wolffd@0
|
41
|
wolffd@0
|
42 fb = feval(f, b, varargin{:});
|
wolffd@0
|
43 num_evals = num_evals + 1;
|
wolffd@0
|
44
|
wolffd@0
|
45 % Assume that we know going from a to b is downhill initially
|
wolffd@0
|
46 % (usually because gradf(a) < 0).
|
wolffd@0
|
47 if (fb > fa)
|
wolffd@0
|
48 % Minimum must lie between a and b: do golden section until we find point
|
wolffd@0
|
49 % low enough to be middle of bracket
|
wolffd@0
|
50 c = b;
|
wolffd@0
|
51 b = a + (c-a)/phi;
|
wolffd@0
|
52 fb = feval(f, b, varargin{:});
|
wolffd@0
|
53 num_evals = num_evals + 1;
|
wolffd@0
|
54 while (fb > fa)
|
wolffd@0
|
55 c = b;
|
wolffd@0
|
56 b = a + (c-a)/phi;
|
wolffd@0
|
57 fb = feval(f, b, varargin{:});
|
wolffd@0
|
58 num_evals = num_evals + 1;
|
wolffd@0
|
59 end
|
wolffd@0
|
60 else
|
wolffd@0
|
61 % There is a valid bracket upper bound greater than b
|
wolffd@0
|
62 c = b + phi*(b-a);
|
wolffd@0
|
63 fc = feval(f, c, varargin{:});
|
wolffd@0
|
64 num_evals = num_evals + 1;
|
wolffd@0
|
65 bracket_found = 0;
|
wolffd@0
|
66
|
wolffd@0
|
67 while (fb > fc)
|
wolffd@0
|
68 % Do a quadratic interpolation (i.e. to minimum of quadratic)
|
wolffd@0
|
69 r = (b-a).*(fb-fc);
|
wolffd@0
|
70 q = (b-c).*(fb-fa);
|
wolffd@0
|
71 u = b - ((b-c)*q - (b-a)*r)/(2.0*(sign(q-r)*max([abs(q-r), TINY])));
|
wolffd@0
|
72 ulimit = b + max_step*(c-b);
|
wolffd@0
|
73
|
wolffd@0
|
74 if ((b-u)'*(u-c) > 0.0)
|
wolffd@0
|
75 % Interpolant lies between b and c
|
wolffd@0
|
76 fu = feval(f, u, varargin{:});
|
wolffd@0
|
77 num_evals = num_evals + 1;
|
wolffd@0
|
78 if (fu < fc)
|
wolffd@0
|
79 % Have a minimum between b and c
|
wolffd@0
|
80 br_min = b;
|
wolffd@0
|
81 br_mid = u;
|
wolffd@0
|
82 br_max = c;
|
wolffd@0
|
83 return;
|
wolffd@0
|
84 elseif (fu > fb)
|
wolffd@0
|
85 % Have a minimum between a and u
|
wolffd@0
|
86 br_min = a;
|
wolffd@0
|
87 br_mid = c;
|
wolffd@0
|
88 br_max = u;
|
wolffd@0
|
89 return;
|
wolffd@0
|
90 end
|
wolffd@0
|
91 % Quadratic interpolation didn't give a bracket, so take a golden step
|
wolffd@0
|
92 u = c + phi*(c-b);
|
wolffd@0
|
93 elseif ((c-u)'*(u-ulimit) > 0.0)
|
wolffd@0
|
94 % Interpolant lies between c and limit
|
wolffd@0
|
95 fu = feval(f, u, varargin{:});
|
wolffd@0
|
96 num_evals = num_evals + 1;
|
wolffd@0
|
97 if (fu < fc)
|
wolffd@0
|
98 % Move bracket along, and then take a golden section step
|
wolffd@0
|
99 b = c;
|
wolffd@0
|
100 c = u;
|
wolffd@0
|
101 u = c + phi*(c-b);
|
wolffd@0
|
102 else
|
wolffd@0
|
103 bracket_found = 1;
|
wolffd@0
|
104 end
|
wolffd@0
|
105 elseif ((u-ulimit)'*(ulimit-c) >= 0.0)
|
wolffd@0
|
106 % Limit parabolic u to maximum value
|
wolffd@0
|
107 u = ulimit;
|
wolffd@0
|
108 else
|
wolffd@0
|
109 % Reject parabolic u and use golden section step
|
wolffd@0
|
110 u = c + phi*(c-b);
|
wolffd@0
|
111 end
|
wolffd@0
|
112 if ~bracket_found
|
wolffd@0
|
113 fu = feval(f, u, varargin{:});
|
wolffd@0
|
114 num_evals = num_evals + 1;
|
wolffd@0
|
115 end
|
wolffd@0
|
116 a = b; b = c; c = u;
|
wolffd@0
|
117 fa = fb; fb = fc; fc = fu;
|
wolffd@0
|
118 end % while loop
|
wolffd@0
|
119 end % bracket found
|
wolffd@0
|
120 br_mid = b;
|
wolffd@0
|
121 if (a < c)
|
wolffd@0
|
122 br_min = a;
|
wolffd@0
|
123 br_max = c;
|
wolffd@0
|
124 else
|
wolffd@0
|
125 br_min = c;
|
wolffd@0
|
126 br_max = a;
|
wolffd@0
|
127 end
|