Inferring sparse Gaussian graphical models with latent structure

Our concern is selecting the concentration matrix's nonzero coefficients for a sparse Gaussian graphical model in a high-dimensional setting. This corresponds to estimating the graph of conditional dependencies between the variables. We describe a novel framework taking into account a latent structure on the concentration matrix. This latent structure is used to drive a penalty matrix and thus to recover a graphical model with a constrained topology. Our method uses an $\ell_1$ penalized likelihood criterion. Inference of the graph of conditional dependencies between the variates and of the hidden variables is performed simultaneously in an iterative \textsc{em}-like algorithm. The performances of our method is illustrated on synthetic as well as real data, the latter concerning breast cancer.


Introduction
Estimating the concentration matrix (namely the inverse of the covariance matrix) of a Gaussian vector in a sparse, high-dimensional setting has received much attention recently. Graphical models provide a convenient setting for modelling multivariate dependence patterns. In this framework, an undirected graph is matched to the Gaussian random vector, where each vertex corresponds to one coordinate of the vector, and an edge is not present between two vertices if the corresponding random variables are independent, conditional on the remaining variables. Now, conditional independence between two coordinates of the Gaussian random vector corresponds exactly to a zero entry in the concentration matrix. Thus, detecting nonzero elements in the concentration matrix is equivalent to reconstructing the Gaussian graphical model (GGM, see e.g. Lauritzen 1996).
We focus here on the crucial problem of selecting the concentration matrix's nonzero coefficients. In other words, we focus on variable selection rather than estimation. Application areas include gene-regulation graph inference in Biology (using gene expression microarray data), as well as spectroscopy, climate studies, functional magnetic resonance imaging, etc. We provide a very novel approach driving the graph selection according to an unobserved modular structure on the vertices.
The idea of covariance selection first appeared in the work of Dempster (1972). In the so-called "large p, small n" setting (namely when the number of observations is smaller than the dimension of the observed response), the need for covariance selection is huge, as the empirical covariance matrix is no longer regular.
In Drton and Perlman (2007), a classification of the different methods for model selection/estimation in GGMs into three group types is suggested: constraint-based methods, performing statistical tests; Bayesian approaches; and scorebased methods, maximizing a model-based criterion. The multiple testing problem has been taken into account in Drton and Perlman (2007;. The authors perform GGM covariance selection by multiple testing of hypotheses about vanishing partial correlation coefficients. Such procedures may also be implemented using the PC-algorithm (Kalisch and Bühlmann 2007). Starting from a complete, undirected graph, the PC-algorithm deletes edges recursively, according to conditional independence decisions. However, the statistical procedure in Drton and Perlman (2007) relies on asymptotic considerations, a regime never attained in real situations.
Another attempt in this vein is to consider limited-order partial correlations (Wille andBühlmann 2006, Castelo andRoverato 2006). In Wille and Bühlmann (2006) the authors consider only zero and first-order conditional dependencies. They argue that for sparse graphical models, these low-order dependencies still reflect reasonably well the full-order conditional dependency structure. Moreover, these dependencies may be well estimated even with a small number of observations. In Castelo and Roverato (2006), the authors introduce a nonrejection rate to reduce the multiple testing and computational problems to which these approaches give rise.
A Bayesian framework is proposed in Dobra et al. (2004), Jones et al. (2005) and their method was applied for evaluating patterns of association in large-scale gene expression data. The approach is based on dependency networks, namely a collection of conditional distributions {P(X i |X \i )} (where X \i stands for the set of all variables but X i ). However, such conditional distributions will not in general result in a coherent joint distribution. This point is further discussed below concerning a similar problem appearing in Meinshausen and Bühlmann (2006). Moreover, constructing priors on the set of concentration matrices is not a trivial task (it mainly relies on Wishart priors). The use of MCMC procedures limits the range of applications to moderate-sized networks.
Before focusing on score-based methods, let us first introduce regularization procedures. In the context of linear regression, the Lasso (least absolute shrinkage and selection operator) technique was introduced by Tibshirani (1996). This procedure performs model selection and parameter estimation at the same time. The idea is that ordinary least-squares criterion may be improved in a sparse context, using an 1 -norm penalty. The 1 -norm penalty shrinks the estimates to zero while preserving the convexity of the optimization problem. Note that the 1 -norm penalization is also known as 'basis pursuit' in signal processing (Chen et al. 2001).
It is well known that if the ultimate goal is parameter estimation, model selection and estimation should be done in a single step. Performing the model selection prior to parameter estimation in the selected model will, in fact, result in a non-robust procedure. However, our primary focus here is on model selection, as we want to infer sparse networks. We therefore concentrate on model selection rather than on estimation performances.
The Lars algorithm (Efron et al. 2004) is one of the most popular techniques for solving the Lasso problem. It gives the path of solutions obtained when varying the penalty parameter (the penalty parameter is used as a scaling factor of the 1 -norm penalty). The larger the penalty parameter, the sparser the Lasso solution.
Using convex optimization techniques (see for instance Minoux 1986), the Lasso problem may be stated as a primal problem, whose dual formulation may be solved more easily. This approach is taken in Osborne et al. (2000b). The authors obtain an iterative algorithm, the "homotopy method" (Osborne et al. 2000a). Other very efficient approaches are based on focusing on each coordinate iteratively. Indeed, for each coordinate, the Lasso problem is solved very simply (assuming the other coordinates are fixed) by soft-thresholding (Donoho and Johnstone 1995). Thus, different 'coordinate optimization' procedures have been proposed in the literature. Following the work of Fu (1998), a cyclic procedure is proposed in Friedman et al. (2007), where optimization with respect to each coordinate is done iteratively; whereas Wu and Lange (2008) propose a greedy approach, computing the solution for each coordinate and choosing that which provides the largest decrease in a surrogate objective function. Note that these approaches rely on the underlying assumption that the predictors for the regression problem are uncorrelated.
Let us now come back to covariance (or concentration) matrix inference in GGMs, using maximization of a model-based criterion. Meinshausen and Bühlmann (2006) were the first authors to apply Lasso techniques for inferring a covariance matrix in a GGM. Their approach is to solve p different Lasso regression problems, where p is the dimension of the observed vector. The main drawback of such a procedure is that a symmetrization step is required to obtain the final network. It might, for instance, be the case that the estimator of the regression coefficient for X i on X j is zero, whereas the estimator for X j on X i is not zero. Meinshausen and Bühlmann propose to use either an "AND" or an "OR" final step procedure to recover an undirected correlation graph. However, these two procedures might result in different estimates and there is no way of choosing between them. Moreover, as previously stated, a set of conditional distributions does not necessarily cohere into a joint distribution. Using a set of possibly non-coherent conditional distributions corresponds to a pseudo-likelihood approach. This aspect was not underlined in Meinshausen and Bühlmann (2006), and we clarify this point in Section 4.
Subsequently, two other articles, Banerjee et al. (2008) and Yuan and Lin (2007) independently provided an improvement of the initial work of Meinshausen and Bühlmann (2006). In both works, the problem is seen as a penalized maximum-likelihood (PML) problem. Instead of considering p different regression problems, these two articles focus on the likelihood of the Gaussian vector, penalizing the entries of the concentration matrix with an 1 -norm penalty. They explain how the PML estimation may be solved as a "Lasso-like" problem. The major issue with PML strategies in the context of the concentration matrix estimation of GGMs is to obtain a positive definite estimate. However, the approach for solving the problem in Yuan and Lin (2007) is not suited to high-dimensional settings, in contrast to the approach proposed in Banerjee et al. (2008). In Yuan and Lin (2007), a non-negative garrote-type estimator is used, and asymptotic properties (as n tends to infinity while p is held fixed) are given. In Banerjee et al. (2008), two different algorithms are proposed for solving problems in a high-dimensional setting. The first approach relies on a block-coordinate descent algorithm. The second is a semi-definite programming algorithm, based on Nesterov's method, which is computationally intensive.
The next improvement in this vein comes with Friedman et al. (2008). Relying on coordinate descent techniques, previously described in Friedman et al. (2007), the authors revisit Banerjee et al.'s first approach and propose an efficient algorithm to solve the PML estimation problem, under the positive definite constraint. In fact, they use the block coordinate descent approach proposed by Banerjee et al. (2008) and combine it with a second coordinate descent method. Our method will make use of this approach.
To conclude this part, we remark that a completely different shrinkage estimate was proposed by Schäfer and Strimmer (2005) in the same context of largescale covariance matrix estimation. The approach consists in using a weighted average of two different estimators, the first being unconstrained (thus having small bias but large variance), the second being low-dimensional (and thus exhibiting small variance but large bias). Now let us motivate the use of hidden structures in networks. Modularity is a property observed in real (biological) networks (see for instance Ihmels et al. 2002). Heterogeneity in the node behaviors is an important property of these data. For example, so-called 'hubs' are highly-connected nodes, showing a different behavior from the rest of the graph nodes. An interesting model capturing these network features is a mixture model for random graphs (see for instance Daudin et al. 2008). This model has been rediscovered many times in the literature, and a non-exhaustive bibliography should include Frank and Harary (1982), Snijders and Nowicki (1997), Nowicki and Snijders (2001), Tallberg (2005), Daudin et al. (2008), Mariadassou andRobin (2007), Zanghi et al. (2008). To state it simply, this model assumes that each node belongs to some unobserved group. Conditional on the node groups, the (weighted) edges are independent and identically-distributed (i.i.d.) random variables, whose distribution depends on the groups of the nodes to be connected. As we are interested in GGMs, weighted edges correspond to entries of the concentration matrix.
In this work, we aim at estimating a hidden structure, namely node groups, while discovering the network. This hidden structure should help us in choosing adaptive penalty parameters. Indeed, we wish to penalize the elements of the concentration matrix, according to the unobserved clusters to which the nodes belong. For instance, if two nodes belong to the same unobserved group, we wish to lower the penalty parameter acting on the corresponding entry in the concentration matrix. Conversely, if we increase the penalty parameters on the entries corresponding to nodes belonging to different groups, we shrink the estimated coefficient to zero. Our approach is completely new and improves inference of sparse modular networks.
Another adaptive Lasso procedure is given in Zou (2006), whose idea is to lower the bias of the large coefficients by adapting the penalty parameter of each coefficient so that it automatically scales with the inferred value. It is known that the non-adaptive Lasso procedure may result in inconsistent parameter estimation. An illustration of the conflict between optimal prediction and consistent variable selection for the Lasso procedure is given in Meinshausen and Bühlmann (2006). They proved that the optimal penalty parameter for prediction gives inconsistent variable selection results, motivating the use of another penalty parameter to ensure the control of the probability of falsely connecting two or more distinct connectivity components of the graph. Like them, we also focus on optimal selection rather than on optimal prediction. The adaptivity of our procedure is not used for lowering the bias of large coefficients, but instead for constraining the prediction to fit the underlying structure of the graph.
Model. Let us now briefly describe the general approach of our work. The model will be presented in detail in Section 2. Let X = (X 1 , . . . , X p ) be a Gaussian random vector in R p , with zero mean and positive definite covariance matrix Σ, namely X ∼ N (0 p , Σ). We observe independent and identicallydistributed (i.i.d) vectors (X 1 , . . . , X n ) with the same distribution as X. The matrix K = Σ −1 is the concentration matrix of the model. Let S be the empirical covariance matrix. The log-likelihood of the observations is given by where c is a constant term. The 1 -penalized estimator proposed by Banerjee et al. (2008) is given by where K 0 stands for positive definiteness, ρ > 0 is a penalty parameter and K 1 = ij |K ij |.
A natural generalization of this approach is to have different penalty parameters for different entries K ij . Namely, where ρ(K) = (ρ ij (K ij )) i,j∈P is a matrix of penalty functions acting on each entry. As a general rule, using as many penalty functions as there are entries in the concentration matrix to be estimated is not meaningful.
Here, we propose to take into account a hidden structure on the correlations between the coordinates random variables X k . Thus, we consider latent i.i.d. random variables Z 1 , . . . , Z p with values in a finite set {1, . . . , Q}. Each variable Z i describes the state of X i , and we wish to adapt the penalty function ρ ij with respect to the states of X i , X j . More precisely, we wish to use a criterion of the form where ρ Z (K) = (ρ ZiZj (K ij )) i,j∈P is a matrix of random penalty functions whose entries depend on the latent structure Z = Z 1 , . . . , Z p . However, the hidden structure is not supposed to be known, thus we cannot rely on the previous criteria. Intuitively, following the principle of Expectation-Maximization (em) algorithm of Dempster et al. (1977), the idea will be to replace the unobserved value ρ Z (K) with its conditional expectation E(ρ Z (K)|X 1 , . . . , X n ; K (m) ) under some model with parameter K (m) , and iterate the following steps (e) Compute E(ρ Z (K)|X 1 , . . . , X n ; K (m) ) (m) Update K (m+1) = argmax K 0 E(ρ Z (K)|X 1 , . . . , X n ; K (m) ).
One of our aims is to provide a very simple framework for such an analysis. Note that the 1 -norm used here acts on diagonal elements of the matrix K. It is counter-intuitive to penalize diagonal elements of the concentration matrix, as these do not reflect sparsity in the correlation structure. However, from a technical point of view, this strategy ensures that the procedure will select a positive definite estimator (see Remark 2). This point was not emphasized in the previous procedures using 1 penalized likelihood of GGMs.
Road-map. In Section 2 we present the model and the penalized maximumlikelihood criterion on which we base our inference procedure, described in Section 3. This procedure relies on a variational em algorithm, combined with a Lasso-like procedure. Section 4 explains how Meinshausen and Bühlmann's approach may be interpreted as a penalized pseudo-likelihood method. Finally, Section 5 illustrates the performance of the method on synthetic data, for which an R-package, SIMoNe (Statistical Inference for Modular Network), can be downloaded from the first author's website. We also test our algorithm on a real data set provided by Hess et al. (2006) and concerning n = 133 patients with breast cancer treated using chemotherapy. According to Hess et al. (2006) and Natowicz et al. (2008), the patient response to the treatment can be classified either as a pathologic complete response (pCR), or as a residual disease (not-pCR). The prediction of the patient response is achieved accurately by studying the expression levels of a limited number of genes (p = 26). Our algorithm is applied on each class of patients (pcR and not-pCR). Two distinct gene-regulatory networks are thus inferred, showing a very different structure according to the selected class of patients.

A latent structure model for network inference
In this section we present a framework for modelling heterogeneity among dependencies between the variables. To this end, let us first recall classical notations from Gaussian Graphical Models (see Lauritzen 1996, for elementary results about GGMs).

Gaussian graphical models: general settings
Let P = {1, . . . , p} be a set of fixed vertices, X = (X 1 , . . . , X p ) a random vector describing a signal over this set and a sample (X 1 , . . . , X n ) of size n with the same distribution as X.
The vector X is assumed to be Gaussian with positive definite covariance matrix Σ = (Σ ij ) (i,j)∈P 2 . No loss of generality is involved when centering X, so we may assume that X ∼ N (0 p , Σ). GGMs are based on a classical result, originally emphasized by Dempster (1972), claiming that any couple of entries (X i , X j ) with i = j are independent conditional on all other variables indexed by P\{i, j}, if and only if the entry (Σ −1 ) ij is zero. The inverse of the covariance matrix K = (K ij ) (i,j)∈P 2 = Σ −1 , known as the concentration matrix, thus describes the conditional independence structure of X. Moreover, each entry K ij , i = j is directly linked to the partial correlation coefficient r ij|P\{i,j} between variables X i and X j . In fact, we have r ij|P\{i,j} = −K ij / K ii K jj , and also K ii = Var(X i |X P\i ) −1 . Hence, after a simple rescaling, the matrix K can be interpreted as the adjacency matrix of an undirected weighted graph G representing the partial correlation structure between variables X 1 , . . . , X p . This graph has no self-loop, with a random set of edges composed by all pairs (i, j) such that K ij = 0. Note that we are seeking only pairs of vertices (i, j) such that i < j, since there is no self-loop, and since K ij = K ji . Inferring nonzero entries of K is equivalent to inferring G, and is therefore a highly relevant issue in this framework.

Providing the network with a latent structure
Let us now extend the modeling by providing the network with an internal latent structure.
The model proposed in Daudin et al. (2008) attempts a better fit of data, as it places the network G in the mixture framework, in order to take account of the heterogeneity among vertices. The same general mixture model is adopted here: vertices of P are distributed among a set Q = {1, . . . , Q} of hidden clusters that model the latent structure of the network. For any vertex i, the indicator variable Z iq is equal to 1 if i ∈ q and 0 otherwise, hence describing which cluster the vertex i belongs to. A vertex is assumed to belong to one cluster only, thus the random vector Z i = (Z i1 , . . . , Z iQ ) obviously follows a multinomial distribution. Namely, where α = (α 1 , . . . , α Q ) is a vector of cluster proportions, so that q α q = 1.
The concentration matrix structure. We shall now extend the clustering of vertices from P to the concentration matrix K. Accordingly, both the existence and the weight of edges, described by the off-diagonal elements of K, will depend on the cluster each vertex belongs to. Conditional on the events i ∈ q and j ∈ where q, are clusters chosen from Q, each K ij (i = j) is a random variable whose probability density function is denoted by f q , that is, It will be remarked that in this formulation the variables K ij are assumed to be independent, conditional on the clusters the vertices belong to. Moreover, we are considering only undirected graphs, so we may assume that f q = f q . For technical reasons (see Remark 2), we also assume a prior distribution on diagonal elements of K, namely Our suggestion is to adopt Laplace distributions; hence where λ q , λ 0 > 0 are scaling parameters and λ q = λ q . Below, the parameter λ 0 will be fixed and not estimated. The reason for choosing a Laplace distribution is that it is reminiscent of the 1 -norm, itself linked to Lasso-techniques for which appropriate tools are available. In fact, when considering the general penalized least-square problem, the penalty term can be seen as a log-prior density on the vector of parameters. In the case of Lasso, the prior distribution corresponding to the 1 -norm is actually the Laplace distribution (see, e.g. Hastie et al. 2001).
The affiliation model. The affiliation model is a special case of network structure (to be investigated below), where there are many different clusters, but where the focus is restricted to two types of edges: edges between nodes of the same cluster, and edges between nodes from different clusters. In the affiliation model the densities f q in (4) are of only two kinds; that is, for all q, ∈ Q, let if q = , the intra-cluster density of edges, f q = f out (·; λ out ) if q = , the inter-cluster density of edges. (5)

The complete likelihood
Having described the modeling of the network, we now focus on the inference issue. We denote as X the n×p matrix that contains the data-set {X 1 , X 2 , . . . , X n } row-wisely organized, i.e., (X k ) is the kth row of X. Furthermore, we denote as Z = {Z iq } i∈P,q∈Q the set of all latent indicator variables for vertices. For the sake of simplicity, the number of clusters Q and the parameters α = (α q ) q∈Q and λ = {λ q } q, ∈Q are assumed to be known for the moment.
The data experiments X are the only observations available, and from these we should like to be able to infer the graph G of conditional dependencies or, equivalently, nonzero entries of K. As the matrix K has been given a prior distribution, our aim is to maximize the posterior probability of K, given the data X, or equivalently, the logarithm of the joint distribution. The estimate is thus defined as follows: where K 0 stands for positive-definiteness.
To solve this problem, we place ourselves in the classical complete-data framework. The distribution of K is only known conditionally on the latent structure described by Z. We denote as Z the set of all possible clusterings over nodes from P. The marginalization over the latent clusters Z leads to where the so-called complete-data likelihood L c (X, K, Z) = P(X, K, Z) is the function we shall develop using an em-like strategy hereafter. For this purpose, a closed form of L c is required.
Proposition 1. The following relation holds for the complete-data likelihood L c .
Proof. Using the Bayes rule, L c divides into three terms: where we make use of the fact that log P(X|K, Z) = log P(X|K).
The first term is the likelihood associated with a size-n sample of a multivariate Gaussian distribution, since X ∼ N (0 p , Σ). Routine computations lead to log P(X|K) = n 2 log det(K) − n 2 Tr(SK) − np 2 log(2π).
As regards the second term, using the expression (4), we have From (2), we have log P(Z) = i,q Z iq log α q , and the result follows.

Inference strategy by alternate optimization
In the classical em framework developed by Dempster et al. (1977), where X is the available data, inferring the unknown parameters K spread over a latent structure Z would make use of the following conditional expectation: where K (m) is the estimation of K from the previous step of the algorithm.
The usual em strategy would be to alternate an E-step computing the conditional expectation (8) with an M-step maximizing this quantity over the parameter of interest K. Unfortunately, no closed form of Q K|K (m) can be formulated in the present case. The technical difficulty lies in the complex dependency structure contained in the model. Indeed, P(Z|K) cannot be factorized, as argued in Daudin et al. (2008). This makes the direct calculation of Q K|K (m) impossible. To tackle this problem we use a variational approach (see, e.g., Jaakkola 2000, for elementary results on variational methods). In this framework, the conditional distribution of the latent variables P(Z|K (m) ) is approximated by a more convenient distribution denoted by R m (Z), which is chosen carefully in order to be tractable. Hence, our em-like algorithm deals with the following approximation of the conditional expectation (8) In the following section we develop a variational argument in order to choose an approximation R m (Z) of P(Z|K (m) ). This enables us to compute the conditional expectation (9) and proceed to the maximization step.

Variational estimation of the latent structure ( E-step)
In this part, K is assumed to be known, and we are looking for an approximate distribution R(·) of the latent variables. The variational approach consists in maximizing a lower bound J of the log-likelihood log P(X, K), defined as follows: where D KL is the Küllback-Leibler divergence. This measures the difference between the probability distribution P(·|K) in the underlying model and its approximation R(·). An intuitively straightforward choice for R(·) is a completely factorized distribution (see Mariadassou andRobin 2007, Zanghi et al. 2008) where h τ i is the density of the multinomial probability distribution M(1; τ i ), and τ i = (τ i1 , . . . , τ iQ ) is a random vector containing the variational parameters to optimize. The complete set of parameters τ = {τ iq } i∈P,q∈Q is what we are seeking to obtain via the variational inference. In the case in hand the variational approach intuitively operates as follows: each τ iq must be seen as an approximation of the probability that vertex i belongs to cluster q, conditional on the data, that is, τ iq estimates P(Z iq = 1|K), under the constraint q τ iq = 1. In the ideal case where P(Z|K) can be factorized as i P(Z i |K) and the parameters τ iq are chosen as τ iq = P(Z iq = 1|K), the Küllback-Leibler divergence is null and the bound J reaches the log-likelihood.
The following proposition gives the form of the lower bound J to be maximized in order to estimate τ .
Proposition 2. Let us assume that R τ can be factorized as in (11), and let us denote J τ (X, K) := J (X, K, R τ (Z)). Then J τ satisfies the following expression where c does not depend on τ and ρ τ (K) = (ρ τ iτ j (K ij )) i,j∈P 2 is defined similarly as (7), replacing Z iq by τ iq .
Proof. Starting from (10), classical results on variational methods show that where H(R τ (·)) is the entropy of the distribution R τ (·) and Q τ (K) is the approximation of the complete log-likelihood conditional expectation, computed under the distribution R τ . Namely, In the special case of factorized distribution (11), the entropy is Moreover, Equation (12) follows via Proposition 1, by using that E Rτ (Z iq ) = τ iq and The optimal approximate distribution R τ is then derived by direct maximization of J τ . The following proposition gives the estimate τ that solves the problem.
Proposition 3. Let α and λ be known. The following fixed-point relationship holds for the optimal variational parameters τ = arg max τ J τ where ∝ means that there is a scaling factor such that for any i ∈ P, we have q τ iq = 1. Proof. This is just an adaptation to the Laplace case of Mariadassou and Robin (2007, Proposition 3).
The initial value of τ is chosen using a classification algorithm such as spectral clustering (see for instance Ng et al. 2002). As a consequence, the initial values for τ iq lie in {0, 1}. We then use an iterative procedure setting τ (m+1) = g( τ (m) ), where g is the function (implicitly defined above) for which τ is a fixed point. Note that we cannot ensure uniqueness of the fixed point for g, nor convergence of this iterative procedure. In practice, we can always use a maximal number of iterations, and if convergence has not occurred, we keep the initial value of τ given by the clustering method. In appendix A.2 we explain that at least in the affiliation model (5), if the current values K (m) ij of the precision matrix are small enough, and if the penalty parameters λ −1 in and λ −1 out are well-chosen, then uniqueness of the fixed point is ensured. However, such a result does not hold in the general case, which is one of the drawbacks of of the variational approach in this context.
Estimation of α and λ. The parameters α and λ have been previously considered as known to keep the statement as clear as possible.
Two different strategies may be used with respect to these parameters. The first approach is to fix their values. Fixing the value of α comes down to choosing a priori the proportions of the groups, which is quite a common strategy in mixture models. As for the choice of λ, this is equivalent to choosing the penalty parameter in the classical lasso. Concerning general parameters λ, a number of values need to be determined, which might be a problem. However in the particular affiliation model (5), only 2 parameters have to be fixed: a parameter λ in that corresponds to a light penalty, since many intra-cluster edges are expected, and another parameter λ out that fits with a heavier penalty, since we do not expect many inter-cluster edges. This is typically the kind of strategy that will be used for numerical applications (see Section 5). More generally, the matrix penalty can be tuned to obtain a desired quantity of inferred edges, or to constrain the topology of the graph, e.g. graphs with hubs.
The second strategy is to make use of the current inferred graph to estimate the parameters. The basic idea is to include this estimation in the variational method. Unfortunately, the maximization of J τ given in equation (12) with respect to τ , λ and α at the same time is not possible. To tackle this problem, we use an alternate strategy. The parameter τ is computed with the fixedpoint relationship (14) for fixed values of λ and α. Then we maximize J τ with respect to λ and α, once R τ is fixed (that is, once τ is fixed), as in the following proposition. We successively iterate these two steps until stabilization.
Proposition 4. For fixed values of τ , the parametersα,λ maximizing J τ are given by Proof. Once terms that do not depend on the parameters of interest have been removed from J τ , the problem becomeŝ Null-differentiation with respect to α q (under the constraint q α q = 1) and λ q leads straightforwardly to the result.

A Lasso-like method to estimate the concentration matrix (the M-step)
Now that we are able to compute the approximate conditional expectation Q τ (K) defined by (13), we wish to infer the concentration matrix K, assuming τ is known. This is the aim of the m-step of our em-like strategy, that deals with the maximization problem arg max K 0 Q τ (K). Using Proposition 1 and the equality E Rτ (Z iq Z j ) = τ iq τ j , it is a simple matter to rewrite the problem as follows Hence, our m-step can be seen as a penalized maximum likelihood estimation problem, exactly like in Friedman et al. (2008), Banerjee et al. (2008). The likelihood considered here is P(X|K), that is, the likelihood which corresponds to the n realizations of the Gaussian vector X for a given concentration matrix K. The difference of our approach lies in the complexity of the penalty term, and in slight discrepancies as regards some constant factors. Remark 1. Since we are using a penalty term 1/λ 0 on matrix K's diagonal elements, the solution to (15) satisfies when λ −1 0 < n|S ii |/2 for any i ∈ P. Indeed, the sub-gradient equation is n/2(K −1 ii − S ii ) + sgn(K ii )/λ 0 = 0, and K ii ≥ 0 since it is the inverse of a conditional variance.
Let us now look at the solution of the m-step: the following proposition gives an equivalent formulation of (15) that is more likely to be solved. The result draws its inspiration from Banerjee et al. (2008).
Proposition 5. The maximization problem (15) over the concentration matrix K is equivalent to the following, dealing with the covariance matrix Σ Σ = argmax where · is the term-by-term division and Remark 2. By penalizing the diagonal terms of the concentration matrix K in the initial problem, the set of matrices Σ over which we maximize our criterion contains, for instance, the matrix S + 2/(nλ 0 )I, (where I stands for the identity matrix). Thus, provided that the value of penalty parameter 1/λ 0 is set sufficiently high, this set contains positive definite matrices. This ensures that our estimator is always invertible. Obviously, when S is invertible, which is usually true for n greater or equal than p, penalizing the diagonal terms becomes futile.
In this case 1/λ 0 is set to zero.
Proof. The penalty term in (15) can be written as follows where is the term-by-term product. The set {T q } q, ∈Q contains p × p symmetric matrices, defined, for each couple (q, ), by Let us now use the fact that A 1 = max U ∞ ≤1 Tr(AU), for a given matrix A. The optimization problem (15) can now be written as since the trace operator is linear. The dual version of the above expression is obtained by swapping max and min. The maximization is solved by differentiating with respect to K. To do this, we recall that in our specific case the matrices T are symmetrical, and thus Tr ((T K)U) = Tr (K(T U)). Then, applying the usual rules for the derivative of the trace operator, null-differentiation with respect to K yields The dual problem therefore becomes or in other words, max log det(Σ).
Finally, we need to write the constraint as a function of Σ rather than the set {U q }. In fact, we simply need to show that which is straightforward (see Appendix A.1 for details).
To solve (17) and thus obtain the estimate Σ, we successively use two coordinate descent methods. The first corresponds to a block-wise strategy suggested by Banerjee et al.. The second one is used to solve the resulting Lasso problem and was suggested by Friedman et al. (2007).
Let us first explain the block-wise strategy. For this purpose, we introduce the following notation for Σ, S and the penalty matrix P τ Σ = Σ 11 σ 12 σ 12 Σ 22 , S = S 11 s 12 s 12 S 22 , P τ = P 11 p 12 p 12 P 22 , where Σ 11 , S 11 and P 11 are (p−1)×(p−1) matrices, σ 12 , s 12 and p 12 are (p−1) length column vectors and Σ 22 , S 22 and P 22 are real numbers. We have already remarked (Remark 1) that the solution to (17) We have det( Σ) = det( Σ 11 )( Σ 22 − σ 12 Σ −1 11 σ 12 ). The full matrix Σ is approximated in the following way: first, if required when p is greater than n, we initialize the procedure with S + 2/(nλ 0 )I, where λ 0 > 0 is chosen so as to make S + 2/(nλ 0 )I invertible; secondly, we permute the columns (and thus the rows) of Σ and iteratively solve problems like (20) until convergence of the procedure. This convergence is ensured by the following lemma.
Lemma 1. The procedure which starts with a positive definite matrix and iteratively updates the columns and rows of this matrix according to the solutions of (20) converges to the solution Σ of (17).
Proof. The proof relies on Banerjee et al. (2008, Theorem 3) and Tseng (2001, Theorem 4.1). Convergence of block-coordinate descent methods is a well-documented topic in convex optimization literature. Here, we have to bear in mind that using 1 -norm penalty leads to non-differentiable functions. Thus, we rely on a result by Tseng (2001, Theorem 4.1), which in our case ensures the convergence of the procedure, provided there is at most one solution to each minimization problem (20). This point is proved in Banerjee et al. (2008, Theorem 3).
Then, starting from a result given in Banerjee et al. (2008), an interpretation of (20) as an 1 -penalized problem is given in Friedman et al. (2008). This 1penalized problem is reminiscent of the Lasso and may thus be solved using a coordinate descent strategy (Friedman et al. 2007). The following proposition enunciates a result similar to those obtained in Banerjee et al. (2008, equation (6)) and Friedman et al. (2008, equation (2.4)), although with a more general penalty term and a factor 1 2 that differs. Since none of these articles gives an explicit proof for this result, it is fitting that we provide our own proof here.
where solution σ 12 to (20) and β to (21) are linked through Proof. Problem (20) can be written as follows, by splitting the constraint: Let us introduce L the so-called Lagrangian, with vectors of Lagrange coefficients denoted by β 1 = (β 1 i ) i≤p−1 , β 2 = (β 2 i ) i≤p−1 with nonnegative entries. Also, let β = β 2 − β 1 . The Lagrange version of the above problem is where, in the present case, L is given by The coefficients β 1 i and β 2 i maximizing L(β) are null when the constraints are satisfied, and for each index i, at least one coefficient among Meanwhile, consider the dual problem of (23), swapping min and max: the solution that minimizes the dual problem with respect to y satisfies the nullgradient hypothesis. We obtain 2 Σ −1 11 y − β = 0, that is y = 1 2 Σ 11 β (which proves equation (22)). Introducing this result in the dual of (23), we get also equivalent to min β 1 4 β Σ 11 β − s 12 β + p 12 β 1 .
Expressing this quantity by using the Euclidean norm achieves the proof.
Hence, the column σ 12 of the estimated covariance matrix Σ is computed by solving the Lasso problem (21) using another coordinate descent method. (21) is computed by updating the jth coordinate of β via
Proof. The proof of this lemma is postponed to Appendix A.3.
Finally, the estimate of the matrix of concentration K is recovered by inverting Σ, which can be done at low computational cost (see appendix A.4 for details). Hence, we solve the initial maximization problem (15) that defines the m-step of our algorithm.
Implementation of the full em algorithm is outlined in Algorithm 1.

Choice of penalty parameters
As previously stated, the penalty parameters λ may be estimated in the e-step of the algorithm (see subsection 3.1). However, this choice is not necessarily optimal for the estimation of K, and other choices might in practice lead to a better solution. A good strategy is to keep the estimated value of λ in the e-step that leads to the estimation of τ , and to impose another value of λ during the m-step. In this part, we indicate a possible choice for the penalty parameters to use in the m-step, ensuring a small error on the connectivity components of the estimated graph.
Let us first introduce some notation. For any node i ∈ P, let C i denote the connectivity component of node i in the true underlying conditional dependency graph, and C i the corresponding component resulting from the estimate K of this graph structure. The following proposition is based on Meinshausen and Bühlmann (2006, Theorem 2) and Banerjee et al. (2008, Theorem 2).
Proposition 7. Fix some ε > 0 and choose the penalty parameters λ such that, for all q, ∈ Q, where 1 − F n−2 is the c.d.f. of Student's t-distribution with n − 2 degrees of freedom. Then Proof. Here we simply indicate the main differences between the proof of Banerjee et al. (2008, Theorem 2) and what is valid in our context. Note that according to (15), the estimator K must satisfy the following sub-gradient equation where ν ij ∈ sgn( K ij ). Following the proof of Banerjee et al. (2008, Theorem 2), we easily get Performing some computations involving the correlation between variables X i and X j , we also obtain which entails the conclusion.
Remark 3. Following Banerjee et al. (2008), note that in order to ensure (25), it is enough to choose the penalty parameter λ such that, for all q, ∈ Q, degrees of freedom, i.e. F n−2 (t n−2 (u)) = u. Remark 4. Inequality (25) does not take into account that different penalty parameters are used for different hidden classes q, ∈ Q. An adaptation of the preceding strategy is to use current values Z (m) obtained from the probabilities τ (m) of the hidden classes and to choose the current penalty parameters λ (m) accordingly. More precisely, let us set, for instance Then, when for all q, ∈ Q, the current estimate K (m) of the dependency graph will approximately satisfy (26). Moreover, in order to ensure (27), it is enough to choose, for all q, ∈ Q, (28) Typically, the kind of values obtained with (28) will lead to large penalties and, consequently, to very sparse graphs: practically, more informative networks can be obtained by replacing the term ε/2p 2 in (28) by greater values. In any cases, (28) should be seen as a starting value.

Link with Meinshausen and Bühlmann's approach
We should also like to fill the gap between, on the one hand solving (15) and, on the other, the approach proposed in Meinshausen and Bühlmann (2006), where p independent penalized regression problems are solved using the Lasso. In fact, we shall show that Meinshausen and Bühlmann's approach is equivalent to maximizing the penalized pseudo log-likelihood corresponding to the sizen sample of the multivariate Gaussian vector X on the set of non symmetric matrices. Let us denote as L this pseudo-likelihood, defined by where X k P\i is the kth realization of the Gaussian vector X, once the ith coordinate has been removed. In this section, the 1 -norm of matrices is restricted to off-diagonal elements only, that is, Proposition 8. Consider the solution K pseudo to the penalized pseudo-likelihood problem (whose diagonal is fixed) and the solution K MB given in Meinshausen and Bühlmann (2006) to the p different regression problems, using the matrix penalty 2P/n. The two solutions have exactly the same null entries.
Proof. Denote by K \i\i and S \i\i , respectively, the matrices K and S once their ith row and ith column have been removed. Moreover, K i\i and S i\i are the ith rows of the matrices with the ith term removed. After some routine computations, and using classical results for Gaussian multivariate vectors (see Appendix A.5), it can be shown that where c does not depend on K. Thus, if we forget the symmetry constraint on K, maximizing the pseudo-likelihood (30) with respect to the non-diagonal entries of K is equivalent to p independent maximization problems with respect to each column K i\i . Consider, for instance, the last column of K, that is, for i = p, and the relative term in (30). This term can be written as where we use the block-wise notation defined above (19). The term C does not depend on K i\i , which is the current column of the concentration matrix to infer. Namely, c = −K 2 22 s 12 S −1 11 s 12 . Consider now the penalized version of the log-likelihood (29): we wish to solve p penalized problems of minimization as defined above, which can be written as follows min Meinshausen and Bühlmann wish to solve p Lasso-problems, for instance for the last variable p, where X p is the pth column of X and X \p is the matrix of data the pth column has been removed (note that we adapted the penalization term corresponding to the framework developed here). The minimum is reached in (31) for null-differentiation, and we get 2S 11 β + 2K 22 s 12 + 2K 22 n p 12 ν = 0, where ν ∈ sign(β). The same for (32), and we get where γ ∈ sign(α). Now, just note that n −1 X \p X \p = S 11 and n −1 X p X \p = s 12 , and problems (31) and (32) are equivalent, provided that α = −β/K 22 . Thus, the columns of the concentration matrix (with a removed diagonal term) inferred from the penalized maximum pseudo-likelihood problem (29), and those inferred with Meinshausen and Bühlmann's approach, share exactly the same null-entries, that is, the same network of conditional dependencies.

Numerical experiments
In this section we present numerical experiments on both synthetic data, to investigate how well the proposed selection procedure behaves, and real data, to demonstrate the practical use of GGM covariance selection with latent structure. In the remainder of this section we focus on an affiliation model (5), the choice of the penalty being made in line with Section 3.3. More precisely, we fix the ratio λ in /λ out = 1.2 and either let the value 1/λ in vary when considering precision/recall curves for synthetic data, or fix this parameter according to (28) when dealing with real data.

Synthetic data
We perform numerical experiments to assess the performance of our approach (SIMoNe, Statistical Inference for Modular Network) and compare it to already existing methods for GGM covariance selection: GLasso (Friedman et al. 2008) and GeneNet (Schäfer and Strimmer 2005).
Data synthesis in our framework requires the simulation of a structured sparse inverse covariance matrix. To this aim, we first simulate a graph with an affiliation structure. We consider a simple binary affiliation model where two types of edges exist: edges between nodes of the same class and edges between nodes of different classes. The binary incidence matrix of the graph is transformed by randomly flipping the sign of some elements in order to simulate both positively and negatively correlated variables. Positive definiteness of this matrix is ensured by adding a large enough constant to the diagonal. The matrix is then further normalized to have a diagonal of ones. A Gaussian sample of size n with zero mean and the above covariance matrix is then simulated 50 times. The results we present below are averaged over the 50 samples. At the end of this section we discuss the performances of our method when there is no latent structure on the data. We simulate sparse graphs with p = 200 and n from 100 to 2000 (n/p ∈ {1/2, 2, 3, 6, 10}). We use a probability of intra-cluster connection of 0.125, a probability of inter-cluster connection of 0.0025, Q = 3 groups and equal group proportions α i = 1/3. With these settings, the theoretical expected number of edges is about 862 and the total number of potential edges is 19900. A sample graph is given in Figure 1. The running times of GLasso and SIMoNe are of the same order. For the settings described above the running time varies from a few seconds to a few minutes, according to the penalty parameter.
We focus the experiments on the ability to recover existing edges of the network, that is the nonzero entries of the concentration matrix. This is a binary decision problem where the compared algorithms are considered as classifiers. The decision made by a binary classifier can be summarized using four numbers: True Positives (T P ), False Positive (F P ), True Negatives (T N ) and False Negatives (F N ). We have chosen to draw precision/recall curves to display this information and compare how well the methods perform ( Figure 2).
Precision (T P/(T P + F P )) is the ratio of the number of true nonzero elements to the total number of nonzero elements in the estimated concentration matrix K. Recall that (T P/(T P + F N )) is the ratio of true nonzero elements in K to all nonzero entries of the real concentration matrix K. In a sparse context where the number of actual positives (T P + F N ) is small compared to the number of actual negatives (F P + T N ), precision/recall curves give a more informative picture of an algorithm's performance than classical Receiver Operator Characteristic (ROC) curves. Indeed, ROC curves plot the False Positive Rate (F P R = F P/(F P + T N )) against the True Positive Rate (T P R = T P/(T P + F N )). When the number of total positives is small compared to the number of total negatives, small variations of F P and T P will result in small variations of F P R and large variations of T P R, which is not relevant for comparing performances. In a statistical framework, the recall is equivalent to the power and the precision is equivalent to one minus the False Discovery Proportion.
Additionally to the GLasso (Friedman et al. 2008) and GeneNet (Schäfer and Strimmer 2005) we consider two other procedures: • When n is greater than p, a straightforward way to obtain an estimate of the inverse covariance matrix is to invert the empirical covariance matrix. Although this approach is unlikely to perform well in a selection context (since it is designed for estimation purposes), it is worth comparing it to its competitors in order to assess the scale of improvement. We call this procedure InvCor. • When the latent structure Z of the concentration matrix is known, our method can be applied without its E-step and produce a relevant selection of the nonzero entries of the concentration matrix. This approach represents the upper limit of our method, since it makes use of an usually unavailable source of information. This procedure is denoted perfect SIMoNe.
In some problems the latent structure of the graph is partially known and this information can be used in the E-step to improve the estimation of the latent structure. For example, when inferring gene regulation networks, a subset of identified genes may be known to belong to the same functional module.
The approach of Meinshausen and Bühlmann (2006) was also tested. The principle of this approach, and the performances obtained are close to those of GLasso, but it was always slightly outperformed. We have therefore decided, for the sake of brevity, to report only the four previously described procedures.
For the methods based on penalization (GLasso, SIMoNe and Perfect SIMoNe), the precision/recall curves are plotted by varying the penalty parameter (namely 1/λ in in our case). The penalty parameter varies from close to zero to a maximum value which forces all off-diagonal elements of K to be null (see Appendix A.6). The GeneNet and InvCor methods are plotted by sorting the elements of K according to their absolute values, and choosing different thresholds to find nonzero entries.
Even when n is really greater than p (Figures 2 (a-b)) Invcor is always dominated by the other methods from a selection point of view. This simple check shows that even in a favorable context with abundant data, penalization procedures improve the selection of nonzero entries of the concentration matrix, in comparison with methods based on estimation of these entries.
Although GeneNet and GLasso can provide different results on a given run, both methods perform similarly on average (50 runs for our experiment). The only parameter we change in this experimental setting is the n/p ratio.
Perfect SIMoNe's curves dominate all other curves for any n/p ratio. This clearly shows that the knowledge of the structure provides a valuable information for selecting the nonzero entries of the concentration matrix. When the structure is hidden, the main problem of our approach is then to find a reliable estimate of this structure from the initial data.
Perfect SIMoNe and SIMoNe perform equivalently when n = 10p and when the ratio n/p decreases, Perfect SIMoNe tends to outperform SIMoNe more clearly. This means that SIMoNe is able to recover the latent structure when there is enough data, but does not find a substantial structure when n drops below p.
When p > n, the empirical covariance matrix ceases to be invertible. Thus, Figures 2 (e-f) do not display the InvCor results. Although it is possible to show that both GLasso and SIMoNe increase the number of inferred true nonzero elements with the number of iterations in all settings, precision/recall curves show the relative poor performances for all tested algorithms when p ≥ n.
Notice that when p > n, the estimated latent structure is not reliable. Nevertheless, the performance of SIMoNe remains comparable to that of GLasso. We can therefore see that assuming the existence of a latent structure when there is none does not impair the selection of nonzero entries of the matrix K.

Breast Cancer data
We tested our algorithm on a gene expression data set provided by Hess et al. (2006) and concerning 133 patients with stage I − III breast cancer. The patients were treated with chemotherapy prior to surgery. Patient response to the treatment is classified as either a pathologic complete response (pCR) or a residual disease (not-pCR). Hess et al. (2006) and Natowicz et al. (2008) developed and tested a multigene predictor for treatment response on this data set. They focused on a set of 26 genes having a high predictive value (see Table 1). We thus consider a total of n = 133 cases containing p = 26 gene expression levels.
When dealing with gene regulatory networks, we typically observe n independent microarray experiments, each giving the expression levels of the same p genes. If the same experimental conditions are used for all microarrays, these may be considered as a sample of the same experiment. In the application in question, cases from the pCR class (34 cases) and from the not-pCR class (99 cases) clearly do not have the same distribution. We apply our algorithm on each class of patients. Two distinct gene regulatory networks are thus inferred. Figure 3 plots the resulting networks obtained for three different penalizations. The penalization parameters were heuristically chosen from the number of expected nonzero entries. We used Q = 2 latent clusters, and it is interesting to note that when assuming more than two clusters, the algorithm systematically produces exactly two non-empty clusters.
The inferred networks exhibit very different structures according to the class of patients. This in itself is interesting and suggests that gene regulation differs with respect to the presence or absence of a pCR.
The network obtained with not-pCR cases displays a two-star pattern. Each star connects to a unique gene, either SCUBE2 or IGFBP4. Almost all the most significant connections imply SCUBE2. This star pattern suggests that further studies of this particular gene would be of interest for understanding residual disease.
The network estimated with the pCR cases has a different two-cluster structure. In particular, it groups IGFBP4 and SCUBE2 in the same cluster with a direct significant link. This again indicates a completely different relationship between the genes in pCR versus non-pCR. Insulin-like growth factor binding protein 4 Table 1 The key genes that composed the inferred networks.

A.2. Fixed-point study
Let us first introduce some notation. For any i, j ∈ P and any q, ∈ Q, consider the random variables Let u : R pQ → R pQ be defined by its coordinate functions u = (u iq ) i∈P,q∈Q in the following way and let g = (g iq ) i∈P,q∈Q : R pQ → R pQ satisfy ∀a ∈ R pQ , g iq (a) = u iq (a) u i (a) .
According to Proposition 3, the optimal parameter τ is a fixed-point of g. Now, let Θ = a = (a iq ) i∈P,q∈Q ∈ R pQ ; ∀i ∈ P, q ∈ Q, a iq ∈ [0, 1] and q a iq = 1 .
We wish to study the fixed-points of g in Θ. First, let us note that as Θ is a compact state space and as the function g satisfies g : Θ → Θ and is continuous, the existence of a fixed-point of g follows from Brouwer's Theorem. We now restrict our attention to a smaller set than the whole state space Θ. For any ε > 0, let Note that we do not claim that g : Θ ε → Θ ε . However, the existence of a fixedpoint of g is ensured in Θ and if we assume α q > 0 for any q ∈ Q (which is a reasonable assumption if the number of classes Q is not too large), it can easily be seen that any fixed-point satisfies a iq > 0, for any i ∈ P and any q ∈ Q. Thus for sufficiently small ε > 0, the fixed-points of g belong to Θ ε .
In order to study the behaviour of g in the vicinity of a fixed-point, we need to look at some kind of contraction property for g. To this end we introduce a distance d on Θ ε that will make use of the form of the state space Θ ε . For all a, b ∈ Θ ε , denote by a i = (a iq ) q∈Q ∈ R Q and b It is well known that d 0 is a distance in [ε, 1 − ε] Q , and it is easy to check that the resulting d is also a distance in Θ ε . Now, fix a, b ∈ R pQ and consider the distance d(g(a), g(b)). It is easily checked that d(g(a), g(b)) = max i∈P d 0 (g i (a), g i (b)) = max i∈P d 0 (u i (a), u i (b)) = max i∈P d 0 (ū i (a),ū i (b)), whereū = (ū i ) i∈P = (ū iq ) i∈P,q∈Q is defined in the following way In the following, fix ε > 0 and a, b ∈ Θ ε and denote by With these notations, we have We only consider the affiliation model described in (5). Thus, there are only two different values for λ q , namely λ in and λ out for intra and extra cluster connectivity.
Before proving the lemma, let us explain the consequences of this result. Consider the function h K defined on (0, +∞) by This function first decreases from +∞ to the value 1 + log 2|K| on the interval (0, |K|) and then increases from 1 + log 2|K| to +∞ on (|K|, +∞). At any step of the algorithm, if the current values K (m) ij of the concentration matrix are small enough, namely smaller than 1/(2e) 0.184 then the functions h K (m) ij take all the values between 1 + log 2|K| < 0 and +∞. Thus, there is room for choosing λ in , λ out such that (34) is satisfied. In such a case, the fixed-point we are looking for is unique and the iterative procedure setting τ (s+1) = g( τ (s) ) converges.
Proof. Using that for any j ∈ P and any ∈ Q, we have c j 1 b j ≤ a j ≤ c j 2 b j and L ijq > 0, we get (35) In the case of the affiliation model, for fixed i, j ∈ P and q ∈ Q, the set of random variables {L ijq } ∈Q is reduced to only two random values, namely For the sake of simplicity, we assume Q = 2 groups (our arguments may be easily generalized to 3 groups or more). Now, denoting L max ij = max(L in ij , L out ij ) and L min ij = min(L in ij , L out ij ), it can easily be seen that (for ε < 1/2), almost surely. Note that if we have Q ≥ 3 groups, explicit bounds can also be obtained (their expression is only slightly more complicated). Coming back to (35), we get This leads to d 0 (ū i (a),ū i (b)) = log max q∈Qūiq (a)/ū iq (b) min q∈Qūiq (a)/ū iq (b) Finally, recall that d(g(a), g(b)) = max i d 0 (ū i (a),ū i (b)), leading to Now, using the inverse triangle inequality, and the fact that c i 1 ≤ 1 ≤ c i 2 , we get for any i ∈ P, Since a and b belong to Θ ε , we get that c i 1 , c i 2 ∈ [ε, ε −1 ] and thus In particular, recalling (33), we have Coming back to (36), we get d(g(a), g(b)) ≤ (1 + ε −1 )2(p − 1) max Now, under assumption (34) the multiplicative random factor (1 + ε −1 )2(p − 1) max j =i L max ij is strictly smaller than 1.
A.3. Proof of Lemma 2 ( Lasso with pathwise coordinate optimization) The following is partly based on Friedman et al. (2007). There are various algorithms for solving the Lasso problem. When there is just one predictor, the Lasso solution is simply given by soft-thresholding (Donoho and Johnstone 1995). The approach used here is based on iterative soft-thresholding with a "partial residual" as a response variable.
The usual formulation of the Lasso problem is the minimization with respect to β of the quantity where (y i ) i=1,...,n is a vector of response and (x ij ) i=1,...,n;j=1,...,p a matrix of predictors such that i x ij = 0, with no loss of generality. Using a coordinatedescent approach, we simply write the problem (38) in the form and minimizing this function with respect to β j will lead to the solution whereỹ (j) i = k =j x ik β k (ρ), the normalizing term N 2 j satisfies N 2 j = n i=1 x 2 ij and the function S(x, ρ) = sgn(x)(|x| − ρ) + is the soft-thresholding operator. This leads to an iterative procedure, repeated on each coordinate of β until stabilization of the full vector. Note that as each coordinate-wise solution is unique, results from Tseng (2001, Theorem 4.1) imply that the procedure converges. Now, we want to apply this approach to solve the problem (21), which can be written From the previous lines, the solution for jth entry of β is Then, using the symmetry of the matrices, it is easy to see that