超越 R 的 optim 函数答案

【问题标题】：Moving beyond R's optim function超越 R 的 optim 函数
【发布时间】：2011-04-15 00:19:05
【问题描述】：

我正在尝试使用 R 来估计具有手动规范的多项 logit 模型。我发现了一些软件包可以让您估计 MNL 模型 here 或 here。

我发现了一些关于“滚动”你自己的 MLE 函数here 的其他著作。然而，根据我的研究，所有这些函数和包都依赖于内部的optim 函数。

在我的基准测试中，optim 是瓶颈。使用具有约 16000 个观察值和 7 个参数的模拟数据集，R 在我的机器上大约需要 90 秒。 Biogeme 中的等效模型大约需要 10 秒。一位在Ox 中编写自己的代码的同事报告说，同一模型大约需要 4 秒。

有没有人有编写自己的 MLE 函数的经验，或者可以向我指出一些超出默认 optim 函数（没有双关语）的优化方向？

如果有人想要 R 代码重新创建模型，请告诉我 - 我很乐意提供。我没有提供它，因为它与优化optim 函数和保留空间的问题没有直接关系......

编辑：感谢大家的想法。基于下面无数的 cmets，对于更复杂的模型，我们能够在与 Biogeme 相同的范围内获得 R，而对于我们运行的几个更小/更简单的模型，R 实际上更快。我认为这个问题的长期解决方案将涉及编写一个独立的最大化函数，该函数依赖于 fortran 或 C 库，但我当然愿意接受其他方法。

【问题讨论】：

魔鬼在细节中。您可能会弄乱optim 参数（请参阅文档中有关control 的部分）。您可以将默认参数与您的同事代码或 Biogeme 使用的默认参数进行比较。它们是否不同，如果是，那为什么？
@Marek - Biogeme 依赖于用 C 编写的自定义最大化例程，这与 Ox 的故事类似。这对我来说是一个新领域，但我开始了解使用的不同方法。
据我了解，nlm() 和 R 中可能的其他优化例程已经用 C 语言编写。我宁愿建议你直接寻找对内部函数的访问，这样你就可以摆脱开销而不是重新发明轮子

标签： optimization r

【解决方案1】：

我是 R 包 optimParallel 的作者，这可能对您的情况有所帮助。该软件包提供了optim() 的基于梯度的优化方法的并行版本。该包的主要功能是optimParallel()，其用法和输出与optim()相同。使用optimParallel() 可以显着减少优化时间，如下图所示（p 是参数个数）。请参阅https://cran.r-project.org/package=optimParallel 和http://arxiv.org/abs/1804.11058 了解更多信息。

【讨论】：

【解决方案2】：

FWIW，我在 C-ish 中使用 OPTIF9 完成了这项工作。你很难比这更快。有很多方法可以让速度变慢，比如运行像 R 这样的解释器。

补充：从 cmets 中，很明显 OPTIF9 被用作优化引擎。这意味着很可能大部分时间都花在了评估 R 中的目标函数上。虽然 C 函数可能在下面用于某些操作，但仍然存在解释器开销。有一种快速的方法可以确定 R 中的哪些代码行和函数调用大部分时间负责，那就是用 Escape 键暂停它并检查堆栈。如果一个语句花费 X% 的时间，那么它在 X% 的时间里在堆栈上。您可能会发现有些操作不属于 C 语言，但应该存在。当您找到一种方法来并行化 R 执行时，您通过这种方式获得的任何加速因素都将被保留。

【讨论】：

R 的一些底层例程是用 C 或 Fortran 编写的。所以理论上，R应该能够在速度上接近C。现在只需要找到那个好的实现即可。不知道是否使用了 OPTIF9，但对于额外的包来说这将是一个好主意 ;-)
Dennis-Schnabel 算法 'optif9' 在 R 中的 nlm 中实现。（快速浏览一下 R 源代码会发现它是对旧的 fortran 子例程的 c 重写）。
@eyjo：@Joris：很好，这也是我所期望的。然后大部分时间将花在 R 中用户编码的目标函数上。如果可以将其翻译成编译器语言，它应该与 C 处于相同的范围。
是的，你会遇到这个问题。正如 Chase 所说，Ox 或多或少是一个定制的 c 编译器，它在 4 秒内完成了他的模型，而在 R 中则需要 90 秒。但是你无法击败在 R 中使用统计数据的灵活性，并且可以轻松扩展它使用您的自定义 c 调用。如果你用包含 Rh、Rinternals.h 和 R_ext/Applic.h 头文件的 c 语言编写目标函数，并从那里调用相关的优化例程，为 R 编译它，我认为你应该非常接近并且它在 R 中无缝工作.
@eyjo：对。对于其他代码，目标函数将被调用 100 次，因为它是 optif9 的内部循环。因此，如果我无法编译它，我只会将该函数置于一个临时的无限循环中并调整它的日光。

【解决方案3】：

已经用 nlm() 函数试过了吗？不知道它是否更快，但它确实提高了速度。还要检查选项。 optim 使用慢速算法作为默认值。通过使用准牛顿算法（method="BFGS"）而不是默认算法，您可以获得 > 5 倍的加速。如果您不太关心最后一位数字，您还可以将 nlm() 的容差级别设置得更高，以获得额外的速度。

f <- function(x) sum((x-1:length(x))^2)

a <- 1:5

system.time(replicate(500,
     optim(a,f)
))
   user  system elapsed 
   0.78    0.00    0.79 

system.time(replicate(500,
     optim(a,f,method="BFGS")
))
   user  system elapsed 
   0.11    0.00    0.11 

system.time(replicate(500,
     nlm(f,a)
))
   user  system elapsed 
   0.10    0.00    0.09 

system.time(replicate(500,
      nlm(f,a,steptol=1e-4,gradtol=1e-4)
))
   user  system elapsed 
   0.03    0.00    0.03

【讨论】：

Joris，你能解释一下不知道是不是更快，但确实提高了速度？我今天只喝了两杯咖啡。您的意思是“更快，但不确定多少”？
Errrr，键盘这一侧也需要更多的咖啡...确实，在我测试过的功能上它更快，有时略微，有时相当多。它还取决于控制设置，如代码所示。
感谢关于 nlm() 的提示。如上所述改变步长和渐变后，它最终成为最快的。

【解决方案4】：

您考虑过CRAN Task View for Optimization 上的材料吗？

【讨论】：