那些遇到过的问题

在 R 中，如何进行涉及求解微分方程的非线性最小二乘优化？

Zhe*_*lei 5 optimization r nonlinear-optimization

使用可重现的示例进行更新以说明我的问题

我最初的问题是“在 R 中实现信任区域反射优化算法”。但是，在生成可重现示例的过程中（感谢 @Ben 的建议），我意识到我的问题是在 Matlab 中，一个函数lsqnonlin很好（意味着不需要选择一个好的起始值，速度快）对于我拥有的大多数情况来说已经足够了，而在 R 中，没有这样一个一应俱全的功能。不同的优化算法在不同的情况下效果很好。不同的算法达到不同的解决方案。这背后的原因可能不是R中的优化算法不如Matlab中的信任区域反射算法，也可能与R如何处理自动微分有关。这个问题其实来自两年前的工作中断。当时，包optimx的作者之一 John C Nash 教授 已经建议 Matlab 进行了大量的自动微分工作，这可能是 Matlab lsqnonlin 比 R 中的优化函数/算法执行得更好的原因。我无法用我的知识弄清楚。

下面的示例显示了我遇到的一些问题（更多可重现的示例即将推出）。要运行示例，首先运行install_github("KineticEval","zhenglei-gao"). 您需要安装包mkin及其依赖项，并且可能还需要为不同的优化算法安装一堆其他包。

基本上我试图解决非线性最小二乘曲线拟合问题，如 Matlab 函数lsqnonlin的文档 ( http://www.mathworks.de/de/help/optim/ug/lsqnonlin.html ) 中所述。在我的例子中，曲线是由一组微分方程建模的。我将通过示例进行更多解释。我尝试过的优化算法包括：

  • Marq from nls.lm, the Levenburg-Marquardt
  • 端口来自 nlm.inb
  • L-BGFS-B 来自 optim
  • 来自 optimx
  • solnp包Rsolnp

我也尝试了其他一些，但没有在这里显示。

我的问题总结

  • R 中是否有可靠的函数/算法可以使用，例如lsqnonlin在 Matlab 中可以解决我的非线性最小二乘问题类型？（我找不到一个。）
  • 对于一个简单的情况，不同的优化达到不同的解决方案的原因是什么？
  • 什么lsqnonlin比 R 中的函数更胜一筹？信任区域反射算法或其他原因？
  • 有没有更好的方法来解决我的问题类型，以不同的方式表达？也许有一个非常简单的解决方案，但我只是不知道。

示例 1：一个简单的案例

我先给出R代码，稍后解释。 示例 1 的拟合图

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))
Run Code Online (Sandbox Code Playgroud)

最后一行的输出是：

L-BFGS-B     Marq     Port      spg    solnp 
5735.744 4714.500 5780.446 5728.361 4714.499
Run Code Online (Sandbox Code Playgroud)

除了“Marq”和“solnp”，其他算法都没有达到最优。此外，'spg'方法（还有'bobyqa'等其他方法）对于这样一个简单的情况需要太多的函数评估。此外，如果我更改起始值并使用k_Parent=0.0058（该参数的最佳值）而不是随机选择的值0.1，“Marq”将无法再找到最佳值！（下面提供了代码）。我也有“solnp”找不到最佳数据的数据集。但是，如果我lsqnonlin在Matlab中使用，对于这种简单的情况，我没有遇到任何困难。

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))
Run Code Online (Sandbox Code Playgroud)

现在最后一行的输出是：

L-BFGS-B     Marq     Port      spg    solnp 
5653.132 4866.961 5653.070 5635.372 4714.499
Run Code Online (Sandbox Code Playgroud)

我将在这里解释我正在优化的内容。如果您已运行上述脚本并查看曲线，我们将使用具有一阶反应的两室模型来描述曲线。表达模型的微分方程为ex1$diffs：

                                                             Parent 
                                    "d_Parent = - k_Parent * Parent" 
                                                               Metab 
"d_Metab = - k_Metab * Metab + k_Parent * f_Parent_to_Metab * Parent"
Run Code Online (Sandbox Code Playgroud)

对于这个简单的情况，我们可以从微分方程推导出描述两条曲线的方程。待优化参数有$M_0,k_p, k_m, c=\mbox{FF_parent_to_Met} $约束条件$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}
$$
Run Code Online (Sandbox Code Playgroud)

因此，我们可以在不求解微分方程的情况下拟合曲线。 

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
Run Code Online (Sandbox Code Playgroud)

这样做，经过的时间几乎为0，达到了最佳状态。然而，我们并不总是有这种简单的情况。该模型可能很复杂，需要求解微分方程。参见示例 2

示例 2，复杂模型

我很久以前就研究过这个数据集，但没有时间自己完成以下脚本的运行。（您可能需要数小时才能完成跑步。） 拟合图示例

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))
Run Code Online (Sandbox Code Playgroud)

这是一个示例，您将在其中看到如下警告消息：

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
Run Code Online (Sandbox Code Playgroud)

原来的问题

Matlab Optimization Toolbox 求解器中使用的许多方法都基于信任区域。根据 CRAN 任务视图页面，只有包trust，trustOptim，minqa 实现了基于信任区域的方法。但是，trust并且trustOptim需要梯度和粗麻布。bobyqa在minqa似乎不是我要找的那个。根据我的个人经验，与我在 R 中尝试的算法相比，Matlab 中的 trust-region-reflective 算法通常性能更好。 所以我试图在 R 中找到该算法的类似实现。

我在这里问了一个相关的问题：R function to search for a function

Matthew Plourde 提供的答案lsqnonlin在 Matlab 中给出了一个具有相同函数名称的函数，但没有实现信任区域反射算法。我编辑了旧的并在这里提出了一个新问题，因为我认为 Matthew Plourde 的回答通常对正在寻找功能的 R 用户非常有帮助。

我再次进行了搜索，但没有运气。是否还有一些函数/包实现了类似的 matlab 函数。如果没有，我想知道是否允许我将 Matlab 函数直接翻译成 R 并将其用于我自己的目的。

Joh*_*nke 1

一般来说，当只查看问题的标题时，我建议仅使用该FME包。但这不是您问题的重点，成功可能取决于您建立模型的方式。

对于您在示例中显示的问题类型（使用多个转换产品拟合退化数据），我创建了该mkin包作为 FME 此类问题的便捷包装器。那么让我们看看 mkin 0.9-29 在这些情况下的表现如何。使用mkin，只能使用FME提供的算法：

实施例1

library(mkin)

ex1_data_wide = data.frame(
  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),
  Parent = 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),
  Metab = 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))

ex1_data = mkin_wide_to_long(ex1_data_wide, time = "time")

ex1_model = mkinmod(Parent = list(type = "SFO", to = "Metab"),
                    Metab = list(type = "SFO"))

algs = c("L-BFGS-B", "Marq", "Port")

times_ex1 <- list()
fits_ex1 <- list()
for (alg in algs) {
  times_ex1[[alg]] <- system.time(fits_ex1[[alg]] <- mkinfit(ex1_model, ex1_data,
                                                             method.modFit = alg))
}

times_ex1
unlist(lapply(fits_ex1, function(x) x$ssr))
Run Code Online (Sandbox Code Playgroud)

因此，nls.lm 中的 Levenberg-Marquardt 和 Port 算法都可以找到最小值，而 LM 速度要快得多：

$`L-BFGS-B`
       User      System verstrichen 
      2.036       0.000       2.051 

$Marq
       User      System verstrichen 
      0.716       0.000       0.714 

$Port
       User      System verstrichen 
      2.032       0.000       2.030 

L-BFGS-B     Marq     Port 
5742.312 4714.498 4714.498
Run Code Online (Sandbox Code Playgroud)

当我告诉 mkin 使用形成分数而不仅仅是比率时

ex1_model = mkinmod(Parent = list(type = "SFO", to = "Metab"),
                    Metab = list(type = "SFO"), use_of_ff = "max")
Run Code Online (Sandbox Code Playgroud)

并使用你的起始值，

