R Series: Optimization

We have learnt various R profiling techniques in the previous article. In this twentieth article in the R series, we shall learn how to optimise R code.

We will use version 4.2.2 installed on Parabola GNU/Linux-libre (x86-64) for the code snippets.

$ R --version
R version 4.2.2 (2022-10-31) -- “Innocent and Trusting”
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under the terms of the
GNU General Public License versions 2 or 3.
For more information about these matters see
https://www.gnu.org/licenses/.

Recursive versus non-recursive functions

Consider a recursive function to compute the nth Fibonacci sequence:

recurFib <- function(n) {
if (n == 0) {0}
else if (n == 1 || n == 2) {1}
else { recurFib(n - 1) + recurFib(n - 2) }
}

You can verify the output sequence numbers, as shown below:

> recurFib(0)
[1] 0
> recurFib(1)
[1] 1
> recurFib(2)
[1] 1
> recurFib(3)
[1] 2
> recurFib(4)
[1] 3
> recurFib(5)
[1] 5

A non-recursive implementation is given by the following fibo() function:

> fibo <- function(n) {
+ phi <- (sqrt(5) + 1)/2
+ fib <- (phi^(n+1) - (1-phi)^(n+1)) / (2*phi - 1)
+ round(fib)
+ }

In this example, we are starting the Fibonacci number with 1, as illustrated below:

> fibo(0)
[1] 1
> fibo(1)
[1] 1
> fibo(2)
[1] 2
> fibo(3)
[1] 3
> fibo(4)
[1] 5
> fibo(5)
[1] 8

The system.time() function can be used to obtain the computation time taken by both these functions.

> recurFib(35)
[1] 9227465
> system.time(recurFib(35))
user system elapsed
6.831 0.000 6.834
> fibo(34)
[1] 9227465
> system.time(fibo(34))
user system elapsed
0 0 0

The recursive function is definitely slower than the non-recursive implementation. Moreover, R uses faster native implementation for computing power of values (^), and hence its usage is recommended.

Preallocate

You do not need to allocate memory for any object in R. In the following example, the initialisation takes time.

> v1 <- c()
> system.time(for (i in 1:10000000) {v1[i] <- i;})
user system elapsed
2.143 0.120 2.263

The length can be specified for a vector and memory can be pre-allocated, so that initialisation can be sped up, as illustrated below:

> v2 <- c(NA)
> length(v2) <- 10000000
> system.time(for (i in 1:10000000) {v2[i] <- i;})
user system elapsed
0.299 0.020 0.319

Matrix multiplication

Consider two matrices, m1 and m2:

> m1 <-matrix(1:10000, nrow=100)
> m2 <-matrix(10001:20000, nrow=100)

You can write your own version of matrix multiplication as follows:

> matrixMul <- function(mat1, mat2) {
rows <- nrow(mat1)
cols <- ncol(mat2)
output <- matrix(nrow = rows, ncol = cols)
for (i in 1:rows) {
for(j in 1:cols) {
vec1 <- mat1[i,]
vec2 <- mat1[,j]
mult1 <- vec1 * vec2
output[i,j] <- sum(mult1)
}
}
return(output)
}

The time taken for multiplying the two matrices, m1 and m2, using the above function is given below:

> system.time(matrixMul(m1,m2))
user system elapsed
0.039 0.000 0.039

The recommended practice is to perform matrix multiplication using the ‘%*%’ operator. The result for the multiplication of matrices m1 and m2 is quite instantaneous.

> system.time(m1*m2)
user system elapsed
0 0 0

Bytecode compiler

The cmpfun() function provides an interface to the byte code compiler for R. It takes the following arguments:

The ‘recurFib’ function can be byte-compiled to improve the performance, as shown below:

> library(compiler)
> compiled.recurFib <- cmpfun(recurFib)
> system.time(compiled.recurFib(34))
user system elapsed
4.296 0.000 4.298

The compiled function, along with the reference to the bytecode, is shown below:

> compiled.recurFib
function(n) {
if (n == 0) {0}
else if (n == 1 || n == 2) {1}
else { recurFib(n - 1) + recurFib(n - 2) }
}
<bytecode: 0x55eff06fcba8>

You can also view the disassembly of the compiled function using the disassemble() function, as follows:

> disassemble(compiled.recurFib)
list(.Code, list(12L, GETVAR.OP, 1L, LDCONST.OP,3L,EQ.OP, 4L,
BRIFNOT.OP, 5L, 13L, LDCONST.OP, 3L,RETURN.OP, GETVAR.OP,
1L, LDCONST.OP, 8L, EQ.OP, 9L, OR1ST.OP,7L,30L, GETVAR.OP,
1L, LDCONST.OP, 10L, EQ.OP, 11L, OR2ND.OP, 7L, BRIFNOT.OP,
12L, 36L, LDCONST.OP, 8L, RETURN.OP, GETFUN.OP, 14L, MAKEPROM.OP,
17L, CALL.OP, 15L, GETFUN.OP, 14L, MAKEPROM.OP, 19L, CALL.OP,
18L, ADD.OP, 20L, RETURN.OP), list({
if (n == 0) {
0
}
else if (n == 1 || n == 2) {
1
}
else {
recurFib(n - 1) + recurFib(n - 2)
}
}, n, structure(c(2L, 3L, 4L, 44L, 3L, 44L, 2L, 4L), srcfile = <environment>, class = “srcref”),
0, n == 0, if (n == 0) {
0
...

The cmpfile() function can be used to compile a large source .R file. It accepts the following arguments:

Bigmemory package

The bigmemory package is used to handle large, Big Data matrices using shared memory and memory-mapped files. You can install the package as follows:

> install.packages(“bigmemory”)
Installing package into ‘/home/guest/R/x86_64-pc-linux-gnu-library/4.1’
(as ‘lib’ is unspecified)
--- Please select a CRAN mirror for use in this session ---
...
Creating a generic function for ‘typeof’ from package ‘base’ in package ‘bigmemory’
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded from temporary location
** checking absolute paths in shared objects and dynamic libraries
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (bigmemory)

The ‘m’ matrix is created using the ‘matrix’ function and written to a .csv file using the ‘write.table()’ function, as shown below:

> library(bigmemory)
> m <- matrix(rnorm(30*15000, mean=0, sd=4), 30, 15000)
> write.table(m, file = “/tmp/m1.csv”)
> system.time(write.table(m, file = “/tmp/m1.csv”))
user system elapsed
0.318 0.001 0.318

You can instead use the big.matrix() function from the bigmemory package to improve the write performance, as demonstrated below:

> system.time(big.matrix(m, nrow=30, ncol=15000, backingfile=”m2.csv”, type=”double”, descriptorfile=”m2.desc”))
user system elapsed
0.006 0.000 0.006

The contents of the m2.csv file are given below for reference:

$ cat m2.desc
new(“big.matrix.descriptor”, description = list(sharedType = “FileBacked”,
filename = “m2.csv”, dirname = “/home/guest/20-optimization/”,
totalRows = 30L, totalCols = 15000L, rowOffset = c(0, 30),
colOffset = c(0, 15000), nrow = 30, ncol = 15000, rowNames = NULL,
colNames = NULL, type = “double”, separated = FALSE))

The R code optimizer

The rco package takes in R code as input, and applies various optimization techniques to improve the code. You can install the package using the following command:

> install.packages(“rco”)
Installing package into ‘/home/guest/R/x86_64-pc-linux-gnu-library/4.1’
(as ‘lib’ is unspecified)
--- Please select a CRAN mirror for use in this session ---
...
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (rco)

The formula to convert Fahrenheit to Celsius is given below:

celsius = (fahrenheit - 32) / 1.8

The code snippets are stored in a ‘code’ variable, as shown below:

> library(rco)
> code <- paste(
+ “f <- 100”,
+ “paren <- (f - 32)”,
+ “v <- paren / 1.8”,
+ “c <- v”,
+ “if (c > 30) {“,
+ ‘ print(“It is hot!”)’,
+ “} else {“,
+ ‘ print(“It is pleasant!”)’,
+ “}”,
+ sep=”\n”
+ )

You can invoke the optimize_code() function on the code and view the optimised code, as illustrated below:

> pass1 <- optimize_text(code, iterations = 1)
Optimization number 1

> pass1
[1] “f <- 100\nparen <- (100 - 32)\nv <- paren / 1.8\nc <- v\nif (c > 30) {\n print(\”It is hot!\”)\n} else {\n print(\”It is pleasant!\”)\n}”
> pass2 <- optimize_text(pass1, iterations = 1)
Optimization number 1

> pass2
[1] “f <- 100\nparen <- 68\nv <- 68 / 1.8\nc <- v\nif (c > 30) {\n print(\”It is hot!\”)\n} else {\n print(\”It is pleasant!\”)\n}”

The ‘optimize_text’ function accepts the following arguments:

You are encouraged to read the ‘Help’ pages of the above R functions to learn more about their arguments, options and usage.