Working with large data in R

Strategies for working with large data with limited resources

Adapted from a tutorial I gave at a lab group meeting. The full workshop has the runnable code and five more workflows.

When doing data analysis, there is an awkward size of dataset, which is small enough that it feels unnecessary to move over to a supercompute cluseter, but where analyzing it locally on your computer can be a serious time-sink if you use the wrong choice of tool. This small demonstration servers to illustrate what tools should be reached for first in this situation.

The task

The example is a table of fake somatic mutation calls meant to replicate the real PCAWG dataset. 46,115,588 rows, 15 columns, 2,778 tumor samples, about 6.5 GB stored as a tab-separated file and 1.1 GB as Parquet. For each sample I want the subclonal mutation fraction (sMF), the share of that sample’s mutations that are subclonal. Every mutation belongs to a cluster; the cluster with the lowest ID is the clonal one, everything else is subclonal.

It is a simple task. Simple grouping and counting within groups, with a very simple calculation. What makes it a useful test is that the grouped minimum has to be broadcast back across every row before the counting starts, so an engine needs a window function, a join, or two passes. That is exactly the shape where implementations stop resembling each other.

How the timing works

To benchmark these methods, we will time each strategy within its own fresh R process, and the clock is split into four phases:

• Startup, R booting.
• Packages, attaching the libraries the method needs.
• Load, getting the data into whatever form the method computes on. For the in-memory tools that means reading the Parquet file into a data frame. For the engines it means opening a connection or building a lazy scan, which is nearly free, because they do the actual reading during the next phase.
• Compute, producing the 2,778-row answer.

The workshop’s own harness fuses reading and computing into one number, so the tables here come from a re-run with the phases instrumented separately. The script and its raw output are alongside this post.

Memory is sampled the same way, from outside the process. R’s own profilers (profmem, bench::mark()’s mem_alloc, peakRAM) all hook R’s allocator, and DuckDB, Polars, and Arrow are compiled libraries that call malloc themselves. To those profilers they are close to invisible. In the earlier edition of this workshop one DuckDB query had profmem reporting a peak of 0.344 GB where the sampled resident set size was 9.779 GB. The harness polls the process’s RSS from a separate R process instead.

Everything below ran on a Windows machine with 8 physical cores and 16 logical, with every engine capped to 8 threads.

Start with the file format

Before any library changes, the single largest speedup available is not storing six and a half gigabytes of numbers as text.

Format	Write (s)	Size (GB)	× smaller than TSV
TSV (text)	2.59	6.40	1.00
TSV + gzip	23.03	1.70	3.76
Parquet (snappy)	11.44	1.09	5.85
Parquet (zstd)	11.80	1.01	6.32
Arrow IPC	5.29	1.92	3.34
fst	1.72	2.10	3.04
qs2	7.28	0.97	6.60

Six times smaller is the headline, but that isn’t the only advantage of the format. Parquet is binary, columnar, and compressed, stored column by column in independently compressed chunks that each carry min/max statistics. That layout buys three separate things:

• Compression. Values within a column are homogeneous, so they compress far better than a row of mixed types.
• Projection pushdown. Reading 2 of 15 columns reads roughly two-fifteenths of the bytes. A row-oriented format has to read every row in full and discard most of it.
• Predicate pushdown. Row-group statistics let a reader skip whole chunks that cannot contain matching rows, without decompressing them.

The last two belong to the reader, not the file, which is why they show up in DuckDB and Polars and not in base R. All three rows below are Polars reading the same 1.1 GB Parquet file, varying only how much of it the reader is allowed to skip:

Polars query against the Parquet file	Optimization	Elapsed (s)	Peak RSS (GB)
All 15 columns, all rows	none	10.75	13.25
2 columns, all rows	projection pushdown	4.59	10.13
2 columns, filtered rows	projection + predicate	3.80	9.88

The choice of reader matters just as much. Every row below reads the same 6.5 GB TSV, whole, into memory:

Reading the 6.5 GB TSV	Elapsed (s)	Peak RSS (GB)	GB/s
base R read.table	286.64	19.65	0.02
data.table::fread	3.67	15.00	1.78
polars scan_csv	4.59	19.16	1.42
DuckDB read_csv	5.47	17.47	1.19
vroom (materialized)	15.09	21.86	0.43
readr::read_delim	15.19	21.86	0.43

Nearly five minutes against under four seconds.

The same data in a binary format, read whole into memory, is faster again, and the gain costs nothing at read time because the work was done once when the file was written:

Reading the same data back	From	Elapsed (s)	Peak RSS (GB)
arrow::read_parquet	Parquet, 1.09 GB	3.89	12.37
DuckDB to a data frame	Parquet, 1.09 GB	3.69	20.20
fst::read_fst	fst, 2.10 GB	3.23	13.93
arrow::read_feather	Arrow IPC, 1.92 GB	0.79	17.57

Switching a TSV to Parquet is one line of code. Almost nothing else in this post has that ratio of benefit to effort.

Partitioning can help in some cases

A Hive-partitioned dataset splits rows into a directory tree named by the partitioning column, so a query filtering on that column can skip whole directories without opening them. The temptation is to partition everything. Three physical layouts of identical data, against three different questions:

Layout	One sample	One chromosome	Full scan (sMF)
Single file	2.40	2.38	1.13
By chromosome	1.30	0.12	1.75
By samplename	0.17	1.38	1.64

Each layout wins its own column and loses the others. Partitioning by chromosome does nothing for a query grouped by sample, which is exactly the sMF task, and both partitioned layouts are slower than the single file on a full scan because the engine now has to open 22 or 2,778 files instead of one.

Push it further and it stops being a tradeoff and becomes a mistake. Partitioning by sample and chromosome produces 89,303 files averaging 26 KB, takes 182 seconds to write, and makes a full scan roughly five times slower than the single file. Partition on a column you filter by constantly, that has modest cardinality, and accept that you are optimizing for one access pattern at the expense of the others.

The main result

Read the file and compute sMF, sorted by total time. Every figure is the mean of nine runs, with the order rotated so each method takes a turn going first.

Method	Startup	Packages	Load	Compute	Total	SD	Peak RSS (GB)
DuckDB over Parquet	0.18	0.17	0.04	0.12	0.52	0.02	0.14
duckplyr	0.19	0.38	0.13	0.20	0.90	0.03	0.17
polars (streaming scan)	0.17	0.09	0.00	0.69	0.95	0.03	2.87
arrow (Acero)	0.16	0.31	0.01	1.05	1.54	0.03	0.46
polars (eager)	0.30	0.10	0.34	0.85	1.59	0.04	10.31
collapse	0.28	0.28	3.96	2.75	7.27	0.06	7.88
dplyr	0.28	0.36	3.99	2.72	7.34	0.10	7.93
data.table	0.19	0.25	8.60	0.88	9.93	0.16	14.05
base R (aggregate)	0.27	0.20	4.05	34.15	38.68	0.87	8.29

The load column explains the ranking. The engines sit between 0.00 and 0.13 seconds there, because at that point they have done nothing but open a handle. The in-memory tools spend 4 seconds turning the Parquet file into an R object before they can start, and data.table spends 8.6.

data.table shows why the split is worth reporting. Its compute is 0.88 seconds, three times quicker than dplyr or collapse and the fastest of anything working on an R object in memory. Its total is 9.93 seconds, because it spent 8.60 seconds loading, most of that converting the Arrow table into a data.table, which is the slowest load in the table.

Which number matters depends on your situation. If the data is already in memory, and often it is, because you are working interactively and read it in an hour ago, then 0.88 seconds is the honest figure and data.table is an excellent choice. Syntax you already know is worth real time too, and no benchmark measures that. What the total column tells you is the cost of getting there from a cold start, and the memory it takes to hold: 14.05 GB against duckplyr’s 0.17 GB.

That memory column is the real limiter, since on a laptop it can make working with a dataset this large simply impossible.

The two Polars rows show the same effect inside one library. polars (eager) and polars (streaming scan) are the same library computing the same answer. The eager version calls pl$read_parquet(), which materializes all 46 million rows before doing anything, so it pays 0.34 seconds of load and peaks at 10.3 GB. The streaming version calls pl$scan_parquet(), pays almost nothing, and lets the optimizer push the aggregation into the file scan, peaking at 2.87 GB.

One caveat on reading this table: the engines’ compute column includes the file read, because that is inseparable from the query by design, so their load and compute are not comparable line by line with the in-memory rows.

Rotating the order turned out to matter less than I expected. Across all nine passes no method’s startup differs by more than 0.05 seconds between leading and following, so position in the sequence carries no penalty worth reporting.

That is not what the single-pass version of this table showed, where DuckDB posted a 0.77 second startup from going first. The difference is what “first” meant: that run was the first to touch DuckDB’s library since the machine booted, so it came off disk rather than out of the operating system’s cache. It is a cold-start cost, not an ordering effect, and it does not recur. Averaging removes it either way, which is the useful property here.

The standard deviations are small apart from base R, whose 34 seconds of compute leave more room to vary.

Key takeaways

Change your storage format when it makes sense. Writing the table once as Parquet costs a line of code and makes every read after it cheaper, and it is the only change here that helps whichever tool you go on to use.

Prefer scan_ to read_, and open_dataset() to reading a file in. When an engine here has a lazy entry point that lets it push work into the file, using that for simple tasks will save lots of time.

If you use dplyr and want it faster without rewriting, use duckplyr. It placed third here on syntax that is already comfortable.

If your data fits comfortably in memory and you already use data.table, that is still a reasonable place to stand. Its aggregation was the second fastest thing in the table. Where this approach can still fail in a restricted environment is if you hit your memory cap.

The general rule is to push the work to the data rather than pulling the data to the work. The follow-up post covers some cases where it breaks down.