Julien Chiquet
email: [julien.chiquet@genopole.cnrs.fr

Yves Grandvalet
email: yves.grandvalet@utc.fr

Christophe Ambroise
email: christophe.ambroise]@genopole.cnrs.fr

Inferring Multiple Graphical Structures

Keywords: Network inference, Gaussian graphical model, Multiple sample setup, Coop-LASSO, Intertwined-LASSO

Gaussian Graphical Models provide a convenient framework for representing dependencies between variables. Recently, this tool has received a high interest for the discovery of biological networks. The literature focuses on the case where a single network is inferred from a set of measurements, but, as wetlab data is typically scarce, several assays, where the experimental conditions affect interactions, are usually merged to infer a single network. In this paper, we propose two approaches for estimating multiple related graphs, by rendering the closeness assumption into an empirical prior or group penalties. We provide quantitative results demonstrating the benefits of the proposed approaches. The methods presented in this paper are embeded in the R package simone from version 1.0-0 and later.

Motivations

Systems biology provides a large amount of data sets that aim to understand the complex relationships existing between the molecular entities that drive any biological process. Depending on the molecule of interest, various networks can be inferred, e.g., gene-to-gene regulation networks or protein-protein interaction networks. The basic idea is to consider that if two molecules interact, a statistical dependency between their expression should be observed.

A convenient model of multivariate dependence patterns is Gaussian Graphical Modeling (GGM). In this framework, a multidimensional Gaussian variable is characterized by the so-called concentration matrix, where conditional independence between pairs of variables is characterized by a zero entry. This matrix may be represented by an undirected graph, where each vertex represents a variable, and an edge connects two vertices if the corresponding pair of random variables are dependent, conditional on the remaining variables.

Merging different experimental conditions from wetlab data is a common practice in GGM-based inference methods [START_REF] Toh | Inference of a genetic network by a combined approach of cluster analysis and graphical gaussian modeling[END_REF]Horimoto 2002, Schäfer and[START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF]. This process enlarges the number of observations available for inferring interactions. However, GGMs assume that the observed data form an independent and identically distributed (i.i.d.) sample. In the aforementioned paradigm, assuming that the merged data is drawn from a single Gaussian component is obviously wrong, and is likely to have detrimental side effects in the estimation process.

In this paper, we propose to remedy this problem by estimating multiple GGMs, each of whom matching different modalities of the same set of variables, which correspond here to the different experimental conditions. As the distributions of these modes have strong commonalities, we propose to estimate these graphs jointly. Considering several tasks at a time has been attracting much attention in the machine learning literature, where the generic problem is usually referred to as "multi-task learning" [START_REF] Caruana | Multitask learning[END_REF]. [START_REF] Efron | The future of indirect evidence[END_REF] used the terms "indirect evidence" and "learning from the experience of others" for similar ideas. The principle is to learn an inductive bias, whose role is to stabilize the estimation process, hopefully enabling more accurate predictions in small sample size regimes [START_REF] Baxter | A model of inductive bias learning[END_REF]. The techniques comprise the empirical and hierarchical Bayes methodologies and regularization schemes (see Argyriou et al. 2008, for example), which may be interpreted as approximations to the latter. Here, we will mainly follow the penalty-based approach to leverage the inference of several related graphs towards a common pattern.

A typical example of this problem arises when inferring gene interactions from data measured on slightly different stem cells, such as wild and mutant. It is reasonable to assume that most, though not all, interactions will be common to both types of cells. The line of attack we propose alleviates the difficulties arising from the scarcity of data in each experimental condition by coupling the estimation problems. Our first proposal biases the estimation of the concentration matrices towards a common value. Our second proposition focuses on the similarities in the sparsity pattern that are more directly related to the graph itself. We propose the Cooperative-LASSO, which builds on the Group-LASSO [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] to favor solutions with a common sparsity pattern, but encodes a further preference towards solutions with similar sign patterns, thus preserving the type of co-regulation (positive or negative) across assays.

To our knowledge, the present work is the first to exploit the multi-task learning framework for learning GGMs. However, coupling the estimation of several networks has recently been investigated for Markov random fields, in the context of time-varying networks. [START_REF] Kolar | Estimating time-varying networks[END_REF] propose two specific constraints, one for smooth variations over time, the other one for abrupt changes. Their penalties are closer to the Fused-LASSO and total variations penalties than to the group penalties proposed here.

Network Inference with GGM

In the GGM framework, we aim to infer the graph of conditional dependencies among the p variables of a vector X from independent observations (X 1 , . . . , X n ). We assume that X is a p-dimensional Gaussian random variable X ∼ N (0 p , Σ). Let K = Σ -1 be the concentration matrix of the model; the non-zero entries of K ij indicate a conditional dependency between the variables X i and X j , and thus define the graph G of conditional dependencies of X.

The GGM approach produces the graph G from an inferred K. The latter cannot be obtained by maximum likelihood estimation that would typically return a full matrix, and hence a useless fully connected graph. To produce sparse networks, [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF] propose to penalize the entries of K by an 1norm. [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] latter addressed the very same problem with an elegant algorithm named the graphical-LASSO. Their well-motivated approach produces a sparse, symmetric and positive-definite estimate of the concentration matrix. However, a cruder though more direct approach has been reported to be more accurate in terms of edge detection [START_REF] Villers | Assessing the validity domains of graphical Gaussian models in order to infer relationships among components of complex biological systems[END_REF][START_REF] Rocha | A path following algorithm for sparse pseudolikelihood inverse covariance estimation[END_REF]. This approach, proposed by [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF] and known as neighborhood selection, determines G via an iterative estimation of the neighborhood of its nodes. For this purpose, it considers p independent 1 -penalized regression problems. Let X be the n × p matrix of stacked observations, whose kth row contains (X k ) . The vertices adjacent to vertex i are estimated by the non-zero elements of β solving min

