Optimización

Sobre Optimización

Para entender en general las ideas de esta entrada basta con tener algunos conocimientos de Calculo Diferencial, de una y varias variables. Si no se cuenta con ideas o conocimientos de ello, no implica el no entender los algoritmos de optimización que se comentan o por lo menos visualizar como funcionan.

Parto de ideas sencillas hasta dar unos ejemplos de optimización estocástica. Espero poder trasmitir la idea de los métodos y que puedan servirles de ayuda para investigar más respecto al tema.

Omito hablar sobre el proceso ” gradiente descendiente estocástico” pero dejo referencias donde pueden consultar en qué consiste este método y por que es que se hace amplio uso en técnicas de Machine Learning.

 La idea de optimizar  seguro es familiar cuando vemos como nuestras madres resuelven la siguiente pregunta: ¿cómo hacer para comer carne con solo cierta cantidad de dinero?

Aunque parece broma la pregunta anterior es una optimización, tenemos que encontrar las cantidades de bienes que podemos comprar considerando que solo tenemos una cantidad de dinero y debemos de incluir en la lista de bienes la carne. Obviamente esto no dice que nuestras madres hagan uso de las técnicas que comento esta entrada, pero espero que de una idea burda de lo que es optimizar.

Otro ejemplo es cómo maximizar las ganancias de una empresa o negocio, considerando R cantidad de ingresos y C cantidad de costos. Claramente este problema implica conocer qué genera nuestros ingresos y costos, para plantear un modelo idealizado sobre la empresa.

En general la idea de optimizar es encontrar los máximos o mínimos de una función. Para tener una idea de qué significa hago dos gráficas como ejemplo en R project.

#Creamos una función
#Ejemplo 1-Mínimo de una función

Func_1<-function(x){
 x[1]^2 
      }

resultado<-optimize(Func_1,c(-5,5),maximum=FALSE,tol=1e-8) 

#Vemos los resultados
print(resultado$minimum)
print(resultado$objective)

#Generamos el rango de valores para x y y=f(x)
x<-seq(-5,5,length.out=100)
y<-Func_1(expand.grid(x))

#Graficamos identificando el mínimo

plot(x,y,xlab="x",ylab="f(x)",col="4")
points(resultado$minimum,resultado$objective, col="red",pch=19)
rect(resultado$minimum-0.3, resultado$objective-0.7, resultado$minimum+0.3,resultado$objective+0.7)

Ejemplo_opt_mínimo

La gráfica muestra en color rojo el mínimo de la función, lo obtenemos por medio de la función optimize. La cual necesita una función objetivo, un rango de valores para optimizar y un nivel de tolerancia en el cálculo. Esto último es el número que utilizará para realizar su aproximación. Recomiendo revisar la información de la función por medio del comando help.

Cuando pedimos el resultado por medio de print(resultado$minimun)en el código, se observa que el valor es 0, lo cual es evidente al ver la gráfica. Hago otro ejemplo con una función diferente.

#Ejemplo 2
#Creamos una función
Func_2<-function(x){
 x^2+6*sin(x*pi*2)
 }
 
#Optimizamos
resultado2_1<-optimize(Func_2,c(-5,5),lower="-1",upper="3",maximum=FALSE,tol=.Machine$double.eps^0.25) 
resultado2<-optimize(Func_2,c(-5,5),maximum=FALSE,tol=1e-8)
#Vemos los resultados
print(resultado2$minimum)
print(resultado2$objetive)

#Rango de calores para nuestras gráficas 
x2<-seq(-5,5,length.out=100)
y2<-Func_2(expand.grid(x2)) 

#Graficamos
plot(x2,as.matrix(y2)[,1],xlab="x",ylab="f(x)",col="4")
points(resultado2$minimum,resultado2$objective, col="red",pch=19)
rect(resultado2$minimum-0.3, resultado2$objective-0.7, resultado2$minimum+0.3,resultado2$objective+0.7)

points(resultado2_1$minimum,resultado2_1$objective,col="green",pch=19)
rect(resultado2_1$minimum-0.3,resultado2_1$objective-0.7,resultado2_1$minimum+0.3,resultado2_1$objective+0.7)

Ejemplo_opt_mínimo_2

La diferencia de esta última gráfica es que se tienen dos puntos, rojo y verde, que a primera vista parece que el verde es el “mínimo” de toda la gráfica. Pero lo que realmente hice fue acotar el cálculo del optimización, por lo cual si observamos en la función  optimize hice distintas cotas (lower, upper) para que se ejecutara la función. Entonces lo que se ve son mínimos locales y dependen de la restricción que se haga a los cálculos.

¿Qué es lo que hace la función optimize? Da el punto máximo o mínimo según las restricción para el cálculo o en otro caso regresa el máximo o mínimo global.

La diferencia en local y global, es que el primero se restringe a una pequeña región, ejemplo a las x que están entre -2 y -1, el punto rojo de la gráfica es un mínimo local. Por otro lado, global se refiere a que para todo el rango de las x’s que se revisan no se puede encontrar otro valor de f(x) menor. En el ejemplo el punto verde de la segunda gráfica es un mínimo global y el punto rojo de la primera también.

Los dos ejemplos muestran gráficamente cual es la idea de obtener un mínimo que depende de una sola variable x. Pero en la mayoría de problemas de optimización no son funciones de una sola variable, sino de varias variables. Ejemplo, cuando  se busca maximizar ganancias o minimizar costos, es difícil que dependa de un solo producto que se vende o el costo de un solo insumo.

Entonces, la idea de optimizar se puede generalizar para varias variables. Pero en este caso se tienen varios métodos para hacer el cálculo y cada uno tienen ventajas y desventajas, depende del tipo de problema o función que  se desea optimizar.

¿Qué tiene que ver todo esto con Machine Learning?

 En general todos los algoritmos de Machine Learning tienen tres piezas, datos de entrada, datos de salida y una función entre los datos. Pero por medio de un conjunto de datos de entrenamiento y un algoritmo de aprendizaje se trata de obtener la mejor función entre los datos.

Entonces, la optimización está ligada a los algoritmos, ya que siempre en los algoritmos  se tiene una función objetivo y un conjunto de parámetros que se busca estimar. Por lo cual es fundamental la optimización al aplicar los algoritmos de Machine Learning.

Por ejemplo, en la regresión lineal. Sobre la cual en otra entrada se mencionó que la función lm permite estimarla. Se toman datos de altura y peso de personas, los datos se pueden descargar desde de Github.

Al usar la función lm  no se muestra lo que hace o como estima los valores de la regresión  a y b, donde la ecuación es:

Peso=a+b*Altura

Lo que se puede hacer es estimar los valores de los parámetros por medio de una optimización y comparar con los valores que se estiman en R project por medio de la función lm.

#Ejemplo de estimación de Parámetros
#Definimos la función de la regresión

Altura.a.peso <- function(altura, a, b)
{
 return(a + b * altura)
}

#Función de error cuadrado
error.cuadrado <- function(altura.peso, a, b)
{
 predic <- with(altura.peso, Altura.a.peso(Altura, a, b))
 errors <- with(altura.peso, Peso - predic)
 return(sum(errors ^ 2))
}

#Optimización

optim(c(0, 0),function (x){error.cuadrado(altura.peso, x[1], x[2])})

El código anterior hace dos cosas. La primera es definir la función de la cual se busca estimar sus parámetros a y b, pero esa función solo define rectas y para elegir la que mejor describe la relación entre  los datos se necesita otra función, error.cuadrado.

La función error.cuadrado como se observa depende de la función de Altura.a.peso. Así que se tiene solo una función compuesta con otra y  se busca en qué valores se minimiza esa función.

La función error.cuadrado permite encontrar la mejor recta entre los valores de Peso contra Altura. Además indica que se encontró la recta que tienen el menor error cuadrado medio, que es medido entre el Peso menos el valor que estimamos mediante la función Altura.a.peso.

Comparo los valores que se encuentran con la estimación por optimización y por la función lm.

#Comparación entre las estimaciones

 >lm(Peso~Altura,data=altura.peso)

Call:
lm(formula = Peso ~ Altura, data = altura.peso)

Coefficients:
(Intercept) Altura 
 -350.737 7.717 


> optim(c(0, 0),function (x){error.cuadrado(altura.peso, x[1], x[2])})
$par
[1] -350.786736 7.718158

$value
[1] 1492936

$counts
function gradient 
 111 NA 

$convergence
[1] 0

$message
NULL

Se nota que los valores de la regresión por medio de lm y los valores estimados por la función optim son prácticamente los mismos. Entonces se obtiene mediante optimización de una función objetivo error.cuadrado  el valor de sus parámetros a, b.

Se observa que la optimización fue hecha para una función de dos variables, el método que se usa para optimizar es Nelder-Mead, más adelante hablo de él con algo de detalle.

En conclusión, lo que parecía que operaba como una caja-negra cuando se usa la función lm, ahora se puede entender o tener una visión de cómo posiblemente es el proceso de estimación de los parámetros.

Aprovechando el ejemplo de la regresión lineal se puede modificar la función de error para estimar lo que se llama Rigde Regression. Es otro tipo de regresión en la cual para estimar los parámetros se afecta su error cuadrado o se penalizan los errores. El objetivo de esta técnica es evitar la sobre estimación de los datos y disminuir la varianza del ajuste lineal.

La Ridge Regression implica agregar un nuevo parámetro, lambda; el cual para hacer una buena elección se requiere hacer uso de la Regularización y Cross-Validation[1]. Se haciendo la modificación a la estimación para revisar el comportamiento al optimizar.

Haciendo dos estimaciones de Ridge Regression y se compara la estimación con respecto a las que se realizo anteriormente.

#Calculo de Ridge Regression

ridge.error <- function(altura.peso, a, b, lambda)
{
 predic <- with(altura.peso, Altura.a.peso(Altura, a, b))
 errores <- with(altura.peso, Peso - predic)
 return(sum(errores ^ 2) + lambda * (a ^ 2 + b ^ 2))
}

lambda=1
optim(c(0, 0),function (x){ridge.error(altura.peso, x[1], x[2], lambda)})

$par
[1] -340.434108 7.562524
lambda=0.75
optim(c(0, 0),function (x){ridge.error(altura.peso, x[1], x[2], lambda)})

$par
[1] -343.13351 7.60324

El ejemplo muestra el uso de la optimización al estimar los parámetros de un modelo. En la Rigde Regression se definió el valor de lambda, lo cual hace que se tenga ese parámetro “libre”  en el modelo y por ello es necesario para un buen uso de esta técnicas elegir la mejor lambda por medio de regularización y cross-validation.

Entonces al usar el paquete ridge de R project para estimar la ridge regression, ya no parece del todo extraño como hace el cálculo el software. Podría usar la función linearRidge del paquete ridge y comparar por medio de optimización, regularización y cross-validacón la elección del lambda.

Algo que no hice pero que puede resultar ser un buen ejercicio es revisar el artículo de investigación de Erika Cule y estudiar el método de elección de lambda contra otro método para elegirla.

Ejemplo_opt_regresion

La gráfica muestra como cambia ligeramente la pendiente de las rectas. Haciendo solamente una modificación en la estimación de los errores, regresión lineal y ridge regression.

Espero que este ejemplo haga notar que los algoritmos de Machine Learning  que parecen cajas negras, en el fondo están combinando cosas, y una de esas es un proceso de optimización para estimar los parámetros.

Este tema un trabajo serio al momento de procesar volúmenes masivos de datos, quienes tengan interés pueden revisar las investigaciones del grupo de optimización de la universidad Carnegie-Mellon.

Cabe mencionar que en las investigaciones sobre optimización su relevancia puede implicar una disminución en el costo computacional o en la mejora de algoritmos.

Ejemplos de técnicas de Optimización

Lo siguiente que comparto en esta entrada son técnicas de optimización convexa y dos técnicas de optimización estocástica. Para algunas solo pondré ejemplos de cómo se estima y como interpretar el resultado siguiendo unos ejemplos sencillos publicados por Jason Brownlee en su página. En otros pondré ejemplos de como usarlas de manera más detallada.

La explicación de cada uno de los métodos no es exhaustiva, tanto por que está fuera de mi alcance como de la intensión de esta entrada.

Dejo referencias en los comentarios y al final de la entrada algunas ligas a notas, reportes de investigación y libros para consultar más sobre el tema.

Método 1.-Búsqueda de la sección dorada

No repito el código, pero la primera gráfica que se hizo en esta entrada hacía uso de la función optimize y muestra un punto rojo que es el mínimo de la función. El método por medio del cual estima ese punto es conocido como  búsqueda de la sección dorada.

Su uso se limita a función que dependen de una sola variable y se busca el punto óptimo entre dos puntos que se definen del rango de las x’s. El nombre provienen de la relación entre los tres primeros puntos que toma el método, ya que esto tienen una proporción con el número conocido como razón dorada o número áureo.

Entre los limitantes de dicha técnica es que solo funciona para funciones de una sola variable, de preferencia que sea una función convexa y que se pueda identificar que se tienen un único óptimo entre la región que define para realizar el cálculo.

Función Rosenbrock

Para tres de los siguientes ejemplos se hace uso de la función Rosenbrock. Es una función no convexa que se usa para probar los métodos de optimización o rendimiento de los algoritmos. En general se puede usar de dos variables o se puede generalizar para n-variables. Forma parte de un conjunto de funciones que son usadas para hacer pruebas con técnicas de optimización.

La ecuación en general de la función para dos variables es:

f(x,y)=(a-x)^2+b(y-x^2)^2

La visualización de la función para un par de valores a=1 y b=100.

#Función Rosenbrock
#La función es f(x,y)=(1-x)^2+100*(y-x^2)^2

nx<-21
ny<-21

x<-seq(-5,5,length=nx)
y<-seq(-5,5,length=ny)

#Se genera la salida de los valores en z

z<-outer(x,y,function(x,y)(1-x)^2+100*(y-x^2)^2)

#Escala de colores
hgt <- 0.25 * (z[-nx,-ny] + z[-1,-ny] + z[-nx,-1] + z[-1,-1])
hgt <- (hgt - min(hgt))/ (max(hgt) - min(hgt))

#Se construye la gráfica en perspectiva con cierto ángulo para visualizarla
persp(x, y, z, col =cm.colors(10)[floor(9* hgt + 1)],theta = 35,phi=30)

Función_Rosenbrock

Método 2.-Nelder-Mead

