{'fraction_non_alphanumeric': 0.05026882432396642, 'fraction_numerical': 0.017903964827714065, 'mean_word_length': 4.043815121560246, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fused Multiple Graphical Lasso

In this paper, we consider the problem of estimating multiple graphical models simultaneously using the fused lasso penalty, which encourages adjacent graphs to share similar structures. A motivating example is the analysis of brain networks of Alzheimer's disease using neuroimaging data. Specifically, we may wish to estimate a brain network for the normal controls (NC), a brain network for the patients with mild cognitive impairment (MCI), and a brain network for Alzheimer's patients (AD). We expect the two brain networks for NC and MCI to share common structures but not to be identical to each other; similarly for the two brain networks for MCI and AD. The proposed formulation can be solved using a second-order method. Our key technical contribution is to establish the necessary and sufficient condition for the graphs to be decomposable. Based on this key property, a simple screening rule is presented, which decomposes the large graphs into small subgraphs and allows an efficient estimation of multiple independent (small) subgraphs, dramatically reducing the computational cost. We perform experiments on both synthetic and real data; our results demonstrate the effectiveness and efficiency of the proposed approach.


Introduction
Undirected graphical models explore the relationships among a set of random variables through their joint distribution. The estimation of undirected graphical models has applications in many domains, such as computer vision, biology, and medicine [11,17,44]. One instance is the analysis of gene expression data. As shown in many biological studies, genes tend to work in groups based on their biological functions, and there exist some regulatory relationships between genes [5]. Such biological knowledge can be represented as a graph, where nodes are the genes, and edges describe the regulatory relationships. Graphical models provide a useful tool for modeling these relationships, and can be used to explore gene activities. One of the most widely used graphical models is the Gaussian graphical model (GGM), which assumes the variables to be Gaussian distributed [2,47]. In the framework of GGM, the problem of learning a graph is equivalent to estimating the inverse of the covariance matrix (precision matrix), since the nonzero off-diagonal elements of the precision matrix represent edges in the graph [2,47].
In recent years many research efforts have focused on estimating the precision matrix and the corresponding graphical model (see, for example [2,10,16,17,22,23,25,26,29,30,33,47]. Meinshausen and Bühlmann [30] estimated edges for each node in the graph by fitting a lasso problem [36] using the remaining variables as predictors. Yuan and Lin [47] and Banerjee et al. [2] proposed a penalized maximum likelihood model using ℓ 1 regularization to estimate the sparse precision matrix. Numerous methods have been developed for solving this model. For example, d'Aspremont et al. [8] and Lu [25,26] studied Nesterov's smooth gradient methods [32] for solving this problem or its dual. Banerjee et al. [2] and Friedman et al. [10] proposed block coordinate ascent methods for solving the dual problem. The latter method [10] is widely referred to as Graphical lasso (GLasso). Mazumder and Hastie [29] proposed a new algorithm called DP-GLasso, each step of which is a box-constrained QP problem. Scheinberg and Rish [35] proposed a coordinate descent method for solving this model in a greedy approach. Yuan [48] and Scheinberg et al. [34] applied alternating direction method of multipliers (ADMM) [4] to this problem. Li and Toh [22] and Yuan and Lin [47] proposed to solve this problem using interior point methods. Wang et al. [40], Hsieh et al. [16], Olsen et al. [33], and Dinh et al. [9] studied Newton method for solving this model. The main challenge of estimating a sparse precision matrix for the problems with a large number of nodes (variables) is its intensive computation. Witten et al. [42] and Mazumder and Hastie [28] independently derived a necessary and sufficient condition for the solution of a single graphical lasso to be block diagonal (subject to some rearrangement of variables). This can be used as a simple screening test to identify the associated blocks, and the original problem can thus be decomposed into a group of smaller sized but independent problems corresponding to these blocks. When the number of blocks is large, it can achieve massive computational gain. However, these formulations assume that observations are independently drawn from a single Gaussian distribution. In many applications the observations may be drawn from multiple Gaussian distributions; in this case, multiple graphical models need to be estimated.
There are some recent works on the estimation of multiple precision matrices [7,11,12,13,19,20,31,49]. Guo et al. [11] proposed a method to jointly estimate multiple graphical models using a hierarchical penalty. However, their model is not convex. Honorio and Samaras [13] proposed a convex formulation to estimate multiple graphical models using the ℓ 1,∞ regularizer. Hara and Washio [12] introduced a method to learn common substructures among multiple graphical models. Danaher et al. [7] estimated multiple precision matrices simultaneously using a pairwise fused penalty and grouping penalty. ADMM was used to solve the problem, but it requires computing multiple eigen decompositions at each iteration. Mohan et al. [31] proposed to estimate multiple precision matrices based on the assumption that the network differences are generated from node perturbations. Compared with single graphical model learning, learning multiple precision matrices jointly is even more challenging to solve. Recently, a necessary and sufficient condition for multiple graphs to be decomposable was proposed in [7]. However, such necessary and sufficient condition was restricted to two graphs only when the fused penalty is used. It is not clear whether this screening rule can be extended to the more general case with more than two graphs, which is the case in brain network modeling.
There are two types of fused penalties that can be used for estimating multiple (more than two) graphs: (a) pairwise fused or (b) sequential fused [37]. In this paper we set out to address the sequential fused case first, because we work on practical applications that can be more appropriately formulated using the sequential formulation. Specifically, we consider the problem of estimating multiple graphical models by maximizing a penalized log likelihood with ℓ 1 and sequential fused regularization. The ℓ 1 regularization yields a sparse solution, and the fused regularization encourages adjacent graphs to be similar. The graphs considered in this paper have a natural order, which is common in many applications. A motivating example is the modeling of brain networks for Alzheimer's disease using neuroimaging data such as Positron emission tomography (PET). In this case, we want to estimate graphical models for three groups: normal controls (NC), patients of mild cognitive impairment (MCI), and Alzheimer's patients (AD). These networks are expected to share some common connections, but they are not identical. Furthermore, the networks are expected to evolve over time, in the order of disease progression from NC to MCI to AD. Estimating the graphical models separately fails to exploit the common structures among them. It is thus desirable to jointly estimate the three networks (graphs).
Our key technical contribution is to establish the necessary and sufficient condition for the solution of the fused multiple graphical lasso (FMGL) to be block diagonal. The duality theory and several other tools in linear programming are used to drive the necessary and sufficient condition. Based on this crucial property of FMGL, we develop a screening rule which enables the efficient estimation of large multiple precision matrices for FMGL. The proposed screening rule can be combined with any algorithms to reduce computational cost. We employ a second-order method [16,21,38] to solve the fused multiple graphical lasso, where each step is solved by the spectral projected gradient method [27,43]. In addition, we propose a shrinking scheme to identify the variables to be updated in each step of the second-order method, which reduces the computation cost of each step. We conduct experiments on both synthetic and real data; our results demonstrate the effectiveness and efficiency of the proposed approach.

Notation
In this paper, ℜ stands for the set of all real numbers, ℜ n denotes the n-dimensional Euclidean space, and the set of all m × n matrices with real entries is denoted by ℜ m×n . All matrices are presented in bold format. The space of symmetric matrices is denoted by S n . If X ∈ S n is positive semidefinite (resp. definite), we write X 0 (resp. X ≻ 0). Also, we write X Y to mean X − Y 0. The cone of positive semidefinite matrices in S n is denoted by S n + . Given matrices X and Y in ℜ m×n , the standard inner product is defined by X, Y := tr(XY T ), where tr(·) denotes the trace of a matrix. X • Y and X ⊗ Y means the Hadamard and Kronecker product of X and Y, respectively. We denote the identity matrix by I, whose dimension should be clear from the context. The determinant and the minimal eigenvalue of a real symmetric matrix X are denoted by det(X) and λ min (X), respectively. Given a matrix X ∈ ℜ n×n , diag(X) denotes the vector formed by the diagonal of X, that is, diag(X) i = X ii for i = 1, . . . , n. Diag(X) is the diagonal matrix which shares the same diagonal as X. vec(X) is the vectorization of X. In addition, X > 0 means that all entries of X are positive.
The rest of the paper is organized as follows. We introduce the fused multiple graphical lasso formulation in Section 2. The screening rule is presented in Section 3. The proposed second-order method is presented in Section 4. The experimental results are shown in Section 5. We conclude the paper in Section 6.

Fused multiple graphical lasso
Assume we are given K data sets, x (k) ∈ R n k ×p , k = 1, . . . , K with K ≥ 2, where n k is the number of samples, and p is the number of features. The p features are common for all K data sets, and all K k=1 n k samples are independent. Furthermore, the samples within each data set x (k) are identically distributed with a p-variate Gaussian distribution with zero mean and positive definite covariance matrix Σ (k) , and there are many conditionally independent pairs of features, i.e., the precision matrix Θ (k) = (Σ (k) ) −1 should be sparse. For notational simplicity, we assume that n 1 = · · · = n K = n. Denote the sample covariance matrix for each data set x (k) as S (k) with S (k) = 1 n (x (k) ) T x (k) , and Θ = (Θ (1) , . . . , Θ (K) ). Then the negative log likelihood for the data takes the form of K k=1 − log det(Θ (k) ) + tr(S (k) Θ (k) ) . (1) Clearly, minimizing (1) leads to the maximum likelihood estimate (MLE)Θ (k) = (S (k) ) −1 . However, the MLE fails when S (k) is singular. Furthermore, the MLE is usually dense. The ℓ 1 regularization has been employed to induce sparsity, resulting in the sparse inverse covariance estimation [2,10,46]. In this paper, we employ both the ℓ 1 regularization and the fused regularization for simultaneously estimating multiple graphs. The ℓ 1 regularization leads to a sparse solution, and the fused penalty encourages Θ (k) to be similar to its neighbors. Mathematically, we solve the following formulation: where λ 1 > 0 and λ 2 > 0 are positive regularization parameters. This model is referred to as the fused multiple graphical lasso (FMGL).
To ensure the existence of a solution for problem (2), we assume throughout this paper that diag(S (k) ) > 0, k = 1, . . . , K. Recall that S (k) is a sample covariance matrix, and hence diag(S (k) ) ≥ 0. The diagonal entries may be not, however, strictly positive. But we can always add a small perturbation (say 10 −8 ) to ensure the above assumption holds.
The following theorem shows that under this assumption the FMGL (2) has a unique solution.
To prove Theorem 2.1, we first establish a technical lemma which regards the existence of a solution for a standard graphical lasso problem.
Lemma 2.2. Let S ∈ S p + and Λ ∈ S p be such that Diag(S) + Λ > 0 and diag(Λ) ≥ 0. Consider the problem min Then the following statements hold: (a) Problem (3) has a unique optimal solution; We first claim that the feasible region of problem (4) is nonempty, or equivalently, there exists U ∈ U such that λ min (S + Λ •Ū) > 0. Indeed, one can observe that where the second equality follows from the Lagrangian duality since its associated Slater condition is satisfied. Let Ω := {X ∈ S p : tr(X) = 1, X 0}. By the assumption Diag(S) + Λ > 0, we see that Λ ij > 0 for all i = j and S ii + Λ ii > 0 for every i. Since Ω ⊂ S p + , we have tr(SX) ≥ 0 for all X ∈ Ω. If there exists some k = l such that X kl > 0, then i =j Λ ij |X ij | > 0 and hence, Otherwise, one has X ij = 0 for all i = j, which, together with the facts that S ii + Λ ii > 0 for all i and tr(X) = 1, implies that for all X ∈ Ω, Hence, (6) again holds. Combining (5) with (6), one can see that max U∈U λ min (S+Λ•U) > 0. Therefore, problem (4) has at least a feasible solution. We next show that problem (4) has an optimal solution. LetŪ be a feasible point of (4), and One can observe that {S + Λ • U : U ∈ U} is compact. Using this fact, it is not hard to see that log det(S + Λ • U) → −∞ as U ∈ U and λ min (S + Λ • U) ↓ 0. Thus there exists some δ > 0 such that which implies that Hence,Ω is a compact set. In addition, one can observe that problem (4) is equivalent to The latter problem clearly has an optimal solution and so is problem (4). Finally we show that X * = (S + Λ • U * ) −1 is the unique optimal solution of (3), where U * is an optimal solution of (4). Since S + Λ • U * ≻ 0, we have X * ≻ 0. By the definitions of U and X * , and the first-order optimality conditions of (4) at U * , one can have where ∂(·) stands for the subdifferential of the associated convex function. For convenience, let f (X) denote the objective function of (3). Then we have which, together with X * = (S + Λ • U * ) −1 , implies that 0 ∈ ∂f (X * ). Hence, X * is an optimal solution of (3) and moreover it is unique due to the strict convexity of − log det(·). (b) By statement (a), problem (3) has a finite optimal value f * . Hence, the above sub-level set L is nonempty. We can observe that for any X ∈ L, where f * := inf{f (X) : X ≻ 0}. By the assumption Diag(S) + Λ > 0, one has Diag(S) + Λ/2 > 0. This together with statement (a) yields f * ∈ ℜ. Notice that Λ ij > 0 for all i = j. This relation and (7) imply that X ij is bounded for all X ∈ L and i = j. In addition, it is well-known that det(X) ≤ X 11 X 22 · · · X pp for all X 0. Using this relation, the definition of f (·), and the boundedness of X ij for all X ∈ L and i = j, we have that for every X ∈ L, for some δ > 0. In addition, notice from the assumption that S ii + Λ ii > 0 for all i, and hence for all i. This relation together with (8) implies that for every X ∈ L and all i, and hence X ii is bounded for all i and X ∈ L. We thus conclude that L is bounded. In view of this result and the definition of f , it is not hard to see that there exists some ν > 0 such that λ min (X) ≥ ν for all X ∈ L. Hence, one has By the continuity of f on {X : X νI}, it follows that L is closed. Hence, L is compact.
We are now ready to prove Theorem 2.1.
In view of Lemma 2.2, we know that Ω k is compact for all k, which implies thatΩ is also compact. Notice that h is continuous and strictly convex onΩ. Hence, problem (9) has a unique optimal solution and so is problem (2).

The screening rule for fused multiple graphical lasso
Due to the presence of the log determinant, it is challenging to solve the formulations involving the penalized log-likelihood efficiently. The existing methods for single graphical lasso are not scalable to the problems with a large amount of features because of the high computational complexity. Recent studies have shown that the graphical model may contain many connected components, which are disjoint with each other, due to the sparsity of the graphical model, i.e., the corresponding precision matrix has a block diagonal structure (subject to some rearrangement of features). To reduce the computational complexity, it is advantageous to first identify the block structure and then compute the diagonal blocks of the precision matrix instead of the whole matrix. Danaher et al. [7] developed a similar necessary and sufficient condition for fused graphical lasso with two graphs, thus the block structure can be identified. However, it remains a challenge to derive the necessary and sufficient condition for the solution of fused multiple graphical lasso to be block diagonal for K > 2 graphs. In this section, we first present a theorem demonstrating that FMGL can be decomposable once its solution has a block diagonal structure. Then we derive a necessary and sufficient condition for the solution of FMGL to be block diagonal for arbitrary number of graphs.
Let C 1 , . . . , C L be a partition of the p features into L non-overlapping sets, with C l ∩C l ′ = ∅, ∀l = l ′ and L l=1 C l = {1, . . . , p}. We say that the solution Θ of FMGL (2) is block diagonal with L known blocks consisting of features in the sets C l , l = 1, . . . , L if there exists a permutation matrix U ∈ ℜ p×p such that each estimation precision matrix takes the form of For simplicity of presentation, we assume throughout this paper that U = I.
The following decomposition result for problem (2) is straightforward. Its proof is thus omitted.
where Θ (k) l and S (k) l are the |C l | × |C l | symmetric submatrices of Θ (k) and S (k) corresponding to the l-th diagonal block, respectively, for k = 1, . . . , K, and Θ l = (Θ The above theorem demonstrates that if a large-scale FMGL problem has a block diagonal solution, it can then be decomposed into a group of smaller sized FMGL problems. The computational cost for the latter problems can be much cheaper. Now one natural question is how to efficiently identify the block diagonal structure of the FMGL solution before solving the problem. We address this question in the remaining part of this section.
The following theorem provides a necessary and sufficient condition for the solution of FMGL to be block diagonal with L blocks C l , l = 1, . . . , L, which is a key for developing efficient decomposition scheme for solving FMGL. Since its proof requires some substantial development of other technical results, we shall postpone the proof until the end of this section.
Theorem 3.2. The FMGL (2) has a block diagonal solution Θ (k) , k = 1, . . . , K with L known blocks C l , l = 1, . . . , L if and only if S (k) , k = 1, . . . , K satisfy the following inequalities: One immediate consequence of Theorem 3.2 is that the conditions (12) can be used as a screening rule to identify the block diagonal structure of the FMGL solution. The steps about this rule are described as follows.
2. Identify the connected components of the adjacency matrix E (for example, it can be done by calling the Matlab function "graphconncomp").
In view of Theorem 3.2, it is not hard to observe that the resulting connected components are the partition of the p features into nonoverlapping sets. It then follows from Theorem 3.1 that a large-scale FMGL problem can be decomposed into a group of smaller sized FMGL problems restricted to the features in each connected component. The computational cost for the latter problems can be much cheaper. Therefore, this approach may enable us to solve large-scale FMGL problems very efficiently.
In the remainder of this section we provide a proof for Theorem 3.2. Before proceeding, we establish several technical lemmas as follows.
(ii) It is easy to see that 0 ∈ P I,J and moreover there exist n linearly independent active inequalities at 0. Hence, 0 is an extreme point of P I,J . On the other hand, suppose y is an arbitrary extreme point of P I,J . Then there exist n linearly independent active inequalities at y, which together with the definition of P I,J immediately implies y = 0. Therefore, 0 is the unique extreme point of P I,J .
(iii) Suppose that P I,J is unbounded. By statement (ii), we know that P I,J has a unique extreme point. Using Minkowski's resolution theorem (e.g., see [3]), we conclude that ext(P I,J ) = ∅. Let d ∈ ext(P I,J ) be arbitrarily chosen. Then d = 0. It follows from (13) that d satisfies the inequalities and moreover, the number of independent active inequalities at d is n−1. If all entries of d are nonzero, then d must satisfy d J − d J+1 = 0 and dJ − dJ +1 = 0 (with a total number n − 1), which implies d 1 = d 2 = · · · = d n and thus d ∈ Q. We now assume that d has at least one zero entry. Then, there exist positive integers k, One can immediately observe that We next divide the rest of proof into four cases.
Case (a): m 1 = 1 and n k = n. In view of (16), one can observe that d mi−1 − d mi = 0 and d ni−1 − d ni−1+1 = 0 for i = 2, . . . , k. We then see from (15) that except the active inequalities given in (17), all other possible active inequalities at d are (with a total number k i=2 (m i − n i−1 − 2)). Notice that the total number of independent active inequalities given in (17) is k i=1 (n i − m i + 1). Hence, the number of independent active inequalities at d is at most Recall that the number of independent active inequalities at d is n−1. Hence, we have n−k+1 ≥ n−1, which implies k ≤ 2. Due to d = 0, we observe that k = 1 holds for this case. Also, we know that k > 0. Hence, k = 2. We then see that all possible active inequalities described in (18) must be active at d, which together with k = 2 immediately implies that d ∈ Q.
Case (b): m 1 = 1 and n k < n. Using (16), we observe that d mi−1 − d mi = 0 for i = 2, . . . , k and d ni − d ni+1 = 0 for i = 1, . . . , k. In view of these relations and a similar argument as in case (a), one can see that the number of independent active inequalities at d is at most Similarly as in case (a), we can conclude from the above relation that k = 1 and d ∈ Q.
Case (c): m 1 > 1 and n k = n. By (16), one can observe that d mi−1 − d mi = 0 for i = 1, . . . , k and d ni − d ni+1 = 0 for i = 1, . . . , k − 1. Using these relations and a similar argument as in case (a), we see that the number of independent active inequalities at d is at most Similarly as in case (a), we can conclude from the above relation that k = 1 and d ∈ Q.
Case (d): m 1 > 1 and n k < n. From (16), one can observe that d mi−1 − d mi = 0 for i = 1, . . . , k and d ni − d ni+1 = 0 for i = 1, . . . , k. By virtue of these relations and a similar argument as in case (a), one can see that the number of independent active inequalities at d is at most Recall that k ≥ 1 and the number of independent active inequalities at d is n − 1. Hence, this case cannot occur.
Lemma 3.5. Let x ∈ ℜ n , λ 1 , λ 2 ≥ 0 be given, and let Then, f (y) ≤ 0 for all y ∈ ℜ n if and only if x satisfies the following inequalities: Proof. Let P I,J be defined in (13)  If P I,J is bounded, we have P I,J = {0} and f (y) = 0 for y ∈ P I,J . Suppose that P I,J is unbounded. By Lemma 3.3 and Minkowski's resolution theorem, P I,J equals the finitely generated cone by ext(P I,J ).
It then follows that f (y) ≤ 0 for all y ∈ P I,J if and only if f (d) ≤ 0 for all d ∈ ext(P I,J ). Using these facts and Lemma 3.4, we see that f (y) ≤ 0 for all y ∈ ℜ n if and only if f (d) ≤ 0 for all d ∈ Q, where Q is defined in (14). By the definitions of Q and f , we further observe that f (y) ≤ 0 for all y ∈ ℜ n if and only if f (d) ≤ 0 for all which together with the definition of f immediately implies that the conclusion of this lemma holds.

Proof. The linear system (19) has a solution if and only if the linear programming
has an optimal solution. The Lagrangian dual of (20) is By the Lagrangian duality theory, problem (20) has an optimal solution if and only if its dual problem (21) has optimal value 0, which is equivalent to f (y) ≤ 0 for all y ∈ ℜ n . The conclusion of this lemma then immediately follows from Lemma 3.5.
We are now ready to prove Theorem 3.2.

Second-order method
The screening rule proposed in Section 3 is capable of partitioning all features into a group of smaller sized blocks. Accordingly, a large-scale FMGL (2) can be decomposed into a number of smaller sized FMGL problems. For each block l, we need to compute its individual estimated precision matrix Θ l . In this section, we discuss how to solve those single block FMGL problems efficiently. For simplicity of presentation, we assume throughout this section that the FMGL (2) has only one block, that is, L = 1.
) be the approximate solution obtained at the t-th Newton iteration.
In the second-order method, we approximate the objective function F (Θ) at the current iterate Θ t by a "quadratic" model Q t (Θ): where q k is the quadratic approximation of f k at Θ (k) t , that is, t . Suppose thatΘ t+1 is the optimal solution of (28). Then we obtain the Newton search direction We shall mention that the subproblem (28) can be suitably solved by the non-monotone spectral projected gradient (NSPG) method (see, for example, [43,27]). It was shown by Lu and Zhang [27] that the NSPG method is locally linearly convergent. Numerous computational studies have demonstrated that the NSPG method is very efficient though its global convergence rate is so far unknown. When applied to (28), the NSPG method requires solving the proximal subproblems in the form of for some G = (G (1) , . . . , G (K) ) and α > 0. By the definition of P (Θ), it is not hard to see that problem (30) can be decomposed into a set of independent and smaller sized problems for all i ≥ j, j = 1, . . . , p, where (α 1 , α 2 ) = α(λ 1 , λ 2 ). The problem (31) is known as the fused lasso signal approximator, which can be solved very efficiently and exactly [6,24]. In addition, they are independent from each other and thus can be solved in parallel. Given the current search direction D = (D (1) , . . . , D (K) ) that is computed above, we need to find the suitable step length β ∈ (0, 1] to ensure a sufficient reduction in the objective function of (2) and positive definiteness of the next iterate Θ (k) t+1 = Θ (k) t + βD (k) , k = 1, . . . , K. In the context of the standard (single) graphical lasso, Hsieh et al. [16] have shown that a step length satisfying the above requirements always exists. We can similarly prove that the desired step length also exists for the FMGL (2).

Lemma 4.2.
For Θ t in the t-th iteration, define the fixed set J f ixed as t,ij satisfy the inequalities below}.
Then, the solution of the following optimization problem is D (1) = · · · = D (K) = 0 : Proof. Consider problem (34), which can be reformulated to where H we only consider the variables in the set J f ixed . According to Lemma 3.6, it is easy to see that D J f ixed = 0 satisfies the optimality condition of the following problem t vec(D (k) ) ≥ 0, the optimal solution of (34) is given by D (1) = · · · = D (K) = 0. Lemma 4.2 provides a shrinking scheme to partition the variables into the free set J f ree and the fixed set J f ixed . With shrinking, each Newton step of the proposed second-order method falls into a block coordinate gradient descent framework [38]. Lemma 4.2 shows that when the variables in the free set J f ree are fixed, no update is needed for the variables in the fixed set J f ixed . Minimization of (28) restricted to the free set can therefore guarantee the convergence to the unique optimal solution [16,38]. In addition, it has been shown that local quadratic convergence rate can be achieved when the exact Hessian is used (see, for example, [16,21]).
The resulting second-order method for solving the fused multiple graphical lasso is summarized in Algorithm 1.

Experimental results
In this section, we evaluate the proposed algorithm and screening rule on synthetic datasets and two real datasets: ADHD-200 2 and FDG-PET images 3 . The experiments are performed on a PC with quad-core Intel 2.67GHz CPU and 9GB memory.

Efficiency
We conduct experiments to demonstrate the effectiveness of the proposed screening rule and the efficiency of our method FMGL. The following algorithms are included in our comparisons: • FMGL: the proposed second-order method in Algorithm 1.
Both FMGL and ADMM are written in Matlab. Since both methods involve solving (30) which involves a double loop, we implement the sub-routine for solving (30) in C for a fair comparison.
The synthetic covariance matrices are generated as follows. We first generate K block diagonal ground truth precision matrices Θ (k) with L blocks, and each block Θ to be about 10p/L so that the total number of nonzeros in the K precision matrices is 10Kp. Given the precision matrices, we draw 5p samples from each Gaussian distribution to compute the sample covariance matrices. The fused penalty parameter λ 2 is fixed to 0.1, and the ℓ 1 regularization parameter λ 1 is selected so that the total number of nonzeros in the solution is about 10Kp. We terminate the NSPG in the FMGL when the relative error FMGL is terminated when the relative error of the objective value is smaller than 1e-5, and ADMM stops until it achieves an objective value equal to or smaller than that of FMGL. The results presented in Table 1 show that FMGL is consistently faster than ADMM. FMGL converges much more quickly than ADMM. Moreover, the screening rule can achieve great computational gain. The speedup with the screening rule is about 10 and 20 times for L = 5 and 10 respectively.

Stability
We conduct experiments to demonstrate the effectiveness of FMGL. The synthetic sparse precision matrices are generated in the following way: we set the first precision matrix Θ (1) as 0.25I p×p , where p = 100. When adding an edge (i, j) in the graph, we add σ to θ (1) ii and θ (1) jj , and subtract σ from θ (1) ij and θ (1) ji to keep the positive definiteness of Θ (1) , where σ is uniformly drawn from [0.1, 0.3]. When deleting an edge (i, j) from the graph, we reverse the above steps with σ = θ (1) ij . We randomly assign Table 1: Comparison of the proposed FMGL and ADMM with and without screening in terms of average computational time (seconds). FMGL-S and ADMM-S are FMGL and ADMM with screening respectively. p stands for the dimension, K is the number of graphs, L is the number of blocks, and λ 1 is the ℓ 1 regularization parameter. The fused penalty parameter λ 2 is fixed to 0.1. Θ 0 represents the total number of nonzero entries in ground truth precision matrices Θ (k) , k = 1, . . . , K, and Θ * 0 is the number of nonzeros in the solution.
Data and parameter setting Computational time (iteration numbers) 200 edges for Θ (1) . Θ (2) is obtained by adding 25 edges and deleting 25 different edges from Θ (1) . 2) in the same way. For each precision matrix, we randomly draw n samples from the Gaussian distribution with the corresponding precision matrix, where n varies from 40 to 200 with a step of 20. We perform 500 replications for each n. For each n, λ 2 is fixed to 0.08, and λ 1 is adjusted to make sure that the edge number is about 200. The accuracy n d /n g is used to measure the performance of FMGL and GLasso, where n d is the number of true edges detected by FGML and GLasso, and n g is the number of true edges. The results are shown in Figure 1. We can see from the figure that FMGL achieves higher accuracies, demonstrating the effectiveness of FMGL for learning multiple graphical models simultaneously.

ADHD-200
Attention Deficit Hyperactivity Disorder (ADHD) affects at least 5-10% of school-age children with annual costs exceeding 36 billion/year in the United States. The ADHD-200 project has released resting-state functional magnetic resonance images (fMRI) of 491 typically developing children and 285 ADHD children, aiming to encourage the research on ADHD. The data used in this experiment is preprocessed using the NIAK pipeline, and downloaded from neurobureau 4 . More details about the preprocessing strategy can be found in the same website. The dataset we choose includes 116 typically developing children (TDC), 29 ADHD-Combined (ADHD-C), and 49 ADHD-Inattentive (ADHD-I). There are 231 time series and 2834 brain regions for each subject. We want to estimate the graphs of the three groups simultaneously. The sample covariance matrix is computed using all data from the same group. Since the number of brain regions p is 2834, obtaining the precision matrices is computationally intensive. We use this data to test the effectiveness of the proposed screening rule. λ 1 and λ 2 are set to 0.6 and 0.015. The convergence criterion is 1e-5. The comparison of FMGL and ADMM in terms of the objective value curve is shown in Figure 2. The result shows that FMGL converges much faster than ADMM. The computational times of FMGL and ADMM are 1557.08 and 8306.35 seconds respectively. However, utilizing the screening, the computational times of FMGL-S and ADMM-S are 18.08 and 119.23 seconds respectively, demonstrating the superiority of the screening rule. The obtained solution has 1443 blocks. The largest one including 634 nodes is shown in Figure 3. The block structures of the FMGL solution are the same as those identified by the screening rule. The screening rule can be used to analyze the rough structures of the graphs. The cost of identifying blocks using the screening rule is negligible compared to that of estimating the graphs. For highdimensional data such as ADHD-200, it is practical to use the screening rule to identify the block structure before estimating the large graphs. We use the screening rule to identify block structures on ADHD-200 data with varying λ 1 and λ 2 . The size distribution is shown in Figure 4. We can observe that the number of blocks increases, and the size of blocks deceases when the regularization parameter value increases.

