24

我前段时间问过这个问题。我不确定是否应该将此作为答案或新问题发布。我没有答案，但我通过nls.lm在 R 中应用 Levenberg-Marquardt 算法“解决”了这个问题，当解决方案在边界处时，我运行 trust-region-reflective 算法（TRR，在 R 中实现）步骤远离它。现在我有新的问题。

根据我的经验，这样做程序会达到最佳状态，并且对起始值不那么敏感。但这只是一种实用的方法，可以避开我在使用nls.lmR 中遇到的问题以及其他优化函数。我想知道为什么nls.lm在具有边界约束的优化问题上表现得这样，以及nls.lm在实践中使用时如何处理边界约束.

下面我举了一个例子来说明这两个问题nls.lm。

  1. 它对起始值很敏感。
  2. 当某个参数到达边界时它会停止。

一个可重现的例子：焦点数据集 D

library(devtools)
install_github("KineticEval","zhenglei-gao")
library(KineticEval)
data(FOCUS2006D)
km <- mkinmod.full(parent=list(type="SFO",M0 = list(ini   = 0.1,fixed = 0,lower = 0.0,upper =Inf),to="m1"),m1=list(type="SFO"),data=FOCUS2006D)
system.time(Fit.TRR <- KinEval(km,evalMethod = 'NLLS',optimMethod = 'TRR'))
system.time(Fit.LM <- KinEval(km,evalMethod = 'NLLS',optimMethod = 'LM',ctr=kingui.control(runTRR=FALSE)))
compare_multi_kinmod(km,rbind(Fit.TRR$par,Fit.LM$par))
dev.print(jpeg,"LMvsTRR.jpeg",width=480)

LM 贴合与 TRR 贴合

描述模型/系统的微分方程是：

"d_parent = - k_parent * parent"                                                 
"d_m1 = - k_m1 * m1 + k_parent * f_parent_to_m1 * parent"

左图是带有初始值的模型，中间是使用“TRR”的拟合模型（类似于Matlablsqnonlin函数中的算法），右图是使用“LM”的拟合模型nls.lm。查看拟合参数（Fit.LM$par），您会发现一个拟合参数（f_parent_to_m1）位于边界处1。如果我将一个参数的起始值M0_parent从 0.1 更改为 100，那么使用nls.lm和得到相同的结果lsqnonlin。我有很多这样的情况。

newpars <- rbind(Fit.TRR$par,Fit.LM$par)
rownames(newpars)<- c("TRR(lsqnonlin)","LM(nls.lm)")
newpars

                   M0_parent   k_parent        k_m1 f_parent_to_m1
TRR(lsqnonlin)  99.59848 0.09869773 0.005260654       0.514476
LM(nls.lm)      84.79150 0.06352110 0.014783294       1.000000

除了上述问题，经常会出现由返回的 Hessiannls.lm是不可逆的（特别是当某些参数在边界上时），所以我无法得到协方差矩阵的估计。另一方面，“TRR”算法（在 Matlab 中）几乎总是通过计算解点处的雅可比来给出估计。我认为这很有用，但我也确信 R 优化算法（我尝试过的那些）没有这样做是有原因的。我想知道我是否错了，通过使用 Matlab 计算协方差矩阵的方法来获得参数估计的标准误差。

最后一点，我在上一篇文章中声称，Matlablsqnonlin在几乎所有情况下都优于 R 的优化功能。我错了。从上面的示例中可以看出，如果在 R 中也实现，则 Matlab 中使用的“信任区域反射”算法实际上更慢（有时慢得多）。但是，它仍然比 R 的基本优化算法更稳定并达到更好的解决方案。

4

2 回答 2

1

首先，我不是 Matlab 和优化方面的专家，也从未使用过 R。

我不确定我明白你的实际问题是什么，但也许我可以对你的困惑有所了解：

LM 略微增强了 Gauß-Newton 方法 - 对于具有多个局部最小值的问题，它对初始状态非常敏感。包括边界通常会产生更多的最小值。

TRR 类似于 LM，但更健壮。它具有更好的“跳出”不良局部最小值的能力。与 LM 相比，它的表现更好但性能更差是完全可行的。实际上解释为什么是非常困难的。您需要详细研究算法并查看它们在这种情况下的行为。

我无法解释 Matlab 和 R 实现之间的区别，但是 TRR 有几个扩展可能 Matlab 使用而 R 没有。您交替使用 LM 和 TRR 的方法是否比单独使用 TRR 更好？

于 2013-11-11T10:26:06.800 回答
0

使用 mkin 包，您可以使用“端口”算法（据我从其文档中得知，这也是一种 TRR 算法）或使用 nls.lm 的“Marq”算法找到参数的背景。然后您可以使用“正常”起始值或“坏”起始值。

library(mkin)
packageVersion("mkin")

如果您的系统上有可用的编译器（例如，您在 Debian/Ubuntu 上安装了 r-base-dev，或在 Windows 上安装了 Rtools），最新的 mkin 版本可以大大加快这一过程，因为它们会从自动生成的 C 代码编译模型。

这定义了模型：

m <- mkinmod(parent = mkinsub("SFO", "m1"),
             m1 = mkinsub("SFO"),
             use_of_ff = "max")

您可以检查微分方程是否正确：

cat(m$diffs, sep = "\n")

然后我们适合四个变体，端口和 LM，有或没有固定为 0.1 的 M0：

f.Port = mkinfit(m, FOCUS_2006_D)
f.Port.M0 = mkinfit(m, FOCUS_2006_D, state.ini = c(parent = 0.1, m1 = 0))
f.LM = mkinfit(m, FOCUS_2006_D, method.modFit = "Marq")
f.LM.M0 = mkinfit(m, FOCUS_2006_D, state.ini = c(parent = 0.1, m1 = 0),
               method.modFit = "Marq")

然后我们看一下结果：

results <- sapply(list(Port = f.Port, Port.M0 = f.Port.M0, LM = f.LM, LM.M0 = f.LM.M0),
  function(x) round(summary(x)$bpar[, "Estimate"], 5))

哪个是

                   Port  Port.M0       LM    LM.M0
parent_0       99.59848 99.59848 99.59848 39.52278
k_parent        0.09870  0.09870  0.09870  0.00000
k_m1            0.00526  0.00526  0.00526  0.00000
f_parent_to_m1  0.51448  0.51448  0.51448  1.00000

所以我们可以看到，即使起始值不好，Port 算法也能找到最佳解决方案（据我所知）。使用自动生成 C 代码可以缓解使用更复杂模型时可能遇到的速度问题。

于 2016-10-19T15:10:48.913 回答