Introdução

Esta nota técnica apresenta alguns detalhes técnicos e implementações de algorítmos em R (gradiente descendente e método de Newton) que talvez possam facilitar um melhor entendimento do exemplo apresentado uso da função Penalty no texto de Bryan e Shibberu, passado em classe.

Descrição do problema e resultado teórico usando condições KKT

A solução desse problema é óbvia, \(x^*=1\).

Pergunta:Como podemos verificar teoricamente que \(x^*=1\) é de fato a solução ótima global?

Resposta: podemos verificar se as condições KKT se verificam para esse caso. Para ficar no formato de um problema de maximização (para facilitar a comparação com o formato das condições apresentada no Hoy et al.), o mesmo problema seria dado por

em que:

As condições KKT serão dadas, no caso, por:

Observe que não há nesse caso restrição de domínio \(x\ge 0\) usada no no modelo geral apresentado no Hoy et al. 2001.

Note que a solução \(x^*=1\) leva a \(\lambda^*=4\), pela primeira condição, uma solução que atende às 4 condições:

Conclusão: Como a função objetivo \(f(x)=-x^4\) e a restrição \(g(x)=1-x\) são côncavas, o atendimento às condições KKT é suficiente e necessário para concluirmos que há de fato um máximo global (ou seja o mínimo global de \(f(x)=x^4\) condicionado pela restrição \(x\ge 1\)).

Achando a solução pelo método Penalty

A idéia geral do método é transformar o problema de otimização condicionada num problema de otimização não condicionada, introduzindo uma penalidade que que cresce na medida que cresce o desvio no atendimento à restrição (ou restrições).

Para o problema específico

A formulação não-condicionada envolvendo a função Penalty seria:

Em que \(\phi(k,x)\) é uma função que adiciona uma penalidade crescente para valores de \(x<1\), que pode ser acentuada pelo parâmetro \(k\). A especificação proposta por Bryan e Shibberu (com algumas modificações na formulação original, para facilitar o entendimento), com o apoio de uma função indicadora é dada por:

x<-seq(0,2,0.01)
k<-10
y<-ifelse(x<1,k*(1-x)^3,0)
plot(x,y,type="l",col="red")
title("Função penalty utilizada")
abline(h=0,lty=3)

Note que a função é contínua e convexa (verifique como exercício), aumentando a com menores valores de \(x\), a partir de \(x=1\). O valor de \(k\) pode ser utilizado para “regular” o efeito da penalidade aplicada.

Com essa especificação, a função original \(f(x)\) (preto para \(x<1\)) e a nova função \(\hat f(x)\) (vermelho), considerando \(k=10\), são representadas na figura a seguir.

plot(x,x^4,type="l",col="black")
title("Função original e com penalty (k=10)")
lines(x,x^4+y,col="red")
abline(h=0,lty=3)

Observe e verifique que a nova função \(\hat f(x)\) é estritamente convexa, com mínimo global estabelecido num ponto estacionário. Teoricamente, para valores crescentes de \(k\), a solução convergirá para a solução desejada no problema original.

Resolvendo o problema pelo método do Gradiente Descendente (GD)

Para um dado \(k\) podemos resolver o problema de otimização utilizando um algoritmo numérico, como o gradiente descendente (lembre que no caso o gradiente tem uma só dimensão, que é a de maior crescimento da função). No caso, partimos de um valor inicial \(x_0\) e obtemos uma recursão

em que para valores apropriados de \(\delta\) tenderá a convergir para a solução ótima global (que existirá em função da convexidade da função). Para implementarmos o método, necessitamos a derivada da função, que será dada por:

Para interrompermos o processo recursivo, testamos a cada passo o valor absoluto da derivada (norma do gradiente no caso unidimensional), verificando se ela está suficientemente próxima de zero, dentro de limites de tolerância pré-estabelecidos. Frequentemente há necessidade de adequação dos parâmetros \(k\) e \(\delta\) para bons resultados.

## Implementação do gradiente descendente
## partindo de x=2 com k=1000000
##
options(digits=15)
phi<-function(x){ifelse(x<1,(1-x)^3,0)}
mfunc<-function(k,x){x^4+k*phi(x)}
dmfunc<-function(k,x){4*x^3-ifelse(x<1,k*3*(1-x)^2,0)}
k<-1000000
x<-2 #ponto inicial
tol<-0.0000001
passos<-0
while(abs(dmfunc(k,x))>tol){
  passos<-passos+1
  x<-x-0.0001*dmfunc(k,x)
  #cat(x,dmfunc(k,x),"\n")
}
cat("iterações:",passos,"x=",x,"norma do gradiente=",abs(dmfunc(k,x)),"\n")
## iterações: 954 x= 0.998847295436723 norma do gradiente= 5.32514294881992e-08

Com essa simples implementação, com pouco mais de 900 iterações, atingimos um valor muito próximo de \(x=1\), que seria o resultado teoricamente correto. A partir de um certo ponto, há maior dificuldade de se melhorar os resultados em função de problemas de precisão. Em muitos casos, o uso de outros métodos, como o Método de Newton, pode apresentar resultados adequados, com menos iterações.

Resolvendo o problema pelo método de Newton

O método de Newton, para problemas de otimização, é aplicado pela recursão

exigindo também a derivada segunda, no caso univariado, ou o hessiano no caso multivariado.

em que para valores apropriados de \(\delta\) tenderá a convergir para a solução ótima global (que existirá em função da convexidade da função). Para implementarmos o método, necessitamos a derivada da função, que será dada por:

Para interrompermos o processo recursivo, testamos a cada passo o valor da derivada, verificando se ela está suficientemente próxima de zero, dentro de limites de tolerância pré-estabelecidos. Frequentemente há necessidade de adequação do parâmetro \(k\) para bons resultados.

## Implementação com método de Newton
## partindo de x=2 com k=1000000
##
phi<-function(x){ifelse(x<1,(1-x)^3,0)}
mfunc<-function(k,x){x^4+k*phi(x)}
dmfunc<-function(k,x){4*x^3-ifelse(x<1,k*3*(1-x)^2,0)}
d2mfunc<-function(k,x){12*x^2+ifelse(x<1,k*6*(1-x),0)}
k<-10000000
x<-2 #ponto inicial
tol<-0.0000001
passos<-0
while(abs(dmfunc(k,x))>tol){
  passos<-passos+1
  x<-x-dmfunc(k,x)/d2mfunc(k,x)
  #cat(x,dmfunc(k,x),"\n")
}
cat("iterações:",passos,"x=",x,"norma do gradiente=",abs(dmfunc(k,x)),"\n")
## iterações: 14 x= 0.999635051500615 norma do gradiente= 3.31845662060459e-11

Com essa simples implementação, com apenas 14 iterações, atingimos um valor muito próximo de \(x=1\), que seria o resultado teoricamente correto. Tal como na situação anterior, a partir de um certo ponto, há maior dificuldade de se melhorar os resultados em função de problemas de precisão. Em muitos problemas reais há necessidade de combinação de vários métodos de otimização para atingir o resultado desejado, assim como ajustes nos parâmetros que controlam os algoritmos.

Caso multivariado

No exemplo envolvendo o caso multivariado, a solução do problema pode seguir o mesmo caminho utilizado para o caso anterior. Inicialmente define-se as funções Penalty, para cada restrição, da forma indicada. A seguir, pode-se implementar a solução tanto pelo uso do Gradiente Descendente ou o Método de Newton, em suas versões multivariadas, que foram os métodos vistos no curso. Nesse caso será necessário a especificação do gradiente e o Hessiano (no caso do método de Newton), como fizemos no exemplo utilizando o WxMaxima.

Solução envolvendo método da Barreira (GD, Newton, Comb, optimize e Monte Carlo)

O fundamento e motivação são similares ao da função Penalty. Nesse caso é necessario partir de um ponto inicial possível (algo que pode ser complicado de encontrar em alguns problemas). Apesar dessa dificuldade, esse método é muito eficiente, sendo a base dos algorítmos interiores que solucionam problemas de programação linear envolvendo milhares de variáveis.

No exemplo anterior, resolvido pelo método da função Penalty, uma função logaritmica usualmente é utilizada na implementação do método da Barreira, de forma que o problema

é transformado numa versão não-condicionada, especificada por:

Note que partindo de um ponto \(x_0\) no domínio à direita de 1, a função estabelece uma penalidade que vai tendendo ao infinito na medida que \(x\rightarrow 1\). Isso é ilustrado no gráfico a seguir com a função original (preto) e a função com barreira (vermelho)

x<-seq(1.00001,2,0.00001)
b<-function(k,x){x^4-k*log(x-1)}
k<-2 # 
plot(x,b(k,x),type="l",col="red",ylim=c(0,30))
title("Função original e com barreira")
lines(x,x^4)
abline(h=0,lty=3)
abline(v=1,lty=3)

Na medida que o valor de \(k\) vai sendo reduzido, o problema com a barreira converge para a solução do problema original.

Para solucionar o problema podemos usar os mesmos procedimentos anteriores, usando a derivada primeira e segunda, no método do gradiente descendente e método de Newton, caracterizadas a seguir:

Note que \(b''(x)>0\) e que a função é estritamente convexa.

## Implementação com gradiente descendente
## partindo de x=2 com k=0.001
##
b<-function(k,x){x^4-k*log(x-1)}
db<-function(k,x){4*x^3-k*1/(x-1)}
k<-0.001
x<-2 #ponto inicial
tol<-0.0000001
passos<-0
while(abs(db(k,x))>tol){
  passos<-passos+1
  x<-x-0.0001*db(k,x)
  #cat(x,dmfunc(k,x),"\n")
  if (passos>10000) break
}
cat("iterações:",passos,"x=",x,"norma do gradiente=",abs(db(k,x)),"\n")
## iterações: 968 x= 1.00024981273861 norma do gradiente= 7.35234451099132e-08
## Implementação com método de Newton
## partindo de x=1.5 com k=0.001
##
b<-function(k,x){x^4-k*log(x-1)}
db<-function(k,x){4*x^3-k*1/(x-1)}
d2b<-function(k,x){12*x^2+k*1/(x-1)^2}
k<-0.001
x<-1.5 #ponto inicial
tol<-0.0000001
passos<-0
while(abs(db(k,x))>tol){
  passos<-passos+1
  x<-x-db(k,x)/d2b(k,x)
  #cat(x,dmfunc(k,x),"\n")
  if (passos>10000) break
}
cat("iterações:",passos,"x=",x,"norma do gradiente=",abs(db(k,x)),"\n")
## iterações: 6 x= 1.00024981273402 norma do gradiente= 1.61204383175573e-12

É interessante observar, na implementação com o Método de Newton, que iniciando no ponto \(x=2\) o método apresenta problemas de convergência (algo que não é incomum nesse método). É garantido, contudo, que nas proximidades da solução, quando a função é convexa, num problema de minimização, o método convergirá rapidamente, algo que aconteceu nesse caso em apenas 6 iterações, iniciando em \(x=1.5\). Em bons programas de otimização, as iterações iniciais são realizadas por um algoritmo mais seguro (eventualmente mais lento), passando-se a utilizar num segundo momento um algoritmo mais rápido (ex. Método de Newton) nas proximidades da solução.

A seguir vai uma possível implementação que usa o método do gradiente descendente até um certo valor da norma do gradiente (valor da derivada nesse caso), definido por tol1, e a partir desse valor passa a usar o Método de Newton.

## Implementação com método Combinado
## partindo de x=2 com k=0.00015
##
b<-function(k,x){x^4-k*log(x-1)}
db<-function(k,x){4*x^3-k*1/(x-1)}
d2b<-function(k,x){12*x^2+k*1/(x-1)^2}
k<-0.00015
x<-3 #ponto inicial
tol1<-0.01
tol2<-0.0000001
passos<-0
while(abs(db(k,x))>tol2){
  passos<-passos+1
  if (db(k,x)>tol1) x<-x-0.0001*db(k,x) #Gradiente Descendente
    else x<-x-db(k,x)/d2b(k,x) # Newton
  #cat(x,dmfunc(k,x),"\n")
  if (passos>10000) break
}
cat("iterações:",passos,"x=",x,"norma do gradiente=",abs(db(k,x)),"\n")
## iterações: 1115 x= 1.00003749578204 norma do gradiente= 1.43887568526679e-10

Se reduzirmos muito o valor de \(k\), para ganhar mais precisão no resultado, passamos a ter problemas numéricos que dificultam a convergência. Esses problemas podem ser solucionados usando operações envolvendo um maior numero de dígitos, usando o package Rmpfr, que usamos no problema envolvendo a estimativa de \(\pi\) pelo método geométrico para aumentar o número de digitos nas operações (veja aqui).

Solução usando funções do R para otimização

No R há várias funções e packages para otimização, usando métodos diversos. Um deles, para o caso unidimensional é a função optimize que é exemplificada a seguir com tolerância e valor de \(k\) idêntico aos da última estimativa usando o método combinado, pesquisando-se a solução no intervalo [0,3]. O resultado obtido foi o mesmo anterior.

k<-0.00015
optimize(b,c(0,3),tol=0.0000001,k=k)
## $minimum
## [1] 1.00003748996467
## 
## $objective
## [1] 1.00167868388232
minx<-1.00003748996467
cat("x=",minx,"norma do gradiente=",abs(db(k,minx)),"\n")
## x= 1.00003748996467 norma do gradiente= 0.00062082529103602

Note que a norma do gradiente (valor absoluto da derivada, no caso) no ponto mínimo estimado foi superior à obtida anteriormente.

Solução por simulação Monte Carlo

Finalmente, uma opção sempre possível, inclusive em mais dimensões, é a utilização de um método robusto (mas mais demandante de recursos computacionais) usando, por exemplo, simulação Monte Carlo (aplicado sem maiores refinamentos), como o implementado a seguir. Note contudo, que o tempo necessário para uma solução que atinja os mesmos níveis de tolerância anteriores exige um elevado número de simulações.

## Implementação com método Monte Carlo
## partindo com k=0.00015
##
b<-function(k,x){x^4-k*log(x-1)}
k<-0.00015
tol<-0.0000001
passos<-0
mindb<-5
dbcalculado<-mindb
pmt<-proc.time() # iniciando cronometragem do tempo
while(dbcalculado>tol){
  passos<-passos+1
  x<-runif(1,1,3)
  dbcalculado<-abs(db(k,x))
  if(dbcalculado<mindb) { # salva o menor valor obtido
    mindb=dbcalculado
    minx<-x
  }
  if (passos>20000000) break
}
cat("iterações:",passos,"x=",minx,"norma do gradiente=",abs(db(k,minx)),"\n")
## iterações: 20000001 x= 1.00003747036681 norma do gradiente= 0.00271371255518993
proc.time()-pmt #mostrando o tempo de processamento em seg
##    user  system elapsed 
##  476.84    0.25  488.80

Também nesse caso a norma do gradiente (valor absoluto da derivada, no caso) no ponto máximo estimado foi um pouco superior à obtida no método combinado e relativamente similar à obtida pelo pela função optimize. A solução demorou quase 4 min (o tempo é medido em segundos no resultado apresentado) para ser obtida, considerando 20 milhões de simulações, mostrando que esse é um problema importante desse método com relação aos anteriores, em que a solução é encontrada em segundos. Esse problema será amplificado em otimização envolvendo mais variáveis.