Under stringent model type and variable distribution assumptions, differentiable score-based causal discovery methods learn a directed acyclic graph (DAG) from observational data by evaluating candidate graphs over an average score function.
Despite great success in low-dimensional linear systems, it has been observed that these approaches overly exploit easier-to-fit samples, thus inevitably learning spurious edges.
Worse still,

the common homogeneity assumption can be easily violated, due to the widespread existence of heterogeneous data in the real world, resulting in performance vulnerability when noise distributions vary.

We propose a simple yet effective model-agnostic framework to boost causal discovery performance by dynamically learning the adaptive weights for the Reweighted Score function, ReScore for short, where the weights tailor quantitatively to the importance degree of each sample.
Intuitively, we leverage the bilevel optimization scheme to alternately train a standard DAG learner and reweight samples — that is, upweight the samples the learner fails to fit and downweight the samples that the learner easily extracts the spurious information from.
Extensive experiments on both synthetic and real-world datasets are carried out to validate the effectiveness of ReScore.
We observe consistent and significant boosts in structure learning performance.
Furthermore, we visualize that ReScore concurrently mitigates the influence of spurious edges and generalizes to heterogeneous data.
Finally, we perform the theoretical analysis to guarantee the structure identifiability and the weight adaptive properties of ReScore in linear systems.
Our codes are available at <https://github.com/anzhang314/ReScore>.

§ INTRODUCTION

Learning causal structure from purely observational data (causal discovery) is a fundamental but daunting task [Chickering et al., 2004, Shen et al., 2020].
It strives to identify causal relationships between variables and encode the conditional independence as a directed acyclic graph (DAG).
Differentiable score-based optimization is a crucial enabler of causal discovery [Vowels et al., 2021].
Specifically, it is formulated as a continuous constraint optimization problem by minimizing the average score function and a smooth acyclicity constraint.
To ensure the structure is fully or partially identifiable (see Section <ref>), researchers impose stringent restrictions on model parametric family (linear, additive) and common assumptions of variable distributions (data homogeneity) [Peters et al., 2014, Ng et al., 2019].
Following this scheme, recent follow-on studies [Kalainathan et al., 2018, Ng et al., 2019, Zhu et al., 2020, Khemakhem et al., 2021, Yu et al., 2021] extend the formulation to general nonlinear problems by utilizing a variety of deep learning models.

However, upon careful inspections, we spot and justify two unsatisfactory behaviors of the current differentiable score-based methods:

* Differentiable score-based causal discovery is error-prone to learning spurious edges or reverse causal directions between variables, which derails the structure learning accuracy [He et al., 2021, Ng et al., 2022].
We substantiate our claim with an illustrative example as shown in Figure <ref> (see another example in Appendix <ref>).
We find that even the fundamental  chain structure in a linear system is easily misidentified by the state-of-the-art method, NOTEARS [Zheng et al., 2018].

* Despite being appealing in synthetic data, differentiable score-based methods suffer from severe performance degradation when encountering heterogeneous data [Huang et al., 2020, Huang et al., 2019].
Considering Figure <ref> again, NOTEARS is susceptible to learning redundant causations when the distributions of noise variables vary.

A simple example of basic chain structure that NOTEARS would learn spurious edges while ReScore can help to mitigate the bad influence.

Taking a closer look at this dominant scheme (optimizing the DAG learner via an average score function under strict assumptions), we ascribe these undesirable behaviors to its inherent limitations:

* The collected datasets naturally include an overwhelming number of easy samples and a small number of informative samples that might contain crucial causation information [Shrivastava et al., 2016].
Averagely scoring the samples deprives the discovery process of differentiating sample importance, thus easy samples dominate the learning of DAG.
As a result, prevailing score-based techniques fail to learn true causal relationship but instead yield the easier-to-fit spurious edges.

* Noise distribution shifts are inevitable and common in real-world training, as the observations are typically collected at different periods, environments, locations, and so forth [Arjovsky et al., 2019].
As a result, the strong assumption of noise homogeneity for differentiable DAG learner is easily violated in real-world data [Peters et al., 2016].
A line of works [Ghassami et al., 2018, Wang et al., 2022] dedicated to heterogeneous data can successfully address this issue.
However, they often require explicit domain annotations (ideal partition according to heterogeneity underlying the data) for each sample, which are prohibitively expensive and hard to obtain [Creager et al., 2021], thus further limiting their applicability.

To reshape the optimization scheme and resolve these limitations, we propose to adaptively reweight the samples, which de facto concurrently mitigates the influence of spurious edges and generalizes to heterogeneous data.
The core idea is to discover and upweight a set of less-fitted samples that offer additional insight into depicting the causal edges, compared to the samples easily fitted via spurious edges.
Focusing more on less-fitted samples enables the DAG learner to effectively generalize to heterogeneous data, especially in real-world scenarios whose samples typically come from disadvantaged domains.
However, due to the difficulty of accessing domain annotations, distinguishing such disadvantaged but informative samples and adaptively assigning their weights are challenging.

Towards this end, we present a simple yet effective model-agnostic optimization framework, coined ReScore, which automatically learns to reweight the samples and optimize the differentiable DAG learner, without any knowledge of domain annotations.
Specifically, we frame the adaptive weights learning and the differentiable DAG learning as a bilevel optimization problem, where the outer-level problem is solved subject to the optimal value of the inner-level problem:

* In the inner loop, the DAG learner is first fixed and evaluated by the reweighted score function to quantify the reliance on easier-to-fit samples, and then the instance-wise weights are adaptively optimized to induce the DAG learner to the worst-case.

* In the outer loop, upon the reweighted observation data where the weights are determined by the inner loop, any differential score-based causal discovery method can be applied to optimize the DAG learner and refine the causal structure.

Benefiting from this optimization scheme, our ReScore has three desirable properties.
First, it is a model-agnostic technique that can empower any differentiable score-based causal discovery method.
Moreover, we theoretically reveal that the structure identifiability is inherited by ReScore from the original causal discovery method in linear systems (Theorem <ref>).
Second, ReScore jointly mitigates the negative effect of spurious edge learning and performance drop in heterogeneous data via auto-learnable adaptive weights.
Theoretical analysis in Section <ref> (Theorem <ref>) validates the oracle adaptive properties of weights.
Third, ReScore boosts the causal discovery performance by a large margin.
Surprisingly, it performs competitively or even outperforms CD-NOD [Huang et al., 2020] and DICD [Wang et al., 2022], which require domain annotation, on heterogeneous synthetic data and real-world data (Section <ref>).

§ DIFFERENTIABLE CAUSAL DISCOVERY

We begin by introducing the task formulation of causal discovery and the identifiability issue.
We then present the differentiable score-based scheme to optimize the DAG learner.

Task Formulation.
Causal discovery aims to infer the Structural Causal Model (SCM) [Pearl, 2000, Pearl et al., 2016] from the observational data, which best describes the data generating procedure.
Formally, let $\Mat{X}\in\Space{R}^{n\times d}$ be a matrix of observational data, which consists of $n$ independent and identically distributed (i.i.d.) random vectors $X =(X_1, \dots, X_d) \in \Space{R}^d$.
Given $\Mat{X}$, we aim to learn a SCM $(P_X, \Set{G})$, which encodes a causal directed acyclic graph (DAG) with a structural equation model (SEM) to reveal the data generation from the distribution of variables $X$.
Specifically, we denote the DAG by $\Set{G}=(V(\Set{G}),E(\Set{G}))$, where $V(\Set{G})$ is the variable set and $E(\gG)$ collects the causal directed edges between variables.
We present the joint distribution over $X$ as $P_X$, which is Markov $\Set{G}$.
The probability distribution function of $P_{X}$ is factored as $p(x)= \prod_{i=1}^d P(x_i|x_{pa(i)})$, where $pa(i) = \{j \in V(\gG): X_j \rightarrow X_i \in E(\gG)\}$ is the set of parents of variable $X_i$ in $\gG$ and $P(x_i|x_{pa(i)})$ is the conditional probability density function of variable $X_i$ given $X_{pa(i)}$.
As a result, the SEM can be formulated as a collection of $d$ structural equations:

\begin{gather}\label{eq:SEM}
X_i = f_i(X_{pa(i)},N_i), \quad i = 1,\cdots,d
\end{gather}

where $f_i:\Space{R}^{|X_{pa(i)}|} \rightarrow \Space{R}$ can be any linear or nonlinear function, and $N = (N_1,\dots,N_d)$ are jointly independent noise variables.

Identifiability Issue.
In general, without further assumption on the SEM (eq:SEM), it is not possible to uniquely learn the DAG $\gG$ by only using the observations of $P_X$.
This is the identifiability issue in causal discovery [Lachapelle et al., 2020].
Nonetheless, with the assumption of the SEM, the DAG $\Set{G}$ is said to be identifiable over $P_X$, if no other SEM can encode the same distribution $P_X$ with a different DAG under the same assumption.
To guarantee the identifiability, most prior studies restrict the form of the structural equations to be additive to noises, additive noise models (ANM).
Assuming ANM, as long as the structural equations are linear with non-Gaussian errors [Shimizu et al., 2006, Loh & Bühlmann, 2014], linear Gaussian model with equal noise variances [Peters & Bühlmann, 2014], or nonlinear structural equation model with mild conditions [Hoyer et al., 2008, Zhang & Hyvarinen, 2009, Peters et al., 2014],
then the DAG $\gG$ is identifiable.

Solution to Causal Discovery.
Prevailing causal discovery approaches roughly fall into two lines: constraint- and score-based methods [Spirtes & Zhang, 2016, Glymour et al., 2019].
Specifically, constraint-based methods [Spirtes et al., 1995, Spirtes & Glymour, 1991, Colombo et al., 2012] determine up to the Markov equivalence class of causal graphs, based on conditional independent tests under certain assumptions.
Score-based methods [Vowels et al., 2021] evaluate the candidate graphs with a predefined score function and search the DAG space for the optimal graph.
Here we focus on the score-based line.

Score-based Causal Discovery.
With a slight abuse of notation, $\Set{G}$ refers to a directed graph in the rest of the paper.
Formally, the score-based scheme casts the task of DAG learning as a combinatorial optimization problem:

\begin{align}\label{eq:score-comb}
\min_{\Set{G}} S(\Set{G};\Mat{X}) = \Lapl(\Set{G};\Mat{X}) + \lambda\Set{R}_{\text{sparse}}(\Set{G})
\quad\textrm{s.t.}\quad\Set{G}\in \text{DAG},
\end{align}

Here this problem consists of two ingredients: the combinatorial acyclicity constraint $\Set{G}\in \text{DAG}$ and the score function $S(\Set{G};\Mat{X})$.
The score function composes two terms:
(1) the goodness-of-fit measure $\Lapl(\Set{G};\Mat{X}) = \frac{1}{n}\sum_{i=1}^{n}l(\Mat{x}_i, f(\Mat{x}_i))$, where $l(\Mat{x}_i, f(\Mat{x}_i))$ represents the loss of fitting observation $\Mat{x}_i$;
(2) the sparsity regularization $\Set{R}_{\text{sparse}}(\Set{G})$ stipulating that the total number of edges in $\Set{G}$ should be penalized;
And $\lambda$ is a hyperparameter controlling the regularization strengths.
Next, we will elaborate on the previous implementations of these two major ingredients.

To implement $S(\Set{G};\Mat{X})$, various approaches have been proposed, such as penalized least-squares loss [Zheng et al., 2020, Zheng et al., 2018, Ng et al., 2019], Evidence Lower Bound (ELBO) [Yu et al., 2019], log-likelihood with complexity regularizers [Kalainathan et al., 2018, Van de Geer & Bühlmann, 2013, Ng et al., 2020], Maximum Mean Discrepancy (MMD) [Goudet et al., 2018], Bayesian Information Criterion (BIC) [Geiger & Heckerman, 1994, Zhu et al., 2020], Bayesian Dirichlet equivalence uniform (BDeu) score [Heckerman et al., 1995], Bayesian Gaussian equivalent (BGe) score [Kuipers et al., 2014], and others [Huang et al., 2018, Bach & Jordan, 2002, Sokolova et al., 2014].

As $\Set{G}\in \text{DAG}$ enforces $\Set{G}$ to be acyclic, it becomes the main obstacle to the score-based scheme.
Prior studies propose various approaches to search in the acyclic space, such as greedy search [Chickering, 2002, Hauser & Bühlmann, 2012], hill-climbing [Gámez et al., 2011, Tsamardinos et al., 2006], dynamic programming [Silander & Myllymäki, 2006, Koivisto & Sood, 2004], A* [Yuan & Malone, 2013], integer linear programming [Jaakkola et al., 2010, Cussens, 2011].

Differentiable Score-based Optimization.
Different from the aforementioned search approaches, NOTEARS [Zheng et al., 2018] reframes the combinatorial optimization problem as a continuous constrained optimization problem:

\begin{gather}\label{eq:constrained-opti}
\min_{\Set{G}}~S(\Set{G};\Mat{X}) \quad \text{s.t.}\quad H(\Set{G}) = 0,
\end{gather}

where $H(\Set{G})=0$ is a differentiable equality DAG constraint.

As for the DAG constraint $H(\Set{G}) = 0$, the prior effort [Zheng et al., 2018] turns to depict the “DAGness” of $\Set{G}$'s adjacency matrix $\Set{A}(\Set{G}) \in \{0,1\}^{d\times d}$.
Specifically, $[\Set{A}(\Set{G})]_{ij} =1$ if the causal edge $X_j \rightarrow X_i$ exists in $E(\Set{G})$, otherwise $[\Set{A}(\Set{G})]_{ij} =0$.
Prevailing implementations of DAGness constraints are $H(\gG) = \Tr(e^{\gA \odot \gA})-d$ [Zheng et al., 2018], $H(\gG) = \Tr[(I+\alpha \gA\odot\gA)^d]-d$ [Yu et al., 2019], and others [Wei et al., 2020, Kyono et al., 2020, Bello et al., 2022, Zhu et al., 2021].
As a result, this optimization problem in eq:constrained-opti can be further formulated via the augmented Lagrangian method as:

\begin{gather} \label{eq:score-cont}
\min_{\gG}~S(\gG;\rmX) + \mathcal{P}_{\text{DAG}}(\gG),
\end{gather}

where $\mathcal{P}_{\text{DAG}}(\gG) = \alpha H(\gG) + \frac{\rho}{2}|H(\gG)|^2$ is the penalty term enforcing the DAGness on $\gG$, and $\rho >0$ is a penalty parameter and $\alpha$ is the Lagrange multiplier.

§ METHODOLOGY OF RESCORE

On the basis of differentiable score-based causal discovery methods, we first devise our ReScore and then present its desirable properties.

§.§ Bilevel Formulation of ReScore

Aiming to learn the causal structure accurately in practical scenarios, we focus on the observational data that is heterogeneous and contains a large proportion of easy samples.
Standard differentiable score-based causal discovery methods apply the average score function on all samples equally, which inherently rely on easy samples to obtain high average goodness-of-fit.
As a result, the DAG learner is error-prone to constructing easier-to-fit spurious edges based on the easy samples, while ignoring the causal relationship information maintained in hard samples.
Assuming the oracle importance of each sample is known at hand, we can assign distinct weights to different samples and formulate the reweighted score function $S_\ervw(\gG;\rmX)$, instead of the average score function:

\begin{align} \label{eq:score-weighted}
S_\ervw(\gG;\rmX) = \Ls_\ervw(\gG;\rmX) + \lambda\mathcal{R}_{\text{sparse}}(\gG)  = \sum_{i=1}^{n} w_i l(\Mat{x}_i, f(\Mat{x}_i)) + \lambda\mathcal{R}_{\text{sparse}}(\gG),
\end{align}

where $\Mat{w} = (w_1,\dots,w_n)$ is a sample reweighting vector with length $n$, wherein $w_i$ indicates the importance of the $i$-th observed sample $\Mat{x}_{i}$.

However, the oracle sample importance is usually unavailable in real-world scenarios. The problem, hence, comes to how to automatically learn appropriate the sample reweighting vector $\Mat{w}$.
Intuitively, samples easily fitted with spurious edges should contribute less to the DAG learning, while samples that do not hold spurious edges but contain critical information about causal edges should be more importance.
We therefore use a simple heuristic of downweighting the easier-to-fit but less informative samples, and upweighting the less-fitted but more informative samples.
This inspires us to learn to allocate weights adaptively, with the aim of maximizing the influence of less well-fitted samples and failing the DAG learner.
Formally, we cast the overall framework of reweighting samples to boost causal discovery as the following bilevel optimization problem:

\begin{align} \label{maxweight}
\min_{\Set{G}}~& S_{\ervw^*}(\gG;\rmX) + \mathcal{P}_{DAG}(\gG), \nonumber\\
\textrm{s.t.}~& \Mat{w}^* \in \argmax_{\Mat{w} \in \sC(\tau)} \, S_{\ervw}(\gG;\rmX),
\end{align}

where $ \sC(\tau) := \{\Mat{w}: 0 < \frac{\tau}{n} \leq w_1, \dots, w_n \leq \frac{1}{\tau n}, \sum_{i=1}^{n} w_i = 1\}$ for the cutoff threshold $ \tau \in (0, 1) $.
The deviation of the weight distribution from the uniform distribution is bound by the hyperparameter $\tau$.
Clearly, maxweight consists of two objectives, where the inner-level objective (optimize $\Mat{w}$ by maximizing the reweighted score function) is nested within the outer-level objective (optimize $\Set{G}$ by minimizing the differentiable score-based loss).
Solving the outer-level problem should be subject to the optimal value of the inner-level problem.

Now we introduce how to solve this bilevel optimization problem.
In the inner loop, we first fix the DAG learner which evaluates the error of each observed sample $\Mat{x}_{i}$, $\forall i\in\{1,\cdots,n\}$, and then maximize the reweighted score function to learn the weight $w_i^*$ correspondingly.
In the outer loop, upon the reweighted observations whose weights are determined in the inner loop, we minimize the reweighted score function to optimize the DAG learner.
By alternately training the inner and outer loops, the importance of each sample is adaptively estimated based on the DAG learner's error, and in turn gradually guides the DAG learner to perform better on the informative samples.
It is worth highlighting that this ReScore scheme can be applied to any differentiable score-based causal discovery method listed in Section <ref>.
The procedure of training ReScore is outlined in Algorithm <ref>.

Furthermore, our ReScore has the following desirable advantages:

* As shown in Section <ref>, under mild conditions, our ReScore inherits the identifiability property of the original differentiable score-based causal discovery method.

* ReScore is able to generate adaptive weights to observations through the bilevel optimize, so as to distinguish more information samples and fulfill their potentials to guide the DAG learning. This is consistent with our theoretical analysis in Section <ref> and empirical results in Section <ref>.

* ReScore is widely applicable to various types of data and models. In other words, it is model-agnostic and can effectively handle heterogeneous data without knowing the domain annotations in advance. Detailed ReScore performance can be found in Section <ref>.

§.§ Theoretical Analysis on Identifiability

The graph identifiability issue is the primary challenge hindering the development of structure learning.
As an optimization framework, the most desired property of ReScore is the capacity to ensure graph identifiability and substantially boost the performance of the differentiable score-based DAG learner.
We develop Theorem <ref> that guarantees the DAG identifiability when using ReScore.

Rendering a DAG theoretically identifiable requires three standard steps [Peters et al., 2014, Zheng et al., 2020, Ng et al., 2022]:
(1) assuming the particular restricted family of functions and data distributions of SEM in eq:SEM;
(2) theoretically proving the identifiability of SEM;
and (3) developing an optimization algorithm with a predefined score function and showing that learned DAG asymptotically converges to the ground-truth DAG.
Clearly, ReScore naturally inherits the original identifiability of a specific SEM as stated in Section <ref>.
Consequently, the key concern lies on the third step — whether the DAG learned by our new optimization framework with the reweighted score function $S_\ervw(\gG;\rmX)$ can asymptotically converge to the ground-truth DAG.
To address this, we present the following theorem.
Specifically, it demonstrates that, by guaranteeing the equivalence of optimization problems (eq:score-comb and maxweight) in linear systems, the bounded weights will not affect the consistency results in identifiability analysis.
See detailed proof in Appendix <ref>.
Suppose the SEM in eq:SEM is linear and the size of observational data $\Mat{X}$ is $n$. As the data size increases, $ n \to \infty $,

\begin{align*}
\argmin_{\gG} \big\{ S_\Mat{w}(\gG; \Mat{X}) + \mathcal{P}_{DAG}(\gG) \big\} - \argmin_{\gG} \big\{ S(\gG; \Mat{X}) + \mathcal{P}_{DAG}(\gG) \big\} \xrightarrow[]{a.s.} \bf{0}
\end{align*}

in the following cases:

a. Using the least-squares loss $\Lapl(\Set{G};\Mat{X}) = \frac{1}{2n}\Vert \Mat{X} - f(\Mat{X}) \Vert_F^2$;

b. Using the negative log-likelihood loss with standard Gaussian noise.

Remark: The identifiability property of ReScore with two most common score functions, namely least-square loss and negative log-likelihood loss, is proved in Theorem <ref>.
Similar conclusions can be easily derived for other loss functions, which we will explore in future work.

§.§ Oracle Property of Adaptive Weights

Our ReScore suggests assigning varying degrees of importance to different observational samples. At its core is the simple yet effective heuristic: the less-fitted samples are more important than the easier-to-fit samples, as they do not hold spurious edges but contain critical information about the causal edges.
Hence, mining hard-to-learn causation information is promising to help DAG learners mitigate the negative influence of spurious edges.
The following theorem shows the adaptiveness property of ReScore, instead of equally treating all samples, ReScore tends to upweight the importance of hard but informative samples while downweighting the reliance on easier-to-fit samples.

Suppose that in the optimization phase, the $i$-th observation has a larger error than the j-th observation in the sense that $l(\Mat{x}_i, f(\Mat{x}_i)) > l(\Mat{x}_j, f(\Mat{x}_j))$, where $i, j \in \{1, \dots, n\}$. Then,

$$w^*_i \czb{\geq} w^*_j,$$

where $w_i^*, w_i^*$ are the optimal weights in maxweight. The equality holds if and only if $ w^*_i = w^*_j = \frac{\tau}{n} $ or $ w^*_i = w^*_j = \frac{1}{\tau n} $.

See Appendix <ref> for the detailed proof.
It is simple to infer that, following the inner loop that maximizes the reweighted score function $S_{\ervw}(\gG;\rmX)$, the observations are ranked by learned adaptive weights $\Mat{w}^*$.
That is, one observation equipped with a higher weight will have a greater impact on the subsequent outer loop to dominate the DAG learning.

§ EXPERIMENTS

We aim to answer the following research questions:

* RQ1: As a model-agnostic framework, can ReScore widely strengthen the differentiable score-based causal discovery baselines?

* RQ2: How does ReScore perform when noise distribution varies? Can ReScore effectively learn the adaptive weights that successfully identify the important samples?

To answer the first question (RQ1), we implement various backbone models

including NOTEARS [Zheng et al., 2018] and GOLEM [Ng et al., 2020] in linear systems, and NOTEARS-MLP [Zheng et al., 2020], and GraN-DAG [Lachapelle et al., 2020] in nonlinear settings.
To answer the second question (RQ2), we compare GOLEM+ReScore, NOTEARS-MLP+ReScore to a SOTA baseline CD-NOD [Huang et al., 2020] and a recently proposed approach DICD [Wang et al., 2022], which both require the ground-truth domain annotation.
For a comprehensive comparison, extensive experiments are conducted on both homogeneous and heterogeneous synthetic datasets as well as a real-world benchmark dataset, Sachs [Sachs et al., 2005].
In Sachs, GES [Chickering, 2002], a benchmark discrete score-based causal discovery method, is also considered.
A detailed description of the employed baselines can be found in Appendix <ref>.

Results for ER graphs of 10 nodes on linear and nonlinear synthetic datasets.

4c|ER2                                                                                4cER4

1cTPR  $\mathbf{\uparrow}$     1cFDR $\mathbf{\downarrow}$     1cSHD $\mathbf{\downarrow}$     1c|SID $\mathbf{\downarrow}$     1cTPR  $\mathbf{\uparrow}$     FDR $\mathbf{\downarrow}$     SHD $\mathbf{\downarrow}$     SID $\mathbf{\downarrow}$
Random          0.08$\pm$0.07     0.93$\pm$0.18     33.2$\pm$7.3     95.6$\pm$12.2     0.09$\pm$0.17     0.93$\pm$0.09     52.3$\pm$16.7     80.3$\pm$17.7

NOTEARS         0.85$\pm$0.09
0.07$\pm$0.07     5.8$\pm$2.2      20.8$\pm$5.2     0.79$\pm$0.11     0.09$\pm$0.05      10.0$\pm$5.2      25.8$\pm$9.9

+ ReScore       $\textbf{0.89}\text{\scriptsize{$\pm$0.07}}^{\color{+}+5\%}$       $0.08\text{\scriptsize{$\pm$0.09}}^{\color{-}-12 \%}$      $\textbf{4.6}\text{\scriptsize{$\pm$2.3}}^{\color{+}+26 \%}$     $12.8\text{\scriptsize{$\pm$7.0}}^{\color{+}+63 \%}$      $\textbf{0.85}\text{\scriptsize{$\pm$0.04}}^{\color{+}+8\%}$     $\textbf{0.05}\text{\scriptsize{$\pm$0.04}}^{\color{+}+57\%}$      $\textbf{7.2}\text{\scriptsize{$\pm$1.9}}^{\color{+}+39 \%}$      $\textbf{24.2}\text{\scriptsize{$\pm$8.4}}^{\color{+}+7 \%} $

GOLEM           0.87$\pm$0.06
0.22$\pm$0.11     6.5$\pm$3.4      13.0$\pm$6.7     0.63$\pm$0.03     0.16$\pm$0.03      17.2$\pm$1.3      48.0$\pm$13.3

+ ReScore       $0.88\text{\scriptsize{$\pm$0.06}}^{\color{+}+1\%}$       $0.21\text{\scriptsize{$\pm$0.11}}^{\color{+}+2 \%}$      $6.0\text{\scriptsize{$\pm$3.4}}^{\color{+}+8 \%}$     $\textbf{12.4}\text{\scriptsize{$\pm$6.3}}^{\color{+}+5 \%}$      $0.66\text{\scriptsize{$\pm$0.04}}^{\color{+}+5\%}$     $0.17\text{\scriptsize{$\pm$0.01}}^{\color{-}-5\%}$      $16.2\text{\scriptsize{$\pm$1.0}}^{\color{+}+6 \%}$      $46.7\text{\scriptsize{$\pm$13.3}}^{\color{+}+3 \%} $
NOTEARS-MLP     0.76$\pm$0.17
0.14$\pm$0.09     7.0$\pm$3.5      17.9$\pm$10.0     0.83$\pm$0.05     0.21$\pm$0.04      10.9$\pm$1.9      28.6$\pm$12.0

+ ReScore       $0.73\text{\scriptsize{$\pm$0.07}}^{\color{-}-4\%}$       $0.10\text{\scriptsize{$\pm$0.09}}^{\color{+}+37 \%}$      $6.8\text{\scriptsize{$\pm$2.9}}^{\color{+}+3 \%}$     $20.3\text{\scriptsize{$\pm$9.7}}^{\color{-}-11 \%}$      $0.94\text{\scriptsize{$\pm$0.06}}^{\color{+}+14\%}$     $0.15\text{\scriptsize{$\pm$0.06}}^{\color{+}+44\%}$      $6.80\text{\scriptsize{$\pm$2.7}}^{\color{+}+60 \%}$      $8.80\text{\scriptsize{$\pm$12.4}}^{\color{+}+225 \%} $

