7 ene. 2016

El algoritmo del gradiente descendente

El objetivo en esta ocasión es estudiar el algoritmo del gradiente descendente (AGD) y presentar una breve aplicación al problema de mínimos cuadrados empleando el lenguaje de programación Python.

Para comenzar, debemos saber que AGD es un algoritmo iterativo de optimización que permite encontrar valores mínimos de funciones convexas y diferenciables en todo su dominio. Cuando se desea encontrar el mínimo de una función convexa a partir de un $p$ de su dominio, se toman pasos proporcionales al sentido opuesto del gradiente de la función en el punto $p$. Por otro lado, si se desea encontrar el máximo de una función convexa, entonces se toman pasos proporcionales al mismo sentido del gradiente y en este caso se denomina algoritmo de gradiente ascendente.

Figura 1. El gradiente descendente construye un sucesión de puntos sobre cada una de las curvas de nivel definidas por el gradiente de la función.

El AGD consiste en considerar una función del tipo $f:\Omega\subseteq \mathbb{R}^{m}\to \mathbb{R}$ convexa, diferenciable y que alcanzan su mínimo en $\Omega$, bajo estas condiciones 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:

Sea $f:\Omega \subseteq \mathbb{R}^{m}\to \mathbb{R}$ una función convexa de clase $C^{1}(\Omega)$ con 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 esto, se comienza en un punto $p_{t}$ y nos desplazamos por una cantidad $-\alpha_{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}-\alpha_{t}\nabla f(p_{t})\end{equation} donde el parámetro $\alpha_{t}$ se selecciona de tal manera que $p_{t+1}\in \Omega$ y $f(p_{t})\geq f(p_{t+1})$.

Que una función $f:\Omega \subseteq \mathbb{R}^{m}\to \mathbb{R}$ sea convenza quiere decir que satisface para todo $w, z\in \Omega$ y $\alpha \in [0, 1]$ la siguiente condición: \begin{equation} f(\alpha w + (1 - \alpha) z) \leq \alpha f(w)+(1-\alpha)f(z). \end{equation} Note que está condición es válida para todo $m\geq 1$.


Note en la figura a) que si la función no es convexa, el gradiente puede encontrar diferentes mínimos locales con diferentes puntos de inicio $p_{_0}$.

El parámetro $\alpha_{t}$ se selecciona para maximizar la cantidad a la que decrece la función $f$ en cada paso. Es decir, $\alpha_{t}$ se escoge para minimizar la aplicación: $$\phi_{t}(\alpha)=f(p_{t}-\alpha \nabla f(p_{t}))$$ Esto es: $$\alpha_{t}=\arg \min_{\alpha \geq 0} f(p_{t}-\alpha \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} En resumen el AGD se describe de la siguiente manera:

AGD
Inicialización: $p = p_{_0}$
while $||\nabla f(p)||_{2}>\epsilon$ do:
$$p \leftarrow p -\alpha\cdot \nabla f(p)$$

AGD aplicado a la 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 + \gamma\end{equation} donde $v\in \mathbb{R}^{m}$, $\gamma\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}=(\gamma^\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, es decir, 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: $$\gamma_{t+1}=\gamma_{t}+\frac{\alpha_{t}}{k}\sum_{i=1}^{k}(y_{i}-h_{\theta_{t}}(x_{i})),$$ $$\beta_{t+1}=\beta_{t}+\frac{\alpha_{t}}{k}\sum_{i=1}^{k}(y_{i}-h_{\theta_{t}}(x_{i}))x_{i}.$$ Implementando lo anterior en Python se tiene algo así:

Lo primero que haremos en importar las librerías necesarias para nuestro algoritmo.

import numpy as np
from sklearn.datasets.samples_generator import make_regression
import pylab
from scipy import stats

Definimos ahora el método de gradiente descendente en la siguiente función:

def gradient_descent(alpha: float, x: list, y: list, epsilon: float,
                     max_iter = 1000)-> float:
    converged = False
    iter = 0
    k = x.shape[0]
    # Inicializamos los parametros alpha y gamma
    alpha = np.random.random(x.shape[1])
    gamma = np.random.random(x.shape[1])
    # La función de error.
    temp_error = sum([(y[i] - gamma - beta * x[i]) ** 2 for i in range(k)])

    # A continuación se establece la iteración principal del AGD.

    while not converged:
        # El gradiente del error.
        grad0 = 1 / k * sum([(y[i] - gamma - beta * x[i]) for i in range(k)])
        grad1 = 1 / k * sum(
            [(y[i] - gamma - beta * x[i]) * x[i] for i in range(k)])
        # Se guarda temporalmente los parámetros gamma y alpha.
        temp_gamma = gamma + alpha * grad0
        temp_beta = beta + alpha * grad1

        # Se actualizan los parámetros gamma y alpha.
        gamma = temp_gamma
        beta = temp_beta

        error_square = sum([(y[i] - gamma - beta * x[i])**2 for i in range(k)])

        # Condiciones de parada.
        if abs(temp_error - error_square) <= epsilon:
            print('The iteration converges: ', iter)
            converged = True

        temp_error = error_square
        iter += 1

        if iter == max_iter:
            print('Maximum iterations exceeded')
            converged = True
    return gamma, beta

El lector puede notar que en la líneas 30 se estableció el criterio de parada $|e(\theta_t)-e(\theta_{t+1})|<\epsilon$ como $\verb|abs(temp_error - error_square) <= epsilon|$ donde $\verb|temp_error|=e(\theta_t)$ y $\verb|error_square|=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 algoritmo principal se tiene:

if __name__ == '__main__':
    kwargs = dict(n_samples=200, n_features=1, n_informative=1,
                  random_state=0, noise=35)

    x, y = make_regression(**kwargs)
    print('Test set size: ', x.shape, y.shape)

    alpha = 0.01
    epsilon = 0.01

    gamma, beta = gradient_descent(alpha, x, y, epsilon, max_iter=1000)
    print('gamma= ', gamma, ', beta =', beta)

    slope, intercept, r, p, slope_std_error = stats.linregress(x[:, 0], y)
    print('Intercept =', intercept, ',slope =', slope)

    y_predict = 0
    for i in range(x.shape[0]):
        y_predict = gamma + beta * x

    pylab.plot(x, y, 'o')
    pylab.plot(x, y_predict, 'k-')
    pylab.savefig('agd.png')
    pylab.show()
    print('Done!')

Cuando se ejecuta este algoritmo los resultados obtenidos son:

»»» Test set size:  (200, 1) (200,)
»»» The iteration converges:  749
»»» gamma=  [-4.85886785] , beta = [97.10010555]
»»» Intercept = -4.886236897512811 ,slope = 97.14302194296974
»»» Done!

La curva encontrada se representa en la siguiente figura:

Figura 3. Aplicación afín de parámetros $\gamma = -4.86$ y $\beta = 97.10$.

Conclusiones

  • Hoy aprendimos que el AGD es un algoritmo iterativo empleado principalmente para resolver problemas de optimización.
  • La principal desventaja del AGD se encuentra en el ajuste adecuado de la tasa de aprendizaje $\alpha$. Si $\alpha$ toma un valor muy pequeño es necesario una gran número de iteraciones para que el proceso converga, si por otro lado $\alpha$ es muy grande, entonces puede ocurrir que el proceso no converga.
  • Finalmente, el código lo encuentran en Github (aqui).

Referencias

0 comentarios: