En R, ¿cómo hacer una optimización no lineal de mínimos cuadrados que implique resolver ecuaciones diferenciales?

Question

Feb 04, 2013, 02:00 PM

En R, ¿cómo hacer una optimización no lineal de mínimos cuadrados que implique resolver ecuaciones diferenciales?

Actualización con ejemplos reproducibles para ilustrar mi problema.

Mi pregunta original era "Implementación del algoritmo de optimización de confianza de región-reflexión en R". Sin embargo, en la forma de producir un ejemplo reproducible (gracias a @Ben por su consejo), me doy cuenta de que mi problema es que en Matlab, una funciónlsqnonlin es bueno (lo que significa que no es necesario elegir un buen valor de inicio, rápido) lo suficiente para la mayoría de los casos que tengo, mientras que en R, no existe una función única para todos. Diferente algoritmo de optimización funciona bien en diferentes casos. Diferentes algoritmos alcanzan diferentes soluciones. La razón detrás de esto puede no ser que los algoritmos de optimización en R sean inferiores al algoritmo de reflexión de la región de confianza en Matlab, también podría estar relacionado con la forma en que R maneja la diferenciación automática. Este problema viene en realidad del trabajo interrumpido hace dos años. En ese momento, el profesor John C Nash, uno de los autores del paquete.optimx ha sugerido que ha habido bastante trabajo para Matlab para la diferenciación automática, lo cual podría ser la razón por la cual el lsqnonlin de Matlab funciona mejor que las funciones / algoritmos de optimización en R. No puedo resolverlo con mi conocimiento.

El ejemplo a continuación muestra algunos problemas que he encontrado (vienen más ejemplos reproducibles). Para ejecutar los ejemplos, primero ejecuteinstall_github("KineticEval","zhenglei-gao"). Necesitas instalar el paquetemkin y sus dependencias, y también puede necesitar instalar un montón de otros paquetes para diferentes algoritmos de optimización.

