如何在 r 中并行计算重复循环？

Question

我正在尝试找到方程组的根。这是我正在使用的 R 代码：

x1 <- 0
x2 <- 0
counter <- 0
eps <- 0.000001
repeat {
       x1_old<-x1
       x2_old<-x2
       counter <- counter + 1
       res<-uniroot(fun_x1,c(0,5),tol = 0.000001)
       x1<-res$root

       res_o<-uniroot(fun_x2,c(0,5),tol = 0.000001)
       x2 <- res_o$root

       print(c(counter,x1,x2,x1_old,x2_old))
       if (abs(sum(c(x1,x2)-c(x1_old,x2_old))) < eps)
          break
     }

这里的fun_x1 和fun_x2 是涉及x1 和x2 的两个等式。此代码需要一段时间才能找到根。我想知道有没有办法在 R 中并行计算这个 repeat 函数？

函数fun_x1 和fun_x2 是嵌套积分。例如，fun_x1 的简化版本是，

fun_x1<-function(x1)
{
  s<-7

  f123_uv<-function(u)
  {
    f123_inner<-function(v)
    {
      prob_23_v<-(exp(-(integrate(fun1,0,v-u)$value*x1+integrate(fun2,0,v-u)$value*x2)))*fun1(v-u)*x1
    }         
  }

  p_123<-integrate(Vectorize(f123_uv),0,s)$value
  return(p_123)
}

Answer 1

由于提供的示例函数不完整（fun1 未定义），我使用了一对微不足道的函数，但使用 sleep 调用来模拟一些繁重的计算：

s <- 0.1
fun_x1 <- function(x1) {
  Sys.sleep(s)
  2 + 0.5 * x2 -x1
}
fun_x2 <- function(x2) {
  Sys.sleep(s)
  3 + 0.25 * x1 -x2
}

作为基线，我们将调用您的代码：

eps <- 0.000001

t1 <- Sys.time()
x1 <- 0
x2 <- 0
counter <- 0
repeat {
  x1_old<-x1
  x2_old<-x2
  counter <- counter + 1
  res<-uniroot(fun_x1,c(0,5),tol = 0.000001)
  x1<-res$root

  res_o<-uniroot(fun_x2,c(0,5),tol = 0.000001)
  x2 <- res_o$root

  if (abs(sum(c(x1,x2)-c(x1_old,x2_old))) < eps) {
    print(c(counter,x1,x2,x1_old,x2_old))
    break
  }
}
#> [1] 10  4  4  4  4
t2 <- Sys.time()
print(t2 -t1)
#> Time difference of 8.089114 secs

这里需要在 8 秒内进行 10 次迭代才能找到公共根。但是，这不能并行化，因为每个步骤都取决于上一步的结果。我们可以通过首先找到两个单独的根然后更新x1 和x2 来解决这个问题。问题是这种方式收敛速度较慢：

t1 <- Sys.time()
x1 <- 0
x2 <- 0
counter <- 0
repeat {
  x1_old<-x1
  x2_old<-x2
  counter <- counter + 1
  res<-uniroot(fun_x1,c(0,5),tol = 0.000001)
  res_o<-uniroot(fun_x2,c(0,5),tol = 0.000001)

  x1<-res$root
  x2 <- res_o$root

  if (abs(sum(c(x1,x2)-c(x1_old,x2_old))) < eps) {
    print(c(counter,x1,x2,x1_old,x2_old))
    break
  }
}
#> [1] 16.000000  4.000000  4.000000  3.999999  4.000000
t2 <- Sys.time()
print(t2 -t1)
#> Time difference of 12.91926 secs

对于我的示例函数，它现在需要在 13 秒内进行 16 次迭代。但是，这种形式可以并行化，因为我们可以使用 future 包并行计算两个根：

library(future)
plan("multiprocess")

t1 <- Sys.time()
x1 <- 0
x2 <- 0
counter <- 0
repeat {
  x1_old<-x1
  x2_old<-x2
  counter <- counter + 1
  res %<-% uniroot(fun_x1,c(0,5),tol = 0.000001)
  res_o <- uniroot(fun_x2,c(0,5),tol = 0.000001)

  x1 <- res$root
  x2 <- res_o$root

  if (abs(sum(c(x1,x2)-c(x1_old,x2_old))) < eps) {
    print(c(counter,x1,x2,x1_old,x2_old))
    break
  }
}
#> [1] 16.000000  4.000000  4.000000  3.999999  4.000000
t2 <- Sys.time()
print(t2 -t1)
#> Time difference of 7.139439 secs

仍然需要 16 次迭代，但这次是在 7 秒内完成。这几乎比以前的版本快两倍，即几乎没有开销。然而，原始版本几乎一样快，因为它收敛得更快。 如果并行执行的加速值得较慢的收敛，您将不得不尝试使用您的实际功能。

顺便说一句，你有没有检查过没有更好的算法来找到这个公共根？