Structure Learning with Adaptive Random Neighborhood Informed MCMC
Authors: Xitong Liang
Alberto Caron
Samuel Livingstone
Jim Griffin
Venue: NeurIPS 2023 poster (2023)
Keywords: Bayesian Networks
structure MCMC on graphs
Structure Learning
Random neighborhood samplers
Locally informed Metropolis-Hastings schemes

Abstract: In this paper, we introduce a novel MCMC sampler, PARNI-DAG, for a fully-Bayesian approach to the problem of structure learning under observational data. Under the assumption of causal sufficiency, the algorithm allows for approximate sampling directly from the posterior distribution on Directed Acyclic Graphs (DAGs). PARNI-DAG performs efficient sampling of DAGs via locally informed, adaptive random neighborhood proposal that results in better mixing properties. In addition, to ensure better scalability with the number of nodes, we couple PARNI-DAG with a pre-tuning procedure of the sampler's parameters that exploits a skeleton graph derived through some constraint-based or scoring-based algorithms. Thanks to these novel features, PARNI-DAG quickly converges to high-probability regions and is less likely to get stuck in local modes in the presence of high correlation between nodes in high-dimensional settings. After introducing the technical novelties in PARNI-DAG, we empirically demonstrate its mixing efficiency and accuracy in learning DAG structures on a variety of experiments.

Structure Learning with Adaptive Random Neighborhood Informed MCMC
Alberto Caron'
Xitong Liang"
The Alan Turing Institute London, UK acaronSturing ac. uk
Department of Statistical Sciences
University College London
London, UK xitong liang. 18@ucl. ac.
Samuel Livingstone Department of Statistical Sciences
Jim Griffin
Department of Statistical Sciences University College London
University College London London, UK samuel samuel. hiv.ngl..iiinggttte@uuul... ac. uk
London, UK
j griffin@ucl. ac uk
Abstract
In this paper, we introduce novel MCMC sampler, PARNI-DAG, for fully- Bayesian approach to the problem of structure learning under observational data. Under the assumption of causal sufficiency, the algorithm allows for approxi- mate sampling directly from the posterior distribution on Directed Acyclic Graphs (DAGs). PARNI-DAG performs efficient sampling DAGs via locally informed, adaptive random neighborhood proposal that results in better mixing properties. In addition, to ensure better scalability with the number of nodes, we couple PARNI- DAG with a pre-tuning procedure of the sampler's parameters that exploits a skeleton graph derived through some constraint-based Thanks to these novel features, PARNI-DAG quickly converges to high-probabilin regions and is less likely to get stuck in local modes in the presence of high corre- lation between nodes in high-dimensional settings. After introducing the technical novelties in PARNI-DAG, we empirically demonstrate its mixing efficiency and accuracy in learning DAG structures variety of experiments.
scoring-based algorithms.
Introduction
Structure Learning and Causal Discovery are concerned with reconstructing a graphical model underlying set of random variables from data, in the form of Directed Acyclic Graph (DAG), provided that causal identifiability assumptions hold [Pearl, 2009, Drton and Maathuis, 2017, Glymour al., 2019]. Causal relations can typically inferred if randomized experiments are carried out [Pearl, 2009]. However, in many relevant applied fields (e.g. biology. genomics, ecology, etc.), only observational data are accessible, as performing interventions usually costly or simply unfeasible. Structure learning from observational data is challenging problem, with the space of possible DAGs growing super-exponentially the number of nodes [Robinson, 1977].
Related work. Early contributions on structure learning algorithms include constraint-based methods [Spirtes al., 2000, Pearl, 2009], such as PC and FCI, that utilize conditional independence
"These authors contributed equally to this work.
'Code to implement the PARNI-DAG proposal and replicate experimental sections available attps://gitaaie com/ ctw./iitooniiinggeeeeeiPllaa wwtph/heeePAAR-sheeeeee wetaameeeeeeeeeeeeeeee SamiiinSttucture-eeaaaee
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

testing to output a Markov Equivalence Class (MEC). Another stream of literature focuses on score- based methods [Geiger and Heckerman, 1994], such as Greedy Equivalent Search (GES) [Chickering. 1996, 2002], which define and consequently maximize score function associated with different MECs. Alternatively, functional causal models, such as LiNGAM [Shimizu et al., 2006] and Additive frrrncttoely. Noise Models [Hoyer et al., 2008bJ, directly output DAGis Recent advancements falling in either of the classes above have considered methods relying on deep learning techniques [Zheng et al., 2018, Yu et al., 2019], but are specifically suited for large sample regimes.
MCMC-based Bayesian approaches, that account for graph uncertainty by learning a posterior distri- bution over graphs, include Madigan et al. [1995] and Giudici and Castelo [2003], who introduced the MCMC Model Composition (MC*) scheme and the Reveeee (ADR) scheme respectively. Grzegorzzyk and Husmeier [2008] improve ADR's 's mixing by proposing new arc reversal move. Another stream of contributions on structure learning MCMC have considered working with score functions in smaller spaces rather than the DAG one. Friedman and Koller [2003] first recasted the problem by operating in the space of "orders" instead of DAGs, so that candidate orders, consistent with multiple DAGs, are proposed. Although this generates considerable benefits in terms of scalabil- ity, it introduces sampling bias, as some DAGs are over-represented [Ellis and Wong, 2008]. This bias addressed in Kuipers and Moffa [2017], who proposed working the space of "ordered partitions" instead, which introduces scalability issues, making partition MCMC applicable to graphs with a small number of nodes. Further works have focused on improving mixing efficiency and scalability of these score-based MCMC [Ninimāki et al., 2016, Viinikka and Koivisto, 2020, Kuipers et al., 2022]. Recent advances have also considered approximate Variational Inference methods [Lorch et al., 2021, Cundy et al., 2021, Geffner et al., 2022]. More closely related to our work, van den Boom et al. [2022] considered the Docally-ifformed proposal on non-decomposable Gaussian Graphical Models (GGM), but with neighborhoods defined by subset of the space around the current state of the Markov chain, i.e. neighborhoods of the MC type, with size growing quadratically in the number of nodes. However, their WWA algorithm does not target the posterior on DAGs directly, unlike PARNI-DAG.
Tangential to this work is also the rich literature on discrete spaces MCMC samplers for Bayesian variable selection. Noteworthy contributions include Zanella [2020], that up-weights a random- walk kernel via local neighborhood information and constructs the locally- -nneoll-i proposal, and Griffin et al. [2021], that develops an adaptive MCMC algorithm addressing high-dimensional sparse settings. The locally-informed proposal can be viewed as a discrete analog to the gradient MCMC Langevin algorithm [Roberts and Rosenthal, 19981), where the gradient of the target distribution does not exist. The adaptive MCMC instead is capable of learning from the shape of the target distribution and exploiting the most important variables. Liang et al. [2022] generalise the locally-informed proposal and the adaptive MCMC to be used within random neighborhood schemes, and introduced the PARNI proposal for Bayesian variable selection.
proposals (e.g. the Metropolis-adjuse
Contribution. This work presents a novel MCMC sampler for structure learning, PARNI-DAG, that builds on top of PARNI [Liang et al.. 2022], but is equipped with new features for efficient sampling in the space of DAGs, PARNI-DAS constructs a random neighborhood of possible DAGs via probabilities proportional to function of the Posterior Edge Probabilities (PEPs), that is, the probabilities of including a direct edge between two nodes. Then, since full enumeration of all the possible DAGs within the neighborhood is virtually impossible in high-dimensions, it proposes new candidate DAG via a "point-wise update", i.e. by considering a sequence of intermediate proposed DAGs belonging to smaller subset of the neighborhood. Finally, to guarantee better scalability with the number of nodes, we propose procedure to pre-tune PARNI-DAG's parameters and warm-start the chain by utilising a skeleton graph derived through any constraint score based algorithm (e.g. PC, GES,... or the iterative search space expansion proposed in Kuipers et al. [2022].
Motivations. As we will discuss in later sections, PARNI-DAG generates clear advantages over structure MCMC methods such as MC and ADR in terms of speed of convergence and mixing, thanks to its new improved proposal, while targeting the same posterior distribution over DAGs. As for the comparison to score-based MCMC methods, it is known that order MCMC [Friedman and Koller, 2003] has better scalability than standard structure MCMC methods, but, operating in the smaller space of order, it introduces sampling bias via a non-uniform prior being assigned to the different DAGs compatible with single order [Ellis and Wong, 2008]. Compared to order MCMC, PARNI-DAC does not incur in any sampling bias as it targets the posterior on DAGs directly, while
2

---

at the same time converging faster than ADR. Kuipers and Moffa [2017] propose a variant called partition MCMC, which operates in the space of ordered partitions instead of orders. Partition MCMC is unbiased in terms of sampling, but it is extremely slow and has high computational complexity (can only be used on very few nodes) [Kuipers et al., 2022]. PARNI-DAG, on the other hand, specifically addresses high-dimensional settings with many nodes.
As solution to scale order and partition MCMC up to high-dimensional scenarios, the hybrid MCMC scheme of Kuipers et al. [2022] proposes to restrict the initial search space by pruning it via skeleton graph (e.g. PC algorithm derived skeleton) H. whose max parent set size is m per node, so that complexity reduces to ((n2)).. Although this increases MCMC speed dramatically in large graphs, it is likely to introduce bias as some true edges might be deleted in H. To tackle this issue they also propose an iterative procedure (Iterative MCMC) to pe-populate H with additional potential parent nodes, at the expenses of an increased computational time. We couple PARNI-DAG with a similar procedure, but we use a previously derived skeleton H (e.g. the Iterative MCMC one) not to restrict the space, and possibly cancel out relevant edges, but to warm-start the chain by pre-tuning some of the MCMC, as we will explain in details in later sections.
We show in the experiments section how PARNI-DAG brings about improvements over the classical structure MCMC ADR [Giudici and Castelo, 2003, Grzegorczyk and Husmeier, 2008] and score- based MCMC methods [Friedman and Koller, 2003, Kuipers and Moffa, 2017, Kuipers et al., 2022] in terms of DAG learning accuracy and MCMC mixing, particularly in settings with high number of nodes, as it able to reach high probability regions in very few steps, while addressing the presence of highly correlated neighborhoods of edges where classical samplers might get trapped. Code to implement PARNI-DAG and fully reproduce the experiments is provided'.
2 Problem Setup
Consider a graph 9 = (V, E), made of nodes V and edges E, suitable to represent probabilistic depen- dencies (edges) between random variables (nodes) through probabilistic graphical model [Koller and Friedman, 2009]. Given a collection of n = [V] continuous random variables (X1,.., X..). a Bayesian Network B (G, 2) associated with a DAG G. is probabilistic graphical model utilized to represent a factorization of the joint probability p(.I...., ) into the conditional distributions p(X,,..., X,, IT= P (X, , (P((X))). where ø E & are the parameters and Pa X,) all the parents of node X,- Suppose that all nodes (X1,-.. X,,) are re-scaled to have zero mean and unit variance. The goal in structure lcarning is to learn the BN factorization of conditional probabilities, given sample of N observations on D = {X1,4, ..., xx.) J. In Bayesian structure learning, we specif- ically want to Iearn posterior distribution on possible DAGs, i.e. p(G[D) x p(D(G))(G) In this work, we define DAG indicator variable 7 E r {0, 1) "=x", where each 7 implies a unique graph gy: Mr = 1 indicates the presence of an edge from X; to Xj. while Yij = 0 indicates absence of it. The linear functional model associated with DAG G, can be specified as
x= W and W,y is a matrix of weights. Within a fully Bayesian approach the inference problem presented in (1), we consider the following prior specification
where B(E) = 0, Var(E) = diag(o?....,, ......,,))..
(I)
(W>)ol of. Ki ~ Normal go
(W+) of Tk = ~ do. p(az) o0²²
where do is the Dirac measure at O. Thanks to their conjugate form, the coefficients (Wa) Jijj and the heteroskedastic error variances {a can be integrated out analytically, Ieading to a marginal likelihood that depends only on y, denoted by p(D)7). The last building block we need for posterior (?>) = P() inference is acyclicity constraints imposed in the prior distribution for the model indicator y, which is specified by
h
p(7) <
(T-T.
x I[G is a DAG]
where the hyperparameter h (0. represents the prior edge probability and d. |E,1 is the number of edges implied by 7. We are interested in Bayesian inference on 7 with posterior distribution (7) x p(D)))((7).
"Code is provided in the supplementary materials. Github link to be added upon acceptance.
3

---

Notice that the likelihood, and thus posterior distribution, could also be replaced by scoring function typically used in score- cosre-based MCMC methods [Friedman and Koller, 2003, Kuipers and Moffa, 2017, Kuipers et al., 2022], such as the BGe score [Geiger and Heckerman, 2002], or the BDe score, which would make PARNI-DAG suitable also for discrete-valued nodes.
3 The novel PARNI-DAG proposal
3.1 Point-wise implementation of Adaptive Random neighborhood Informed proposal
Recently, Liang et al. [2022] introduced the Point-wise Adaptive Random Neighborhood Informed proposal (PARNI) in the context of Bayesian variable selection problems. PARNI is characterized by MMetropolis-Hssting (MH) proposal [Metropolis et al., 1953, Hastings, 1970] which couples the Adaptively scaled Individual Adaption scheme [Griffin et al., 2021] with the locally-informed proposal [Zanella, 2020J. We begin by firstly describing the generic PARNI proposal for DAG learning problems, then we will introduce the modifications needed to make it efficient and that result in the PARNI-DAG proposal.
The PARNI proposal falls into the class of random neighborhood informed proposals, which are characterized by two steps: i) random sampling of neighbothood N, and then ii) proposal of new DAG within this neighborhood W according to informed proposal [Zanella, 2020]. The random neighborhoods are drawn based on an auxiliary neighborhood indicator variable k E K. with conditional distribution pkl]r), such that the neighborhood is constructed a function of ((; k)) S r. Let K. r = {0. 13axn, then in PARNI each k- E value indicates whether the corresponding position of is included in the neighborhood. The conditional distribution of takes a product form (k)(k) = 1k. Pn (kkz)11s)) characterized by the set of tuning parameters
= [ where Mkj E (e, <) for a small e E (0, 1/2), and each pO (kis (Ks)) is given by
Pty (kig = 1/y = 0) min {
1, 1 - Mij
(ky, = 1//s = 1) = min {
1, 1 nkj
7i
(3)
Tkj
The methods for adapting 77 will be discussed in Section 3.2 specitically. The neighborhood is thus the constructed as
((, (-" ET 1 - The vii.j) s.t. ki- 0]
The neighborhood (n,, R)) contains 2d models where de denotes the number of is in The full enumeration over the whole neighborhood is computationally expensive if dN is big, and infeasible when de is beyond 30. Liang et al. [2022] considered a point-wise implementation of the algorithm, which dramatically reduces the computational complexity from ((2d)) to ((2d)). The idea of how the point-wise implementation works is the following. Let variable K = { represent the collection of positions where ky 1. Instead of working with the full neighborhood (1 k). construct sequence of smaller neighborhoods (N(()r), (,))))))) K/ C _N(r, k), and new DAG 'h' drawn from these neighborhoods {W(r)r), K== is sequentially proposed according to the sub-proposals Kg (7(r- - 1), at each time Given the intermediate DAGs 7 70) + (1)
+ (R) = , each of their corresponding sub-proposal is defined by
Jg.k. (r(r 1),)(r)) 9 (K.
1))) I(((r) N 1), K,)}
(5)
Z,-
where g [0, 00) + [0, oo) is continuous function and 2, is the normalising constant. The construction of and of the sub-neighborhoods (N(x(r), K.) will be discussed later in Section 3.3.
The choice of function g(-) is crucial to the performance of the PARNI proposal. If g(x) the sub-proposal K. in (5) reverts back to a random walk proposal which proposes a new DAG from N(>(r- 1), K,) with uniform probability. would then be hard for it to explore important neighborhoods of edges in high-dimensions as it is likely to get stuck in neighborhoods with highly correlated nodes. In the locally-informed proposal of Zanella [2020], function g(-) chosen to be a non-decreasing and continuous balancing function g [0, oo) + [0, 20). which satisfies g(x) = (1/)) Zanella [2020] shows that the locally-informed proposal with balancing function g(-) is asymptotically optimal compared to other choices of g(-) in terms of Peskun ordering. In this

---

