24 Improving performance

24.1 Check for existing solutions

  1. Q: What are faster alternatives to lm? Which are specifically designed to work with larger datasets?

    A: Within the Cran task view for HighPerformanceComputing we can find for example the speedglm package and its speedlm() function. We might not gain any performance improvements on small datasets:

    However on bigger datasets it can make a difference:

    For further speed improvements, you might consider switching your linear algebra libraries as stated in ?speedglm::speedlm

    The functions of class ‘speedlm’ may speed up the fitting of LMs to large data sets. High performances can be obtained especially if R is linked against an optimized BLAS, such as ATLAS.

    Note that there are many other opportunities mentioned in the task view, also some that make it possible to handle data which is not in memory.

    When it comes to pure speed a quick google search on r fastest lm provides a stackoverflow thread where someone already solved this problem for us…

  2. Q: What package implements a version of match() that’s faster for repeated lookups? How much faster is it?

    A: Again google gives a good recommendation for the search term r faster match:

    On my laptop fastmatch::fmatch() is around 25 times as fast as match().

  3. Q: List four functions (not just those in base R) that convert a string into a date time object. What are their strengths and weaknesses?

    A: At least these functions will do the trick: as.POSIXct(), as.POSIXlt(), strftime(), strptime(), lubridate::ymd_hms(). There might also be some in the timeseries packages xts or zoo and in anytime. An update on this will follow.

  1. Q: Which packages provide the ability to compute a rolling mean?

    A: Again google r rolling mean provides us with enough information and guides our attention on solutions in the following packages:
    • zoo
    • TTR
    • RcppRoll
    • caTools

    Note that an exhaustive example on how to create a rolling mean function is provided in the textbook.

  2. Q: What are the alternatives to optim()?

    A: Depending on the use case a lot of different options might be considered. For a general overview we would suggest the corresponding taskview on Optimization.

24.2 Do as little as possible

  1. Q: What’s the difference between rowSums() and .rowSums()?

    A: .rowSums() is defined as

    this means, that the internal .rowSums() function is called via .Internal().

    .Internal performs a call to an internal code which is built in to the R interpreter.

    The internal .rowSums() is a complete different function than the “normal” rowSums() function.

    Of course (since they have the same name) in this case these functions are heavily related with each other: If we look into the source code of rowSums(), we see that it is a wrapper around the internal .rowSums(). Just some input checkings, conversions and the special cases (complex numbers) are added:

  2. Q: Make a faster version of chisq.test() that only computes the chi-square test statistic when the input is two numeric vectors with no missing values. You can try simplifying chisq.test() or by coding from the mathematical definition.

    A: Since chisq.test() has a relatively long source code, we try a new implementation from scratch:

    We check if our new implementation returns the same results

    Finally we benchmark this implementation against a compiled version of itself and the original stats::chisq.test():

  3. Q: Can you make a faster version of table() for the case of an input of two integer vectors with no missing values? Can you use it to speed up your chi-square test?

    A:

24.3 Vectorise

  1. Q: The density functions, e.g., dnorm(), have a common interface. Which arguments are vectorised over? What does rnorm(10, mean = 10:1) do?

    A: We can see the interface of these functions via ?dnorm:

    They are vectorised over their numeric arguments, which is always the first argument (x, q, p, n), mean and sd. Note that it’s dangerous to supply a vector to n in the rnorm() function, since the behaviour will change, when n has length 1 (like in the second part of this question).

    rnorm(10, mean = 10:1) generates ten random numbers from different normal distributions. The normal distributions differ in their means. The first has mean 10, the second has mean 9, the third mean 8 etc.

  2. Q: Compare the speed of apply(x, 1, sum) with rowSums(x) for varying sizes of x.

    A: We compare regarding different sizes for square matrices:

    The graph is a good indicator to notice, that apply() is not “vectorised for performance”.

  3. Q: How can you use crossprod() to compute a weighted sum? How much faster is it than the naive sum(x * w)?

    A: We can just give the vectors to crossprod() which converts them to row- and column-vectors and then multiplies these. The result is the dot product which is also a weighted sum.

    A benchmark of both alternatives for different vector lengths indicates, that the crossprod() variant is about 2.5 times faster than sum():

24.4 Old Exercises

  1. Q: How many different ways can you compute a 1d density estimate in R?

    A: According to Deng and Wickham (2011) density estimation is implemented in over 20 R packages.

  1. Q: When the input is a list, we can make a more efficient as.data.frame() by using special knowledge. A data frame is a list with class data.frame and row.names attribute. row.names is either a character vector or vector of sequential integers, stored in a special format created by .set_row_names(). This leads to an alternative as.data.frame():

    What impact does this function have on read_delim()? What are the downsides of this function?

  2. Q: Line profile the following function with torture = TRUE. What is surprising? Read the source code of rm() to figure out what’s going on.

  1. Q: How do the results change if you compare mean() and mean.default() on 10,000 observations, rather than on 100?

    A: We start with 100 observations as shown in the textbook:

    In case of 10,000 observations we can observe that using mean.default() preserves only a small advantage over the use of mean():

    When using even more observations - like in the next lines - it seems that mean.default() doesn’t preserve anymore any advantage at all:

  2. Q: The following code provides an alternative implementation of rowSums(). Why is it faster for this input?

    A:

  3. Q: Imagine you want to compute the bootstrap distribution of a sample correlation using cor_df() and the data in the example below. Given that you want to run this many times, how can you make this code faster? (Hint: the function has three components that you can speed up.)

    Is there a way to vectorise this procedure?

    A: The three components (mentioned in the questions hint) are:

    1. sampling of indices
    2. subsetting the data frame/conversion to matrix (or vector input)
    3. the cor() function itself.

    Since a run of lineprof like shown in the textbook suggests that as.matrix() within the cor() function is the biggest bottleneck, we start with that:

    Remember the outgoing function:

    First we want to optimise the second line (without attention to the cor() function itself). Therefore we exclude the first line from our optimisation approaches and define i within the global environment:

    Then we define our approaches, check that their behaviour is correct and do the first benchmark:

    According to the resulting medians, lower and upper quartiles of our benchmark all three new versions seem to provide more or less the same speed benefit (note that the maximum and mean can vary a lot for these approaches). Since the second version is most similar to the code we started, we implement this line into a second version of cor_df() (if this sounds too arbitrary, note that in the final solution we will come back to the vector input version anyway) and do a benchmark to get the overall speedup:

    Now we can focus on a speedup for the random generation of indices. (Note that a run of linepfrof suggests to optimize cbind(). However, after rewriting cor() to a version that only works with vector input, this step will be unnecessary anyway). We could try differnt approaches for the sequence generation within sample() (like seq(n), seq.int(n), seq_len(n), a:n) and a direct call of sample.int(). In the following, we will see, that sample.int() is always faster (since we don’t include the generation of the sequence into our benchmark). When we look into sample.int() we see that it calls two different internal sample versions depending on the input. Since in our usecase always the second version will be called, we also provide this version in our benchmark:

    The sample.int() versions give clearly the biggest improvement. Since the internal version doesn’t provide any clear improvement, but restricts the general scope of our function, we choose to implement sample.int() in a third version of cor_df() and benchmark our actual achievements:

    As a last step, we try to speedup the calculation of the pearson correlation coefficient. Since quite a lot of functionality is build into the stats::cor() function this seems like a reasonable approach. We try this by working with another cor() function from the WGCNA package and an own implementation which should give a small improvement, because we use sum(x) / length(x) instead of mean(x) for internal calculations:

    In our final benchmark, we also include compiled verions of all our attempts:

    Our final solution benefits most from the switch from data frames to vectors. Working with sample.int gives only little improvement. Reimplementing and compiling a new correlation function adds only minimal speedup.

    To trust our final result we include a last check for similar return values:

    Vectorisation of this problem seems rather difficult, since attempts of using matrix calculus, always depend on building and handling big matrices in the first place.

    We can for example rewrite our correlation function to work with matrices and build a new (vectorised) version of cor_df() on top of that

    However this still doesn’t provide any improvement over the use of lapply():

    Further improvements can be achieved using parallelisation (for example via parallel::parLapply())