Descenso de Gradiente y Variaciones Sobre el Tema

Mariano Rivera

agosto 2018

Una excelente compilación de métodos tipo descenso puede encontrase en el blog
[http://ruder.io/optimizing-gradient-descent/index.html#fn:6]

Descenso de Gradiente Simple (Descenso de Gradiente, GD)

Sea el problema de optimización
argminx    f(x) \underset{x}{\arg\min} \;\; f(x)

para ff suave (2 veces continuamente diferenciable).

Luego, sea xtx^t el punto actual. Entonces el vector
pt= f(xt) p^t = - \nabla\,f(x^t)
es una dirección de descenso (esto es, satisface p f(xt)<0p^\top \nabla\,f(x^t) <0).

Notación. Para simplificar nuestra notación definimos:
 ft=def f(xt). \nabla\,f^t \overset{def}{=} \nabla\,f(x^t).

Por lo que el punto xt+1x^{t+1} obtenido mediante la formula de actualización
xt+1=xtα  ft \boxed{ x^{t+1} = x^{t} - \alpha \, \nabla\,f^t }

donde α\alpha es el tamaño de paso y para una α\alpha sufientemente pequeña se garantiza:
f(xtα  ft)<ft f(x^{t} - \alpha \, \nabla\,f^t) < f^{t}
si  ft0\nabla\,f^t \ne 0.

Ejemplo GD

Calcular la raiz cuadrada de 2.

Notémos que x=2x=\sqrt{2} es solución de la ecuación no lineal:
f(x)=defx22=0. f(x) \overset{def}{=} x^2 -2 = 0.
Luego

# solo por usar sympy

import sympy as sym
from sympy.abc import x
sym.init_printing()
sym.integrate(x**2-2, x)

x332x\frac{x^{3}}{3} - 2 x

entonces definimos la función gg tal que f=gf= g'. Y resolvemos:
argminx    g(x)=x332x \underset{x}{\arg\min} \;\; g(x) = \frac{x^3}{3} -2 x

import matplotlib.pyplot as plt
import numpy as np

x = np.arange(-3,3,.1)  # from -3 to 3 in steps of .1
plt.plot(x,x**3/3-2*x)
[<matplotlib.lines.Line2D at 0x11d67cd68>]

png

Note que el mínimo esta precisamente en 2&ThinSpace;&ThinSpace;\sqrt{2}\,\, (y la otra raiz está en el máximo, en 2-\sqrt{2}).

Usamos descenso de gradiente para culcular el mínimo.

x = 0        # valor inicial de x
alpha = 0.2   # tamaño (pequeño) de paso

for t in range(40):
    x = x - alpha*(x**2-2)
    print('x1({0}) = {1}'.format(t,x))
    
x1(0) = 0.4
x1(1) = 0.768
x1(2) = 1.0500352
x1(3) = 1.229520415752192
x1(4) = 1.3271763252019033
x1(5) = 1.3748969255666177
x1(6) = 1.3968286143801103
x1(7) = 1.4066025787898986
x1(8) = 1.41089641585822
x1(9) = 1.4127706766019057
x1(10) = 1.4135864796686644
x1(11) = 1.413941132568255
x1(12) = 1.414095227294575
x1(13) = 1.414162164923116
x1(14) = 1.414191239183109
x1(15) = 1.4142038669866575
x1(16) = 1.4142093515066543
x1(17) = 1.41421173352888
x1(18) = 1.414212768078728
x1(19) = 1.4142132173993485
x1(20) = 1.4142134125459451
x1(21) = 1.4142134973009757
x1(22) = 1.4142135341113242
x1(23) = 1.414213550098596
x1(24) = 1.414213557042101
x1(25) = 1.4142135600577665
x1(26) = 1.4142135613675142
x1(27) = 1.4142135619363567
x1(28) = 1.4142135621834133
x1(29) = 1.4142135622907135
x1(30) = 1.4142135623373155
x1(31) = 1.4142135623575556
x1(32) = 1.414213562366346
x1(33) = 1.414213562370164
x1(34) = 1.414213562371822
x1(35) = 1.414213562372542
x1(36) = 1.414213562372855
x1(37) = 1.4142135623729908
x1(38) = 1.4142135623730498
x1(39) = 1.4142135623730754

que es muy cercano a

np.sqrt(2)

1.41421356237309511.4142135623730951

Ejemplo de regresión

Datos

from sklearn import linear_model, datasets

n_samples = 500
X, y = datasets.make_regression(n_samples=n_samples,
                                n_features=1,
                                n_informative=2, 
                                noise=5,
                                random_state=0) #2)
n_outliers=100
X[:n_outliers], y[:n_outliers] = datasets.make_regression(n_samples=n_outliers,
                                n_features=1,
                                n_informative=2, 
                                noise=2,
                                random_state=61)
y=np.expand_dims(y,axis=1)
plt.scatter(X[:],y[:], marker='.')

<matplotlib.collections.PathCollection at 0x120012eb8>

png

Agregamos in 1 a los datos de tal forma que nos queden en coordenadas homogéneas:
x^i=[xi,1] \hat{x}_i = [x_i, 1]
y en vector de corficientes para regresión lineal esta dado por
θ=[θ1,θ0] \theta = [\theta_1, \theta_0]^\top

Luego, la función objetivo la definimos como

f(θ)=12κ&ThinSpace;Ni=0N11exp(κ[θx^iyi]2) f(\theta) = \frac{1}{2 \kappa \, N} \sum_{i=0}^{N-1} 1- \exp\left( - \kappa [\theta^\top \hat x_i -y_i ]^2\right)

con derivadas parciales:
θ1f(θ)=1Ni=0N1xi[θx^iyi]exp(κ[θx^iyi]2) \frac{\partial}{\partial \theta_1}f(\theta) = - \frac{1}{N} \sum_{i=0}^{N-1} x_i [\theta^\top \hat x_i -y_i ] \exp\left( - \kappa [ \theta^\top \hat x_i -y_i]^2\right)
y
θ0f(θ)=1Ni=0N1[θx^iyi]exp(κ[θx^iyi]2) \frac{\partial}{\partial \theta_0}f(\theta) = - \frac{1}{N} \sum_{i=0}^{N-1} [\theta^\top \hat x_i -y_i] \exp\left( - \kappa [\theta^\top \hat x_i -y_i ]^2\right)

Entonces el gradiente de la función objetivo es

#-------------------------------------------------------------
def grad_quadratic(theta, f_params):
    '''
    Gradiente de la funcion de costo 
           sum_i (theta@x[i]-y[i])**2
    '''    
    X = f_params['X']
    y = f_params['y']

    err=theta[0]*X+theta[1]-y
    partial0=err
    partial1=X*partial0
    gradient= np.concatenate((partial1, partial0), axis=1)
    return np.sum(gradient, axis=1)
#-------------------------------------------------------------
def grad_exp(theta, f_params):
    '''
    Gradiente de la funcion de costo 
           sum_i 1-exp(-k(theta@x[i]-y[i])**2)
    '''
    kappa= f_params['kappa']
    X    = f_params['X']
    y    = f_params['y']
    err=theta[0]*X+theta[1]-y
    partial0=err*np.exp(-kappa*err**2)
    partial1=X*partial0
    gradient= np.concatenate((partial1, partial0), axis=1)
    return np.mean(gradient, axis=0)
#-------------------------------------------------------------
 

Implementación descenso de gradiente simple (GD)

def GD(theta=[], grad=None, gd_params={}, f_params={}):
    '''
    Descenso de gradiente
    
    Parámetros
    -----------
    theta     :   condicion inicial
    grad      :   función que calcula el gradiente
    gd_params :   lista de parametros para el algoritmo de descenso, 
                     nIter = gd_params[0] número de iteraciones
                     alpha = gd_params[1] tamaño de paso alpha

    f_params  :   lista de parametros para la funcion objetivo
                     kappa = f_params['kappa'] parametro de escala (rechazo de outliers)
                     X     = f_params['X'] Variable independiente
                     y     = f_params['y'] Variable dependiente                   
                  
    Regresa
    -----------
    Theta     :   trayectoria de los parametros
                     Theta[-1] es el valor alcanzado en la ultima iteracion
    '''
    
    nIter = gd_params['nIter'] 
    alpha = gd_params['alpha']
    Theta=[]
    for t in range(nIter):
        p = grad(theta,f_params=f_params)
        theta = theta - alpha*p
        Theta.append(theta)
    return np.array(Theta)
        

Descenso de Gradiente Estocástico (SGD)

Ahora considermos una variante del problema general, uno cuya función objetivo, o costo, se pueda denotar como la suma de muchos (si, muchos) pequeños costos. Esto es

argminx&ThickSpace;&ThickSpace;f(x)=def1ΩiΩfi(x) \underset{x}{\arg\min} \;\; f(x) \overset{def}{=} \frac{1}{\sharp \Omega} \sum_{i \in \Omega} f_i(x)

En este caso, la dirección de descenso de gradiente esta dado por
pt=&ThinSpace;ft&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;=1ΩiΩfit p^t = - \nabla\,f^t \\ \;\;\;\;\;\; = - \frac{1}{\sharp \Omega} \sum_{i \in \Omega} \nabla f_i^t

Note que, ptp^t puede ser interpretado como un valor esperado (promedio sobre toda la población). Esto es:

pt=E{fit}&ThickSpace;&ThickSpace;&ThickSpace;para&ThickSpace;&ThickSpace;i=1,2,,Ω; p^t = - \mathbb{E} \{\nabla f_i^t \} \;\;\; \text{para} \;\; i=1,2,\ldots, \sharp \Omega;

donde Ω\sharp \Omega denota la cardinalidad (número de elementos) en el conjunto Ω\Omega.

El término estocástico viene por el hecho de que, si en cada iteración, en vez de tomar la suma sobre toda la población, solo lo hacemnos sobre una muestra

StΩ S^t \subset \Omega
Luego

p~t=1StiStfit \tilde p^t = - \frac{1} {\sharp S^t} \sum_{i \in S^t} \nabla f_i^t

Esto equivale a calcular el gradiente como el promedio de los gradientes de la muestra.

Ventajas del SGD:

Desventajas del SGD:

[1] Herbert Robbins and Sutton Monro, A Stochastic Approximation Method, Ann. Math. Statist., Vol 22(3), 400-407 (1951).

Implementación descenso de gradiente estocástico (SGD)

def SGD(theta=[], grad=None, gd_params=[], f_params=[]):
    '''
    Descenso de gradiente estocástico
    
    Parámetros
    -----------
    theta     :   condicion inicial
    grad      :   funcion que calcula el gradiente
    
    gd_params :   lista de parametros para el algoritmo de descenso, 
                      nIter = gd_params['nIter'] número de iteraciones
                      alpha = gd_params['alpha'] tamaño de paso alpha
                      batch_size = gd_params['batch_size'] tamaño de la muestra
                      
    f_params  :   lista de parametros para la funcion objetivo, 
                      kappa = f_params['kappa'] parametro de escala (rechazo de outliers)
                      X     = f_params['X'] Variable independiente
                      y     = f_params['y'] Variable dependiente                   
                  
    Regresa
    -----------
    Theta     :   trayectoria de los parametros
                     Theta[-1] es el valor alcanzado en la ultima iteracion
    '''
    (high,dim) = f_params['X'].shape
    batch_size = gd_params['batch_size']
    
    nIter      = gd_params['nIter']
    alpha      = gd_params['alpha']
        
    Theta=[]
    for t in range(nIter):
        # Set of sampled indices
        smpIdx = np.random.randint(low=0, high=high, size=batch_size, dtype='int32')
        # sample 
        smpX = f_params['X'][smpIdx]
        smpy = f_params['y'][smpIdx]
        # parametros de la funcion objetivo
        smpf_params ={'kappa' : f_params['kappa'], 
                      'X'     : smpX , 
                      'y'     : smpy}
        
        p = grad(theta,f_params=smpf_params)
        theta = theta - alpha*p
        Theta.append(theta)
        
    return np.array(Theta)

Descenso de Gradiente con Momento (Inercia)

Descenso de gradiente es un método muy robusto y se mantiene aproximandose constantemente hacia un mínimo local. Este puede ser a su vez un problema:

Para reducir el efecto de los problemas mencionados, se ha propuesto incluir inercia. Esto es, mientras GD puede comprenderse con la analogía de un caminante que siempre dá un paso cosntante en la dirección que localmente tienen el mayor descenso. Descenso de Gardiente con Momento (MGD) sería el equivalente a una partícula masiva bajo el efecto de la gravedad. En tal caso, la partícula se acelera conforme acumule varios pasos en descenso. En el caso de la partícula, la acelaración puede incrementarse hasta alcanzar su velocidad límite: aquella en la que la fuerza de gravedad se equipare con la fuerza que ejerce la fricción de la superficie y del aire. La velocidad límite es la explicación al porqué un proyectil lanzado por un arma de fuego verticalmente no regresa con la misma velocidad con que salió del fusil, o el porqué la caida de los paracaidistas no se acelera durante todo el descenso; de hecho los paracaidistas maniobran para modificar su arrastre y así incrementar o decrementar la velocidad de caida; y así poder mantenerse acercarse y permanecer junto a otros miembros del grupo. Igualmente, podemos imponer a nuestra partícula en su descenso por la superficie de la función costo que además de momento (momentum), tenga una velocidad límite. Esto lo dejaremos para después, por lo pronto nos centraremos en el momento o inercia.

La modificación para incluir momento en GD consiste en lo siguiente

  1. Sea gt=ftg^t = \nabla f^t el gradiente en el punto actual; luego

  2. calcular la dirección de descenso con pt=gt+ηpt1p^t = g^t + \eta p^{t-1}; y

  3. actualizar el punto con xt+1=xtαptx^{t+1} = x^t - \alpha p^t.

La siguiente figura muestra gráficamente la actualización mediante MGD

Analizando el cálculo de la dirección de descenso observamos que:
pt=gt+ηpt1&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;=gt+ηgt1+ηpt2]&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;&ThickSpace;=gt+ηgt1+η2gt2+η3pt3&ThickSpace;&ThickSpace;=i=0tηigti p^t = g^t + \eta p^{t-1} \\ \;\;\;\;\;\;\;\; = g^t + \eta g^{t-1} + \eta p^{t-2}] \\ \;\;\;\;\;\;\;\;\;\;\;\; = g^t + \eta g^{t-1} + \eta^2 g^{t-2} + \eta^3 p^{t-3} \\ \;\; = \sum_{i=0}^{t} \eta^{i} g^{t-i}
donde asumimos que p(1)=0p^{(-1)}=0.

Notamos que se están integrando todos los gradientes de la trayectoria, pesando mas los gradientes mas recientes.

Para conservar el efecto de los gradientes recientes y suavizar la trayectoria, se recomienda η\eta cercana y menor que 1 , de hecho η=0.9\eta=0.9 es un valor comúnmente utilizado.

eta =0.9
t=10
etapow = [eta**i for i in range(t)]
print('Efecto de t (peso) en gradientes pasados')
print('\n'.join('{0} : {1:2.2E}'.format(*k) for k in enumerate(etapow)))

plt.plot(etapow)
plt.show()
Efecto de t (peso) en gradientes pasados
0 : 1.00E+00
1 : 9.00E-01
2 : 8.10E-01
3 : 7.29E-01
4 : 6.56E-01
5 : 5.90E-01
6 : 5.31E-01
7 : 4.78E-01
8 : 4.30E-01
9 : 3.87E-01

png

La integración de los gradientes pasados permite accumular la componente de los gradientes que apuntan en la misma dirección y cancelar las componentes normales a dicha trayectoria.

Implementación descenso de gradiente con momento (MGD)

def MGD(theta=[], grad=None, gd_params={}, f_params={}):
    '''
    Descenso de gradiente con momento (inercia)
    
    Parámetros
    -----------
    theta     :   condicion inicial
    grad      :   funcion que calcula el gradiente
    gd_params :   lista de parametros para el algoritmo de descenso, 
                      nIter = gd_params['nIter'] número de iteraciones
                      alpha = gd_params['alpha'] tamaño de paso alpha
                      eta   = gd_params['eta']  parametro de inercia (0,1]
    f_params  :   lista de parametros para la funcion objetivo, 
                      kappa = f_params['kappa'] parametro de escala (rechazo de outliers)
                      X     = f_params['X'] Variable independiente
                      y     = f_params['y'] Variable dependiente                   

    Regresa
    -----------
    Theta     :   trayectoria de los parametros
                     Theta[-1] es el valor alcanzado en la ultima iteracion
    '''
    nIter = gd_params['nIter']
    alpha = gd_params['alpha'] 
    eta   = gd_params['eta']
    p_old = np.zeros(theta.shape)
    Theta=[]
    for t in range(nIter):
        g = grad(theta, f_params=f_params)
        p = g + eta*p_old
        theta = theta - alpha*p
        p_old=p
        Theta.append(theta)
    return np.array(Theta)

