Vectorizing nested ifelse

Question

I'm trying to fasten my function in R. It contains of three ifelse statements where one of it is nested. For the single one I conducted vectorization which reduced my computation time. Unfortunately I don't see how I can vectorize the nested one. Every way I apply it returns an error. Furthemore if there is any another quirk I can use to speed it up?

cont.run <- function(reps=10000, n=10000, d=0.005, l=10 ,s=0.1) {
  r <- rep(0, reps)
  theta <- rep(0, n)
  for (t in 1:reps) {
    epsilon <- rnorm(1, 0, d)
    Zt = sum(ifelse(epsilon > theta, 1, 
                ifelse(epsilon < -theta, -1, 0)))
    r[t] <- Zt / (l * n)
    theta <- ifelse(runif(n) < s, abs(r[t]), theta)
  }
  return(mean(r))
}

system.time(cont.run())

I got:

cont.run <- function(reps=10000, n=10000, d=0.005, l=10 ,s=0.1) {
  r <- rep(0, reps)
  theta <- rep(0, n)
  for (t in 1:reps) {
    epsilon <- rnorm(1, 0, d)
    Zt = rep(NA, length(theta))
    Zt = sum(Zt[epsilon > theta, 1])
    Zt = sum(Zt[epsilon < -theta, -1])
    r[t] <- Zt / (l * n)
    theta = rep(theta, length(s))
    theta[runif(n)  < s] = abs(r[t])  
  }
  return(mean(r))
}

system.time(cont.run())

user_1286528 · Answer

Here's a little bit improved code.
Main change is that we don't use double ifelse, but instead perform two sums on TRUE vectors (sum(epsilon > theta) - sum(epsilon < -theta)) (we don't care about zeroes here). I added a couple of other improvements (eg., replaced rep with numeric, moved some operations outside the for loop).

contRun <- function(reps = 1e4, n = 1e4, d = 5e-3, l = 10, s = 0.1) {
    # Replace rep with numeric
    r <- numeric(reps)
    theta <- numeric(n)    
    # Define before loop
    ln <- l * n
    # Don't use t as it's a function in base R
    for (i in 1:reps) {
        epsilon <- rnorm(1, 0, d)
        # Sum two TRUE vectors
        r[i] <- (sum(epsilon > theta) - sum(epsilon < -theta)) / ln
        # Define before ifelse
        absr <- abs(r[i])
        theta <- ifelse(runif(n) < s, absr, theta)
    }
    return(mean(r))
}

library(microbenchmark)
microbenchmark(cont.run(), contRun())

Unit: seconds                       
       expr       min        lq      mean    median        uq       max neval
 cont.run() 13.652324 13.749841 13.769848 13.766342 13.791573 13.853786   100
  contRun()  6.533654  6.559969  6.581068  6.577265  6.596459  6.770318   100

PS. For this kind of computing you might one to set seed (set.seed() before the for loop) to make sure that you can reproduce your results.

user_5861244 · Answer

Furthemore if there is any another quirk I can use to speed it up?

In addition to PoGibas answer, you can avoid calling ifelse and get a faster function as follows

contRun <- function(reps = 1e4, n = 1e4, d = 5e-3, l = 10, s = 0.1) {
  # Replace rep with numeric
  r <- numeric(reps)
  theta <- numeric(n)    
  # Define before loop
  ln <- l * n
  # Don't use t as it's a function in base R
  for (i in 1:reps) {
    epsilon <- rnorm(1, 0, d)
    # Sum two TRUE vectors
    r[i] <- (sum(epsilon > theta) - sum(epsilon < -theta)) / ln
    # Define before ifelse
    absr <- abs(r[i])
    theta <- ifelse(runif(n) < s, absr, theta)
  }
  mean(r)
}

contRun2 <- function(reps = 1e4, n = 1e4, d = 5e-3, l = 10, s = 0.1) {
  r <- numeric(reps)
  theta <- numeric(n)    
  ln <- l * n
  for (i in 1:reps) {
    epsilon <- rnorm(1, 0, d)
    r[i] <- (sum(epsilon > theta) - sum(epsilon < -theta)) / ln
    absr <- abs(r[i])
    # avoid ifelse
    theta[runif(n) < s] <- absr
  }
  mean(r)
}

contRun3 <- function(reps = 1e4, n = 1e4, d = 5e-3, l = 10, s = 0.1) {
  r <- numeric(reps)
  theta <- numeric(n)    
  ln <- l * n
  for (i in 1:reps) {
    epsilon <- rnorm(1, 0, d)
    r[i] <- (sum(epsilon > theta) - sum(epsilon < -theta)) / ln
    absr <- abs(r[i])
    # replace runif
    theta[sample(c(T, F), prob = c(s, 1 - s), size = n, replace = TRUE)] <- absr
  }
  mean(r)
}

# gives the same
set.seed(1)
o1 <- contRun()
set.seed(1)
o2 <- contRun2()
set.seed(1)
o3 <- contRun3()

all.equal(o1, o2)
#R [1] TRUE
all.equal(o1, o3) # likely will not match
#R [1] [1] "Mean relative difference: 0.1508537"

# but distribution is the same
set.seed(1)
c1 <- replicate(10000, contRun2(reps = 100, n = 100))
c2 <- replicate(10000, contRun3(reps = 100, n = 100))
par(mfcol = c(1, 2), mar = c(5, 4, 2, .5))
hist(c1, breaks = seq(-.015, .015, length.out = 26))
hist(c2, breaks = seq(-.015, .015, length.out = 26))

# the latter is faster
microbenchmark::microbenchmark(
  contRun  = {set.seed(1); contRun ()}, 
  contRun2 = {set.seed(1); contRun2()},
  contRun3 = {set.seed(1); contRun3()},
  times = 5)
#R Unit: seconds
#R      expr      min       lq     mean   median       uq      max neval
#R   contRun 7.121264 7.371242 7.388159 7.384997 7.443940 7.619352     5
#R  contRun2 3.811267 3.887971 3.892523 3.892158 3.921148 3.950070     5
#R  contRun3 1.920594 1.920754 1.998829 1.999755 2.009035 2.144005     5

The only bottleneck now is runif in contRun2. Replacing it with sample yields quite an improvement.

Vectorizing nested ifelse

2 Answers