andrey.score_matching

andrey.score_matching(data, *, alpha=0.001, knowledge=None)

Learn a DAG with SCORE, a causal order from the Jacobian of the score (additive noise).

SCORE [1] models each variable as a nonlinear function of its parents plus Gaussian noise, x_j = f_j(parents) + e_j. Under that model the score of the joint density, grad log p(x), marks the leaves: the derivative of its entry j in x_j is the same at every point exactly when j has no children [1]. SCORE estimates those derivatives at the rows by a kernel estimator built on Stein’s identity [2], takes the variable whose derivative varies least as the next leaf, repeats on the others, and prunes the complete DAG of the order by the smooth-term tests of CAM [3]. It returns a DAG and a causal order.

Parameters:

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

• alpha (float) – Significance level of the pruning test, in (0, 1): a parent stays when its smooth term’s p-value is at most alpha. The causal order does not read it. 0.001 is the value of the paper’s experiments.

• knowledge (Knowledge | None) – Background knowledge, a Knowledge: required and forbidden edges, the pairs allowed to be adjacent, and temporal tiers, all honored. A variable is the next leaf only when no remaining variable is in a later tier or is the effect of a required edge out of it. Forbidden edges, edges into an earlier tier, and pairs outside the allowed set are in no pruning fit and never in the result; the score estimate still reads every remaining variable. The pruning never removes a required edge.

Returns:

Wraps a dag GraphStructure; metadata["algorithm"] is "score_matching". ordering is the causal order of the leaf search, as a tuple of column indices, causes first. metadata["tie_break"] is the tie-break rule the fit ran under, andrey.config.tie_break.

Return type:

StructureOutput

Raises:

• TypeError – If an argument has the wrong type, such as a string for a number; the message names the method, the argument, what it accepts, and the value. Also if knowledge is neither None nor a Knowledge.

• ValueError – If an argument lies outside its range or is not one of its choices; the message names the method, the argument, what it accepts, and the value. Also if data is not a usable 2-D matrix: not numeric, holding NaN or inf, with fewer than twice as many rows as columns, or with a constant column, a column of fewer than 3 distinct values, or a column that repeats another up to scale and shift, each named in the message; or if more than half of the pairs of rows are equal on the remaining variables at a step, which leaves the kernel no bandwidth. Also if knowledge names a column data lacks or contradicts itself, as Knowledge lists.

Warns:

DataTypeWarning – When a column holds at most 10 distinct whole numbers and so may hold category codes, which the fit reads as numbers. Each message warns once per process.

Notes

When to use. SCORE suits i.i.d. continuous observational data whose effects are nonlinear in every parent, with additive Gaussian noise and no unmeasured common causes [1]; under that model the DAG is identified [4], and it returns the DAG with a causal order. Its order search fits no regression: each step is one kernel estimate over all the remaining variables, so its cost grows with the number of variables only linearly, and with the cube of the rows. The paper shows it on 10 to 50 variables and 1,000 rows. A link that is linear in a parent leaves that parent’s derivative constant too, so it can be taken for a leaf: with linear effects and Gaussian noise use pc() or ges(), which return the Markov equivalence class. For additive noise that is not Gaussian use resit(), and for the same model searched by regressions, cam().

The steps. The paper’s two phases [1]:

1. Causal order. The score’s Jacobian has the diagonal entries d^2 log p(x) / d x_j^2. They are estimated at every row of the remaining variables by the paper’s extension of the Stein gradient estimator [2] to second derivatives (its Eq. 12), with the Gaussian kernel exp(-||x - y||^2 / (2 s^2)), the bandwidth s the median distance between the rows, and one ridge term for both of its regressions. The variable whose entry has the smallest variance over the rows is the next leaf; it is removed, and the estimate, bandwidth included, is made again on the others. The causal order is the leaves in reverse.

2. Pruning. Each variable is fit on all the variables before it in the order, an additive model of penalized thin plate regression splines with ten basis functions per parent and smoothing chosen by GCV [5], with ceil(n_samples / (3 n_parents)) functions per term, at least 3, when there are fewer than 30 rows per term, and at most one per distinct value. A parent stays when its smooth term’s p-value is at most alpha, the pruning of cam() [3].

The choices the paper leaves open. It states no value of the ridge term: it is 0.001, the value in the code published with the paper, added to the kernel matrix as the paper’s Eq. 12 writes it, a matrix with ones on its diagonal. The median distance is taken over the pairs of distinct rows. The variances are compared as estimated, as the paper’s Algorithm 1 states, not divided by the entries’ means.

One kernel reads all the remaining columns, so the estimate depends on their units, which the model does not: a leaf’s entry is constant in any units. The columns are therefore centered and divided by their standard deviations before the order search, so rescaling or shifting a column changes no step; the pruning reads each variable through an affine map to a fixed interval. The paper’s experiments run on the columns as simulated. The fit draws no random numbers and runs on one BLAS thread.

Ties. The variances are compared rounded to 9 decimal places, as hc() compares its score changes per row, so two variances that differ only in their rounding tie. A tie goes to the variable with the lower rank computed from the data, as for hc(); the fit runs on the columns in that rank order, so the result does not depend on the column order unless two columns are identical. The column-order guide says what the rank keeps under a new column order.

Missing values. A NaN in data raises ValueError: every fit reads every row of its variables, so a row with a missing value has no place in it. Drop the incomplete rows first, which suits values missing completely at random, or use pc() or fci(), whose tests delete rows per test.

Limitations. The numbers below were measured on the paper’s simulation with standardized columns, Gaussian noise, and 1,000 rows, on one core, 20 graphs per setting: Erdos-Renyi graphs with as many edges as variables or four times as many, each link one draw of a Gaussian process, at 10, 20, and 50 variables. With as many edges as variables the order put 81%, 73%, and 65% of the true edges in the true direction at 10, 20, and 50 variables, and the DAG was 2.8, 7.8, and 25.0 edges from the true one (SHD), where the empty graph is 10, 20, and 50 away; with four times as many edges the order put 97%, 90%, and 85% in the true direction, and the SHD was 22.1, 51.9, and 142.7 against the empty graph’s 40, 80, and 200. On the graphs with as many edges as variables cam() was more accurate on the same data: its order put 96% to 98% of the true edges in the true direction and its SHD was 0.7, 2.1, and 5.5. At 50 variables with four times as many edges, over 40 graphs, the order put 0.007 to 0.010 less of the true edges in the true direction than the code published with the paper did (86% against 87%). The paper’s tables are on data as simulated, and the order search reads standardized columns, where it orders fewer edges of a sparse graph correctly: in a separate check on ten graphs of 10 variables with as many edges, the same search put 99% of the true edges in the true direction when it read the columns as simulated and 80% on standardized columns. On standardized columns the error of the search’s estimate for a variable grows with about the square of that variable’s variance-to-noise ratio (its variance over the variance of its noise term). Where links have a strong linear part, variance builds up along the causal order, the sinks have the largest ratios, and the search ranks them last: on additive-noise data with links z + sin z (1,000 rows, 10 and 20 variables, as many edges as variables), where a sink’s ratio had a median of 10 to 13 against 1.6 to 1.7 in the paper’s simulation, the order put 41% to 53% of the true edges in the true direction, a random order’s share, and 16,000 rows (52%), another ridge term, or another bandwidth did not change that. Use cam() on such data (76% to 94% on the same datasets). One kernel bandwidth serves every variable, and links it does not resolve leave a cause’s entry as flat as a leaf’s: on x0 -> x2 <- x1, x2 -> x3 with x2 = 2 x0 / (1 + x0^2) + 2 / (1 + x1^2), x3 = x2^2 / 2, and noise of standard deviation 0.7 or 1, over 400 and 800 rows, the search took x0 and x1 for the leaves in 24 of 24 fits. The leaf rule is exact only for Gaussian noise: under other noise a leaf’s entry changes with its residual, and nogam() and resit() drop that assumption. A fit took 1.0 s at 10 variables, 3.0 s at 20, and 28 to 32 s at 50, within 0.16 GB. Each step factors a matrix with one row and one column per sample, so the order search’s time grows with the cube of n_samples and its memory with the square. The pruning starts from the complete DAG of the order, so its last fits have nearly n_variables smooth terms each: it needs at least twice as many rows as columns, and it takes most of a fit’s time as the variables grow. Binary columns, and any column with fewer than 3 distinct values, are not supported: the pruning fits a smooth term of each parent, and such a column raises ValueError before the fit. The numbers are a reference, not a hard limit: they change with the data generator, the density, the method’s settings, and the machine.

References

[1] (1,2,3,4)

Rolland, P., Cevher, V., Kleindessner, M., Russell, C., Janzing, D., Scholkopf, B., and Locatello, F. (2022). Score matching enables causal discovery of nonlinear additive noise models. ICML, PMLR 162, 18741-18753.

[2] (1,2)

Li, Y., and Turner, R. E. (2018). Gradient estimators for implicit models. ICLR.

[3] (1,2)

Buhlmann, P., Peters, J., and Ernest, J. (2014). CAM: Causal additive models, high-dimensional order search and penalized regression. Annals of Statistics 42(6), 2526-2556.

[4]

Peters, J., Mooij, J. M., Janzing, D., and Scholkopf, B. (2014). Causal discovery with continuous additive noise models. Journal of Machine Learning Research 15, 2009-2053.

[5]

Wood, S. N. (2003). Thin plate regression splines. Journal of the Royal Statistical Society B 65(1), 95-114.

Examples

>>> import numpy as np
>>> from andrey import score_matching
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(300)
>>> y = np.sin(2 * x) + 0.3 * rng.standard_normal(300)
>>> out = score_matching(np.column_stack([x, y]))
>>> out.structure.kind, out.ordering
('dag', (0, 1))