Monday, November 28, 2011

Another aspect of speeding up loops in R

Any frequent reader of R-bloggers will have come across several posts concerning the optimization of code - in particular, the avoidance of loops.

Here's another aspect of the same issue. If you have experience programming in other languages besides R, this is probably a no-brainer, but for laymen, like myself, the following example was a total surprise. Basically, every time you redefine the size of an object in R, you are also redefining the allotted memory - and this takes some time. It's not necessarily a lot of time, but if you are having to do it during every iteration of a loop, it can really slow things down.

The following example shows three versions of a loop that creates random numbers and stores those numbers in a results object. The first example (Ex. 1) demonstrates the wrong approach, which is to concatenate the results onto the results object ("x") , thereby continually changing the size of x after each loop. The second approach (Ex. 2) is about 150x faster - x is defined as an empty matrix containing NAs, which is gradually filled (by row) during each loop. The third example (Ex. 3) shows another possibility if one does not know what the size of the results from each loop will be. An empty list is created of length equaling the number of loops. The elements of the list are then gradually filled with the loop results. Again, this is at least 150x faster than Ex. 1 (and I'm actually surprised to see that it may even be faster than Ex.2).

the example...
NROW=5000
NCOL=100
 
#Ex. 1 - Creation of a results matrix where its memory
#allocation must continually be redefined
t1 <- Sys.time()
x <- c()
for(i in seq(NROW)){
 x <- rbind(x, runif(NCOL))
}
T1 <- Sys.time() - t1
 
 
#Ex. 2 - Creation of a results matrix where its memory
#allocation is defined only once, at the beginning of the loop.
t2 <- Sys.time()
x <- matrix(NA, nrow=NROW, ncol=NCOL)
for(i in seq(NROW)){
 x[i,] <- runif(NCOL)
}
T2 <- Sys.time() - t2
 
 
#Ex. 3 - Creation of a results object as an empty list of length NROW. 
#Much faster than Ex. 1 even though the size of the list is
#not known at the start of the loop.
t3 <- Sys.time()
x <- vector(mode="list", NROW)
for(i in seq(NROW)){
 x[[i]] <- runif(NCOL)
}
T3 <- Sys.time() - t3
 
png("speeding_up_loops.png")
barplot(c(T1, T2, T3), names.arg = c("Concatenate result", "Fill empty matrix", "Fill empty list"),ylab="Time in seconds")
dev.off()
 
T1;T2;T3
 
#On my PC:
#Time difference of 54.906 secs  #T1
#Time difference of 0.2820001 secs  #T2
#Time difference of 0.1880000 secs  #T3
Created by Pretty R at inside-R.org

8 comments:

  1. Yeah, read that in the R inferno. Never grow object.

    ReplyDelete
  2. Thanks! Very helpful!!

    ReplyDelete
  3. By using NA instead of NA_real_ [or easier to remember as.double(NA)], your allocated 'x' matrix is of type logical. Then at the first iteration, x[i,] <- runif(NCOL) will coerce 'x' to a numeric matrix. That takes time and unnecessary allocations and garbage collections.

    ReplyDelete
  4. # In your second function, I think "T3 <- Sys.time() - t2" should be "T2 <- Sys.time() - t2"
    # Have you thought about using the rbenchmark package to do timings? I find it easier for this stuff.
    # Also, I know you didn't include an lapply version for comparison as you're specifically using a
    # For loop to show the difference but it's still nice to see it, so I include an example belo :)
    # Nice post!

    library(rbenchmark)

    f2 <- function(NROW, NCOL) {
    x <- matrix(NA, nrow=NROW, ncol=NCOL)
    for(i in seq(NROW)){
    x[i,] <- runif(NCOL)
    }
    return(x)
    }

    f3 <- function(NROW, NCOL) {
    x <- vector(mode="list", NROW)
    for(i in seq(NROW)){
    x[[i]] <- runif(NCOL)
    }
    return(x)
    }


    f4 <- function(NROW, NCOL) {
    x <- vector(mode="list", NROW)
    x = lapply(seq_len(NROW), function(r) runif(NCOL))
    return(x)
    }

    NROW=5000
    NCOL=100
    benchmark(f2(NROW, NCOL), f3(NROW, NCOL), f4(NROW, NCOL),
    columns = c("test", "replications", "elapsed", "relative"),
    order = "relative",
    replications = 100)

    # test replications elapsed relative
    #3 f4(NROW, NCOL) 100 1.999 1.000000
    #2 f3(NROW, NCOL) 100 2.806 1.403702
    #1 f2(NROW, NCOL) 100 3.826 1.91395

    ReplyDelete
  5. NaN will do, because it is numeric, cf. str(NaN) or storage.mode(NaN).

    ReplyDelete
  6. Thanks tonybreyal,
    Yeah, I intentionally did not use lapply as I was focusing on loops. I recomend that people should look to R-bloggers for more great examples of optimization. Thanks for catching my typo - had it fixed in my version at home, but have updated the example above now too.
    Cheers, Marc

    ReplyDelete
  7. You should plot time on a log scale in your graph. A time difference of 0.2820001 secs for T2 is quite different from Time difference of 0.1880000 secs for T3.

    For your example, it doesnt really matter. But if you had say, NROW= 5000000 and NCOL= 10000, it would be a pretty substantial difference.

    Just saying,

    sebastien

    ReplyDelete