更新可重现的例子来说明我的问题
我最初的问题是“在 R 中实现信任区域反射优化算法”。但是，在生成可重现示例的过程中（感谢@Ben 的建议），我意识到我的问题是在 Matlab 中，一个函数lsqnonlin对于我拥有的大多数情况来说，它是好的（意味着不需要选择一个好的起始值，速度很快），而在 R 中，没有这样一个万能的功能。不同的优化算法在不同的情况下效果很好。不同的算法得出不同的解决方案。这背后的原因可能不是R中的优化算法不如Matlab中的trust-region-reflective算法，也可能与R如何处理自动微分有关。这个问题其实源于两年前的工作中断。当时，包optimx的作者之一 John C Nash 教授有人建议，Matlab 已经为自动微分做了大量工作，这可能是 Matlab lsqnonlin 比 R 中的优化函数/算法执行得更好的原因。我无法用我的知识弄清楚。
下面的示例显示了我遇到的一些问题（更多可重现的示例即将出现）。要运行示例，首先运行install_github("KineticEval","zhenglei-gao"). 您需要安装包mkin及其依赖项，并且可能还需要为不同的优化算法安装一堆其他包。
lsqnonlin基本上，我正在尝试解决 Matlab 函数文档（http://www.mathworks.de/de/help/optim/ug/lsqnonlin.html）中描述的非线性最小二乘曲线拟合问题。在我的例子中，曲线是由一组微分方程建模的。我将通过示例进行更多解释。我尝试过的优化算法包括：
- 来自
nls.lm列文堡-马夸特的马克 - 港口从
nlm.inb - L-BGFS-B 从
optim - spg 来自
optimx solnp包Rsolnp
我也尝试了其他一些，但没有在这里展示。
我的问题总结
- R 中是否有可靠的函数/算法可供使用，例如
lsqnonlin在 Matlab 中可以解决我的非线性最小二乘问题类型？（我找不到一个。） - 对于一个简单的案例，不同的优化达到不同的解决方案的原因是什么？
- 是什么让
lsqnonlinR 中的功能更胜一筹？信任区域反射算法或其他原因？ - 有没有更好的方法来解决我的问题类型，以不同的方式制定？也许有一个非常简单的解决方案，但我只是不知道。
示例 1：一个简单的案例
我会先给出R代码，稍后再解释。
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))
最后一行的输出是：
L-BFGS-B Marq Port spg solnp
5735.744 4714.500 5780.446 5728.361 4714.499
除了“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))
现在最后一行的输出是：
L-BFGS-B Marq Port spg solnp
5653.132 4866.961 5653.070 5635.372 4714.499
我将在这里解释我正在优化的内容。如果您已运行上述脚本并查看曲线，我们将使用具有一阶反应的二室模型来描述曲线。表达模型的微分方程在ex1$diffs：
Parent
"d_Parent = - k_Parent * Parent"
Metab
"d_Metab = - k_Metab * Metab + k_Parent * f_Parent_to_Metab * Parent"
对于这个简单的例子，我们可以从微分方程推导出描述两条曲线的方程。要优化的参数是$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}
$$
因此，我们可以在不求解微分方程的情况下拟合曲线。
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
这样做，经过的时间几乎为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))
这是一个示例，您将看到如下警告消息：
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
原来的问题
Matlab 优化工具箱求解器中使用的许多方法都基于信任区域。根据 CRAN 任务视图页面，只有包trust、trustOptim、minqa 实现了基于信任区域的方法。但是，trust并且trustOptim需要梯度和粗麻布。bobyqa在minqa似乎不是我要找的那个。根据我的个人经验，与我在 R 中尝试的算法相比，Matlab 中的信任区域反射算法通常表现更好。所以我试图在 R 中找到该算法的类似实现。
我在这里问了一个相关的问题：R function to search for a function
Matthew Plourde 提供的答案lsqnonlin在 Matlab 中给出了一个具有相同函数名称的函数，但没有实现 trust-region-reflective 算法。我编辑了旧的并在这里提出了一个新问题，因为我认为 Matthew Plourde 的回答通常对正在寻找功能的 R 用户非常有帮助。
我再次搜索并没有运气。是否还有一些功能/包可以实现类似的 matlab 功能。如果没有，我想知道是否允许我将 Matlab 函数直接翻译成 R 并用于我自己的目的。