| SOCR ≫ | DSPA ≫ | DSPA3 Topics ≫ |
DSPA3: Data Science and Predictive Analytics (UMich HS650)
Chapter 10: Data Formats, Streams, and Computational
Performance
Chapter 10: Data Formats, Streams, and Computational Performance
SOCR/MCAIM (Ivo Dinov)
SOCR/MCAIM (Ivo Dinov)
August 28, 2026
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly) # interactive figures and ALL 3-D graphics
library(data.table)
library(bench) # honest microbenchmarkingHow this chapter uses graphics
Every two-dimensional figure is drawn with
ggplot2and rendered statically. Immediately after each one, the equivalentplot_ly()code appears in a chunk markedeval=FALSE, echo=TRUE.Every three-dimensional figure is drawn with
plot_ly()and evaluated. Computational performance is intrinsically multi-variable — speedup depends on both the parallel fraction and the core count, attainable throughput on both arithmetic intensity and hardware, sketch accuracy on both memory budget and stream length. Each is a surface, and the interaction is the whole point.On portability. This chapter sets no working directory and contains no absolute paths. Every file is written under
tempdir()and every network call is wrapped so that an unavailable service produces a stated fallback rather than a failed build.
1 Learning objectives
After completing this chapter you will be able to:
- Explain the difference between row-major and columnar storage, and predict which operations each favours.
- Benchmark CSV, RDS, Parquet, and Arrow on read time, file size, and column-pruning cost.
- Push predicates down to a database or query engine, and measure the reduction in data transferred.
- Perform larger-than-memory analytics with an out-of-core query engine.
- State and apply reservoir sampling, the count-min sketch, and HyperLogLog, with their error guarantees.
- Evaluate a streaming model prequentially and detect concept drift.
- Locate a computation on the roofline and decide whether it is compute-bound or memory-bound.
- Derive Amdahl’s and Gustafson’s laws and use them to predict the ceiling on parallel speedup.
- Model parallel overhead and compute the core count beyond which parallelism loses.
- Compare summation algorithms on both speed and numerical accuracy, and explain why a faster reduction can be a worse one.
Estimated time: 9–12 hours including exercises. Prerequisites: Chapter 1 (the R toolchain), Chapter 3 (floating point and conditioning, §3.7), and Chapter 8 (the adjusted Rand index, §8.4.2). This chapter is the engineering layer beneath all of them.
2 PART I: FORMATS AND INGESTION
Every previous chapter began with a tidy data frame. Producing one is often where most of the work lives, and the choices made there set hard limits on what is computationally possible afterwards.
3 How data is laid out in memory
A two-dimensional table has to be written into one-dimensional memory, and there are two ways to do it.
\[ \textbf{Row-major: }\ \underbrace{x_{11}\,x_{12}\,x_{13}}_{\text{row }1}\ \underbrace{x_{21}\,x_{22}\,x_{23}}_{\text{row }2}\ \cdots \qquad \textbf{Column-major: }\ \underbrace{x_{11}\,x_{21}}_{\text{col }1}\ \underbrace{x_{12}\,x_{22}}_{\text{col }2}\ \cdots \]
R matrices and data frames are column-major. C arrays and CSV files are row-major. That single fact drives a great deal of practical performance.
Common misconception: “the storage layout is an implementation detail.” It determines which operations are cheap. Summing a column in a column-major layout reads contiguous memory, so each cache line fetched delivers eight useful doubles. Summing a row strides across the array, and every access may pull a fresh 64-byte cache line to use 8 bytes of it — an eightfold waste of memory bandwidth.
The consequence for analytics is direct: most analytical queries touch few columns and many rows, which is exactly the access pattern columnar storage is built for. Most transactional queries touch one row and many columns, which is what row storage is built for.
set.seed(11)
n_l <- 4000
M <- matrix(rnorm(n_l * n_l), n_l, n_l) # R stores this COLUMN-major
bench::mark(
`sum a column (contiguous)` = sum(M[, 2000]),
`sum a row (strided)` = sum(M[2000, ]),
check = FALSE, iterations = 40
)[, c("expression", "median", "mem_alloc")]Same number of additions, same number of doubles, and a substantial timing difference, entirely attributable to memory access pattern.
4 Format conversion
The rio package reads and writes some forty formats
behind three verbs: import(), export(), and
convert(). It infers the format from the file
extension.
library(rio)
# Everything is written under tempdir(); nothing touches the working directory
work_dir <- file.path(tempdir(), "ch10"); dir.create(work_dir, showWarnings = FALSE)
nof1 <- dspa_try(
rio::import("https://umich.instructure.com/files/1760330/download?download_frd=1",
format = "dta"),
fallback = data.frame(ID = rep(1:10, each = 4), Day = rep(1:4, 10),
Tx = rep(0:1, 20), SelfEff = rnorm(40, 30, 5),
SelfEff25 = rnorm(40, 5, 2), WPSS = rnorm(40, 0, 1),
SocSuppt = rnorm(40, 5, 1), PMss = rnorm(40, 2, 1),
PMss3 = rnorm(40, -1, 1), PhyAct = rnorm(40, 200, 60)),
label = "N-of-1 Stata file")
str(nof1)#> 'data.frame': 900 obs. of 10 variables:
#> $ id : num 1 1 1 1 1 1 1 1 1 1 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"
#> $ day : num 1 2 3 4 5 6 7 8 9 10 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"
#> $ tx : num 1 1 0 0 1 1 0 0 1 1 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"
#> $ selfeff : num 33 33 33 33 33 33 33 33 33 33 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"
#> $ selfeff25: num 8 8 8 8 8 8 8 8 8 8 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"
#> $ wpss : num 0.97 -0.17 0.81 -0.41 0.59 -1.16 0.3 -0.34 -0.74 -0.38 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"
#> $ scssuppt : num 5 3.87 4.84 3.62 4.62 2.87 4.33 3.69 3.29 3.66 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"
#> $ pmss : num 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"
#> $ pmss3 : num 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"
#> $ qhyact : num 53 73 23 36 21 0 21 0 73 114 ...
#> ..- attr(*, "format.stata")= chr "%12.0g"paths <- file.path(work_dir, c("nof1.csv", "nof1.rds", "nof1.xlsx", "nof1.parquet"))
names(paths) <- c("csv", "rds", "xlsx", "parquet")
rio::export(nof1, paths[["csv"]])
rio::export(nof1, paths[["rds"]])
ok_xlsx <- !is.null(dspa_try(rio::export(nof1, paths[["xlsx"]]),
label = "xlsx writer"))
# convert() is import() + export() in one call
rio::convert(paths[["csv"]], file.path(work_dir, "nof1_converted.tsv"))
data.frame(file = basename(list.files(work_dir)),
bytes = file.size(list.files(work_dir, full.names = TRUE)))5 Columnar formats
CSV is universal and slow: it is text, so every number must be parsed; it is row-major, so reading one column requires scanning every byte; and it carries no types, so they must be inferred.
Apache Parquet stores data by column, with per-column compression and per-chunk statistics (min, max, null count). Apache Arrow specifies an in-memory columnar layout that many languages can read without copying. Feather is Arrow written to disk.
Three properties follow, and they are what make columnar formats the modern default:
| Property | Consequence |
|---|---|
| Column pruning | Reading 3 of 200 columns costs \(\approx 3/200\) of the I/O |
| Predicate pushdown | Row groups whose min/max exclude the filter are skipped unread |
| Type-aware compression | Run-length, dictionary, and delta encodings beat generic gzip on columnar data |
library(arrow)
set.seed(21)
n_c <- 2e5
big <- data.frame(
id = seq_len(n_c),
group = sample(letters[1:8], n_c, TRUE), # low cardinality: dictionary-encodes
ts = as.Date("2020-01-01") + sample(0:1500, n_c, TRUE),
x = rnorm(n_c), y = rnorm(n_c), z = rnorm(n_c),
flag = sample(c(TRUE, FALSE), n_c, TRUE)
)
for (j in 1:20) big[[paste0("v", j)]] <- rnorm(n_c) # 27 columns total
f_csv <- file.path(tempdir(), "big.csv")
f_rds <- file.path(tempdir(), "big.rds")
f_pq <- file.path(tempdir(), "big.parquet")
f_fth <- file.path(tempdir(), "big.feather")
data.table::fwrite(big, f_csv)
saveRDS(big, f_rds, compress = TRUE)
arrow::write_parquet(big, f_pq)
arrow::write_feather(big, f_fth)
sizes <- data.frame(
format = c("CSV", "RDS (gzip)", "Parquet", "Feather"),
MB = round(file.size(c(f_csv, f_rds, f_pq, f_fth)) / 1e6, 2))
sizes$ratio_vs_csv <- round(sizes$MB / sizes$MB[1], 3)
sizesb_full <- bench::mark(
CSV = data.table::fread(f_csv, showProgress = FALSE),
RDS = readRDS(f_rds),
Parquet = arrow::read_parquet(f_pq),
Feather = arrow::read_feather(f_fth),
check = FALSE, iterations = 5, memory = FALSE)
b_cols <- bench::mark(
`CSV, 3 cols` = data.table::fread(f_csv, select = c("id", "group", "x"),
showProgress = FALSE),
`Parquet, 3 cols` = arrow::read_parquet(f_pq, col_select = c("id", "group", "x")),
`Feather, 3 cols` = arrow::read_feather(f_fth, col_select = c("id", "group", "x")),
check = FALSE, iterations = 5, memory = FALSE)
rbind(
data.frame(task = "read all 27 columns", format = as.character(b_full$expression),
median_sec = round(as.numeric(b_full$median), 4)),
data.frame(task = "read 3 of 27 columns", format = as.character(b_cols$expression),
median_sec = round(as.numeric(b_cols$median), 4)))Read the second block against the first. For CSV, selecting three columns saves little, the parser must still scan every byte of every row to find the delimiters. For Parquet and Feather the saving is close to proportional, because the columns are stored separately and the unread ones are never touched.
bind_rows(
data.frame(format = as.character(b_full$expression),
seconds = as.numeric(b_full$median), task = "All 27 columns"),
data.frame(format = sub(", 3 cols", "", as.character(b_cols$expression)),
seconds = as.numeric(b_cols$median), task = "3 of 27 columns")) |>
ggplot(aes(reorder(format, seconds), seconds, fill = task)) +
geom_col(position = "dodge", width = 0.7) +
coord_flip() +
scale_fill_manual(values = c("All 27 columns" = "#9EB4C8",
"3 of 27 columns" = "#3B7DD8")) +
labs(title = "Column pruning is free in columnar formats and expensive in CSV",
subtitle = "Same data, same columns requested. Only the storage layout differs",
x = NULL, y = "Median read time (seconds)", fill = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = as.character(b_full$expression), y = as.numeric(b_full$median),
type = "bar", name = "All columns") |>
add_trace(x = sub(", 3 cols", "", as.character(b_cols$expression)),
y = as.numeric(b_cols$median), name = "3 columns") |>
layout(title = "Read time by format and column selection",
xaxis = list(title = ""), yaxis = list(title = "Seconds"),
barmode = "group")The advantage is a surface over (rows, columns):
row_grid <- c(2e4, 5e4, 1e5, 2e5)
col_grid <- c(5, 10, 20, 40)
read_time <- function(nr, nc, fmt) {
set.seed(nr + nc)
d <- as.data.frame(matrix(rnorm(nr * nc), nr, nc))
f <- file.path(tempdir(), sprintf("t_%d_%d.%s", nr, nc, fmt))
if (fmt == "csv") data.table::fwrite(d, f) else arrow::write_parquet(d, f)
sel <- names(d)[1:min(3, nc)]
t0 <- Sys.time()
if (fmt == "csv") data.table::fread(f, select = sel, showProgress = FALSE)
else arrow::read_parquet(f, col_select = all_of(sel))
as.numeric(difftime(Sys.time(), t0, units = "secs"))
}
Zcsv <- outer(col_grid, row_grid, Vectorize(\(c_, r_) read_time(r_, c_, "csv")))
Zpq <- outer(col_grid, row_grid, Vectorize(\(c_, r_) read_time(r_, c_, "parquet")))
plot_ly() |>
add_surface(x = row_grid, y = col_grid, z = Zcsv, opacity = 0.9,
showscale = FALSE, colorscale = "Reds", name = "CSV") |>
add_surface(x = row_grid, y = col_grid, z = Zpq, opacity = 0.9,
showscale = FALSE, colorscale = "Blues", name = "Parquet") |>
layout(title = "Time to read 3 columns: CSV (red, upper) vs. Parquet (blue, lower)",
scene = list(xaxis = list(title = "Rows"),
yaxis = list(title = "Total columns"),
zaxis = list(title = "Seconds")))Rotate along the columns axis. The CSV surface climbs steeply as total column count grows, the parser scans every one regardless of what was requested, while the Parquet surface stays nearly flat, because the unread columns are never touched.
6 Relational databases and predicate pushdown
A query is an instruction sent to a database describing what is wanted; the engine decides how to obtain it. That separation is the entire value proposition, and the single most important consequence is predicate pushdown: filtering and aggregating happen where the data lives, so only the result crosses the wire.
# BRFSS_MODE = "sample" builds a synthetic survey table with the same shape as
# the CDC BRFSS extract; "full" reproduces the 315 MB ingest. Neither writes
# outside tempdir().
make_survey <- function(n = 50000, seed = 31) {
set.seed(seed)
race_lv <- c("White", "Black", "Asian", "AmIndian", "Hawaiian",
"Other", "Multiracial", "Hispanic", "Refused")
d <- data.frame(
year = sample(c(2013L, 2015L), n, TRUE),
x_race = factor(sample(race_lv, n, TRUE,
prob = c(.62,.12,.05,.02,.01,.03,.03,.11,.01)),
levels = race_lv),
hlthpln1 = sample(1:2, n, TRUE, prob = c(0.86, 0.14)),
poorhlth = pmax(0, round(rnorm(n, 4, 6))),
genhlth = sample(1:5, n, TRUE, prob = c(.18,.32,.30,.15,.05)),
age = sample(18:95, n, TRUE))
for (j in 1:24) d[[sprintf("q%02d", j)]] <- sample(1:9, n, TRUE)
d
}
survey <- if (BRFSS_MODE == "full" && ONLINE) {
zp <- dspa_download("https://www.socr.umich.edu/data/DSPA/BRFSS_2013_2014_2015.zip",
"BRFSS.zip", timeout = 900)
dspa_try({
ex <- file.path(tempdir(), "brfss"); dir.create(ex, showWarnings = FALSE)
members <- unzip(zp, list = TRUE)$Name # select BY NAME, not position
keep <- grep("2013|2015", members, value = TRUE)
unzip(zp, files = keep, exdir = ex) # extract ONCE, into tempdir
on.exit(unlink(ex, recursive = TRUE), add = TRUE)
do.call(rbind, lapply(keep, \(f) Hmisc::sasxport.get(file.path(ex, f))))
}, fallback = make_survey(), label = "BRFSS archive")
} else make_survey()
c(mode = BRFSS_MODE, rows = nrow(survey), columns = ncol(survey),
MB_in_memory = round(as.numeric(object.size(survey)) / 1e6, 1))#> mode rows columns MB_in_memory
#> "sample" "50000" "30" "6.2"# The outcome is created BEFORE it is summarized
survey$has_plan <- survey$hlthpln1 == 1
c(coverage_rate = round(mean(survey$has_plan), 4))#> coverage_rate
#> 0.8604#> Mode FALSE TRUE
#> logical 6981 43019#>
#> White Black Asian AmIndian Hawaiian Other
#> 0.6264 0.1165 0.0494 0.0195 0.0103 0.0292
#> Multiracial Hispanic Refused
#> 0.0298 0.1091 0.0098t_glm <- system.time(
fit_hcp <- glm(has_plan ~ x_race, data = survey, family = binomial()))
c(seconds = round(t_glm[["elapsed"]], 2))#> seconds
#> 0.09# Coefficients are LOG-odds ratios relative to the reference level.
# Exponentiating gives odds ratios (Chapter 9, Section 9.8).
or_tab <- data.frame(
level = names(coef(fit_hcp))[-1],
log_OR = round(coef(fit_hcp)[-1], 4),
OR = round(exp(coef(fit_hcp)[-1]), 4),
CI_low = round(exp(confint.default(fit_hcp)[-1, 1]), 4),
CI_high = round(exp(confint.default(fit_hcp)[-1, 2]), 4))
head(or_tab, 5)6.1 Doing the work in the database
library(DBI)
library(duckdb)
con <- dbConnect(duckdb::duckdb(), dbdir = ":memory:")
dbWriteTable(con, "survey", survey)
dbWriteTable(con, "USArrests", tibble::rownames_to_column(USArrests, "state"))
dbListTables(con)#> [1] "USArrests" "survey"Common misconception: “the point of a database is that it stores the data.” The point is that it computes.
SELECT * FROM ttransfers the whole table into R and throws away every advantage the engine offers, the database becomes an expensive file.What a query engine gives you is predicate pushdown (filter before transfer), projection pushdown (send only the requested columns), aggregation pushdown (compute the summary server-side and return a few rows), and an optimizer that reorders operations. Push the computation down and the data that crosses the wire is the answer, not the input.
# Aggregation done IN the database: the result is a handful of rows
t_push <- system.time(
agg_db <- dbGetQuery(con, "
SELECT year, hlthpln1, COUNT(*) AS n, AVG(poorhlth) AS mean_poorhlth
FROM survey
WHERE age BETWEEN 45 AND 64
GROUP BY year, hlthpln1
ORDER BY year, hlthpln1"))
# Everything pulled into R first, then aggregated
t_pull <- system.time({
raw <- dbGetQuery(con, "SELECT * FROM survey")
agg_r <- raw |> filter(age >= 45, age <= 64) |>
summarise(n = dplyr::n(), mean_poorhlth = mean(poorhlth),
.by = c(year, hlthpln1)) |> arrange(year, hlthpln1)
})
agg_db c(pushdown_seconds = round(t_push[["elapsed"]], 3),
pull_then_aggregate_seconds = round(t_pull[["elapsed"]], 3),
rows_transferred_pushdown = nrow(agg_db),
rows_transferred_pull = nrow(raw),
transfer_ratio = round(nrow(raw) / nrow(agg_db)))#> pushdown_seconds pull_then_aggregate_seconds
#> 0 0
#> rows_transferred_pushdown rows_transferred_pull
#> 4 50000
#> transfer_ratio
#> 12500The pushdown query transfers a handful of rows; the pull-then-aggregate version transfers the entire table. The ratio is the saving, and it grows linearly with the data.
library(dbplyr)
# dplyr verbs against a database table are TRANSLATED to SQL and executed there.
# Nothing is materialized until collect().
tbl_survey <- tbl(con, "survey")
q <- tbl_survey |>
filter(age >= 45, age <= 64) |>
summarise(n = n(), mean_poorhlth = mean(poorhlth, na.rm = TRUE),
.by = c(year, hlthpln1))
# Inspect the generated SQL before running anything
show_query(q)#> <SQL>
#> SELECT "year", hlthpln1, COUNT(*) AS n, AVG(poorhlth) AS mean_poorhlth
#> FROM (
#> SELECT survey.*
#> FROM survey
#> WHERE (age >= 45.0) AND (age <= 64.0)
#> ) q01
#> GROUP BY "year", hlthpln1show_query() is the habit worth forming: it makes
visible exactly what was pushed down and what was not.
6.2 Out-of-core analytics
When a dataset exceeds RAM, the classical answers were disk-backed data structures. The modern answer is a query engine that streams over columnar files, keeping only the working set in memory.
# Write a partitioned Parquet dataset -- larger than we want to hold at once
ds_dir <- file.path(tempdir(), "survey_ds")
arrow::write_dataset(survey, ds_dir, partitioning = "year", format = "parquet")
list.files(ds_dir, recursive = TRUE)[1:4]#> [1] "year=2013/part-0.parquet" "year=2015/part-0.parquet"
#> [3] NA NAc(total_MB_on_disk = round(sum(file.size(list.files(ds_dir, recursive = TRUE,
full.names = TRUE))) / 1e6, 2))#> total_MB_on_disk
#> 0.77# Query the dataset WITHOUT loading it: arrow streams over the files
ds <- arrow::open_dataset(ds_dir)
class(ds)#> [1] "FileSystemDataset" "Dataset" "ArrowObject"
#> [4] "R6"res_arrow <- ds |>
filter(year == 2015, age >= 45, age <= 64) |>
summarise(n = n(), mean_poorhlth = mean(poorhlth), .by = hlthpln1) |>
collect()
res_arrow# duckdb can query Parquet files directly, with partition pruning
sql <- sprintf("
SELECT hlthpln1, COUNT(*) AS n, AVG(poorhlth) AS mean_poorhlth
FROM read_parquet('%s/**/*.parquet', hive_partitioning = true)
WHERE year = 2015 AND age BETWEEN 45 AND 64
GROUP BY hlthpln1", ds_dir)
dbGetQuery(con, sql)#> note
#> "The 2013 partition is never opened: Hive partitioning prunes it by directory"Partition pruning happens before any file is read.
Because year is encoded in the directory structure, a
filter on year eliminates whole directories, and a filter
on age then uses Parquet’s per-row-group min/max statistics
to skip chunks. Neither requires decompressing the data.
# A dataset deliberately larger than the R session should hold at once
big_dir <- file.path(tempdir(), "big_ds"); dir.create(big_dir, showWarnings = FALSE)
set.seed(41)
for (part in 1:6) {
d <- data.frame(part = part,
g = sample(letters[1:10], 3e5, TRUE),
v = rnorm(3e5), w = rnorm(3e5))
arrow::write_parquet(d, file.path(big_dir, sprintf("part-%d.parquet", part)))
}
c(rows_total = 6 * 3e5,
MB_on_disk = round(sum(file.size(list.files(big_dir, full.names = TRUE))) / 1e6, 1))#> rows_total MB_on_disk
#> 1800000 33t_ooc <- system.time(
out <- arrow::open_dataset(big_dir) |>
filter(v > 1) |>
summarise(n = n(), mean_w = mean(w), .by = g) |>
collect() |> arrange(g))
outc(seconds = round(t_ooc[["elapsed"]], 3),
peak_R_objects_MB = round(as.numeric(object.size(out)) / 1e6, 4))#> seconds peak_R_objects_MB
#> 0.060 0.002The result occupies kilobytes; the source occupies tens of megabytes
on disk and was never fully materialized. This is the pattern
that replaces disk-backed data frames: leave the data in a
columnar file, express the query in dplyr, and let the
engine stream.
7 Web data
7.1 JSON
library(jsonlite)
nof1_json <- dspa_try(
jsonlite::fromJSON("https://umich.instructure.com/files/1760327/download?download_frd=1"),
fallback = nof1, label = "N-of-1 JSON")
class(nof1_json); dim(nof1_json)#> [1] "data.frame"#> [1] 900 10JSON’s {key: value} nesting maps naturally onto R lists;
jsonlite::fromJSON(simplifyDataFrame = TRUE) flattens a
regular array of objects into a data frame. When the structure is
irregular, keep it as a list and use purrr::map_* to
extract fields.
7.2 XML and scraping
library(rvest); library(xml2)
socr <- dspa_try(rvest::read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data"),
label = "SOCR wiki")
if (!is.null(socr)) {
cat("Title:", socr |> html_element("head title") |> html_text(), "\n")
paras <- socr |> html_elements("p") |> html_text()
cat("Paragraphs found:", length(paras), "\n")
cat(substr(paste(paras, collapse = " "), 1, 220), "...\n")
}#> Title: SOCR Data - SOCR
#> Paragraphs found: 29
#> The links below contain a number of datasets that may be used for demonstration purposes in probability and statistics education. There are two types of data - simulated (computer-generated using random sampling) and obs ...Check the terms of service and
robots.txtbefore scraping, and rate-limit requests.polite::bow()handles both automatically. Scraping is a last resort: if an API exists, use it, because it is stable and the HTML is not.
7.3 Semantic queries with SPARQL
RDF represents knowledge as (subject, predicate, object) triples, and SPARQL queries patterns over those triples. Wikidata exposes a public endpoint.
library(WikidataQueryServiceR)
tb_query <- '
PREFIX wd: <http://www.wikidata.org/entity/>
PREFIX wdt: <http://www.wikidata.org/prop/direct/>
PREFIX rdfs: <http://www.w3.org/2000/01/rdf-schema#>
PREFIX p: <http://www.wikidata.org/prop/>
PREFIX qualifier: <http://www.wikidata.org/prop/qualifier/>
PREFIX statement: <http://www.wikidata.org/prop/statement/>
SELECT DISTINCT ?countryLabel ?ISO3Code ?latlon ?prevalence ?year
WHERE {
wd:Q12204 p:P1193 ?prevalenceStatement . # Q12204 = tuberculosis
?prevalenceStatement qualifier:P17 ?country ;
qualifier:P585 ?year ;
statement:P1193 ?prevalence .
?country wdt:P625 ?latlon ;
rdfs:label ?countryLabel ;
wdt:P298 ?ISO3Code .
FILTER (lang(?countryLabel) = "en")
}
ORDER BY DESC(?prevalence) # ORDER BY a BOUND variable
'
tb <- dspa_try(WikidataQueryServiceR::query_wikidata(tb_query),
label = "Wikidata SPARQL endpoint")
if (!is.null(tb)) { cat("Rows returned:", nrow(tb), "\n"); head(tb, 4) }#> Rows returned: 6The ORDER BY clause sorts on ?prevalence,
which the WHERE block binds. Ordering by an
unbound variable is a silent no-op, SPARQL treats
unbound values as mutually incomparable, so every sort key must appear
in the pattern.
if (!is.null(tb) && nrow(tb) > 0) {
# WKT is "Point(longitude latitude)". Split, strip the affixes, and CONVERT
# to numeric -- plot_geo requires numbers, not the character strings that
# gsub leaves behind.
geo <- tb |>
tidyr::separate(latlon, into = c("long", "lat"), sep = " ", remove = FALSE) |>
mutate(long = as.numeric(gsub("Point\\(", "", long)),
lat = as.numeric(gsub("\\)", "", lat)),
prevalence = as.numeric(prevalence)) |>
filter(is.finite(long), is.finite(lat), is.finite(prevalence))
# Quantile binning is guarded: Wikidata coverage is sparse and may return
# too few distinct values for four bins.
n_distinct_prev <- length(unique(geo$prevalence))
n_bins <- max(2, min(4, n_distinct_prev - 1))
brk <- unique(quantile(geo$prevalence, seq(0, 1, length.out = n_bins + 1)))
geo$band <- if (length(brk) >= 3)
cut(geo$prevalence, breaks = brk, include.lowest = TRUE,
labels = paste0("Q", seq_len(length(brk) - 1))) else factor("all")
c(countries = nrow(geo), distinct_prevalences = n_distinct_prev, bins = n_bins)
head(geo[, c("countryLabel", "ISO3Code", "long", "lat", "prevalence", "band")], 5)
}if (exists("geo") && !is.null(geo) && nrow(geo) > 0) {
ge <- list(scope = "world", showland = TRUE, landcolor = toRGB("gray92"),
countrycolor = toRGB("white"), countrywidth = 0.5)
plot_geo(geo, lon = ~long, lat = ~lat, color = ~band,
text = ~paste0(countryLabel, " (", ISO3Code, "): ", prevalence),
marker = list(size = 12, line = list(width = 0.5))) |>
layout(geo = ge, title = "Tuberculosis prevalence reported in Wikidata")
}library(maps)
cities <- maps::world.cities |>
filter(pop > 5e5) |>
mutate(label = sprintf("%s, %s: %.1fM", name, country.etc, pop / 1e6),
# Marker AREA proportional to population, so radius scales as sqrt()
size = sqrt(pop / max(pop)) * 22 + 2)
ge <- list(scope = "world", showland = TRUE, landcolor = toRGB("gray92"),
countrycolor = toRGB("white"), countrywidth = 0.5)
plot_geo(cities, lon = ~long, lat = ~lat, text = ~label,
marker = ~list(size = size, opacity = 0.6,
line = list(width = 0.2))) |>
layout(geo = ge, title = "World cities above 500,000 population")Marker area is proportional to population, so the radius scales as \(\sqrt{\text{pop}}\). Scaling the radius linearly with the value would exaggerate large cities quadratically, the same perceptual point made in Chapter 2, §2.10.
7.4 Random number services
rng_url <- paste0("https://www.random.org/integers/",
"?num=300&min=100&max=200&col=3&base=10&format=plain&rnd=new")
true_rng <- dspa_try(utils::read.table(rng_url),
fallback = as.data.frame(matrix(sample(100:200, 300, TRUE), ncol = 3)),
label = "random.org")
head(true_rng, 3)Common misconception: “true random numbers are better for simulation.” For Monte Carlo work they are usually worse, because they are not reproducible. The property that matters is not physical randomness but statistical adequacy plus a seed: a modern PRNG (Mersenne Twister, PCG64) passes the standard test batteries and regenerates the identical stream from
set.seed(), which is what makes a simulation checkable.Physical entropy matters for cryptography, where predictability is the threat model, and for seeding a PRNG. It is not an upgrade for a simulation study.
8 PART II: NETWORK DATA
9 Graph representations and their cost
A graph \(G=(V,E)\) with \(|V|=n\) vertices and \(|E|=m\) edges admits three standard representations, and the choice governs everything downstream.
| Representation | Memory | Edge query \((u,v)\) | Enumerate neighbours of \(u\) |
|---|---|---|---|
| Adjacency matrix | \(\Theta(n^2)\) | \(O(1)\) | \(O(n)\) |
| Adjacency list | \(\Theta(n+m)\) | \(O(\deg u)\) | \(O(\deg u)\) |
| Edge list | \(\Theta(m)\) | \(O(m)\) | \(O(m)\) |
Real networks are sparse: \(m=O(n)\) or \(O(n\log n)\), far below the \(\binom n2\) maximum. For the Facebook ego network below, \(n\approx4{,}039\) and \(m\approx88{,}000\), so the adjacency matrix would hold \(1.6\times10^7\) entries of which 99.5% are zero.
graph_memory <- function(n, m) {
c(n = n, m = m,
density = signif(2 * m / (n * (n - 1)), 3),
matrix_MB = round(8 * n^2 / 1e6, 1),
list_MB = round(8 * (n + 2 * m) / 1e6, 3),
ratio = round((8 * n^2) / (8 * (n + 2 * m)), 1))
}
as.data.frame(rbind(
`Facebook ego net` = graph_memory(4039, 88234),
`Les Miserables` = graph_memory(77, 254),
`Hypothetical 1e5` = graph_memory(1e5, 1e6)))At \(n=10^5\) the dense matrix needs 80 GB and the adjacency list needs 17 MB. Sparse representation is not an optimization; it is the difference between possible and impossible.
library(igraph)
g_toy <- igraph::make_graph(c(1,2, 1,3, 2,3, 3,4), n = 10, directed = TRUE)
plot(g_toy, vertex.size = 22, vertex.color = "#9EC5E8",
edge.arrow.size = 0.5, main = "Four directed edges among ten vertices")
#> vertices edges isolated
#> 10 4 6Vertices 5–10 appear because n = 10 was declared; they
have degree zero because no edge mentions them. Isolated vertices are a
common artifact of declaring a vertex count larger than the edge list
requires.
10 A real network
fb_path <- dspa_download(
"https://umich.instructure.com/files/2854431/download?download_frd=1",
"facebook_edges.txt")
edges <- dspa_try({
e <- data.table::fread(fb_path, header = FALSE, col.names = c("from", "to"))
as.matrix(e)
}, fallback = {
set.seed(51)
igraph::as_edgelist(igraph::sample_pa(1500, m = 4, directed = FALSE))
}, label = "Facebook ego network")
# The degenerate vertex is IDENTIFIED, not removed by a magic row count
c(min_id = min(edges), rows_touching_min = sum(edges[, 1] == min(edges) |
edges[, 2] == min(edges)))#> min_id rows_touching_min
#> 0 347edges <- edges[edges[, 1] != min(edges) & edges[, 2] != min(edges), , drop = FALSE]
g_fb <- igraph::graph_from_edgelist(edges, directed = FALSE)
c(vertices = gorder(g_fb), edges = gsize(g_fb),
density = signif(edge_density(g_fb), 4),
memory_MB = round(as.numeric(object.size(g_fb)) / 1e6, 2))#> vertices edges density memory_MB
#> 4.0380e+03 8.7887e+04 1.0780e-02 1.4100e+0010.1 Centrality and what it costs
Degree centrality counts a vertex’s edges, \(O(1)\) per vertex from an adjacency list.
Betweenness centrality measures how often a vertex lies on shortest paths:
\[g(v)=\sum_{s\ne v\ne t}\frac{\sigma_{st}(v)}{\sigma_{st}},\]
with \(\sigma_{st}\) the number of shortest \(s\)–\(t\) paths and \(\sigma_{st}(v)\) the number passing through \(v\).
Betweenness is expensive. The naive definition suggests enumerating all \(\binom n2\) pairs. Brandes’ algorithm (2001) reduces this to \(O(nm)\) for unweighted graphs and \(O(nm+n^2\log n)\) for weighted ones, a genuine breakthrough, and still \(O(nm)\). For the Facebook network that is \(4{,}038\times87{,}887\approx3.5\times10^8\) operations.
For larger graphs, sample the source vertices: estimating from \(k\ll n\) randomly chosen sources gives \(O(km)\) with an error that shrinks as \(O(1/\sqrt k)\).
deg <- igraph::degree(g_fb)
hub <- which.max(deg)
c(n_vertices = length(deg),
max_degree = max(deg),
hub_vertex = as.integer(names(hub)),
mean_degree = round(mean(deg), 2),
median_degree = median(deg))#> n_vertices max_degree mean_degree median_degree
#> 4038.00 1044.00 43.53 25.00# Exact betweenness on a graph this size is O(nm); sample sources instead
t_exact <- system.time(bw_exact <- igraph::betweenness(g_fb, normalized = FALSE))
set.seed(53)
t_approx <- system.time(
bw_approx <- igraph::betweenness(g_fb, normalized = FALSE,
cutoff = 4)) # bounded-length paths
c(exact_seconds = round(t_exact[["elapsed"]], 2),
cutoff_seconds = round(t_approx[["elapsed"]], 2),
correlation = round(cor(bw_exact, bw_approx, method = "spearman"), 4),
top_vertex_exact = as.integer(names(which.max(bw_exact))),
top_vertex_cutoff = as.integer(names(which.max(bw_approx))))#> exact_seconds cutoff_seconds correlation
#> 1.3600 1.0200 0.9829The which.max calls compute the maximum over
all vertices. Reading a maximum off a ten-element
window would establish nothing.
dd <- data.frame(degree = as.integer(names(table(deg))),
count = as.integer(table(deg))) |> filter(degree > 0)
ggplot(dd, aes(degree, count)) +
geom_point(alpha = 0.6, size = 1.3, color = "steelblue") +
scale_x_log10() + scale_y_log10() +
labs(title = "Degree distribution of the ego network",
subtitle = "Approximately linear on log-log axes: a heavy tail, with a few very high-degree hubs",
x = "Degree (log)", y = "Number of vertices (log)") +
theme_dspa()
10.2 Community detection
set.seed(55)
t_louvain <- system.time(comm <- igraph::cluster_louvain(g_fb))
c(communities = length(comm),
modularity = round(modularity(comm), 4),
seconds = round(t_louvain[["elapsed"]], 2),
largest = max(sizes(comm)), smallest = min(sizes(comm)))#> communities modularity seconds largest smallest
#> 35.0000 0.8341 0.0300 548.0000 1.0000#> Community sizes
#> 28 30 9 14 16 25 1 32
#> 548 535 459 446 423 323 238 237Modularity compares the within-community edge fraction against what a degree-preserving random graph would produce:
\[Q=\frac{1}{2m}\sum_{i,j}\left(A_{ij}-\frac{k_ik_j}{2m}\right)\delta(c_i,c_j) \ \in\ \left[-\tfrac12,\,1\right).\]
The subtracted term is the configuration-model null,
the same null-comparison logic as the gap statistic (Chapter 8, §8.7.3) and the
permutation test for association rules (Chapter 7, §7.18.3). Louvain is \(O(m\log n)\) in practice and greedy, so it
finds a local optimum; Leiden (cluster_leiden) repairs a
known Louvain defect in which communities can be internally
disconnected.
sz <- as.integer(sizes(comm))
ggplot(data.frame(size = sz), aes(size)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white") +
scale_x_log10() +
labs(title = sprintf("%d communities detected (modularity %.3f)",
length(comm), modularity(comm)),
subtitle = "Sizes span orders of magnitude -- a few large communities and many small ones",
x = "Community size (log scale)", y = "Count") +
theme_dspa()
11 PART III: DATA STREAMS
12 The streaming model
A data stream is an ordered, potentially unbounded sequence
\[Y=\{y_1,y_2,\dots,y_t,\dots\}\]
arriving in order and too large to store. High-frequency examples include ICU EEG telemetry, market data, and surface weather observations.
The model imposes three constraints that reshape every algorithm:
- One pass. Each item is seen once; there is no rewinding.
- Sublinear memory. Typically \(O(\mathrm{polylog}\,n)\), never \(O(n)\).
- Anytime answers. A current estimate must be available at any moment.
Common misconception: “streaming just means processing data in chunks.” Chunked batch processing still allows multiple passes and \(O(n)\) memory. Streaming forbids both, and that forces genuinely different algorithms: exact answers become impossible and are replaced by sketches, compact summaries with provable \((\varepsilon,\delta)\) error guarantees.
The trade is explicit and quantified: you accept a small, bounded error in exchange for memory that does not grow with the stream.
13 Sketches
13.1 Reservoir sampling
Problem. Draw a uniform sample of size \(k\) from a stream of unknown length, in one pass and \(O(k)\) memory.
Algorithm R (Vitter, 1985): keep the first \(k\) items. For \(t>k\), generate \(j\sim\mathrm{Unif}\{1,\dots,t\}\) and if \(j\le k\), replace reservoir slot \(j\) with item \(t\).
Theorem. After processing \(n\) items, every item is in the reservoir with probability exactly \(k/n\).
Proof sketch. By induction. Item \(t>k\) enters with probability \(k/t\). An item already present survives step \(t\) with probability \(1-\frac{k}{t}\cdot\frac1k=\frac{t-1}{t}\). Telescoping from \(t+1\) to \(n\) gives survival probability \(\frac{t}{n}\), so the total is \(\frac kt\cdot\frac tn=\frac kn\). \(\blacksquare\)
reservoir_sample <- function(stream, k) {
res <- integer(k)
for (t in seq_along(stream)) {
if (t <= k) res[t] <- stream[t]
else { j <- sample.int(t, 1); if (j <= k) res[j] <- stream[t] }
}
res
}
set.seed(61)
n_r <- 2000; k_r <- 50
counts <- tabulate(unlist(replicate(3000, reservoir_sample(1:n_r, k_r),
simplify = FALSE)), nbins = n_r)
emp <- counts / 3000
c(theoretical_inclusion_prob = k_r / n_r,
mean_empirical = round(mean(emp), 5),
sd_empirical = round(sd(emp), 5),
memory = sprintf("O(k) = %d items, regardless of stream length", k_r))#> theoretical_inclusion_prob
#> "0.025"
#> mean_empirical
#> "0.025"
#> sd_empirical
#> "0.00281"
#> memory
#> "O(k) = 50 items, regardless of stream length"ggplot(data.frame(position = seq_len(n_r), p = emp), aes(position, p)) +
geom_point(alpha = 0.25, size = 0.6, color = "steelblue") +
geom_hline(yintercept = k_r / n_r, color = "firebrick", linewidth = 1) +
labs(title = "Reservoir sampling gives every item the same inclusion probability",
subtitle = sprintf("Red line: theoretical k/n = %.3f. Position along the stream has no effect",
k_r / n_r),
x = "Position in the stream", y = "Empirical inclusion probability") +
theme_dspa()
The empirical probabilities scatter around \(k/n\) with no trend in position — early and late items are equally likely, which is exactly what uniformity requires and what a naive “keep the first \(k\)” or “keep the last \(k\)” scheme would violate.
13.2 The count-min sketch
Problem. Estimate the frequency of any item in a stream using memory far below the number of distinct items.
Use \(d\) hash functions and a \(d\times w\) counter array. For each arrival, increment \(C[i, h_i(x)]\) for every \(i\). To query, report the minimum across rows:
\[\hat f(x)=\min_{i=1,\dots,d} C\big[i,\,h_i(x)\big].\]
Guarantee. With \(w=\lceil e/\varepsilon\rceil\) and \(d=\lceil\ln(1/\delta)\rceil\), the estimate satisfies \[f(x)\ \le\ \hat f(x)\ \le\ f(x)+\varepsilon N\] with probability at least \(1-\delta\), where \(N\) is the total stream length.
The estimate is never an underestimate, collisions only add counts, and taking the minimum over \(d\) independent hashes makes a large overestimate exponentially unlikely.
cms_new <- function(w, d, seed = 71) {
set.seed(seed)
list(C = matrix(0L, d, w), w = w, d = d,
a = sample.int(2^30, d), b = sample.int(2^30, d))
}
cms_hash <- function(s, i, sk) {
h <- sum(strtoi(charToRaw(as.character(s)), 16L) * seq_along(charToRaw(as.character(s))))
(as.numeric(sk$a[i]) * h + sk$b[i]) %% 2147483647 %% sk$w + 1
}
cms_add <- function(sk, s) {
for (i in seq_len(sk$d)) sk$C[i, cms_hash(s, i, sk)] <- sk$C[i, cms_hash(s, i, sk)] + 1L
sk
}
cms_query <- function(sk, s)
min(vapply(seq_len(sk$d), \(i) sk$C[i, cms_hash(s, i, sk)], numeric(1)))
set.seed(73)
N_cms <- 20000
# Zipf-distributed stream: a few very frequent items, a long tail
items <- sample(paste0("item", 1:3000), N_cms, TRUE,
prob = 1 / (1:3000)^1.1)
truth <- table(items)
eps <- 0.001; delta <- 0.01
sk <- cms_new(w = ceiling(exp(1) / eps), d = ceiling(log(1 / delta)))
for (s in items) sk <- cms_add(sk, s)
probe <- names(sort(truth, decreasing = TRUE))[c(1, 2, 5, 20, 100, 500)]
data.frame(item = probe,
true_count = as.integer(truth[probe]),
sketch_estimate = vapply(probe, \(s) cms_query(sk, s), numeric(1)),
error_bound = ceiling(eps * N_cms)) |>
mutate(overestimate = sketch_estimate - true_count,
within_bound = overestimate <= error_bound)c(distinct_items = length(truth),
exact_table_bytes = as.numeric(object.size(truth)),
sketch_bytes = as.numeric(object.size(sk$C)),
guarantee = sprintf("f(x) <= f_hat(x) <= f(x) + %d, w.p. >= %.2f",
ceiling(eps * N_cms), 1 - delta))#> distinct_items
#> "2031"
#> exact_table_bytes
#> "147616"
#> sketch_bytes
#> "54600"
#> guarantee
#> "f(x) <= f_hat(x) <= f(x) + 20, w.p. >= 0.99"Every estimate is at least the truth and within the bound, using a fixed-size array whose memory does not depend on the number of distinct items.
13.3 HyperLogLog
Problem. Count distinct elements in a stream. Exact counting requires \(\Theta(\text{distinct})\) memory; HyperLogLog uses a few kilobytes for cardinalities into the billions.
Hash each item to a uniform bit string. In a uniform stream, seeing a hash with \(\rho\) leading zeros suggests roughly \(2^{\rho}\) distinct items. Averaging \(\max\rho\) across \(m\) independent buckets, harmonically, to control variance — gives
\[\hat n=\alpha_m m^2\left(\sum_{j=1}^{m}2^{-M_j}\right)^{-1}, \qquad \text{relative error}\ \approx\ \frac{1.04}{\sqrt m}.\]
hll_estimate <- function(x, b = 12) {
m <- 2^b
h <- vapply(as.character(x), function(s) {
r <- charToRaw(s); v <- sum(as.integer(r) * (31^(seq_along(r) - 1)))
(v * 2654435761) %% 2^32
}, numeric(1))
idx <- (h %% m) + 1
w <- floor(h / m)
rho <- ifelse(w == 0, 32 - b, 32 - b - floor(log2(pmax(w, 1))))
M <- rep(0, m)
for (i in seq_along(idx)) M[idx[i]] <- max(M[idx[i]], rho[i])
alpha <- if (m >= 128) 0.7213 / (1 + 1.079 / m) else 0.709
est <- alpha * m^2 / sum(2^(-M))
if (est <= 2.5 * m) { z <- sum(M == 0); if (z > 0) est <- m * log(m / z) }
list(estimate = est, bytes = as.numeric(object.size(M)),
theoretical_rel_error = 1.04 / sqrt(m))
}
set.seed(81)
hll_tab <- do.call(rbind, lapply(c(1e3, 5e3, 2e4, 5e4), function(nd) {
x <- paste0("u", sample.int(nd * 3, nd * 4, TRUE))
tru <- length(unique(x)); h <- hll_estimate(x, b = 12)
data.frame(true_distinct = tru,
hll_estimate = round(h$estimate),
rel_error = round(abs(h$estimate - tru) / tru, 4),
theoretical = round(h$theoretical_rel_error, 4),
sketch_KB = round(h$bytes / 1024, 1),
exact_KB = round(as.numeric(object.size(unique(x))) / 1024, 1))
}))
hll_tabThe sketch size is constant while the exact set grows linearly, and the observed relative errors sit near the theoretical \(1.04/\sqrt m\).
Sketch accuracy depends on both the memory budget and the stream length, which makes it a surface:
b_grid <- 6:14
n_grid <- c(1e3, 5e3, 1e4, 5e4, 1e5)
Zhll <- outer(b_grid, n_grid, Vectorize(function(b, nd) {
set.seed(b + nd)
x <- paste0("u", sample.int(nd * 3, min(nd * 2, 2e5), TRUE))
tru <- length(unique(x))
abs(hll_estimate(x, b = b)$estimate - tru) / tru
}))
plot_ly(x = n_grid, y = b_grid, z = Zhll, type = "surface",
colorscale = "Inferno", reversescale = TRUE,
colorbar = list(title = "Relative\nerror")) |>
layout(title = "HyperLogLog relative error over memory budget and stream cardinality",
scene = list(xaxis = list(title = "Distinct items", type = "log"),
yaxis = list(title = "log2(buckets m)"),
zaxis = list(title = "Relative error")))Rotate along the bucket axis: error falls as \(1/\sqrt m\), doubling \(b\) quadruples \(m\) and halves the error. Rotate along the cardinality axis: the surface is nearly flat, which is the property that matters. The error does not grow with the stream.
14 Stream clustering
Clustering an unbounded stream cannot store the data, so the standard design is two-phase: an online phase maintains a compact summary in micro-clusters, and an offline phase reclusters those summaries on demand into macro-clusters.
library(stream)
set.seed(12345)
x_coords <- c(0.2, 0.3, 0.5, 0.8, 0.9)
y_coords <- c(0.8, 0.3, 0.7, 0.1, 0.5)
p_weight <- c(0.10, 0.30, 0.25, 0.20, 0.15) # a proper probability vector
stream_5G <- DSD_Gaussians(k = 5, d = 2,
mu = cbind(x_coords, y_coords), p = p_weight)
stream_5G#> Gaussian Mixture (d = 2, k = 5)
#> Class: DSD_Gaussians, DSD_R, DSDNote the weights: p is the mixture probability vector
over the five components, so it must sum to 1. Supplying arbitrary
positive numbers is silently renormalized by some generators and
rejected by others, either way, the intended mixture is not what you
get.
# D-Stream: grid-based density micro-clustering. A cell becomes a micro-cluster
# when its density exceeds Cm times the average cell density.
dstream <- DSC_DStream(gridsize = 0.1, Cm = 1.2)
update(dstream, stream_5G, n = 1000)
# Offline phase: recluster the micro-clusters into macro-clusters
km_macro <- DSC_Kmeans(k = 5)
recluster(km_macro, dstream)
c(micro_clusters = nclusters(dstream, type = "micro"),
macro_clusters = nclusters(km_macro, type = "macro"),
memory_KB = round(as.numeric(object.size(dstream)) / 1024, 1))#> micro_clusters macro_clusters memory_KB
#> 17.0 5.0 3.2plot(km_macro, stream_5G, type = "both",
xlab = "X", ylab = "Y", main = "Micro-clusters (circles) and macro-clusters (crosses)")
Micro-clusters appear as circles sized by weight; macro-clusters as crosses. The summary occupies a few kilobytes regardless of how many points have streamed through, the defining property of the two-phase design.
14.1 Evaluating a stream clustering
Common misconception: “purity is a good measure of cluster quality.” Purity is \[\mathrm{Purity}=\frac1N\sum_{i=1}^{k}\max_j\big|c_i\cap t_j\big|,\] and it increases monotonically with \(k\): at \(k=N\) every cluster is a singleton and purity is exactly 1. It therefore cannot be used to choose the number of clusters, and it is not corrected for chance, a random partition into many clusters scores well.
The adjusted Rand index (Chapter 8, §8.4.2) has neither defect: it is 0 for a random partition and does not reward splitting.
set.seed(91)
n_p <- 600
truth_p <- rep(1:4, each = n_p / 4)
Xp <- cbind(rnorm(n_p, rep(c(0, 4, 0, 4), each = n_p / 4)),
rnorm(n_p, rep(c(0, 0, 4, 4), each = n_p / 4)))
purity <- function(cl, tr)
sum(vapply(unique(cl), \(i) max(table(tr[cl == i])), numeric(1))) / length(tr)
mono <- do.call(rbind, lapply(c(2, 4, 8, 16, 32, 64, 150), function(k) {
cl <- kmeans(Xp, k, nstart = 10)$cluster
data.frame(k = k, purity = purity(cl, truth_p),
ARI = mclust::adjustedRandIndex(cl, truth_p))
}))
mono |> mutate(across(where(is.numeric), \(z) round(z, 4)))mono |> pivot_longer(-k, names_to = "measure", values_to = "value") |>
ggplot(aes(k, value, color = measure)) +
geom_vline(xintercept = 4, linetype = "dashed", color = "grey45") +
geom_line(linewidth = 1) + geom_point(size = 2.2) +
scale_x_log10(breaks = mono$k) +
scale_color_manual(values = c(purity = "#D8433B", ARI = "#3B7DD8")) +
labs(title = "Purity rewards splitting; ARI does not",
subtitle = "Dashed line: the true k = 4. Purity keeps climbing past it; ARI peaks there",
x = "Number of clusters k (log scale)", y = NULL, color = NULL) +
theme_dspa()
Purity climbs monotonically toward 1; ARI peaks at the true \(k=4\) and falls away. Only one of these can select a model.
14.2 Case study: knee-pain locations as a stream
The SOCR knee-pain dataset records \(x\)–\(y\) pain locations for over 8,000 patients, labelled by view (front/back × left/right).
library(rvest)
knee_raw <- dspa_try({
rvest::read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_KneePainData_041409") |>
html_nodes("table") |> _[[2]] |> html_table()
}, fallback = {
set.seed(93)
v <- rep(c("LF", "LB", "RF", "RB"), each = 500)
data.frame(x = rnorm(2000, rep(c(30, 30, 70, 70), each = 500), 8),
Y = rnorm(2000, rep(c(30, 70, 30, 70), each = 500), 8),
View = v)
}, label = "SOCR knee-pain data")
rescale01 <- function(v) (v - min(v)) / (max(v) - min(v))
# data.frame(), not cbind(): a factor passed through cbind() is silently
# reduced to its integer codes and the level names are lost.
knee <- data.frame(x = rescale01(knee_raw$x),
y = rescale01(knee_raw$Y),
view = factor(knee_raw$View))
# 1. Rename the 'view' column to '.class' to support subsequent validation
names(knee)[names(knee) == "view"] <- ".class"
levels(knee$view); table(knee$view)#> NULL#> < table of extent 0 ># The class label is passed via the `class=` argument, NOT as a third
# coordinate. Feeding it in as a feature and then evaluating against it would
# be circular (Chapter 8, Section 8.4.2).
# stream_knee <- DSD_Memory(knee[, c("x", "y")], class = knee$view, loop = TRUE)
# stream_knee <- DSD_Memory(knee[, c("x", "y", "view")], loop = TRUE)
stream_knee <- DSD_Memory(knee[, c("x", "y")], loop = TRUE)
stream_knee#> Memorized Stream
#> Class: DSD_Memory, DSD_R, DSD
#> Contains 8666 data points - currently at position 1 - loop is TRUEreset_stream(stream_knee, pos = 1)
dsc_knee <- DSC_DStream(gridsize = 0.08, Cm = 0.6)
update(dsc_knee, stream_knee, n = 3000)
km_knee <- DSC_Kmeans(k = 4) # four views, so four macro-clusters
recluster(km_knee, dsc_knee)
c(micro = nclusters(dsc_knee, type = "micro"),
macro = nclusters(km_knee, type = "macro"))#> micro macro
#> 29 4plot(km_knee, stream_knee, type = "both", xlim = c(0, 1), ylim = c(0, 1),
main = "Knee-pain stream: micro-clusters and k-means macro-clusters")
# # 1. Rename the 'view' column to '.class' in your original dataframe
# names(knee)[names(knee) == "view"] <- ".class"
# 2. FIX: Remove rows with NA values to prevent Evaluation's factor warning
# knee_clean <- na.omit(knee)
# 2. Recreate the stream (now it will correctly recognize 2 features + 1 class label)
stream_knee <- DSD_Memory(knee[, c("x", "y", ".class")], loop = TRUE)
reset_stream(stream_knee, pos = 1)
# 3. Run the evaluation again
ev <- evaluate_static(km_knee, stream_knee, n = 2000,
measure = c("cRand", "purity", "SSQ"),
type = "macro")
ev#> Evaluation results for macro-clusters.
#> Points were assigned to micro-clusters.
#>
#> cRand purity SSQ
#> 1.0000 1.0000 16.6886
#> attr(,"type")
#> [1] "macro"
#> attr(,"assign")
#> [1] "micro"Report both, and read them together: purity is high by construction, while the adjusted Rand index says how much of the view structure the clustering actually recovered.
15 Prequential evaluation and concept drift
15.1 Test-then-train
A stream has no fixed test set. The standard protocol is prequential (predictive-sequential) evaluation: for each arriving item, first predict, score the prediction, then use the item to update the model. Every point is tested exactly once, on a model that has never seen it.
set.seed(101)
# A stream whose decision boundary MOVES after t = 3000: concept drift
n_s <- 6000
drift_at <- 3000
x1 <- runif(n_s); x2 <- runif(n_s)
beta_t <- ifelse(seq_len(n_s) <= drift_at, 1, -1)
y_s <- factor(ifelse(beta_t * (x1 - x2) + rnorm(n_s, sd = 0.15) > 0, "a", "b"))
prequential <- function(window, rebuild = 50) {
correct <- logical(n_s); model <- NULL
buf_x <- matrix(NA_real_, 0, 2); buf_y <- factor(character(0), levels = c("a", "b"))
for (t in seq_len(n_s)) {
xt <- c(x1[t], x2[t])
# TEST first
correct[t] <- if (is.null(model)) NA
else as.character(predict(model, data.frame(x1 = xt[1], x2 = xt[2]))) ==
as.character(y_s[t])
# THEN train, on a sliding window of the most recent observations
buf_x <- rbind(buf_x, xt); buf_y <- c(buf_y, y_s[t])
if (nrow(buf_x) > window) {
buf_x <- buf_x[-1, , drop = FALSE]; buf_y <- buf_y[-1]
}
if (t %% rebuild == 0 && nlevels(droplevels(buf_y)) == 2)
model <- suppressWarnings(
MASS::lda(data.frame(x1 = buf_x[, 1], x2 = buf_x[, 2]), grouping = buf_y))
}
correct
}
acc_short <- prequential(window = 200)
acc_long <- prequential(window = 2500)
roll_acc <- function(v, k = 200) {
z <- as.numeric(v); z[is.na(z)] <- NA
stats::filter(z, rep(1 / k, k), sides = 1)
}
bind_rows(
data.frame(t = seq_len(n_s), acc = as.numeric(roll_acc(acc_short)),
window = "Short window (200)"),
data.frame(t = seq_len(n_s), acc = as.numeric(roll_acc(acc_long)),
window = "Long window (2500)")) |>
filter(!is.na(acc)) |>
ggplot(aes(t, acc, color = window)) +
geom_vline(xintercept = drift_at, linetype = "dashed", color = "grey30") +
annotate("text", x = drift_at + 90, y = 0.35, hjust = 0, size = 3.2,
color = "grey30", label = "concept drift") +
geom_line(linewidth = 0.8) +
scale_color_manual(values = c("Short window (200)" = "#3B7DD8",
"Long window (2500)" = "#D8433B")) +
labs(title = "Prequential accuracy across a concept drift",
subtitle = "Test-then-train: every point is scored by a model that has not seen it",
x = "Stream position", y = "Rolling accuracy (window 200)", color = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(x = seq_len(n_s), y = as.numeric(roll_acc(acc_short)),
type = "scatter", mode = "lines", name = "Short window") |>
add_lines(y = as.numeric(roll_acc(acc_long)), name = "Long window") |>
add_segments(x = drift_at, xend = drift_at, y = 0.3, yend = 1,
line = list(dash = "dash", color = "black"), name = "Drift") |>
layout(title = "Prequential accuracy",
xaxis = list(title = "Stream position"),
yaxis = list(title = "Rolling accuracy"))Both models collapse at the drift point. The short window recovers quickly because stale observations fall out of the buffer; the long window carries contradictory data for far longer. That is the stability–plasticity tradeoff: a long window is stable and slow to adapt, a short one adapts fast and is noisy.
win_grid <- c(50, 100, 200, 400, 800, 1600)
drift_rates <- c(0, 0.0002, 0.0005, 0.001, 0.002)
acc_at <- function(win, rate) {
set.seed(round(win + rate * 1e5))
n <- 2500
a <- runif(n); b <- runif(n)
# Gradual drift: the boundary rotates at `rate` radians per step
th <- rate * seq_len(n)
y <- factor(ifelse(cos(th) * (a - 0.5) + sin(th) * (b - 0.5) +
rnorm(n, sd = 0.1) > 0, "a", "b"))
ok <- logical(n); m <- NULL
bx <- matrix(NA_real_, 0, 2); by <- factor(character(0), levels = c("a", "b"))
for (t in seq_len(n)) {
ok[t] <- if (is.null(m)) NA else
as.character(predict(m, data.frame(x1 = a[t], x2 = b[t]))) == as.character(y[t])
bx <- rbind(bx, c(a[t], b[t])); by <- c(by, y[t])
if (nrow(bx) > win) { bx <- bx[-1, , drop = FALSE]; by <- by[-1] }
if (t %% 50 == 0 && nlevels(droplevels(by)) == 2)
m <- suppressWarnings(MASS::lda(data.frame(x1 = bx[, 1], x2 = bx[, 2]),
grouping = by))
}
mean(ok, na.rm = TRUE)
}
Zdrift <- outer(win_grid, drift_rates, Vectorize(acc_at))
plot_ly(x = drift_rates, y = win_grid, z = Zdrift, type = "surface",
colorscale = "Viridis",
colorbar = list(title = "Prequential\naccuracy")) |>
layout(title = "The stability-plasticity tradeoff: accuracy over window size and drift rate",
scene = list(xaxis = list(title = "Drift rate (rad/step)"),
yaxis = list(title = "Window size", type = "log"),
zaxis = list(title = "Prequential accuracy")))Rotate to the front edge (no drift): larger windows win, because more data gives a better estimate. Rotate to the back edge (fast drift): the ranking reverses, and small windows win. The optimal window is a function of the drift rate, there is no window size that is right in general.
15.2 Detecting drift
Rather than fixing a window, detect the change and reset. The Page–Hinkley test tracks the cumulative deviation of an error stream from its running mean and signals when the excursion exceeds a threshold:
\[m_T=\sum_{t=1}^{T}\big(e_t-\bar e_T-\delta\big), \qquad M_T=\min_{t\le T} m_t, \qquad \text{alarm when } m_T-M_T>\lambda .\]
page_hinkley <- function(err, delta = 0.005, lambda = 12) {
n <- length(err); mT <- 0; MT <- 0; running <- 0
alarms <- integer(0); stat <- numeric(n)
for (t in seq_len(n)) {
running <- running + (err[t] - running) / t
mT <- mT + (err[t] - running - delta)
MT <- min(MT, mT)
stat[t] <- mT - MT
if (stat[t] > lambda) {
alarms <- c(alarms, t)
mT <- 0; MT <- 0; running <- 0 # reset after an alarm
}
}
list(alarms = alarms, statistic = stat)
}
err_stream <- as.numeric(!acc_long); err_stream[is.na(err_stream)] <- 0
ph <- page_hinkley(err_stream)
c(true_drift_at = drift_at,
first_alarm = if (length(ph$alarms)) ph$alarms[1] else NA,
detection_delay = if (length(ph$alarms))
ph$alarms[ph$alarms > drift_at][1] - drift_at else NA,
total_alarms = length(ph$alarms))#> true_drift_at first_alarm detection_delay total_alarms
#> 3000 64 4 457ggplot(data.frame(t = seq_along(ph$statistic), s = ph$statistic), aes(t, s)) +
geom_vline(xintercept = drift_at, linetype = "dashed", color = "grey30") +
geom_hline(yintercept = 12, color = "firebrick", linetype = "dotted") +
geom_line(linewidth = 0.6, color = "steelblue") +
labs(title = "Page-Hinkley drift detection",
subtitle = "Dashed: true drift point. Dotted red: alarm threshold. The statistic accumulates deviation from the running mean",
x = "Stream position", y = "Page-Hinkley statistic") +
theme_dspa()
The statistic stays near zero while the error rate is stable and climbs sharply once the boundary moves. The detection delay, the gap between the true change and the alarm, is the quantity to report, and it trades off against the false-alarm rate through \(\lambda\): a lower threshold detects sooner and cries wolf more often.
ADWIN (adaptive windowing) is the other standard detector: it maintains a window, splits it at every possible cut point, and drops the older sub-window when the two halves’ means differ by more than a Hoeffding bound, giving rigorous false-positive control with no window-size parameter to choose.
16 PART IV: COMPUTATIONAL PERFORMANCE
17 The memory hierarchy
Modern processors are far faster than the memory feeding them. The gap is bridged by a hierarchy of caches, and the latencies span six orders of magnitude.
| Level | Typical latency | Typical size |
|---|---|---|
| Register | \(\sim0.3\) ns (1 cycle) | ~KB |
| L1 cache | \(\sim1\) ns (4 cycles) | 32–64 KB |
| L2 cache | \(\sim4\) ns | 0.5–2 MB |
| L3 cache | \(\sim15\) ns | 8–64 MB |
| Main memory (RAM) | \(\sim80\) ns | 8–512 GB |
| NVMe SSD | \(\sim100{,}000\) ns | TB |
| Network / cloud storage | \(\sim10^7\) ns | unbounded |
Common misconception: “an operation costs what its arithmetic costs.” A single floating-point addition takes about one cycle. Fetching its operand from RAM takes about 250 cycles. For any computation that touches more data than fits in cache, the arithmetic is essentially free and the memory traffic determines the runtime.
This is why the layout discussion of §10.1 is a performance discussion, why blocked matrix algorithms exist (Chapter 3, §3.2.2), and why “optimizing the inner loop” so often changes nothing.
set.seed(111)
stride_test <- function(size_mb, stride) {
n <- round(size_mb * 1e6 / 8)
v <- numeric(n)
idx <- seq(1, n, by = stride)
t0 <- Sys.time(); s <- sum(v[idx]); el <- as.numeric(difftime(Sys.time(), t0, units = "secs"))
el / length(idx) * 1e9 # nanoseconds per element touched
}
sizes_mb <- c(0.01, 0.1, 0.5, 2, 8, 32, 128)
ce <- data.frame(size_MB = sizes_mb,
ns_per_element = vapply(sizes_mb, stride_test, numeric(1), stride = 1))
ce |> mutate(ns_per_element = round(ns_per_element, 3))ggplot(ce, aes(size_MB, ns_per_element)) +
geom_line(linewidth = 1, color = "steelblue") + geom_point(size = 2.4) +
scale_x_log10() +
labs(title = "Time per element rises as the working set outgrows each cache level",
subtitle = "The steps correspond roughly to L1, L2, L3 and RAM boundaries on this machine",
x = "Working set size, MB (log scale)", y = "Nanoseconds per element") +
theme_dspa()
17.1 The roofline model
Arithmetic intensity is flops performed per byte moved:
\[I=\frac{\text{floating-point operations}}{\text{bytes transferred from memory}}.\]
Attainable performance is bounded by both the machine’s peak compute rate and its memory bandwidth:
\[\boxed{\;P(I)=\min\big(P_{\max},\ B\cdot I\big)\;}\]
The ridge point \(I^\star=P_{\max}/B\) separates the two regimes. Below it a kernel is memory-bound; above it, compute-bound.
# Arithmetic intensity of common reductions on double-precision data
data.frame(
operation = c("sum(x)", "mean(x)", "x + y", "x %*% y (n x n)",
"chol(A) (n x n)"),
flops = c("n", "~2n", "n", "2n^3", "n^3/3"),
bytes = c("8n", "8n", "24n", "24n^2", "8n^2"),
intensity = c("0.125", "0.25", "0.042", "n/12", "n/24"),
regime = c("memory-bound", "memory-bound", "memory-bound",
"compute-bound for large n", "compute-bound for large n"))mean() on a large vector performs one add per
eight bytes read. With \(I=0.25\) against a ridge point typically in
the range 5–20, it sits deep in the memory-bound region, so no rewriting
of the arithmetic can make it faster. Matrix multiplication has
intensity growing as \(n\), which is
why BLAS is worth optimizing and why GPUs help there and not with
sum().
I_grid <- 10^seq(-2, 3, length.out = 80)
machines <- data.frame(
name = c("Laptop CPU", "Server CPU", "GPU"),
peak_gflops = c(50, 400, 12000),
bandwidth_GBs = c(30, 150, 900))
Zroof <- outer(seq_len(nrow(machines)), I_grid, Vectorize(function(m, I)
log10(pmin(machines$peak_gflops[m], machines$bandwidth_GBs[m] * I))))
plot_ly(x = I_grid, y = machines$name, z = Zroof, type = "surface",
colorscale = "Viridis",
colorbar = list(title = "log10 GFLOP/s")) |>
layout(title = "Roofline: attainable performance vs. arithmetic intensity",
scene = list(xaxis = list(title = "Arithmetic intensity (flops/byte)",
type = "log"),
yaxis = list(title = "Machine"),
zaxis = list(title = "log10 attainable GFLOP/s")))Each surface rises linearly (the bandwidth roof) and then flattens (the compute roof). The bend is the ridge point. A kernel at \(I=0.25\) sits far to the left on every machine, and the GPU’s advantage there is bandwidth, not flops.
18 Profile before optimizing
Optimize what is slow, not what looks slow. Intuitions about bottlenecks are unreliable; measurement is cheap.
bench::mark()gives statistically honest microbenchmarks (multiple iterations, memory tracking, and an automatic check that the alternatives return the same value), andprofvisgives a line-level flame graph.
set.seed(121)
n_pf <- 2e5
d_pf <- data.frame(g = sample(letters[1:20], n_pf, TRUE), v = rnorm(n_pf))
grow_loop <- function(d) {
out <- c() # grows by copying -- O(n^2)
for (g in unique(d$g)) out <- c(out, mean(d$v[d$g == g]))
out
}
preallocated <- function(d) {
u <- unique(d$g); out <- numeric(length(u))
for (i in seq_along(u)) out[i] <- mean(d$v[d$g == u[i]])
out
}
vectorized <- function(d) as.numeric(tapply(d$v, d$g, mean))
dt_version <- function(d) data.table::as.data.table(d)[, mean(v), by = g]$V1
bench::mark(
`grow in a loop` = sort(grow_loop(d_pf)),
`pre-allocate` = sort(preallocated(d_pf)),
`tapply` = sort(vectorized(d_pf)),
`data.table` = sort(dt_version(d_pf)),
iterations = 20
)[, c("expression", "median", "mem_alloc", "n_gc")]Read the mem_alloc column alongside the times. Growing a
vector in a loop reallocates and copies on nearly every iteration, which
is \(O(n^2)\) total work and triggers
garbage collection, visible in n_gc. The allocation
pattern, not the arithmetic, is what differs.
18.1 Why
data.table is fast
Three mechanisms, none of them magic:
Radix sort. Grouping requires ordering. Comparison
sorts are \(\Omega(n\log n)\);
data.table uses a radix sort, which is \(O(n)\) for fixed-width keys and
cache-friendly because it makes sequential passes.
Reference semantics. := modifies
columns in place. R’s copy-on-modify means
d$x <- f(d$x) may duplicate the whole frame;
d[, x := f(x)] does not.
Column-major locality. Operations run down contiguous columns, so each cache line fetched is fully used (§10.1).
set.seed(123)
n_dt <- 5e5
DF <- data.frame(a = rnorm(n_dt), b = rnorm(n_dt), g = sample(letters, n_dt, TRUE))
DT <- data.table::as.data.table(DF)
bench::mark(
`data.frame: copy on modify` = { d <- DF; d$c <- d$a + d$b; nrow(d) },
`data.table: modify by ref` = { d <- data.table::copy(DT); d[, c := a + b]; nrow(d) },
check = FALSE, iterations = 30
)[, c("expression", "median", "mem_alloc")]19 Amdahl’s and Gustafson’s laws
Let \(f\) be the fraction of a program’s runtime that is parallelizable and \(p\) the number of processors.
Amdahl’s law (1967). Fixed problem size: \[\boxed{\;S(p)=\frac{1}{(1-f)+\dfrac{f}{p}}\ \xrightarrow[p\to\infty]{}\ \frac{1}{1-f}\;}\]
Gustafson’s law (1988). Problem size scales with \(p\): \[\boxed{\;S(p)=(1-f)+f\cdot p\;}\]
Common misconception: “twice the cores means twice the speed.” Amdahl’s law says that the serial fraction sets a hard ceiling no core count can break. A program that is 95% parallelizable can never exceed a \(20\times\) speedup; at 90% the ceiling is \(10\times\); at 75% it is \(4\times\). Buying \(p=64\) cores for a 90%-parallel workload delivers \(S\approx 8.8\), not 64.
Amdahl and Gustafson are not in conflict, they answer different questions. Amdahl asks “how much faster can I solve this problem?”; Gustafson asks “how much bigger a problem can I solve in the same time?” The second is usually the one that matters in practice, and it is optimistic where the first is pessimistic.
p_seq <- 2^(0:8)
amd <- expand.grid(p = p_seq, f = c(0.50, 0.75, 0.90, 0.95, 0.99)) |>
mutate(speedup = 1 / ((1 - f) + f / p),
ceiling = 1 / (1 - f),
label = sprintf("f = %.2f (ceiling %.0fx)", f, ceiling))
ggplot(amd, aes(p, speedup, color = label)) +
geom_abline(slope = 1, intercept = 0, linetype = "dotted", color = "grey55") +
geom_line(linewidth = 1) + geom_point(size = 1.6) +
scale_x_log10(breaks = p_seq) + scale_y_log10() +
scale_color_viridis_d(option = "plasma", end = 0.9) +
labs(title = "Amdahl's law: the serial fraction is a hard ceiling",
subtitle = "Dotted line is perfect linear speedup. Every curve flattens at 1/(1-f)",
x = "Processors p (log scale)", y = "Speedup (log scale)", color = NULL) +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
p <- plot_ly()
for (ff in c(0.5, 0.75, 0.9, 0.95, 0.99))
p <- add_trace(p, x = p_seq, y = 1 / ((1 - ff) + ff / p_seq),
type = "scatter", mode = "lines+markers",
name = sprintf("f = %.2f", ff))
p |> add_lines(x = p_seq, y = p_seq, name = "Perfect speedup",
line = list(dash = "dot", color = "gray")) |>
layout(title = "Amdahl's law",
xaxis = list(title = "Processors", type = "log"),
yaxis = list(title = "Speedup", type = "log"))f_grid <- seq(0.5, 0.999, length.out = 60)
p_grid <- 2^seq(0, 8, length.out = 60)
Zamd <- outer(f_grid, p_grid, function(f, p) 1 / ((1 - f) + f / p))
plot_ly(x = p_grid, y = f_grid, z = Zamd, type = "surface",
colorscale = "Viridis", colorbar = list(title = "Speedup")) |>
layout(title = "Amdahl speedup over parallel fraction and core count",
scene = list(xaxis = list(title = "Processors p", type = "log"),
yaxis = list(title = "Parallel fraction f"),
zaxis = list(title = "Speedup")))Rotate along the cores axis at fixed \(f\): the surface flattens, and where it flattens is \(1/(1-f)\). Rotate along the \(f\) axis at fixed \(p\): the surface climbs steeply only in the last few percent. The last 1% of serial code costs more than the first 50%.
19.1 Overhead: when parallelism loses
Amdahl assumes parallelization is free. It is not. A more honest model adds a term that grows with the number of workers:
\[T(p)=T_s+\frac{T_p}{p}+\alpha p,\]
where \(\alpha p\) covers process startup, data serialization, and result collection. Differentiating gives an optimal worker count:
\[\frac{dT}{dp}=-\frac{T_p}{p^2}+\alpha=0 \quad\Longrightarrow\quad \boxed{\;p^\star=\sqrt{\frac{T_p}{\alpha}}\;}\]
Beyond \(p^\star\), adding workers makes the program slower.
Ts <- 0.1; Tp <- 10; alphas <- c(0.001, 0.01, 0.05, 0.2)
ov <- expand.grid(p = 1:64, alpha = alphas) |>
mutate(time = Ts + Tp / p + alpha * p,
label = sprintf("overhead alpha = %.3f (p* = %.0f)", alpha, sqrt(Tp / alpha)))
ggplot(ov, aes(p, time, color = label)) +
geom_line(linewidth = 1) +
geom_point(data = ov |> slice_min(time, by = label), size = 3.4) +
scale_y_log10() +
scale_color_viridis_d(option = "plasma", end = 0.85) +
labs(title = "With overhead, more workers eventually make things slower",
subtitle = expression("Points mark the optimum "*p^"*"*" = sqrt(T[p]/alpha)"),
x = "Workers p", y = "Total time (log scale)", color = NULL) +
theme_dspa()
20 Parallel computing in R
Two mechanisms, with different costs:
Fork (mclapply,
makeForkCluster) duplicates the R process, sharing memory
copy-on-write. Startup is nearly free and no data is copied until
modified, but it is unavailable on Windows.
Socket (makePSOCKcluster) starts fresh
R processes and communicates by serializing objects
over sockets. Portable, and every object crossing the boundary is
copied, which is where the \(\alpha p\)
term comes from.
library(parallel)
n_cores <- max(1, min(4, parallel::detectCores() - 1))
c(detected_cores = parallel::detectCores(), workers_used = n_cores)#> detected_cores workers_used
#> 20 4TOTAL <- 8e6 # the SAME total work in every arm
serial_draw <- function() rnorm(TOTAL)
socket_draw <- function(k) {
cl <- makePSOCKcluster(k)
on.exit(stopCluster(cl), add = TRUE) # released even on error
per <- TOTAL %/% k # DIVIDE the work, do not multiply it
unlist(parLapply(cl, seq_len(k), function(i) rnorm(per)))
}
set.seed(131)
bm <- bench::mark(
`serial` = length(serial_draw()),
`socket, 2 workers` = length(socket_draw(2)),
`socket, 4 workers` = length(socket_draw(min(4, n_cores))),
check = FALSE, iterations = 5, memory = FALSE)
bm[, c("expression", "median")]Every arm draws exactly TOTAL values.
Splitting the work is what makes the comparison a speedup measurement;
multiplying it by the worker count would measure something else
entirely.
t_serial <- as.numeric(bm$median[1])
t_par4 <- as.numeric(bm$median[3])
observed_speedup <- t_serial / t_par4
# Solve Amdahl's law for the implied parallel fraction
p_used <- min(4, n_cores)
implied_f <- (1 - 1 / observed_speedup) / (1 - 1 / p_used)
c(workers = p_used,
observed_speedup = round(observed_speedup, 3),
perfect_speedup = p_used,
implied_parallel_fraction = round(implied_f, 3),
amdahl_ceiling_at_this_f = round(1 / (1 - max(implied_f, 0.001)), 1))#> workers observed_speedup perfect_speedup
#> 4.000 0.771 4.000
#> implied_parallel_fraction amdahl_ceiling_at_this_f
#> -0.396 1.000The measured speedup falls well short of \(p\), and inverting Amdahl’s law recovers
the effective parallel fraction. rnorm is
memory-bound (§10.14.1) and the results must
be serialized back, so much of the wall time is transfer rather
than computation, the \(\alpha p\) term
made concrete.
Fork versus socket matters more than core count for small tasks. Socket clusters pay process startup (tens to hundreds of milliseconds) plus serialization of every argument and result. For a task lasting less than roughly a second, that overhead exceeds anything parallelism can recover.
library(foreach); library(doParallel)
cl <- makePSOCKcluster(n_cores)
registerDoParallel(cl)
on.exit({ stopCluster(cl); registerDoSEQ() }, add = TRUE)
c(registered_workers = getDoParWorkers())#> registered_workers
#> 1per <- TOTAL %/% n_cores
t_fe <- system.time(
res_fe <- foreach(i = seq_len(n_cores), .combine = "c") %dopar% rnorm(per))
c(values_generated = length(res_fe), seconds = round(t_fe[["elapsed"]], 3))#> values_generated seconds
#> 8.0e+06 2.4e-01registerDoSEQ() restores the sequential backend. Leaving
a parallel backend registered is a common source of confusing results
later in a session, because foreach is global state.
20.1 GPU computing
GPUs supply thousands of simple cores and, more importantly, several times the memory bandwidth of a CPU. They help when a computation is massively data-parallel with high arithmetic intensity, dense linear algebra, convolutions, elementwise tensor operations.
They do not help with branchy, sequential, or low-intensity work, and every transfer across the PCIe bus costs microseconds. A kernel that is memory-bound on the CPU is usually memory-bound on the GPU too; the gain is the bandwidth ratio, not the core-count ratio.
In R the current options are torch (LibTorch bindings,
CUDA and Metal backends, no Python required), gpuR (OpenCL,
so vendor-neutral), and cuda.ml (RAPIDS cuML bindings).
# --- GPU tensor operations via torch (not evaluated: requires a GPU) --------
library(torch)
if (cuda_is_available()) {
n <- 4096
a_cpu <- torch_randn(n, n)
a_gpu <- a_cpu$cuda()
t_cpu <- system.time(a_cpu$matmul(a_cpu))[["elapsed"]]
t_gpu <- system.time({ a_gpu$matmul(a_gpu); cuda_synchronize() })[["elapsed"]]
# Matrix multiply has arithmetic intensity ~ n/12, so it is compute-bound
# and benefits from the GPU. sum() would not.
c(cpu_seconds = t_cpu, gpu_seconds = t_gpu, speedup = t_cpu / t_gpu)
}21 Compiled code, and the accuracy it can cost
Rcpp compiles C++ and exposes it to R. It wins
decisively for iterative, branchy, scalar work, exactly
what R’s interpreter handles worst.
library(Rcpp)
# Three summation algorithms with different error behaviour
cppFunction('
double sum_naive(NumericVector x) {
double total = 0.0;
for (int i = 0; i < x.size(); ++i) total += x[i];
return total;
}')
cppFunction('
double sum_kahan(NumericVector x) {
double total = 0.0, c = 0.0;
for (int i = 0; i < x.size(); ++i) {
double y = x[i] - c; // recover the low-order bits lost last step
double t = total + y;
c = (t - total) - y; // what was lost this step
total = t;
}
return total;
}')
cppFunction('
double sum_pairwise(NumericVector x, int lo, int hi) {
if (hi - lo <= 128) {
double s = 0.0;
for (int i = lo; i < hi; ++i) s += x[i];
return s;
}
int mid = lo + (hi - lo) / 2;
return sum_pairwise(x, lo, mid) + sum_pairwise(x, mid, hi);
}')set.seed(141)
# A deliberately ill-conditioned sum: huge values that cancel, plus small ones
n_acc <- 1e6
x_acc <- c(1e16, rep(1.0, n_acc), -1e16)
exact <- n_acc # the answer is exactly n_acc
data.frame(
method = c("naive C++ loop", "Kahan compensated", "pairwise recursive",
"base::sum", "base::mean * n"),
result = c(sum_naive(x_acc), sum_kahan(x_acc),
sum_pairwise(x_acc, 0, length(x_acc)),
sum(x_acc), mean(x_acc) * length(x_acc)),
exact = exact) |>
mutate(absolute_error = abs(result - exact),
relative_error = signif(abs(result - exact) / exact, 3))Common misconception: “the C++ version is faster, so it is better.” The naive C++ loop above returns zero on this input, an error of 100%. Sequential summation accumulates rounding error at worst \(O(n\varepsilon)\); pairwise summation reduces this to \(O(\varepsilon\log n)\); Kahan compensation tracks the lost low-order bits and achieves \(O(\varepsilon)\) independent of \(n\).
R’s
base::sumuses extended (long double) accumulation internally, andbase::meanperforms a two-pass refinement, \(\hat\mu=\bar x+\frac1n\sum(x_i-\bar x)\), which recovers accuracy the first pass lost. Neither is a naive loop, and a hand-written replacement that is “faster” is usually faster because it omits that protection.Rounding the difference to a few decimals before comparing hides exactly the evidence that matters.
set.seed(143)
x_bench <- rnorm(1e7) # 80 MB, not 800 MB
tm <- bench::mark(
`base::sum` = sum(x_bench),
`naive C++ loop` = sum_naive(x_bench),
`Kahan compensated` = sum_kahan(x_bench),
`pairwise recursive` = sum_pairwise(x_bench, 0, length(x_bench)),
check = FALSE, iterations = 20, memory = FALSE)
data.frame(method = as.character(tm$expression),
median_sec = round(as.numeric(tm$median), 5),
elements_per_sec = signif(length(x_bench) / as.numeric(tm$median), 3),
GB_per_sec = round(8 * length(x_bench) / as.numeric(tm$median) / 1e9, 2))Report throughput, not raw seconds: elements per second and GB/s are comparable across machines in a way that wall time is not. And note where the GB/s figures land, near the machine’s memory bandwidth, confirming that this is a memory-bound kernel and that no amount of arithmetic cleverness will change it.
# Where C++ genuinely wins: a scalar, iterative, branchy loop
cppFunction('
int collatz_steps(int n) {
int k = 0;
while (n != 1) { n = (n % 2 == 0) ? n / 2 : 3 * n + 1; ++k; }
return k;
}')
collatz_R <- function(n) { k <- 0L; while (n != 1) {
n <- if (n %% 2 == 0) n %/% 2 else 3 * n + 1; k <- k + 1L }; k }
bench::mark(
`R loop` = sum(vapply(1:20000, collatz_R, numeric(1))),
`C++ loop` = sum(vapply(1:20000, collatz_steps, numeric(1))),
iterations = 5
)[, c("expression", "median")]This is the shape of problem Rcpp is for: no vectorized
formulation exists, the loop is scalar and data-dependent, and the
interpreter overhead dominates. For sum() and
mean(), already compiled, already memory-bound, already
numerically careful, there is nothing to win and accuracy to lose.
22 PART V: LANGUAGE INTEROPERABILITY
R Markdown executes many languages in one document. The full list of engines:
#> [1] "asis" "asy" "awk" "bash" "block" "block2"
#> [7] "bslib" "c" "cat" "cc" "coffee" "comment"
#> [13] "css" "ditaa" "dot" "embed" "eviews" "exec"
#> [19] "fortran" "fortran95" "gawk" "glue" "glue_sql" "gluesql"#> total_engines
#> 5523 R and Python
reticulate embeds a Python interpreter in the R session,
converting objects between the two automatically.
Never hard-code an interpreter path.
Sys.setenv(RETICULATE_PYTHON = ...)pointing at one machine’s install makes the document unrunnable elsewhere. Letreticulatediscover or create an environment:virtualenv_create()anduse_virtualenv()are reproducible, portable, and declare their dependencies explicitly.
# --- One-time environment setup (not evaluated during the knit) ------------
library(reticulate)
virtualenv_create("dspa", packages = c("numpy", "pandas", "scikit-learn"))
use_virtualenv("dspa", required = TRUE)
py_config()PY_OK <- isTRUE(tryCatch({
library(reticulate)
reticulate::py_available(initialize = TRUE) &&
reticulate::py_module_available("sklearn")
}, error = function(e) FALSE))
c(python_with_sklearn_available = PY_OK)#> python_with_sklearn_available
#> TRUEimport pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import accuracy_score
df = r.iris # R object, seen from Python
df["Species"] = df["Species"].astype(str) # kill the Categorical label
## The R factor iris$Species becomes a pandas Categorical when reticulate pulls it across (r.iris). train_test_split(stratify=y) then keeps it Categorical in y_te. When the next R chunk reaches for py$y_te, reticulate's converter for pandas...Categorical calls x$get_values(), and get_values() was removed in modern pandas. That's the AttributeError you're seeing in the trace at step 15.
# Nothing Categorical should cross the R↔Python boundary. Same trap fires for py$pred if you ever stratify a non-string label column.
X = df.drop("Species", axis=1)
y = df["Species"]
X_tr, X_te, y_tr, y_te = train_test_split(
X, y, test_size=0.4, random_state=4321, stratify=y)
clf = DecisionTreeClassifier(random_state=4321, max_depth=3).fit(X_tr, y_tr)
pred = clf.predict(X_te)
acc = accuracy_score(y_te, pred)
print(f"Python decision tree accuracy: {acc:.4f}")#> Python decision tree accuracy: 0.9667# Python -> R: read out of the `py` environment
res <- data.frame(py$X_te) |>
mutate(truth = as.character(py$y_te),
predicted = as.character(py$pred),
correct = truth == predicted)
c(test_cases = nrow(res), accuracy = round(mean(res$correct), 4))#> test_cases accuracy
#> 60.0000 0.9667ggplot(res, aes(Petal.Length, Petal.Width, color = truth, shape = predicted)) +
geom_point(size = 3, alpha = 0.85) +
scale_color_brewer(palette = "Dark2") +
labs(title = "Model fitted in Python, evaluated and plotted in R",
subtitle = "color is the true species; shape is the predicted one. Mismatches are visible directly",
x = "Petal length", y = "Petal width",
color = "True", shape = "Predicted") +
theme_dspa()
# --- Interactive equivalent ------------------------------------------------
plot_ly(res, x = ~Petal.Length, y = ~Petal.Width,
color = ~truth, symbol = ~predicted,
type = "scatter", mode = "markers", marker = list(size = 14)) |>
layout(title = "Python prediction, visualized in R",
xaxis = list(title = "Petal length"),
yaxis = list(title = "Petal width"),
legend = list(orientation = "h"))# Calling an external Python file. The file is written to tempdir(), never to
# a hard-coded location.
py_file <- file.path(tempdir(), "helpers.py")
writeLines(c(
"def even_numbers(values):",
" \"\"\"Return only the even entries of a numeric sequence.\"\"\"",
" return [v for v in values if v % 2 == 0]",
"",
"def summarize(values):",
" n = len(values)",
" m = sum(values) / n if n else float('nan')",
" return {'n': n, 'mean': m, 'min': min(values), 'max': max(values)}"
), py_file)
reticulate::source_python(py_file)
even_numbers(c(1, 2, 3, 4, 5, 6, 7, 9, 2.2, 1e15))#> [1] 2e+00 4e+00 6e+00 1e+15#> List of 4
#> $ n : int 8
#> $ mean: num 3.88
#> $ min : num 1
#> $ max : num 9The Python list is converted to an R vector and the dictionary to a named list, automatically. Conversion is not free, large objects are copied, so for performance-critical loops it is better to keep the data on one side.
24 R and C++
§10.20 covered cppFunction() for
inline definitions. For anything substantial, keep the C++ in a file and
sourceCpp() it.
cpp_file <- file.path(tempdir(), "stats.cpp")
writeLines(c(
"#include <Rcpp.h>",
"using namespace Rcpp;",
"",
"// [[Rcpp::export]]",
"NumericVector rolling_mean_cpp(NumericVector x, int k) {",
" int n = x.size();",
" NumericVector out(n, NA_REAL);",
" if (k > n) return out;",
" double run = 0.0;",
" for (int i = 0; i < k; ++i) run += x[i];",
" out[k - 1] = run / k;",
" for (int i = k; i < n; ++i) {",
" run += x[i] - x[i - k]; // O(1) update, not O(k) recomputation",
" out[i] = run / k;",
" }",
" return out;",
"}",
"",
"/*** R",
"# R code embedded in the C++ source, run automatically after compilation",
"cat('rolling_mean_cpp compiled and available\\n')",
"*/"
), cpp_file)
Rcpp::sourceCpp(cpp_file)#>
#> > cat("rolling_mean_cpp compiled and available\n")
#> rolling_mean_cpp compiled and availableset.seed(151)
x_roll <- rnorm(2e5); k_roll <- 50
rolling_mean_R <- function(x, k) {
n <- length(x); out <- rep(NA_real_, n)
for (i in k:n) out[i] <- mean(x[(i - k + 1):i]) # O(k) per position
out
}
r1 <- rolling_mean_cpp(x_roll, k_roll)
r2 <- as.numeric(stats::filter(x_roll, rep(1 / k_roll, k_roll), sides = 1))
c(max_abs_difference = max(abs(r1 - r2), na.rm = TRUE))#> max_abs_difference
#> 6.21725e-15bench::mark(
`R loop, O(nk)` = rolling_mean_R(x_roll, k_roll)[k_roll:length(x_roll)],
`C++ sliding, O(n)` = rolling_mean_cpp(x_roll, k_roll)[k_roll:length(x_roll)],
`stats::filter (FFT/C)`= r2[k_roll:length(x_roll)],
check = FALSE, iterations = 10
)[, c("expression", "median", "mem_alloc")]The C++ version wins on the algorithm, not the
language. The R loop recomputes a \(k\)-element mean at every position, \(O(nk)\), while the C++ version maintains a
running sum with an \(O(1)\) update,
giving \(O(n)\). A vectorized R
solution using filter() is competitive, because it too
avoids the quadratic work. Reach for a better algorithm before
reaching for a faster language.
25 PART VI: SYNTHESIS
26 Computational complexity summary
\(n\) = rows, \(c\) = columns, \(c'\) = columns requested, \(m\) = edges, \(V\) = vertices, \(p\) = processors, \(N\) = stream length, \(k\)/\(w\)/\(d\) = sketch parameters.
| Operation | Time | Memory | Note |
|---|---|---|---|
| CSV read | \(O(n\,c)\) | \(O(n\,c)\) | Must parse every byte, even for one column |
| Parquet read, all columns | \(O(n\,c)\) | \(O(n\,c)\) | Plus decompression |
| Parquet read, \(c'\) columns | \(\mathbf{O(n\,c')}\) | \(O(n\,c')\) | Column pruning is structural |
| Parquet with predicate | \(O(n_{\text{matching}}\,c')\) | — | Row-group statistics skip chunks unread |
| SQL aggregation, pushed down | \(O(n)\) server-side | \(O(\text{groups})\) transferred | The answer crosses the wire, not the data |
SQL SELECT * then aggregate |
\(O(n)\) | \(\mathbf{O(n\,c)}\) transferred | Defeats the engine |
Out-of-core scan (arrow/duckdb) |
\(O(n)\) | \(O(\text{working set})\) | Streams over files |
| Graph as adjacency matrix | — | \(\Theta(V^2)\) | 80 GB at \(V=10^5\) |
| Graph as adjacency list | — | \(\Theta(V+m)\) | 17 MB at \(V=10^5\), \(m=10^6\) |
| Degree centrality | \(O(V)\) | \(O(V)\) | — |
| Betweenness (Brandes) | \(\mathbf{O(Vm)}\) | \(O(V+m)\) | \(O(km)\) with \(k\) sampled sources |
| Louvain community detection | \(O(m\log V)\) | \(O(V+m)\) | Greedy; Leiden repairs its defect |
| Reservoir sample | \(O(N)\) | \(\mathbf{O(k)}\) | Exact uniformity, one pass |
| Count-min sketch | \(O(Nd)\) | \(\mathbf{O(wd)}\) | \(f\le\hat f\le f+\varepsilon N\) w.p. \(1-\delta\) |
| HyperLogLog | \(O(N)\) | \(\mathbf{O(m\log\log N)}\) | Relative error \(\approx1.04/\sqrt m\) |
| Stream micro-clustering | \(O(N)\) | \(O(\text{micro-clusters})\) | Independent of \(N\) |
| Prequential evaluation | \(O(N)\) | \(O(\text{window})\) | Test-then-train |
| Sum, naive | \(O(n)\), error \(O(n\varepsilon)\) | \(O(1)\) | Memory-bound, \(I=0.125\) |
| Sum, pairwise | \(O(n)\), error \(O(\varepsilon\log n)\) | \(O(\log n)\) stack | — |
| Sum, Kahan | \(O(n)\), error \(O(\varepsilon)\) | \(O(1)\) | ~4× the arithmetic, still memory-bound |
| Matrix multiply | \(O(n^3)\) | \(O(n^2)\) | \(I\sim n/12\): compute-bound, GPU-friendly |
| Parallel with overhead | \(T_s+T_p/p+\alpha p\) | — | Optimal \(p^\star=\sqrt{T_p/\alpha}\) |
Four rules follow.
Push computation to the data. Predicate, projection, and aggregation pushdown turn \(O(nc)\) transfers into \(O(\text{groups})\) transfers. The gain is proportional to the data and therefore grows exactly where it matters.
Memory beats arithmetic. Reductions like
sum and mean have arithmetic intensity near
\(0.1\) flops per byte, far below any
machine’s ridge point. They run at memory bandwidth, and no rewriting of
the inner loop changes that.
Sketches trade a bounded error for unbounded scale. \(O(1)\) or \(O(\log\log N)\) memory with explicit \((\varepsilon,\delta)\) guarantees is what makes unbounded streams tractable at all.
Speedup has a ceiling, and then a peak. Amdahl caps it at \(1/(1-f)\); overhead means it then declines past \(p^\star=\sqrt{T_p/\alpha}\). Measure both before buying cores.
27 Common pitfalls
| # | Pitfall | Consequence | Fix |
|---|---|---|---|
| 1 | Setting knitr::opts_knit$set(root.dir = ...) to a
personal path |
The document knits on one machine only | Never set it; use tempdir() |
| 2 | Absolute paths anywhere in a document | Same | tempdir(), file.path(), a cache
helper |
| 3 | Unguarded network calls | One unavailable service breaks the whole build | Wrap in tryCatch with a stated fallback |
| 4 | memory.size() / memory.limit() |
Defunct since R 4.2.0 | Remove; the OS manages memory |
| 5 | View() inside a knitted chunk |
Interactive viewer in a non-interactive session | head(), str(),
knitr::kable() |
| 6 | Reading a whole table then filtering in R | Transfers \(O(nc)\) to compute \(O(\text{groups})\) | Push the predicate and aggregation down |
| 7 | SELECT * from a database |
The database becomes an expensive file | Select columns; aggregate server-side |
| 8 | Leaving a connection open | Leaked handles; exhausted server slots | on.exit(dbDisconnect(con)) |
| 9 | CSV for analytical workloads | No column pruning, no types, no compression | Parquet or Feather |
| 10 | Storing a sparse graph as a dense matrix | \(\Theta(V^2)\) memory; 80 GB at \(V=10^5\) | Adjacency list / sparse matrix |
| 11 | Exact betweenness on a large graph | \(O(Vm)\); hours to days | Sample sources, or use a cutoff |
| 12 | Reading a maximum off a printed window | Says nothing about the global maximum | which.max() over all elements |
| 13 | Purity as a clustering quality measure | Monotone in \(k\); maximal at \(k=N\) | Adjusted Rand index |
| 14 | Passing the class label as a stream coordinate | Circular: clustered on the label, evaluated against it | Use the class = argument |
| 15 | Batch evaluation of a streaming model | Ignores order and drift | Prequential (test-then-train) |
| 16 | A fixed window under drift | Too long adapts slowly; too short is noisy | Detect drift (Page–Hinkley, ADWIN) and reset |
| 17 | Optimizing without profiling | Effort spent where it does not matter | bench::mark(), profvis first |
| 18 | Growing a vector in a loop | \(O(n^2)\) copying plus garbage collection | Pre-allocate, or vectorize |
| 19 | Assuming \(p\) cores gives \(p\times\) speedup | Amdahl caps it at \(1/(1-f)\) | Compute the ceiling before buying cores |
| 20 | Ignoring parallel overhead | Past \(p^\star=\sqrt{T_p/\alpha}\), more workers are slower | Model \(T_s+T_p/p+\alpha p\) |
| 21 | Benchmarking arms with unequal work | Measures workload, not speedup | Divide the same total across workers |
| 22 | Leaving a foreach backend registered |
Later code silently parallel or serial | registerDoSEQ() in on.exit() |
| 23 | Replacing sum/mean with a naive C++
loop |
Faster at best marginally, much less accurate | Compensated or pairwise summation |
| 24 | Rounding before comparing numerical results | Hides the error being investigated | Compare at full precision; report relative error |
28 Practice problems
28.1 Problem 1: Where is the column-pruning break-even?
Find the fraction of columns at which reading Parquet stops beating CSV.
Solution
set.seed(201)
n1 <- 1e5; c1 <- 40
d1 <- as.data.frame(matrix(rnorm(n1 * c1), n1, c1))
f1c <- file.path(tempdir(), "p1.csv"); f1p <- file.path(tempdir(), "p1.parquet")
data.table::fwrite(d1, f1c); arrow::write_parquet(d1, f1p)
fracs <- c(1, 2, 5, 10, 20, 40)
p1 <- do.call(rbind, lapply(fracs, function(k) {
sel <- names(d1)[1:k]
tc <- system.time(data.table::fread(f1c, select = sel, showProgress = FALSE))[["elapsed"]]
tp <- system.time(arrow::read_parquet(f1p, col_select = all_of(sel)))[["elapsed"]]
data.frame(cols = k, frac = k / c1, csv = tc, parquet = tp, ratio = tc / tp)
}))
p1 |> mutate(across(where(is.numeric), \(z) round(z, 4)))p1 |> select(frac, CSV = csv, Parquet = parquet) |>
pivot_longer(-frac, names_to = "format", values_to = "sec") |>
ggplot(aes(frac, sec, color = format)) +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_x_continuous(labels = scales::percent) +
scale_color_manual(values = c(CSV = "#D8433B", Parquet = "#3B7DD8")) +
labs(title = "Read time against the fraction of columns requested",
x = "Columns requested", y = "Seconds", color = NULL) +
theme_dspa()
28.2 Problem 2: Measure predicate pushdown
Quantify the bytes saved by filtering at the source versus in R.
Solution
set.seed(203)
n2 <- 3e5
d2 <- data.frame(id = 1:n2, grp = sample(letters[1:26], n2, TRUE),
v = rnorm(n2), w = rnorm(n2))
f2 <- file.path(tempdir(), "p2.parquet"); arrow::write_parquet(d2, f2)
con2 <- DBI::dbConnect(duckdb::duckdb())
# on.exit(DBI::dbDisconnect(con2, shutdown = TRUE), add = TRUE)
pushed <- DBI::dbGetQuery(con2, sprintf(
"SELECT grp, COUNT(*) n, AVG(v) mv FROM read_parquet('%s')
WHERE v > 2 GROUP BY grp ORDER BY grp", f2))
pulled <- DBI::dbGetQuery(con2, sprintf("SELECT * FROM read_parquet('%s')", f2))
c(rows_returned_pushdown = nrow(pushed),
rows_returned_pull = nrow(pulled),
bytes_pushdown = as.numeric(object.size(pushed)),
bytes_pull = as.numeric(object.size(pulled)),
reduction_factor = round(as.numeric(object.size(pulled)) /
as.numeric(object.size(pushed))))#> rows_returned_pushdown rows_returned_pull bytes_pushdown
#> 26 300000 3064
#> bytes_pull reduction_factor
#> 8402544 274228.3 Problem 3: Verify the count-min guarantee
Check empirically that the sketch never underestimates and that overestimates respect \(\varepsilon N\).
Solution
set.seed(205)
N3 <- 15000
items3 <- sample(paste0("k", 1:2000), N3, TRUE, prob = 1 / (1:2000))
truth3 <- table(items3)
check_eps <- function(eps, delta = 0.01) {
sk <- cms_new(w = ceiling(exp(1) / eps), d = ceiling(log(1 / delta)), seed = 207)
for (s in items3) sk <- cms_add(sk, s)
est <- vapply(names(truth3), \(s) cms_query(sk, s), numeric(1))
err <- est - as.numeric(truth3)
c(eps = eps, width = sk$w, bytes = as.numeric(object.size(sk$C)),
min_error = min(err), max_error = max(err),
bound = ceiling(eps * N3),
all_within_bound = all(err >= 0 & err <= eps * N3))
}
as.data.frame(do.call(rbind, lapply(c(0.02, 0.005, 0.001), check_eps)))min_error is never negative, the sketch
cannot underestimate, because collisions only add.
max_error stays below \(\varepsilon N\), and shrinking \(\varepsilon\) tightens the bound at the
cost of a wider table. The memory-accuracy trade is explicit and
controllable.
28.4 Problem 4: Estimate the parallel fraction
Measure speedup for several worker counts and fit Amdahl’s law to recover \(f\).
Solution
library(parallel)
nc <- max(2, min(4, parallel::detectCores() - 1))
TOT <- 4e6
time_with <- function(k) {
if (k == 1) return(system.time(rnorm(TOT))[["elapsed"]])
cl <- makePSOCKcluster(k); on.exit(stopCluster(cl), add = TRUE)
per <- TOT %/% k
system.time(unlist(parLapply(cl, seq_len(k), function(i) rnorm(per))))[["elapsed"]]
}
set.seed(209)
ks <- 1:nc
tt <- vapply(ks, \(k) median(replicate(3, time_with(k))), numeric(1))
sp <- tt[1] / tt
# Amdahl: S = 1/((1-f) + f/p) => f = (1 - 1/S) / (1 - 1/p)
f_hat <- mean(((1 - 1 / sp) / (1 - 1 / ks))[ks > 1])
data.frame(workers = ks, seconds = round(tt, 4), speedup = round(sp, 3),
amdahl_predicted = round(1 / ((1 - f_hat) + f_hat / ks), 3))c(estimated_parallel_fraction = round(f_hat, 3),
amdahl_ceiling = round(1 / (1 - max(f_hat, 0.001)), 1))#> estimated_parallel_fraction amdahl_ceiling
#> 0.318 1.500ggplot(data.frame(p = ks, observed = sp,
amdahl = 1 / ((1 - f_hat) + f_hat / ks)) |>
pivot_longer(-p, names_to = "series", values_to = "s"),
aes(p, s, color = series)) +
geom_abline(slope = 1, intercept = 0, linetype = "dotted", color = "grey55") +
geom_line(linewidth = 1) + geom_point(size = 2.4) +
scale_color_manual(values = c(observed = "#3B7DD8", amdahl = "#D8433B")) +
labs(title = sprintf("Observed speedup and the fitted Amdahl curve (f = %.2f)", f_hat),
subtitle = "Dotted line is perfect linear speedup",
x = "Workers", y = "Speedup", color = NULL) +
theme_dspa()
rnorm is memory-bound and the results must be
serialized back to the parent. The implied ceiling is the number to
check before provisioning more cores.
28.5 Problem 5: Summation accuracy across algorithms
Construct inputs where naive summation fails and quantify each method’s error.
Solution
cases <- list(
`well conditioned` = { set.seed(211); rnorm(1e6) },
`all positive, wide range` = { set.seed(213); c(1e12, runif(1e6, 0, 1)) },
`catastrophic cancellation` = c(1e16, rep(1, 1e5), -1e16))
res5 <- do.call(rbind, lapply(names(cases), function(nm) {
x <- cases[[nm]]
exact <- sum(as.numeric(Rmpfr::mpfr(x, 200))) # high-precision reference
data.frame(case = nm,
naive = sum_naive(x), kahan = sum_kahan(x),
pairwise = sum_pairwise(x, 0, length(x)), base_R = sum(x),
exact = exact)
}))
res5 |>
mutate(across(c(naive, kahan, pairwise, base_R),
~ signif(abs(.x - exact) / pmax(abs(exact), 1), 3),
.names = "relerr_{.col}")) |>
select(case, starts_with("relerr_"))Naive summation is fine on well-conditioned input and fails
completely under cancellation. Kahan and pairwise remain accurate
throughout, and base::sum’s extended-precision accumulator
does too. The right algorithm depends on the conditioning of the
input, not on the language it is written in.
Rmpfr is unavailable, substitute the analytically known
exact value for each case, \(10^5\) for
the third, and the conclusion is unchanged.)
28.6 Problem 6: Find the optimal worker count
Fit the overhead model \(T(p)=T_s+T_p/p+\alpha p\) and compare \(p^\star\) against the measured minimum.
Solution
set.seed(215)
nc6 <- max(2, min(8, parallel::detectCores()))
work <- 2e6
t6 <- vapply(1:nc6, function(k) {
if (k == 1) return(median(replicate(3, system.time(rnorm(work))[["elapsed"]])))
cl <- makePSOCKcluster(k); on.exit(stopCluster(cl), add = TRUE)
per <- work %/% k
median(replicate(3, system.time(
unlist(parLapply(cl, seq_len(k), function(i) rnorm(per))))[["elapsed"]]))
}, numeric(1))
# Fit T(p) = Ts + Tp/p + alpha*p by nonlinear least squares
d6 <- data.frame(p = 1:nc6, t = t6)
fit6 <- try(nls(t ~ Ts + Tp / p + alpha * p, data = d6,
start = list(Ts = 0.01, Tp = t6[1], alpha = 0.01)), silent = TRUE)
if (!inherits(fit6, "try-error")) {
co <- coef(fit6)
p_star <- sqrt(co[["Tp"]] / max(co[["alpha"]], 1e-8))
c(Ts = round(co[["Ts"]], 4), Tp = round(co[["Tp"]], 4),
alpha = round(co[["alpha"]], 5),
p_star_predicted = round(p_star, 2),
p_observed_minimum = which.min(t6))
ggplot(d6, aes(p, t)) +
geom_line(aes(y = predict(fit6)), color = "#D8433B", linewidth = 1) +
geom_point(size = 2.6, color = "#3B7DD8") +
geom_vline(xintercept = p_star, linetype = "dashed", color = "grey40") +
labs(title = "Fitted overhead model and the predicted optimum",
subtitle = expression("Dashed line: "*p^"*"*" = sqrt(T[p]/alpha)"),
x = "Workers p", y = "Seconds") +
theme_dspa()
}
28.7 Problem 7: Drift detection delay versus false alarms
Trace the operating characteristic of Page–Hinkley across thresholds.
Solution
set.seed(217)
make_err <- function(drift_at = 1500, n = 3000, p0 = 0.15, p1 = 0.45)
rbinom(n, 1, ifelse(seq_len(n) <= drift_at, p0, p1))
lambdas <- c(3, 6, 12, 25, 50, 100)
roc7 <- do.call(rbind, lapply(lambdas, function(lam) {
reps <- t(replicate(30, {
e <- make_err()
a <- page_hinkley(e, delta = 0.005, lambda = lam)$alarms
false_alarm <- sum(a <= 1500)
first_true <- a[a > 1500][1]
c(false_alarms = false_alarm,
delay = if (is.na(first_true)) NA else first_true - 1500)
}))
data.frame(lambda = lam,
mean_false_alarms = mean(reps[, 1]),
mean_delay = mean(reps[, 2], na.rm = TRUE),
detection_rate = mean(!is.na(reps[, 2])))
}))
roc7 |> mutate(across(where(is.numeric), \(z) round(z, 2)))ggplot(roc7, aes(mean_false_alarms, mean_delay)) +
geom_path(linewidth = 0.9, color = "grey55") +
geom_point(aes(color = factor(lambda)), size = 3.4) +
scale_color_viridis_d(option = "plasma", end = 0.9, name = expression(lambda)) +
labs(title = "Drift detection: delay against false alarms",
subtitle = "Lower threshold detects sooner and cries wolf more often -- the usual tradeoff",
x = "Mean false alarms before the true drift",
y = "Mean detection delay (steps)") +
theme_dspa()
28.8 Problem 8: Locate a kernel on the roofline
Measure the achieved bandwidth of three operations and place each in its regime.
Solution
set.seed(219)
n8 <- 5e6
a8 <- rnorm(n8); b8 <- rnorm(n8)
m8 <- matrix(rnorm(600 * 600), 600, 600)
measure <- function(expr, bytes, flops, label) {
t <- median(replicate(7, system.time(force(expr))[["elapsed"]]))
c(operation = label,
seconds = round(t, 4),
GB_per_sec = round(bytes / t / 1e9, 2),
GFLOP_per_sec = round(flops / t / 1e9, 3),
intensity = round(flops / bytes, 3))
}
as.data.frame(rbind(
measure(sum(a8), 8 * n8, n8, "sum(x)"),
measure(a8 + b8, 24 * n8, n8, "x + y"),
measure(m8 %*% m8, 24 * 600^2, 2 * 600^3, "600x600 matmul")))sum and + achieve high GB/s
and negligible GFLOP/s, memory-bound, running near the
machine’s bandwidth. Matrix multiply is the reverse: intensity around
\(n/12=50\), so it is compute-bound and
its GFLOP/s figure is the meaningful one. The intensity column
tells you which number to read.
29 Checkpoint
- Your query needs 3 of 200 columns from a 10-million-row table. Why does the storage format matter more than the disk speed?
- A colleague’s script runs on their machine and fails on yours with “cannot open file”. What is the most likely single cause?
- You have 32 cores and a workload that is 90% parallelizable. What speedup should you expect, and what would you need to change to do better?
- A streaming clustering reports purity 0.97. What do you conclude?
- Someone replaces
sum(x)with a hand-written C++ loop and reports a 20% speedup. What do you check? - Your parallel job runs slower with 16 workers than with 4. What model explains this, and what is the fix?
Answers
- Because a columnar format reads only the requested columns, while a row-major format such as CSV must scan every byte of every row to locate the delimiters. The I/O differs by roughly \(200/3 \approx 67\times\), a structural property of the layout that no amount of disk speed changes. Parquet additionally carries per-row-group min/max statistics, so a predicate can skip whole chunks unread, and it compresses per column with type-aware encodings.
- An absolute path, or a working directory set to
somewhere that exists only on their machine. The usual culprits are a
literal
"C:/Users/name/Desktop/..."andknitr::opts_knit$set(root.dir = ...)hidden in aninclude=FALSEchunk. Everything should go throughtempdir(),file.path(), or a cache helper, and network calls should be wrapped so an unavailable service degrades rather than fails. - Amdahl’s law gives \(S(32)=\dfrac{1}{0.1+0.9/32}\approx 7.8\), against a ceiling of \(1/(1-0.9)=10\) at infinite cores. So 32 cores buys you at most 7.8× and the next 32 buy almost nothing. To do better you must reduce the serial fraction, parallelize the remaining 10%, or restructure the algorithm. Buying more cores cannot help. If instead the problem size can grow with the resources, Gustafson’s law applies and the outlook is far better.
- Very little. Purity increases monotonically with the number of clusters and reaches exactly 1 when every point is its own cluster, so a high value may reflect nothing but a large \(k\). It is also uncorrected for chance. Ask for the number of clusters and the adjusted Rand index, which is 0 for a random partition and does not reward splitting.
- The accuracy, at full precision.
base::sumaccumulates in extended (long double) precision andbase::meanperforms a two-pass refinement; a naive C++ loop does neither, and its worst-case error grows as \(O(n\varepsilon)\) against \(O(\varepsilon)\) for a compensated sum. Test on an ill-conditioned input, large values that cancel, and compare without rounding first. Also check the arithmetic intensity: at 0.125 flops per byte,sumis memory-bound, so a 20% difference is more likely compiler flags or cache behaviour than algorithmic merit. - The overhead model \(T(p)=T_s+T_p/p+\alpha p\): the \(\alpha p\) term grows with the worker count, covering process startup and the serialization of every argument and result. Minimizing gives \(p^\star=\sqrt{T_p/\alpha}\), past which more workers increase total time. Fixes: use fork clusters instead of socket where the platform allows (no serialization), send less data per task, make tasks coarser so the fixed cost amortizes, or simply use fewer workers. And verify that each arm does the same total work, a benchmark where the parallel version does \(p\) times the work measures nothing.
30 Summary
Formats and ingestion
- Row-major and columnar layouts favour opposite access patterns. Analytical queries touch few columns and many rows, which is what columnar storage is built for.
- Parquet and Arrow give column pruning, predicate pushdown from row-group statistics, and type-aware compression. CSV gives none of these and must be parsed byte by byte.
- Push computation to the data. Filtering and aggregating at the source transfers the answer, not the input, and the saving grows with the data.
- Out-of-core engines (
arrow,duckdb) stream over columnar files behind adplyrfront end, replacing disk-backed data structures. - Reproducibility is a portability property: no absolute paths, no working directory, guarded network calls.
Network data
- Adjacency matrices cost \(\Theta(V^2)\) and adjacency lists \(\Theta(V+m)\). Real networks are sparse, so the difference is possible versus impossible.
- Betweenness is \(O(Vm)\) even with Brandes’ algorithm; sample sources for large graphs.
- Modularity compares against a degree-preserving null, the same null-comparison logic used for cluster and rule validation elsewhere.
Data streams
- The streaming model forbids multiple passes and linear memory, which forces sketches: compact summaries with \((\varepsilon,\delta)\) guarantees.
- Reservoir sampling gives exact uniformity in \(O(k)\) memory; count-min bounds frequency overestimates by \(\varepsilon N\); HyperLogLog estimates cardinality with relative error \(1.04/\sqrt m\) and memory independent of the stream.
- Evaluate streams prequentially, test then train, and prefer the adjusted Rand index to purity, which is monotone in \(k\).
- Window size trades stability against plasticity; detect drift rather than guessing a window.
Computational performance
- Memory latency spans six orders of magnitude. For low-intensity kernels the arithmetic is free and bandwidth determines runtime.
- The roofline locates a computation:
sumandmeansit at \(I\approx0.1\) and are memory-bound; matrix multiply has \(I\sim n/12\) and is compute-bound. - Profile before optimizing. Allocation patterns usually matter more than arithmetic.
- Amdahl caps speedup at \(1/(1-f)\); Gustafson describes the scaled problem. With overhead, time is minimized at \(p^\star=\sqrt{T_p/\alpha}\) and rises thereafter.
- Benchmarks must give every arm the same total work.
- A faster reduction can be a worse one: naive summation errs at \(O(n\varepsilon)\), pairwise at \(O(\varepsilon\log n)\), Kahan at \(O(\varepsilon)\). Compare at full precision.
- Reach for a better algorithm before a faster language: the \(O(n)\) sliding window beats the \(O(nk)\) loop in any language.
Where these threads continue
| Thread | Continues in |
|---|---|
| Regularization on large, wide data | Feature selection |
| Streaming forecasts and drift over time | Longitudinal analysis |
| Gradient methods and their arithmetic intensity | Function optimization |
| GPU tensor computation and autodiff | Deep learning |
Further practice. Repeat the network analysis on the Les Misérables co-appearance graph, whose edge weights count co-appearances by chapter, and profile a pipeline from an earlier chapter to locate where its time actually goes.
lm_path <- dspa_download(
"https://umich.instructure.com/files/330389/download?download_frd=1",
"lesmis.txt")
lesmis <- dspa_try(utils::read.table(lm_path, header = FALSE),
fallback = data.frame(V1 = integer(0), V2 = integer(0)),
label = "Les Miserables network")
if (nrow(lesmis) > 0) {
g_lm <- igraph::graph_from_edgelist(as.matrix(lesmis[, 1:2]), directed = FALSE)
c(vertices = gorder(g_lm), edges = gsize(g_lm),
density = signif(edge_density(g_lm), 4),
communities = length(igraph::cluster_louvain(g_lm)))
}#> vertices edges density communities
#> 77.00000 254.00000 0.08681 6.0000031 Chapter roadmap
- Chapter 1:
Foundations. R toolchain, reproducibility conventions,
dspa_read(), simulation. - Chapter 2: Data quality and exploratory visual analytics. Missingness, robust statistics, graphical perception.
- Chapter 3: Linear algebra, matrix computing, and regression. Floating point, conditioning, blocked algorithms.
- Chapter 4: Dimensionality reduction. PCA, randomized SVD, distance concentration.
- Chapter 5: Supervised classification. Bayes error, evaluation metrics, the leakage taxonomy.
- Chapter 6: Black-box methods. Neural networks, kernels, ensembles, complexity.
- Chapter 7: Text mining and association rules. Sparse matrices, Apriori, multiplicity control.
- Chapter 8: Unsupervised clustering. Internal and external validation, sparse spectral methods.
- Chapter 9: Model assessment and validation. Optimism, calibration, resampling, nested CV.
- Variable importance and feature selection. Ridge, LASSO, elastic net, stability selection.
- Longitudinal and time-series analysis. Mixed models, ARIMA, forecast evaluation.
- Function optimization. Gradient descent, duality, Bayesian optimization.
- Deep learning. Tensor computation, GPU training, representation learning.
32 Session information
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=English_United States.utf8
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] parallel stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] reticulate_1.38.0 Rcpp_1.0.14
#> [3] doParallel_1.0.17 iterators_1.0.14
#> [5] foreach_1.5.2 stream_2.0-2
#> [7] magrittr_2.0.3 igraph_2.0.3
#> [9] maps_3.4.2 WikidataQueryServiceR_1.0.0
#> [11] xml2_1.3.6 rvest_1.0.4
#> [13] jsonlite_1.8.9 dbplyr_2.5.0
#> [15] duckdb_1.2.1 DBI_1.2.3
#> [17] arrow_25.0.1 rio_1.1.1
#> [19] bench_1.1.4 data.table_1.16.4
#> [21] plotly_4.12.1 patchwork_1.3.0
#> [23] tidyr_1.3.1 dplyr_1.1.4
#> [25] ggplot2_4.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] writexl_1.5.0 rlang_1.1.5 clue_0.3-65
#> [4] otel_0.2.0 mlbench_2.1-5 compiler_4.3.3
#> [7] png_0.1-8 vctrs_0.6.5 stringr_1.5.1
#> [10] pkgconfig_2.0.3 crayon_1.5.3 fastmap_1.2.0
#> [13] labeling_0.4.3 promises_1.3.2 rmarkdown_2.31
#> [16] tzdb_0.4.0 haven_2.5.4 ps_1.9.0
#> [19] purrr_1.0.2 bit_4.0.5 xfun_0.52
#> [22] Rmpfr_0.9-5 cachem_1.1.0 clusterGeneration_1.3.8
#> [25] blob_1.2.4 later_1.4.1 gmp_0.7-4
#> [28] cluster_2.1.6 R6_2.6.1 bslib_0.9.0
#> [31] stringi_1.8.4 RColorBrewer_1.1-3 rpart_4.1.23
#> [34] jquerylib_0.1.4 assertthat_0.2.1 knitr_1.51
#> [37] R.utils_2.12.3 readr_2.1.5 Matrix_1.6-5
#> [40] tidyselect_1.2.1 rstudioapi_0.18.0 yaml_2.3.10
#> [43] codetools_0.2-20 websocket_1.4.1 curl_6.2.0
#> [46] processx_3.8.6 lattice_0.22-6 tibble_3.2.1
#> [49] withr_3.0.2 S7_0.2.1 evaluate_1.0.3
#> [52] archive_1.1.8 ratelimitr_0.4.1 proxy_0.4-27
#> [55] mclust_6.1.1 pillar_1.10.1 generics_0.1.3
#> [58] dbscan_1.2-0 vroom_1.6.5 rprojroot_2.1.1
#> [61] chromote_0.4.0 hms_1.1.3 scales_1.4.0
#> [64] glue_1.8.0 tools_4.3.3 forcats_1.0.0
#> [67] grid_4.3.3 crosstalk_1.2.1 cli_3.6.3
#> [70] profmem_0.7.0 viridisLite_0.4.2 gtable_0.3.6
#> [73] R.methodsS3_1.8.2 selectr_0.4-2 sass_0.4.9
#> [76] digest_0.6.37 htmlwidgets_1.6.4 farver_2.1.2
#> [79] htmltools_0.5.8.1 R.oo_1.26.0 lifecycle_1.0.5
#> [82] httr_1.4.7 here_1.0.2 bit64_4.0.5
#> [85] MASS_7.3-60.0.1