Performance, Profiling and Benchmarking

system.time(expr) calls proc.time(), evaluates expr, calls proc.time() again and returns the difference. The result is an object of class proc_time: a numeric vector of length five holding the user and system CPU time of the R process (user.self, sys.self), the elapsed ("real", wall-clock) time (elapsed), and the user and system time of child processes that R has waited for (user.child, sys.child). Its print method folds the child times into the process times and shows three numbers, called user, system and elapsed. The help page defines the first two as the CPU time charged to the process for executing user instructions and for work done by the system on its behalf, and says that child times are not available on Windows and are NA there. The clock resolution is system specific: on Unix-alikes times are rounded down to milliseconds, and on Windows they are typically available to 10 ms (R documentation: proc.time). Two consequences follow. Anything shorter than a few milliseconds cannot be timed once with system.time(), and elapsed can be much smaller than user when a multi-threaded BLAS or other library keeps several cores busy at once.

The difference between CPU time and elapsed time is the first thing to understand about any timing. Sys.sleep() waits without computing: it costs elapsed time and no CPU time. A computation costs both. The script below shows the structure of the result, the sleeping against the computing, that the expression runs in the calling environment (so an assignment inside system.time() is visible afterwards), and the gcFirst argument, whose default TRUE runs a garbage collection just before the evaluation so that its cost is not charged to the expression that happens to trigger it.

Source: R documentation: system.time

If the expression fails, system.time() still reports how far it got, in a message that starts with Timing stopped at: and gives the same three numbers (see the section on common errors below). It is a good tool for a first look and for a long job, because it adds no overhead and no dependency. It is a poor tool for comparing two fast expressions, which is the reason for the next two sections.

A single timing is one sample from a noisy process. The computer is doing other things, the CPU changes its clock speed, the memory caches are warm or cold, the garbage collector runs when it happens to, and the first time anything is used it can cost more than later: a package namespace is loaded, a function is byte-compiled by the JIT (a small closure on its second call, a closure with a loop on its first; see the section on the compiler), a file is read from disk instead of the cache. The remedy is the usual one in measurement: repeat, and report a statistic that does not move with the outliers, the minimum or the median rather than the mean. The script below repeats a computation ten times and shows that the numbers differ, shows that the first load of a package costs far more than the second, and times the parts of a computation separately, because the total of a section of a script hides which of its lines is slow.

Source: R documentation: system.time

Writing the loop by hand (for (i in 1:100) system.time(...) and averaging) is what the next two tools automate, with the repeats, the warm-up, the choice of statistic, the check that the expressions agree, and, in bench, the memory.

bench::mark() is the benchmarking function of the bench package (version 1.1.4 here; bench package documentation: mark). It takes any number of expressions, named or not, evaluates each at least twice, once to measure the memory allocation and store the result, and then as many times as needed to fill min_time seconds, and returns a tibble with one row per expression. Its arguments, with their defaults, are these.

Source: bench package documentation: mark

The columns of the result are expression (class bench_expr), min and median (bench_time), itr/sec (iterations per second), mem_alloc (bench_bytes), gc/sec, n_itr (the number of timed iterations after GC filtering), n_gc (the garbage collections over all iterations, a pseudo-measure of the pressure on the collector), total_time, and the list columns result, memory, time (the raw bench_time vector of each run) and gc. The classes print readably (1.2ms, 7.63MB) and are numeric underneath. The check argument is the feature that matters most for correctness: an optimisation that returns a different answer is a bug, and mark() stops with an error if the expressions disagree. The script below is the comparison of the previous section made properly: the days a loan is overdue, for 120000 loans, computed by a preallocated loop and by vectorised arithmetic. The relative column is rounded to the nearest power of ten because the exact ratio changes from run to run (it was between 1000 and 2000 in the runs that were made), and the script passes filter_gc = FALSE to mark() and to summary() because the loop allocates enough small objects to collect garbage in every iteration, which would otherwise make the function print a warning.

Source: bench package documentation: summary.bench_mark

bench::press() runs mark() across a grid of parameters and binds the results into one tibble, which is how one finds out how a cost grows (bench package documentation: press). Named arguments define the grid, and the one unnamed expression, usually a braced block, runs once per row of the grid with the parameters available as variables. The script below grows a vector with c() and fills a preallocated one, for three sizes, and prints the memory column. That column shows the difference with no noise: doubling the length roughly quadruples the memory of the growing version, because every iteration allocates a new vector of the size so far, and only doubles the memory of the preallocated one. The function bench::system_time() is the single-shot counterpart, a system.time() that returns bench_time values called process (CPU) and real (elapsed).

Source: bench package documentation: press

The script warms both functions up before the first measurement. The first call of a function allocates memory for its byte code, which mark() would otherwise count in the memory of the first expression it measures (without the warm-up the preallocated version at 1000 elements reported a larger figure than at 2000). With a ggplot2 installation, ggplot2::autoplot() of a bench_mark result draws the distribution of the timings of each expression; that is mentioned and not run here.

The microbenchmark package (version 1.5.0, whose help page names Olaf Mersmann as its author, on CRAN; microbenchmark on CRAN) is a community package that predates bench and is still widely used in older answers and tutorials. microbenchmark(..., times = 100L) evaluates each expression times times, in an order chosen by control$order ("random", the default, "inorder" or "block"), after control$warmup warm-up iterations (default 2), using the nanosecond timer of the operating system and subtracting the estimated cost of the timing calls themselves. It returns a data frame of class microbenchmark with one row per evaluation, whose summary() has the columns expr, min, lq, mean, median, uq, max and neval, and whose unit argument ("relative" among others) chooses the unit. Its differences from bench are worth knowing. The check argument defaults to NULL, which skips the check, and takes a function or one of the strings "equal", "equivalent" and "identical". There is no memory column. And its own help page says it is meant for micro-benchmarks of small pieces of code, and recommends the R profiler for hot spots in larger chunks. The script below compares vapply() over 100000 numbers with the vectorised expression, shows the check refusing expressions that disagree, and uses the control list. Its relative median is printed rounded down to a power of ten (the ratio was about 350).

Source: microbenchmark on CRAN

The help page also says that a timing which turns out negative after the overhead is subtracted is replaced by zero with a warning, and that if every evaluation does so, the function stops with an error: the symptom of an expression too fast to measure, such as NULL. The data.table documentation, in its guide to benchmarking, adds a caution that applies to every tool in this section: repeating a benchmark hundreds of times is the right design for atomic calculations, but it is not a good picture of data-processing tasks, where a handful of runs on data large enough to take seconds is more representative, because call overhead and garbage collection dominate the tiny ones (data.table package documentation: Benchmarking data.table).

The guidelines that follow are the ones that remove most wrong conclusions.

The usage of Rprof() is Rprof(filename = "Rprof.out", append = FALSE, interval = 0.02, memory.profiling = FALSE, gc.profiling = FALSE, line.profiling = FALSE, filter.callframes = FALSE, numfiles = 100L, bufsize = 10000L, event = c("default", "cpu", "elapsed")). Enabling profiling disables any profiling that was running, so Rprof(NULL) (or "") stops it; calling Rprof() again with append = TRUE adds to the same file (R documentation: Rprof).

Source: R documentation: Rprof

Two details of what is recorded matter when reading a profile. Only functions that put a context on the call stack appear: some primitives, those of type "special", do not, which is why <-, [ and $ never show up as time even when they take much of it (their cost is charged to the function that called them, and to <GC> when gc.profiling is on and they allocate). And since the profile is of R code, the time spent in compiled code that R calls, .Call() or .Fortran() for instance, is charged to the R function that made the call; the profiler cannot look inside it. The section on profiling compiled code below says what to use for that.

summaryRprof() returns a list: by.self, a data frame sorted by self time, by.total, sorted by total time, sample.interval and sampling.time. Both data frames have the columns self.time, self.pct, total.time and total.pct, with times in seconds and percentages of the total sampled time. Self time is the time when the function was the last one on the stack, the one doing the work. Total time is the time when it was anywhere on the stack, so it includes everything it called. A function with a high total and a low self time is a router, calling something that is slow; a function with a high self time is itself the work. The row names are the function names inside double quotes (the quotes are part of the name). The script below profiles the ISBN check from the scenario code, is_valid_isbn13(), over 600000 strings, and prints the structure of the result and the two tables. The seconds and percentages are replaced with <n>; the function names and their order are what was printed, with one concession: by.self is printed for four rows picked by name, because the first few rows of that table are near ties whose order changes from run to run.

Source: R documentation: summaryRprof

The profile reads top-down in by.total. The whole call, is_valid_isbn13(), is 100% of the time. vapply() is most of it, and FUN, the anonymous function that is_valid_isbn13() passes to vapply(), takes three quarters (R names the function that vapply() calls FUN, which is the name of its argument; an anonymous function that is not an argument shows as <Anonymous>). Below FUN the time is in strsplit() and as.integer(), the work of splitting one ISBN into digits and converting them, and a smaller part, outside FUN, is in gsub() and toupper(), called by normalise_isbn(). The reading is that the function is slow because it splits and converts one string at a time, and the remedy is to do that for all the elements at once. It is also clear from by.self that vapply() itself, the loop, is a small part of the cost: replacing vapply() with a for loop would have bought almost nothing.

The profiler can also be run from the shell on a file that Rprof() wrote: R CMD Rprof prints the same two tables as text. The script below writes Rprof.out in the working directory (the default name) after profiling a small job, and the commands below were run in that directory. Its options (--lines, --total, --self, --linesonly, --min%total=, --min%self=) are listed by R CMD Rprof --help. The output is shown with <n> for the numbers, and the header line with the total run time is likewise edited.

Source: R documentation: Rprof

Source: Writing R Extensions: Profiling R code for speed

Source: Writing R Extensions: Profiling R code for speed

The file that Rprof() writes is plain text, and knowing its shape explains several of the behaviours above. The first line is sample.interval= followed by the interval in microseconds, and every following line is one sample, the call stack with the innermost call first, each function name in double quotes (those are the quotes that appear in the row names of the summary tables). append = TRUE adds new samples to an existing file, which is how a profile of several separate stretches of a program is collected into one summary: the sampling time of the combined file is larger than that of the first session alone. The script below shows the first line, the shape of the others, the effect of append, and that Rprof(NULL) when nothing is running does nothing.

Source: R documentation: Rprof

Two arguments change what the samples contain. gc.profiling = TRUE records, for each sample, whether the garbage collector was running, and the summary then has a <GC> row. A function that allocates all the time, like the growing vector of the earlier sections, spends a visible share of its time there, and the row is the profile’s way of showing allocation as time; without the argument the collector’s time is charged to whatever function it interrupted. filter.callframes = TRUE (R 4.0.0) removes the frames that lazy evaluation and eval() put between a caller and the function that finally evaluates its argument. The help page gives try(EXPR) as the example: its call tree has tryCatch(), tryCatchList(), tryCatchOne() and doTryCatch() between try() and EXPR, and the filter drops those four. The two scripts show each: with gc.profiling the <GC> row appears and without it does not, and with the filter tryCatchList has disappeared from a profile of a function called inside try() while the function itself remains.

Source: R documentation: Rprof

Source: R documentation: Rprof

Reading a profile is a matter of asking the same three questions each time.

Three cautions apply. Profiling has a small cost of its own, and the output file can grow large for a long run at the default interval. A short profile is unreliable: a run of a few hundred milliseconds has only a few dozen samples, so repeat the work (the scripts here loop over a few hundred thousand elements for that reason) or lengthen the data. And the garbage collector runs from time to time and its time is charged to whichever function was running when it started; Writing R Extensions suggests comparing a profile taken straight after a gc() with one taken without (Writing R Extensions: Profiling R code for speed). R-admin adds a fourth, for installations that use a multi-threaded BLAS: profiling uses the SIGPROF signal, which may cause problems with threaded libraries, and the advice is to disable profiling when one is used (R Installation and Administration: Linear algebra).

With line.profiling = TRUE the profiler also records the source line that is being executed, and summaryRprof(lines = "show") adds a by.line table, sorted by source location, whose row names are file#line. lines = "both" mixes function names and locations in one table, and lines = "hide" (the default) shows only function names; basenames is the number of components of the path to display. The condition is that the code was parsed with source references: the help page names keep.source = TRUE in source() and KeepSource: yes in a package’s DESCRIPTION as ways to get them (R documentation: summaryRprof). This matters because Rscript does not keep them. The option keep.source is FALSE unless the session is interactive, so a function defined in a script has no source reference and a line profile of it is empty; the next script shows the check, and the one after it works around it by writing the function to a file and sourcing it with keep.source = TRUE.

Source: R documentation: capabilities

The function in the next script computes a late fee for each loan with a for loop (the file is written under the session’s temporary directory, which keeps its name hotspot.R and has no path in the output). The profile points at line 5, as.integer(end - loan_date[[i]]), which is about four fifths of the samples (75 to 83% in the runs that were made): the subtraction of two dates dispatches to -.Date, which builds a difftime. Line 4, the if that handles unreturned loans, is the second most expensive; line 6, the arithmetic, is so cheap that some runs did not sample it at all, which is why the script picks the two rows by name. Without a line profile the by-function table would have said -.Date and difftime; with it, it names the line of the program that calls them, which is what one edits.

Source: R documentation: summaryRprof

Before R 4.4.0 the profiler measured one thing: CPU time on Unix-alikes, elapsed time on Windows. R 4.4.0 added the event argument and, with it, profiling by elapsed (real, wall-clock) time on Unix (R NEWS: "Rprof() gains an event argument and support for elapsed (real) time profiling on Unix", R NEWS: Changes in R 4.4.0). The help page explains the difference. With "cpu", the default on Unix-alikes and typically the preferred event for finding performance bottlenecks, the profiler counts CPU time of the R process, so it excludes the time R spends waiting for input or for processes started by system(); it can go slower than the elapsed clock when the process often waits, and faster when concurrent threads (OpenMP, a threaded BLAS) compute. With "elapsed", the default and the only supported event on Windows, it counts the time the process is running and executing an R command, including a readline() that waits for input. The help page also says that elapsed profiling on Unix is new and considered experimental: to avoid missing a sample, R tries to use the real-time FIFO scheduling policy with the maximum priority for the internal thread that takes each sample, and if that is not allowed it falls back to the maximum priority of the current policy, which matters only when profiling a machine under high load. "default" is the platform’s choice.

The choice changes what the profile says about code that waits, and the script below is built to show it. wait_then_work() sleeps for a second and then sorts a million random numbers. Sys.sleep() stands in for a slow download or a database query. The CPU profile never sees the sleep, only the sorting (a fraction of a second of CPU in a function that takes over a second), and the elapsed profile shows Sys.sleep() as nine tenths of the time. For a program that is slow because it waits on a network, a disk or another process, the CPU profile is blind to the cause and the elapsed profile is the right one. For a computation, the CPU profile is the cleaner. The script ends with the two errors the arguments can raise, with their real messages: an unknown event is an error from match.arg(), and an interval below the platform minimum is a warning that says what was used instead.

Source: R documentation: Rprof

profvis, described below, follows the same rule: its timing argument defaults to elapsed time where that is possible, that is on Windows or on R 4.4.0 or later, and "cpu" selects CPU time (profvis package documentation: profvis).

Memory use is part of speed, because allocation, copying and garbage collection are time, and the subject has its own page: Memory and Copy Semantics. In brief, there are three profiling tools. Rprof(memory.profiling = TRUE) writes the total allocation in small vectors, large vectors and cons cells at each sample and the number of calls to the internal duplicate function that copies objects, and summaryRprof(memory = "both") adds the change in memory to the timing tables; its caveat is that memory is attributed to the function running at the end of the sampling interval, and a collection can make the memory fall (Writing R Extensions: Memory statistics from Rprof). Rprofmem() is a non-sampling profiler that writes a stack trace each time a large vector is allocated (above a threshold given in bytes) or a new page of the R heap is allocated; it is what bench::mark() uses for its mem_alloc column (Writing R Extensions: Tracking memory allocations). tracemem() marks one object and prints a message whenever it is copied (Writing R Extensions: Tracing copies of an object), and gc() reports the memory in use and triggers a collection (R documentation: gc). The last two need R to have been built with memory profiling, which capabilities("profmem") (shown in the script above) answers. Writing R Extensions notes that memory profiling is difficult mainly because of garbage collection: allocation happens at well-defined times and release at whenever the collector runs.

profvis (version 0.4.0 here, profvis package documentation) draws the data of Rprof() as an interactive page. profvis(expr) runs the expression with line profiling and memory profiling on, and returns an htmlwidget of class profvis that RStudio and Positron show in a viewer and that htmlwidgets::saveWidget() writes to a self-contained HTML file. The page has the source code on top, with the time and memory attributed to each line, and a flame graph below it, where the horizontal axis is time and the vertical axis is the call stack, so the bottom block is the outermost call and every block above it is a function called by the one below; the width of a block is the time spent in it. A Data tab gives the same information as a table that expands the call stack. The documentation lists what the graph cannot show: some internal functions do not appear (<-, [ and $ among them, because they are not on the call stack, though their side effects show up as <GC> blocks), calls through an anonymous function appear as <Anonymous> unless the function is saved in a variable, white blocks are code for which profvis does not have the source (base R and packages), and nothing in C, C++ or Fortran is visible, only the R function that called it (profvis package documentation: Getting started with profvis, profvis package documentation: Frequently asked questions). Memory figures on a line are the allocation between the previous sample and this one, which may have happened on an earlier line.

The arguments of profvis() that matter are expr, the expression to profile (turned into the body of a function, so variables created inside it are not visible outside), interval (default 0.01, and values below 0.005 will probably not give accurate timings), prof_output (a file or directory in which to keep the Rprof() data; the default is a temporary file removed on exit), prof_input (the path to an Rprof() data file produced elsewhere, in which case there is no expr), timing (the event, described above), simplify (the equivalent of filter.callframes), rerun (run again until a sample is recorded, for expressions that return before the first sample), torture (garbage-collect after every so many allocations, which makes memory figures more exact) and split ("h" or "v", the orientation of the divider). A terminal has no viewer, so the script below does what a script, a test or a CI job can do: it creates the widget, checks its class, saves it to a temporary HTML file, and builds a second one from an Rprof() file with prof_input. The flame graph itself is an interactive page and is described, not shown. For the first widget, which profiles the loop of the benchmarking section for 36000 loans, the data behind it has late_loop at the bottom of the stack, -.Date above it with difftime and as.POSIXct.Date above that, and [[.Date with .Date as the other tower: the Date arithmetic inside the loop is what the graph is made of. Under Rscript the code has no source references (keep.source is FALSE, as the script above showed), so the flame graph has function names only and the code panel has no per-line figures; start R interactively, or run options(keep.source = TRUE) earlier in the script, which was checked to bring them back, to get them.

Source: profvis package documentation: profvis

In RStudio the viewer’s save icon writes an .Rprofvis file that RStudio opens and that, renamed to .html, opens in a browser (profvis package documentation: Frequently asked questions).

R’s profiler stops at the boundary of compiled code, so a function whose self time is large and whose work is in C, C++ or Fortran (a .Call() into a package, a BLAS routine) needs a different tool. Writing R Extensions devotes a section to it, with the hints "gleaned from various R users" for Linux (perf, oprofile and operf, sprof for a shared object, and Google’s gperftools), macOS and Windows, with the warning that this is highly system specific and that the chapter is based on reports that may not be current (Writing R Extensions: Profiling compiled code). The tools work best when R and the packages were built with debugging symbols. For the common case, a BLAS call, perf record R -f script.R followed by perf report --sort=dso shows how the time divides between R itself, the C library, the BLAS shared object and the packages' shared objects, which is how one finds out that the time was spent in libRblas or libopenblas and not in the interpreter.

A vectorised function runs its loop in compiled code; a loop in R pays the interpreter’s cost for every element. The effect is the largest single win in everyday R, and Vectorisation and Iteration explains the functions that replace the common loops. The measurement here is of a kind that is easy to miss: a loop that looks up each loan’s book in a table. The loop calls which(keys == isbn[[i]]) for every loan, which scans the 16 keys every time and allocates a logical vector each time; titles[match(isbn, keys)] builds a hash table of the keys once and looks every loan up in it. The difference is about two orders of magnitude for 60000 loans. The relative column is printed with one significant digit (the ratio was between 100 and 125 in the runs that were made).

Source: R documentation: match

Vectorising is not always possible. A loop whose step depends on the result of the previous step, like the running stock of a shelf below, has no vectorised form; the sections on the compiler and on C++ are for it. Even then the vectorised parts of the loop body are worth keeping out of the loop: compute everything that does not depend on the previous step in one pass before the loop.

Growing a result one element at a time (out <- c(out, value)) allocates a new vector and copies the old one on every iteration, so the total work grows with the square of the length. The section on bench::press() showed it in memory (doubling the length quadrupled the allocation of the growing version). Vectorisation and Iteration covers the idiom in the section called "Preallocation and the growing-vector trap" and the assignment-beyond-the-end behaviour that R 3.4.0 added; the rule is numeric(n), integer(n), character(n), logical(n) or vector("list", n) before the loop, x[[i]] <- value inside it, and, when the length is unknown, a list of the pieces combined once at the end.

The same applies to data frames and is much more costly: rbind(df, row) inside a loop copies the whole data frame and checks the columns every time. The script below builds a 2000-row data frame three ways: rbind() in a loop, a list of one-row data frames combined by one do.call(rbind, parts), and the columns built as vectors and wrapped in data.frame() once. The time ratios are printed rounded down to a power of ten, and the memory ratios to the nearest. The loop allocates a thousand times the memory of the vector version, and the list version about ten times.

Source: R documentation: cbind and rbind

R’s functions behave as if arguments were copied (call by value), but a copy is made only when an object is modified while something else refers to it; Writing R Extensions explains the design with a package example in which the single line x$call <- x$printcall copies the whole object x, and shows that data frames are slower than matrices and that "this is an example of why" (Writing R Extensions: Profiling R code for memory use). The consequences for speed are three. Modifying an argument inside a function copies it, however large. Setting one cell of a data frame column in a loop, df$x[i] <- value, goes through $<-.data.frame and copies the column, and more, at every iteration, whereas the same assignment on a plain vector is done in place. And a function that takes a large object only to read a small piece of it costs nothing for the object. The script below measures the second and the first. For 10000 cells the data-frame version allocates ten thousand times the memory of the vector version, and the function that modifies its million-element argument allocates the 7.63 MB of a copy of it, while the function that only reads it allocates nothing. The remedy for the data frame is to do the loop on vectors and build the data frame once, or to work with a matrix or a list; the remedy for the argument is to make the modification a new, smaller object, or to accept the copy knowingly. Memory and Copy Semantics is the place for the rules (what counts as "something else refers to it", reference counting since R 4.0.0, and when R can modify in place) and for tracemem() to see a copy happen.

Source: Writing R Extensions: Tracing copies of an object

A very common piece of advice is to replace loops by sapply() or lapply() because they are faster. They are not, by much: the apply functions are loops written in C that call an R function once per element, and the call is the cost. The profvis documentation says the same thing, and gives it as the example of a belief that is not true. The script below applies a small fee function to 120000 loan durations with a preallocated for loop, vapply(), sapply(), lapply() plus unlist(), and the vectorised pmax() form, and prints two comparisons (the four loops are within a factor of three of each other, and the vectorised form is more than ten times faster than the best of them; the factor was a hundred in the runs made) and the memory ratios, which are the same on every run. What the apply functions offer is not speed but clarity and safety, and vapply() is the safe one: it declares the type and length of each result, FUN.VALUE, and an element that does not conform is an error, not a surprise. sapply() guesses: it simplifies the result to a vector, a matrix or a list depending on what the function returns, and on an empty input it returns an empty list whatever the function would have returned. The script after it shows the two errors that vapply() raises, with their real messages, and the empty-input case.

Source: R documentation: lapply

Source: R documentation: lapply

Matrix operations are where R hands its work to the BLAS and LAPACK (see the next section), and the way they are written still matters, because the library is fast only for the operation it is asked to do. Three rules cover most of it. Do not form an inverse that is only going to be multiplied: solve(a, b) solves the linear system directly, and solve(a) %*% b computes the whole inverse first, which is more work and less accurate. Use crossprod(x) for t(x) %*% x and crossprod(x, y) for t(x) %*% y: it does not build the transpose, and for a single argument it computes only half of the symmetric result (R documentation: crossprod). Prefer a decomposition (qr(), chol()) when the same matrix is used for several right-hand sides, and keep the work in matrices: a data frame converted to a matrix once is faster to compute on than the data frame, as the section on copies showed. The script below solves the normal equations of a least-squares problem with 5000 rows and 500 columns written three ways, checking that they agree, and prints which is fastest, that the naive form is at least one and a half times slower, and the memory ratios. The crossprod form allocates about a seventeenth of the memory of the naive form, which is the cleanest measure of the difference.

Source: R documentation: solve

Matrices and Arrays has the matrix functions themselves, and its section on BLAS and LAPACK the matprod option, which selects how %*%, crossprod() and tcrossprod() treat NaN and Inf: the default checks for them and falls back to a plain loop so that they propagate correctly, and the setting "blas" skips the check and needs care.

R does not implement its own dense linear algebra. The Installation and Administration manual states it at the start of the appendix on linear algebra: the routines make use of BLAS (Basic Linear Algebra Subprograms) and, for most, of LAPACK; the sources contain reference implementations in Fortran, and they can be replaced by external libraries, usually tuned for specific CPUs (R Installation and Administration: Linear algebra). %*%, crossprod(), solve(), chol(), qr(LAPACK = TRUE), eigen() and svd() end up there, and a tuned BLAS multiplies large matrices many times faster than the reference one, partly through CPU-specific code and partly through threads. It is the one speed-up on this page that needs no change to the program.

The choices, as the manual describes them, are these.

Source: R Installation and Administration: BLAS

The ways to change the BLAS follow from where it is linked. If R was built with the shared BLAS (--enable-BLAS-shlib, the default on all platforms but AIX), every reference goes through libRblas, so the library can be swapped without rebuilding R or the packages: by symlinking an external shared BLAS over libRblas, or, on Debian-based systems, by building R with --with-blas and using update-alternatives; the manual points to the Debian wiki page on linear algebra libraries and says that ATLAS, MKL and OpenBLAS alternatives are available (R Installation and Administration: Shared BLAS). Since 2025 this is less free: the BLAS gained routines (dgemmtr, zgemmtr) that the bundled LAPACK 3.12.1 uses, so a replacement BLAS must contain them (OpenBLAS does from version 0.3.29) or R must be built with --with-2025blas, and the NEWS of R 4.5.0 says that 3.12.1 is expected to be the last bundled version, after which external libraries will be required (R NEWS: Changes in R 4.5.0). Binary distributions make the choice for the user: the NEWS of R 4.5.0 says Windows builds use the internal BLAS by default, the manual says some CRAN builds for macOS can be switched to Accelerate, and the Linux image used for this page runs on OpenBLAS.

R can say what it is using. extSoftVersion() returns a named character vector whose BLAS entry names the library file, La_library() returns the file that holds LAPACK and La_version() its version, and sessionInfo() prints a "Matrix products" line (the matprod option), a BLAS: line and a LAPACK: line with the version (R documentation: extSoftVersion, R documentation: La_library, R documentation: sessionInfo). The values are file paths and machine specific, so the script below asks questions about them (is the entry present, is it empty, does the LAPACK file name mention OpenBLAS) and shows the sessionInfo() lines with the paths replaced by <path>, and the answer to the last question, TRUE, is a fact about the machine that ran the examples and not about R: it is an OpenBLAS build, and on a build with the reference BLAS the answer is FALSE. The same script shows the effect of the threads. It multiplies two 1500 by 1500 matrices and checks that the CPU time (user.self) is larger than the elapsed time, which is only possible if several cores worked at once, and then times the same multiplication in two fresh R processes started by callr with OPENBLAS_NUM_THREADS set to 1 and to 4 (the variable is read when the library loads, so setting it in a running session is too late) and prints whether one thread is slower. On a computer with one core the last answer would be FALSE.

Source: R Installation and Administration: Linear algebra

The manual lists the caveats of an external BLAS, and all of them matter for a decision. The results are not identical to the reference implementation; some differences are benign, such as the signs of singular vectors, but an optimised library can be less accurate, and its LAPACK can be older. R relies on IEEE 754 behaviour, and a library that assumes a term with a zero factor is zero can break it, since x * 0 can be NaN. A multi-threaded BLAS can use more CPU time, and occasionally much more elapsed time, than a single-threaded one, notably on a busy machine and on hyper-threaded CPUs, because the threads contend for the CPU caches. BLAS and LAPACK calls inside threaded code (OpenMP sections in a package such as mgcv) may be unsafe with external libraries, whose thread-safety is not guaranteed. A threaded BLAS and the profiler disagree (the SIGPROF signal, noted in the section on profiling). And a threaded BLAS inside a process that is itself one of several parallel workers oversubscribes the machine: set OPENBLAS_NUM_THREADS=1 (or OMP_NUM_THREADS=1) in the workers, a point Parallel and Asynchronous Computing returns to. The links to the projects are OpenBLAS and FlexiBLAS; for MKL and Accelerate the manual’s sections are the reference (R Installation and Administration: Intel MKL, R Installation and Administration: Accelerate). The Fortran Reference page BLAS, LAPACK and Numerical Libraries covers the libraries from the Fortran side.

R evaluates functions by walking their syntax tree, or, if they are compiled, by running a compact byte-code form of them on a virtual machine. The compiler package, part of R since 2.13.0, compiles closures (cmpfun()), expressions (compile()) and files (cmpfile()), and disassemble() prints what it produced. Since R 3.4.0 a just-in-time compiler is on by default: the NEWS says that functions are compiled on first or second use and that top-level loops are compiled and then run (R NEWS: Changes in R 3.4.0); and since R 3.5.0 all packages are byte-compiled on installation by default (R NEWS: Changes in R 3.5.0). Vectorisation and Iteration has the section "Byte-compilation" with the first look at all this; this section is about checking what the compiler did and how it is controlled.

The levels of enableJIT(level) are these (the help page, R documentation: compile); the function returns the previous level, a negative argument returns the current one, and level 3 needs the compiler option optimize to be 2 or 3.

Source: R documentation: compile

The level of a new R process is set by the environment variable R_ENABLE_JIT, and compilePKGS(TRUE) (or _R_COMPILE_PKGS_ set to a positive integer) makes the compiler compile packages when they are installed, but the help page says not to enable it outside installation, because it makes any serialized function be compiled, with a time and space overhead; R_COMPILE_PKGS is the variable that tells R CMD INSTALL whether to compile. A package can opt out in its DESCRIPTION with ByteCompile: no, and an installation can with R CMD INSTALL --no-byte-compile (Writing R Extensions: The DESCRIPTION file, R Installation and Administration: Byte-compilation). The compiler options are optimize (an integer from 0 to 3; the out-of-the-box default is 2), suppressAll, suppressUndefined and suppressNoSuperAssignVar, read by getCompilerOption() and set by setCompilerOptions(). Two of its details are visible in everyday use. The compiler does not compile code that contains explicit calls to browser(), so that one can single-step from there (R NEWS, 3.4.0). And the help page remarks that the "byte code" is, when R is compiled with gcc, actually threaded code rather than byte code.

How to tell what happened? Printing a compiled function ends with a <bytecode: ...> line that holds a memory address, so the script defines is_compiled(), which asks whether such a line is printed, and uses it to show the three facts about the default JIT: the functions of an installed package are compiled; a closure with a loop is compiled by its first call; a small closure by its second. It ends with cmpfun() and disassemble(), whose result is a list of the form (.Code, instructions, constants). For x + 1 the instruction sequence is GETVAR.OP (push x), LDCONST.OP (push the constant), ADD.OP and RETURN.OP, and the operands (1L, 2L, 0L) are positions in the constant pool that is the third element: x + 1, x and 1.

Source: R documentation: compile

The effect of compilation can be measured, if the JIT is kept from compiling the function that is meant to be the interpreted one. The next script switches the JIT off with enableJIT(0L), which returns the previous level (3), defines a loop, compiles a copy of it by hand with cmpfun(), compares the two with bench::mark(), and restores the level. With the JIT back on, the "interpreted" function would be compiled by its first call, so the comparison has to happen with the JIT off. The compiled copy is about ten times faster on a loop of a million scalar additions (the ratio was between 8.6 and 9.3 in the runs made, and is printed rounded to the nearest power of ten). That is the entire gain of the compiler and it is large for scalar loops, but it does nothing for work that R already does in C: the loop of the first benchmark on this page was compiled by the JIT and was still a thousand times slower than the vectorised form. The script ends with R_ENABLE_JIT, the variable that sets the level for a new R process, shown through callr: a child process started with 0 reports level 0, and one started with 1 reports level 1.

Source: R documentation: compile

Two more entry points of the package complete the picture, and both are rarely needed. compile() turns an expression into a byte-code object (type "bytecode") that eval() runs, and cmpfile() compiles a file of R code into an .Rc file, writing a line that names the file, which loadcmp() later loads and runs. They exist for tools that distribute compiled code; ordinary programs and packages do not need them. The script also reads and sets the compiler options: setCompilerOptions() returns a named list of the previous values, so the old setting can be restored.

Source: R documentation: compile

Calling cmpfun() by hand is therefore rarely needed: the JIT compiles what is called, and installation compiles packages. It remains useful in two cases, the one above where the interpreted time must be measured, and when the JIT is switched off on purpose (R_ENABLE_JIT=0) to debug a problem that might be the compiler’s. Compiled code changes the format of tracebacks and some messages, which is why the NEWS of R 3.5.0 mentioned it among the user-visible changes.

data.table (version 1.18.6.1 here) is a package for tables that are large or whose operations are heavy: grouping, joins, sorting and reading and writing files. It is fast because of how it works: it modifies columns by reference instead of copying, it uses its own radix sort and optimised grouping (GForce), it can build indexes and keys for binary-search joins, and parts of it (fread(), fwrite(), forder(), grouping) use OpenMP threads. data.table covers the package, including its section on printing and threads; this page shows only where it changes the measurements. The script below, run on the shelf loans repeated 10000 times (600000 rows), makes three comparisons. A grouped mean over the eight members: here tapply() is as fast as data.table (the script prints that only aggregate() is at least three times slower), which is the honest result for few groups. The same mean over 100000 groups, where tapply() and split() followed by vapply() are at least ten times slower than data.table (the factor was 20 to 35 in the runs made, printed rounded down to a power of ten). Reading a 600000-row CSV file, where fread() is more than ten times faster than read.csv() (about 16 in the runs made). The R FAQ gives the reason for the third part: read.table() reads everything as character data first and then works out what the columns are, which takes time and memory for a large data set, and giving colClasses helps substantially (R FAQ 7.28: Why is read.table() so inefficient?). fread() detects the types itself and in parallel.

Source: data.table package documentation: setDTthreads and getDTthreads

The number of threads is the thing to know before benchmarking: getDTthreads() reports it, setDTthreads(n) changes it for the session and returns the previous value (setDTthreads(0) is all logical CPUs), and R_DATATABLE_NUM_THREADS and R_DATATABLE_NUM_PROCS_PERCENT set it at load time. The script ends by setting one thread, reading it back and restoring the setting. The default is documented as half of the logical CPUs, and the script does not print the number because it is a fact about the machine. The package’s own guide to benchmarking gives rules that carry over to any tool: clear the file cache before timing fread(); benchmark index creation and the indexed query separately; benchmark set* functions only on the first run, because they update their input by reference and later runs see the processed table; prepare the inputs in the class each tool expects so that the coercion is not timed; prefer a few runs on large data over many runs on small data; keep the number of threads in mind; inside a loop prefer set() to := and setDT() to data.table(); and where the benchmark is meant to be published, time the atomic steps separately as well as the end-to-end pipeline (data.table package documentation: Benchmarking data.table). Attaching data.table masks the %notin% that base R gained in 4.6.0 and prints a message about it; the script loads the package with suppressPackageStartupMessages().

When the profile points at a loop that cannot be vectorised, because each iteration depends on the previous one, and the byte-code compiler has already done what it can, the next step is to write that loop in a compiled language. Rcpp (version 1.1.2 here, Rcpp project site, Rcpp on CRAN) is the usual way for C++: Rcpp::cppFunction("code") compiles a function from a string and makes it callable from R, Rcpp::sourceCpp("file.cpp") does the same for a file whose exported functions carry a // [[Rcpp::export]] comment, and a package that contains such code lists Rcpp under LinkingTo and Imports in its DESCRIPTION. cpp11 is an alternative C++ interface to R. The compilation uses the tool chain of the machine: gcc and g++ 13.3 here, Rtools on Windows, the Xcode command line tools on macOS. The shelf example is a stock level that rises with returns, falls with loans and stays between zero and a capacity: level = min(capacity, max(0, level + change)), applied one day after another over a million days. The clamping makes it impossible to write with cumsum(), so the R version is a loop. The C++ version is the same loop on a NumericVector. The script checks that both return identical results before it compares them; the C++ function is about a hundred times faster (printed to the nearest power of ten; the ratio was between 90 and 105 in the runs made), at the price of a compiler dependency, a compilation step on every new session for cppFunction() (a package compiles once, at installation), and code in a second language to maintain. Because cppFunction() compiles the text it is given, the script’s compilation takes a few seconds on first use; it prints nothing when it succeeds.

Source: Rcpp on CRAN

The C and C++ sides of the story, .Call(), PROTECT(), registering routines, src/Makevars and the rules for portable and compliant code, are the subject of Interfacing C, C++ and Fortran; the languages themselves have their references, C++ Reference, C Reference and Fortran Reference. The order matters: a compiled loop is the most expensive remedy to write, test and keep working across compilers and R versions, so it comes after the cheaper ones.