-
Notifications
You must be signed in to change notification settings - Fork 26
/
anneal.m
210 lines (196 loc) · 7.76 KB
/
anneal.m
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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
function [minimum,fval] = anneal(loss, parent, options)
% ANNEAL Minimizes a function with the method of simulated annealing
% (Kirkpatrick et al., 1983)
%
% ANNEAL takes three input parameters, in this order:
%
% LOSS is a function handle (anonymous function or inline) with a loss
% function, which may be of any type, and needn't be continuous. It does,
% however, need to return a single value.
%
% PARENT is a vector with initial guess parameters. You must input an
% initial guess.
%
% OPTIONS is a structure with settings for the simulated annealing. If no
% OPTIONS structure is provided, ANNEAL uses a default structure. OPTIONS
% can contain any or all of the following fields (missing fields are
% filled with default values):
%
% Verbosity: Controls output to the screen.
% 0 suppresses all output
% 1 gives final report only [default]
% 2 gives temperature changes and final report
% Generator: Generates a new solution from an old one.
% Any function handle that takes a solution as input and
% gives a valid solution (i.e. some point in the solution
% space) as output.
% The default function generates a row vector which
% slightly differs from the input vector in one element:
% @(x) (x+(randperm(length(x))==length(x))*randn/100)
% Other examples of possible solution generators:
% @(x) (rand(3,1)): Picks a random point in the unit cube
% @(x) (ceil([9 5].*rand(2,1))): Picks a point in a 9-by-5
% discrete grid
% Note that if you use the default generator, ANNEAL only
% works on row vectors. For loss functions that operate on
% column vectors, use this generator instead of the
% default:
% @(x) (x(:)'+(randperm(length(x))==length(x))*randn/100)'
% InitTemp: The initial temperature, can be any positive number.
% Default is 1.
% StopTemp: Temperature at which to stop, can be any positive number
% smaller than InitTemp.
% Default is 1e-8.
% StopVal: Value at which to stop immediately, can be any output of
% LOSS that is sufficiently low for you.
% Default is -Inf.
% CoolSched: Generates a new temperature from the previous one.
% Any function handle that takes a scalar as input and
% returns a smaller but positive scalar as output.
% Default is @(T) (.8*T)
% MaxConsRej: Maximum number of consecutive rejections, can be any
% positive number.
% Default is 1000.
% MaxTries: Maximum number of tries within one temperature, can be
% any positive number.
% Default is 300.
% MaxSuccess: Maximum number of successes within one temperature, can
% be any positive number.
% Default is 20.
%
%
% Usage:
% [MINIMUM,FVAL] = ANNEAL(LOSS,NEWSOL,[OPTIONS]);
% MINIMUM is the solution which generated the smallest encountered
% value when input into LOSS.
% FVAL is the value of the LOSS function evaluated at MINIMUM.
% OPTIONS = ANNEAL();
% OPTIONS is the default options structure.
%
%
% Example:
% The so-called "six-hump camelback" function has several local minima
% in the range -3<=x<=3 and -2<=y<=2. It has two global minima, namely
% f(-0.0898,0.7126) = f(0.0898,-0.7126) = -1.0316. We can define and
% minimise it as follows:
% camel = @(x,y)(4-2.1*x.^2+x.^4/3).*x.^2+x.*y+4*(y.^2-1).*y.^2;
% loss = @(p)camel(p(1),p(2));
% [x f] = ANNEAL(loss,[0 0])
% We get output:
% Initial temperature: 1
% Final temperature: 3.21388e-007
% Consecutive rejections: 1027
% Number of function calls: 6220
% Total final loss: -1.03163
% x =
% -0.0899 0.7127
% f =
% -1.0316
% Which reasonably approximates the analytical global minimum (note
% that due to randomness, your results will likely not be exactly the
% same).
% Reference:
% Kirkpatrick, S., Gelatt, C.D., & Vecchi, M.P. (1983). Optimization by
% Simulated Annealing. _Science, 220_, 671-680.
% joachim.vandekerckhove@psy.kuleuven.be
% $Revision: v5 $ $Date: 2006/04/26 12:54:04 $
def = struct(...
'CoolSched',@(T) (.8*T),...
'Generator',@(x) (x+(randperm(length(x))==length(x))*randn/100),...
'InitTemp',1,...
'MaxConsRej',1000,...
'MaxSuccess',20,...
'MaxTries',300,...
'StopTemp',1e-8,...
'StopVal',-Inf,...
'Verbosity',1);
% Check input
if ~nargin %user wants default options, give it and stop
minimum = def;
return
elseif nargin<2, %user gave only objective function, throw error
error('MATLAB:anneal:noParent','You need to input a first guess.');
elseif nargin<3, %user gave no options structure, use default
options=def;
else %user gave all input, check if options structure is complete
if ~isstruct(options)
error('MATLAB:anneal:badOptions',...
'Input argument ''options'' is not a structure')
end
fs = {'CoolSched','Generator','InitTemp','MaxConsRej',...
'MaxSuccess','MaxTries','StopTemp','StopVal','Verbosity'};
for nm=1:length(fs)
if ~isfield(options,fs{nm}), options.(fs{nm}) = def.(fs{nm}); end
end
end
% main settings
newsol = options.Generator; % neighborhood space function
Tinit = options.InitTemp; % initial temp
minT = options.StopTemp; % stopping temp
cool = options.CoolSched; % annealing schedule
minF = options.StopVal;
max_consec_rejections = options.MaxConsRej;
max_try = options.MaxTries;
max_success = options.MaxSuccess;
report = options.Verbosity;
k = 1; % boltzmann constant
% counters etc
itry = 0;
success = 0;
finished = 0;
consec = 0;
T = Tinit;
initenergy = loss(parent);
oldenergy = initenergy;
total = 0;
if report==2, fprintf(1,'\n T = %7.5f, loss = %10.5f\n',T,oldenergy); end
while ~finished;
itry = itry+1; % just an iteration counter
current = parent;
% % Stop / decrement T criteria
if itry >= max_try || success >= max_success;
if T < minT || consec >= max_consec_rejections;
finished = 1;
total = total + itry;
break;
else
T = cool(T); % decrease T according to cooling schedule
if report==2, % output
fprintf(1,' T = %7.5f, loss = %10.5f\n',T,oldenergy);
end
total = total + itry;
itry = 1;
success = 1;
end
end
newparam = newsol(current);
newenergy = loss(newparam);
if (newenergy < minF),
parent = newparam;
oldenergy = newenergy;
break
end
if (oldenergy-newenergy > 1e-6)
parent = newparam;
oldenergy = newenergy;
success = success+1;
consec = 0;
else
if (rand < exp( (oldenergy-newenergy)/(k*T) ));
parent = newparam;
oldenergy = newenergy;
success = success+1;
else
consec = consec+1;
end
end
end
minimum = parent;
fval = oldenergy;
if report;
fprintf(1, '\n Initial temperature: \t%g\n', Tinit);
fprintf(1, ' Final temperature: \t%g\n', T);
fprintf(1, ' Consecutive rejections: \t%i\n', consec);
fprintf(1, ' Number of function calls:\t%i\n', total);
fprintf(1, ' Total final loss: \t%g\n', fval);
end