work, we consider the Hastings' choice (i.e. g(x) = mi,, satisfying g(z) = 2g(1/ar) as it has been shown be the most stable choice for most problems (see the discussion in Supplement B.1.3 of Zanella [2020]).
The informed proposal in this context not only helps speed up convergence, but it also avoids proposing invalid DAGs. This is because the acyclicity checks for candidate models 7' E WN(1rr= 1), K2) can be run before computing their posterior model probabilities, and thus one can assign zero mass to their prior (resulting in zero mass posterior) if y is not a valid DAG. The resulting proposal kernel is then given by
n II k. (1)
ag.k (7.x())))
(6)
In order to construct a m-reversible chain and calculate the MH acceptance probability, we define the new collection of variables K" {K2 1. for the reverse moves, where K" contains the same element as K but with reverse order (i.e. K² = KR-r+1). The MH acceptance probability is then defined as
a(T, min {1 =(1'),
x()))),(()))))
Proposition 3 of Liang et al. [2022] shows that using the above K' can simplify the calculation of the MH acceptance probability in (7).
In the following sections, will introduce the necessary modifications and developments to effi- ciently sample DAGs, the combination of which defines the novel PARNI-DAG proposal. The full details and algorithmic pseudo code of PARNI-DAG are also provided in the supplementary material (Appendix D).
3.2 Warm-start of r) combined with hybrid scheme
Griffin et al. [2021] and Liang et al. [2022] both suggested using the Posterior Inclusion Probabilities in Bayesian variable selection context to define the Mιÿ hyperparameters described above. Perfectly analogous, we will make use of the Posterior Edge Probabilities (PEPs) instead, i.e. π(k; 1). The distribution in (3) with the choice of parameters described, encourages to 'flip' edge values based on variables' marginal importance. For example, suppose the posterior for 7ij is < = 1) = 0.9, and the chain is currently at 7ij = The edge should be then assigned higher probability to be included in the DAG (i.e. Jay should be switched from 0 to 1). Thus, we propose k; 1 with probability min 0. = 1, so that the new neighborhood features this possible move and hy can be effectively switched from 0 to in the next moves. Another advantage of using the conditional distribution featuring the PEPs in (3) is that when all variables are independent this proposal will have acceptance probability 1 and achieve optimal asymptotic variance as showed in Griffin et al. [2021].
The PEPs, <(%y 1). however cannot be directly calculated for most problems. Griffin et al. [2021] considered the Rao- coo-BBackwellieed estimates ((%y 1/y- where Y denote excluding the edge indicator Tij. and v is updated on the fly given the current output of the chain. The Rao-Blackwellisationn is however infeasible in the context of DAG sampling because it entails computing n² model probabilities at each iteration of the MCMC. This would result an impractical computational scheme when working with many nodes. Alternatively, one could take the ergodic sum over the output (y(1)] ]i=1 as
(1) = Σ 0 =1].
(8)
but the ergodic sum would be too slow to converge. This updating scheme would also lead to other issues such as poor exploration and mis-specification of the important edges. It is worth mentioning that the 17's do not have to be the exact PEPs, as we can still draw valid samples from the target posterior distribution π as long as the MH acceptance probability preserves r-rrversibility. However, choosing 77's that are close to the true PEPs will cause each component to change and to target the right proportion of their marginal probabilities (higher probability of including important edges).
Our solution then to approximate the PEPs म.y before actually running the PARNI-DAG's MCMC, and use this approximation to warm-start the chain. While running the algorithm at time t, we update
5

---

the tuning parameters My adaptively based on
ش. +1 (1-
where {ói] is a decreasing sequence of weights which control the trade-off between the provided warm-start PEPs Hij and the ergodic sum of the output Rly. The choice of [ is not unique and we give a full specification of our choice in the supplementary material (Appendix B). In what follows, we will briefly describe the methods used to calculate the warm-start approximation Tzy.
We argue that the iterative procedure used to restrict the initial search space in Kuipers et al. [2022] can be particularly useful in our case efficiently approximate the PEPs before running the chain. Following the same notation as in Kuipers et al. [2022], let H. be the adjacency matrix underlying a skeleton graph obtained from any DAG learning algorithm (e.g. PC, GES, etc.). For node X,. the
adjacency matrix H defines the permissible parent set h = {X, Hij 1}. All possible parent combinations of node X2 are included in the collection of permissible parent set which is defined by
h [m]] m C J.
(10)
The collection of permissible parent set h can be then extended to a bigger collection h² as described in Section 4.3 of Kuipers et al. [2022]. The extended collection of permissible parent set ht, essentially includes one additional node from outside the permissible parent set h as an additional parent, otherwise the candidate parents of X, restricted to the set h After deriving hy via this procedure as in Kuipers al. [2022]. we then approximate ßij through the following procedure:
1. We remove the acyclicity constraint from the prior (2) and compute the unconstrained posterior distribution x² (1), approximating the marginal probability x² ("iy) under h
2. We approximate the joint probability of a pair edges by the product
# ((rg 79x) R x =(1g)) x ="(s))-
(11)
3. We impose the simplest acyclicity constraint which avoids creating two edges with opposite directions, i.e. let ( = 1, Tsx = 1) Approximate
# 1, Tja = 0)
(12)
# ( = 1, 7yz 0) + #( = 0, 7js = # 0, =)
The full details for the calculations of Fiy are given in the supplementary material (Appendix C).
3.3 Introducing the "reversal" neighborhood
In the original PARNI proposal for Bayesian variable selection, each K, in the intermediate proposals corresponds to exactly one position (i.j). such that ky 1 and that each intermediate neighborhood (((((r - 1), K₂) only contains two DAGs: 7(r- - 1) itself and the DAG obtained by flipping K,, in 7(r (((r-11 1). This sequential neighborhood construction is clearly not appropriate for the problem of DAG sampling and might result in very slow mixing. In addition, it is natural to assume that a pair of edges i <>> > j are highly correlated and thus equally likely be included in a DAG. It is then necessary to specify a reversal move that achieves good mixing between pairs of edges. The neighborhood construction in the original PARNI proposal does not support a reversal move in one single sub-prope aub-proposal. To elaborate on this, suppose that currently Vij = 1 and 7js = 0. to propose reversal move according to the original PARNI proposal, we should first propose to remove Yi from the DAG, then include "/j in the following tub--roposals. These flips might be running into issues if the posterior model probability for excluding both edges, Vis = Tip+ = 0, is low, so that these moves rarely proposed and the chain gets stuck at Yk- = 1.
Therefore, we construct a bigger neighborhood containing the edge reversal move of DAG. This modification allows for an edge reversal move within a single sub cub-proposal Kg.K, without going through the low posterior probability models and avoids the chain getting stuck. To clarify this, suppose l: is sampled at each iteration of the chain. For all those pairs which satisfy k; = 1 and kj; = in k. we construct the reversal neighborhoods for edges Ki and 75t- For other positions such that kij- = 1, but kyę = O, we construct the original neighborhood which only contains the current model and model with Tay flipped. This reversal neighborhood specification is included into the main PARNI-DAG proposal given in the supplementary material (Appendix D).

---

3.4 New adaptive scheme for befter computational efficiency
Although the point-wise implementation reduces the per-iteration costs from exponential to linear in the neighborhood size (equivalent to the number of sub-net nelb--eighboroodds in PARNI), the PARNI proposal is still computationally expensive when the sampled DAGs are not very sparse. For example, considering the gsim100 simulated data of Suter et al. [2023], featuring 161 true edges, we would be often required to compute 320 posterior DAG probabilities on average a single MCMC iteration, and the computation thus becomes 320 times more expensive than the standard ADR proposal in this case. To reduce the computational cost of the PARNI-DAG proposal even more, we make use of a "neighborhood thinning' parameter w to control the number of neighborhoods to be evaluated.
In the original PARNI proposal, the thinning parameter w represents the random walk jumping probability, so it does not take the posterior inclusion probabilities into account and does not affect the number of neighborhoods evaluated. For the new PARNI-DAG proposal instead, with probability 1 - w, we will skip evaluation of some of the neighborhoods. Within each iteration, for each point r in the sequence, given the current sub-proposal 7(r - 1) and sub-neighborhood indicator K,, PARNI-DAG:
1. with probability w. proposes 7()) ~ qK, (T(r - 1); )) as in (5).
2. with probability 1 wi does evaluate the posterior probabilities of other DAGs in N(1(r
1), K,), and directly propose 7() = 7(r 1).
Varying w vary the number of posterior DAG probability evaluations. The larger the value of w,
the more the DAG probabilities evaluated. If w = 1, all neighborhoods are required to be evaluated to propose new DAG.
The parameter w is updated adaptively. Here, we consider the Robins-Monro adaptation scheme to update w on the fly. At time f, the new (f+1)) parameter is updated via
logit, (f+1) logit, (NG - M)
(13)
where Nú is number of sub-neighborhoods evaluated at time t and N is a pre-specified expected number of evaluations. We set th = 0-7.7 meaning that the adaptation of diminishes at the rate of O(--0.7). From the empirical results, we find that N 10 evaluated neighborhoods on average seems to achieve optimal balance between mixing performance and computational cost, Therefore, we choose to set N = 10 for most of the numerical studies in the next section and in the supplementary material.
4 Experiments
4.1 Convergence and mixing efficiency
Protein data. We first consider the real-world protein-signallinn signalling dataset [Sachs al., 2005], found also in Cundy et al. [2021], to test PARNI-DAG's mixing. The dataset consists of n = -1 nodes and N = 853 observations and the ground-trutt DAG (17 edges) is provided by expert knowledge. We compare the performance of PARNI-DAG to ADR, Order MCMC and Partition MCMC. We provide the full details on the implementation of the compared models in the supplementary materials (Appendix E). For a fair comparison of the models, we re-define the score functions for Order and Partition MCMC such that all schemes target the same posterior distribution as specified in Section 2. Notice that for this dataset, it is feasible for all the four MCMC schemes to use the full skeleton, without the need to restrict the initial space [Kuipers et al., 2022]. For ADR and PARNI-DAG, we use prior parameters g 10 and h 1/11.
We compare trace plots of log posterior DAG probabilities for ADR, PARNI-DAG and Partition MCMC in Figure 1. Since we could not extract the full score trace of Order MCMC using the BiDAG package [Suter et al., 2023], Order MCMC is not included in Figure! We ran ADR, PARNI-DAG and Partition MCMC for 480,000 5,0,000,000 60,000 and 20,000 iterations respectively which take about the same CPU time?. and thin the output from ADR and PARNI-DAG so that 20,000 measures are kept. From Figure 1, all these three algorithms mix equally well, the only difference is that the Partition
Using Intel 17 2.80 GHz processor
7

---

Add-Ge
-25200
-2730
algorithm
MSE
posterio modes probatiity
25750
ADR
8.84
-25720
PARNI-DAG
19.44
-25730
Order MCMC
145.17
-25750
Partition MCMC
121.83
co
Table 1: Protein dataset: time-normalised median MSE (x10-$) esti- mating posterior edge probabilities. (Lower is better.
-25730
-25740
-25750
Iterations
Figure 1: Protein dataset: trace plots of log posterior model probabil- ities. All algorithms were ran for the same CPU-time and thinned to 20,000 measures.
MCMC algorithm explores more DAGs with lower probabilities. We then compare the median mean squared errors (MSE) on estimating PEPs compared to a ground-truth estimate. The MSE assesses both the convergence of the chain and the bias in the Monte Carlo estimates. The ground-truth PEPs were obtained by running Partition MCMC for approximately hours. We ran all algorithms for 20 replications. Each individual chain was for minutes. The MSE estimates presented in Table 1. The structure MCMC schemes (ADR and PARNI-DAG) generally outperform Order MCMC and Partition MCMC and return estimates with lower MSE. Notice that ADR in this specific also outperforms PARNI-DAS DAG, since the dataset is very low-dimensional (1 nodes), with PARNI-DAC trailing very far behind. We will demonstrate in the next experiments how PARNI-DAG instead performs better in higher-dimensiona. settings.
gsim100 data. We also study the performance of ADR, PARNI-DAG, Order MCMC and Partition MCMC on a more complex graph, the gsim100 simulated dataset, found in the BiDAG package [Suter al., 2023]. This features randomly generated DAG with 100 nodes and N 100 observations, with 161 true edges. For ADR and PARNI-DAG we set prior parameters to 10 and 1/100. Since 20, is feasible consider two skeletons to restrict the initial space: PC derived skeleton Hpc and the skeleton derived with the iterative procedure of Kuipers et [2022] Hiteer featuring more edges than Hpc.
the full skeleton as in the protein dataset. So
We examine the trace plot of log posterior DAG probabilities from ADR, PARNI-DAG and Partition MCMC using both skeletons Hyc and Harr respectively (Figure 2). We ran ADR, PARNI-DAG and Partition MCMC for 80,000, 20,000 and 20,000 iterations respectively. Notice that, under both Hpc and Haee. PARNI-DAG converges very quickly to the probabbbbaiilittcc region, whereas ADR takes longer. Partition MCMC under Hre cannot reach the same gigh-probability region as targets posterior distribution restricted HPc. which biased if Hec excludes positive edges. Under the more populated Huer instead it converges, but still slower than PARNI-DAG. When comparing median MSEs PEPs this dataset, considered 20 replications with each individual chain being run for 15 minutes. As for Order and Partition MCMC. Hpc seems to lead to slightly better performance than Hirr. This might be due the fact that in this high-dimensional example Hater includes considerably more DAGs than Hpc, so Order MCMC, and even more so Partition MCMC, take longer to explore the space. As for ADR and PARNI-DAG, the choice of H does not have big impact, as in their case this is used only to warm-start the chain and not to restrict the space. PARNI-DAG here outperforms ADR and is competitive Order MCMC as well, without introducing biases stemming from search space restriction and sampling.
4.2 DAG learning accuracy
Lastly, we compare performance of PARNI-DAG with other main methods in recovering the underly- ing true DAG on four different real-world graph structures. These graph structures are taken from the
8

---

Hyz
skeleton
algorithm
MSE
log posteror crodes probebiity
ADR
2,91
-14020
PARNI-DAG
1.83
Hyc
Order MCMC
1.79
-20000
Partition MCMC
12.75
ADR
2.37
PARNI-DAG
1.72
20000
Hiuer
Order MCMC
1.84
Partition MCMC
31.16
5003
12000
15200
70030
Teratons
15000
15000
20000
Table 2: table
gsim100 dataset: time-normalised median MSE 10 ?) on estimating posterior edge probabilities. (Lower is better.)
Algorithm
PAR
POPTSOMMCC
Figure 2: gsim100 dataset: trace plots of log posterior model probabilities. All algorithms were ran for the same CPU-time and thinned 20,000 measures.
R package bnlearn' and include: i) ecoli70 dataset, with = 46 nodes and 70 edges; ii) magic-niab dataset, with n = 44 nodes and 66 edges; iii) magic-irri dataset, with n = 64 nodes and 102 edges; iv) arth/50 dataset, with
107 nodes and 150 edges. Given the graph structure and its associated parameters, we simulate N = 100 Li.dd observations for each of the four cases. The models we compare include: PC algorithm; ii) GES; iii) LiNGAM; iv) DiBS, differentiable variational method introduced by Lorch et al. [2021]: v) DiBS+, a version DiBS with the weighted particle mixture; vi) DECI, another non- non-linear additive noise model that uses variational inference introduced by Geffner et al. [2022]; vii) Order MCMC, on a restricted space defined by the skeleton of the PC algorithm: viii) Iterative MCMC, the version of Order MCMC run on search space defined by starting with the PC skeleton, and iteratively expanding the possible parent sets until the score cannot be improved [Kuipers et al., 2022]; ix) Partition MCMC, on the same space defined by the iterative procedure in Iterative MCMC; x) ADR scheme and xi) PARNI-DAG, both with pre- tuning of the 17's using the Iterative-derieed skeleton. All the score-based methods (GES, Order, Partition, Iterative) assigned an edge-penalty parameter equal to 21)g(n)). as suggested in [Kuipers et al., 2022], to adjust to potentially sparse graphs. Performance is measured via the Structural Hamming Distance between the true DAG and the estimated DAG, averaged over 20 replications of the experiment, for each dataset. Results are collected in Table 3. Notice how PARNI-DAG is either the best or the second best method across all four experiments. In particular, it statistically matches the performance of the resulting best MCMC scheme the first three experiments featuring a lower number nodes, while it significantly outperforms all the methods in the case with high number of nodes (arth/ /50), as expected. All variational methods (DiBS, DiBS+ and DECI) struggle to detect edges in these sample setups they require stronger signal (or larger sample, as featured in the experiments in Lorch et al. [2021], Geffner et al. [2022]) to perform well. Additionally, they are significantly slow converge and run into numerical instability issues the last large-n dataset (arth/50).
Limitations
Although the PARNI-DAG proposal has demonstrated significant improvements the state-of- the-art algorithms considered in the previous section, it still comes with a few limitations that can be addressed. The first limitation lies in the linearity assumption. In fact, in order to derive the closed-form marginal likelihood featured in the locally informed proposal, PARNI-DAG is currently constrained to the linear model assumption described in (1). This can potentially result in model misspecification in presence of non- non-linee linear relationships, which is threat in cases where this hampers edges identifiability. The recent work of Liang et al. [2023] has shed light on how efficiently extend the PARNI proposal to more general non-lineat models, and can possibly be adjusted to Bayesian structure learning settings that employ Additive Noise Models (ANMs) [Hoyer et al., 2008a], where a closed-fom marginal likelihood does not exist. The second limitation of the PARNI-DAG proposal is the computational cost associated with the likelihood evaluation, which scales at least linearly in the
More information sttp://bnlearn. com/.

---

ecoli70
magic-niab
magic-irri
arth150
PC
62.65 + 1.64
66.95 + 0.82
103.75 + 0.91
126.70 + 2.10
GES
47.95 = 1.72
65.50 = 0.42
97.60 = 0.99
129.70 + 3.47
LINGAM
107.4 = 3.96
71.60 +: 1.34
111.25 +: 2.34
DiBS
71.05 + 0.17
65.80 + 0.24
102.35 + 0.17
DIBS+
71.25 + 0.28
64.70 + 0.36
100.45 + 0.40
DECI
70.95 = 0.02
66.95 = 0.02
102.00 = 0.00
Order MCMC
42.85 + 1.30
62.05 1 1.52
92.90 +: 1.66
110.00 +: 1.96
Iterative MCMC
39.55 + 2.15
65.00 + 2.43
96.15 + 2.11
119.45 + 3.43
Partition MCMC
48.85 + 2.76
66.20 + 1.62
103.35 + 1.81
140.75 + 4.04
ADR
39.40 = 1.27
64.25 = 0.64
95.20 = 1.16
143.15 + 1.62
PARNI-DAG
37.55 = 1.59
64.14 + 0.55
91.85 = 0.97
97.60 + 2.16
Table SHD averaged over 20 replications with associated 95% confidence intervals of the various models compared on the four datasets (n = 100). All algorithms were ran for approximately the same CPU-time and the significantly best performing models are given in bold. (Lower is better)
number of datapoints. The issue arises as the likelihood features both in the informed proposals and in the MH step. A viable solution to this problem consists in the possibility of coupling PARNI-DAG with optimal sub-sampling MCMC procedures, such as the ones presented in Korattikara al. [2014] and Maclaurin and Adams [2014] (FireflyMC). These limitations associated with the PARNI-DAG proposal provide us with interesting future research directions, aimed at improving its scalability and capability to handle more complex data.
6 Discussion
In this work, we proposed a novel MCMC sampler, PARNI-DAG, for fully Bayesian approach to the problem of structure learning. PARNI-DAG samples directly from the space of DAGs and introduces improvements in terms of computational complexity and MCMC mixing. stemming from the nature of its adaptive random neighborhood informed proposal. This proposal facilitates moves to higher probability regions, and is particularly useful when dealing with high-dimensional settings with high number of nodes, demonstrated in the experimental section. In future work, will consider adding global moves to the PARNI-DAG proposal, which re-generates a new DAG model every few iterations, specifically tackle settings with highly correlated structures.
In addition, we are also interested in studying PARNI's theoretical mixing time bounds. In fact, although the mixing time bound for random walk proposals has been largely studied, similar results for the class of locally informed proposals are relatively under-developed, due to the complexities arising in the proposal distributions. Most results mixing time in discrete sample spaces focus on the problem of Bayesian variable selection. It has been shown by Yang et al. [2016] that a random mixing under
walk proposal (specifically the add-delete- swap proposal) can achieve polynomial-tinge mild conditions on the posterior distribution. Under similar conditions, Zhou et al. [2022] has shown that the mixing time of the Locally Informed and Thresholded proposal (LIT) does not depend on the number of covariates. The mixing time bound of MCMC samplers in Bayesian structure learning settings is notoriously harder study, due the higher complexity of the sample space of DAGs. The recent work of Zhou and Chang [2023] has shed light on the analysis of mixing time bounds for Bayesian structure learning as a generalisation of the results in Yang al. [2016]. Zhou and Chang [2023] has proven that the mixing time of the Random Walk Greedy Equivalent Search (RW-GES) proposal is at most linear in the number of covariates and the number of datapoints. Moreover, they also presented the necessary conditions for posterior consistency in Bayesian structure learning. It has not been formally proven yet that the informed proposal achieves faster theoretical mixing time compared to the random walk proposal in the context of Bayesian structure learning, although the empirical experiments we conducted suggest so. Similarly, it is highly likely then that the PARNI-DAG proposal can also achieve dimension-fre mixing, such as the LIT proposal. We leave the topic of theoretical mixing time bounds of the PARNI proposal in different applications (c.g., Bayesian variable selection, Bayesian structure learning) to future research.
10

---

Acknowledgements
SL supported by EPSRC research grant EP/VO55380/L. AC's research is funded by the Defence Science and Technology Laboratory (Dstl) which is an executive agency of the UK Ministry of Defence providing world class expertise and delivering cutting-edge science and technology for the benefit of the nation and allies. The research supports the Autonomous Resilient Cyber Defence (ARCD) project within the Dstl Cyber Defence Enhancement programme.
References
David Maxwell Chickering. Learning Bayesian networks NP-complete. Learning from data: Artificial intelligence and statisties V, pages 121-130, 1996.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of machine learning research, 3005v4,507-554, 2002.
Chris Cundy, Aditya Grover, and Stefano Ermon. BDC nets: Scalable variational approaches for Bayesian causal discovery. Advances in Neural Information Processing Systems, 34:7095-7110, 2021.
Mathias Drton and Marloes H. Maathuis. Structure Learning in Graphical Modeling. Annual Review of Statistics and Application, 4(1):365-339) 2017.
Byron Ellis and Wing Hung Wong. Learning causal Bayesian network structures from experimental data. Journal of the American Statistical Association, 103(482):778-788) 2008.
Nir Friedman and Daphne Koller. Being Bayesian about network structure. A Bayesian approach structure discovery in Bayesian networks. Machine learning, 50:95-125, 2003.
Tomas Geffner, Javier Antoran, Adam Foster, Wenbo Gong, Chao Ma, Emre Kiciman, Amit Sharma, Angus Lamb, Martin Kukla, Nick Pawlowski, et al. Deep end-to-end causal inference. arXiv preprint 20Xii:22022.22995, 2022.
Dan Geiger and David Heckerman. Learning Gaussian Networks. In Proceedings of the Tenth International Conference on Uncertainty in Artificial Intelligence, page 235-243, 1994.
Dan Geiger and David Heckerman. Parameter priors for directed acyelic graphical models and the characterization of several probability distributions. The Annals of Statistics, 1015):1412-1440,
2002.
Paolo Giudici and Robert Castelo. Improving Markov chain Monte Carlo model search for data mining. Machine learning, 50: 127-158, 2003.
Clark Glymour, Kun Zhang, and Peter Spirtes. Review of Causal Discovery Methods Based on Graphical Models. Frontiers in Genetics, 10, 2019.
JE Griffin, KG Łatuszyáski, and MFJ Steel. In search of lost mixing time: adaptive Markov chain Monte Carlo schemes for Bayesian variable selection with very large p. Biometrika, 108(1):53-69, 2021.
Marco Grzegorezyk and Dirk Husmeier. Improving the structure MCMC sampler for bayesian networks by introducing new edge reversal move. Machine Learning, 7 71(2-3):265, 2008
WK Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1):97, 1970.
Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21, 2008a.
Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Koller, D. Schuurmans, Bengio, and I.. Bottou, editors. Advances in Neural Information Processing Systems, volume 21, 2008b.

---

Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
Anoop Korattikara, Yutian Chen, and Max Welling. Austerity in meme land: Cutting the metropolis- hastings budget. In International conference on machine learning, pages 181-189, PMLR, 2014.
Jack Kuipers and Giusi Moffa. Partition MCMC for Inference Acyelic Digraphs. Journal of the American Statistical Association, 112651))282-299) 2017.
Jack Kuipers, Polina Suter, and Giusi Moffa. Efficient Sampling and Structure Learning of Bayesian Networks. Journal of Computational and Graphical Statistics, (1)):639-650) 2022.
Xitong Liang, Samuel Livingstone, and Jim Griffin. Adaptive random neighbourhood informed Markov chain Monte Carlo for high-dimensional Bayesian variable selection. Statistics and Computing, 32(5):84, 2022.
Xitong Liang, Samuel Livingstone, and Jim Griffin. Adaptive meme for bayesian variable selection in generalised linear models and survival models. Entropy, 25(9):1310, 2023.
Lars Lorch, Jonas Rothfuss, Bernhard Schölkopf, and Andreas Krause. DiBS: Differentiable Bayesian structure learning. Advances Neural Information Processing Systems, 34:24111-24122 2021,
Dougal Maclaurin and Ryan P Adams. Firefly monte carlo: Exact meme with subsets of data. In 30th Conference on Uncertainty in Artificial Intelligence, UAI 2014, pages 543-552. AUAI Press, 2014.
David Madigan, Jeremy York, and Denis Allard. Bayesian graphical models for discrete data. International Statistical Review/Revue Internationale de Statistique, pages 215-232, 1995.
Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller. Equation of state calculations by fast computing machines. The journal of chemical physics, 1186::08771002, 1953.
Teppo Nüinimäki, Pekka Parviainen, and Mikko Koivisto. Structure discovery in Bayesian networks by sampling partial orders. The Journal of Machine Learning Research, 17(1):2002-2088, 2016.
Judea Pearl. Causality. Cambridge university press, 2009.
Gareth O Roberts and Jeffrey S Rosenthal. Optimal scaling of discrete approximations to Langevin diffusions. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 60(1): 255-268, 1998.
Robert W Robinson. Counting unlabeled acyclic digraphs. In Combinatorial Mathematics V: Proceedings of the Fifth Australian Conference, pages 28-43, 1977.
Karen Sachs, Omar Perez, Dana Pe' er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein- signaling networks derived from multiparameter single-cell data. Science, 308(5721):523-529,
2005.
Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, and Antti Kerminen. A Linear Non-Gaussian Acyclic Model for Causal Discovery. Journal of Machine Learning Research, 2032),2003-2000, 2006.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
Polina Suter, Jack Kuipers, Giusi Moffa, and Niko Beerenwinkel. Bayesian Structure Learning and Sampling of Bayesian Networks with the R Package BiDAG. Journal of Statistical Software, 105 (9):1-31, 2023.
Willem van den Boom, Alexandros Beskos, and Maria De lorio. The G- G-sharrrt weighted proposal al- gorithm: Efficient posterior computation for Gaussian graphical models. Journal of Computational and Graphical Statistics, 1114):1215-1224, 2022.
12

---

Jussi Viinikka and Mikko Koivisto. Layering MCMC for Structure Learning in Bayesian Networks. In Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI), volume 124, pages 839-848, 2020.
Mayering-MCMC
Yun Yang. Martin J Wainwright, and Michael Jordan. On the computational complexity of high- dimensional baysian variable selection. The Annals of Statistics, (446))2447-25322 2016.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG Structure Learning with Graph Neural Networks. In Proceedings of the 36th International Conference on Machine Learning, volume 97, pages 7154-7163, 2019.
G. Zanella. Informed proposals for local MCMC in discrete spaces. Journal of the American Statistical Association, 115533)):52-865, 2020.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. DAGS with NO TEARS: Continuous Optimization for Structure Learning. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, page 1492-9500 2018
Quan Zhou and Hyunwoong Chang. Complexity analysis of bayesian learning of high-dimensione dag models and their equivalence classes. The Annals of Statistics,
2023.
Quan Zhou, Jun Yang, Dootika Vats, Gareth O Roberts, and Jeffrey S Rosenthal Dimension-fie mixing for high-dimensional Bayesian variable selection. Journal of the Royal Statistical Society Series B, 84(5):7511-1784 2022.
13


