andrey.grasp

andrey.grasp(data, *, score_func=None, lambda_value=1.0, depth=3, seed=None, categorical=None, max_categories=10, knowledge=None)

Learn a CPDAG with GRaSP (Greedy Relaxation of the Sparsest Permutation).

GRaSP searches variable orders (permutations) using linear-Gaussian BIC (Bayesian information criterion). It converts the best-scoring order to a CPDAG, which represents an equivalence class of graphs. It searches the same space as boss() using chains of up to depth tucks [1], the moves described in Notes.

Parameters:

• data (ArrayLike) – Observed data, one row per sample. Coerced to float64. A DataFrame’s column names become the structure’s labels.

• score_func (str | None) – Local score objective: "local_score_BIC_from_cov", the covariance-based linear-Gaussian BIC; "local_score_BDeu" and "local_score_BIC_categorical" are the BDeu score and the multinomial BIC for categorical data coded by whole numbers, and "local_score_degen_gauss" and "local_score_cond_gauss" are the degenerate and conditional Gaussian BICs for mixed data, as for ges(). None reads the column types from categorical, or infers them, and picks "local_score_BDeu" when every column is categorical or "local_score_cond_gauss" when some are, each with a DataTypeWarning; otherwise it picks "local_score_BIC_from_cov". Pass "local_score_BIC_from_cov" to read whole numbers, such as counts, as numbers.

• lambda_value (float) – Multiplier on the BIC penalty of log(n_samples) per parameter: per parent for the linear-Gaussian BIC, per free parameter of the conditional probability table for the multinomial BIC. Larger values give sparser graphs; 1.0 is the standard BIC. Must be finite and non-negative. BDeu ignores it.

• depth (int) – Maximum number of tucks in one chain of moves, at least 1. The depth-first search follows a chain of tucks that leave the score unchanged until one improves it or the chain reaches depth tucks. 1 tries single tucks only; larger depths search more chains, which can reach a better score, and run longer.

• seed (int | None) – Seed for the random starting order and the order in which the search visits variables. None uses the ANDREY_SEED environment variable, or 0 when it is unset, so repeated calls return the same graph. Different seeds can reach different graphs.

• categorical (Sequence[str | int] | None) – The categorical columns of data, for a mixed-data score: column labels (strings) or 0-based positions. None infers them: a DataFrame’s category, bool, and text columns, and every column of at most max_categories distinct whole numbers, named in a DataTypeWarning. [] marks none. The other scores take no column types and reject it. See Data types in the getting-started guide.

• max_categories (int) – The most distinct whole numbers a column may hold to be inferred as categorical, at least 2. With a continuous score, columns within it get a one-time DataTypeWarning that they may be category codes. With score_func=None, the columns it infers pick the default score.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. Each variable’s parents exclude the forbidden ones and include the required ones, and every causal order the search visits puts the cause of each required edge, and each earlier tier, first. The result directs every edge the knowledge decides, then Meek’s rules direct the edges that follow, as Knowledge describes.

Returns:

Wraps a cpdag GraphStructure; metadata["algorithm"] is "GRaSP" and metadata["score"] the score of the DAG read off the final order; for the BIC, the sum over variables of n * log(residual variance) + lambda_value * log(n) * n_parents, for BDeu, -2 times the log of the structure prior times the marginal likelihood, and for the multinomial BIC, -2 times the log-likelihood plus the penalty; lower is better. Directed marks (->) are oriented edges, undirected (--) the ones the equivalence class leaves open. Under a mixed-data setting, metadata["categorical"] lists the categorical columns it read, by name or 0-based position.

Return type:

StructureOutput

Raises:

• NotImplementedError – If score_func names no available score.

• TypeError – If depth is not an int, categorical is not a list of column names and positions, or max_categories is not an int. Also if knowledge is neither None nor a Knowledge.

• ValueError – If depth is less than 1, lambda_value is negative or non-finite, or if data or categorical does not suit the score, as for ges(). With the BIC score, also if a column repeats another up to scale and shift given more than two rows or, with more rows than columns, is a linear function of other columns. One message names every group of copies and says to keep one column per group. Also if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

Warns:

DataTypeWarning – When score_func is None and the columns are all categorical or mixed, so the search scores with BDeu or the conditional Gaussian BIC; when a mixed-data score infers the categorical columns, or finds none and so runs as the linear-Gaussian form; or when a continuous setting meets columns that look like category codes. Each message warns once per process.

See also

andrey.boss

Searches the same space of orders with different moves; faster on the published benchmarks at 400 variables.

andrey.ges

Searches over equivalence classes for the same kind of data.

andrey.fci

Allows unmeasured common causes and selection bias.

andrey.direct_lingam

Orients every edge when the noise is non-Gaussian.

Notes

When to use it. GRaSP suits the same data as ges(): continuous data with linear effects, Gaussian noise, and no unmeasured common causes. Its authors report accurate search even on dense graphs [1]. When every column is categorical, the default score is "local_score_BDeu", and on mixed data "local_score_cond_gauss"; "local_score_BIC_categorical" and "local_score_degen_gauss" also score those data.