for (alg in algs) {
  times_ex1[[alg]] <- system.time(fits_ex1[[alg]] <- mkinfit(ex1_model, ex1_data,
    state.ini = c(195, 0),
    parms.ini = c(f_Parent_to_Metab = 0.1, k_Parent = 0.0058, k_Metab = 0.1),
    method.modFit = alg))
}
Run Code Online (Sandbox Code Playgroud)

所有三种算法都找到相同的解决方案，甚至更快。但是，如果我随后在调用 mkinfit ( transform_rates = FALSE, transform_fractions = FALSE) 时关闭比率和分数的转换，我会得到

L-BFGS-B     Marq     Port 
5653.132 4714.498 5653.070
Run Code Online (Sandbox Code Playgroud)

所以它似乎与参数内部转换的方式有关（当你给出边界时，FME 也会这样做）。在 mkin 中，我进行了显式的内部参数转换，因此使用默认设置的优化参数不需要限制。

实施例2

library(mkin)
library(KineticEval) # for the dataset BCS2
data(BCS2)

ex2_data = mkin_wide_to_long(BCS2, time = "time")

ex2_model = mkinmod(Parent = list(type = "SFO", to = paste0("Met", 1:5)),
                    Met1 = list(type = "SFO", to = c("Met3", "Met4")),
                    Met2 = list(type = "SFO", to = "Met3"),
                    Met3 = list(type = "SFO"),
                    Met4 = list(type = "SFO", to = "Met5"),
                    Met5 = list(type = "SFO"))

times_ex2 <- list()
fits_ex2 <- list()

for (alg in algs) {
  times_ex2[[alg]] <- system.time(fits_ex2[[alg]] <- mkinfit(ex2_model, ex2_data,
    method.modFit = alg))
}   

times_ex2
unlist(lapply(fits_ex2, function(x) x$ssr))
Run Code Online (Sandbox Code Playgroud)

同样，LM 是最快的，但最低的最小值是由 Port 找到的：

$`L-BFGS-B`
       User      System verstrichen 
     75.728       0.004      75.653 

$Marq
       User      System verstrichen 
      6.440       0.004       6.436 

$Port
       User      System verstrichen 
     51.200       0.028      51.180 

L-BFGS-B     Marq     Port 
485.3099 572.9635 478.4379
Run Code Online (Sandbox Code Playgroud)

我过去总是推荐 LM，但最近我也发现它有时会陷入局部最小值，具体取决于不明确参数的起始值。一个例子是 Schaefer 07 数据，在 mkin 包中mkinfit 的最后一个单元测试中处理，称为test.mkinfit.schaefer07_complex_example。

希望这有用，亲切的问候，

约翰内斯

PS：当我注意到您在 github 上的 KineticEval 包中添加了信任区域反射优化的纯 R 实现作为函数 lsqnonlin() 时，我发现了这个问题，并且我正在对信任区域反射进行搜索。

归档时间：

查看次数：

3945 次

最近记录：

11 年，2 月 前
R函数搜索函数 3
更多相关链接
相关归档
是否可以在Emacs ESS中获得类似于Rstudio中可用的R的代码完成？ 34
如何在R中舍入到整数？ 23
data.table join然后将列添加到现有data.frame而不重新复制 19
R:*_join(dplyr)的标准评估 13
如何加快随机森林的训练？ 11
Kniter kable(Sweave)中的LaTeX数学表达式 11
scipy.optimize minus：目标函数中的两个输出变量？ 5
列表搜索优化 4
执行2个elif块之间的语句 1
使用python scipy.optimize.minimize进行优化 0
难疑归档
C++中的" - >"运算符是什么？ 8590
如果__name__ =="__ main__":怎么办？ 5545
什么是反思,为什么它有用？ 2011
为什么Python 3中的"1000000000000000在范围内(1000000000000001)"如此之快？ 1890
如何在没有换行或空格的情况下打印？ 1760
禁用Chrome缓存以进行网站开发 1563
如何检查Bash中是否设置了变量？ 1417
如何迭代Bash中变量定义的一系列数字？ 1409
丢弃Git中的本地提交 1304
如何删除选择框的所有选项,然后添加一个选项并使用jQuery选择它？ 1028