El método es adecuado para funciones no lineal de varias variables y se busca el óptimo sin restricciones en el dominio. La herramienta que usa para buscar el óptimo en su algoritmo es una estructura geométrica conocida como simplex. El método para una función de dos variables parte de la elección de 3 puntos iniciales de los cuales bajo un proceso de varias etapas se va modificando el simplex hasta aproximarse al punto óptimo de la función. El proceso de transformación del simplex es la parte principal del algoritmo.

Es adecuado para funciones de varias variables y al no requerir la derivada para su algoritmo puede ser usada para funciones no diferenciables, no continuas o con ruido. Un resumen de sus ventajas y desventajas se pueden consultar en, Nelder-Mead_algorithm.

#Creamos una función

#Definición de una función Rosenbrock, el optimo de ubica en (1,1)

rosenbrock <- function(v){ 
 #Valores de a=1, b=100
(1 - v[1])^2 + 100 * (v[2] - v[1]*v[1])^2
}

# Identificación del mínimo local por el método Nelder-Mead 
resultado3 <- optim(c(runif(1,-3,3), runif(1,-3,3)),rosenbrock,NULL,method="Nelder-Mead",control=c( maxit=100,reltol=1e-8,alpha=1.0, beta=0.5,gamma=2.0))

#Revisión de resultados
print(resultado3$par) 
print(resultado3$value) 
print(resultado3$counts)

#Gráfica de cortorno
x3<- seq(-3, 3, length.out=100)
y3<- seq(-3, 3, length.out=100)
z3<- rosenbrock(expand.grid(x3, y3))
contour(x3, y3, matrix(log10(z3), length(x3)), xlab="Valores de X",ylab="Valores de Y")
#Punto óptimo
points(resultado3$par[1], resultado3$par[2], col="red", pch=19)
#Úbicación del punto óptimo
rect(resultado3$par[1]-0.2, resultado3$par[2]-0.2, resultado3$par[1]+0.2, resultado3$par[2]+0.2, lwd=2)

Ejemplo_opt_Nelder-Mead2

Método 3.-Gradiente Descendente o Algoritmo de Descenso Rápido.

El métodos hace uso del gradiente de una función y cuando busca el mínimo hace uso del negativo del gradiente para aproximarse al óptimo.

Funciona bien si la función es convexa, ya que el método estimará el mínimo global. El método es muy sensible a la elección del punto inicial para hacer la estimación y como depende de un parámetro para ir haciendo los avances, es susceptible a malas elecciones de él.

Como referencias dejo las siguientes ligas: Aplicación del método a procesamiento de imágenes y curso de Andrew Ng en línea.

Usamos el método para calcular el mínimo de una función que tienen su mínimo global en el punto (0,0), el cual se visualiza en la imagen.

Ejemplo_función_min_global

#Creamos una función
#Se define una función, el óptimo se ubica en (0,0)

Func_3<- function(x) {
x[1]^2 + x[2]^2
}

#Definición de la derivadas
derivada <- function(x) {
c(2*x[1], 2*x[2])
}

#Definición del método del gradiente
gradient_descent <- function(func, derv, start, step=0.05, tol=1e-8) {
 pt1 <- start
 grdnt <- derv(pt1)
 pt2 <- c(pt1[1] - step*grdnt[1], pt1[2] - step*grdnt[2])
 
 while (abs(func(pt1)-func(pt2)) > tol) {
 pt1 <- pt2
 grdnt <- derv(pt1)
 pt2 <- c(pt1[1] - step*grdnt[1], pt1[2] - step*grdnt[2])
 print(func(pt2)) #regresa el valor en el punto
 }
 pt2 # Regresa el valor óptimizado
 }

#Método del Gradiente Descendente
resultados4 <- gradient_descent(Func_3, derivada,c(runif(1,-3,3), runif(1,-3,3)), 0.05,1e-8)

#Resumen de resultados
print(resultados4) 
print(Func_3(resultados4)) 

#Gráfica de contorno
x4 <- seq(-3, 3, length.out=100)
y4 <- seq(-3, 3, length.out=100)
z4 <- Func_3(expand.grid(x4, y4))
contour(x4, y4, matrix(z4, length(x4)), xlab="x",ylab="y")
#Dibujamos el punto óptimo
points(resultados4[1], resultados4[2], col="red", pch=19)
#Dibujamos el cuadrado del óptimo
rect(resultados4[1]-0.2, resultados4[2]-0.2, resultados4[1]+0.2, resultados4[2]+0.2, lwd=2)

Ejemplo_gradiente_descendente

Método 4.-Gradiente Conjugado

Es adecuado para funciones no lineales de varias variables. Le idea del método es construir una base ortogonal de vectores y hacer uso de ella para hacer la búsqueda del óptimo de un modo más eficiente. Aspectos teóricos de manera breve puede consultar en las notas en español del Prof. Maximenko y una buena referencia para ver más detalles recomiendo las notas del Prof.Shewchuk.

Igual que otros métodos su eficiencia se ve mermada ante una mala elección del punto inicial en  funciones no  convexas, pero una de sus ventajas es que no hace uso de la matriz Hessiana en su algoritmo lo cual permite hacer uso del método para funciones de muchas variables de manera eficiente.

Para hacer el ejemplo hacemos uso de la función Rosenbrock.

#Definimos la función Rosenbrock

rosenbrock <- function(v) { 
(1 - v[1])^2 + 100 * (v[2] - v[1]*v[1])^2
}

#Definición de la derivada de la función
derivada2 <- function(v) {
c(-400 * v[1] * (v[2] - v[1]*v[1]) - 2 * (1 - v[1]), 
200 * (v[2] - v[1]*v[1]))
}


Ejemplo_opt_gradente_conjugado2

Método 5.-BFGS

Este método hace uso tanto del gradiente como de una aproximación a la inversa de la matriz Hessiana de la función, esto para hacer una aproximación al cálculo de la segunda derivada. Por ser un método de aproximación a la segunda derivada se dice que es un método Quasi-Newtoniano. El problema con el método es el costo computacional para funciones de muchas variables, ya que se almacena una matriz cuadrada de datos tan grande como el cuadrado de la cantidad de variables.

Es adecuado para funciones no lineales de varias variables y la búsqueda del óptimo sin restricciones. Si bien el método hace que la convergencia al óptimo sea rápida por ser de aproximación a la segunda derivada, el costo de procesamiento es bastante alto.

Se pueden consultas las notas de Galen Andrew o de Peter Blomger, para ver detalles del método y el algoritmo.

#Método BFGS

#Definición la función Rosenbrock, se óptimiza en (1,1)
rosenbrock <- function(v) { 
(1 - v[1])^2 + 100 * (v[2] - v[1]*v[1])^2
}
 
#Definición del Gradiente
derivada4 <- function(v) {
c(-400 * v[1] * (v[2] - v[1]*v[1]) - 2 * (1 - v[1]), 200 * (v[2] - v[1]*v[1]))
}
 
#Calculo del método BFGS 
resultados6 <- optim(c(runif(1,-3,3), runif(1,-3,3)),rosenbrock, derivada4,method="BFGS",control=c(maxit=100,reltol=1e-8))
 
#Resumen de resultados
print(resultados6$par)
print(resultados6$value) 
print(resultados6$counts)
 
#Gráfica de contorno
x5 <- seq(-3, 3, length.out=100)
y5 <- seq(-3, 3, length.out=100)
z5 <- rosenbrock(expand.grid(x5, y5))
contour(x5, y5, matrix(log10(z5), length(x5)), xlab="x", ylab="y")

points(resultados6$par[1], resultados6$par[2], col="red", pch=19)
# draw a square around the optima to highlight it
rect(resultados6$par[1]-0.2, resultados6$par[2]-0.2, resultados6$par[1]+0.2, resultados6$par[2]+0.2, lwd=2)

Ejemplo_BFGS

Los métodos anteriores son considerados de optimización convexa, para ver detalles teóricos recomiendo consultar un texto especializado en el tema. Una buena referencia es el libro de Stephen P. Boyd y Lieven Vandenberghe el cual se puede descargar libremente desde su página. El segundo autor tienen dos cursos en línea con formato de presentaciones; optimización convexa y métodos de optimización para sistemas de larga escala.

Optimización Estocástica (o ¿computación evolutiva?)

Los siguientes ejemplos que comparto son abordados con técnicas de optimización estocástica, las cuales también son estudiadas en computación evolutiva. Por supuesto que no doy explicaciones exhaustivas, pero dejo referencias al final de la entrada.

Los métodos estocásticos que generalmente se usan son: algoritmos genéticos (genetic algorithm)  y recorridos simulados (Simulated anneling).

En el ejemplo considero el problema siguiente, si se  tiene un grupo de viajeros y todos parten de lugares distintos pero se dirigen al mismo sitio, si todos están de acuerdo en reunirse en un punto y se esperaran para compartir transporte terrestre, esto implica que entre los posibles vuelos tienen que elegir aquellos que mejor les convenga en costo, tiempo de espera, horarios de llegada y salida. Así que se tienen varios factores a considerar en la optimización.

Como cualquier problema de optimización se requiere tener la función objetivo, en este ejemplo se le llamará función de costo.

Los datos pueden ser descargados desde la página de Toby Segaran, cuenta con con una versión del código en python. Hice una adaptación en R y considerando las restricciones del problema y de los datos.

Los pasos para abordar el problema son los siguientes:

  1. Se representarán las soluciones como cadenas de números.
  2. Se presentaran las soluciones como un tabla.
  3. Se definirán una función objetivo o función de costo.
  4. Se estimaran soluciones por medio del método búsqueda aleatoria, recorrido simulados y algoritmos genéticos.

Paso 1

Se representa las solución como una cadena de 12 dígitos que pueden estar entre 1 y 9. Ejemplo, solución=[1,3,4,5,8,10,2,3,7,9,3,5], representa que la persona uno toma el primero vuelo y regresa en su tercer vuelo posible. La persona dos toma el cuarto vuelo y regresa en su quinto vuelo posible, así sucesivamente para las 6 personas.

#Carga de los datos desde el directorio
viajes.path<-file.path("data","schedule.txt")

viajes<-read.table(viajes.path,sep=",",stringsAsFactor=FALSE)

colnames(viajes)<-c('Origen','Destino','Hora de Salida','Hora de arrivo','Precio')

#Generamos lista de viajeros

gente=list(c('Seymour','Franny','Zooey','Walt','Buddy','Les'),c('BOS','DAL','CAK','MIA','ORD','OMA'))

#Destino al que viajan
Destino=c('LGA')

#Definimos una función para medir el tiempo, será necesaria en otras funciones
#Librería para manipulas fechas
library(lubridate)

#Función que nos regresa los minutos
minutos<-function(horario){
 min<-hour(hm(horario))*60+minute(hm(horario))
 return (min)}

Paso 2

Teniendo los datos cargados, se procede a representar por una pequeña tabla la solución al problema. Para ello se define la función impri.horario.

#Impresión de soluciones

impri.horarios<-function(r){

#r es la cadena con datos de los viajeros y lugars
 for( i in 1:(length(r)/2)){
 nombre=gente[[1]][i]
 origen=gente[[2]][i]
 Destino=c('LGA')
 
 salida=viajes[which(viajes$Origen==origen),][r[i*2-1],]
 regreso=viajes[which(viajes$Origen==Destino),][r[i*2],]
 
 a<-c(nombre,origen,paste(salida[,3],"-",salida[,4],sep=""," ",salida[,5]),paste(regreso[,3],"-",regreso[,4],sep=""," ",regreso[,5]))
 print(a,sep="\t")
 
 } 
 
} 

Paso 3

Se define una función objetivo o la función de costo, es un ejemplo solamente, por lo cual la función puede ser modificada para una situación real de optimización. Pero cabe mencionar que la definición de ella es la parte más importante para aplicar estos algoritmos.

#Creamos una función a optimizar
#Función Objetivo o Función de Costo

costo.tiempo<-function(r){
 precio.total=0
 ultimo.arrivo=0
 salida.temprana=24*60
 
 for(i in 1:(length(r)/2)){
 #nombre=gente[[1]][i]
 origen=gente[[2]][i]
 Destino=c('LGA')
 salida=viajes[which(viajes$Origen==origen),][r[i*2-1],]
 regreso=viajes[which(viajes$Origen==Destino),][r[i*2],]
 precio.total=precio.total+salida[,5]+regreso[,5]
 
 if(ultimo.arrivo<minutos(salida[,4])){ultimo.arrivo=minutos(salida[,4])}
 if(salida.temprana>minutos(regreso[,3])){salida.temprana=minutos(regreso[,3])}
 
 }
 tiempo_espera=0 
 for( i in 1:(length(r)/2)){
 origen=gente[[2]][i]
 Destino=c('LGA')
 salida=viajes[which(viajes$Origen==origen),][r[i*2-1],]
 regreso=viajes[which(viajes$Origen==Destino),][r[i*2],]
 
 tiempo_espera=tiempo_espera+(ultimo.arrivo-minutos(salida[,4]))
 tiempo_espera=tiempo_espera+(minutos(regreso[,3])-salida.temprana)
 }
 if(ultimo.arrivo>salida.temprana){precio.total=precio.total+50} 
 
 return(precio.total+tiempo_espera)
 }

Si se hace un prueba con esta función se nota que devuelve un número entero, el cual representa el valor de la solución o configuración que asignamos a los viajantes.

#Ejemplo del uso de la función de costo
#Definimos una posible solución
 vector
 [1] 6 7 4 2 6 7 4 2 6 2 6 8
> costo.tiempo(vector)
[1] 3806

Paso 4

Para comparar métodos uso el método de búsqueda aleatoria. Consiste básicamente en tomar posibles soluciones y comparar entre las que se van encontrado, así con un proceso de comparación se elige la solución que minimiza  la función objetivo. Como es de esperar no es un método eficiente, pero permite acercarse a las ideas de la optimización estocástica.

#Búsqueda Aleatoria

optimizacion.aleatoria<-function(funcost){
 #dominio va ser igual al número máximo que puede tomar cada entrada
 best=999999999
 bestr=list()
 for( i in (0:1000)){
 r=sample(10,size=12,replace=TRUE) 
 cost=funcost(r)
 
 if(cost<best){
 best=cost
 bestr=r}
 }

 return(bestr)}

Cabe mencionar que el método fue adaptada para el problema, el algoritmo es estándar. Se ponemos a prueba el método y se evalúa la función de costo.

#Resultados del método Búsqueda Aleatoria

> optimizacion.aleatoria(costo.tiempo)
 [1] 8 9 6 7 2 8 8 7 6 8 7 8

> solucionba=c(8, 9, 6, 7, 2, 8, 8, 7, 6, 8, 7, 8)

> costo.tiempo(solucionba)
[1] 3689

> impri.horarios(solucionba)
[1] "Seymour" "BOS" "17:11-18:30 108" "18:25-20:34 205"
[1] "Franny" "DAL" "13:54-18:02 294" "15:07-17:21 129"
[1] "Zooey" "CAK" "8:27-10:45 139" "16:35-18:56 144"
[1] "Walt" "MIA" "17:07-20:04 291" "15:07-17:21 129"
[1] "Buddy" "ORD" "14:22-16:32 126" "16:35-18:56 144"
[1] "Les" "OMA" "15:03-16:42 135" "16:35-18:56 144"

Ahora usamos el algoritmo de recorrido simulados, tienen su origen en un problema de física. Por lo cual algunos de sus parámetros prevalecen en uso por analogía del fenómeno físico.

La idea es iniciar con una solución aleatoria y con un parámetro de temperatura alto, este último representa en la búsqueda la guía de un proceso de “enfriamiento” donde conforme decrece se van comparando los valores de la función objetivo. De los cuales se elige una posible solución al óptimo global que se espera obtener al tener un nivel de temperatura bajo.

El método es una adaptación al método Metropolis-Hastings, por lo cual en el algoritmo aparece una exponencial para comparar un valorar aleatorio y continuar en las comparaciones de los valores de la función objetivo.

#Algoritmo recorrido simulado
annealingoptimize<-function(costf,T=10000.0,cool=0.95,step=1){
 #Creamos una solución aleatória
 vec=sample(9,size=12,replace=TRUE)

 while(T>0.1){
 #Elejimos un índice
 i=sample(8,size=1)
 #Definimos la dirección del algortimos
 x=c(-1,0,1)
 dir=sample(x,1) 
 vecb=vec
 vecb[i]=vecb[i]+dir
 
 if(vecb[i]<0){vecb[i]=0}
 if(vecb[i]>9){vecb[i]=9}
 
 ea=costf(vec)
 eb=costf(vecb)
 
 
 if(eb<=ea){
 vec=vecb
 }
 else
 { 
 a=runif(1,min=0,max=1)
 
 if(a<exp((-eb-ea)/T))
 {
 vec=vecb
 
 } 
 else
 {
 vec=vec}
 
 
 }
 
 T=T*cool
 
 }
 return(vec)
 
 }

Hacemos una prueba de la optimización por este método. Los resultados fueron los siguientes.

#Método de Recorrido Simulado

annealingoptimize(costo.tiempo,T=10000.0,cool=0.95,step=1)
 [1] 8 7 4 7 8 3 7 3 1 5 8 7

> solucionrs=c(8,7,4,7,8,3,7,3,1,5,8,7)

#Evaluamos nuestra función de costo
> costo.tiempo(soucionrs)
[1] 4573
#Imprimimos nuestros vuelos

> impri.horarios(soucionrs)
[1] "Seymour" "BOS" "17:11-18:30 108" "15:07-17:21 129"
[1] "Franny" "DAL" "10:30-14:57 290" "15:07-17:21 129"
[1] "Zooey" "CAK" "17:08-19:08 262" "9:31-11:43 210" 
[1] "Walt" "MIA" "15:34-18:11 326" "9:31-11:43 210" 
[1] "Buddy" "ORD" "6:05-8:32 174" "12:31-14:02 234"
[1] "Les" "OMA" "16:51-19:09 147" "15:07-17:21 129"

Por último hago un pequeño código de un algoritmo genético, este fue adaptado para el problema que se está resolviendo. Lo que se debe de mencionar son las etapas que conforman dicho algoritmo: Mutación, Cruza, Generación de Población, Elitismo y procesamiento.

Cada etapa tienen varias opciones para ser realizada, la mutación puede ser solo una variación o afectación a la cadena de entrada. La cruza puede tener varios métodos para implementarse, y la idea base es “pegar” dos cadenas para obtener una nueva. La Población, tal cual como su nombre es, genera una cantidad de cadenas obtenidas de manera aleatoria y por “elitismo” elegir las mejores, es decir por un criterio se discrimina entre las cadenas.

Después de tener una Población inicial y haber elegido entre ellas,  se usa la mutación y la cruza para ir “mejorando” la población y por ultimo se  aproxima al óptimo.

El código que muestro es una versión simple y hasta cierto punto fea o mala, ya que hay pasos que pueden hacerse de manera más simplificada o compacta. Pero preferí mostrar todos los pasos por simples que fueran, pensando dar visibilidad en lugar de elegancia.

#Algoritmo Genético

genetic.optimize<-function(costf,popsize=50,step=1,mutprod=0.2,elite=0.2,maxiter=100){
 
 #Mutación
 mutacion<-function(vec){
 i=sample(12,size=1)
 a=runif(1,min=0,max=1)
 if((a<0.5)){
 if(i==1){return(c(vec[1]-step,vec[2:12]))}
 else if (i==12){return(c(vec[1:11],vec[12]-step))} 
 else
 {
 return(c(vec[1:(i-1)],vec[i]-step,vec[(i+1):12]))
 }
 }
 else if(vec[i]<9) 
 { 
 if(i==1){return(c(vec[1]+step,vec[2:12]))}
 else if (i==12){return(c(vec[1:11],vec[12]+step))}
 else 
 {
 return(c(vec[1:(i-1)],vec[i]+step,vec[(i+1):12]))
 }
 
 }
 
 }

 #Cruce entre dos cadenas de soluciones
 crossover<-function(r1,r2){
 i=sample(12,size=1) 
 if(i==1){return(c(r1[1],r2[2:12]))}
 else if (i==12){return(c(r1[1:11],r2[12]))}
 else 
 {return(c(r1[1:(i-1)],r2[i:12]))}
 }
 
 #Creación de la población inicial
 A=matrix(nrow=popsize,ncol=12)
 
 for(i in (1:popsize)){
 A[i,]=rbind(sample(9,size=12,replace=TRUE))
 }

 #Ganadores por generación
 topelite=as.integer(elite*popsize)
 
 #Loop principal
 for( i in (1:maxiter)){
 score=sapply(1:nrow(A),function(p){costf(A[p,])})
 B=cbind(A,score)
 Bdata<-data.frame(B)
 
 Bdata<-Bdata[with(Bdata,order(score)),]
 A=as.matrix(Bdata[1:topelite,1:12])
 while (length(A[,1])<popsize){
 a=runif(1,max=1,min=0)
 if(a<mutprod){
 c3=sample(topelite, size=1)
 A=rbind(A,mutacion(A[nrow=c3,1:12]))
 }
 else
 {
 c1=sample(topelite, size=1)
 c2=sample(topelite, size=1)
 A=rbind(A,crossover(A[nrow=c1,1:12],A[nrow=c2,1:12]))
 }
 } 
 
 }
 print(A[1,1:12])
 
 }

Se  compara el resultado y  se obtenemos lo siguiente.

#Método de Algoritmo Genético

 sol5
 6 5 4 5 6 7 4 5 5 7 5 5

costo.tiempo(sol5)
 2898

impri.horarios(sol5)
[1] "Seymour" "BOS" "13:40-15:37 138" "12:31-14:02 234"
[1] "Franny" "DAL" "10:30-14:57 290" "12:31-14:02 234"
[1] "Zooey" "CAK" "13:40-15:38 137" "15:07-17:21 129"
[1] "Walt" "MIA" "11:28-14:40 248" "12:31-14:02 234"
[1] "Buddy" "ORD" "12:44-14:17 134" "15:07-17:21 129"
[1] "Les" "OMA" "12:18-14:56 172" "12:31-14:02 234"

Es muy curioso como los resultados de los algoritmo muestran comportamientos muy distintos en cuando a las elecciones de vuelos y precios. Claro que esto es dependiente de la función de costo definida. El costo por medio de la solución de algoritmos genéticos es mucho mejor que las otras dos. Es casi 1000 unidades menos que la mejor de las otras dos.

Espero esta entrada de una idea de lo que es en general la optimización y el tipo de algoritmos que existen. Por su puesto que hay muchos temas que ni siquiera menciono, pero espero puedan consultarlos y que lo visto en esta entrada les ayude a tener una visión global del tema.

Comentario: El algoritmo que actualmente se emplea para hacer Machine Learning con grandes volúmenes de datos se llama Gradiente Descendiente Estocástico, se ha verificado que es adecuado para la mayoría de algoritmo de clasificación y predicción de ML, lo cual lo ha vuelto popular y se ha implementado en varias librerías, para ver  una breve, pero buena explicación el método recomiendo el vídeo de 12 minutos de Andrew Ng: https://class.coursera.org/ml-003/lecture/107

Referencias.

Libros:

http://www.amazon.com/Optimization-Machine-Learning-Information-Processing/dp/026201646X

http://www.princeton.edu/~sbubeck/Bubeck14.pdf

http://www.amazon.es/Stochastic-Optimization-Scientific-Computation-Schneider/dp/3540345590

Centros de investigación:

http://www.cs.manchester.ac.uk/mlo/

http://research.microsoft.com/en-us/groups/mlo/

Notas, artículos y cursos:

http://stoprog.org/SPTutorial/SPTutorial.php

http://www.people.fas.harvard.edu/~plam/teaching/methods/mcmc/mcmc_print.pdf

http://online.stanford.edu/course/convex-optimization-winter-2014

http://www.econ.uiuc.edu/~roger/research/conopt/coptr.pdf

Anuncios

Un comentario sobre “Optimización

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s