Inferring large graphs using ℓ1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _1$$\end{document}-penalized likelihood

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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _1$$\end{document}-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.


Introduction
Revealing the true structure of a complex system is paramount in many fields to identify system regulators, predict its behav-  (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 (Barabási and 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 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 and Sood 2004;Tsamardinos et al. 2006;Grzegorczyk and 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 (Bühlmann and 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 and Wong 2008;Guyon et al. 2010;Fu and 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 and 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).
A large number of methods have been proposed for estimating DAGs, including for instance score-based methods (Bayesian score, 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 high-dimensional case, but requires a test of conditional independences that quickly becomes computationally intractable (Kalisch and Bühlmann 2007). Recent works stress the limitations of the absence of cycles in DAGs in the study of complex systems (De Smet and 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 Modeling 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 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).
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 and Werhli 2007), or a partial knowledge of the network can inform the inference process efficiently, e.g., in the semisupervised framework of Mordelet and Vert (2008). For a known order among the variables in the graph, Shojaie and Michailidis (2010) present results for the estimation of highdimensional graphs based on independent linear regressions using an adaptive Lasso scheme. When the order of the variables is unknown, van de Geer and Bühlmann (2013) studied the convergence of the 0 -penalized likelihood. However, the 0 -regularized approaches (Silander and Myllymäki 2006;Hauser and Bü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 and 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 unifying framework for pseudolikelihood-based graphical modeling 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 and Zhou (2013) and Shojaie and 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 twofold: 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 Sect. 4 our inference algorithm. Section 5 reports numerical experiments both on toy problems and realistic data sets.

DAG's modeling and estimation
This work considers the framework of an unknown DAG G 0 = (V, E), consisting of vertices V = {1, . . . , p} and a set of edges E ⊆ V × V . The p nodes are associated to random variables X 1 , . . . , X p . A natural approach, developed by Meinshausen and Bühlmann (2006) to solve the network 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 = i) through the Gaussian Structural Equation Model: with ε j ∼ N (0, σ 2 j ) (σ 2 j 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 nonzero coefficients of G 0 , i.e., (G 0 ) j i 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 Eq. (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: , j≤ p is the p × p matrix compatible with the graph G 0 and ε := (ε 1 , . . . , ε p ) is the n × p matrix of noise vectors.
The negative log-likelihood of the model is then (Rau et al. 2013): where I denotes the p × p identity matrix and X k the vector of length p corresponding to the k-th observation of X 1 , . . . , X p . To recover the structure of the DAG G 0 and make the estimated graph sparse enough, we focus on a penalized maximum likelihood procedure (Bickel and Li 2006): where (.) is the negative log-likelihood of Eq. (3), pen(.) is a penalization function, λ is a trade-off parameter between penalization and fit to the data, and G DAG 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. (2011Peters et al. ( , 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 Eq. (4) is an attractive option to infer sparse graphs. From a computational point of view, the main difficulty when solving the optimization problem in Eq. (4) lies in exploring the set of DAGs G DAG . (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 G DAG , a dynamic programming method has been proposed in Silander and Myllymä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 and Bü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 G 0 ). 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 Sect. 2.2). Given Eq. (3) and omitting constant terms, the 1 -penalized likelihood estimator we consider is: where for any matrix M :

A new formulation for the estimator
We propose here a new formulation of the minimization problem of Eq. (5). It naturally decouples the initial problem into two steps of the minimization 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 1 from Bühlmann (2013) then gives an equivalent condition for a matrix to be compatible with a DAG.

Proposition 1 (Bü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:
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 nonzero entries of G, as illustrated in Example 1.
Example 1 Consider the DAG G given in Fig. 1 (top). The corresponding matrix G can then be written as the strictly  Fig. 1 An example of DAG G (top) and the action of P and T on G: P is associated to a complete graph 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 ) lower triangular matrix T by permutation of its rows and columns using P: Looking at the nonzero 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 Fig. 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 ).
Using Proposition 1, the estimator in (5) leads to the following equivalent optimization problem: where the optimization space of constraints C is defined as C = P p (R) × T p (R), with P p (R) the set of permutation matrices and T p (R) the set of strictly lower triangular matrices. Note that a similar formulation has already been proposed by van de Geer and Bü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 (Sect. 3) and a computationally effficient algorithm to solve Problem (6) (Sect. 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 Eq. (6), both in prediction (Eq. 7) and estimation (Eq. 8). Following the works of van de Geer and Bü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.

Estimating the true order of variables
For a known ordering among the variables of the graph (Shojaie and Michailidis 2010), an unrealistic 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 G 0 : Π 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 variablesP given by Eq. (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 ∈ P p (R), we denote by G 0 (P) the matrix defined as: From a graphical point of view, while P / ∈ Π 0 , the graph G 0 (P) associated to G 0 (P) is obtained from G 0 by permuting some of its nodes (see Example 2), otherwise, if P ∈ Π 0 , G 0 (P) = G 0 . We also denote by ε(P) := X − X G 0 (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 Sect. 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 Sect. 3.3).

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

Assumptions on the model
For a square matrix M ∈ M p× p (R) and a subset S of 1, p 2 , we denote by M S ∈ M p× 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.
Hypotheses H 1 There exists σ 2 > 0 such that H 2 There exists σ 2 0 , independent of p and n, such that max 1≤ j≤ p Graph G 0 (top) and the permuted graph G 0 (P) (bottom) associated to the permutation P H 3 There exists λ * > 0 such that the minimal eigenvalue of the covariance matrix Σ of X satisfies H 4 There exists g max < ∞ such that the maximal weight of the DAG G 0 is bounded max 1≤i, j≤ p H 5 The number of nodes p satisfies where the minimum is taken over the set of p × p matrices satisfying M S C 1 ≤ 3 M S 1 , with S ⊂ 1, p 2 and |S| ≤ t.
Assumption H 1 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 results are still valid even if σ 2 = 1, by small modifications in the constant terms as long as they are all equal.
Assumption H 5 deserves a special attention since it bounds the high-dimensional setting. The considered problem is obviously non-trivial and requires a sufficient amount of information. A more detailled discussion about assumptions H 3 and H 5 is proposed in Sect. 3.4.
Assumption H 6 is a natural extension of the Restricted Eigenvalue condition of Bickel et al. (2009) to our multitask setting. More precisely, denoting H 6 is equivalent to assuming that the Gram matrixXX T n is non-degenerate on a restricted cone (Lounici et al. 2009;Bühlmann and 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 and Wolf 2015).
The last assumption H 7 is an identifiability condition needed to ensure that the estimated permutationP is in Π 0 . This assumption was introduced by van de Geer and Bühlmann (2013) as the "omega-min" condition. In a sense, it separates the set of compatible permutations from its complement in a finite sample scenario.

Main result
The result we establish in this section is double-edged: (a) with large probability, we ensure that the estimatedP belongs to Π 0 , and (b) we provide convergence inequalities both in prediction and estimation for the graph estimated from the minimization problem (6). This result clearly states the desirable theoretical properties of the derived estimator, assuming reasonable conditions on the complex system embedding the data.
Then, with probability greater than 1 − 5/ p, any solutionĜ =PTP T of the minimization problem (6) satisfies thatP ∈ Π 0 . Moreover, with at least the same probability, the following inequalities hold: The proof of this result is deferred in Section C of the Supplementary Materials.
Theorem 1 states that with probability at least 1 − 5/ p, we choose a compatible order of variables over the set of permutations. Inequalities (7) and (8) give non-asymptotic upper bounds on the loss under conditions depending on s, p and n (see Sect. 3.4). They also ensure that the estimated T is close to the true T 0 with large probability.

Discussion on the high-dimensional scenario
Sparsity of the graph Assumption H 7 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Ĝ choosing η ≤ α 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 H 5 . Taking at least p 2 log( p) = O(n) ensures a good estimation of the graph. Note, however, that universal conditions cannot be overcome and the ultra-high dimension settings (e.g., Wainwright (2009); Verzelen (2012)) is an insurmountable limit.
Minimal eigenvalue condition Assumption H 3 ensures that the minimal eigenvalue of the covariance matrix Σ of X is not too small. In the high-dimensional scenario, this could be hard to verify, λ min decreasing while n, p growing to infinity (Hogben 2007). A natural bound for λ min is: with g max and s as in H 4 and Theorem 1. Assumption H 3 can thus be relaxed by allowing λ min to decrease with 1/ p √ s. The price to pay for this relaxation is a data dimensionality reduction p 3 log( p) = O(n), which automatically implies: with Eq. (9) (for more details, see Section A of the Supplementary Material Proof details).

Inference algorithm 4.1 Global algorithm overview
In this section, we present GADAG (Genetic Algorithm for learning Directed Acyclic Graphs), a computational procedure devoted to solve Eq. (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 score function, then the roles of P and T are reversed and the procedure is repeated iteratively until convergence for some criterion (Csiszár and Tusná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 Fig. 3). As we show in the following, population-based meta-heuristics algorithms are a natural and efficient choice for exploring the space of permutation matrices (Sect. 4.3). The nested optimization problem can be resolved using a steepest descent approach (Sect. 4.2).

Graph structure learning when the variable order is fixed
Assume first that the variable ordering P ∈ P p (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: The minimization problem given by Eq. (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 and Bü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 and 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 Eq. (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 Fig. 3 and the "Evaluate the new individuals" item in Algorithm 2). Moreover, its mechanistic Algorithm 1: Graph structure learning -minimization of Eq. (10) Input: λ, L , > 0. Initialization: T 0 the null squared p × p matrix, k = 0 and e = +∞. while e > do with Eq. (13); Using Eq. (11), compute the current matrix the unique solution of (10).

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 of the LARS algorithm with matrix arguments.
Let (T k ) k≥0 the sequence of matrices defined for all i, j ∈ 1, p 2 as: where for all k ≥ 0, U k = T k − ∇ 1 n X (I −PT k P T ) 2 F L , L is the Lipschitz constant of the gradient function ∇ 1 n X (I − PT k P T ) 2 F and sign() is the sign of any element. Then, a solution of (10) is given by performing Algorithm 1, where: -the gradient of 1 n X The detailed calculations are deferred to Section B of the Supplementary Material Proof details.

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 P p (R) space of permutation matrices in dimension p and to evaluating the quality of permutation candidates P ∈ P p (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;Dré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 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 Sect. 5.1. The details of the different operators are given in the following.
As we show in Example 3, any P ∈ P p (R) is uniquely defined by a permutation vector of 1, p . Hence, we use as a the search space S p the set of permutations of 1, p , which is a well-suited formulation for GAs.
Generate P k+1 as a random selection of N individuals from P k ; Pick an even subset P xo of P k+1 (each individual of P k+1 selected with probability p xo ); Perform crossover on P xo by randomly pairing the individuals; Mutate each obtained individual with probability p m ; Evaluate the new individuals P m by running Algorithm 1; Replace P xo by P m in P k+1 ; Compute the Shannon entropy H and the difference in the average fitness e J = max 0≤i≤imax J (P k+1 ) −J (P k−i ) ; Increase k: k ← k + 1; end Example 3 Consider the permutation matrix ( p = 5): Then, P is represented by the 5 3 4 1 2 vector, looking at the ranks of non-null values of P column by column. The nodes are ranked according to their topological ordering.
Note that our problem closely resembles the classical Traveling Salesman Problem (TSP), which has 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. 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 p i ∈ S p , the fitness function is defined as: with P i constructed from p i as in Example 3 and T * i the solution of Eq. (10) with P = P i . 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 algorithm GADAG.
Selection operator The selection operator (or survival step) consists in generating a population of N individuals from the N existing individuals by random sampling (with replacement, hence some individuals are duplicated and others are deleted). It aims at improving the average quality of the population by giving to the best potential solutions a higher probability to be copied in the intermediate population. We have chosen to use the classical proportional selection of Holland (1975): each individuals is selected with a probability inversely proportional to its fitness value of Eq. (14).
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 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 p xo ; 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 p 1 and p 2 , a random set of crossover points are selected, which we denote Ω. It consists in a permutation of k elements taken from 1, p , with k uniformly drawn between 0 and p. A first child C 1 between p 1 and p 2 is then generated by: 1. fixing the crossover points of p 1 , 2. completing C 1 with the missing numbers in the order they appear in p 2 .
Example 4 Consider the two following parents: p 1 4 3 10 7 5 9 1 2 6 8 p 2 6 1 9 4 10 2 8 3 7 5 Assume that the crossover points randomly chosen are 4, 9, 2 and 8 (in bold red above). Then, the child C 1 is defined by inheriting those points from p 1 and filling the other points in the order they appear in p 2 : 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 C 1 4 * * * * 9 * 2 * 8 p 2 6 1 9 4 10 2 8 3 7 5 / / / / ⇓ 4 6 1 10 3 9 7 2 5 8 From a graphical point of view, a crossover between p 1 and p 2 , which encode two complete graphs G P 1 and G P 2 , constructs two new graphs. One of them, G C 1 is composed of the sub-graph of G P 1 induced by the set of crossover points Ω and the sub-graph of G P 1 induced by the complementary set Ω C of Ω in 1, p (see Fig. 4). The second child graph G C 2 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 neighboring 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 p m .
Example 5 A possible mutation for the first child of Example 4 is to swap the genes "1" and "10" (in bold red below): M 1 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 ∈ 1, p as: where N i, j is the number of times when i appears in position j. H j = 0 if all the individuals "agree" 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 distribution of the different nodes at a given position and the population is in this case at a maximum of heterogeneity or disorder for this position. The algorithm stops if the population entropy value H = N j=1 H j drops below a threshold since H = 0 if all the individuals are identical. A second criterion can terminate GADAG if difference in the average fitness (denotedJ thereafter) of the population between a given number of consecutive iterations, does not change by more than a predefined threshold.

Numerical experiments
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 analyze 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 Sect. 5.4.1. In Sect. 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 Sect. 5.3 for the simulated high-dimensional toy datasets and in Sect. 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 Sect. 5.3 correspond to a Windows 7 laptop computer with 8 threads on a 4-core hyperthreaded 2.50GHz processor, with 4GB of RAM.

Algorithm parameters
Running the procedure of Algorithm 2 requires to define parameters of the outer loop, which generates our population of P's, and of the nested loop to find the optimal T * . The evaluation of the Lipschitz gradient constant L, used to find the optimal graph structure T * , is known as a hard established problem in optimization. Some authors propose to choose an estimate of L from a set of possible values (Jones et al. 1993;Sergeyev and Kvasov 2006), to estimate local Lipschitz constants (Sergeyev 1995), or to set it a priori to a fixed value (Evtushenko et al. 2009;Horst and Pardalos 1995). Here, observing Eq. (13), a major bound for L is given by: We found that setting L to this bound worked well in practice in all our scenarios. Five parameters need to be tuned to run the Genetic Algorithm: the crossover rate p xo , the mutation rate p m , the constant of the stopping criteria H and J and the size of the population N . For the first four parameters, we observed that their value had a limited effect on the efficiency, hence we chose commonly used values in the literature (see Table 1). The size of the population has a more complex effect and has been investigated in several prospective papers (e.g., Schaffer et al. 1989;Alander 1992;Piszcz and Soule 2006;Ridge 2007) but without providing a definitive answer to the problem. In our simulation study, we chose as a rule-of thumb N = 5 p, which was found as a good compromise between computational cost and space exploration on several experiments.
The complete parameter settings used in our experiments are reported in Table 1.

Performance metrics
A classical performance measure for graph inference methods consists in comparing predicted interactions with the known edges in the true graph G 0 using precision versus recall (P/R) curves. We denote by TP, FP, FN and TN, the true positive (correctly predicted) edges, the false positive (inferred by mistake) edges, the false negative (missed) edges, and the true negative (correctly non-predicted) edges. The recall, defined as T P T P+F N , measures the power (or sensitivity) of reconstruction of nonzero elements of the true matrix G (or equivalently of the true network) for one method, whereas the precision, equal to T P T P+F P , measures the accuracy of the reconstruction. The closer to one the precision and the recall the better. P/R curves represent the evolution of those quantities when varying the sparsity of the methods. GADAG is based on penalized optimization: it seeks linear dependencies between the variables with a controlled level of parsimony (λ in Eq. (5)). For λ varying from 0 (complete graph) to +∞ (empty graph), it thus produces a list of edges successively introduced in the model. This list of edges defines the precision versus recall curve. As a summary performance measurement, we also computed the classical area under the P/R curve (AUPR).

Exploratory analysis on toy datasets
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 nonzero parameters of the matrix G 0 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 Sect. 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 1000 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 summarized in Table 2. We can first remark a crude decrease in 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 in dimensionality. For large graphs ( p = 500 and 1, 000), even if n/ p ≤ 10%, it succeeds in recovering them, which makes it a competitive algorithm with regard 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 300 s for hundred of nodes and 1500 s for thousand of nodes.

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 organismstranscriptional 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) versus 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 analyzing a real data 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 very challenging task.

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 Boot-Lasso (Allouche et al. 2013) was one of the key components of the best performing approach of one of the DREAM5 sub-challenges (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 and Wiener 2002) for our results. The Lasso is a 1 -norm penalization technique for solving linear regression. Following the works of Bach (2008), Boot-Lasso uses 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) modeling. 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 maximizing 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), although 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 independence 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 Sect. 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).

Numerical results
The P/R curves for the five DREAM problems are shown in Fig. 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 Fig. 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 ( Fig. 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 precision 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.

Conclusion and discussion
In this paper, we proposed a hybrid genetic/convex algorithm for inferring large graphs based on a particular decomposition of the 1 -penalized maximum likelihood criterion. We obtained two convergence inequalities that ensure that the  graph estimator converges to the true graph under assumptions that mainly control the model structure: graph size (balance between sparsity, number of nodes and maximal degree) and signal-to-noise ratio. From an algorithmic point of view, the estimation task is split into two sub-problems: node ordering estimation and graph structure learning. The first one is a non-trivial problem since we optimize over a discrete non-convex large dimensional set. It led us to use a heuristic approach we specifically tailored to achieve the optimization task. The second one is a more common problem, related to the Lasso one, for which we proposed a sound procedure with theoretical guarantees. The potential of such an approach clearly appeared in the numerical experiments, for which the behavior of our algorithm seemed to be very competitive when compared to the state-of-the-art. 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 be to consider other optimization schemes for the node ordering with more established convergence proofs (e.g., simulated annealing (Granville et al. 1994)).
Second, other potential extensions involve algorithmic considerations in order to improve the calculation time, including a finer calibration of the algorithm parameters, an initialization step for the gradient descent, and, in general, improving the interactions 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 the noise variances of all graph nodes to be equal to ensure graph identifiability (that is no equivalence class of graphs). Such a hypothesis is of course restrictive and likely to be violated for real datasets. In order to infer networks for any noise variances, one solution consists in incorporating interventional data on the model. These data are obtained from perturbations of the biological system (e.g., gene knockouts or over-expressions) and make the equivalence class of graphs smaller (Hauser and Bühlmann 2012). The use of additional data, informative yet very costly interventional data could be combined with observational on the MLE estimator. It was recently proposed by Hauser and Bühlmann (2015) for a BIC-score penalized MLE, or by Rau et al. (2013) for learning Gaussian Bayesian networks in the case of GRN inference. A modification of our hybrid algorithm GADAG could then lead to a more accurate identification of the true graph. Lastly, the cyclic structure framework could also be considered by using a Markov equivalence characterization (Richardson 1997) to relax the strictly triangular assumption on our matrix T using Equation of Proposition 1. It would pave the way for totally new theoretical developments, and a more realistic modeling of a complex system.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.