β∈R p-1 1 n X i -X \i β 2 2 + λ β 1 , (1)
where X i is the ith column of X and X \i is X deprived of its ith column: the ith variable is "explained" by the remaining ones. As the neighborhood of the p variables are selected separately, a post-symmetrization must be applied to manage inconsistencies between edge selections; Meinshausen and Bühlmann suggest AND or OR rules, which are both asymptotically consistent (as n goes to infinity). Solving the p regression problems (1) may be interpreted as inferring the concentration matrix in a penalized, pseudo maximum likelihood framework, where the joint distribution of X is approximated by the product of the p distributions of each variable conditional on the other ones [START_REF] Rocha | A path following algorithm for sparse pseudolikelihood inverse covariance estimation[END_REF][START_REF] Ambroise | Inferring sparse Gaussian graphical models with latent structure[END_REF][START_REF] Ravikumar | High-dimensional Ising model selection using 1 -regularized logistic regression[END_REF], that is

L(K|X) = p i=1 n k=1 log P(X k i |X k \i ; K i ) ,
where X k \i is the kth realization of the vector X deprived of the ith coordinate. Considering the Gaussian assumption on the generation of the data X, the pseudo-log-likelihood admits a compact and simple expression (see derivation in Appendix A.1):

L(K|X) = n 2 log det(D) - n 2 Tr D -1 2 KSKD -1 2 - np 2 log(2π) , (2)
where S = n -1 X X is the empirical covariance matrix, and D is a p×p diagonal matrix with elements D ii = K ii . In the sequel, it will be convenient to use the sufficiency of S for K, and, by a slight abuse of notations, write L(K|X) = L(K|S).

Following [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF], an 1 penalty may be added to obtain a sparse estimate of K. Nevertheless, our approach to maximizing the pseudo-loglikelihood is much simpler than the optimization of the log-likelihood proposed by [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF]. Indeed, as stated formally in the following proposition, maximizing the penalized pseudo-log-likelihood on the set of arbitrary matrices (not constrained to be either symmetric or positive definite) boils down to solving p independent LASSO problems of size p-1. Furthermore, for the purpose of discovering the graph structure, additional computational savings are achieved by remarking that D needs not to be estimated. We thus avoid the iterative scheme of [START_REF] Rocha | A path following algorithm for sparse pseudolikelihood inverse covariance estimation[END_REF] alternating optimization with respect to D and to the off-diagonal elements of D -1 K.

Proposition 1. Consider the following reordering of the rows and columns of K and S:

K \i\i K i\i K i\i K ii , S \i\i S i\i S i\i S ii , (3)
where K \i\i is matrix K deprived of its ith column and its ith line, and where K i\i is the ith column of K deprived of its ith element. The problem

max {Kij :i =j} L(K|S) -λ K 1 , (4)
where K 1 is the componentwise 1 -norm, can be solved column-wisely by considering p LASSO problems in form min

β∈R p-1 1 2 S 1/2 \i\i β + S -1/2 \i\i S i\i 2 2 + λ n β 1 , (5)
where the optimal β is the maximizer of (4) with respect to K -1 ii K i\i as defined in (3). Hence, Problem (4) may be decomposed into the p Problems (5) of size p -1 generated by the p possible permutations in (3).

The full solution to Problem (4), in {K ij : i = j}, also requires K ii (see [START_REF] Rocha | A path following algorithm for sparse pseudolikelihood inverse covariance estimation[END_REF]. However, since our interest is to unveil the graph structure, the sparsity pattern of the penalized maximum likelihood estimate of K is sufficient, and the latter is directly recovered from β.

Proof. See appendix A.2.

From the definition of the covariance matrix S, it is clear that Problem (1) is a slight reparameterization of Problem (5). Hence, the graph produced by the approach of [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF] is identical to the one obtained by maximizing the penalized pseudo likelihood [START_REF] Ambroise | Inferring sparse Gaussian graphical models with latent structure[END_REF], Proposition 8).

Inferring Multiple GGMs

In transcriptomics, it is a common practice to conduct several assays where the experimental conditions differ, resulting in T samples measuring the expression of the same molecules. From a statistical viewpoint, we have T samples belonging to different sub-populations, hence with different distributions. Assuming that each sample was drawn independently from a Gaussian distribution X (t) ∼ N (0 p , Σ (t) ), the T samples may be processed separately by following the approach described in Section 2. The objective function is expressed compactly as a sum: max

{K (t) ij :i =j} T t=1 T t=1 L(K (t) |S (t) ) -λ K (t) 1 . ( 6
)
Note that it is sensible to apply the same penalty parameter λ for all samples provided that the T samples have similar sizes and originate from similar distributions, in particular regarding scaling and sparseness. Problem (6) ignores the relationships between regulation networks. When the tasks are known to have strong commonalities, the multi-task learning framework is well adapted, especially for small sample sizes, where sharing information may considerably improve estimation accuracy. To couple the estimation problems, we have to break the separability in K (1) , . . . , K (T ) in Problem (6). This may be achieved either modifying the data-fitting term or the penalizer. These two options result respectively in the graphical Intertwined-LASSO and the graphical Cooperative-LASSO presented below.

Intertwined Estimation

In the Maximum A Posteriori framework, the estimation of a concentration matrix can be biased towards a specific value, say S -1 0 . From a practical viewpoint, this is usually done by considering a conjugate prior on K , that is, a Wishart distribution W(S -1 0 , n). The MAP estimate is then computed as if we had observed additional observations of empirical covariance matrix S 0 .

Here, we would like to bias each estimation problem towards the same concentration matrix, whose value is unknown. An empirical Bayes solution would be to set S 0 = S, where S is the weighted average of the T empirical covariance matrices. As in the maximum likelihood framework, this approach would lead to a full concentration matrix. Hence, we will consider here a penalized criterion, which does not exactly fit the penalized maximum likelihood nor the MAP frameworks, but that will perform the desired coupling between the estimates of K (1) , . . . , K (T ) while pursuing the original sparseness goal.

Formally, let n 1 , . . . , n T be the sizes of the respective samples, whose empirical covariance matrices are denoted by S (1) , . . . , S (T ) . Also denote n = n t , we consider the following problem: max

{K (t) ij :i =j} T t=1 T t=1 L(K (t) | S(t) ) -λ K (t) 1 , (7)
where S(t) = αS (t) +(1-α) S and S = n -1 T t=1 n t S (t) . As this criterion amounts to consider that we observed a blend of the actual data for task t and data from the other tasks, we will refer to this approach as intertwined estimation. The idea is reminiscent of the compromise between linear discriminant analysis and its quadratic counterpart performed by the regularized discriminant analysis of [START_REF] Friedman | Regularized discriminant analysis[END_REF]. Although the tools are similar, the primary goals differ: [START_REF] Friedman | Regularized discriminant analysis[END_REF] aims at getting a control on the number of effective parameters, we want to bias empirical distributions towards a common model.

