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_r_kernel_local -p 8888:8888 -d palmoreck/jupyterlab_r_kernel:1.1.0password para jupyterlab: qwerty
Detener el contenedor de docker:
docker stop jupyterlab_r_kernel_localDocumentación de la imagen de docker palmoreck/jupyterlab_r_kernel: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.
Ver Coordinate descent.
Consideramos con la norma $1$ el método steepest descent:
$$\Delta x_{\text{nsd}} = \text{argmin} \{ \nabla f_o(x)^Tv : ||v||_1 \leq 1, \nabla f_o(x)^Tv < 0 \} $$se prueba que $\Delta x_{\text{sd}} = - \frac{\partial f_o(x)}{\partial x_i} e_i$ con $e_i$ $i$-ésimo vector canónico y el índice $i$ es la entrada del vector $\nabla f_o(x)$ de máxima magnitud: $i$ tal que $\left |(\nabla f_o(x))_i \right | = ||\nabla f_o(x)||_\infty$.
Ver 4.2.Algoritmos_para_optimizacion_sin_restricciones. Entonces:
Algoritmo de descenso por coordenadas
Dado un punto inicial $x$ en $\text{dom}f_o$
Repetir el siguiente bloque para $k=0,1,2,...$
- Obtener el índice $i$ de $\nabla f_o(x)$ con máximo valor absoluto: $\left |(\nabla f_o(x))_i \right | = ||\nabla f_o(x)||_\infty$.
- Calcular la dirección de descenso más pronunciado bajo la norma $1$: $\Delta x_{\text{sd}} = - \frac{\partial f_o(x)}{\partial x_i} e_i$ con $e_i$ $i$-ésimo vector canónico.
- 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]:
install.packages("latex2exp",lib="/usr/local/lib/R/site-library/",
repos="https://cran.itam.mx/")
In [1]:
options(repr.plot.width=6, repr.plot.height=6) #esta línea sólo se ejecuta para jupyterlab con R
In [2]:
library(ggplot2)
library(latex2exp)
In [3]:
inc_index<-function(vec,index,h){
'
Auxiliary function for gradient and Hessian computation.
Args:
vec (double): vector
index (int): index.
h (float): quantity that vec[index] will be increased.
Returns:
vec (double): vector with vec[index] increased by h.
'
vec[index]<-vec[index]+h
vec
}
In [4]:
gradient_approximation<-function(f,x,h=1e-8){
'
Numerical approximation of gradient for function f using forward differences.
Args:
f (expression): definition of function f.
x (double): vector 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<-length(x)
gf<-vector("double",n)
for(i in 1:n){
gf[i]=(f(inc_index(x,i,h))-f(x))
}
gf/h
}
In [5]:
Hessian_approximation<-function(f,x,h=1e-6){
'
Numerical approximation of Hessian for function f using forward differences.
Args:
f (expression): definition of function f.
x (double): vector that holds values where Hessian will be computed.
h (float): step size for forward differences, tipically h=1e-6
Returns:
Hf (double): matrix of numerical approximation to Hessian of f.
'
n<-length(x)
Hf<-matrix(rep(0,n^2),nrow=n,ncol=n)
f_x<-f(x)
for(i in 1:n){
x_inc_in_i<-inc_index(x,i,h)
f_x_inc_in_i<-f(x_inc_in_i)
for(j in i:n){
dif<-f(inc_index(x_inc_in_i,j,h))-f_x_inc_in_i-f(inc_index(x,j,h))+f_x
Hf[i,j]<-dif
if(j!=i)
Hf[j,i]<-dif
}
}
Hf/h^2
}
In [6]:
line_search_by_backtracking<-function(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 (expression): definition of function f.
dir_desc (double): vector of descent direction.
x (double): vector 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
}
}
t
}
In [7]:
Euclidian_norm<-function(vec){
sqrt(sum(vec*vec))
}
In [8]:
compute_error<-function(x_obj,x_approx){
'
Relative or absolute error between x_obj and x_approx.
'
if (Euclidian_norm(x_obj) > .Machine$double.eps){
Err<-Euclidian_norm(x_obj-x_approx)/Euclidian_norm(x_obj)
}else
Err<-Euclidian_norm(x_obj-x_approx)
Err
}
Coordinate descent
In [9]:
coordinate_descent<-function(f, x_0, tol,
tol_backtracking,
x_ast, p_ast, maxiter){
'
Method of coordinate descent to numerically approximate solution of min f.
Args:
f (expression): definition of function f.
x_0 (double): vector of initial point for coordinate 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 (double): vector solution of min f, now its required that user knows the solution...
p_ast (double): vector value of f(x_ast), now its required that user knows the solution...
maxiter (int): maximum number of iterations
Returns:
x (double): vector approximation of x_ast.
iteration (int): number of iterations.
Err_plot (double): vector array of absolute error between p_ast and f(x) with x approximation.
of x_ast. Useful for plotting.
x_plot (double): vector array that containts in columns vector of approximations. Last column
contains x, approximation of solution. Useful for plotting.
'
iteration <- 1
x <- x_0
feval <- f(x)
gfeval <- gradient_approximation(f,x)
normgf <- Euclidian_norm(gfeval)
Err_plot_aux <- vector("double",maxiter)
Err_plot_aux[iteration] <- compute_error(p_ast,feval)
Err <- compute_error(x_ast,x)
n <- length(x)
x_plot <- matrix(0,nrow=n,ncol=maxiter)
x_plot[,iteration] <- x
cat(sprintf("I\tNormagf\t\tError x_ast\tError p_ast\tline search\n"))
cat(sprintf("%d\t%.2e\t%0.2e\t%0.2e\t%s\n",iteration,normgf,Err,Err_plot_aux[iteration],"---"))
iteration<-iteration + 1
while(normgf>tol && iteration <= maxiter){
ind_maximo <- which.max(abs(gfeval))
e_canonical <- vector("integer",n)
e_canonical[ind_maximo] <- gfeval[ind_maximo]
dir_desc <- -e_canonical
der_direct <- sum(gfeval*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 <- Euclidian_norm(gfeval)
Err_plot_aux[iteration] <- compute_error(p_ast,feval)
x_plot[,iteration] <- x
Err <- compute_error(x_ast,x)
cat(sprintf("%d\t%.2e\t%0.2e\t%0.2e\t%0.2e\n",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 <- iteration + 1
} #while
cat(sprintf("Error of x with respect to x_ast: %.2e\n",Err))
cat(sprintf("Approximate solution:"))
print(x)
cond <- Err_plot_aux > .Machine$double.eps*10**(-2)
Err_plot <- Err_plot_aux[cond]
if (iteration == maxiter && t < tol_backtracking){
print("Backtracking value less than tol_backtracking, check approximation")
iteration<-iter_salida
x_plot <- x_plot[,1:iteration]
iteration <- iteration -1
}
else{
if (iteration == maxiter) print("Reached maximum number of iterations, check approximation")
iteration <- iteration - 1
x_plot <- x_plot[,1:iteration]
}
list(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 por coordenadas con $x_0=(5,5,1,0)^T$:
In [10]:
fo <-function(x)(x[1]-2)**2 + (2-x[2])**2 + x[3]**2 + x[4]**4
In [11]:
x_ast <- c(2,2,0,0)
x_0 <- c(5,5,1,0)
tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 5
p_ast <- fo(x_ast)
In [12]:
l<-coordinate_descent(fo, x_0, tol, tol_backtracking, x_ast, p_ast, maxiter)
In [13]:
x <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
x_plot <- l[[4]]
In [14]:
gg <- ggplot()
In [15]:
gg +
geom_line(aes(x=1:total_of_iterations,y=Err_plot)) +
xlab('Iterations') + ylab(TeX('Error relativo entre f_o(x^k) y p^*'))
2) $$\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 por coordenadas con $x_0=(0.5,0.5)^T$.
Descenso por coordenadas para caso $C=10$
In [16]:
cte<-10
In [17]:
fo <- function(x) 1/2*(x[1]**2+cte*x[2]**2)
In [18]:
x_ast <- c(0,0)
x_0 <- c(0.5,0.5)
tol <- 1e-7
tol_backtracking <- 1e-14
maxiter <- 50
p_ast <- fo(x_ast)
l<-coordinate_descent(fo, x_0, tol, tol_backtracking, x_ast, p_ast, maxiter)
x <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
x_plot <- l[[4]]
In [19]:
total_of_iterations
In [20]:
gg <- ggplot()
In [21]:
gg +
geom_line(aes(x=1:total_of_iterations,y=Err_plot)) +
xlab('Iterations') + ylab(TeX('Error relativo entre f_o(x^k) y p^*'))
In [22]:
gg +
geom_point(aes(x=x_plot[1,],y=x_plot[2,]),size=2) +
annotate(geom='text', x=0.5, y=0.47,
label=TeX("x^{(0)}", output='character'), parse=TRUE) +
xlab('x') + ylab('y') +
ggtitle(TeX('Iter del método de descenso por coordenadas para $f_o$ con $C=10$'))
In [23]:
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 [24]:
fo <- function(x) exp(x[1]+3*x[2]-0.1) + exp(x[1]-3*x[2]-0.1) + exp(-x[1]-0.1)
In [25]:
x_ast <- c(-3.4654e-01,-7.6725e-06)
In [26]:
x_ast
In [27]:
x_0 <- c(0,0)
tol <- 1e-8
tol_backtracking <- 1e-14
maxiter <- 50
p_ast <- fo(x_ast)
l<-coordinate_descent(fo, x_0, tol, tol_backtracking, x_ast, p_ast, maxiter)
x <- l[[1]]
total_of_iterations <- l[[2]]
Err_plot <- l[[3]]
x_plot <- l[[4]]
In [28]:
total_of_iterations
In [29]:
gg <- ggplot()
In [30]:
gg +
geom_line(aes(x=1:total_of_iterations,y=Err_plot)) +
xlab('Iterations') + ylab(TeX('Error relativo entre f_o(x^k) y p^*'))
In [31]:
gg +
geom_point(aes(x=x_plot[1,],y=x_plot[2,]),size=2) +
annotate(geom='text', x=0.5, y=0.47,
label=TeX("x^{(0)}", output='character'), parse=TRUE) +
xlab('x') + ylab('y') +
ggtitle(TeX('Iter del método de descenso por coordenadas para $f_o$'))
In [32]:
x_plot
Comentario: el método de descenso por coordenadas es lento en situaciones como la situación en el lado izquierdo y rápido como el lado derecho del siguiente dibujo:
Imagen tomada del libro de J. Nocedal & S. J. Wright, Numerical Optimization, Springer.
Referencias: