Simulated Annealing

Introdução

Encontar ótimos globais nem sempre é uma tarefa simples, especialmente quando lidamos com muitas dimensões. Geralmente utilizamos métodos que envolvem derivadas para solucionar esses problemas, como o método do gradiente por exemplo. Acontece que algumas funções podem ser muito difíceis de derivar analiticamente. Além disso, certas funções tem muitos ótimos locais, e encotrar todos os ótimos locais já pode ser algo muito custoso computacionalmente, por isso as vezes utilizamos técnicas de busca aproximativa.

Podemos citar dois tipos de abordagem para o problema de maximizar ou de minimizar uma função objetivo: uso de algoritmos exatos e uso de algoritmos aproximativos (meta-heurísticas). Onde os algoritmos exatos fornecem resultados que garantem a solução ótima, mas as vezes tem formulação demasiadamente complicada e nem sempre conseguem produzir uma boa solução viável rápidamente, ou seja, envolve um gasto de tempo proibitivo para chegar na solução ótima.

Uma heurística é uma estratégia para encontrar boas soluções de um problema, que muitas vezes é de otimização combinatória. Alguns algoritmos utilizam não apenas uma heurística, mas um conjunto de heurísticas, e normalmente heurísticas que melhoram outras heurísticas, chamamos isso de meta-heurística. Esse tipo de técnica apesar de não garantir valores ótimos globais, é capaz de fornecer soluções muito boas em um tempo de convergencia adequado e com pouco esforço computacional.

Os algoritmos de subida da montanha (Hill Climbing), busca tabu (tabu search) e arrefecimento simulado (simulated annealing), são exemplos de técnicas para otimizar funções de maneira independente de diferenciação, se baseando em observar valores na vizinhança de um chute inicial e aceitar o novo valor se houver melhoramento na função objetivo, e repetir esse procedimento interativamente gerando valores na vizinhança dos novos pontos, discutiremos quais são as vantagens e desvantagens desses métodos ao decorrer deste trabalho e mostraremos alguns exemplos de implementação no R.

Em estatística, um problema comum é o de estimar parâmetros de uma população, para isso geralmente procuramos maximizar a função de verossimilhança da amostra. Veremos que em alguns casos, não é facil obter uma estimação satisfatória com o método de Newton-Raphson, neste caso aplicaremos o algoritmo de arrefecimento simulado (simulated annealing) para encontrar uma estimativa do parâmetro.

Arrefecimento Simulado

Nesta seção iremos apresentar o algoritmo de arrefecimento simulado e discutir qual a motivação de utiliza-lo, especialmente em estatística.

Arrefecimento Simulado (Simulated Anealing) é um algorítimo de otimização estocástica para busca local que, diferentemente de outros algoritimos dessa classe, tem maior facilidade de escapar de pontos ótimos locais.

Antes de falarmos sobre busca local, precisamos definir o conceito de Vizinhança.

Definição Uma função \(N:S \rightarrow 2^S\) que atribui a cada solução \(x\) de \(S\) um conjunto de soluções $N(x) S $ é chamada de vizinhança de x.

Exemplo Em um espaço de buscas contínuo \(S\) a vizinhança de \(x\), \(N(x)\), seria a bola aberta de centro \(x\) e raio \(\epsilon\). i.e: \(N(x) = \{x^{\prime} \in S : | \ \ || x^{\prime} - x|| < \epsilon\}\)

Busca Local

Busca local é um conjunto de técnicas de otmização que visam resolver problemas do tipo $,  x S $, onde inicia-se com uma solução \(x_0\) e a cada iteração substitui-se a solução atual \(x_k\) por uma solução na vizinhança que melhora a solução atual de forma que $f(x) < f(x^{}),  x N(s^{}) $. A busca acaba quando todos os candidatos vizinhos são piores que a solução atual, o que implica que um ótimo local foi alcançado.

Os principais algoritmos dessa classe são:

Subida da montanha ou Hill Climbimg

O conceito desse algoritmo é de fato uma analogia a subir uma montanha. Ao dar um passo em alguma direção e perceber que desceu, pode significar que aquele caminho é o de descida da montanha, o mesmo raciocínio se aplica ao subir a montanha. A figura a seguir ilustra esse caso:

Busca Tabu ou Tabu Search

O problema do algoritmo anterior é que ele cai rapidamente em ótimos globais, e isso pode ser muito ruim, já que dessa forma ele desconsidera outras soluções potencialmente melhores. Para melhorar isso adaptamos ligeiramente a técnica de subida da montanha da seguinte maneira: encontramos um ótimo local mas continuamos ``caminhando na montanha“, guardando na memória os valores de cada passo até encontrar uma solução melhor ou atingir o limite de passos. Esse procedimento é chamado de busca tabu, esse nome é devido ao fato de ser um tabu enfrentar uma piora na solução apostando que ela pode melhorar no futuro. Chamamos os pontos armazenados na memória como elementos tabu. A figura a seguir ilustra um caso bem sucedido de aplicação de busca tabu:

Arrefecimento Simulado ou Simulated Annealing

Esta técnica é o foco deste trabalho, conforme já foi dito o arrefecimento simulado foi baseado no porcesso de fundição de metais para fazer ligas. A mecânica de funcionamento dele é muito parecida com a técnica de subida de montanha ou com a busca tabu, a diferença é que para escapar dos ótimos locais o algoritmo pode aceitar um valor que piora a função objetivo com uma certa probabilidade.

Proposto inicialmente por Scoot Kirkpatrick et al. (1983) em Optimization by Simulated Annealing, com a proposta de modificar o algoritmo de Metropolis (Gibbs, 1953). O objetivo deste método é minimizar funções.

O processo inicia com um valor de temperatura T elevado, e a cada interação geram-se novas soluções até que o equilíbrio térmico seja atingido para aquela temperatura. Então a temperatura é rebaixada e o processo segue até que se chege ao congelamentou seja, quando a temperatura chega próxima a 0.

A tanto o número de interações em cada temperatura, quanto o tamanho do passo de rebaixamento da temperatura são parâmetros que devem ser ajustados empiricamente.

A analogia com o processo físico, indica que os estados possíveis do metal são o conjunto de soluções no espaço de busca. A energia em cada estado é o valor da função objetivo.

A cada interação um novo estado é gerado a partir de uma modificação aleatória no estado atual. Usaremos \(\Delta_\theta = f(\theta_{i+1})-f(\theta_i)\) para indicar a diferença entre o novo estado e o estado anterior, em que \(\theta\) é o parâmetro a ser estimado. Note que, se \((f(\theta_{i+1})-f(\theta_i)) \geq 0\) significa que aumentamos o valor da função objetivo, e no caso, como gostariamos de minimizar, então esta solução não é interessante, mas para evitar cair em um ótimo local podemos aceitar essa nova estimativa com a espectativa de chegar em soluções ainda melhores no futuro. A probabilidade \(p\) de aceitar essa nova solução é definida por:

\[\begin{equation} p = e^{-(\Delta_\theta)/kT} \end{equation}\]

  p = function(t){
    p = exp(-1/t)
  }

  t = seq(0,50,0.5)
  
  plot(t, p(t),type = "l", xlab = "Temperatura (t)", ylab = "exp(-1/t)")

Perceba que quanto maior a temperatura T maior a chance de aceitarmos a nova solução e quando T tende para infinito p tende para 1. Logo, não é difícil ver que, quando a temperatura aumenta a probabilidade de aceitar novos estados é maior.

Utilizando esse fato, podemos controlar a chance de aceitação de piores valores por meio da temperatura T. Apresentamos a seguir, o algoritmo em alto nível:

Algoritmo: Arrefecimento Simulado

1. Escolha um ponto inicial \(\theta_0\). Se isso não for possível, então escolha aleatoriamente.

2. Calcule \(f(\theta_0)\), o valor inicial da função.

3. Inicie a partir do ponto atual, escolhendo um ponto aleatório na vizinhança, dentro de uma esfera \(n\)-dimensional. Onde n é a dimensão do problema. Isso esfecifica uma direção aleatória.

4. Tome um passo de tamanho \(\Delta_\theta\).

5. No novo ponto \(\theta_1\), calcule o valor \(f(\theta_1)\) e aceite o passo com probabilidade 1 se \(\Delta_\theta = (f(\theta_1)-f(\theta_0)) \leq 0\). Se \(\Delta_\theta > 0\), aceite o passo com propabilidade

\[ p = e^{-(\Delta_\theta)/T} \]

6. Se o passo anterior não aceita o novo valor, então vá para o passo 3. Caso contrário use a nova solução no passo 3 e repita.

Estudo de casos

Para exemplificar o algoritmo fizemos alguns estudos de casos, em que minimizamos algumas funções de verossimilhança no R com o algoritmo de arrefecimento simulado.

Considere a densidade da cauchy:

\[ f(x) = \frac{\beta}{\pi\{\beta^2+(x-\alpha)^2\}} \]

para \[-\infty < x < \infty, \beta \geq 0, -\infty < \alpha < \infty \]

Para estimar \(\alpha\) e \(\beta\), de uma amostra aleatória \(x_1, x_2, ..., x_n\),podemos maximizar a função de verossimilhança, ou a função log-verossimilhança, que é:

\[ l = n \, log (\beta) - {\sum_{i=1}^{n} \{\beta^2 + (x_i-\alpha)^2\}} - n \, log (\pi) \]

Com \(\beta\) fixo, a função de máxima verossimilhança estima \(\alpha\) maximizando:

\[ l(\alpha) = \sum_{i=1}^{n}{log \, \{\beta^2 + (x_i-\alpha)^2\}} \]

que é equivalente a minimizar:

\[ l(\alpha) = -\sum_{i=1}^{n}{log \, \{\beta^2 + (x_i-\alpha)^2\}} \]

Suponha a amostra (-4.20, -2.85, -2.30, 1.02, 0.70. 0.98, 2.72, 3.50) com \(\beta = 1\) a seguir podemos visualizar o plot da função log-verossimilhança desta amostra:

cauchy.l = function(alpha){
  x = c(-4.2,-2.85,-2.3,-1.02,0.7,0.98,2.72,3.5)
  return(sum(log(0.01+(x-alpha)^2)))
}
x = seq(-6,6,.1)
y = sapply(as.list(x),cauchy.l)
plot(x,y, type = "l", ylab = "Cauchy.l (beta = 0.1, alpha)", xlab = "alpha")

sa <- function(fun, N, rho, t, theta){
  y = NULL
  theta.old = theta
  
  repeat{
    s = 0
    for(k in 1:N){
      f.old = fun(theta.old)
      theta.new = runif(1,-1,1) + theta.old
      f.new = fun(theta.new)
      
      if(f.new<f.old){
        theta.star = theta.new
        s = s+1
      }else{
       u = runif(1)
       b = (u <= exp(-(f.new-f.old)/t))
       theta.star = theta.new * b + theta.old * (1-b)
       s = s+b
      }
      theta.old = theta.star
      y = c(y, theta.star)
    }
    
    write(y, file = "anneal.result", ncol = 1, append = T)
    r = y[length(y)]
    y = NULL
    if(s==0) break
    t = rho*t
    if(t<.1) break
  }
  return(r)
}

Agora iremos rodar o arrefecimento simulado para estimar o parâmetro \(\alpha\). com N = 50 interações para cada temperatura t, o parâmetro rho = 0.95 é o quanto a temperatura t decai quando fazemos t = rho*t. Utilizamos a temperatura t = 10 e um chute inicial theta = 4.

r = sa(cauchy.l,50,.95,10,4)

r

## [1] 0.7426755

plot(x,y, type = "l", ylab = "Cauchy.l (beta = 0.1, alpha)", xlab = "alpha")

abline(v=r, col = "red")

Valor estimado de alfa com o arrefecimento simulado

Podemos visualisar o histórico de valores encontrados para entender melhor o como o algoritmo, visita os candidatos a solução:

Neste caso, o valor encontrado é aparentemente muito bom. Poderiamos fazer o mesmo para outras funções e comparar os resultados obtidos com o valor dos estimadores. Por exemplo, considere uma amostra \(x_1, x_2, ..., x_n\) de uma população normal, com variância igual a 4 e \(\mu\) desconhecido. Seja a amostra (12,15,9,10,17,12,11,18,15,13), a função log-verossimilhança é dada por

\[ l(\mu) = -5 log (8 \pi) - \frac{1}{8}(\sum_{i=1}^{10}x_i^2-2 \mu \sum_{i=1}^{10}x_i+10 \mu^2) \]

Inserindo a função de máxima verossimilhança com sinal negativo para a amostra exemplo no R:

  normal.l = function(mu){
  x = c(12,15,9,10,17,12,11,18,15,13)
  return((sum(x^2)-2*mu*sum(x)+10*mu^2))
}

Agora podemos fazer o plot e verificar os resultados de estimação de \(\mu\) com o arrefecimento simulado em respeito ao estimador média amostral e visualizar as soluções gráficamente.

x = seq(11,15,.1)
y = sapply(as.list(x),normal.l)
plot(x,y, type = "l", ylab = "Normal (sigma2 = 4, mu)", xlab = "mu")

Agora podemos tentar minimizar essa função que é simples via arrefecimento similado para verificar se os resultados obtidos são satisfatórios:

x = c(12,15,9,10,17,12,11,18,15,13)
normal.l = function(mu){
  x = c(12,15,9,10,17,12,11,18,15,13)
  return((sum(x^2)-2*mu*sum(x)+10*mu^2))
 }
u = seq(-6,40,1)
v = sapply(as.list(u),normal.l)
sol <- sa(normal.l,50,0.85,30,5)
##Solução analítica
mean(x)

## [1] 13.2

##Solução encontrada via Simulated Annealing
sol

## [1] 13.37643

plot(u,v, type = "l", ylab = expression("-1*log- verossimilhança  Dist. Normal"), xlab = "alpha")
abline(v= sol)
abline(v = mean(x), lty = 2, col = "red")

Considerações Finais

Vimos nas sessões anteriores que existem diversos problemas no mundo real que encontram-se sem solução. Vimos tambem alguns algoritmos de metaheuristicas, Entretanto, estes foram uma pequan quantidade em relação ao total de algoritmos que podemos construir.

O método de arrefecimento simulado deve ser empregado com cautela, já que se for possível obter um estimador com boas características pelas técnicas mais convencionais, é preferivel o estimador em detrimento do simulated annealing. Já que os resultados podem não ser tão precisos quando os parâmetros do algoritmo não são bem calibrados.

Tendo em vista a enorme gama de possibilidades, um passo futuro seria estudar e implementar novas metaheuristicas, tendo como objetivo investigar sua eficiência em relação aos já implementados nesse trabalho.