7 ene. 2016

El gradiente descendente y la regresión lineal - Introducción

El objetivo de esta entrada es introducir al lector en algunos aspectos teóricos el método numérico del gradiente descendente, la regresión lineal y presentar una breve implementación en Python.
Para comenzar hablaremos del algoritmo del gradiente descendente (AGD) y luego mostraremos un ejemplo con el método de la regresión lineal. El AGD busca los puntos $p\in \Omega$ donde funciones del tipo $f:\Omega\subseteq \mathbb{R}^{m}\to \mathbb{R}$ alcanzan su mínimo. La idea de este método surge de la siguiente situación: Si $f$ es una función diferenciable en todo su dominio $\Omega$, entonces la derivada de $f$ en un punto $p\in \Omega$ en dirección de un vector unitario $v\in \mathbb{R}^{m}$ se define como: \begin{equation}\label{eq:01}df_{p}(v)=\nabla f(p)\cdot v.\end{equation} Observe que la magnitud de la ecuación $\eqref{eq:01}$ es: $$|df_{p}(v)|=||\nabla f(p)||||v||\cos \theta = ||\nabla f(p)||\cos \theta$$ Dicha magnitud es máxima cuando $\theta = 2n \pi, n\in \mathbb{Z}$. Es decir, para que $|df_{p}(v)|$ sea máxima los vectores $\nabla f(p)$ y $v$ deben ser paralelos. De esta manera la función $f$ crece más rápidamente en la dirección del vector $\nabla f(p)$ y decrece más rápidamente en la dirección del vector $-\nabla f(p)$. Nótese además que cualquier dirección $v\in \mathbb{R}^m$ ortogonal a $\nabla f(p)$ es una dirección de cambio nulo. Esto sugiere que la dirección negativa del gradiente $-\nabla f(p)$ es una buena dirección de búsqueda para encontrar el minimizador de la función $f$. Es importante recordar que dada una hipersuperfice de nivel $S_{k}=\{p\in \Omega: f(p)=k\}$, si $p\in S_{k}$ entonces $\nabla f(p)\perp t_{p}$, donde $t_{p}$ es cualquier vector tangente a $S_{k}$ en $p$. Por lo tanto surge de manera natural el siguiente algoritmo:

Algoritmo del gradiente descendente

Sea $f:\Omega \subseteq \mathbb{R}^{m}\to \mathbb{R}$, si $f$ tiene un mínimo en $p$, para encontrar a $p$ se construye una sucesión de puntos $\{p_t\}$ tal que $p_{t}$ converge a $p$. Para resolver esto, comenzamos en $p_{t}$ y nos desplazamos por una cantidad $-\lambda_{t}\nabla f(p_{t})$ para encontrar el punto $p_{t+1}$ más cercano a $p$. Es decir: \begin{equation}\label{eq:02} p_{t+1}=p_{t}-\lambda_{t}\nabla f(p_{t})\end{equation} donde el parámetro $\lambda_{t}$ se selecciona de tal manera que $p_{t+1}\in \Omega$ y $f(p_{t})\geq f(p_{t+1})$.
El parámetro $\lambda_{t}$ se selecciona para maximizar la cantidad a la que decrece la función $f$ en cada paso. Es decir, $\lambda_{t}$ se escoge para minimizar la aplicación: $$\phi_{t}(\lambda)=f(p_{t}-\lambda \nabla f(p_{t}))$$ Esto es: $$\lambda_{t}=\arg \min_{\lambda \geq 0} f(p_{t}-\lambda \nabla f(p_{t})). $$
Al momento de realizar una implementación computacional de este método, es necesario establecer un criterio de parada para cada iteración del algoritmo. Para esto, el lector puede definir un umbral $\epsilon>0$ y detener el algoritmo cuando los errores relativos en los valores de $f(p)$ cumplan: \begin{equation}\label{eq:03}\frac{|f(p_{t+1})-f(p_{t})|}{\max\{1, |f(p_{t})|\}}<\epsilon.\end{equation}
El criterio de parada se selecciona de acuerdo a las necesidades y creencias del lector.

Ejemplo del AGD aplicado al método de regresión lineal

