5. Optimização Monte Carlo

Este capítulo é o equivalente para problemas de otimização do que o Capítulo 3 é para problemas de integração. Distinguimos entre dois usos separados de variáveis aleatórias geradas por computador para resolver problemas de otimização.

O primeiro uso, conforme vai ser visto na Seção 5.3, é produzir técnicas de busca estocástica para atingir o máximo ou mínimo de uma função, elaborando técnicas de exploração aleatória na superfície dessa função que evitem ficar presos em máximos ou mínimos locais e são suficientemente atraídos pelo máximo ou mínimo global. O segundo uso, descrito na Seção 5.4, está mais próximo do Capítulo 3 em que a simulação é usada para aproximar a função a ser otimizada.



5.1 Introdução

Os problemas de otimização podem ser vistos principalmente como um de dois tipos: ou precisamos encontrar o extremo de uma função \(h(\theta)\) sobre um domínio ou encontrar as soluções para uma equação implícita \(g(\theta) = 0\) sobre um domínio \(\Theta\). Ambos os problemas são intercambiáveis até certo ponto, pois o segundo é um problema de minimização para uma função como \(h(\theta) = g^2(\theta)\), na dimensão um; enquanto o primeiro é equivalente a resolver \(\partial h(\theta)/\partial\theta = 0\), assumindo que a função \(h\) pode ser diferenciada.

Portanto, focamos apenas no problema de maximização \[\begin{equation} \tag{5.1} \max_{\theta \in \Theta} h(\theta) \end{equation}\] uma vez que um problema de minimização pode ser tratado como um problema de maximização ao substituir \(-h\) ou \(1/h\) por \(h\).

Embora a escolha da transformada seja inócua do ponto de vista matemático, uma vez que o argumento do mínimo não muda, seu impacto no método de aproximação usado para encontrar esse argumento está longe de ser inócuo e a seleção da transformada deve ser considerada com cuidado.

Semelhante ao problema de integração tratado no Capítulo 3, o problema de otimização (5.1) pode ser processado por meios numéricos ou estocásticos. Na perspectiva numérica, o desempenho é altamente dependente das propriedades analíticas da função de destino, como convexidade, limitação e suavidade; enquanto as propriedades de \(h\) desempenham um papel menor em abordagens baseadas em simulação. Portanto, se \(h\) for muito complexo para permitir um estudo analítico ou se o domínio \(\Theta\) for muito irregular, o método de escolha é a abordagem estocástica.

O capítulo inteiro aborda a questão de encontrar extremos usando técnicas estocásticas, mas vamos apontar aqui que os problemas de otimização são genericamente mais difíceis de resolver do que os problemas de integração. Ou seja, os primeiros são muito mais locais do que os últimos, pelo menos do ponto de vista computacional (já que, do ponto de vista matemático, a informação sobre os extremos de uma função é tão global quanto a de suas integrais). Em outras palavras, é mais difícil apontar um único ponto extremo em um domínio do que a média de uma função regular no mesmo domínio. Observe que, de acordo com o restante do livro, este capítulo trata de problemas de otimização relacionados à estatística, com domínios quase sempre contínuos. Isso significa que problemas combinatórios difíceis que levam à otimização em um conjunto finito (se grande), como o problema do caixeiro viajante (ver, por exemplo, Spall, 2003, Robert e Casella, 2004) não são considerados aqui.


5.2 Métodos de otimização numérica

No R, existem diversas funções incorporadas para resolver problemas de otimização. O mais simples é optimize ou optimise, que processa alvos unidimensionais.


Exemplo 5.1

Ao considerar a maximização da função de verossimilhança Cauchy \(C(\theta,1)\), \[ \ell(\theta; x_1,\cdots,x_n) = \prod_{i=1}^n \dfrac{1}{1+(x_i-\theta)^2}, \] a sequência de máximos, ou seja, dos estimadores de máxima verossimilhança está convergindo para \(\theta^*=10\) quando \(n\) vai para \(\infty\).

Isso é refletido pela Figura 5.1 (esquerda), que corresponde ao código:

onde a log-verossimilhança é maximizada sequencialmente à medida que a amostra aumenta. No entanto, ao olhar diretamente para a verossimilhança, usar a otimização eventualmente produz uma sequência divergente, pois a verossimilhança fica muito pequena em torno de n = 300 observações, embora as sequências de MLEs sejam as mesmas até esse ponto para ambas as funções. Quando substituímos a verossimilhança suave por uma versão perturbada, como em 

set.seed(8643)
mi = xm = rcauchy(500, location = 10)
f = function(y){-sum(log(1+(x-y)^2))}
for (i in 1:500){ 
  x = xm[1:i]
  mi[i] = optimise(f, interval=c(9.5,11), maximum=TRUE)$max}
par(mfrow=c(1,2))
plot(1:500,mi, type = "l", col = "brown", xlab = "Iterations", ylab = expression(hat(theta)))
grid()
f = function(y){-sin(y*100)^2-sum(log(1+(x-y)^2))}
for (i in 1:500){ 
  x = xm[1:i]
  mi[i] = optimise(f, interval=c(9.5,11), maximum=TRUE)$max}
plot(1:500,mi, type = "l", col = "brown", xlab = "Iterations", ylab = expression(hat(theta)))
grid()

Figura. 5.1. (esquerda) Sequência de estimadores de máxima verossimilhança obtida de 500 simulações da distribuição Cauchy \(C(0,10)\) obtida pela aplicação de optimize para a log-verossimilhança; (à direita) as mesmas sequências ao usar uma verossimilhança perturbada.



5.3 Pesquisa estocástica


5.3.1 Uma solução básica

Uma maneira natural, embora rudimentar, de usar simulação para obter uma aproximação para a solução de (5.1) é simular pontos sobre \(\Theta\) de acordo com uma distribuição arbitrária \(f\), positiva em todos os lugares até que um valor suficientemente alto de \(h(\theta)\) seja observado.

Esta solução pode ser muito ineficiente se \(f\) não for escolhida em conexão com \(h\), mas, dado um número infinito de simulações e alguns requisitos de regularidade no problema, incluindo a compacidade do domínio \(\Theta\), é obrigado a convergir (Spall, 2003).

Por exemplo, se \(\Theta\) for limitado, podemos simular a partir de uma distribuição uniforme em \(\Theta\), \(U_1,\cdots,U_m\sim U_\Theta\) e usar \(h_m = \max\big(h(u_1),\cdots,h(u_m)\big)\) como uma aproximação para a solução de (5.1).


Exemplo 5.3

Lembre-se da função simples e regular, mas altamente variável \[ h(x) = \big( \cos(50 x) + \sin(20 x)\big)^2 \] definida em \([0,1]\), vista pela primeira vez no Exemplo 3.3.

Uma chamada para optimise fornece uma identificação do máximo em \(x^* = 0.379\) com um valor de \(h(x^*) = 3.8325\). Se quisermos avaliar a variabilidade de um amostrador uniforme, podemos usar várias sequências uniformes cujo resultado é mostrado na Figura 5.3.

Embora o valor inicial \(h(u_1)\) da sequência seja altamente variável, o intervalo reduz muito rapidamente e, após \(10^3\) iterações, a pior sequência entre as \(10^3\) sequências paralelas está dentro de 0.24 do máximo verdadeiro. 

set.seed(6578)
h=h=function(x){(cos(50*x)+sin(20*x))^2}
rangom=h(matrix(runif(10^6),ncol=10^3))
monitor=t(apply(rangom,1,cummax))
plot(monitor[1,],type="l",col="white",xlab="Iterations",ylab=expression(h(theta)))
polygon(c(1:10^3,10^3:1),c(apply(monitor,2,max), rev(apply(monitor,2,min))),col="grey")
abline(h=optimise(h,int=c(0,1),maximum=T)$ob, col="brown",lwd=3)
grid()

Figura. 5.3. Vaiação das \(10^3\) sequências de máximos sucessivos encontrados por amostragem uniforme aleatória em \(10^3\) iterações. O valor máximo verdadeiro é identificado pela linha marrom na parte superior do gráfico.


