<span xmlns:dct="http://purl.org/dc/terms/" property="dct:title">Guia3 p3</span> por <a xmlns:cc="http://creativecommons.org/ns#" href="http://nbviewer.ipython.org/urls/raw.githubusercontent.com/elsuizo/Redes_neuronales_Fuzzy/master/guia1.ipynb?create=1" property="cc:attributionName" rel="cc:attributionURL">Martin Noblía Marcos Panizzo Damián Presti</a> se distribuye bajo una Licencia Creative Commons Atribución-CompartirIgual 4.0 Internacional.
Minimizar la función:
$$F(x_{1},x_{2})=(x_{2} - 2 \, x_{1} + e^{-x_{1}})^{2} + (x_{1} - 2 \, x_{2} + e^{-x_{2}})^{2}$$i) Utilizar el método gradiente descendente que parte de un $x_{0}$ al azar entre $[-1,1]$
ii) Utilizar el método steepest descent que utiliza información sobre el Gradiente y el Hessiano de la función $F$ para calcular el parámetro $\alpha$(learning rate)
iii) ¿Que algoritmo converge más rápido a la solución? ¿Cuales son los pros y contras de cada método?
Las definiciones básicas que necesitamos recordar son:
En otras palabras si nos movemos del punto en cualquier dirección una distancia pequeña la función tomará valores mayores
Mínimo Global: El punto $\mathbf{x}^{*}$ es un único mínimo global de $F(\mathbf{x})$ si $F(\mathbf{\mathbf{x}^{*}})< F(\mathbf{\mathbf{x}^{*}}+ \Delta \mathbf{x})$ para todo $\Delta \mathbf{x} \neq 0$
Mínimo débil: El punto $\mathbf{x}^{*}$ es un mínimo débil de $F(\mathbf{x})$ si no es un mínimo fuerte y existe un escalar $\delta > 0$ tal que $F(\mathbf{\mathbf{x}^{*}}) < F(\mathbf{\mathbf{x}^{*}}+ \Delta \mathbf{x})$ para todo $\Delta \mathbf{x}$ tal que $\delta > \| \Delta \mathbf{x} \| > 0$
Se demuestra a traves del desarrollo de series de Taylor que en un entorno reducido los terminos de orden superior no son importantes y por ello podemos resumir lo anterior en las dos condiciones siguientes:
El punto $\mathbf{x}^{*}$ es un candidato a mínimo si:
$$\nabla F (\mathbf{x}) |_{\mathbf{x}= \mathbf{x}^{*} } = 0$$En este caso suponemos que se cumple la primer condición y del análisis de la serie de Taylor alrededor del punto candidato queda:
$$\Delta \mathbf{x} \, \nabla^{2} F(\mathbf{x})|_{\mathbf{x} = \mathbf{x}^{*}} \Delta \mathbf{x} > 0 $$Para que lo anterior se cumpla para un $\Delta \mathbf{x} \neq 0 $ se requiere que la matriz Hessiana sea definida positiva
Los algoritmos de optimización que vamos a utilizar para hallar el/los minimo/s de $F(x)$ se basan en comenzar con algún iterante inicial $\mathbf{x}_{0}$ y actualizar de acuerdo a alguna ecuación de la forma:
$$\mathbf{x}_{k+1}=\mathbf{x}_{k} + \alpha_{k} \mathbf{p}_{k}$$o $$\Delta \mathbf{x}_{k}=(\mathbf{x}_{k+1}-\mathbf{x}_{k})=\alpha_{k}\mathbf{p}_{k}$$
Donde el vector $\mathbf{p}_{k}$ representa una dirección de búsqueda y $\alpha_{k}$ es lo que se conoce como learning rate. Las diferentes estrategias para elegir a estos parámetros definen a cada algoritmo
En este caso se elige a $\mathbf{p}_{k}$ con la dirección opuesta a la del gradiente o sea $\mathbf{p}_{k}= -\mathbf{g}_{k} = - \nabla F (\mathbf{x}) |_{\mathbf{x}= \mathbf{x}_{k} } $
En este caso se minimiza con respecto a $\alpha_{k}$ en cada iteración, o sea elegir $\alpha_{k}$ para minimizar $F(\mathbf{x}_{k}+\alpha_{k} \mathbf{p}_{k})$ Para hacer esto en una función arbitraria se requiere:
$$\alpha_{k}=\frac{-\mathbf{g}_{k}\mathbf{p}_{k}}{\mathbf{p}_{k} \, A_{k} \, \mathbf{p}_{k}}$$Donde la Matriz $A_{k} = \nabla^{2} F (\mathbf{x}) |_{\mathbf{x}= \mathbf{x}_{k} } $ (Que en el caso que $F$ sea una forma cuádratica es constante)
A continuación las implementaciones con el maravilloso y libre lenguaje de Programación Python :P
In [21]:
# Imports and setups
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
%matplotlib inline
sym.init_printing() # Imprime la salida en LaTeX
plt.rcParams['figure.figsize'] = 8,6 #parámetros de tamaño
In [22]:
x_1, x_2 = sym.symbols('x_1 x_2')
In [23]:
# Funcion simbolica a minimizar
f = (x_2 - 2 * x_1 + sym.exp(-x_1))** 2 + (x_1 - 2 * x_2 + sym.exp(-x_2))** 2
f
Out[23]:
In [24]:
# Funcion numerica a minimizar
f_numeric = sym.lambdify((x_1,x_2),f,'numpy')
In [25]:
# Definimos el gradiente simbolico
def grad(f):
return sym.Matrix([f.diff(x_1),f.diff(x_2)])
In [26]:
#definimos el hessiano simbolico
def hess(f):
return sym.Matrix([[f.diff(x_1,2),f.diff(x_1,x_2)],[f.diff(x_2,x_1),f.diff(x_2,2)]])
In [27]:
g = grad(f)
h = hess(f)
h
Out[27]:
In [28]:
g_numeric = sym.lambdify((x_1,x_2),g,'numpy') # Gradiente numerico
h_numeric = sym.lambdify((x_1,x_2),h,'numpy') # Hessiano numerico
In [29]:
%%timeit
# gradiente descendente
ite1_max = 100
x_sol = np.empty((2,ite1_max))
x_0 = np.random.uniform(-1,1,size=2)
alpha = .03 # Learning rate
x_sol[:,0] = x_0 # Condicion inicial
tol1 = 1e-7 # tolerancia al error
delta1 = 1
i = 0
while (delta1 > tol1) and (i < (ite1_max - 1)):
p = np.asarray(g_numeric(x_sol[0,i],x_sol[1,i])).reshape(2,)
x_sol[:,i+1] = x_sol[:,i] - alpha * p
delta1 = np.linalg.norm(x_sol[:,i+1] - x_sol[:,i])
i += 1
In [30]:
print 'cantidad de iteraciones: --> ', i
In [31]:
# Plots
X = np.arange(-2, 2, 0.1)
Y = np.arange(-2, 2, 0.1)
X, Y = np.meshgrid(X, Y)
Z = f_numeric(X, Y)
CS = plt.contour(X, Y, Z)
plt.clabel(CS, inline=1, fontsize=10)
plt.title('Curvas de Nivel', fontsize=20)
plt.plot(x_sol[0,0:i],x_sol[1,0:i],'--o')
plt.plot(x_0[0],x_0[1],'ro') # Punto inicial
plt.plot(x_sol[0,i],x_sol[1,i],'go') #Punto final
Out[31]:
In [32]:
# Verificamos la condición de primer orden
g = np.asarray(g_numeric(x_sol[0,i],x_sol[1,i])).reshape(2,)
In [33]:
# Es casi cero
g
Out[33]:
In [34]:
# Verificamos la condicion de segundo orden(definida positiva)
H = np.asarray(h_numeric(x_sol[0,i],x_sol[1,i])).reshape(2,2)
In [35]:
H
Out[35]:
In [36]:
# Verificamos que es un minimo fuerte ya que los autovalores del Hessiano en x* son positivos
np.linalg.eigvals(H)
Out[36]:
In [37]:
%%timeit
# Steepest descent
ite2_max = 100
x_sol2 = np.empty((2,ite2_max))
#x_02 = np.random.uniform(-1,1,size=2)
k = 0
x_sol2[:,0] = x_0
tol = 1e-7
delta2 = 1
while (delta2 > tol) and (k < (ite2_max-1)):
grad_k = np.asarray(g_numeric(x_sol2[0,k],x_sol2[1,k])).reshape(2,)
p_k = -grad_k
hess_k = np.asarray(h_numeric(x_sol2[0,k],x_sol2[1,k])).reshape(2,2)
alpha_k = -(np.dot(grad_k,p_k) ) / (np.dot(p_k,hess_k).dot(p_k))
x_sol2[:,k+1] = x_sol2[:,k] + alpha_k * p_k
delta2 = np.linalg.norm(x_sol2[:,k] - x_sol2[:,k-1])
k += 1
In [38]:
print 'Cantidad de iteraciones: --> ', k
In [39]:
X = np.arange(-2, 2, 0.1)
Y = np.arange(-2, 2, 0.1)
X, Y = np.meshgrid(X, Y)
Z = f_numeric(X, Y)
CS = plt.contour(X, Y, Z)
plt.clabel(CS, inline=1, fontsize=10)
plt.title('Curvas de Nivel', fontsize=20)
plt.plot(x_sol2[0,0:k],x_sol2[1,0:k],'--o')
plt.plot(x_0[0],x_0[1],'ro') # Punto inicial
plt.plot(x_sol2[0,k],x_sol2[1,k],'go') #Punto final
Out[39]:
Graficamos la función en el espacio para observar ese hermoso valle que posee alrededor del/los minimos
In [40]:
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
ax = fig.gca(projection='3d')
X = np.arange(-2, 2, 0.1)
Y = np.arange(-2, 2, 0.1)
X, Y = np.meshgrid(X, Y)
Z = f_numeric(X, Y)
#ax.scatter(x_sol2[0],x_sol2[1],'ro')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.YlGnBu_r,
linewidth=0, antialiased=False)
#ax.set_zlim(-1.01, 1.01)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
ax.view_init(elev=30., azim=30.)
fig.colorbar(surf, shrink=0.5, aspect=5)
Out[40]:
Como vemos para este caso el algoritmo que converge más rápido(por muy poco en cuanto a tiempo) es el de Steepest descent y además realiza la tarea en menos iteraciones(con la misma tolerancia). Los pros de cada método son:
Gradiente descendente:
Donde $\lambda_{max}$ es el autovalor más grande de la representación matricial de $F$ y esta relacionado con la curvatura de la función
Steepest descent:
Desventajas:
Gradiente descendente:
Steepest descent:
Si la función $F$ no es una forma cuadrática el cálculo de $\alpha_{k}$ se vuelve más costoso computacionalmente ya que por cada iteración debemos evaluar el Hessiano de $F$ ($\nabla^{2}F$)
Si la función $F$ posee valles como por ejemplo el famoso caso de la función de Rosenbrock( $F(x_{1},x_{2})=(1-x_{2})^{2}+100(x_{1}-x_{2}^{2})^{2}$ ) el método tardará mucho en converger