El método de steepest descent, es un método de optimización, que utiliza el negativo gradiente de x \(-\nabla f(x)\) de la función objetivo \(f(x)\) con el propósito de encontrar el máximo cambio en dirección a un mínimo global o local.
Esto se logra al valuar la función en un punto aleatorio, el cual sera el punto de inicio del algoritmo. A partir de ese punto, utilizando el gradiente y un valor \(\alpha\) es posible acercarse paso a paso al valor mínimo mas cercano utilizando la siguiente función:
\(x^{(k+1)} = x^{(k)} - \alpha \nabla f(x^{(k)})\)
a modo de lograr encontrar:
\(\min f(x^{(k)} - \alpha f(x^{(k)}))\)
En la siguiente imagen se puede ver la superficie de la función objetivo \(f(x) = e^{(0.10(x_2-x_1^{2})^{2} + 0.5(1-x_1)^{2})}\) que se quiere optimizar en el rango de \(-1.5 \leqslant x \leqslant 1.5\) y \(-1.5 \leqslant y \leqslant 1.5\)
Es posible interactuar con las graficas utilizando el mouse
#Libraries that will be used in the project
library(Deriv)
library(plotly)
library(dplyr)
library(matlib)
library(neldermead)
#Function to evaluate the Steepest descent function
f1 <- function (x, y) {
exp((0.1*(y-x^2)^2+0.05*(1-x)^2))
}
#Plot values
minValue <- -1.5
maxValue <- 1.5
x <- y <- seq(minValue, maxValue, length= 200)
#We use outer to create a matrix so the surface can be plotted
graphResults <- outer(x, y, f1)
surfacePlot <- plot_ly(x = seq(minValue,maxValue, length = 200),
y = seq(minValue,maxValue, length = 200),
z = graphResults,
type = "surface",
color=seq(0, 100),
colorbar=list(title='Colorbar'),
colorscale='Viridis',
reversescale = T
)
surfacePlotEs posible interactuar con las graficas utilizando el mouse
Al implementar el método de steepest descent es posible ver: en negro el punto de inicio, en rojo los diferentes “pasos” que tomo el algoritmo y en azul el mínimo local que fue encontrado.
Al mover el cursor sobre los diferentes puntos se puede apreciar: los diferentes puntos calculados por la función de steepest descent y el resultado de la función valuada en cada uno de estos puntos. También es posible observar la progresión del algoritmo viendo las lineas que conectan los diferentes puntos encontrados:
#Initial values for the Steepest descent
iteration <- 1
alfa <- 0.05
step <- 0.03
startingPoint <- x0y0 <- c(-0.3,1)
#Partial derivatives
dfdx <- Deriv(f1, "x")
dfdy <- Deriv(f1, "y")
minimumCondition <- TRUE
#Starting point
points <- c(x0y0[1],x0y0[2], f1(x0y0[1],x0y0[2]))
#While will execute until we have found the minimum local
while(minimumCondition){
#We calculate the building blocks of the steepest descent function
#Alfa and step
alfa <- alfa + step
#Initial value
z0 <- f1(x0y0[1],x0y0[2])
#Derivatives
dfdxValued <- dfdx(x0y0[1],x0y0[2])
dfdyValued <- dfdy(x0y0[1],x0y0[2])
#Step forward toward the minimum value of the function
dfdxValued <- x0y0[1] - alfa*dfdxValued
dfdyValued <- x0y0[2] - alfa*dfdyValued
minimumCondition <- (z0 > f1(dfdxValued,dfdyValued))
if(minimumCondition){
x0y0[1] <- dfdxValued
x0y0[2] <- dfdyValued
}
#After the calculations have been completed we add the results to a matrix so
#this new points can be added to the surface in order to plot the trayectori
#toward the minimum point in the function
points <- c(points, x0y0[1],x0y0[2],f1(x0y0[1],x0y0[2]))
iteration <- iteration + 1
}
#Results
mPoints <- matrix(points, ncol = 3, byrow = TRUE)
mPoints <- as.data.frame(mPoints)
colnames(mPoints) <- c("x", "y", "z")
#Plot of the surface with markers and the trayectory.
#The minimum is plotted in BLUE
surfacePlot <- plot_ly(x = seq(minValue,maxValue, length = 200),
y = seq(minValue,maxValue, length = 200),
z = graphResults,
type = "surface",
color=seq(1, dim(mPoints)[1]),
colorbar=list(title='Colorbar'),
colorscale='Viridis',
reversescale = T,
text = "Not the minimum"
) %>% add_trace(data = mPoints,
x = mPoints[[1]],
y = mPoints[[2]],
z = mPoints[[3]],
mode = "markers+lines",
type = "scatter3d",
marker = list(
size = 5,
color = "red",
symbol = 104
)
) %>% add_trace(data = mPoints,
x = mPoints[[dim(mPoints)[1],1]],
y = mPoints[[dim(mPoints)[1],2]],
z = mPoints[[dim(mPoints)[1],3]],
mode = "markers",
type = "scatter3d",
marker = list(
size = 8,
color = "blue",
symbol = 104
),
text = "Minimum") %>% add_trace(data = mPoints,
x = mPoints[[1,1]],
y = mPoints[[1,2]],
z = mPoints[[1,3]],
mode = "markers",
type = "scatter3d",
marker = list(
size = 8,
color = "black",
symbol = 104
),
text = "Starting point")
surfacePlotEs posible interactuar con las graficas utilizando el mouse
Los resultados de este algoritmo en base al punto de inicio son:
## [1] "Punto de startingPoint -0.3, 1"
## [2] "Mínimo de la función 1.00000211082261"
## [3] "Coordenadas: 0.994190440811671, 0.986357276317702"
## [4] "Alfa: 3.52999999999999"
## [5] "Iteraciones: 117"Para comprobar que la implementación del algoritmo fue adecuado, se decidió tomar otros puntos de inicio como se puede apreciar en el summary de las siguientes gráficas:
Es posible interactuar con las graficas utilizando el mouse
## [1] "Punto de startingPoint 0.6, -0.8"
## [2] "Mínimo de la función 1.00000226730365"
## [3] "Coordenadas: 0.993276602451698, 0.986864817578448"
## [4] "Alfa: 3.55999999999999"
## [5] "Iteraciones: 118"Es posible interactuar con las graficas utilizando el mouse
## [1] "Punto de startingPoint -0.6, 0.8"
## [2] "Mínimo de la función 1.00000183441067"
## [3] "Coordenadas: 0.993966096393236, 0.988342891433304"
## [4] "Alfa: 3.61999999999999"
## [5] "Iteraciones: 120"Es posible interactuar con las graficas utilizando el mouse
## [1] "Punto de startingPoint -1, -1"
## [2] "Mínimo de la función 1.00000250620662"
## [3] "Coordenadas: 0.992937814976849, 0.986278783601224"
## [4] "Alfa: 3.55999999999999"
## [5] "Iteraciones: 118"El objetivo de este método es encontrar el mínimo local, o global, de una función utilizando una función de \(f(x)\) de \(x \in \mathbf{R^{n}}\) utilizando la expansión de Taylor para aproximar la función objetivo. La expansión de la serie de Taylor es una suma infinita de términos los cuales son derivadas de la misma función. Esta suma infinita de \(n\) términos es un polinomio de grado \(n\).
En el método de Newton se expande la función a 2 términos:
\(f(x + a) \approx f(x) + \nabla f(x) h + \frac{1}{2} h^{'}\nabla^{2}f(x) h\)
al igualar la función a 0 y resolver la ecuación linear por h es posible encontrar el siguiente punto en dirección al mínimo local o global.
Es posible observar que el método de Newton requiere menos pasos para encontrar dicho mínimo si comparamos con el método anterior. Debido a que este algoritmo necesita del calculo de la matriz Hessiana en cada iteración del algoritmo, este algoritmo puede resultar muy costoso, computacionalmente, al momento de tratar de utilizarlo para encontrar mínimos locales o globales.
A continuación se pueden observar los resultados del algoritmo valuado en los puntos anteriores.
#Initial values for the Newton’s method including the function at the initial point
iteration <- 1
startingPoint <- x0y0 <- c(-0.3,1)
minimumCondition <- TRUE
z0 <- f1(x0y0[1],x0y0[2])
points <- c(x0y0[1], x0y0[2], f1(x0y0[1],x0y0[2]))
#Calculations
##Partial derivatives (Gradient)
dfdx <- Deriv(f1, "x")
dfdy <- Deriv(f1, "y")
gradientVector <- c(dfdx,dfdy)
##Partial derivatives (Hessian)
dfdxdx <- Deriv(dfdx,"x")
dfdydy <- Deriv(dfdy,"y")
dfdxdy <- Deriv(dfdx,"y")
dfdydx <- Deriv(dfdy,"x")
hessianMatrix <- c(dfdxdx, dfdxdy, dfdydx, dfdydy)
while(minimumCondition){
z0 <- f1(x0y0[1],x0y0[2])
##Gradient at the initial point
gradientVectorVal <- as.vector(c(gradientVector[[1]](x0y0[1],x0y0[2]),
gradientVector[[2]](x0y0[1],x0y0[2])))
##Hessian matrix valued at the initial point
hessianMatrixVal <- matrix(c(hessianMatrix[[1]](x0y0[1],x0y0[2]), hessianMatrix[[2]](x0y0[1],x0y0[2]),
hessianMatrix[[3]](x0y0[1],x0y0[2]), hessianMatrix[[4]](x0y0[1],x0y0[2])),
ncol = 2,
byrow = TRUE)
invHessianMatrixVal <- inv(hessianMatrixVal)
#Calculation of the next point toward the minimum
x0y0 <- x0y0 - invHessianMatrixVal%*%gradientVectorVal
minimumCondition <- (f1(x0y0[1],x0y0[2]) < z0)
#Storing results
points <- c(points, x0y0[1], x0y0[2], f1(x0y0[1],x0y0[2]))
iteration <- iteration + 1
}
#Results
mPoints <- matrix(points, ncol = 3, byrow = TRUE)
mPoints <- as.data.frame(mPoints)
names(mPoints) <- c("x","y","z")
surfacePlot <- plot_ly(x = seq(minValue,maxValue, length = 200),
y = seq(minValue,maxValue, length = 200),
z = graphResults,
type = "surface",
color=seq(1, dim(mPoints)[1]),
colorbar=list(title='Colorbar'),
colorscale='Viridis',
reversescale = T,
text = "Not the minimum"
) %>% add_trace(data = mPoints,
x = mPoints[[1]],
y = mPoints[[2]],
z = mPoints[[3]],
mode = "markers+lines",
type = "scatter3d",
marker = list(
size = 5,
color = "red",
symbol = 104
)
) %>% add_trace(data = mPoints,
x = mPoints[[dim(mPoints)[1],1]],
y = mPoints[[dim(mPoints)[1],2]],
z = mPoints[[dim(mPoints)[1],3]],
mode = "markers",
type = "scatter3d",
marker = list(
size = 8,
color = "blue",
symbol = 104
),
text = "Minimum") %>% add_trace(data = mPoints,
x = mPoints[[1,1]],
y = mPoints[[1,2]],
z = mPoints[[1,3]],
mode = "markers",
type = "scatter3d",
marker = list(
size = 8,
color = "black",
symbol = 104
),
text = "Starting point")
surfacePlotEs posible interactuar con las graficas utilizando el mouse
## [1] "Punto de startingPoint -0.3, 1" "Mínimo de la función 1"
## [3] "Coordenadas: 1, 1" "Iteraciones: 11"Es posible interactuar con las graficas utilizando el mouse
## [1] "Punto de startingPoint -0.3, 0.8"
## [2] "Mínimo de la función 1.14443378067596"
## [3] "Coordenadas: -1.00023907206444, 0.516457040121432"
## [4] "Iteraciones: 2"## [1] "Punto de startingPoint -1, -1"
## [2] "Mínimo de la función 1"
## [3] "Coordenadas: 0.999999999999999, 0.999999999999997"
## [4] "Iteraciones: 11"## [1] "Punto de startingPoint -0.3, 1" "Mínimo de la función 1"
## [3] "Coordenadas: 1, 1" "Iteraciones: 11"Los métodos de búsqueda directa (direct search) son algoritmos que no hacen necesitan hacer el calculo de derivadas para calcular los mínimos de una función. Estos métodos valúan la función en una cantidad finita de puntos, según los resultados obtenidos en cada uno de estos puntos, el algoritmo toma la decisión de cual sera su siguiente paso.
Una de las razones por la que se desea implementar métodos que no necesiten el uso de derivadas para encontrar el mínimo local es que las funciones que se desean optimizar deben cumplir con las condiciones First Order Necessary Conditions (FONC) para problemas sin restricciones. \(x_*\) es un mínimo local de \(f\) y \(f\) es continuamente diferenciable si: \(\nabla f(x_*) = 0\)
Esto puede llegar a ser un problema si la primer derivada no es posible de encontrar o si el costo computacional de calcular las derivadas es demasiado alto.
Uno de los intercambios mas grandes que hay al hacer uso de estos métodos es que: al tratar de lograr de reducir el costo computacional la eficiencia de los algoritmos puede resultar no ser la optima, ademas de la informacion que no se logra obtener acerca de la dirección del máximo cambio con el Gradiente y la curvatura de la curva en el punto valuado con la matriz Hessiana. Los problemas que presentan los algoritmos de búsqueda directa también se pueden presentar al momento de detener el algoritmo debido a que han encontrado un punto optimo. Debido a la falta de informacion respecto al punto que se esta valuando, es difícil conocer si realmente hay alguna mejor solución o no.
Algunos de las características mas relevantes de los algoritmos que no utilizan derivadas para la búsqueda de puntos óptimos son:
Al igual que los algoritmos que utilizan derivadas, estos métodos deben implementar métodos para alejarse de la estacionalidad y lograr salir de ella. Para eso implementan métodos como Positive Bases o Spanning sets.
Deben tener control sobre la geometría de la muestra que se utiliza, asegurándose que cualquier indicio que un punto sea estacionario se considere como tal y así se pueda salir de ese punto.
El algoritmo debe detenerse en el momento que el step que se esta tomando, en búsqueda del mínimo optimo, sea cero o este muy cercano a cero.
Es una función de la base de R y contiene varios algoritmos que no necesitan del calculo de derivadas para la búsqueda de los puntos óptimos.
El primer algoritmo que se utiliza es el que la función tiene por Default el cual es: Nelder and Mead.
Es un método heuristico de búsqueda capaz de converger en puntos no estacionarios. Utiliza el método simplex como base utilizando segmentos de rectas, triángulos en planos, tetraedros en 3D y así consecutivamente.
Este método utiliza la desviación estándar de la función del simplex que se esta resolviendo para determinar el step donde debe de detenerse y así detener el ciclo. En la grafica es posible observar que la forma en la que el algoritmo encuentra el minimo no sigue un patron regular sino que se ve, hasta cierto grado, un patron mas errático.
#Starting points
i <- 1
vals <- c()
startingPoint <- c(-0.3,1)
#Same function as before but using a vector instead of using two coordenates
f2 <- function(x){
z <- exp((0.1*(x[2]-x[1]^2)^2+0.05*(1-x[1])^2))
vals <<- c(vals, x[1], x[2], z)
i <<- i + 1
exp((0.1*(x[2]-x[1]^2)^2+0.05*(1-x[1])^2))
}
#Optim function, from the base of R
directSearchSol <- optim(fn = f2, startingPoint)
#Graph of the different results obtained from the function
mPoints <- matrix(vals, byrow = TRUE, ncol = 3)
mPoints <- as.data.frame(mPoints)
names(mPoints) <- c("x","y","z")
surfacePlot <- plot_ly(x = seq(minValue,maxValue, length = 200),
y = seq(minValue,maxValue, length = 200),
z = graphResults,
type = "surface",
color=seq(1, dim(mPoints)[1]),
colorbar=list(title='Color bar'),
colorscale='Viridis',
reversescale = T,
text = "Not the minimum"
) %>% add_trace(data = mPoints,
x = mPoints[[1]],
y = mPoints[[2]],
z = mPoints[[3]],
mode = "markers+lines",
type = "scatter3d",
marker = list(
size = 5,
color = "red",
symbol = 104
)
) %>% add_trace(data = mPoints,
x = mPoints[[dim(mPoints)[1],1]],
y = mPoints[[dim(mPoints)[1],2]],
z = mPoints[[dim(mPoints)[1],3]],
mode = "markers",
type = "scatter3d",
marker = list(
size = 8,
color = "blue",
symbol = 104
),
text = "Minimum") %>% add_trace(data = mPoints,
x = mPoints[[1,1]],
y = mPoints[[1,2]],
z = mPoints[[1,3]],
mode = "markers",
type = "scatter3d",
marker = list(
size = 8,
color = "black",
symbol = 104
),
text = "Starting point")
surfacePlot## [1] "Punto de inicio: ( -0.3, 1 )"
## [2] "Mínimo de la función 1.00000000785118"
## [3] "Coordenadas: 1.0002971300155, 1.00040896053752"
## [4] "Iteraciones: 69"## [1] "Punto de inicio: ( 1, -0.3 )"
## [2] "Mínimo de la función 1.00000000261786"
## [3] "Coordenadas: 1.00022881581729, 1.00045816455671"
## [4] "Iteraciones: 65"## [1] "Punto de inicio: ( -1, -1 )"
## [2] "Mínimo de la función 1.00000000580395"
## [3] "Coordenadas: 1.00028355063041, 1.00070074446339"
## [4] "Iteraciones: 83"