Gradient descent

exemplos


Henrique Aparecido Laureano [Lattes, LEG GitLab, GitLab]

Março de 2016

De forma intuitiva o algoritmo de gradiente descendente é um método para encontrar o mínimo de uma função de forma iterativa


Batch


\[ f(x, y) = x^{2} + 2 y^{2} \]

library(animation)
grad.desc()[1:3]
$par
            x             y 
-0.0675851986  0.0009735557 

$value
          x 
0.004569655 

$iter
[1] 36

\[ f(x, y) = \sin(0.5x^{2} - 0.25y^{2} + 3) \cos(2x + 1 - e^{y}) \]

Aumentando o tamanho do passo, argumento gamma

O default é 0.05

fn <- function(x, y) sin(.5*x**2 - .25*y**2 + 3) * cos(2*x + 1 - exp(y))
grad.desc(fn, init = c(-1, .5), gamma = .3, tol = 1e-04)[1:3]
Warning in grad.desc(fn, init = c(-1, 0.5), gamma = 0.3, tol = 1e-04):
Maximum number of iterations reached!
$par
         x          y 
2.26956698 0.07494521 

$value
        x 
0.1617829 

$iter
[1] 50

Referência: Demonstration of the Gradient Descent Algorithm


Construindo um exemplo

x0 <- rep(1, 5) ; x1 <- 1:5
x <- as.matrix(cbind(x0, x1))
 
y <- as.matrix(c(3, 7, 5, 11, 14))
m <- nrow(y)

x.scaled <- x
x.scaled[ , 2] <- (x[ , 2] - mean(x[ , 2])) / sd(x[ , 2])
# essa padronização aumenta a velocidade de convergência do algoritmo

solve(t(x) %*% x) %*% t(x) %*% y
   [,1]
x0  0.2
x1  2.6
lm(y ~ x[ , 2])$coefficients
(Intercept)      x[, 2] 
        0.2         2.6 
solve(t(x.scaled) %*% x.scaled) %*% t(x.scaled) %*% y
       [,1]
x0 8.000000
x1 4.110961
lm(y ~ x.scaled[ , 2])$coefficients
  (Intercept) x.scaled[, 2] 
     8.000000      4.110961 

Hipótese linear

\[ h_{\theta} (x) = \theta_{0} + \theta_{1} x \]

A ideia é minimizar a soma de quadrados dos erros, que no contexto de machine learning é representada pela função custo

\[ J(\theta_{0}, \theta_{1}) = \frac{1}{2m} \sum_{i=1}^{m} (h_{\theta} (x^{(i)}) - y^{(i)})^{2} \]

Para \(j = 1\) e \(j = 0\) o algoritmo do gradiente descendente é

Repetir até convergência:

\[ \theta_{j} := \theta_{j} - \alpha \frac{\partial}{\partial \theta_{j}} J(\theta_{0}, \theta_{1}) \]

Onde

\[ \frac{\partial}{\partial \theta_{j}} J(\theta) = \frac{1}{m} \sum_{i=1}^{m} (h_{\theta} (x^{(i)}) - y^{(i)}) x_{j}^{(i)} \]

Dessa forma o algoritmo do gradiente descendente fica

Repetir até convergência:

\[ \theta_{j} := \theta_{j} - \alpha \frac{1}{m} \sum_{i=1}^{m} (h_{\theta} (x^{(i)}) - y^{(i)}) x_{j}^{(i)} \]


Esse código não é eficiente, mas é intuitivo

# definindo a derivada da função custo em forma matricial
grad <- function(x, y, theta){
  gradient <- (1/m)*( t(x) %*% ( (x %*% t(theta)) - y ) )
  return(t(gradient))}
 
# definindo o algoritmo
grad.descent <- function(x, maxit){
  theta <- matrix(c(0, 0), nrow = 1)
  alpha = .05 # taxa de aprendizado padrão
  for(i in 1:maxit){
    theta <- theta - alpha*grad(x, y, theta)}
  return(theta)}

# objetivo: 0.2 e 2.6
grad.descent(x, 100)
            x0       x1
[1,] 0.4067232 2.542741
grad.descent(x, 500)
            x0      x1
[1,] 0.2069313 2.59808
grad.descent(x, 1000)
            x0       x1
[1,] 0.2000994 2.599972
# objetivo: 8 e 4.110961
grad.descent(x.scaled, 100)
           x0       x1
[1,] 7.952636 4.041608
grad.descent(x.scaled, 500) # convergência obtida mais rapidamente
     x0       x1
[1,]  8 4.110961

Custo e intuição de convergência

Tipicamente nós iteramos o algoritmo até que a mudança na função custo (como os valores de \(\beta_{0}\) e \(\beta_{1}\) atualizados) seja extremamente pequena, i.e., \(c\). \(c\) pode ser definido como o critério de convergência. Se \(c\) não é atingido após um dado número de iterações, você pode aumentar as iterações ou mudar a taxa de aprendizado \(\alpha\) para acelerar a convergência

beta <- grad.descent(x, 1000)
 
# definindo a hipótese linear
h <- function(xi, b0, b1) b0 + b1*xi
 
# definindo a função custo
cost <- t(mat.or.vec(1, m))
for(i in 1:m) cost[i, 1] <- (1/(2*m))*(h(x[i, 2], beta[1, 1], beta[1, 2]) - y[i, ])**2
 
(totalCost <- colSums(cost))
[1] 1.24
cost1000 <- totalCost
 
# aumentando o número de iterações para 1001
beta <- grad.descent(x, 1001)

cost <- t(mat.or.vec(1, m))
for(i in 1:m) cost[i, 1] <- (1/(2*m))*(h(x[i, 2], beta[1, 1], beta[1, 2]) - y[i, ])**2

cost1001 <- colSums(cost)
 
# essa diferença atende à seus critérios de convergência?
cost1000 - cost1001
[1] 1.515144e-11

Referência: Regression via Gradient Descent in R


Mais um exemplo

x <- runif(1000, -5, 5) ; y <- x + rnorm(1000) + 3
lm(y ~ x)$coefficients
(Intercept)           x 
  2.9653794   0.9966173 
plot(y ~ x, main = "Regressão linear") ; abline(lm(y ~ x), col = 2, lwd = 2)


Outra maneira de programar o algoritmo de gradiente descendente

# função custo
cost <- function(X, y, theta) sum((X %*% theta - y)**2) / (2*length(y))

# taxa de aprendizado e número limite de iterações
alpha <- .01 ; num_iters <- 1000

# armazenando o histórico
cost_history <- double(num_iters)
theta_history <- list(num_iters)

# inicializando os coeficientes
theta <- matrix(c(0, 0), nrow = 2)

# inserindo uma coluna de 1's para o intercepto
X <- cbind(1, matrix(x))

# gradiente descendente
for(i in 1:num_iters){
  error <- X %*% theta - y
  delta <- t(X) %*% error / length(y)
  theta <- theta - alpha * delta
  cost_history[i] <- cost(X, y, theta)
  theta_history[[i]] <- theta}

theta
          [,1]
[1,] 2.9652511
[2,] 0.9966184
# objetivo: 2.9653794 e 0.9966173

plot(y ~ x, main = "Regressão linear por gradiente descendente")
for(i in c(1, 3, 6, 10, 14, seq(20, num_iters, 10))) abline(theta_history[[i]], col = "blue")
abline(theta, col = 2, lwd = 3)


E como a função custo decai?

plot(cost_history, type = "l"
     , col = 2, lwd = 2, main = "Função custo", xlab = "Iterações", ylab = "Custo")


Referência: Linear regression by gradient descent


Stochastic


A diferença no algoritmo de gradiente descendente estocástico é que aqui é usada apenas uma observação aleatória da base em cada iteração

Quando a base de dados é muito grande esse algoritmo é útil, contudo, as estimativas obtidas com ele não são tão boas quanto a solução ótima global, já que apenas uma observação é usada por iteração

x <- runif(1000000, 1, 100) ; y <- 5 + 4*x
lm(y ~ x)$coefficients
(Intercept)           x 
          5           4 
# gds: gradiente descendente estocástico
gds <- function(x, y, b0, n, alpha){
  beta <- as.matrix(cbind(rep(b0[1], n), rep(b0[2], n)))
  for(i in 2:n){
    m <- length(x)
    sample_num <- sample.int(m, 1)
    xx <- x[sample_num]
    yy <- y[sample_num]
    beta[i, 1] <- beta[i-1, 1] - alpha*(beta[i-1, 1] + beta[i-1, 2]*xx - yy)
    beta[i, 2] <- beta[i-1, 2] - alpha*(xx*(beta[i-1, 1] + beta[i-1, 2]*xx - yy))}
  return(beta)}

b0 <- c(0, 0)

beta <- gds(x, y, b0, 1000000, .00005)

beta[1e06, ]
[1] 4.999974 4.000000
plot(beta[1:1000000, 1], type = "l"
     , col = 2, lwd = 2, xlab = "Iteração", ylab = expression(beta), main = "Coeficientes")
lines(beta[1:1000000, 2], col = 2, lwd = 2, lty = 2)
legend(85e04, 1, bty = "n", lty = c(1, 2), lwd = 2, col = 2
       , legend = c(expression(beta[0]), expression(beta[1])))


Referência: Stochastic gradient descent to find least square in linear regression


Boosting


Aqui \(h(x_{i})\) é atualizado com o resíduo de um modelo ajustado, onde esse resíduo é equivalente ao negativo do gradiente

É escolhido um chute inicial, \(h(x_{i})^{(0)}\), é calculado o \(- \partial J(y_{i}, h(x)^{(k)}) / \partial h(x_{i})^{(k)}\)

E é ajustado um modelo de regressão \(g(x_{i})^{(k)}\) fundamentado no negativo do gradiente

\[ h(x_{i})^{(k+1)} = h(x_{i})^{(k)} + \alpha g(x_{i})^{(k)}, \quad k = 0, 1, \cdots \]

até atingir convergência


library("mboost")

Exemplo: predição de gordura corporal

data("bodyfat", package = "TH.data")

lm(DEXfat ~ hipcirc + kneebreadth + anthro3a, bodyfat)$coefficients
(Intercept)     hipcirc kneebreadth    anthro3a 
-75.2347840   0.5115264   1.9019904   8.9096375 
model.boost <- glmboost(DEXfat ~ hipcirc + kneebreadth + anthro3a, bodyfat)
coefficients(model.boost, off2int = TRUE)
(Intercept)     hipcirc kneebreadth    anthro3a 
-75.2073365   0.5114861   1.9005386   8.9071301 
par(mfrow = c(1, 2), mar = c(5, 4, 4, 5))
plot(model.boost, off2int = TRUE, col = 2, lwd = 2, main = "Covariáveis")
plot(model.boost, col = 2, lwd = 2
     , ylim = range(coefficients(model.boost, which = model.boost$basemodel[[1]]$Xnames[-1]))
     , main = "Covariáveis") ; layout(1)

plot(AIC(model.boost), col = 2, lwd = 2, main = "AIC pelo número de iterações")


Referência: Model-based Boosting in R, A Hands-on Tutorial Using the R Package mboost