GraN-DAG        0.88$\pm$0.06
0.02$\pm$0.03     2.7$\pm$1.6      8.70$\pm$4.8     0.98$\pm$0.02     0.12$\pm$0.03      5.4$\pm$1.1      3.70$\pm$4.71

+ ReScore       $\textbf{0.90}\text{\scriptsize{$\pm$0.05}}^{\color{+}+2\%}$       $\textbf{0.01}\text{\scriptsize{$\pm$0.03}}^{\color{+}+35 \%}$      $\textbf{2.4}\text{\scriptsize{$\pm$1.1}}^{\color{+}+13 \%}$     $\textbf{7.20}\text{\scriptsize{$\pm$3.0}}^{\color{+}+21 \%}$      $\textbf{0.99}\text{\scriptsize{$\pm$0.01}}^{\color{+}+1\%}$     $\textbf{0.11}\text{\scriptsize{$\pm$0.01}}^{\color{+}+12\%}$      $\textbf{4.80}\text{\scriptsize{$\pm$0.6}}^{\color{+}+13 \%}$      $\textbf{0.50}\text{\scriptsize{$\pm$0.81}}^{\color{+}+640 \%} $

TPR $\lambda$

FDR $\lambda$

SHD $\lambda$

SID $\lambda$

Performance comparison between NOTEARS-MLP and ReScore on ER4 graphs of 10 nodes on nonlinear synthetic datasets. The hyperparameter $\lambda$ defined in eq:score-comb refers to the graph sparsity. See more results in Appendix <ref>

Evaluation Metrics.
To evaluate the quality of structure learning, four metrics are reported: True Positive Rate (TPR), False Discovery Rate (FDR), Structural Hamming Distance (SHD), and Structural Intervention Distance (SID) [Peters & Bühlmann, 2015], averaged over ten random trails.

§.§ Overall Performance Comparison on Synthetic Data (RQ1)

The generating data differs along three dimensions: number of nodes, the degree of edge sparsity, and the type of graph.
Two well-known graph sampling models, namely Erdos-Renyi (ER) and scale-free (SF) [Barabási & Albert, 1999] are considered with $kd$ expected edges (denoted as ER$k$ or SF$k$) and $d=\{10,20,50\}$ nodes.
Specifically, in linear settings, similar to [Zheng et al., 2018, Gao et al., 2021], the coefficients are assigned following $U(-2,-0.5) \cup U(0.5,2)$ with additive standard Gaussian noise.
In nonlinear settings, following [Zheng et al., 2020], the ground-truth SEM in eq:SEM is generated under the Gaussian process (GP) with radial basis function kernel of bandwidth one, where $f_i(\cdot)$ is additive noise models with $N_i$ as an i.i.d. random variable following standard normal distribution.
Both of these settings are known to be fully identifiable [Peters & Bühlmann, 2014, Peters et al., 2014].
For each graph, 10 data sets of 2,000 samples are generated and the mean and standard deviations of the metrics are reported for a fair comparison.

Tables <ref>, <ref> and Tables in Appendix <ref> report the empirical results on both linear and nonlinear synthetic data.
The error bars depict the standard deviation across datasets over ten trails.
The red and blue percentages separately refer to the increase and decrease of ReScore relative to the original score-based methods in each metric.
The best performing methods are bold.
We find that:

* ReScore consistently and significantly strengthens the score-based methods for structure learning across all datasets.
In particular, it achieves substantial gains over the state-of-the-art baselines by around 3% to 60% in terms of SHD, revealing a lower number of missing, falsely detected, and reversed edges.
We attribute the improvements to the dynamically learnable adaptive weights, which boost the quality of score-based DAG learners.
With a closer look at the TPR and FDR, ReScore typically lowers FDR by eliminating spurious edges and enhances TPR by actively identifying more correct edges.
This clearly demonstrates that ReScore effectively filters and upweights the more informative samples to better extract the causal relationship.
Figure <ref> also illustrates the clear trend that ReScore is excelling over NOTEARS-MLP as the sparsity penalty climbs.
Additionally, as Table <ref> indicates, ReScore only adds a negligible amount of computational complexity as compared to the backbone score-based DAG learners.

* Score-based causal discovery baselines suffer from a severe performance drop on high-dimensional dense graph data.
Despite the advances, beyond linear, NOTEARS-MLP and GraN-DAG fail to scale to more than 50 nodes in SF4 and ER4 graphs, mainly due to difficulties in enforcing acyclicity in high-dimensional dense graph data [Varando, 2020, Lippe et al., 2022].
Specifically, the TPR of GraN-DAG and NOTEARS-MLP in SF4 of 50 nodes is lower than 0.2, which indicates that they are not even able to accurately detect 40 edges out of 200 ground-truth edges.
ReScore, as an optimization framework, relies heavily on the performance of the score-based backbone model.
When the backbone model fails to infer DAG on its own as the number of nodes and edge density increase, adding ReScore will not be able to enhance the performance.

§.§ Performance on Heterogeneous Data (RQ2)

§.§.§ Evaluation on Synthetic Heterogeneous Data

Weights on linear data at the 1st and last epochs

Weights on nonlinear data at the 1st and last epochs

Illustration of adaptive weights learned by ReScore sample loss on both linear and nonlinear synthetic data.
For each dataset, the left and right plots refer to the distribution of adaptive weights at the first and last epochs in the outer loop, respectively (the value of $\Mat{w}^*$, when $k_1=0$ and $k_1=K_{outer}$ in Algorithm <ref>, respectively).
The disadvantaged but more informative samples are represented by the red dots.
The dominant and easy samples, in contrast, are in blue.

It is commonplace to encounter heterogeneous data in real-world applications, of which the underlying causal generating process remain stable but the noise distribution may vary.
Specific DAG learners designed for heterogeneous data are prone to assume strict conditions and require the knowledge of group annotation for each sample.
Group annotations, however, are extremely costly and challenging to obtain.
We conjecture that a robust DAG learner is able to successfully handle heterogeneous data without the information of group annotation.

Synthetic heterogeneous data in both linear and nonlinear settings ($n=1000$, $d = 20$, ER2) containing two distinct groups are also considered.

10% of observations come from the disadvantaged group, where half of the noise variables $N_i$ defined in eq:SEM follow $\mathcal{N}(0,1)$ and the remaining half of noise variables follow $\mathcal{N}(0,0.1)$.
90% of the observations, in contrast, are generated from the dominant group where the scales of noise variables are flipped.

To evaluate whether ReScore can handle heterogeneous data without requiring the group annotation by automatically identifying and upweighting informative samples, we compare baseline+ReScore to CD-NOD and DICD, two SOTA causal discovery approaches that rely on group annotations and are developed for heterogeneous data.
Additionally, a non-adaptive reweighting method called baseline+IPS is taken into account, in which sample weights are inversely proportional to group sizes.
Specifically, we divide the whole observations into two subgroups.
Obviously, a single sample from the disadvantaged group is undoubtedly more informative than a sample from the dominant group, as it offers additional insight to depict the causal edges.

Results on heterogeneous data.

Linear     1cTPR$\mathbf{\uparrow}$     1cFDR$\mathbf{\downarrow}$     1c|SHD$\mathbf{\downarrow}$     Nonlinear     1cTPR$\mathbf{\uparrow}$     1cFDR$\mathbf{\downarrow}$     1cSHD$\mathbf{\downarrow}$
GOLEM           0.79                                 0.33                                 18.7                                 NOTEARS-MLP     0.62                                 0.36                                 25.8

+  IPS      0.65                                 0.19                                 18.6                                 + IPS       0.35                                 0.21                                 28.7

+  ReScore      0.81                                 0.24                                 16.4                                 + ReScore       0.63                                 0.32                                 23.8
CD-NOD          0.51                                 0.17                                 24.1                                 CN-NOD          0.60                                 0.29                                 26.0

DICD            0.82                                 0.28                                 16.7                                 DICD            0.50                                 0.24                                 23.5

As Figure <ref> shows, dots of different colours are mixed and scattered at the beginning of the training.
After multiple iterations of training in inner and outer loops, the red dots from the disadvantaged group are gradually identified and assigned to relatively larger weights as compared to those blue dots with the same measure-of-fitness.
This illustrates the effectiveness of ReScore and further offers insight into the reason for its performance improvements when handling heterogeneous data.
Overall, all figures show clear positive trends, the underrepresented samples tend to learn bigger weights.
These results validate the property of adaptive weights in Theorem <ref>.

Table <ref> indicates that ReScore drives impressive performance breakthroughs in heterogeneous data, achieving competitive or even lower SHD without group annotations compared to CD-NOD and DICD recognized as the lower bound.
Specifically, both GOLEM and NOTEARS-MLP are struggling from notorious performance drop when homogeneity assumption is invalidated, and posing hurdle from being scaled up to real-world large-scale applications.
We ascribe this hurdle to blindly scoring the observational samples evenly, rather than distilling the crucial group information from distribution shift of noise variables.
To better highlight the significance of the adaptive property, we also take Baseline+IPS into account, which views the ratio of group size as the propensity score and exploits its inverse to re-weight each sample's loss.
Baseline+IPS suffers from severe performance drops in terms of TPR, revealing the limitation of fixed weights.
In stark contrast, benefiting from adaptive weights, ReScore can even extract group information from heterogeneous data that accomplish more profound causation understanding, leading to higher DAG learning quality.
This validates that ReScore endows the backbone score-based DAG learner with better robustness against the heterogeneous data and alleviates the negative influence of spurious edges.

§.§.§ Evaluations on Real Heterogeneous Data.

The performance comparison on Sachs dataset.

1lTPR  $\mathbf{\uparrow}$     1lFDR $\mathbf{\downarrow}$     1lSHD $\mathbf{\downarrow}$     1lSID $\mathbf{\downarrow}$     1l#Predicted Edges
Random          0.076                       0.899                       23                          63                                   22
GOLEM           0.176                       0.026                       15                          53                                   4

+ ReScore       0.294                       0.063                       14                 49                                   6
NOTEARS-MLP     0.412                       0.632                       16                          45                                   19

+ ReScore       0.412                       0.500                       13                 43                                   14
GraN-DAG        0.294                       0.643                       16                          60                                   14

+ ReScore       0.353                       0.600                       15                 58                                   15

GES        0.294                       0.853                       31                          54                                   34

+ ReScore       0.588                       0.722                       28                 50                                   36

CD-NOD        0.588                       0.444                       15                          -                                   18

Sachs [Sachs et al., 2005] contains the measurement of multiple phosphorylated protein and phospholipid components simultaneously in a large number of individual primary human immune system cells.
In Sachs, nine different perturbation conditions are applied to sets of individual cells, each of which administers certain reagents to the cells.
With the annotations of perturbation conditions, we consider the Sachs as real-world heterogeneous data [Mooij et al., 2020].
We train baselines on 7,466 samples, where the ground-truth graph (11 nodes and 17 edges) is widely accepted by the biological community.

As Table <ref> illustrates, ReScore steadily and prominently boosts all baselines, including both differentiable and discrete score-based causal discovery approaches SHD and SID metrics.
This clearly shows the effectiveness of ReScore to better mitigate the reliance on easier-to-fit samples.
With a closer look at the TPR and FDR, baseline+ReScore surpasses the state-of-the-art corresponding baseline by a large margin in most cases, indicating that ReScore can help successfully predict more correct edges and fewer false edges.
Remarkably, compared to CD-NOD, which is designed for heterogeneous data and utilizes the annotations as prior knowledge, GES+ReScore obtains competitive TPR without using ground-truth annotations.
Moreover, GraN-DAG+ReScore can reach the same SHD as CD-NOD when 15 and 18 edges are predicted, respectively.
These findings validate the potential of ReScore as a promising research direction for enhancing the generalization and accuracy of DAG learning methods when dealing with real-world data.

§ CONCLUSION

Today's differentiable score-based causal discovery approaches are still far from being able to accurately detect the causal structures, despite their great success on synthetic linear data.
In this paper, we proposed ReScore, a simple-yet-effective model-agnostic optimization framework that simultaneously eliminates spurious edge learning and generalizes to heterogeneous data by utilizing learnable adaptive weights.
Grounded by theoretical proof and empirical visualization studies, ReScore successfully identifies the informative samples and yields a consistent and significant boost in DAG learning.
Extensive experiments verify that the remarkable improvement of ReScore on a variety of synthetic and real-world datasets indeed comes from adaptive weights.

Two aspects of ReScore's limitations will be covered in subsequent works.
First, the performance of ReScore is highly related to the causal discovery backbone models, which leads to minor improvements when the backbone methods fail.
Second, having empirically explored the sensitivity to pure noise samples in <ref>, we will theoretically analyze and further enhance the robustness of ReScore against these noises.

It is expected to substantially improve the DAG learning quality, as well as distinguish true informative samples from pure noise samples.
We believe that ReScore provides a promising research direction to diagnose the performance degradation for nonlinear and heterogeneous data in the structure learning challenge and will inspire more works in the future.

§.§.§ Acknowledgments

This research is supported by the Sea-NExT Joint Lab, the National Natural Science Foundation of China (9227010114), and CCCD Key Lab of the Ministry of Culture and Tourism.

[Arjovsky et al., 2019]
Martin Arjovsky, Léon Bottou, Ishaan Gulrajani, and David Lopez-Paz.
Invariant risk minimization.
arXiv preprint arXiv:1907.02893, 2019.

[Bach & Jordan, 2002]
Francis R. Bach and Michael I. Jordan.
Learning graphical models with mercer kernels.
In NIPS, 2002.

[Barabási & Albert, 1999]
Albert-László Barabási and Réka Albert.
Emergence of scaling in random networks.
science, 2860 (5439):0 509–512, 1999.

[Bello et al., 2022]
Kevin Bello, Bryon Aragam, and Pradeep Ravikumar.
Dagma: Learning dags via m-matrices and a log-determinant acyclicity
arXiv preprint arXiv:2209.08037, 2022.

[Cai et al., 2020]
Ruichu Cai, Jincheng Ye, Jie Qiao, Huiyuan Fu, and Zhifeng Hao.
Fom: Fourth-order moment based causal direction identification on the
heteroscedastic data.
Neural Networks, 124:0 193–201, 2020.

[Chickering, 2002]
David Maxwell Chickering.
Optimal structure identification with greedy search.
Journal of machine learning research, 30
(Nov):0 507–554, 2002.

[Chickering et al., 2004]
David Maxwell Chickering, David Heckerman, and Christopher Meek.
Large-sample learning of bayesian networks is np-hard.
Journal of machine learning research, 5:0 1287–1330,

[Colombo et al., 2012]
Diego Colombo, Marloes H Maathuis, Markus Kalisch, and Thomas S Richardson.
Learning high-dimensional directed acyclic graphs with latent and
selection variables.
The Annals of Statistics, pp.  294–321, 2012.

[Creager et al., 2021]
Elliot Creager, Jörn-Henrik Jacobsen, and Richard Zemel.
Environment inference for invariant learning.
In ICML, 2021.

[Cussens, 2011]
James Cussens.
Bayesian network learning with cutting planes.
In UAI, 2011.

[Gámez et al., 2011]
José A Gámez, Juan L Mateo, and José M Puerta.
Learning bayesian networks by hill climbing: efficient methods based
on progressive restriction of the neighborhood.
Data Mining and Knowledge Discovery, 220 (1):0
106–148, 2011.

[Gao et al., 2021]
Yinghua Gao, Li Shen, and Shu-Tao Xia.
DAG-GAN: causal structure learning with generative adversarial
In ICASSP, 2021.

[Geiger & Heckerman, 1994]
Dan Geiger and David Heckerman.
Learning gaussian networks.
In UAI, 1994.

[Ghassami et al., 2017]
AmirEmad Ghassami, Saber Salehkaleybar, Negar Kiyavash, and Kun Zhang.
Learning causal structures using regression invariance.
In NIPS, 2017.

[Ghassami et al., 2018]
AmirEmad Ghassami, Negar Kiyavash, Biwei Huang, and Kun Zhang.
Multi-domain causal structure learning in linear systems.
In NeurIPS, 2018.

[Glymour et al., 2019]
Clark Glymour, Kun Zhang, and Peter Spirtes.
Review of causal discovery methods based on graphical models.
Frontiers in genetics, 10:0 524, 2019.

[Goudet et al., 2018]
Olivier Goudet, Diviyan Kalainathan, Philippe Caillou, Isabelle Guyon, David
Lopez-Paz, and Michele Sebag.
Learning functional causal models with generative neural networks.
In Explainable and interpretable models in computer vision and
machine learning, pp.  39–80. 2018.

[Hauser & Bühlmann, 2012]
Alain Hauser and Peter Bühlmann.
Characterization and greedy learning of interventional markov
equivalence classes of directed acyclic graphs.
Journal of machine learning research, 13:0 2409–2464,

[He et al., 2021]
Yue He, Peng Cui, Zheyan Shen, Renzhe Xu, Furui Liu, and Yong Jiang.
DARING: differentiable causal discovery with residual independence.
In KDD, 2021.

[Heckerman et al., 1995]
David Heckerman, Dan Geiger, and David M Chickering.
Learning bayesian networks: The combination of knowledge and
statistical data.
Machine learning, 200 (3):0 197–243, 1995.

[Hoyer et al., 2008]
Patrik O. Hoyer, Dominik Janzing, Joris M. Mooij, Jonas Peters, and Bernhard
Nonlinear causal discovery with additive noise models.
In NIPS, 2008.

[Huang et al., 2018]
Biwei Huang, Kun Zhang, Yizhu Lin, Bernhard Schölkopf, and Clark Glymour.
Generalized score functions for causal discovery.
In KDD, 2018.

[Huang et al., 2019]
Biwei Huang, Kun Zhang, Mingming Gong, and Clark Glymour.
Causal discovery and forecasting in nonstationary environments with
state-space models.
In ICML, 2019.

[Huang et al., 2020]
Biwei Huang, Kun Zhang, Jiji Zhang, Joseph D. Ramsey, Ruben Sanchez-Romero,
Clark Glymour, and Bernhard Schölkopf.
Causal discovery from heterogeneous/nonstationary data.
Journal of machine learning research, 21:0
89:1–89:53, 2020.

[Jaakkola et al., 2010]
Tommi S. Jaakkola, David A. Sontag, Amir Globerson, and Marina Meila.
Learning bayesian network structure using LP relaxations.
In AISTATS, 2010.

[Kalainathan et al., 2018]
Diviyan Kalainathan, Olivier Goudet, Isabelle Guyon, David Lopez-Paz, and
Michèle Sebag.
Structural agnostic modeling: Adversarial learning of causal graphs.
arXiv preprint arXiv:1803.04929, 2018.

[Khemakhem et al., 2021]
Ilyes Khemakhem, Ricardo Pio Monti, Robert Leech, and Aapo Hyvärinen.
Causal autoregressive flows.
In AISTATS, 2021.

[Koivisto & Sood, 2004]
Mikko Koivisto and Kismat Sood.
Exact bayesian structure discovery in bayesian networks.
Journal of machine learning research, 5:0 549–573,

[Kuipers et al., 2014]
Jack Kuipers, Giusi Moffa, and David Heckerman.
Addendum on the scoring of gaussian directed acyclic graphical
The Annals of Statistics, 420 (4):0
1689–1691, 2014.

[Kyono et al., 2020]
Trent Kyono, Yao Zhang, and Mihaela van der Schaar.
CASTLE: regularization via auxiliary causal graph discovery.
In NeurIPS, 2020.

[Lachapelle et al., 2020]
Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon
Gradient-based neural DAG learning.
In ICLR, 2020.

[Lippe et al., 2022]
Phillip Lippe, Taco Cohen, and Efstratios Gavves.
Efficient neural causal discovery without acyclicity constraints.
In ICLR, 2022.

[Loh & Bühlmann, 2014]
Po-Ling Loh and Peter Bühlmann.
High-dimensional learning of linear causal networks via inverse
covariance estimation.
Journal of machine learning research, 150
(1):0 3065–3105, 2014.

[Mooij et al., 2020]
Joris M. Mooij, Sara Magliacane, and Tom Claassen.
Joint causal inference from multiple contexts.
Journal of machine learning research, 21:0
99:1–99:108, 2020.

[Ng et al., 2019]
Ignavier Ng, Zhuangyan Fang, Shengyu Zhu, and Zhitang Chen.
Masked gradient-based causal structure learning.
CoRR, abs/1910.08527, 2019a.

[Ng et al., 2019]
Ignavier Ng, Shengyu Zhu, Zhitang Chen, and Zhuangyan Fang.
A graph autoencoder approach to causal structure learning.
CoRR, abs/1911.07420, 2019b.

[Ng et al., 2020]
Ignavier Ng, AmirEmad Ghassami, and Kun Zhang.
On the role of sparsity and DAG constraints for learning linear
In NeurIPS, 2020.

[Ng et al., 2022]
Ignavier Ng, Sébastien Lachapelle, Nan Rosemary Ke, Simon Lacoste-Julien,
and Kun Zhang.
On the convergence of continuous constrained optimization for
structure learning.
In International Conference on Artificial Intelligence and
Statistics, pp.  8176–8198, 2022.

[Pearl, 2000]
Judea Pearl.
Causality: Models, Reasoning, and Inference.

[Pearl et al., 2016]
Judea Pearl, Madelyn Glymour, and Nicholas P Jewell.
Causal inference in statistics: A primer.
John Wiley & Sons, 2016.

[Peters & Bühlmann, 2014]
Jonas Peters and Peter Bühlmann.
Identifiability of gaussian structural equation models with equal
error variances.
Biometrika, 1010 (1):0 219–228, 2014.

[Peters & Bühlmann, 2015]
Jonas Peters and Peter Bühlmann.
Structural intervention distance for evaluating causal graphs.
Neural computation, 270 (3):0 771–799, 2015.

[Peters et al., 2014]
Jonas Peters, Joris M. Mooij, Dominik Janzing, and Bernhard Schölkopf.
Causal discovery with continuous additive noise models.
Journal of machine learning research, 150
(1):0 2009–2053, 2014.

[Peters et al., 2016]
Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen.
Causal inference by using invariant prediction: identification and
confidence intervals.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 780 (5):0 947–1012, 2016.

[Sachs et al., 2005]
Karen Sachs, Omar Perez, Dana Pe'er, Douglas A Lauffenburger, and Garry P
Causal protein-signaling networks derived from multiparameter
single-cell data.
Science, 3080 (5721):0 523–529, 2005.

[Shen et al., 2020]
Xinpeng Shen, Sisi Ma, Prashanthi Vemuri, and Gyorgy Simon.
Challenges and opportunities with causal discovery algorithms:
application to alzheimer’s pathophysiology.
Scientific reports, 100 (1):0 1–12, 2020.

[Shimizu et al., 2006]
Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, and Antti J. Kerminen.
A linear non-gaussian acyclic model for causal discovery.
Journal of machine learning research, 7:0 2003–2030,

[Shrivastava et al., 2016]
Abhinav Shrivastava, Abhinav Gupta, and Ross B. Girshick.
Training region-based object detectors with online hard example
In CVPR, pp.  761–769. IEEE Computer Society, 2016.

[Silander & Myllymäki, 2006]
Tomi Silander and Petri Myllymäki.
A simple approach for finding the globally optimal bayesian network
In UAI, 2006.

[Sokolova et al., 2014]
Elena Sokolova, Perry Groot, Tom Claassen, and Tom Heskes.
Causal discovery from databases with discrete and continuous
In European Workshop on Probabilistic Graphical Models, pp.    442–457, 2014.

[Spirtes & Glymour, 1991]
Peter Spirtes and Clark Glymour.
An algorithm for fast recovery of sparse causal graphs.
Social science computer review, 90 (1):0
62–72, 1991.

[Spirtes & Zhang, 2016]
Peter Spirtes and Kun Zhang.
Causal discovery and inference: concepts and recent methodological
In Applied informatics, volume 3, pp.  1–28, 2016.

[Spirtes et al., 1995]
Peter Spirtes, Christopher Meek, and Thomas S. Richardson.
Causal inference in the presence of latent variables and selection
In UAI, 1995.

[Tsamardinos et al., 2006]
Ioannis Tsamardinos, Laura E Brown, and Constantin F Aliferis.
The max-min hill-climbing bayesian network structure learning
Machine learning, 650 (1):0 31–78, 2006.

[Van de Geer & Bühlmann, 2013]
Sara Van de Geer and Peter Bühlmann.
$l_0$-penalized maximum likelihood for sparse directed acyclic
The Annals of Statistics, 410 (2):0 536–567,

[Varando, 2020]
Gherardo Varando.
Learning dags without imposing acyclicity.
CoRR, abs/2006.03005, 2020.

[Vowels et al., 2021]
Matthew J Vowels, Necati Cihan Camgoz, and Richard Bowden.
D’ya like dags? a survey on structure learning and causal
ACM Computing Surveys (CSUR), 2021.

[Wang et al., 2022]
Yu Wang, An Zhang, Xiang Wang, Xiangnan He, and Tat-Seng Chua.
Differentiable invariant causal discovery.
CoRR, abs/2205.15638, 2022.

[Wang et al., 2020]
Yuhao Wang, Vlado Menkovski, Hao Wang, Xin Du, and Mykola Pechenizkiy.
Causal discovery from incomplete data: A deep learning approach.
CoRR, abs/2001.05343, 2020.

[Wei et al., 2020]
Dennis Wei, Tian Gao, and Yue Yu.
Dags with no fears: A closer look at continuous optimization for
learning bayesian networks.
In NeurIPS, 2020.

[Yang et al., 2021]
Mengyue Yang, Furui Liu, Zhitang Chen, Xinwei Shen, Jianye Hao, and Jun Wang.
Causalvae: Disentangled representation learning via neural structural
causal models.
In CVPR, 2021.

[Yu et al., 2019]
Yue Yu, Jie Chen, Tian Gao, and Mo Yu.
DAG-GNN: DAG structure learning with graph neural networks.
In ICML, 2019.

[Yu et al., 2021]
Yue Yu, Tian Gao, Naiyu Yin, and Qiang Ji.
Dags with no curl: An efficient DAG structure learning approach.
In ICML, 2021.

[Yuan & Malone, 2013]
Changhe Yuan and Brandon M. Malone.
Learning optimal bayesian networks: A shortest path perspective.
Journal of Artificial Intelligence Research, 48:0
23–65, 2013.

[Zhang & Hyvarinen, 2009]
Kun Zhang and Aapo Hyvarinen.
On the identifiability of the post-nonlinear causal model.
UAI, 2009.

[Zheng et al., 2018]
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing.
Dags with NO TEARS: continuous optimization for structure
In NeurIPS, 2018.

[Zheng et al., 2020]
Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing.
Learning sparse nonparametric dags.
In AISTATS, 2020.

[Zhu et al., 2021]
Rong Zhu, Andreas Pfadler, Ziniu Wu, Yuxing Han, Xiaoke Yang, Feng Ye, Zhenping
Qian, Jingren Zhou, and Bin Cui.
Efficient and scalable structure learning for bayesian networks:
Algorithms and applications.
In ICDE, 2021.

[Zhu et al., 2020]
Shengyu Zhu, Ignavier Ng, and Zhitang Chen.
Causal discovery with reinforcement learning.
In ICLR, 2020.

§ RELATED WORK

Differentiable score-based causal discovery methods.
Learning the directed acyclic graph (DAG) from purely observational data is challenging, owing mainly to the intractable combinatorial nature of acyclic graph space.
A recent breakthrough, NOTEARS [Zheng et al., 2018], formulates the discrete DAG constraint into a continuous equality constraint, resulting in a differentiable score-based optimization problem.
Recent subsequent works extends the formulation to deal with nonlinear problems by using a variety of deep learning models, such as neural networks (NOTEARS+ [Zheng et al., 2020], GraN-DAG [Lachapelle et al., 2020], CASTLE [Kyono et al., 2020], MCSL [Ng et al., 2019], DARING [He et al., 2021]), generative autoencoder (CGNN [Goudet et al., 2018], CausalVAE [Yang et al., 2021], ICL [Wang et al., 2020], DAG-GAN [Gao et al., 2021]), graph neural network (DAG-GNN [Gao et al., 2021], GAE [Ng et al., 2019]), generative adversarial network (SAM [Kalainathan et al., 2018], ICL [Wang et al., 2020]), and reinforcement learning (RL-BIC [Zhu et al., 2020]).

Multi-domain causal structure learning.
Most multi-domain causal structure learning methods are constraint-based and have diverse definition of domains.
In our paper, the multi-domain or multi-group refers to heterogeneous data whose underlying causal generating process remain stable but the distributions of noise variables may vary.
In literature, our definition of multi-domain is consistent with MC [Ghassami et al., 2018], CD-NOD [Huang et al., 2020], LRE [Ghassami et al., 2017], DICD [Wang et al., 2022], and others [Peters et al., 2016].
In addition to the strict restriction of knowing the domain annotation in advance, the majority of structure learning models dedicated to heterogeneous data exhibit limited applicability, due to linear case assumption [Ghassami et al., 2018, Ghassami et al., 2017], causal direction identification only [Huang et al., 2019, Cai et al., 2020], and time-consuming [Huang et al., 2020].

§ ALGORITHM OF RESCORE

ReScore Algorithm for Differentiable Score-based Causal Discovery

Input: observational data $\mathcal{D}$: {$\Mat{x}_i: i = 1, 2, ..., n $}, DAG learner parameters $\theta_{\gG}$, reweighting model parameters $\theta_w$, cutoff threshold $\tau$, epoch to start reweighting $K_{reweight}$, maximum epoch in the inner loop $K_{inner}$, maximum epoch in the outer loop $K_{outer}$
Initialize: initialize $\theta_w$ to uniformly output $\frac{1}{n}$, $k_1=0$, $k_2=0$
$k_1 \leq K_{outer} $
Fix reweighting model parameters $\theta_w$
Calculate $\Mat{w}^*$ by applying threshold $[{\tau\over{n}}, {1\over{n\tau}}]$
Optimize $\theta_{\gG}$ by minimizing $S_{\ervw^*}(\gG;\rmX) + \mathcal{P}_{DAG}(\gG)$   #  Outer optimization in maxweight
$k_1 \geq k_{reweight}$
$k_2 \leq K_{inner}$
Fix the DAG learner's parameters $\theta_{\gG}$
Get $\Mat{w}$ from $\theta_w$ by applying threshold $[{\tau\over{n}}, {1\over{n\tau}}]$
Optimize $\theta_{w}$ by maximizing $S_{\ervw}(\gG;\rmX)$   # Inner optimization in maxweight
$k_2 \leftarrow k_2+1$
$k_1 \leftarrow k_1+1$
$k_2 \leftarrow 0$
return predicted $\gG$ from DAG learner

§ IN-DEPTH ANALYSIS OF RESCORE

§.§ Proof of Theorem <ref>

Let $ B = (\beta_1, \dots, \beta_d) \in \mathbb{R}^{d \times d} $ be the weighted adjacent matrix of a SEM, the linear SEM can be written in the matrix form:

\begin{equation} \label{eq:matrixSEM}
X = X B + N
\end{equation}

where $ \mathbb{E}( N | X ) = \overrightarrow{\mathbf{0}} $, $ \Var( N | X ) = diag(\sigma_1^2, \dots, \sigma_d^2) $, and $ B_{ii} = 0 $ since $ X_i $ cannot be the parent of itself. Let $ \Mat{X} \in \mathbb{R}^{n \times d} $ be the observational data and $ \Mat{N} \in \mathbb{R}^{n \times d} $ be the corresponding errors, then

\Mat{X} = \Mat{X} B + \Mat{N}.

The original and reweighted functions for optimization are

\begin{align*}
S(B; \Mat{X}) + \mathcal{P}_{DAG}(B) &= \Ls(B; \Mat{X}) + \lambda \mathcal{R}_{sparse}(B) + \mathcal{P}(B), \\
S_\ervw(B; \Mat{X}) + \mathcal{P}_{DAG}(B) &= \Ls_\ervw(B; \Mat{X}) + \lambda \mathcal{R}_{sparse}(B) + \mathcal{P}_{DAG}(B).
\end{align*}

Comparing the above functions, only the first goodness-of-fit term are different, we will only consider this term.

For the least-squares loss case, the optimization problem is

\begin{align*}
& \min_{B} \Ls_\Mat{w}(B; \rmX) = \min_{B} \sum_{i=1}^{n} w_i l(\Mat{x}_i, \Mat{x}_i B) , \\
& s.t. ~ B_{ii} = 0, \quad i = 1, \dots, d.
\end{align*}

Let $ W = diag(w_1, \dots, w_n) $ be the $n$-dimensional matrix, and rewrite the loss function as

\begin{align*}
\Ls_\Mat{w}(B; \rmX) &= \sum_{i=1}^{n} w_i \Vert \Mat{x}_i - \Mat{x}_i B \Vert_2^2 \\
&= \sum_{i=1}^{n} w_i (\Mat{x}_i - \Mat{x}_i B) (\Mat{x}_i - \Mat{x}_i B)^\top \\
&= \sum_{i=1}^{n} \sum_{j=1}^{d} w_i (\rmX_{ij} - \Mat{x}_i \beta_j)^2 \\
&= \sum_{j=1}^{d} (\Mat{x}^j - \rmX \beta_j)^\top W (\Mat{x}^j - \rmX \beta_j),
\end{align*}

where $ \Mat{x}^j $ is the $j$-th column in matrix $ \rmX $. Let $ D_j $ be the $d$-dimensional identify matrix by setting $j$-th element as 0, for $ j = 1, \dots, d $. The above optimization is able to be written without the restriction:

\begin{align*}
\min_{B} \tilde{\Ls}_\ervw(B; \rmX) &= \min_{B} \sum_{j=1}^{d} (\Mat{x}^j - \rmX D_j \beta_j)^\top W (\Mat{x}^j - \rmX D_j \beta_j) \\
&= \min_{B} \sum_{j=1}^{d} \big( (\Mat{x}^j)^\top W \Mat{x}^j - 2 (\Mat{x}^j)^\top W \rmX D_j \beta_j + \beta_j^\top D_j^\top \rmX^\top W \rmX D_j \beta_j \big).
\end{align*}

The partial derivative of the loss function with respect to $\beta_j$ is

\begin{align*}
\frac{\partial \tilde{\Ls}_\ervw(B; \rmX)}{\partial \beta_j} &= \frac{\partial \bigg[ \sum_{j=1}^{d} \big( (\Mat{x}^j)^\top W X_j - 2 (\Mat{x}^j)^\top W \rmX D_j \beta_j + \beta_j^\top D_j^\top \rmX^\top W \rmX D_j \beta_j \big) \bigg]}{\partial \beta_j} \\
&= \frac{\partial \big( (\Mat{x}^j)^\top W \Mat{x}^j - 2 (\Mat{x}^j)^\top W \rmX D_j \beta_j + \beta_j^\top D_j^\top \rmX^\top W \rmX D_j \beta_j \big)}{\partial \beta_j} \\
&= - 2 D_j^\top \rmX^\top W \Mat{x}^j  + 2 D_j^\top \rmX^\top W \rmX D_j \beta_j.
\end{align*}

Setting the partial derivative to zero produces the optimal parameter:

\begin{align} \label{eq:weightedbeta}
\hat{\beta}_j &= D_j^\top (\rmX^\top W \rmX)^{-1} D_j D_j^\top \rmX^\top W \Mat{x}^j \nonumber \\
&= D_j^\top (\rmX^\top W \rmX)^{-1} D_j D_j^\top \rmX^\top W (\rmX D_j \beta_j + \rmN^j) \nonumber \\
&= D_j \beta_j + D_j (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j,
\end{align}

where $ \rmN^j \in \mathbb{R}^n $ is the $j$-th column in matrix $ \rmN $. In the above equation, the second equality holds because $ \Mat{x}^j = \rmX D_j \beta_j + \rmN^j $. Similarly, one can easily obtain that the optimum parameter for ordinary mean-squared loss is

\begin{align} \label{eq:ordinarybeta}
\tilde{\beta}_j &= D_j \beta_j + D_j (\rmX^\top \rmX)^{-1} \rmX^\top \rmN^j.
\end{align}

It is obvious that the difference between eq:weightedbeta and eq:ordinarybeta is the second term. Compute the mean and variance matrix of the second term in eq:weightedbeta, we can get

\begin{align*}
\mathbb{E} \big[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j \big] &= \mathbb{E} \bigg(\mathbb{E} \big[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j | \rmX \big] \bigg) \\
&= \mathbb{E} \bigg( (\rmX^\top W \rmX)^{-1} \rmX^\top W \cdot \mathbb{E} \big[ \rmN^j | \rmX \big] \bigg) \\
&= \overrightarrow{\mathbf{0}}, \\
% \Var \big[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j \big] &= \mathbb{E} \bigg(\Var \big[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j | \rmX \big] \bigg) \\
% &= \mathbb{E} \bigg( (\rmX^\top W \rmX)^{-1} \rmX^\top W \cdot \Var \big[ \rmN^j | \rmX \big] \cdot W \rmX (\rmX^\top W \rmX)^{-1} \bigg) \\
% &= \sigma_j^2 \mathbb{E} \big[(\rmX^\top W \rmX)^{-1} \rmX^\top W^2 \rmX (\rmX^\top W \rmX)^{-1} \big].
\end{align*}

\begin{align*}
&\Var \big[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j \big] \\
=~& \mathbb{E} \bigg( (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j (\rmN^j)^\top  W^\top \rmX (\rmX^\top W \rmX)^{-1} \bigg) \\
& \qquad \qquad - \bigg( \mathbb{E} \big[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j \big] \bigg) \bigg(\mathbb{E} \big[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j \big]\bigg)^\top\\
=~& \mathbb{E} \bigg( \mathbb{E} \bigg[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \rmN^j (\rmN^j)^\top W \rmX (\rmX^\top W \rmX)^{-1} \big| \rmX \bigg] \bigg) \\
=~& \mathbb{E} \bigg[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \cdot \mathbb{E}\big[ \rmN^j (\rmN^j)^\top | \rmX \big] \cdot W \rmX (\rmX^\top W \rmX)^{-1}  \bigg] \\
=~& \mathbb{E} \bigg[ (\rmX^\top W \rmX)^{-1} \rmX^\top W \cdot \mathbb{E}\big[ \rmN^j (\rmN^j)^\top | \rmX \big] \cdot W \rmX (\rmX^\top W \rmX)^{-1}  \bigg] \\
=~& \sigma_j^2 \mathbb{E} \big[(\rmX^\top W \rmX)^{-1} \rmX^\top W^2 \rmX (\rmX^\top W \rmX)^{-1} \big].
\end{align*}

The last equality holds because $\mathbb{E}(NN^\top|X)=\Var(N|X) + \mathbb{E}(N|X)[\mathbb{E}(N|X)]^\top = diag(\sigma_1^2, \dots, \sigma_d^2) $.

Since $ \Mat{w} \in \sC(\tau) $, it is easy to know that the variance matrix is finite. By the Kolmogorov's strong law of large numbers, the second term converges to zero, thus

\begin{align*}
\hat{\beta}_j \czb{\xrightarrow[]{a.s.}} D_j \beta_j,
\end{align*}

which is same as the ordinary case. Since noise $ N = (N_1, \dots, N_d) $ are jointly independent, the previous process can be apply to the other $ j \in \{1, \dots, d\} $. Let $\hat{B} = (\hat{\beta}_1, \dots, \hat{\beta}_d)$ and $\tilde{B} = (\tilde{\beta}_1, \dots, \tilde{\beta}_d)$, then

\begin{align*}
\hat{B} - \tilde{B} \xrightarrow[]{a.s.} \mathbf{0}.
\end{align*}

Therefore, the convergence has been shown for `case a.'

Since the noise follows a Gaussian distribution, i.e.

$$ X - XB = N = (N_1, \dots, N_d) \sim \mathcal{N} \big( \overrightarrow{\mathbf{0}}, diag(\sigma_1^2, \dots, \sigma_d^2) \big), $$

the loss function (negative log-likelihood function) is

\begin{align} \label{eq:likelihoodloss}
\Ls_\ervw(B; \rmX) &= - \sum_{i=1}^{n} w_i \sum_{j=1}^{d} \bigg[ \log\bigg( \frac{1}{\sigma_j \sqrt{2\pi}} \bigg) - \frac{(\rmX_{ij} - \Mat{x}_i \beta_j)^2}{2 \sigma_j^2} \bigg] \nonumber \\
&= \sum_{j=1}^{d} \sum_{i=1}^{n} w_i \log \big( \sigma_j \sqrt{2\pi} \big) + \sum_{j=1}^{d} \sum_{i=1}^{n} \frac{w_i}{2 \sigma_j^2} (\rmX_{ij} - \Mat{x}_i \beta_j)^2 \nonumber \\
&= \sum_{j=1}^{d} \sum_{i=1}^{n} w_i \log \big( \sigma_j \sqrt{2\pi} \big) + \sum_{j=1}^{d} \frac{1}{2 \sigma_j^2} (\Mat{x}^j - \rmX \beta_j)^\top W (\Mat{x}^j - \rmX \beta_j).
\end{align}

To minimize the loss function above w.r.t. $ B $, it is equivalent to minimize the second term in eq:likelihoodloss:

\min_{B} \Ls_\ervw(B; \rmX) \quad \Longleftrightarrow \quad \min_{B} \sum_{j=1}^{d} \frac{1}{2 \sigma_j^2} (\Mat{x}^j - \rmX \beta_j)^\top W (\Mat{x}^j - \rmX \beta_j).

It can be seen that the RHS above is similar to the loss function in `case a.' except the coefficients $ \frac{1}{2 \sigma_j^2}, j = 1, \dots, d $. Therefore, one can use same approaches to get the equivalence result for `case b.'

Consequently, the proofs of the two special cases have been done.

§.§ Proof of Theorem <ref>

We will show the theorem by contradiction. Without loss of generality, let $i = 1$, $j = 2$, and suppose $w^*_1 \czb{<} w^*_2$. Since $\ervw^* \in \mathbb{C}(\tau)$, one can find a small constant $\varepsilon \in \big( 0, min\{w^*_1-\frac{\tau}{n}, \frac{1}{\tau n}-w^*_2\} \big) $, such that

\begin{align} \label{eq:newweight}
\ervw^{**} = (w^{*}_1 + \varepsilon, w^*_2 - \varepsilon, w^*_3 \dots, w^*_n) \in \czb{\mathbb{C}(\tau)}.
\end{align}

\begin{align*}
&S_{\ervw^*}(\gG;\rmX) - S_{\ervw^{**}}(\gG;\rmX) \\
=~& \big[ w^*_1 \cdot l(\Mat{x}_1, f(\Mat{x}_1)) + w^*_2 \cdot l(\Mat{x}_2, f(\Mat{x}_2)) \big] - \big[(w^*_1+\varepsilon) \cdot l(\Mat{x}_1, f(\Mat{x}_1)) + (w^*_2-\varepsilon) \cdot l(\Mat{x}_2, f(\Mat{x}_2))\big] \\
=~& \varepsilon \cdot [l(\Mat{x}_2, f(\Mat{x}_2)) - l(\Mat{x}_1, f(\Mat{x}_1))] < 0,
\end{align*}

which contradicts $ \Mat{w}^* \in \argmax_{\ervw} \, S_{\ervw}(\gG;\rmX) $. Thus, by contradiction, we can get $w^*_1 \geq w^*_2$ as stated in the theorem.

When $\frac{\tau}{n} < w^*_1 = w^*_2 < \frac{1}{\tau n}$, we can also find a small $\varepsilon \in \big( 0, min\{w^*_1-\frac{\tau}{n}, \frac{1}{\tau n}-w^*_2\} \big)$ such that eq:newweight holds. Similarly, we can get $S_{\ervw^*}(\gG;\rmX) < S_{\ervw^{**}}(\gG;\rmX)$, and $ w^*_1 = w^*_2 = \frac{\tau}{n} $ or $ w^*_1 = w^*_2 = \frac{1}{\tau n} $ by contradiction.

§ SUPPLEMENTARY EXPERIMENTS

§.§ Baselines

We select seven state-of-the-art causal discovery methods as baselines for comparison:

* NOTEARS [Zheng et al., 2018] is a breakthrough work that firstly recasts the combinatorial graph search problem as a continuous optimization problem in linear settings.
NOTEARS estimates the true causal graph by minimizing the reconstruction loss with the continuous acyclicity constraint.

* NOTEARS-MLP [Zheng et al., 2020] is an extension of NOTEARS for nonlinear settings, approximating the generative SEM model by MLP while only applying the continuous acyclicity constraint to the first layer of the MLP.

* GraN-DAG [Lachapelle et al., 2020] adapts the continuous constrained optimization formulation to allow for nonlinear relationships between variables using neural networks and makes use of a final pruning step to remove spurious edges, thus achieving good results in nonlinear settings.

* GOLEM [Ng et al., 2020] improves on the least squares score function [Zheng et al., 2018] by proposing a score function that directly maximizes the data likelihood.
They show the likelihood-based score function with soft sparsity regularization is sufficient to asymptotically learn a DAG equivalent to the ground-truth DAG.

* DICD [Wang et al., 2022] aims to discover the environment-invariant causation while removing the environment-dependent correlation based on ground truth domain annotation.

* CD-NOD [Huang et al., 2020] is a constrained-based causal discovery method that is designed for heterogeneous data, datasets from different environments.
CD-NOD utilizes the independent changes across environments to predict the causal orientations and proposes constrained-based and kernel-based methods to find the causal structure.

* GES [Chickering, 2002] is a score-based search algorithm that searches over the space of equivalence classes of Bayesian network structures.

§.§ Experimental Settings

For NOTEARS, we follow the original linear implementation. For GOLEM, we adopt the GOLEM-NV setting from the original repo. For NOTEARS-MLP, we follow the original non-linear implementation which consist a Multilayer Perceptron (MLP) comprising of two hidden layers with ten neurons each and ReLU activation functions (except for the Sachs dataset, which uses only one hidden layer, inherent the settings from [Zheng et al., 2020]). For GraN-DAG, we employ the pns, training, and cam-pruning stages from the original code and tune three pipeline stages together for best performance. The ReScore adaptive weights learning model for all nonlinear baselines consists of two hidden layer and ReLU activation, and for linear baselines the layer size is reduced to one. All Experiments are conducted on a single Tesla V100 GPU. Detailed hyperparameter search space for different methods is shown in Table <ref>.

Hyperparameter search spaces for each algorithm.

1cHyperparameter space
c]@l@NOTEARS /
NOTEARS+ReScore               c]@l@$\lambda \sim \{0.002,0.005,0.01,0.015,0.02,0.03,0.09,0.1,0.25\}$
Gumbel softmax temperature $\sim \{0.1,1,5,10,20,30,40,50,100\}$
Cut-off threshold $\tau \sim \{0.01,0.1,0.3,0.5,0.7,0.9,0.99 \}$
Constraint convergence tolerance $\sim \{ 10^{-6}, 10^{-8}, 10^{-10}\}$
Log(learning rate of ReScore) $\sim U[-1,-5]$
c]@l@GOLEM /
GOLEM+ReScore                      c]@l@$\lambda \sim \{0.002,0.005,0.01,0.015,0.02,0.03,0.09,0.1,0.25\}$
Gumbel softmax temperature $\sim \{0.1,1,5,10,20,30,40,50,100\}$
Cut-off threshold $\tau \sim \{0.01,0.1,0.3,0.5,0.7,0.9,0.99 \}$
Log(learning rate of ReScore) $\sim U[-1,-5]$
c]@l@NOTEARS-MLP /
NOTEARS-MLP+ReScore             c]@l@$\lambda \sim \{0.002,0.005,0.01,0.015,0.02,0.03,0.09,0.1,0.25\}$
Gumbel softmax temperature $\sim \{0.1,1,5,10,20,30,40,50,100\}$
# hidden units of ReScore $\sim \{1,10,20,50,80,100\}$
# hidden layers of ReScore $\sim \{1,2,3,4\}$
# hidden units of NOTEARS-MLP $\sim \{1,10,20,50,80,100\}$
# hidden layers of ReScore $\sim \{1,2,3\}$
Cut-off threshold $\tau \sim \{0.01,0.1,0.3,0.5,0.7,0.9,0.99 \}$
Constraint convergence tolerance $\sim \{ 10^{-6}, 10^{-8}, 10^{-10}\}$
Log(learning rate of ReScore) $\sim U[-1,-5]$
c]@l@GraN-DAG /
GraN-DAG+ReScore           c]@l@$\lambda \sim \{0.002,0.005,0.01,0.015,0.02,0.03,0.09,0.1,0.25\}$
Gumbel softmax temperature $\sim \{0.1,1,5,10,20,30,40,50,100\}$
# hidden units of ReScore $\sim \{1,10,20,50,80,100\}$
# hidden layers of ReScore $\sim \{1,2,3,4\}$
Log(learning rate of ReScore) $\sim U[-1,-5]$
PNS threshold $\sim \{ 0.5,0.75,1,2\}$
Log(Pruning cutoff) $\sim \{ 0.001,0.005,0.01,0.03,0.1,0.2,0.3\}$

c]@l@GES /
GES+ReScore           c]@l@$\lambda \sim \{0.002,0.005,0.01,0.015,0.02,0.03,0.09,0.1,0.25\}$
Gumbel softmax temperature $\sim \{0.1,1,5,10,20,30,40,50,100\}$
# hidden units of ReScore $\sim \{1,10,20,50,80,100\}$
# hidden layers of ReScore $\sim \{1,2,3,4\}$
Cut-off threshold $\tau \sim \{0.01,0.1,0.3,0.5,0.7,0.9,0.99 \}$

§.§ Study on ReScore

§.§.§ Illustrative Examples of ReScore

To fully comprehend the benefits of reweighting, two research hypotheses need to be verified.
First, we have to determine the validity of the fundamental understanding of ReScore, which states that real-world datasets inevitably include samples of varying importance.
In other words, there are many informative samples that come from disadvantaged groups in real-world scenarios.
Additionally, we must confirm that the adaptive weights learned by ReScore are the faithful reflection of sample importance, less-fitted samples typically come from disadvantaged groups, which are more important than those well-fitted samples.

Real-world Sachs [Sachs et al., 2005] dataset naturally contains nine groups, where each group corresponds with a different experimental condition.
We first rank the importance of each group in Sachs by using the average weights for each group learned by ReScore as the criterion.
Then we eliminate 500 randomly selected samples in one specific group, perform NOTEARS-MLP, and show its DAG accuracy inferred from the remaining samples.
Note that the sample size in each group, which ranges from 700 to 900, is fairly balanced.

Performance comparison for removing samples in different groups

Group Index      3        5      7      1      2      6       4       8       0
Avg. ranking                     578.4     2856.7     3368.1     3877.0     3949.4     4549.4     4573.2     4590.6      4910.1

SHD w/o group     16     16     17     16     16     17     17      19       19

TPR w/o group     0.529     0.412     0.412     0.412     0.412     0.412     0.412      0.353      0.294

Results. Table <ref> clearly shows a declining trend SHD and TPR metrics as the significance of deleting groups grows.
Specifically, removing samples from disadvantaged groups such as Groups 8 and 0, which have the highest average weights, will significantly influence the DAG learning quality.
In contrast, the SHD and TPR of NOTEARS-MLP can even be maintained or slightly decreased by excluding the samples from groups with relatively low average weights.
This illustrates that samples of different importance are naturally present in real-world datasets, and ReScore is capable of successfully extracting this importance.

§.§.§ Sensitivity to Pure Noise Samples

A basic assumption of ReScore is that no pure noise outliers are involved in the training process.
Otherwise, the DAG learner might get overwhelmed by arbitrarily up-weighting less well-fitted samples, in this case, pure noise data.
The good news is that the constraint of the cutoff threshold $\tau \in \sC(\tau) = \{\Mat{w}: 0 < \frac{\tau}{n} \leq w_1, \dots, w_n \leq \frac{1}{\tau n}, \sum_{i=1}^{n} w_i = 1\}$ prevents over-exploitation of pure noise samples, which further strengthens ReScore's ability to withstand outliers.
To evaluate the robustness of ReScore against pure noise samples, the following experiments are conducted.

We produce $p_{corrupt}$ percentage pure noise samples in nonlinear settings ($n=2000$, $d = 20$, ER2), where those noise samples are generated from a different structural causal model.
We try out a broad range of $p_{corrupt} = \{0, 0.01, 0.02, 0.05, 0.08, 0.1, 0.2, 0.3, 0.5\}$.

SHD for $p_{corrupt}$ percentage noise samples.

0        0.01     0.02     0.05     0.08     0.1      0.2      0.3     0.5
NOTEARS-MLP                       14.9     15.2     15.3     18.9     19.8     21.3     23.9     23.8     28.3

+ ReScore ($\tau \rightarrow 0$)     13.8     14.2     15.0     18.3     19.5     20.7     24.0      24.4       29.3

+ ReScore (Optimal $\tau$)             13.7     14.1     15.0     18.1     19.2     19.9     21.9       24.0      28.9

Imp. %                           +8%      +7%      +2%      +4%      +3%      +7%      +8%       -1%       -2%

Results. Table <ref> reports the comparison of performance in NOTEARS-MLP and two ReScore methods (no cut-off threshold and optimal $\tau$) when encountering pure noise data.
The best-performing methods are bold; Imp.% measures the relative improvements of ReScore (Optimal $\tau$) over the backbone NOTEARS-MLP.
We observe that ReScore (Optimal $\tau$) consistently yields remarkable improvements compared with NOTEARS-MLP in the case that less than 20% of samples are corrupted.
These results demonstrate the robustness of ReScore when handling data that contains a small proportion of pure noise data.
Surprisingly, when the cutoff threshold $\tau$ is set to be close to 0, the ReScore can still achieve relative gains over the baseline when less than 10% of the samples are pure noise.
Although it is more sensitive to noise samples than the optimum cutoff threshold $\tau$.
These surprising findings support the effectiveness of adaptive weights and show the potential of ReScore.

§.§.§ Effect of Hyperparameter $\tau$.

We investigate the effect of cut-off threshold $\tau$ on the performance of ReScore.
Intuitively, ReScore relies on the hyperparameter $\tau$ to control the balance between hard sample mining and robustness towards extremely noisy samples.
On one hand, setting the threshold closer to 0 results in no weight-clipping and leaves the model susceptible to noises, which results in sub-optimal performance.
On the other hand, setting the threshold closer to 1 disables the reweighting scheme and eventually reduces ReScore performance to its backbone model.

We conduct experiments under different settings of $\tau$  using $n=2000$ samples generated from GP model on ER4 graphs with $d=20$ nodes.
The weight distribution under best performing threshold $tau=0.9$ and the trend of SHD to $\tau$ is shown in Figure <ref>.
One can observe that ReScore obtains its best performance at $\tau=0.9$, while a smaller or bigger threshold results in sub-optimal performance.
Furthermore, we find that in different settings, the optimal threshold $\tau$ usually falls in the range of $[0.7,0.99]$. This indicates that ReScore performs best when adaptive reweighting is conducted within a restricted range.

weight distribution with $\tau=0.9$

SHD threshold $\tau$

Study of varying $\tau$ in ReScore model.

§.§.§ Sensitivity to Neural Network Complexity.

We also investigated the effect of number of hidden units in our adaptive weights learning model for ReScore. We plot the TPR, FDR, SHD, SID with varing number of hidden units ranging from 10 to 100 units in nonlinear settings, using $n=600$ and $n=2,000$ samples generated from GP model on ER4 graph with $d=10$ nodes.
Detailed results could be found in Figure <ref>.
One can first observe our model is stable when increasing the neurons, illustrating the insensitivity of ReScore the number of neurons in adaptive weights learning model.
On the other hand, more observational samples to estimate the parameters could help the ReScore achieve higher performance, indicating rich samples bring benefit.

TPR neurons

FDR neurons

SHD neurons

SID neurons

Performance with varying neurons in ReScore model.

§.§.§ Training Costs.

In terms of time complexity, as shown in Table <ref>, we report the time for each baseline and ReScore on Sachs. Compared with backbone methods, ReScore adds very little computing cost to training.

Training cost on Sachs (seconds per iteration/in total).

NOTEARS      0.74 / 2.97

+ ReScore     3.8 / 15.3

NOTEARS-MLP      0.87 / 3.48

+ ReScore     4.3 / 17.0

GOLEM      13.4 / 53.5

+ ReScore     14.6 / 58.2

GraN-DAG      4.9 / 197.3

+ ReScore     5.5 / 221.6

§.§ More Experimental Results for RQ1

More experimental results on both the linear and nonlinear synthetic data are reported in Figures <ref> - <ref> and Tables <ref> - <ref>.
The error bars depict the standard deviation across datasets over ten trails.
The red and blue percentages separately refer to the increase and decrease of ReScore relative to the original score-based causal discovery methods in each metric.
The best performing methods per task are bold.

TPR $\lambda$

FDR $\lambda$

SHD $\lambda$

SID $\lambda$

Performance comparison between NOTEARS-MLP and ReScore on ER2 graphs of 10 nodes on nonlinear synthetic datasets. The hyperparameter $\lambda$ defined in eq:score-comb refers to the graph sparsity.

TPR $\lambda$

FDR $\lambda$

SHD $\lambda$

SID $\lambda$

Performance comparison between NOTEARS-MLP and ReScore on ER2 graphs of 20 nodes on nonlinear synthetic datasets. The hyperparameter $\lambda$ defined in eq:score-comb refers to the graph sparsity.

TPR $\lambda$

FDR $\lambda$

SHD $\lambda$

SID $\lambda$

Performance comparison between NOTEARS-MLP and ReScore on ER4 graphs of 20 nodes on nonlinear synthetic datasets. The hyperparameter $\lambda$ defined in eq:score-comb refers to the graph sparsity.

Results for ER graphs of 20 nodes on linear and nonlinear synthetic datasets.

1cTPR  $\mathbf{\uparrow}$     1cFDR $\mathbf{\downarrow}$     1cSHD $\mathbf{\downarrow}$     1c|SID $\mathbf{\downarrow}$     1cTPR  $\mathbf{\uparrow}$     FDR $\mathbf{\downarrow}$     SHD $\mathbf{\downarrow}$     SID $\mathbf{\downarrow}$
Random          0.11$\pm$0.09     0.89$\pm$0.08    56.8$\pm$8.7     292.3$\pm$45.7     0.07$\pm$0.03     0.90$\pm$0.08     86.9$\pm$7.0     387.5$\pm$52.3

NOTEARS         0.85$\pm$0.08
0.09$\pm$0.03     9.2$\pm$3.8      55.4$\pm$31.1     0.74$\pm$0.02     0.23$\pm$0.03      39.4$\pm$7.9      185.8$\pm$38.1

+ ReScore       $\textbf{0.87}\text{\scriptsize{$\pm$0.07}}^{\color{+}+2\%}$       $0.11\text{\scriptsize{$\pm$0.05}}^{\color{-}-17 \%}$      $\textbf{8.8}\text{\scriptsize{$\pm$3.5}}^{\color{+}+5 \%}$     $\textbf{50.6}\text{\scriptsize{$\pm$26.3}}^{\color{+}+9 \%}$      $\textbf{0.79}\text{\scriptsize{$\pm$0.05}}^{\color{+}+7\%}$     $0.28\text{\scriptsize{$\pm$0.05}}^{\color{-}-17\%}$      $\textbf{36.8}\text{\scriptsize{$\pm$7.9}}^{\color{+}+7 \%}$      $\textbf{180.8}\text{\scriptsize{$\pm$43.5}}^{\color{+}+3 \%} $

GOLEM           0.75$\pm$0.07
0.20$\pm$0.11     17.0$\pm$6.1      78.2$\pm$22.6     0.46$\pm$0.06     0.50$\pm$0.06      73.6$\pm$7.9      249.8$\pm$7.8

+ ReScore       $0.76\text{\scriptsize{$\pm$0.06}}^{\color{+}+2\%}$       $0.20\text{\scriptsize{$\pm$0.10}}^{\color{+}+1 \%}$      $15.8\text{\scriptsize{$\pm$5.8}}^{\color{+}+8 \%}$     $77.0\text{\scriptsize{$\pm$21.5}}^{\color{+}+2 \%}$      $0.48\text{\scriptsize{$\pm$0.06}}^{\color{+}+3\%}$     $0.43\text{\scriptsize{$\pm$0.06}}^{\color{+}+16\%}$      $70.2\text{\scriptsize{$\pm$8.3}}^{\color{+}+5 \%}$      $246.2\text{\scriptsize{$\pm$11.4}}^{\color{+}+1 \%} $
NOTEARS-MLP     0.70$\pm$0.12     0.13$\pm$0.07     14.9$\pm$5.4     98.4$\pm$22.5     0.44$\pm$0.09     0.26$\pm$0.10     55.0$\pm$9.2     176.3$\pm$33.3

+ ReScore       $0.73\text{\scriptsize{$\pm$0.09}}^{\color{+}+3\%}$     $0.11\text{\scriptsize{$\pm$0.05}}^{\color{+}+7 \%}$     $13.7\text{\scriptsize{$\pm$5.1}}^{\color{+}+8 \%}$     $88.8\text{\scriptsize{$\pm$23.3}}^{\color{+}+11 \%}$     $0.41\text{\scriptsize{$\pm$0.07}}^{\color{-}-6 \%}$     $0.17\text{\scriptsize{$\pm$0.08}}^{\color{+}+54 \%}$     $\textbf{51.6}\text{\scriptsize{$\pm$6.4}}^{\color{+}+7 \%}$     $179.9\text{\scriptsize{$\pm$33.7}}^{\color{-}-2 \%}$

GraN-DAG        0.81$\pm$0.15     0.08$\pm$0.08     9.3$\pm$5.4     53.4$\pm$24.4     0.20$\pm$0.07     0.18$\pm$0.08     57.4$\pm$4.6     131.5$\pm$21.4

+ ReScore       $\textbf{0.81}\text{\scriptsize{$\pm$0.14}}^{\color{+}+0\%}$     $\textbf{0.05}\text{\scriptsize{$\pm$0.04}}^{\color{+}+64\%}$     $\textbf{8.5}\text{\scriptsize{$\pm$5.7}}^{\color{+}+9\%}$     $\textbf{51.0}\text{\scriptsize{$\pm$24.6}}^{\color{+}+5\%}$     $0.21\text{\scriptsize{$\pm$0.07}}^{\color{+}+5\%}$     $\textbf{0.17}\text{\scriptsize{$\pm$0.09}}^{\color{+}+8\%}$     $56.2\text{\scriptsize{$\pm$4.6}}^{\color{+}+2\%}$     $\textbf{125.4}\text{\scriptsize{$\pm$23.3}}^{\color{+}+5\%}$

Results for ER graphs of 50 nodes on linear and nonlinear synthetic datasets.

1cTPR  $\mathbf{\uparrow}$     1cFDR $\mathbf{\downarrow}$     1cSHD $\mathbf{\downarrow}$     1c|SID $\mathbf{\downarrow}$     1cTPR  $\mathbf{\uparrow}$     FDR $\mathbf{\downarrow}$     SHD $\mathbf{\downarrow}$     SID $\mathbf{\downarrow}$
Random          0.04$\pm$0.02     0.90$\pm$0.03     397.3$\pm$12.7     1082.0$\pm$182.2     0.09$\pm$0.08     0.92$\pm$0.08     998.2$\pm$45.9     3399.1$\pm$489.2

NOTEARS         0.79$\pm$0.06     0.09$\pm$0.03     27.6$\pm$7.7     427.0$\pm$186.1     0.51$\pm$0.12     0.27$\pm$0.10     133.4$\pm$29.5     1643.8$\pm$172.2

+ ReScore       $\textbf{0.88}\text{\scriptsize{$\pm$0.06}}^{\color{+}+11\%}$     $0.15\text{\scriptsize{$\pm$0.04}}^{\color{-}-39\%}$     $\textbf{26.2}\text{\scriptsize{$\pm$7.6}}^{\color{+}+5\%}$     $\textbf{266.0}\text{\scriptsize{$\pm$146.4}}^{\color{+}+61\%}$     $\textbf{0.52}\text{\scriptsize{$\pm$0.21}}^{\color{+}+3\%}$     $0.29\text{\scriptsize{$\pm$0.07}}^{\color{-}-7\%}$     $\textbf{130.2}\text{\scriptsize{$\pm$37.4}}^{\color{+}+2\%}$     $\textbf{1453.6}\text{\scriptsize{$\pm$336.5}}^{\color{+}+13\%}$

GOLEM           0.80$\pm$0.09
0.35$\pm$0.09     68.6$\pm$19.7      433.5$\pm$215.6     0.31$\pm$0.11     0.68$\pm$0.06      150.6$\pm$25.1       1775.4$\pm$161.6

+ ReScore       $0.82\text{\scriptsize{$\pm$0.15}}^{\color{+}+3\%}$       $0.33\text{\scriptsize{$\pm$0.14}}^{\color{+}+5 \%}$      $63.4\text{\scriptsize{$\pm$27.9}}^{\color{+}+8\%}$     $430.2\text{\scriptsize{$\pm$155.5}}^{\color{+}+1 \%}$      $0.39\text{\scriptsize{$\pm$0.06}}^{\color{+}+24\%}$     $0.66\text{\scriptsize{$\pm$0.06}}^{\color{+}+3\%}$      $146.3\text{\scriptsize{$\pm$26.3}}^{\color{+}+3 \%}$      $1643.6\text{\scriptsize{$\pm$114.8}}^{\color{+}+8 \%} $
NOTEARS-MLP     0.32$\pm$0.04     0.13$\pm$0.08     69.5$\pm$4.7     884.4$\pm$172.8     0.17$\pm$0.02     0.06$\pm$0.04     167.0$\pm$4.1     1607.6$\pm$97.0

+ ReScore       $0.51\text{\scriptsize{$\pm$0.08}}^{\color{+}+59\%}$     $\textbf{0.10}\text{\scriptsize{$\pm$0.07}}^{\color{+}+30\%}$     $53.5\text{\scriptsize{$\pm$8.7}}^{\color{+}+30\%}$     $628.1\text{\scriptsize{$\pm$120.6}}^{\color{+}+41\%}$     $0.26\text{\scriptsize{$\pm$0.04}}^{\color{+}+52\%}$     $0.11\text{\scriptsize{$\pm$0.05}}^{\color{-}-51\%}$     $154.4\text{\scriptsize{$\pm$6.4}}^{\color{+}+8\%}$     $1437.7\text{\scriptsize{$\pm$111.1}}^{\color{+}+12\%}$

GraN-DAG        0.52$\pm$0.09     0.15$\pm$0.05     51.6$\pm$9.3     632.8$\pm$140.3     0.32$\pm$0.04     0.08$\pm$0.16     141.6$\pm$8.2     1379.0$\pm$91.3

+ ReScore       $\textbf{0.53}\text{\scriptsize{$\pm$0.06}}^{\color{+}+3\%}$     $0.11\text{\scriptsize{$\pm$0.02}}^{\color{+}+36\%}$     $\textbf{46.0}\text{\scriptsize{$\pm$6.0}}^{\color{+}+12\%}$     $\textbf{581.0}\text{\scriptsize{$\pm$104.7}}^{\color{+}+9\%}$     $0.31\text{\scriptsize{$\pm$0.03}}^{\color{-}-4\%}$     $0.06\text{\scriptsize{$\pm$0.04}}^{\color{+}+32\%}$     $\textbf{138.8}\text{\scriptsize{$\pm$7.5}}^{\color{+}+2\%}$     $\textbf{1351.0}\text{\scriptsize{$\pm$98.2}}^{\color{+}+2\%}$

Results for SF graphs of 10 nodes on linear and nonlinear synthetic datasets.

1cTPR  $\mathbf{\uparrow}$     1cFDR $\mathbf{\downarrow}$     1cSHD $\mathbf{\downarrow}$     1c|SID $\mathbf{\downarrow}$     1cTPR  $\mathbf{\uparrow}$     FDR $\mathbf{\downarrow}$     SHD $\mathbf{\downarrow}$     SID $\mathbf{\downarrow}$
Random          0.05$\pm$0.03     0.91$\pm$0.09     32.2$\pm$7.97     35.1$\pm$7.3     0.13$\pm$0.01     0.93$\pm$0.15     57.2$\pm$10.3     79.1$\pm$8.7

NOTEARS         0.98$\pm$0.02     0.02$\pm$0.03     0.8$\pm$0.5     1.0$\pm$2.0     0.95$\pm$0.03     0.03$\pm$0.02     12.2$\pm$1.2     6.2$\pm$5.3

+ ReScore       $\textbf{0.99}\text{\scriptsize{$\pm$0.02}}^{\color{+}+1\%}$     $0.04\text{\scriptsize{$\pm$0.04}}^{\color{-}-45\%}$     $\textbf{0.4}\text{\scriptsize{$\pm$0.7}}^{\color{+}+100\%}$     $\textbf{1.0}\text{\scriptsize{$\pm$0.9}}^{\color{+}+0\%}$     $\textbf{0.97}\text{\scriptsize{$\pm$0.03}}^{\color{+}+2\%}$     $0.03\text{\scriptsize{$\pm$0.03}}^{\color{+}+27\%}$     $10.2\text{\scriptsize{$\pm$1.5}}^{\color{+}+20\%}$     $\textbf{3.0}\text{\scriptsize{$\pm$1.9}}^{\color{+}+107\%}$

GOLEM           0.96$\pm$0.07     0.07$\pm$0.12     1.8$\pm$3.1     1.2$\pm$2.4     0.85$\pm$0.03     0.12$\pm$0.08     7.0$\pm$2.3     12.8$\pm$7.9

+ ReScore       $0.97\text{\scriptsize{$\pm$0.07}}^{\color{+}+1\%}$     $0.07\text{\scriptsize{$\pm$0.12}}^{\color{+}+3\%}$     $1.4\text{\scriptsize{$\pm$2.9}}^{\color{+}+29\%}$     $1.2\text{\scriptsize{$\pm$2.4}}^{\color{+}+0\%}$     $0.87\text{\scriptsize{$\pm$0.06}}^{\color{+}+3\%}$     $\textbf{0.10}\text{\scriptsize{$\pm$0.08}}^{\color{+}+17\%}$     $\textbf{5.8}\text{\scriptsize{$\pm$2.9}}^{\color{+}+21\%}$     $9.8\text{\scriptsize{$\pm$8.2}}^{\color{+}+31\%}$
NOTEARS-MLP     0.84$\pm$0.17     0.25$\pm$0.12     6.7$\pm$3.4     8.1$\pm$7.3     0.73$\pm$0.14     0.23$\pm$0.05     12.0$\pm$3.9     19.4$\pm$7.4

+ ReScore       $0.82\text{\scriptsize{$\pm$0.22}}^{\color{-}-2\%}$     $0.17\text{\scriptsize{$\pm$0.08}}^{\color{+}+45\%}$     $5.8\text{\scriptsize{$\pm$3.3}}^{\color{+}+16\%}$     $\textbf{6.0}\text{\scriptsize{$\pm$3.8}}^{\color{+}+35\%}$     $\textbf{0.88}\text{\scriptsize{$\pm$0.09}}^{\color{+}+20\%}$     $0.27\text{\scriptsize{$\pm$0.07}}^{\color{-}-16\%}$     $11.0\text{\scriptsize{$\pm$3.4}}^{\color{+}+9\%}$     $12.8\text{\scriptsize{$\pm$9.3}}^{\color{+}+52\%}$

GraN-DAG        0.69$\pm$0.20     0.05$\pm$0.05     5.9$\pm$3.0     12.0$\pm$8.2     0.82$\pm$0.11     0.11$\pm$0.08     8.7$\pm$1.8     8.4$\pm$4.1

+ ReScore       $0.72\text{\scriptsize{$\pm$0.17}}^{\color{+}+4\%}$     $\textbf{0.04}\text{\scriptsize{$\pm$0.03}}^{\color{+}+28\%}$     $\textbf{5.3}\text{\scriptsize{$\pm$2.8}}^{\color{+}+11\%}$     $10.5\text{\scriptsize{$\pm$8.7}}^{\color{+}+14\%}$     $0.86\text{\scriptsize{$\pm$0.12}}^{\color{+}+5\%}$     $0.12\text{\scriptsize{$\pm$0.08}}^{\color{-}-12\%}$     $\textbf{8.1}\text{\scriptsize{$\pm$2.0}}^{\color{+}+7\%}$     $\textbf{7.0}\text{\scriptsize{$\pm$6.7}}^{\color{+}+20\%}$

Results for SF graphs of 20 nodes on linear and nonlinear synthetic datasets.

1cTPR  $\mathbf{\uparrow}$     1cFDR $\mathbf{\downarrow}$     1cSHD $\mathbf{\downarrow}$     1c|SID $\mathbf{\downarrow}$     1cTPR  $\mathbf{\uparrow}$     FDR $\mathbf{\downarrow}$     SHD $\mathbf{\downarrow}$     SID $\mathbf{\downarrow}$
Random          0.11$\pm$0.10     0.89$\pm$0.03     43.2$\pm$5.4     96.8$\pm$10.4     0.09$\pm$0.05     0.88$\pm$0.05     108.2$\pm$12.9     155.6$\pm$37.2

NOTEARS         0.90$\pm$0.06     0.02$\pm$0.01     4.0$\pm$1.9     19.8$\pm$12.8     0.90$\pm$0.05     0.12$\pm$0.06     45.2$\pm$7.0     28.6$\pm$20.2

+ ReScore       $0.95\text{\scriptsize{$\pm$0.04}}^{\color{+}+6\%}$     $0.06\text{\scriptsize{$\pm$0.04}}^{\color{-}-70\%}$     $\textbf{3.6}\text{\scriptsize{$\pm$1.8}}^{\color{+}+11\%}$     $\textbf{9.8}\text{\scriptsize{$\pm$8.1}}^{\color{+}+102\%}$     $\textbf{0.93}\text{\scriptsize{$\pm$0.03}}^{\color{+}+3\%}$     $\textbf{0.02}\text{\scriptsize{$\pm$0.07}}^{\color{+}+624\%}$     $45.0\text{\scriptsize{$\pm$6.8}}^{\color{+}+0\%}$     $\textbf{25.6}\text{\scriptsize{$\pm$12.1}}^{\color{+}+12\%}$

GOLEM           0.96$\pm$0.03     0.19$\pm$0.06     9.0$\pm$3.2     10.4$\pm$7.0     0.83$\pm$0.05     0.35$\pm$0.09     42.8$\pm$13.0     41.4$\pm$14.8

+ ReScore       $\textbf{0.96}\text{\scriptsize{$\pm$0.02}}^{\color{+}+0\%}$     $0.18\text{\scriptsize{$\pm$0.06}}^{\color{+}+4\%}$     $8.6\text{\scriptsize{$\pm$3.1}}^{\color{+}+5\%}$     $10.4\text{\scriptsize{$\pm$7.0}}^{\color{+}+0\%}$     $0.85\text{\scriptsize{$\pm$0.43}}^{\color{+}+2\%}$     $0.34\text{\scriptsize{$\pm$0.09}}^{\color{+}+5\%}$     $\textbf{39.8}\text{\scriptsize{$\pm$14.0}}^{\color{+}+8\%}$     $37.6\text{\scriptsize{$\pm$12.8}}^{\color{+}+10\%}$
NOTEARS-MLP     0.42$\pm$0.13     0.23$\pm$0.13     25.5$\pm$4.5     49.9$\pm$7.4     0.20$\pm$0.03     0.22$\pm$0.12     58.9$\pm$3.1     115.6$\pm$25.0

+ ReScore       $0.41\text{\scriptsize{$\pm$0.13}}^{\color{-}-2\%}$     $\textbf{0.10}\text{\scriptsize{$\pm$0.10}}^{\color{+}+121\%}$     $\textbf{23.5}\text{\scriptsize{$\pm$4.5}}^{\color{+}+9\%}$     $\textbf{47.6}\text{\scriptsize{$\pm$9.4}}^{\color{+}+5\%}$     $\textbf{0.21}\text{\scriptsize{$\pm$0.04}}^{\color{+}+3\%}$     $\textbf{0.09}\text{\scriptsize{$\pm$0.09}}^{\color{+}+131\%}$     $\textbf{56.4}\text{\scriptsize{$\pm$2.2}}^{\color{+}+4\%}$     $\textbf{109.0}\text{\scriptsize{$\pm$21.8}}^{\color{+}+6\%}$

GraN-DAG        0.03$\pm$0.15     0.24$\pm$0.17     27.1$\pm$4.15     77.0$\pm$28.0     0.20$\pm$0.06     0.18$\pm$0.12     56.8$\pm$4.5     133.4$\pm$21.0

+ ReScore       $0.03\text{\scriptsize{$\pm$0.15}}^{\color{-}-6\%}$     $0.15\text{\scriptsize{$\pm$0.10}}^{\color{+}+63\%}$     $25.7\text{\scriptsize{$\pm$4.4}}^{\color{+}+5\%}$     $72.8\text{\scriptsize{$\pm$26.0}}^{\color{+}+6\%}$     $0.21\text{\scriptsize{$\pm$0.07}}^{\color{+}+2\%}$     $0.17\text{\scriptsize{$\pm$0.08}}^{\color{+}+8\%}$     $\textbf{56.4}\text{\scriptsize{$\pm$4.6}}^{\color{+}+1\%}$     $125.4\text{\scriptsize{$\pm$23.3}}^{\color{+}+6\%}$

Results for SF graphs of 50 nodes on linear and nonlinear synthetic datasets.

1cTPR  $\mathbf{\uparrow}$     1cFDR $\mathbf{\downarrow}$     1cSHD $\mathbf{\downarrow}$     1c|SID $\mathbf{\downarrow}$     1cTPR  $\mathbf{\uparrow}$     FDR $\mathbf{\downarrow}$     SHD $\mathbf{\downarrow}$     SID $\mathbf{\downarrow}$
Random          0.10$\pm$0.08     0.89$\pm$0.07     334.2$\pm$16.9     1093.3$\pm$145.4     0.12$\pm$0.11     0.89$\pm$0.04     1023.5$\pm$49.5     1903.9$\pm$194.3

NOTEARS         0.82$\pm$0.03     0.07$\pm$0.05     23.6$\pm$6.2     135.4$\pm$47.5     0.71$\pm$0.18     0.25$\pm$0.07     97.6$\pm$36.9     276.2$\pm$131.0

+ ReScore       $\textbf{0.94}\text{\scriptsize{$\pm$0.03}}^{\color{+}+15\%}$     $0.15\text{\scriptsize{$\pm$0.06}}^{\color{-}-55\%}$     $\textbf{21.6}\text{\scriptsize{$\pm$9.1}}^{\color{+}+9\%}$     $\textbf{61.2}\text{\scriptsize{$\pm$22.6}}^{\color{+}+121\%}$     $\textbf{0.73}\text{\scriptsize{$\pm$0.05}}^{\color{+}+3\%}$     $\textbf{0.10}\text{\scriptsize{$\pm$0.03}}^{\color{+}+138\%}$     $\textbf{67.6}\text{\scriptsize{$\pm$12.3}}^{\color{+}+44\%}$     $\textbf{275.2}\text{\scriptsize{$\pm$55.0}}^{\color{+}+0\%}$

GOLEM           0.77$\pm$0.07     0.19$\pm$0.11     38.6$\pm$16.7     161.6$\pm$53.2     0.62$\pm$0.17     0.21$\pm$0.09     114.2$\pm$37.5     384.0$\pm$107.4

+ ReScore       $0.79\text{\scriptsize{$\pm$0.09}}^{\color{+}+2\%}$     $0.24\text{\scriptsize{$\pm$0.12}}^{\color{-}-20\%}$     $32.2\text{\scriptsize{$\pm$11.1}}^{\color{+}+20\%}$     $143.4\text{\scriptsize{$\pm$63.0}}^{\color{+}+13\%}$     $0.68\text{\scriptsize{$\pm$0.17}}^{\color{+}+9\%}$     $0.21\text{\scriptsize{$\pm$0.09}}^{\color{+}+1\%}$     $113.7\text{\scriptsize{$\pm$37.5}}^{\color{+}+0\%}$     $366.4\text{\scriptsize{$\pm$107.0}}^{\color{+}+5\%}$
NOTEARS-MLP     0.22$\pm$0.04     0.04$\pm$0.04     75.8$\pm$4.0     266.8$\pm$46.0     0.11$\pm$0.02     0.03$\pm$0.02     168.8$\pm$3.8     461.6$\pm$54.9

+ ReScore       $\textbf{0.23}\text{\scriptsize{$\pm$0.05}}^{\color{+}+4\%}$     $0.07\text{\scriptsize{$\pm$0.07}}^{\color{-}-47\%}$     $\textbf{75.6}\text{\scriptsize{$\pm$4.3}}^{\color{+}+0\%}$     $267.2\text{\scriptsize{$\pm$36.6}}^{\color{-}-0\%}$     $\textbf{0.13}\text{\scriptsize{$\pm$0.04}}^{\color{+}+10\%}$     $0.07\text{\scriptsize{$\pm$0.06}}^{\color{-}-52\%}$     $\textbf{167.7}\text{\scriptsize{$\pm$7.0}}^{\color{+}+1\%}$     $\textbf{453.4}\text{\scriptsize{$\pm$57.7}}^{\color{+}+2\%}$

GraN-DAG        0.19$\pm$0.03     0.28$\pm$0.05     80.2$\pm$3.5     380.8$\pm$56.1     0.11$\pm$0.03     0.25$\pm$0.11     171.4$\pm$6.3     549.6$\pm$84.9

+ ReScore       $0.20\text{\scriptsize{$\pm$0.03}}^{\color{+}+5\%}$     $\textbf{0.24}\text{\scriptsize{$\pm$0.05}}^{\color{+}+17\%}$     $79.8\text{\scriptsize{$\pm$0.3}}^{\color{+}+1\%}$     $349.2\text{\scriptsize{$\pm$49.6}}^{\color{+}+9\%}$     $0.11\text{\scriptsize{$\pm$0.02}}^{\color{+}+0\%}$     $0.24\text{\scriptsize{$\pm$0.10}}^{\color{+}+5\%}$     $170.8\text{\scriptsize{$\pm$4.0}}^{\color{+}+0\%}$     $548.0\text{\scriptsize{$\pm$91.4}}^{\color{+}+0\%}$