When not to. When unmeasured common causes are plausible, use fci(). When the noise is non-Gaussian and every edge should be oriented, use direct_lingam().

GRaSP or BOSS. Both search the same space of orders with different moves. Either can reach the better score on a given dataset. On the published benchmarks, GRaSP took far longer than boss() at 400 variables.

Choosing seed. Both the starting order and the order of visiting variables are random, so different seeds can reach different graphs; comparing a few shows whether a result depends on the seed.

How it works. GRaSP searches variable orders for the best DAG with the lowest total score [1]. Each variable uses a grow-shrink tree [2] to cache its best parent set for any set of variables allowed to precede it. The search starts from a shuffled order and changes it with tucks. A tuck on an edge x -> y moves y just before x, along with its ancestors that sit between x and y. A depth-first search keeps the first chain of up to depth tucks that improves the total score. It repeats until no chain improves it. The final parent sets form a DAG, which is reduced to its CPDAG.

Tucks must improve the score by more than 1e-6, the fixed absolute tolerance used by boss(). A tuck within this tolerance counts as leaving the score unchanged. Under the linear-Gaussian BIC each comparison is decided as full local scores decide it, as for boss().

Column order. Each column gets a tie rank computed from its own data, as for hc(), and the search runs on the columns sorted by it: the seed shuffles the sorted columns into the starting order and the order of visiting variables, and ties between candidate parents go to the lowest rank. The column-order guide says what the rank keeps under a new column order and on another machine.

Categorical data. "local_score_BDeu" reads every column as categorical, each distinct value one category, and scores a variable given its parents with the BDeu score of ges(): a Dirichlet prior with an equivalent sample size of 1 and a structure prior of one expected parent per variable, on the scale of the BIC. "local_score_BIC_categorical" reads the columns the same way and scores a variable with the multinomial BIC of ges(): -2 times the log-likelihood of its counts in each configuration of its parents, plus lambda_value * log(n) for each free parameter. Under either score each column gets its tie rank as a categorical column, and each candidate parent set is scored by full local scores.

Mixed data. "local_score_degen_gauss" [3] and "local_score_cond_gauss" [4] score data with continuous and categorical columns, as ges() describes, and read the categorical columns as pc() does. Each candidate parent set is scored by full local scores, without the Schur update the BIC search uses.

Limitations. Measured on simulated linear-Gaussian data from sparse random graphs (mean degree 2), from 20 to 400 variables with ten rows per variable, at lambda_value=4.0: its median errors (SHD) stayed below the empty graph’s at every size, 0 against 400 at 400 variables. The median fit took 0.037 s at 20 variables, 40 s at 200, and 1,040 s at 400, on 16 cores; larger graphs were not measured. On the 223-variable ANDES network with 1,000 simulated rows, at lambda_value=1.0, its median errors passed the empty graph’s (471 against 338). With fewer rows than columns the BIC score is unreliable. A parent set of n_samples - 1 or more variables fits its child exactly, so the search never selects one, nor any parent set that leaves its child a residual variance of at most 1e-12 times its variance. Against exact_search() on small simulated categorical networks, with BDeu or the multinomial BIC, GRaSP stopped at a DAG that scores worse than the optimum on some draws; comparing a few seeds shows whether a result depends on the seed. Neither categorical score is benchmarked with GRaSP yet. The mixed-data scores are experimental and not benchmarked yet; their tests run GRaSP on 5 variables with 3,000 rows.

References

[1] (1,2,3)

Lam, W.-Y., Andrews, B., and Ramsey, J. (2022). Greedy Relaxations of the Sparsest Permutation Algorithm. UAI, PMLR 180:1052-1062. https://proceedings.mlr.press/v180/lam22a.html

[2]

Andrews, B., Ramsey, J., Sanchez Romero, R., Camchong, J., and Kummerfeld, E. (2023). Fast Scalable and Accurate Discovery of DAGs Using the Best Order Score Search and Grow Shrink Trees. NeurIPS 36. https://doi.org/10.52202/075280-2794

[3]

Andrews, B., Ramsey, J., and Cooper, G. F. (2019). Learning high-dimensional directed acyclic graphs with mixed data-types. Proceedings of Machine Learning Research 104, 4-21.

[4]

Andrews, B., Ramsey, J., and Cooper, G. F. (2018). Scoring Bayesian networks of mixed variables. International Journal of Data Science and Analytics 6(1), 3-18.

Examples

>>> import numpy as np
>>> import pandas as pd
>>> from andrey import grasp
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(500)
>>> y = rng.standard_normal(500)
>>> z = x + y + 0.5 * rng.standard_normal(500)
>>> grasp(pd.DataFrame({"x": x, "y": y, "z": z})).structure.oriented_edges()
[('x', 'z', 'directed'), ('y', 'z', 'directed')]