Descenso acelerado de Nesterov (NAG)

El algoritmo de descenso con momento no ha demostrado que si bien el momento ayuda a accelerar la convergencia y reducir el riesgo de quedar atrapado en mínimos locales, también tienen una tendencia a sobrepasar los valles. Por ello, mediante una pequeña modificación se reducirá el effecto de sobrepaso (overshooting). Esta modificación es el algoritmo de Descenso Acelerado de Neasterov (NAG) [2].

El algoritmo NAG puede considerarse del tipo de dos pasos Predictor-Corrector. En el paso predictor se extrapola linealmente la trayectoria actual. Luego, en el punto predicho, se evalua el gradiente y se hace la correción de la trayectoria. Con ello se logra un aproximación de segundo orden de la trayectoria con un costo computacional similar al de Descenso de Gardiente con Momento (MGD). Al mantener el efecto de inercia pero reducir los sobrepasos, NAG también se incrementa la razón de convergencia

Las fórmulas de actualización del algoritmo NAG estan dadas por los siguientes puntos

  1. Sea x~t=xtαpt1\tilde x^{t} = x^t - \alpha p^{t-1} el punto predicho en primer orden (usando la dirección anterior); luego

  2. Calcular el gradiente en el punto predicho: gt=f(x~t)g^t = \nabla f(\tilde x^{t}); entonces

  3. Calcular la dirección de descenso con pt=gt+ηpt1p^t = g^t + \eta p^{t-1}; y

  4. Actualizar el punto con xt+1=xtαptx^{t+1} = x^t - \alpha p^t

