Gradiente Descendente

Para estimar parâmetros de uma RNA com um único neurônio

Author

Jessica Kubrusly

Suponha conhecidos a matriz \(X\) de dimensão \(n \times m\), com os valores observados para as covariáveis, e o vetor \(Y\) com os valores observados para a variável alvo. O Método do Gradiente Descendentes busca minimizar uma função de custo \(J\), que mede o erro no modelo, a partir de um método iterativo. O método iterativo do gradiente descendente para encontrar o ponto de mínimo da função de custo \(J\) é definido por:

  1. Atribua valores para o hiperparâmetro \(\delta\), \(\varepsilon\) e \(K\).

  2. Escolha um ponto inicial qualquer \((\Theta_0,\mathbf{w}_0)\) e faça \(k=0\).

  3. Calcule \(\nabla J(\Theta_k,\mathbf{w}_k)\).

  4. Atualize o ponto \((\Theta_{k+1},\mathbf{w}_{k+1}) = (\Theta_k,\mathbf{w}_k) - \delta \nabla J(\Theta_k,\mathbf{w}_k)\).

  5. Verifique os critérios de parada:

    • \(||\nabla J(\Theta_k,\mathbf{w}_k)|| < \varepsilon\)

    • \(||(\Theta_{k+1},\mathbf{w}_{k+1}) - (\Theta_k,\mathbf{w}_k)|| < \varepsilon\)

    • \(|J(\Theta_{k+1},\mathbf{w}_{k+1}) - J(\Theta_k,\mathbf{w}_k)| < \varepsilon\)

    • \(k > K\)

      Se os critérios de parada não foram satisfeitos, faça \(k = k + 1\) e volte para o Passo 2.

Para realizar o método iterativo precisamos do gradiente de \(J\). Vejamos como é o vetor gradiente para \(J\) definida acima. \[ \nabla J(\mathbf{w},\Theta)= \left( \frac{\partial J}{\partial w_1}(\mathbf{w},\Theta), \frac{\partial J}{\partial w_2}(\mathbf{w},\Theta), \ldots, \frac{\partial J}{\partial w_m}(\mathbf{w},\Theta), \frac{\partial J}{\partial \Theta}(\mathbf{w},\Theta) \right) \]

Problemas de regressão

Para o caso de um problema de regressão, uma função de custo adequada é a soma dos erros ao quadrado (SSE). Nesse caso, a função \(J\) é definida a seguir
\[ J(\mathbf{w},\Theta)= \sum_{i=1}^n \left( y_i - g({w_1x_{ì,1}+w_2x_{i,2} + \ldots + w_m x_{i,m}+\Theta }) \right)^2 \] sendo \(g\) a função de ativação. No R ela pode ser definida assim:

J = function(theta,w){
  soma = 0
  n = nrow(X)
  for(i in 1:n){
    linha_x = X[i,]
    yi = Y[i]
    soma = soma + (yi - g(sum(w*linha_x) + theta))^2
  }
  return(soma)
}

O vetor gradiente de \(J\) nesse caso é definido pelas coordenadas: \[ \frac{\partial J}{\partial w_j}(\mathbf{w},\Theta) = \sum_{i=1}^n -2\left( y_i - g({w_1x_{ì,1}+w_2x_{i,2} + \ldots + w_m x_{i,m}+\Theta}) \right) g'({w_1x_{ì,1}+w_2x_{i,2} + \ldots + w_m x_{i,m}+\Theta}) x_{i,j} \] \[ \frac{\partial J}{\partial \Theta}(\mathbf{w},\Theta) = \sum_{i=1}^n -2\left( y_i - g({w_1x_{ì,1}+w_2x_{i,2} + \ldots + w_m x_{i,m}+\Theta}) \right) g'({w_1x_{ì,1}+w_2x_{i,2} + \ldots + w_m x_{i,m}+\Theta}) \]

E no R ele pode ser implementado por:

gradiente_J = function(theta,w){
  
  n = nrow(X)
  m = ncol(X)
  grad = NULL
  #derivada em relacao a theta
  soma = 0
  for(i in 1:n){
      linha_x = X[i,]
      yi = Y[i]
      soma = soma - 2*(yi - g(sum(w*linha_x) + theta))*dg(sum(w*linha_x) + theta)
  }
  grad[1] = soma
  #derivadas em relacao a wj
  for(j in 1:m){
    soma = 0
    for(i in 1:n){
      linha_x = X[i,]
      yi = Y[i]
      soma = soma - 2*(yi - g(sum(w*linha_x) + theta))*dg(sum(w*linha_x) + theta)*X[i,j]
    }
    grad[j+1] = soma
  }
  return(grad)
}

Função de Ativação Logística

A função de ativação logística e sua derivada são definidas por:

\[ g(x) = \dfrac{1}{1+e^{-x}} \qquad \text{ e } \qquad g'(x) = \dfrac{e^{-x}}{(1+e^{-x})^2} \]

No R podemos fazer uma função que recebe um número \(x\) e retorna \(g(x)\):

g = function(x){
  1/(1+exp(-x))
}

E também a função uma função que recebe \(x\) e retorna \(g(x)\):

dg = function(x){
  exp(-x)/(1+exp(-x))^2
}

Aqui vamos supor que a matriz de dados \(X\) e o vetor \(Y\) estão definidos globalmente, isto é, temos acesso a eles dentro da função. Agora que já temos como avaliar \(\nabla J\), podemos implementar o processo iterativa para buscar os valore para os parâmetros que minimizam \(J\).

gradiente_descendente = function(theta0 = rnorm(1), w0 = rnorm(ncol(X)), k=1, d=0.01, e=0.01, K=1000){
  
  gradJ = gradiente_J(theta0,w0)
  novo = c(theta0,w0) - d*gradJ
  theta1 = novo[1]
  w1 = novo[-1]
  if(sqrt(sum(gradJ^2))<e &
     sqrt(sum((novo - c(theta0,w0))^2))<e &
     sqrt(sum((J(theta0,w0) - J(theta1,w1))^2))<e){
    print(paste("Convergiu em k =",k))
    return(novo)
  }
  if(k > K){
    print(paste("Não convergiu com k =",k))
    return(novo)
  }
  
  return(gradiente_descendente(theta0 = theta1, w0 = w1, k = k+1,d,e,K))
}

Vejamos os parâmetros estimados para a base de dados de seguros trabalhada nas aulas anteriores. Depois vamos comparar com os valores dos parâmetros fornecidos pela função neuralnet.

Primeiro, carregamos os pacotes necessários.

library(tidyverse)
library(neuralnet)

Depois vamos carregar a matriz de treino com os valores de \(X\) e \(Y\). Aqui vamos usar já a matriz pré-processada, isto é, não há dados faltantes e as variáveis já foram padronizadas.

base = readRDS("base_treino_final_log.RDS")
Y = base[,"charges"]
X = base[,c("age","sexmale","bmi","children","smokeryes","regionnorthwest","regionsoutheast","regionsouthwest","marriedyes")]
pesos_gd = gradiente_descendente()
[1] "Convergiu em k = 298"
theta_gd = pesos_gd[1]
w_gd = pesos_gd[-1]

Vamos agora encontrar os pesos a partir do pacote neuralnet para comparação.

modelo_completo_log = neuralnet(
  formula = charges ~ age + sexmale + bmi + children + smokeryes + regionnorthwest + regionsoutheast +  regionsouthwest + marriedyes,
  data = base,
  hidden = 0,
  linear.output = FALSE)
pesos_nn = modelo_completo_log$weights[[1]][[1]][,1]
theta_nn = pesos_nn[1]
w_nn = pesos_nn[-1]

Será que chegamos em valores parecido?

J(theta_gd,w_gd);J(theta_nn,w_nn)
       1 
4.996842 
      1 
4.99717 
theta_gd;theta_nn
[1] -1.997767
[1] -1.988397
w_gd;w_nn
[1]  0.36263366 -0.02204345  0.33392079  0.05798238  2.08896674  0.01212005
[7] -0.03859576 -0.05691824 -0.01678537
[1]  0.362309538 -0.022591548  0.334424944  0.058045852  2.087482595
[6]  0.003040218 -0.045943508 -0.064378338 -0.020544182

Função de Ativação Tangente Hiperbólica

A função de ativação tangente hiperbólica e sua derivada são definidas por:

\[ g(x) = \dfrac{e^{2x}-1}{e^{2x}+1} \qquad \text{ e } \qquad g'(x) = \dfrac{4e^{2x}}{(e^{2x}+1)^2} \]

No R as funções g e dg devem ser atualizadas:

g = function(x){
  (exp(2*x)-1)/(exp(2*x)+1)
}
dg = function(x){
  (4*exp(2*x))/((exp(2*x)+1)^2)
}

Já a função que retorna o gradiente de \(J\) não precisa ser atualizada. Vamos agora testar na base de seguros. Mas se a função de ativação é a tangente hiperbólica será preciso usar a base onde a variável \(Y\) foi transformada para valores entre -1 e 1.

base = readRDS("base_treino_final_tan.RDS")
Y = base[,"charges"]
X = base[,c("age","sexmale","bmi","children","smokeryes","regionnorthwest","regionsoutheast","regionsouthwest","marriedyes")]

Agora vamos calcular o valor dos pesos que minimizam \(J\), usando uma taxa de aprendizado menor (usei \(d=0,01\) e o método não convergiu).

pesos_gd = gradiente_descendente(d = 0.001)
[1] "Convergiu em k = 442"
theta_gd = pesos_gd[1]
w_gd = pesos_gd[-1]

Em seguida calculamos os pesos pelo pacote neuralnet

modelo_completo_tan = neuralnet(
  formula = charges ~ age + sexmale + bmi + children + smokeryes + regionnorthwest + regionsoutheast +  regionsouthwest + marriedyes,
  data = base,
  hidden = 0,
  act.fct = "tanh",
  linear.output = FALSE)
pesos_nn = modelo_completo_tan$weights[[1]][[1]][,1]
theta_nn = pesos_nn[1]
w_nn = pesos_nn[-1]

E por fim comparamos os resultados.

J(theta_gd,w_gd);J(theta_nn,w_nn)
       1 
19.98717 
       1 
19.98717 
theta_gd;theta_nn
[1] -1.001901
[1] -1.001106
w_gd;w_nn
[1]  0.181421058 -0.010559870  0.166839902  0.028959058  1.044957881
[6]  0.008905517 -0.016770068 -0.025758976 -0.007696399
[1]  0.181376698 -0.010688033  0.166852674  0.028962066  1.044763912
[6]  0.008340786 -0.017295544 -0.026316809 -0.007956693

Problemas de classificação

Para o caso de um problema de classificação, a função de custo adequada é a entropia cruzada (CE). Nesse caso, a função \(J\) é definida a seguir.

\[ J(\mathbf{w},\Theta)= -\dfrac{1}{n} \sum_{i=1}^n y_i \ln(g({w_1x_{ì,1}+w_2x_{i,2} + \ldots + w_m x_{i,m}+\Theta })) + (1-y_i)\ln(1-g({w_1x_{ì,1}+w_2x_{i,2} + \ldots + w_m x_{i,m}+\Theta })) \]

sendo \(g\) a função logística. N o R essa função é implementada por:

clip = function(p,eps=0.001){
  pmax(pmin(p,1-eps),eps)
}
J = function(theta,w){
  soma = 0
  n = nrow(X)
  for(i in 1:n){
    linha_x = X[i,]
    yi = Y[i]
    p = g(sum(w*linha_x) + theta)
    p = clip(p)
    soma = soma + yi*log(p) + (1-yi)*log(1 - p)
  }
    soma = -soma/n
  return(soma)
}

A função clip evita \(p=0\) ou \(p=1\), logo, evita \(\ln(0)\).

Quando combinamos a função de custo CE com a função de ativação logistica as derivadas parciais de \(J\) simplificam e ficam:

\[ \frac{\partial J}{\partial w_j}(\mathbf{w},\Theta) = -\dfrac{1}{n}\sum_{i=1}^n \left( y_i - g({w_1x_{i,1}+w_2x_{i,2} + \ldots + w_m x_{i,m}+\Theta})\right) x_{i,j} \]

\[ \frac{\partial J}{\partial \Theta}(\mathbf{w},\Theta) = -\dfrac{1}{n}\sum_{i=1}^n \left( y_i - g({w_1x_{i,1}+w_2x_{i,2} + \ldots + w_m x_{i,m}+\Theta})\right) \]

No R o gradiente passa a ser calculado da seguinte forma:

gradiente_J = function(theta,w){
  n = nrow(X)
  m = ncol(X)
  grad = NULL
  #derivada em relacao a theta
  soma = 0
  for(i in 1:n){
    linha_x = X[i,]
    yi = Y[i]
    p = g(sum(w*linha_x) + theta)
    p = clip(p)
    soma = soma + (yi - p)
  }
  grad[1] = -soma/n
  #derivadas em relacao a wj
  for(j in 1:m){
    soma = 0
    for(i in 1:n){
      linha_x = X[i,]
      yi = Y[i]
      p = g(sum(w*linha_x) + theta)
      p = clip(p)
      soma = soma + (yi - p)*X[i,j]
    }
    grad[j+1] = -soma/n
  }
  return(grad)
}

Para testarmos os resultados vamos precisar da base onde a variável alvo é categórica: cliente de alto custo e cliente normal.

base = readRDS("base_treino_final_log_cla.RDS")
colnames(base)
 [1] "age"             "sexmale"         "bmi"             "children"       
 [5] "smokeryes"       "regionnorthwest" "regionsoutheast" "regionsouthwest"
 [9] "custonormal"     "marriedyes"     
Y = base[,"custonormal"]
table(Y)
Y
  0   1 
123 432 
X = base[,c("age","sexmale","bmi","children","smokeryes","regionnorthwest","regionsoutheast","regionsouthwest","marriedyes")]
colnames(X)
[1] "age"             "sexmale"         "bmi"             "children"       
[5] "smokeryes"       "regionnorthwest" "regionsoutheast" "regionsouthwest"
[9] "marriedyes"     
dim(X)
[1] 555   9

Vamos chamar o método com uma taxa de aprendizado agora maior, caso contrário a convergência demora muito para ocorrer.

pesos_gd = gradiente_descendente(d = 0.1, e=0.1)
[1] "Convergiu em k = 153"
theta_gd = pesos_gd[1]
w_gd = pesos_gd[-1]

Vamos agora encontrar os pesos a partir do pacote neuralnet para comparação.

modelo_completo_log_cla = neuralnet(
  formula = custonormal ~ age + sexmale + bmi + children + smokeryes + regionnorthwest + regionsoutheast +  regionsouthwest + marriedyes, 
  data = base, 
  hidden = 0,
  linear.output = FALSE,
  err.fct = "ce")
pesos_nn = modelo_completo_log_cla$weights[[1]][[1]][,1]
theta_nn = pesos_nn[1]
w_nn = pesos_nn[-1]

Será que chegamos em valores parecido?

J(theta_gd,w_gd);J(theta_nn,w_nn)
        1 
0.4511685 
        1 
0.5066871 
theta_gd;theta_nn
[1] 1.866671
[1] 2.822368
w_gd;w_nn
[1] -0.14769352 -0.87046305 -0.13183466  0.12259710 -0.57173540 -0.28481822
[7] -0.30981415 -0.03060572  0.40839782
[1] -0.72076646 -0.09985034 -0.50818692  0.01321182 -4.50361647  0.19484267
[7] -0.19231708  0.21353938  0.16311746