Usar uma distribuição uniforme sobre o domínio \(\Theta\) somene é relevante quando \(\Theta\) tem uma forma regular. Caso contrário, é mais eficiente simular a partir de uma distribuição uniforme sobre um domínio contendo \(\Theta\) e descartar as simulações fora de \(\Theta\).

Obviamente, esta solução cega, cega no sentido de que não leva \(h\) em conta, torna-se rapidamente impraticável à medida que a dimensão ou a complexidade do problema aumenta. Por exemplo, em uma configuração bayesiana, quando o tamanho de \(n\) de uma amostra iid \(X_1,\cdots,X_n\) de \(f(x|\theta)\) cresce, a distribuição a posteriori associada \(\pi(\theta|x_1,\cdots,x_n)\) fica cada vez mais concentrada em torno de sua moda, o que é cada vez mais difícil de se aproximar dessa forma. Isso, no entanto, nem sempre é o caso, como mostra o exemplo a seguir.


Exemplo 5.4 (Continuação do Exemplo 5.1)

Exemplo 5.4. Usando o mesmo modelo de Cauchy do Exemplo 5.1, podemos monitorar a discrepância entre as soluções fornecidas por optimise e por amostragem uniforme em [-5;5] à medida que o tamanho da amostra \(n\) aumenta de \(n = 1\) para \(n = 5001\)

A Figura 5.4 (topo) mostra o valor do argumento verdadeiro \(\theta^*\), conforme fornecido por optimise, contra o valor da amostra uniforme com a maior verossimilhança e não apresenta um aumento significativo do erro próximo a zero, pois há mais pontos nessa vizinhança. Da mesma forma, a Figura 5.4 (abaixo) mostra que o erro relativo entre o máximo verdadeiro e a aproximação na amostra uniforme não exibe uma tendência crescente.

Fig. 5.4. Comparação de uma maximização numérica e estocástica de uma verossimilhança de Cauchy em termos do tamanho da amostra via (topo) respectivas localizações das avaliações numéricas e estocásticas dos argumentos, plotadas ao longo da diagonal; (abaixo) erro relativo da avaliação estocástica em relação à avaliação numérica em função do tamanho da amostra.


Portanto, é mais proveitoso projetar o experimento de simulação em próxima conexão com \(h\), bem como com o domínio \(\Theta\). Intuitivamente, faz sentido aumentar a probabilidade de simular em regiões onde \(h\) é grande e diminuir essa probabilidade em regiões onde é pequeno. Isso significa criar uma distribuição de probabilidade relacionada com \(h\) de maneira não linear, mas com modos idênticos ou quase idênticos.


Exemplo 5.6

Considere minimizar a função construída artificialmente em \(\mathbb{R}^2\) \[\begin{equation} \begin{array}{rcl} h(x,y) & = & \big( x \sin(20 y)+ y \oc(20 x)\big)^2 \cosh\big( \sin(10 x)x\big) \\ & & \qquad + \big( x \cos(10 y)- y\sin(10 x)\big)^2 \cosh\big( \cos(20 y)y\big), \end{array} \end{equation}\] cujo mínimo global é 0, obtido em \((x,y) = (0,0)\). 

h=function(x,y){(x*sin(20*y)+y*sin(20*x))^2*cosh(sin(10*x)*x)+(x*cos(10*y)-y*sin(10*x))^2*cosh(cos(20*y)*y)}
x=y=seq(-3,3,le=435) #defines a grid for persp
z=outer(x,y,h)
par(bg="wheat",mar=c(1,1,1,1)) #bg stands for background
persp(x,y,z,theta=155,phi=30,col="green4", ltheta=-120,shade=.75,border=NA,box=FALSE)

Figura. 5.6. Representação via persp da função \(h(x,y)\) do Exemplo 5.6 em \([-3,3]^2\).

Uma vez que esta função tem muitos mínimos locais, conforme mostrado na Figura 5.6, ela não satisfaz as condições sob as quais os métodos de minimização padrão são garantidos para fornecer o mínimo global.