Básicamente, estoy tratando de resolver problemas de ajuste de curvas de mínimos cuadrados no lineales como se describe en la función Matlablsqnonlin la documentación (http://www.mathworks.de/de/help/optim/ug/lsqnonlin.html). Las curvas en mi caso están modeladas por un conjunto de ecuaciones diferenciales. Explicaré un poco más con los ejemplos. Algoritmos de optimización que he probado incluyendo:

Marq denls.lm, el Levenburg-MarquardtPuerto denlm.inbL-BGFS-B deoptimspg deoptimxsolnp del paqueteRsolnp

También he intentado algunos otros, pero no se muestra aquí.

Un resumen de mis preguntas.¿Hay en R una función / algoritmo confiable para usar, comolsqnonlin ¿En Matlab que puede resolver mi tipo de problemas de mínimos cuadrados no lineales? (No pude encontrar uno.)¿Cuál es la razón por la que, en un caso simple, diferentes optimizaciones alcanzan diferentes soluciones?Lo que hacelsqnonlin Superior a las funciones en R? ¿El algoritmo de reflexión de región de confianza u otras razones?¿Hay mejores maneras de resolver mi tipo de problemas, formular de manera diferente? Tal vez haya una solución mucho más simple pero simplemente no lo sé.Ejemplo 1: Un caso simple

Daré los códigos R primero y luego explicaré.

ex1 <- mkinmod.full(
  Parent = list(type = "SFO", to = "Metab", sink = TRUE,
                k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
                M0 = list(ini = 195, fixed = 0,lower = 0,upper = Inf),
                FF = list(ini = c(.1),fixed = c(0),lower = c(0),upper = c(1)),
                time=c(0.0,2.8,   6.2,  12.0,  29.2,  66.8,  99.8, 127.5, 154.4, 229.9, 272.3, 288.1, 322.9),
                residue = c( 157.3, 206.3, 181.4, 223.0, 163.2, 144.7,85.0, 76.5, 76.4, 51.5, 45.5, 47.3, 42.7),
                weight = c( 1,  1,   1, 1, 1,   1,  1,     1,     1,     1,     1,     1,     1)),
  Metab = list(type = "SFO",
               k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
              M0 = list(ini = 0, fixed = 1,lower = 0,upper = Inf),
                    residue =c( 0.0,  0.0,  0.0,  1.6,  4.0, 12.3, 13.5, 12.7, 11.4, 11.6, 10.9,  9.5,  7.6),
               weight = c( 1,  1,   1, 1, 1,   1,  1,     1,     1,     1,     1,     1,     1))
  )
ex1$diffs
Fit <- NULL
alglist <- c("L-BFGS-B","Marq", "Port","spg","solnp")
for(i in 1:5) {
  Fit[[i]] <- mkinfit.full(ex1,plot = TRUE, quiet= TRUE,ctr = kingui.control(method = alglist[i],submethod = 'Port',maxIter = 100,tolerance = 1E-06, odesolver = 'lsoda'))
  }
names(Fit) <- alglist
kinplot(Fit[[2]])
(lapply(Fit, function(x) x$par))
unlist(lapply(Fit, function(x) x$ssr))

La salida de la última línea es:

L-BFGS-B     Marq     Port      spg    solnp 
5735.744 4714.500 5780.446 5728.361 4714.499

Excepto por "Marq" y "solnp", los otros algoritmos no alcanzaron el método óptimo. Además, el método 'spg' (también otros métodos como 'bobyqa') necesita demasiadas evaluaciones de funciones para un caso tan simple. Además, si cambio el valor de inicio y hagok_Parent=0.0058 (el valor óptimo para ese parámetro) en lugar del aleatorio elegido0.1, "Marq" ya no puede encontrar el óptimo! (Código provisto abajo). También he tenido conjuntos de datos donde "solnp" no encuentra el óptimo. Sin embargo, si usolsqnonlin En Matlab, no he encontrado ninguna dificultad para casos tan simples.

ex1_a <- mkinmod.full(
  Parent = list(type = "SFO", to = "Metab", sink = TRUE,
                k = list(ini = 0.0058,fixed = 0,lower = 0,upper = Inf),
                M0 = list(ini = 195, fixed = 0,lower = 0,upper = Inf),
                FF = list(ini = c(.1),fixed = c(0),lower = c(0),upper = c(1)),
                time=c(0.0,2.8,   6.2,  12.0,  29.2,  66.8,  99.8, 127.5, 154.4, 229.9, 272.3, 288.1, 322.9),
                residue = c( 157.3, 206.3, 181.4, 223.0, 163.2, 144.7,85.0, 76.5, 76.4, 51.5, 45.5, 47.3, 42.7),
                weight = c( 1,  1,   1, 1, 1,   1,  1,     1,     1,     1,     1,     1,     1)),
  Metab = list(type = "SFO",
               k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
              M0 = list(ini = 0, fixed = 1,lower = 0,upper = Inf),
                    residue =c( 0.0,  0.0,  0.0,  1.6,  4.0, 12.3, 13.5, 12.7, 11.4, 11.6, 10.9,  9.5,  7.6),
               weight = c( 1,  1,   1, 1, 1,   1,  1,     1,     1,     1,     1,     1,     1))
  )

Fit_a <- NULL
alglist <- c("L-BFGS-B","Marq", "Port","spg","solnp")
for(i in 1:5) {
  Fit_a[[i]] <- mkinfit.full(ex1_a,plot = TRUE, quiet= TRUE,ctr = kingui.control(method = alglist[i],submethod = 'Port',maxIter = 100,tolerance = 1E-06, odesolver = 'lsoda'))
  }
names(Fit_a) <- alglist
lapply(Fit_a, function(x) x$par)
unlist(lapply(Fit_a, function(x) x$ssr))

Ahora la salida de la última línea es:

L-BFGS-B     Marq     Port      spg    solnp 
5653.132 4866.961 5653.070 5635.372 4714.499

Explicaré lo que estoy optimizando aquí. Si ha ejecutado el script anterior y ve las curvas, utilizamos un modelo de dos compartimentos con reacciones de primer orden para describir las curvas. Las ecuaciones diferenciales para expresar el modelo están enex1$diffs:

                                                             Parent 
                                    "d_Parent = - k_Parent * Parent" 
                                                               Metab 
"d_Metab = - k_Metab * Metab + k_Parent * f_Parent_to_Metab * Parent"

Para este caso simple, de las ecuaciones diferenciales podemos derivar las ecuaciones para describir las dos curvas. Los parámetros a optimizar son $M_0,k_p, k_m, c=\mbox{FF_parent_to_Met} $ con las restricciones $M_0>0,k_p>0, k_m>0, 1> c >0$ .

$
\begin{split}
            y_{1j}&= M_0e^{-k_pt_i}+\epsilon_{1j}\\
            y_{2j} &= cM_0k_p\frac{e^{-k_mt_i}-e^{-k_pt_i}}{k_p-k_m}+\epsilon_{2j}
            \end{split}
$

Por lo tanto podemos ajustar la curva sin resolver ecuaciones diferenciales.

BCS1.l <- mkin_wide_to_long(BCS1)
BCS1.l <- na.omit(BCS1.l)
indi <- c(rep(1,sum(BCS1.l$name=='Parent')),rep(0,sum(BCS1.l$name=='Metab')))
sysequ.indi <- function(t,indi,M0,kp,km,C)
  {
    y <- indi*M0*exp(-kp*t)+(1-indi)*C*M0*kp/(kp-km)*(exp(-km*t)-exp(-kp*t));
    y
  }
M00 <- 100
kp0 <- 0.1
km0 <- 0.01
C0 <- 0.1
library(nlme)
result1 <- gnls(value ~ sysequ.indi(time,indi,M0,kp,km,C),data=BCS1.l,start=list(M0=M00,kp=kp0,km=km0,C=C0),control=gnlsControl())
#result3 <- gnls(value ~ sysequ.indi(time,indi,M0,kp,km,C),data=BCS1.l,start=list(M0=M00,kp=kp0,km=km0,C=C0),weights = varIdent(form=~1|name))
## Coefficients:
##          M0           kp           km            C 
## 1.946170e+02 5.800074e-03 8.404269e-03 2.208788e-01

De esta manera, el tiempo transcurrido es casi 0, y se alcanza el óptimo. Sin embargo, no siempre tenemos este caso simple. El modelo puede ser complejo y es necesario resolver las ecuaciones diferenciales. Ver ejemplo 2

Ejemplo 2, un modelo complejo.

Trabajé en este conjunto de datos hace mucho tiempo y no tengo tiempo para terminar de ejecutar el siguiente script. (Es posible que necesites horas para terminar de correr).

data(BCS2)
ex2 <- mkinmod.full(Parent= list(type = "SFO",to = c( "Met1", "Met2","Met4", "Met5"),
                                 k = list(ini = 0.1,fixed = 0,lower = 0,upper = Inf),
                                 M0 = list(ini = 100,fixed = 0,lower = 0,upper = Inf),
                                 FF = list(ini = c(.1,.1,.1,.1),fixed = c(0,0,0,0),lower = c(0,0,0,0),upper = c(1,1,1,1))),
                    Met1 = list(type = "SFO",to = c("Met3", "Met4")),
                    Met2 = list(type = "SFO",to = c("Met3")),
                    Met3 = list(type = "SFO" ),
                    Met4 = list(type = "SFO", to = c("Met5")),
                    Met5 = list(type = "SFO"),
                    data=BCS2)
ex2$diffs
Fit2 <- NULL
alglist <- c("L-BFGS-B","Marq", "Port","spg","solnp")
for(i in 1:5) {
  Fit2[[i]] <- mkinfit.full(ex2,plot = TRUE, quiet= TRUE,ctr = kingui.control(method = alglist[i],submethod = 'Port',maxIter = 100,tolerance = 1E-06, odesolver = 'lsoda'))
  }
names(Fit) <- alglist
(lapply(Fit, function(x) x$par))
unlist(lapply(Fit, function(x) x$ssr))

Este es un ejemplo donde verá mensajes de advertencia como:

DLSODA-  At T (=R1) and step size H (=R2), the    
  corrector convergence failed repeatedly     
  or with ABS(H) = HMIN   
In above message, R = 
[1] 0.000000e+00 2.289412e-09

La pregunta original

Muchos de los métodos utilizados en los solucionadores de Matlab Optimization Toolbox se basan en regiones de confianza. De acuerdo con la página de Vista de Tareas CRAN, solo paquetesconfianza, trustOptim, minqa tener implementados los métodos basados en la región de confianza. Sin embargo,trust ytrustOptim requiere gradiente y arpillera.bobyqa enminqa Parece que no es el que estoy buscando. Desde mi experiencia personal, el algoritmo de reflexión de la región de confianza en Matlab a menudo funciona mejor en comparación con los algoritmos que probé en R. Así que traté de encontrar una implementación similar de este algoritmo en R.

He hecho una pregunta relacionada aquí:Función R para buscar una función

La respuesta proporcionada por Mateo Plourde da una función.lsqnonlin con el mismo nombre de función en Matlab pero no tiene implementado el algoritmo de reflexión de región de confianza. Edité el anterior y hice una nueva pregunta aquí porque creo que la respuesta de Matthew Plourde es en general muy útil para los usuarios de R que buscan una función.

Hice una búsqueda de nuevo y no tengo suerte. ¿Todavía hay algunas funciones / paquetes que implementan funciones similares de matlab? Si no, me pregunto si está permitido que traduzca la función de Matlab directamente a R y la use para mi propio propósito.