Notas para contenedor de docker:
Comando de docker para ejecución de la nota de forma local:
nota: cambiar <ruta a mi directorio> por la ruta de directorio que se desea mapear a /datos dentro del contenedor de docker.
docker run --rm -v <ruta a mi directorio>:/datos --name jupyterlab_numerical -p 8888:8888 -d palmoreck/jupyterlab_numerical:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_numericalDocumentación de la imagen de docker palmoreck/jupyterlab_numerical:1.1.0 en liga.
En esta nota se consideran resolver problemas de la forma:
$$\min f_o(x)$$con $f_o:\mathbb{R}^n \rightarrow \mathbb{R}$ fuertemente convexa y $f_o \in \mathcal{C}^2(\text{dom}f_o)$ para buscar óptimos locales. Además se asume que los puntos iniciales $x^{(0)}$ de los métodos iterativos están en $\text{dom}f_o$ y los conjuntos $f_o(x^{(0)})$-subnivel son conjuntos cerrados. Ver 1.4.Polinomios_de_Taylor_y_diferenciacion_numerica y 4.1.Optimizacion_numerica_y_machine_learning para definiciones utilizadas en esta nota.
También se asume que existe un punto óptimo $x^*$ por lo que el problema tiene solución y el valor óptimo se denota por $p^* = f_o(x^*) = \inf f_o(x)$
Las suposición que una función $f$ sea convexa asegura que una condición necesaria y suficiente para que $x^*$ sea óptimo es: $\nabla f(x^*) = 0$ la cual es en general es un conjunto de $n$ ecuaciones no lineales en $n$ variables y que resuelve el problema de optimización planteado al inicio.
Algoritmo de descenso en gradiente
Dado un punto inicial $x$ en $\text{dom}f_o$
Repetir el siguiente bloque para $k=0,1,2,...$
- Calcular la dirección de descenso $\Delta x = - \nabla f_o(x)$.
- Búsqueda de línea. Elegir un tamaño de paso $t > 0$.
- Hacer la actualización: $x = x + t\Delta x$.
hasta convergencia (satisfacer criterio de paro).
Nota: los pasos en el algoritmo anterior representan una guía para la implementación. Al describirse los pasos de un algoritmo no implica que se tengan que implementar uno a continuación del otro como se describe. Si una implementación respeta la lógica y al mismo método, entonces pueden seguirse los pasos de una forma distinta.
In [1]:
import math
import matplotlib.pyplot as plt
import numpy as np
In [2]:
def inc_index(vec,index,h):
'''
Auxiliary function for gradient and Hessian computation.
Args:
vec (array): numpy array.
index (int): index.
h (float): quantity that vec[index] will be increased.
Returns:
vec (array): numpy array vec with vec[index] increased by h.
'''
vec[index] +=h
return vec
In [3]:
def dec_index(vec,index,h=1):
'''
Auxiliary function for gradient and Hessian computation.
Args:
vec (array): numpy array.
index (int): index.
h (float): quantity that vec[index] will be decreased.
Returns:
vec (array): numpy array vec with vec[index] decreased by h.
'''
vec[index] -=h
return vec
In [4]:
def gradient_approximation(f,x,h=1e-8):
'''
Numerical approximation of gradient for function f using forward differences.
Args:
f (lambda expression): definition of function f.
x (array): numpy array that holds values where gradient will be computed.
h (float): step size for forward differences, tipically h=1e-8
Returns:
gf (array): numerical approximation to gradient of f.
'''
n = x.size
gf = np.zeros(n)
f_x = f(x)
for i in np.arange(n):
inc_index(x,i,h)
gf[i] = f(x) - f_x
dec_index(x,i,h)
return gf/h
In [5]:
def Hessian_approximation(f,x,h=1e-6):
'''
Numerical approximation of Hessian for function f using forward differences.
Args:
f (lambda expression): definition of function f.
x (array): numpy array that holds values where Hessian will be computed.
h (float): step size for forward differences, tipically h=1e-6
Returns:
Hf (array): numerical approximation to Hessian of f.
'''
n = x.size
Hf = np.zeros((n,n))
f_x = f(x)
for i in np.arange(n):
inc_index(x,i,h)
f_x_inc_in_i = f(x)
for j in np.arange(i,n):
inc_index(x,j,h)
f_x_inc_in_i_j = f(x)
dec_index(x,i,h)
f_x_inc_in_j = f(x)
dif = f_x_inc_in_i_j-f_x_inc_in_i-f_x_inc_in_j+f_x
Hf[i,j] = dif
if j != i:
Hf[j,i] = dif
dec_index(x,j,h)
inc_index(x,i,h)
dec_index(x,i,h)
return Hf/h**2
In [8]:
def line_search_by_backtracking(f,dir_desc,x,
der_direct, alpha=.15, beta=.5):
"""
Line search that sufficiently decreases f restricted to a ray in the direction dir_desc.
Args:
alpha (float): parameter in line search with backtracking, tipically .15
beta (float): parameter in line search with backtracking, tipically .5
f (lambda expression): definition of function f.
dir_desc (array): descent direction.
x (array): numpy array that holds values where line search will be performed.
der_direct (float): directional derivative of f.
Returns:
t (float): positive number for stepsize along dir_desc that sufficiently decreases f.
"""
t=1
if alpha > 1/2:
print('alpha must be less than or equal to 1/2')
t=-1
if beta>1:
print('beta must be less than 1')
t=-1;
if t!=-1:
eval1 = f(x+t*dir_desc)
eval2 = f(x) + alpha*t*der_direct
while eval1 > eval2:
t=beta*t
eval1=f(x+t*dir_desc)
eval2=f(x)+alpha*t*der_direct
return t
In [9]:
def compute_error(x_obj,x_approx):
'''
Relative or absolute error between x_obj and x_approx.
'''
if np.linalg.norm(x_obj) > np.nextafter(0,1):
Err=np.linalg.norm(x_obj-x_approx)/np.linalg.norm(x_obj)
else:
Err=np.linalg.norm(x_obj-x_approx)
return Err
In [10]:
def gradient_descent(f, x_0, tol,
tol_backtracking, x_ast=None, p_ast=None, maxiter=30):
'''
Method of gradient descent to numerically approximate solution of min f.
Args:
f (lambda expression): definition of function f.
x_0 (numpy ndarray): initial point for gradient descent method.
tol (float): tolerance that will halt method. Controls norm of gradient of f.
tol_backtracking (float): tolerance that will halt method. Controls value of line search by backtracking.
x_ast (numpy ndarray): solution of min f, now it's required that user knows the solution...
p_ast (float): value of f(x_ast), now it's required that user knows the solution...
maxiter (int): maximum number of iterations
Returns:
x (numpy ndarray): numpy array, approximation of x_ast.
iteration (int): number of iterations.
Err_plot (numpy ndarray): numpy array of absolute error between p_ast and f(x) with x approximation
of x_ast. Useful for plotting.
x_plot (numpy ndarray): numpy array that containts in columns vector of approximations. Last column
contains x, approximation of solution. Useful for plotting.
'''
iteration = 0
x = x_0
feval = f(x)
gfeval = gradient_approximation(f,x)
normgf = np.linalg.norm(gfeval)
Err_plot_aux = np.zeros(maxiter)
Err_plot_aux[iteration]=compute_error(p_ast,feval)
Err = compute_error(x_ast,x)
n = x.size
x_plot = np.zeros((n,maxiter))
x_plot[:,iteration] = x
print('I\tNormagf\t\tError x_ast\tError p_ast\tline search')
print('{}\t{:0.2e}\t{:0.2e}\t{:0.2e}\t{}'.format(iteration,normgf,Err,Err_plot_aux[iteration],"---"))
iteration+=1
while(normgf>tol and iteration < maxiter):
dir_desc = -gfeval
der_direct = gfeval.dot(dir_desc)
t = line_search_by_backtracking(f,dir_desc,x,der_direct)
x = x + t*dir_desc
feval = f(x)
gfeval = gradient_approximation(f,x)
normgf = np.linalg.norm(gfeval)
Err_plot_aux[iteration]=compute_error(p_ast,feval)
x_plot[:,iteration] = x
Err = compute_error(x_ast,x)
print('{}\t{:0.2e}\t{:0.2e}\t{:0.2e}\t{:0.2e}'.format(iteration,normgf,Err,
Err_plot_aux[iteration],t))
if t<tol_backtracking: #if t is less than tol_backtracking then we need to check the reason
iter_salida=iteration
iteration = maxiter - 1
iteration+=1
print('{} {:0.2e}'.format("Error of x with respect to x_ast:",Err))
print('{} {}'.format("Approximate solution:", x))
cond = Err_plot_aux > np.finfo(float).eps*10**(-2)
Err_plot = Err_plot_aux[cond]
if iteration == maxiter and t < tol_backtracking:
print("Backtracking value less than tol_backtracking, check approximation")
iteration=iter_salida
else:
if iteration == maxiter:
print("Reached maximum of iterations, check approximation")
x_plot = x_plot[:,:iteration]
return [x,iteration,Err_plot,x_plot]
1) $$\min \quad (x_1-2)^2+(2-x_2)^2+x_3^2+x_4^4$$
La solución del problema anterior es $x^* = (2,2,0,0)^T$.
Método de descenso gradiente con $x_0=(5,5,1,0)^T$:
In [11]:
fo = lambda x: (x[0]-2)**2 + (2-x[1])**2 + x[2]**2 + x[3]**4
x_ast = np.array([2,2,0,0],dtype=float)
x_0 = np.array([5,5,1,0],dtype=float)
tol=1e-8
tol_backtracking=1e-14
maxiter=5
p_ast=fo(x_ast)
[x,total_of_iterations,Err_plot,x_plot]=gradient_descent(fo, x_0, tol, tol_backtracking, x_ast, p_ast, maxiter)
In [12]:
x
Out[12]:
In [13]:
total_of_iterations
Out[13]:
In [14]:
Err_plot
Out[14]:
In [15]:
x_plot
Out[15]:
In [16]:
plt.plot(np.arange(Err_plot.size),Err_plot,'o-')
plt.ylabel('Error relativo entre: $f_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Iterations',size=12)
plt.grid()
plt.show()
Comentario: en esta gráfica se observa una convergencia lineal del método de descenso en gradiente para $f_o$.
2) El método de descenso en gradiente es altamente sensible a la forma de las curvas de nivel de la función objetivo $f_o$. Para observar esto considérese el problema: $$\min \frac{1}{2}\left(x_1^2+Cx_2^2 \right)$$
La solución del problema anterior es $x^*=(0,0)^T$.
Método de descenso en gradiente con $x_0=(0.5,0.5)^T$.
Descenso en gradiente para caso $C=10$
In [17]:
cte=10
In [18]:
fo = lambda x: 1/2*(x[0]**2+cte*x[1]**2)
In [19]:
x_ast=np.array([0,0],dtype=float)
x_0 = np.array([0.5,0.5],dtype=float)
tol=1e-8
tol_backtracking=1e-14
maxiter=50
p_ast=fo(x_ast)
[x,total_of_iterations,Err_plot,x_plot]=gradient_descent(fo, x_0, tol, tol_backtracking, x_ast, p_ast, maxiter)
In [20]:
x
Out[20]:
In [21]:
total_of_iterations
Out[21]:
In [22]:
Err_plot.shape
Out[22]:
In [23]:
plt.yscale('log') #escala logarítmica en el eje vertical
plt.plot(np.arange(Err_plot.size),Err_plot,'.-')
plt.ylabel('Log Error relativo entre: $f_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Iterations',size=12)
plt.grid()
plt.show()
In [24]:
x_plot[:,0]
Out[24]:
In [25]:
x_plot.shape
Out[25]:
In [26]:
x_plot[:,-1]
Out[26]:
In [27]:
plt.figure(figsize=(8,4))
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.ylabel('y')
plt.xlabel('x')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12,xytext=(0.5, 0.47))
plt.title('Iteraciones del método de descenso en gradiente para $f_o$ con $C=10$')
plt.grid()
plt.show()
In [28]:
plt.plot(np.arange(Err_plot.size-40),Err_plot[40:],'.-')
plt.ylabel('Error relativo entre: $f_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Últimas 10 iteraciones del método descenso en gradiente para $f_o$ con $C=10$')
plt.grid()
plt.show()
Comentario: en esta gráfica se observa una convergencia lineal del método de descenso en gradiente para $f_o$.
Curvas de nivel de $f_o$
In [29]:
z=lambda x_mesh,y_mesh,cte: 1/2*(x_mesh**2+cte*y_mesh**2)
In [30]:
density=1e-1
xl=-1
yl=-1
xr=1
yr=1
x_p=np.arange(xl,xr,density)
y_p=np.arange(yl,yr,density)
x_mesh,y_mesh = np.meshgrid(x_p,y_p)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.contour(x_p,y_p,z(x_mesh,y_mesh,cte))
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12)
plt.grid()
plt.title("Curvas de nivel de $f_o$, $C=10$")
plt.show()
Curvas de nivel $C=1$
In [31]:
cte=1
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.contour(x_p,y_p,z(x_mesh,y_mesh,cte))
plt.grid()
plt.title("Curvas de nivel de $f_o$, $C=1$")
plt.show()
Curvas de nivel $C=0.1$
In [32]:
cte=.1
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.contour(x_p,y_p,z(x_mesh,y_mesh,cte))
plt.grid()
plt.title("Curvas de nivel de $f_o$, $C=0.1$")
plt.show()
Se puede probar que las curvas de nivel para $f_o$ bajo las suposiciones que están al inicio de esta nota tienen forma de elipsoides cerca del óptimo $x^*$ y son aproximadamente iguales al elipsoide definido como:
$$\mathcal{E}_{x^*} = \{z \in \mathbb{R}^2 : ||z-x^*||_{\nabla^2 f_o(x^*)}^2 \leq 1\}$$*La norma $||z||_P$ se nombra norma cuadrática y está definida para matrices simétricas positivas definidas $P$. La definición es: $$||z||_P = \sqrt{z^TPz}.$$
Obsérvese que $f_o(x)=\frac{1}{2}\left(x_1^2+Cx_2^2 \right)$ tiene una Hessiana igual a: $\nabla^2f(x)=\left [\begin{array}{cc} 1 & 0\\ 0 & C\\ \end{array} \right] $ $\forall x \in \mathbb{R}^2$.
A continuación se presenta el elipsoide $\mathcal{E}_{x^{(k)}}=\{z \in \mathbb{R}^2 : ||z-x^{(k)}||_{\nabla^2 f_o(x^{(k)})}^2 \leq 1\}$ para la última iteración $x^{(k)}$ de los resultados del método de descenso en gradiente anteriores usando $C=10$:
In [33]:
cte=10
In [34]:
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.contour(x_p,y_p,z(x_mesh,y_mesh,cte))
plt.grid()
plt.title("Curvas de nivel de $f_o$, $C=10$")
plt.show()
In [35]:
x_k = x_plot[:,-1].reshape(2,1)
Hf_x_k = Hessian_approximation(fo,x_k)
norm_cuad = lambda z : np.sum((z-x_k)*(Hf_x_k@(z-x_k)),axis=0)
density_p=int(2.5*10**5)
z_p=np.random.uniform(-2,2,(2,density_p))
lim_sup=2
lim_inf=1.95
ind=(norm_cuad(z_p)<=lim_sup) & (lim_inf<=norm_cuad(z_p))
z_p_subset = z_p[:,ind]
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.scatter(z_p_subset[0,:],z_p_subset[1,:],marker='.')
plt.title('Puntos en el plano que cumplen $||x|| \leq 1$ para norma cuadrática usando la segunda derivada de $f_o$ evaluada en la última iteración')
plt.grid()
plt.show()
Número de condición de la matriz Hessiana
In [36]:
Hf_x_k
Out[36]:
In [37]:
np.linalg.cond(Hf_x_k)
Out[37]:
La Hessiana es medianamente bien condicionada
Descenso en gradiente para caso $C=1$
In [38]:
cte=1
In [39]:
x_ast=np.array([0,0],dtype=float)
x_0 = np.array([0.5,0.5],dtype=float)
tol=1e-8
tol_backtracking=1e-14
maxiter=50
p_ast=fo(x_ast)
[x,total_of_iterations,Err_plot,x_plot]=gradient_descent(fo, x_0, tol, tol_backtracking, x_ast, p_ast, maxiter)
In [40]:
x
Out[40]:
In [41]:
total_of_iterations
Out[41]:
In [42]:
plt.plot(np.arange(Err_plot.size),Err_plot,'.-')
plt.ylabel('Error relativo entre: $f_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Iterations',size=12)
plt.grid()
plt.show()
Comentario: en esta gráfica se observa una convergencia lineal del método de descenso en gradiente para $f_o$ con $C=1$.
In [43]:
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.contour(x_p,y_p,z(x_mesh,y_mesh,cte))
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12)
plt.title("Curvas de nivel de $f_o$, $C=1$")
plt.grid()
plt.show()
In [44]:
x_k = x_plot[:,-1].reshape(2,1)
Hf_x_k = Hessian_approximation(fo,x_k)
norm_cuad = lambda z : np.sum((z-x_k)*(Hf_x_k@(z-x_k)),axis=0)
density_p=int(2.5*10**4)
z_p=np.random.uniform(-1,1,(2,density_p))
lim_sup=1
lim_inf=.95
ind=(norm_cuad(z_p)<=lim_sup) & (lim_inf<=norm_cuad(z_p))
z_p_subset = z_p[:,ind]
plt.xlim(-1.1, 1.1)
plt.ylim(-1.1, 1.1)
plt.scatter(z_p_subset[0,:],z_p_subset[1,:],marker='.')
plt.title('Puntos en el plano que cumplen $||x|| \leq 1$ para norma cuadrática usando la segunda derivada de $f_o$ evaluada en la última iteración')
plt.grid()
plt.show()
In [45]:
np.linalg.cond(Hf_x_k)
Out[45]:
La matriz Hessiana es bien condicionada
Descenso en gradiente para caso $C=0.1$
In [46]:
cte=.1
In [47]:
x_ast=np.array([0,0],dtype=float)
x_0 = np.array([0.5,0.5],dtype=float)
tol=1e-8
tol_backtracking=1e-14
maxiter=70
p_ast=fo(x_ast)
[x,total_of_iterations,Err_plot,x_plot]=gradient_descent(fo, x_0, tol, tol_backtracking, x_ast, p_ast, maxiter)
In [48]:
x
Out[48]:
In [49]:
total_of_iterations
Out[49]:
In [50]:
Err_plot.shape
Out[50]:
In [51]:
plt.yscale('log') #escala logarítmica en el eje vertical
plt.plot(np.arange(Err_plot.size),Err_plot,'.-')
plt.ylabel('Log Error relativo entre: $f_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Iterations',size=12)
plt.grid()
plt.show()
In [52]:
x_plot[:,-1]
Out[52]:
In [53]:
plt.figure(figsize=(8,4))
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.ylabel('y')
plt.xlabel('x')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12,xytext=(0.5, 0.47))
plt.title('Iteraciones del método de descenso en gradiente para $f_o$ con $C=10$')
plt.grid()
plt.show()
In [54]:
plt.plot(np.arange(Err_plot.size-65),Err_plot[65:],'.-')
plt.ylabel('Error relativo entre: $f_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Últimas 5 iteraciones del método descenso en gradiente para $f_o$ con $C=0.1$')
plt.grid()
plt.show()
Comentario: en esta gráfica se observa una convergencia lineal del método de descenso en gradiente para $f_o$.
In [55]:
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.contour(x_p,y_p,z(x_mesh,y_mesh,cte))
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12)
plt.title("Curvas de nivel de $f_o$, $C=0.1$")
plt.grid()
plt.show()
In [56]:
x_k = x_plot[:,-1].reshape(2,1)
Hf_x_k = Hessian_approximation(fo,x_k)
norm_cuad = lambda z : np.sum((z-x_k)*(Hf_x_k@(z-x_k)),axis=0)
density_p=int(2.5*10**4)
z_p=np.random.uniform(-1,1,(2,density_p))
lim_sup=.5
lim_inf=.45
ind=(norm_cuad(z_p)<=lim_sup) & (lim_inf<=norm_cuad(z_p))
z_p_subset = z_p[:,ind]
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.scatter(z_p_subset[0,:],z_p_subset[1,:],marker='.')
plt.title('Puntos en el plano que cumplen $||x|| \leq 1$ para norma cuadrática usando la segunda derivada de $f_o$ evaluada en la última iteración')
plt.grid()
plt.show()
In [57]:
np.linalg.cond(Hf_x_k)
Out[57]:
La matriz Hessiana es medianamente bien condicionada
Comentarios:
El número de condición de la matriz Hessiana de la función $f_o(x) = \frac{1}{2}\left ( x_1^2 + Cx_2^2 \right)$ es igual a $\max\left\{C, \frac{1}{C}\right\}$ lo cual se obtiene del cociente del eigenvalor máximo entre el eigenvalor mínimo de $\nabla ^2f_o(x)$: $\frac{\max\{1,C\}}{\min\{1,C\}}$ ya que $\nabla ^2 f_o(x)$ es una matriz diagonal y positiva definida pues $C>0$: $\nabla^2f_o(x)=\left [\begin{array}{cc} 1 & 0\\ 0 & C\\ \end{array} \right] $ $\forall x \in \mathbb{R}^2$.
Por el punto anterior el número de condición de la matriz Hessiana es grande si $C >> 1$ o $C <<1$.
El número de condición de la matriz Hessiana se observa que influye mucho en la convergencia del método de descenso en gradiente. Aún siendo la Hessiana medianamente bien condicionadas (número de condición igual a 10 con $C=10$ o $C=0.1$) el número de iteraciones cambió drásticamente (un número de iteraciones mayor a $50$). Esto se cumple para este método y el método más general: el método de descenso más pronunciado: steepest descent method bajo la norma cuadrática con matrices $P$ que no aproximen bien a la Hessiana de $f_o$. Para matrices mal condicionadas el método de gradiente es muy lento por lo que no se recomienda en la práctica.
La simplicidad del método de descenso en gradiente es su mayor ventaja pero su enorme dependencia en el número de condición en la Hessiana o los conjuntos subnivel es su desventaja.
La elección de valores en los parámetros de búsqueda de línea por backtracking $\alpha$ y $\beta$ tienen un efecto visible pero no dramático en la convergencia. Típicamente se utilizan $\alpha=0.01, 0.15$ y $\beta = 0.5$.
Ejercicio: realizar mismas gráficas y análisis de convergencia y curvas de nivel para la función objetivo $f_o(x) = \frac{1}{2}\left ( Cx_1^2 + x_2^2 \right)$ para $C=10^{-2}, 10^2$.
Algoritmo de descenso en gradiente (visto como un caso particular del método steepest descent para matrices simétricas definidas positivas)
Dado un punto inicial $x$ en $\text{dom}f_o$
Repetir el siguiente bloque para $k=0,1,2,...$
- Definir $P$ matriz simétrica definida positiva.
- Transformar $\hat{x}$ con: $x = P^{-1/2} \hat{x}$ con $\hat{x}$ variable del problema transformado*.
- Calcular la dirección de descenso en gradiente $\Delta x = - \nabla f_o(x)$.
- Búsqueda de línea. Elegir un tamaño de paso $t > 0$.
- Hacer la actualización: $x = x + t\Delta x$.
hasta convergencia (satisfacer criterio de paro).
*El problema transformado bajo $P$ es:
$$min \hat{f}_o(\hat{x})$$con $\hat{f}_o (\hat{x}) = f_o(P^{-1/2}\hat{x}) = f_o(x)$ y $x$ variable "original": $\hat{x} = P^{1/2}x$.
Entonces aplicamos descenso en gradiente a $f_o$ pero transformando los datos de entrada de $f_o$ por la matriz raíz cuadrada simétrica* $P^{-1/2}$.
*$P^{1/2}$ se nombra raíz cuadrada simétrica o symmetric squareroot y está definida para matrices $P$ simétricas semidefinidas positivas como $P^{1/2}=Qdiag(\lambda_1^{1/2},\dots,\lambda_n^{1/2})Q^T$ con $Q$ y $diag(\lambda_1^{1/2},\dots,\lambda_n^{1/2})$ obtenidas con la descomposición espectral de $P$, ver 3.3.d.SVD.
1) $$\min \quad e^{(x_1+3x_2-0.1)}+e^{x_1-3x_2-0.1}+e^{-x_1-0.1}$$
a)Usando $P_1=\left [\begin{array}{cc} 1 & 0\\ 0 & 1\\ \end{array} \right] $ por lo que no hay transformación de coordenadas.
In [58]:
P=np.eye(2)
In [59]:
P
Out[59]:
In [60]:
def coord_transform_original_variable(x):
eigenvalues, eigenvectors = np.linalg.eig(P)
return eigenvectors@((eigenvalues**(-1/2))*(np.transpose(eigenvectors)@x))
In [61]:
def coord_transform(x):
eigenvalues, eigenvectors = np.linalg.eig(P)
return eigenvectors@((eigenvalues**(1/2))*(np.transpose(eigenvectors)@x))
In [62]:
fo = lambda x: math.exp(x[0]+3*x[1]-0.1)+math.exp(x[0]-3*x[1]-0.1)+math.exp(-x[0]-0.1)
In [63]:
fo_transf = lambda x: fo(coord_transform_original_variable(x))
In [64]:
x_ast=np.array([-0.3465735941251887,0], dtype=float)
In [65]:
p_ast=fo(x_ast)
In [66]:
x_ast = coord_transform(x_ast)
In [67]:
x_ast
Out[67]:
In [68]:
x_0 = np.array([0,0],dtype=float)
tol=1e-8
tol_backtracking=1e-14
maxiter=50
[x,total_of_iterations,Err_plot,x_plot]=gradient_descent(fo_transf, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
In [69]:
x
Out[69]:
In [70]:
total_of_iterations
Out[70]:
In [71]:
Err_plot.shape
Out[71]:
In [72]:
plt.yscale('log') #escala logarítmica en el eje vertical
plt.plot(np.arange(Err_plot.size),Err_plot,'.-')
plt.ylabel('Log Error relativo entre: $f_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Iterations',size=12)
plt.grid()
plt.show()
In [73]:
x_plot.shape
Out[73]:
In [74]:
x_plot[:,0]
Out[74]:
In [75]:
x_plot[:,-1]
Out[75]:
In [76]:
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.ylabel('y')
plt.xlabel('x')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12)
plt.title('Iteraciones del método de descenso en gradiente para $f_o$')
plt.grid()
plt.show()
In [77]:
plt.plot(np.arange(Err_plot.size-22),Err_plot[22:],'.-')
plt.ylabel('Error relativo entre: $f_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Últimas 7 iteraciones del método descenso en gradiente para $f_o$')
plt.grid()
plt.show()
Comentario: en esta gráfica se observa una convergencia lineal del método de descenso en gradiente para $f_o$.
Curvas de nivel de $f_o$
In [78]:
z=lambda x_mesh,y_mesh: np.exp(x_mesh+3*y_mesh-0.1)+np.exp(x_mesh-3*y_mesh-0.1)+np.exp(-x_mesh-0.1)
In [79]:
def coord_transform_original_variable_2(x):
eigenvalues, eigenvectors = np.linalg.eig(P)
return eigenvectors@((eigenvalues**(-1/2)).reshape(x.shape[0],1)*(eigenvectors@x))
In [80]:
density=1e-1
xl=-.5
yl=-.5
xr=.5
yr=.5
x_p=np.arange(xl,xr,density)
y_p=np.arange(yl,yr,density)
x_y_p_t = coord_transform_original_variable_2(np.row_stack((x_p,y_p)))
x_p_t = x_y_p_t[0,:]
y_p_t = x_y_p_t[1,:]
x_mesh,y_mesh = np.meshgrid(x_p_t,y_p_t)
plt.xlim(-.5, .5)
plt.ylim(-.5, .5)
plt.contour(x_p,y_p,z(x_mesh,y_mesh))
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12)
plt.grid()
plt.title("Curvas de nivel de $f_o$")
plt.show()
In [81]:
x_k = x_plot[:,-1]
Hf_x_k = Hessian_approximation(fo_transf,x_k)
x_k = x_k.reshape(2,1)
norm_cuad = lambda z : np.sum((z-x_k)*(Hf_x_k@(z-x_k)),axis=0)
density_p=int(2.5*10**4)
z_p=np.random.uniform(-1,1,(2,density_p))
lim_sup=.9
lim_inf=.85
ind=(norm_cuad(z_p)<=lim_sup) & (lim_inf<=norm_cuad(z_p))
z_p_subset = z_p[:,ind]
plt.xlim(-.5, .5)
plt.ylim(-.5, .5)
plt.scatter(z_p_subset[0,:],z_p_subset[1,:],marker='.')
plt.title('Puntos en el plano que cumplen $||x|| \leq 1$ para norma cuadrática usando la segunda derivada de $f_o$ evaluada en la última iteración')
plt.grid()
plt.show()
In [82]:
np.linalg.cond(Hf_x_k)
Out[82]:
La matriz Hessiana está bien condicionada.
En los siguientes ejemplos continuamos usando el método de descenso en gradiente pero con diferentes matrices de transformación para el método de steepest descent bajo la norma cuadrática y ver su influencia en las curvas de nivel y en el número de condición de la Hessiana.
3) Usando $P_1=\left [\begin{array}{cc} 2 & 0\\ 0 & 8\\ \end{array} \right] $ para transformar a $f_0$.
In [83]:
P=np.array([[2,0],
[0,8]])
In [84]:
P
Out[84]:
In [85]:
x_ast=np.array([-0.3465735941251887,0], dtype=float)
In [86]:
p_ast=fo(x_ast)
In [87]:
x_ast = coord_transform(x_ast)
In [88]:
x_ast
Out[88]:
In [89]:
x_0 = np.array([0,0],dtype=float)
tol=1e-8
tol_backtracking=1e-14
maxiter=50
[x,total_of_iterations,Err_plot,x_plot]=gradient_descent(fo_transf, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
In [90]:
x
Out[90]:
In [91]:
total_of_iterations
Out[91]:
In [92]:
Err_plot.shape
Out[92]:
In [93]:
plt.yscale('log') #escala logarítmica en el eje vertical
plt.plot(np.arange(Err_plot.size),Err_plot,'.-')
plt.ylabel('Log Error relativo entre: $\hat{f}_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Iterations',size=12)
plt.grid()
plt.show()
In [94]:
x_plot.shape
Out[94]:
In [95]:
x_plot[:,0]
Out[95]:
In [96]:
x_plot[:,-1]
Out[96]:
In [97]:
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.ylabel('y')
plt.xlabel('x')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12)
plt.title('Iteraciones del método de descenso en gradiente para $\hat{f}_o$')
plt.grid()
plt.show()
In [98]:
plt.plot(np.arange(Err_plot.size-10),Err_plot[10:],'.-')
plt.ylabel('Error relativo entre: $\hat{f}_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Últimas 7 iteraciones del método descenso en gradiente para $\hat{f}_o$')
plt.grid()
plt.show()
Comentario: en esta gráfica se observa una convergencia lineal del método de descenso en gradiente para $f_o$ transformada.
Curvas de nivel de $f_o$ transformada
In [99]:
density=1e-1
xl=-1.5
yl=-1
xr=.5
yr=1
x_p=np.arange(xl,xr,density)
y_p=np.arange(yl,yr,density)
x_y_p_t = coord_transform_original_variable_2(np.row_stack((x_p,y_p)))
x_p_t = x_y_p_t[0,:]
y_p_t = x_y_p_t[1,:]
x_mesh,y_mesh = np.meshgrid(x_p_t,y_p_t)
plt.xlim(-1.5, .5)
plt.ylim(-1, 1)
plt.contour(x_p,y_p,z(x_mesh,y_mesh))
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12)
plt.grid()
plt.title("Curvas de nivel de $\hat{f}_o$")
plt.show()
In [100]:
x_k = x_plot[:,-1]
Hf_x_k = Hessian_approximation(fo_transf,x_k)
x_k = x_k.reshape(2,1)
norm_cuad = lambda z : np.sum((z-x_k)*(Hf_x_k@(z-x_k)),axis=0)
density_p=int(2.5*10**4)
z_p=np.random.uniform(-2,1,(2,density_p))
lim_sup=.9
lim_inf=.85
ind=(norm_cuad(z_p)<=lim_sup) & (lim_inf<=norm_cuad(z_p))
z_p_subset = z_p[:,ind]
plt.xlim(-1.5, .5)
plt.ylim(-1, 1)
plt.scatter(z_p_subset[0,:],z_p_subset[1,:],marker='.')
plt.title('Puntos en el plano que cumplen $||x|| \leq 1$ para norma cuadrática usando la segunda derivada de $\hat{f}_o$ evaluada en la última iteración')
plt.grid()
plt.show()
In [101]:
np.linalg.cond(Hf_x_k)
Out[101]:
Obsérvese que es menor el número de condición con la transformación de $f_o$.
c) Usando $P_2=\left [\begin{array}{cc} 8 & 0\\ 0 & 2\\ \end{array} \right] $ para transformar a $f_o$.
In [102]:
P=np.array([[8,0],
[0,2]])
In [103]:
P
Out[103]:
In [104]:
x_ast=np.array([-0.3465735941251887,0], dtype=float)
In [105]:
p_ast=fo(x_ast)
In [106]:
x_ast = coord_transform(x_ast)
In [107]:
x_ast
Out[107]:
In [108]:
x_0 = np.array([0,0],dtype=float)
tol=1e-8
tol_backtracking=1e-14
maxiter=60
[x,total_of_iterations,Err_plot,x_plot]=gradient_descent(fo_transf, x_0,tol, tol_backtracking, x_ast, p_ast, maxiter)
In [109]:
x
Out[109]:
In [110]:
total_of_iterations
Out[110]:
In [111]:
Err_plot.shape
Out[111]:
In [112]:
plt.yscale('log') #escala logarítmica en el eje vertical
plt.plot(np.arange(Err_plot.size),Err_plot,'.-')
plt.ylabel('Log Error relativo entre: $\hat{f}_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Iterations',size=12)
plt.grid()
plt.show()
In [113]:
x_plot.shape
Out[113]:
In [114]:
x_plot[:,0]
Out[114]:
In [115]:
x_plot[:,-1]
Out[115]:
In [116]:
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.ylabel('y')
plt.xlabel('x')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12)
plt.title('Iteraciones del método de descenso en gradiente para $\hat{f}_o$')
plt.grid()
plt.show()
In [117]:
plt.plot(np.arange(Err_plot.size-53),Err_plot[53:],'.-')
plt.ylabel('Error relativo entre: $\hat{f}_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Últimas 7 iteraciones del método descenso en gradiente para $\hat{f}_o$')
plt.grid()
plt.show()
Comentario: en esta gráfica se observa una convergencia lineal del método de descenso en gradiente para $f_o$ transformada.
Curvas de nivel de $f_o$ transformada
In [118]:
density=1e-1
xl=-1.5
yl=-1
xr=.5
yr=1
x_p=np.arange(xl,xr,density)
y_p=np.arange(yl,yr,density)
x_y_p_t = coord_transform_original_variable_2(np.row_stack((x_p,y_p)))
x_p_t = x_y_p_t[0,:]
y_p_t = x_y_p_t[1,:]
x_mesh,y_mesh = np.meshgrid(x_p_t,y_p_t)
plt.xlim(-1.5, .5)
plt.ylim(-1, 1)
plt.contour(x_p,y_p,z(x_mesh,y_mesh))
plt.plot(x_plot[0,:],x_plot[1,:],'-*')
plt.annotate('$x^{(0)}$',(x_plot[0,0],x_plot[1,0]),fontsize=12)
plt.grid()
plt.title("Curvas de nivel de $\hat{f}_o$")
plt.show()
In [119]:
x_k = x_plot[:,-1]
Hf_x_k = Hessian_approximation(fo_transf,x_k)
x_k = x_k.reshape(2,1)
norm_cuad = lambda z : np.sum((z-x_k)*(Hf_x_k@(z-x_k)),axis=0)
density_p=int(2.5*10**4)
z_p=np.random.uniform(-2,1,(2,density_p))
lim_sup=.9
lim_inf=.85
ind=(norm_cuad(z_p)<=lim_sup) & (lim_inf<=norm_cuad(z_p))
z_p_subset = z_p[:,ind]
plt.xlim(-1.5, .5)
plt.ylim(-1, 1)
plt.scatter(z_p_subset[0,:],z_p_subset[1,:],marker='.')
plt.title('Puntos en el plano que cumplen $||x|| \leq 1$ para norma cuadrática usando la segunda derivada de $\hat{f}_o$ evaluada en la última iteración')
plt.grid()
plt.show()
In [120]:
np.linalg.cond(Hf_x_k)
Out[120]:
Obsérvese que es mayor el número de condición con la transformación de $f_o$.
Comentarios:
Como se observa en los ejemplos anteriores la elección de la norma en el método de descenso más pronunciado tiene un efecto fuerte en la tasa de convergencia.
Siempre existe una matriz $P$ para la cual el método de descenso más pronunciado tiene una convergencia buena. El reto está en encontrar tal matriz. La idea es identificar una matriz $P$ para la cual el problema transformado tenga un número de condición moderado.
1) $$\min \quad e^{(x_1+3x_2-0.1)}+e^{x_1-3x_2-0.1}+e^{-x_1-0.1}$$
óptimo:
$x^* = \left [\begin{array}{c} -0.3465735941251887 \\ 0 \end{array} \right] $
In [122]:
!pip3 install --user -q cvxpy
In [123]:
import cvxpy as cp
In [124]:
x1 = cp.Variable()
x2 = cp.Variable()
In [125]:
obj = cp.Minimize(cp.exp(x1+3*x2-0.1)+cp.exp(x1-3*x2-0.1)+cp.exp(-x1-0.1))
In [126]:
prob = cp.Problem(obj)
In [127]:
prob.solve()
Out[127]:
In [128]:
print("status:", prob.status)
print("optimal value", prob.value)
print("optimal var", x1.value, x2.value)
Referencias: