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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
function [x, options] = linemin(f, pt, dir, fpt, options, ...
varargin)
%LINEMIN One dimensional minimization.
%
% Description
% [X, OPTIONS] = LINEMIN(F, PT, DIR, FPT, OPTIONS) uses Brent's
% algorithm to find the minimum of the function F(X) along the line DIR
% through the point PT. The function value at the starting point is
% FPT. The point at which F has a local minimum is returned as X. The
% function value at that point is returned in OPTIONS(8).
%
% LINEMIN(F, PT, DIR, FPT, OPTIONS, P1, P2, ...) allows additional
% arguments to be passed to F().
%
% The optional parameters have the following interpretations.
%
% OPTIONS(1) is set to 1 to display error values.
%
% OPTIONS(2) is a measure of the absolute precision required for the
% value of X at the solution.
%
% OPTIONS(3) is a measure of the precision required of the objective
% function at the solution. Both this and the previous condition must
% be satisfied for termination.
%
% OPTIONS(14) is the maximum number of iterations; default 100.
%
% See also
% CONJGRAD, MINBRACK, QUASINEW
%
% Copyright (c) Ian T Nabney (1996-2001)
% Set up the options.
if(options(14))
niters = options(14);
else
niters = 100;
end
options(10) = 0; % Initialise count of function evaluations
display = options(1);
% Check function string
f = fcnchk(f, length(varargin));
% Value of golden section (1 + sqrt(5))/2.0
phi = 1.6180339887499;
cphi = 1 - 1/phi;
TOL = sqrt(eps); % Maximal fractional precision
TINY = 1.0e-10; % Can't use fractional precision when minimum is at 0
% Bracket the minimum
[br_min, br_mid, br_max, num_evals] = feval('minbrack', 'linef', ...
0.0, 1.0, fpt, f, pt, dir, varargin{:});
options(10) = options(10) + num_evals; % Increment number of fn. evals
% No gradient evals in minbrack
% Use Brent's algorithm to find minimum
% Initialise the points and function values
w = br_mid; % Where second from minimum is
v = br_mid; % Previous value of w
x = v; % Where current minimum is
e = 0.0; % Distance moved on step before last
fx = feval('linef', x, f, pt, dir, varargin{:});
options(10) = options(10) + 1;
fv = fx; fw = fx;
for n = 1:niters
xm = 0.5.*(br_min+br_max); % Middle of bracket
% Make sure that tolerance is big enough
tol1 = TOL * (max(abs(x))) + TINY;
% Decide termination on absolute precision required by options(2)
if (max(abs(x - xm)) <= options(2) & br_max-br_min < 4*options(2))
options(8) = fx;
return;
end
% Check if step before last was big enough to try a parabolic step.
% Note that this will fail on first iteration, which must be a golden
% section step.
if (max(abs(e)) > tol1)
% Construct a trial parabolic fit through x, v and w
r = (fx - fv) .* (x - w);
q = (fx - fw) .* (x - v);
p = (x - v).*q - (x - w).*r;
q = 2.0 .* (q - r);
if (q > 0.0) p = -p; end
q = abs(q);
% Test if the parabolic fit is OK
if (abs(p) >= abs(0.5*q*e) | p <= q*(br_min-x) | p >= q*(br_max-x))
% No it isn't, so take a golden section step
if (x >= xm)
e = br_min-x;
else
e = br_max-x;
end
d = cphi*e;
else
% Yes it is, so take the parabolic step
e = d;
d = p/q;
u = x+d;
if (u-br_min < 2*tol1 | br_max-u < 2*tol1)
d = sign(xm-x)*tol1;
end
end
else
% Step before last not big enough, so take a golden section step
if (x >= xm)
e = br_min - x;
else
e = br_max - x;
end
d = cphi*e;
end
% Make sure that step is big enough
if (abs(d) >= tol1)
u = x+d;
else
u = x + sign(d)*tol1;
end
% Evaluate function at u
fu = feval('linef', u, f, pt, dir, varargin{:});
options(10) = options(10) + 1;
% Reorganise bracket
if (fu <= fx)
if (u >= x)
br_min = x;
else
br_max = x;
end
v = w; w = x; x = u;
fv = fw; fw = fx; fx = fu;
else
if (u < x)
br_min = u;
else
br_max = u;
end
if (fu <= fw | w == x)
v = w; w = u;
fv = fw; fw = fu;
elseif (fu <= fv | v == x | v == w)
v = u;
fv = fu;
end
end
if (display == 1)
fprintf(1, 'Cycle %4d Error %11.6f\n', n, fx);
end
end
options(8) = fx;
|
