wolffd@0
|
1 function [x, options] = linemin(f, pt, dir, fpt, options, ...
|
wolffd@0
|
2 varargin)
|
wolffd@0
|
3 %LINEMIN One dimensional minimization.
|
wolffd@0
|
4 %
|
wolffd@0
|
5 % Description
|
wolffd@0
|
6 % [X, OPTIONS] = LINEMIN(F, PT, DIR, FPT, OPTIONS) uses Brent's
|
wolffd@0
|
7 % algorithm to find the minimum of the function F(X) along the line DIR
|
wolffd@0
|
8 % through the point PT. The function value at the starting point is
|
wolffd@0
|
9 % FPT. The point at which F has a local minimum is returned as X. The
|
wolffd@0
|
10 % function value at that point is returned in OPTIONS(8).
|
wolffd@0
|
11 %
|
wolffd@0
|
12 % LINEMIN(F, PT, DIR, FPT, OPTIONS, P1, P2, ...) allows additional
|
wolffd@0
|
13 % arguments to be passed to F().
|
wolffd@0
|
14 %
|
wolffd@0
|
15 % The optional parameters have the following interpretations.
|
wolffd@0
|
16 %
|
wolffd@0
|
17 % OPTIONS(1) is set to 1 to display error values.
|
wolffd@0
|
18 %
|
wolffd@0
|
19 % OPTIONS(2) is a measure of the absolute precision required for the
|
wolffd@0
|
20 % value of X at the solution.
|
wolffd@0
|
21 %
|
wolffd@0
|
22 % OPTIONS(3) is a measure of the precision required of the objective
|
wolffd@0
|
23 % function at the solution. Both this and the previous condition must
|
wolffd@0
|
24 % be satisfied for termination.
|
wolffd@0
|
25 %
|
wolffd@0
|
26 % OPTIONS(14) is the maximum number of iterations; default 100.
|
wolffd@0
|
27 %
|
wolffd@0
|
28 % See also
|
wolffd@0
|
29 % CONJGRAD, MINBRACK, QUASINEW
|
wolffd@0
|
30 %
|
wolffd@0
|
31
|
wolffd@0
|
32 % Copyright (c) Ian T Nabney (1996-2001)
|
wolffd@0
|
33
|
wolffd@0
|
34 % Set up the options.
|
wolffd@0
|
35 if(options(14))
|
wolffd@0
|
36 niters = options(14);
|
wolffd@0
|
37 else
|
wolffd@0
|
38 niters = 100;
|
wolffd@0
|
39 end
|
wolffd@0
|
40 options(10) = 0; % Initialise count of function evaluations
|
wolffd@0
|
41
|
wolffd@0
|
42 display = options(1);
|
wolffd@0
|
43
|
wolffd@0
|
44 % Check function string
|
wolffd@0
|
45 f = fcnchk(f, length(varargin));
|
wolffd@0
|
46
|
wolffd@0
|
47 % Value of golden section (1 + sqrt(5))/2.0
|
wolffd@0
|
48 phi = 1.6180339887499;
|
wolffd@0
|
49 cphi = 1 - 1/phi;
|
wolffd@0
|
50 TOL = sqrt(eps); % Maximal fractional precision
|
wolffd@0
|
51 TINY = 1.0e-10; % Can't use fractional precision when minimum is at 0
|
wolffd@0
|
52
|
wolffd@0
|
53 % Bracket the minimum
|
wolffd@0
|
54 [br_min, br_mid, br_max, num_evals] = feval('minbrack', 'linef', ...
|
wolffd@0
|
55 0.0, 1.0, fpt, f, pt, dir, varargin{:});
|
wolffd@0
|
56 options(10) = options(10) + num_evals; % Increment number of fn. evals
|
wolffd@0
|
57 % No gradient evals in minbrack
|
wolffd@0
|
58
|
wolffd@0
|
59 % Use Brent's algorithm to find minimum
|
wolffd@0
|
60 % Initialise the points and function values
|
wolffd@0
|
61 w = br_mid; % Where second from minimum is
|
wolffd@0
|
62 v = br_mid; % Previous value of w
|
wolffd@0
|
63 x = v; % Where current minimum is
|
wolffd@0
|
64 e = 0.0; % Distance moved on step before last
|
wolffd@0
|
65 fx = feval('linef', x, f, pt, dir, varargin{:});
|
wolffd@0
|
66 options(10) = options(10) + 1;
|
wolffd@0
|
67 fv = fx; fw = fx;
|
wolffd@0
|
68
|
wolffd@0
|
69 for n = 1:niters
|
wolffd@0
|
70 xm = 0.5.*(br_min+br_max); % Middle of bracket
|
wolffd@0
|
71 % Make sure that tolerance is big enough
|
wolffd@0
|
72 tol1 = TOL * (max(abs(x))) + TINY;
|
wolffd@0
|
73 % Decide termination on absolute precision required by options(2)
|
wolffd@0
|
74 if (max(abs(x - xm)) <= options(2) & br_max-br_min < 4*options(2))
|
wolffd@0
|
75 options(8) = fx;
|
wolffd@0
|
76 return;
|
wolffd@0
|
77 end
|
wolffd@0
|
78 % Check if step before last was big enough to try a parabolic step.
|
wolffd@0
|
79 % Note that this will fail on first iteration, which must be a golden
|
wolffd@0
|
80 % section step.
|
wolffd@0
|
81 if (max(abs(e)) > tol1)
|
wolffd@0
|
82 % Construct a trial parabolic fit through x, v and w
|
wolffd@0
|
83 r = (fx - fv) .* (x - w);
|
wolffd@0
|
84 q = (fx - fw) .* (x - v);
|
wolffd@0
|
85 p = (x - v).*q - (x - w).*r;
|
wolffd@0
|
86 q = 2.0 .* (q - r);
|
wolffd@0
|
87 if (q > 0.0) p = -p; end
|
wolffd@0
|
88 q = abs(q);
|
wolffd@0
|
89 % Test if the parabolic fit is OK
|
wolffd@0
|
90 if (abs(p) >= abs(0.5*q*e) | p <= q*(br_min-x) | p >= q*(br_max-x))
|
wolffd@0
|
91 % No it isn't, so take a golden section step
|
wolffd@0
|
92 if (x >= xm)
|
wolffd@0
|
93 e = br_min-x;
|
wolffd@0
|
94 else
|
wolffd@0
|
95 e = br_max-x;
|
wolffd@0
|
96 end
|
wolffd@0
|
97 d = cphi*e;
|
wolffd@0
|
98 else
|
wolffd@0
|
99 % Yes it is, so take the parabolic step
|
wolffd@0
|
100 e = d;
|
wolffd@0
|
101 d = p/q;
|
wolffd@0
|
102 u = x+d;
|
wolffd@0
|
103 if (u-br_min < 2*tol1 | br_max-u < 2*tol1)
|
wolffd@0
|
104 d = sign(xm-x)*tol1;
|
wolffd@0
|
105 end
|
wolffd@0
|
106 end
|
wolffd@0
|
107 else
|
wolffd@0
|
108 % Step before last not big enough, so take a golden section step
|
wolffd@0
|
109 if (x >= xm)
|
wolffd@0
|
110 e = br_min - x;
|
wolffd@0
|
111 else
|
wolffd@0
|
112 e = br_max - x;
|
wolffd@0
|
113 end
|
wolffd@0
|
114 d = cphi*e;
|
wolffd@0
|
115 end
|
wolffd@0
|
116 % Make sure that step is big enough
|
wolffd@0
|
117 if (abs(d) >= tol1)
|
wolffd@0
|
118 u = x+d;
|
wolffd@0
|
119 else
|
wolffd@0
|
120 u = x + sign(d)*tol1;
|
wolffd@0
|
121 end
|
wolffd@0
|
122 % Evaluate function at u
|
wolffd@0
|
123 fu = feval('linef', u, f, pt, dir, varargin{:});
|
wolffd@0
|
124 options(10) = options(10) + 1;
|
wolffd@0
|
125 % Reorganise bracket
|
wolffd@0
|
126 if (fu <= fx)
|
wolffd@0
|
127 if (u >= x)
|
wolffd@0
|
128 br_min = x;
|
wolffd@0
|
129 else
|
wolffd@0
|
130 br_max = x;
|
wolffd@0
|
131 end
|
wolffd@0
|
132 v = w; w = x; x = u;
|
wolffd@0
|
133 fv = fw; fw = fx; fx = fu;
|
wolffd@0
|
134 else
|
wolffd@0
|
135 if (u < x)
|
wolffd@0
|
136 br_min = u;
|
wolffd@0
|
137 else
|
wolffd@0
|
138 br_max = u;
|
wolffd@0
|
139 end
|
wolffd@0
|
140 if (fu <= fw | w == x)
|
wolffd@0
|
141 v = w; w = u;
|
wolffd@0
|
142 fv = fw; fw = fu;
|
wolffd@0
|
143 elseif (fu <= fv | v == x | v == w)
|
wolffd@0
|
144 v = u;
|
wolffd@0
|
145 fv = fu;
|
wolffd@0
|
146 end
|
wolffd@0
|
147 end
|
wolffd@0
|
148 if (display == 1)
|
wolffd@0
|
149 fprintf(1, 'Cycle %4d Error %11.6f\n', n, fx);
|
wolffd@0
|
150 end
|
wolffd@0
|
151 end
|
wolffd@0
|
152 options(8) = fx;
|