Simulated Annealing

Home * Programming * Algorithms * Simulated Annealing

Glowing train wheels 1 Glowing train wheels ¹

train wheel production, Bochumer Verein, Bochum, Germany, ExtraSchicht 2010, The Industrial Heritage Trail, image by Rainer Halama, June 19, 2010, CC BY-SA 3.0, Wikimedia Commons, Glühen from Wikipedia.de (German)↩︎

Simulated Annealing, (SA)
a Monte Carlo based algorithm for combinatorial optimization problems inspired by statistical mechanics in thermodynamics with the statistical ensemble of the probability distribution over all possible states of a system described by a Markov chain, where its stationary distribution converts to an optimal distribution during a cooling process after reaching the equilibrium. Thus, the annealing algorithm simulates a nonstationary finite state Markov chain whose state space is the domain of the cost function called energy to be minimized ¹.

Contents

History

The annealing algorithm is an adaptation of the Metropolis–Hastings algorithm to generate sample states of a thermodynamic system, invented by Marshall Rosenbluth and published by Nicholas Metropolis et al. in 1953 ² ³ , later generalized by W. Keith Hastings at University of Toronto ⁴. According to Roy Glauber and Emilio Segrè, the original algorithm was invented by Enrico Fermi and reinvented by Stanislaw Ulam ⁵.

SA was independently described by Scott Kirkpatrick, C. Daniel Gelatt and Mario P. Vecchi in 1983 ⁶, at that time affiliated with IBM Thomas J. Watson Research Center, Yorktown Heights, and by Vlado Černý from Comenius University, Bratislava in 1985 ⁷.

Quotes

In the 2003 conference proceedings Celebrating the 50th Anniversary of the Metropolis Algorithm ⁸ ⁹, Marshall Rosenbluth describes the algorithm in the following beautifully concise and clear manner ¹⁰:

A simple way to do this [sampling configurations with the Boltzmann weight], as emerged after discussions with Teller, would be to make a trial move: if it decreased the energy of the system, allow it; if it increased the energy, allow it with probability exp(−ΔE/kT) as determined by a comparison with a random number. Each step, after an initial annealing period, is counted as a member of the ensemble, and the appropriate ensemble average of any quantity determined.

Applications

SA has multiple applications in discrete NP-hard optimization problems such as the Travelling salesman problem, in machine learning, in training of neural networks, and in the domain of computer games and computer chess in automated tuning as elaborated by Peter Mysliwietz in his Ph.D. thesis ¹¹ to optimize the evaluation weight vector in Zugzwang. In its variant of temporal difference learning to adjust pattern weights in Morph, Robert Levinson at al. used simulated annealing as metaheuristic to set its own learning rate for each pattern, the more frequently a pattern is updated, the slower becomes its learning rate ¹² ¹³ ¹⁴.

Algorithm

Description

The control flow of the algorithm is determined by two nested loops, the outer loop over decreasing temperature simulates the cooling, and an inner loop times n Monte Carlo iterations. Each time a randomly picked neighbor state inside the inner loop provides a better energy or fitness than the current state, the neighbor becomes the new current and even new optimum if fitter than fittest so far. Otherwise, if the neighbor fitness does not exceed current, it might still become current depending on the positive fitness or energy difference ΔE, and absolute temperature T, with a probability p according to the Boltzmann factor:

where k the Boltzmann constant, and e base of the exponential function whose negative exponent ensures the [0, 1] probability interval. Accepting worse solutions is a primary feature of SA, and important to stop greedy exploitation a local optimum but to explore other areas - higher temperatures favor exploration, while decreasing temperatures make the algorithm to behave greedier in favoring exploitation of the hopefully global optimum.

Animation

Simulated annealing - searching for a maximum. ¹⁵
With the high temprature, the numerous local maxima are left quickly through the strong noise movement -
but the global maximum is reliably found because of cooling temperature is no longer sufficient to leave it.

Pseudo Code

The C like pseudo code is based on Peter Mysliwietz' description as given in his Ph.D. thesis ¹⁶. Several neighbor functions used to modify the weight vector were tried, where one randomly chosen element changed randomly performed well. The fitness function inside the inner loop is of course the most time consuming part. For Zugzwang, Mysliwietz used a database of 500 test-positions with a search depth of one ply, which took about three minutes on a T 800 Transputer per iteration - the higher the hit rate of found expert moves, the fitter. The whole optimization used a tHight to tLow ratio of 100, a reduction factor r of 0.95, and n=40 inner iterations.

/**
 * simulatedAnnealing
 * @author Peter Mysliwietz, slightly modified
 * @param tHigh is the start temperature
 * @param  tLow is the minimal end temperature
 * @param     r is the temperature reduction factor < 1.0   
 * @param     n number of iterations for each temperature     
 * @return best weight vector
 */
vector simulatedAnnealing(double tHigh, double tLow, double r, int n) {
   vector currentWeights = randomWeights();
   vector bestWeights = currentWeights;
   double fittest = fitness(currentWeights);

   for (double t = tHigh; t > tLow; t *= r) {
      for (int i = 0; i < n; ++i) {
         vector neighborWeights = neighbor(currentWeights);
         if ( fitness(neighborWeights ) > fitness(currentWeights) ) {
            currentWeights = neighborWeights;         
            if ( fitness(neighborWeights ) > fittest ) {
               fittest = fitness(neighborWeights);
               bestWeights = neighborWeights;
            }
         } else if (accept( fitness(currentWeights) - fitness(neighborWeights ), t) ) {
            currentWeights = neighborWeights;
         }
      } /* for i */
   } /* for t */
   return bestWeights;
}

/**
 * accept
 * @param d is the energy difference >= 0
 * @param t is the current temperature
 * @return true with probability of Boltzmann factor e^(-d/kt) 
 */
bool accept(double d, double t ) {
   const double k = 1.38064852e−23; /* joule / kelvin */
   double p = exp(-d / (k*t) ); 
   double r = rand() / (RAND_MAX + 1.0);
   return r < p;
}

¹⁷ ¹⁸