在 R 中使用数值迭代器改进收敛算法答案

【问题标题】：Improving Convergence Algorithms with Numerical Iterator in R在 R 中使用数值迭代器改进收敛算法
【发布时间】：2019-06-26 23:47:25
【问题描述】：

我正在执行迭代计算以检查 R 中 y 与 x 的变化情况。我的目标是估计 x 截距。现在每次迭代的计算量都很大，因此实现这一点所需的迭代越少越好。

这是y 与x 对比的图像我通过定义一个充分捕捉问题的渐近函数创建了一个工作示例

y <- (-1/x)+0.05

绘制时产生的结果

x <- 1:100
y <- (-1/x)+0.05
DT <- data.frame(cbind(x=x,y=y))
ggplot(DT, aes(x, y)) + geom_point() + geom_hline(yintercept=0, color="red")

我开发了两种迭代算法来近似 x 截距。

解决方案 1 ：x 最初很小，是步进 1...n 倍。步长的大小是预先定义的，从大开始（增加 10 倍）。在每一步之后计算y.i。如果 abs(y.i) < y[i-1] 则重复该大步，除非 y.i 更改符号，这表明该步超出了 x 截距。如果算法超调，则我们回溯并采取更小的步骤（增加 2 倍）。随着每次超调，从 10*,2*,1.1*,1.05*,1.01*,1.005*,1.001* 开始，步长越来越小。

x.i <- x <- runif(1,0.0001,0.001)
y.i <- y <- (-1/x.i)+0.05
i <- 2
while(abs(y.i)>0.0001){
  x.i <- x[i-1]*10
  y.i <- (-1/x.i)+0.05
  if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
    x <- c(x,x.i); y <- c(y,y.i)
  } else {
    x.i <- x[i-1]*2
    y.i <- (-1/x.i)+0.05
    if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
      x <- c(x,x.i); y <- c(y,y.i)
    } else {
      x.i <- x[i-1]*1.1
      y.i <- (-1/x.i)+0.05
      if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
        x <- c(x,x.i); y <- c(y,y.i)
      } else {
        x.i <- x[i-1]*1.05
        y.i <- (-1/x.i)+0.05
        if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
          x <- c(x,x.i); y <- c(y,y.i)
        } else {
          x.i <- x[i-1]*1.01
          y.i <- (-1/x.i)+0.05
          if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
            x <- c(x,x.i); y <- c(y,y.i)
          } else {
            x.i <- x[i-1]*1.005
            y.i <- (-1/x.i)+0.05
            if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
              x <- c(x,x.i); y <- c(y,y.i)
            } else {
              x.i <- x[i-1]*1.001
              y.i <- (-1/x.i)+0.05
            }
          }
        }
      }
    }
  }
  i <- i+1
}

解决方案 2：该算法基于 Newton-Raphson 方法的思想，其中步骤基于 y 的变化率。变化越大，采取的步骤就越小。

x.i <- x <- runif(1,0.0001,0.001)
y.i <- y <- (-1/x.i)+0.05
i <- 2
d.i <- d <- NULL
while(abs(y.i)>0.0001){
  if(is.null(d.i)){
    x.i <- x[i-1]*10
    y.i <- (-1/x.i)+0.05
    d.i <- (y.i-y[i-1])/(x.i-x[i-1])
    x <- c(x,x.i); y <- c(y,y.i); d <- c(d,d.i)
  } else {
    x.i <- x.i-(y.i/d.i)
    y.i <- (-1/x.i)+0.05
    d.i <- (y.i-y[i-1])/(x.i-x[i-1])
    x <- c(x,x.i); y <- c(y,y.i); d <- c(d,d.i)
  }
  i <- i+1
}

比较

解决方案 1 所需的迭代次数始终少于解决方案 2（如果不是 1/3，则为 1/2）。
解决方案 2 更优雅，不需要任意减小步长。
我可以设想解决方案 1 卡住的几种情况（例如，即使在最小的步骤中，循环也不会收敛到足够小的 y.i 值）

问题

在这种情况下，是否有更好（迭代次数更少）逼近 x 截距的方法？
谁能给我指出一些解决此类问题的文献（最好是写给像我这样的初学者可以理解的）？
欢迎对代表此类问题/算法的命名法或关键词提出任何建议。
我提出的解决方案可以改进吗？
欢迎任何有关如何使更广泛的社区或具有潜在解决方案的专家更容易访问标题/问题的建议。

【问题讨论】：

为什么要重新发明轮子？你不能只为你的数据拟合一个模型并从中推断出属性吗？
您可以使用approxfun 或拟合模型，如@RomanLuštrik 建议的那样，将y 作为x 的函数，然后使用uniroot 找到根。
@RomanLuštrik 每个点在计算上都是昂贵的，所以我宁愿不计算一大堆，然后为它拟合一个模型。此外，即使我确实有一个模型，我可能仍然需要一种迭代方法来计算 x 截距，因为可能无法求解 y=0。（例如，rms 包中的 rcs 有很多 pmax() 函数来指定节点位置，这使得解决 x 变得困难
@JustGettinStarted 如果我理解正确，您只是在数字上寻找根，对吗？如果是这样，this 将是一个很好的资源。
正是在这种情况下，昂贵的函数评估，导数不可用，Dekker、Müller、Brent 等是否开发了他们改进的 regula-falsi 和二等分方法的组合。 Brent 的实现应该在其中一个标准库中可用。

标签： r algorithm iterator numerical-methods convergence

【解决方案1】：

根据@Lyngbakr 和@LutzL 的阅读建议和建议，被称为Brent's Method 的寻根算法被证明是有效的，并且比我实现的Newton-Raphson（解决方案2）要快得多。该算法由uniroot在R中实现。

f <- function(x){(-1/x)+0.05}
uniroot(f, interval=c(0,100), maxiter=100)

【讨论】：