1

I am trying to find the roots of a system of equations. Here is the R code I am using:

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
     }

Here fun_x1 and fun_x2 are the two equations both involving x1 and x2. This code takes a while to find the roots. I would like to know is there any way to parallel compute this repeat function in R?

The functions fun_x1 and fun_x2 are nested integrals. For example, a simplified version of fun_x1 is,

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)
}
ssaha
  • 459
  • 2
  • 10
  • It seems each iteration of your loop is dependent on the values from the previous iteration. How would you make such a loop parallel? Doing parallel operations makes more sense when the computations aren't dependent on each other. – MrFlick Jul 31 '18 at 16:13
  • Yes, my current iteration values depend on the previous iteration. That's why I am having difficulties to parallel compute the process. @MrFlick – ssaha Jul 31 '18 at 16:51
  • 2
    How many steps does it take? If there are many steps then remove the print statement, since IO in a tight loop can be expensive. If there are only a few steps needed, then it might be possible to parallelise the two uniroot calls, e.g. using the future package. – Ralf Stubner Aug 01 '18 at 08:38
  • @RalfStubner: It only takes 6-7 steps to converge. I have never used the future package. I don't know how I can parallel compute a dependent loop. Would you mind providing me with a detailed answer of how to program such loop? Thanks! – ssaha Aug 03 '18 at 17:44
  • Can you provide (maybe simplified) definitions for `fun_x1` and `fun_x2`? – Ralf Stubner Aug 03 '18 at 17:52
  • Parallel processing involves running the different iterations of the loop simultaneously. If each iteration relies on values generated by the previous iteration it isn't possible to parallelize your loop, at least without totally rethinking your algorithm. To do parallel processing, you need to be able to divide your program into chunks that can be run independently and then put back together at the end. The second answer to this question might help clarify this: https://stackoverflow.com/questions/1050222/what-is-the-difference-between-concurrency-and-parallelism?rq=1 – divibisan Aug 07 '18 at 00:38
  • @RalfStubner: The functions are nested integrals. I updated my question and provided a simplified version of `fun_x1`. Please see the question. – ssaha Aug 07 '18 at 00:58

1 Answers1

0

Since the provided sample function is not complete (fun1 is not defined), I am using a pair of trivial functions but with a sleep call to simulate some heavy computation:

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
}

As a baseline, we than call your code:

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

Here it takes 10 iterations in 8s to find the common root. However, this cannot be parallelized since each step depends on the result of the previous step. We can untangle this at one point by first finding both individual roots and then updating x1 and x2. Problem is that convergence is slower this way:

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

For my sample functions it now takes 16 iterations in almost 13s. However, this form can be parallelized since we can compute both roots in parallel using the future package:

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

It still takes 16 iterations, but this time they are done in 7s. This is almost a factor of two faster than the previous version, i.e. there is little overhead. However, the original version is nearly as fast, since it converges faster. You will have to try with your real functions if the speed-up from parallel execution is worth the slower convergence.

BTW, have you checked that there are no better algorithms for finding this common root?

Ralf Stubner
  • 26,263
  • 3
  • 40
  • 75