Benchmarks

Every graph in CausalStructures precomputes and stores each node's parents, children, spouses, and neighbors directly, so parents, children, spouses, and neighbors are effectively $O(1)$ lookups rather than scans over the whole edge set. Many higher-level algorithms, including d_separated/m_separated, minimal_separator, and adjustment-set search, are built on these primitive operations, so their performance depends directly on the efficiency of these basic queries.

using CausalStructures
using BenchmarkTools
using Random

Queries stay fast as graphs grow

Single-hop queries such as children are just slice lookups, so their runtime is essentially unchanged whether the graph has 100 nodes or 10,000:

for n in (100, 1_000, 10_000)
    dag = generate_graph(Random.Xoshiro(1), n; p = 5 / n, class = DAG)
    x = Symbol("V", n ÷ 2)
    t = @benchmark children($dag, $x) samples = 200 evals = 1
    println(rpad("n=$n", 10), "children: ", median(t))
end

n=100     children: TrialEstimate(50.000 ns)
n=1000    children: TrialEstimate(60.000 ns)
n=10000   children: TrialEstimate(50.000 ns)

These primitives keep more involved queries fast too, on graphs with hundreds of nodes and across graph classes:

dag = generate_graph(Random.Xoshiro(1), 500; p = 0.1, class = DAG)
mag = generate_graph(Random.Xoshiro(1), 500; p = 0.01, class = MAG, latents = 20)
x, y = :V10, :V290

rows = [
    ("d_separated (DAG)", () -> d_separated(dag, x, y)),
    ("m_separated (MAG)", () -> m_separated(mag, x, y)),
    ("topological_sort", () -> topological_sort(dag)),
    ("minimal_separator", () -> minimal_separator(dag, x, y)),
    ("markov_blanket", () -> markov_blanket(dag, x)),
]
for (label, f) in rows
    t = @benchmark $f() samples = 200 evals = 1
    println(rpad(label, 20), median(t))
end

d_separated (DAG)   TrialEstimate(56.060 μs)
m_separated (MAG)   TrialEstimate(2.665 μs)
topological_sort    TrialEstimate(45.630 μs)
minimal_separator   TrialEstimate(58.118 μs)
markov_blanket      TrialEstimate(4.308 μs)

Searching for adjustment sets

all_backdoor_sets, all_adjustment_sets, all_iv_sets, and all_frontdoor_sets search subsets of the candidate universe up to max_size, checking each one against the relevant validity criterion. Consequently, the runtime grows exponentially with max_size:

admg = generate_graph(Random.Xoshiro(1), 35; p = 0.25, class = ADMG, latents = 8)
x, y = :V3, :V33

for ms in (2, 3, 4)
    t = @benchmark all_adjustment_sets($admg, $x, $y; minimal = false, max_size = $ms) samples =
        20 evals = 1
    r = all_adjustment_sets(admg, x, y; minimal = false, max_size = ms)
    println("max_size=$ms  ", median(t), "  (", length(r), " sets)")
end

max_size=2  TrialEstimate(2.004 ms)  (0 sets)
max_size=3  TrialEstimate(23.334 ms)  (0 sets)
max_size=4  TrialEstimate(205.259 ms)  (0 sets)

Expensive algorithms

Here we will show the performance of some of the most expensive algorithms.

Exact uniform DAG sampling

generate_graph uses an Erdős–Rényi model: it samples each edge independently, which is cheap regardless of n but is not uniform over the space of DAGs. uniform_dag instead draws exactly uniformly from all labelled DAGs on n nodes using the recursive enumeration algorithm of Kuipers and Moffa (2015).

for n in (10, 20, 40)
    t = @benchmark uniform_dag(Random.Xoshiro(1), $n) samples = 20 evals = 1
    println("n=$n  ", median(t))
end

n=10  TrialEstimate(110.787 μs)
n=20  TrialEstimate(677.750 μs)
n=40  TrialEstimate(5.521 ms)

Markov equivalence class enumeration

count_dags and enumerate_dags operate on a CPDAG, PDAG, or MPDAG by considering every DAG in its Markov equivalence class. The size of that class grows combinatorially with the number of undirected edges. For example, a clique on $k$ nodes has $k!$ consistent orientations, so the class can grow very large very quickly:

prettyresult(t) =
    string(BenchmarkTools.prettytime(t.time), " / ", BenchmarkTools.prettymemory(t.memory))

for k in (5, 7, 9)
    names = [Symbol("V$i") for i = 1:k]
    clique_edges = [undirected(names[i], names[j]) for i = 1:k for j = (i+1):k]
    pdag = PDAG(clique_edges...)
    c = count_dags(pdag)
    tc = @benchmark count_dags($pdag) samples = 5 evals = 1
    te = @benchmark enumerate_dags($pdag) samples = 5 evals = 1
    println(rpad("k=$k", 6), "DAGs=$c")
    println(rpad("  count_dags:", 20), prettyresult(median(tc)))
    println(rpad("  enumerate_dags:", 20), prettyresult(median(te)))
end

k=5   DAGs=120
  count_dags:       569.262 μs / 728.84 KiB
  enumerate_dags:   611.601 μs / 831.38 KiB
k=7   DAGs=5040
  count_dags:       34.579 ms / 38.81 MiB
  enumerate_dags:   36.435 ms / 45.58 MiB
k=9   DAGs=362880
  count_dags:       3.406 s / 3.32 GiB
  enumerate_dags:   5.396 s / 4.04 GiB

enumerate_dags must materialize every DAG rather than simply count them, so it is slower and uses more memory than count_dags, as the numbers above show.

MAG equivalence class enumeration

enumerate_mags is the PAG/MAG counterpart to enumerate_dags. While enumerate_dags is fairly efficient via Chickering's recursive pruning, the algorithm for enumerate_mags is simply just brute-forcing every tail/arrowhead assignment for each circle endpoint in the PAG (2^k candidates for k circle endpoints), and checks if it's a valid PAG. Thus, it's considerably more expensive than count_dags/enumerate_dags:

for k in (3, 4, 5)
    names = [Symbol("V$i") for i = 1:k]
    circle_edges = [partial(names[i], names[j]) for i = 1:k for j = (i+1):k]
    pag = PAG(circle_edges...)
    m = length(enumerate_mags(pag))
    t = @benchmark enumerate_mags($pag) samples = 5 evals = 1
    println("k=$k  MAGs=$m  ", median(t))
end

k=3  MAGs=23  TrialEstimate(704.464 μs)
k=4  MAGs=242  TrialEstimate(60.814 ms)
k=5  MAGs=4457  TrialEstimate(17.915 s)

Comparison to CausalInference.jl

CausalInference.jl is another Julia package implementing some of the same identification criteria (DAGs only): the generalized adjustment criterion, the backdoor criterion, and the frontdoor criterion.

On the tested graphs, CausalStructures enumerates all valid adjustment, backdoor, and frontdoor sets substantially faster than CausalInference.jl. See benchmark/README.md for the full setup and numbers.