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] 50Referê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.6lm(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.110961lm(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.542741grad.descent(x, 500) x0 x1
[1,] 0.2069313 2.59808grad.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.041608grad.descent(x.scaled, 500) # convergência obtida mais rapidamente x0 x1
[1,] 8 4.110961Custo 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.24cost1000 <- 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-11Referê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.000000plot(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