High Performance Computing in R using Rcpp

Abstract

R is inherently a high-performance scripting language since its core is written in C. Writing a high-performance R program, however, may not be as straightforward as it is for some other general purpose languages. This notenook demonstrates quickly the key concept in writing a good R program and a fact check about the decomposition of the language source code. It also demonstrates in more details with several working examples and benchmarks on how we can boost the performance when the underlying task is not easy to optimize in native R code: by using Rcpp.

---

The notebook is a full re-writing of its original version for a lightning talk at Data Science Conference 2015 at Taipei. The coding examples in the original notebook are for Linux/macOS users only. The revised notebook here is platform-independent.

---

1 Motivation

The performance of a R program can vary a lot, depending on whether it is written by an experienced R programmer or a newbie who don’t know how to properly write R. Unfortunately, unlike other general purpose langauges, it is much more likely to write improper R codes which run embarrassingly slow. Two important and common factors that affect the performance of R code:

1. Vectorization
2. Memory Copy

In general, we should use vectorization as much as possible and avoid memory copy as much as possible. For example the following code is TERRIBLE:

numbers <- 1:1e5

pow_all_1 <- function(numbers) {
  out <- c()
  for ( n in numbers ) {
    out <- c(out, n^2)  # Memory copy is happening here.
  }
  out
}

system.time(out_1 <- pow_all_1(numbers))

   user  system elapsed
  8.593   0.328   8.923 

We are just trying to calculate the base-2 power of a sequence of numbers. It will be very straghtforward to just loop over the numbers and compute the result for a general purpose language.

But NOT for R.

The above code is just WRONG because it violates both our two principles: it didn’t vectorize when it could, and it didn’t avoid memory copy. Let’s first fix the second issue about the memory copy:

pow_all_2 <- function(numbers) {
  out <- numeric(length(numbers))  # Pre-allocate memory.
  for ( i in 1:length(numbers) ) {
    out[i] <- numbers[i]^2
  }
  out
}

system.time(out_2 <- pow_all_2(numbers))

   user  system elapsed
  0.015   0.000   0.010 

stopifnot(all.equal(out_1, out_2))

The output vector size is pre-determined for our second attempt, which avoid copies happening again and again during the unnecessary concatenation in the loop.

Still the code is terrible, because it didn’t vectorize when it actually could:

pow_all_3 <- function(numbers) {
  numbers^2
}

system.time(out_3 <- pow_all_3(numbers))

   user  system elapsed
      0       0       0 

stopifnot(all.equal(out_1, out_3))

The difference may not be obvious at the first glance because the task is too trivial. Now for a serious benchmarking at microsecond-level:

library(microbenchmark)

numbers <- 1:1e4  # Make it smaller for faster benchmarking.

ben <- microbenchmark(
  pow_all_1(numbers),
  pow_all_2(numbers),
  pow_all_3(numbers),
  times=10
)

print(ben)

Unit: microseconds
               expr     min      lq     mean   median      uq      max neval
 pow_all_1(numbers) 85837.5 86857.6 92445.84 87417.35 88266.0 138826.5    10
 pow_all_2(numbers)   624.8   627.5   632.62   628.45   637.7    645.7    10
 pow_all_3(numbers)    15.9    17.0   100.34    17.35    24.0    833.2    10

A properly written R program can be of several orders of magnitude faster.

Even if we already carefully write our R code, it may still occur that we want even faster speed. This is especially true when there is no trivial way of vectorization for the implementation desired. This is the main topic of this notebook: we can still be faster, by using R’s C interface.

Rcpp is the most dominant way of achieving that. But before our deep dive, let’s understand what is R a bit more.

2 Decompose the R Source Code

Let’s examine the source code of the R language.

# Download the source code of R lang.
# This may take a while.

TEMPDIR <- tempdir()
target_file <- "https://cran.r-project.org/src/base/R-3/R-3.6.1.tar.gz"
downloaded_file <- file.path(TEMPDIR, basename(target_file))

if ( !file.exists(downloaded_file) ) {
  download.file(target_file, downloaded_file)
}

untar(downloaded_file, exdir=TEMPDIR)  # This may generate warnings under Windows.

We’d like to know the file type distribution in the source code.

# List all files in the source and count the extensions.

library(data.table)

src_dir <- file.path(tools::file_path_sans_ext(downloaded_file, compression=TRUE), "src")
src_files <- list.files(src_dir, recursive=TRUE, full.names=TRUE)
src_exts <- tools::file_ext(src_files)
src_exts <- src_exts[src_exts != ""]

ext_counts <- as.data.table(sort(table(src_exts), decreasing=TRUE), keep.rownames=TRUE)
setnames(ext_counts, c("ext", "count"))
ext_counts[, ext:=factor(ext, levels=rev(ext))]

library(wordcloud2)

# Scale down to avoid over-sizing.
wordcloud2(sqrt(table(src_exts)))

library(ggplot2)

top_n <- 10
ggplot(ext_counts[1:top_n], aes(x=ext, y=count)) +
  geom_bar(stat="identity") +
  labs(title=sprintf("Top %s File Extensions in R Source Code", top_n)) +
  coord_flip()

Here are the most common file types found in the source code of R.

• .Rd: Document file
• .R: R source file
• .c: C source file
• .mo: Binary data file
• .po: Language translation file
• .h: Header file for C
• .afm: Font file
• .in: Template configuration file for some macro preprocessor
• .win: Same above but specifically for Windows
• .f: Fortran source

If we further limit to only programming language related files:

lang_ext <- c("c", "h", "f", "R")
lang_ext_counts <- ext_counts[ext %in% lang_ext]
lang_ext_counts[, pct:=count / sum(count)]

lang_ext_counts[, ext:=factor(ext, levels=lang_ext)]

ggplot(lang_ext_counts, aes(x=ext, y=count)) +
  geom_bar(stat="identity") +
  labs(y="Number of Files", x="File Extension",
       title="Language File Distribution in R Source Code") +
  geom_text(aes(label=scales::percent(pct), y=count),
            vjust=1, size=5, color="white")

The fact is that, R is heavily written in C. Primitive functions are writtin in C. Most of them are vectorized so it is very fast to apply the function directly to a vector. For example the power function we just examined previously is a primitive function:

`^`  # A primitive.

function (e1, e2)  .Primitive("^")

(1:10)^2  # A vectorized primitive function call.

 [1]   1   4   9  16  25  36  49  64  81 100

Instead of counting files, we can further count the number of lines in those language files to arrive at the distribution at file line level:

lang_src_files <- src_files[tools::file_ext(src_files) %in% lang_ext]
lang_src_file_lens <- sapply(lang_src_files, R.utils::countLines)

lang_len <- data.table(lang=gsub("^.*\\.", "", names(lang_src_file_lens)), len=lang_src_file_lens)
lang_len_counts <- lang_len[, .(tot_lines=sum(len)), by="lang"]
lang_len_counts[, pct:=tot_lines / sum(tot_lines)]

lang_len_counts[, ext:=factor(lang, levels=lang_ext)]

ggplot(lang_len_counts, aes(x=lang, y=tot_lines)) +
  geom_bar(stat="identity") +
  labs(y="Number of Lines", x="File Extension",
       title="Language File Line Distribution in R Source Code") +
  geom_text(aes(label=scales::percent(pct), y=tot_lines),
            vjust=1.25, size=5, color="white")

The take-away is that, R simply cannot be slow for most of the fundamentally demanding computating tasks. Its core is written in C and C is super fast.

3 Working Examples of Using Rcpp

There will be cases where the native R code is not suitable for the desired implementation. This is usually because the built-in primitive functions are not available for a particular kind of algorithm. In such scenario we can use Rcpp to re-write only the most demanding part of our implementation in C++ and port it to our R program seamlessly. In this section we are going to illustrate several such use cases.

3.1 N-Gram Generation

N-gram generation is a common task for natural language processing and understanding. Despite its simplicity, surprisingly, it is not easy for native R code to fulfill the task in an efficient manner.

Let’s assume we’d like to count the bigrams of the hex dump of a (potentially binary) file. The byte code information could be useful for a downstream machine learning task such as file category classification. But that is not our concern for now. We just want to implement a function that can extract ngrams given a long string, and count their frequencies.

# Prepare hex dump of a file.
# We will use the Rmd file itself as an example.
infile <- "rcpp.Rmd"
con <- file(infile, "rb")
hex <- readBin(con, what="raw", n=file.info(infile)$size)
hex_c <- as.character(hex)
close(con)

str(hex_c)  # A character vector of hexadecimals.

 chr [1:31050] "2d" "2d" "2d" "0a" "74" "69" "74" "6c" "65" "3a" "20" "22" "48" "69" "67" "68" "20" "50" "65" "72" ...

Native R

Now we implement the minimum ngram counter given a character vector, assuming each element is a token:

ngram_r <- function(x, n) {
    len <- length(x) - (n - 1)
    out <- character(len)  # Pre-allocate size.
    for ( i in 1:len ) {
        out[i] <- x[i]
        if ( n > 1 )
            for ( j in 1:(n-1) )
                out[i] <- paste0(out[i], x[i+j])
    }
    sort(table(out), decreasing=TRUE)
}

out_r <- ngram_r(hex_c, 2)

print(length(out_r))  # Total ngrams.

[1] 1460

print(out_r[1]) # The most frequent ngram.

2020
 851 

The above function, though implementation is straightforward, is terribly slow because it didn’t utilize vectorization. And there doesn’t seem to have a readily available primitive function for this very purpose. This is the usually frustrated experience when people familiar with other general purpose languages first come to use R.

Rcpp

For a function requiring arbitraily large explicit loop and there is no obvious way to vectorization, we have Rcpp come to the rescue.¹

For exactly the same minimum implementation we can re-write it in C++ using the Rcpp API:

library(Rcpp)
print(installed.packages()["Rcpp", "Version"])

[1] "1.0.2"

cppFunction("
  CharacterVector ngramRcpp (CharacterVector hexvector, int ngram) {
    int len = hexvector.size() - (ngram - 1);
    CharacterVector out(len);
    for (int i = 0; i < len; i++) {
      out(i) = hexvector[i];
      for (int j = 1; j < ngram; j++) {
        out(i) += hexvector[i+j];
      }
    }
    return out;
}")

ngramRcpp  # A primitive function.

function (hexvector, ngram)
.Call(<pointer: 0x7fe6ac4864a0>, hexvector, ngram)

ngram_rcpp <- function(x, n) {
  ngrams <- ngramRcpp(x, 2)
  sort(table(ngrams), decreasing=TRUE)
}

out_rcpp <- ngram_rcpp(hex_c, 2)

all.equal(out_rcpp, out_r, check.attributes=FALSE)

[1] TRUE

Notice that we only implement the ngram extraction loop but leave the frequency count for native R. This is because the table function is not a bottleneck of our task.

In this example we use the CharacterVector class implemented in Rcpp which mimics the character vector we usually use in native R code. There are many such high-level classes we can find in the Rcpp library. These classes will save our development time when writing the C++ code since we won’t limit ourselves to only C++ standard library.

print(microbenchmark(
  ngram_r(hex_c, 2),
  ngram_rcpp(hex_c, 2),
  times=10
))

Unit: milliseconds
                 expr     min      lq     mean   median      uq     max neval
    ngram_r(hex_c, 2) 59.4615 59.7028 60.16685 59.94985 60.2995 61.8870    10
 ngram_rcpp(hex_c, 2)  5.0451  5.1562  6.02980  5.34610  7.3100  8.7338    10

The Rcpp is much faster than the native R implementation. And it will be even faster if the data is getting larger.

# Now take the entire R source code as one big file.
con <- file(downloaded_file, "rb")
big_hex <- readBin(con, what="raw", n=file.info(downloaded_file)$size)
big_hex_c <- as.character(big_hex)
close(con)

str(big_hex_c)

 chr [1:30463021] "1f" "8b" "08" "08" "99" "f6" "1e" "5d" "02" "03" "52" "2d" "33" "2e" "36" "2e" "31" "2e" "74" ...

This time we benchmark with different file size, we’ll see that theoretically Rcpp can be hundreds of times faster than native R.

# Benchmark the native R code with varying file sizes.
# Since it may take too long for the redundant tasks we'll use parallel computing.
# Another more clever way is to time on-the-fly when we parse the file and hence parse only once.
# In that case we need to modify the original function.
# Here we go for the brutal but simple way.
library(parallel)

cl <- makeCluster(detectCores() / 2)
clusterExport(cl, c("big_hex_c"))

n_size <- 10
sizes <- seq(1e3, length(big_hex_c), length.out=n_size)

time_ngram <- function(func, s, n) {
  st <- proc.time()
  func(big_hex_c[1:s], n)
  et <- proc.time()
  (et - st)["elapsed"]
}

# This could take a while.
time_r <- parLapplyLB(cl, sizes, time_ngram, func=ngram_r, n=2)

stopCluster(cl)


# Benchmark the Rcpp code with varying file sizes.
time_rcpp <- numeric(n_size)
for ( i in 1:length(sizes) ) {
  time_rcpp[i] <- time_ngram(ngram_rcpp, sizes[i], 2)
}

out_time <- data.table(r=unlist(time_r), rcpp=unlist(time_rcpp), size=sizes)
out_time <- melt(out_time, id.vars="size", variable.name="impl", value.name="time")
ggplot(out_time, aes(x=size, y=time, group=impl, color=impl)) +
  geom_line() +
  geom_point() +
  labs(x="File Size (Bytes)", y="Seconds",
       title="Speed on Bigram Generation of Varying Input Length")

Other Packages

For completeness, in this section we include the benchmark for some other well-developed packages for text processing. These packages are usually using R’s native C interface or Rcpp to implement the critical parts of the computation for their intended tasks. We will expect them to be also much more efficient than a naive implementation in native R code.

The first package we explore is quanteda:

library(quanteda)
print(installed.packages()["quanteda", "Version"])

[1] "1.5.1"

# For quanteda we need to convert text to its `tokens` API.
# This may introduce some overhead so we will create two functions for fair comparison latter.
hex_s <- paste(hex_c, collapse=" ")
hex_t <- tokens(hex_s)

ngram_quanteda_1 <- function(x, n) {
  # Tokenization included.
  ngrams <- tokens_ngrams(tokens(x), n=n)[[1]]
  sort(table(ngrams), decreasing=TRUE)
}

ngram_quanteda_2 <- function(tok, n) {
  # Assume pre-tokenized.
  ngrams <- tokens_ngrams(tok, n=n)[[1]]
  sort(table(ngrams), decreasing=TRUE)
}

out_quanteda <- ngram_quanteda_1(hex_s, n=2)

all.equal(out_quanteda, out_r, check.attributes=FALSE)

[1] TRUE

text2vec is also a high-performance package implemented with Rcpp:

library(text2vec)
print(installed.packages()["text2vec", "Version"])

[1] "0.5.1"

# Again to use text2vec we need to follow its API:
# Converting text into `itoken` iterator.
# Similar to the case of quanteda, we use two functions where one excludes the conversion overhead.
hex_i <- itoken(hex_s, progressbar=FALSE)

ngram_text2vec_1 <- function(x, n=n) {
  it <- itoken(x, progressbar=FALSE)
  create_vocabulary(it, ngram=c(n, n))
}

ngram_text2vec_2 <- function(it, n=n) {
  create_vocabulary(it, ngram=c(n, n))
}

out_text2vec <- ngram_text2vec_1(hex_s, n=2)

# Tidy the result to align with our previous functions.
out_text2vec <- setNames(out_text2vec$term_count, out_text2vec$term)

all.equal(as.table(sort(out_text2vec, decreasing=TRUE)), out_r, check.attributes=FALSE)

[1] TRUE

Lastly, unlike the previous two packages both aim at a larger scope of natural language processing tasks, the package ngram dedicates only at ngram generation task, written mainly in C:

library(ngram)
print(installed.packages()["ngram", "Version"])

[1] "3.0.4"

ngram_ngram <- function(x, n) {
  ngram(x, n=n)
}

out_ngram <- ngram(hex_s)

# Tidy the result.
out_ngram <- get.phrasetable(out_ngram)
out_ngram <- setNames(out_ngram$freq, out_ngram$ngrams)

all.equal(as.table(out_ngram), out_r, check.attributes=FALSE)

[1] TRUE

Now we benchmark all the above implementations using the large file:

benchmark_all <- function(x) {
  
  x_s <- paste(x, collapse=" ")
  x_t <- tokens(x_s)
  x_i <- itoken(x_s, progressbar=FALSE)
  
  microbenchmark(
    ngram_r(x, 2),
    ngram_rcpp(x, 2),
    ngram_quanteda_1(x_s, 2),
    ngram_quanteda_2(x_t, 2),
    ngram_text2vec_1(x_s, 2),
    ngram_text2vec_2(x_i, 2),
    ngram_ngram(x_s, 2),
    times=3
  )
}

x <- big_hex_c[1:as.integer(length(big_hex_c) / 10)]  # Smaller.
ben_all <- benchmark_all(x)
print(ben_all)

Unit: seconds
                     expr      min       lq     mean   median       uq      max neval
            ngram_r(x, 2) 7.534342 7.548868 7.688973 7.563393 7.766289 7.969184     3
         ngram_rcpp(x, 2) 1.157837 1.261585 1.353972 1.365333 1.452040 1.538746     3
 ngram_quanteda_1(x_s, 2) 3.369109 3.544664 4.206217 3.720219 4.624771 5.529324     3
 ngram_quanteda_2(x_t, 2) 1.936233 1.944721 1.960462 1.953209 1.972577 1.991945     3
 ngram_text2vec_1(x_s, 2) 1.405728 1.438129 1.466014 1.470530 1.496158 1.521786     3
 ngram_text2vec_2(x_i, 2) 1.313886 1.629596 1.965398 1.945306 2.291154 2.637002     3
      ngram_ngram(x_s, 2) 3.727586 3.781383 3.876381 3.835181 3.950779 4.066378     3

3.2 Moving-Window Calculation

Task like n-gram generation is a type of moving-window operation. Such task is in general harder to vectorize and hence will usually result in performance issue in R. In this section we generalize the discussion to a generic moving window operation. For illustration we will use a sample log file from a Apache web server publicly available from Loghub (Zhu et al. (2019)). The log will look like the following text lines:

# https://github.com/logpai/loghub/tree/master/Apache
target_file <- "https://raw.githubusercontent.com/logpai/loghub/master/Apache/Apache_2k.log"
downloaded_file <- file.path(TEMPDIR, basename(target_file))

if ( !file.exists(downloaded_file) ) {
  download.file(target_file, downloaded_file)
}

logs <- readLines(downloaded_file, warn=FALSE)
print(head(logs))

[1] "[Sun Dec 04 04:47:44 2005] [notice] workerEnv.init() ok /etc/httpd/conf/workers2.properties"
[2] "[Sun Dec 04 04:47:44 2005] [error] mod_jk child workerEnv in error state 6"
[3] "[Sun Dec 04 04:51:08 2005] [notice] jk2_init() Found child 6725 in scoreboard slot 10"
[4] "[Sun Dec 04 04:51:09 2005] [notice] jk2_init() Found child 6726 in scoreboard slot 8"
[5] "[Sun Dec 04 04:51:09 2005] [notice] jk2_init() Found child 6728 in scoreboard slot 6"
[6] "[Sun Dec 04 04:51:14 2005] [notice] workerEnv.init() ok /etc/httpd/conf/workers2.properties"

print(length(logs))

[1] 2000

Suppose the task is to count number of records by log type (notice or error) within the last 5 minutes for every 1 minute. This is an overlapping moving window operation where the actual window size (number of records involved in each time-fixed window) must be further determined by the original timestamp. The idea is rather simple but the actual implementation presents several challanges from the angle of performance optimization in native R code.

Let’s tidy up the log lines first to have a cleaner input to our task:

apache_logs <- as.data.table(stringr::str_match(logs, "^\\[([a-zA-Z0-9: ]+)\\] \\[([a-z]+)\\]")[,-1])
setnames(apache_logs, c("ts", "type"))
apache_logs <- apache_logs[, ts:=as.POSIXct(ts, format="%a %b %d %H:%M:%S %y", tz="UTC")]
apache_logs <- apache_logs[, sec:=as.integer(ts)]

head(apache_logs)

Native R

Again let’s try the naive approach using plain R code first.

mw_count_r <- function(DT, wsize, inter, select=NULL) {
    require(data.table)
    sts <- min(DT$sec)
    ninter <- ceiling(diff(range(DT$sec)) / inter)
    if ( !is.null(select) )
        DT <- DT[type %in% select]
    res <- list()
    for ( ca in unique(DT$type) ) {
        sts_ <- sts
        cnt <- integer(ninter)  # Pre-allocate memory.
        for ( i in 1:ninter ) {
            cnt[i] <- nrow(DT[type == ca][sec <= sts_ & sec > sts_ - wsize])
            sts_ <- sts_ + inter
        }
        res[[ca]] <- cnt
    }
    res
}

system.time(mw_res_r <- mw_count_r(apache_logs, 300, 60))

   user  system elapsed
 38.547   0.078   6.459 

The results are two numeric series, which can be plotted as below.

mv_res <- as.data.table(mw_res_r)
mv_res <- mv_res[, i:=.I]
mv_res <- melt(mv_res, id.vars="i", variable.name="log_type", value.name="count")

ggplot(mv_res, aes(x=i, y=count, group=log_type, color=log_type)) +
  geom_line() +
  labs(x="5-Min Moving Window with 1-Min Interval",
       y="Number of Logs")

One should notice that this is a very small dataset but our function already exhibits bottlenecked. Before we proceed to implement the same idea but in Rcpp, let’s struggle a bit more:

mw_count_r_2 <- function(DT, wsize, inter, select=NULL) {
    require(data.table)
    all_secs <- data.table(sec=min(DT$sec):max(DT$sec))
    if ( !is.null(select) )
        DT <- DT[type %in% select]
    res <- lapply(split(DT, DT$type), 
                  function(x) {
                      x <- merge(x, all_secs, by="sec", all=TRUE)
                      sec_unit_cnt <- x[, .(cnt=sum(!is.na(type))), by="sec"]
                      sec_unit_cumcnt <- c(rep(0, wsize), cumsum(sec_unit_cnt$cnt))
                      as.integer(tail(sec_unit_cumcnt, -wsize) - head(sec_unit_cumcnt, -wsize))
                  })
    lapply(res, function(x) x[seq(1, length(x), inter)])
}

system.time(mw_res_r_2 <- mw_count_r_2(apache_logs, 300, 60))

   user  system elapsed
  0.484   0.000   0.248 

identical(mw_res_r$notice, mw_res_r_2$notice) && 
  identical(mw_res_r$error, mw_res_r_2$error)

[1] TRUE

Apparently we can improve drastically with native R code if we do some creative twisting of the original task. Here we re-work the original problem by forcely expanding the input data.frame to have at least one row for every second in the covered period, and we do a cumulative sum by second to arrive at a second-level moving window and then subset its minute-level results. In doing so we effectively remove the need to do the moving window loop (which R is not good at) and transform the task to use the built-in vectorized primitive cumsum function.

The improved solution still has several flaws:

1. It uses more memory than it actually needs by expanding the data.frame from a sparse representation (a second without log will not present) into a dense one; this can be a bottleneck for really huge datasets
2. It computes things more than we need: moving window operation is performed at second level but we only need a minute level (This can be solved if we refactor the code a bit, so less of an issue and also its performance impact is quite limited)
3. It is much harder to read and understand the intent of the code

And, as a matter of fact, it is still not super fast.

Let’s see if we can achieve even faster computing time by using Rcpp.

Rcpp

cppFunction("
  NumericVector MWCountRcpp(NumericVector ts, int wsize, int inter, int sts, int nt) {
    NumericVector out(nt);
    for (int i = 0; i < nt; i++) {
      NumericVector cnts = ts[(ts <= sts) & (ts > sts - wsize)];
      out(i) = cnts.size();
      sts += inter;
    }
    return out;
  }")

mw_count_rcpp <- function(DT, wsize, inter, select=NULL) {
    require(data.table)
    sts <- min(DT$sec)
    ninter <- ceiling(diff(range(DT$sec)) / inter)
    if ( !is.null(select) )
        DT <- DT[type %in% select]
    lapply(split(DT, DT$type), 
           function(x) as.integer(MWCountRcpp(x$sec, wsize, inter, sts, ninter)))
}

system.time(mw_res_rcpp <- mw_count_rcpp(apache_logs, 300, 60))

   user  system elapsed
  0.141   0.016   0.085 

identical(mw_res_r$notice, mw_res_rcpp$notice) && 
  identical(mw_res_r$error, mw_res_rcpp$error)

[1] TRUE

print(microbenchmark(
  mw_count_r(apache_logs, 300, 60),
  mw_count_r_2(apache_logs, 300, 60),
  mw_count_rcpp(apache_logs, 300, 60),
  times=3
))

Unit: milliseconds
                                expr       min         lq      mean    median        uq       max neval
    mw_count_r(apache_logs, 300, 60) 7081.7981 7213.51185 7525.0865 7345.2256 7746.7306 8148.2357     3
  mw_count_r_2(apache_logs, 300, 60)  257.9965  263.29445  267.7426  268.5924  272.6157  276.6390     3
 mw_count_rcpp(apache_logs, 300, 60)   84.1526   86.41185   90.0316   88.6711   92.9711   97.2711     3

As one can see, using Rcpp is brutal and simple and extremely fast. It also supports vector slicing as in R so writing Rcpp is indeed more like writing R and less like writing C++.

In this toy example our Rcpp solution is more than 3 times faster than our creative twisted approach. And it cerrtainly can be even faster when it comes to larger application.

3.3 Time-Dependent Feature Generation

Another common use case similar to a moving window operation is when we need to calculate derived features based on a look-back time window. This is usually for a machine learning dataset where an entity can have time dependent behavior statistics such as number of a certain activity in the past 30 days at any given point when a model needs to make a prediction about the future.

For this scenario we will create a MOOC-like student-course dataset as our working example:

# Assume user is encoded by integer and class is encoded by a single English letter.
create_mooc_data <- function(nlog, nuser, seed=528491) {
  set.seed(seed)
  dat <- data.table(user=sample(1:nuser, nlog, replace=TRUE, 
                                prob=rbeta(1:nuser, .1, 2)),
                    course=sample(letters, nlog, replace=TRUE, prob=1:26),
                    ts=runif(nlog, as.integer(as.POSIXct("2015-01-01")),
                             as.integer(as.POSIXct("2015-06-30"))))
  dat[, t:=as.POSIXct(ts, origin="1970-01-01")]
  setorder(dat, user, course, ts)
  dat
}

mdata <- create_mooc_data(nlog=1e6, nuser=1e3)
head(mdata)

To interpret the dataset, a row represents a pair of studient-class interaction at a given timestamp.

Now assume we’d like to extract a feature: How many unique courses in total did a student register at the time of his/her last activity record for each course? For example, if student 1’s last activity for course a is today, then how many courses in total (including course a) did she also have at least one interaction record? We may want to use this information as one of the feature to build a model to predict whether a user will drop out from a registered class.

For this use case we will again benchmark with 3 different approaches:

1. The naive R approach
2. The improved R approach using data.table API
3. The naive approach but written in Rcpp

Native R

First the naive approach in native R code:

naive_r_func <- function(DT) {
  require(data.table)
  DT <- copy(DT)
  DT[, `:=`(tsmin=min(ts), tsmax=max(ts)), by="user,course"]
  DT <- unique(DT[, list(user, course, tsmin, tsmax)])
  res <- integer(uniqueN(DT$user) + uniqueN(DT$course))
  cnt <- 1
  for ( i in unique(DT$user) ) {
    tmpdat <- DT[user == i]
    for ( j in tmpdat$course ) {
      res[cnt] <- sum(tmpdat$tsmin <= tmpdat[course == j, tsmax])
      cnt <- cnt + 1
    }
  }
  cbind(DT[, list(user, course)], res)
}

system.time(fe_res_r1 <- naive_r_func(mdata))

   user  system elapsed
 99.000   0.110  16.598 

The idea is simple: we loop over each student and inner-loop over each courses to calculate the required metric. Such huge nested loop is doomed to fail for large application.

The result will look something like:

fe_res_r1[sample(.N, 10)]

So for example user 491 has interacted with 18 courses (or 17 other courses) in total when she last interacted with course t in the records.

To stick with R, let’s improve the performance by using data.table’s special group by each i functionality:

better_r_func <- function(DT) {
    require(data.table)
    DT <- copy(DT)
    DT[, `:=`(tsmin=min(ts), tsmax=max(ts)), by="user,course"]
    tmpdat1 <- unique(DT[, list(user, course, tsmin)])
    tmpdat2 <- unique(DT[, list(user, course, tsmax)])
    setkey(tmpdat1, user)
    setkey(tmpdat2, user)
    cbind(tmpdat1[, list(user, course)], 
          res=tmpdat1[tmpdat2, list(res=sum(tsmin <= tsmax)), by=.EACHI][, res])
}

system.time(fe_res_r2 <- better_r_func(mdata))

   user  system elapsed
  0.594   0.000   0.124 

setkey(fe_res_r1, user)
identical(fe_res_r1, fe_res_r2)

[1] TRUE

The solution is indeed very fast. The only drawback is that the code is hard to read and understand, as it always will be when we try to re-work a problem from its original representation.

Rcpp

Can we beat the performance of a pure data.table approach by simply implement the nested loop using Rcpp?

cppFunction("
  NumericVector rcppFunc(NumericVector first, NumericVector last) {
    int len = first.size();
    NumericVector out(len);
    for (int i = 0; i < len; i++) {
      out(i) = 0;
      for (int j = 0; j < len; j++) {
        if (first(j) <= last(i)) {
          out(i) += 1;
        }
      }
    }
    return out;
  }")

rcpp_func <- function(DT) {
  require(data.table)
  DT <- copy(DT)
  DT[, `:=`(tsmin=min(ts), tsmax=max(ts)), by="user,course"]
  DT <- unique(DT[, list(user, course, tsmin, tsmax)])
  cbind(DT[, list(user, course)],
        res=as.integer(DT[, list(res=rcppFunc(.SD$tsmin, .SD$tsmax)), by="user"][, res]))
}

system.time(fe_res_cpp <- rcpp_func(mdata))

   user  system elapsed
  0.516   0.110   0.165 

setkey(fe_res_cpp, user)
identical(fe_res_cpp, fe_res_r1)

[1] TRUE

print(microbenchmark(better_r_func(mdata),
                     rcpp_func(mdata), 
                     times=10))

Unit: milliseconds
                 expr     min      lq      mean   median       uq      max neval
 better_r_func(mdata) 90.9778 94.6369  98.63532  98.3861 102.4294 108.4087    10
     rcpp_func(mdata) 83.8232 85.3324 131.85706 137.3354 160.0626 174.7365    10

It turns out that, Rcpp is still the fastest! But the difference is not significant. We shouldn’t be surprised if we know that data.table itself is indeed written in C.

4 Final Wrap-Up

In this notebook we’ve demonstrated the key concept to write a high-performance R program: to vectorize and to avoid memory copy. They will work 90% of the time to make sure our R program is fast enough. But there are use cases where the concept, especially vectorization, is hard to apply. To overcome such problem we can:

1. Re-work the original problem with an alternative algorithm that can be vectorized
2. Survey any existing package that already handles such case efficiently
3. Implement the core of our algorithm in Rcpp

Obviously the last option is the most flexible one. And we’ve also demonstrated 3 different common use cases where Rcpp can gain tremendous improvement in computing time, with not really much effort to implement.

The key to using Rcpp with ease is to identify only the bottleneck in the algorithm (which is usually a very small part of it) and implement only the bottleneck part in Rcpp.

5 References

Bengtsson, Henrik. 2019. R.utils: Various Programming Utilities. https://CRAN.R-project.org/package=R.utils.

Benoit, Kenneth, Kohei Watanabe, Haiyan Wang, Paul Nulty, Adam Obeng, Stefan Müller, and Akitaka Matsuo. 2018. “Quanteda: An R Package for the Quantitative Analysis of Textual Data.” Journal of Open Source Software 3 (30): 774. https://doi.org/10.21105/joss.00774.

Dowle, Matt, and Arun Srinivasan. 2019. Data.table: Extension of ‘Data.frame‘. https://CRAN.R-project.org/package=data.table.

Eddelbuettel, Dirk. 2013. Seamless R and C++ Integration with Rcpp. New York: Springer. https://doi.org/10.1007/978-1-4614-6868-4.

Eddelbuettel, Dirk, and James Joseph Balamuta. 2017. “Extending extitR with extitC++: A Brief Introduction to extitRcpp.” PeerJ Preprints 5 (August): e3188v1. https://doi.org/10.7287/peerj.preprints.3188v1.

Eddelbuettel, Dirk, and Romain François. 2011. “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software 40 (8): 1–18. https://doi.org/10.18637/jss.v040.i08.

Lang, Dawei, and Guan-tin Chien. 2018. Wordcloud2: Create Word Cloud by ’Htmlwidget’. https://CRAN.R-project.org/package=wordcloud2.

Mersmann, Olaf. 2018. Microbenchmark: Accurate Timing Functions. https://CRAN.R-project.org/package=microbenchmark.

Schmidt, Drew, and Christian Heckendorf. 2017. “ngram: Fast N-Gram Tokenization.” https://cran.r-project.org/package=ngram.

Selivanov, Dmitriy, and Qing Wang. 2018. Text2vec: Modern Text Mining Framework for R. https://CRAN.R-project.org/package=text2vec.

Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

———. 2019. Stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.

Zhu, Jieming, Shilin He, Jinyang Liu, Pinjia He, Qi Xie, Zibin Zheng, and Michael R Lyu. 2019. “Tools and Benchmarks for Automated Log Parsing.” In Proceedings of the 41st International Conference on Software Engineering: Software Engineering in Practice, 121–30. IEEE Press.

---

1. And also surprisingly, we don’t really have many such cases in R programming where we really need to resort to its C interface. The languange is doing a very good job in vectorizing the most utilized functions for the majority of programming operations, especially in the data analytical field.↩