Reviews:
[{"review_id": "shFsE6Gklz", "reviewer": "Reviewer_xD8f", "summary": "This paper presents a Markov chain Monte Carlo (MCMC) sampler adapted from previous PARNI sampler, called PARNI-DAG, which is designed for Bayesian structure learning under observational data. The authors assume causal sufficiency and target the posterior distribution on Directed Acyclic Graphs (DAGs). The proposed PARNI-DAG mainly relies on the PARNI sampler and modify it to suite the purpose of structure learning, including 1. warm-start of neighbourhood sampling probability 2. reversal neighbourhood step and 3. adaptive neighbourhood skipping probability.\n\nThe main contributions of the paper are\n\n1. A warm-start procedure of the sampler's parameters that exploits skeleton graphs derived through constraint-based or scoring-based algorithms, ensuring better scalability and mixing property with the number of nodes.\n\n2. A reverse step to avoid getting trapped in the local mode.\n\n3. An adaptive skipping probability such that not all intermediate neighbourhood sampling steps are executed. \n\nEmpirically, the author demonstrates the advantage of PARNI-DAG using real-world examples, which shows advantages compared to the previous MCMC methods. ", "questions": "All the questions are mentioned in the weakness section. ", "limitations": "The author does not explicitly write a limitation section, I have made suggestions to include one like computational challenges and limitations of using linear model. ", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "## Originality\nThe originality of the paper lies in the development of the PARNI-DAG algorithm, which combines the Point-wise Adaptive Random Neighborhood Informed (PARNI) proposal with new features specifically designed for structure learning in the space of Directed Acyclic Graphs (DAGs). The proposed algorithm addresses the challenges of mixing and convergence in high-dimensional settings, which is not adequately addressed by existing MCMC methods for structure learning. This work does not provides a completely new sampling algorithm, instead, modify the existing approaches to suite the purpose of structure learning. \n\n## Clarity\nThe paper is logically structured with clear presentation of the modifications made to the PARNI proposal. The authors have made efforts to provide intuitive explanations and motivate their choices in the development of the PARNI-DAG algorithm. The paper is easy to follow, and the appendices provide additional details on the derivations and calculations.\n\n## Significance\nThe significance of the paper lies in its potential impact on the field of structure learning and causal discovery. The PARNI-DAG algorithm aims to address the challenges of mixing and convergence in high-dimensional settings. The algorithm's improved performance over existing MCMC methods for structure learning should make it a reasonable contribution to the field but there are some limitations, which I will elaborate in the following.", "weaknesses": "## Empirical experiments:\nWhile the experimental results demonstrate the advantages of PARNI-DAG over existing MCMC methods, the experiments primarily focus on the comparison with MCMC based approach. However, in the literature review, the author also mentioned several structure learning approach. Although the main claim of this paper is the improvement over existing MCMC, it is still beneficial to include a comparison to state-of-the-art Bayesian structure learning approach like [1,2,3].\n\n[1] Cundy, C., Grover, A., & Ermon, S. (2021). Bcd nets: Scalable variational approaches for bayesian causal discovery. Advances in Neural Information Processing Systems, 34, 7095-7110.\n[2] Lorch, L., Rothfuss, J., Schölkopf, B., & Krause, A. (2021). Dibs: Differentiable bayesian structure learning. Advances in Neural Information Processing Systems, 34, 24111-24123.\n[3] Geffner, T., Antoran, J., Foster, A., Gong, W., Ma, C., Kiciman, E., ... & Zhang, C. (2022). Deep end-to-end causal inference. arXiv preprint arXiv:2202.02195.\n\n## Limitation with linear model\nThe proposed PARNI-DAG mainly targets at the linear model with Gaussian assumptions because the marginal probability are needed. For general nonlinear model, such integration is not tractable, which rendering the following PARNI-DAG proposal invalid. However, nonlinearity is everywhere in the real-world applications. To make the paper stronger, the author should discuss the implication of using linear model or demonstrates that linearity assumption does not harm the performances too much. This can be achieved by comparing some of the previous mentioned baselines [1,2,3]\n\n## Computational complexity\nAlthough the proposed PARNI-DAG method uses many modifications to reduce the computation cost, this approach relies on the MH step to correct the bias. It is known that MH step scales linearly with the number of datapoints, which can be a huge computational bottleneck for large dataset. I suggests the author should explicitly discuss the computational complexity or limitations, and potential approach to remove this constraints. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the time spent carefully reviewing the paper and their positive comments about its originality and clarity. We proceed by addressing the points raised in the weaknesses section.\n\n1)\tEmpirical Experiments\n\nWe have implemented DiBS and DiBS+ (Lorch, L. et al., 2021), and DECI (Geffner, T. et al., 2022) models on the experimental setups of Section 4.2. However, we have noticed that, perhaps not so surprisingly, these variational methods struggle to detect edges in these low-sample regimes, as they require a stronger signal (or larger samples, as featured in their experiments) to perform well. In the pdf file attached to the global response, we have added a table reporting these methods’ average SHD and number of edges detected in the first three datasets (ecoli, magic-niab, magic-irri), as in the last dataset (arth150), characterised by large number of nodes, they are significantly slow and run into numerical instability issues.\n\n2)\tLimitations of Linear Model\n\nWe have added a new ‘limitations’ section in the manuscript, where we discuss both implications of the linearity assumption and computational complexity of PARNI-DAG. Linearity is naturally restrictive, particularly when the objective of the inference is the joint posterior $p(G, \\theta | D)$ rather than the marginal posterior $p(G | D)$ that we focus on, where model misspecification is a threat only when it hampers edges identifiability - meaning that some edges could potentially still be detected in presence of non-linear relationships. Unfortunately, the PARNI-DAG proposal is not directly extendable to Additive Noise Models (ANM) (Hoyer, P. et al., 2008) where the closed-form marginal likelihood does not exist. As part of a separate, but related, research project we have developed a new PARNI proposal coupled with (approximate) Laplace approximation (Rossell, D. et al, 2021) suitable in contexts where the marginal likelihood is not closed form. We have shown that it works well empirically in Bayesian variable selection tasks in generalised linear models and survival models. This provides a solid base to potentially extend the PARNI-DAG sampler to ANM with approximate marginal likelihood, although we leave this for future research. \n\n3)\tComputational Complexity \n\nThe other limitation of the PARNI-DAG proposal pointed out by the reviewer, and discussed in the newly added limitation section, is the computational cost associated with the likelihood evaluation, which scales at least linearly in the number of datapoints. The issue arises as the likelihood features both in the informed proposals and in the MH step. A potential solution to this problem, that we discuss in the new section, regards the possibility of coupling PARNI-DAG with sub-sampling MCMC procedures, such as the ones presented in Korattikara, A. et al. (2014) and Maclaurin, D. & Adams, R. (2015) (FireflyMC). We are actively investigating this topic as a part of a separate ongoing research project.\n\nWe look forward to hearing back from you, \\\nThe Authors\n\n\n[1] Hoyer, P., Janzing, D., Mooij, J.M., Peters, J. and Schölkopf, B., 2008. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21.\n\n[2] Rossell, D., Abril, O. and Bhattacharya, A., 2021. Approximate Laplace approximations for scalable model selection. Journal of the Royal Statistical Society Series B: Statistical Methodology, 83(4), pp.853-879.\n\n[3] Korattikara, A., Chen, Y. and Welling, M., 2014, January. Austerity in MCMC land: Cutting the Metropolis-Hastings budget. In International conference on machine learning (pp. 181-189). PMLR.\n\n[4] Mclaurin, D. and Adams, R., P., 2015, Firefly Monte Carlo: Exact MCMC with Subsets of Data. Proceedings of the 24th International Conference on Artificial Intelligence (pp. 4289–4295). IJCAI'15.\n\n"}, {"review_id": "Mfticji8ub", "reviewer": "Reviewer_1zPb", "summary": "This paper proposes PARNI-DAG, a new MCMC method for sampling Directed Acyclic Graphs (DAGs) that can be used for the problem of structure learning under observational data. PARNI-DAG builds on top of PARNI, and similarly uses locally informed, adaptive random neighborhood with an efficient point-wise implementation, but introduces additional improvements include: pre-tune sampler parameters using a skeleton graph derived from other methods, augmenting the search neighborhood with an edge-reversal move, and neighborhood thinning to improve computational efficiency. Experiments on some toy datasets demonstrate advantage over existing baselines.", "questions": "See weaknesses", "limitations": "The authors did not discuss limitations", "rating": 4, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 1, "strengths": "The proposed method seems to be an effective adaptation of PARNI to do Bayesian learning on DAGs, and outperfoms a few (classical) baseline methods on a few toy benchmarks.", "weaknesses": "I am not an expert in this field, but it seems to me this paper is largely adapting the existing PARNI method for Bayesian variable selection for DAG learning, and combining various techniques that are already present in the current literature on top of PARNI. It would be helpful if the authors can make a more compelling case for the contributions of the paper over existing work.\n\nFor example, in L62-L72 contribution, it seems up until L69 it is just describing what PARNI already has. The procedure for pre-tune sampler parameters and do warm start also seem like a trivial application of exsiting method. In Section 3 titled \"The novel PARNI-DAG\" proposal, the entire Section 3.1 seems to be just a recap of wht PARNI already has. Section 3.2 is mostly a recap of Kuipers et al. [2022]. Section 3.3 introduces the reversal neighborhood but this is also done in for example partition MCMC. Section 3.4 also seems like a simple adaptation of the thinning procedure that is already present in PARNI. In fact a large part of Section 3 seems to belong to background review rather than a description of a novel method. It would be helpful if the authors can reorganze this way and also clearly state what exactly is the novel contribution of PARNI-DAG.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the time spent reviewing our submission. In what follows, we address concerns about novelty and originality of the work.\n\nThe modifications made to the original PARNI proposal, in order to adjust the sampler from a pure variable selection setting to the more complex structure learning one, are several, and non-trivial. Applying a vanilla PARNI proposal to a DAG learning problem results in sub-optimal performance and particularly poses serious concerns about the sampler’s scalability, due to the complexity of DAGs space. For this reason, the adjustments we introduce are necessary to make the PARNI sampler available to the structure learning community to use.\n\nWe proceed by addressing some of the specific points raised by the reviewer:\n\n•\tIn Section 3.2, the procedure we introduce is fundamentally different from the one described in Kuipers, J. et al. (2022). While they make use of the PC algorithm to restrict the starting search space, we instead utilise it to compute warm-start estimates (instead of the computationally intractable Rao-Blackwellised estimates used in original PARNI proposal) for the posterior edge probabilities $\\pi (\\gamma_{ij} = 1)$. While the former procedure introduces bias due to true positive edges deletion, ours does not, as the sampler can potentially still revert back and target the true posterior distribution.\n\n•\tIn Section 3.3, we describe the reversal neighbourhood construction, introduced to improve the sampler’s mixing. Although the idea of a reversal neighbourhood features also in the classical structure MCMC sampler, designing the reversal neighbourhood for the PARNI-DAG proposal is a relevant modification - the way the reversal neighbourhood is constructed in a vanilla PARNI proposal is trivial and would result in extremely slow mixing.\n\n•\tSection 3.4 introduces the parameter adaptation scheme that controls the neighbourhood size and can thus significantly reduce the computational cost. The novel PARNI-DAG’s adaptation scheme comes with a completely different objective from the classic PARNI proposal, as it attempts to control the neighbourhood size directly. \n\nWithout these adjustments, vanilla PARNI proposal would not be efficaciously applicable to DAG learning problems, as it would suffer from extremely slow mixing issues.\n\nWe look forward to hearing back from you,\\\nThe Authors\n"}, {"review_id": "P60b5FZbjk", "reviewer": "Reviewer_PZBP", "summary": "By adding various elaborations to the state-of-the-art Markov chain Monte Carlo inference algorithm, this paper achieves an efficient and effective Bayesian network inference algorithm. Bayesian networks are one of the main tools of machine learning with a long history. In general, their learning is known to be a computationally hard problem, but efficient inference algorithms based on Markov chain Monte Carlo have been actively studied. Section 2 describes the problem setup of a linear functional Bayesian network with a normally distributed weight matrix as a concrete example. Section 3.1 introduces the baseline algorithm as a state-of-the-art inference algorithm for Bayesian networks, applying the method of Liang et al. [2022]  for Bayesian variable selection problems to Bayesian networks. In Sections 3.2 through 3.4, the authors' further elaborations are carefully described step by step. Section 3.2 shows how to effectively adjust the neighborhood of the proposals, inspired by Liang et al. [2022]  for Bayesian variable selection problems. Section 3.3 introduces an effective device for (nested) sequential sampling of neighborhoods to eliminate problems that can arise in the DAG estimation problem, but not in Bayesian variable selection. Section 3.4 presents a method for properly restricting the neighborhoods, which tend to be enormous.", "questions": "Thank you for sharing this very interesting paper. This paper has a lot of devices to help the reader understand. I really enjoyed reading about many of the parts. On the other hands, for some parts, I also worry that perhaps I am underestimating the value of this paper due to my own lack of understanding. Therefore, I would like to present some arguments below to improve my own understanding. Authors do not need to respond to the inconsequential ones, but if there is anything that authors infer that my understanding is lacking, it would be very helpful if you could respond.\n\n- I find the author's argument against mixing times in the proposed algorithm somewhat unclear.\n\nIn this paper, there are several mentions of mixing (e.g., Line 13, 288, 365). As the author states \"PARNI-DAG quickly converges to high-probability regions (Line 11),” I interpret these claims as \"the proposed algorithm achieves more rapid mixing times or shorter arrival time expectations to regions of high probability (note: these two claims are equivalent (Theorem 1.4 of [*]))”.\n\n[*] Yuval Peres and Perla Sousi. Mixing times are hitting times of large sets. Journal of Theoretical Probability, 28(2):488–519, 2015. \n\nHowever, I think that in general (especially when it depends on input data) analyzing or empirically observing MCMC mixing times is not an easy task.\nTo my knowledge, there are not many examples of MCMC mixing times being analyzed for many problems in general. The exceptions are the following few problems:\n- Bayesian variable selection:  Yun Yang, Martin J. Wainwright, and Michael I. Jordan. On the computational complexity of high- \ndimensional Bayesian variable selection. The Annals of Statistics, 44(6):2497 – 2532, 2016. \n- Bayesian community detection: Bumeng Zhuo and Chao Gao. Mixing time of Metropolis-Hastings for Bayesian community detection. Journal of Machine Learning Research, 22:10:1–10:89, 2021,\n- Structure learning with directed acyclic graph: [+] Quan Zhou and Hyunwoong Chang. Complexity analysis of Bayesian learning of high-dimensional  DAG models and their equivalence classes, Annals of Statistics, 2023 (or arXiv:2101.04084).\n\nFortunately, the DAG inference problems addressed in this paper seem to be able to guarantee polynomial-time mixing times for certain algorithms (Theorem 6 of the above literature [+]). \nIn this light, what observations can be made about the mixing time of the proposed MCMC algorithm?\nCan existing results for mixing times be easily applied to the proposed algorithm? Or is theoretical analysis of the proposed algorithm for mixing time a future issue? Also, as we can observe from the experimental results, can the proposed algorithm speed up the mixing time by an order of magnitude that depends on the problem size? Or is it an improvement of constant orders?\nFortunately, theoretical analysis of MCMC mixing time for DAGs seems to be in progress, so we hope that the author's mention of these issues (separating what is known from what is not known) will provide important insights for the reader.", "limitations": "As discussed in the above Questions, I am not sure the (theoretical) guarantee of mixing time of the proposed MCMC algorithm. ", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- This paper is a solid proposal for possible reasonable improvements in the Bayesian network inference problem, with broad coverage of the latest developments in the surrounding fields of Bayesian machine learning.\n- The text is very detailed so that a wide variety of readers (from beginners to experts) can follow the history and the latest developments in the field.\n- The code is provided in such a way that the proposed algorithm can be easily followed up by subsequent research, which is a very significant contribution to the field.", "weaknesses": "The candidate weaknesses listed below are based on questions I had during my initial peer review. As my misconceptions are resolved, they may cease to be weaknesses.\n\n- The paper is somewhat unclear in its claims about the mixing time analysis (or empirical observation) of the proposed MCMC algorithm, although there are several mentions of it (e.g., Line 13, 288, 365).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are thankful for the reviewer’s positive comments about the clarity of exposition and the developments brought about by our work. We proceed by answering the reviewer’s clarifying questions about mixing time here below.\n\nIn some parts of the paper, we use the wording ‘mixing time’ to describe the empirical mixing performance observed in the output of simulated Markov chains, rather than the formal theoretical mixing time property of the MCMC algorithm. These numerical results (Section 4) suggest that the PARNI-DAG proposal is very likely to have a faster theoretical mixing time compared to the add-delete-reverse proposal, even though this has not been formally proved. Investigating the mixing time bounds of the PARNI(-DAG) proposal is a very interesting topic, but currently beyond the scope of this work. Nonetheless, to acknowledge this, we have included a brief discussion in the new version of the paper.\n\nAlthough the theoretical mixing time bound for random walk proposals has been largely studied, similar results for the class of locally informed proposals (including PARNI and PARNI-DAG) are relatively under-developed, due to the complexities arising in the proposal distributions. Most results on mixing time in discrete sample spaces focus on the problem of Bayesian variable selection. It has been shown in [1] that a random walk proposal (specifically the add-delete-swap proposal) can achieve polynomial-time mixing under mild conditions on the posterior distribution. Another theoretical result in [2] shows that the informed proposals can achieve dimension-free mixing time bounds under the same mild conditions. Under the conditions that posterior mass concentrates on a small set of models and the chain starts at a model close enough to the underlying “true” model, it has been shown in [2] that the mixing time of the Locally Informed and Thresholded proposal (LIT) does not depend on the number of covariates. \n\nNotoriously, the mixing time of MCMC samplers in Bayesian structure learning settings is harder to study, due to the higher complexity of the DAGs sample space. Recent work [3] has shed light on the analysis of mixing time bounds for Bayesian structure learning as a generalisation of the result from [1]. In [3], the authors showed that the mixing time of the Random Walk Greedy Equivalent Search (RW-GES) proposal is at most linear in the number of covariates and the number of datapoints. Moreover, they also presented the necessary conditions for posterior consistency in Bayesian structure learning. Proving that the informed proposal can achieve faster mixing time than the random walk proposal (as [2] did for the case of Bayesian variable selection) in Bayesian structure learning settings is still an ongoing area of research.\n\nWe conclude by mentioning that, as part of a separate ongoing research, we are currently studying the mixing time of the PARNI proposal in comparison to other state-of-the-art schemes on various applications. From the empirical results in [4], the PARNI proposal has faster empirical mixing time compared to the LIT proposal, and it is highly likely that the PARNI proposal can also achieve dimension-free mixing like the LIT proposal on Bayesian variable selection. Considering the theoretical results in [2], it appears feasible to generalise them to the Bayesian structure learning setting, but we leave this for future research.\nWe are happy to include in the new version of the paper a discussion about the theoretical mixing time of the PARNI-DAG proposal and how it relates to the work mentioned above. The discussion also mentions possible methods (e.g., spectral gap, canonical path analysis and drift-and-minorization methods) that can be employed to find the theoretical mixing time bound of PARNI-DAG on structure learning problems.\n\nWe look forward to hearing back from you, \\\nThe Authors\n\n\n[1] Yun Yang, Martin J. Wainwright, and Michael I. Jordan. On the computational complexity of high- dimensional Bayesian variable selection. The Annals of Statistics, 44(6):2497 – 2532, 2016.\n\n[2] Zhou, Q., Yang, J., Vats, D., Roberts, G.O. and Rosenthal, J.S., 2022. Dimension-free mixing for high-dimensional Bayesian variable selection. Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(5), pp.1751-1784.\n\n[3] Quan Zhou and Hyunwoong Chang. Complexity analysis of Bayesian learning of high-dimensional DAG models and their equivalence classes, Annals of Statistics, 2023 (orarXiv:2101.04084).\n\n[4] Liang, X., Livingstone, S. and Griffin, J., 2022. Adaptive random neighbourhood informed Markov chain Monte Carlo for high-dimensional Bayesian variable Selection. Statistics and Computing, 32(5), p.84.\n"}]

Author rebuttal:
We are thankful to the reviewers for the time spent going through our submission and for their insightful comments, as we believe these have significantly contributed to improving the paper. 

In particular, we have addressed Reviewer xD8f’s request for additional experimental comparison by implementing the variational methods DiBS, DiBS+ and DECI (Lorch, L. et al., 2021; Geffner, T. et al., 2022) on the experimental setups of Section 4.2, and showed that these are not particularly well-suited for these types of low data regimes (a table of the results can be found in the pdf attached this global response). We have also added a new section in the paper discussing the limitations of the linearity assumption and the computational complexity of PARNI-DAG together with potential ways to alleviate this (sub-sampling MCMC solutions), which we are currently investigating.

We have addressed Reviewer 1zPb’s concerns about the originality of the work, by stressing how PARNI-DAG brings about several relevant adjustments to the vanilla PARNI proposal that makes the methods applicable to the higher complexity of DAG spaces. We have stressed also that vanilla PARNI would yield sub-optimal performance and would pose serious scalability concerns in structure learning settings. 

Finally, in response to Reviewer PZBP’s questions, we have included in the new version of the paper a brief discussion on the theoretical mixing time of the PARNI and PARNI-DAG proposals. As mentioned in the last paragraph of the individual response to Reviewer PZBP, as part of a separate ongoing project, we are currently studying mixing time bounds of the classic PARNI proposal in Bayesian variable selection settings, and plan to potentially extend this study to Bayesian structure learning settings in the future, although this is notoriously a harder problem.

We hope we have adequately responded to all the concerns regarding the paper, and we look forward to engaging with the reviewers during the interactive discussion period, to respond to any outstanding queries.

Kindest Regards, \
The Authors


Meta-review: While we didn't reach full consensus among the reviewers, two of the three reviewers were satisfied with the author responses and recommend acceptance.  The authors are strongly encouraged to make the changes the proposed in the replies to these reviewers.  There was good back-and-forth between the third reviewer and the authors, although the third reviewer never replied to the last author reply.  This meta-reviewer found that author reply satisfactory and so is recommending acceptance.  Authors are encouraged to take into careful consideration the comments of that reviewer as well in preparing the final version of the paper.
Decision: Accept (poster)