Por outro lado, a distribuição em \(\mathbb{R}^2\) com densidade proporcional a \(\exp\big(-h(x,y)\big)\) pode ser simulada, mesmo que não seja uma distribuição padrão, usando um algoritmo Accept-Reject baseado na distribuição uniforme; uma vez que \(h\) é positivo, o que está anulando o propósito de simular a partir de \(h\) em vez da distribuição uniforme sobre \(\Theta\) ou técnicas MCMC mais avançadas introduzidas posteriormente no Capítulo 6.



5.3.2 Método de gradiente estocástico

Dado que gerar simulações diretas a partir da função alvo \(h\), definida na seção anterior, costuma ser uma grande dificuldade; uma abordagem estocástica diferente para a maximização de \(h\) é explorar a superfície de \(h\) de maneira local, isto é, definindo uma sequência \(\{\theta_n\}\) movendo-se de \(\theta_n\) para \(\theta_{n+1}\) em uma etapa dependente, em vez de independentemente como no algoritmo básico de busca estocástica.

A dependência de \(\theta_{n+1}\) em \(\theta_n\) é frequentemente escolhida para ser linear, no sentido de que é representada como \[\begin{equation} \tag{5.3} \theta_{n+1}=\theta_n+\epsilon_n, \end{equation}\] onde \(\epsilon_n\) é a perturbação local do valor atual.

Em termos matemáticos, isso torna a sequência \(\{\theta_n\}\) uma Cadeia de Markov. Embora haja uma conexão entre esses métodos e os algoritmos MCMC, a ser visto nos Capítulos 6 e 7, a propriedade de Markov, no entanto, é menos importante no presente cenário simplesmente porque a matemática que justifica a convergência para um máximo global é muito avançada para ser considerada aqui; veja, por exemplo, Hàjek (1988) ou Haario and Sacksman (1991).

Ao considerar a implementação da ideia de atualizaç_ão local em (5.3), a perturbação \(\epsilon_n\) pode ser simulada a partir de uma distribuição arbitrária como uma distribuição \(Np(0, \sigma^2\pmb{I}_p)\) se \(\Theta\in\mathbb{R}^p\). No entanto, dado que estamos especificamente procurando o máximo de \(h\), usar alguma informação sobre \(h\) na construção da distribuição da perturbação certamente aumentará a eficiência do método. Em particular, faz sentido favorecer movimentos crescentes em \(h\) sobre movimentos decrescentes em \(h\), mesmo que o último não seja impossível se você quiser evitar máximos locais. Uma abordagem natural é usar o gradiente de \(h\), \(\nabla h\), se disponível.

Na otimização numérica, o método do gradiente é uma abordagem numérica determinística para o problema de otimização (5.1) relacionado ao método de Newton-Raphson já apresentado na Seção 5.2. Ele produz uma sequência \(\{\theta_n\}\) definida por \[\begin{equation} \tag{5.4} \theta_{n+1} = \theta_n+\alpha_n \nabla h(\theta_n), \qquad \alpha_n > 0, \end{equation}\] que converge para a solução exata de (5.1), \(\theta^*\), quando o domínio \(\Theta\in\mathbb{R}^d\) e a função \((-h)\) é convexa, assumindo assim que existe um único máximo para várias escolhas da sequência decrescente \(\{\alpha_n\}\), ver Thisted (1988). Para problemas menos regulares, é mais provável que a sequência gradiente fique presa em um extremo local da função \(h\).

O método de gradiente estocástico aproveita este método para construir a perturbação em (5.3). Por exemplo, a proposta de diferenças finitas é construir um substituto numérico para o verdadeiro gradiente \[\begin{equation} \nabla h(\theta_n) \approx \dfrac{h(\theta_n+\beta_n\zeta_n)-h(\theta_n-\beta_n\zeta_n)}{2\beta_n}\zeta_n \, = \, \dfrac{\Delta h(\theta_n,\beta_n\zeta_n)}{2\beta_n}\zeta_n, \end{equation}\] onde \(\{\beta_n\}\) é uma segunda sequência decrescente e \(\zeta_n\) é distribuído uniformemente sobre a esfera unitária \(||\zeta|| = 1\).

Em contraste com a abordagem determinista, a atualizar _{n+1}=_n+ h(_n,_n_n)_n, \end{equation} não prossegue ao longo da inclinação mais íngreme de \(h\) em \(\theta_n\), pois cada vez que ele escolhe a direção é aleatoriamente, mas essa propriedade é geralmente uma vantagem no sentido de que pode evitar ficar preso em máximos locais ou em pontos de sela de \(h\).

Se \(\{\theta_n\}\) definido por (5.5) converge ou não para o argumento de (5.1) dependerá muito da escolha das sequências \(\{\alpha_n\}\) e \(\{\beta_n\}\). Por exemplo, \(\alpha_n\) precisa diminuir lentamente o suficiente para 0 para a série correspondente \(\sum_n \alpha_n\) divergir, enquanto \(\beta_n\) deve diminuir ainda mais lentamente para a série \(\sum_n \big(\alpha_n/\beta_n\big)^2\) convergir (Spall, 2003).


Exemplo 5.7 (Continuação do Exemplo 5.6)

Aplicamos a construção iterativa (5.5) para a função multimodal \(h(x,y)\) com diferentes sequências de \(\alpha_n\)’s e \(\beta_n\)’s para verificar seu impacto.

Uma regra de parada natural para o algoritmo é verificar a estabilização na sequência \(\theta_n\) , levando à implementação no R onde as sequências \(\alpha_n\)’s e \(\beta_n\)’s devem ser inseridas. Ao executar este código R, na verdade, tivemos que incluir um loop de segurança interno para proteger contra avaliações divergentes do gradiente que ocorrem de tempos em tempos e levam o programa a abortar. 

start=c(.65,.8)
theta=matrix(start,ncol=2)
dif=iter=1
# scenario 1
alpha=matrix(1/log(2),ncol=1)
beta=matrix(1/log(2)^0.1,ncol=1)
#while (dif>10^-5){
# zeta=rnorm(2)
# zeta=zeta/sqrt(t(zeta)%*%zeta)
# grad=0.5*alpha[iter]*zeta*(outer(theta[iter,1]+beta[iter]*zeta[1],theta[iter,2]+beta[iter]*zeta[1],h)-
#                          outer(theta[iter,1]-beta[iter]*zeta[2],theta[iter,2]-beta[iter]*zeta[2],h))/beta[iter]
# theta=rbind(theta,theta[iter,]+grad)
# dif=sqrt(t(grad)%*%grad)
# iter=iter+1
# alpha=rbind(alpha,1/log(iter+1))
# beta=rbind(beta,1/log(iter+1)^0.1)}
#scale=sqrt(t(grad)%*%grad)
scale = 10
#while (scale>0){
# zeta=rnorm(2);zeta=zeta/sqrt(t(zeta)%*%zeta)
# grad=alpha[iter]*zeta*(h(theta[iter,]+beta[iter]*zeta)-
# h(theta[iter,]-beta[iter]*zeta))/beta[iter]
# grad=alpha[iter]*zeta*(outer(theta[iter,1]+beta[iter]*zeta[1],theta[iter,2]+beta[iter]*zeta[1],h)-
#                          outer(theta[iter,1]-beta[iter]*zeta[2],theta[iter,2]-beta[iter]*zeta[2],h))/beta[iter]
# theta=rbind(theta,theta[iter,]+grad)
# iter=iter+1
# alpha=rbind(alpha,1/log(iter+1))
# beta=rbind(beta,1/log(iter+1)^0.1)
# scale=sqrt(t(grad)%*%grad)}



5.5 Exercícios

5.3 Quando é de nido em R2 pela restrição x2(1 + sin(3y) cos(8x)) + y2(2 + cos(5x) cos(8y)) 1 ; proponha uma simulação uniforme simples em um domínio maior e avalie o desempenho desse método por meio do número médio de pontos rejeitados.

5.5 Reproduza a análise do Exemplo 5.7 acima sobre o impacto da dinâmica das sequências \(\alpha_j\) e \(\beta_j\) na convergência do método das diferenças finitas no ajuste da versossimilhança de mistura do Exemplo 5.2.