R 中的优化:Levenberg-Marquardt 在 minpack.lm 中使用 nls.lm:将“maxiter”重置为 1024

Optimization in R: Levenberg-Marquardt using nls.lm in minpack.lm: resetting `maxiter' to 1024

我正在尝试学习如何使用 R 库 minpack.lm 中的 nls.lm,方法是使用 Rosenbrock 函数查看算法是否收敛于 f(x,y) 处的全局最小值= (1,1)。无论是否使用解析雅可比行列式,我都会这样做。在这两种情况下,我都会收到一条警告,告诉我算法已决定将调用 nls.lm 时指定的最大迭代次数还原为 1024:

Warning messages:
1: In nls.lm(par = initpar, fn = objective_rosenbrock, jac = gradient_rosenbrock,  :
  resetting `maxiter' to 1024!
2: In nls.lm(par = initpar, fn = objective_rosenbrock, jac = gradient_rosenbrock,  :
  lmder: info = -1. Number of iterations has reached `maxiter' == 1024.

鉴于我最初猜测的 (-1.2, 1.0),该算法从未完全达到 (1,1)。我在 GitHub 上找到了该库的源代码,以下代码行与此处相关:

https://github.com/cran/minpack.lm/blob/master/src/nls_lm.c

OS->maxiter = INTEGER_VALUE(getListElement(control, "maxiter"));
if(OS->maxiter > 1024) {
  OS->maxiter = 1024;
  warning("resetting `maxiter' to 1024!");
}

为什么将最大迭代次数限制为 1024 有什么逻辑吗?有位和 2^10 的东西?我想将该库用于不同的应用程序,但这种迭代上限可能会阻止这种情况。任何见解将不胜感激。

Git blame says that this code limiting the max iterations was introduced in version 1.1-0, in 2008. The NEWS file for the package 只能追溯到版本 1.1-6。除了你指向的那个(它只是一个 CRAN 镜像;它不包含任何 comments/commit messages/etc。来自开发人员的可能会给我们一些线索。)

除了联系维护者之外,我认为很难弄清楚这个限制的基本原理是什么。

虽然我确实有一些猜测。

代码中实际使用 maxiter 的地方只有 here and here - in R code, not Fortran or C code, so it seems extremely unlikely that we are dealing with something like a 10-bit unsigned integer type (which seems an unlikely choice in any case). I think the limitation is there because we also have a buffer defined for holding trace information here:

  double rsstrace[1024];

如您所见,它被硬编码为 1024 的长度。如果我们试图将 1025 次迭代的跟踪信息填充到此数组中,可能会发生不好的事情...

我的建议:

  • 将代码中“1024”的所有实例更改为更大的值,看看会发生什么。只有四个:
$ find . -type f -exec grep -Hn 1024 {} \;
./src/nls_lm.c:141:    if(OS->maxiter > 1024) {
./src/nls_lm.c:142:      OS->maxiter = 1024;
./src/nls_lm.c:143:      warning("resetting `maxiter' to 1024!");
./src/minpack_lm.h:20:  double rsstrace[1024];
  • 最好在 src/minpack_lm.h#define MAXITER 2048(或其他)并使用它来代替数值。
  • 联系维护者 (maintainer("minpack.lm")) 并询问他们这个问题。