In order to avoid multiple hyper-parameter tuning, the results shown in the experimental section were obtained with α arbitrarily set to 1/2. More refined strategies are left for future work.

Graphical Cooperative-LASSO

The second approach consists in devising penalties that encourage similar sparsity patterns across tasks, such as the Group-LASSO [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], which has already inspired some multi-task learning strategies [START_REF] Argyriou | Convex multi-task feature learning[END_REF], but was never considered for learning graph models. We shortly describe how Group-LASSO may be used for inferring multiple graphs before introducing a slightly more complex penalty that was inspired by the application to biological interactions, but should be relevant in many other applications.

As in the single task case, sparsity of the concentration matrices is obtained via an 1 penalization of their entries. An additional constraint imposes the similarity between the two concentration matrices. Each interaction is considered as a group.

The Group-LASSO is a mixed norm that encourages sparse solutions with respect to groups, where groups form a pre-defined partition of variables. In the GGM framework, by grouping the partial correlations between variables across the T tasks, such a penalty will favor graphs G 1 , . . . , G T with common edges. The learning problem is then max

{K (t) ij :i =j} T t=1 T t=1 L(K (t) |S (t) ) -λ i =j T t=1 K (t) ij 2 1/2 . ( 8
)
Though this formalization expresses some of our expectations regarding the commonalities between tasks, it is not really satisfying here since we aim at inferring the support of the solution (that is, the set of non-zero entries of K (t) ).

To enable the inference of different networks (t, u), we must have some (i, j) such that

K (t) ij = 0 and K (u) ij = 0.
This event occurs with probability zero with the Group-LASSO, whose variables enter or leave the support group-wise [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. However, we may cure this problem by considering a regularization term that better suits our needs. Namely, when the graphs represent the regulation networks of the same set of molecules across experimental conditions, we expect a stronger similarity pattern than the one expressed in (8). Specifically, the co-regulation encompasses up-regulation and down-regulation and the type of regulation is not likely to be inverted across assays: in terms of partial correlations, sign swaps are very unlikely. This additional constraint is formalized in the following learning problem: max

{K (t) ij :i =j} T t=1 T t=1 L(K (t) |S (t) ) -λ i =j T t=1 K (t) ij 2 + 1/2 + T t=1 -K (t) ij 2 + 1/2 , (9)
where (u) + = max(0, u).

Figures 1 and 2 illustrate the role of each penalty on a problem with T = 2 tasks and p = 2 variables. They represent several views of the unit balls

2 i=1 2 t=1 β (t) i 2 1/2
≤ 1 , and

2 i=1 2 t=1 β (t) i 2 + 1/2 + 2 t=1 -β (t) i 2 + 1/2 ≤ 1 ,
that is, the admissible set for a penalty for a problem with two tasks and two features.

These plots also provide some insight on the sparsity pattern that originate from the penalty, since sparsity is related to the singularities at the boundary of the admissible set [START_REF] Nikolova | Local strong homogeneity of a regularized estimator[END_REF]. In particular, the first column illustrates that, when

β (2) 2 is null, β (1) 2
may also be exactly zero, while the second column shows that this event is improbable when β

(2) 2 differs from zero. The second row illustrates the same type of relationship between β

(2) 1 and β

(1) 1 that are expected due to the symmetries of the unit ball.

Figure 2 corresponds to a Cooperative-LASSO penalty. These plots should be compared with their Group-LASSO counterpart in Figure 1. We see that there are additional discontinuities in the unit ball resulting in new vertices on the 3-D plots. As before, we have that, when

β (2) 2 is null, β (1) 2
may also be exactly zero, but in addition, we may also have β

(1) 1 or β

(2) 1 exactly null. Accordingly, in the second and third row, we see that we may have

β (1) 2 null when β (2)
2 is non-zero. These new edges will result in some new zeroes when the Group-LASSO would have allowed a solution with opposite signs between tasks.

The second main striking difference with Group-LASSO is the loss of the axial symmetry of the Cooperative-LASSO when some variables are non-zero. These plots illustrate that the decoupling of the positive and negative parts of the regression coefficients in the penalty favors solutions where these coefficients are of same sign across tasks. The penalties are identical in the positive and negative orthant, but the Cooperative-LASSO penalization is more stringent elsewhere, when there are some sign mismatches between tasks. The most extreme situation occurs when there is no sign agreement across tasks for all variables. In the setup represented here, with only two tasks, the effective penalty then reduces to the LASSO.

β (2) 2 = 0 β (2) 2 = 0.1 β (2) 2 = 0.3 β (2) 1 = 0 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (2) 1 = 0.1 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (2) 1 = 0.3 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 Fig 1.
Admissible set for the Group-LASSO penalty for a problem with two tasks and two features. Top row: cuts of the unit ball through (β

(1) 1 , β (2) 1 , β (1)
2 ) for various values of β

(2) 2 , where (β

(1) 1 , β (2)
1 ) span the horizontal plane, and β

(1) 2 is on the vertical axis; bottom rows:

cuts through (β (1) 1 , β (1)
2 ) for various values of (β

(2) 1 and β (2) 2 ). β (2) 2 = 0 β (2) 2 = 0.1 β (2) 2 = 0.3 β (2) 1 = 0 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (2) 1 = 0.1 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (2) 1 = 0.3 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 β (1) 2 1 1 -1 -1 β (1) 1 Fig 2.
Admissible set for the Cooperative-LASSO penalty for a problem with two tasks and two features. Top row: cuts of the unit ball through (β

(1) 1 , β (2) 1 , β (1)
2 ) for various values of β

(2) 2 , where (β

(1) 1 , β (2)
1 ) span the horizontal plane, and β

(1) 2 is on the vertical axis; bottom rows:

cuts through (β (1) 1 , β (1)
2 ) for various values of (β

(2) 1 and β (2) 2 ).

Algorithms

In this section, we describe the strategy proposed for solving the three optimization problems introduced above, based upon the proposal of Osborne et al. (2000a) for solving the LASSO. This part also draws its inspiration from Osborne et al. (2000b), [START_REF] Kim | Blockwise sparse regression[END_REF], [START_REF] Roth | The group-lasso for generalized linear models: uniqueness of solutions and efficent algorithms[END_REF].

Problem Decomposition

The multiple independent tasks Problem (6) can be solved by considering either T single tasks like (4) of size (p -1) × p (each one possibly decomposed in p LASSO sub-problems of size p-1), or a single large problem of size T ×(p-1)×p, which can be decomposed into p LASSO sub-problems of size (p-1)×T , through Proposition 1. This line of attack is not computationally efficient here, but it will become advantageous when considering the penalties presented in Section 3.2. It is introduced at this point to provide a unified conceptual view of all algorithms.

Consider the (p T )×(p T ) block-diagonal matrix C composed by the empirical covariance matrices of each tasks

C =    S (1) 0 . . . 0 S (T )    ,
and define

C \i\i =     S (1) \i\i 0 . . . 0 S (T ) \i\i     , C i\i =     S (1) i\i . . . S (T ) i\i     .
(10)

The (p-1) T ×(p-1) T matrix C \i\i is the matrix C where we removed each line and each column pertaining to variable i. We define C, C\i\i and Ci\i similarly, with S (t) being replaced by S(t) for each t = 1, . . . , T in the above definitions.

Let β (t) ∈ R (p-1) denote the vector estimating

K (t)
i\i , defined from K (t) as in (3), and let β ∈ R T ×(p-1) be the vector of the concatenated estimates

β = (β (1) , • • • , β (T )
). The optimization of ( 6) is achieved by solving p subproblems in form:

min β∈R T ×(p-1) 1 2 C 1/2 \i\i β + C -1/2 \i\i C i\i 2 2 + λ T t=1 1 n t β 1 . (11)
Note that we do not need to perform the costly matrix operations that are expressed in the the first term of the objective function of Problem (11). In practice, we compute

f (β; C) = 1 2 β C \i\i β + β C i\i ,
which only differs from the squared 2 norm in (11) by a constant that is irrelevant for the optimization process.

Accordingly, Problems ( 7), ( 8) and ( 9) can be decomposed into p minimization sub-problems whose objective functions may be decomposed as

L k (β) = f (β) + λg k (β) , (12)
where, with a slight abuse of notation, f (β) is either f (β; C) for Problem (7) or f (β; C) for Problems ( 8) and ( 9), and where g k (β) stands for the penalty functions respectively defined below:

• for the graphical Intertwined LASSO

g 1 (β) = T t=1 1 n t β (t) 1 ;
• for the graphical Group-LASSO

g 2 (β) = p-1 i=1 β [1:T ] i 2
,

where β

[1:T ] i = β (1) i , . . . , β (T ) i
∈ R T is the vector of the ith component across tasks;

• for the graphical Cooperative-LASSO

g 3 (β) = p-1 i=1 β [1:T ] i + 2 + -β [1:T ] i + 2
.

Since f is convex with respect to β, and all penalties are norms, all these objective functions are convex and thus easily amenable to optimization. They are non-differentiable at zero, due to the penalty terms, which all favor zero coefficients. Bearing in mind the typical problems in biological data, where graphs have a few tens or hundreds nodes, and where connectivity is very weak1 , we need convex optimization tools that are efficient for medium-size problems with extremely sparse solutions. We thus chose a greedy strategy that aims at solving a series of small-size sub-problems, and will offer a simple monitoring of convergence.

Solving the Sub-Problems

The minimizers β of the objective functions ( 12) are assumed to have many zero coefficients. The approach developed for the LASSO by Osborne et al. (2000a) takes advantage of this sparsity by solving a series of small linear systems, whose size is incrementally increased/decreased, similarly to a column generation algorithm. The master problem is the original problem, but solved only with respect to the subset of variables currently identified as non-zero β coefficients. The subproblem of identifying new non-zero variables simply consists in detecting the violations of the first-order optimality conditions with respect to all variables. When there are no more such violations, the current solution is optimal.

The objective functions L k (β) are convex and smooth except at some locations with zero coefficients. Thus, the minimizer is such that the null vector 0 ∈ R p-1 is an element of the subdifferential ∂ β L k (β). In our problems, the subdifferential is given by

∂ β L k (β) = ∇ β f (β) + λ∂ β g k (β) , (13)
where ∇ β f (β) = C \i\i β + C i\i and where the form of ∂ β g k (β) differs for the three problems and will be detailed below.

The algorithm is started from a sparse initial guess, that is, β = 0 or, if available, the solution obtained on a more constrained problem with a larger penalization parameter λ. Then, one converges to the global solution iteratively, by managing the index A of the non-zero coefficients of β and solving the master problem over A, where the problem is continuously differentiable. The management of A requires two steps: the first one removes from A the coefficients that have been zeroed when solving the previous master problem, ensuring its differentiability at the next iteration, and the second one examines the candidate non-zero coefficients that could enter A. In this process, summarized in Algorithm 1, the size of the bigger master problems is typically of the order of magnitude of the number of non-zero entries in the solution. Solving the master problem with respect to the non-zero coefficients β A can be formalized as

solving min h L k (β A + h), where h ∈ R |A| is optimal if 0 ∈ ∂ h L k (β A + h).
Algorithm 1: General optimization algorithm

// 0. INITIALIZATION β ← 0 A ← ∅ while 0 / ∈ ∂ β L(β) do
// 1. MASTER PROBLEM: OPTIMIZATION WITH RESPECT TO β A Find a (approximate) solution h to the smooth problem

∇ h f (β A + h) + λ∂ h g k (β A + h) = 0 . // where ∂ h g k = {∇ h g k } β A ← β A + h
// 2. IDENTIFY NEWLY ZEROED VARIABLES while ∃i ∈ A : β i = 0 and min

θ∈∂ β i g k ∂f (β) ∂β i + λθ = 0 do A ← A\{i} // 3. IDENTIFY NEW NON-ZERO VARIABLES // Select i ∈ A c such that an infinitesimal change of β i provides the highest reduction of L k i ← arg max j∈A c v j , where v j = min θ∈∂ β j g k ∂f (β) ∂β j + λθ if v i = 0 then A ← A ∪ {i} else
Stop and return β, which is optimal

Implementation Details

We provide below the implementation details that are specific to each optimization problem. Specificity of each problem relies on ∂ β g k (β), denoted θ herein.

Intertwined LASSO -This LASSO problem is solved as proposed by Osborne et al. (2000a), except that we consider here the Lagrangian formulation with λ fixed.

The components of θ in the subdifferential (13) read

if β i = 0 then θ i ∈ [-1, 1] , else θ i = sign(β i ) .
Solving the master problem on A requires an estimate of θ A at β A + h . It is computed based on a local approximation, where the components of sign(β A +h) are replaced by sign(β A ). 2 This leads to the following descent direction h:

h = -β A -C-1 \i\i (A, A)( Ci\i (A) + λθ A ) ,
where, in order to avoid double subscripts, we use the notation M(A, A) for the square submatrix of M formed by the rows and columns indexed by A, and v(A) for the subvector formed by the columns of v indexed by A.

Then, before updating β A , one checks whether the local approximation used to compute h is consistent with the sign of the new solution. If not the case, one looks for the largest step size ρ in direction h such that β + A = β A + ρh is sign-consistent with β A . This amounts to zero a coefficient, say β i , and i is removed from A if |∂f (β + )/∂β i | < λ , otherwise, the corresponding θ i is set to -sign(∂f (β + )/∂β i ) . In any case, a new direction h is computed as above, and β A is updated until the optimality conditions are reached within A.

Finally, the global optimum is attained if the first-order optimality conditions are met for all the components of β, that is, if β verifies

0 ∈ C \i\i β + C i\i + λθ ,
where θ is such

θ A = sign( β A ) and θ A c ∞ ≤ 1 .
Graphical Group-LASSO -In this problem, the subdifferential ( 13) is conditioned on the norm of β

[1:T ] i
, the vector of the ith component across tasks.

Let θ

[1:T ] i = θ (1) i , . . . , θ (T ) i ∈ R T be defined similarly to β [1:T ] i , we have that, if β [1:T ] i 2 = 0 then θ [1:T ] i 2 ≤ 1 , else θ [1:T ] i = β [1:T ] i -1 2 β [1:T ] i ,
where, here and in what follows, 0/0 is defined by continuation as 0/0 = 0. As the subgradient w.r.t. β is non-zero, the management of the null variables is done here by subsets 2 When A is updated and that β i = 0, the corresponding θ i is set to -sign(∂f (β)/∂β i ) .

of T variables, according to

∇ β [1:T ] i f (β) 2
, instead of the one by one basis of the LASSO. Hence, we only need to index the groups i ∈ {1, . . . , p -1} in A.

Here also, solving the master problem on A requires an estimate of θ

[1:T ] A at β [1:T ] A +h . Provided that β [1:T ] i 2 = 0 for all i ∈ A, θ [1:T ] A is differentiable w.r.t. β [1:T ] A
. It will thus be approximated by a first-order Taylor expansion, resulting in a Newton-Raphson or quasi-Newton step. Here, we used quasi-Newton with BFGS updates. Note that, whenever β [1:T ] i = 0, that is, when a new group of variables has just been activated or is about to be deactivated, the corresponding θ

[1:T ] i is set so that ∇ β [1:T ] i f (β) + λθ [1:T ] i 2 (14)
is minimum (that is, with θ

[1:T ] i proportional to ∇ β [1:T ] i f
). The updates of A are also based on the minimal value of ( 14).

Graphical Cooperative-LASSO -As the Group-LASSO, the Cooperative-LASSO considers a group structure, but its implementation differs considerably from the former in the management of A. Though several variables are usually activated or deactivated at the same time, they typically correspond to subsets of β

[1:T ] i
, and these subsets are context-dependent; they are not defined beforehand. As a result, the index of non-zero β

A + = i ∈ {1, . . . , p -1} : β [1:T ] i + 2 > 0 ,
and A -= i ∈ {1, . . . , p -1} : -β

[1:T ] i + 2 > 0 .
Let T denote the index of non-zero entries of β

[1:T ] i
, with complement T c ; the subdifferential at the current solution is such that:

• if i ∈ A c + ∩ A c -, then max θ [1:T ] i + 2
, -θ

[1:T ] i + 2
≤ 1 ;

• if i ∈ A c + ∩ A -then θ T i = -β T i + -1 2 β T i , θ T c i : θ T c i + 2 ≤ 1 and -θ T c i + 2
= 0 ;

• i ∈ A + ∩ A c -, then θ T i = β T i + -1 2 β T i , θ T c i : -θ T c i + 2 ≤ 1 and θ T c i + 2 = 0 ; • if i ∈ A + ∩ A -, then θ (t) i = sign β (t) i β [1:T ] i + -1 2 β (t) i
, t = 1, . . . , T .

Once A + and A -are determined, the master problem is solved as for the Group-LASSO, with BFGS updates, with box constraints to ensure sign feasible β

[1:T ] i for i such that i ∈ A c + ∩ A -or i ∈ A + ∩ A c -.
When a new variable has just been activated or is about to be deactivated, the corresponding θ

[1:T ] i is set so that ∇ β [1:T ] i f (β) + λθ [1:T ] i + 2 + -∇ β [1:T ] i f (β) -λθ [1:T ] i + 2 (15)
is minimum. The updates of A are also based on the minimal value of (15).

Experiments

In most real-life applications, the major part of the inferred graphs are unknown, with little available information on the presence/absence of edges. We essentially face an unsupervised learning problem, where there is no objective criterion allowing to compare different solutions. As a result, setting the hyper-parameters is particularly troublesome, alike, say, choosing the number of components in a mixture model, and it is a common practice to visualize several networks corresponding to a series of penalties.

Regarding the first issue, we chose to present here synthetic and well-known real data that allow for an objective quantitative evaluation. Regarding the second issue, the problem of choosing penalty parameters can be guided by theoretical results that provide a bound on the rate of false edge discovery [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF][START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF][START_REF] Ambroise | Inferring sparse Gaussian graphical models with latent structure[END_REF], or by more traditional information criteria targeting the estimation of K (Yuan andLin 2007, Rocha et al. 2008). However, these proposals tend to behave poorly, and it is a usual practice to compare the performance of learning algorithms by providing a series of results, such as precision-recall plots or ROC-curves, letting the choice of penalty parameters as a mostly open question for future research. Although the shortcomings of this type of comparison are well-known [START_REF] Drummond | Cost curves: An improved method for visualizing classifier performance[END_REF]Holte 2006, Bengio et al. 2005), we will use precision vs. recall plots, which can be considered as valuable exploratory tools.

