Optimización

Esta es la función de Rosenbrock, también conocida como función plátano o —en algunos contextos— como el coco:

Es una de esas funciones contra la que tienen que demostrar su valía los algoritmos de optimización que los matemáticos discurren por ahí. La función ilustra uno de los problemas habituales de la optimización: las variables se confabulan para que las ideas simples no funcionen: los gradientes no apuntan hacia el mínimo, este se encuentra en un valle estrecho, etc. Y que conste que las he visto peores en la práctica.

He escrito esta entrada como una introducción a lo que se cuenta aquí, aquí y aquí sobre el asunto de la relación entre la optimización (como parte del proceso de ajuste de modelos) y la generalización (o su capacidad para aprender sobre el mundo y no solo sobre los datos de entrenamiento). En los enlaces, el lector encontrará planteadas una serie de cuestiones sobre cómo y por qué generalizan los (o cierto tipo de) modelos en lugar de, simplemente, no hacerlo.

Vale, admito que no funciona siempre. Pero una manera de distinguir a un matemático de un ingeniero es por una casi imperceptible pausa que los primeros realizan antes de pronunciar optimización. Un matemático nunca conjuga el verbo optimizar en vano.

[Una vez, hace tiempo, movido por una mezcla de paternalismo y maldad, delegué un subproblema que incluía el fatídico optim de R en una ingeniera. Aún le debe doler el asunto.]

Una de los proyectos en los que estoy trabajando últimamente está relacionado con un problema de optimización no lineal: tengo un modelo (o una familia de modelos) no lineales con una serie de parámetros, unos datos y se trata de lo que no mercería más explicación: encontrar los que minimizan cierta función de error.

Tengo implementadas dos vías:

• La nls, que usa un optimizador numérico genérico para encontrar esos mínimos. (Nótese que uso nls y no nls porque esa función me queda muy corta).
• La stan, donde especifico el modelo, introduzco una serie de prioris más o menos informativas según lo que sepa de mi problema y estimo la distribución a posteriori de mis parámetros.

Ambas tienen sus ventajas y desventajas. La una es rápida y la otra no; la una me da poca información sobre los parámetros y la otra, mucha; una me permite introducir mucha información a priori y la otra casi nada, etc.

A veces tomas un artículo de vaya uno a saber qué disciplina, sismología, p.e., y no dejas de pensar: los métodos estadísticos que usa esta gente son de hace 50 años. Luego cabe preguntarse: ¿pasará lo mismo en estadística con respecto a otras disciplinas?

Por razones que no vienen al caso, me he visto en la tesitura de tener que encontrar mínimos de funciones que podrían cuasicatalogarse como de mínimos cuadrados no lineales. Y por algún motivo, pareciere que no hubiese en el mundo un algoritmo de ajuste que no fuese IRLS. Que tiene una gran tradición en estadística; es, de hecho, la base de la optimización propuesta por Nelder y McCullagh en 1972.

In illo tempore me llamaba mucho la atención encontrar métodos de ciencia de datos basados en factorización de matrices cuando la matriz a factorizar tenía nulos. Ocurre, por ejemplo, en sistemas de recomendación (cuando un usuario no ha visto o no nos ha dicho si le gusta determinada película).

Y claro, con un nulo en la cosa, te comes los apuntes de álgebra lineal con papas.

¿Cómo se hace? Si buscas $latex U$ y $latex V$ tales que $latex Y = UV^\prime$:

Dependiendo de con quién hables, la optimización (de funciones) es un problema fácil o difícil.

Si hablas con matemáticos y gente de la escuela de optim y derivados (BFGS y todas esas cosas), te contarán una historia de terror.

Si hablas con otro tipo de gente, la de los que opinan que el gradiente es un tobogán que te conduce amenamente al óptimo, el de la optimización no alcanza siquiera talla de problema.

Alguien que igual no me lee (porque está de vacaciones) está aprendiendo a punta de palo a manejar la función optim de R con una función objetivo de las enrevesadas. Creo que ahora entiende por qué a los matemáticos nos sobrecoge la palabra optimizar y tratamos de no mencionarla en vano.

Por eso, además, y aunque no estilo envolver meros enlaces de terceros en entradas hechas y derechas, voy a hacer una excepción con An Interactive Tutorial on Numerical Optimization, que no tiene desperdicio.

Durante la charla de Carlos Ortega del pasado jueves sobre el paquete caret y sus concomitancias, se planteó el asunto de la optimización de los parámetros de un modelo usando rejillas (grids) de búsqueda.

Cuando un determinado algoritmo depende de, p.e., cuatro parámetros, se puede definir una rejilla como en

gbmGrid <-  expand.grid(interaction.depth = c(1, 5, 9),
      n.trees = (1:30)*50,
      shrinkage = 0.1,
      n.minobsinnode = 20)

y caret se encarga de ajustar el modelo bajo todas esas combinaciones de parámetros (90 en el ejemplo) para ver cuál de ellas es, con las debidas salvedades, óptima.

Tengo por delante otro proyecto que tiene mucho de análisis exploratorio de datos. Sospecho que más de un árbol construiré. Los árboles son como la Wikipedia: prácticamente nunca el último pero casi siempre el primer recurso.

Esta vez, además, por entretenerme un poco, probaré el paquete [evtree](http://cran.r-project.org/web/packages/evtree/index.html). Aunque no porque espere sorprendentes mejoras con respecto a los tradicionales, ctree y rpart.

¿Qué tiene aquél que los diferencie de los otros dos? Que la optimización es global. Tanto ctree como rpart utilizan algoritmos recursivos: al definir un nuevo corte del espacio, el algoritmo solo tiene en cuenta la región definida por los cortes anteriores. La optimización es local. evtree utiliza un algoritmo global de la familia de los evolucionarios (¡qué tufillo a lentorro!). Los detalles están aquí.

Avisé en mi entrada del otro día: no me preguntéis por qué (imponer restricciones en un problema de mínimos cuadrados).

Pero cuanto más pienso sobre ello, menos claro lo tengo. ¿Por qué restricciones?

Primero, el contexto. O el casi contexto. Porque no es exactamente así. Pero sí parecido. Supongamos que queremos predecir algo y construimos, p.e., 4 modelos. Se nos ocurre (y hay buenas razones para ello) combinar los predictores.

Uno puede pensar en usar la media de las predicciones. O la mediana. O tratar de usar un peso revelado por los datos.

Sí, había restricciones. No me preguntéis por qué, pero los coeficientes tenían que ser positivos y sumar uno. Es decir, buscaba la combinación convexa de cuatro vectores que más se aproximase a y en alguna métrica razonable. Y lo resolví así:

# prepare constrained optimization

y <- dat.clean$actual
x <- t(dat.clean[,2:5])

# target function: L2 first, then other metrics

L2 <- function(coef){
  sum(abs((y - colSums(x * coef)))^1.5)
}

# restrictions: coefs > 0, sum(coefs) ~ 1

ui <- rbind(diag(4), c(-1,-1,-1,-1), c(1,1,1,1))
ci <- c(0,0,0,0,-1.000001,0.999999)

theta <- rep(0.25, 4)

best.coef <- constrOptim(theta, L2,
  grad = NULL, ui = ui, ci = ci)

coefs <- best.coef$par

Objetos aparte de x e y, hay:

Esta entrada tiene un notable contenido publicitario: el pasado viernes se inauguró Martina Cocina (nota: la página es provisional), una cafetería-restaurante en el que ostento una participación del 50%.

Se trata de un lugar al que, por supuesto, mis lectores están invitados a descubrir en número 11 de la castiza plaza de Cascorro de Madrid.

Pero como estoy en un foro al que me obligo a no traer temas extemporáneos, confesaré a mis lectores que tengo dos proyectos en mente que guardan relación tanto con Martina Cocina como con los asuntos que esperan encontrar en estas páginas. El primero es que (si me dejan) trataré de aplicar lo que Dan Arieli nos enseña en Predictably Irrational a, principalmente aunque no exclusivamente, el diseño de la carta y los precios de los productos para ayudaros a vosotros, los clientes a elegir, ejem, más convenientemente (¿para quién?).