arXiv:1507.02018v3 [] 6 Oct 2017

Inferring large graphs using ℓ1-penalized likelihood
Magali Champion
Laboratoire MAP5, Université Paris Descartes,
Sorbonne Paris Cité, France
Victor Picheny
INRA, UR875 Applied Mathematics and Computer Science Unit,
Castanet-Tolosan, France
and
Matthieu Vignes
Institute of Fundamental Sciences, Massey University,
Palmerston North, New Zealand.

Abstract

1

We address the issue of recovering the structure of
large sparse directed acyclic graphs from noisy observations of the system. We propose a novel procedure based on a specific formulation of the ℓ1 norm regularized maximum likelihood, which decomposes the graph estimation into two optimization sub-problems: topological structure and node
order learning. We provide convergence inequalities
for the graph estimator, as well as an algorithm to
solve the induced optimization problem, in the form
of a convex program embedded in a genetic algorithm. We apply our method to various data sets
(including data from the DREAM4 challenge) and
show that it compares favorably to state-of-the-art
methods. This algorithm is available on CRAN as
the R package GADAG.
Keywords: Directed Acyclic Graphs Lasso Convex
program Optimization

Revealing the true structure of a complex system is
paramount in many fields to identify system regulators, predict its behavior or decide where interventions are needed to disentangle direct relationships
(Newman, 2003; Souma et al., 2006; Verma et al.,
2014). This problem can often be seen as a graph inference problem. Given observational data, we aim
at predicting the presence (or absence) of edges between elements of the system, which form the vertices of a graph. Edges specify the relational structure of the system depicted as a graph or network.
As a motivating problem, the reconstruction of Gene
Regulatory Networks (GRN), which model activation and inhibition relationships between genes, is
one of the main challenges in modern computational
biology (Barabási & Oltvai, 2004).
A popular approach consists in assuming that the
data are generated by a Directed Acyclic Graph
(DAG) (Pearl, 2009). DAGs are made of a collection of vertices, which stand for variables, and
directed edges to model the dependency structure
among the variables, avoiding self-loops and cycles.
However, inferring a DAG is a rather challenging
1

Introduction

cycles requires temporal data, e.g. in the context
of dynamic Bayesian networks (Perrin et al., 2003;
Dondelinger et al., 2013), data which are difficult
and very expensive to collect despite efforts in this
direction (Sachs et al., 2009).

problem. Firstly, the number of nodes p of the graph
may be so large that exploring relevant DAG topologies is simply infeasible, since the number of possible
DAG structures is super-exponential in p (Robinson,
1973; Koivisto & Sood, 2004; Tsamardinos et al.,
2006; Grzegorczyk & Husmeier, 2008). Another dimension flaw occurs when p, even being reasonable,
is larger than the number of observations, and parameter estimation is jeopardized. High-dimensional
statistical techniques are then needed to overcome
this issue (Bühlmann & van de Geer, 2011; Giraud,
2015). Secondly, even if the ratio between p and
the sample size n is not impeding model estimation, the nature of the data can be an additional
obstacle (Ellis & Wong, 2008; Guyon et al., 2010;
Fu & Zhou, 2013). The available observational data
are in general not sufficient to identify the true underlying DAG, and can only determine an equivalence class of DAGs (Verma & Pearl, 1991). This
approach relies on the assumption that the joint distribution is Markov and faithful with respect to the
true graph (Spirtes et al., 2000).

In this work, we focus on Gaussian structural
equation models associated with maximum likelihood estimators (MLE). In the last years, the ℓ0 regularization of the MLE drew the attention of
a large number of works since it leads to infer
sparse graphs. In DAGs, a not necessarily topological ordering of the nodes can always be defined according to edge distribution (Kahn, 1962).
Identifying this ordering is known to be a challenging problem (Cook, 1985). Additional data,
like gene knock-out data or more general perturbations data (Maathuis et al., 2010; Shojaie et al.,
2014) can give information in that way. More generally, biological prior knowledge, retrieved from
specific data bases, can be used to assist the network reconstruction algorithm (Husmeier & Werhli,
2007), or a partial knowledge of the network can
inform the inference process efficiently, e.g. in
the semi-supervised framework of Mordelet & Vert
(2008). For a known order among the variables in
the graph, Shojaie & Michailidis (2010) present results for the estimation of high-dimensional graphs
based on independent linear regressions using an
adaptive Lasso scheme. When the order of the variables is unknown, van de Geer & Bühlmann (2013)
studied the convergence of the ℓ0 -penalized likelihood. However, the ℓ0 -regularized approaches
(Silander & Myllymäki, 2006; Hauser & Bühlmann,
2012) remain impractical for estimating graphs with
more than 20 vertices, either due to an exhaustive exploration of the set of DAGs or overfitting
(Chen & Chen, 2008). Quite recently, Aragam et al.
(2015) explored a penalized least-squares estimator
for a variety of concave penalization terms. They
obtain theoretical guarantees of sparsity bounds and
consistency results in model selection. Their theoretical results would greatly benefit an implementation
of the methods and an empirical study to demonstrate the effectiveness of the approach. The uni-

A large number of methods have been proposed
for estimating DAGs, including for instance scorebased methods (Bayesian score, Friedman & Koller
2003 or Bayesian Information Criterion, Schwarz
1978), complex space sampling (Zhou, 2011) or the
PC algorithm (Spirtes et al., 2000). The latter has
been proved to be uniformly consistant in the highdimensional case, but requires a test of conditional
independences that quickly becomes computationally intractable (Kalisch & Bühlmann, 2007). Recent works stress the limitations of the absence of
cycles in DAGs in the study of complex systems
(De Smet & Marchal, 2010). Wright (1921) already
described more general directed dependencies between variables when introducing genetic path analysis, but the data were then very limited. Structural Equation Modelling later introduced the notion of noise measurement (Hoyle, 1995) and Pearl
(2009) extended them beyond linearity. Moreover,
the directed cyclic graph (Spirtes, 1995) framework
received little attention as compared to its acyclic
counterpart. Finally, the actual discovery of causal
2

work inference problem is to consider that each variable X i (1 ≤ i ≤ p) of the DAG can be represented
as a linear function of all other variables X j (j 6= i)
through the Gaussian Structural Equation Model:

fying framework for pseudolikelihood-based graphical modelling of Kahre et al. (2015) extends classical
regularization methods. The authors obtain theoretical convergence results, and offer an algorithm with
an associated implementation. Their simulation results are quite promising, in particular in terms of
computational time.
Our objective consists in overcoming this drastic
dimensional limitation, and find inference strategies
for graphs with up to several hundred nodes. Such
strategies must ensure a high level of sparsity, be
supported by computationally affordable algorithms,
while preserving sound theoretical bases. Here, we
propose to use the ℓ1 -regularization, similarly to
Fu & Zhou (2013) and Shojaie & Michailidis (2010),
to penalize the MLE. From a computational point
of view, this regularization makes the criterion to
maximize partially convex while ensuring sparse estimates. Our contribution is two-fold: firstly, we provide convergence inequalities that guarantee good
theoretical performances of our proposed estimator
in the sparse high-dimensional setting. Secondly, we
provide an efficient algorithm to infer the true unknown DAG, in the form of a convex program embedded in a genetic algorithm.
The next section covers the model definition and
the associated penalized MLE problem. Section 3
details the convergence inequalities, and Section 4
our inference algorithm. Section 5 reports numerical
experiments both on toy problems and realistic data
sets.

∀j ∈ J1, pK, X j =

p
X
(G0 )ji X i + εj ,

(1)

i=1

with εj ∼ N (0, σj2 ) (σj2 known) a Gaussian residual
error term. The set of edges E, which is assumed
to be of size s (s ≤ p(p − 1)/2), corresponds to the
non-zero coefficients of G0 , i.e. (G0 )ji encodes the
relationship from variable X i to variable X j .
Assume that we observe an n-sample consisting
of n i.i.d. realizations (X 1 , ..., X p ) of Equation (1),
distributed according to a N (0, Σ) law where Σ is
non-singular. We denote by X := (X 1 , ..., X p ) the
n × p data matrix. The relations between the variables can then be represented in its matrix form:
X = XG0 + ε,

(2)

where G0 = ((G0 )ji )1≤i,j≤p is the p × p matrix compatible with the graph G0 and ε := (ε1 , ..., εp ) is the
n × p matrix of noise vectors.
The negative log-likelihood of the model is then
(Rau et al., 2013):
ℓ(G) =

p
X
np
log σj
log(2π) + n
2
j=1

p
n X
X

1
j 2
X
(I
−
G)
, (3)
+
k
σ2
k=1 j=1 j

where I denotes the p × p identity matrix and Xk
the vector of length p corresponding to the k-th observation of X 1 , ..., X p .
To recover the structure of the DAG G0 and make
2.1 DAG’s modelling and estimation the estimated graph sparse enough, we focus on a penalized maximum likelihood procedure (Bickel & Li,
This work considers the framework of an un- 2006):
known DAG G0 = (V, E), consisting of vertices
(4)
Ĝ = argmin {ℓ(G) + λ pen(G)},
V = {1, ..., p} and a set of edges E ⊆ V × V .
G∈GDAG
The p nodes are associated to random variables
X 1 , ..., X p . A natural approach, developped by where ℓ(.) is the negative log-likelihood of Equation
Meinshausen & Bühlmann (2006) to solve the net- (3), pen(.) is a penalization function, λ is a trade-off

2

The ℓ1-penalized likelihood
for estimating DAGs

3

where, for any
M := (Mij )1≤i,j≤p , we denote
Pmatrix
by kMkF = i,j (Mij )2 the Frobenius norm and by
P
kMk1 = i,j |Mij | the ℓ1 -norm.

parameter between penalization and fit to the data,
and GDAG is the set of p×p matrices compatible with
a DAG over p nodes.
In the setting of Gaussian structural equation
models with equal noise variance, Peters et al.
(2011, 2014) showed that the true DAG was identifiable for respectively discrete and continuous data.
In a nutshell, it implies that the true DAG could be
inferred, not just the Markov equivalence class of the
underlying DAG - a partially directed graph exactly
encoding the conditional dependency structure. Using an ℓ0 -norm regularization in Equation (4) is an
attractive option to infer sparse graphs. From a computational point of view, the main difficulty when
solving the optimization problem in Equation (4)
lies in exploring the set of DAGs GDAG . (Chickering,
1996) showed it to be an NP-hard problem: an ℓ0 regularization does not set a favorable framework
for this task. To avoid the whole exploration of
GDAG , a dynamic programming method has been
proposed in Silander & Myllymäki (2006), using a
particular decomposition of the ℓ0 -penalized maximum likelihood. The greedy equivalent search algorithm of Chickering (2002) is a hill climbing alternative method. (Hauser & Bühlmann, 2012) rather
restricted the search space to the smaller space of
equivalence classes, and they provide an efficient
algorithm without enumerating all the equivalent
DAGs. They showed that they are asymptotically
optimal under a faithfulness assumption (i.e. independences in the distribution are those read from
G0 ). However, none of the approaches above can be
used on high-dimensional data to estimate graphs
with a large number of nodes. In this context, we
focus on the ℓ1 -norm convex regularization instead
of ℓ0 for its sparse, high-dimensional and computational properties.
The ℓ1 -regularization clearly improves the computation of (4). It allows us to write a convex formulation of the problem (see Section 2.2). Given
Equation (3) and omitting constant terms, the ℓ1 penalized likelihood estimator we consider is:


1
2
Ĝ = argmin
kX(I − G)kF + λ kGk1 , (5)
n
G∈GDAG

2.2

A new formulation for the estimator

We propose here a new formulation of the minimization problem of Equation (5). It naturally decouples
the initial problem into two steps of the minimisation procedure: node ordering and graph topology
search. A key property is that any DAG leads to a
topological ordering of its vertices, denoted ≤, where
a directed path from node X i to node X j is equivalent to X j ≤ X i (Kahn, 1962; Cormen et al., 2001)
(see Example 1 below for more explanations). This
ordering is not unique in general. Proposition 2.1
from Bühlmann (2013) then gives an equivalent condition for a matrix to be compatible with a DAG.
Proposition 2.1 (Bühlmann, 2013) A matrix G
is compatible with a DAG G if and only if there exists
a permutation matrix P and a strictly lower triangular matrix T such that:
G = PTPT.
Graphically, the permutation matrix sets an ordering of the nodes of the graph and is associated to a
complete graph. The strictly lower triangular matrix
T sets the graph structure, i.e. the non-zero entries
of G, as illustrated in Example 1.
Example 1 Consider the DAG G given in Figure
1 (left). The corresponding matrix G can then be
written as the strictly lower-triangular matrix T by
permutation of its rows and columns using P :

0
2

G=
0
0
0
4

0
0
0
0
0

0
1
0
0
0

7
6
0
0
0


5
4

T
0
 = PTP ,
3
0


0
0

T =
3
5
4

0
0
0
0
1

0
0
0
7
6


0
0 0
0
0 0


0
0 0
and
P
=


0
0 0
1
2 0


0
0
1
0
0

0
0
0
1
0

1
0
0
0
0



0
1

0
.
0
0

Looking at the non-zero values of P column by column, P defines a node hierarchy X 5 ≤ X 3 ≤ X 4 ≤
X 1 ≤ X 2 compatible with the topological orderings
of G. Graphically, P is associated to the complete
graph represented in Figure 1 (bottom). The dashed
edges then correspond to the lower zero entries of
T . Note that since X 3 is not connected with X 1 ,
X 4 and X 5 , four topological ordering are possible
(X 5 ≤ X 4 ≤ X 1 ≤ X 3 , X 5 ≤ X 4 ≤ X 3 ≤ X 1 ,
X 5 ≤ X 3 ≤ X 4 ≤ X 1 and X 3 ≤ X 5 ≤ X 4 ≤ X 1 ).
X2

1
X3

2

6
4
7
3

X5

1

X2

7
X

4

3

X1
5

4

argmin
(P,T )∈C



1
X(I − P T P T )
n

2
F



+ λ kT k1 , (6)

where the optimization space of constraints C
is defined as C = Pp (R) × Tp (R), with Pp (R)
the set of permutation matrices and Tp (R) the
set of strictly lower-triangular matrices.
Note
that a similar formulation has already been proposed by van de Geer & Bühlmann (2013) to ensure
good theoretical properties for the ℓ0 -penalized loglikelihood estimation. However, it has never been
exploited from a computational point of view to recover the graph structure optimizing problem (5).
In the following two sections, we propose a theoretical analysis of the proposed estimator (Section 3)
and a computationally effficient algorithm to solve
Problem (6) (Section 4).

Convergence inequalities for
the DAG estimation

The main result of this section deals with convergence rates: in Theorem 1, we provide upper
bound for error associated with the ℓ1 -penalized
maximum likelihood estimator considered in Equation (6), both in prediction (Equation 7) and estimation (Equation 8). Following the works of
van de Geer & Bühlmann (2013) on the ℓ0 -penalized
maximum likelihood estimator and of Bickel et al.
(2009) on the Lasso and the Dantzig Selector, we obtain two convergence results under some mild sparsity assumptions, when the number of variables is
large but upper bounded by a function ϕ(n) of the
sample size n.

2

6

(P̂ , T̂ ) =

3

5

X4

X3

X1

Using Proposition 2.1, the estimator in (5) leads
to the following equivalent optimization problem:

X5

Figure 1: An example of DAG G (top) and the action
of P and T on G: P is associated to a complete graph 3.1
that orders the nodes of the graph (bottom) and T
sets the weights on the edges. The dashed edges
correspond to null weight edges (a zero entry in T ). For a
graph
5

Estimating the true order of
variables
known ordering among the variables of the
(Shojaie & Michailidis, 2010), an unrealistic

a wrong permutation.
In Figure 2, we represent the permuted graph
G0 (P ) (bottom) associated to the graph G0 (top). The
latter is obtained from G0 after permutation of its
nodes using P P0T , where P0 (corresponding to the
matrix P in Example 1) defines a right order of variables.

assumption in many applications, the DAG inference problem can be cast in a convex optimization
problem. To provide convergence inequalities of the
proposed estimator in the most general case of an
unknown order we consider here, we first focus on
the problem of estimating the true variable order.
Let us denote by Π0 the set of permutation matrices
compatible with the true DAG G0 :

Π0 = P ∈ Pp (R), P T G0 P ∈ Tp (R) .

1
X3

Π0 contains one or more permutation matrice(s) (see
Example 1). We will have to make a decision as to
whether the estimated order of variables P̂ given by
Equation (6) is in Π0 or not.
To answer this question, we investigate the effect
of learning an erroneous order of variables P ∈
/ Π0 .
We introduce the following additional notations: for
any permutation matrix P ∈ Pp (R), we denote by
G0 (P ) the matrix defined as:

X2

0
1
0
0
0

3

X5

X1
2

6

7

4
5

X4
3

X3

Figure 2: The graph G0 (top) and the permuted
graph G0 (P ) (bottom) associated to the permutation P .

3.2

Assumptions on the model

For a square matrix M ∈ Mp×p (R) and a subset S of
J1, pK2 , we denote by MS ∈ Mp×p (R) the matrix that
has the same elements as M on S and zero on the
complementary set S C of S. We now introduce the
assumptions we used to obtain statistical properties
of our estimator.

Example 2 Let
0
0
0
0
1

5

1

with T0 = P0T G0 P0 a lower triangular decomposition of G0 . From a graphical point of view, while
P ∈
/ Π0 , the graph G0 (P ) associated to G0 (P ) is obtained from G0 by permuting some of its nodes (see
Example 2), otherwise, if P ∈ Π0 , G0 (P ) = G0 . We
also denote by ε(P ) := X − XG0 (P ) the associated
residual term. We denote by Ω(P ) the covariance
matrix of ε(P ) and ωj (P ) := Var(εj (P )) the associated noise variances of each node.
With these notations and checking that the assumptions presented in Section 3.2 hold, we ensure
that, with large probability, we choose a right order
of variables and the estimated graph converges to
the true graph when n and p grow to infinity (see
Section 3.3).

0
0
1
0
0

X1

4

X5

G0 (P ) = P T0 P T ,


0
0

P =
0
1
0

2
6

7
X4

X2


1
0

0
/ Π0
∈
0
0

Hypotheses Nothing
6

results are still valid even if σ 2 6= 1, by small modifications in the constant terms as long as they are
∀j ∈ J1, pK, Var(εj ) = σ 2 .
all equal.
Assumption H5 deserves a special attention since
H2 There exists σ02 , independent of p and n, such it bounds the high dimensional setting. The considthat
ered problem is obviously non-trivial and requires a
max Var(X j ) ≤ σ02 .
sufficient amount of information. A more detailled
1≤j≤p
discussion about assumptions H3 and H5 is proH3 There exists λ∗ > 0 such that the minimal eigenposed in Section 3.4.
value of the covariance matrix Σ of X satisfies
Assumption H6 is a natural extension of the Re∗
stricted Eigenvalue condition of Bickel et al. (2009)
λmin ≥ λ > 0.
to our multi-task setting. More precisely, denoting
H4 There exists gmax < ∞ such that the maximal


0
X❆
weight of the DAG G0 is bounded
❆❆
 n × p,
❆❆
X̃ = 
j
max (G0 )i ≤ gmax .
0
X

H1 There exists σ 2 > 0 such that

1≤i,j≤p

p2

H5 The number of nodes p satisfies

H6 is equivalent to assuming that the Gram maT
is non-degenerate on a restricted cone
trix X̃ X̃
n
(Lounici et al., 2009; Bühlmann & van de Geer,
2011). Notice that this condition is very classical in
the literature. It yields good practical performance
even for small sample sizes, and some recent works
discuss an accurate population eigenvalue estimation even in a large dimension setting (Mestre, 2008;
El Karoui, 2008; Liu et al., 2014; Ledoit & Wolf,
2015).
The last assumption H7 is an identifiability condition needed to ensure that the estimated permutation P̂ is in Π0 . This assumption was introduced
by van de Geer & Bühlmann (2013) as the “omegamin” condition. In a sense, it separates the set of
compatible permutations from its complement in a
finite sample scenario.

p log p = O(n).
H6 There exists κ(t) > 0 with 1 ≤ t ≤ p2 such that:


kXMkF
min √
≥ κ(t),
n kMS kF
where the minimum is taken over the set of p × p
matrices satisfying kMS C k1 ≤ 3 kMS k1 , with S ⊂
J1, pK2 and |S| ≤ t.
n
× √1s such that, for
H7 There exists 0 < η ≤ C p log
p
all permutations P ∈
/ Π0 ,
p

2 1
1X
|ωj (P )|2 − 1 > .
p j=1
η
Assumption H1 states that the noise variances
are the same among all variables. This assumption is clearly hard to test in practice but makes
the problem identifiable and ensures that we can
recover the true DAG. Otherwise, minimizing (5)
only leads to the identification of one element of the
Markov equivalence class of the true DAG (partially
directed graph). To simplify the theoretical results
and proofs, until the end of this work, we assume
that the noise variances σ 2 are equal to one. Our

3.3

Main result

The result we establish in this section is doubleedged: (a) with large probability, we ensure that the
estimated P̂ belongs to Π0 , and (b) we provide convergence inequalities both in prediction and estimation for the graph estimated from the minimisation
problem (6). This result clearly states the desirable
7

Note however that universal conditions cannot be
theoretical properties of the derived estimator, assuming reasonable conditions on the complex system overcome and the ultra-high dimension settings (e.g.
Wainwright (2009); Verzelen (2012)) is an insurembedding the data.
mountable limit.
Theorem 1 Assume thatP
H1−7 are satisfied, with s
⊂ J1, p2 K in H6 such that i,j 1(G0 )j 6=0 ≤ s (G0 is s- Minimal eigenvalue condition Assumption H
i
3
q
log
p
sparse). Let λ = 2C s1/2 n . Then, with probabil- ensures that the minimal eigenvalue of the covariance matrix Σ of X is not too small. In the highity greater than 1 − 5/p, any solution Ĝ = P̂ T̂ P̂ T of
dimensional scenario, this could be hard to verthe minimization problem (6) satisfies that P̂ ∈ Π0 .
ify, λmin decreasing while n, p growing to infinity
Moreover, with at least the same probability, the fol(Hogben, 2007). A natural bound for λmin is:
lowing inequalities hold:
1
2
16C 2 3/2 log p
1
√ ,
λmin ≥
(9)
2
X Ĝ − XG0 ≤ 2 s
.
(7)
p max (1, gmax ) (1 + s)
n
κ (s)
n
F
with gmax and s as in H4 and Theorem 1.
r
log
p
16C
Assumption H3 can thus
√ be relaxed by allowing
.
(8)
s5/2
Ĝ − G0 ≤ 2
λmin to decrease with 1/p s. The price to pay for
κ (s)
n
1
this relaxation is a data dimensionality reduction
The proof of this result is deferred in Section C of p3 log(p) = O(n), which automatically implies:
the Supplementary Materials.
r
r
Theorem 1 states that with probability at least
log(p)
2p log(p)
3λmin
−2
− 3σ0
> 0,
1 − 5/p, we choose a compatible order of variables
4
n
n
over the set of permutations. Inequalities (7) and (8)
give non-asymptotic upper bounds on the loss under with Equation (9) (for more details, see Section A of
conditions depending on s, p and n (see Section 3.4). the Supplementary Material Proof details).
They also ensure that the estimated T̂ is close to the
true T0 with large probability.

4

3.4

Discussion
on
the
dimensional scenario

high- 4.1

Inference algorithm

Global algorithm overview

In this section, we present GADAG (Genetic Algorithm for learning Directed Acyclic Graphs), a
computational procedure devoted to solve Equation (6) and available as a R package on CRAN
at https://cran.r-project.org/package=GADAG.
Although decomposing the original problem made it
more natural to handle, this problem is still a very
challenging task from an optimization point of view,
due to the different nature of the variables P and T ,
the non-convexity of the cost function and the high
dimension of the search space.
An intuitive approach consists in using an alternating scheme: one of the variables P or T is fixed
and the other one is sought so as to optimize the

Sparsity of the graph Assumption H7 and Theorem 1 naturally require a trade-off between signal
sparsity, dimensionality and sample size. In the ultra sparse regime (where the sparsity s of the true
graph is bounded by s∗ > 0), Theorem 1 provides
convergence inequalities for Ĝ choosing η ≤ √αs∗ with
p log(p) = αn in Assumption H7 .
In the standard sparsity scenario, if s is at least
of the order of p, then η should be of the order of
√
α/ p, which is unrealistic as p → +∞. This case
thus requires a stronger dimensional assumption H5 .
Taking at least p2 log(p) = O(n) ensures a good estimation of the graph.
8

The minimization problem given by Equation (10)
looks like a well-studied problem in machine learning, as it is closely related to the ℓ1 -constrained
quadratic program, known as the Lasso in the statistics literature (Tibshirani, 1996). Indeed, the ℓ1 regularization leads to variable selection and convex constraints that make the optimization problem
solvable. We note here that this allows us to always provide a locally optimal solution, i.e. optimal weight estimates given a hierarchy between the
nodes.
A large number of efficient algorithms are available for computing the entire path of solutions as λ
is varied, e.g. the LARS algorithm of Efron et al.
(2004) and its alternatives. For example, in the context of the estimation of sparse undirected graphical models, Meinshausen & Bühlmann (2006) fit a
Lasso model to each variable, using the others as
predictors, and define some rules for model symmetrization as they do not work on DAGs. The
graphical Lasso (or glasso, Friedman et al. 2007)
algorithm directly relies on the estimation of the
inverse of a structure covariance matrix assumed
to be sparse. Improvements were proposed for
example by Duchi et al. (2008) (improved stopping criterion) and Witten et al. (2011) (estimation of a block-diagonal matrix). Other authors
propose to solve the optimization problem using
an adaptation of classical optimization methods,
such as interior point (Yuan & Lin, 2007) or block
coordinate descent methods (Banerjee et al., 2008;
Friedman et al., 2007).
We propose here an original convex optimization
algorithm to find the solution in Equation (10) in a
form similar to a steepest descent algorithm. Our
proposed algorithm is much quicker than a glasso
approach, a desirable feature as it will run at each
iteration of the global algorithm (see the “Search of
an optimal T ∗ ” box in Figure 3 and the “Evaluate
the new individuals” item in Algorithm 2). Moreover, its mechanistic components (see Section B of
the Supplementary Material Proof details) allowed
us to derive the theoretical results of Theorem 1.
The proposed scheme can be seen as an adaptation

score function, then the roles of P and T are reversed and the procedure is repeated iteratively until
convergence for some criterion (Csiszár & Tusnády,
1984). However, the structure of our problem does
not allow us to use such a scheme: looking for an optimal T given a fixed P makes sense, but changing
P for a fixed T does not.
In our inference algorithm GADAG, an outer loop
is used to perform the global search among the DAGs
space, which is driven by the choice of P , while a
nested loop is used to find an optimal T for each
given fixed P (see Figure 3). As we show in the following, population-based metaheuristics algorithms
are a natural and efficient choice for exploring the
space of permutation matrices (Section 4.3). The
nested optimization problem can be resolved using
a steepest descent approach (Section 4.2).

Choice of P
Search of an optimal T ∗
Evaluate the likelihood
Problem solved?

NO

YES
END
Figure 3:
GADAG.

4.2

Overview of our hybrid algorithm

Graph structure learning when
the variable order is fixed

Assume first that the variable ordering P ∈ Pp (R)
is fixed. The problem of inferring a graph is then
reduced to estimating the graph structure, which can
be solved by finding a solution of:


1
2
T
X(I − P T P ) F + λ kT k1 . (10)
min
T ∈Tp (R)
n
9

of the LARS algorithm with matrix arguments.
Let (Tk )k≥0 the sequence of matrices defined for
all i, j ∈ J1, pK2 as:



λ
j
j
j
(Tk+1 )i = sign (Uk )i max 0, (Uk )i −
, (11)
L
∇



1
n

2

kX(I−P Tk P T )kF

Algorithm 1: Graph structure learning minimization of Equation (10)
Input: λ, L, ǫ > 0.
Initialization: T0 the null squared p × p
matrix, k = 0 and e = +∞.
while e > ǫ do


2
1
∇ n
kX(I−P Tk P T )kF
Compute Uk = Tk −
L
with Equation (13);
Using Equation (11), compute
the current

matrix Tk+1 = (Tk+1 )ji i,j ;
Project Tk+1 on Tp (R) with Equation
(12): Tk+1 ← ProjTp (R) (Tk+1 );
Compute e = kTk+1 − Tk kF ;
Increase k: k ← k + 1;
end
Output: Tk ∈ Tp (R) the unique solution of
(10).



,
where for all k ≥ 0, Uk = Tk −
L
L
is the Lipschitz constant
 of the gradient function
2
1
T
∇ n X(I − P Tk P ) F and sign() is the sign of
any element. Then, a solution of (10) is given by
performing Algorithm 1, where:
• the projection ProjTp (R) (T ) of any p × p realvalued matrix T = ((Tk )ji )i,j on the set Tp (R) is
given by
(

j
0
if i < j,
ProjTp (R) (Tk ) =
(12)
j
i
(Tk )i otherwise.
• the gradient of
∇



1
n

X(I − P Tk P T )

2
F

is


1
2
T
X(I − P Tk P ) F =
n
2
− (XP )T (X − XP Tk P T )P. (13)
n

The detailed calculations are deferred to Section B
of the Supplementary Material Proof details.

4.3

A Genetic Algorithm for a
global exploration of the permutation matrices space incorporating network topologies

As the optimal T can be calculated for any P using Algorithm 1 and with a very good approximation accuracy according to Theorem 1, the optimization task (6) comes down to exploring the Pp (R)
space of permutation matrices in dimension p and
to evaluating the quality of permutation candidates
P ∈ Pp (R). We first note that the number of permutation matrices is p!, which rules out any exact

enumeration method, even for relatively small p. We
propose instead to use a meta-heuristic approach,
which has proven to be successful for many discrete optimization problems like wire-routing, transportation problems or traveling salesman problem
(Michalewicz, 1994; Dréo et al., 2006).
Among the different meta-heuristics (Simulated
annealing, Tabu search, Ant Colony,...) we focused
on Genetic Algorithms (GA) because, despite limited convergence results (Cerf, 1998; Michalewicz,
1994), they were found much more efficient in problems related to ours than alternatives with more established convergence proofs (e.g. Granville et al.
(1994) for simulated annealing), while allowing the
use of parallel computation.
GAs mimic the process of natural evolution, and
use a vocabulary derived from natural genetics: populations (a set of potential solutions of the optimization problem), individuals (a particular solution) and genes (the components of a potential solution). Each generation/iteration of the algorithm
will improve the constituting elements of the population. In short, a population made of N potential
solutions of the optimization problem samples the

10

search space. This population is sequentially modified, with the aim of achieving a balance between exploiting the best solutions and exploring the search
space, until some termination condition is met.
We use here a classical Genetic Algorithm, as described in Michalewicz (1994) for instance, which is
based on three main operators at each iteration: selection, crossover and mutation. The population is
reduced by selection; selection shrinks the population diversity based on the individual fitness values.
The crossover allows the mixing of good properties
of the population to create new composite individuals. Mutations change one (or a few in more general
GAs) components of the individuals to allow random
space exploration. The complete sketch of algorithm
GADAG is given in Algorithm 2. A discussion on
parameters to set in Algorithm 2 is found in Section
5.1. The details of the different operators are given
in the following.
As we show in Example 3, any P ∈ Pp (R) is
uniquely defined by a permutation vector of J1, pK.
Hence, we use as a the search space Sp the set of
permutations of J1, pK, which is a well-suited formulation for GAs.

However, unlike the TSP, the problem is not circular (in the TSP, the last city is connected to the
first one), and the permutation here defines a hierarchy between nodes rather than a path, which
makes the use of TSP-designed operators a potentially poor solution. As we show in the following, we
carefully chose these operators in order to respect
the nature of the problem at hand. In particular,
we emphasize two of their desirable features: their
efficiency in exploring the search space and the interpretable aspect they offer in terms of modifications
on a given network or the blend of two different networks (crossover).
Fitness function Given a potential solution pi ∈
Sp , the fitness function is defined as:
Ji = J(pi ) =

1
X(I − Pi Ti∗ PiT )
n

2
F

+ λ kTi∗ k1 ,

(14)
with Pi constructed from pi as in Example 3 and
Ti∗ the solution of Equation (10) with P = Pi .
As mentioned earlier, at each step of the proposed
GA, the evaluation of the fitness function thus requires running the nested loop of our global algoExample 3 Consider the permutation matrix (p = rithm GADAG.
5):


0 0 0 1 0
Selection operator The selection operator (or
 0 0 0 0 1


survival step) consists in generating a population

P =
 0 1 0 0 0 .
of N individuals from the N existing individuals by
 0 0 1 0 0
random sampling (with replacement, hence some in1 0 0 0 0
dividuals are duplicated and others are deleted). It
aims at improving the average quality of the popula4
Then, P is represented by the 5 3 1 2 vector,
tion by giving to the best potential solutions a higher
looking at the ranks of non-null values of P column
probability to be copied in the intermediate popuby column. The nodes are ranked according to their
lation. We have chosen to use the classical proportopological ordering.
tional selection of Holland (1975): each individuals
Note that our problem closely resembles the clas- is selected with a probability inversely proportional
sical Traveling Salesman Problem (TSP), which has to its fitness value of Equation (14).
been succesfully addressed by means of genetic algorithms (Grefenstette et al., 1985; Davis, 1991).
Identically to the TSP, we optimize over the space
of permutations, which induces specific constraints
for defining the crossover and mutation operators.

Crossover operator A crossover operator generates a new set of potential solutions (children)
from existing solutions (parents). Crossover aims at
achieving at the same time (i) a good exploration of

11

the search space by mixing the characteristics of the
parents to create potentially new ones while (ii) preserving some of the desirable characteristics of the
parents. By desirable features, we mean features of
the network which lead to good fitness values, and
which in turn are favored by selection over the generations. The crossover population (set of parents)
is obtained by selecting each individual of the population with a probability pxo ; the parents are then
paired randomly.
We have chosen the order-based crossover, originally proposed for the TSP (Michalewicz, 1994,
Chapter 10), which is defined as follows. Given two
parents p1 and p2 , a random set of crossover points
are selected, which we denote Ω. It consists in a
permutation of k elements taken from J1, pK, with k
uniformly drawn between 0 and p. A first child C1
between p1 and p2 is then generated by:
1. fixing the crossover points of p1 ,

GP1 and GP2 , constructs two new graphs. One of
them, GC1 is composed of the sub-graph of GP1 induced by the set of crossover points Ω and the subgraph of GP1 induced by the complementary set ΩC
of Ω in J1, pK (see Figure 4). The second child graph
GC2 is obtained in an identical manner by reversing
the roles played by the two parental graphs.
Mutation operator Mutation operators usually
correspond to the smallest possible change in an individual (unary operator). We thus define it as an
alteration of two neighbouring genes (see Example
5). Graphically, a mutation consists in switching the
arrowhead of an edge between two nodes. Mutation
is applied to each child with probability pm .
Example 5 A possible mutation for the first child
of Example 4 is to swap the genes “1” and “10” (in
bold red below):
M1

2. completing C1 with the missing numbers in the
order they appear in p2 .
Example 4 Consider the two following parents:
p1

4 3 10 7 5 9 1 2 6 8

p2

6 1 9 4 10 2 8 3 7 5

4 6 1 10 3 9 7 2 5 8

Stopping criterion Two quantities are monitored along the iterations: the heterogeneity of the
population and the value of the objective function.
For the first indicator, we use the Shannon entropy, defined for each rank position j ∈ J1, pK as:



p
X
Ni,j
Ni,j
Assume that the crossover points randomly chosen
,
log
Hj = −
N
N
are 4, 9, 2 and 8 (in bold red above). Then, the child
i=1
C1 is defined by inheriting those points from p1 and
filling the other points in the order they appear in p2 : where Ni,j is the number of times when i appears
in position j. Hj = 0 if all the individuals “agree”
p2 6 1 /9 /4 10 /2 /8 3 7 5
on the position of a node and the population is perfectly homogeneous at this node. On the contrary,
it is maximum when we observe a uniform distribuC1 4 * * * * 9 * 2 * 8
tion of the different nodes at a given position and the
population is in this case at a maximum of hetero⇓
geneity or disorder for this position. The algorithm
P
4 6 1 10 3 9 7 2 5 8
stops if the population entropy value H = N
j=1 Hj
drops below a threshold since H = 0 if all the inFrom a graphical point of view, a crossover be- dividuals are identical. A second criterion can tertween p1 and p2 , which encode two complete graphs minate GADAG if difference in the average fitness
12

10

4

10

4
9

9

5

9

5

5

1

2
6

1

2
6

3

1

2
6

8
7

10

4

8
7

3

8
7

3

Figure 4: Graphical representation of crossover between two 10-node graphs. The two parental graphs are
represented on the left. The third graph, on the right, is obtained by combining the blue and red part of
its parents using the crossover operator.
(denoted J¯ thereafter) of the population between
a given number of consecutive iterations, does not
change by more than a predefined threshold.

5

Numerical experiments

Algorithm 2: GADAG overview
Input: pxo , pm , ǫH > 0, ǫJ > 0, kmax > 0,
imax > 0, λ, L.
Initialization: Generate the initial
population P0 with N permutations of J1, pK,
k = 0 and eJ = +∞.
while H > ǫH & eJ > ǫJ & k ≤ kmax do
Generate Pk+1 as a random selection of
N individuals from Pk ;
Pick an even subset Pxo of Pk+1 (each
individual of Pk+1 selected with
probability pxo );
Perform crossover on Pxo by randomly
pairing the individuals;
Mutate each obtained individual with
probability pm ;
Evaluate the new individuals Pm by
running Algorithm 1;
Replace Pxo by Pm in Pk+1 ;
Compute the Shannon entropy H and the
difference in the average fitness

¯ k+1 ) − J(P
¯ k−i ) ;
eJ = max0≤i≤imax J(P
Increase k: k ← k + 1;
end

This section is dedicated to experimental studies to
assess practical performances of our method through
two kinds of datasets. In a first phase, the aim of
these applications is to show that GADAG has a
sound behavior on simulated toy data with a variety of different settings. In a second phase, we
demonstrate the ability of our algorithm to analyse
datasets that mimic the activity of a complex biological system, and we compare it to other state-of-the
art methods. The competing methods are presented
in Section 5.4.1. In Section 5.1, we present the calibration of the Genetic Algorithm parameters. Section 5.2 introduces the measures we used to assess
the merits of the methods. Experimental results are
then detailed in Section 5.3 for the simulated highdimensional toy datasets and in Section 5.4.2 for the
dataset with features encountered in real situations.
All experiments have been performed on R
(R Core Team, 2017) using the package GADAG
(Champion et al., 2017). The computational times
reported in Section 5.3 correspond to a Windows 7
laptop computer with 8 threads on a 4-core hyper- threaded 2.50GHz processor, with 4GB of RAM.
13

5.1

Algorithm parameters

Table 1: Algorithm parameter settings
Parameter
Value

Running the procedure of Algorithm 2 requires to
define parameters of the outer loop, which generpxo
0.25
ates our population of P ’s, and of the nested loop
p
0.5
m
to find the optimal T ∗ . The evaluation of the LipsN
5×p
chitz gradient constant L, used to find the optimal
2
L
XT X F
n
graph structure T ∗ , is known as a hard established
max. nb. of eval. 104
problem in optimization. Some authors propose to
ǫH
10−6
choose an estimate of L from a set of possible values
ǫJ
10−4
(Jones et al., 1993; Sergeyev & Kvasov, 2006), to estimate local Lipschitz constants (Sergeyev, 1995), or
to set it a priori to a fixed value (Evtushenko et al.,
2009; Horst & Pardalos, 1995). Here, observing by TP, FP, FN and TN, the true positive (correctly
Equation (13), a major bound for L is given by:
predicted) edges, the false positive (inferred by mistake) edges, the false negative (missed) edges, and
2
L≤
XT X F .
the true negative (correctly non-predicted) edges.
n
P
The recall, defined as T PT+F
, measures the power
N
We found that setting L to this bound worked well (or sensitivity) of reconstruction of non-zero elein practice in all our scenarios.
ments of the true matrix G (or equivalently of the
Five parameters need to be tuned to run the Ge- true network) for one method, whereas the precinetic Algorithm: the crossover rate pxo , the muta- sion, equal to T P , measures the accuracy of the
T P +F P
tion rate pm , the constant of the stopping criteria reconstruction. The closer to one the precision and
ǫH and ǫJ and the size of the population N. For the recall the better.
the first four parameters, we observed that their
P/R curves represent the evolution of those quanvalue had a limited effect on the efficiency, hence tities when varying the sparsity of the methods.
we chose commonly used values in the literature GADAG is based on penalized optimization: it seeks
(see Table 1). The size of the population has a linear dependencies between the variables with a
more complex effect and has been investigated in controlled level of parsimony (λ in Equation (5)).
several prospective papers (e.g. Schaffer et al. 1989; For λ varying from 0 (complete graph) to +∞
Alander 1992; Piszcz & Soule 2006; Ridge 2007) but (empty graph), it thus produces a list of edges sucwithout providing a definitive answer to the prob- cessively introduced in the model. This list of edges
lem. In our simulation study, we chose as a rule-of defines the precision versus recall curve. As a sumthumb N = 5p, which was found as a good compro- mary performance measurement, we also computed
mise between computational cost and space explo- the classical area under the P/R curve (AUPR).
ration on several experiments.
The complete parameter settings used in our experiments are reported in Table 1.
5.3 Exploratory analysis on toy

datasets
5.2

Performance metrics

A classical performance measure for graph inference
methods consists in comparing predicted interactions with the known edges in the true graph G0 using precision versus recall (P/R) curves. We denote

We first considered simulated data from randomly
generated networks with different characteristics in
order to assess the capabilities and limits of our algorithm. Given a number of nodes p, a random set
of s edges were generated, and the non-zero param-

14

eters of the matrix G0 associated to the corresponding DAG were uniformly sampled between 0 and 1.
Using this graph, we generated N observations following the hypotheses of Gaussian, homoscedastic
and centred error. We then ran GADAG on this
dataset to recover the graph. Note that other assumptions presented in Section 3.2 may not be fulfilled here, but we aimed at evaluating the robustness
of GADAG for recovering DAGs in such scenario.
In our experiments, we varied the number of nodes
p, of edges s and of available observations n. We
chose four different settings p = 50, 100, 500 and
1, 000 with n/p varying from 100% to 10% and s/p
from 100% to 400%. Unless otherwise stated, all
experiments were replicated 50 times each and results were averaged over these replicates. Averaged
computational times correspond to one iteration of
GADAG, for a fixed parameter of penalization λ.
Results, in terms of area under the P/R curves and
computational time are sumarized in Table 2. We
can first remark a crude decrease of performance results when the number of samples is very small (p =
50 and 100, n/p = 10%, so respectively 5 and 10
samples). In that case, GADAG is incapable of recovering any signal (AUPR < 10%). When the sample size is of the order of p (n/p = 100%, first row of
Table 2 a), GADAG works well, although it is clearly
a favorable case, far from the high-dimensional one.
With half of the samples (n/p = 50%), performance
remains satisfactory (AUPR around 50%, or more),
which is critical since this situation corresponds to
realistic biological studies, where subsets of genes
(i.e. nodes) are preselected beforehand. Interestingly, p = 500 and 1, 000 work better than smaller
values of p since the number of samples is larger to
estimate the graph that generated the data. Indeed,
for a given number of samples, e.g. n = 50, performance results slightly decrease from 65% (p = 50)
and 55% (p = 100) to 50% (p = 500). GADAG is
thus not considerably affected by a relative increase
of dimensionality. For large graphs (p = 500 and
1, 000), even if n/p ≤ 10%, it suceeds in recovering
them, which makes it a competitive algorithm with
regards to other state-of-the-art approaches.

An interesting remark is that the number of edges
s does not significantly change numerical results (see
each row of Table 2 a), although GADAG succeeds
slightly better in estimating sparser graphs. This
may be due to the particular structure of our algorithm, which looks for topological ordering between
nodes (genetic algorithm) and then makes the inferred graph sparse.
Concerning computational time, we can finally
note that growing the dimension p clearly makes
the problem harder to solve : each call to GADAG
requires more than 300s for hundred of nodes and
1, 500s for thousand of nodes.

5.4

DREAM data analysis

The second type of datasets we used mimic activations and regulations that occur in gene regulatory networks. It is provided by the DREAM4
challenge on “In Silico Network Challenge”. Note
that although plausibly simulated, DREAM4 data
sets are not real biological data sets. However, the
used network structures (five in total) were extracted
from E. coli and S. cerevisae -two biological model
organisms- trancriptional networks. These networks
contain cycles, but self-loops were discarded. The
gene expression observations were not simulated by
an equal noise Gaussian multivariate model, stochastic differential equations were used to mimic the kinetic laws of intricate and intertwined gene regulations. In addition to the biological noise simulated
from the stochastic differential equations, technical
noises were added to reproduce actual gene measurement noise. All data sets were generated by the GNW
software (Marbach et al., 2009).
Working with simulated networks, we are able to
quantitatively and objectively assess the merit of
competing methods in terms of true predictions (true
positive TP and true negative TN) vs. incorrect predictions (false positive FP and false negative FN)
edges. While the analysis of a real data set is certainly the final goal of a methodology motivated by a
real problem like ours, there are only imprecise ways
of validating a method when analysing a real data

15

Table 2: Performance results of GADAG on a toy dataset with different characteristics (number of nodes
p, number of edges s and sample size n) in terms of area under the Precision vs. Recall curve (a) and
computational time, in seconds (b). All results are averaged over 50 replicates (∗ 5 replicates only as the
running time was 1 day per network).
(a)
p=50

n/p

p=100

p=1,000∗

p=500

s/p

100%

200%

400%

100%

200%

400%

100%

200%

400%

100%

200%

400%

100%
50%
10%

0.72
0.53
0.02

0.65
0.46
0.03

0.65
0.49
0.05

0.69
0.56
0.04

0.63
0.53
0.06

0.64
0.53
0.09

0.79
0.75
0.51

0.79
0.76
0.53

0.70
0.66
0.41

0.86
0.82
0.65

0.80
0.75
0.60

0.65
0.63
0.49

(b)
p=50

n/p

p=100

s/p

100%

200%

400%

100%

200%

400%

100%

200%

400%

100%

200%

400%

100%
50%
10%

0.55
0.48
0.31

0.43
0.48
0.30

0.41
0.41
0.36

8.69
8.63
6.52

9.87
8.90
5.81

7.51
6.62
5.41

253
250
183

258
258
183

214
172
165

1,640
1,520
1,590

1,500
1,590
1,550

1,230
1,270
1,290

set. Well known systems are often small and even
if knowledge has accumulated on them, these can
be noisy and difficult to gather to obtain a fair picture of what can adequately be considered as sets of
true positive and true negative sets of edges. Even if
the data generation process of the DREAM4 In Silico Network Challenge is completely understood, no
existing method is able to predict all regulatory relationships, but at the price of including many false
positive predictions. The DREAM4 datasets we considered have p = 100 nodes and only n = 100 observations making it a a very challenging task.
5.4.1

p=1,000∗

p=500

Comparison to state-of-the art

We compare GADAG to five state-of-the-art inference methods. Among them, the Genie3 method
(Huynh-Thu et al., 2010), based on random forests,
was the best performer of one of the DREAM4
sub-challenges, while the BootLasso (Allouche et al.,
2013) was one of the key components of the best

performing approach of one of the DREAM5 subchallenges (Allouche et al., 2013). The two methods decompose the prediction of the network into
p feature selection sub-problems. In each of the p
sub-problems, one of the node is predicted from the
other ones using random forests (Breiman, 2001)
for Genie3 or a bootstrapped version of the Lasso
(Bach, 2008) for BootLasso. For the random forest approach, parents of each node were detected as
most significant explanatory variables according to
a variance reduction criterion in a regression tree
framework. The process was repeated on a randomized set of trees, which made up the so-called
random forest. This method allowed us to derive
a ranking of the importance of all variables for the
target by averaging the scores over all the trees. We
used the R package randomforest (Liaw & Wiener,
2002) for our results. The Lasso is a ℓ1 -norm penalization technique for solving linear regression. Following the works of Bach (2008), BootLasso uses

16

bootstrapped estimates of the active regression set
based on a Lasso penalty: only those variables that
are selected in every bootstrap are kept in the model.
In both cases, actual coefficient values are estimated
from a straightforward least square procedure. Note
that we slightly relax the condition for a variable
to be included in the model, a variable was selected
at a given penalty level if more than 80% of bootstrapped samples led to selecting it in the model
(Allouche et al., 2013).
We also compare our algorithm to three classical methods for Bayesian Networks (BNs) modelling.
BNs are graphical models (Pearl, 2009) defined by
a DAG and parameters that set quantitative relationships between variables. Algorithms devoted to
structure and parameter learning in BNs either aim
at maximising a score that reflects the fit of the data
to the learnt structure, or test for independencies
between variables. They are often used as references in a gene regulatory network inference context (Tsamardinos et al., 2006), athough mainly for
moderate size networks. The first compared algorithm we used is the PC-algorithm (Spirtes et al.,
2000), a popular constraint-based method that drastically reduces the number of conditional independence tests. It first builds the skeleton of the graph
by removing edges from a complete undirected graph
before determining the orientation of the edges,
when possible. A large number of implementations
of the PC-algorithm exists. The numerical results
presented here were obtained using the pcAlg function of the R-package pcalg, based on standard
correlation estimates for conditional indepence testing. We also ran ARACNE (Margolin et al., 2006),
an improved version of minimum-weight spanning
tree that uses the information inequality to eliminate the majority of indirect relationships between
variables. We used the ARACNE function of the
R-package bnlearn. We finally compare GADAG
to the Greedy Equivalence Search (GSE) algorithm
(Chickering, 2002), implemented in the R-package
pacalg, which heuristically searches in the space
of equivalent classes the model with the highest
Bayesian score.

To compare our algorithm with these competing
methods, we used the P/R curves presented in Section 5.2. As GADAG, BootLasso leads to a sparse
inferred graph while controlling the level of parsimony, which builds the P/R curve. Genie3 produces
as an output a ranked list of regulatory interactions,
which corresponds to the edges of the inferred graph.
Edges are the successively introduced with decreasing confidence scores to produce the random forest
P/R curve. For the PC and the GSE algorithms,
inherent parameters regulating the sparsity of the
produced graphs helped us to define such curves.
Note that the implementation we used for running
ARACNE was only able to produce a final network
prediction (interaction ranking is not available).
5.4.2

Numerical results

The P/R curves for the five DREAM problems are
shown in Figure 5. Each curve corresponds to one
of the five networks used in the challenge. In general, for all the problems the five methods are able to
achieve a precision equal to one (that is, to include
only true edges), but these correspond to overly
sparse graphs (very small recall). Conversely, a recall equal to 1 can only be reached by adding a large
number of FP edges, whatever the method we consider, even if some fail earlier than others. The main
differences between the methods appear on the leftmost part of the P/R curves, especially those of Figure 5 B, C and D: while the precision of BootLasso,
Genie3, PCalg and GSE drops rapidly with a slow
increases in recall above 20% recall, it remains higher
for GADAG. Hence, its first predicted edges are at
least as accurate than those of the four other methods and it produces a larger set of reliable edges. For
graphs of lesser sparsity, none of the five methods is
really able to identify clearly reliable edges. Large
number of FP edges are produced to achieve a recall
higher than 60%.
For Networks 1 and to a lesser extent 5 (Figure
5 A and E), GADAG recovers with more difficulty
the first true edges than other methods, with a high
level of FP edges at the beginning of the curve (low

17

clearly appeared in the numerical experiments, for
which the behavior of our algorithm seemed to be
very competitive when compared to the state-of-theart.
Nevertheless, we see many opportunities for further improvements. First, convergence proof for
the algorithm, although a challenging task, is worth
investigating, for instance using the works of Cerf
(1998) on genetic algorithms. An alternative would
Table 3: Area under the Precision vs. Recall curve be to consider other optimization schemes for the
node ordering with more established convergence
for all networks and methods (except ARACNE).
proofs (e.g. simulated annealing (Granville et al.,
1994)).
Method
Net 1 Net 2 Net 3 Net 4 Net 5 Second, other potential extensions involve algoGADAG
0.182 0.236 0.348 0.317 0.267 rithmic considerations in order to improve the calcuGenie3
0.154 0.155 0.231 0.208 0.197 lation time, including a finer calibration of the algoBootLasso 0.118 0.061 0.171 0.147 0.169 rithm parameters, an initialization step for the graPCalg
0.116 0.089 0.171 0.149 0.130 dient descent, and, in general, improving the interacGSE
0.101 0.089 0.170 0.153 0.133 tions between the nested and outer loops. Tackling
very large datasets from several thousands of nodes
may also require a particular treatment, for instance
by adding local search operators to GADAG.
Finally, we would like to emphasize the graph
identifiability
problem: in our settings, we assume
6 Conclusion and discussion
the noise variances of all graph nodes to be equal to
In this paper, we proposed a hybrid genetic/convex ensure graph identifiability (that is no equivalence
algorithm for inferring large graphs based on a par- class of graphs). Such a hypothesis is of course reticular decomposition of the ℓ1 -penalized maximum strictive and likely to be violated for real datasets.
likelihood criterion. We obtained two convergence In order to infer networks for any noise variances,
inequalities that ensure that the graph estimator one solution consists in incorporating interventional
converges to the true graph under assumptions that data on the model. These data are obtained from
mainly control the model structure: graph size (bal- perturbations of the biological system (e.g. gene
ance between sparsity, number of nodes and max- knockouts or over-expressions) and make the equivimal degree) and signal-to-noise ratio. From an alence class of graphs smaller (Hauser & Bühlmann,
algorithmic point of view, the estimation task is 2012). The use of additional data, informative yet
split into two subproblems: node ordering estima- very costly interventional data could be combined
tion and graph structure learning. The first one is with observational on the MLE estimator. It was rea non-trivial problem since we optimize over a dis- cently proposed by Hauser & Bühlmann (2015) for
crete non-convex large dimensional set. It led us to a BIC-score penalized MLE, or by Rau et al. (2013)
use a heuristic approach we specifically tailored to for learning Gaussian Bayesian networks in the case
achieve the optimization task. The second one is a of GRN inference. A modification of our hybrid almore common problem, related to the Lasso one, for gorithm GADAG could then lead to a more accuwhich we proposed a sound procedure with theoret- rate identification of the true graph. Lastly, the
ical guarantees. The potential of such an approach cyclic structure framework could also be considprecision and low recall). However, as soon as the
recall exceeds the 10%, resp. 15%, for graph A, resp.
for graph E, GADAG performance is again superior
to that other methods.
Table 3 gives the areas under the P/R curves for
all methods and networks. For this indicator, GA
significantly outperforms the state-of-the-art methods for all networks.

18

0.6
0.4

Precision

0.8

GADAG
ARACNE
BootLasso
Genie3
PCalg
GSE

0.0

0.0

0.2

0.4

0.6

0.8

GADAG
ARACNE
BootLasso
Genie3
PCalg
GSE

0.2

Precision

1.0

B

1.0

A

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

Recall

0.8

1.0

1.0

D
0.6
0.2

Precision

0.8

GADAG
ARACNE
BootLasso
Genie3
PCalg
GSE

0.0

0.0

0.2

0.4

0.6

0.8

GADAG
ARACNE
BootLasso
Genie3
PCalg
GSE

0.4

1.0

C
Precision

0.6
Recall

0.0

0.2

0.4

0.6

0.8

1.0

0.0

Recall

0.2

0.4

0.6

0.8

1.0

Recall

1.0

E
0.6
0.4
0.0

0.2

Precision

0.8

GADAG
ARACNE
BootLasso
Genie3
PCalg
GSE

0.0

0.2

0.4

0.6

0.8

1.0

Recall

Figure 5: P/R curves for the five Dream networks and the five compared methods.
ered by using a Markov equivalence characterization
(Richardson, 1997) to relax the strictly triangular assumption on our matrix T using Equation of Proposition 2.1. It would pave the way for totally new
theoretical developments, and a more realistic modelling of a complex system.

(2013). A panel of learning methods for the
reconstruction of gene regulatory networks in a
systems genetics. In Gene Network Inference,
Verification of Methods for Systems Genetics
Data. Springer, Heidelberg, alberto de la fuente
ed., pp. 9–32.

Aragam, B., Amini, A. & Zhou, Q. (2015).
Learning directed acyclic graphs with penalized
References
neighbourhood regression. Preprint on arxiv at
https://arxiv.org/pdf/1511.08963.pdf.
Alander, J. (1992). On optimal population size
of genetic algorithms. In Proceedings of the IEEE Bach, F. (2008). Bolasso: model consistent lasso
Computer Systems and Software Engineering.
estimation through the bootstrap. In Proceedings
of the Twenty-fifth International Conference on
Allouche,
A.,
Cierco-Ayrolles,
C.,
Machine Learning.
de Givry, S., Guillermin, G., Mangin,
B., Schiex, T., Vandel, J. & Vignes, M. Banerjee, O., El Ghaoui, L. & d’Aspremont,
19

A. (2008). Model selection through sparse maxi- Cook, S. (1985). A taxonomy of problems with fast
parallel algorithms. Inform Control 64, 2–22.
mum likelihood estimation for multivariate Gaussian or binary data. J Mach Learn Res 9, 485–516.
Cormen, T. H., Leiserson, C. E., Rivest,
R. L. & Stein, C. (2001). Introduction to algoBarabási, A. & Oltvai, Z. (2004). Network biolrithms. MIT Press, Cambridge, MA; McGraw-Hill
ogy: understanding the cell’s functional organizaBook Co., Boston, MA, 2nd ed.
tion. Nat Rev Genet 5, 101–113.
Bickel, P. & Li, B. (2006). Regularization in Csiszár, I. & Tusnády, G. (1984). Information geometry and alternating minimization prostatistics. Test 15, 271–344.
cedures. Statist. Decisions , 205–237Recent results in estimation theory and related topics.
Bickel, P. J., Ritov, Y. & Tsybakov, A. B.
(2009). Simultaneous analysis of lasso and Dantzig
Davis, L., ed. (1991). Handbook of Genetic Algoselector. Ann Stat 37, 1705–1732.
rithms. Van Nostrand Reinhold.
Breiman, L. (2001). Random forests. Mach Learn
De Smet, R. & Marchal, K. (2010). Advantages
45, 532.
and limitations of current network inference methods. Nature Reviews Microbiology 8, 717–729.
Bühlmann, P. (2013). Causal statistical inference
in high dimensions. Math Method Oper Res 77, Dondelinger, F., Lèbre, S. & Husmeier, D.
357–370.
(2013). Non-homogeneous dynamic bayesian networks with bayesian regularization for inferring
Bühlmann, P. & van de Geer, S. (2011). Statisgene regulatory networks with gradually timetics for high-dimensional data. Springer Series in
varying structure. Machine Learning 90, 191–230.
Statistics. Springer, Heidelberg. Methods, theory
and applications.

Dréo, J., Pétrowski, A., Siarry, P. & Taillard, E. a. (2006). Metaheuristics for hard
Cerf, R. (1998). Asymptotic convergence of geoptimization. Springer-Verlag, Berlin. Methods
netic algorithms. Adv Appl Probab 30, 521–550.
and case studies, Translated from the 2003 French
original by Amitava Chatterjee.
Champion, M., Picheny, V. & Vignes, M.
(2017). GADAG: A Genetic Algorithm for learn- Duchi, J., Gould, S. & Koller, D. (2008).
ing Directed Acyclic Graphs. R package version
Projected subgradient methods for learning sparse
0.99.0.
gaussians. In Proceedings of the 24th Conference
on Uncertainty in Artificial Intelligence.
Chen, J. & Chen, Z. (2008). Extended bayesian
information criteria for model selection with large Efron, B., Hastie, T., Johnstone, I. & Tibmodel spaces. Biometrika , 759–771.
shirani, R. (2004). Least angle regression. Ann
Stat 32, 407–499.
Chickering, D. M. (1996). Learning Bayesian
networks is NP-complete. In Learning from data El Karoui, N. (2008). Spectrum estimation for
(Fort Lauderdale, FL, 1995), vol. 112 of Lect.
large dimensional covariance matrices using ranNotes Stat. Springer, New York, pp. 121–130.
dom matrix theory. Ann Stat. 36, 2757–2790.
Chickering, D. M. (2002). Optimal structure Ellis, B. & Wong, W. (2008). Learning causal
Bayesian network structures from experimental
identification with greedy search. J Mach Learn
data. J Am Stat Assoc 103, 778–789.
Res 3, 507–554.
20

Evtushenko,
Y.,
Malkova,
V.
& Hauser, A. & Bühlmann, P. (2015). Jointly interventional and observational data: estimation
Stanevichyus, A. (2009).
Parallel global
of interventional Markov equivalence classes of dioptimization of functions of several variables.
rected acyclic graphs. J R Stat Soc. Series B 77,
Comput. Math. Math. Phys. 49, 246–260.
291–318.
Friedman, J., Hastie, T. & Tibshirani, R.
(2007). Sparse inverse covariance estimation with
Hogben, L., ed. (2007). Handbook of linear algebra.
the graphical lasso. Biostatistics 9, 432–441.
Discrete Mathematics and its Applications (Boca
Raton). Chapman & Hall/CRC, Boca Raton, FL.
Friedman, N. & Koller, D. (2003). Being
Associate editors: Richard Brualdi, Anne GreenBayesian about network structure: a Bayesian
baum and Roy Mathias.
approach to structure discovery in Bayesian networks. Mach Learn 50, 95–125.
Holland, J. H. (1975). Adaptation in natural and
Fu, F. & Zhou, Q. (2013). Learning sparse causal
artificial systems. University of Michigan Press,
Gaussian networks with experimental intervenAnn Arbor, Mich. An introductory analysis with
tion: regularization and coordinate descent. J Am
applications to biology, control, and artificial inStat Assoc 108, 288–300.
telligence.
Giraud, C. (2015).
Introduction to highdimensional statistics, vol. 139 of Monographs on Horst, R. & Pardalos, P. M., eds. (1995).
Handbook of global optimization. Nonconvex optiStatistics and Applied Probability. CRC Press,
mization and its applications. Dordrecht, Boston:
Boca Raton, FL.
Kluwer Academic Publishers.
Granville, V., Krivanek, M. & Rasson, J.-P.
(1994). Simulated annealing: A proof of conver- Hoyle, R., ed. (1995). Structural equation modelgence. IEEE T Pattern Anal 16, 652–656.
ing. Thousand Oaks [u.a.]: Sage Publ.
Grefenstette, J., Gopal, R., Rosmaita, B. &
Van Gucht, D. (1985). Genetic algorithms for Husmeier, D. & Werhli, A. (2007). Bayesian
integration of biological prior knowledge into the
the traveling salesman problem. In Proceedings
reconstruction of gene regulatory networks with
of the first International Conference on Genetic
bayesian networks. Comput Syst Bioinformatics
Algorithms and their Applications.
6, 85–95.
Grzegorczyk, M. & Husmeier, D. (2008). Improving the structure MCMC sampler for bayesian Huynh-Thu, V. A., Irrthum, A., Wehenkel,
networks by introducing a new edge reversal move.
L. & Geurts, P. (2010). Inferring regulatory
Mach Learn 71.
networks from expression data using tree-based
methods. PLoS ONE 5, e12776.
Guyon, I., Aliferis, C. & Cooper, G. (2010).
Causation and Prediction Challenge, vol. 2. MiJones, D. R., Perttunen, C. D. & Stuckman,
crotome Publishing. Challenges in Machine LearnB. E. (1993). Lipschitzian optimization without
ing.
the lipschitz constant. J Optimiz Theory App 79,
157–181.
Hauser, A. & Bühlmann, P. (2012). Characterization and greedy learning of interventional
Markov equivalence classes of directed acyclic Kahn, A. B. (1962). Topological sorting of large
graphs. J Mach Learn Res 13, 2409–2464.
networks. Commun. ACM 5, 558–562.
21

Kahre, K., Oh, S. & Rajaratnam, B. (2015). Meinshausen, N. & Bühlmann, P. (2006). Highdimensional graphs and variable selection with the
A convex pseudolikelihood framework for high dilasso. Ann Stat 34, 1436–1462.
mensional partial correlation estimation with convergence guarantees. JRSS B 77, 803–825.
Mestre, X. (2008). Improved estimation of eigenKalisch, M. & Bühlmann, P. (2007). Estimating
values and eigenvectors of covariance matrices ushigh-dimensional directed acyclic graphs with the
ing their sample estimates. IEEE T Inform ThePC-algorithm. J Mach Learn Res 8, 613–636.
ory 54, 5113–5129.
Koivisto, M. & Sood, K. (2004). Exact bayesian
Michalewicz, Z. (1994). Genetic algorithms +
structure discovery in Bayesian networks. J Mach
data structures = evolution programs. SpringerLearn Res 5, 549–573.
Verlag, Berlin, 2nd ed.
Ledoit, O. & Wolf, M. (2015). Spectrum estimation: A unified framework for covariance ma- Mordelet, F. & Vert, J. (2008). Sirene: Supervised inference of regulatory networks. Bioinfortrix estimation and pca in large dimensions. J
matics 24, i76–i82.
Multivariate Anal 139, 360384.
Liaw, A. & Wiener, M. (2002). Classification and Newman, M. (2003). The structure and function
of complex networks. SIAM Rev 45, 157–256.
regression by randomforest. R News 2, 18–22.
Liu, H., Wang, L. & Zhao, T. (2014). Sparse Pearl, J. (2009). Causality. Cambridge University
covariance matrix estimation with eigenvalue conPress, Cambridge, 2nd ed. Models, reasoning, and
straints. J Comput Graph Stat 23, 439–459.
inference.
Lounici, K., Pontil, M., Tsybakov, A. & Perrin, B., Ralaivola, L., Mazurie, A., Botvan de Geer, S. (2009). Taking advantage of
tani, S., Mallet, J. & d’Alché Buc, F.
sparsity in multi-task learning. In Proceedings of
(2003). Gene networks inference using dynamic
the 22nd Conference on Learning Theory.
bayesian networks. Bioinformatics Supp2, ii138–
148.
Maathuis, M., Colombo, D., Kalisch, M. &
Bühlmann, P. (2010). Predicting causal effects
Peters, J., Mooij, J., Janzing, D. &
in large-scale systems from observational data.
Schölkopf, B. (2014). Causal discovery with
Nature Methods 7, 247–248.
continuous additive noise models. J Mach Learn
Res 15, 2009–2053.
Marbach, D., Schaffter, T., Mattiussi, C.
& Floreano, D. (2009). Generating realistic in
silico gene networks for performance assessment of Peters, J., Mooij, J. M., Janzing, D. &
Schölkopf, B. (2011). Identifiability of causal
reverse engineering methods. J Comput Bio 16,
graphs using functional models. In 27th Confer229–239.
ence on Uncertainty in Artificial Intelligence (UAI
Margolin, A., Nemenman, I., Basso, K., Wig2011).
gins, C., Stolovitzky, G., Dalla Favera,
R. & Califano, A. (2006). Aracne: An algo- Piszcz, A. & Soule, T. (2006). Genetic programming: Optimal population sizes for varying comrithm for the reconstruction of gene regulatory
plexity problem. In Proceedings of the Genetic and
networks in a mammalian cellular context. BMC
Evolutionary Computation Conference.
Bioinformatics 7, S7.
22

R Core Team (2017). R: A Language and Envi- Shojaie, A., Jauhiainen, A., Kallitsis, M.
& Michailidis, G. (2014). Inferring regulatory
ronment for Statistical Computing. R Foundation
networks by combining perturbation screens and
for Statistical Computing, Vienna, Austria.
steady state gene expression profiles. PLoS ONE
Rau, A., Jaffrézic, F. & Nuel, G. (2013). Joint
9, e82393.
estimation of causal effects from observational and
intervention gene expression data. BMC Syst Biol Shojaie, A. & Michailidis, G. (2010). Penalized
likelihood methods for estimation of sparse high7, 111.
dimensional directed acyclic graphs. Biometrika
97, 519–538.
Richardson, T. (1997). A characterization of
markov equivalence for directed cyclic graphs. Int
Silander, T. & Myllymäki, T. (2006). A simple
J Approx Reason 17, 107–162.
approach for finding the globally optimal bayesian
network structure. In Proceedings of the TwentyRidge, E. (2007). Design of Experiments for the
second Conference on Uncertainty in Artificial InTuning of Optimisation Algorithm. Ph.D. thesis,
telligence.
The University of York, Department of Computer
Science.

Souma, W., Fujiwara, Y. & Aoyama, H.
(2006). Heterogeneous economic networks. In
Robinson, R. W. (1973). Counting labeled acyclic
The complex networks of economic interactions,
digraphs , 239–273.
vol. 567 of Lecture Notes in Econom. and Math.
Systems. Springer, Berlin, pp. 79–92.
Sachs, K., Itani, S., Fitzegerald, J., Wille,
L., Schoeberl, B., Dahleh, M. & Nolan, G. Spirtes, P. (1995). Directed cyclic graphical rep(2009). Learning cyclic signaling pathway strucresentations of feedback models. In Proceedings
tures while minimizing data requirements. In Proof the 11th Annual Conference on Uncertainty in
ceedings of the Pacific Symposium on BiocomputArtificial Intelligence.
ing.
Spirtes, P., Glymour, C. & Scheines, R.
(2000). Causation, prediction, and search. AdapSchaffer, J. D., Caruana, R., Eshelman,
tive Computation and Machine Learning. CamL. J. & Das, R. (1989). A study of control pabridge, MA: MIT Press, 2nd ed. With additional
rameters affecting online performance of genetic
material by David Heckerman, Christopher Meek,
algorithms for function optimization. In ProceedGregory F. Cooper and Thomas Richardson, A
ings of the Third International Conference on GeBradford Book.
netic Algorithms.
Schwarz, G. E. (1978). Estimating the dimension
of a model. The Annals of Statistics 6, 461–464.

Tibshirani, R. (1996). Regression shrinkage and
selection via the lasso. J R Stat Soc Series. B 58,
267–288.

Sergeyev, Y. (1995). An information global op- Tsamardinos, I., Brown, L. & Aliferis, C.
timization algorithm with local tuning. SIAM J
(2006). The max-min hill-climbing bayesian netOptimiz 5, 858–870.
work structure learning algorithm. Mach Learn
65, 31–78.
Sergeyev, Y. & Kvasov, D. (2006). Global
search based on efficient diagonal partitions and van de Geer, S. & Bühlmann, P. (2013). ℓ0 penalized maximum likelihood for sparse directed
a set of lipschitz constants. SIAM J Optimiz 16,
acyclic graphs. Ann Stat 41, 536–567.
910–937.
23

Verma, T., Araújo, N. & Herrmann, H.
(2014). Revealing the structure of the world airline network. Scientific Reports 5.
Verma, T. & Pearl, J. (1991). Equivalence and
synthesis of causal models. In Proceedings of the
6th Annual Conference on Uncertainty in Artificial Intelligence.
Verzelen, N. (2012). Minimax risks for sparse regressions: Ultra-high dimensional phenomenons.
Electron J Stats 6, 38–90.
Wainwright, M. (2009). Information-theoretic
limits on sparsity recovery in the high-dimensional
and noisy setting. IEEE T Info Theory 55, 5728–
5741.
Witten, D., Freidman, J. & Simon, N. (2011).
New insights and faster computations for the
graphical lasso. J COmput Graph Stat 20, 892–
900.
Wright, S. (1921). Corelation and causation. J
Agric Res , 558–585.
Yuan, M. & Lin, Y. (2007). Model selection
and estimation in the Gaussian graphical model.
Biometrika 94, 19–35.
Zhou, Q. (2011). Multi-domain sampling with applications to structural inference of Bayesian networks. J Am Stat Assoc 106, 1317–1330.

24

