用 R 中的数值迭代器改进收敛算法
Improving Convergence Algorithms with Numerical Iterator in R
我正在执行迭代计算以检查 y
在 R 中如何随 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 截距的方法?
- 谁能给我指出一些解决此类问题的文献(最好是为像我这样的初学者写的)?
- 欢迎就代表 problem/algorithm 的 class 的术语或关键字提出任何建议。
- 我提出的解决方案可以改进吗?
- 欢迎就如何使更广泛的社区或具有潜在解决方案的专家更容易访问 title/question 提出任何建议。
根据@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)
我正在执行迭代计算以检查 y
在 R 中如何随 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 截距的方法?
- 谁能给我指出一些解决此类问题的文献(最好是为像我这样的初学者写的)?
- 欢迎就代表 problem/algorithm 的 class 的术语或关键字提出任何建议。
- 我提出的解决方案可以改进吗?
- 欢迎就如何使更广泛的社区或具有潜在解决方案的专家更容易访问 title/question 提出任何建议。
根据@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)