Precision is the ratio of the number of true selected edges to the total number of selected edges in the inferred graphs. Recall is the ratio of the number of true selected edges in the inferred graphs to the total number of edges in the true graphs. In a statistical framework, the recall is equivalent to the power and the precision is equivalent to one minus the false discovery proportion.

Synthetic Data

Simulation Protocol

To generate T samples stemming from a similar graph, we first draw an "ancestor" graph with p nodes and k edges according to the Erdős-Rényi model. Here, we consider a simple setting with T = 4 and a network with p = 20 nodes and k = 20 edges, as illustrated in Figure 3. Then, T children graphs are produced by random addition and deletion of δ edges in the ancestor graph. The T concentration matrices are built from the normalized graph Laplacians, whose offdiagonal elements are slightly deflated to produce strictly diagonally dominant matrices. To allow for positively and negatively correlated variables, we generate a strictly triangular matrix of random signs drawn from a Rademacher distribution. This matrix is symmetrized, complemented with ones on the diagonal, and its component-wise multiplication with the deflated Laplacians produces the ground-truth for the concentration matrices K (1) , . . . , K (T ) . Each K (t) is finally used to generate n Gaussian vectors with zero mean and covariance K (t) -1 .

Experimental Setup

The precision-recall plots are computed by considering the cumulative number of true and false edge detections among the T = 4 children networks. Let E (t) be the set of edges for children network t, precision and recall are respectively formaly defined as:

T t=1 (i,j)∈E (t) 1( K (t) ij ) T t=1 i>j 1( K (t) ij ) and T t=1 (i,j)∈E (t) 1( K (t) ij ) T t=1 E (t)
, where K (t) ij is the estimated partial correlation between variables i and j for network t, 1(u) returns 1 if u = 0 and 0 otherwise, and |E| is the cardinal of E.

To ensure representativeness, the precision-recall figures are averaged over 100 random draws of the ancestor graph, the averaging being performed for fixed values of the penalization parameter λ. That is, each point in the (precision,recall) plane is the average of the 100 points obtained for each random draw of the ancestor graph, for a given estimation method and for a given value of the penalization parameter. We compare our proposals, namely the Graphical Intertwined (α = 1/2), Cooperative and Group LASSO to two baselines: the original neighborhood selection of [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF], either applied separately to each graph (annotated "independent"), or computed on the data set merging the data originating from all graphs (annotated "pooled").

Results

Figures 4, 5 and 6 display precision-recall plots for nine prototypical situations. From Figure 4 to 6, the sample size increases, and from top to bottom, the differentiation between networks increases. First, note that the independent strategy is not influenced by the level of perturbation, yet only by the sub-sample size, as expected.

The top graph in Figure 4 represents the small-sample low-perturbation situation, where merging data sets is a good strategy, leveraging the independent analysis. The latter performs poorly, and our multi-task approaches dominate the pooled strategy, the Cooperative-LASSO being superior to the Group-LASSO, which has the advantage on the Intertwined LASSO. The medium/large-sized-sample low-perturbation (top graph in Figures 5 and6), small/medium-sized-sample medium-perturbation (middle graph in Figures 4 and5) and small-sized-sample large-perturbation (bottom graph in Figure 4) are qualitatively similar, with less differences between all competitors. For the large-sample medium-perturbation (middle graph in Figure 6) and mediumsized-sample low-perturbation (bottom graph in Figure 5) cases, all methods perform similarly. There is a slight advantage for the multi-task strategies for high recalls, and a slight advantage for the independent analysis for low recalls (that is, for high penalization, where there is less effective degrees of freedom to determine). The bottom graph in Figure 6 represents the large-sample highperturbation situation, where merging data is a bad strategy, since the networks differ significantly and there is enough data to estimate each network independently. The independent strategy works best, closely followed by the Intertwined LASSO. The Cooperative and Group-LASSO behave equally well for high recalls (low penalization parameters), but for highly penalized solutions (low recalls), they eventually become slightly worse than the pooled estimation.

These experiments show that our proposals are valuable, especially in the most common situation where data are scarce. Among the baselines, the usual pooled sample strategy is good in the small-sample low-perturbation, and the opposite independent strategy is better in the large-sample high-perturbation case. The intertwined LASSO is very robust, in the sense that it always performs favorably compared to the best baseline method over the whole spectrum of situations. Furthermore, except for the large-sample high-perturbation case, the Group-LASSO performs even better, and the Cooperative-LASSO improves further the supremacy of the multiple graph inference approach.

Protein Signaling Network

Only a few real data sets come with a reliable and exhaustive ground-truth allowing quantitative assessments. We make use of a multivariate flow cytometry 1.0 q q q q q q q q q q q q q q q q q q q q CoopLasso GroupLasso Intertwined Independent Pooled 1.0 q q q q q q q q q q q q q q q q q q q CoopLasso GroupLasso Intertwined Independent Pooled 1.0 q q q q q q q q q q q q q q q q q q q q CoopLasso GroupLasso Intertwined Independent Pooled 1.0 q q q q q q q q q q q q q q q q q q q CoopLasso GroupLasso Intertwined Independent Pooled q q q q q q q q q q q q q q q qq q q CoopLasso GroupLasso Intertwined Independent Pooled q q q q q q q q q q q q q q q q q q q CoopLasso GroupLasso Intertwined Independent Pooled 1.0 qq q q q q q q q q q q q q q q q q CoopLasso GroupLasso Intertwined Independent Pooled 1.0 qq q q q q q q q q q q q q q q q q CoopLasso GroupLasso Intertwined Independent Pooled 1.0 qq q q q q q q q q q q q q q q q q CoopLasso GroupLasso Intertwined Independent Pooled data set pertaining to a well-studied human T-cell signaling pathway [START_REF] Sachs | Causal protein-signaling networks derived from multiparameter single-cell data[END_REF]. The latter involves 11 signaling molecules (phospholipids and phosphorylated proteins) and 20 interactions described in the literature. The signaling network is perturbed by activating or inhibiting the production of a given molecule. Fourteen assays have been conducted, aiming to reveal different part of the network. Here, we used only four assays (inhibition of PKC, activation of PKC, inhibition of AKT, activation of PKA).

Recall

Recall

Recall

Recall

Recall

Recall

Recall

Recall

Recall

Graphs inferred using only one assay at a time show that each assay really focus on different part of the network (see Figure 7). When considering a strategy based on inference from multiple assays, the first false positive inferred by the Intertwined Graphical LASSO occurs when 11 true interactions out of 20 are detected (see Figure 8). This edge, between p38 and Jnk, is in fact due to an indirect connection via unmeasured MAP kinase kinases [START_REF] Sachs | Causal protein-signaling networks derived from multiparameter single-cell data[END_REF], which is a typical problem of confounding arising in this context. Considering partial correlations within the subset of available variables, the edge is correctly detected, but it is a false positive with respect to the biological ground truth. Furthermore, in larger biological networks, the absence of edge in the ground truth pathway often merely means that there is yet no evidence that the co-regulation exists. As a result, most real data evaluation of graph inference methods are based on qualitative subjective assessments by experts.

This caveat being, the various inference algorithms behave here as in the synthetic experiments: all inference methods perform about equally well for large samples (each assay consists here of about 1000 repeated measurements).

Figure 9 displays the results obtained for small sample sizes. Here also, the precision-recall plots are averaged over 100 independent random draws of samples of size n t , that is n = 4n t observations over the four considered assays. It is worth noticing that the large sample size limit is almost obtained for n t = 20. As for the synthetic experiments, the averaging is performed for fixed values of the penalization parameter λ. In this situation, our proposals dominate the best baseline strategy, which is pooled estimation. Again, the Intertwined LASSO is very robust, but the Group-LASSO, and to a greater extent the Cooperative-LASSO perform better in the small-sample-size regime.

Conclusion

This paper presents the first methods dedicated to the inference of multiple graphs in the Gaussian Graphical Model framework. In this setup, the two baseline approaches consist in either handling the inference problems separately or as a single one by merging the available data sets. Our proposals, motivated by bioinformatics applications, were devised to describe the dependencies between pairs of variables in analogous operating conditions, such as measurements recorded in different assays. This situation occurs routinely with omics data.

Our approaches are based on the neighborhood selection of [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF]. The first one, the Intertwined Graphical LASSO, relaxes the uniqueness constraint that is implicit when the tasks are processed as a single one, merely biasing the results towards a common answer. Our second approach, the Graphical Cooperative-LASSO, is based on a group-penalty that favors similar graphs, with homogeneous dependencies between the same pairs of variables. Homogeneity is quantified here by the magnitude and sign of partial correlations. The Cooperative-LASSO contrasts the Group-LASSO in being able to infer differing graph structures across tasks. Our experimental results show that our proposals are valuable and robust, consistently performing at least as well as the best of the two baseline solutions.

The algorithms developed in this paper are made available within the Rpackage simone from version 1.0-0 and later. This package also embeds extension of the multi-task framework to time-course data, that is, when transcriptomics data are collected by considering the same individual across time. This implementation builds on the 1 -penalized VAR(1) model described in [START_REF] Charbonnier | Weighted-lasso for structured network inference from time course data[END_REF]. 1.0 q q q q q q q q q q q q q q q q q Intertwined CoopLasso GroupLasso Independent Pooled 1.0 q q q q q q q q q q q q q q q q q q q q Intertwined CoopLasso GroupLasso Independent Pooled 1.0 q q q q q q q q q q q q q q q q q q q q Intertwined CoopLasso GroupLasso Independent Pooled As future work, we will provide a theoretical analysis of the Cooperative-LASSO regarding uniqueness of the solution and selection consistency, or sparsistence, that corresponds here to the asymptotic convergence of the set of detected edges towards the set of true edges.

Recall

Recall

Recall

Using the above identities in (17), we obtain

σ i = (1 -Σ i\i K i\i )/K ii + Σ i\i K i\i /K ii = 1/K ii , µ i = -K i\i X \i /K ii .
where µ i = (µ 1 i , . . . , µ n i ) . Using these notations and the corresponding blockwise notations for S (S ii = n -1 X i X i , S i\i = n -1 X \i X i and S \i\i = n -1 X \i X \i ), Equation ( 16) reads

L(K|X) = - n 2 p i=1 log σ i - p i=1 1 2σ i (X i -µ i ) (X i -µ i ) - np 2 log(2π) = n 2 p i=1 log K ii - np 2 log(2π) - n 2 p i=1 K ii S ii + 2 K ii S i\i K i\i + 1 K 2 ii K i\i S \i\i K i\i (18) = n 2 log det D - np 2 log(2π) - n 2 p i=1 1 K ii (K i SK i ) ,
where K i is the ith column of K and D is the diagonal matrix such that D ii = K ii . Finally, we use that p i=1

1 Kii (K i SK i ) = Tr(D -1 2 KSKD -1 2 ) to conclude the proof.

A.2. Blockwise Optimization of the pseudo-log-likelihood

Proof of Proposition 1. From (18), we have

L(K|S) = - n 2 p i=1 2S i\i K i\i + 1 K ii K i\i S \i\i K i\i + c, (19)
where c does not depend on K ij with j = i. Thus, if we discard the symmetry constraint on K, maximizing the pseudo-likelihood (19) with respect to the nondiagonal entries of K amounts to solve p independent maximization problems with respect to K i\i , i = 1, . . . , p. The summands of ( 19) can be rewritten as

- n 2K ii 2K ii S i\i K i\i + K i\i S \i\i K i\i = - nK ii 2 K -1 ii S 1/2 \i\i K i\i + S -1/2 \i\i S i\i 2 2 + c ,
where c = n/2 K ii S i\i S \i\i S i\i does not depend on K ij with j = i. Adding an 1 penalty term on K i\i and defining β = K -1 ii K i\i leads to the objective function of Problem (5), which concludes the proof.

A.3. Subdifferential for the Cooperative-LASSO

By definition, for a convex function g, the subdifferential is

