# Robustness of Algorithms for Causal Structure Learning to Hyperparameter
Choice

Damian Machlanski<EMAIL_ADDRESS>
Department of Computer Science and Electronic Engineering
University of Essex Spyridon Samothrakis<EMAIL_ADDRESS>
Institute for Analytics and Data Science
University of Essex Paul Clarke<EMAIL_ADDRESS>
Institute for Social and Economic Research
University of Essex

###### Abstract

Hyperparameters play a critical role in machine learning. Hyperparameter
tuning can make the difference between state-of-the-art and poor prediction
performance for any algorithm, but it is particularly challenging for
structure learning due to its unsupervised nature. As a result, hyperparameter
tuning is often neglected in favour of using the default values provided by a
particular implementation of an algorithm. While there have been numerous
studies on performance evaluation of causal discovery algorithms, how
hyperparameters affect individual algorithms, as well as the choice of the
best algorithm for a specific problem, has not been studied in depth before.
This work addresses this gap by investigating the influence of hyperparameters
on causal structure learning tasks. Specifically, we perform an empirical
evaluation of hyperparameter selection for some seminal learning algorithms on
datasets of varying levels of complexity. We find that, while the choice of
algorithm remains crucial to obtaining state-of-the-art performance,
hyperparameter selection in ensemble settings strongly influences the choice
of algorithm, in that a poor choice of hyperparameters can lead to analysts
using algorithms which do not give state-of-the-art performance for their
data.

Keywords: Hyperparameters, model selection, causal discovery, structure
learning, performance evaluation, misspecification, robustness.

##  1 Introduction

Uncovering causal graphs is an immensely useful tool in data-driven decision-
making as it helps understand the underlying data generating process. A large
number of causal structure learning algorithms incorporate Machine Learning
(ML) methods. These, in turn, heavily rely on hyperparameters for accurate
predictions (Bergstra et al., 2011). In addition, there has been growing
evidence that correctly specified hyperparameters can close the performance
gap between State-of-the-Art (SotA) and other methods (Paine et al., 2020;
Zhang et al., 2021a; Machlanski et al., 2023; Tönshoff et al., 2023). Are
hyperparameters as important in structure recovery?

Hyperparameter optimisation is extremely challenging in structure learning as
the true graphs are inaccessible outside of simulated environments. This
inability to reliably tune could be one of the reasons behind the struggle to
apply some of the algorithms to real data problems (Kaiser & Sipos, 2021), or
why hyperparameters are often completely neglected in this area. On the one
hand, benchmarks and evaluation frameworks (e.g. Raghu et al. (2018); Tu et
al. (2019)) usually focus on finding a learning algorithm that works best
under specific circumstances but without considering hyperparameters as part
of the problem. On the other hand, studies that address hyperparameter tuning
(e.g. Strobl (2021); Biza et al. (2022)) consider individual algorithms but
not the impact of tuning (or the lack of it) on selecting the best algorithm
for the available data. Understanding how hyperparameters affect algorithm
choice, as well as individual methods, is clearly missing but can be a crucial
next step towards more stable causal discovery in real data applications. To
make matters worse, the evaluation metrics used for tuning can be imperfect
and sometimes favour specific learning methods (Curth & van der Schaar, 2023).
This brings us to the core questions of this paper: Do different algorithms
perform similarly given access to a hyperparameter oracle? How robust are they
against misspecified hyperparameters?

In this work, we set out to address these questions and investigate the impact
hyperparameters have on graph recovery performance of individual algorithms,
as well as on the best algorithm choice (see Figure 1). We start by showing
how a single hyperparameter plays a crucial role in the simplest graph problem
(two variables). More extensive experiments strengthen this observation and
confirm it as a more general phenomenon. The experimental setup involves many
seminal structure learning algorithms tested against real and simulated
datasets.

Contributions. This paper offers the following contributions:

  * •

Compare algorithms’ performances and their winning percentages across
hyperparameters.

  * •

Compare algorithms’ performances under well-specified and misspecified
hyperparameters.

  * •

Compare algorithms’ winning percentages under well-specified and misspecified
hyperparameters.

Related work. This work connects with the existing literature mainly through
the topics of performance evaluation and hyperparameter analysis. The
performance of structure learning algorithms has been evaluated from a number
of different perspectives, such as mixed data types (Raghu et al., 2018) or
time series data (Assaad et al., 2022). In an attempt to strengthen the
evaluation, there have been efforts to develop testing environments that
closely resemble real-life datasets. Some examples include simulators based on
gene regulatory networks (Van den Bulcke et al., 2006) or neuropathic pain
pathology (Tu et al., 2019). Grünbaum et al. (2023) take evaluation further by
proposing to test algorithms on the parts of real-life datasets that are known
a priori. Furthermore, to improve reproducibility, Rios et al. (2021)
developed a benchmarking platform that covers a wide range of learning methods
and data scenarios. The importance of hyperparameters and their impact on
performance have been mostly studied in other areas outside of structure
learning. In the offline Reinforcement Learning (RL) setting, Paine et al.
(2020) reported, among other aspects, that robustness to hyperparameter
choices is an important issue and that careful tuning can deliver close to
optimal policies. Zhang et al. (2021a), on the other hand, make a case for
hyperparameter tuning in model-based RL. Furthermore, according to Machlanski
et al. (2023) and Tönshoff et al. (2023), hyperparameters alone can be
responsible for reaching, or even surpassing, SotA performance levels.
Hyperparameters in structure learning have mostly been discussed in the
context of tuning. One common approach is to select hyperparameters that
result in stable structure predictions across random data samples (Liu et al.,
2010; Sun et al., 2013; Strobl, 2021). Another strand of work performs out-of-
sample validation for tuning purposes based on predictive accuracy of models
fitted in accordance with the recovered graph structure (Biza et al., 2022),
or assigned scores developed specifically for structure recovery tuning
(Chobtham & Constantinou, 2023). Metrics based on regression error have been
also considered (Marx & Vreeken, 2019), though in the context of two
variables.

Structure. In Section 2 we briefly discuss the basics of structure learning.
Section 3 demonstrates the importance of hyperparameters in structure learning
via a bivariate example, further motivating more extensive numerical
experiments presented in Section 4. Section 5 concludes the paper and offers
potential future work directions.

Figure 1: Summary of the main idea of the paper. We explore various structure
learning algorithms and investigate how hyperparameters affect their
performance. Notation: $h_{1}$ and $h_{2}$ are different hyperparameter
values; $X$ denotes i.i.d. data provided to algorithms; $\hat{A}$ is recovered
adjacency matrix while $\mathcal{G}(\hat{A})$ is a causal graph based on
$\hat{A}$; SHD is structural Hamming distance (lower is better). Note how
recovered graphs differ between different hyperparameters of the same
algorithm (green edges are correct; red incorrect). sim_mean are
hyperparameters that achieved the best average performance across all
simulations.

##  2 Structure Learning

We briefly describe the notation and data assumptions used throughout the
paper as well as important details of structure learning methods necessary to
read the technical parts of the document. For a more detailed review, see
recent literature (e.g. Eberhardt (2017); Glymour et al. (2019)).

###  2.1 Graphs

Let $\mathcal{G}=(V,\mathcal{E})$ be a graph with nodes/vertices
$V=\\{1,\ldots,p\\}$ and edges $\mathcal{E}\subseteq V^{2}$. Edges are pairs
of nodes $(j,k)\in\mathbf{V}$ where $(v,v)\notin\mathcal{E}$ to exclude self-
cycles. Nodes $j,k$ are adjacent in $\mathcal{G}$ if either
$(j,k)\in\mathcal{E}$ or $(k,j)\in\mathcal{E}$. An edge is undirected if
$(j,k)\in\mathcal{E}$ and $(k,j)\in\mathcal{E}$, whereas it is directed if
only one pair appears in $\mathcal{E}$; if this pair is $(j,k)$ then $j$ is
called a child of parent $k$. The set of parents of $j$ in $\mathcal{G}$ is
denoted by $\textbf{PA}^{\mathcal{G}}_{j}$. We call $\mathcal{G}$ undirected
if all its edges are undirected; conversely, $\mathcal{G}$ consisting only of
directed edges is directed. A mixed graph consists of both directed and
undirected edges. The skeleton of any directed or mixed graph $\mathcal{G}$ is
an equivalent graph with all directed edges replaced by undirected ones. A
fully connected graph $\mathcal{G}$ is one where all pairs of nodes are
adjacent. A (directed) path is a sequence of nodes connected by (directed)
edges. A partially directed acyclic graph (PDAG) is a mixed graph such that
there is no pair $(j,k)$ such that there are directed paths from $j$ to $k$
and vice versa. Then, $\mathcal{G}$ is a directed acyclic graph (DAG) if it is
a PDAG and is directed. Two graphs are Markov equivalent, or belong to the
same equivalence class, when they involve the same sets of d-separations
(Pearl, 2000). A completed PDAG (CPDAG) can encode such a class of graphs, in
which undirected edges mean that the graphs within the class may contain a
directed edge in either direction; directed edges denote agreement in edge
direction in subsumed graphs.

###  2.2 Assumptions

Now consider a vector of random variables $\mathbf{X}=(X_{1},\ldots,X_{p})$
generated according to an unknown data generating process (DGP) leading to
joint distribution $\mathcal{L}(\mathbf{X})$. The node $j\in\mathbf{V}$
represents random variable $X_{j}$ and the edge between nodes $j$ and $k$ in
$\mathcal{E}$ is directed if and only if $X_{k}$ is used in the DGP to
generate $X_{j}$. We further assume the following.

######  Assumption 1 (Sufficiency).

There are no hidden confounders.

######  Assumption 2 (Markov condition).

Two variables are independent in $\mathcal{L}(\mathbf{X})$ if they are
d-separated in $\mathcal{G}$.

######  Assumption 3 (Faithfulness).

Two variables are d-separated in $\mathcal{G}$ if they are independent in
$\mathcal{L}(\mathbf{X})$.

###  2.3 Learning Methods

The goal of causal structure learning is to infer (or identify) graph
$\mathcal{G}$ given the distribution $\mathcal{L}(\mathbf{X})$. If it is
possible to do this, we say $\mathcal{G}$ is identifiable from $\mathcal{L}$.
Traditional methods, such as PC (Spirtes & Glymour, 1991) or GES (Chickering,
2002), were often built around assumptions a-c) above, which are in most cases
not enough to identify a unique DAG solution, only the class of CPDAGs. If
one, however, assumes the DGP is a Structural Causal Model (SCM) then
identification is possible, for example, if the DAG is a linear SCM with
additive non-Gaussian noise:
$X_{j}=f_{j}(X_{\textbf{PA}^{\mathcal{G}}_{j}})+\epsilon_{j}$ (Shimizu et al.,
2006). Subsequent research has extended this result to nonlinear SCMs with
additive noise (Hoyer et al., 2008), linear models with Gaussian noise terms
of equal variances (Peters & Bühlmann, 2014), and additive models of the form
$X_{j}=\sum_{k\in\textbf{PA}^{\mathcal{G}}_{j}}f(X_{k})+\epsilon_{j}$
(Bühlmann et al., 2014). More recent approaches also involve neural network-
based algorithms that specifically restate the learning task as a continuous
optimisation problem (Zheng et al., 2018; 2020).

##  3 Hyperparameters in Structure Learning

###  3.1 Bivariate Example

An illustrative example, strongly inspired by Marx & Vreeken (2019), is the
classic cause-effect pairs challenge (Guyon et al., 2019) that consists of two
(synthetically generated here) variables $X$ and $Y$, with the goal of
establishing the existence and direction of the causal link between them
($X\rightarrow Y$, $X\leftarrow Y$, no link) given only observed data. One
possible solution is to fit two regressors $y=f(x)$ and $x=g(y)$, and predict
the causal direction based on the lower prediction error of the two models
($\epsilon_{f}=[y-\hat{f}(x)]^{2}$ and $\epsilon_{g}=[x-\hat{g}(y)]^{2}$; no
link if the errors are comparable). As shown in Figure 2, changing the
hyperparameter that controls the number of regression parameters can result in
a different causal direction being predicted. This is precisely what
constitutes the problem, since the true DGP and the correct hyperparameter
value are unknown.

######  Observation 1.

Incorrect hyperparameters can cause prediction mistakes.

######  Observation 2.

There might be more than one correct and incorrect hyperparameter choice.

Figure 2: The number of allowed regression parameters, a hyperparameter,
clearly affects prediction error of the two models and can determine predicted
causal direction (see decision colours at the bottom). The algorithm predicts
$X\rightarrow Y$ if $\epsilon_{f}<\epsilon_{g}$, $X\leftarrow Y$ if
$\epsilon_{f}>\epsilon_{g}$, and is inconclusive otherwise. Note that the true
causal direction is unknown in practice.

The problem grows in complexity as the number of graph nodes and edges
increases. This is because each edge is a potentially different function to
approximate that will require a different hyperparameter value (function
complexity) to obtain the correct answer. In addition, many algorithms provide
multiple hyperparameters to tune, making even more room for further mistakes,
effectively increasing the chance of hyperparameter misspecification. In fact,
even the bivariate example can involve more hyperparameters by, for instance,
introducing a threshold such that the algorithm predicts ‘no link’ if
$|\epsilon_{f}-\epsilon_{g}|<threshold$. To this end, we make the following
claim and set up the following key definition:

######  Claim 1.

The existence and direction of an edge in the predicted graph strongly depends
on algorithm’s hyperparameters.

######  Definition 1 (Hyperparameter Misspecification).

Mistakes in predicted graph structure arising from incorrect hyperparameters.

###  3.2 General Form of the Problem

Let $\mathbf{X}\in\mathbb{R}^{n\times p}$ represent i.i.d. data of $n$
observations and $p$ features. Furthermore, let $A\in\\{0,1\\}^{p\times p}$ be
a binary adjacency matrix of directed graph $\mathcal{G}(A)$ such that
$a_{jk}=1$ if $(j,k)\in\mathcal{E}$ and $a_{jk}=0$ otherwise. A weighted
adjacency matrix $W\in\mathbb{R}^{p\times p}$ is defined such that
$A(W)_{jk}=1$ if $w_{jk}\neq 0$ and zero otherwise, which results in a
weighted graph $\mathcal{G}(W)$. Let us also define distance $d(A,B)$ between
two adjacency matrices $A$ and $B$ such that $d(A,B)=0$ if and only if
$\mathcal{G}(A)=\mathcal{G}(B)$ are the same graph. From now on, let us denote
$A$ as the true adjacency matrix and $\hat{A}$ its estimate obtained by a
programme $P$ from i.i.d. data $X$ using programme options $O$ so that

$\hat{A}=P(O,X),$ (1)

where $O$ generally involves the user specifying algorithm $S$ and
hyperparameter values $H$. Therefore, given $K$ user-specified candidate
programmes $P=\\{P_{1},\ldots,P_{K}\\}$ and

$\hat{A}_{k}=P(S_{k},H_{k},X),$ (2)

where $S_{k}$ and $H_{k}$ are the algorithm and hyperparameter choices
associated with program candidate $k$, then the best program is

$k^{*}=\operatorname*{argmin}_{k\in\\{1,\ldots,K\\}}d(A,\hat{A}_{k}),$ (3)

Note that $k^{*}$ is generally not identifiable unless $A$ is known.
Furthermore, identification of $\mathcal{G}$ does not guarantee
$\hat{A}_{k^{*}}=A$ as $\hat{A}_{k^{*}}$ depends on algorithms $S$ and their
ability to identify $A$, as well as the choice of their hyperparameters $H$.

More generally, when considering different algorithms $S$ and hyperparameters
$H$, Equation 3 is the standard model selection problem, whereas if the choice
of algorithms is fixed to a specific value, leaving hyperparameters as the
only variable, the task reduces to hyperparameter tuning. In practical terms,
obtaining the distance $d(A,\hat{A})$ is not feasible, as the true matrix $A$
and its corresponding graph $\mathcal{G}(A)$ are inaccessible outside of
simulated environments. As algorithms can vary substantially in design, the
appropriate way to compare them requires the use of distance measures that
incorporate the ground truth $A$. This renders model selection impractical in
structure learning problems. Tuning hyperparameters of a single algorithm
might be feasible by comparing its relative scores across explored
hyperparameter values.

###  3.3 Common Hyperparameters

Despite differences in algorithms, many of them share similar hyperparameters.
Commonly used ones are briefly described here.

Significance level of independence tests. Refers to the p-value of
independence tests and the desired level of certainty. Decreasing the value
(increasing certainty) will usually result in fewer predicted edges. Often
named as alpha and incorporated by traditional and pairwise algorithms.

Sparsity. A penalty term that encourages sparser solutions. Higher values
result in fewer predicted edges. Similar in mechanism to L1 regularisation
which discards less relevant features. Often employed by regression-based
solutions, especially if they perform some form of feature selection.

Model complexity. A penalty that encourages simpler models to avoid
overfitting (L2 regularisation). As shown in the example in Section 3.1, its
influence on the final prediction is complicated. Usually applies to solutions
that model the assumed form of SEMs.

Post-prunning threshold. Many SEM-based methods output weighted adjacency
matrices $W$ that need to be converted to the binary form of $A$. This is
usually done by applying a threshold below which all edges are set to zero.
That is, $a_{ij}=1$ if $w_{ij}>w\\_thresh$; $0$ otherwise..

Note that alpha and w_threshold are algorithm agnostic and can be transferred
between methods, whereas the other two hyperparameters may differ in value
between algorithms.

##  4 Experiments

Since our analysis in Section 3 is based merely on a simple and artificial
example, our next step is to study the influence of hyperparameters more
rigorously in a more general setting. We devise a set of experiments
consisting of diverse data sets of various size and difficulty (Section 4.1),
processed by a representative set of structure learning algorithms (Section
4.2). Different hyperparameter selection scenarios are also detailed (Section
4.4).

The experimental framework is implemented through Benchpress (Rios et al.,
2021), a benchmarking platform to evaluate structure learning algorithms.
Performances of all algorithms are collected from Benchpress, followed by some
mild post-processing of the results to suit our analysis of hyperparameters.
All numerical experiments can be fully replicated using the code and data that
are available online at: https://github.com/dmachlanski/benchpress-dm.

###  4.1 Graphs and Data

####  4.1.1 Simulations

We follow recent literature in structure learning when it comes to simulating
different DGPs. The simulation procedure starts with generating a random DAG
$\mathcal{G}$ with $p$ nodes and $d$ edges, built according to a random graph
model. The resulting graph is binary, with $A\in\\{0,1\\}^{p\times p}$. Next,
i.i.d. data $\mathrm{X}\in\mathbb{R}^{n\times p}$ are sampled from a simulated
SEM of choice, with $n$ being the sample size. Each individual combination of
settings is repeated for $10$ seeds and forms a single experiment. In our
experiments, we explore $p\in\\{10,20,50\\}$, $d\in\\{1p,4p\\}$, and
$n\in\\{200,1000\\}$. Included random graph models are Erdös-Rényi (ER) (Erdös
& Rényi, 1959) and Barabási-Albert (Barabási & Albert, 1999), with the latter
also known as scale-free (SF). We also explore $n=10,000$ but only for sparse
ER graphs with $p=50$ nodes due to computational limitations. As for explored
SEMs, we include the following:

Linear Non-Gaussian (gumbel). $X=XW^{T}+z\in\mathbb{R}^{p}$, with
$W\in\mathbb{R}^{p\times p}$ as edge weights assigned independently from
$U(\left[-2,-0.5\right]\cup\left[0.5,2\right])$ and based on $A$. Noise $z$
follows the Gumbel distribution $z\sim$ Gumbel($0,I_{p\times p}$).

Nonlinear Gaussian (gp).
$X_{j}=f_{j}(X_{\textbf{PA}^{\mathcal{G}}_{j}})+z_{j}$ for all
$j\in\left[p\right]$ in the topological order of $\mathcal{G}$. Noise $z_{j}$
follows Gaussian distribution $z_{j}\sim\mathcal{N}(0,1)$, $j=1,\ldots,p$.
Where functions $f_{j}$ represent a draw from a Gaussian process with a unit
bandwidth RBF kernel.

Note that both settings have been shown to be identifiable. That is, linear
non-Gaussian additive models (Shimizu et al., 2006) and nonlinear additive
models (Hoyer et al., 2008).

####  4.1.2 Real Datasets

We also tested structure learning algorithms against real or semi-real
datasets. The most popular ones in the literature are protein signaling and
SynTReN.

Protein signaling comes from Sachs et al. (2005) which measures protein and
phospholipid expression levels in human cells. The ground truth causal graph
has been established and accepted by the experts in the field. We use the
second dataset that is already logged and standardised and consists of $n=902$
observations, $p=11$ nodes and $d=17$ edges.

SynTReN is a generator of synthetic transcriptional regulatory networks and
related gene expression data that simulate a real experiment (Van den Bulcke
et al.,
2006)111http://bioinformatics.intec.ugent.be/kmarchal/SynTReN/index.html. We
use the same data as in Lachapelle et al. (2019), which consist of $10$ random
seeds, $n=500$ samples and $p=20$ nodes.

###  4.2 Structure Learning Algorithms

We consider in our setup the following algorithms. Due to high computational
demands, we only focus on well-established and seminal algorithms that, in our
view, effectively represent different classes of solutions. More details about
the hyperparameters involved can be found in Appendix A.1 and A.2.

  * •

PC (Spirtes & Glymour, 1991). Peter and Clark algorithm. Constraint-based
approach that starts with a fully-connected undirected graph and removes edges
based on conditional independence tests. Next, it attempts to orient as many
of the remaining edges as possible. The result is a CPDAG.
Hyperparameters: alpha (significance level for conditional independence
tests).

  * •

FCI (Spirtes et al., 1993). Fast Causal Inference. Constraint-based. An
important generalisation of PC to unknown confounding variables.
Hyperparameters: alpha (significance level for conditional independence
tests).

  * •

FGES (Ramsey et al., 2017). Fast Greedy Equivalence Search. Optimised and
parallelised version of the original score-based GES algorithm (Chickering,
2002). It starts with an empty graph and adds an edge that yields maximum
score improvement until no significant score gain is achieved. Then it removes
edges in the same greedy manner until a plateau.
Hyperparameters: penaltyDiscount (sparsity penalty).

  * •

LiNGAM (Shimizu et al., 2006). Linear Non-Gaussian Acyclic Model. Assumes
linear SEMs and non-Gaussian noise that enters additively:
$X_{j}=\sum_{k\in\textbf{PA}^{\mathcal{G}}_{j}}w_{jk}X_{k}+\epsilon_{j}$.
Hyperparameters: max_iter (FastICA (Hyvarinen, 1999)), thresh (post-prunning
threshold).

  * •

ANM (Hoyer et al., 2008). Additive Noise Model. Assumes nonlinear SEMs and
additive noise: $X_{j}=f_{j}(X_{\textbf{PA}^{\mathcal{G}}_{j}})+\epsilon_{j}$.
Hyperparameters: alpha (significance level for the independence test).

  * •

CAM (Bühlmann et al., 2014). Causal Additive Models. Assumes a generalised
additive noise model with additive noise and functions:
$X_{j}=\sum_{k\in\textbf{PA}^{\mathcal{G}}_{j}}f(X_{k})+\epsilon_{j}$.
Hyperparameters: cutoff (variable selection threshold).

  * •

NOTEARS (Zheng et al., 2018). Score-based continuous DAG optimisation with a
smooth acyclicity regularisation term. Assumes linear SEMs with additive
noise.
Hyperparameters: lambda1 (sparsity term), max_iter (optimisation steps) and
w_threshold (post-prunning threshold).

  * •

NOTEARS MLP (Zheng et al., 2020). Nonlinear extension of NOTEARS by
incorporating the Multi-Layer Perceptron (MLP). Assumes nonlinear SEMs with
additive noise.
Hyperparameters: lambda1 (sparsity term), lambda2 (regularisation strength),
w_threshold (post-prunning threshold), hidden_units (number of units in the
hidden layer).

Many traditional algorithms, such as PC, FCI and FGES, make the standard set
of assumptions that involve sufficiency, faithfulness and Markov condition.
These, however, are often not enough to identify a unique DAG as a solution,
which is a major drawback of these methods (they output CPDAGs). Making
assumptions about distributions and functional forms of the data generating
process seems to be critical to overcome this issue (all methods above except
for PC, FCI and FGES output DAGs).

###  4.3 Evaluation

In order to compare algorithms’ performances, we employ the commonly used
structural Hamming distance (SHD) metric, which is provided via Benchpress
(Rios et al., 2021, appendix A.1.). For the convenience of the reader, we
briefly describe it here as well.

Let us define $E$ and $E^{\prime}$ as a set of edges of the true and predicted
DAG respectively. Then, for $e\in E^{\prime}$, true positives (TP) and false
positives (FP) are assigned as follows:

$TP(e)=\begin{cases}1&\text{ if }e\in E\text{ and correctly oriented}\\\
0.5&\text{ if }e\in E\text{ and incorrectly oriented}\\\ 0&\text{ otherwise
}\end{cases}$ (4)

$FP(e)=\begin{cases}1&\text{ if }e\notin E\\\ 0.5&\text{ if }e\in E\text{ and
incorrectly oriented}\\\ 0&\text{ otherwise }\end{cases}$ (5)

where TP and FP are sums of all TP(e) and FP(e) scores respectively. The
structural Hamming distance (SHD) aggregates the number additions, removals
and reversals in predicted edges so they match the true ones ($E=E^{\prime}$).
It can be defined as:

$SHD=|E|-TP+FP$ (6)

Note the SHD defined as above allows to evaluate mixed graphs, that is,
compare DAGs to CPDAGs. If, for instance, a predicted undirected edge exists
in $E$ but is supposed to be directed, it will result in $TP=0.5$ and
$FP=0.5$, ultimately leading to $SHD=1$. This shows that the need to orient an
undirected edge is treated equally as the need to add, remove or reverse an
edge so $E=E^{\prime}$. Such evaluation puts algorithms outputting CPDAGs at a
disadvantage compared to DAG-only methods. We justify it on the grounds that
the main focus of this study is DAG recovery, hence any predicted undirected
edge is treated as any other mistake. The ability to evaluate mixed graphs is
an important feature for this study as it allows us to compare algorithms
outputting CPDAGs and DAGs.

###  4.4 Hyperparameters

All incorporated learning algorithms have at least one hyperparameter. We
collect performances of algorithms across all hyperparameter combinations
(exhaustive grid search; see Appendix A.1). To better understand the influence
of hyperparameters on structure recovery performance, we experiment with four
different hyperparameter selection strategies described below.

BEST. To simulate the choice of the best hyperparameters (as if we had access
to hyperparameter oracle), we pick hyperparameter values that achieved the
lowest SHD in that particular data setting. Each data setting can have a
different set of the best hyperparameters.

WORST. Identified similarly to ‘best’ except the criterion here is the highest
SHD.

DEFAULT. Default hyperparameter values recommended by the authors of an
algorithm. See Appendix A.2 for specific values.

SIM_MEAN. An alternative to ‘default’. We found a single set of hyperparameter
values per algorithm that achieved the lowest average SHD across all
simulations. In other words, these are simulation-derived default values. See
Appendix A.1 for identified values.

###  4.5 Results

Presented results employ the following naming convention with respect to the
data generating process: graph_p (number of nodes), graph_d (edge density),
graph_type (graph models; ER or SF), data_n (sample size), data_sem (SEM type;
gumbel or gp). Error bars are standard errors unless stated otherwise. Note
only the most important results are presented in the main content of the
paper. The rest of the results, that do not change conclusions, are in
Appendix B.

####  4.5.1 Performance vs. Hyperparameter Quality

Figure 3: SHD performances depending on the quality of selected
hyperparameters. Numbers are averaged across data seeds as well as sample
sizes and graph types due to negligible differences.

As shown in Figure 3, achieved performance is clearly affected by both
algorithm and hyperparameter choices. Even assuming access to hyperparameter
oracle (‘best’), selecting different algorithms will have a significant impact
on the result. Furthermore, fixed hyperparameters seem to be a viable strategy
as they are relatively close in SHD to the best cases. The differences between
simulation-derived and paper-default values (‘sim_mean’ and ‘default’
respectively) are negligible in most cases. The worst hyperparameters, on the
other hand, can result in performances substantially worse than the fixed
ones. This shows the risks of hyperparameter misspecification, the degree of
which clearly varies across algorithms (different robustness). SHDs against
dense graphs are generally high even with the best hyperparameters, suggesting
that those cases should be avoided with current tools. For the sake of
clarity, consider the upper left subfigure ($graph\\_p=10$, $graph\\_d=1$,
$data\\_sem=gumbel$) as an example. First, we can see that SHD varies
significantly across algorithms even if they have access to a hyperparameter
oracle as the blue bars range from around 2 SHD (FGES) to as much as over 15
(ANM). Next, we can observe that the orange and green bars that refer to
‘default’ or ‘sim_mean’ hyperparameters are in most cases not that far from
blue bars (‘best’ hyperparameters) within the same algorithm, showing that a
fixed set of hyperparameter values can work surprisingly well. Looking at the
red bars that represent the ‘worst’ hyperparameter values, it is clear that
SHDs associated with them are much higher than with any other bar colour. This
shows how misspecified hyperparameters can degrade algorithm’s performance.
The differences between best and worst cases (blue and red respectively),
which we read as robustness to misspecification, can also vary between
algorithms. For instance, ANM is fairly robust due to its small worst-best
difference in SHD, whereas NOTEARS, despite its low SHD in the best case, has
a substantially higher SHD in the worst case (high worst-best difference as a
result), showing its poor robustness to misspecification.

####  4.5.2 Performance Distribution Across Hyperparameters

As per Figure 7, all algorithms perform similarly when averaged across all
simulations and hyperparameters. But as simulation circumstances change,
algorithms respond differently. Specifically, not only do the best and worst
hyperparameter cases vary across algorithms, but also the proportion of
relatively good performances. A high proportion of low SHD values, for
example, could be read as a proxy for robustness to hyperparameter
misspecification. Figure 9 confirms that increased sample size generally
helps, even with relatively large graphs ($p=50$), although some algorithms
require more data to notice significant benefits (see LiNGAM in gumbel and
NOTEARS_MLP in gp). Positive effects can be noticed with respect to improved
best hyperparameter cases and an increased proportion of good performances.
This case also confirms that relatively big and sparse graphs can be recovered
with high accuracy given the right hyperparameters.

(a) all simulations

(b) data_n

(a) data_sem

(b) graph_d

(c) graph_p

(d) graph_type

Figure 7: Distributions of SHD performances across all hyperparameters.

Figure 8: across all hyperparameters Figure 9: Performances across all
hyperparameters for ER graphs with $p=50$ and $d=1$. ANM was excluded due to
long execution time against $10,000$ samples.

####  4.5.3 Winning Algorithms vs. Hyperparameter Quality

Previous analysis revealed that the best algorithm choice may depend on
specific DGP properties as well as hyperparameter choices. To make it clearer,
we collect winning algorithms across different DGP properties and
hyperparameter selection strategies. An algorithm with the lowest SHD in a
given setting wins. The cumulative winning percentages are presented in Figure
10. The results confirm that no one algorithm wins in all settings (i.e, there
are no free lunches). Specific DGP properties may favour certain methods. When
looking at how winning odds change depending on different hyperparameter
choices, it is clear that the best algorithm selection depends not only on DGP
properties, but also on the type of available hyperparameters. For instance,
in order to minimise the consequences of potential hyperparameter
misspecification, it might be best to choose the algorithm that wins most
frequently in a specific setting with the worst hyperparameters. In addition,
it is clear that sample size and the type of the simulated graph do not impact
winning odds across hyperparameters in a significant way.

####  4.5.4 Semi-Synthetic and Real Data

We put our simulation-derived findings to a test by performing structure
recovery on SynTReN (Figure 11) and Sachs (Figure 12) datasets. All numbers
are compared to SotA performances retrieved from Lachapelle et al. (2019),
which are $33.7\pm 3.7$ and $12$ SHD for SynTReN and Sachs respectively. Both
cases generally confirm that fixed hyperparameters (sim_mean and default) can
work almost as well as the best hyperparameters, and that even the best
hyperparameters may not be enough to reach the best possible performance as
some algorithms perform better than other under those conditions. It is also
clear from both cases that hyperparameters play an important role and, in
fact, can decide whether an algorithm reaches or beats SotA. For instance,
against SynTReN, both NOTEARS methods and LiNGAM seem to be good options under
the best and fixed hyperparameters. But under the worst hyperparameters,
NOTEARS methods can be extremely inaccurate, making ANM the safest choice in
this case. In the Sachs dataset, this is no longer the case with ANM, showing
that the best algorithm pick indeed strongly depends on DGP properties. All
algorithms except ANM can, in fact, beat SotA on Sachs data. However, when it
comes to robustness to hyperparameter misspecification and safety of use,
NOTEARS methods appear to be the most risky, with LiNGAM being extraordinarily
robust, as it beats SotA even with its worst hyperparameters.

(a) all simulations

(b) $graph\\_p$

(c) $data\\_n$

(d) $data\\_sem$

(e) $graph\\_d$

(f) $graph\\_type$

Figure 10: Winning percentage depending on hyperparameter quality and
simulation properties.

##  5 Conclusion

In this work, we have successfully shown that hyperparameters play an
important role in causal structure learning. However, the way hyperparameters
influence the methods is somewhat different than recent results from the ML
literature. More specifically, Machlanski et al. (2023); Tönshoff et al.
(2023) found that many learners can reach SotA performance levels with the
right hyperparameters, reducing the importance of model selection. But in this
study, we observe that algorithms still differ significantly in performance
even with access to a hyperparameter oracle. However, reliable tuning is not
always available in structure learning, leading to hyperparameter
misspecification and prediction errors. This is where hyperparameters become
important as we showed that different learning algorithms vary in robustness
to hyperparameter misspecification, and that strong performance under the
right hyperparameters does not necessarily translate to misspecification
robustness. As a consequence, an algorithm that is the best pick under correct
hyperparameters, might be a suboptimal choice when its hyperparameters are
misspecified; another algorithm with better misspecification robustness might
be safer to use, especially in those cases where minimising the consequences
of potential misspecification is a priority. Thus, overall, the best algorithm
choice may depend not only on the properties of the data generating process,
but also on the quality of selected hyperparameter values. In terms of
secondary findings, default hyperparameters seem to perform surprisingly well
in many cases, and hence may constitute a viable alternative to tuning, though
the risk of potential misspecification should be carefully considered when
doing so. Another interesting observation is that relatively large sparse
graphs ($50$ nodes) can be recovered with high accuracy, subject to large
sample size and the right hyperparameters.

(a) quality of selected hyperparameters

(b) across all hyperparameters

Figure 11: Performances against the SynTReN dataset. Numbers are averaged
across data seeds.

(a) quality of selected hyperparameters

(b) across all hyperparameters

Figure 12: Performances against the Sachs dataset.

###  5.1 Limitations

This study is naturally limited by our choice of explored algorithms,
hyperparameters and simulation properties. It is worth noting, however, that
we do not intend to identify the best possible learning algorithm or
hyperparameters. On the contrary, the objective of this work is to show that
the appropriate algorithm choices are nuanced, as also recently shown in the
treatment effect estimation domain (Curth & van der Schaar, 2023), and that
hyperparameters should be part of that subtle decision-making process, further
confirming the importance of hyperparameters reported in the literature (Paine
et al., 2020; Zhang et al., 2021a; Machlanski et al., 2023; Tönshoff et al.,
2023). More extensive search spaces are unlikely to negate such conclusions.

###  5.2 Future work

Provided our observation that increasing the sample size generally increases
the proportion of good performances across all hyperparameters of an
algorithm, a natural next step would be to examine whether the hyperparameter
misspecification problem vanishes as the sample size approaches infinity. As
hyperparameter tuning is one possible strategy of mitigating misspecification
with finite sample sizes, further work into reliable tuning schemes might help
stabilise performance and facilitate better adoption of algorithms in real-
world datasets. Another issue is the barrier of using tuning metrics, as
obtaining them often requires algorithm modifications that cannot be done
without expert knowledge of the method. Making the metrics easier to use could
facilitate further research on the topic, as well as provide means for good
tuning practices. Finally, choosing an algorithm appropriate for the problem
at hand has a considerable impact on the final performance. Although some
hyperparameter tuning metrics are general enough to perform algorithm
selection as well (e.g. Liu et al. (2010); Biza et al. (2022)), doing so on
real-life datasets is still a challenge. One promising direction could be
validation that incorporates domain knowledge (Grünbaum et al., 2023).

#### Acknowledgments

All three authors (DM, SS and PC) were supported by the ESRC Research Centre
on Micro-Social Change (MiSoC) - ES/S012486/1. This research was supported in
part through computational resources provided by the Business and Local
Government Data Research Centre BLG DRC (ES/S007156/1) funded by the Economic
and Social Research Council (ESRC).

## References

  * Assaad et al. (2022) Charles K Assaad, Emilie Devijver, and Eric Gaussier.  Survey and evaluation of causal discovery methods for time series.  _Journal of Artificial Intelligence Research_ , 73:767–819, 2022.
  * Barabási & Albert (1999) Albert-László Barabási and Réka Albert.  Emergence of scaling in random networks.  _Science_ , 286(5439):509–512, 1999.
  * Bergstra et al. (2011) James Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl.  Algorithms for Hyper-Parameter Optimization.  In _Advances in Neural Information Processing Systems_ , volume 24. Curran Associates, Inc., 2011.
  * Biza et al. (2022) Konstantina Biza, Ioannis Tsamardinos, and Sofia Triantafillou.  Out-of-Sample Tuning for Causal Discovery.  _IEEE Transactions on Neural Networks and Learning Systems_ , pp. 1–11, 2022.  ISSN 2162-2388.  doi: 10.1109/TNNLS.2022.3185842.
  * Bühlmann et al. (2014) Peter Bühlmann, Jonas Peters, and Jan Ernest.  CAM: Causal additive models, high-dimensional order search and penalized regression.  _The Annals of Statistics_ , 42(6):2526–2556, 2014.  ISSN 0090-5364.
  * Chickering (2002) David Maxwell Chickering.  Optimal Structure Identification With Greedy Search.  _Journal of Machine Learning Research_ , 3(Nov):507–554, 2002.  ISSN ISSN 1533-7928.
  * Chobtham & Constantinou (2023) Kiattikun Chobtham and Anthony C Constantinou.  Tuning structure learning algorithms with out-of-sample and resampling strategies.  _arXiv preprint arXiv:2306.13932_ , June 2023.
  * Curth & van der Schaar (2023) Alicia Curth and Mihaela van der Schaar.  In Search of Insights, Not Magic Bullets: Towards Demystification of the Model Selection Dilemma in Heterogeneous Treatment Effect Estimation.  _arXiv preprint arXiv:2302.02923_ , June 2023.
  * Eberhardt (2017) Frederick Eberhardt.  Introduction to the foundations of causal discovery.  _International Journal of Data Science and Analytics_ , 3(2):81–91, March 2017.  ISSN 2364-4168.  doi: 10.1007/s41060-016-0038-6.
  * Erdös & Rényi (1959) P. Erdös and A. Rényi.  On random graphs.  _Publ. math. debrecen_ , 6(290-297):18, 1959.
  * Glymour et al. (2019) Clark Glymour, Kun Zhang, and Peter Spirtes.  Review of causal discovery methods based on graphical models.  _Frontiers in Genetics_ , 10, 2019.  ISSN 1664-8021.  doi: 10.3389/fgene.2019.00524.
  * Grünbaum et al. (2023) Daniel Grünbaum, Maike L. Stern, and Elmar W. Lang.  Quantitative probing: Validating causal models with quantitative domain knowledge.  _Journal of Causal Inference_ , 11(1), January 2023.  ISSN 2193-3685.  doi: 10.1515/jci-2022-0060.
  * Guyon et al. (2019) Isabelle Guyon, Alexander Statnikov, and Berna Bakir Batu.  _Cause effect pairs in machine learning_.  Springer, 2019.
  * Hoyer et al. (2008) Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf.  Nonlinear causal discovery with additive noise models.  In _Advances in Neural Information Processing Systems_ , volume 21. Curran Associates, Inc., 2008.
  * Hyvarinen (1999) A. Hyvarinen.  Fast and robust fixed-point algorithms for independent component analysis.  _IEEE Transactions on Neural Networks_ , 10(3):626–634, 1999.  doi: 10.1109/72.761722.
  * Kaiser & Sipos (2021) Marcus Kaiser and Maksim Sipos.  Unsuitability of NOTEARS for Causal Graph Discovery.  _arXiv:2104.05441 [cs, math, stat]_ , April 2021.
  * Kalainathan et al. (2020) Diviyan Kalainathan, Olivier Goudet, and Ritik Dutta.  Causal Discovery Toolbox: Uncovering causal relationships in Python.  _Journal of Machine Learning Research_ , 21(37):1–5, 2020.  ISSN 1533-7928.
  * Kalisch et al. (2012) Markus Kalisch, Martin Mächler, Diego Colombo, Marloes H. Maathuis, and Peter Bühlmann.  Causal inference using graphical models with the R package pcalg.  _Journal of Statistical Software_ , 47(11):1–26, 2012.  doi: 10.18637/jss.v047.i11.
  * Lachapelle et al. (2019) Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien.  Gradient-Based Neural DAG Learning.  In _International Conference on Learning Representations_ , September 2019.
  * Liu et al. (2010) Han Liu, Kathryn Roeder, and Larry Wasserman.  Stability Approach to Regularization Selection (StARS) for High Dimensional Graphical Models.  In _Advances in Neural Information Processing Systems_ , volume 23. Curran Associates, Inc., 2010.
  * Machlanski et al. (2023) Damian Machlanski, Spyridon Samothrakis, and Paul Clarke.  Hyperparameter Tuning and Model Evaluation in Causal Effect Estimation.  _arXiv preprint arXiv:2303.01412_ , March 2023.
  * Marx & Vreeken (2019) Alexander Marx and Jilles Vreeken.  Identifiability of Cause and Effect using Regularized Regression.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, KDD ’19, pp. 852–861, New York, NY, USA, July 2019. Association for Computing Machinery.  ISBN 978-1-4503-6201-6.  doi: 10.1145/3292500.3330854.
  * Paine et al. (2020) Tom Le Paine, Cosmin Paduraru, Andrea Michi, Caglar Gulcehre, Konrad Zolna, Alexander Novikov, Ziyu Wang, and Nando de Freitas.  Hyperparameter Selection for Offline Reinforcement Learning.  _arXiv preprint arXiv:2007.09055_ , July 2020.
  * Pearl (2000) Judea Pearl.  _Causality: Models, Reasoning, and Inference_.  Cambridge University Press, March 2000.  ISBN 978-0-521-77362-1.
  * Peters & Bühlmann (2014) J. Peters and P. Bühlmann.  Identifiability of Gaussian structural equation models with equal error variances.  _Biometrika_ , 101(1):219–228, March 2014.  ISSN 0006-3444.  doi: 10.1093/biomet/ast043.
  * Raghu et al. (2018) Vineet K. Raghu, Allen Poon, and Panayiotis V. Benos.  Evaluation of Causal Structure Learning Methods on Mixed Data Types.  In _Proceedings of 2018 ACM SIGKDD Workshop on Causal Discovery_ , pp. 48–65. PMLR, August 2018.
  * Ramsey et al. (2017) Joseph Ramsey, Madelyn Glymour, Ruben Sanchez-Romero, and Clark Glymour.  A million variables and more: The Fast Greedy Equivalence Search algorithm for learning high-dimensional graphical causal models, with an application to functional magnetic resonance images.  _International Journal of Data Science and Analytics_ , 3(2):121–129, March 2017.  ISSN 2364-4168.  doi: 10.1007/s41060-016-0032-z.
  * Ramsey et al. (2018) Joseph D Ramsey, Kun Zhang, Madelyn Glymour, Ruben Sanchez Romero, Biwei Huang, Imme Ebert-Uphoff, Savini Samarasinghe, Elizabeth A Barnes, and Clark Glymour.  Tetrad—a toolbox for causal discovery.  In _8th international workshop on climate informatics_ , pp. 29, 2018.
  * Rios et al. (2021) Felix L Rios, Giusi Moffa, and Jack Kuipers.  Benchpress: A Scalable and Versatile Workflow for Benchmarking Structure Learning Algorithms.  _arXiv preprint arXiv:2107.03863_ , May 2021.
  * Sachs et al. (2005) Karen Sachs, Omar Perez, Dana Pe’er, Douglas A. Lauffenburger, and Garry P. Nolan.  Causal protein-signaling networks derived from multiparameter single-cell data.  _Science_ , 308(5721):523–529, 2005.  doi: 10.1126/science.1105809.
  * Shimizu et al. (2006) Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, and Antti Kerminen.  A Linear Non-Gaussian Acyclic Model for Causal Discovery.  _Journal of Machine Learning Research_ , 7(Oct):2003–2030, 2006.  ISSN ISSN 1533-7928.
  * Spirtes & Glymour (1991) Peter Spirtes and Clark Glymour.  An Algorithm for Fast Recovery of Sparse Causal Graphs.  _Social Science Computer Review_ , 9(1):62–72, April 1991.  ISSN 0894-4393.  doi: 10.1177/089443939100900106.
  * Spirtes et al. (1993) Peter Spirtes, Clark Glymour, and Richard Scheines.  Discovery algorithms for causally sufficient structures.  In _Causation, Prediction, and Search_ , pp. 103–162. Springer New York, New York, NY, 1993.  ISBN 978-1-4612-2748-9.  doi: 10.1007/978-1-4612-2748-9_5.
  * Strobl (2021) Eric V. Strobl.  Automated hyperparameter selection for the PC algorithm.  _Pattern Recognition Letters_ , 151:288–293, November 2021.  ISSN 0167-8655.  doi: 10.1016/j.patrec.2021.09.009.
  * Sun et al. (2013) Wei Sun, Junhui Wang, and Yixin Fang.  Consistent Selection of Tuning Parameters via Variable Selection Stability.  _Journal of Machine Learning Research_ , 14(107):3419–3440, 2013.  ISSN 1533-7928.
  * Tönshoff et al. (2023) Jan Tönshoff, Martin Ritzert, Eran Rosenbluth, and Martin Grohe.  Where Did the Gap Go? Reassessing the Long-Range Graph Benchmark.  _arXiv preprint arXiv:2309.00367_ , September 2023.
  * Tu et al. (2019) Ruibo Tu, Kun Zhang, Bo Bertilson, Hedvig Kjellstrom, and Cheng Zhang.  Neuropathic Pain Diagnosis Simulator for Causal Discovery Algorithm Evaluation.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Van den Bulcke et al. (2006) Tim Van den Bulcke, Koenraad Van Leemput, Bart Naudts, Piet van Remortel, Hongwu Ma, Alain Verschoren, Bart De Moor, and Kathleen Marchal.  SynTReN: A generator of synthetic gene expression data for design and analysis of structure learning algorithms.  _BMC bioinformatics_ , 7:43, January 2006.  ISSN 1471-2105.  doi: 10.1186/1471-2105-7-43.
  * Zhang et al. (2021a) Baohe Zhang, Raghu Rajan, Luis Pineda, Nathan Lambert, André Biedenkapp, Kurtland Chua, Frank Hutter, and Roberto Calandra.  On the Importance of Hyperparameter Optimization for Model-based Reinforcement Learning.  In _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , pp. 4015–4023. PMLR, March 2021a.
  * Zhang et al. (2021b) Keli Zhang, Shengyu Zhu, Marcus Kalander, Ignavier Ng, Junjian Ye, Zhitang Chen, and Lujia Pan.  gcastle: A python toolbox for causal discovery.  _arXiv preprint arXiv:2111.15155_ , 2021b.
  * Zheng et al. (2018) X. Zheng, Bryon Aragam, P. Ravikumar, and E. Xing.  DAGs with NO TEARS: Continuous Optimization for Structure Learning.  In _NeurIPS_ , 2018.
  * Zheng et al. (2020) Xun Zheng, Chen Dan, Bryon Aragam, P. Ravikumar, and E. Xing.  Learning Sparse Nonparametric DAGs.  In _AISTATS_ , 2020.

##  Appendix A Experimental Details

###  A.1 Hyperparameters

Table 1: Hyperparameter search spaces defined per algorithm. algorithm | hyperparameters and values
---|---
ANM | $alpha\in\\{0.001^{*},0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.5\\}$
CAM | $cutoff\in\\{0.001^{*},0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.5\\}$
|  $score=$ nonlinear
|  $selmethod=$ gamboost
|  $prunmethod=$ gam
FCI | $alpha\in\\{0.001^{*},0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.5\\}$
|  $test=$ fisher-z-test
FGES | $penaltyDiscount\in\\{0.001,0.01,0.1,0.25,0.5,0.75,1,1.5^{*}\\}$
|  $score=$ sem-bic
LiNGAM | $thresh\in\\{0.001,0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.5^{*}\\}$
| $max\\_iter\in\\{100^{*},1000\\}$
NOTEARS | $lambda1\in\\{0.001,0.002,0.005,0.01,0.02,0.05,0.1,0.2^{*},0.3,0.5\\}$
| $max\\_iter\in\\{100^{*},1000\\}$
| $w\\_threshold\in\\{0.001,0.002,0.005,0.01,0.02,0.05,0.1,0.2^{*},0.3,0.5\\}$
|  $loss\\_type=$ l2
| $h\\_tol=1e-8$
| $rho\\_max=1e+16$
NOTEARS MLP | $lambda1\in\\{0.001,0.01^{*},0.1\\}$
| $lambda2\in\\{0.001,0.01,0.1^{*}\\}$
| $w\\_threshold\in\\{0.1,0.3,0.5^{*}\\}$
| $hidden\\_layers=1$
| $hidden\\_units\in\\{8,16^{*},32\\}$
| $max\\_iter=100$
| $h\\_tol=1e-8$
| $rho\\_max=1e+16$
PC | $alpha\in\\{0.001,0.002^{*},0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.5\\}$
|  $indepTest=$ gaussCItest
${}^{*}=$ found to perform best on average across all simulations (sim_mean).

###  A.2 Summary of Algorithms

Table 2: Summary of incorporated algorithms and their sources. Recommended default hyperparameters have been derived from respective papers as much as possible. If necessary, they have been further supplemented with defaults suggested within respective packages. algorithm | default hyperparameters | package | paper
---|---|---|---
ANM | $alpha=0.05$ | gCastle | (Hoyer et al., 2008)
CAM | $cutoff=0.001$ | cdt | (Bühlmann et al., 2014)
|  $score=$ nonlinear |  |
|  $selmethod=$ gamboost |  |
|  $prunmethod=$ gam |  |
FCI | $alpha=0.01$ | tetrad | (Spirtes et al., 1993)
FGES | $penaltyDiscount=2.0$ | tetrad | (Ramsey et al., 2017)
LiNGAM | $thresh=0.3$ | gCastle | (Shimizu et al., 2006)
| $max\\_iter=1000$ |  |
NOTEARS | $lambda1=0.1$ | gCastle | (Zheng et al., 2018)
| $max\\_iter=100$ |  |
| $w\\_threshold=0.3$ |  |
|  $loss\\_type=$ l2 |  |
| $h\\_tol=1e-8$ |  |
| $rho\\_max=1e+16$ |  |
NOTEARS MLP | $lambda1=0.01$ | gCastle | (Zheng et al., 2020)
| $lambda2=0.01$ |  |
| $w\\_threshold=0.3$ |  |
| $hidden\\_layers=1$ |  |
| $hidden\\_units=10$ |  |
| $max\\_iter=100$ |  |
| $h\\_tol=1e-8$ |  |
| $rho\\_max=1e+16$ |  |
PC | $alpha=0.01$ | pcalg | (Spirtes & Glymour, 1991)

###  A.3 Summary of Packages

Table 3: Summary of incorporated algorithm packages. package | paper | link
---|---|---
gCastle | (Zhang et al., 2021b) | https://github.com/huawei-noah/trustworthyAI/tree/master/gcastle
cdt | (Kalainathan et al., 2020) | https://fentechsolutions.github.io/CausalDiscoveryToolbox
tetrad | (Ramsey et al., 2018) | https://cmu-phil.github.io/tetrad/manual/
pcalg | (Kalisch et al., 2012) | https://cran.r-project.org/package=pcalg

##  Appendix B Supplemental Results

The following results complement the ones presented in the main content of the
paper. Although they do not change the overall conclusions of the paper, they
offer additional analysis that may facilitate a deeper understanding of the
problem and the final outcomes.

###  B.1 Best Performances

Figures 13 and 14 focus on performances achieved specifically with the best
(oracle) hyperparameters. These correspond to the best performances presented
in Figure 3 (blue bars). Some important observations: a) algorithms differ in
performance even with access to a hyperparameter oracle, b) number of graph
nodes and edge density can significantly impact performance, and c) no
algorithm performs best under all conditions (linear vs. nonlinear).

(a) linear SEM with Gumbel noise

(b) nonlinear (GP) SEM with Gaussian noise

Figure 13: Best performances against ER graphs. Error bars are standard
errors.

(a) linear SEM with Gumbel noise

(b) nonlinear (GP) SEM with Gaussian noise

Figure 14: Best performances against SF graphs. Error bars are standard
errors.

###  B.2 Large Graphs with Large Sample Size

Previous results showed that $p=50$ graphs are much more challenging than
smaller ones. Figure 15a demonstrates that with enough data samples, even such
larger graphs are possible to solve accurately. As presented in the subfigure
(b), increased sample size can also help with robustness. Note also how the
degree of the benefits vary between algorithms.

(a) best performances

(b) worst performances

Figure 15: Performances for ER1 $p=50$ graphs.

###  B.3 Performance vs. Hyperparameter Quality

Figure 16 complements Figure 3 from the main content. Instead of showing raw
SHD values achieved with different hyperparameter values, it shows differences
in SHD between best and fixed hyperparameters, and between fixed and worst
ones. This point of view makes the following observations more prominent: a)
there is very little performance difference between the best and fixed
hyperparameters, b) the worst hyperparameters substantially degrade achieved
performance, and c) robustness to misspecified hyperparameters (orange bars)
can vary significantly between algorithms.

Figure 16: Change in performance as a cost of changing the quality of selected
hyperparameters. Colours: blue – difference in $SHD$ achieved with best
(oracle) and fixed hyperparameters; orange – difference in $SHD$ achieved with
fixed and worst hyperparameters.

###  B.4 Winning Algorithms vs. Simulation Properties

Figure 17 complements Figure 10 from the main content but removes
hyperparameters from the picture in order to analyse how DGP properties alone
affect winning odds of the algorithms. The results presented here involve the
use of the best hyperparameter values. From this perspective, it is clear that
no algorithm is the best under all conditions, and that SEM types involved and
edge density have the strongest impact on the winning odds.

Figure 17: Percentage of wins per algorithm.

###  B.5 Winning Algorithms vs. Hyperparameter Quality

Figure 18 offers a point of view alternative to Figure 10. More specifically,
Figure 10 shows how winning odds change as specific DGP properties change as
well, while aggregating the numbers across all other DGP properties. Whereas
Figure 18 aggregates the numbers over the sample size and simulated graph type
only. For example, Figure 10d shows how winning odds change as edge density
changes, and the numbers are averaged across all other DGP properties (SEM
types, sample size, etc.). On the other hand, the upper-left subfigure in
Figure 18 presents the numbers that correspond to $data\\_sem=gumbel$,
$graph\\_d=1$, $graph\\_p=10$ and the averaging is done across only data_n and
graph_type specifically as there is no significant difference in winning odds
across these two specific dimensions, as clearly demonstrated by Figure 10.

Figure 18: Percentages of winning algorithms under different data, graph and
hyperparameter conditions.

##  Appendix C A Guide To Algorithm Selection

###  C.1 General Recommendations

  * •

Current algorithms seem to work reasonably well for sparse graphs of up-to 20
nodes. Bigger graphs (50 nodes) are also possible to solve, but more data
might be required ($10,000$ samples) to achieve good accuracy. The accuracy of
recovered structures drops dramatically for dense graphs.
Recommendation: Stick to sparse graphs with up-to 20 nodes (moderate amount of
data) or up to 50 nodes (a lot of data). Avoid dense graphs.

  * •

No single algorithm is the best option for all problems. Some perform the best
under very specific conditions.
Recommendation: Choose an algorithm that is the most likely to accurately
solve the problem at hand based on assumptions derived from data.

  * •

The best choice of an algorithm may depend not only on graph and data
properties, but also on the availability of quality hyperparameters. This is
because algorithms vary in robustness to misspecified hyperparameters.
Recommendation: When selecting the best algorithm for the problem at hand,
take into account the type of hyperparameters that are available and
algorithm’s robustness to misspecified hyperparameters.

###  C.2 Hyperparameter Selection Strategies

Optimal hyperparameters are almost never available in structure recovery
problems due to inaccessible ground truth. Some methods provide scores that
can be used to decide whether a set of hyperparameters is better than others
for the same algorithm. However, this strategy cannot be used to compare
different algorithms to each other, as they are likely to use different score
metrics (while some use none at all). In addition, in order to use those
scores, an algorithm’s internal code must be modified in most cases, creating
a substantial barrier to practitioners.

Thus, default hyperparameters might be a reasonable selection strategy as they
often work almost as well as the optimal ones. This is especially the case if
the recommended defaults have been derived from data problems similar in
nature to the problem at hand.

If, however, the problem to solve is believed to be fairly unique, blindly
using default hyperparameters without considering other factors might not be
safe. In this case, the best course of action (assuming hyperparameter tuning
is not an option) might be to still use default hyperparameters but choose the
algorithm that is the most robust to hyperparameter misspecification under
specific graph and data conditions that are believed to be applicable to the
problem at hand.

###  C.3 How to Select Algorithms

Figure 19 summarises the best algorithm choices based on Figure 18, and which
takes into consideration data and graph properties as well as the types of
hyperparameters available.

Figure 19: Recommended algorithm choices based on the number of graph nodes,
SEM types and edge density in the graph. The final choice depends also on the
type of hyperparameters available – see ‘quality’ colours. In case there is no
clear winner, multiple choices are provided in the order of higher winning
percentage.