Stable Differentiable Causal Discovery

Achille Nazaret * 1 2 Justin Hong * 1 2 Elham Azizi 1 2 3 David Blei 1

Abstract

methods are intractable, even for modest numbers of variables (Chickering, 1996). Yet datasets in fields like biology
routinely involve thousands of variables (Dixit et al., 2016).

Inferring causal relationships as directed acyclic
graphs (DAGs) is an important but challenging
problem. Differentiable Causal Discovery (DCD)
is a promising approach to this problem, framing the search as a continuous optimization. But
existing DCD methods are numerically unstable,with poor performance beyond tens of variables.
In this paper, we propose Stable Differentiable
Causal Discovery (SDCD), a new method that improves previous DCD methods in two ways: (1)It employs an alternative constraint for acyclicity;
this constraint is more stable, both theoreticallyand empirically, and fast to compute. (2) It uses
a training procedure tailored for sparse causal
graphs, which are common in real-world scenarios. We first derive SDCD and prove its stability
and correctness. We then evaluate it with observational and interventional data and in both small
and large scale settings. We find SDCD outperforms existing methods in convergence speed andaccuracy, and can scale to thousands of variables.

To address this problem, Zheng et al. (2018) introduced Differentiable Causal Discovery (DCD), which formulates the
DAG search as a continuous optimization over the space ofgraph adjacency matrices. A key element of this strategy isan acyclicity constraint, in the form of a penalty, that guides
an otherwise unconstrained search toward acyclic graphs.
This optimization-oriented formulation often scales better
than previous methods, and it has opened opportunities toharness neural networks (Lachapelle et al., 2019), incorporate interventional data (Brouillard et al., 2020), and use
matrix approximation techniques (Lopez et al., 2022). But,
while promising, existing DCD methods still struggle toscale consistently beyond tens of variables, or they rely on
approximations that limit their applicability (see Section 4).
In this paper, we study the problems of DCD and improveon it, so that it can scale more easily and apply to many
types of causal discovery problems. We trace the issues
with DCD to the instability of its objective function; inparticular, properties of the acyclicity constraint it uses tofind a DAG solution. We formalize this notion of stability,
show that previous DCD methods are unstable, and thenformulate a method, Stable Differentiable Causal Discovery(SDCD), that is stable and scalable.

1. Introduction
Inferring cause-and-effect relations between variables is afundamental challenge in many scientific fields, including
biology (Sachs et al., 2005), climate science (Zhang et al.,2011), and economics (Hoover, 2006). Mathematically, a setof causal relations can be represented with a directed acyclic
graph (DAG) where nodes are variables, and directed edges
indicate direct causal effects. The goal of causal discovery
is to recover the graph from the observed data.

The main contributions of this work are:
• We develop a theoretical analysis of the acyclicity constraints used in DCD, and their numerical instabilities.
• We motivate an alternative acyclicity constraint withsuperior stability, both theoretically and empirically.
• We propose the SDCD method for efficient DCD. It
leverages the stable constraint within a two-stage optimization procedure designed for training robustness.
We prove that SDCD does not compromise accuracy.

The challenge of causal discovery is that searching for thetrue DAG underlying the data is an NP-hard problem. Exact*

Equal contribution 1 Department of Computer Science,
Columbia University, New York, USA 2 Irving Institute for Cancer
Dynamics, Columbia University, New York, USA 3 Department of
Biomedical Engineering, Columbia University, New York, USA.
Correspondence to: Elham Azizi <ea2690@columbia.edu>.

• We empirically study SDCD, and show that it efficiently solves problems with thousands of variables.

2. Stable Acyclicity Constraint

Published at the 2 nd Differentiable Almost Everything Workshop
at the 41 st International Conference on Machine Learning, Vienna, Austria. July 2024. Copyright 2024 by the author(s).

In this section, we demonstrate how most existing acyclicity
constraints can lead to unstable numerical behaviors during
1Stable Differentiable Causal Discovery

optimization, especially with large numbers of variables d.
We then motivate an alternative constraint, which we show
to be theoretically and empirically more stable. We briefly
introduce differentiable discovery and review its standard
notations in Appendix B.

Name

ak

fa

ha

∇h⊤
a

hexp
hlog
hinv
hbinom

1/k!
1/k
1

d

exp(x) − 1
1
log 1−xk(1 + x)d − 1

Tr exp(A) − d− log det(I − A)
Tr(I − A)−1
Tr(I + A)d − dexp(A)(I − A)−1(I − A)−2d(I + A)d−1––|λd (A)|⊤v d u⊤
d /vd ud

hρ

1
1−x

2.1. Power Series Trace Constraints
We first introduce a family of constraints. It generalizes
existing constraints and reveals their similarities.

Table 1: (Top) Existing PST constraints with their power
series and gradients. (Bottom) The spectral acyclicity constraint, which is not PST.

Definition 2.1 (The Power Series Trace Family). For any∗
non-negative coefficients (ak )k∈N∗ ∈ RN
≥0 , consider the∞P
ak xk .
power series fa (x) =Definition 2.3. An acyclicity constraint h is stable if thesethree criteria hold for almost every A ∈ Rd×d
≥0 :

k=1

- E-stable. h(sA) = Os→∞ (s)- V-stable. h(A) ̸= 0 ⇒ h(εA) = Ωε→0+ (ε)- D-stable. h(A) and ∇h(A) are well-defined.

Then, for any matrix A ∈ Rd×d
≥0 with non-negative entries,
we define the Power Series Trace (PST) function
ha (A) = Tr [fa (A)] =

∞
X

k=1

E-stability ensures that h does not explode to infinity; Vstability ensures h does not vanish rapidly to 0; D-stability
ensures that h and its gradient exist.

 
ak Tr Ak .

The quantity ha (A) is closely related to thein the cyclesgraph represented by A. In ha (A), each Tr Ak equals thetotal weight of all length-k cycles in A – where the weightof a cycle is the product of its edge weights (Bapat, 2010).
The next theorem generalizes the result of Wei et al. (2020)
to show that most ha can be used to characterize acyclicity.

These three criteria are all important for maximizing Sα,β (θ)
under the constraint h(Aθ ) = 0. D-stability and E-stability
ensure the constraint remains well-defined and with boundedvalues throughout the optimization procedure. The Vstability is related to the nature of constrained optimization.
We now show that PST constraints are unstable as d grows.
Theorem 2.4 (PST instability). For d ≥ 2, any PST constraint h is both E-unstable and V-unstable. More precisely,d- E-unstable. ∃A ∈ Rd×d
≥0 , h(sA) = Ωs→∞ (s )d
- V-unstable. ∃A ∈ Rd×d
≥0 , h(εA) = Oε→0+ (ε )Theorem 2.2 (PST constraint). For any sequence∗(ak )k∈N∗ ∈ RN
≥0 , if we have ak > 0 for all k ∈ J1, dK, then,for any matrix A ∈ Rd×d
≥0 , we have ha (A) = 0 ⇔ A is acyclic,ha (A) ≥ 0,
∇ha (A) = ha′ (A⊤ ) with a′k = (k + 1)ak+1 .

Also, any PST constraint for which fa has a finite radius ofconvergence is D-unstable (e.g., hlog , hinv ).

We say that ha is a PST constraint.

Theorem 2.4 is proved in Appendix C.2. It shows that the
instability of the PST constraints worsens exponentially in d.
Fig. 1 empirically corroborates the theorem with two typesof adjacency matrices encountered during DCD: a cycle
and some uniformly random noise. It shows that all PST
constraints escalate to infinity or vanish to zero as the scaleof noise ε changes (Fig. 1 left) or as the number of variables
d increases (Fig. 1 right), reflecting their E-instability andV-instability. In addition, the D-instability of hlog and hinv
appears even in small ϵ or d (vertical lines). We encounterall three instabilities during causal discovery experiments(Section 4.1), leading existing approaches to fail.

The proof is in Appendix C.1. In particular, sequencesof strictly positive ak satisfy the conditions for any d, so
several standard power series are PST constraints.
1
For example, the sequence aexp
= k!
recovers the penaltykhexp (A) = Tr(exp(A)) − d originally proposed in Zheng
et al. (2018).

Table 1 shows that other constraints introduced in the literature are also PST (Ng et al., 2020; Bello et al., 2022; Yu
et al., 2019; Zheng et al., 2018).
2.2. Limitations of PST constraints
In this section, we provide the criteria necessary for constraints to exhibit stable optimization behavior. We provethat PST constraints do not satisfy these criteria and show
empirically that optimization with these constraints can be
slow or fail. As a solution, we suggest an alternative, nonPST acyclicity constraint and demonstrate its stability.

2.3. The Spectral Acyclicity Constraint
To overcome the limits of the PST constraint family, we
propose to use another type of constraint, one based on the
spectrum of A, which was first used in Lee et al. (2019).
This constraint draws from a characterization of DAG ma2

Stable Differentiable Causal Discovery:

random weights in [0,ε]length-d/2 cycle of weights ε

To further understand the impact of the constraints’ stability
on optimization, we empirically study the optimization pathof the augmented Lagrangian and the penalty method witheach constraint in Appendix Fig. 5. The instabilities of PST
constraints effectively slow their convergence and require
increasing γ and µ to excessively large values. In contrast,
the optimization paths with hρ take the least number of
iterations to converge, especially with the penalty method.

:

We are ready to perform DCD with the stable hρ .

3. Stable Differentiable Causal Discovery

Figure 1: Constraint behaviors when evaluated on uniform
random matrices in [0, ϵ]d×d (dashed) or a cycle of length
d/2 with weights ε (solid). The y-axis is the constraint’svalue, and the x-axis is (left) the weights’ scale ε (right) the
number of variables d. Only the proposed hρ (orange) remains stable; others vanish to zero exponentially or escalate
to infinity (as soon as d > 10). The vertical dotted lines
indicate the constraint escaped its domain of definition. All
these scenarios were encountered during DCD experiments.

With the stable acyclicity constraint hρ in hand, we now
introduce Stable Differentiable Causal Discovery (SDCD).
SDCD efficiently learns causal graphs in two stages.
3.1. The SDCD method
To solve the DCD optimization problem (8) with the spectral
acyclicity constraint hρ , SDCD optimizes the following
objective with gradient-based optimization:

trices from graph theory - that A is acyclic if and only if allits eigenvalues are zero (Cvetković et al., 1980).

θ̂ = arg max Sα,β (θ) − γ · hρ (Aθ ),We write λ1 (A) to λd (A) ∈ C, the d eigenvalues of A,
sorted from smallest to highest complex magnitude The
spectral radiushρ (A) = |λd (A)|,(1)

where hρ is used as a penalty with coefficient γ. SDCD
proceeds in two stages.

(2)θStage 1: Edge Preselection. First, SDCD solves Eq. (2)without the constraint, by setting γ = 0.

is the largest eigenvalue magnitude of A.

θ̂1 = arg max Sα1 ,β1 (θ)

The next theorem shows that the spectral radius can be used
as an acyclicity constraint.
Theorem 2.5 (Cvetković et al. (1980); Lee et al. (2019)).
The spectral radius Eq. (1) is an acyclicity constraint.

(3)θ

∀j,Aθ,jj =0This stage amounts to solving simultaneously d independent
prediction problems of each variable given the others (theconstraint Aθ,jj = 0 prevents self-loops). The goal is to
identify nonpredictive edges and remove them in stage 2,
akin to feature selection.

hρ (A) = 0 ⇔ A is a DAG.
We refer to it as the spectral acyclicity constraint. It isdifferentiable almost everywhere, with gradient

SDCD selects the removed edges as R̂1 = {(j, l) ∈ J1, dK2 |Aθ̂1 ,jl < τ1 } where τ1 is a threshold.

⊤∇hρ (A) = vd u⊤
d /vd ud ,Stage 2: Differentiable Causal Discovery. Next, SDCD
re-solves Eq. (2), this time with the constraint and withmasking the removed edges from stage 1,

where ud , vd are respectively the right and left eigenvectors
associated with λd (A) (Magnus, 1985).
Theorem 2.5 is proved in Appendix C.3. It implies that hρ
is D-stable. Next, we prove hρ is E-and-V-stable.
Theorem 2.6. hρ is stable.

θ̂2 =arg max

Sα2 ,β2 (θ) − γ2 hρ (Aθ ).

(4)

θ∀(j,l)∈R̂1 ,Aθ,jl =0

The term γ2 is initialized at 0 and is increased by a constant,γδ , after each epoch. Like other DCD methods, SDCD
forms the final graph ĜSDCD by selecting the edges in Aθ̂2
with weight above a threshold τ2 . The details of the algorithm can be found in Appendix E.2.

We refer to Appendix C.4 for the proof.
Remark 2.7. As a corollary of Theorem 2.6, hρ is not another PST constraint (since it is stable).
We complete Fig. 1 with the empirical behavior of hρ . As
theoretically expected, hρ retains non-extreme values and issuitable for constraint-based optimization.

Compared to other methods, SDCD innovates in two ways:(1) by using the constraint hρ with the penalty method and3

Stable Differentiable Causal Discovery# of edges
SDCD

DCDI
NOTEARSNOBEARS
DCDFG

DAGMA
SCORE

Number of true edges = d x s

s=2

s=4

= Edge Density = 2s / (d-1)

SHD

s=6

d=20

2k

200

1k

100

10

20

30
d

40

50

0

d=30

0

sortnregress
GIES

50

200 300 400 500dd=40

Figure 2: SHD on observational data across simulations
with increasing numbers of variables d. SDCD achieves thebest SHDs. It is the only method scaling above 200 variables
with nontrivial SHD. Missing data points imply the method
failed to run. Error bars indicate std on 30 random datasetsfor d ≤ 50 and five for d > 50 (175 total). Lower is better.

Fraction of variables intervened on

Figure 3: SHD across simulations with an increasing
proportion of variables intervened on, varying the total
number of variables d (columns) and average edges pervariable s (rows). SDCD is the only method to consistently
improve with interventional data and has the best SHDs forsparse graphs (edge density δ ≤ 45%). Each boxplot over
5 random datasets (45 datasets total).

(2) by using a two-stage optimization that preselects edges
in stage 1 and optimize the DCD objective only on those instage 2. Without explicit masking, stage 2 would be similar
to the barrier or penalty method (Ng et al., 2020; Bello et al.,2022), with stage 1 only providing a warm start.
3.2. Theoretical guarantees

for d > 50 NaNs appear right at initialization when hexp
overflows due to E-instability. DAGMA fails to converge
under 6 hours for as few as 30 variables, which we find iscaused by the learned adjacency matrix escaping the domainof definition of hlog due to D-instability.

In Appendix D, we show that theoretically, stage 1 does
not remove true causal parents, and so, stage 2 returns an
optimal graph. In the next section, we study SDCD’s empirical performance to examine the impact of finite data,nonconvex optimization, and relaxations.

To show that SDCD’s performance is robust to a comprehensive set of scenarios, we provide additional metrics forthese experiments in Appendix F.5 and Figs. 7 and 10.

4. Empirical studies
We compare SDCD to state-of-the-art baselines on multiple
datasets. We find that SDCD achieves significantly betterscores in both observational and interventional settings, particularly excelling at recovering sparser graphs. SDCD is
the only method to scale to thousands of variables without
sacrificing accuracy. We provide details on the evaluation
setup and additional experiments in Appendix F.

The runtimes associated with Fig. 2 are presented in Appendix Fig. 6. SDCD-GPU runs under 15 minutes for all
values of d in Fig. 2 experiments (e.g., d = 500). In Appendix Fig. 8, we further demonstrate that SDCD can scale
up to 4,000 variables under 2h45.
4.2. Interventional Data Experiments

4.1. Observational Data Experiments

Next, we compare SDCD, DCDI, DCD-FG, and GIES overdatasets with an increasing proportion of intervened variables. We characterize the edge density of a graph, δ, as theratio of true edges to the maximum number of edges possible in a DAG. We show the results for SDCD and DCDI
in Figure 3 and all methods in the Appendix (DCD-FG
and GIES performed consistently worse). As expected, themethods generally improve with more interventional data,
although SDCD is the only method to do so consistently.
We find that SDCD performs the best in most scenarios,
particularly on sparser graphs.

We evaluate all eight methods on a wide range of numberof variables d, with a fixed average number of edges pervariable s = 4, and repeat the experiments over 30 random
datasets. Fig. 2 reports the results and detailed tables areprovided in Appendix G with additional baselines.
SDCD outperforms the other methods in accuracy at every
scale and speed. DCDI is competitive on small d but crashes
for d > 40 – as discussed in 2.2, for d = 50, NaNs appear
during training when hexp underflows due to V-instability;4

Stable Differentiable Causal Discovery

References

Lambert, S. A., Jolma, A., Campitelli, L. F., Das, P. K., Yin,
Y., Albu, M., Chen, X., Taipale, J., Hughes, T. R., andWeirauch, M. T. The human transcription factors. Cell,
172(4):650–665, 2018.

Bapat, R. B. Graphs and matrices, volume 27. Springer,
2010.
Bello, K., Aragam, B., and Ravikumar, P. DAGMA:
Learning DAGs via M-matrices and a Log-Determinant
Acyclicity Characterization. In Neural Information Processing Systems, 2022.

Lee, H.-C., Danieletto, M., Miotto, R., Cherng, S. T., andDudley, J. T. Scaling structural learning with NO-BEARS
to infer causal transcriptome networks. In Pacific Symposium on Biocomputing, 2019.

Brouillard, P., Lachapelle, S., Lacoste, A., Lacoste-Julien,S., and Drouin, A. Differentiable causal discovery frominterventional data. In Neural Information Processing
Systems, 2020.

Lippe, P., Cohen, T., and Gavves, E. Efficient neural causal
discovery without acyclicity constraints. arXiv preprint
arXiv:2107.10483, 2021.
Loh, P.-L. and Bühlmann, P. High-dimensional learning of
linear causal networks via inverse covariance estimation.
The Journal of Machine Learning Research, 15(1):3065–
3105, 2014.

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

Lopez, R., Hütter, J.-C., Pritchard, J., and Regev, A. Largescale differentiable causal discovery of factor graphs. In
Neural Information Processing Systems, 2022.

Chickering, D. M. Learning Bayesian networks is NPcomplete. Learning from Data: Artificial Intelligenceand Statistics V, pp. 121–130, 1996.
Chickering, D. M. Optimal structure identification withgreedy search. Journal of Machine Learning Research, 3(Nov):507–554, 2002.

Lorch, L., Sussex, S., Rothfuss, J., Krause, A., andSchölkopf, B. Amortized inference for causal structurelearning. Advances in Neural Information Processing
Systems, 35:13104–13118, 2022.

Cvetković, D. M., Doob, M., and Sachs, H. Spectra ofGraphs: Theory and Application. Academic Press, 1980.

Magnus, J. R. On differentiating eigenvalues and eigenvectors. Econometric theory, 1(2):179–191, 1985.

Deng, C., Bello, K., Aragam, B., and Ravikumar, P. K.
Optimizing notears objectives via topological swaps. In
International Conference on Machine Learning, pp. 7563–
7595. PMLR, 2023.

Neapolitan, R. E. et al. Learning Bayesian Networks. Prentice Hall, 2004.
Ng, I., Ghassami, A., and Zhang, K. On the role of sparsity
and DAG constraints for learning linear DAGs. In Neural
Information Processing Systems, 2020.

Dixit, A., Parnas, O., Li, B., Chen, J., Fulco, C. P., JerbyArnon, L., Marjanovic, N. D., Dionne, D., Burks, T., Raychowdhury, R., et al. Perturb-seq: dissecting molecular
circuits with scalable single-cell rna profiling of pooled
genetic screens. Cell, 167(7):1853–1866, 2016.

Ng, I., Lachapelle, S., Ke, N. R., Lacoste-Julien, S., andZhang, K. On the convergence of continuous constrained
optimization for structure learning. In International Conference on Artificial Intelligence and Statistics, 2022.
Ng, I., Huang, B., and Zhang, K. Structure learning with
continuous optimization: A sober look and beyond. In
Causal Learning and Reasoning, pp. 71–105. PMLR,
2024.

Eberhardt, F. and Scheines, R. Interventions and causal
inference. Philosophy of Science, 74(5):981–995, 2007.
Glymour, C., Zhang, K., and Spirtes, P. Review of causal
discovery methods based on graphical models. Frontiersin genetics, 10:524, 2019.

Peters, J. and Bühlmann, P. Identifiability of gaussian
structural equation models with equal error variances.
Biometrika, 101(1):219–228, 2014.

Hauser, A. and Bühlmann, P. Characterization and greedy
learning of interventional markov equivalence classes of
directed acyclic graphs. Journal of Machine Learning
Research, 13(1):2409–2464, 2012.

Peters, J., Mooij, J. M., Janzing, D., and Schölkopf, B.
Causal discovery with continuous additive noise models.
Journal of Machine Learning Research, 15:2009–2053,
2014.

Hoover, K. D. Causality in Economics and Econometrics,
pp. 1–13. Palgrave Macmillan UK, 2006.

Reisach, A. G., Seiler, C., and Weichwald, S. Beware of thesimulated DAG! causal discovery benchmarks may beeasy to game. In Neural Information Processing Systems,
2021.

Horn, R. A. and Johnson, C. R. Matrix analysis. Cambridge
university press, 2012.
Lachapelle, S., Brouillard, P., Deleu, T., and Lacoste-Julien,
S. Gradient-based neural DAG learning. In International
Conference on Learning Representations, 2019.

Replogle, J. M., Saunders, R. A., Pogson, A. N., Hussmann,
J. A., Lenail, A., Guna, A., Mascibroda, L., Wagner,
5

Stable Differentiable Causal Discovery

E. J., Adelman, K., Lithwick-Yanai, G., et al. Mapping information-rich genotype-phenotype landscapeswith genome-scale perturb-seq. Cell, 185(14):2559–2575,
2022.

large-scale human crisis. In Proceedings of the National
Academy of Sciences, 2011.
Zheng, X., Aragam, B., Ravikumar, P. K., and Xing, E. P.
DAGs with NO TEARS: Continuous optimization forstructure learning. In Neural Information Processing
Systems, 2018.

Rolland, P., Cevher, V., Kleindessner, M., Russell, C., Janzing, D., Schölkopf, B., and Locatello, F. Score matchingenables causal discovery of nonlinear additive noise models. In Chaudhuri, K., Jegelka, S., Song, L., Szepesvari,C., Niu, G., and Sabato, S. (eds.), Proceedings of the 39th
International Conference on Machine Learning, volume
162 of Proceedings of Machine Learning Research, pp.
18741–18753. PMLR, 2022.

Zheng, X., Dan, C., Aragam, B., Ravikumar, P., and Xing,
E. Learning sparse nonparametric dags. In International
Conference on Artificial Intelligence and Statistics, 2020.

Sachs, K., Perez, O., Pe’er, D., Lauffenburger, D. A., andNolan, G. P. Causal protein-signaling networks derived
from multiparameter single-cell data. Science, pp. 523–
529, 2005.
Schölkopf, B., Locatello, F., Bauer, S., Ke, N. R., Kalchbrenner, N., Goyal, A., and Bengio, Y. Toward causal
representation learning. Proceedings of the IEEE, pp.
612–634, 2021.
Spirtes, P., Glymour, C. N., and Scheines, R. Causation,
prediction, and search. MIT press, 2000.
Triantafillou, S. and Tsamardinos, I. Constraint-based
causal discovery from multiple interventions over overlapping variable sets. Journal of Machine Learning Research,
16(1):2147–2205, 2015.
Tsamardinos, I., Brown, L. E., and Aliferis, C. F. The maxmin hill-climbing bayesian network structure learningalgorithm. Machine learning, 65:31–78, 2006.
Wang, Z. and Chan, L. Learning bayesian networks frommarkov random fields: an efficient algorithm for linear
models. ACM Transactions on Knowledge Discovery
from Data (TKDD), 6(3):1–31, 2012.
Wei, D., Gao, T., and Yu, Y. Dags with no fears: A closer
look at continuous optimization for learning bayesian
networks. Advances in Neural Information Processing
Systems, 33:3895–3906, 2020.
Yang, K., Katcoff, A., and Uhler, C. Characterizing andlearning equivalence classes of causal dags under interventions. In International Conference on Machine Learning, 2018.
Yu, Y., Chen, J., Gao, T., and Yu, M. DAG-GNN: DAG
structure learning with graph neural networks. In International Conference on Machine Learning, 2019.
Yu, Y., Gao, T., Yin, N., and Ji, Q. Dags with no curl:
An efficient dag structure learning approach. In International Conference on Machine Learning, pp. 12156–
12166. Pmlr, 2021.
Zhang, D. D., Lee, H. F., Wang, C., Li, B., Pei, Q., Zhang, J.,
and An, Y. The causality analysis of climate change and6Supplementary Materials: Stable Differentiable Causal Discovery

A. Related Work
Causal discovery methods mainly fall into two categories: constraint-based methods and score-based methods (Glymour
et al., 2019).
Constraint-based methods identify causal relationships by testing for conditional independence among variables in the data.
For example, the PC algorithm (Spirtes et al., 2000) finds the graphs which conform to all the independencies present in thedata. COmbINE is an extension to support interventional data (Triantafillou & Tsamardinos, 2015).
On the other hand, score-based methods design a score S(G) that is maximized by the true graph G∗ , and they aim to findits maximizer Ĝ = arg maxG∈DAG S(G). Existing score-based methods differ in their choice of S and of the optimization
method to maximize it. GES (Chickering, 2002) and GDS (Peters & Bühlmann, 2014) optimize the BIC score of a Gaussian
linear model by greedily adding or removing edges. GIES modifies GES to support interventional data (Hauser & Bühlmann,2012), and CAM supports non-linear additive models (Bühlmann et al., 2014).
Differentiable Causal Discovery (DCD), which our work extends, is a type of score-based approach that reformulates thesearch of the score maximizer into a continuous optimization problem. It uses a numerical criterion to distinguish acyclic
graphs from cyclic ones (called the acyclicity constraint). It was initially introduced in (Zheng et al., 2018) as NO-TEARS,
which uses linear models, augmented Lagrangian optimization, and a constraint based on the adjacency matrix exponential.
Other works extend the methodology to incorporate polynomial regression (Lee et al., 2019), neural networks (Lachapelle
et al., 2019; Zheng et al., 2020), and support for interventional data (Brouillard et al., 2020).
Alternative acyclicity constraints (Lee et al., 2019; Ng et al., 2020; Bello et al., 2022) have been proposed, as well as
optimization schemes different than augmented Lagrangian (Ng et al., 2020, 2022; Bello et al., 2022). Deng et al. (2023)
introduced a hybrid approach combining gradient optimization with combinatorial optimization, while Lippe et al. (2021)
explored removing the acyclicity constraint entirely.
Despite a rich literature, DCD has difficulty scaling to a large number of variables, exhibiting long training times andnumerical instability. Wei et al. (2020); Ng et al. (2024) identified those limitations, and Lee et al. (2019); Lopez et al.
(2022) addressed these problems with elegant approximations, but they resulted in poor accuracy.
Here, we provide a theoretical understanding of some algorithmic issues with DCD and then use that understanding todevelop a better DCD method. Our proposed method builds on ideas that have been partially studied in previous work,
including the acyclicity constraint of Lee et al. (2019) and the penalty method of Ng et al. (2020). Lee et al. (2019) uses the
spectral acyclicity constraint for computational efficiency but otherwise does not expand on its advantages. In this work,
we provide a novel analysis of the constraint and incorporate it into a new strategy that is more accurate and scalable thanexisting DCD algorithms. Table 8 recapitulates research in DCD and how this work fits in.

7

Stable Differentiable Causal Discovery

B. Background and Notations
We review the Differentiable Causal Discovery (DCD) approach and define the notations used in the paper.
B.1. Background on Causal Discovery
Causal discovery is the task of learning cause-and-effect relationships among a set of variables. In this work, we consider
variables that can be intervened on such that they no longer are affected by their causal parents. These interventions arecalled structural or perfect interventions (Eberhardt & Scheines, 2007).
Causal Graphical Models. Causal graphical models (CGMs) provide a mathematical framework for reasoning aboutcausal relationships between variables. Consider a CGM over d variables and K possible interventions on it.

There are three components:
1. A directed acyclic graph (DAG), G∗ = (V, E), where each node, j ∈ V , represents a variable xj , and each edge,(j, k) ∈ E, indicates a direct causal relationship from xj to xk .
2. A list of conditional distributions, p∗j (xj | xpaG∗ ; 0), which specify the distribution of each xj given its causal parentsj
xpaG∗ without intervention (the 0 indicates no intervention).
j

3. A list of interventions, I = {I0 , I1 , ..., IK }, where I0 = ∅ (no intervention) and the others Ik ⊂ V define the target
variables of intervention k. Alongside is a list of interventional distributions, p∗m (xm ; k), for each k > 0 and m ∈ Ik ,which define the distributions over xm after intervention k.
The joint distribution under intervention k writes:p∗ (x; k) =Yp∗j (xj | xpaG∗ ; 0)jYp∗j (xj ; k).

(5)

j∈Ik

j∈V \Ik

Note p∗(x; 0) is the joint on observational data.
Data. We observe n data points of the d variables X = {(xi1 , ..., xid )}ni=1 with labels T = {ti }ni=1 where ti ∈ {0, ..., K}
indicates which intervention was applied to xi (0 indicating no intervention). For example, in genomics, Perturb-seq screens(Replogle et al., 2022) measure the expression of d genes across n cells, where each cell can be edited once to change theexpression of one of its genes.
Causal discovery with score-based methods. The goal of causal discovery is to infer the graph G∗ from the data (X, T ).
In particular, score-based methods assign a score S(G) to every possible graph G, where the score function is designed sothat it is maximized on the true graph G∗ . Score-based methods aim to find the maximizerĜ = arg max S(G).

(6)

G∈DAG

Fix a model class that defines the conditional (and interventional) distributions for each possible DAG G as {p(· | G; θ, k)}θ .
We can define the score S(G) to be the maximum log-likelihood that can be achieved under graph G with some regularization
over the number of edges |G| (Chickering, 2002):
#
n
1X
i
i
log p(x | G; θ, t ) − λ|G|.
Smle (G) = sup
n i=1
θ
"

(7)

In the limit of infinite samples (n → ∞), and under a few assumptions, any maximizer Ĝ of Eq. (7) is close to the true G∗(Brouillard et al., 2020). More precisely, Ĝ and G∗ are I-Markov-equivalent: they share the same skeleton, v-structures,and other restrictive properties at the intervened variables in I (Yang et al., 2018).
8

Stable Differentiable Causal Discovery

B.2. Differentiable Causal Discovery
The main challenge to a score-based method for causal discovery is how to search over the large space of DAGs. Differentiable Causal Discovery (DCD) reformulates this combinatorial search into a continuous optimization problem over the
space of all graphs, including cyclic ones (Zheng et al., 2018). It introduces three key components.
Model Class with Implicit Graph. First, DCD definesa model class with no apparent underlying graph, where eachQvariable is conditioned on all others as p(· | θ, k) ∝ pj (xj |x−j ; θ, k). Instead, θ defines the graph G implicitly, such thatif θ renders x−j 7→ pj (xj |x−j ; θ, 0) invariant to some xℓ ⊂ x−j , then there is no edge from l to j. The induced adjacency
matrix is denoted Aθ . When θ induces an acyclic Aθ , then p(· | θ, k) defines a valid CGM.
Acyclicity Function. Second, DCD introduces a differentiable function h(Aθ ) that quantifies how “cyclic” Aθ is. h(Aθ )is high when Aθ contains cycles with large edge weights, it is low when Aθ contains cycles with small weights, andh(Aθ ) = 0 when Aθ contains no cycles.

Optimization. Finally, DCD reformulates Eq. (7) into a constrained optimization problem only over θ.
θ̂ =arg maxSα,β (θ).

(8)θ

s.t. h(Aθ ) = 0.
It uses Aθ in place of G, uses the constraint h(Aθ ) = 0 in place of G ∈ DAG, and uses a new objective Sα,β :n1Xlog p(xi ; θ, ti ) − α∥Aθ ∥1 − β∥θ∥22 ,Sα,β (θ) =
n i=1

