在 R 中使用“nls.lm”时如何处理边界约束

Question

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

根据我的经验，这样做程序会达到最佳状态，并且对起始值不太敏感。但这只是一种实用的方法，可以避开我在使用 nls.lm 以及 R 中的其他优化函数时遇到的问题。我想知道为什么 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)

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

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

左图是带有初值的模型，中间是使用“TRR”拟合的模型（类似于Matlab中的算法lsqnonlin函数），右图是使用“ 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

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

最后一点，我在previous post 中声称，Matlab lsqnonlin 在几乎所有情况下都优于 R 的优化函数。我错了。从上面的例子可以看出，如果在 R 中也实现了，Matlab 中使用的“Trust-Region-Reflective”算法实际上更慢（有时慢得多）。但是，它仍然比 R 的基本优化算法更稳定并达到更好的解决方案。

Answer 1

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

我不确定我明白你的实际问题是什么，但也许我可以解释一下你的困惑：

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

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

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

Answer 2

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

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 代码可以缓解使用更复杂模型时可能遇到的速度问题。