Permutation-based search (BOSS / GRaSP)

Download this notebook

BOSS and GRaSP search over orderings of the variables. For an ordering, a BIC score chooses each variable’s parents from the variables before it; the search keeps the ordering with the best score. Both return a CPDAG, like GES, and suit linear-Gaussian data.

Start

Choose a CPU runtime and run this cell first. In Colab it installs andrey-core==0.9.0 from PyPI; local kernels use the installed environment.

import os
import sys

IN_COLAB = "google.colab" in sys.modules

if IN_COLAB:
    os.environ["ANDREY_DEVICE"] = "cpu"
    !pip install --quiet "andrey-core[viz]==0.9.0" pandas

import numpy as np
import pandas as pd
from IPython.display import display

import andrey
import andrey.metrics as metrics
import andrey.viz
from andrey.data import benchmarks

Data

linear_gauss_er draws an Erdos-Renyi graph with linear mechanisms and Gaussian noise: here 10 variables and 2,000 samples.

seed = 7
dataset = benchmarks.scm("linear_gauss_er", d=10).sample(n=2000, seed=seed, scale="standardize")
X, truth = dataset.data, dataset.observed_graph
print(f"{X.shape[0]} observations, {X.shape[1]} variables, seed {seed}")
display(andrey.viz.draw(truth))

2000 observations, 10 variables, seed 7

Fit BOSS and GRaSP

BOSS starts from an ordering computed from the data and visits the variables in a random order; GRaSP starts from a random ordering. seed fixes the random choices, so a run repeats exactly; without seed, the ANDREY_SEED environment variable sets it (default 0). metadata["score"] is the BIC objective of the final ordering, which the search minimizes. A larger lambda_value weighs the BIC penalty more and gives a sparser graph.

results = {
    "BOSS": andrey.boss(X, seed=0),
    "GRaSP": andrey.grasp(X, seed=0),
}
for name, result in results.items():
    print(f"{name}: {result.structure.kind}, score {result.metadata['score']:.1f}")
    display(andrey.viz.draw(result.structure))

BOSS: cpdag, score -22398.8

GRaSP: cpdag, score -22406.2

Compare with GES over ten datasets

GES searches over CPDAGs with a BIC score as well. Structural Hamming distance (SHD) counts errors against the truth’s CPDAG; lower is better. The empty graph’s SHD is the number of true edges, so an estimate at or above it has not improved on returning no edges. One dataset can favor any method, so the table covers ten seeds and their median.

methods = {
    "GES": lambda X: andrey.ges(X),
    "BOSS": lambda X: andrey.boss(X, seed=0),
    "GRaSP": lambda X: andrey.grasp(X, seed=0),
}
rows = {}
for s in range(10):
    draw = benchmarks.scm("linear_gauss_er", d=10).sample(n=2000, seed=s, scale="standardize")
    d = draw.observed_graph.n_nodes
    empty = andrey.GraphStructure.from_numpy(np.zeros((d, d), dtype=np.uint8), kind="cpdag")
    rows[s] = {"Empty graph": metrics.score(empty, draw.observed_graph)["shd"]}
    for name, fit in methods.items():
        rows[s][name] = metrics.score(fit(draw.data), draw.observed_graph)["shd"]
shd = pd.DataFrame(rows)
shd.columns.name = "seed"
shd["median"] = shd.median(axis=1)
shd

seed	0	1	2	3	4	5	6	7	8	9	median
Empty graph	15	15	15	15	15	15	15	15	15	15	15.0
GES	21	10	0	10	17	7	11	17	0	0	10.0
BOSS	6	8	0	10	1	10	13	11	0	14	9.0
GRaSP	6	0	0	1	1	0	1	6	0	0	0.5

A seed where GES reaches the empty graph’s SHD is one where it has not improved on returning no edges. The release cards show the same on larger graphs: GES stays close to the empty graph’s SHD, while BOSS and GRaSP come close to zero.

Use another dataset

Change seed, or the d and n passed to benchmarks.scm and sample, and rerun. The benchmarks page also compares both methods with causal-learn.