(9)

Sα,β is a relaxed version of Smle (Eq. (7)) where the discrete |G| is changed into an L1 regularization of Aθ (for α ≥ 0)
and an L2 regularization of θ is included (for β ≥ 0). The supθ in Smle (Eq. (7)) is now removed, as it merges with thearg maxθ of Eq. (8).
With Eq. (8) in hand, different methods for DCD solve the constrained optimization in different ways. Some approaches use
the augmented Lagrangian method (Zheng et al., 2018; Lachapelle et al., 2019; Brouillard et al., 2020; Lopez et al., 2022),
some use the barrier method (Bello et al., 2022), and others use h as a regularizing penalty (Ng et al., 2020).
In all these approaches, the choices of h and the optimization method dictate the optimization behavior and the ultimate
quality of the inferred graph. In the next section, we highlight the importance of h.

9Stable Differentiable Causal Discovery

C. Theoretical results
C.1. Proof of Theorem 2.2
Before proving Theorem 2.2, we precisely define an acyclic matrix and prove a few lemmas.
Definition C.1 (Cyclic and acyclic matrices). Take a matrix A ∈ Rd×d
≥0 .
We say that A has a cycle of length k if and only if:∃(i0 , ..., ik ) ∈ J1, kKk+1 , such that,i0 = ik∀ℓ ∈ J1, kK, Aiℓ−1 ,iℓ > 0(10)

We say that A is cyclic if it contains at least one cycle. We note that if A contains a cycle of length k for k ∈ N∗ (the setof strictly positive integers), then A also contains a cycle of length k ′ for k ′ ∈ J1, dK (this follows from the pigeon holeprinciple).
We say that A is acyclic if it does not contain any cycle (or equivalently if it does not contain any cycle of length k ≤ d).
Lemma C.2. For any matrix A ∈ Rd×d
≥0 , • Tr Ak ≥ 0 for any k • A has a cycle of length k if and only if Tr Ak > 0• Ad = 0 if and only if A is acyclic.

Proof. Fix a matrix A ∈ Rd×d
≥0 . We have, 
Tr Ak =kY

X

Aℓ−1,ℓ .

(11)(i0 ,...,ik )∈J1,dKk+1 ℓ=1i0 =ik =i 
Each addend is non-negative so Tr Ak ≥ 0. Furthermore, the total sum is strictly positive if and only if at least one addend
is strictly positive. This happens if and only if A has a cycle of length k by definition.
Similarly, we have
(Ad )i,j =XkY

Aℓ−1,ℓ .

(12)

d+1 ℓ=1(i0 ,...,id )∈J1,dK
i0 =id =j

If Adi,j > 0, then one addend is strictly positive and so there exists (i0 , ..., id ) ∈ J1, dKd+1 such that i0 = id = j andQkℓ=1 Aℓ−1,ℓ > 0. By the pigeon-hole principle, two iℓ are identical, which provides a cycle. Reciprocally, if (i0 , ..., ik )is a cycle of length k, then by repeating (i0 , ..., ik , i1 , i2 ...id mod k ) until having a path of length d + 1, we have that
(Ad )i0 ,id mod k > 0. Hence, Ad = 0 if and only if A is acyclic.
We recall Theorem 2.2.
∗

Theorem 2.2 (PST constraint). For any sequence (ak )k∈N∗ ∈ RN
≥0 , if we have ak > 0 for all k ∈ J1, dK, then, for any
d×d
matrix A ∈ R≥0 , we have ha (A) = 0 ⇔ A is acyclic,ha (A) ≥ 0,
∇ha (A) = ha′ (A⊤ ) with a′k = (k + 1)ak+1 .

We say that ha is a PST constraint.

10

Stable Differentiable Causal Discovery∗

N
Proof. Fix a matrix A ∈ Rd×d
≥0 and a sequence (ak )k∈N∗ ∈ R≥0 such that ak > 0 for any k ∈ J1, dK.

By definition, we have,ha (A) = Tr

"+∞
X

k=1

=

+∞
X

ak A

k

#

 
ak Tr Ak

k=1

(13)

(14)
(15)

 
1. By Lemma C.2, Tr Ak ≥ 0 and so ha (A) ≥ 0. This proves the second property.
 
2. Then, ha (A) = 0 if and only if Tr Ak = 0 for all k for which ak > 0. Since ak > 0 for any k ∈ J1, dK we conclude
that if ha (A) = 0, then A does not contain cycles of length k ≤ d, so A is acyclic by Definition C.1. Reciprocally, if A
is acyclic, it does not contain cycles of any length, so ha (A) = 0.
3. Finally, if we write ra the radius of convergence of fa , then ha converge
P+∞ absolutely over the set of matrices withhρ (A) < ra so it is differentiable with gradient given by: ∇ha (A) = k=1 ak k(A⊤ )k−1 .

This concludes the proof.

C.2. Proof of Theorem 2
We recall Theorem 2.4.
Theorem 2.4 (PST instability). For d ≥ 2, any PST constraint h is both E-unstable and V-unstable. More precisely, dE-unstable. ∃A ∈ Rd×d
≥0 , h(sA) = Ωs→∞ (s )d×d
- V-unstable. ∃A ∈ R≥0 , h(εA) = Oε→0+ (εd )Also, any PST constraint for which fa has a finite radius of convergence is D-unstable (e.g., hlog , hinv ).
Proof. Take a PST constraint ha for some (ak )k ∈ Rd×d
≥0 with ak > 0 for k ∈ J1, dK.
We will show the E-unstable and V-unstable results using a particular adjacency matrix C.
Define C as the adjacency matrix of the cycle 1 → 2 → ... → d → 1 with edges weights of 1. That is:
0 1 0 ... ... 0
0 1
0 . . . 0


 ..
.. 
..
..
..
.

.
.
.
.


C=

..
.
.
.
.

. 0
.
.




..
0
. 1
1 0
...
0
d



We have C = Id and Tr Ck

We obtain for any w ∈ R≥0 ,



=



d
0

(16)

if k = 0 mod d.
if k ̸= 0 mod dha (wC) = d+∞Xaℓd wℓd .

(17)

ℓ=1

• In particular, we have for any s ≥ 0, ha (wC) = dad sd = Ωs→+∞ sd (since ad > 0). This proves the E-instability.
11

Stable Differentiable Causal Discovery• Define u = min(1, ra /2) where ra is the radius of convergence of fa .
Then, for any ε ∈ [0, u2 ],ha (εC) = d+∞Xaℓd εℓd

(18)

ℓ=1

d

=ε d
d

≤ε d
≤ εd d

+∞
X

ℓ=1
+∞
X

ℓ=1
+∞
X

aℓd ε

(ℓ−1)d

aℓd u

!

(19)

!

(20)

!

(21)

2(ℓ−1)daℓd uℓd + adℓ=1

≤ εd d (fa (uC) + ad ) .

(22)

d

(23)

= Oε→0+ (ε )Where we obtain Eq. (21) by noting that 2(ℓ − 1) ≥ ℓ and u ≤ 1. Finally, since u < ra , fa (uC) is finite. Hence theresult.
The D-instability result follows from the definition of the radius of convergence.
C.3. Proof of Theorem 3
We recall Theorem 2.5.
Theorem 2.5 (Cvetković et al. (1980); Lee et al. (2019)). The spectral radius Eq. (1) is an acyclicity constraint.
hρ (A) = 0 ⇔ A is a DAG.
We refer to it as the spectral acyclicity constraint. It is differentiable almost everywhere, with gradient⊤∇hρ (A) = vd u⊤
d /vd ud ,where ud , vd are respectively the right and left eigenvectors associated with λd (A) (Magnus, 1985).
The two properties stated in Theorem 2.5 are standard results.
Proof.

• We provide proof for the statement hρ (A) = 0 ⇔ A is a DAG for the sake of completeness.
  Pd ⇒ If hρ (A) = 0 then all eigenvalues λj (A) are zeros. But since Tr Ak = j=1 λj (A)k , we have Tr Ak = 0 for
any k ≥ 1 and by Lemma C.2, A is acyclic.
⇐ Assume A is acyclic, then Ad = 0 by Lemma C.2. But then all eigenvalues are 0 (as for eigenvalue λj (A) andassociated eigenvector vj (A), we have Ad vj (A) = λj (A)d vj (A) = 0.
Hence, hρ is a valid acyclicity constraint.
• Magnus (1985) shows that hρ is differentiable at every A that has mutually distinct eigenvalues, with the formulaprovided in Theorem 2.6. The set of matrices with all distinct eigenvalues is dense in the set of matrices (Horn &Johnson, 2012)[Theorem 2.4.7.1], which proves the result

12Stable Differentiable Causal Discovery

C.4. Proof of Theorem 4
We recall Theorem 2.6.
Theorem 2.6. hρ is stable.
Proof. We prove each stability criterion.
• E-stable: For any s > 0 and matrix A, hρ (sA) = |s|hρ (A) = Os→+∞ (s).
• V-stable:For any ε > 0 and matrix A such that hρ (A) > 0, hρ (εA) = |ε|hρ (A) = Ωε→0+ (ε).
• D-stable: Every matrix has eigenvalues (C is algebraically closed), so hρ is well defined everywhere. In addition,
Theorem 2.5 proved that hρ was differentiable almost everywhere.
Hence, hρ is a stable constraint.

13

Stable Differentiable Causal Discovery

D. Theoretical guarantees of SDCD
We provide correctness guarantees for the two stages of SDCD. We show that, theoretically, stage 1 does not remove truecausal parents, and so stage 2 returns an optimal graph.
With infinite data, Eq. (3)’s unrelaxed version writes,θ̃1 = arg max
θKX∀j,Aθ,jj =0 k=0
Ik ̸∋jπk∗ E [log pj (xj |x−j ; θ, 0)]−λ|Aθ |,(24)

p (x;k)where πk is the proportion of data coming from intervention k. The next theorem characterizes the graph G̃1 = Aθ̃1 interms of Markov boundaries in the true graph G∗ . A Markov boundary for j is a minimal set of variables that render j
independent of all the others. In a causal graph, each j has a unique Markov boundary, consisting of j’s parents, j’s children,and j’s children’s parents (Neapolitan et al., 2004).
1
Theorem D.1. Under regularity assumptions detailed below, the candidate parents paG̃of j selected by stage 1 are
j
∗

G
G
G
1
precisely the Markov boundary of j in the true graph G∗ , That is, paG̃j = paj ∪ chj ∪ pachG∗ .
∗∗j

The assumptions are similar to the ones detailed in Brouillard et al. (2020) to guarantee that differentiable causal discovery
can identify causal graphs. The assumptions are:• π0 > 0 – we observe some observational data,• ∃θ, s.t. ∀k, p∗ (· ; k) = p(· ; θ, k) – the model class can express the true model p∗ ,• The observational distribution p∗ (x; 0) is faithful to the graph G∗ (that is any edge in G∗ indeed result in a nonzero
cause-and-effect relation in the distribution p∗ (x; 0). See Neapolitan et al. (2004) for more details.
• The true distributions p∗ (x; k) and any distribution of the model class p(x; θ, k) have strictly positive density p∗ (x; k) >0, p(x; θ, k) > 0. This avoids technical difficulty when forming conditional distributions (e.g., p∗ (xj |xT ; k)).
• The expectations Ep∗ (x;k) [log p∗ (x; k)] are well defined (they are finite). This enables us to consider the likelihood
expectations in the first place.
• The regularization strength λ is strictly positive and small enough (see the proof for how small).
The assumptions are reasonable: p∗ should be in the model class {pθ }, the expectations should be well defined, and“faithfulness” should hold (that is, G∗ doesn’t have superfluous edges).
Theorem D.1 makes two guarantees: (1) stage 1 does not remove causal parents and (2) stage 1 returns only a subset of theedges, not all of them. For instance, if G∗ is sparse such that each node has at most k parents, then only O(dk 2 ) edges are
returned, which is essentially linear in d for small k (Appendix D.2).
Theorem D.1 implies that Brouillard et al. (2020, Theorem 1) still applies, and we deduce that stage 2 remains optimal under
the stated assumptions (see Appendix D.1).
Proof. Fix j ∈ J1, jK.
For clarity of notations, we rewrite ?? as1paG̃
j = arg max supS⊂J1,dK\{j}

θ

K
X

k=0
j̸∈Ik

πk ∗ Ep (x;k)[log pj (xj |xS ; θ, k)] − λ|S|,(25)

θ
where the condition S = paA
j is fully captured by the notation pj (xj |xS ; θ, k).

Then, defineψ(T ) = supθXEp∗ (x;k) [πk log pj (xj |xS ; θ, k)] − λ|S|.

k
j̸∈Ik

14

(26)

Stable Differentiable Causal Discovery∗∗
Further, define B = boG
j to be the Markov boundary of node j in the true causal graph G .

We will show that ψ(B) > ψ(T ) for any other T ⊂ J1, dK\ {j}.
We compute,ψ(B) − ψ(T ) = supθXπk Ep∗ (x;k) [log pj (xj |xB ; θ, k)] − sup
θ

k
j̸∈Ik

X

πk Ep∗ (x;k) [log pj (xj ; xT ; θ, k)]− λ|B| + λ|T |X= − inf
πk Ep∗ (x−j ;k) DKL p∗j (xj |x−j ; k) ∥ pj (xj |xB ; θ, k)
θ

+ inf
θ

k
j̸∈Ik

X

k
j̸∈Ik

(28)



πk Ep∗ (x−j ;k) DKL p∗j (xj |x−j ; k) ∥ pj (xj ; xT ; θ, k)− λ(|B| − |T |)
X


= inf
πk Ep∗ (x−j ;k) DKL p∗j (xj |x−j ; k) ∥ pj (xj ; xT ; θ, k)
θ

(27)

k
j̸∈Ik

(29)

k
j̸∈Ik

+ λ(|T | − |B|).
The line 28 comes from Ep∗ (x;k) [log pj (xj |xB ; θ, k)] = −Ep∗ (x−j ;k) DKL p∗j (xj |x−j ; k) ∥ pj (xj |xB ; θ, k) +Ep∗ (x;k) log p∗j (xj |x−j ; k) where we added and substracted the log p(k) term (the Ep∗ (x;k) is decomposed intoEp∗ (x−j ;k) Ep∗ (xj ;k) , where the second expectation is in the KL divergence). We use the assumption of strictly positive
density here to define the conditional p∗j (xj |x−j ; k) without technical difficulties.
The line 29 comes from the assumption of sufficient model class capacity and the definition of the Markov boundary. Indeed,
we first have p∗j (xj |x−j ; k) = p∗j (xj |xB ; k) by definition of the Markov boundary B, and since the model class is expressive
enough, there exists θ such that DKL p∗j (xj |x−j ; k) ∥ pj (xj |xB ; θ, k) = 0.

We further have:ψ(B) − ψ(T ) ≥ π0 inf Ep∗ (x−j ;0) DKL p∗j (xj |xB ; 0) ∥ pj (xj |xT ; θ, 0) + λ(|T | − |B|)
θ




p∗j (xj |xT ; 0)∗∗
= π0 Ep∗ (x−j ;0) DKL pj (xj |xB ; 0) ∥ pj (xj |xT ; 0) + π0 inf Ep∗ (x;0) log
θpj (xj |xT ; θ, 0)+ λ(|T | − |B|)
≥ π0 Ep∗ (x−j ;0) DKL p∗j (xj |xB ; 0) ∥ p∗j (xj |xT ; 0) +λ(|T | − |B|).
|
{z
}

(30)
(31)

(32)

η(T )

i


p∗ (xj |xT ;0)where line 32 follows from Ep∗ (x;0) log pjj(xj |xT ;θ,0) = Ep∗ (xT ) DKL p∗j (xj |xT )∥pj (xj |xT ; θ, 0) ≥ 0.
o

n
η(T )
|
T
⊂
J1,
dK\{j}
and
η(T
)
>
0
and
|B|
>
|T
|
∪
{1}
and fix any λ ∈]0, u[.
Let’s finally define u = min|B|−|T |hLet’s assume now that ψ(T ) ≥ ψ(B) for some T ⊂ J1, dK\{j}, and show that we obtain contradictions.

First, we would have λ(|B| − |T |) ≥ η(T ). In particular we deduce that |B| ≥ |T | (since η(T ) ≥ 0).
Now, two possibilities:
1. If η(T ) > 0, then |B| > |T | and by definition of λ, λ > λ which is absurd.

2. If η(T ) = 0, then π0 Ep∗ (x−j ;0) DKL p∗j (xj |xB ; 0) ∥ p∗j (xj |xT ; 0)
= 0.
This implies that
DKL p∗j (xj |xB ; 0) ∥ p∗j (xj |xT ; 0) = 0 for all (x−j ); since p∗ (x−j ; 0) has positive density and π0 > 0. Hence, theconditional p∗j (xj |xB ; 0) and p∗j (xj |xT ; 0) are identical. Since B was the Markov boundary of xj , that makes T also a
15

Stable Differentiable Causal Discovery

Markov blanket of xj . But then, by minimality of the Markov boundary in a faithful graph, we have B ⊂ T . Rememberthat we had deduced |B| ≥ |T |. So B = T .
This ends the proof, where λ ∈]0, u[.
D.1. Proof for Stage 2
Since stage 1 does not remove any true causal parents, theorem 1 of Brouillard et al. (2020) remains valid.
D.2. Lemma: Asymptotic Bound on number of edges returned in Stage 1
∗

∗

∗

∗

∗

G
G
G
G
We denote the Markov boundary of j in G∗ by boG
j , and recall that boj = paj ∪ chj ∪ pachG∗ \{j}.
jThe following lemma upper-bounds the theoretical number of edges returned by stage 1 when each node has at most kparents.
Lemma D.2. Assume G∗ is sparse such that each node has at most k parents. Then, the total size of all the Markov
boundaries is upper-bounded by dk(k + 2) = O(dk 2 ).
Proof. First, note that if each node has at most k parents, then |E| ≤ dk. Finally,
X
X
∗
∗
|boG
|boG
j |
j |=

(33)

j∈V

j∈V

≤

X

∗

∗

∗

G
G
|paG
j | + |chj | + |pachG∗ |

(34)

j

j∈V

≤ |E| + |E| +

X X

∗

|paG
k |

(35)

j∈V k∈chG∗j≤ 2kd +

X X

∗

|paG
k |

(36)

k∈V j∈paG∗j≤ 2kd + dk 2

16

(37)

Stable Differentiable Causal Discovery

E. Methods
Stage 1
Edge Preselection

True edge

Stage 2Differentiable Causal Discovery

Preselected
edge

Inferrededge

Figure 4: Visual representation of the SDCD method.

E.1. Model Details
In SDCD, the conditional distributions, pj (xj |x−j ; θ, k), are modeled as Gaussian distributions where the mean and varianceare learned by a neural network that takes in all of the other x−j as input. The initial layer of the network applies d
independent linear transformations followed by a sigmoid nonlinearity to the input and outputs d hidden states of size 10.
Each of the d hidden states corresponds to the features then used to predict each variable. Each hidden state is fed into
two linear layers: one to predict the mean parameter of the conditional and one to predict the variance parameter of theconditional. For the variance, a softplus operation is applied to the output of the linear layer to constrain the variance to bestrictly positive.
E.2. Algorithm Details
E.2.1. S PECTRAL ACYCLICITY C ONSTRAINT E STIMATION⊤As described in Theorem 2.5, the gradient of the spectral acyclicity constraint can be computed as hρ (A) = vd u⊤
d /vd ud ,where ud , vd are the right and left eigenvectors of A respectively. Using the power iteration method, which involves a fixed
number of matrix-vector multiplications, ud , vd can be estimated in O(d2 ). Specifically, the updates are as follows:

(i+1)

ud
(1)

(i)

:=

(i)

A⊤ ud
Avd
(i+1)
, vd
:=
∥ud ∥2∥vd ∥2(1)

where ud , vd are initialized as ud , vd := [ √1d , . . . , √1d ] at the very first epoch of SDCD. In our implementation, we use
15 iterations to estimate the spectral acyclicity constraint value.
Importantly, we re-use the estimates of ud and vd from one epoch to another, as we don’t expect A (and its eigenvectors) tochange drastically.
Hence, at each epoch, we initialize ud , vd using their last epoch’s value and perform 15 power iterations.
E.2.2. SDCD A LGORITHM
The SDCD algorithm follows a two-stage procedure. In the first stage, the coefficient of the spectral acyclicity constraint, γ,
is fixed at zero. We use an Adam optimizer with a learning rate, η1 , specific to stage 1 to perform minibatch gradient-based
optimization. The coefficients corresponding to the L1 and L2 penalties, α1 and β1 , respectively, are fixed throughout
training. The stage 1 training loss is written as:
L1 (X, θ, α1 , β1 ) = Sα1 ,β1 (θ)
n1X
log p(xi ; θ, ti ) − α1 ∥Aθ ∥1 − β1 ∥θ∥22 .
=
n i=1
17

Stable Differentiable Causal Discovery

To prevent the model from learning implicit self-loops, the weights corresponding to the predicted variable are masked out
for every hidden state output by the initial neural network layer. Thus, the prediction of each variable is prevented frombeing a function of the same variable.
In interventional regimes, the log-likelihood terms corresponding to the prediction of intervened variables are zeroed out.
The intervened variables do not have to be modeled as we assume perfect interventions.
Stage 1 is run for a fixed number of epochs. By default, stage 1 also has an early stopping mechanism that uses thereconstruction loss of a held-out validation set of data (sampled uniformly at random from the training set) as the early
stopping metric. If the validation reconstruction loss does not achieve a new minimum after a given number of epochs, thestage 1 training loop is exited.
At the end of stage 1, the learned input layer weights are used to compute a set of removed edges, R̂, for stage 2. Let
W ∈ Rd×d×10 represent the input layer weights. Then, each value of the implicitly defined weighted adjacency matrix is
computed as the L2 vector norm for the corresponding set of weights (i.e., Aθ,i,j := ∥Wi,j,: ∥2 ). This weighted adjacency
matrix is discretized with a fixed threshold, τ1 , such that each edge, (i, j), is removed if it falls below the threshold (i.e.,Aθ,i,j < τ1 ).
In stage 2, the spectral acyclicity constraint is introduced. Like stage 1, we use an Adam optimizer with learning rate, η2 ,and perform minibatch gradient-based optimization. Once again, the L1 and L2 coefficients, α2 , β2 , are fixed throughout
training. Rather than a fixed γ, SDCD takes an increment value, γ + ∈ R+ , determining the rate at which γ increases every
epoch. The training loss for stage 2 is as follows:

L2 (X, θ, R̂, α2 , β2 , γ) = Sα2 ,β2 (θ) − γhρ (Aθ )n1X=
log p(xi ; θ, ti ) − α2 ∥Aθ ∥1 − β2 ∥θ∥22 − γhρ (Aθ ).
n i=1
The same masking strategy as in stage 1 is used to prevent self-loops in Aθ . However, the input layer weights corresponding
to edges (i, j) ∈ R̂ are also masked.
Like before, the reconstruction loss terms corresponding to intervened variables are removed from the loss.
To reduce the sensitivity of stage 2 to the choice of γ + and to prevent the acyclicity constraint term from dominating theloss, the linear increment schedule is frozen when Aθ achieves a DAG at the final threshold, τ2 . In practice, the DAG checkis performed every 20 epochs. If the adjacency matrix returns to being cyclic throughout training, the γ increment schedule
restarts to increase from where it left off.
The early stopping metric is computed similarly to stage 1, but in stage 2, the early stopping can only kick in when γ hasbeen frozen. If the γ schedule is resumed due to Aθ reintroducing a cycle, the early stopping is reset.
Lastly, once stage 2 is complete, Aθ is computed and thresholded according to a fixed threshold, τ2 . All values exceeding
the threshold (i.e., Aθ , i, j ≥ τ2 ) are considered edges in the final graph prediction.
The thresholded adjacency matrix may contain cycles if stage 2 runs to completion without hitting early stopping. To ensure
a DAG, we follow a greedy edge selection procedure detailed in Algorithm 2.
Pseudocode for a simplified SDCD algorithm (excludes γ freezing and early stopping) is provided in Algorithm 1.
E.2.3. M OTIVATION FOR S TAGE 1
Dedicating stage 1 to removing unlikely edges is motivated by the hypothesis that real-life causal graphs are sparse. For
example, individual genes in biological systems are typically regulated by a few other genes rather than all other genes(Lambert et al., 2018). A similar hypothesis underlies work in sparse mechanism shift (Schölkopf et al., 2021). Hence,stage 1 will likely remove many false edges and facilitate stage 2. Alternative approaches for variable selection (e.g.,
markov boundaries (Loh & Bühlmann, 2014; Wang & Chan, 2012), skeletons (Tsamardinos et al., 2006), preliminaryneighborhood selection (Bühlmann et al., 2014; Lachapelle et al., 2019)) can be motivated for the same reasons. Here, we
found a simple modification to the objective function can effectively serve this purpose. In practice, we find that stage 1
improves convergence speed and accuracy (Section 4, Table 6). Notably, we find that stage 1 improves the stability of the18

Stable Differentiable Causal Discovery

training in stage 2, even when PST constraints are used in place of the spectral one (Table 7). For this reason, stage 1 may
also serve as a beneficial preprocessing step for other causal discovery methods. In Theorem D.1, we prove that stage 1 does
not remove true causal parents.
E.2.4. T IME AND S PACE C OMPLEXITY
The time complexity of each iteration of SDCD is O(d2 ). The forward pass in stage 1 can be computed in O(d2 ) time.
On the other hand, each of the d prediction problems can be computed independently. This allows for parallelizing the dproblems, each taking O(d) time. Stage 2 also takes O(d2 ) time as the spectral acyclicity constraint and the forward pass
both take O(d2 ) time to compute. Thus, the time complexity of each iteration in both stages is O(d2 ).
If the sparsity pattern of the underlying causal graph is known beforehand such that each variable has at most k parents, we
can further tighten the time complexity of SDCD. By Appendix D.2, we know the size of the set of remaining edges afterstage 1 is O(dk 2 ). Using sparse matrix multiplication, the spectral acyclicity constraint can be done in O(dk 2 ), which iseffectively linear in d if k ≪ d. However, this improvement only becomes significant when d > 10, 000 (from experiments
not reported in this paper).
The space complexity of the algorithm is O(d2 ), as the number of parameters in the input layer scale quadratically in thenumber of features.
Algorithm 1 SDCD
Require: α1 , α2 ∈ R+ , β1 , β2 ∈ R+ , γ + ∈ R+ ,
τ1 , τ2 ∈ R+ , η1 , η2 ∈ R+ , E1 , E2 ∈ Z+
(0)Aθ ← ⃗0G×G(0)θ−Aθ ← RandomGaussianInit()e←0while e < E1 doθ(e+1) := AdamUpdate(θ(e) , ∇L1 (X, θ(e) , α1 , β1 ), η1 )e←e+1
end while(E )R̂ := Threshold(Aθ 1 , τ2 )(E )Aθ 1 ← ⃗0G×G(E1 )θ−Aθ ← RandomGaussianInit()
γ←0while e < E1 + E2 doθ(e+1) := AdamUpdate(θ(e) , ∇L1 (X, θ(e) , R̂, α2 , β2 , γ), η2 )γ ← γ + γ+e←e+1
end while(E )output DAGTrim(Aθ 2 , τ2 )Algorithm 2 DAGTrim
Require: Aθ ∈ RD×D , τ ∈ R+
E ← ∅ {Initialize the set of final edges.}C ← [(i, j) ∈ J1, dK2 | (Aθ,i,j > τ ] {Candidate edges above threshold τ .}Sort C by decreasing Aθ,i,j .
for each (i, j) ∈ C do
if the graph with edges E ∪ {(i, j)} is still acyclic thenE ← E ∪ {(i, j)} {We add the edge if it does not create a cycle.}end ifend for19

Stable Differentiable Causal Discovery

F. Empirical Studies Details
F.1. Simulation Details
We simulate observational and interventional data for a wide range of d (number of variables), varying the graph densitywith s (the average number of parents per node), and varying the number of variables that are intervened on. The simulations
proceed as done in Brouillard et al. (2020); Bello et al. (2022), by sampling a random graph, modeling its conditionals withrandom neural networks, setting its interventional distribution to Gaussian, and drawing samples from the obtained model.
In detail, we generated simulated data according to the following procedure:• Draw a random undirected graph from the Erdős-Rényi distribution.
• Convert the undirected graph into a DAG G∗ by setting the direction of each edge i → j if π(i) < π(j), where π is arandom permutation of the nodes.
• Form d possible sets of interventions that target one variable at a time: Ij = {j} and I0 = ∅.
G∗

• Draw a set of random fully connected neural networks MLP(j) : R| paj | → R100 → 1, each one with one 100dimensional hidden layer. Each neural network parametrizes the mean of the observational conditional distributions:p∗j (xj | xpaG∗ ; 0) ∼ N µ = MLP(j) (xpaG∗ ), σ = 0.5 .
jj• For intervention distribution k ≥ 1, perform a hard intervention on variable k and setp∗j (xk ; k) ∼ N (0, 0.1).
• Draw the data according to the model, with 10,000 observational samples and 500 extra interventional samples pertarget variable.
• Standardize the data.
We consider several values of d to simulate different scenarios.
F.2. Choice of Hyperparameters
Consistent with prior work (e.g., DAGMA, NOTEARS), we do not conduct hyperparameter optimization for the experiments.
Instead, we fix a single set of parameters for all experiments.
We fixed the hyperparameters as follows: α1 := 1e−2, β1 := 2e−4, η1 := 2e−3, τ1 := 0.2, α2 := 5e−4, β2 :=
5e−3, η2 := 1e−3, γ + := 0.005, τ2 := 0.1. We found that these selections worked well empirically across multiple
simulated datasets and were used in all experiments without simulation-specific fine-tuning.
Each stage was run for 2000 epochs with a batch size of 256, and the validation loss was computed over a held-out fractionof the training dataset (20% of the data) every 20 epochs for early stopping. In stage 2, the DAG check of the implicit
adjacency matrix was performed every 20 epochs before the validation loss computation.
F.3. Baseline Methods
For interventional data. For datasets with interventional data, we compare SDCD against DCDI (Brouillard et al., 2020),
DCD-FG (Lopez et al., 2022), and GIES (Hauser & Bühlmann, 2012).
For observational data. When the dataset contains only observational data, we include the interventional methods andfurther compare against NO-TEARS (Zheng et al., 2018), NO-BEARS (Lee et al., 2019), DAGMA (Bello et al., 2022), andSCORE (Rolland et al., 2022). In addition, we report sortnregress (Reisach et al., 2021), a trivial baseline that should beoutperformed (see Robustness Checks). We further included NOCURL (Yu et al., 2021) and AVICI (Lorch et al., 2022) inAppendix G.
Implementation details. Here, we provide details on the baseline methods and cite which implementations were used forthe experiments. For DCDI and DCDFG, we used the implementations from https://github.com/Genentech/20

Stable Differentiable Causal Discovery

dcdfg, using the default parameters for optimization. For DCDFG, we used 10 modules in our benchmarks, as reported in
the paper experiments. For GIES, we used the Python implementation from https://github.com/juangamella/gies, using the default parameters. For DAGMA, we used the original implementation from https://github.
com/kevinsbello/dagma with the default parameters. For NOTEARS, we used the implementation from https://github.com/xunzheng/notears, and for NOBEARS, we used the implementation from https://github.
com/howchihlee/BNGPU. For NOTEARS and NOBEARS, we found the default thresholds for determining the final
adjacency matrix performed poorly or did not return a DAG, so for each of these baselines, we followed the same proceduredescribed in Lopez et al. (2022): we find the threshold that returns the largest possible DAG via binary search. sortnregress(Reisach et al., 2021) is a trivial baseline meant to ensure that the causal graph cannot be easily inferred from the variancepattern across the variables. For this baseline, we used the implementation in https://github.com/Scriddie/
Varsortability. For all methods, the training time on CPU is measured on an AMD 3960x with 4-core per method; on
GPU on an AMD 3960x with 16-core and an Nvidia A5000.
F.4. Metrics
We evaluate performance using the structural Hamming distance (SHD) between the true DAG G∗ and each method’s
output graph. SHD is the standard metric reported in causal discovery. It quantifies the minimum number of edge additions,deletions, and reversals needed to transform one graph into the other. A lower SHD indicates better reconstruction of G∗ .
F.5. Robustness Checks
Previous works detailed common issues with the SHD metric (Tsamardinos et al., 2006) and data simulation processes(Reisach et al., 2021). Below, we discuss three categories of issues that commonly arise when evaluating causal discovery
methods and address each issue with a diagnostic metric.
Sparsity. Particularly when the true causal graph is sparse, SHD may favor sparser predictions since, in the extreme case,the empty graph achieves an SHD equal to the number of true edges. To show the relative performance of the benchmarked
methods with respect to this trivial solution, we indicate the number of true edges for each simulated setting in Fig. 2 andFig. 3. We find that most methods outperform this baseline. Additionally, we report the F1 score and the recall of thepredictions (see Figs. 7 and 10), two metrics that suffer when a method predicts many false negatives. We find that SDCD
still outperforms other methods with these metrics.
Identifiability. In settings with incomplete or no interventional data, the true causal graph may be unidentifiable, meaning
multiple I-Markov equivalent graphs can maximize the score (Brouillard et al., 2020). Therefore, graphs in the same
Markov equivalence class as the true causal graph may have positive SHD values despite being optimal with respect to theavailable data. As proposed in Peters et al. (2014), we also compute an adapted version of the SHD to compare the Markov
equivalence class of the methods’ results against the true Markov equivalence class instead of the graphs themselves. This
metric, called SHD-CPDAG, is computed as the SHD between the completed partially directed acyclic graph (CPDAG) of
the predicted graph and the CPDAG of the true graph. Unlike the regular SHD metric, this metric is zero if two graphs arein the same equivalence class. We report it alongside SHD for our experiments in Fig. 7 to better represent the results in
scenarios with an unidentifiable causal graph. We find very similar results.
Simulation issues. As discussed in Reisach et al. (2021), certain simulation processes used for causal discovery
benchmarking exhibit an issue where the order of the variables, when sorted by sample variance, reflects the true causal
ordering of the graph. As a result, methods that exploit this phenomenon to accurately infer the causal graph may be
misrepresented. To ensure that our simulation process does not suffer from this issue and that the methods are being properly
evaluated, we take two complementary steps recommended in Reisach et al. (2021): (i) we standardize the data before
being input into any of the evaluated methods so that no artificial sample variance information can be exploited, and (2)we include the trivial baseline, sortnregress (Reisach et al., 2021), which is designed to exploit sample variance artifacts
from a flawed simulation, and should be outperformed by an effective, scale-invariant algorithm. We find that sortnregress
performs poorly, which confirms that our normalization scheme removes simulation artifacts, and we find that SDCD and its
competing methods beat sortnregress by a wide margin.

21

Stable Differentiable Causal Discovery

F.6. Additional Experiments
Ablation Experiments. We performed ablation studies to judge the impact of each innovation implemented in SDCD. We
evaluated modifications of SDCD where (1) only stage 2 is performed without stage 1 and where (2) in stage 2, the spectral
constraint is substituted for alternative PST acyclicity constraints. As the results show in Tables 6 and 7, both stages are
essential to the success of SDCD.
Experiment against the best baseline. To further validate the results against the strongest baseline, we evaluate SDCD
on the simulated data generated in Brouillard et al. (2020) (DCDI) and compare our results against their reported SHD
values. In Supplementary Table 2 we report the results of SDCD on the simulated data presented in Brouillard et al. (2020)
alongside their original DCDI results. SDCD outperforms DCDI on all its sparse datasets (s = 1). Only for datasets whered = 10, s = 4, does SDCD perform worse than DCDI. However, we find the edge density (δ = 88.9%) of these graphs to
be unrepresentative of realistic scenarios.

22Stable Differentiable Causal Discovery

G. Supplementary Figures and Tables

Threshold to DAG

With Penalty

With Augmented Lagrangian

0.04
0.03

hexp
hlog
hinv

0.02
0.01
0.00

0k

2k

4k

Epoch

h
h+

6k

0k

50k

100k

Epoch

150k

Figure 5: The effect of constraints on the learned graph throughout training. The training with penalty hρ + (dashed purple,
exactly hρ with a hard mask on the diagonal as to prevent self-loops, as implemented in SDCD) converges the fastest towarda DAG. (left) training with h as a regularization penalty. (right) training with h as an augmented Lagrangian constraint.
Threshold to DAG is the smallest η at which all edges with weight > η form a DAG.

s

d

δ

MethodSDCD

DCDI-G

DCDI-DSF

1

10

22.2%

L
NL-Add
NL-NN

0.7±1.2
0.6±0.7
0.7±0.7

1.3±1.9
5.2±7.5
2.3±3.6

0.9±1.3
4.2±5.6
7.0±10.7

20

10.5%

L
NL-Add
NL-NN

1.4±3.4
4.1±3.0
3.0±2.5

5.4±4.5
21.8±30.1
13.9±20.3

3.6±2.7
4.3±1.98.3±4.1

10

88.9%

L
NL-Add
NL-NN

5.2±3.5
4.8±2.1
7.3±3.0

3.3±2.1
4.3±2.4
2.4±1.6

3.7±2.3
5.5±2.4
1.6±1.6

20

42.1%

L
NL-Add
NL-NN

18.8±10.5
18.0±7.3
14.9±1.923.7±5.635.2±13.2
16.8±8.7

16.6±6.4
26.7±16.9
11.8±2.1

4

Table 2: Means and standard deviations of SHD scores over simulations from Brouillard et al. (2020). The “Method”
column refers to the model used to simulate the causal relationships. “L” refers to linear model, “NL-Add” refers tononlinear, additive model, and “NL-NN” refers to nonlinear, non-additive (neural network) model. We refer to Brouillard
et al. (2020) for the simulation details. The results are reported alongside the values presented in the original paper. s refers
to the expected number of edges per node, d denotes the number of nodes, and the edge density, δ, is computed as thefraction of d(d−1), the maximum number of edges for a DAG. The lowest average SHD values are set in bold.
2

23

Stable Differentiable Causal Discovery

SDCD
DCDI

Training Time

1.5h

NOTEARSNOBEARS

DCDFG
GIES

1h

30m
0m

10

20

30

40

d

DAGMA

5h
4h
3h
2h
1h
0

50 50

SDCD-GPU

200 300 400 500dFigure 6: Training runtimes across simulations from Fig. 2. SDCD on GPU (dashed) scales to 500 variables in under 20
minutes. Error bars indicate std on 5 random datasets for d < 50 and 3 random datasets for d ≥ 50.
SDCD

SDCD-GPU

DCDI

DCDFG

GIES

DAGMA

NOTEARSNOBEARS

SCORE

sortnregress

AVICI

NOCURL

14.7 ±5.535.7 ±6.253.8 ±11.964.0 ±13.7
69.9 ±12.3
92.7 ±9.1
225.3 ±13.7
350.0 ±12.5
466.3 ±62.4
621.7 ±10.7

NT
NT
NT
NT
68.3 ±13.3
89.7 ±11.0
228.0 ±18.3
360.0 ±nan471.7 ±68.0
621.0 ±10.524.3 ±3.931.7 ±6.5
55.5 ±10.7102.4 ±21.6
NA
NA
NA
NA
NA
NA

24.6 ±6.0
108.0 ±14.3258.8 ±32.6
426.6 ±73.7
660.8 ±126.1
1807.3 ±788.2
5657.3 ±2982.6
7284.7 ±5072.3
3779.7 ±507.3
7252.7 ±3284.6

27.8 ±3.9123.4 ±12.3
NA
NA
NA
NA
NA
NA
NA
NA

25.3 ±6.2
62.0 ±12.0
89.4 ±10.9115.3 ±13.4
NA
NA
NA
NA
NA
NA

35.3 ±1.975.4 ±4.0
113.1 ±4.8147.9 ±6.0
183.4 ±7.4327.3 ±7.5619.0 ±4.2
NA
NA
NA

33.5 ±3.074.2 ±3.0113.7 ±3.3151.1 ±3.4192.0 ±3.5389.0 ±3.6770.0 ±7.81149.0 ±14.0
1534.7 ±3.1
1915.7 ±18.8

14.4 ±4.0118.6 ±8.5275.9 ±21.0454.3 ±52.8
619.4 ±59.7
NA
NA
NA
NA
NA

28.6 ±4.781.6 ±6.3
134.6 ±8.4172.9 ±12.2
216.6 ±12.4
421.3 ±12.0
824.0 ±19.0
1190.7 ±26.3
1585.0 ±59.3
1974.3 ±34.6

23.52
60.93
97.13
135.83170.83366.50
NT
NT
NT
NT

33.23
82.37
134.60179.93240.93513.07
NT
NT
NT
NT

d
10
20
30
40
50
100
200
300
400
500

Table 3: Detailed results of SHD means and standard deviations from Fig. 2. SDCD-GPU was only run for d ≥ 50. Allother NA values correspond to failed runs (possibly from timeout after 6h, e.g., GIES, or from training error, e.g., DCDI).
NT corresponds to the method not having been tested on that particular example.

1200

F1 Score

SHD-CPDAG

0.75

100
0

0.5

0.25

10

20

30
d

40

50

0

10

20

30
d

40

50

# of edges
SDCD
DCDI
NOTEARSNOBEARS
DCDFG
sortnregress
GIES
DAGMA
SCORE

Figure 7: F1 and SHD-CPDAG metrics across simulations from Fig. 2, observational data with increasing numbers ofvariables d. Missing data points imply the method failed to run. Error bars indicate std on 30 random datasets.

24

Stable Differentiable Causal Discovery# of edges

150m
120m
90m
60m
30m
0

Training Time

SDCD-GPU

SHD

16000
12000
8000
4000
0

1000 2000 3000 4000
d1000 2000 3000 4000dFigure 8: SDCD on GPU (dashed) scales to 4000 variables under 3 hours while maintaining competitive SHD. Error bars
indicate std on 3 random datasets for d = 1000, 2000 and 2 random datasets for d = 3000, 4000.

d

SDCD-GPU

1000
2000
3000
4000

1438.7 ±59.2
3356.7 ±70.0
5172.5 ±89.8
7567.0 ±343.7

Table 4: Detailed results of SHD means and standard deviations from Fig. 8.

25

Stable Differentiable Causal Discovery

In addition to SHD, we computed precision and recall metrics over the predicted edges with respect to the true edges forboth observational and interventional scenarios. The precision is the fraction of true edges among all the predicted edges.
The recall is the fraction of true edges that have been correctly predicted.

SDCD
DCDI

NOTEARSNOBEARS

DCDFG
sortnregress

GIES
DAGMA

SCORE

1.0
Precision

0.8
0.6
0.4
0.2
0.0

10

20

30
d

40

50

50

200 300 400 500dFigure 9: Precision across simulations from Fig. 2, observational data with increasing numbers of variables d. The
SDCD(-CPU) and SDCD-GPU lines overlap, indicating consistent results. Missing data points imply the method failed torun. Error bars indicate std on 30 random datasets for d < 50 and 5 random datasets for d ≥ 50.

SDCD
DCDI

NOTEARSNOBEARS

DCDFG
sortnregress

GIES
DAGMA

SCORE

1.0
Recall

0.8
0.6
0.4
0.2
0.0

10

20

30
d

40

50

50

200 300 400 500dFigure 10: Recall across simulations from Fig. 2, observational data with increasing numbers of variables d. The SDCD(CPU) and SDCD-GPU lines overlap, indicating consistent results. Missing data points imply the method failed to run. Error
bars indicate std on 30 random datasets for d < 50 and 5 random datasets for d ≥ 50.

26

Stable Differentiable Causal Discovery

1.0

1.0

0.8

0.8
RecallPrecision

SDCD-GPU

0.6
0.4
0.2

0.6
0.4
0.2

0.0

0.0
1000

2000

d

3000

4000

1000

2000

d

3000

4000

Figure 11: Precision and recall across simulations from Fig. 8, observational data with increasing numbers of variables d.
Error bars indicate std on 3 random datasets for d = 1000, 2000 and 2 random datasets for d = 3000, 4000.

27

Stable Differentiable Causal Discovery

s=2

Number of true edges = d x s

s=4

s=6

= Edge Density = 2s / (d-1)

s=8

d=20
d=30
d=40
Figure 12: SHD across simulations with an increasing proportion of variables intervened on, varying the total number ofvariables d (columns) and average edges per variable s (rows). Extended version of Fig. 3 with DCDFG and GIES ands = 8. Boxplots over 5 random datasets.

28

Stable Differentiable Causal Discovery

s

2

4

6

8

d

δ

20

21%

30

13%

40

10%

20

42%

30

27%

40

20%

20

63%

30

41%

40

30%

20

84%

30

55%

40

41%

Fraction of Variables
Intervened on

SDCD

DCDI

DCDFG

GIES

0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00

18.0 ±6.5
13.4 ±5.7
9.4 ±5.4
9.0 ±2.1
9.0 ±2.8
24.6 ±9.626.8 ±13.1
21.8 ±9.0
16.2 ±3.418.2 ±2.2
44.2 ±4.1
33.4 ±3.633.0 ±3.527.6 ±3.427.0 ±7.524.8 ±4.6
27.4 ±8.625.4 ±8.5
19.8 ±4.8
18.8 ±3.356.0 ±32.960.4 ±40.956.2 ±26.643.2 ±17.8
31.7 ±6.1
109.2 ±27.6123.8 ±62.1
105.6 ±38.8
76.0 ±24.4
63.0 ±26.9112.3 ±25.8
99.0 ±7.8109.0 ±31.476.0 ±14.2
98.0 ±19.7
183.0 ±102.1
204.3 ±51.4
206.7 ±84.3
207.0 ±119.7
283.7 ±15.2
465.0 ±131.8
372.7 ±84.3
479.7 ±103.5
469.3 ±126.6
333.7 ±135.5

80.3 ±26.1
70.0 ±36.8
69.7 ±41.962.7 ±22.0
65.7 ±18.5
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA

0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00

36.0 ±6.432.2 ±6.6
34.6 ±10.625.8 ±8.822.4 ±7.1
54.0 ±9.8
43.8 ±10.3
39.2 ±8.635.0 ±9.929.0 ±6.5
69.0 ±11.7
56.8 ±15.4
50.4 ±13.0
41.4 ±10.734.4 ±11.3

31.2 ±4.1
33.0 ±3.630.4 ±7.6
29.6 ±8.2
25.8 ±6.4
57.6 ±7.967.0 ±4.0
72.4 ±15.572.2 ±10.760.3 ±5.7
99.0 ±30.7
107.0 ±39.2
105.6 ±41.6
97.8 ±33.993.5 ±16.3

104.0 ±3.5105.7 ±4.2
107.7 ±9.5102.0 ±6.1
107.0 ±10.4
234.7 ±42.9269.3 ±36.7
262.0 ±15.7
232.7 ±24.0
236.0 ±42.8
460.0 ±47.8438.3 ±73.1
457.7 ±21.2
426.3 ±52.6458.7 ±55.3

110.7 ±25.3
110.7 ±20.0
101.3 ±22.3
105.7 ±22.8105.3 ±12.2
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA

0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00

51.2 ±7.544.0 ±6.542.2 ±8.6
38.8 ±8.3
34.0 ±5.185.4 ±6.579.8 ±12.5
69.4 ±14.8
67.0 ±12.2
57.4 ±11.3
95.4 ±27.7
75.6 ±23.6
63.6 ±17.0
57.4 ±19.6
47.2 ±12.2

56.6 ±10.437.8 ±10.2
29.0 ±10.825.2 ±10.0
23.8 ±7.775.8 ±4.469.2 ±15.4
80.6 ±17.2
86.2 ±23.3
82.0 ±5.3
107.8 ±36.3
130.2 ±43.3
146.0 ±51.7
128.3 ±16.3
114.5 ±31.8112.0 ±12.3
92.3 ±17.0
97.0 ±11.5
94.0 ±11.5
93.3 ±1.5256.7 ±30.6
260.7 ±12.4
257.3 ±18.6
245.3 ±23.7
259.7 ±19.5
460.3 ±63.3
401.0 ±112.8
370.0 ±87.1
387.7 ±81.1
469.0 ±28.2

117.7 ±11.9117.3 ±11.6
112.7 ±7.6
91.3 ±16.6
86.0 ±7.9
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA

0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00

53.6 ±10.2
51.0 ±8.1
47.0 ±5.340.8 ±6.7
43.0 ±9.0111.8 ±9.8
90.8 ±7.5
89.6 ±13.2
77.6 ±9.3
71.0 ±4.6150.4 ±16.9127.0 ±12.4
113.4 ±20.2
104.4 ±21.3
92.0 ±17.6

82.8 ±14.558.2 ±12.6
41.4 ±13.826.0 ±3.8
19.8 ±4.993.4 ±13.375.6 ±8.2
78.8 ±12.2
81.2 ±9.4
81.6 ±12.5151.0 ±23.1
188.4 ±27.4
200.0 ±33.9
193.0 ±14.1
190.0 ±nan117.7 ±11.4
108.3 ±19.7
100.7 ±7.673.7 ±8.1
81.7 ±22.4
255.0 ±27.1
242.7 ±14.6
255.0 ±22.1
226.3 ±15.1
222.3 ±13.0
450.0 ±44.5452.0 ±7.0424.3 ±43.0
426.7 ±34.6
422.0 ±18.1

111.7 ±16.5
96.7 ±15.8
91.0 ±23.5
86.0 ±31.2
79.7 ±10.8
NA
NA
NA
NA
NA
NA
NA
NA
NA
NA

Table 5: Detailed results of SHD means and standard deviations from Fig. 12. GIES failed to run on d ≥ 30.

29

Stable Differentiable Causal Discovery

s=2

Number of true edges = d x s

s=4

s=6

= Edge Density = 2s / (d-1)

s=8

d=20
d=30
d=40
Figure 13: Precision across simulations from Fig. 12 increasing proportion of variables intervened on, varying the total
number of variables d (columns) and average edges per variable s (rows). Boxplots over 5 random datasets.

30

Stable Differentiable Causal Discovery

s=2

Number of true edges = d x s

s=4

s=6

= Edge Density = 2s / (d-1)

s=8

d=20
d=30
d=40
Figure 14: Recall across simulations from Fig. 12 increasing proportion of variables intervened on, varying the total numberof variables d (columns) and average edges per variable s (rows). Boxplots over 5 random datasets.

31

Stable Differentiable Causal Discovery

Name

d=10

d=20

d=30

d=40

SDCD
SDCD-warm
SDCD-warm-nomask
SDCD-no-s1
SDCD-no-s1-2

14.7
14.7
19.3
19.3
16.3

40.3
40.7
69.7
68.3
56.7

54.3
55.0
156.0
155.3
95.0

69.0
68.7
272.7
272.3
135.0

DCDI

24.0

35.7

56.7

87.0

Table 6: Ablation study for SDCD stage 1. We observe that the described version of SDCD performs the best out of all
variations. SDCD-warm performs competitively but generally provides little benefit. SDCD-warm-nomask performs muchworse than SDCD, demonstrating that enforcing the mask during stage 2 is important. We report mean SHD over threerandom seeds of observational data (no interventions) with a fixed number of edges per variable, s = 4, for a range of
numbers of variables, d. SDCD-warm refers to starting stage 2 of SDCD, where the input layer is ported over from stage 1
instead of re-learned. SDCD-warm-nomask performs the same warmstart as SDCD-warm but does not enforce the mask instage 2. SDCD-no-s1 only performs stage 2. SDCD-no-s1-2 only does stage 2, but sets (α2 , β2 ) to the default values fromstage 1 (α1 , β1 ). We report these values alongside DCDI. The lowest SHD values are bolded for each value of d.

Name

d = 10d = 20d = 30d = 40d = 50SDCD
SDCD-exp
SDCD-log

13.33
11.60
11.20

33.47
44.33
52.00

54.07
69.07
87.47

70.80
85.07
117.8776.60
89.93
116.00

Table 7: Ablation study for SDCD stage 2 (choice of the constraint). We observe that SDCD performs the best out of thethree variations. Additionally, the variations using the PST constraints do not crash for any of the runs, even for those withd = 50. We attribute this improved stability (as compared to DCDI and DAGMA) to stage 1 since there are fewer non-zeroparameters contributing to the value of the constraint. We report mean SHD over five random seeds of observational data(no interventions) with a fixed number of edges per variable, s = 4, for a range of numbers of variables, d. SDCD-exp is
SDCD except using the hexp constraint in place of the hρ , and SDCD-log uses the hlog constraint.

MethodStableTraining

Scalable
Constraint

Can UseInterventions

Expressive
Model Class

SDCD

✓

O(d2 )

✓

✓

✓
✓
✕
✕
✕

✓
✕
✓
✕
✕

DCDI
DCDFG
DAGMA
NO-TEARS
NO-BEARS

✕
✕
✕
✕
✓

3

O(d )
O(md)
O(d3 )O(d3 )O(d2 )

Table 8: Comparison of Differentiable Causal Discovery methods including our proposed SDCD method. Expressive model
class refers to the capability to approximate any causal graph with non-linear structural equations.

32