FDG-PET
In this experiment, we use FDG-PET images from 74 Alzhei-mer's disease (AD), 172 mild cognitive impairment (MCI), and 81 normal control (NC) subjects downloaded from the Alzheimer's disease neuroimaging initiative (ADNI) database. The different regions of the whole brain volume can be represented by 116 anatomical volumes of interest (AVOI), defined by Automated Anatomical Labeling (AAL) [39]. Then we extracted data from each of the 116 AVOIs, and derived the average of each AVOI for each subject. The 116 AVOIs can be categorized into 10 groups: prefrontal lobe, other parts of the frontal lobe, parietal lobe, occipital lobe, thalamus, insula, temporal lobe, corpus striatum, cerebellum, and vermis. More details about the categories can be found in [39,41]. We remove two small groups (thalamus and insula) containing only 4 AVOIs in our experiments.
To examine whether FMGL can effectively utilize the information of common structures, we randomly select g percent samples from each group, where g varies from 20 to 100 with a step size of 10. For each g, λ 2 is fixed to 0.1, and λ 1 is adjusted to make sure the number of edges in each group is about the same. We perform 500 replications for each g. The edges with probability larger than 0.85 are considered as stable edges. The results showing the numbers of stable edges are summarized in Figure 5. We can observe that FMGL is more stable than GLasso. When the sample size is too small (say 20%), there are only 20 stable edges in the graph of NC obtained by GLasso. But the graph of NC obtained by FMGL still has about 140 stable edges, illustrating the superiority of FMGL in stability.  The brain connectivity models obtained by FMGL are shown in Figure 6. We can see that the number of connections within the prefrontal lobe significantly increases, and the number of connections within the temporal lobe significantly decreases from NC to AD, which are supported by previous literatures [1,14]. The connections between the prefrontal and occipital lobes increase from NC to AD, and connections within cerebellum decrease. We can also find that the adjacent graphs are similar, indicating that FMGL can identify the common structures, but also keep the meaningful differences.

Conclusion
In this paper, we consider simultaneously estimating multiple graphical models by maximizing a fused penalized log likelihood. We have derived a set of necessary and sufficient conditions for the FMGL solution to be block diagonal for an arbitrary number of graphs. A screening rule has been developed to enable the efficient estimation of large multiple graphs. The second-order method is employed to solve the fused multiple graphical lasso. The global convergence of the proposed method is guaranteed, and the convergence rate is local quadratic. A shrinking scheme is proposed to identify the variables to be updated during the Newton iterations, thus reduces the computation. Numerical experiments on synthetic and real data demonstrate the efficiency and effectiveness of the proposed method and the screening rule. We plan to explore the convergence properties of the second-order method using the inexact Newton direction. Due to the shrinking scheme, the proposed second-order method is suitable for warm-start techniques. A good initial solution can further speedup the computation. As part of the future work, we plan to explore how to efficiently find a good initial solution to further improve the efficiency of the proposed method. One possibility is to use divide-and-conquer techniques [15].