La siguiente figura muestra gráficamente la actualización mediante NAG

[2] Nesterov, Y., A method for unconstrained convex minimization problem with the rate of convergence o(1/k2). Doklady ANSSSR (translated as Soviet.Math.Docl.), vol. 269, pp. 543-547 (1993)

Implementación del algoritmo NAG

def NAG(theta=[], grad=None, gd_params={}, f_params={}):
    '''
    Descenso acelerado de Nesterov
    
    Parámetros
    -----------
    theta     :   condicion inicial
    grad      :   funcion que calcula el gradiente
    gd_params :   lista de parametros para el algoritmo de descenso, 
                      nIter = gd_params['nIter'] número de iteraciones
                      alpha = gd_params['alpha'] tamaño de paso alpha
                      eta   = gd_params['eta']  parametro de inercia (0,1]
    f_params  :   lista de parametros para la funcion objetivo, 
                      kappa = f_params['kappa'] parametro de escala (rechazo de outliers)
                      X     = f_params['X'] Variable independiente
                      y     = f_params['y'] Variable dependiente                   

    Regresa
    -----------
    Theta     :   trayectoria de los parametros
                     Theta[-1] es el valor alcanzado en la ultima iteracion
    '''
    nIter = gd_params['nIter']
    alpha = gd_params['alpha'] 
    eta   = gd_params['eta']
    p     = np.zeros(theta.shape)
    Theta=[]
    
    for t in range(nIter):
        pre_theta = theta - 2.0*alpha*p
        g = grad(pre_theta, f_params=f_params)
        p = g + eta*p
        theta = theta - alpha*p
        Theta.append(theta)
    return np.array(Theta)

Descenso de Gradiente Adaptable (ADAGRAD)

Considere el caso en que la magnitud del gradiente de la función es relativamente grande (aun estamos lejos del óptimo). Dado que la magitud del gradiente es una cantidad a la que contribuyen todas las derivadas parciales, es posible que a pesar de tener una magnitud grande, existan coordenadas donde la función cambia poco (derivadas parciales de norma pequeña). Lo que convienen hacer es que se den pasos grandes en las componentes del gradiente que son grandes (en valor absoluto) y pasos cortos en las componentes pequeñas.

Los algortimos MGD como NAG se enfocan en estimar mejor la dirección de descenso. Sin embargo, el desempeño de dichos algoritmos se ve comprometido por la correcta selección del tamaño de paso α\alpha. Para ello, el algoritmo ADAGRAD [3] busca estimar el correcto tamaño de paso en cada iteración estimando un tamaño de paso para cada una de las variables en forma independiente.

Las fórmulas de actualización del algoritmo NAG estan dadas por los siguientes puntos

  1. Sea
    git=defxift g^t_i \overset{def}{=} \frac{\partial }{\partial x_i}f^t
    la ii-ésima derivada parcial de la función en el punto actual xtx^{t}; esto es git=[ft]ig^t_i = [\nabla f^t]_i, entonces

  2. Calcular la suma de los cuadrados de las parciales hasta la iteración actual:
    Git=defk=0t(gik)2; G^t_i \overset{def}{=} \sum_{k=0}^t (g^k_i)^2;
    y

  3. Actualiza el punto con
    xt+1=xtαGit+ϵpt x^{t+1} = x^t - \frac{\alpha}{\sqrt{G_i^t + \epsilon}} p^t
    donde ϵ\epsilon es una constante pequeña que evita la división por cero.

[3] Duchi, J., Hazan, E., & Singer, Y., Adaptive Subgradient Methods for Online Learning and Stochastic Optimization. Journal of Machine Learning Research, 12, pp 2121–2159 (2011). see http://jmlr.org/papers/v12/duchi11a.html

Descenso de Gradiente Adaptable (ADADELTA)

Note que el término GitG^t_i en ADAGRAD se mantienen acumulando las derivadas parciales desde el inicio de las iteraciones. Esto puede reducir en forma temprana el tamaño de paso para algunos parámetros. Por ello, es mejor ir reduciendo paulatinamente la contribución de términos pasados.

Para ello, es conveniente realizar ajuste similar a la intergación de las direcciones pasadas:

  1. Sea
    git=defxift g^t_i \overset{def}{=} \frac{\partial }{\partial x_i}f^t
    la ii-ésima derivada parcial de la función en el punto actual xtx^{t}; esto es git=[ft]ig^t_i = [\nabla f^t]_i, entonces

  2. integrar los cuadrados de las parciales hasta la iteración actual:
    Git=η(git)2+(1η)Git1 G^t_i = \eta (g^t_i)^2 + (1-\eta) G^{t-1}_i
    con η\eta grandes (típicamente η=0.9\eta =0.9); y

  3. actualizar el punto, donde cada elemento tienen su propio paso
    xit+1=xitαGit+ϵpit x^{t+1}_i = x^t_i - \frac{\alpha}{\sqrt{G_i^t} + \epsilon} p^t_i

[4] Bengio, Y., Boulanger-Lewandowski, N., & Pascanu, R., Advances in Optimizing Recurrent Networks (2012)

Implementación del algoritmo ADADELTA

def ADADELTA(theta=[], grad=None, gd_params={}, f_params={}):
    '''
    Descenso de Gradiente Adaptable (ADADELTA) 
    
    Parámetros
    -----------
    theta     :   condicion inicial
    grad      :   funcion que calcula el gradiente
    gd_params :   lista de parametros para el algoritmo de descenso, 
                      nIter    = gd_params['nIter'] número de iteraciones
                      alphaADA = gd_params['alphaADADELTA'] tamaño de paso alpha
                      eta      = gd_params['eta']  parametro adaptación del alpha
    f_params  :   lista de parametros para la funcion objetivo, 
                      kappa = f_params['kappa'] parametro de escala (rechazo de outliers)
                      X     = f_params['X'] Variable independiente
                      y     = f_params['y'] Variable dependiente                   

    Regresa
    -----------
    Theta     :   trayectoria de los parametros
                     Theta[-1] es el valor alcanzado en la ultima iteracion
    '''
    epsilon= 1e-8
    nIter    = gd_params['nIter']
    alpha    = gd_params['alphaADADELTA'] 
    eta      = gd_params['eta']
    G        = np.zeros(theta.shape)
    g        = np.zeros(theta.shape) 
    Theta=[]
    for t in range(nIter):
        g = grad(theta, f_params=f_params)
        G = eta*g**2 + (1-eta)*G
        p = 1.0/(np.sqrt(G)+epsilon)*g
        theta = theta - alpha * p
        Theta.append(theta)
    return np.array(Theta)

Momentum Adaptable (ADAM)

El algoritmo ADAM [5]_ Calcula la dirección de descenso usando momentum (similar a MGD) y utiliza una estrategia similar para calcular adaptar el tamaño de paso. Es decir, utiliza momuntum para actualizar el paso, lo que evita cambios bruscos en el paso. Esto lo hace muy estable para su uso en estrategia tipo Gradiente Estocástico (SGD) donde las muestras pueden provocar cambios grandes en la magnitud del gradiente, además calcula un paso global en vez de usar un paso para cada variable. También adecuado en estrategias de entrenamiento tipo estocásticas o por lotes, como en el caso de Redes Neuronales Profundas (Deep Learning)

Una mejora importante es la correción del sezgo (bias) en la estimación de los momentos.

Generalmente las razones de aprendizaje (momentum) son cercanas a 1, típicamente: η1=0.9\eta_1=0.9 y η2=0.999\eta_2=0.999.

[5] D. P. Kingma and J. L. Ba. Adam: a Method for Stochastic Optimization. In procc. ICLR 2015, 1–13 (2015)

Una iteración del algoritmo ADAM se resume en los siguientes pasos

Sea
gt=defft g_t \overset{def}{=} \nabla f_t
la ii-ésima derivada parcial de la función en el punto actual xtx^{t}. Entonces

  1. Calcular la dirección de descenso con momentum
    pt=η1pt1+(1η1)gt p_t = \eta_1 p_{t-1} + (1-\eta_1) g_t
    donde ptp_t conserva la escala de pt+1p_{t+1} y el gradiente
    gtg^t.

  2. Luego, actualizar, a la vez con momentum, el factor adaptativo del descenso vtv^{t} :
    vt=η2vt1+(1η2)Gt. v_t = \eta_2 v_{t-1} + (1-\eta_2) G_t.
    donde G=[gi2,g22,,gn2]\overline G = [g_i^2, g_2^2, \ldots, g_n^2]^\top es el vector de elementos a cuadrado del gradiente.

  3. Escalar la dirección de descenso y el momentum (reducción del sezgo):
    p^t=11(η1)tptv^t=11(η2)tvt \hat p_t = \frac{1}{1- (\eta_1)^t} p_t \\ \hat v_t = \frac{1}{1- (\eta_2)^t} v_t

  4. Actualizar el punto con la fórmula de paso adaptable
    xt+1=xtαv^t+ϵp^t x_{t+1} = x_t - \frac{\alpha}{\sqrt{\hat v_t} + \epsilon} \hat p_t

Importancia de la reducción del sezgo

El paso 1 del algoritmo ADAM pretende mejorar el cálculo del gradiente, calculado general mente sobre solo una muestra, promediandolo con los gradientes recientemente calculados:

pt=η1pt1+(1η1)gt=η12pt2+η1(1η1)gt1+(1η1)gt=η13pt3+η12(1η1)gt2+η1(1η1)gt1+(1η1)gt p_t = \eta_1 p_{t-1} + (1-\eta_1) g_t \\ = \eta_1^2 p_{t-2} + \eta_1 (1-\eta_1) g_{t-1} + (1-\eta_1) g_t \\ = \eta_1^3 p_{t-3} + \eta_1^2 (1-\eta_1) g_{t-2} + \eta_1 (1-\eta_1) g_{t-1} + (1-\eta_1) g_t \\ \ldots

Si asumimos p0=0p_0=0 como condicion inicial de las iteraciones, nos lleva a
pt=(1η1)i=0tη1tigi p_t = (1-\eta_1) \sum_{i=0}^{t} \eta_1^{t-i} g_{i}

Tomando valores esperados, y usando E{xy}=E{x}E{y}\mathbb{E}\{xy\} = \mathbb{E}\{x\} \mathbb{E}\{y\}, si xx y yy son independientes:

E{pt}=E{gt}(1η1)i=0tη1ti \mathbb{E}\{p_t\} = \mathbb{E}\{g_t\} (1-\eta_1) \sum_{i=0}^{t} \eta_1^{t-i}
Luego usamos la propiedad de la serie geométrica
k=0n1xk=1xn1x \sum_{k=0}^{n-1} x^k = \frac{1-x^n}{1-x}

lo que resulta en
E{pt}=E{gt}(1η1t) \mathbb{E}\{p_t\} = \mathbb{E}\{g_t\} (1-\eta_1^{t})

por loque para corregur este facor se hace
p^t=pt1ηt \hat p_t = \frac{p_t}{1-\eta^t}

Implementación del algoritmo ADAM

def ADAM(theta=[], grad=None, gd_params={}, f_params={}):
    '''
    Descenso de Gradiente Adaptable con Momentum(A DAM) 
    
    Parámetros
    -----------
    theta     :   condicion inicial
    grad      :   funcion que calcula el gradiente
    gd_params :   lista de parametros para el algoritmo de descenso, 
                      nIter    = gd_params['nIter'] número de iteraciones
                      alphaADA = gd_params['alphaADAM'] tamaño de paso alpha
                      eta1     = gd_params['eta1'] factor de momentum para la direccion 
                                 de descenso (0,1)
                      eta2     = gd_params['eta2'] factor de momentum para la el 
                                 tamaño de paso (0,1)
    f_params  :   lista de parametros para la funcion objetivo, 
                      kappa = f_params['kappa'] parametro de escala (rechazo de outliers)
                      X     = f_params['X'] Variable independiente
                      y     = f_params['y'] Variable dependiente                   

    Regresa
    -----------
    Theta     :   trayectoria de los parametros
                     Theta[-1] es el valor alcanzado en la ultima iteracion
    '''
    epsilon= 1e-8
    nIter    = gd_params['nIter']
    alpha    = gd_params['alphaADAM'] 
    eta1     = gd_params['eta1']
    eta2     = gd_params['eta2']
    p        = np.zeros(theta.shape)
    v        = 0.0
    Theta    = []
    eta1_t = eta1
    eta2_t = eta2
    for t in range(nIter):
        g  = grad(theta, f_params=f_params)
        p  = eta1*p + (1.0-eta1)*g
        v  = eta2*v + (1.0-eta2)*(g**2)
        #p = p/(1.-eta1_t)
        #v = v/(1.-eta2_t)
        theta = theta - alpha * p / (np.sqrt(v)+epsilon)
        eta1_t *= eta1
        eta2_t *= eta2
        Theta.append(theta)
    return np.array(Theta)

Evaluación del los algoritmos de descenso de gradiente

# condición inicial
theta=10*np.random.normal(size=2)
#theta= [-0.61752689 -0.76804482]
# parámetros del algoritmo
gd_params = {'alpha'          : 0.95, 
             'alphaADADELTA'  : 0.7,
             'alphaADAM'      : 0.95,
             'nIter'          : 300,
             'batch_size'     : 100,
             'eta'            : 0.9,
             'eta1'           : 0.9,
             'eta2'           : 0.999}

# parámetros de la función objetivo
f_params={'kappa' : 0.01, 
          'X'     : X , 
          'y'     : y}

Descenso de Gradiente

ThetaGD = GD(theta=theta, grad=grad_exp, 
             gd_params=gd_params, f_params=f_params)
print('Inicio:', theta,'-> Fin:', ThetaGD[-1,:])
Inicio: [  6.92950223 -10.81700592] -> Fin: [45.67020876 -0.27659597]

Descenso de Gradiente Estocástico

ThetaSGD = SGD(theta=theta, grad=grad_exp, 
               gd_params=gd_params, f_params=f_params)
print('Inicio:', theta,'-> Fin:', ThetaSGD[-1,:])
Inicio: [  6.92950223 -10.81700592] -> Fin: [45.89263223 -0.2931078 ]

Descenso de Gradiente con Inercia

ThetaMGD = MGD(theta=theta, grad=grad_exp, 
               gd_params=gd_params, f_params=f_params)
print('Inicio:', theta,'-> Fin:', ThetaMGD[-1,:])
Inicio: [  6.92950223 -10.81700592] -> Fin: [45.67020783 -0.27659608]

Nesterov

ThetaNAG = NAG(theta=theta, grad=grad_exp, 
               gd_params=gd_params, f_params=f_params)
print('Inicio:', theta,'-> Fin:', ThetaMGD[-1,:])
Inicio: [  6.92950223 -10.81700592] -> Fin: [45.67020783 -0.27659608]

ADADELTA

ThetaADADELTA = ADADELTA(theta=theta, grad=grad_exp, 
                         gd_params=gd_params, f_params=f_params)
print('Inicio:', theta,'-> Fin:', ThetaADADELTA[-1,:])
Inicio: [  6.92950223 -10.81700592] -> Fin: [46.02421512 -0.62897892]

ADAM

ThetaADAM = ADAM(theta=theta, grad=grad_exp, 
                 gd_params=gd_params, f_params=f_params)
print('Inicio:', theta,'-> Fin:', ThetaADAM[-1,:])
Inicio: [  6.92950223 -10.81700592] -> Fin: [45.67020748 -0.27659608]
Tmax=100
plt.figure(figsize=(10,10))

plt.subplot(211)
plt.plot(ThetaNAG[:Tmax], '.')
plt.title('NAG')

plt.subplot(212)
plt.plot(ThetaADAM[:Tmax], '.')
plt.title('ADAM')

plt.show()

png

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

mpl.rcParams['legend.fontsize'] = 14

fig = plt.figure(figsize=(15,15))
ax = fig.gca(projection='3d')
nIter=np.expand_dims(np.arange(ThetaGD.shape[0]),1) 
Tmax=200
ax.plot(ThetaGD[:Tmax,0],  ThetaGD [:Tmax,1], nIter[:Tmax,0], label='GD')
ax.plot(ThetaSGD[:Tmax,0], ThetaSGD[:Tmax,1], nIter[:Tmax,0], label='SGD')
ax.plot(ThetaMGD[:Tmax,0], ThetaMGD[:Tmax,1], nIter[:Tmax,0], label='MGD')
ax.plot(ThetaNAG[:Tmax,0], ThetaNAG[:Tmax,1], nIter[:Tmax,0], label='NAG')
ax.plot(ThetaADADELTA[:Tmax,0], ThetaADADELTA[:Tmax,1], nIter[:Tmax,0], label='ADADELTA')
ax.plot(ThetaADAM[:Tmax,0], ThetaADAM[:Tmax,1], nIter[:Tmax,0], label='ADAM')
ax.legend()
ax.set_title(r'Trayectorias los parámetros calculados con distintos algoritmos')
ax.set_xlabel(r'$\theta_1$')
ax.set_ylabel(r'$\theta_0$')
ax.set_zlabel('Iteración')
plt.show()

png