1 Introducción

En esta publicación de Rpubs se muestra como implementar el algoritmo del gradiente descendente en R.

El ejemplo está inspirado en el video que se muestra abajo. Por favor vea el video antes de seguir leyendo este contenido.

2 Función de interés

La función de interés mostrada en el video es \(f(x,y)\) y el objetivo es encontrar las coordenadas \(x\) e \(y\) que minimizan a \(f(x,y)\).

\[ f(x,y)=(x^2+y-11)^2+(x+y^2-7)^2 \]

Para dibujar la función podemos utilizar el siguiente código de R.

fun1 <- function(x, y) (x^2+y-11)^2 + (x+y^2-7)^2

a <- 5
x <- seq(from=-a, to=a, length.out=50)
y <- seq(from=-a, to=a, length.out=50)
z <- outer(x, y, fun1)

image(x=x, y=y, z=z)
contour(x=x, y=y, z=z, add=TRUE, nlevels=40, col=gray(0.5))

3 Gradiente descendente

Para iniciar el algoritmo se parte de un valor \(\boldsymbol{\theta}_0\) y luego se actualiza utilizando la siguiente expresión.

\[ \boldsymbol{\theta}_{i+1} = \boldsymbol{\theta}_i - \alpha \, \Delta f \quad \text{para} \quad i=0,1,2, \ldots, \] donde \(\alpha\) es la tasa de aprendizaje (learning rate) y \(\Delta\) representa el vector Score con las primeras derivadas de la función \(f\) de interés.

4 Implementación

Implementar el algoritmo es muy fácil, primero debemos definir el \(\theta_0\), la tasa de aprendizaje \(\alpha\) y un número fijo de iteraciones de \(n=9\)

theta_0 <- c(0, 0)
alfa <- 0.01
n <- 9

Ahora vamos a re-escribir la función de interés con un único parámetro vectorial, esto se necesita para poder obtener el gradiente (derivada) en un punto de forma computacional. Si usted quiere profundizar en el tema de cálculo de derivadas en forma computacional con R visite este enlace.

fun2 <- function(w) (w[1]^2+w[2]-11)^2 + (w[1]+w[2]^2-7)^2

Compare la función fun1 y fun2, son en realidad la misma función.

Luego vamos a usar la ecuación del algoritmo para actulizar los valores, el código se muestra a continuación.

library(numDeriv)                           # Package to obtain derivates

thetas <- matrix(NA, ncol=2, nrow=n+1)      # To store results
colnames(thetas) <- c('x', 'y')
thetas[1, ] <- theta_0                      # Initial value

for (i in 2:(n+1))
  thetas[i, ] <- thetas[i-1, ] - alfa * grad(func=fun2, x=thetas[i-1, ])

En el objeto thetas se encuentran almacenados los valores en cada una de las iteraciones, para explorar su contenido escribimos lo siguiente.

thetas
##               x         y
##  [1,] 0.0000000 0.0000000
##  [2,] 0.1400000 0.2200000
##  [3,] 0.3364902 0.4951501
##  [4,] 0.6047242 0.8301042
##  [5,] 0.9560018 1.2156580
##  [6,] 1.3865301 1.6151022
##  [7,] 1.8605039 1.9584806
##  [8,] 2.3018487 2.1722241
##  [9,] 2.6263942 2.2410364
## [10,] 2.8089320 2.2201119

Podemos superponer estos valores a la figura inicial y el resultado sería el mostrado a continuacion.

En la figura podemos observar la evolución del vector \(\boldsymbol{\theta}\), desde el valor de inicio hasta el óptimo.

5 Reto 1

  1. Modifique la tasa de aprendizaje \(\alpha\) para volver más lento el proceso de búsqueda.
  2. Busque un valor de \(\alpha\) que perjudique el proceso de búsqueda haciéndolo “enloquecer”.
  3. Vuelva a repetir todo el ejercicio pero sin usar un valor de \(n\). Haga que el algoritmo busque y busque hasta que la diferencia entre las componentes de dos iteraciones sucesivas sean ambas menores que 0.0001.

6 Reto 2

El objetivo es explicar la media de la variable \(Y\) en función de los valores de la variable \(X\) por medio del modelo \(Y \sim N(\beta_0 + \beta_1 X, \sigma^2=3)\). Suponga que se tiene como información disponible los vectores x e y que se muestran a continuación.

x <- c(1, 3, 4, 6, 9)
y <- c(-5, 6, 6, 17, 23)
  1. Construya la función de log-verimilitud \(f(\beta_0, \beta_1)\).
  2. Dibuje \(f(\beta_0, \beta_1)\).
  3. Obtenga los valores de \(\beta_0\) y \(\beta_1\) que maximizan \(f(\beta_0, \beta_1)\) usando el algoritmo gradiente descendiente.
  4. Haga una figura similar a la de la evolución del vector \(\boldsymbol{\theta}\) del ejemplo.