Diffusion Models for Causal Discovery via Topological Ordering
Authors: Pedro Sanchez
Xiao Liu
Alison Q O'Neil
Sotirios A. Tsaftaris
Venue: ICLR 2023 poster (2023)
Keywords: Diffusion Models
Causal Discovery
Topological Ordering
Score-based Methods

Abstract: Discovering causal relations from observational data becomes possible with additional assumptions such as considering the functional relations to be constrained as nonlinear with additive noise (ANM). Even with strong assumptions, causal discovery involves an expensive search problem over the space of directed acyclic graphs (DAGs). \emph{Topological ordering} approaches reduce the optimisation space of causal discovery by searching over a permutation rather than graph space.
For ANMs, the \emph{Hessian} of the data log-likelihood can be used for finding leaf nodes in a causal graph, allowing its topological ordering. However, existing computational methods for obtaining the Hessian still do not scale as the number of variables and the number of samples are increased. Therefore, inspired by recent innovations in diffusion probabilistic models (DPMs), we propose \emph{DiffAN}, a topological ordering algorithm that leverages DPMs for learning a Hessian function. We introduce theory for updating the learned Hessian without re-training the neural network, and we show that computing with a subset of samples gives an accurate approximation of the ordering, which allows scaling to datasets with more samples and variables. We show empirically that our method scales exceptionally well to datasets with up to $500$ nodes and up to $10^5$ samples while still performing on par over small datasets with state-of-the-art causal discovery methods.
Implementation is available at \url{https://github.com/vios-s/DiffAN} .

Published as a conference paper at ICLR 2023
DIFFUSION MODELS FOR CAUSAL DISCOVERY VIA TOPOLOGICAL ORDERING
Pedro Sanchez Xiao Liu', Alison Q O'Neil Sotirios A. Tsaftaris" 1.3
The University of Edinburgh
Cannn Medical Research Europe
The Alan Turing Institute
ABSTRACT
Discovering causal relations from observational data becomes possible with addi- tional assumptions such as considering the functional relations to constrained as nonlinear with additive noise (ANM). Even with strong assumptions, causal discovery involves an expensive search problem over the space of directed acyclic graphs (DAGs). Topological ordering approaches reduce the optimisation space of causal discovery by searching over permutation rather than graph space. For ANMs, the Hessian of the data log-likelihooe can used for finding leaf nodes in causal graph, allowing its topological ordering. However, existing compu- tational methods for obtaining the Hessian still do not scale as the number variables and the number of samples are increased. Therefore, inspired by re- cent innovations in diffusion probabilistic models (DPMs), we propose DijfAAN topological ordering algorithm that leverages DPMs for learning Hessian fune- tion. We introduce theory for updating the learned Hessian without re-training the neural network, and we show that computing with a subset of samples gives ac- curate approximation of the ordering, which allows scaling to datasets with more samples and variables. We show empirically that our method scales exceptionally well to datasets with up to 500 nodes and up to 10 samples while still performing on par over small datasets with sfate-of-the-art causal discovery methods.
40000
method
INTRODUCTION
Runtime (seconds)
DifIAN
30000
SCORE
Understanding the causal structure of problem important for areas such as economics, biology (Sachs et al. 2005) and healthcare (Sanchez et al. (0022) es-
20000
10000
pecially when reasoning about the ef-
fect of interventions. When interven- tional data from randomised trials are not available, causal discovery meth- ods (Glymour et al., 2019) may be em- ployed to discover the causal structure of a problem solely from observational data. Causal structure is typically mod- elled as a directed acyclic graph (DAG) G in which each node is associated with a random variable and each edge repre- sents causal mechanism i.e. how one variable influences another.
10a
10
10
10
Sample Size
Figure Plot showing run time in seconds for differ- ent sample sizes, for discovery of causal graphs with 500 nodes. Most causal discovery methods have pro- hibitive run time and memory cost for datasets with many samples; the previous state-of-the-ar SCORE algorithm RRolland al. 2022) which is included in this graph can- not be computed beyond 2000 samples in a machine with 64GB of RAM. By contrast,
method DiffAN has reasonable run time even for numbers of samples two or- ders of magnitude larger than capable by most existing methods.
However, learning such model from data NP-hard (Chickering 1996). Tra- ditional methods search the DAG space
ppcdrossanched @ ed.ac.uk
Implementation is available at http:://ggccbb..cm/////////////e 0

---

Published as a conference paper at ICLR 2023
by testing for conditional independence between variables (Spirtes al. 1993) or by optimising some goodness of fit measure (Chickering 2002). Unfortunately, solving the search problem with a greedy combinatorial optimisation method can be expensive and does not scale to high- high-dimensional
problems.
In line with previous work (Teyssier & Koller 2005 Park & Klabjan 2017 Bühlmann al. 2014
Solus et al. 2021 Wang et al. [202]] Rolland et al.] (0022). we can speed up the combinatorial search problem over the space of DAGs by rephrasing it as a lopological ordering task, ordering from leaf nodes to root nodes. The search space over DAGs with d nodes and (d² d)/2 possible edges is much larger than the space of permutations over variables. Once a topological ordering of the nodes is found, the potential causal relations between later (cause) and earlier (effect) nodes can be pruned with a feature selection algorithm (e.g. [Bühlmann et al. (2014)) to yield a graph which is naturally directed and acyclic without further optimisation.
Recently, Rolland et al.) (2022) proposed the SCORE algorithm for topological ordering- SCORE uses the Hessian of the data log-likelihood, V? logg(x))BB verifying which elements of V pog) diagonal are constant across all data points, leaf nodes can be iteratively identified and removed. [Rolland et al. (2022) estimate the Hessian point-wise with second-order Stein gra- dient estimator (Li & Turner 2018) over a radial basis function (RBF) kernel. However, point-wise estimation with kernels scales poorly to datasets with large number of samples because it requires inverting n x n kernel matrix.
Here, we enable scalable causal discovery by utilising neural networks (NNs) trained with denois- ing diffusion instead of Rolland et al.'s kernel-basee based estimation. We use the ordering procedure, based on Rolland et al. 2022). which requires re-computing the score's Jacobian at each iteration. Training NNs at each iteration would not be feasible. Therefore, we derive a theoretical analysis that allows updating the learned score without re-training. In addition, the NN is trained over the entire dataset (n samples) but only subsample is used for finding leaf nodes. Thus, once the score model learned, we can use it to order the graph with constant complexity on n, enabling causal discovery for large datasets in nigh-dimensional settings. Interestingly, our algorithm does not require archi- tectural constraints on the neural network, as previous causal discovery methods based on neural networks (Lachapelle et al., 2020; Zheng et al., 2020 Yu et al., 2019; Ng al., (022). Our training procedure does not learn the causal mechanism directly, but the score of the data distribution.
Contributions. In summary, we propose DiffAN, an identifiable algorithm leveraging a diffusion probabilistic model for topological ordering that enables causal discovery assuming an additive noise model: (i) To the best of our knowledge, we present the first causal discovery algorithm based on denoising diffusion training which allows scaling to datasets with up to 500 variables and 10 sam- ples. The score estimated with the diffusion model is used to find and remove leaf nodes iteratively; (ii) We estimate the second-order derivatives (score's ' Jacobian or Hessian) of a data distribution us- ing neural networks with diffusion training via backpropagation; (iii) The proposed deciduous score (Section 3) allows efficient causal discovery without re-training the score model at each iteration. When a leaf node removed, the score of the new distribution can be estimated from the original score (before leaf removal) and its Jacobian.
2 PRELIMINARIES
2.1 PROBLEM DEFINITION
We consider the problem of discovering the causal structure between d variables, given a probability distribution p(x) from which a d-dimensional random vector x = (XI....,Xx) Xa) can be sampled. We assume that the true causal structure is described by DAG & containing d nodes. Each node represents a random variable x; and edges represent the presence of causal relations between them. In other words, we can say that G defines structural causal model (SCM) consisting of collection of assignments x, == (P())))), Ez), where Pa(x;) are the parents of x; in G, and is noise term independent of x, also called exogenous noise, 6; are i.i.d. from a smooth distribution p²_ The SCM entails a unique distribution p(x) = I[__ p(x; Pa(x;)) over the variables x (Peters et al. 2017). The observational input data are X E 2n where n is number of samples. The target output is an adjacency matrix A E Rdd
2

---

Published as a conference paper at ICLR 2023
The topological ordering (also called causal ordering causal list) of DAG G is defined non-unique permutation и of nodes such that a given node always appears first in the list than its descendants. Formally, < node in G (Appendix B Peters et al. (201)))
E Dxg(x;) where Dxg(x;) are the descendants of the ith
2.2 NONLINEAR ADDITIVE NOISE MODELS
Learning unique A from X with observational data requires additional assumptions. A common class of methods called additive noise models (ANM) (Shimizu et al. 2006 Hoyer et al., 2008 Peters 2014 Bühlmann al. 2014 explores asymmetries the data by imposing functional assumptions on the data generation process. In most cases, they assume that assignments take the form x; (P((xx:)) + 4 with c p². Here we focus on the case described by Peters et al. (014)) where fi is nonlinear. We use the notation ft for fi (((()())) because the arguments of Ji will always be Pa(x,) throughout this paper. We highlight that f; does not depend on
Identifiability. We assume that the SCM follows additive noise model (ANM) which is known to be identifiable from observational data (Hoyer et al. 2008 Peters et al.) 2014). We also assume common cause of at least two observed variables. In addition, corollary 33 from Peters et al. (2014) states that the true topological ordering of the DAG, in our setting, is identifiable from p(x) generated by ANM without requiring causal minimality assumptions.
causal sufficiency, i.. there are hidden variables that
Finding Leaves with the Score. Rolland al. (2022) propose that the score of an ANM with distribution p(x) can be used to find leaves Before presenting how to find the leaves, we derive,
following Lemma in Rolland al. (2022) an analytical expression for the score which can be
written as
Vx, logp(x)= = V. log II
p(x: Pa(x:))
1=1
V
log p(x: Pa(xx))
(1)
log (x; %--fi)
b Using G; =x fi
0 log (x₃ fz)
aʃ. a p' (x; Si)
0x,
0x,
iECh(x))
Where Ch(x;) denotes the children of Xj- We now proceed, based on Rolland et al. (0022) to derive condition which can be used to find leaf nodes.
Lemma 1. Given a nonlinear ANM with a noise distribution p' and leaf node j: assume that 8 log dx
a, where a constant, then
Var H1gg = =0
See proof in Appendix A.I
Remark. Lemma[eenbbler findingg leaf nodes based on the diagonal the the logllikelssood's Hessian.
Rolland et al. (2022), using a similar conclusion, propose a topological ordering algorithm that iteratively finds and removes leaf nodes from the dataset. At each iteration Rolland al. (2022) re-compute the Hessian with kernel-based estimation method. In this paper, we develop a more efficient algorithm for learning the Hessian high-dimension and for a large number of samples. Note that Rolland et al. (2022) prove that Equation 2 can identify leaves in nonlinear ANMs with Gaussian noise. We derive formulation which, instead, requires the second-order derivative of
the noise distribution to be constant. Indeed, the condition
8 logg'' = a is true for p" following dzd Gaussian distribution which is consistent with Rolland et (2022), but could potentially be true for other distributions as well.
We refer to nodes without children in DAG 9 leaves.
3

---

Published as a conference paper at ICLR 2023
2.3 DIFFUSION MODELS APPROXIMATE THE SCORE
The process of learning to denoise (Vincent 2011 can approximate that of matching the score [Hyvärinen] 2005). A diffusion process gradually adds noise to data distribution over time. Dif- fusion probabilistic models (DPMs) Sohl-Dickstein et al. 2015); Ho et al. (2020); Song et al. (202T) learn to reverse the diffusion process, starting with noise and recovering the data distr- bution. The diffusion process gradually adds Gaussian noise, with a time-dependent variance at, to a sample Xo Piaa (x) from the data distribution. Thus, the noisy variable x,, with [0, is learned to correspond to versions of Xo perturbed by Gaussian noise following t E
p (x; xo) N (xL: (xxi Varxo; (1 ax) I), where a IT (1 By). Bg is the variance sched- uled between [Smin: Smax] and I is the identity matrix. DPMs (Ho et al. 2020) are learned with weighted sum of denoising score matching objectives at different perturbation scales with
0t = arg min Ex..t. (A(E) [x,, -c
(3)
where Xx (o vI cre, with Xn ~ p(x) being sample from the data distribution, u (O, and N~ (0, I) is the noise. X(t) is a loss weighting term following Ho et al. 2200))
Remark. Throughoul this paper, we leverage the fact that the trained model eg approximates the score Vx, p(g)(x) of the data (Song Ermon 2019).
THE DECIDUOUS SCORE
Discovering the complete topological ordering with the distribution's Hessian (Rolland et al. 2022) is done by finding the leaf node (Equation 2). appending the leaf node x; to the ordering list и and removing the data column corresponding to X₁ from before the next iteration d- 1 times. Rolland al. (022) estimate the score's Jacobian (Hessian) at each iteration.
Instead, we explore an alternative approach that does not require estimation of new score after each leaf removal. In particular, we describe how to adjust the score of a distribution after each leaf removal, terming this a "deciduous score'] We obtain an analytical expression for the deciduous score and derive a way of computing it. based on the original score before leaf removal. In this sec- tion, we only consider that p(x) follows a distribution described by an ANM, we pose no additional assumptions over the noise distribution.
Definition 1. Considering a DAG & which entails a distribution p(x) = II p(x: Pa((x)). Let p(x-i) p(xi (a(s)) be p(x) without the random variable corresponding to the leaf node X5. The deciduous score V log E gd-1 is the score of the distribution p(x-1).
p(x)
Lemma 2. Given ANM which entails a distribution p(x), we can use Equation / to find an- alytical expression for an additive residue Di between the distribution' score log l(x)(x) and its deciduous score V log p(x-1) such that
A Vlogp(x) V logp(x-)).
(4)
In particular, A is vector {ō₂ denoted as
Vj E [1,
d (1 where the residue w.E.t a node x, can be
55 = Vx, logg(x) Vx, (ogp(x
Of. Olog pf (x, fi) 0x, Òx
(5)
x/x, $ Pa(xi). 83 ==:0. =0
Proof. Observing Equation [ the score Vx, log ((g)(x) only depends on the following random vari- ables (i) P.(x,). (ii) Ch(x,). and (ii) Pa(Ch(xx)). We consider x to be a leaf node, therefore Vx, log p(x) only depends on if x, E Pa(xı). If x, Pa(x)), the only term depending on Vx, log p(x) dependent on x, one of the terms inside the summation.
An analogy deciduous which seasonally shed leaves during autumn.
4

---

Published as a conference paper at ICLR 2023
Find nodes
variables
C
D
0.7
0.1
03
49
2.5
0,4
0.7
2.1
0.5
0.1
0.3
samples
GT edges Spurious edges
X M
repeat d times
Figure 2: Topological ordering with diffusion models by iteratively finding leaf nodes. At each iteration, one leaf node is found using Equation 9 In the subsequent iteration, the previous leaves are removed, reducing the search space. After topological ordering (as illustrated on the right side), the presence of edges (causal mechanisms) between variables can be inferred such that parents of each variable are selected from the preceding variables the ordered list. Spurious edges can be pruned with feature selection post-processing step (Bühlmann al. 2014 Lachapelle et al. 2020 Rolland et al. 2022).
However, we wish to estimate the deciduous score V logp(x-1) without direct access to the function fi. to its derivative, nor the distribution p¹. Therefore, we now derive an expression for A using solely the score and the Hessian of log p(x).
Theorem 1. Consider an ANM of distribution p(x) with score V log)(x) and the score's Jaco- bian H (log p(x)). The additive residue A₁ necessary for computing the deciduous score (as in Proposition 2) can be estimated with
Vx, log
A = (logp((x))
(6)
(1,((og
See proof in Appendix A.2
4 CAUSAL DISCOVERY WITH DIFFUSION MODELS
DPMs approximate the score of the data distribution (Song & Ermon 2019). In this section, we explore how to use DPMs to perform leaf discovery and compute the deciduous score, based on Theorem for iteratively finding and removing leaf nodes without re-training the score.
4.1 APPROXIMATING THE SCORE'S ' JACOBIAN VIA DIFFUSION TRAINING
The score's Jacobian can approximated by learning score 6g with denoising diffusion training of neural networks and back- propppppggg (Rumelhart et al. 1986 from the output to the input variables. can be written, for an input data point æ E Rd as
H., p(x) ss Vz,eo(z,t).
where Ee (z, t) means the ith output Eø is backpropagated to the jth input. The diagonal of the Hessian in Equation can, then, be used for finding leaf nodes as in Equation
In a two variable setting, it is sufficient for causal discovery to (i) train a diffusion model (Equation 3b; (ii) approximate the score's Jacobian via backpropagation (Equation 7 (iii) compute variance of
'The Jacobian of neural network can be efficiently computed with auto-differentiation libraries such as
funetorch (Horace Hc 2021.

---

Published as a conference paper at ICLR 2023
the diagonal across all data points: (iv) identify the variable with lowest variance effect (Equation 2 We illustrate in Appendix the Hessian of a two variable SCM computed with a diffusion model.
4.2 TOPOLOGICAL ORDERING
When a DAG contains more than two nodes, the process of finding leaf nodes (i.e. the topological order) needs to be done iteratively as illustrated in Figure 2 The naive (greedy) approach would be to remove the leaf node from the dataset, recompute the score, and compute the variance of the new distribution's 's Hessian identify the next leaf node (Rolland al.) 2022). Since we employ diffusion models to estimate the score, this equates to re-training the model each time after a leaf is removed.
We hence propose a method to compute the deciduous score V Plogg((x)) using Theorem to remove leaves from the initial score without re-training the neural network. In particular, assuming that leaf X₁ is found, the residue A can be approximated with with
c0(2, ()
(a = Vico(æ,
(8)
V1,co(1,
where (a () output corresponding the leaf node. Note that the term Vieg(æ; () is vector of size d and the other term is a scalar. During topological ordering, we compute A., which the summation of A over all leaves already discovered and appended to Naturally, we only compute A w.r.t. nodes Xj & π because Xj π have already been ordered and are not taken into account anymore.
In practice, we observe that training Eg on X but using a subsample B E Igd of size k randomly sampled from X increases speed without compromising performance (see (ection|4.3) In addition, analysing Equation 5) the absolute value of the residue 6 decreases if the values of are set to zero once the leaf node discovered. Therefore, we apply a mask M, E 1 over leaves discovered in the previous iterations and compute only the Jacobian of the outputs corresponding to x-1. Mg is updated after each itcration based on the ordered nodes т. We then find leaf node according to
leaf arg min Var V (score(M, @ B,t))),
x
where is DPM trained with Equation See Appendix B.3 for the choice of t. This topological means that we only consider ordering procedure is formally described in Algorithm score(- (- π) the outputs for nodes xy &
Algorithm 1: Topological Ordering with DiffAN
Input: X Raxd trained diffusion model es. ordering batch size ic [] A- = 0exd M, = 1 Exd score = Eg
while ||| * d do
B b X
// Randomly sample batch of k elements // Mask removed leaves // Sum of Equation g over
B -- B o M,
A GetA, (score, B)
score = score(-a) As
// Update score with residue // Equation
leaf GetLeaf B)
π [leaf, T] M.inaf
// Append leaf to ordered list
// Set discovered leaf to zero
end
Output: Topological order π
'The diffusion model itself is approximation of the score, therefore its gradients approximations of the derivatives.

---

Published as a conference paper at ICLR 2023
4.3 COMPUTATIONAL COMPLEXITY AND PRACTICAL CONSIDERATIONS
We now study the complexity of topological ordering with DiffAN w.r.t. the number of samples and number of variables d in dataset. In addition, we discuss what are the complexities of a greedy version as well as approximation which only utilises masking.
Complexity on n. Our method separates learning the score E from computing the variance of the Hessian's 's diagonal across data points, in contrast to Rolland et al. (2022). We use all samples in X for learning the score function with diffusion training (Equation 3 does not involve expensive constrained optimisation techniques" and we train the model for a fixed number of epochs (which is linear with n) or until reaching the early stopping criteria. We use a MLP that grows in width with d but it does not significantly affect complexity. Therefore, Moreover, Algorithm i computed over a batch B with size k < instead of the entire dataset X, as described in Equation 9 Note that the number of samples
consider training be (n))
B can be arbitrarily small and constant for different datasets. In Section 5.2 we verify that the accuracy of causal discovery initially improves as is increased but eventually tapers off.
Complexity on d. Once is trained, a topological ordering can be obtained by running Va (2, t) d times. Moreover, computing the Jacobian of the score requires back propagaaggg the gradients d - times, where i is the number of nodes already ordered in a given iteration. Finally, computing the deciduous score's residue (Equation K means computing gradient of the nodes. Resulting in a complexity of O(d (d final topological ordering complexity is therefore Onn + dd).
i) with i varying from 0 to d which can be described by O(d)). The
DiSTAN Masking, We verify empirically that the masking procedure described in Section 4.2 can significantly reduce the deciduous score's 's residue absolute value while maintaining causal discov- ery capabilities. In DiffAN Masking, do not re-train the Eg compute the deciduous score. This ordering algorithm is an approximation but has shown to work well in practice while showing remarkable scalability. DiffAN Masking has Onn + d²) ordering complexity.
5 EXPERIMENTS
In our experiments, we train a NN with a DPM objective to perform topological ordering and follow this with a pruning post-processing step (Bühlmann et al. 2014). The performance is evaluated on synthetic and real data and compared state-of-the-uue causal discovery methods from observational data which are either ordering- based or gradient-hasee methods, NN architecture. We use 4- layer multilayer perceptron (MLP) with LeakyRel.U and layer normalisation. Metrics. We use the structural Hamming distance (SHD). Structural Intervention Distance (SID) (Peters & Bühlmann) 2015), Order Divergence (Rolland et al. 2022) and run time in seconds. See Appendix [D3 for details of each metric. Baselines. We use CAM (Bühlmann et al., 2014). GranDAG (Lachapelle (et al. 2020) and SCORE (Rolland et al.) 2022). We apply the pruning procedure Bühlmann et al. (2014) all methods. See detailed results in the Appendix [ Experiments with real data from Sachs (Sachs et al. 2005) and SynTReN (Van den Bulcke et al. 2006) datasets are in the Appendix
5.1 SYNTHETIC DATA
In this experiment, we consider causal relationships with f. being a function sampled from a Gaus- sian Process (GP) with radial basis function kernel of bandwidth one. We generate data from additive noise models which follow Gaussian, Exponential or Laplace distributions with noise scales in the
intervals [0..4, 4,0.40000 0.8). [0.8, 1.2], [1. 11). which are known to be identifiable (Peters et al. 2014). The causal graph is generated using the Erdös-Rényi (ER) (Erdós al.) 1960)) and Scale Free (SF) (Bol- lobás et af., (203) models. For a fixed number of nodes d, we vary the sparsity of the sampled graph by setting the average number of edges be either d or 5d. We use the notation [d][graph typel[sparsity] for indicating experiments over different synthetic datasets. We show that DiffAN performs on par with baselines while being extremely fast, see Figure 3 We also explore the role of overfitting in Appendix E.2 the difference between DiffAN with masking only and the greedy
Suuc as the Augmented Lagrangian method (Zheng et al. 2018 Lachapelle al.
2020
7

---

Published as a conference paper at ICLR 2023
100
method
80
GranDAG
400
SHD
60
SCORE
runtime
DifIAN
40
200
20
20ER1
20SFI
20ERS
20SF5
20ER1
20SF1
20ER5
20SF5
Experiment
Experiment
Figure 3: SHD (left) and rum time in seconds (right) for experiments of synthetic data graphs for graphs with 20 nodes. The variation in the violinplots come from 3 different seeds over dataset generated from 3 different noise type and different noise scales. Therefore, we run a total of 27 experiments for each method and synthetic datasets type
version in Appendix E.4 how we choose t during ordering in Appendix E.3 and we give results stratified by experiment in Appendix E.5
5.2 SCALING WITH DIFFAN MASKING
We now verify how DiffAN scales to bigger datasets, in terms of the number of samples n. Here, we use DiffAN masking because computing the residue with DiffAN would be too expensive for very big d. We evaluate only the topological ordering, ignoring the final pruning step. Therefore, the performance will be measured solely with the Order Divergence metric.
Scaling to large datasets. We evaluate how DiffAN compares to SCORE (Rolland et al. 2022). the previous state-of-the-arr, in terms of run time (in seconds) and and the performance (order diver- gence) over datasets with d 500 and different sample sizes 102 ..... 10ah the error bars results across dataset (different samples of ER and SF graphs). As illustrated in Figure DiffAN the more tractable option the size of the dataset increases. SCORE relies inverting very large n x n matrix which is expensive in memory and computing for large n. Running SCORE for d 500 and 2000 is intractable in machine with 64Gb of RAM. Figure 4 (left) shows that, since DiffAN can learn from bigger datasets and therefore achieve better results as sample size increases.
Ordering batch size. An important aspect of our method, discussed Section that allows scalability in terms of the separation between learning the score function Eø and computing the Hessian variance across a batch of size k, with k < n. Therefore, we show empirically, as illustrated Figure ( (right), that decreasing is does not strongly impact performance for datasets with d E 10, 20, 50.
6 RELATED WORKS
Ordering- based Causal Discovery. The observation that causal DAG can be partially represented with a topological ordering goes back to Verma & Pearl (1990). Searching the topological ordering space instead of searching over the space of DAGs has been done with greedy Markov Chain Monte
Carlo (MCMC) (Friedman & Koller 2003), greedy hill-climbing search (Teyssier & Koller 2005). search (Park & Klabjan 2017). restricted maximum likelihood estimators (Bühlmann et al.
2014). sparsest permutation Raskutti Uhler 2018 Lam et al. 2022 [Solus et al., (021), and reinforcement learning (Wang al. 2021). In linear additive models, Ghoshal & Honorio (2018);
Chen et al. (2019) propose an approach, under some assumptions on the noise variances, to discover the causal graph by sequentially identifying leaves based on estimation of the precision matrix.
Hessian the Log-likelihodd. Estimating H (loog (x)) is the most expensive task the ordering algorithm. Our baseline (Rolland et al. 2022) propose extension of & Turner (2018) which
8

---

Published as a conference paper at ICLR 2023
1400
method DifIANM SCORE
Order Divergence
Order Divergence
1200
10
20
50
1000
10²
10
10
10"
10
Sample Size
Ordering Batch Size
Figure 4: Accuracy of DiffAN as the dataset size is scaled up for datasets with 500 variables and increasing numbers of data samples n (left) and as the batch size for computing the Hessian variance is changed (right). We show 95% confidence intervals over 6 datasets which have different graph structures sampled different graph types (ER/SF).
utilises the Stein's 's identity over a RBF kernel (Schölkopf & Smola 2002). Rolland al.js method cannot obtain gradient estimates at positions out of the training samples. Therefore, evaluating the Hessian over subsample of the training dataset is not possible. Other promising kernel-based approaches rely spectral decomposition (Shi et al., 2018) solve this issue and can be promising future directions. Most importantly, computing the kernel matrix is expensive for memory and computation
There are, however, methods (Achlioptas et al. 2001; Halko et al. 2011; Si jet al. 2017) that help scaling kernel techniques, which were not considered in the present work. Other approaches also possible with deep likelihood methods such normalizing flows (Durkan iet al. [2019] Dinh et al. (017)) and further compute the Hessian via backpropagation. This would require two backpropagation passes giving O(d) complexity and be less scalable than denoising diffusion. Indeed, preliminary experiments proved impractical in our high- digh-dimensional settings.
We use DPMs because they can efficiently approximate the Hessian with a single backpropagation pass and while allowing Hessian evaluation on a subsample of the training dataset. It has been shown (ong Ermon 2019) that denoising diffusion can better capture the score than simple denoising
Vincent 2011) because noise at multiple scales explore regions of low data density.
CONCLUSION
We have presented a scalable method using DPMs for causal discovery. Since DPMs approximate the score of the data distribution, they can be used efficiently compute the log-likelihood's ' Hessian by backpropagating each element the output with respect to each element the input. The decid- hous score allows adjusting the score remove the contribution of the leaf most recently removed, avoiding re-training the NN. Our empirical results show that neural networks can be efficiently used for topological ordering in high-dimensional graphs (up to 500 nodes) and with large datasets (up to 10t samples).
Our deciduous score can be used with other Hessian estimation techniques as long as obtaining the score and its full Jacobian is possible from a trained model, e.g. sliced score matching (Song al. (020)) and approximate backpropagation (Kingma & Cun 2010p. Updating the score is more practical than re-training in most settings with neural networks. Therefore, our theoretical result enables the community to efficiently apply new score estimation methods to topological ordering. Moreover, DPMs have been previously used generative diffusion models in the context of causal es- timation (Sanchez & Tsaftaris 2022). In this work, we have not explored the generative aspect such as Geffner et al. (2022) does with normalising flows. Finally, another promising direction involves constraining the NN architecture as in Lachapelle et al. (2020) with constrained optimisation losses (Zheng et al. 2018).

---

Published as a conference paper at ICLR 2023
8 ACKNOWLEDGEMENT
This work was supported by the University of Edinburgh, the Royal Academy of Engineering and Canon Medical Research Europe via P. Sanchez's PhD studentship. S.A. Tsaftaris acknowledges the support of Canon Medical and the Royal Academy of Engineering and the Research Chairs and Senior Research Fellowships scheme (grant RCSRF1819 (825).
REFERENCES
Dimitris Achlioptas, Frank Mcsherry, and Bernhard Schölkopf. Sampling Tech- niques for Kernel Methods. In T Dietterich, S Becker, and Z Ghahramani (eds.), Advances in Neural Information Processing Systems, volume 14. MIT Press, 2001.
URL https://prooced ings neurips. cc.ppaper/2001//111 VTcST666of11667772222442266----3 pdd
Béla Bollobás, Christian Borgs, Jennifer Chayes, and Oliver Riordan. Directed Scale-Free Graphs. In Proceedings of the Fourteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 003, pp. 132-139, USA, 2003, Society for Industrial and Applied Mathematics.
Peter Bühhmaan, Jonas Peters, and Jan Ernest. CAM:
Causal Additive Models, High-dimensinn Order Search And Penalized Regression. The Annals of Statistics, 42(6):252622566, 2014.
Wenyu Chen, Mathias Drton, and Y Samuel Wang. On causal discovery with an equal vauia-vrriance assumption. Biometrika, 106(4):973-98, 12 2019.
David Maxwell Chickering. Learning Bayesian Networks is NP-Complete. In Doug Fisher and Hans-J. Lenz (eds.), Learning from Data: Artificial Intelligence and Statistics V, pp. 121-130. Springer New York, New York, NY, 1996.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3(Nov):507-554, 2002.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using Real NVP. In- ternational Conference on Learning Representations, 12 2017.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural Spline Flows. In Proc. Advances in Neural Information Processing Systems (NeurIPS), 2019.
Paul Erdös, Alfréd Rényi, and others. On the evolution of random graphs. Publ, Math. Inst. Hung. Acad. Sei, 51): 1960.
Nir Friedman and Daphne Koller. Being Bayesian About Network Structure. A Bayesian Approach to Structure Discovery in Bayesian Networks. Machine Learning, 50(1))95-125, 2003.
Tomas Geffner, Javier Antoran, Adam Foster, Wenbo Gong, Chao Ma, Emre Kiciman, Amit Sharma, Angus Lamb, Martin Kukla, Nick Pawlowski, Miltiadis Allamanis, and Cheng Zhang. Deep End- to-end Causal Inference. 2022.
Asish Ghoshal and Jean Honorio. Learning linear structural equation models in polynomial time and sample complexity. In Amos Storkey and Fernando Perez-Cruz (eds.), Proceedings of the Twenty- First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pp. 1466-1475. PMLR, 8 2018.
Clark Glymour, Kun Zhang, and Peter Spirtes. Review of Causal Discovery Methods Based on Graphical Models. Frontiers in Genetics, 10(JUN):244, 12 2019.
Halko, Martinsson, and A Tropp. Finding Structure with Randomness: Probabilistic Al- gorithms for Constructing Approximate Matrix Decompositions. SIAM Review, 53(2):217-288, 2011. doi: 10.1137/090718800 URL htttptt/t///to.... 1t.1.//00000011111000007773100 777000/088880//3//00000000000000
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising Diffusson Probabilistic Models, In Advances Neural Information Processing Systems, 2020.
10

---

Published as a conference paper at ICLR 2023
Richard Zou Horace He. functorch: JAX-like composable function transforms for PyTorch. http::/github..omm/pp htw.g//hhc...ooooooooo//tttttttt 2021.
Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2008.
Aapo Hyvärinen. Estimation of Non-Normalized Statistical Models by Score Matching. Journal of Machine Learning Research, 6:695-709, 2005.
Durk P Kingma and Yann Cun. Regularized estimation of image statisties by Score Matching. In J Lafferty, C Williams, J Shawe- Taylor, R Zemel, and A Culotta (cds.), Advances in Neural Information Processing Systems, volume 23. Curran Associates, Inc., 2010.
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien Gradient- Based Neural DAG Learning. In International Conference on Learning Representations, 2020.
Wai- Yin Lam, Bryan Andrews, and Joseph Ramsey. Greedy Relaxations of the Sparsest Permutation Algorithm. In The 38th Conference on Uncertainty in Artificial Intelligence, 2022.
Yingzhen Li and Richard E Turner. Gradient Estimators for Implicit Models. In International Conference on Learning Representations, 2018.
Ignavier Ng, Shengyu Zhu, Zhuangyan Fang. Haoyang Li, Zhitang Chen, and Jun Wang. Masked Cradient-Based Causal Structure Learning. In Proceedings of the SIAM International Conference on Data Mining (SDM), pp. 424-432. 2022.
Young Woong Park and Diego Klabjan. Bayesian Network Learning via Topological Order. Journal of Machine Learning Research, (899))11332 2017.
Jonas Peters and Peter Bühlmann. Structural intervention distance for evaluating causal graphs. Neural Computation, 27(3):771-799, 2015.
Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal Discovery with Continuous Additive Noise Models. Journal of Machine Learning Research, 15(58):2009-2053 2014.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference. MIT Press, 2017.
Garvesh Raskutti and Caroline Uhler. Learning directed acyclic graph models based sparsest permutations. Stat, 7(1):e183, 1 2018.
Paul Rolland, Volkan Cevher, Matthaus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score Matching Enables Causal Discovery of Nonlinear Additive Noise Models. In Proceedings of the 39th International Conference on Machine Learn ing, volume 162 of Proceedings of Machine Learning Research, pp. 18741-18753, 2022.
David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams. Learning representations by back- propagating errors. Nature, (23688))) 1986.
Karen Sachs, Omar Perez, Dana Pe'er, Douglas Lauffenburger, and Garry P Nolan. Causal Srotein-Signalle Networks Derived from Multiparameter Single-Cell Data. Science, 308(5721): 523-529, 2005.
Pedro Sanchez and Sotirios A Tsaftaris. Diffuusinn Causal Models for Counterfactual Estimation. In Proceedings of the First Conference on Causal Learning and Reasoning, volume 177 of Proceed- ings of Machine Learning Research, pp. 647-668. PMLR, 8 2022.
Pedro Sanchez, Jeremy P Voisey, Tian Xia, Hannah Watson, Alison Q O'Neil, and Sotirios A Tsaftaris. Causal machine learning for healthcare and precision medicine. Royal Society Open Science, 98). 8 2022.
Bernhard Schölkopf and Alexander Smola. Learning with kernels. Optimization, and Beyond. MIT press, 1(2), 2002.
11

---

Published as a conference paper at ICLR 2023
Jiaxin Shi, Shengyang Sun, and Jun Zhu. A Spectral Approach Gradient Estimation for Implicit Distributions. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp- 18a.
4653. PMIR, 11 2018.
ngs hww.mm.m..........m..mmmm....m rallr... waress/v0///0 htt.
html
Shohei Shimizu, Patrik O Hoyer, Aapo Hyvrinen, and Antti Kerminen. A Linear Non-Gaussian Acyclic Model for Causal Discovery. Journal of Machine Learning Research, 7(72):2003-20300 2006.
Si Si, Cho-Jui Hsieh, and Inderjit S Dhillon. Memory Efficient Kernel Approximation. Journal of Machine Learning Research, 18(20)1-32,, 2017. URL
15
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep Unsuper- vised Learning using Nonequilibrium Thermodynamics. Proc. of 32nd International Conference on Machine Learning. 3:2246 2255, 12 2015.
L Solus, Y Wang, and C Uhler. Consistency guarantees for greedy permutation-bastt causal infer- ence algorithms. Biometrika, 108(4):795-814, 8 2021.
Yang Song and Stefano Ermon. Generative Modeling by Estimating Gradients of the Data Distribu- tion. Advances in Neural Information Processing Systems, 32, 2019.
Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon. Sliced Score Matching: A Scalable Ap- proach to Density and Score Estimation. In Ryan P Adams and Vibhav Gogate (eds.), Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, volume 115 of Proceedings of Ma- chine Learning Research, pp. 574-584. PMLR, 9 2020.
Yang Song Song, Jascha Sohl- Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-Based Generative Modeling Through Stochastic Differential Equations. In ICLR, 2021.
Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search. 81, 1993.
Mare Teyssier and Daphne Koller. Ordering-Basee Search: A Simple and Effective Algorithm for Learning Bayesian Networks. In Proceedings of the Twenty-Firs, Conference on Uncertainty in Artificial Intelligence, UAI'05, pp. 584-590, Arlington, Virginia, USA, 2005. AUAI Press.
Tim Van den Bulcke, Koenraad Van Leemput, Bart Naudts, Piet van Remortel, Hongwu Ma, Alain Verschoren, Bart De Moor, and Kathleen Marchal. SynTReN: a generator of synthetic gene expression data for design and analysis of structure learning algorithms. BMC Bioinformatics, 7 (1):43, 2006.
Thomas Verma and Judea Pearl. Causal Networks: Semantics and Expressiveness. In Ross D SHACHTER, Tod S LEVITT, Laveen N KANAL, and John F LEMMER (eds.), Uncertainty in Artificial Intelligence, volume of Machine Intelligence and Pattern Recognition, pp. 69-76. North-Hollane 1990.
Pascal Vincent. A connection between score matching and denoising autoencoders. Neural Compu- 12 2011.
tation,
Xiaoqiang Wang, Yali Du, Shengyu Zhu, Liangjun Ke, Zhitang Chen, Jianye Hao, and Jun Wang- Brdering-Based Causal Discovery with Reinforcement Learning. In Zhi-Hua Zhou (ed.), Pro- ceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI, pp. International Joint Conferences Artificial Intelligence Organization, 8 2021.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN; DAG structure learning with graph neural networks. In 36th International Conference on Machine Learning, ICML 2019, volume 2019- June, pp. 12395-12409 International Machine Learning Society (IMLS), 2019.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P Xing. DAGs with NO TEARS: Con- tinuous Optimization for Structure Learning. In n Advances in Neural Information Processing Systems, 2018
12

---

Published as a conference paper at ICLR 2023
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. Learning Sparse Non- parametric DAGs. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pp. 3414- 3425. PMLR, 9 2020.
13

---

Published as a conference paper at ICLR 2023
A PROOFS
We re-write Equation here for improved readability:
Ölogp' (x; f3)
0f: Ologp' p' (x, fi)
Vx, log
(10)
Õx, (ECH(X))
8x,
Ōr
A.1 PROOF LEMMA
Proof. We start by showing the "=" direction by deriving Equation [ the first term of the equation is present, then taking derivative results
w.r.t. X5. If x, is leaf, only
H(,(log((x)) 8 log Df (xi fi).
dfi
Ox
dx:
(11)
8 log p' (xl fi)
Ox
Therefore, only j is a leaf and d'logp' = a, Var (log p(x)) ((x)) 0. The remaining of the proof dy follows Rolland et al. (2022) (which was done for Gaussian noise only), we prove by contradiction that + is also true. In particular, if we consider that is not a leaf and Hjj log (og)(x) = C, with e being constant,
write
x, log (() = cx₃ + g(x-)).
(12)
Replacing Equation 12 in to Equation
we have
ex + -(x_)) 2 Olog p' (x, - fj)
0f: Olog p' (x, - fi)
(13)
dr
0x,
E.C(()))
Let X. E Ch((;) such that X, & Pa(Ch(x,)), X, always exist since xj is not a leaf, and it suffices pick child of x, appearing last position in some topological order. If we isolate the terms depending on X on the RHS of Equation 13
we have
ex, + 2 log pd (x; fj)
Of log (xi fi)
Ofe Olog (Xe (xe-fe) fo) Ox,
Σ
g(x_x)
dx,
Ōx
dx
(14)
Deriving both sides w.r.t. Xc, since the LHS of Equation
14 does not depend on Xc, we can write
0fe d Ologrf (xe Sc)
9(x-;)
Ox.
Öx,
dx
(15)
Se θ llogr
(N_- 1-)= 0g(x
0x
Ox²
Ox.
0h. does not depend on xj neither, implying that fo is linear in Xj. Since g does not depend Xj, 8xy contradicting the non-linearity assumption.
A.2 PROOF THEOREM
Proof. Using Equation I
we will derive expressions for each of the elements in Equation 6 and
show that it is equivalent to in Equation 5 First, note that the score of a leaf node x can be denoted as:
2 log (x, - fz)
Vx, log p(x)
0x,
i log pf (xi ff) d (xi fi)
(16)
dr.
dxi
a logpp (xi fi)
dr
14

---

Published as a conference paper at ICLR 2023
Second, replacing Equation into each element of (logg((x)) E R write
H1,(logp(x))
xx logp(x)
0xj
8² log p$ (xt - ft)
(x1 - fi)
(17)
dx²
dx3
0 log pf (xt fi). dft
Ox²
dxj
If = in Equation 17
we have
82 logp' of (xi -fi)
dfi
(log(o)((x))
Ox2
dx:
(18)
8 log pt (xi
fi)
dx
Finally, replacing Equations
16
and into the Equation for single node xy, j * I, we have
6, H (l,s(og p(x)) V log (xx)
Hi, (log p(x))
log (x
ft)
dft
Blog p' f1)
Ox
dxy
dR
(19)
Bog D' (x1- /t)
Ou
dfi & log pe (xi ft)
dx,
dr
The last line in Equation 19 is the same as Equation from Lemma proving that 69 can be written using the first and second order derivative of the likelikel.
D
B SCORE OF NONLINEAR ANM WITH GAUSSIAN NOISE
SCM entails a distribution
p(x) II p(x: Pa(x:)),
(20)
over the variables (Peters et al. 2017 By assuming that the noise variables ~ (0, and
inserting the ANM function, Equation 20 can be written
1
X; [((Pa(x)))
2
logp(x) -
log (2(o7)))
(21)
2
The score of p(x) can hence be written as
Vx, logp(x) Xj [s(Pa(x,))
05, (Pa(x;)) [((Pa(x:))
Σ
(22)
of
:Echildren(s)
C VISUALISATION OF THE SCORE'S JACOBIAN FOR TWO VARIABLES
Considering two variables problem where the causal mechanisms are B = (A) and = with EA, EB NO 1) and f., being a two-layer MLP with randomly initialised weights. Note that, in Figure 5 the variance of 8² log ((A.)) while not 0 predicted by Equation is smaller than 8 log p(A.B) allowing discovery of the causal direction.
15

---

Published as a conference paper at ICLR 2023
@ log pl B) 8A2
8f log, pi A. B) 3
A
3.0
1.5
0.0
B
-1.5
-3.0
-3.0
-1.5
0.0
1.5
3.0
Var: 0.12
Var: 0.07
Figure 5: Visualisation of the diagonal of the score's Jacobian estimated with diffusion model as in Equation 7 for two-variable SCM where A + B.
Table 1: MLP architecture. The hyperparameters of each Linear layer depend on such that big max(1024, 5 d) and small max 128, 3 d).
maz 1024,
Layer
Hyperparameters Ch: (d small)
Linear
LeakyReLU
LayerNorm
Dropout
Prob: 0.2
Linear
Ch: (small, big)
LeakyRel.U
LayerNorm
Linear
Ch: (big, big)
LeakyReLU
Linear
Ch: (big, big)
LeakyReLU
Linear
Ch: (big, d)
D EXPERIMENTS DETAILS
D.1 HYPERPARAMETERS OF DPM TRAINING
We describe the hyperparameters for the diffusion training. We number of time steps T= 100, 3e a linearly scheduled between 8min 0.0001 and 3max = 0.02. The model trained
according Equation which follows Ho al. 2020). During sampling, is sampled from Uniform distribution.
D.2 NEURAL ARCHITECTURE
The neural network follows simple MIP with 5 Linear layers, LeakyReI.U activation function, Layer Normalization and Dropout the first layer. The full architecture is detailed in Table T
D.3
METRICS
For each method, we compute the
SHD. Structural Hamming distance between the output and the true causal graph, which counts the number missing, falsely detected, or reversed edges.
SID. Structural Intervention Distance is based graphical criterion only and quantifies the closeness between two DAGs in terms of their corresponding causal inference statements PPtters)
& Bühlmann 2015).
Order Divergence. Rolland et al. (2022) propose this quantity for measuring how well the topolog- ical order is estimated, For ordering म, and target adjacency matrix A, we define the topological
16

---

Published as a conference paper at ICLR 2023
order divergence Dtop (T, A)
d Σ (=1 1:8,33
Dtop (TT (T, A)
Σ Aij.
(23)
If is correct topological order for A, then (x, 0. Otherwise, Drop( A) counts the number edges that cannot be recovered due to the choice topological order. Therefore, provides a lower bound on the SHD of the final algorithm (irrespective of the pruning method).
E OTHER RESULTS
E.1
REAL DATA
We consider real datasets: (i) Sachs: A protein signaling network based expression levels of proteins and phospholipids (Sachs et al. 2005). We consider only the observational data (n 853 samples) since our method targets discovery of causal mechanisms when only observational data is available. The ground truth causal graph given by [Sachs et al. (2005) has 11 nodes and 17 edges. (ii) RynTReN: We also evaluate the models on a pseudo-real dataset sampled from SynTReN gener- ator (Van den Bulcke et al. 2006). Results, in Table 2 show that our method is competitive against other state-of-the-ee causal discovery baselines on real datasets.
Table 2: SHD and SID results over real datasets.
Sachs
SynTReN
SHD
SID
SHD
SID
CAM
12
55
40.5
152.3
GraN-DAG
13
47
34.0
161.7
SCORE
12
45
36.2
193.4
DiffAN (ours)
13
56
39.7
173.5
E.2 OVERFITTING
The data used for topological ordering (inference) is subset of the training data. Therefore, it obvious if overfitting would be an issue with our algorithm. Therefore, we run an experiment where fix the number of epochs another set in order to verify if overfitting is an issue. On average across all 20 nodes datasets, the early stopping strategy output an ordering diverge of 9.5 whilst overfitting at 11.1 showing that the method does not benefit from overfitting.
2000 considered high for set of and use early stopping for
E.3 OPTIMAL / FOR SCORE ESTIMATION
As noted by Vincent 2011]) b, the best approximation of the score by learned denoising function is when the training lignal-tonoiis (SNR) ratio is low. In diffusion model training, = O corresponds the coefficient with lowest SNR. However, we found empirically that the best score estimate varies somehow randomly across different values of t. Therefore, we run the the leaf finding function (Equation 9 N times for values of evenly spaced in the [o, 2]] interval and choose the best leaf' based on majority vote. We show in Figure that majority voting is better approach than choosing a constant value for
E.4 ABLATIONS OF DIFFAN MASKING AND GREEDY
We now verify how DiffAN masking and DiffAN greedy compare against the original version de- tailed the main text which computes the deciduous score. Here, we use the same datasets decribed in Section 5. which comprises 4 (20ER1, 20ER5, 20SF1, 20SF5) synthetic dataset types with 27 variations over seeds, noise type and noise scale.
DiffAN Greedy. A greedy version of the algorithm re-trains the eø after each leaf removal iteration. In this case, the deciduous score is not computed, decreasing the complexity w.r.t. but increasing w.r.t. n. DiffAN greedy has O(nd²) ordering complexity.
17

---

Published as a conference paper at ICLR 2023
10ER5
majority voting
10SF5
20ER5
20SF5
50ER5
50SF5
10
20
30
40
50
Topological Order Divergence
Figure 6: The distribution of order divergence measured for different values of is highly variable. Therefore, we show that we obtain better approximation with majority voting.
We observe, in Figure that the greedy version performs the best but it is the slowest, as seen in Figure DiffAN masking
400
method
300
DifFAN
DifIANM BilfANMmmakkig
DifTAN DiffAN_greedy
SID
200
100
20ER1
20SF1
20ER5
20SF5
Experiment
Figure 7: SID metric for different versions DiffAN.
E.5 DETAILED RESULTS
We present the numerical results for the violinplots in Section 5.1 in Tables 3 and The results are presented in mean with statistics acquired over experiments with 3 seeds.
18

---

Published as a conference paper at ICLR 2023
200
150
runtime
100
50
20ER1
20SF1
20ER5
20SF5
Experiment
Figure 8: Runtime in seconds for different versions of DiffAN.
Table 3: Erdös-Rényi (ER) graphs
Exp Name
noisetype
method
shd
sid
runtime (seconds)
20ERI
exp
CAM
6.503.21
35.5027.54
112.17
DiffAN
7.67
38.8321.94
30.071.00
DiffAN Greedy
6.67 1.77.68
34.1718.37
121.661.19
GranDAG
12.835.60
65. 65.1721.99
341.6974.82
SCORE
4.172.86
17.5011.26
50.4941.87
gauss
CAM
1.331.03
6.506.69
83. 83.2711.75
DifFAN
8.174.02
41.6724.21
30.982.28 28
DifFAN Greedy
3.332.16
18.3311.40
129. 229.10g.64
GranDAG
10.50
44.5037.90
304.0717.39
SCORE
8.67
41.6721.54
24.77a.20
laplace
CAM
0.83g.08
4.006.96
90.8515.49
DiffAN
2.332.07
8.674.68
31.01
DiffAN Greedy
3.001.67
14.334.76
175.3019.71
GranDAG
10.33
41.1724.98
302.0016.23
SCORE
4.173.43
24.17
26.62
20ER5
exp
CAM
60.679.7T
240.0033.02
105.6527.49
DiffAN
67.50a.32
278.8322.88
30.332.12
DiffAN Greedy
63.006.07
266.3350.25
186.2 286.213.28
GranDAG
80.177.99
304.1726.90
310.241a.71
SCORE
55.837.35
190.6732.91
18.41
gauss
CAM
64.675.85
214.6711.11
130.8226.11
DiffAN
68.50
289.0042.31
29. 29.581.99
DiffAN Greedy
63.176.18
228.83
178. 178.177.17
GranDAG
81.006.72
273.0046.72
319.1718.81
SCORE
63.50
233.8339.0
21 21.296.63
laplace
CAM
68.007.04
228. 228.5027.80
157. 7022.502..5.
DiffAN
68.837.86
248.67 30.36
29.521.73
DiffAN Greedy
67.33g.55
238.00g5.19
176. 176.71480
GranDAG
82.67
271.3335.19
305.20 205.2913.93
SCORE
66.337.09
218.6725.81
23.01 4.33
19

---

Published as a conference paper at ICLR 2023
Table 4: Scale Free (SF) graphs.
Exp Name
noisetype
method
shd
sid
runtime (seconds)
20SFI
exp
CAM
7.172.88
32.67
85.0827.71
DiffAN
8.672.80
28.6720.99
29. 10.2311mm 36
DiffAN Greedy
9.172.32
35.3317.20
183,80 183.801.s1
GranDAG
15.502.59
68.6724.19
998.884.96
SCORE
5.835.04
23.5023.65
26. 26.131.39
gauss
CAM
1.832.23
9.17
84.4619.07
DiffAN
9.836.24
56.6741.50
30.19h.9
DiffAN Greedy
7.335.32
46.8338.42
175. 175.63 63,.8288
GranDAG
12.172.71
46.6716.22
299. 299.2511.14
SCORE
8.174.15
41.6728.75
26. 26.291885
laplace
CAM
2.833.92
1.83
1521.84
DifFAN
6.003.74
24.1716.18
29. 29.682.00
DiffAN Greedy
4.833.66
17.1717.47
177.058.82
GranDAG
9.503.27
35.6720.84
295.122.88
SCORE
5.673.50
22.5015.27
26.391.79
20SF5
exp
CAM
47.83s.78
228. 128.8363.88
113.8213.92
DiffAN
11.83
243.5024.68
30.032.17
DiffAN Greedy
43.506.80
236.8323.17
173.133.52
GranDAG
60.67
275.0022.56
284.345.33
SCORE
38.50
180.3357.44
19. 19.792.92
gauss
CAM
46.83
199. 199.1753.13
130.4120.55
DiffAN
50.834.62
259.5047.54
28.581.16
DiffAN Greedy
45.176.15
224.1741.80
174.903.69
GranDAG
61.674.08
241.6743.11
292.777.29
SCORE
44.505.24
217.8350.05
19. 19.772.51
laplace
CAM
49.50g.01
191.3327.43
160.5323.16
DiffAN
52.005.66
230.6750.66
29.08 29.051.6
DiffAN Greedy
47.00g.91
191.3326.16
173.423.52
GranDAG
65.007.18
262.5044.25
295.407.72
SCORE
46.6710.03
193.8334.14
43.8427.14
20


Reviews:
[{"review_id": "g5jc9M9zFx", "reviewer": "Reviewer_1RL4", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper considers non-linear causal discovery under the same setting as Rolland et al., (2022). The authors proposed to replace the scoring function in the SCORE method proposed by Rolland et al., (2022) with diffusion probabilistic models (DPMs). This allows them to update the learned Hessian without re-training the neural network and perform ordering over a batch, and thus they can large-scale data at a relatively faster speed. However, in total, the novelty of the approach and the profound impact are both limited. In addition, the main methodology and numerical results require additional justification.", "summary_of_the_review": "Though the paper considers an important problem in causal discovery for the non-linear additive noise model, the novelty of the approach and the profound impact are both limited. In addition, the main methodology and numerical results require additional investigation.", "strength_and_weaknesses": "*Strength*\n\n- The paper considers an important problem in causal discovery and leveraged a popular machine-learning tool.\n\n*Weaknesses*\n\n- The proposed DiffAN method highly depends on the SCORE method proposed by Rolland et al., (2022), by replacing the scoring part with diffusion probabilistic models (DPMs). More justification on novelty is needed. For instance, can the authors comment on if their method can be extended to other causal discoveries for non-linear models? In addition, please comment on why we cannot decouple DPMs from the SCORE method and consider a new causal discovery method built upon DPMs. I am also curious whether DPMs are really necessary for this approach. Why not consider other deep learning approaches since we are not handling images/texts?\n\n- Although firstly introducing DPMs into causal discovery in this paper, the authors did not well explain why simply dropping the leaf in DPMs is sufficient for updating the score function without retraining. For instance, according to Algorithm 1, the scoring step is batch-dependent. If I understand correctly, when we update the score by pruning DPMs, it is not matching to the distribution of the next batch of data. In particular, the data $x$ is changing over batches in Equations 4 and 8. More theoretical justification is needed besides Figure 4.\n\n- The numerical studies did not show the best performance of the proposed method. Indeed, Figure 3 suggests that both CAM and SCORE are better than the proposed DiffAN method. In addition, the runtime of SCORE in Figure 3 is still reasonable. Similarly, both GraN-DAG and SCORE perform better than the DiffAN method in real data  I would suggest running more studies that vote for the authors' method. Finally,\nwhy not consider batches in SCORE so we can improve the run time as well? \n \n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "Mwq2w8lgIxF", "reviewer": "Reviewer_YaE5", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The structure of the Hessian of the data log likelihood carries information about topological ordering.  Since this information is learned by score based generative models, DPMs carry within them information necessary to reconstruct topological orderings.  The paper runs with this insight to obtain a practical algorithm with interesting scaling properties for causal discovery.", "summary_of_the_review": "Please provide a short summary justifying your recommendation of the paper.\nThe paper rigorously explores a topic of interest to the community but does not comment on what would seem to me to be a more direct approach to solving the same problem (i.e. using a KDE).  The authors should comment on such alternatives to computing / approximating these Hessians that do not rely on diffusion models.\n\nIf the authors do not provide an explanation for why other approaches to computing the Hessian are ill-suited, or discuss how they might be applied here instead, I will significantly lower my score (to a reject recommendation).\n", "strength_and_weaknesses": "Strenghts:\n* The paper presents a timely analysis of how causal relationships may be uncovered using diffusion models.  Given that both causality and diffusion models are of interest to many in the community, the possible connections are certainly on the minds of many (including my own!).  \n*The paper runs with this direction and does a very nice job rigorously establishing some of what is possible with the Hessians of the log densities learned within diffusion probabilistic models.\n*The paper is clearly written and the was enjoyable to read.\n\nWeaknesses:\n* Why not do the same analysis on a Gaussian kernel density estimate (KDE) of the data distribution? In my opinion, the paper is unpublishable without explicitly acknowledging this alternative route to the same ends, and how it would compare to the proposed method involving diffusion models.  In particular, in the limit of a large enough network / good enough optimization the Hessians computed coincide with what would obtain with the KDE.  So why go through the trouble of the diffusion model?  And is there anything special about diffusion models in this paper other than the fact that they provide computationally tractable Hessians (which is also possible with energy based and likelihood based models)?\n* There may be good reasons to use a neural network here (e.g. amortization of Hessian computation, or the generalization properties may lead to better estimates of the Hessians of the population distribution as compared to the KDE). But these possibilities should be made explicit.  Without acknowledging that the approach involves what may be viewed as an approximation of something that may in principle be computed analytically, I view the paper as disingenuous.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "LC89-Jyam6", "reviewer": "Reviewer_vVJW", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes a more efficient approach for finding a topological ordering based on identifying low variance entries along the diagonal of the Hessian of the log density. Unlike the previous approach, this method does not require iteratively reestimating the entire Hessian to identify each leaf node. The approach is implemented using diffusion models, since they approximate the score function. The resulting method performs comparably to the previous state of the art but is more efficient in high dimensions.", "summary_of_the_review": "The paper proposes an improvement to causal discovery based on estimating the variance of the Hessian log density diagonal that generalizes the causal discovery method and improves efficiency. The paper makes two core technical contributions and empirical results demonstrate the gain in efficiency. The clarity of the paper could be improved to emphasize the specific contributions made and the role (and necessity) of diffusion models in the proposed approach.", "strength_and_weaknesses": "Strengths:\n- New theoretical insights are provided, which generalize the results from Rolland et al. (2022) \n- A more efficient algorithm is provided \n- Empirical results demonstrate the gain in efficiency \n\nWeaknesses:\n- Too much emphasis on the diffusion models aspect makes the contribution confusing (see below)\n- More explicit connections between the technical results and the correctness of the algorithm would improve clarity ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "-Kmsaprj-K", "reviewer": "Reviewer_F1bu", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors proposed DiffAN, a topological ordering algorithm that leverages DPMs for computing the Hessian. The algorithm updates the learned Hessian without re-training the neural network and performs ordering over a batch, which allows scaling to datasets with more samples and variables. In the experiments, they show that the method scales up well and achieves SOTA results. ", "summary_of_the_review": "The proposed method is novel and interesting.  The author should discuss the reason why the method cannot outperform the score-based method. \n\n\n", "strength_and_weaknesses": "Strength:\nIntegrating causal discovery and diffusion models is novel and interesting. The proposed method using diffusion methods to compute the Hessian of distribution scores is interesting.\n\n\nWeakness:\nThe method only achieves comparable results on the data sets. The author should give detailed reasons why the method could not beat existing methods given that deep neural nets can provide better representations. ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster