7 ene 2020

El algoritmo del gradiente descendente

El gradiente descendente (GD) es un algoritmo de optimización genérico, capaz de encontrar soluciones óptimas para una amplia gama de problemas. La idea del gradiente descendente es ajustar los parámetros de forma iterativa para minimizar una función.

Concretamente, se tiene de una función diferenciable convexa1, $f:\Omega\subset\mathbb{R}^n \to \mathbb{R}$, el algoritmo GD permite encontrar un $w$ en $\Omega$ tal que $f(w)$ es un mínimo, en otras palabras, GD se utiliza para determinar los elementos del siguiente conjunto2:

\begin{equation}\tag{1} w \in \operatorname*{argmin\,\,}_{ w\in \Omega} f(w). \end{equation}

Para determinar los valores de $w$ que optimiza la función $f(w)$, GD hace uso de una serie de iteracciones que se hacen de acuerdo con la siguiente regla de actualización:

\begin{equation}\tag{2} w_{t+1} = w_{t} -\eta_{t} \nabla f(w_{t}), \end{equation}

que usualmente se inicializa en cero y cada iteración, como se puede observar, se hace en la dirección negativa del gradiente. Recuerde que el gradiente se define como el vector:

\begin{equation}\nonumber \nabla f(w)=\left(\frac{\partial f}{\partial x_{1}}(w),\dots, \frac{\partial f}{\partial x_{n}}(w)\right). \end{equation}

Un hiperparámetro importante en GD es el $\eta_{t} > 0$, denominado como tasa de aprendizaje. Si la tasa de aprendizaje es demasiado pequeña, entonces el algoritmo tendrá que pasar por muchas iteraciones para converger, lo que llevará mucho tiempo. Por otro lado, si la tasa de aprendizaje es demasiado alta, es posible que se salte el mínimo global y termine en otro lado, posiblemente incluso más alto que antes. Esto podría hacer que el algoritmo diverja, con valores cada vez mayores, sin encontrar una buena solución.

¿Cómo funciona este algoritmo?

Este algoritmo se define a partir de las dos características esenciales que tiene el gradiente, las cuales se mencionan a continuación:

  1. El gradiente es perpendicular a las curvas de nivel de $f$, de manera que para cualquier dirección $v\in \mathbb{R}^{n}$ ortogonal a $\nabla f(p)$, es una dirección de cambio nulo. Esto se observa facilmente al parametrizar la curva $S_{k}=\{p\in \Omega: f(p)=k\}$ mediante una función $\alpha:I\subset \mathbb{R}\to S_{k}$ tal que $\alpha(0)=p$, pues al calcular el producto punto de $\nabla f(p)$ con la velocidad de $\alpha$ en $p$ se obtiene que la tasa de cambio es:

    \begin{equation}\nonumber|df_{p}(\alpha'(0))|=|\nabla f(p)\cdot \alpha' (0)| = 0, \end{equation}
    es decir, la tasa de cambio en la dirección de $\alpha'(0)$ es cero.



  1. El gradiente indica la dirección ascendente de la tasa de máximo cambio de $f$ en el punto $p$. La tasa máxima se calcula como $||\nabla f(p)||$. La razón de esto, se aprecia cuando se considera un vector $v\in \mathbb{R}$ tal que $||v||=1$, de manera que para este vector la tasa de cambio es:

    \begin{equation}\nonumber |df_{p}(v)|=||\nabla f(p)||||v|||cos\theta|\leq ||\nabla f(p)|| \end{equation}
    Dicha magnitud es máxima cuando $\theta = 2n\pi$ con $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 de $-\nabla f(p)$, en efecto, si $v=\frac{\nabla f(p)}{||\nabla f(p)||}$, entonces $df_{p}(v)=||\nabla f(p)||$.

La iteración definida en la Ecuación 2 permite construir una sucesión de puntos $\{w_t\}_{t\in [m]_{\mathbb{N}_0}}$ de tal manera que $f(w_{t+1})<f(w_{t})$. Este hecho se puede evidenciar mediante el polinomio de Taylor, en efecto si se considera un punto inicial $w_{0}$, entonces para el primer termino de la expansión de Taylor alrededor de $w_{0}$ se tiene:

\begin{equation}\nonumber f(w_{1})-f(w_{0})\approx \langle w_{1}-w_{0}, \nabla f(w_{0})\rangle=-\eta ||\nabla f(w_{o})||^2 \end{equation}

Por consiguiente:

\begin{equation}\nonumber f(w_{1})-f(w_{0}))=-\eta ||\nabla f(w_{o})||^2 + o(\eta) \end{equation}

de tal manera que para un $\eta$ adecuado se puede garantizar que $f(w_{1})<f(w_{0})$. El razanomiento se puede repetir para $w_{t+1}$ y $w_{t}$, de tal manera que $w_{t+1}$ es una mejora de $w_{t}$.

De acuerdo a las consideraciones anteriores, el algoritmo del gradiente descendente se define de la siguiente manera:


Algoritmo del gradiente descendente (AGD)

Input: $w_0$, $m$, $\eta$, $\nabla f(w)$ 
1.     for $k=0$ to $m$ do: 
2.          $w\leftarrow w - \eta \nabla f(w)$ 
3.     end 
return: $w$

donde $w_0$ es la condición inicial del algoritmo, $m$ es el número máximo de iteraciones, $\eta$ es la tasa de aprendizaje y $\nabla f(w)$ es la función gradiente de $f$. Observe que la regla de actualización se definió apartir de la Ecuación 2. Es importante notar que finalmente $w_{t+1}$ es tal que:

\begin{equation}\nonumber w_{t+1}\in \operatorname*{argmin\,\,}_{ w\in \Omega} \frac{1}{2\eta_{t}}||w-w_{t}+\eta_{t}\nabla f(w_{t})||. \end{equation}

Es fácil comprobar que el problema de optimización anterior se puede reescribir como:

\begin{equation} \operatorname*{argmin\,\,}_{ w\in \Omega} \frac{1}{2\eta_{t}}||w-w_{t}+\eta_{t}\nabla f(w_{t})||=\operatorname*{argmin\,\,}_{ w\in \Omega} \left(f(w_{t}) + \langle \nabla f(w_t), w-w_t \rangle + \frac{1}{2\eta_{t}}||w-w_{t}||\right) \end{equation}

Por lo tanto, el $w_{t+1}$ es obtenido para minimizar la linealización de la función $f$ alrededor del punto $w_{t}$, manteniendolo lo suficientemente aproximado a este el punto $w_{t}$.

¿Cuándo parar las iteraciones?

Es importante aclarar, que a pesar de que el algoritmo ejecute el máximo de iteraciones, el resultado arrojado no necesariamente es una buena aproximación a el elemento minimizador de la función $f$. Por lo tanto es necesario definir un criterio que permita decidir si el resultado obteniendo es adecuado o no.

Como primer criterio de parada del algoritmo que se le puede ocurrir al lector es detener las iteraciones cuando $||\nabla f(x_{t})||=0$, pero esto no es práctico debido a diferentes factores que influyen, como el comportamiento de la coma flotante y la elección adecuada de la tasa de aprendizaje. Usualmente en la práctica se suele definir un parametro de tolerancia $\epsilon>0$ que junto con alguno de los siguientes criterios se usan para detener las iteraciones:

  • Condición sobre el gradiente:

    \begin{equation}\nonumber ||\nabla f(x_{t})||<\epsilon\end{equation}
  • Condición sobre las diferencias sucesivas relativas de la función objetivo:

    \begin{equation}\nonumber \frac{|f(x_{t+1})-f(x_{t})|}{|f(x_{t})|}<\epsilon,\end{equation}
    si el denominador es muy pequeño, es conveniente remplazarlo por $\max\{1, |f(x)|\}$.
  • Condición sobre las diferencias sucesivas relativas de la variable independiente:

    \begin{equation}\nonumber \frac{||x_{t+1}-x_{t}||}{||x_{t}||}<\epsilon,\end{equation}
    si el denominador es muy pequeño, es conveniente remplazarlo por $\max\{1, ||x_{t}||\}$.

Algunas consideraciones sobre las tasas de aprendizaje.

  • La principal desventaja del AGD se encuentra en el ajuste adecuado de la tasa de aprendizaje $\eta$. Si $\eta$ toma un valor muy pequeño, es necesario un gran número de iteraciones para que el proceso converga; si por otro lado $\eta$ es muy grande, entonces puede ocurrir que el proceso no converga.

  • La tasa de aprendizaje $\eta$ es determinada por la minimización exacta de:

    \begin{equation}\nonumber\eta_{t} \in \operatorname*{argmin\,\,}_{\eta>0}f(x_{t}-\eta\nabla f(x_{t})). \end{equation}
    Esto se usa principalmente para problemas de caracter cuadrático y en donde el cálculo de $\eta$ es económico pero una evaluación de gradiente costosa; de lo contrario, no vale la pena el esfuerzo de resolver este subproblema exactamente.

Ejemplo: Gradiente descente para una forma cuadrática

Asuma que $Q$ es simétrica y definida positiva ($x^{\top}Q x>0$ para cualquier $x\neq 0$). Considere la forma cuadrática:

\begin{equation}\nonumber f(x)=\frac{1}{2}x^{\top}Qx - b^{\top}x \end{equation}

el lector puede comprobar que su gradiente es:

\begin{equation}\nonumber \nabla f(x)=Qx-b. \end{equation}

Así la secuencia de $\{x_{t}\}_{t\in [m]}$ que inicia en cualquier $x_{0}$ viene dada por:

\begin{equation}\nonumber x_{t+1}=x_{t}-\eta_{k}(Qx-b) \end{equation}

con $g_{t}:=\nabla f(x_t)$ se define:

\begin{equation}\nonumber \eta_{t}=\frac{g_{t}^{\top}g_{t}}{g_{t}^{\top}Qg_{t}}. \end{equation}

El lector puede validar que con el valor de $\eta_{t}$ definido anteriormente se tiene:

\begin{equation}\nonumber \eta_{t}\in \operatorname*{argmin\,\,}_{\eta>0}f(x_{t}-\eta\nabla f(x_{t})). \end{equation}

¿Cómo implementarlo en python?

Para las consideraciones del ejemplo anterior, es fácil definir una clase en Python para todas las formas cuadraticas, junto con tres métodos principales que permiten evaluar la forma cuadrática, calcular $\eta$ y el gradiente en un punto $x$. Esto sería algo así:

In [1]:
import numpy as np
class QuadraticForm:
def __init__(self, Q, b):
"""
Inputs:
Q: Positive definite symmetric matrix.
b: Rn vector.
"""
self.Q = np.array(Q)
self.b = np.array(b).reshape(-1, 1)
def evaluate(self, x):
"""
Method to evaluate the quadratic form.
Inputs:
x: Rn vector.
Ouput:
Value of the quadratic form in x.
"""
x = np.array(x).reshape(-1, 1)
Q = self.Q
b = self.b
return (1/2) * x.T.dot(Q.dot(x)) - b.T.dot(x)
def eta(self, x):
"""
Method to evaluate eta.
Inputs:
x: Rn vector.
Output:
Value of eta in x.
"""
gradient_x = self.gradient(x)
numerator = gradient_x.T.dot(gradient_x)
denominator = gradient_x.T.dot(Q.dot(gradient_x))
return numerator / denominator
def gradient(self, x):
"""
Method to evaluate the gradient.
Inputs:
x: Rn vector.
Output:
Value of the gradient in x.
"""
x = np.array(x).reshape(-1, 1)
return Q.dot(x) - b

El AGD se implementaría de la siguiente forma:

In [2]:
def gradient_descent(initial_x, eta, epsilon, function):
"""
Gradient descent Algorithm.
Inputs:
initial_x: Rn vector.
eta: learning rate.
epsilon: precision.
function: Quadractic Form.
Output:
Point where the function reaches the minimum.
"""
k = 0
x = np.array(initial_x).reshape(-1, 1)
while True:
gradient_x = function.gradient(x)
if np.linalg.norm(gradient_x) < epsilon:
print('stop: {}'.format(k))
break
if isinstance(eta, (int, float)):
x = x - eta * gradient_x
else:
x = x - eta(x) * gradient_x
k += 1
return x

Para ejemplicar el funcionamiento del código anterior, se considera la matriz simétrica y positiva definida:

\begin{equation}\nonumber Q = \left(\begin{array}{cc} 1 & 0.5 \\ 0.5 & 3 \end{array}\right) \end{equation}

y vector $b$ dado por:

\begin{equation}\nonumber b = \left(\begin{array}{c} 3 \\ 0.5 \end{array}\right) \end{equation}

De esta manera la forma cuadrática en Python y usando código anterior quedaría así:

In [3]:
Q = np.array([[1, 0.5], [0.5, 3]])
b = np.array([3, 0.5]).reshape(-1, 1)
In [4]:
quadratic_form = QuadraticForm(Q, b) 

Para encontrar el $x$ minimizador de la función se tiene dos opciones, una es definir el parámetro $\eta$ manualmente, y la otra es usar el método eta de la clase QuadracticForm. A continuación se hará uso de los dos casos.

Para empezar se ejecuta el algoritmo en algún punto $x$:

In [5]:
initial_x = np.array([105.5, 105.8])

Luego ejecuta el AGD con $\eta = 0.0001$ y $\epsilon=0.0000001$:

In [6]:
gradient_descent(initial_x, eta=0.0001, epsilon=0.0000001, function=quadratic_form)
stop: 230300
Out[6]:
array([[ 3.18181829],
[-0.36363639]])

Como el lector podrá notar el algoritmo tardó 230300 iteraciones para obtener un candidato al mínimo con la precisión deseada. A continuación se ejecuta usando el método eta de la forma cuadrática:

In [7]:
gradient_descent(initial_x, eta=quadratic_form.eta, epsilon=0.0000001, function=quadratic_form)
stop: 15
Out[7]:
array([[ 3.1818182 ],
[-0.36363637]])

En esta ocasión el algoritmo tardó 15 iteraciones. Esto representa una mejora en tiempo de ejecución bastante considerable con respecto a el experimento anterior.

¿Cómo se puede asegurar que este es el mínimo de la función? Para este caso en particular, el mínimo ocurre cuando $Qx = b$, por lo tanto es suficiente con resolver este sistema lineal. Usando los métodos de la librería numpy se puede resolver rápidamente así:

In [8]:
np.linalg.solve(Q, b)
Out[8]:
array([[ 3.18181818],
[-0.36363636]])

Como se puede observar, el resultado es muy aproximado al valor que se obtuvo al ejecutar AGD, por lo que el algoritmo funciona bastante bien.

Conclusiones

  • Se logró aprender 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 $\eta$. Si $\eta$ toma un valor muy pequeño, es necesario un gran número de iteraciones para que el proceso converga; si por otro lado $\eta$ es muy grande, entonces puede ocurrir que el proceso no converga.
  • Finalmente te invito a leer este otro artículo (Jugando con el gradiente descendente y Python) en donde se análiza el comportamiento general de los algoritmos por gradiente descendente.

No olvides comentar y suscribirte al blog para que estés enterando de los posts que voy a ir subiendo semana a semana. También sientete en libertad de seguirme en LinkedIn, Twitter, Github e Instagram.

Notas

[1]. En el contexto de este artículo, se dira que una función $f:\Omega \subset \mathbb{R}^{m}\to \mathbb{R}$ es diferenciable si tiene derivada continua en $\Omega$ y es convexa, si para todo $w, z\in \Omega$ y $\alpha \in [0, 1]$ se cumple la siguiente condición:

\begin{equation}\nonumber f(\alpha w + (1-\alpha)z)\leq \alpha f(w) + (1-\alpha)f(z). \end{equation}

Note que la condición es valida para todo $m\geq 1$.



[2]. Si $y_{\min}$ es el mínimo global de una función $f(x)$, entonces se define el operador $\operatorname*{argmin\,\,}$ como siguiente conjunto: $$\operatorname*{argmin\,\,}_{ x \in X} f(x) = \{x\in X: f(x) = y_\min\}.$$

Referencias

Contacto

  • Participa de la canal de Nerve a través de Discord.
  • Se quieres conocer más acerca de este tema me puedes contactar a través de Classgap.