Simulated Annealing

Sometimes, a prob­lem does­n’t have a known, ef­fi­cient, de­ter­min­is­tic al­go­rithm. This hic­cup is­n’t nec­es­sar­ily a big deal be­cause of­ten we only care about get­ting an an­swer that is close enough. One al­go­rithm that can give ap­prox­i­mate val­ues for hard to solve prob­lems is sim­u­lated an­neal­ing. Before jump­ing into the R code, let’s look at a mo­ti­vat­ing ex­am­ple.

Motivation

The trav­el­ing sales­man prob­lem does­n’t have a poly­no­mial-time de­ter­min­is­tic al­go­rithm and thus is a good can­di­date for us to ex­plore. The prob­lem name it­self nearly gives away the de­tails. Imagine a map with some fixed num­ber of ran­domly placed cities. We now want to hire a sales­per­son—let’s name her Sally—who needs to visit every city, with­out re­vis­it­ing any city, to sell some prod­uct.

When the num­ber of cities is small, we can use brute force tech­niques to try every pos­si­ble per­mu­ta­tion of city or­der­ings. Of course, for a rel­a­tively large amount of cities, try­ing to find the route this way will be in­tractable. We aren’t com­pletely doomed though. Do we care if Sally takes the most op­ti­mal route? For this sit­u­a­tion, not re­ally. While we want Sally to take a fast route, it does­n’t have to be the fastest route. Luckily, it turns out that sim­u­lated an­neal­ing can tractably gen­er­ate fast routes.

How does it work?

The trav­el­ing sales­man prob­lem is a lit­tle com­pli­cated, so let’s set up a toy prob­lem in­stead. For this demon­stra­tion, we are go­ing to find the global min­i­mum of the sine func­tion for the do­main [0, 2π]. For the math­e­mat­i­cally in­clined, you could use tech­niques from cal­cu­lus class to find that the an­swer is -1. Thus, if our al­go­rith­m’s an­swer is not in the same ball­park as -1, then we know some­thing went wrong.

State Spaces

Simulated an­neal­ing searches through a state space it­er­a­tively. An it­er­a­tive al­go­rithm be­gins with an ini­tial guess. This guess is then po­ten­tially up­graded every it­er­a­tion of the al­go­rithm to pro­duce a bet­ter ap­prox­i­ma­tion. For our demon­stra­tion, we’ll end the al­go­rithm af­ter a fixed num­ber of it­er­a­tions.

Taking care of the ini­tial guess is easy–just choose a ran­dom valid state. From the ini­tial guess, we’ll need to find it’s neigh­bor­ing states and check to see if they’re bet­ter. Disappointingly, there is no sin­gle method to de­rive neigh­bor­ing states. We’ll have to come up with a method spe­cific to our prob­lem.

Here, I want to keep things sim­ple. The im­ple­men­ta­tion be­low adds a ran­dom num­ber from the range [-0.1, 0.1] to the cur­rent state to gen­er­ate the suc­ces­sor. I also wrap the new state value around if it ex­ceeds the range [0, 2π].

gen_successor <- function(x) {
  TWO_PI <- 2 * pi
  successor <- x + runif(1, -.1, .1)
  
  if(successor < 0) {
    return(TWO_PI + successor)
  }
  else if(successor > TWO_PI) {
    return(0 + (successor - TWO_PI))
  }
    
  successor
}

Of course, you can ex­per­i­ment on your own. For ex­am­ple, try think­ing of new ways to gen­er­ate neigh­bor­ing states.

Temperature Schedule Go Brrr

This tem­per­a­ture idea is what gives sim­u­lated an­neal­ing its name. The tem­per­a­ture value of the al­go­rithm is anal­o­gous to the role of tem­per­a­ture when an­neal­ing metal. To an­neal metal, the metal is first heated and then slowly cooled over time. This cool­ing process al­lows the metal to have im­proved duc­til­ity and re­duced brit­tle­ness. By en­dur­ing a long cool­ing pe­riod, the metal is im­proved. Similarly, the rate at which the “heat” of the al­go­rithm is de­creased de­ter­mines the ac­cu­racy of the ap­prox­i­ma­tion.

The “heat” of the al­go­rithm af­fects the prob­a­bil­ity that it will make a bad move on pur­pose. Surprisingly, the al­go­rithm some­times needs to make bad moves to reach the global op­ti­mum. Why? There is a chance that it might get stuck on a lo­cal op­ti­mum. By mak­ing bad moves, the al­go­rithm can po­ten­tially stum­ble its way out of lo­cal op­ti­mum (but this is­n’t guar­an­teed).

You can choose any num­ber to be the ini­tial tem­per­a­ture. The tem­per­a­ture sched­ule then de­ter­mines how the tem­per­a­ture value will de­crease every it­er­a­tion. Like choos­ing neigh­bor­ing states, there is no right way to de­crease the tem­per­a­ture. The only rule of thumb is that the slower the tem­per­a­ture is dropped, the bet­ter the ap­prox­i­ma­tion will be. The trade­off is that you’ll have to wait longer to get an an­swer. Thus, choose some­thing that short enough for your at­ten­tion span and long enough to get a good an­swer.

For our toy ex­am­ple, I choose to sub­tract an ar­bi­trary num­ber every it­er­a­tion.

schedule <- function(t) {
  t - .00009
}

Once again, feel free to get cre­ative and try dif­fer­ent ini­tial val­ues and sched­ules for your im­ple­men­ta­tion.

Probabilistically Making Moves

Every it­er­a­tion a neigh­bor­ing state is gen­er­ated. The al­go­rithm will al­ways re­place the cur­rent state with the new state if it’s bet­ter. As we just dis­cussed, the new state may be worse–drag­ging us fur­ther away from an op­ti­mum. Yet, we want to ran­domly make some bad moves with some prob­a­bil­ity. The key to sim­u­lated an­neal­ing is that this prob­a­bil­ity is not con­stant. Rather, it on av­er­age de­creases over time.

The prob­a­bil­ity of choos­ing the new state is de­fined as the fol­low­ing for­mula:

Two pa­ra­me­ters af­fect the prob­a­bil­ity The “change in E” pa­ra­me­ter is the dif­fer­ence be­tween the func­tion ap­plied with the po­ten­tial new state and the func­tion ap­plied with the cur­rent state. The big­ger this dif­fer­ence, the smaller the prob­a­bil­ity. The next pa­ra­me­ter is the tem­per­a­ture value, which we’ve al­ready seen. Obviously, for every it­er­a­tion of the al­go­rithm, the tem­per­a­ture value will de­crease–caus­ing the prob­a­bil­ity to de­crease slowly.

As the tem­per­a­ture gets closer to zero, the prob­a­bil­ity of choos­ing a bad move will also get closer to zero. Towards the end, the al­go­rithm will mostly take greedy moves to­wards bet­ter states. Once the tem­per­a­ture hits zero, then the al­go­rithm is done. We can get our an­swer and go home.

Put Everything Together

We’ve cov­ered all the im­por­tant parts of this al­go­rithm. Let’s put every­thing to­gether and see it in ac­tion. Here is the full code for my so­lu­tion.

library(ggplot2)

simulated_annealing <- function(gen_successor, get_value, schedule) {
  generated_states <- c()
  temp <- 10 # initial temperature
  current_state <- .5 # initial state
  current_state_value <- get_value(current_state)
  
  while(temp > 0) {
    generated_states <- append(generated_states, current_state)
    temp <- schedule(temp)
    next_state <- gen_successor(current_state)
    next_state_value <- get_value(next_state)
    change_in_e = next_state_value - current_state_value
    
    if(change_in_e < 0) {
      current_state <- next_state
      current_state_value <- next_state_value
    }
    else {
      probability <- exp(-change_in_e / temp)
      samp <- runif(1)
      
      if(samp <= probability) {
        current_state <- next_state
        current_state_value <- next_state_value
      }
    }
  }
  
  return(generated_states)
}

gen_successor <- function(x) {
  TWO_PI <- 2 * pi
  successor <- x + runif(1, -.1, .1)
  
  if(successor < 0) {
    return(TWO_PI + successor)
  }
  else if(successor > TWO_PI) {
    return(0 + (successor - TWO_PI))
  }
    
  successor
}

get_value <- function(x) {
  sin(x)
}

schedule <- function(t) {
  t - .00009
}

generated_states <- simulated_annealing(gen_successor, get_value, schedule)
index <- seq(length(generated_states))
sin_values <- sin(generated_states)

data <- data.frame(
  index = index, 
  x = generated_states, 
  sin = sin_values
  )

ggplot(data, aes(x = index, y = sin)) +
  geom_smooth() +
  ylab("Sine value") +
  xlab("Iteration")

I de­cided to cap­ture the state for every it­er­a­tion and graph the re­sult. As you can see, the ini­tial state is nowhere near -1. As the al­go­rithm pro­gresses, it be­gins to slowly con­verge. Towards the end, the prob­a­bil­ity of choos­ing bad moves is near 0, which is why the down­ward slope for later it­er­a­tions is steep.

For this run, I got a fi­nal state of 4.703329. The cor­re­spond­ing sine value is -0.9999590, which is close to -1. These two val­ues match what we ex­pected.