- exponentially-decaying temperature annealing object
copyright 2009-2012 Blair Armstrong, Christine Watson, David Plaut
This file is part of SOS
SOS is free software: you can redistribute it and/or modify
it for academic and non-commercial purposes
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version. For commercial or for-profit
uses, please contact the authors (sos@cnbc.cmu.edu).
SOS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.0001 % - exponentially-decaying temperature annealing object 0002 % 0003 % copyright 2009-2012 Blair Armstrong, Christine Watson, David Plaut 0004 % 0005 % This file is part of SOS 0006 % 0007 % SOS is free software: you can redistribute it and/or modify 0008 % it for academic and non-commercial purposes 0009 % under the terms of the GNU General Public License as published by 0010 % the Free Software Foundation, either version 3 of the License, or 0011 % (at your option) any later version. For commercial or for-profit 0012 % uses, please contact the authors (sos@cnbc.cmu.edu). 0013 % 0014 % SOS is distributed in the hope that it will be useful, 0015 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0016 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0017 % GNU General Public License for more details. 0018 0019 % You should have received a copy of the GNU General Public License 0020 % along with SOS (see COPYING.txt). 0021 % If not, see <http://www.gnu.org/licenses/>. 0022 0023 0024 0025 classdef expAnneal < genericAnneal 0026 %% Provides support for exponentially decaying temperature annealing 0027 % 0028 % Provides support for an exponential anneal function governed by the 0029 % following equation: 0030 % 0031 % temp = N*e^(-lambda*curStep); 0032 % 0033 % The Schedule operates as follows: 0034 % {blockSize} iterations are run with temperature set to Inf (the -1 0035 % value of curStep, which reflects the current temperature step). 0036 % provides an initial sample of the types of cost values encountered 0037 % when making purely random switches. From this sample, the maximum 0038 % delta cost in the block is calculated by subtracting the biggest 0039 % cost from the smallest and multiplying this value by 10. The x10 0040 % multiple is applied because in the context of the logicstic 0041 % function, a temperature 10x greater than deltaCost provides a good 0042 % approximation of random behavior. 0043 % 0044 % Next, to fit the lambda term in the equation 0045 % the user must express how large of a decrease in this initial temp 0046 % is desired on the next step, as a probability. For slower annealing, the 0047 % decrease in temperature should be small, and values should be close 0048 % to 1to 0. For faster descent, values closer to 1 should be used. 0049 % 0050 % Once the values for the variable in the formula have been derived, 0051 % the curStep is set to '0' and the initial temperature is calculated 0052 % using the formula (in this case = N). Thus, the '0' step should 0053 % also basically yield random. 0054 % 0055 % Subsequent changes in temperature occur when the state of the 0056 % optimizer is said to have reached thermal equilibrium and a new 0057 % temperature step begins. This 0058 % approximated in the algorithm as when the previous two 'blocks' of 0059 % cost values are non-significantly different from one another 0060 % according to an independent samples t-test. The value of this 0061 % t-test can be adjusted to make it easier or harder to reject the 0062 % null hypothesis. In the present case, thermal equilibrium is 0063 % reached when we fail to reject the null hypothesis, and thus one 0064 % could also make a Type-II error and incorrectly fail to reject a 0065 % false null hypothesis. This can be avoided by using p-values 0066 % closer to 1. It doesn't take that much to cause a p-value to 0067 % exceed even the standard cutoff of 0.05 though if blockSizse 0068 % becomes large (e.g., 1000), so keep that in mind, otherwise thermal 0069 % equilibrium may never be reached 0070 % 0071 % It is worth noting, however, that this measurement of 0072 % thermal equilibrium is merely intended to serve as a quick and 0073 % readily-interpretable approximation thereof. A failure to reject 0074 % the null hypothesis does not stricly == confirmation of the null 0075 % hypothesis, and there are no explicit checks of the t-test 0076 % assumptions to have a strong basis for believing that the exact 0077 % p-value set is the true p-value for a given data set. Again 0078 % though, for present purposes this approximation has proven to 0079 % adequately satisfy its intended purpose. 0080 % 0081 % 0082 % PROPERTIES: 0083 % blockSize % size of the block of cost values used to determine initial max Delta cost and assess thermal equilibrium 0084 % pval % p-value for the test of equality of cost; p-values greater than this are indicative of thermal equilibrium. 0085 % N % N from exponential decay equation 0086 % lambda % lambda from exponental decay equation 0087 % curStep % current temperatuer step 0088 % pDecrease % proportion decrease in temperature on next step. 0089 % tempLog % log of all temperature values produced by the algorithm, and the iteration they were set at 0090 % prevBlock % history of cost values from previous block 0091 % curBlock % history of cost values from current block 0092 % startDeltaCost % history of delta cost values from the first block for equation calibration. 0093 % 0094 %METHODS: 0095 % expAnneal(varargin) - CONSTRUCTOR 0096 % temp = getTemp(obj) - Returns current temperature 0097 % anneal(obj,curIt,cost,deltaCost) - anneal current temperature using the exponentially decay annealing schedule. 0098 % 0099 %METHODS (Static) 0100 % numSteps(initDeltaCost, finalDeltaCost, pDecrease) % calculates the number of steps that would be necessary to go from initDeltaCost to finaldeltaCost 0101 % maxpDecrease(initDeltaCost, finalDeltaCost,nStep) % calculates the maxDecrease that would be necessary to go from initDeltaCost to finaldeltaCost using at least 3 drops in temperature. 0102 0103 %% PROPERTIES 0104 properties 0105 blockSize % size of the block of cost values used to determine initial max Delta cost and assess thermal equilibrium 0106 pval % p-value for the test of equality of cost; p-values greater than this are indicative of thermal equilibrium. 0107 N % N from exponential decay equation 0108 lambda % lambda from exponental decay equation 0109 curStep % current temperature step 0110 pDecrease % proportion decrease in temperature on next step. 0111 tempLog % log of all temperature values produced by the algorithm, and the iteration they were set at 0112 prevBlock % history of cost values from previous block 0113 curBlock % history of cost values from current block 0114 startDeltaCost % history of delta cost values from the first block for equation calibration. 0115 overrideInitDeltaCost % user-specified override for the calculated initDeltaCost value. 0116 end 0117 0118 0119 methods 0120 %% expAnneal CONSTRUCTOR 0121 function obj = expAnneal(varargin) 0122 % creates and initializes an expAnneal object 0123 % 0124 %PARAMETERS: 0125 %Optional: 0126 % 'blockSize'/integer - size of blocks to use when estimating the initial temperature and thermal equilibrium. Default: 1000 0127 % 'pVal'/numeric - p-value to use when testing for thermal equilibrium. Default = 0.05 0128 % 'pDecrease'/numeric - % decrease in previous temperature when temperature is decreased 0129 % 'schedule'/'exp' - name of annealing schedule 0130 0131 verbosePrint('Creating exponential decay temperature annealing function...',... 0132 'expAnneal_constructor_startObjCreation'); 0133 0134 p = expAnneal.parseConstructorArgs(varargin); 0135 0136 % potentially altered by user input 0137 obj.blockSize = p.Results.blockSize; % must be >=2 0138 obj.pval = p.Results.pval; 0139 obj.pDecrease = p.Results.pDecrease; 0140 0141 0142 if(any(strcmp(p.UsingDefaults,'overrideInitDeltaCost'))) 0143 obj.overrideInitDeltaCost = NaN; 0144 else 0145 obj.overrideInitDeltaCost = p.Results.overrideInitDeltaCost; 0146 end 0147 0148 0149 % variables for expDecay equation 0150 obj.curStep = -1; 0151 obj.N = NaN; 0152 obj.lambda = NaN; 0153 0154 obj.tempLog = []; 0155 obj.temp = Inf; 0156 obj.tempLog = [obj.temp 0]; 0157 0158 obj.curBlock = nan(obj.blockSize,1); 0159 obj.startDeltaCost = nan(obj.blockSize,1); 0160 end 0161 0162 %% temp = getTemp() METHOD 0163 function temp = getTemp(obj) 0164 %Returns current temperature 0165 temp = obj.temp; 0166 end 0167 0168 %% anneal(obj,sosIt,cost) METHOD 0169 function anneal(obj,curIt,cost,deltaCost) 0170 %anneal current temperature using the exponentially decay annealing schedule. 0171 % 0172 % PARAMETERS: 0173 % curIt - current iteration in sosObj 0174 % cost - cost value produced during optimization 0175 % deltaCost - deltaCost value during optimization 0176 0177 if(mod(curIt,obj.blockSize)) == 0 0178 curIndex = obj.blockSize; 0179 else 0180 curIndex = mod(curIt,obj.blockSize); 0181 end 0182 0183 obj.curBlock(curIndex) = cost; 0184 0185 0186 if(curIt <= obj.blockSize) 0187 obj.startDeltaCost(curIt) = deltaCost; 0188 end 0189 0190 % if < 1 stepSize of iterations have been run, temp = inf to 0191 % derive an initial estimation of typical cost values. After 0192 % this block, we can fit the decay equation 0193 if(mod(curIt,obj.blockSize) == 0 && (curIt)/obj.blockSize == 1 ) 0194 0195 % The safest thing to do is to always use the largest 0196 % possible difference (maxCost - minCost in the set) as the 0197 % largest possible delta cost. However, this makes it 0198 % harder to reliably replicate a given run because the 0199 % initial calibration can be influenced by those two data 0200 % points alone. So instead, use the 90th-10th percentile 0201 % as the basis for calibrating the annealing equation. 0202 0203 percentiles = [2.5 97.5]; 0204 scores = prctile(obj.startDeltaCost, percentiles); 0205 0206 DeltaCost95 = scores(2) - scores(1); 0207 0208 %apply override if it has been specified; 0209 0210 if isnan(obj.overrideInitDeltaCost) == 0 0211 DeltaCost95 = obj.overrideInitDeltaCost; 0212 0213 verbosePrint('Applying user-specified override to initDeltaCost', ... 0214 'expAnneal_anneal_calibration'); 0215 end 0216 0217 % fit the decay function. 0218 % for the logistic funciton, if temp is 10x greater than 0219 % the largest cost, even the minCost 0220 % would only have a 0.0525 chance of eliciting a swap. 0221 0222 0223 0224 obj.N = DeltaCost95; 0225 0226 step1Cost = (1-obj.pDecrease)*obj.N; 0227 obj.lambda = -1*log((step1Cost/obj.N))/1; 0228 0229 obj.curStep = 0; 0230 0231 obj.temp = obj.N*exp(-obj.lambda*obj.curStep); 0232 obj.tempLog = [obj.tempLog ; obj.temp curIt]; 0233 0234 obj.prevBlock = obj.curBlock; 0235 obj.curBlock = nan(obj.blockSize,1); 0236 0237 0238 verbosePrint([sprintf('%i) ', curIt), ... 0239 'Annealing Equation calibrated, changing temperature from ', ... 0240 num2str(obj.tempLog(length(obj.tempLog)-1),1), ... 0241 ' to ', num2str(obj.temp)], ... 0242 'expAnneal_anneal_calibration'); 0243 0244 % otherwise, if a new block has been completed, see if thermal 0245 % thermal equilibrium has been reached, i.e., cost is no longer 0246 % descending 0247 elseif (mod(curIt,obj.blockSize) == 0 && (curIt)/obj.blockSize > 1 ) 0248 0249 [H,p] = ttest2(obj.prevBlock,obj.curBlock,obj.pval,'both',... 0250 'unequal'); 0251 0252 verbosePrint([sprintf('%i) ', curIt), ... 0253 'p(thermEquil): ', num2str(p), ... 0254 ' prevBlock m = ', num2str(mean(obj.prevBlock)), ... 0255 ' (se = ', num2str(std(obj.prevBlock)/length(obj.prevBlock)), ... 0256 ') curBlock m = ', num2str(mean(obj.curBlock)), ... 0257 ' (se = ', num2str(std(obj.curBlock)/length(obj.curBlock)), ... 0258 ')'], ... 0259 'expAnneal_anneal_pthermalEquil'); 0260 0261 %we've reached our operational definition of thermal 0262 %equilibrium if H == 0 0263 if H == 0 0264 obj.curStep = obj.curStep + 1; 0265 obj.temp = obj.N*exp(-obj.lambda*obj.curStep); 0266 obj.tempLog = [obj.tempLog ; obj.temp curIt]; 0267 0268 verbosePrint([sprintf('%i) ', curIt), ...' 0269 'Thermal Equilibrium Reached - Dropping temperature from ', ... 0270 num2str(obj.tempLog(length(obj.tempLog)-1,1)), ... 0271 ' to ', num2str(obj.temp)], ... 0272 'expAnneal_anneal_dropTemp'); 0273 end 0274 0275 obj.prevBlock = obj.curBlock; 0276 obj.curBlock = nan(obj.blockSize,1); 0277 0278 end 0279 end % anneal 0280 end 0281 0282 0283 methods (Static) 0284 0285 %% p = expAnnealInputParser() 0286 function p = expAnnealInputParser() 0287 % creates an input parser for the constructor args 0288 % 0289 %PARAMETERS: 0290 % see constructor args 0291 0292 p = inputParser; 0293 0294 p.addParamValue('blockSize',1000, ... 0295 @(blockSize)validateattributes(blockSize, {'numeric'}, ... 0296 {'scalar', 'integer', 'positive', '>', 1})); 0297 p.addParamValue('pval',0.05, ... 0298 @(pval)validateattributes(pval, {'numeric'}, ... 0299 {'scalar', 'positive', '>', 0, '<',1})); 0300 p.addParamValue('pDecrease',0.5,... 0301 @(pDecrease)validateattributes(pDecrease, {'numeric'}, ... 0302 {'scalar', 'positive', '>', 0, '<',1})); 0303 p.addParamValue('schedule','exp', ... 0304 @(schedule)strcmp(schedule,'exp')); 0305 %override value should be overriden later to set it to Nan 0306 p.addParamValue('overrideInitDeltaCost',0,... 0307 @(overrideInitDeltaCost)validateattributes(overrideInitDeltaCost, {'numeric'}, ... 0308 {'scalar'})); 0309 0310 end 0311 0312 0313 %% numSteps(initDeltaCost, finalDeltaCost, pDecrease) STATIC FUNCTION 0314 function numSteps(initDeltaCost, finalDeltaCost, pDecrease) 0315 % calculates the number of steps that would be necessary to go 0316 % from initDeltaCost to finaldeltaCost, given a particular 0317 % pDecrease. 0318 % 0319 % The intended use of this function is to help determine the 0320 % value of pDecrease. To do so, a couple of steps are 0321 % invovled. First, run sos with an exponential anneal function 0322 % for expAnneal.blockSize iterations. During this first block, 0323 % temperature is set to infinitie, so cost values are changed 0324 % randomly. Once blockSize iteractions have been completed, 0325 % run <sosObj>.deltaCostPercentiles to see the deltaCost value 0326 % for the 97.5th - 2.5th percentile, which is the deltaCost 0327 % value used by the exponential anneal object during it's 0328 % calibration of the initial temperature value. The next step 0329 % is determining the finalDeltaCost value. To determine that 0330 % value, run SOS again, this time in greedy mode, and examine 0331 % the distribution of deltaCosts once the algorithm reaches a 0332 % frozen state. For the stochastic annealing to find a better 0333 % solution than what was found by the exponential annealer, it 0334 % must approach that frozen state slowly. For the frozen state 0335 % to have been reached, by definition no swaps occured and 0336 % therefore all delta costs should have an upper bound of zero. 0337 % For a very slow anneal, you could try to reach the point at 0338 % which a large number of swaps would have occured quite slowly, 0339 % which corresponds to the delta cost from a high decile like 90%). 0340 % This could take a while though. Alternatively, you could try 0341 % to reach a lower decile instead (e.g., 10%). A good 0342 % compromise might be to aim for the 50th percentile so that 0343 % you approach that point more slowly. So to sum up to this 0344 % point, input the initdeltaCost for the 97.5th to the 2.5th 0345 % percentile from the first block in an 0346 % exponential anneal optimization and the finalDeltaCost from 0347 % a perceintile of your choice as the from the percentiles 0348 % listed when a greedy anneal was run to a frozen state. 0349 % 0350 % Now, you have to decide what the pDecrease value will be. It 0351 % must be bounded between 0-1, and it should lead to at least 3 0352 % steps for the exponential anneal to really be effective at 0353 % gradually lowering temperature. of course, more steps than 3 0354 % will result in even slower and potentially better annealing. 0355 % (1 step at inf, 1 intermediate step, and 1 step at a 0356 % temperature low enough to cause freezing in a greedy anneal). 0357 % Try a few different values to see how they would correspond 0358 % to this constraint and the amount of time you're willing to 0359 % wait while the algorithm runs. Lowever values should lead to 0360 % more steps, whereas higher value will lead to fewer steps. 0361 % 0362 % 0363 % PARAMETERS 0364 % -initDeltaCost - the initial delta cost for the first 0365 % block of optimizatoin 0366 % -finalDeltaCost - the approximate delta cost value that you 0367 % wish to approach slowly to avoid premature freezing 0368 % -pDecrease - the reduction in the previous temperature to be 0369 % applied to generate the new temperature 0370 % 0371 % 0372 0373 0374 0375 validateattributes(initDeltaCost, {'numeric'}, ... 0376 {'scalar', 'positive', '>', 0}); 0377 validateattributes(finalDeltaCost, {'numeric'}, ... 0378 {'scalar', 'positive', '>', 0}); 0379 validateattributes(pDecrease, {'numeric'}, ... 0380 {'scalar', 'positive', '>', 0, '<',1}); 0381 0382 if(finalDeltaCost >= initDeltaCost) 0383 error('final delta cost must be < init delta cost. Consult the manual if this is not the case for your data.'); 0384 end 0385 0386 0387 0388 % calibrate the expAnneal equation 0389 N = initDeltaCost*1; %#ok<PROP> 0390 step1Cost = (1-pDecrease)*N; %#ok<PROP> 0391 lambda = -1*log((step1Cost/N))/1; %#ok<PROP> 0392 0393 %calculate initial temperature; 0394 initTemp = N*exp(-lambda*0); %#ok<PROP> 0395 0396 %calculate the ifnal temperature. 0397 finalTemp = (initTemp * (finalDeltaCost/(initDeltaCost*10)))*10; 0398 0399 %calculate the number of steps. 0400 nStep = -1*log(finalTemp/N)/lambda; %#ok<PROP> 0401 nStep = upper(nStep); 0402 0403 msg = sprintf('\n%g steps needed to reach a temp low enough to freeze given finalDeltaCost\n', ... 0404 nStep); 0405 0406 verbosePrint(msg, ... 0407 'expAnneal_numSteps_nStep'); 0408 0409 if(nStep < 3) 0410 verbosePrint(['Warning: pDecrease should be lowered until nStep >= 3 ',... 0411 char(10),' at least if exponentially decaying annealing ',... 0412 char(10),' is to have an appreciable effect on annealing'],... 0413 'expAnneal_numSteps_Warn'); 0414 end 0415 0416 end 0417 0418 %% maxpDecrease(initDeltaCost, finalDeltaCost) STATIC FUNCTION 0419 function maxpDecrease(initDeltaCost, finalDeltaCost,nStep) 0420 % calculates the maxDecrease that would be necessary to go 0421 % from initDeltaCost to finaldeltaCost using at least 3 drops 0422 % in temperature. 0423 % 0424 % The intended use of this function is to help determine the 0425 % value of pDecrease. To do so, a couple of steps are 0426 % invovled. First, run sos with an exponential anneal function 0427 % for expAnneal.blockSize iterations. During this first block, 0428 % temperature is set to infinitie, so cost values are changed 0429 % randomly. Once blockSize iteractions have been completed, 0430 % run <sosObj>.deltaCostPercentiles to see the deltaCost value 0431 % for the 97.5th - 2.5th percentile, which is the deltaCost 0432 % value used by the exponential anneal object during it's 0433 % calibration of the initial temperature value. The next step 0434 % is determining the finalDeltaCost value. To determine that 0435 % value, run SOS again, this time in greedy mode, and examine 0436 % the distribution of deltaCosts once the algorithm reaches a 0437 % frozen state. For the stochastic annealing to find a better 0438 % solution than what was found by the exponential annealer, it 0439 % must approach that frozen state slowly. For the frozen state 0440 % to have been reached, by definition no swaps occured and 0441 % therefore all delta costs should have an upper bound of zero. 0442 % For a very slow anneal, you could try to reach the point at 0443 % which a large number of swaps would have occured quite slowly, 0444 % which corresponds to the delta cost from a high decile like 90%). 0445 % This could take a while though. Alternatively, you could try 0446 % to reach a lower decile instead (e.g., 10%). A good 0447 % compromise might be to aim for the 50th percentile so that 0448 % you approach that point more slowly. So to sum up to this 0449 % point, input the initdeltaCost for the 97.5th to the 2.5th 0450 % percentile from the first block in an 0451 % exponential anneal optimization and the finalDeltaCost from 0452 % a perceintile of your choice as the from the percentiles 0453 % listed when a greedy anneal was run to a frozen state. 0454 % 0455 % From this point, the algorithm will determine what the value 0456 % of pDecrease would be if you were to have 3 steps in 0457 % temperature during the annealing. This is probably the 0458 % minimum number that youèd wnat oth ave to be useful 0459 % (otherwise you simply go from a first step of near-random 0460 % behavior to a next step of basically stochastic behavior). 0461 % If desired, you can then experiment with how many steps you 0462 % 0463 % 0464 % PARAMETERS 0465 % -initDeltaCost - the initial delta cost for the first 0466 % block of optimizatoin 0467 % -finalDeltaCost - the approximate delta cost value that you 0468 % wish to approach slowly to avoid premature freezing 0469 % -pDecrease - the reduction in the previous temperature to be 0470 % applied to generate the new temperature 0471 % 0472 % 0473 0474 0475 0476 validateattributes(initDeltaCost, {'numeric'}, ... 0477 {'scalar', 'positive', '>', 0}); 0478 validateattributes(finalDeltaCost, {'numeric'}, ... 0479 {'scalar', 'positive', '>', 0}); 0480 validateattributes(nStep, {'numeric'}, ... 0481 {'integer', 'scalar', 'positive', '>', 0}); 0482 0483 if(finalDeltaCost >= initDeltaCost) 0484 error('final delta cost must be < init delta cost. See the manual for what to do if this is the case'); 0485 end 0486 0487 if(nStep < 3) 0488 error('nStep should be >=3 for exponentially-decaying annealing to have an appreciable effect'); 0489 end 0490 0491 % calibrate the expAnneal equation 0492 N = initDeltaCost*1; %#ok<PROP> 0493 0494 %calculate initial temperature; 0495 initTemp = N*exp(0); %#ok<PROP> 0496 0497 %calculate the ifnal temperature 0498 finalTemp = (initTemp * (finalDeltaCost/(initDeltaCost*10)))*10; 0499 0500 0501 lambda= -1*log(finalTemp/N)/(nStep); %#ok<PROP> 0502 0503 step1Cost = N*exp(-lambda*1); %#ok<PROP> 0504 0505 pDecrease = -1*(step1Cost/N -1); %#ok<PROP> 0506 0507 0508 msg = sprintf('\n%f is max pDecrease to ensure %d steps during exp anneal', ... 0509 pDecrease,nStep); %#ok<PROP> 0510 0511 verbosePrint(msg, ... 0512 'expAnneal_maxpDecrease_maxpDecrease'); 0513 0514 0515 end 0516 0517 0518 end 0519 0520 methods (Static, Access = private) 0521 0522 %% parseConstructorArgs(varargin) 0523 function p = parseConstructorArgs(varargin) 0524 % parsers the constructor arguments 0525 % 0526 %PARAMETERS: 0527 % see constructor args 0528 0529 varargin = varargin{1}; 0530 p = expAnneal.expAnnealInputParser(); 0531 p.parse(varargin{:}); 0532 end 0533 end 0534 end 0535