La técnica de la regresión lineal se clasifica en la categoría del aprendizaje mecánico supervisado y dentro de esta es de tipo aprendizaje mecánico por regresión. El objetivo de la regresión lineal es ajustar una aplicación afín a un conjunto de datos $M=\{(x_{i},y_{i})\}_{1\leq i\leq k}$ con $(x_{i},y_{i})\in \mathbb{R}^m\times\mathbb{R}^n$. Es decir, se busca una aplicación $h_{\theta}:\mathbb{R}^{m}\to \mathbb{R}^{n}$ de la forma: \begin{equation}\label{eq:04}h_{\theta}(v)=v\cdot\beta + \alpha\end{equation} donde $v\in \mathbb{R}^{m}$, $\alpha\in \mathbb{R}^{n}$ son vectores filas, $\beta\in M_{m\times n}(\mathbb{R})$ es una matriz de orden $m\times n$ y $\theta^{\top}=(\alpha^\top, \beta^\top)\in M_{n\times(m+1)}(\mathbb{R})$, de tal forma que la suma de las distancias entre los puntos $h_{\theta}(x_{i})$ e $y_{i}$ sea mínima. Se desea encontrar un parámetro $\theta$ tal que: \begin{equation}\label{eq:05}\theta = \arg\min_{\theta}\frac{1}{2k}\sum_{i=1}^{k}||y_{i}-h_{\theta}(x_{i})||^{2} \end{equation}
Por lo tanto para encontrar los parámetros $\alpha$ y $\beta$ (de tal forma que se cumpla la condición $\eqref{eq:05}$), utilizamos el AGD empleando la ecuación $\eqref{eq:02}$. Para cuando $m=n=1$ las iteraciones que se deben realizar son: $$\alpha_{t+1}=\alpha_{t}+\frac{\lambda_{t}}{k}\sum_{i=1}^{k}(y_{i}-h_{\theta_{t}}(x_{i})),$$ $$\beta_{t+1}=\beta_{t}+\frac{\lambda_{t}}{k}\sum_{i=1}^{k}(y_{i}-h_{\theta_{t}}(x_{i}))x_{i}.$$ Implementando lo anterior en Python se tiene algo así:
#! /usr/bin/env python
# -*- coding: UTF-8 -*-

from __future__ import division, print_function
import numpy as np
import random
import sklearn
from sklearn.datasets.samples_generator import make_regression
import pylab
from scipy import stats

def gradientDescent(Lambda, x, y, epsilon, max_iter = 10000):
    """
    Esta función permite calcular el parametro «theta» para la aplicación afín que mejor se ajusta
    a los datos (x, y).
    :param Lambda: Tasa de aprendizaje.
    :param x: Valores independientes de entrada.
    :param y: Valores dependientes de salida.
    :param epsilon: Umbral de parada.
    :param max_iter: Numero de veces que se repite la iteración. 
    """
    
    converged = False
    iter = 0
    k = x.shape[0] # tamaño de la muestra.
    
    # Iniciamos el parametro theta.
    alpha = np.random.random(x.shape[1])
    beta = np.random.random(x.shape[1])
    
    # Error total, e(theta_t)
    temp_et = sum([(y[i]-alpha - beta*x[i])**2 for i in range(k)])
    
    #Ciclo de iteraciones
    while not converged:
        #Para el conjunto muestral se calcula el gradiente de la funcion e(theta)
        grad0 = 1.0/k * sum([(y[i]-alpha - beta*x[i]) for i in range(k)])
        grad1 = 1.0/k * sum([(y[i]-alpha - beta*x[i])*x[i] for i in range(k)])
    
        #Guardamos temporamente el parametro theta
        temp_alpha = alpha + Lambda*grad0
        temp_beta = beta + Lambda*grad1
        
        #Actualizamos el parametro theta
        alpha = temp_alpha
        beta = temp_beta
        
        # Error cuadrático medio.
        et = sum([(y[i]-alpha - beta*x[i])**2 for i in range(k)])
        
        if abs(temp_et-et) <= epsilon:
            print('La iteración converge: ', iter)
            converged = True
            
        temp_et = et #Actualización del error
        iter += 1 #Actualización del numero de iteraciones
        
        if iter == max_iter:
            print('¡Máximo de iteraciones excedido!')
            converged = True
            
    return alpha, beta       
    pass

if __name__ == '__main__':
    x, y = make_regression(n_samples = 200, n_features=1, n_informative=1, random_state=0, noise=35)
    print('Tamaño del conjunto de prueba: ', x.shape, y.shape)
    
    Lambda = 0.001 # Tasa de aprendizaje
    epsilon = 0.01 # Umbra de convergencia
    
    # Llamamos la función gradientDescent para obtener el parametro theta.
    alpha, beta = gradientDescent(Lambda, x, y, epsilon, max_iter=10000)
    print('alpha = ', alpha, ', beta = ', beta)
    
    #Comprobando los resultados de nuestro algoritmo con scipy linear regression
    slope, intercept, r_value, p_value, slope_std_error = stats.linregress(x[:,0], y)
    print('Intercepto = ', intercept, ', Pendiente = ', slope) 
 
    # plot
    for i in range(x.shape[0]):
        y_predict = alpha + beta*x 
 
    pylab.plot(x, y, 'o')
    pylab.plot(x, y_predict, 'k-')
    pylab.show()
    print('¡Listo!')
El lector puede notar que en la líneas 53-55 se estableció el criterio de parada $|e(\theta_t)-e(\theta_{t+1})|<\epsilon$ como $\verb|abs(temp_et-et) <= epsilon|$ donde $\verb|temp_et|=e(\theta_t)$ y $\verb|et|=e(\theta_{t+1})$. Esta condición es equivalente a la condición de parada dada por los errores relativos en la ecuación $\eqref{eq:03}$.
Para el caso más general de $m$ y $n$, conviene escribir la ecuación $\eqref{eq:04}$ como: \begin{equation} h_{\theta}(v)=\xi\cdot \theta \end{equation} donde $\xi=(1, v)$. Definimos además las matrices $X^\top = (\xi_{1}^\top,\dots\xi_{k}^\top)$ e $Y^\top = (y_{1}^\top,\dots y_{k}^\top)$. Esto nos permite escribir la condición $\eqref{eq:05}$ de forma equivalente como: \begin{equation} \theta = \arg \min_{\theta} \frac{1}{2m}||X\theta -Y||^2. \end{equation} Donde la función $E(\theta)=||X\theta -Y||^2$ es conocida como la función de costo o pérdida y la función $h_{\theta}(x)=\xi\cdot\theta$ conocida como la función hipótesis. Por lo tanto se debe aplicar el AGD a la función $E$. Para esto tenemos que el gradiente de $E$ es: \begin{eqnarray*} \nabla_{\theta}E(\theta)&=&\frac{1}{2m}\nabla_{\theta}[(X\theta -Y)^\top (X\theta -Y)]\\ &=&-\frac{1}{2m}\nabla_{\theta}(Y^\top X\theta)-\frac{1}{2m}\nabla_{\theta}(\theta^\top X^\top Y)+\frac{1}{2m}\nabla_{\theta}(\theta^\top X^\top X\theta)\\ &=& \frac{1}{m} X^\top(X\theta-Y) \end{eqnarray*} Luego se debe realizar la siguiente iteración: \begin{equation} \theta_{t+1}=\theta_{t}-\frac{\lambda_{t}}{m}X^\top(X\theta_{t}-Y) \end{equation} Implementando esto en Python sería algo así:
#! /usr/bin/env python
# -*- coding: UTF-8 -*-

from __future__ import division, print_function
import numpy as np
import random
import sklearn
from sklearn.datasets.samples_generator import make_regression
import pylab
from scipy import stats

def gradientDescent2(Lambda, x, y, max_iter = 10000):
    """
    Esta función permite calcular el parametro «theta» para la aplicación afín que mejor se ajusta
    a los datos (x, y).
    :param Lambda: Tasa de aprendizaje.
    :param x: Valores independientes de entrada.
    :param y: Valores dependientes de salida.
    :param epsilon: Umbral de parada.
    :param
    """
    k = x.shape[0] # número de muestras.
    theta = np.ones(x.shape[1]+1)    
    X = np.c_[ np.ones(x.shape[0]), x]
    Y = y
    for iter in range(0, max_iter):
        hypothesis = np.dot(X, theta)
        loss = hypothesis - Y
        gradient = np.dot(X.T, loss) / k         
        theta = theta - alpha * gradient  # update
    return theta   
    pass

if __name__ == '__main__':
    x, y = make_regression(n_samples = 200, n_features = 1, n_informative = 1, random_state = 0, noise = 35) 
    m, n = np.shape(x)
    alpha = 0.001 # learning rate
    theta = gradientDescent2(alpha, x, y, 10000)
    # plot
    for i in range(x.shape[1]):
        y_predict = theta[0] + theta[1]*x 
    pylab.plot(x, y, 'o')
    pylab.plot(x, y_predict,'k-')
    pylab.show()
    print(theta)
    print('¡Listo!')
Este último algoritmo es más eficiente que el anterior, además el lector, si así lo desea, puede modificarlo para establecer su propio criterio de parada. Para ir finalizando, merece mencionar que el problema de la regresión lineal se puede resolver de una manera más eficiente que los dos métodos anteriores, siempre y cuando en la ecuación $\nabla_{\theta}E(\theta)=\frac{1}{m} X^\top(X\theta-Y)$ la matriz $ X^\top X$ sea invertible, podemos optimizar $E$ empleando el criterio de los puntos críticos de la primera y segunda derivada. En este caso nosotros encontramos que el parámetro $\theta$ que deseamos es: \begin{equation} \theta = (X^\top X)^{-1}X^\top Y \end{equation} Esta última expresión se puede implementar fácilmente en Python así que se sugiere como ejercicio para el lector. Finalmente, espero compartir pronto con ustedes otra agradable entrada de aplicaciones de las matemáticas y la programación.

Referencias