andrey.calm

andrey.calm(data, *, seed=None, **params)

Learn a weighted DAG with CALM, a continuous-optimization structure learner.

CALM (Continuous and Acyclicity-constrained L0-penalized likelihood with an estimated Moral graph) [1] fits a linear-Gaussian DAG with torch. By default it keeps only the candidate edges of an estimated moral graph, then minimizes the Gaussian likelihood plus a penalty per edge under an acyclicity constraint. The optimizer draws random numbers and is seeded on every call; see seed.

Parameters:

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

• seed (int | None) – Seed for the optimizer’s torch random number generator. None uses the ANDREY_SEED environment variable, or 0 when it is unset, so every call is seeded and repeated calls return the same graph on one machine with the same torch version. The draws go to the columns in order of each column’s rank, as in hc() (the column-order guide). metadata["seed"] records the seed used.

• **params (Any) –

  Optimizer settings passed to the solver; a name not listed here raises TypeError.

  • lambda1 (default 0.005): penalty per edge. Larger values give sparser graphs.

  • alpha (default 0.01): significance level of the Fisher-Z tests that estimate the moral graph. Only pairs it keeps can become edges; smaller values keep fewer.

  • use_moral_graph (default True): False makes every pair a candidate edge.

  • standardize (default False): scale each column to unit variance first.

  • w_threshold (default 0.3) and gate_threshold (default 0.5): an edge is kept when its weight magnitude and its gate probability reach these values.

  • max_outer (default 25), subproblem_iter (default 3000), and lr (default 0.01): quadratic-penalty rounds, Adam steps per round, and learning rate.

  • tau (default 0.5), rho_init (default 0.001), rho_mult (default 3.0), rho_max (default 1e16), and h_tol (default 1e-8): gate temperature and the acyclicity-penalty schedule.

  • labels (default None): node names; a DataFrame’s column names take precedence.

Returns:

Wraps a dag GraphStructure; weighted_adjacency[i, j] is the coefficient of edge i -> j. metadata holds "algorithm" ("CALM"), the "seed" used, "n_edges", "objective" (half the sum of the log residual variances plus lambda1 times the edge count; lower is better), and "h_final" (the acyclicity residual of the returned graph, 0 when it is acyclic).

Return type:

StructureOutput

Raises:

• ImportError – If the [torch] extra is not installed.

• TypeError – If params holds a name the solver does not accept.

• ValueError – If the fit’s loss or weights stop being finite, for example when the data’s scale overflows its covariance; standardize=True rescales the columns.

Notes

CALM first estimates the moral graph with IAMB [2]: for each variable it adds, one at a time, the variable most dependent on it given those already added (Fisher-Z at level alpha), then drops any member the others make independent of it. The most dependent variable has the largest absolute partial correlation; equal ones are chosen by each column’s rank, as in hc(). A pair is a candidate edge, in either direction, when each is in the other’s set.

Each candidate edge i -> j has a weight and a gate. During the fit a gate is a Gumbel-Softmax sample [3], sigmoid((logit + noise) / tau), so lambda1 times the sum of the gates approximates a penalty per edge. With B the gated weights and S the sample covariance, the fit minimizes the Gaussian negative log-likelihood with one noise variance per variable [4], 0.5 * sum_j log(((I - B).T @ S @ (I - B))[j, j]) - log|det(I - B)|, plus that penalty and 0.5 * rho * h**2. Here h = tr((I + G / d)**d) - d [5] for the gates G and d variables, and h is zero exactly when the gated graph is acyclic. Adam runs up to subproblem_iter steps per round; a round ends early once, after 50 steps, the loss changes by less than a 1e-7 fraction from one step to the next. Between rounds rho grows by the factor rho_mult, the quadratic penalty method [6], until h of the noise-free gates sigmoid(logit / tau) is at most h_tol, rho reaches rho_max, or max_outer rounds have run.

If the edges that gate_threshold and w_threshold keep still form a cycle, edges are removed in order of increasing weight magnitude until none does.

The defaults are a shorter schedule than the paper’s [1], which starts rho at 1e-5, runs 40,000 steps per round at learning rate 0.001, and keeps every edge whose gate reaches 0.5. rho_init=1e-5, subproblem_iter=40000, lr=0.001, w_threshold=0.0 with a larger max_outer runs that schedule.

Limitations. Measured on simulated linear-Gaussian data from sparse random graphs (mean degree 2), from 10 to 50 variables with ten rows per variable. Its median errors (SHD) stayed below the empty graph’s at every size (15 against 50 at 50 variables). The median fit took 18 s at 10 variables and 59 s at 50, on 4 cores. The paper’s schedule (lr=0.001, max_outer=100, subproblem_iter=40000, rho_init=1e-5, w_threshold=0.0) took 174 s and 472 s for median errors of 6 and 12. Every optimizer step computes a determinant and a matrix power of the full weight matrix, one row and one column per variable, so a step’s cost grows with the cube of the variable count; the defaults allow up to 75,000 steps (max_outer times subproblem_iter). The fit is nonconvex and can stop at a local optimum: on a collider whose two parents are independent, it can return both edges reversed and an edge between the parents, with the defaults and with the paper’s schedule alike.

References

[1] (1,2)

Jin, K., Ng, I., Zhang, K., and Huang, B. (2026). Revisiting differentiable structure learning: inconsistency of l1 penalty and beyond. AAAI 40(27), 22399-22407. https://doi.org/10.1609/aaai.v40i27.39398

[2]

Tsamardinos, I., Aliferis, C. F., and Statnikov, A. (2003). Algorithms for large scale Markov blanket discovery. FLAIRS 2003, 376-380.

[3]

Jang, E., Gu, S., and Poole, B. (2017). Categorical reparameterization with Gumbel-Softmax. ICLR 2017.

[4]

Ng, I., Ghassami, A., and Zhang, K. (2020). On the role of sparsity and DAG constraints for learning linear DAGs. NeurIPS 33, 17943-17954.

[5]

Yu, Y., Chen, J., Gao, T., and Yu, M. (2019). DAG-GNN: DAG structure learning with graph neural networks. ICML 2019, PMLR 97, 7154-7163.

[6]

Ng, I., Lachapelle, S., Ke, N. R., Lacoste-Julien, S., and Zhang, K. (2022). On the convergence of continuous constrained optimization for structure learning. AISTATS 2022, PMLR 151, 8176-8198.

Examples

>>> import numpy as np
>>> from andrey import calm
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(500)
>>> y = x + 0.1 * rng.standard_normal(500)
>>> out = calm(np.column_stack([x, y]), seed=0)
>>> out.structure.kind
'dag'