Performance

In caugi: Causal Graph Interface

#| label: setup
#| include: false
#| purl: false

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  fig.width = 7,
  fig.height = 4,
  warning = FALSE,
  message = FALSE,
  echo = FALSE
)

#| label: bench-setup
#| include: false
#| eval: false

suppressPackageStartupMessages({
  library(caugi)
  library(bench)
  library(data.table)
  library(jsonlite)
  library(igraph)
  library(bnlearn)
  library(dagitty)
  library(ggm)
  library(pcalg)
  library(graph)
})

FIXTURES_DIR <- "fixtures"
RESULTS_DIR <- "results"
dir.create(RESULTS_DIR, showWarnings = FALSE, recursive = TRUE)

spec <- jsonlite::fromJSON(
  file.path(FIXTURES_DIR, "spec.json"),
  simplifyVector = FALSE
)
.skip_table <- if (is.null(spec$skip)) list() else spec$skip

#| label: bench-load-fixture
#| include: false
#| eval: false

load_fixture <- function(fx) {
  edges_path <- file.path(FIXTURES_DIR, fx$edges_file)
  edges <- data.table::fread(
    edges_path,
    col.names = c("from", "to"),
    colClasses = "character"
  )
  nodes <- paste0("V", seq_len(fx$n))

  cg <- caugi::add_edges(
    caugi::caugi(class = "DAG", nodes = nodes),
    from = edges$from,
    edge = rep("-->", nrow(edges)),
    to = edges$to
  )
  cg <- caugi::build(cg)

  ig <- igraph::graph_from_data_frame(
    as.data.frame(edges),
    directed = TRUE,
    vertices = data.frame(name = nodes)
  )

  am <- as.matrix(igraph::as_adjacency_matrix(ig))
  storage.mode(am) <- "integer"

  bn <- bnlearn::empty.graph(nodes)
  bnlearn::amat(bn) <- am

  dg_str <- paste0(
    "dag {\n",
    paste(sprintf("%s -> %s", edges$from, edges$to), collapse = "\n"),
    "\n}"
  )
  dg <- dagitty::dagitty(dg_str)

  # pcalg's searchAM expects the PAG-style 0/1/2/3 amat coding (mark at the
  # column-end of edge {row, col}): for a -> b, amat[a,b] = 2 (arrowhead at b)
  # and amat[b,a] = 3 (tail at a).
  amat_pag <- matrix(0L, nrow(am), ncol(am), dimnames = dimnames(am))
  amat_pag[am == 1L] <- 2L
  amat_pag[t(am) == 1L] <- 3L

  # pcalg::dsep operates on a graphNEL.
  gNEL <- igraph::as_graphnel(ig)

  list(
    cg = cg,
    ig = ig,
    bn = bn,
    dg = dg,
    am = am,
    amat_pag = amat_pag,
    gNEL = gNEL,
    nodes = nodes
  )
}

#| label: bench-helpers
#| include: false
#| eval: false

# bench::mark returns a tibble; we convert each row into a long-format record.
mark_to_rows <- function(bm, operation, fx) {
  pkg <- as.character(bm$expression)
  data.table::data.table(
    language = rep("R", length(pkg)),
    package = pkg,
    operation = rep(operation, length(pkg)),
    fixture_id = rep(fx$id, length(pkg)),
    n = rep(fx$n, length(pkg)),
    p = rep(fx$p, length(pkg)),
    n_edges = rep(fx$n_edges, length(pkg)),
    median_ns = as.numeric(bm$median) * 1e9,
    min_ns = as.numeric(bm$min) * 1e9,
    total_time_ns = as.numeric(bm$total_time) * 1e9,
    n_iter = as.integer(bm$n_itr),
    mem_alloc_bytes = as.numeric(bm$mem_alloc)
  )
}

# A skip rule matches when its package + operation match (with "*" as
# wildcard) and the fixture's n / id fall within optional bounds.
pkg_skip <- function(pkg, op, fx, skip_table = .skip_table) {
  if (length(skip_table) == 0L) {
    return(FALSE)
  }
  for (rule in skip_table) {
    pkg_match <- identical(rule$package, pkg) || identical(rule$package, "*")
    op_match <- identical(rule$operation, op) || identical(rule$operation, "*")
    n_min_ok <- is.null(rule$n_min) || fx$n >= rule$n_min
    n_max_ok <- is.null(rule$n_max) || fx$n <= rule$n_max
    id_ok <- is.null(rule$fixture_id) || identical(rule$fixture_id, fx$id)
    if (pkg_match && op_match && n_min_ok && n_max_ok && id_ok) {
      return(TRUE)
    }
  }
  FALSE
}

# Filter the per-package expression list against the skip table, then call
# bench::mark on the survivors and serialise the result.
run_bench <- function(calls, operation, fx) {
  keep <- !vapply(
    names(calls),
    pkg_skip,
    logical(1),
    op = operation,
    fx = fx
  )
  if (!any(keep)) {
    return(NULL)
  }
  args <- c(
    calls[keep],
    list(check = FALSE, min_iterations = 5L, time_unit = "s")
  )
  bm <- eval(as.call(c(quote(bench::mark), args)))
  mark_to_rows(bm, operation, fx)
}

#| label: data
#| include: false
#| purl: false

benchmarks <- readRDS("data/benchmarks.rds")
metadata <- jsonlite::fromJSON("data/metadata.json", simplifyVector = FALSE)
benchmarks$median_ms <- benchmarks$median_ns / 1e6

This vignette compares the performance of caugi against R packages (igraph, bnlearn, dagitty, ggm, pcalg), the Python package pgmpy, and the Java library Tetrad. We focus on the practical trade-offs that arise from different core data structures and design choices.

The headline is: caugi frontloads computation. That is, caugi spends more effort when constructing a graph and preparing indexes, so that queries are wicked fast. The high performance is also due to the Rust backend.

A note on Tetrad: it is included as a reference point—the de facto gold standard implementation in causal discovery—but the comparison is not strictly apples-to-apples, since Tetrad runs on the JVM (compiled Java) while the R packages and caugi run with R-level dispatch overhead on every call.

Benchmarks below are precomputed, with the last run on r metadata$generated on host r metadata$host$nodename (r metadata$host$sysname).

Versions:

• caugi r metadata$versions$caugi
• igraph r metadata$versions$igraph
• bnlearn r metadata$versions$bnlearn
• dagitty r metadata$versions$dagitty
• ggm r metadata$versions$ggm
• pcalg r metadata$versions$pcalg
• pgmpy r metadata$versions$pgmpy
• Tetrad r metadata$versions$tetrad.

Design Choices

Compressed Sparse Row (CSR) Representation

The core data structure in caugi is a compressed sparse row (CSR) representation of the graph. CSR representations store, for each vertex, a contiguous slice of neighbor IDs with a pointer (offset) array that marks the start and end of each slice. This format is memory-efficient for sparse graphs. The caugi graph object also stores important query information directly, allowing parent, child, and neighbor queries to be done in $\mathcal{O}(1)$ time. This yields a larger memory footprint, but the trade-off is that queries are extremely fast.

Mutation and Lazy Building

The caugi graph objects are expensive to build. This is the performance downside of using caugi. Each time we make a modification to a caugi graph object, we need to rebuild the graph completely since the graph object is immutable by design. This has complexity $\mathcal{O}(|V| + |E|)$, where $V$ is the vertex set and $E$ is the edge set.

However, the graph object will only be rebuilt when the user either calls build() directly or queries the graph. Therefore, you do not need to worry about wasting compute time by iteratively making changes to a caugi graph object, as the graph rebuilds lazily when queried. By doing this, caugi graphs feel mutable, but, in reality, they are not.

By doing it this way, we ensure

• that the graph object is always in a consistent state when queried and
• that queries are as fast as possible,

while keeping the user experience smooth.

Comparison

The comparison covers seven operations across a parameter grid of graph sizes $n$ and edge densities. All packages run on the same generated DAG fixtures so query results are directly comparable.

The grid covers $n \in {100, 1\,000, 10\,000}$ at two density levels, parametrized by the average in-degree $d$ (equivalently the average out-degree, or total edges divided by number of nodes). We sample DAGs with per-edge probability $p = 2d / (n - 1)$, which makes the expected in-degree exactly $d$ regardless of $n$. The grid uses $d = 3$ for the sparser cells and $d = 6$ for the denser ones, so a node has on average $\approx 3$ (or $\approx 6$) parents and $\approx 3$ (or $\approx 6$) children, with edge counts that grow linearly with $n$. The slowest packages (dagitty, ggm, plus bnlearn::dsep and pcalg::dsep) are skipped at $n = 10\,000$ via the skip table in spec.json; the affected lines simply end at $n = 1\,000$.

#| label: helpers
#| include: false
#| purl: false

library(ggplot2)

pkg_levels <- c(
  "caugi",
  "igraph",
  "bnlearn",
  "dagitty",
  "ggm",
  "pcalg",
  "pgmpy",
  "tetrad"
)
benchmarks$package <- factor(benchmarks$package, levels = pkg_levels)

# Stable package → colour mapping. Used in every plot so a package keeps the
# same colour even when it is missing from a panel (e.g. dagitty at n=10000).
# ColorBrewer Set1 + near-black for Tetrad.
pkg_colors <- c(
  caugi = "#E41A1C",
  igraph = "#377EB8",
  bnlearn = "#4DAF4A",
  dagitty = "#984EA3",
  ggm = "#FF7F00",
  pcalg = "#A65628",
  pgmpy = "#F781BF",
  tetrad = "#1B1B1B"
)

plot_op <- function(op) {
  d <- benchmarks[benchmarks$operation == op, ]
  op_manual <- switch(op, d_separated = "d-separation", op)
  op_title <- tools::toTitleCase(gsub("_", " ", op_manual))

  ggplot(d, aes(n, median_ms, color = package)) +
    geom_line() +
    geom_point() +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = pkg_colors, drop = FALSE) +
    facet_wrap(
      ~avg_degree,
      labeller = as_labeller(c(
        "3" = "Sparser (avg degree 3)",
        "6" = "Denser (avg degree 6)"
      ))
    ) +
    labs(
      x = "Number of nodes",
      y = "Median time (ms)",
      color = NULL,
      title = op_title
    ) +
    theme_minimal(base_size = 11) +
    theme(legend.position = "bottom")
}

Relational Queries

Direct neighbor lookups (parents, children) are $O(1)$ for caugi, igraph, bnlearn, and Tetrad---they all maintain adjacency tables. pgmpy is a thin wrapper over NetworkX's predecessor/successor lookups. pcalg and ggm scan a row or column of the adjacency matrix on each call ($O(n)$), with ggm paying noticeably more R-level overhead. dagitty is the outlier: it re-derives neighbors from its string DSL representation each call and pays a linear-in-edges cost.

#| label: bench-parents
#| eval: false
#| echo: true

bench_parents <- function(graphs, fx) {
  v <- fx$test_node
  v_idx <- match(v, graphs$nodes)

  calls <- list(
    caugi = bquote(caugi::parents(graphs$cg, .(v))),
    igraph = bquote(igraph::neighbors(graphs$ig, .(v), mode = "in")),
    bnlearn = bquote(bnlearn::parents(graphs$bn, .(v))),
    dagitty = bquote(dagitty::parents(graphs$dg, .(v))),
    ggm = bquote(ggm::pa(.(v), graphs$am)),
    pcalg = bquote(pcalg::searchAM(graphs$amat_pag, .(v_idx), type = "pa"))
  )

  run_bench(calls, "parents", fx)
}

#| label: plot-parents
#| purl: false

plot_op("parents")

#| label: bench-children
#| eval: false
#| echo: true

bench_children <- function(graphs, fx) {
  v <- fx$test_node
  v_idx <- match(v, graphs$nodes)

  calls <- list(
    caugi = bquote(caugi::children(graphs$cg, .(v))),
    igraph = bquote(igraph::neighbors(graphs$ig, .(v), mode = "out")),
    bnlearn = bquote(bnlearn::children(graphs$bn, .(v))),
    dagitty = bquote(dagitty::children(graphs$dg, .(v))),
    ggm = bquote(ggm::ch(.(v), graphs$am)),
    pcalg = bquote(pcalg::searchAM(graphs$amat_pag, .(v_idx), type = "ch"))
  )

  run_bench(calls, "children", fx)
}

#| label: plot-children
#| purl: false

plot_op("children")

As you can see from the plots, Tetrad is fastest (sub-microsecond), with pgmpy and caugi close behind; all three stay flat as $n$ grows, as do igraph and bnlearn. The cost of pcalg's $O(n)$ row scan, hidden at small $n$, becomes visible by $n = 10000$, where it climbs to a few tenths of a millisecond. ggm and dagitty pay the most for re-deriving neighbors on every call: by $n = 1000$ they are already one to two orders of magnitude slower than the rest, and both are skipped at $n = 10000$ (their lines simply stop). For the front-runners the absolute times are tiny and the gaps between them mostly reflect dispatch overhead differences between Python, R, and Java rather than anything that matters in practice.

Ancestors and Descendants

These require transitive closure (BFS/DFS over the directed graph). Here caugi's CSR pays off---both forward and reverse adjacency are precomputed, so traversal is tight pointer-chasing in Rust.

#| label: bench-ancestors
#| eval: false
#| echo: true

bench_ancestors <- function(graphs, fx) {
  v <- fx$test_node
  v_idx <- match(v, graphs$nodes)

  calls <- list(
    caugi = bquote(caugi::ancestors(graphs$cg, .(v))),
    igraph = bquote(igraph::subcomponent(graphs$ig, .(v), mode = "in")),
    bnlearn = bquote(bnlearn::ancestors(graphs$bn, .(v))),
    dagitty = bquote(dagitty::ancestors(graphs$dg, .(v))),
    pcalg = bquote(pcalg::searchAM(graphs$amat_pag, .(v_idx), type = "an"))
  )

  run_bench(calls, "ancestors", fx)
}

#| label: plot-ancestors
#| purl: false

plot_op("ancestors")

#| label: bench-descendants
#| eval: false
#| echo: true

bench_descendants <- function(graphs, fx) {
  v <- fx$test_node
  v_idx <- match(v, graphs$nodes)

  calls <- list(
    caugi = bquote(caugi::descendants(graphs$cg, .(v))),
    igraph = bquote(igraph::subcomponent(graphs$ig, .(v), mode = "out")),
    bnlearn = bquote(bnlearn::descendants(graphs$bn, .(v))),
    dagitty = bquote(dagitty::descendants(graphs$dg, .(v))),
    pcalg = bquote(pcalg::searchAM(graphs$amat_pag, .(v_idx), type = "de"))
  )

  run_bench(calls, "descendants", fx)
}

#| label: plot-descendants
#| purl: false

plot_op("descendants")

caugi is the most consistent here: sub-millisecond on every fixture for both ancestors and descendants, and barely moving as $n$ grows. igraph::subcomponent() and pgmpy's NetworkX-backed traversal also stay fast across the grid. Descendants on the dense $n = 10000$ fixture is where the design differences blow up: Tetrad takes about ten seconds, while pcalg and bnlearn reach several hundred milliseconds---caugi finishes the same query in roughly a tenth of a millisecond, a five-orders-of-magnitude gap to Tetrad. dagitty is again one to two orders of magnitude slower than the front-runners on the sizes it can handle, paying its DSL-reparsing cost on every call, and is skipped at $n = 10000$.

Markov Blanket

The Markov blanket of a node $v$ is parents(v) ∪ children(v) ∪ parents(children(v)). caugi, bnlearn, and Tetrad have direct accessors backed by adjacency tables; dagitty re-derives it from its DSL representation each call; pgmpy walks NetworkX edges. pcalg and ggm do not expose a Markov-blanket helper and are omitted.

#| label: bench-markov-blanket
#| eval: false
#| echo: true

bench_markov_blanket <- function(graphs, fx) {
  v <- fx$test_node

  calls <- list(
    caugi = bquote(caugi::markov_blanket(graphs$cg, .(v))),
    bnlearn = bquote(bnlearn::mb(graphs$bn, .(v))),
    dagitty = bquote(dagitty::markovBlanket(graphs$dg, .(v)))
  )

  run_bench(calls, "markov_blanket", fx)
}

#| label: plot-markov-blanket
#| purl: false

plot_op("markov_blanket")

caugi, pgmpy, and bnlearn all resolve Markov blankets in microseconds across the entire grid, and Tetrad stays well under a millisecond too. dagitty is the outlier, climbing into the hundreds of milliseconds on the dense $n = 1000$ fixture because it re-derives the blanket from its DSL representation on every call; it is skipped at $n = 10000$.

D-Separation

For each fixture we pick a random (X, Y) pair and compute a minimal d-separator Z with caugi::minimal_separator()---that is, a smallest set of nodes whose conditioning blocks every path between X and Y in the DAG. The same triple (X, Y, Z) is then passed to every package's d-separation routine, so every package answers "is X ⫫ Y | Z?" on identical inputs. Tetrad uses the m-separation generalisation (equivalent on DAGs). pcalg::dsep() follows Lauritzen's moralization-based test on a graphNEL.

#| label: bench-dsep
#| eval: false
#| echo: true

bench_dsep <- function(graphs, fx) {
  if (is.null(fx$dsep)) {
    return(NULL)
  }

  x <- fx$dsep$x
  y <- fx$dsep$y
  z <- unlist(fx$dsep$z)

  calls <- list(
    caugi = bquote(caugi::d_separated(graphs$cg, .(x), .(y), .(z))),
    bnlearn = bquote(bnlearn::dsep(graphs$bn, .(x), .(y), .(z))),
    dagitty = bquote(dagitty::dseparated(graphs$dg, .(x), .(y), .(z))),
    pcalg = bquote(pcalg::dsep(.(x), .(y), .(z), graphs$gNEL))
  )

  run_bench(calls, "d_separated", fx)
}

#| label: plot-dsep
#| purl: false

plot_op("d_separated")

caugi is the clear winner here, staying near a tenth of a millisecond on every fixture, including the dense $n = 10000$ graph. pgmpy is the runner-up, within a few milliseconds throughout. This is the one operation where Tetrad's compiled code does not save it: its m-separation routine explores paths in a way that scales badly on dense graphs, reaching over a second at $n = 1000$ and roughly nine seconds at $n = 10000$---a striking example of an algorithmic choice dominating the language advantage. Among the R packages, bnlearn and pcalg pay to build a moralized graph on each call (pcalg into the hundreds of milliseconds at $n = 1000$, slowest of the three R competitors), and dagitty re-parses its DSL; all three are skipped at $n = 10000$.

Subgraph Extraction

Subgraph extraction is where we explicitly test graph building performance. For caugi we time subgraph(cg, nodes) |> build() together so the CSR rebuild cost is captured. igraph, bnlearn, pgmpy, and Tetrad all return adjacency-backed subgraphs without an analogous index rebuild step. pcalg exposes no native subgraph constructor (users typically call graph::subGraph() on the underlying graphNEL), so it is omitted here.

#| label: bench-subgraph
#| eval: false
#| echo: true

bench_subgraph <- function(graphs, fx) {
  sub <- unlist(fx$subgraph_nodes)

  calls <- list(
    caugi = bquote({
      sg <- caugi::subgraph(graphs$cg, .(sub))
      caugi::build(sg)
    }),
    igraph = bquote(igraph::subgraph(graphs$ig, .(sub))),
    bnlearn = bquote(bnlearn::subgraph(graphs$bn, .(sub)))
  )

  run_bench(calls, "subgraph", fx)
}

#| label: plot-subgraph
#| purl: false

plot_op("subgraph")

caugi and igraph are the only packages that stay in the low-single-digit milliseconds at $n = 10000$ (around 1-2 ms). The two run neck-and-neck: caugi is ahead at $n = 100$ and $n = 1000$, while igraph edges it out on the largest fixture---remarkably close given that caugi pays for a full CSR rebuild on every call. Everyone else is one to three orders of magnitude slower on the largest fixtures, copying nodes and edges into a fresh object each call: pgmpy reaches tens to hundreds of milliseconds, Tetrad several hundred, and bnlearn over two seconds. The takeaway is that even though caugi graphs are nominally expensive to construct, the rebuild after subgraph() is fast enough to beat adjacency-list constructors that rebuild no indexes at all.

Summary

Across the seven operations, caugi either leads or sits in the leading group:

• It matches or beats igraph and bnlearn on the $O(1)$ adjacency lookups (parents, children, Markov blanket).
• It is the fastest R package by a wide margin on transitive-closure operations (ancestors, descendants, d-separation), where the CSR's precomputed forward and reverse adjacency pays off.
• Even on subgraph(), where the CSR has to be rebuilt, it stays in the low-single-digit milliseconds at $n = 10000$ (matched only by igraph), while bnlearn, pgmpy, and Tetrad reach hundreds of milliseconds to seconds.

The most striking pattern emerges at $n = 10000$: caugi is the only package that stays fast across every operation. Tetrad (compiled Java) and pgmpy (NetworkX-backed) beat caugi on the simplest lookups by a small constant factor, which is unsurprising given the R-level dispatch overhead caugi pays on every call. But several competitors either run out of room (dagitty, ggm, and the moralization-based d-separation tests are too slow to include at $n = 10000$) or collapse on specific operations. Tetrad in particular, despite being compiled, takes seconds on dense descendants and d-separation queries because of how those routines scale, not because of the language.

The headline holds: by frontloading work into the build step and querying a CSR-backed structure, caugi keeps query times effectively constant in graph size for the operations users hit most often. The price is paid up front, in the rebuild step that runs lazily on the first query after a mutation.

Reproducing These Benchmarks

The R benchmark code shown in the chunks above is the actual source: each chunk has eval = FALSE so the rendered vignette does not run them, and the wrapper script at tools/benchmark/run_r_bench.R extracts them with knitr::purl() and sources the result into the R session.

The full harness lives in tools/benchmark/ and uses Task to orchestrate the cross-language runners. To reproduce just the R numbers (no Python or Java toolchains needed):

cd tools/benchmark
task fixtures && task r

task r purls the vignette into tools/benchmark/bench_r.R, runs it, prints a per-package, per-operation summary, and writes the raw timings to results/r.csv. Do not edit bench_r.R directly---change the chunks in this vignette instead. See tools/benchmark/README.md for the full pipeline (task all) and per-language details.

#| label: bench-driver
#| include: false
#| eval: false

results <- list()

for (fx in spec$fixtures) {
  message(sprintf("[bench_r] %s (n=%d, edges=%d)", fx$id, fx$n, fx$n_edges))
  graphs <- load_fixture(fx)
  results[[length(results) + 1L]] <- bench_parents(graphs, fx)
  results[[length(results) + 1L]] <- bench_children(graphs, fx)
  results[[length(results) + 1L]] <- bench_ancestors(graphs, fx)
  results[[length(results) + 1L]] <- bench_descendants(graphs, fx)
  results[[length(results) + 1L]] <- bench_markov_blanket(graphs, fx)
  results[[length(results) + 1L]] <- bench_dsep(graphs, fx)
  results[[length(results) + 1L]] <- bench_subgraph(graphs, fx)
}

results <- results[!vapply(results, is.null, logical(1))]
out <- data.table::rbindlist(results)
data.table::fwrite(out, file = file.path(RESULTS_DIR, "r.csv"))

message(sprintf("[bench_r] wrote %d rows to %s/r.csv", nrow(out), RESULTS_DIR))

op_levels <- c(
  "parents",
  "children",
  "ancestors",
  "descendants",
  "markov_blanket",
  "d_separated",
  "subgraph"
)

summary_dt <- out[,
  .(median_ms = median(median_ns) / 1e6),
  by = .(package, operation)
]
summary_dt[, operation := factor(operation, levels = op_levels)]
summary_wide <- data.table::dcast(
  summary_dt,
  package ~ operation,
  value.var = "median_ms"
)

message(
  "\n[bench_r] median ms per (package, operation), aggregated over fixtures:"
)

print(summary_wide, digits = 3, row.names = FALSE)

Session Info

#| label: session-info
#| purl: false

sessionInfo()

caugi documentation built on Oct. 28, 2026, 1:07 a.m.