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.
Consideramos con la norma cuadrática: $P=\nabla ^2 f_o(x^*)$ el método steepest descent:
$$\Delta x_{\text{nsd}} = \text{argmin} \{ \nabla f_o(x)^Tv : ||v||_P \leq 1, \nabla f_o(x)^Tv < 0 \} $$se prueba que $\Delta x_{\text{sd}} = - P^{-1} \nabla f_o(x)$.
Ver 4.2.Algoritmos_para_optimizacion_sin_restricciones. Entonces:
Algoritmo de Newton (visto como un caso particular del método steepest descent para matrices Hessianas simétricas definidas positivas)
Dado un punto inicial $x$ en $\text{dom}f_o$
Repetir el siguiente bloque para $k=0,1,2,...$
- Calcular la Hessiana de $f_o$ y definir $P$ como $P=\nabla ^2 f_o(x^*)$ con $x^*$ óptimo de $f_o$.
- 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 \hat{x} = - \nabla \hat{f}_o(\hat{x}) = - \nabla f_o(P^{-1/2}\hat{x})$.
- Búsqueda de línea. Elegir un tamaño de paso $t > 0$.
- Hacer la actualización: $\hat{x} = \hat{x} + t\Delta \hat{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.
Obs: en el cambio de coordenadas dado por $\hat{x} = P^{1/2}x$ y por tanto $\hat{f}_o({\hat{x}}) = f_o(P^{-1/2}\hat{x})$ se tiene $\nabla ^2 \hat{f}_o(\hat{x}) = P^{-1/2} \nabla f_o^2(x) P^{-1/2}$. Por lo que si $P$ se toma igual a $\nabla^2 f_o (x^*)$ entonces $P^{-1/2} \nabla f^2(x) P^{-1/2} = I$ (recuérdese que $P$ se asume simétrica definida positiva).
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 [6]:
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 [7]:
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 [8]:
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 e^{(x_1+3x_2-0.1)}+e^{x_1-3x_2-0.1}+e^{-x_1-0.1}$$
In [9]:
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 [10]:
x_ast=np.array([-0.3465735941251887,0], dtype=float)
In [11]:
P=Hessian_approximation(fo,x_ast)
In [12]:
P
Out[12]:
In [13]:
def coord_transform_original_variable(x):
eigenvalues, eigenvectors = np.linalg.eig(P)
return eigenvectors@((eigenvalues**(-1/2))*(np.transpose(eigenvectors)@x))
In [14]:
def coord_transform(x):
eigenvalues, eigenvectors = np.linalg.eig(P)
return eigenvectors@((eigenvalues**(1/2))*(np.transpose(eigenvectors)@x))
In [15]:
fo_transf = lambda x: fo(coord_transform_original_variable(x))
In [16]:
p_ast=fo(x_ast)
In [17]:
x_ast = coord_transform(x_ast)
In [18]:
x_ast
Out[18]:
In [19]:
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 [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.shape
Out[24]:
In [25]:
x_plot[:,0]
Out[25]:
In [26]:
x_plot[:,-1]
Out[26]:
In [27]:
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 [28]:
plt.plot(np.arange(Err_plot.size-2),Err_plot[2:],'.-')
plt.ylabel('Error relativo entre: $\hat{f}_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Iteraciones del método descenso en gradiente para $\hat{f}_o$')
plt.grid()
plt.show()
In [29]:
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 [30]:
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)
Curvas de nivel de $f_o$ transformada
In [31]:
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 [32]:
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.6, .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 [33]:
np.linalg.cond(Hf_x_k)
Out[33]:
Obsérvese que la Hessiana con la transformación de $f_o$ está muy bien condicionada.
La dirección de descenso $\Delta x_{\text{nt}} = -\nabla ^2 f_o(x) \nabla f(x)$ para $x \in \text{dom}f_o$ se le nombra paso de Newton para $f_o$ en $x$.
Obs: si $\nabla ^2 f_o(x)$ es simétrica definida positiva en $\text{dom}f_o$ y $x$ no es óptimo entonces $\Delta x_{\text{nt}}$ es dirección de descenso pues: $\nabla f_o(x)^T \Delta x_\text{nt} = -\nabla f_o(x)^T \nabla^2f_o(x)^{-1} \nabla f(x) < 0 $ con $x \in \text{dom}f_o$ y $\nabla f(x) \neq 0$.
El paso de Newton puede obtenerse de diferentes formas:
1) Utilizando la aproximación de segundo orden por Taylor $\hat{f}_o$ de $f_o$ en $x$:
$$\hat{f}_o(x+v) \approx f_o(x) + \nabla f_o(x)^Tv + \frac{1}{2} v^T \nabla ^2 f_o(x) v$$la cual si $\nabla ^2 f_o(x)$ es definida positiva en $\text{dom}f_o$ entonces es una función cuadrática convexa (estricta) en $v$ cuyo mínimo es $\hat{v} = \Delta x_\text{nt}$:
En el dibujo anterior $\hat{f}$ es un modelo cuadrático, $\Delta x_{nt}$ es dirección de descenso de Newton y $x^*$ es el óptimo de $f$. Del punto $(x,f(x))$ nos debemos mover al punto $(x+\Delta x_{nt}, f(x + \Delta x_{nt}))$ para llegar al óptimo y el valor de $f$ decrece: $f(x+\Delta x_{nt}) < f(x)$.
Comentario: encontrar el mínimo de $f_o(x) + \nabla f_o(x)^Tv + \frac{1}{2} v^T \nabla ^2 f_o(x) v$ es un problema de mínimos cuadrados, ver 4.1.Optimizacion_numerica_y_machine_learning.
2) Como un caso particular del método steepest descent (visto al inicio de la nota).
3) Utilizando la aproximación lineal de la condición de optimalidad $\nabla f(x^*) = 0$. Para valores cercanos a $x$ se tiene:
$$\nabla \hat{f}(x+v) \approx \nabla f(x) + \nabla ^2 f_o(x) v = 0$$entonces $\Delta x_{\text{nt}}$ es la solución de esta ecuación y es lo que debe añadirse a $x$ para que la condición de linealidad se satisfaga. Si $x$ es cercano a $x^*$ entonces $x + \Delta x_{\text{nt}}$ es una buena aproximación a $x^*$:
Una cantidad que se utiliza como criterio de paro del algoritmo y resultados de convergencia del método, es el decremento de Newton para $f_o$ en $x$ definido como:
$$\lambda (x) = (\nabla f(x) ^T \nabla ^2 f(x)^{-1} \nabla f(x))^{1/2}.$$Esta cantidad tiene propiedades como son:
$\frac{1}{2} \lambda ^2 (x)$ estima $f(x)-p^*$.
$|| \nabla^2 f(x)^{-1} \nabla f(x)||_{\nabla ^2f(x)} = \left ( \nabla f(x)^T \nabla^2 f(x)^{-1} \nabla ^2 f(x) \nabla ^2 f(x)^{-1} \nabla f(x) \right )^{1/2} = \lambda(x) $ que indica que $\lambda(x)$ es la norma del paso de Newton en la norma cuadrática definida por la Hessiana.
En el método de búsqueda de línea por backtracking $-\lambda (x) ^2$ es la derivada direccional de $f$ en $x$ en la direccióin de $\Delta x_{\text{nt}}$:
Algoritmo de Newton para problemas de optimización sin restricciones
Dado un punto inicial $x$ en $\text{dom}f_o$ y una tolerancia $\epsilon >0$.
Repetir el siguiente bloque para $k=0,1,2,...$
- Calcular la dirección de descenso de Newton $\Delta x_{\text{nt}} = - \nabla ^2 f_o(x)^{-1} \nabla f_o(x)$ y el decremento de Newton al cuadrado: $\lambda^2(x)=\nabla f(x)^T \nabla ^2 f(x)^{-1} \nabla f(x)$.
- Criterio de paro: finalizar el método si $\frac{\lambda^2(x)}{2} \leq \epsilon$.
- Búsqueda de línea. Elegir un tamaño de paso $t > 0$ (usar el cálculo de $\lambda (x)$ del paso anterior).
- Hacer la actualización: $x = x + t\Delta x_{\text{nt}}$.
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 [35]:
def Newtons_method(f, x_0, tol,
tol_backtracking, x_ast=None, p_ast=None, maxiter=30):
'''
Newton's method to numerically approximate solution of min f.
Args:
f (lambda expression): definition of function f.
x_0 (numpy ndarray): initial point for Newton's method.
tol (float): tolerance that will halt method. Controls stopping criteria.
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)
Hfeval = Hessian_approximation(f,x)
normgf = np.linalg.norm(gfeval)
condHf= np.linalg.cond(Hfeval)
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
#Newton's direction and Newton's decrement
dir_Newton = np.linalg.solve(Hfeval, -gfeval)
dec_Newton = -gfeval.dot(dir_Newton)
print('I\tNormgf \tNewton Decrement\tError x_ast\tError p_ast\tline search\tCondHf')
print('{}\t{:0.2e}\t{:0.2e}\t{:0.2e}\t{:0.2e}\t{}\t\t{:0.2e}'.format(iteration,normgf,
dec_Newton,Err,
Err_plot_aux[iteration],"---",
condHf))
stopping_criteria = dec_Newton/2
iteration+=1
while(stopping_criteria>tol and iteration < maxiter):
der_direct = -dec_Newton
t = line_search_by_backtracking(f,dir_Newton,x,der_direct)
x = x + t*dir_Newton
feval = f(x)
gfeval = gradient_approximation(f,x)
Hfeval = Hessian_approximation(f,x)
normgf = np.linalg.norm(gfeval)
condHf= np.linalg.cond(Hfeval)
#Newton's direction and Newton's decrement
dir_Newton = np.linalg.solve(Hfeval, -gfeval)
dec_Newton = -gfeval.dot(dir_Newton)
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}\t{:0.2e}\t{:0.2e}'.format(iteration,normgf,
dec_Newton,Err,
Err_plot_aux[iteration],t,
condHf))
stopping_criteria = dec_Newton/2
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 e^{(x_1+3x_2-0.1)}+e^{x_1-3x_2-0.1}+e^{-x_1-0.1}$$
In [36]:
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 [37]:
x_ast=np.array([-0.3465735941251887,0], dtype=float)
In [38]:
x_ast
Out[38]:
In [39]:
p_ast=fo(x_ast)
In [40]:
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]=Newtons_method(fo, x_0,tol,
tol_backtracking, x_ast, p_ast, maxiter)
In [41]:
x
Out[41]:
In [42]:
total_of_iterations
Out[42]:
In [43]:
Err_plot.shape
Out[43]:
In [44]:
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 [45]:
x_plot.shape
Out[45]:
In [46]:
x_plot[:,0]
Out[46]:
In [47]:
x_plot[:,-1]
Out[47]:
In [48]:
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 Newton para $f_o$')
plt.grid()
plt.show()
In [51]:
plt.plot(np.arange(Err_plot.size-1),Err_plot[1:],'.-')
plt.ylabel('$f_o(x^k)-p^*$',size=12)
plt.xlabel('Iteraciones del método de Newton para $f_o$')
plt.grid()
plt.show()
2) El método de Newton es insensible 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$.
Método de Newton para caso $C=10$
In [52]:
cte=10
In [53]:
fo = lambda x: 1/2*(x[0]**2+cte*x[1]**2)
In [54]:
x_ast=np.array([0,0],dtype=float)
In [55]:
p_ast=fo(x_ast)
In [56]:
x_0 = np.array([0.5,0.5],dtype=float)
tol=1e-8
tol_backtracking=1e-14
maxiter=50
[x,total_of_iterations,Err_plot,x_plot]=Newtons_method(fo, x_0,
tol, tol_backtracking, x_ast, p_ast, maxiter)
In [57]:
x
Out[57]:
In [58]:
total_of_iterations
Out[58]:
In [59]:
Err_plot.shape
Out[59]:
In [60]:
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 [61]:
x_plot[:,0]
Out[61]:
In [62]:
x_plot.shape
Out[62]:
In [63]:
x_plot[:,-1]
Out[63]:
In [64]:
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 Newton para $f_o$ con $C=10$')
plt.grid()
plt.show()
In [65]:
plt.plot(np.arange(Err_plot.size-1),Err_plot[1:],'.-')
plt.ylabel('Error relativo entre: $f_o(x^k)$ y $p^*$',size=12)
plt.xlabel('Iteraciones del método de Newton para $f_o$ con $C=10$')
plt.grid()
plt.show()
Comentarios:
Tanto el paso de Newton como el decremento de Newton son invariantes ante transformaciones afines. Por ejemplo para el paso de Newton: si $T \in \mathbb{R}^{n \times n}$ es no singular y $\hat{f}_o: \mathbb{R}^n \rightarrow \mathbb{R}$ está definida como $\hat{f}_o(\hat{x}) = f_o(T^{-1}x)$ entonces se puede probar que $\Delta \hat{x}_{\text{nt}} = T^{-1} \Delta x_{\text{nt}}$ con $\Delta x_{\text{nt}}$ el paso de Newton para $f_o$ en $x$. Esto indica que los pasos de Newton están relacionados por la misma transformación lineal $T$, de hecho: $x + \Delta x_{\text{nt}} = T(\hat{x} + \Delta \hat{x}_{\text{nt}})$.
Contrástese el punto anterior con el método de descenso en gradiente que es afectado fuertemente por el cambio de coordenadas, ver 4.2.Descenso_en_gradiente.
Al igual que en el método de descenso en gradiente los valores de $\alpha$ y $\beta$ en el método de búsqueda de línea por backtracking tienen un efecto pequeño en el performance del método de Newton. Típicamente se utilizan $\alpha=0.01, .15$ y $\beta = 0.5$.
El método de Newton puede dividirse en dos etapas:
La etapa damped o guarded que se caracteriza por valores de búsqueda de línea por backtracking menores a $1$: $t < 1$.
La etapa pure que se caracteriza por valores de búsqueda de línea por backtracking iguales a $1$: $t =1$ . En esta etapa el método de Newton tiene una convergencia cuadrática. Una vez en esta etapa le toma no más de $6$ iteraciones al método para obtener una aproximación a $p^*$ muy buena (precisiones mayores o iguales a $8$ dígitos).
En la práctica se ha observado que el método de Newton escala bien con el tamaño del problema: ejemplos de su performance para problemas con variables en $\mathbb{R}^{10^4}$ muestran un ligero incremento en el número de iteraciones que en problemas en $\mathbb{R}^{10}$.
La gran desventaja del método de Newton radica en construir, almacenar a la Hessiana y resolver el sistema de ecuaciones lineales asociado que típicamente es costoso. Para estos casos se sugieren:
Referencias: