The boundary

Shu Wan · October 27, 2026 · Andrey 0.10.0

• Release
• Speed

In 0.9.0, the rocket lifted off. Every method we accelerated produced the same graphs, only faster.

Now comes chapter two. How far can that speed take us?

Every method has its limits. Beyond a certain point, it either runs out of time or memory, or produces a graph no better than the empty graph—a graph with no edges.

In 0.10.0, we chart these limits for every supported method across three dimensions—the number of variables, the number of samples, and graph density.

The results push those boundaries further than ever. BOSS scales to 1,600 variables, GRaSP to 800, PC to 400, and GES to 400 with 10,000 samples, all producing graphs with fewer errors than the empty-graph baseline.

TL;DR

• A map for each method: the problems it finishes within the time cap and solves better than the empty graph.

• How far each family goes, at 2,000 samples and a mean degree of 2:

  • Permutation-based: BOSS to 1,600 variables, GRaSP to 800.

  • Constraint-based: PC to 400.

  • Score-based: GES, GIES, and XGES to 200, and to 400 with 10,000 samples; CALM to 400.

  • Linear non-Gaussian: DirectLiNGAM to 400.

• More samples move an accuracy boundary to the right; a denser graph moves it to the left.

• A boundary is a reference, not a hard limit: it moves with the data, the settings, and the machine.

Two ways to stop

It does not finish. The work of a search grows with the number of variables and with how many causes each one has. Past some size, a fit no longer finishes within our time cap of one hour, or it runs out of memory.

It finishes, but the graph is no better than nothing. The empty graph is the simplest answer: it draws no edges, so it misses every true edge and adds no false one. A method that makes more errors than the empty graph has added more mistakes than it found edges. This happens when there are too few samples for the tests or scores to tell an edge from noise, or when the graph is so dense that each variable has more causes than the data can separate.

The second way matters as much as the first. A method that finishes in a second on 3,000 variables is no help if its graph is worse than drawing none. So a problem is inside the boundary only when the method does both: finishes, and beats the empty graph.

Why the boundary is there

Most of the time goes to conditional independence tests. To decide whether two variables share an edge, PC tests them given sets of their neighbors, first one at a time, then in pairs, triples, and so on. The number of pairs grows with the square of the number of variables, and the number of conditioning sets grows exponentially with the number of neighbors. So a denser graph costs far more tests, and a larger one more again. The score-based and permutation-based searches face the same growth in the candidate parent sets they score.

Too many tests also cost accuracy. A test given a large conditioning set has few samples for each degree of freedom, so it loses power, and each wrong answer removes or keeps an edge that later steps build on. This is why a dense graph can stop a method on accuracy before it stops on time.

The boundary of each method

Each method gets two maps. On the size map, the number of variables grows to the right and the number of samples grows upward, with each variable having two causes or effects on average. On the density map, the samples stay at 2,000 and the graph grows denser upward. Here are BOSS’s:

BOSS

How to read it. Each cell is one problem: a number of variables, with a number of samples or a density. A blue cell is one the method finishes with fewer errors than the empty graph: the darker the blue, the fewer errors. A gray cell finishes, but does no better than the empty graph. A mark is where a row stops: a cross at the one-hour time cap, a square at the memory cap, an exclamation mark when the fit fails, and a circle past the accuracy boundary, after two sizes in a row that do not beat the empty graph. The outline is the boundary: in each row, every size up to the first that does not finish or does not beat the empty graph. Hover a cell for its errors as a share of the empty graph’s.

So BOSS is inside its boundary up to 1,600 variables with 500 or 2,000 samples. With 10,000 samples, each fit takes longer, and at 1,600 variables the fits reach the time cap, so its boundary ends at 800. On the density map, a mean degree of 4 ends at 800 on time. At a mean degree of 8, BOSS does no better than the empty graph at 20 variables, so that row has no boundary, and at 800 variables it runs out of memory.

Every method’s maps:

Constraint-based

CD-NOD

PC

Score-based

CALM

Exact search

GES

GIES

Hill climbing

Tabu search

XGES

Permutation-based

BOSS

GRaSP

Linear non-Gaussian

DirectLiNGAM

ICA-LiNGAM

MultiGroupDirectLiNGAM

The maps show two kinds of boundary. Most score-based searches stop on accuracy: with more samples they reach further, GES from 50 variables at 500 samples to 400 at 10,000, while on a graph with a mean degree of 8 GES, XGES, tabu search, and hill climbing do no better than the empty graph even at 20 variables. The permutation-based searches, exact search, and CALM stop on time, as BOSS does. PC stops on accuracy past 400 variables at every sample size. On denser graphs PC beats the empty graph further, but by less, because the empty graph then misses many more edges; hover a cell for its share.

How far each method goes

How far each method goes

Max variables at		500 samples	2,000 samples	10,000 samples
Constraint-based	PC	400	400	400
	CD-NOD	200	100	100
Score-based	GES	50	200	400
	GIES	100	200	400
	XGES	50	200	400
	Hill climbing	20	100	200
	Tabu search	50	100	200
	Exact search	22	22	22
	CALM	400	400	400
Permutation-based	BOSS	1,600	1,600	800
	GRaSP	800	800	800
Linear non-Gaussian	DirectLiNGAM	200	400	200
	ICA-LiNGAM	20	50	100

The largest number of variables inside the boundary, with a mean degree of 2: every size up to it finished and beat the empty graph.

The 0.10.0 release card gives every size, the conditions, and the methods measured on their own axes: FCI and GFCI, the latent-variable, time-series, and pairwise methods.

Check your own problem

Your data are not our simulations, so the surest way to find your own boundary is to measure it. On a simulation the size and density of your problem, compare a method’s errors with the empty graph’s. Below, GES on 50 variables and 2,000 samples makes a third of the empty graph’s errors on a sparse graph, and three times as many on a dense one, as its density map shows.

example.py

import numpy as np

import andrey
from andrey.data import sample_scm
from andrey.metrics import shd

# GES on 50 variables and 2,000 samples: a sparse graph, then a dense one.
for degree in (2, 8):
    data = sample_scm(d=50, n=2000, density=degree, scale="standardize", seed=0)
    errors = shd(andrey.ges(data.data).structure, data.graph)
    empty = shd(np.zeros_like(data.graph.adjacency), data.graph)
    print(f"mean degree {degree}: {errors} errors; the empty graph makes {empty}")

example.py

mean degree 2: 16 errors; the empty graph makes 50
mean degree 8: 657 errors; the empty graph makes 200

And beyond

0.9.0 made each test and score cheaper, and 0.10.0 shows how far that goes. To go much further, a method has to run fewer tests, and to use the sparsity of real graphs. Both are active research:

• Fewer tests. Shiragur, Zhang, and Uhler (2024) learn a coarser causal graph with a polynomial number of tests, one that equals the true graph when the true graph is identifiable. Kim and Lee (2026) deduce the answers of many high-order tests from lower-order ones instead of running them, so a search runs fewer tests and fewer of the weakest ones.

• Sparsity. PC stays feasible on high-dimensional data when the graph is sparse (Kalisch and Bühlmann, 2007): when each variable has few neighbors, the number of tests grows polynomially, not exponentially, with the number of variables. Spread over the cores of one machine, such a search reaches thousands of variables (Le et al.).

Very large causal discovery, sparsity used to reach millions of variables, is a stage of our roadmap. The boundary in this post is where that work starts.

Limits

• A reference, not a hard limit. A boundary holds for what we measured: our simulated graphs and data (Erdos-Renyi graphs with linear relations, standardized), the sizes and densities above, each method’s settings, and our machine. Most methods run at their defaults; BOSS and GRaSP run with a BIC penalty of 4, and GES with a limit of 4 parents. The release card lists every setting. Another graph structure, other data, other settings (a smaller alpha, a max_depth), or more cores can move it either way. Use it to know where to start, then measure on your own data.

• Time. The cap is one hour per fit. A longer cap moves the time boundary to the right; it does not move the accuracy boundary.

• Other axes. FCI and GFCI, whose samples grow with the variables, and the latent-variable, time-series, and pairwise methods are measured on their own axes; the release card lists them.

Try it

pip install -U andrey-core

Help us build it

Andrey is in alpha. The API may change between releases. When you try it on your own data, tell us what breaks or which method you need. We also welcome pull requests. The contributing guide explains how you can set up the development environment.

---

• Getting started

• Release cards

• Changelog

• GitHub

Nonlinear dependence   Ready for real data