方向性证据引导的搜索空间缩减用于精确 DAG 学习

# Directional Evidence Guided Search-Space Reduction for Exact DAG Learning
Source: [https://arxiv.org/html/2610.09136](https://arxiv.org/html/2610.09136)
Upala Junaida IslamEmail:[uislam@asu\.edu](mailto:[email protected])Affiliation:School of Computing and Augmented Intelligence, Arizona State University, Tempe, AZ 85281, USAAbdelmonem ElrefaeyEmail:[aelrefaey@asu\.edu](mailto:[email protected])Affiliation:School of Computing and Augmented Intelligence, Arizona State University, Tempe, AZ 85281, USARong PanEmail:[rong\.pan@asu\.edu](mailto:[email protected])Affiliation:School of Computing and Augmented Intelligence, Arizona State University, Tempe, AZ 85281, USA

###### Abstract

Learning a directed acyclic graph \(DAG\) from observational data is a challenging combinatorial problem due to the exponential growth in the number of candidate parent\-set configurations\. Existing exact score\-based methods often require computationally intensive combinatorial search, whereas constraint\-based methods can become unreliable or computationally demanding as graph size and conditioning\-set complexity increase\. We develop a non\-parametric hybrid framework, referred to as DECO \(Directional Evidence\-guided Configuration Optimization\), that extracts dependency and directional evidence from observation data to construct admissible parent sets prior to exact optimization\. It reduces the optimization search space by eliminating empirically unsupported parent configurations while preserving flexibility for all plausible edge orientations\. Theoretical analysis establishes an exponential reduction in the admissible parent\-set configuration space and quantifies how bounded edge\-level omission affects the probability of retaining the true parent structure\. Experiments on benchmark Bayesian networks and synthetic discrete and continuous DAGs demonstrate substantial search\-space reduction while achieving competitive structure\-recovery performance, with favorable structural Hamming distance across many evaluated settings\. These results show that directional evidence can provide an effective preprocessing mechanism for reducing the computational burden of exact DAG learning without requiring a fixed parametric structural model\.

###### keywords

Causal discovery, Causal structure learning, Directed acyclic graph, Exact Bayesian network learning, Integer programming

## 1Introduction

As a core task in causal inference, causal discovery aims to recover causal relationships among variables by hypothesizing and constructing possible parent\-child links leading to a Directed Acyclic Graph \(DAG\)\. Such a task has been found particularly helpful for researchers in selecting appropriate confounders, avoiding bias, and designing effective interventions by identifying true causes rather than mere correlations\([Cunningham, 2021](https://arxiv.org/html/2610.09136#bib.bib11);[Williamson and Forbes, 2021](https://arxiv.org/html/2610.09136#bib.bib17)\)\. This is found in literature across multiple scientific disciplines, e\.g\., genomics\([Pingault et al\., 2018](https://arxiv.org/html/2610.09136#bib.bib21)\), neuroscience\([Seth et al\., 2015](https://arxiv.org/html/2610.09136#bib.bib22)\), economics\([Henningsen et al\., 2025](https://arxiv.org/html/2610.09136#bib.bib23)\), and the social sciences\([Cowls and Schroeder, 2015](https://arxiv.org/html/2610.09136#bib.bib30)\), among other domains\.

In the realm of DAG learning, we can categorize learning algorithms into constraint\-based, score\-based, and hybrid\. Constraint\-based algorithms use Conditional Independence \(CI\) tests to infer the graph that underpins causal or associative interactions among variables\([Spirtes et al\., 2000](https://arxiv.org/html/2610.09136#bib.bib24);[Spirtes, 2010](https://arxiv.org/html/2610.09136#bib.bib27)\)\. These CI\-test methods are unreliable in finite samples and therefore require massive datasets, which may be unrealistic in many applications\([Zhang et al\., 2017](https://arxiv.org/html/2610.09136#bib.bib5)\)\. Score\-based algorithms, on the other hand, optimize decomposable score functions such as the Bayesian Information Criterion \(BIC\) or Akaike Information Criterion \(AIC\)\([Bouckaert, 1995](https://arxiv.org/html/2610.09136#bib.bib26);[Gámez et al\., 2011](https://arxiv.org/html/2610.09136#bib.bib10)\)\. Among the vast literature of DAG search algorithms, the Integer Programming \(IP\), as an exact optimization approach, is a widely used score\-based method\([Bartlett and Cussens, 2017](https://arxiv.org/html/2610.09136#bib.bib25)\)\. Existing IP\-based DAG learning methods primarily focus on identifying the best parent set for each node in the network\([Manzour et al\., 2021](https://arxiv.org/html/2610.09136#bib.bib28)\)\. Starting without any*a priori*knowledge would make the problem explode in the number of variables, and a complete search becomes impossible\. The limitations of score\- and constraint\-based methods motivate hybrid DAG learning approaches that incorporate features from both classes\. A representative example is Max\-Min Hill\-Climbing \(MMHC\), which first uses a constraint\-based procedure to estimate the graph skeleton, then applies a greedy, score\-based search to orient edges and refine the structure\([Tsamardinos et al\., 2006](https://arxiv.org/html/2610.09136#bib.bib7)\)\. Related methods, such as H2PC, show that restricting the search space through CI information can substantially improve scalability and statistical efficiency\([Gasse et al\., 2015](https://arxiv.org/html/2610.09136#bib.bib6)\)\. However, existing hybrid methods primarily rely on dependency information, while directional information is either underutilized or enforced only locally\. They rarely exploit directional evidence as global structural constraints within exact optimization\. As a result, hybrid methods may remain unstable when CI tests are ambiguous, sample sizes are limited, or nonlinear dependence is present\([Scutari et al\., 2018](https://arxiv.org/html/2610.09136#bib.bib1)\)\. They may fail to distinguish noisy local dependencies from structurally persistent directional patterns\. Incorporating reliable directional evidence as structural constraints can significantly reduce the feasible orientation space while preserving flexibility for uncertain edge directions\.

To this end, we present a unified hybrid framework for learning DAGs over both discrete and continuous causal networks, which we refer to as DECO \(Directional Evidence\-guided Configuration Optimization\)\. Here, we first extract an undirected skeleton from the observational data using a statistical dependency measure; then we use a divergence\-based directional measure to infer reliable orientations and restrict infeasible candidate orientations\. Finally, an existing exact score\-based optimizer searches over the resulting restricted parent\-set space to recover the highest\-scoring DAG\.

Our key contributions can be summarized as follows\.

- •We introduce a non\-parametric directional orientation metric based on maximal conditional divergence\. Our metric exploits pairwise distributional asymmetries in both continuous and discrete data, effectively accounting for local shifts in location and scale, to propose edge directions\.
- •We propose a method for converting this directional evidence into admissible parent\-set restrictions that can be supplied to any exact parent\-set optimizer\. DECO translates divergence scores into targeted restrictions on parent sets\. This yields an exponential reduction in the admissible parent\-set configuration space while preserving mathematical flexibility when edge orientations are uncertain\.
- •We provide a theoretical characterization of the resulting search\-space reduction showing that a linear reduction in the candidate parent\-set size per node yields an exponential reduction in the admissible parent\-set configuration space, and quantify how edge\-level omission risk compounds with node in\-degree and graph size in determining retention of the true parent structure

## 2Methodology

In this section, we present the general framework of DECO for causal discovery\. LetVVdenote the set of nodes, or vertices, with\|V\|=D\\left\|V\\right\|=D, andX∈𝒳N×DX\\in\\mathcal\{X\}^\{N\\times D\}be the observed dataset withNNsamples\. Here,𝒳\\mathcal\{X\}denotes the variable domain, which may be continuous or discrete\. The ultimate goal is to learn the underlying graph structure fromXX\. Our proposed methodology consists of \(i\) extracting structural and directional evidence from the observed data and converting them into an admissible parent\-child construction, and \(ii\) applying existing score\-based exact optimization over the restricted space to recover the structure of the highest\-scoring DAG\. We will go into detail on step one in Sections 2\.1 and 2\.2 and on the optimization procedure in Section 2\.3\.

### 2\.1Directional evidence extraction

The primary objective of the preprocessing phase is to systematically prune the feasible graph search space by eliminating empirically unsupported relations, thereby reducing the computational burden on the downstream global optimization\. Hybrid methods such as MMHC\([Tsamardinos et al\., 2006](https://arxiv.org/html/2610.09136#bib.bib7)\)and H2PC\([Gasse et al\., 2014](https://arxiv.org/html/2610.09136#bib.bib16)\)reduce the optimization space by estimating an undirected skeleton using conditional independence tests, leaving the directional decisions largely to subsequent optimization\.[Rolland et al\. \(2022\)](https://arxiv.org/html/2610.09136#bib.bib34)and[Sanchez et al\. \(2022\)](https://arxiv.org/html/2610.09136#bib.bib33)find topological order using Hessian \(or its suitable approximation\) of the log\-likelihood of the data to iteratively identify and remove leaf nodes and use the topological information as a base before seeking causal relations\. However, such methods strictly require differentiable likelihood functions and are largely restricted to parametric models\. To overcome these limitations, we introduce a non\-parametric framework that operates independently of structural functional forms or differentiability assumptions, and it is applicable to both discrete and continuous data distributions\. Our pipeline extracts directional evidence by first isolating plausible dependencies via a skeleton construction phase, ranking candidate edge orientations through local dominance, and subsequently validating them via asymmetry\-based directional inference\.

First, an undirected graph skeleton is established by evaluating the pairwise Mutual Information \(MI\) across all variables, which quantifies the reduction in uncertainty of one variable given another\([Cover et al\., 1991](https://arxiv.org/html/2610.09136#bib.bib15)\)\. Under standard causal Markov and faithfulness assumptions, statistical independence implies the absence of a structural edge\([Spirtes et al\., 2000](https://arxiv.org/html/2610.09136#bib.bib24);[Pearl, 2009](https://arxiv.org/html/2610.09136#bib.bib29)\)\. While[Chow and Liu \(1968\)](https://arxiv.org/html/2610.09136#bib.bib35)utilized an MI\-maximizing framework to construct tree structures, we apply a conservative, user\-defined thresholdδMI\\delta\_\{\\text\{MI\}\}to eliminate candidate edges with negligible mutual dependence\. To prevent structural fragmentation, a local threshold relaxation is dynamically applied to any node isolated at this screening step, thereby connecting it to a small set of candidate neighbors that share the highest mutual information with the isolated node\. This guarantees that every variable retains at least a limited set of admissible parent candidates for subsequent optimization\. Because mutual information is inherently symmetric, the resulting skeleton delineates the undirected constraint space containing all plausible causal interactions without determining edge orientations\. We conducted a preliminary sensitivity analysis over a range of the threshold,δMI\\delta\_\{\\text\{MI\}\}, on representative synthetic graphs\. The selected thresholds provided a favorable trade\-off between preserving true edges and reducing the search space\. Details are provided in Appendix C\.1\.

Next, we use an asymmetrical measure to obtain the directional preference of the preserved edges\.[Duong and Nguyen \(2022\)](https://arxiv.org/html/2610.09136#bib.bib31)proposed a metric, Conditional Divergence based Causal Inference \(CDCI\), for both discrete and continuous data\. CDCI estimates the average divergence between a variable’s marginal distribution and its standardized conditional distributions, given the other variable, to infer the causal direction between the two variables\. This approach removes any strict assumptions on functional forms between variables and their distribution models, and it supports general divergence metrics\. However, CDCI is formulated for pairwise causal directions only and thus is not suitable for determining the complete network structure\. In addition, by standardizing the conditional distributions, CDCI suppresses location and scale differences between the marginal and conditional distributions that may be informative about causal direction\. Moreover, because CDCI averages divergence over all conditioning values, strong directional effects that occur only for a subset of conditions can be diluted by weaker effects elsewhere\. Motivated by CDCI, we propose a pairwise measure, Maximized Conditional Divergence \(MCD\), which maximizes the divergence between the marginal and conditional empirical distributions for all ordered pairs\. If\(Xu,Xv\)\(X\_\{u\},X\_\{v\}\)are two variables from the observational data, we defineM​C​DMCDofXvX\_\{v\}givenXuX\_\{u\}as

M​C​D​\(v\|u\)\\displaystyle MCD\\left\(v\|u\\right\)=maxx∈𝒳u​d​\(P⁡\(Xv\|Xu=x\),P⁡\(Xv\)\)\\displaystyle\{\}=\\underset\{x\\in\\mathcal\{X\}\_\{u\}\}\{\\max\}~\{d\\left\(P\\left\(X\_\{v\}\|X\_\{u\}=x\\right\),P\\left\(X\_\{v\}\\right\)\\right\)\}\(1\)whereddis any valid probability distance measure, e\.g\., total variation, Hellinger distance, etc\. It is important to note that these distributions are purely empirical, derived solely from observational data\. Access to the actual distributions or even its form is not required to findM​C​DMCD\. The intuition is that when conditioned on a true parent, the conditional distributions ofXvX\_\{v\}will diverge from its marginal distribution\. Therefore, if conditioning onXuX\_\{u\}significantly alters the distribution ofXvX\_\{v\}, we interpret it as strong directional evidence supportingXu→XvX\_\{u\}\\xrightarrow\{\}X\_\{v\}than for the reverse direction\.

M​C​DMCDaims to identify which variable exhibits the strongest distributional shift for any configuration\. This measure is inherently asymmetric, asM​C​D​\(v\|u\)≠M​C​D​\(u\|v\)MCD\\left\(v\|u\\right\)\\neq MCD\\left\(u\|v\\right\)in general\. This aligns with a margin\-based view, where larger divergence differences provide stronger evidence for causal direction\. However,M​C​DMCDvalues are not directly comparable across source nodes due to scale differences\. As a result, an edge with relatively low globalM​C​DMCDmay still be the strongest outgoing candidate for a given source node\. To capture such local dominance, we introduce a row\-normalized metricN​M​C​DNMCD, where

N​M​C​D​\(v\|u\)\\displaystyle NMCD\\left\(v\|u\\right\)=\{M​C​D​\(v\|u\)max𝑤​M​C​D​\(w\|u\)if​M​C​D​\(v\|u\)\>00otherwise\\displaystyle\{\}=\\begin\{cases\}\\frac\{MCD\\left\(v\|u\\right\)\}\{\\underset\{w\}\{\\max\}MCD\\left\(w\|u\\right\)\}&\\text\{if \}MCD\\left\(v\|u\\right\)\>0\\\\ 0&\\text\{otherwise\}\\end\{cases\}\(2\)
which detects the strongest relative shift or causal direction from each source node\. This normalization removes scale differences across source variables\. It also places the outgoing directional strengths from each source node on a common scale, allowing local directional dominance to be compared independently of the absolute magnitude of the directional shift\. Importantly,max𝑤​M​C​D​\(w\|u\)=0\{\\underset\{w\}\{\\max\}MCD\}\\left\(w\|u\\right\)=0if and only ifM​C​D​\(w\|u\)=0MCD\\left\(w\|u\\right\)=0for allwwimplying that the marginal distribution ofXwX\_\{w\}is identical to every conditional distributionP⁡\(Xw\|Xu=x\)P\(X\_\{w\}\|X\_\{u\}=x\)which ultimately makesN​M​C​D​\(v\|u\)=0NMCD\\left\(v\|u\\right\)=0by definition\. To orient an edge, we convertM​C​DMCDandN​M​C​DNMCDinto asymmetrical directional preference scores and compare the two directions of each edge of the skeleton\. We define the absolute directional preference \(SS\) and normalized directional preference \(N​SNS\) in the following way,

Su→v\\displaystyle S\_\{u\\xrightarrow\{\}v\}=M​C​D​\(v\|u\)−M​C​D​\(u\|v\)\\displaystyle\{\}=MCD\\left\(v\|u\\right\)\-MCD\\left\(u\|v\\right\)\(3\)N​Su→v\\displaystyle NS\_\{u\\xrightarrow\{\}v\}=N​M​C​D​\(v\|u\)−N​M​C​D​\(u\|v\)\\displaystyle\{\}=NMCD\\left\(v\|u\\right\)\-NMCD\\left\(u\|v\\right\)\(4\)
For continuous variables, we discretize each variable using a mean/standard\-deviation\-based binning scheme, following[Duong and Nguyen \(2022\)](https://arxiv.org/html/2610.09136#bib.bib31), so that MCD and NMCD can be computed from empirical marginal and conditional distributions in the same nonparametric manner as for discrete data\. For a variable with sample meanμ\\muand standard deviationσ\\sigma, the central range\[μ−k​σ,μ\+k​σ\]\\left\[\\mu\-k\\sigma,\\mu\+k\\sigma\\right\]is partitioned into4​k\+14k\+1equal\-width bins, and the two outermost bin edges are extended so that values beyond the central range are absorbed into the adjacent extreme bin rather than discarded\. We set k = 2, yielding 9 bins per continuous variable\. Constant or near\-constant columns are assigned a single bin\. The discretized data are then processed identically to native discrete data for MI\-based skeleton construction and MCD/NMCD\-based directional evidence extraction\. This avoids introducing additional density\-estimation assumptions and keeps the directional\-evidence stage consistent across data types\.

### 2\.2Admissible parent\-set construction

The directional evidence extracted in Section 2\.1 is now converted into admissible parent\-set constraints\. Our edge orientation strategy is based on two principles: First,N​M​C​DNMCDranks candidate outgoing edges from each source node based on their relative local dominance\. Second, directional asymmetry can be supported when either normalized or absolute asymmetry exhibits a sufficiently strong preference for one direction over the other\. Accordingly, we process candidate outgoing edges from each source node byN​M​C​DNMCD\. Then along the descending order of rank, we evaluate directional asymmetry usingN​SNSandSSscores for an edge and its reverse direction\. If there is a clear directional asymmetry such thatN​Su→v≥δN​SNS\_\{u\\xrightarrow\{\}v\}\\geq\\delta\_\{NS\}orSu→v≥δSS\_\{u\\xrightarrow\{\}v\}\\geq\\delta\_\{S\}, we confirm the direction fromu→vu\\xrightarrow\{\}vas long as it does not create a cycle with respect to the current directed graph\. Here,δN​S\\delta\_\{NS\}andδS\\delta\_\{S\}are the two thresholds calibrated to score scales and set as fixed constants in an experiment\. Empirical analyses over a range of threshold values are provided in Figure[4](https://arxiv.org/html/2610.09136#S3.F4)\.N​SNSandSScapture local and global directional preference, respectively\. We use them as complementary measures, allowing an orientation to be confirmed either by strong local dominance or by a sufficiently large absolute asymmetry\. Each oriented edge is assigned a confidence score, obtained by normalizing the directional\-evidence metric \(N​SNSorSS\) associated with its orientation, so stronger directional asymmetry yields higher confidence\. Confidently oriented edges forbid their reverse directions, and for detached nodes it serves as a tie\-breaking ranking term when selecting fallback parents\. For retained edges, confident orientations strictly enforce parent\-child roles, while ambiguous edges maintain bidirectional candidates, allowing the Integer Programming to orient them flexibly during global optimization\.

The resulting reduction in the admissible parent\-set configuration space is quantified in the following theorem\. We provide a sketch proof of the theorem here; a complete proof appears in Appendix A\.

###### Theorem 1\(Reduction in Parent\-Set Configuration Space\)\.

LetVVbe the set of variables with\|V\|=D\|V\|=D\. IfΩf​u​l​l\\Omega\_\{full\}denotes the number of parent\-set assignments prior to enforcing global acyclicity constraints,Ωfull=2D⁡\(D−1\)\\Omega\_\{\\text\{full\}\}=2^\{D\(D\-1\)\}\. Under the proposed DECO preprocessing, letmv≤D−1m\_\{v\}\\leq D\-1denote the number of admissible parents of nodevvandm¯=1D​∑v∈Vmv\\bar\{m\}=\\frac\{1\}\{D\}\\sum\_\{v\\in V\}\{m\_\{v\}\}be the average candidate parent pool size per node\. Then the number of admissible parent\-set configurations isΩD​E​C​O=2∑v∈Vmv=2D​m¯\\Omega\_\{DECO\}=2^\{\\sum\_\{v\\in V\}m\_\{v\}\}=2^\{D\\bar\{m\}\}\. Finally, if the average admissible parent pool satisfiesm¯≤\(1−α\)​\(D−1\)\\bar\{m\}\\leq\(1\-\\alpha\)\(D\-1\)forα∈\(0,1\]\\alpha\\in\(0,1\],

ΩDECOΩfull≤2−α​D​\(D−1\)\\frac\{\\Omega\_\{\\text\{DECO\}\}\}\{\\Omega\_\{\\text\{full\}\}\}\\leq 2^\{\-\\alpha D\(D\-1\)\}

###### Proof Sketch\.

In the unconstrained model, each node may choose any subset of the remainingD−1D\-1variables, yielding2D−12^\{D\-1\}candidate parent sets\. In the DECO framework, the admissible parent pool of nodevvis𝒫v\\mathcal\{P\}\_\{v\}of sizemvm\_\{v\}\. Before imposing global acyclicity constraints, the parent\-set choices can be counted through the Cartesian product over nodes\. Therefore,

ΩDECO=2∑v∈Vmv\\Omega\_\{\\text\{DECO\}\}=2^\{\\sum\_\{v\\in V\}m\_\{v\}\}Applying the assumed average pruning fractionα\\alphagives the following bound on the ratio of the parent\-search configurations:

ΩDECOΩfull=2∑v∈Vmv2D⁡\(D−1\)≤2−α​D​\(D−1\)\\frac\{\\Omega\_\{\\text\{DECO\}\}\}\{\\Omega\_\{\\text\{full\}\}\}=\\frac\{2^\{\\sum\_\{v\\in V\}m\_\{v\}\}\}\{2^\{D\(D\-1\)\}\}\\leq 2^\{\-\\alpha D\(D\-1\)\}Therefore, a linear reduction in the average admissible parent pool induces an exponential reduction in the candidate parent\-set assignments considered by the exact optimizer\. ∎

This theorem bounds the number of admissible parent\-set configurations that the optimizer considers\. It does not bound the maximum in\-degree of any individual node or the number of valid DAGs\. However, reducing the search space is meaningful only if feasible solutions of interest remain accessible\. The following proposition establishes this property\.

###### Proposition 1\(True\-Parent Retention under Bounded Omission\)\.

For every true directed edgeee, let the probability that DECO preprocessing excludes its true parent relation be at mostε\\varepsilon\. Then for a nodevvwith true in\-degreedvd\_\{v\}, the probability that all true parent relations ofvvremain admissible is

P​r​\(P​aG​\(v\)⊆𝒫v\)≥max⁡\(0,1−dv​ε\)Pr\\left\(Pa\_\{G\}\\left\(v\\right\)\\subseteq\\mathcal\{P\}\_\{v\}\\right\)\\geq\\max\{\\left\(0,1\-d\_\{v\}\\varepsilon\\right\)\}Moreover, for a DAG with\|E\|\\left\|E\\right\|true edges, probability that the true DAG remains feasible after the preprocessing is at leastmax⁡\(0,1−\|E\|​ε\)\\max\{\\left\(0,1\-\\left\|E\\right\|\\varepsilon\\right\)\}\.

Under independent edge\-retention events, the two probabilities strengthen to\(1−ε\)dv\\left\(1\-\\varepsilon\\right\)^\{d\_\{v\}\}and\(1−ε\)\|E\|\\left\(1\-\\varepsilon\\right\)^\{\\left\|E\\right\|\}respectively\.

###### Proof\.

LetAeA\_\{e\}denote the event that true parent relationeeis omitted during preprocessing\. For nodevv, failure to retain alldvd\_\{v\}true parents occurs only if at least one of its true parent relations is omitted\. By the union bound,

P​r​\(P​aG​\(v\)⊈𝒫v\)≤∑e∈P​aG​\(v\)P​r​\(Ae\)≤dv​εPr\\left\(Pa\_\{G\}\\left\(v\\right\)\\nsubseteq\\mathcal\{P\}\_\{v\}\\right\)\\leq\\sum\_\{e\\in Pa\_\{G\}\\left\(v\\right\)\}Pr\\left\(A\_\{e\}\\right\)\\leq d\_\{v\}\\varepsilonwhich gives the bound\.

Under no assumption on dependence among events excluding the true parent relations, the union bound gives the stated lower bounds\. If the edge\-retention events are additionally independent, then

P​r​\(P​aG​\(v\)⊆𝒫v\)=∏e∈P​aG​\(v\)P​r​\(Aec\)≥\(1−ε\)dvPr\\left\(Pa\_\{G\}\\left\(v\\right\)\\subseteq\\mathcal\{P\}\_\{v\}\\right\)=\\prod\_\{e\\in Pa\_\{G\}\\left\(v\\right\)\}Pr\\left\(A\_\{e\}^\{c\}\\right\)\\geq\\left\(1\-\\varepsilon\\right\)^\{d\_\{v\}\}and similarly, probability ofGGremaining feasible be bounded by\(1−ε\)\|E\|\\left\(1\-\\varepsilon\\right\)^\{\\left\|E\\right\|\}\. The independence assumption is not required for the union\-bound result and is stated only to obtain the sharper multiplicative bound\. ∎

Algorithm[1](https://arxiv.org/html/2610.09136#alg1)now summarizes the preprocessing procedure that constructs the admissible parent pools used by the exact optimizer\.

Algorithm 1Directional\-evidence\-guided admissible parent\-set constructionInput:Observational dataset,

X∈𝒳N×DX\\in\\mathcal\{X\}^\{N\\times D\}
Compute the mutual information matrix MI, directional scores

M​C​DMCDusing Equation[1](https://arxiv.org/html/2610.09136#S2.E1), and normalized scores

N​M​C​DNMCDusing Equation[2](https://arxiv.org/html/2610.09136#S2.E2);

Construct the undirected skeleton:

C=\{\{u,v\}:M​I​\(u,v\)≥δM​I\}C=\\\{\\\{u,v\\\}:MI\(u,v\)\\geq\\delta\_\{MI\}\\\};

Locally relax threshold for isolated node after MI thresholding to obtain fallback candidate neighbors;

Initialize set of confirmed directed edges,

Ed​i​r=∅E\_\{dir\}=\\emptysetand set of undirected edges,

Eu​n​d=CE\_\{und\}=C;

For each source node

uu, rank candidate edges

u→vu\\xrightarrow\{\}vby

N​M​C​DNMCD;

for*each candidateu→vu\\xrightarrow\{\}vfrom the ranked order*do

Compute

SSand

N​SNSusing Equation[3](https://arxiv.org/html/2610.09136#S2.E3)and[4](https://arxiv.org/html/2610.09136#S2.E4)respectively;

if*N​Su→v≥δN​S∨Su→v≥δSNS\_\{u\\xrightarrow\{\}v\}\\geq\\delta\_\{NS\}\\lor S\_\{u\\xrightarrow\{\}v\}\\geq\\delta\_\{S\}andu→vu\\xrightarrow\{\}vdoes not create a cycle*then

add

u→vu\\xrightarrow\{\}vto

Ed​i​rE\_\{dir\}and remove

\{u,v\}\\\{u,v\\\}from

Eu​n​dE\_\{und\};

elseKeep

\{u,v\}\\\{u,v\\\}in

Eu​n​dE\_\{und\};

end for

Build the partially directed graph

G′=\(V,Ed​i​r∪Eu​n​d\)G^\{\\prime\}=\\left\(V,E\_\{dir\}\\cup E\_\{und\}\\right\);

Construct an admissible parent set

PvP\_\{v\}for each node

vvsuch that

- •confidently directed edges permit only the identified parent direction
- •ambiguous edges remain feasible in both directions\.

Output:A partially directed graph,

G′G^\{\\prime\}, and admissible parent pool,

𝒫v\\mathcal\{P\}\_\{v\}for every

v∈Vv\\in V

### 2\.3Optimization over restricted parent\-set space

The proposed preprocessing is independent of the optimization framework and can be integrated with any exact parent\-set optimizer\.

In this work, we adopt the classical integer programming formulation for Bayesian network structure learning introduced in[Jaakkola et al\. \(2010\)](https://arxiv.org/html/2610.09136#bib.bib36), further developed via cutting\-plane methods by[Cussens \(2011\)](https://arxiv.org/html/2610.09136#bib.bib8)and[Bartlett and Cussens \(2017\)](https://arxiv.org/html/2610.09136#bib.bib25), as it is the standard benchmark in exact Bayesian network learning\.

Obtained from Algorithm[1](https://arxiv.org/html/2610.09136#alg1),𝒫v\\mathcal\{P\}\_\{v\}is the set of all possible parent nodes of nodevv, and𝒬v\\mathcal\{Q\}\_\{v\}denote the power set of𝒫v\\mathcal\{P\}\_\{v\}\. The empty parent set is always included in𝒬v\\mathcal\{Q\}\_\{v\}, allowing root nodes to be represented within the optimization model\. We define the decision variablezv,𝒬z\_\{v,\\mathcal\{Q\}\}to denote whether𝒬∈𝒬v\\mathcal\{Q\}\\in\\mathcal\{Q\}\_\{v\}is the true set of parent nodes ofvv\(zv,𝒬=1z\_\{v,\\mathcal\{Q\}\}=1\) or not \(zv,𝒬=0z\_\{v,\\mathcal\{Q\}\}=0\)\.

The objective function of IP\-based DAG learning models is defined in terms of decomposable scores, enabling independent evaluation of parent sets\. Since the preprocessing stage already constrains the set of feasible parent configurations, we use the standard BIC score as the local objective without additional penalty terms\. The optimization model selects𝒬\\mathcal\{Q\}for each nodevv\. Formally, we define the score of a node–parent set as

S​c​o​r​e​\(v,𝒬\)=B​I​C​\(v,𝒬\)Score\\left\(v,\\mathcal\{Q\}\\right\)=BIC\\left\(v,\\mathcal\{Q\}\\right\)\(5\)The resulting integer program maximizes the combinatorial score across all nodes and their parent\-configuration search space, subject to the standard convexity and cluster constraints\. This formulation allows the IP to favor structural evidence while permitting corrections when the likelihood strongly supports alternative configurations\. The false\-positive edges introduced during MI\-based skeleton construction can be removed if they yield poor local BN scores, and ambiguous directional decisions can be resolved by global model fit\.

We enforce the constraint∑𝒬∈𝒬vzv,𝒬=1\\sum\_\{\\mathcal\{Q\}\\in\\mathcal\{Q\}\_\{v\}\}z\_\{v,\\mathcal\{Q\}\}=1,∀v\\forall v, to ensure each node has only one subset of parents selected\. A valid DAG should not contain any cycles; however, the initial graph obtained from the IP does not guarantee this requirement\. If the initial graph does contain cycles, we eliminate them by adding cycle\-elimination constraints

∑v∈C∑𝒬∈𝒬v𝒬∩C=∅zv,𝒬≥1\\sum\_\{v\\in C\}\{\\sum\_\{\\begin\{subarray\}\{c\}\\mathcal\{Q\}\\in\\mathcal\{Q\}\_\{v\}\\\\ \\mathcal\{Q\}\\cap C=\\emptyset\\end\{subarray\}\}z\_\{v,\\mathcal\{Q\}\}\}\\geq 1\(6\)and iteratively resolve the IP until a valid DAG is obtained\([Cussens, 2011](https://arxiv.org/html/2610.09136#bib.bib8)\)\. The IP formulation returns one globally optimal DAG within the admissible parent space\. As with other score\-based approaches, multiple DAGs may achieve the same optimal score; in that case, the solver returns a single optimal representative\. The complete MILP formulation over the DECO\-restricted parent sets is provided in Appendix B\.

Figure[1](https://arxiv.org/html/2610.09136#S2.F1)illustrates the complete progression of the proposed DECO framework on the Asia benchmark network described in[Lauritzen and Spiegelhalter \(1988\)](https://arxiv.org/html/2610.09136#bib.bib32)\. The preprocessing stage first constructs an MI\-based dependency skeleton, then estimates directional evidence using the proposed asymmetry measures\. The resulting admissible directional relations are then incorporated into the constrained optimization stage to obtain the final DAG\.

![Refer to caption](https://arxiv.org/html/2610.09136v1/steps.png)Figure 1:DECO pipeline on the Asia network: \(a\) Dependency screening from mutual information, \(b\) Directional evidence extraction and admissible parent\-set construction from the proposed asymmetry measures \(confidence of these directions is implied through edge thickness, and ambiguous edges are dashed in gray\), \(c\) Final DAG obtained after the exact optimization over the restricted search space \(the false positive edge is marked red\)\.Since DECO restricts the admissible parent sets without altering the underlying optimization model, the validity of the cluster constraints and finite termination of the optimization procedure continue to hold\.

## 3Experimental results

In this section, we present the experimental setup and evaluate DECO’s performance at each stage\. Our evaluation follows the three stages of the proposed framework\. We assess the quality of the preprocessing stage by examining \(i\) skeleton construction accuracy and \(ii\) the resulting reduction in admissible parent\-set search space\. Finally, we evaluate the end\-to\-end DAG recovery performance of DECO on benchmark Bayesian networks from Bayesian Network Repository\([Scutari, 2010](https://arxiv.org/html/2610.09136#bib.bib9)\)and synthetic Erdős\-Rényi DAGs\. Erdős\-Rényi is a random directed graph parameterized by the number of ordered nodes and probability they are connected with strictly from lower\-indexed nodes to higher\-indexed nodes to enforce acyclicity\([Erdős and Rényi, 1959](https://arxiv.org/html/2610.09136#bib.bib2)\)\. Additionally, we examine the computational cost of the proposed framework and discuss the trade\-off between search\-space reduction, recovery accuracy, and runtime\.

### 3\.1Skeleton construction

To assess the quality of the preprocessing stage, we first evaluate the dependency\-screening stage by comparing the recovered skeleton against the true graph skeleton\. Skeleton performance is reported using the missing\-edge rate \(omission\) and the extra\-edge rate \(commission\)\. It is important to choose our priority here wisely\. Since the admissible parent sets are constructed directly from the recovered skeleton, a low commission rate will help reduce the search space, but preserving true dependencies at this stage is critical, as omitting true parent relations may exclude the true DAG from the feasible optimization region, as formalized in Proposition 1\. Table[1](https://arxiv.org/html/2610.09136#S3.T1)compares the proposed mutual\-information \(MI\) skeleton construction with Maximum Mean Discrepancy \(MMD\) based dependency screening\([Gretton et al\., 2012](https://arxiv.org/html/2610.09136#bib.bib4)\)and graphical Lasso or GLasso\([Friedman et al\., 2008](https://arxiv.org/html/2610.09136#bib.bib18)\)\.

Table 1:Comparison of skeleton quality measured by omission \(missing\-edge\)=F​NT​P\+F​N=\\frac\{FN\}\{TP\+FN\}and commission \(extra\-edge\)=F​PF​P\+T​N=\\frac\{FP\}\{FP\+TN\}rates\. Lower omission preserves true parent relations, whereas lower commission yields a smaller admissible search space\.Table 1 shows that MI generally provides a favorable balance between omission and commission, with particularly low omission on several benchmark networks while avoiding the very high commission rates observed for MMD\. Since omission errors may eliminate true parent relations and consequently exclude the true DAG from the feasible optimization space \(Proposition 1\), DECO prioritizes preserving true dependencies even at the cost of retaining a small number of additional candidate edges\. Despite this conservative strategy, the preprocessing still substantially reduces the search space, as demonstrated in the following analysis\.

On MI\-skeleton, we examine the trade\-off between preserving true dependencies and reducing the admissible search space with different values ofδM​I\\delta\_\{MI\}\. Figure[2](https://arxiv.org/html/2610.09136#S3.F2)shows the results of Erdős\-Rényi DAGs with fixed size ofn=50n=50and edge probabilityp=0\.05p=0\.05while we vary the number of samples,NN\. From Figure[2](https://arxiv.org/html/2610.09136#S3.F2), one can see that for omission rates up to approximately0\.250\.25, increasing the threshold produces substantial additional search\-space reduction; beyond this region, the marginal gain in reduction becomes small relative to the increase in omission\. Since we want the search space reduction as high as possible without jeopardizing the type II error, we fixδM​I=0\.10\\delta\_\{MI\}=0\.10for this ER\-setting, and it is the same threshold used to produce the results from Table[1](https://arxiv.org/html/2610.09136#S3.T1)\. The qualitative trade\-off remained consistent across the evaluated settings, although the preferred operating region varied with graph size and sample availability; additional results are provided in Appendix C\.1\.

![Refer to caption](https://arxiv.org/html/2610.09136v1/MI_threshold_tradeoff.png)Figure 2:Trade\-off between search\-space reduction and Type II error illustrating the effect of increasing the MI threshold\. Each point corresponds to one threshold value for the ER graph withn=50,p=0\.05n=50,p=0\.05and variableNN\.To further evaluate the effectiveness of the proposed preprocessing stage, we measure the fraction of admissible parent\-set configurations eliminated prior to optimization\. Figure[3](https://arxiv.org/html/2610.09136#S3.F3)reports the average reduction in search space over5050random Erdős\-Rényi DAG realizations of varying graph sizes\. Results are shown separately for discrete and continuous synthetic data\. DECO consistently eliminates a substantial fraction of candidate parent\-set configurations, reducing the search space by approximately50−60%50\-60\\%on discrete datasets and30−55%30\-55\\%on continuous datasets\. These observations empirically support the theoretical results on search\-space reduction\. Together, Table 1, Figure[2](https://arxiv.org/html/2610.09136#S3.F2)and Figure[3](https://arxiv.org/html/2610.09136#S3.F3)demonstrate that the preprocessing stage achieves a favorable balance between preserving true dependencies and reducing the combinatorial search space\.

![Refer to caption](https://arxiv.org/html/2610.09136v1/reduction.png)Figure 3:Average search\-space reduction achieved by DECO on Erdős\-Rényi DAGs of varying size and edge probability=0\.15=0\.15\. The reduction ratio is computed as1−ΩD​E​C​OΩF​u​l​l1\-\\frac\{\\Omega\_\{DECO\}\}\{\\Omega\_\{Full\}\}
### 3\.2Edge orientation

ThresholdsδS\\delta\_\{S\}andδN​S\\delta\_\{NS\}jointly control the degree of confidence required before DECO commits to an edge orientation\. Figure[4](https://arxiv.org/html/2610.09136#S3.F4)examines this trade\-off for the ER\-graph withn=50,p=0\.05n=50,p=0\.05, andN=5,000N=5,000\. We demonstrate the sensitivity in terms of \(a\) correct\-orientation coverage, measured as the fraction of true edges that are correctly oriented; \(b\) ambiguity rate, measured as the fraction of true edges left unoriented; and \(c\) orientation precision, measured as the fraction of committed orientations that are correct\.

![Refer to caption](https://arxiv.org/html/2610.09136v1/heatmap.png)Figure 4:Sensitivity of directional evidence thresholdsδS\\delta\_\{S\}andδN​S\\delta\_\{NS\}\. Increasing the thresholds generally improves the reliability of committed orientations at the expense of directional coverage and increased ambiguity\. The outlined cell indicates the threshold pair used in the main experiments\.Increasing either threshold generally improves the reliability of committed orientations but reduces directional coverage and increases the number of ambiguous edges\. The selected threshold pair lies in an intermediate regime of\(δS,δN​S\)=\(0\.25,0\.5\)\\left\(\\delta\_\{S\},\\delta\_\{NS\}\\right\)=\\left\(0\.25,0\.5\\right\)that achieves high orientation precision while retaining sufficient directional information to meaningfully restrict the subsequent optimization\.

### 3\.3End\-to\-end DAG recovery

For end\-to\-end DAG recovery, we evaluate DECO on both benchmark Bayesian networks and synthetic Erdős\-Rényi DAGs \(see Table 2\)\. For the benchmark networks, observation datasets are generated from the corresponding ground\-truth Bayesian networks using forward sampling\. For the synthetic Erdős\-Rényi directed acyclic graphs \(ER\-DAGs\), we sample a random topological ordering ofnnnodes and add directed edges from earlier to later nodes with probabilityppfor each experiment\. This topological sampling ensures that the generated graph is acyclic\. Then, discrete datasets are generated using conditional probability tables associated with each node, while continuous datasets are generated from linear or non\-linear structural equation models with additive Gaussian noise\.

For each network, we conducted the experiments with sample sizes of1,0001,000,5,0005,000, and10,00010,000, and repeated each experiment over 50 independently generated datasets\. We compare DECO against PC, GES, and Hill\-Climbing \(HC\) on discrete datasets\. For continuous synthetic datasets, we additionally compare against NOTEARS and GOLEM\. We implemented PC and GES usingCausal\-Learn\([Zheng et al\., 2024](https://arxiv.org/html/2610.09136#bib.bib3)\)and HC usingpgmpypackages\([Ankan and Textor, 2024](https://arxiv.org/html/2610.09136#bib.bib20)\)on Python\. GOLEM and NOTEARS were implemented usinggCastlepackage\([Zhang et al\., 2021](https://arxiv.org/html/2610.09136#bib.bib19)\)\. GOBNILP\([Cussens, 2020](https://arxiv.org/html/2610.09136#bib.bib12)\)was run using pygobnilp \[version 1\.0\], with the maximum parent\-set size of 3, and pygobnilp’s default score\-dominance pruning under which a candidate parent set is retained only if its local score exceeds that of all its proper subsets\. GOBNILP’s IP was solved using Gurobi\. This configuration constitutes a standard and non\-directional pruning baseline against which DECO’s performance can be directly compared\. The integer program in DECO was implemented using the PuLP Python library\([Mitchell et al\., 2011](https://arxiv.org/html/2610.09136#bib.bib13)\)and solved with PuLP’s default backend, the open\-source CBC \(COIN\-OR Branch and Cut\) solver\([Forrest and Lougee\-Heimer, 2005](https://arxiv.org/html/2610.09136#bib.bib14)\), with cluster constraints added iteratively as cycles were detected\.

The representative DAGs returned by the implementations were used directly for the metric evaluations\. Performance is assessed using Structural Hamming Distance \(SHD\), False Discovery Rate \(FDR\), True Positive Rate \(TPR\), False Positive Rate \(FPR\), and runtime, and compared against PC, GES, HC, GOBNILP, and continuous\-data baselines where applicable\.

Table 2:Reference networks with their featuresFigure[5](https://arxiv.org/html/2610.09136#S3.F5)compares the SHD of DECO and competing methods on representative datasets\. Since the observed SHD values span over a large range across datasets, the y\-axis is displayed on a logarithmic scale to facilitate visual comparison\. On the benchmark discrete networks \(Figure[5](https://arxiv.org/html/2610.09136#S3.F5)\(a\)\), DECO achieves among the lowest SHD on Asia, Child, Insurance, and Water, while GOBNILP obtains lower SHD on Sachs\. On the larger Hailfinder network, HC achieves lower SHD, whereas DECO remains substantially better than PC and GES\. This behavior suggests that the proposed admissible\-parent construction is highly effective for structure recovery on small and medium\-sized networks, although aggressive pruning may occasionally exclude parent configurations that are beneficial for recovering certain large\-scale structures\.

![Refer to caption](https://arxiv.org/html/2610.09136v1/disc_bar.png)

![Refer to caption](https://arxiv.org/html/2610.09136v1/cont_bar.png)

![Refer to caption](https://arxiv.org/html/2610.09136v1/barlabel.png)

Figure 5:Structural Hamming Distance \(SHD\) comparison of DECO and other competing methods on \(a\) benchmark discrete Bayesian networks and \(b\) nonlinear Gaussian Erdős\-Rényi DAGs, all atN=10,000N=10,000observations\. Error bars denote one standard deviation over repeated runs\. Lower values indicate better structure recovery\. The y\-axis is displayed on a logarithmic scale\.On nonlinear continuous Erdős\-Rényi graphs \(Figure[5](https://arxiv.org/html/2610.09136#S3.F5)\(b\)\), DECO achieves the lowest SHD in two of the three evaluated settings and remains competitive with NOTEARS and GOLEM in the remaining case\. In particular, DECO outperforms all competing methods on the 25\-node \(p=0\.15p=0\.15\) and 50\-node \(p=0\.05p=0\.05\) graphs, while GOLEM achieves lower SHD on the sparser 25\-node \(p=0\.05p=0\.05\) setting\. These results suggest that the proposed directional\-evidence constraints can preserve sufficient flexibility for accurate graph recovery while substantially reducing the search space\. Complete numerical results, including FDR, TPR, and FPR, are provided in Appendix C\.2\. SHD is used as the primary metric in Figure[5](https://arxiv.org/html/2610.09136#S3.F5)because it jointly reflects both missing and spurious edges\.

While DECO achieves favorable SHD across the nonlinear continuous settings, this advantage does not translate into uniform dominance across individual edge\-recovery metrics\. As shown in Appendix C\.2, NOTEARS and GOLEM frequently attain lower FDR and FPR, whereas DECO recovers a larger fraction of true edges in several settings, as reflected by its comparatively higher TPR\. For example, for the nonlinearn=25,p=0\.15,N=10,000n=25,p=0\.15,N=10,000setting, DECO achieves a TPR of0\.680\.68, compared with0\.550\.55for GOLEM and0\.400\.40for NOTEARS, while GOLEM and NOTEARS obtain lower FDR\. This indicates a precision–recall trade\-off: DECO retains greater structural coverage, which can recover additional true relations but may also admit more false\-positive relations\. Consequently, favorable aggregate SHD does not necessarily imply simultaneous superiority in FDR, TPR, and FPR\.

In the linear Gaussian setting, this distinction becomes more pronounced\. Both GOLEM and NOTEARS achieve lower SHD than DECO, consistent with their design for continuous structural equation models under assumptions closely aligned with this data\-generating setting\. Detailed results are provided in Appendix C\.2\.

Across most benchmark and synthetic datasets, increasing the number of observations improves DECO’s recovery performance\. Larger sample sizes yield more reliable estimates of mutual information and directional asymmetry, enabling more accurate admissible\-parent construction and improved optimization outcomes\. Consequently, SHD and FDR generally decrease while TPR tends to improve as the sample size increases\. Although DECO frequently achieves the lowest SHD and FDR, its TPR is not always the highest on larger graphs\. This reflects the conservative nature of the preprocessing stage: by aggressively removing implausible parent relations, DECO reduces false discoveries but may occasionally exclude true dependencies, leading to slightly lower recall\.

Furthermore, we monitored the computational cost of the optimization stage\. On nonlinear Gaussian ER graphs with 50 nodes and 10,000 observations, DECO required approximately146±17146\\pm 17seconds, comparable to HC \(139±206139\\pm 206seconds\) and substantially lower than GOLEM \(1,705±1001,705\\pm 100seconds\)\. GOBNILP achieved the lowest runtime among the exact optimization approaches in this comparison \(10±210\\pm 2seconds\)\. This difference should be interpreted in light of the optimization implementations: GOBNILP employs the highly optimized Gurobi solver together with a fixed maximum parent\-set size of three and built\-in score\-dominance pruning, whereas DECO uses the general\-purpose open\-source CBC solver without imposing a fixed maximum in\-degree\. Thus, the runtime comparison reflects not only the parent\-set reduction strategies but also substantial differences in solver implementation and pruning mechanisms\.

The objective of DECO is therefore not to outperform specialized exact Bayesian\-network solvers in raw optimization time, but to investigate whether data\-driven dependency and directional evidence can reduce the admissible parent\-set space while retaining accurate structural recovery\. In the evaluated setting, DECO remains computationally practical despite using a general\-purpose open\-source MILP solver, while often providing favorable structural recovery relative to the competing methods\. For example, on the 50\-node nonlinear Gaussian ER graphs considered here, DECO achieves lower SHD than GOBNILP despite its higher computational cost\. These results highlight a trade\-off between computational efficiency and structural recovery and suggest that further improvements in the downstream optimization implementation could complement DECO’s preprocessing\-based search\-space reduction\.

## 4Conclusion

In this paper, we propose DECO, a directional evidence\-guided framework for exact DAG learning\. The preprocessing combines dependency screening and distributional directional evidence to construct admissible parent sets before applying a classical exact optimization framework\. Theoretical analysis showed that the proposed preprocessing stage reduces the admissible parent\-set search space exponentially relative to the unconstrained formulation\. Experimental results on benchmark Bayesian networks and simulated DAGs demonstrated that DECO substantially reduces the search space while achieving competitive or superior structure recovery performance across a variety of discrete and nonlinear continuous settings\. Detailed edge\-recovery metrics further reveal a precision–recall trade\-off, particularly for continuous data, where greater true\-edge coverage can be accompanied by higher false\-discovery rates than more conservative competing methods\.

One of the main limitations of the current framework is that the directional evidence is derived from pairwise distributional asymmetry and may not fully capture higher\-order causal interactions among variables\. Consequently, directional decisions currently rely on confidence thresholds that may discard useful information in more complex dependency structures\. Second, although the proposed preprocessing significantly reduces the search space, the exact optimization stage can remain computationally demanding for larger graphs\. In addition, extending the framework to mixed\-type variables and more general structural equation models remains an important direction for future research\.

## Statements and Declarations

### Data availability

The benchmark Bayesian network structures and parameters used in this study \(Asia, Sachs, Child, Insurance, Water, Barley, and Hailfinder\) were obtained from the bnlearn Bayesian Network Repository, maintained by Marco Scutari, and are publicly available at https://www\.bnlearn\.com/bnrepository/\. The observational datasets used in the benchmark experiments were generated from these networks by forward sampling as described in the manuscript\. The synthetic Erdős–Rényi DAGs and corresponding discrete and continuous datasets were generated using the procedures described in Section 3\.3\. No restricted\-access or proprietary datasets were used\.

## References

- A\. Ankan and J\. TextorPgmpy: a python toolkit for bayesian networks\.Journal of Machine Learning Research25\(265\),pp\. 1–8\.External Links:[Link](https://www.jmlr.org/papers/v25/23-0487.html)Cited by:[§3\.3](https://arxiv.org/html/2610.09136#S3.SS3.p2.1)\.
- Bartlett and Cussens \(2017\)M\. Bartlett and J\. CussensInteger Linear Programming for the Bayesian network structure learning problem\.Artificial Intelligence244,pp\. 258–271\.Cited by:[Appendix B](https://arxiv.org/html/2610.09136#A2.p1.1),[§1](https://arxiv.org/html/2610.09136#S1.p2.1),[§2\.3](https://arxiv.org/html/2610.09136#S2.SS3.p2.1)\.
- Bouckaert \(1995\)R\. R\. BouckaertBayesian belief networks: from construction to inference\.Ph\.D\. Thesis,Utrecht University\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p2.1)\.
- Chow and Liu \(1968\)C\. Chow and C\. LiuApproximating discrete probability distributions with dependence trees\.IEEE transactions on Information Theory14\(3\),pp\. 462–467\.Cited by:[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p2.1)\.
- Coveret al\.\(1991\)T\. M\. Cover J\. A\. Thomaset al\.Entropy, relative entropy and mutual information\.Elements of information theory2\(1\),pp\. 12–13\.Cited by:[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p2.1)\.
- Cowls and Schroeder \(2015\)J\. Cowls and R\. SchroederCausation, correlation, and big data in social science research\.Policy & Internet7\(4\),pp\. 447–472\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p1.1)\.
- Cunningham \(2021\)S\. CunninghamCausal inference: the mixtape\.Yale University Press,New Haven\.External Links:ISBN 9780300251685,[Link](https://scunning.com/)Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p1.1)\.
- Cussens \(2011\)J\. CussensBayesian network learning with cutting planes\.UAI 2011, Proceedings of the 27th Conference on Uncertainty in Artificial Intelligence\.External Links:[Link](https://arxiv.org/pdf/2304.02146)Cited by:[Appendix B](https://arxiv.org/html/2610.09136#A2.p1.1),[Appendix B](https://arxiv.org/html/2610.09136#A2.p6.1),[§2\.3](https://arxiv.org/html/2610.09136#S2.SS3.p2.1),[§2\.3](https://arxiv.org/html/2610.09136#S2.SS3.p5.2)\.
- Cussens \(2020\)J\. CussensGOBNILP: learning Bayesian network structure with integer programming\.InInternational Conference on Probabilistic Graphical Models,pp\. 605–608\.Cited by:[§3\.3](https://arxiv.org/html/2610.09136#S3.SS3.p2.1)\.
- Duong and Nguyen \(2022\)B\. Duong and T\. NguyenBivariate causal discovery via conditional divergence\.InConference on Causal Learning and Reasoning,pp\. 236–252\.Cited by:[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p3.1),[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p6.1)\.
- Erdős and Rényi \(1959\)P\. Erdős and A\. RényiOn random graphs i\.Publ\. math\. debrecen6\(290\-297\),pp\. 18\.Cited by:[§3](https://arxiv.org/html/2610.09136#S3.p1.1)\.
- Forrest and Lougee\-Heimer \(2005\)J\. Forrest and R\. Lougee\-HeimerCBC user guide\.InEmerging Theory, Methods, and Applications,pp\. 257–277\.External Links:[Document](https://dx.doi.org/10.1287/educ.1053.0020),[Link](https://doi.org/10.1287/educ.1053.0020)Cited by:[§3\.3](https://arxiv.org/html/2610.09136#S3.SS3.p2.1)\.
- Friedmanet al\.\(2008\)J\. Friedman, T\. Hastie, and R\. TibshiraniSparse inverse covariance estimation with the graphical lasso\.Biostatistics9\(3\),pp\. 432–441\.Cited by:[§3\.1](https://arxiv.org/html/2610.09136#S3.SS1.p1.1)\.
- Gámezet al\.\(2011\)J\. A\. Gámez, J\. L\. Mateo, and J\. M\. PuertaLearning Bayesian networks by hill climbing: efficient methods based on progressive restriction of the neighborhood\.Data Mining and Knowledge Discovery22,pp\. 106–148\.External Links:[Document](https://dx.doi.org/10.1007/s10618-010-0178-6)Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p2.1)\.
- Gasseet al\.\(2014\)M\. Gasse, A\. Aussem, and H\. ElghazelA hybrid algorithm for bayesian network structure learning with application to multi\-label learning\.Expert Systems with Applications41\(15\),pp\. 6755–6772\.Cited by:[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p1.1)\.
- Gasseet al\.\(2015\)M\. Gasse, A\. Aussem, and H\. ElghazelOn the optimality of multi\-label classification under subset zero\-one loss for distributions satisfying the composition property\.InInternational Conference on Machine Learning,pp\. 2531–2539\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p2.1)\.
- Grettonet al\.\(2012\)A\. Gretton, K\. M\. Borgwardt, M\. J\. Rasch, B\. Schölkopf, and A\. SmolaA kernel two\-sample test\.The journal of machine learning research13\(1\),pp\. 723–773\.Cited by:[§3\.1](https://arxiv.org/html/2610.09136#S3.SS1.p1.1)\.
- Henningsenet al\.\(2025\)A\. Henningsen, G\. Low, D\. Wuepper, T\. Dalhaus, H\. Storm, D\. Belay, and S\. HirschEstimating causal effects with observational data: guidelines for agricultural and applied economists\.Journal of Agricultural Economics\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p1.1)\.
- Jaakkolaet al\.\(2010\)T\. Jaakkola, D\. Sontag, A\. Globerson, and M\. MeilaLearning Bayesian network structure using LP relaxations\.InProceedings of the thirteenth international conference on artificial intelligence and statistics,pp\. 358–365\.Cited by:[Appendix B](https://arxiv.org/html/2610.09136#A2.p1.1),[§2\.3](https://arxiv.org/html/2610.09136#S2.SS3.p2.1)\.
- Lauritzen and Spiegelhalter \(1988\)S\. L\. Lauritzen and D\. J\. SpiegelhalterLocal computations with probabilities on graphical structures and their application to expert systems\.Journal of the Royal Statistical Society: Series B \(Methodological\)50\(2\),pp\. 157–194\.Cited by:[§2\.3](https://arxiv.org/html/2610.09136#S2.SS3.p6.1)\.
- Manzouret al\.\(2021\)H\. Manzour, S\. Küçükyavuz, H\. Wu, and A\. ShojaieInteger programming for learning directed acyclic graphs from continuous data\.INFORMS journal on optimization3\(1\),pp\. 46–73\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p2.1)\.
- Mitchellet al\.\(2011\)S\. Mitchell, M\. O’Sullivan, and I\. DunningPuLP: a linear programming toolkit for python\.External Links:[Link](https://optimization-online.org/)Cited by:[§3\.3](https://arxiv.org/html/2610.09136#S3.SS3.p2.1)\.
- Pearl \(2009\)J\. PearlCausality\.Cambridge university press\.Cited by:[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p2.1)\.
- Pingaultet al\.\(2018\)J\. Pingault, P\. F\. O’reilly, T\. Schoeler, G\. B\. Ploubidis, F\. Rijsdijk, and F\. DudbridgeUsing genetic data to strengthen causal inference in observational research\.Nature Reviews Genetics19\(9\),pp\. 566–580\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p1.1)\.
- Rollandet al\.\(2022\)P\. Rolland, V\. Cevher, M\. Kleindessner, C\. Russell, D\. Janzing, B\. Schölkopf, and F\. LocatelloScore matching enables causal discovery of nonlinear additive noise models\.InInternational Conference on Machine Learning,pp\. 18741–18753\.Cited by:[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p1.1)\.
- Sanchezet al\.\(2022\)P\. Sanchez, X\. Liu, A\. Q\. O’Neil, and S\. A\. TsaftarisDiffusion models for causal discovery via topological ordering\.arXiv preprint arXiv:2210\.06201\.Cited by:[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p1.1)\.
- Scutariet al\.\(2018\)M\. Scutari, C\. E\. Graafland, and J\. M\. GutiérrezWho learns better bayesian network structures: constraint\-based, score\-based or hybrid algorithms?\.InInternational conference on probabilistic graphical models,pp\. 416–427\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p2.1)\.
- Scutari \(2010\)M\. ScutariLearning bayesian networks with the bnlearn r package\.Journal of Statistical Software35\(3\),pp\. 1–22\.External Links:[Link](https://www.jstatsoft.org/article/view/v035i03/394)Cited by:[§3](https://arxiv.org/html/2610.09136#S3.p1.1)\.
- Sethet al\.\(2015\)A\. K\. Seth, A\. B\. Barrett, and L\. BarnettGranger causality analysis in neuroscience and neuroimaging\.Journal of Neuroscience35\(8\),pp\. 3293–3297\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p1.1)\.
- Spirteset al\.\(2000\)P\. Spirtes, C\. N\. Glymour, and R\. ScheinesCausation, prediction, and search\.MIT press\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p2.1),[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p2.1)\.
- Spirtes \(2010\)P\. SpirtesIntroduction to causal inference\.\.Journal of Machine Learning Research11\(5\)\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p2.1)\.
- Tsamardinoset al\.\(2006\)I\. Tsamardinos, L\. E\. Brown, and C\. F\. AliferisThe Max\-Min Hill\-Climbing Bayesian network structure learning algorithm\.Machine learning65\(1\),pp\. 31–78\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p2.1),[§2\.1](https://arxiv.org/html/2610.09136#S2.SS1.p1.1)\.
- Williamson and Forbes \(2021\)E\. J\. Williamson and A\. ForbesA tutorial on directed acyclic graphs\.Journal of Clinical Epidemiology139,pp\. 257–265\.External Links:[Document](https://dx.doi.org/10.1016/j.jclinepi.2021.08.012),[Link](https://pmc.ncbi.nlm.nih.gov/articles/PMC8821727/)Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p1.1)\.
- Zhanget al\.\(2017\)H\. Zhang, S\. Zhou, K\. Zhang, and J\. GuanCausal discovery using regression\-based conditional independence tests\.InProceedings of the AAAI conference on artificial intelligence,Vol\.31\.Cited by:[§1](https://arxiv.org/html/2610.09136#S1.p2.1)\.
- Zhanget al\.\(2021\)K\. Zhang, S\. Zhu, M\. Kalander, I\. Ng, J\. Ye, Z\. Chen, and L\. PanGCastle: a python toolbox for causal discovery\.arXiv preprint arXiv:2111\.15155\.Cited by:[§3\.3](https://arxiv.org/html/2610.09136#S3.SS3.p2.1)\.
- Zhenget al\.\(2024\)Y\. Zheng, B\. Huang, W\. Chen, J\. Ramsey, M\. Gong, R\. Cai, S\. Shimizu, P\. Spirtes, and K\. ZhangCausal\-learn: causal discovery in python\.Journal of Machine Learning Research25\(60\),pp\. 1–8\.Cited by:[§3\.3](https://arxiv.org/html/2610.09136#S3.SS3.p2.1)\.

## Appendix AFull proof of Theorem 1

This appendix provides the complete proof of Theorem 1, which characterizes the reduction in admissible parent\-set configurations induced by DECO preprocessing\.

###### Proof\.

In the unconstrained model, each nodevvcan select any subset of the remainingD−1D\-1variables\. By the binomial theorem, the number of candidate configurations for a single node is∑k=0D−1\(D−1k\)=2D−1\\sum\_\{k=0\}^\{D\-1\}\\binom\{D\-1\}\{k\}=2^\{D\-1\}\. Across allDDnodes, the unrestricted parent\-set configuration space is therefore:

Ωfull=\(2D−1\)D=2D⁡\(D−1\)\\Omega\_\{\\text\{full\}\}=\\left\(2^\{D\-1\}\\right\)^\{D\}=2^\{D\(D\-1\)\}In the DECO framework, the candidate parent pool is restricted to𝒫v\\mathcal\{P\}\_\{v\}, containingmvm\_\{v\}allowable parents\. The number of allowable parent sets in𝒬v\\mathcal\{Q\}\_\{v\}becomes∑k=0mv\(mvk\)=2mv\\sum\_\{k=0\}^\{m\_\{v\}\}\\binom\{m\_\{v\}\}\{k\}=2^\{m\_\{v\}\}\. Taking the Cartesian product across all nodes yields an admissible parent\-set configuration space prior to global acyclicity enforcement:

ΩDECO=∏v∈V2mv=2∑v∈Vmv\\Omega\_\{\\text\{DECO\}\}=\\prod\_\{v\\in V\}2^\{m\_\{v\}\}=2^\{\\sum\_\{v\\in V\}m\_\{v\}\}Under the assumptionm¯≤\(1−α\)​\(D−1\)\\bar\{m\}\\leq\\left\(1\-\\alpha\\right\)\\left\(D\-1\\right\), equivalently∑v∈Vmv≤D⁡\(1−α\)​\(D−1\)\\sum\_\{v\\in V\}m\_\{v\}\\leq D\(1\-\\alpha\)\(D\-1\), the ratio of the two configuration spaces satisfies

ΩDECOΩfull=2∑v∈Vmv2D⁡\(D−1\)≤2D​\(1−α\)​\(D−1\)2D⁡\(D−1\)=2−α​D​\(D−1\)\\frac\{\\Omega\_\{\\text\{DECO\}\}\}\{\\Omega\_\{\\text\{full\}\}\}=\\frac\{2^\{\\sum\_\{v\\in V\}m\_\{v\}\}\}\{2^\{D\(D\-1\)\}\}\\leq\\frac\{2^\{D\(1\-\\alpha\)\(D\-1\)\}\}\{2^\{D\(D\-1\)\}\}=2^\{\-\\alpha D\(D\-1\)\}Thus, a linear reduction in the average admissible parent\-pool size yields an exponential reduction in the parent\-set configuration space considered prior to global acyclicity enforcement\. ∎

## Appendix BComplete MILP formulation

For completeness, we provide the integer programming formulation used for downstream exact optimization\. This formulation follows the standard parent\-set\-based Bayesian network structure\-learning framework of[Jaakkola et al\. \(2010\)](https://arxiv.org/html/2610.09136#bib.bib36),[Cussens \(2011\)](https://arxiv.org/html/2610.09136#bib.bib8), and[Bartlett and Cussens \(2017\)](https://arxiv.org/html/2610.09136#bib.bib25); DECO modifies the admissible parent\-set collections supplied to the optimizer\.

For eachv∈Vv\\in Vlet𝒫v\\mathcal\{P\}\_\{v\}and𝒬v\\mathcal\{Q\}\_\{v\}denote the admissible parent pool constructed by DECO and its power set respectively\. For eachQ∈𝒬vQ\\in\\mathcal\{Q\}\_\{v\}, define binary variables:

zv,Q=\{1if​Q​is selected as the parent set of​v0otherwisez\_\{v,Q\}=\\begin\{cases\}1&\\text\{if \}Q\\text\{ is selected as the parent set of \}v\\\\ 0&\\text\{otherwise\}\\end\{cases\}
The restricted parent\-set MILP is,

maxz∑v∈V∑Q∈𝒬vScore⁡\(v,Q\)​zv,Q\\max\_\{z\}\\quad\\sum\_\{v\\in V\}\\sum\_\{Q\\in\\mathcal\{Q\}\_\{v\}\}\\mathrm\{Score\}\(v,Q\)z\_\{v,Q\}subject to

∑Q∈𝒬vzv,Q=1,\\displaystyle\\sum\_\{Q\\in\\mathcal\{Q\}\_\{v\}\}z\_\{v,Q\}=1,∀v∈V\\displaystyle\{\}\\quad\\forall v\\in V∑v∈C∑𝒬∈𝒬v𝒬∩C=∅zv,𝒬≥1,\\displaystyle\\sum\_\{v\\in C\}\{\\sum\_\{\\begin\{subarray\}\{c\}\\mathcal\{Q\}\\in\\mathcal\{Q\}\_\{v\}\\\\ \\mathcal\{Q\}\\cap C=\\emptyset\\end\{subarray\}\}\}z\_\{v,\\mathcal\{Q\}\}\\geq 1,∀C⊆V,\|C\|≥2\\displaystyle\\quad\\forall C\\subseteq V,\\left\|C\\right\|\\geq 2zv,𝒬∈\{0,1\},\\displaystyle z\_\{v,\\mathcal\{Q\}\}\\in\\\{0,1\\\},∀v∈V,𝒬∈𝒬v\.\\displaystyle\\quad\\forall v\\in V,\\mathcal\{Q\}\\in\\mathcal\{Q\}\_\{v\}\.
The first constraint selects exactly one parent set for each node, while the cluster constraints enforce acyclicity\. The violated cluster constraints are added iteratively following the cutting\-plane procedure of[Cussens \(2011\)](https://arxiv.org/html/2610.09136#bib.bib8)\.

## Appendix CAdditional experimental results

### C\.1Threshold sensitivity analyses

The main text \(Figure[2](https://arxiv.org/html/2610.09136#S3.F2)\) illustrates the effect of the MI threshold,δM​I\\delta\_\{MI\}, for a representative ER\-DAG setting\. Here, we extend the sensitivity analysis across different graph sizesnnand sample\-size\-to\-dimension ratiosN/nN/nto assess whether the observed trade\-off between skeleton preservation and search\-space reduction persists across settings\. Figure 6 summarizes the results for discrete ER\-DAGs withp=0\.05p=0\.05\.

![Refer to caption](https://arxiv.org/html/2610.09136v1/MI_threshold_p05.png)Figure 6:Sensitivity of MI\-based dependency screening toδM​I\\delta\_\{MI\}across graph sizesnnand sample size to graph size ratioN/nN/nfor discrete Erdős–Rényi DAGs withp=0\.05p=0\.05\. IncreasingδM​I\\delta\_\{MI\}generally decreases commission and increases search\-space reduction, but eventually increases omission, illustrating the trade\-off between structural preservation and pruning\.
### C\.2DAG recovery results

We compared the SHD of different DAG learning methods, including DECO, in the main text \(Figure[5](https://arxiv.org/html/2610.09136#S3.F5)\)\. Here, we provide the complete numerical results, including Structural Hamming Distance \(SHD\), False Discovery Rate \(FDR\), True Positive Rate \(TPR\), and False Positive Rate \(FPR\), across the evaluated discrete and continuous networks and sample sizes\.

Table 3:Comparison of DAG learning performance on benchmark discrete networksData sizeMethodSHD↓\\downarrowFDR↓\\downarrowTPR↑\\uparrowFPR↓\\downarrowNetwork: Asia \(Nodes: 8, Arcs: 8\)1000PC11\.36±\\pm0\.9750\.88±\\pm0\.0650\.12±\\pm0\.0700\.14±\\pm0\.033GES12\.58±\\pm1\.0020\.89±\\pm0\.0510\.12±\\pm0\.0470\.17±\\pm0\.028HC4\.04±\\pm2\.1910\.35±\\pm0\.2260\.59±\\pm0\.1800\.06±\\pm0\.040GOBNILP2\.82±\\pm1\.5060\.24±\\pm0\.1340\.74±\\pm0\.0980\.04±\\pm0\.026DECO1\.00±\\pm0\.2830\.01±\\pm0\.0300\.88±\\pm0\.0350\.00±\\pm0\.0055000PC12\.24±\\pm0\.7090\.87±\\pm0\.0320\.17±\\pm0\.0600\.18±\\pm0\.026GES12\.54±\\pm0\.6990\.90±\\pm0\.0190\.13±\\pm0\.0250\.19±\\pm0\.028HC6\.64±\\pm2\.1330\.64±\\pm0\.2080\.37±\\pm0\.1900\.12±\\pm0\.043GOBNILP2\.14±\\pm1\.3270\.19±\\pm0\.1330\.78±\\pm0\.1120\.03±\\pm0\.026DECO0\.88±\\pm0\.3820\.01±\\pm0\.0340\.89±\\pm0\.0430\.00±\\pm0\.00610000PC12\.86±\\pm0\.7210\.88±\\pm0\.0310\.15±\\pm0\.0530\.19±\\pm0\.023GES12\.72±\\pm0\.750\.90±\\pm0\.030\.13±\\pm0\.030\.20±\\pm0\.02HC7\.02±\\pm2\.8670\.60±\\pm0\.2320\.44±\\pm0\.2000\.12±\\pm0\.053GOBNILP1\.34±\\pm0\.7100\.14±\\pm0\.0680\.84±\\pm0\.0780\.02±\\pm0\.012DECO0\.58±\\pm0\.5330\.02±\\pm0\.0400\.93±\\pm0\.0620\.00±\\pm0\.007Network: Sachs \(Nodes: 11, Arcs: 17\)1000PC23\.9±\\pm1\.7230\.81±\\pm0\.0600\.19±\\pm0\.0700\.15±\\pm0\.022GES23\.84±\\pm1\.3320\.88±\\pm0\.0460\.10±\\pm0\.0410\.12±\\pm0\.017HC9\.94±\\pm2\.7160\.46±\\pm0\.1510\.45±\\pm0\.1410\.07±\\pm0\.022GOBNILP7\.00±\\pm1\.7090\.37±\\pm0\.1010\.59±\\pm0\.1010\.06±\\pm0\.017DECO9\.62±\\pm2\.2440\.49±\\pm0\.1330\.44±\\pm0\.1220\.08±\\pm\.0195000PC26\.14±\\pm1\.2500\.84±\\pm0\.0500\.23±\\pm0\.0780\.21±\\pm0\.020GES26\.84±\\pm1\.2550\.88±\\pm0\.0460\.13±\\pm0\.0520\.17±\\pm0\.0171HC8\.44±\\pm1\.8460\.47±\\pm0\.0940\.55±\\pm0\.0970\.09±\\pm0\.018GOBNILP5\.02±\\pm1\.8810\.3±\\pm0\.1110\.7±\\pm0\.1110\.05±\\pm0\.020DECO7\.68±\\pm1\.6180\.42±\\pm0\.1010\.55±\\pm0\.0950\.07±\\pm0\.01710000PC27\.38±\\pm1\.3840\.83±\\pm0\.0480\.25±\\pm0\.0770\.22±\\pm0\.016GES27\.62±\\pm0\.9780\.85±\\pm0\.0350\.18±\\pm0\.04890\.19±\\pm0\.014HC9\.24±\\pm2\.1960\.51±\\pm0\.1200\.51±\\pm0\.1160\.10±\\pm0\.024GOBNILP5\.08±\\pm1\.9780\.30±\\pm0\.1160\.70±\\pm0\.1160\.05±\\pm0\.021DECO7\.42±\\pm1\.9610\.40±\\pm0\.1230\.56±\\pm0\.1150\.07±\\pm0\.021Network: Child \(Nodes: 20, Arcs: 25\)1000PC43\.38±\\pm1\.7420\.89±\\pm0\.0260\.14±\\pm0\.0320\.08±\\pm0\.006GES48\.10±\\pm2\.7590\.90±\\pm0\.0340\.15±\\pm0\.0550\.09±\\pm0\.010HC14\.64±\\pm2\.6960\.42±\\pm0\.0880\.55±\\pm0\.0790\.03±\\pm0\.007GOBNILP17\.04±\\pm1\.7770\.56±\\pm0\.0700\.44±\\pm0\.0670\.11±\\pm0\.004DECO11\.30±\\pm1\.6160\.16±\\pm0\.0680\.57±\\pm0\.0580\.01±\\pm0\.0045000PC50\.74±\\pm1\.9680\.91±\\pm0\.0300\.14±\\pm0\.0500\.10±\\pm0\.007GES60\.12±\\pm2\.7900\.89±\\pm0\.0230\.24±\\pm0\.0570\.14±\\pm0\.009HC9\.78±\\pm3\.8750\.33±\\pm0\.1180\.69±\\pm0\.0970\.03±\\pm0\.010GOBNILP12\.66±\\pm1\.0120\.31±\\pm0\.0400\.59±\\pm0\.0390\.11±\\pm0\.003DECO9\.68±\\pm1\.1740\.21±\\pm0\.0350\.66±\\pm0\.0380\.01±\\pm0\.00210000PC53\.72±\\pm1\.9700\.91±\\pm0\.0270\.16±\\pm0\.0500\.11±\\pm0\.010GES65\.88±\\pm2\.2680\.90±\\pm0\.0220\.24±\\pm0\.0600\.15±\\pm0\.010HC10\.86±\\pm2\.4080\.38±\\pm0\.0690\.68±\\pm0\.0580\.03±\\pm0\.006GOBNILP12\.48±\\pm0\.9220\.40±\\pm0\.0370\.50±\\pm0\.0370\.11±\\pm0\.003DECO9\.3±\\pm1\.1360\.24±\\pm0\.0320\.67±\\pm0\.0450\.02±\\pm0\.002Network: Insurance \(Nodes: 27, Arcs: 52\)10000PC101\.34±\\pm2\.9090\.86±\\pm0\.0190\.21±\\pm0\.0280\.11±\\pm0\.005GES129\.00±\\pm3\.5780\.93±\\pm0\.0200\.13±\\pm0\.0410\.15±\\pm0\.007HC28\.96±\\pm5\.3520\.36±\\pm0\.0830\.57±\\pm0\.0650\.03±\\pm0\.007GOBNILP39\.74±\\pm1\.1800\.46±\\pm0\.0350\.27±\\pm0\.0170\.02±\\pm0\.001DECO22\.50±\\pm0\.6710\.21±\\pm0\.0230\.61±\\pm0\.0130\.01±\\pm0\.002Network: Water \(Nodes: 32, Arcs: 66\)10000PC100\.58±\\pm3\.8640\.91±\\pm0\.0120\.02±\\pm0\.0340\.09±\\pm0\.004GES102\.22±\\pm2\.2500\.96±\\pm0\.0170\.03±\\pm0\.0120\.05±\\pm0\.003HC58\.70±\\pm4\.8470\.58±\\pm0\.1020\.22±\\pm0\.0510\.02±\\pm0\.004GOBNILP59\.18±\\pm2\.7550\.42±\\pm0\.0640\.24±\\pm0\.0320\.02±\\pm0\.003DECO53\.66±\\pm3\.1980\.46±\\pm0\.0420\.38±\\pm0\.0440\.02±\\pm0\.004Network: Hailfinder \(Nodes: 56, Arcs: 66\)1000PC98\.34±\\pm2\.0550\.97±\\pm0\.0130\.02±\\pm0\.0090\.02±\\pm0\.001GES162\.14±\\pm6\.6420\.96±\\pm0\.0150\.07±\\pm0\.0240\.04±\\pm0\.002HC13\.06±\\pm3\.2950\.37±\\pm0\.1070\.59±\\pm0\.0870\.03±\\pm0\.008GOBNILP55\.86±\\pm3\.3710\.62±\\pm0\.0390\.31±\\pm0\.0320\.01±\\pm0\.001DECO39\.95±\\pm4\.7530\.35±\\pm0\.0600\.52±\\pm0\.0450\.01±\\pm0\.0015000PC102\.80±\\pm1\.8760\.97±\\pm0\.0120\.02±\\pm0\.0090\.02±\\pm0\.001GES188\.63±\\pm6\.3040\.95±\\pm0\.0180\.10±\\pm0\.0370\.05±\\pm0\.002HC10\.24±\\pm3\.4090\.35±\\pm0\.1050\.70±\\pm0\.0830\.03±\\pm0\.010GOBNILP53\.68±\\pm4\.0170\.60±\\pm0\.0430\.34±\\pm0\.0360\.01±\\pm0\.001DECO38\.21±\\pm4\.2700\.30±\\pm0\.0440\.57±\\pm0\.0370\.01±\\pm0\.00110000PC103\.08±\\pm2\.0480\.97±\\pm0\.0090\.02±\\pm0\.0070\.02±\\pm0\.001GES204\.6±\\pm5\.030\.77±\\pm0\.4250\.12±\\pm0\.0260\.05±\\pm0\.000HC10\.48±\\pm2\.9880\.37±\\pm0\.0850\.70±\\pm0\.0690\.03±\\pm0\.008GOBNILP53\.28±\\pm4\.3360\.60±\\pm0\.0440\.33±\\pm0\.0360\.01±\\pm0\.001DECO37\.10±\\pm6\.0680\.28±\\pm0\.0630\.63±\\pm0\.0950\.00±\\pm0\.001Table 4:Comparison of DAG learning performance on Erdős\-Rényi Graphs and discrete dataData sizeMethodSHD↓\\downarrowFDR↓\\downarrowTPR↑\\uparrowFPR↓\\downarrowNodes: 25, Edge probability: 0\.151000PC19\.38±\\pm5\.7790\.3±\\pm0\.0850\.74±\\pm0\.1140\.03±\\pm0\.007GES16\.28±\\pm5\.9230\.27±\\pm0\.0870\.82±\\pm0\.0880\.02±\\pm0\.009HC19\.82±\\pm6\.0690\.31±\\pm0\.0990\.63±\\pm0\.0920\.02±\\pm0\.008GOBNILP31\.80±\\pm5\.4370\.43±\\pm0\.0830\.30±\\pm0\.0540\.02±\\pm0\.003DECO13\.96±\\pm5\.2690\.23±\\pm0\.0940\.69±\\pm0\.0940\.02±\\pm0\.0075000PC15\.08±\\pm6\.3470\.25±\\pm0\.0700\.89±\\pm0\.0870\.02±\\pm0\.009GES16\.30±\\pm11\.9880\.28±\\pm0\.1200\.91±\\pm0\.0930\.03±\\pm0\.020HC16\.76±\\pm7\.6630\.3±\\pm0\.1210\.74±\\pm0\.0960\.03±\\pm0\.013GOBNILP32\.24±\\pm5\.9650\.41±\\pm0\.0820\.30±\\pm0\.0510\.02±\\pm0\.003DECO10\.54±\\pm4\.0500\.20±\\pm0\.0950\.78±\\pm0\.0770\.01±\\pm0\.00810000PC15\.97±\\pm5\.0860\.25±\\pm0\.0550\.9±\\pm0\.0660\.03±\\pm0\.008GES18\.45±\\pm10\.4370\.31±\\pm0\.0950\.92±\\pm0\.0700\.04±\\pm0\.017HC13\.16±\\pm6\.5900\.23±\\pm0\.1070\.79±\\pm0\.0770\.02±\\pm0\.011GOBNILP34\.24±\\pm5\.2520\.43±\\pm0\.1020\.28±\\pm0\.0520\.02±\\pm0\.004DECO12\.22±\\pm2\.3930\.23±\\pm0\.0800\.75±\\pm0\.0390\.02±\\pm0\.007Nodes: 25, Edge probability: 0\.251000PC53\.11±\\pm8\.1620\.41±\\pm0\.0740\.45±\\pm0\.0620\.05±\\pm0\.010GES53\.04±\\pm11\.7660\.41±\\pm0\.1110\.54±\\pm0\.0980\.06±\\pm0\.019HC47\.26±\\pm10\.4580\.32±\\pm0\.1070\.42±\\pm0\.0090\.03±\\pm0\.001GOBNILP59\.98±\\pm10\.0290\.38±\\pm0\.0730\.20±\\pm0\.0380\.02±\\pm0\.008DECO43\.78±\\pm6\.9270\.24±\\pm0\.0820\.43±\\pm0\.0610\.02±\\pm0\.0075000PC44\.40±\\pm8\.5420\.36±\\pm0\.0680\.62±\\pm0\.0830\.05±\\pm0\.012GES51\.18±\\pm12\.9120\.41±\\pm0\.0780\.71±\\pm0\.0770\.07±\\pm0\.022HC35\.46±\\pm7\.9810\.26±\\pm0\.0940\.58±\\pm0\.0740\.03±\\pm0\.012GOBNILP58\.62±\\pm7\.3890\.35±\\pm0\.1000\.21±\\pm0\.0370\.02±\\pm0\.005DECO43\.47±\\pm6\.3650\.26±\\pm0\.0730\.44±\\pm0\.0520\.02±\\pm0\.00610000PC44\.16±\\pm9\.5010\.36±\\pm0\.0570\.65±\\pm0\.0780\.05±\\pm0\.011GES59\.51±\\pm15\.6630\.46±\\pm0\.0930\.73±\\pm0\.0900\.09±\\pm0\.029HC32\.65±\\pm9\.9670\.25±\\pm0\.1110\.64±\\pm0\.0870\.03±\\pm0\.015GOBNILP58\.58±\\pm7\.7870\.36±\\pm0\.0930\.21±\\pm0\.0330\.02±\\pm0\.004DECO40\.5±\\pm6\.640\.27±\\pm0\.0880\.48±\\pm0\.0880\.02±\\pm0\.007Nodes: 50, Edge probability: 0\.0510000PC28\.70±\\pm5\.2480\.35±\\pm0\.0700\.90±\\pm0\.0590\.01±\\pm0\.002GES15\.73±\\pm4\.3990\.27±\\pm0\.0580\.97±\\pm0\.0410\.01±\\pm0\.002HC17\.23±\\pm5\.9760\.26±\\pm0\.0770\.81±\\pm0\.0580\.01±\\pm0\.003GOBNILP17\.76±\\pm4\.4880\.53±\\pm0\.0800\.48±\\pm0\.0800\.01±\\pm0\.002DECO9\.80±\\pm2\.2800\.14±\\pm0\.0300\.86±\\pm0\.0410\.01±\\pm0\.001Table 5:Comparison of DAG learning performance on Erdős\-Rényi Graphs and continuous \(Non\-linear Gaussian\) dataData sizeMethodSHD↓\\downarrowFDR↓\\downarrowTPR↑\\uparrowFPR↓\\downarrowNodes: 25, Edge probability: 0\.151000PC21\.53±\\pm6\.4920\.34±\\pm0\.0900\.71±\\pm0\.1000\.03±\\pm0\.009GES21\.90±\\pm10\.2390\.34±\\pm0\.1100\.81±\\pm0\.0990\.04±\\pm0\.015HC29\.10±\\pm6\.5190\.46±\\pm0\.0800\.70±\\pm0\.0810\.05±\\pm0\.008NOTEARS25\.27±\\pm3\.8440\.00±\\pm0\.0130\.46±\\pm0\.0430\.00±\\pm0\.001GOLEM21\.30±\\pm4\.5180\.02±\\pm0\.0260\.54±\\pm0\.0490\.00±\\pm0\.001GOBNILP43\.42±\\pm7\.2360\.28±\\pm0\.4140\.01±\\pm0\.0180\.00±\\pm0\.001DECO20\.67±\\pm3\.6820\.40±\\pm0\.0300\.67±\\pm0\.0100\.03±\\pm0\.0075000PC18\.72±\\pm4\.0540\.34±\\pm0\.0720\.76±\\pm0\.0950\.03±\\pm0\.007GES24\.36±\\pm13\.3280\.39±\\pm0\.1490\.85±\\pm0\.1250\.05±\\pm0\.024HC43\.82±\\pm13\.7500\.55±\\pm0\.1110\.71±\\pm0\.1600\.08±\\pm0\.027NOTEARS35\.20±\\pm4\.2840\.09±\\pm0\.0950\.38±\\pm0\.1050\.01±\\pm0\.003GOLEM25\.86±\\pm5\.4550\.02±\\pm0\.0470\.53±\\pm0\.0890\.01±\\pm0\.007GOBNILP32\.72±\\pm6\.6900\.33±\\pm0\.1140\.29±\\pm0\.0770\.01±\\pm0\.004DECO19\.67±\\pm5\.6760\.34±\\pm0\.0890\.65±\\pm0\.0830\.03±\\pm0\.00710000PC19\.40±\\pm7\.0110\.33±\\pm0\.0780\.78±\\pm0\.1130\.03±\\pm0\.010GES27\.85±\\pm14\.1150\.40±\\pm0\.1170\.86±\\pm0\.0900\.05±\\pm0\.024HC52\.25±\\pm10\.3210\.58±\\pm0\.0590\.80±\\pm0\.0910\.10±\\pm0\.018NOTEARS26\.53±\\pm7\.5460\.02±\\pm0\.0310\.40±\\pm0\.0810\.01±\\pm0\.005GOLEM20\.13±\\pm6\.7100\.01±\\pm0\.0050\.55±\\pm0\.0750\.01±\\pm0\.003GOBNILP31\.64±\\pm6\.0390\.39±\\pm0\.0950\.31±\\pm0\.0640\.02±\\pm0\.004DECO18\.00±\\pm3\.1220\.08±\\pm0\.0230\.68±\\pm0\.0510\.03±\\pm0\.001Nodes: 25, Edge probability: 0\.0510000PC14\.7±\\pm3\.8640\.55±\\pm0\.1060\.88±\\pm0\.1030\.03±\\pm0\.008GES6\.54±\\pm3\.7210\.4±\\pm0\.1170\.95±\\pm0\.0850\.02±\\pm0\.008HC14\.98±\\pm6\.5160\.57±\\pm0\.1110\.69±\\pm0\.1260\.03±\\pm0\.011NOTEARS3\.07±\\pm1\.9310\.00±\\pm0\.0000\.82±\\pm0\.0850\.00±\\pm0\.000GOLEM1\.33±\\pm0\.9430\.00±\\pm0\.0020\.9±\\pm0\.0730\.00±\\pm0\.000GOBNILP8\.38±\\pm2\.9390\.47±\\pm0\.1450\.45±\\pm0\.1300\.01±\\pm0\.004DECO2\.87±\\pm2\.0290\.2±\\pm0\.1410\.84±\\pm0\.1180\.00±\\pm0\.003Nodes: 50, Edge probability: 0\.0510000PC35\.67±\\pm2\.8870\.41±\\pm0\.0590\.86±\\pm0\.0350\.02±\\pm0\.006GES20\.00±\\pm3\.6060\.29±\\pm0\.0120\.73±\\pm0\.4100\.01±\\pm0\.001HC86\.33±\\pm33\.0050\.64±\\pm0\.0250\.74±\\pm0\.0590\.04±\\pm0\.014NOTEARS14\.50±\\pm2\.7780\.00±\\pm0\.0010\.73±\\pm0\.0320\.00±\\pm0\.002GOLEM13\.63±\\pm4\.3680\.00±\\pm0\.0010\.77±\\pm0\.0560\.00±\\pm0\.004GOBNILP31\.62±\\pm14\.3430\.35±\\pm0\.1560\.34±\\pm0\.1500\.01±\\pm0\.003DECO11\.69±\\pm4\.1340\.05±\\pm0\.0020\.79±\\pm0\.0560\.00±\\pm0\.002Table 6:Comparison of DAG learning performance on Erdős\-Rényi Graphs and continuous \(Linear Gaussian\) dataData sizeMethodSHD↓\\downarrowFDR↓\\downarrowTPR↑\\uparrowFPR↓\\downarrowNodes: 25, Edge probability: 0\.151000PC32\.33±\\pm10\.1620\.41±\\pm0\.1000\.52±\\pm0\.1400\.03±\\pm0\.008GES33\.35±\\pm30\.0750\.42±\\pm0\.1980\.83±\\pm0\.1460\.06±\\pm0\.053HC78\.25±\\pm13\.1390\.71±\\pm0\.0660\.64±\\pm0\.0970\.13±\\pm0\.033NOTEARS4\.00±\\pm5\.1720\.02±\\pm0\.0380\.94±\\pm0\.0560\.00±\\pm0\.004GOLEM1\.25±\\pm2\.2220\.01±\\pm0\.0220\.99±\\pm0\.0270\.00±\\pm0\.002DECO25±\\pm3\.6110\.12±\\pm0\.0710\.85±\\pm0\.0130\.01±\\pm0\.004