∂g| β 0 = θ : ∀β , g(β) -g(β 0 ) ≥ θ (β -β 0 )
The function g(β) = (β) + 2 + (-β) + 2 has kinks whenever β has at least one zero component and that it has either no positive or no negative component. There are thus three situations where the subdifferential does not reduce to the gradient :

1. (β 0 ) + 2 = 0 and (-β 0 ) + 2 = 0, 2. (β 0 ) + 2 = 0 and (-β 0 ) + 2 = 0, 3. (β 0 ) + 2 = 0 and (-β 0 ) + 2 = 0, i.e. β 0 = 0.

For the first situation, denoting A the index of non-zero entries of β 0 and A c its complement, the subdifferential is defined as

(-β 0 ) + -1 2 β 0 + θ : θ A = 0 and ∀β A c , (β A c ) + 2 ≥ θ A c β A c . (20)
The set of admissible θ is explicitly given by θ : θ A = 0 , (θ A c ) + 2 ≤ 1 and (-θ A c ) + 2 = 0 .

(21)

Proof. We first show that, for any θ in the set defined in ( 21), the inequality in definition (20) always holds. Dropping the subscript A c for readability, we have: To finish the proof, it is sufficient to exhibit some β such that the inequality in definition (20) does not hold when (θ) + 2 > 1 or when (-θ) + 2 > 0.

For (θ) + 2 > 1, we choose β = (θ) + , yielding θ β = (θ) + 2 2 , and (β) + 2 = (θ) + 2 < (θ) + 2 2 , hence (β) + 2 < θ β; for (-θ) + 2 > 0, we choose β = -(-θ) + , yielding θ β = (-θ) + 2 2 > 0, and (β) + 2 = 0, hence (β) + 2 < θ β.

The second situation is treated as the first one, yielding

∂g| β 0 = (β 0 ) + -1 2 β 0 + θ : θ A = 0 , (-θ A c ) + 2 ≤ 1
and (θ A c ) + 2 = 0 .

For the last situation, the subdifferential, defined as

∂g| β 0 = θ : ∀β , (β) + 2 + (-β) + 2 ≥ θ β , (22)
reads ∂g| β 0 = θ : max (θ) + 2 , (-θ) + 2 ≤ 1 ,

Proof. We first show that, for all the elements of ∂g as explicitly defined in ( 23), the inequality in definition ( 22 ≤ max (θ) + 2 , (-θ) + 2 (β) + 2 + (-β) + 2 .

To finish the proof, it is sufficient to exhibit some β such that the inequality in definition ( 22

Fig 3 .

 3 Fig 3. Set of simulated graphs: ancestor (top) and two children (bottom) engendered by a δ = 2 perturbation.

Fig 4 .

 4 Fig 4. Precision-recall curves for the Intertwined, Cooperative, Group and the two baseline LASSO, for inferring four graphs (each with p = 20 nodes, k = 20 edges and a perturbation δ from the ancestor graph) from four samples of size nt = 25.

Fig 5 .

 5 Fig 5. Precision-recall curves for the Intertwined, Cooperative, Group and the two baseline LASSO, for inferring four graphs (each with p = 20 nodes, k = 20 edges and a perturbation δ from the ancestor graph) from four samples of size nt = 50.

Fig 6 .

 6 Fig 6. Precision-recall curves for the Intertwined, Cooperative, Group and the two baseline LASSO, for inferring four graphs (each with p = 20 nodes, k = 20 edges and a perturbation δ from the ancestor graph) from four samples of size nt = 100.

Fig 7 .Fig 8 .

 78 Fig 7. Four graphs inferred from single assay. From left to right, top to bottom, we have respectively graphs inferred from an assay: inhibiting akt, activating pka, inhibiting pkc, activating pkc. Thick black lines represent true positive and thin red lines are false positive.

Fig 9 .

 9 Fig 9. Precision-recall curves for the Intertwined, Cooperative, Group and the two baseline LASSO, for inferring the graphs on four assays of Sachs' data from four samples of size nt = 7, 10 and 20.

  θ β = (θ) + β -(-θ) + β = (θ) + β = (θ) + (β) + -(θ) + (-β) + ≤ (θ) + (β) + ≤ (θ) + 2 (β) + 2 ≤ (β) + 2 .

  ) always holds:θ β = (θ) + β -(-θ) + β ≤ (θ) + (β) + + (-θ) + (-β) + ≤ (θ) + 2 (β) + 2 + (-θ) + 2 (-β) + 2

  ) does not hold for max (θ) + 2 , (-θ) + 2 > 1. Without loss of generality, we assume (θ) + 2 > 1, and choose β = (θ) + , yieldingθ β = (θ) + 2 2 , and (β) + 2 + (-β) + 2 = (β) + 2 = (θ) + 2 < (θ) + 2 2 , hence (β) + 2 + (-β) + 2 < θ β.

Typically, the expected number of vertices in graphs to scale as the number of nodes, that is, we expect order of √ pT non-zero coefficients in each sub-problem of size T × (p -1).

Acknowledgments

Yves Grandvalet was partially supported by the PASCAL2 Network of Excellence, the European ICT FP7 under grant No 247022 -MASH, and the French National Research Agency (ANR) under grant ClasSel ANR-08-EMER-002. Christophe Ambroise and Julien Chiquet were partially supported by the ANR under grants GD2GS ANR-05-MMSA-0013 and NeMo ANR-08-BLAN-0304-01.

Appendix A: Proofs

A.1. Derivation of the pseudo-log-likelihood

We show here that the pseudo-log-likelihood

associated to a sample of size n drawn independently from the multivariate Gaussian vector X ∼ N (0 p , Σ) reads

where S = n -1 X X is the empirical variance-covariance matrix and D is the diagonal matrix such that D ii = K ii , for i = 1, . . . , p.

Proof. Since the joint distribution of X k is Gaussian, the distributions of X k i conditioned on the remaining variables X k \i are also Gaussian. Their parameters (µ k i , σ i ) are given by

where Σ \i\i is matrix Σ deprived of its ith column and its ith line, Σ i\i is the ith column of matrix Σ deprived of its ith element.

As K = Σ -1 , reordering the rows and columns of the matrices yields

where K \i\i is matrix K deprived of its ith column and its ith line, K i\i is the ith column of matrix K deprived of its ith element, and I p-1 is the identity matrix of size p -1. Two of these blockwise equalities are rewritten as follows: