Un artículo de métodos numéricos publicado en Physical Review Letters

Los que trabajamos (investigamos) en métodos numéricos y somos físicos de formación tenemos una espinita clavada, publicar en Physical Review Letters. Siempre uno quiere creerse que hace física computacional y que algún día podrá publicar en PRL, pero ojeando los índices de artículos de esta revista cada semana… y sabiendo que publicar en PRL cada día es más difícil… Pero la esperanza es lo último que se pierde y ver un artículo de métodos numéricos en PRL le da a uno una gran alegría. Nunca hubiera creído que el artículo de los físicos Mark K. Transtrum, Benjamin B. Machta y James P. Sethna, de la Universidad de Cornell, Ithaca, New York, titulado «Why are Nonlinear Fits to Data so Challenging?,» que ya ojeé hace meses en ArXiv, y al que estuve a punto de dedicar una entrada en este blog, haya acabado siendo publicado en Physical Review Letters 104: 060201, 12 Feb. 2010. ¡Enhorabuena!

El artículo estudia el ajuste de los parámetros de un modelo a partir de datos experimentales. Este problema se puede resolver fácilmente con Mathematica, Matlab y cualquier otro software científico al uso. Sin embargo, conforme el número de parámetros crece las dificultades se acumulan y los métodos numéricos convencionales que utilizan dichos programas encuentran dificultades de convergencia. Yo recuerdo una vez que estuve ajustando los parámetros de un modelo de la respiración de monos bajo estrés. Matlab (usando la variante del algoritmo de Levenberg-Marquardt desarrollado por Powell) necesitó varios días de trabajo (hace ya unos 10 años) y el resultado al final no me convenció. Ajustar, ajustaba, pero los valores de los parámetros eran poco «razonables» y no acabó de gustarme el resultado. Seguramente el problema tenía múltiples soluciones y había una solución que se me escapaba que ajustaba los datos y parecía «razonable». La verdad sea dicha, tras cierto tiempo dedicado al asunto, preferí abandonarlo sin publicarlo. Los que sois investigadores sabéis que eso nos pasa muchas veces. Y cuando el tópico es colateral (yo estudio propagación de ondas en medios no lineales dispersivos), uno lo abandona sin remordimiento alguno.

El artículo de Transtrum et al. trata de explicar la razón por la que muchos algoritmos de optimización utilizados para ajustar parámetros de un modelo encuentran soluciones no razonables. Según ellos, es debido a que entran en una región del espacio de parámetros en la que el modelo no responde a cambios en dichos parámetros. La única opción para escapar de dichas regiones «malas» es ajustar a mano los parámetros (es decir, reparametrizar el problema con inteligencia e intuición). Los autores del artículo tratan de explicar este fenómeno aludiendo a la teoría de la interpolación para la aproximación de funciones, que subyace a todos los algoritmos de optimización. Transtrum et al. introducen una noción de curvatura (extrínseca) en el espacio de parámetros del modelo (una idea de geometría diferencial poco explotada en análisis numérico) que permite identificar estas regiones de búsqueda «malas» y la utilizan para mejorar la convergencia del algoritmo de Levenberg-Marquardt, permitiéndole salir de las mismas.

Los lectores de este blog que no hayan oido hablar del algoritmo de Levenberg-Marquardt y que nunca hayan optimizado los parámetros de un modelo a partir de datos experimentales seguramente no sabrán lo importante que es esta tarea numérica en la práctica diaria de los físicos experimentales, biólogos, economistas, etc. El mejor ajuste suele ser de tipo mínimocuadrático, se minimiza cierta función coste definida positiva, normalmente la suma de los cuadrados de las desviaciones entre el modelo y los datos experimentales. Muchos modelos tienen parámetros que dependen de forma no lineal, con lo que este proceso de minimización requiere métodos numéricos para la resolución de sistemas no lineales de ecuaciones. Encontrar la solución óptima es difícil porque el espacio de parámetros (posibles modelos) está repleto de mínimos locales que no corresponden al mejor ajuste posible (el mínimo global). Más aún, la función coste puede ser muy sensible a ciertos valores de los parámetros y muy poco sensible a otros. La idea de Transtrum et al. es que en el espacio de los parámetros hay una hipersuperficie (variedad diferencial) en la que se encuentra la solución óptima que puede pasar desapercibida al algoritmo de optimización porque tiene un grosor muy delgado para que la pueda resolver correctamente, le llaman «hipercinta» (hiper-ribbon). Para poder detectar la presencia de este tipo de «hipercintas» proponen el uso de una noción de curvatura basada en asumir que en cada paso el método numérico recorre una geodésica en el espacio de parámetros a cierta «velocidad» y «aceleración» geodésicas. La curvatura se calcula gracias a estas segundas derivadas y permite acelerar la convergencia de un método numérico (lo aplican sólo al algoritmo de Levenberg-Marquardt, LM).

¿Qué speedup obtienen con el nuevo algoritmo? Para el ejemplo que incluyen en el artículo (un problema de biología de sistemas basado en una red metabólica modelada con leyes de Michaelis-Menten) el algoritmo LM acelerado obtiene la respuesta óptima el 65% de las veces, cuando el algoritmo LM sin acelerar la obtiene sólo el 33%. ¡Un gran logro! Bueno, no tanto, calcular la solución en 2 horas en lugar de en 4 horas no parece un gran avance científico. Pero … me corroe la envidia.

¿Enviaré mi próximo artículo numérico a PRL a ver si cuela? No sé, no sé, me lo pensaré… pero tengo que venderlo bien ya que la respuesta que espero es que Physical Review Letters publica física no métodos numéricos. Bueno… casi siempre.