Network classification with applications to brain connectomics
Jesús D. Arroyo-Relión Daniel Kessler Elizaveta Levina Stephan F. Taylor





Network classification


,




J.D. Arroyo Relión et al.

Johns Hopkins Universitym1 and University of Michiganm2

J. D. Arroyo Relión

Center for imaging science

Johns Hopkins University

Baltimore, Maryland 21218-2682

USA

e1

D. Kessler

E. Levina

Department of Statistics

University of Michigan

Ann Arbor, Michigan 48109–1107

USA

e2

E-mail: 
*e3

S. F. Taylor

Department of Psychiatry

University of Michigan

Ann Arbor, Michigan 48109–1107

USA

e4






While statistical analysis of a single network has received a lot of attention in recent years, with a focus on social networks, analysis of a sample of networks presents its own challenges which require a different set of analytic tools.  Here we study the problem of classification of networks with labeled nodes, motivated by applications in neuroimaging.   Brain networks are constructed from imaging data to represent functional connectivity between regions of the brain, and previous work has shown the potential of such networks to distinguish between various brain disorders, giving rise to a network classification problem.   Existing approaches tend to either treat all edge weights as a long vector, ignoring the network structure, or focus on graph topology as represented by summary measures while ignoring the edge weights.  Our goal is to design a classification method that uses both the individual edge information and the network structure of the data in a computationally efficient way, and that can produce a parsimonious and interpretable representation of differences in brain connectivity patterns between classes. We propose a graph classification method that uses edge weights as predictors but incorporates the network nature of the data via penalties that promote sparsity in the number of nodes, in addition to the usual sparsity penalties that encourage selection of edges.    We implement the method via efficient convex optimization  and provide a detailed analysis of data from two fMRI studies of schizophrenia.  











graph classification, high-dimensional data, variable selection, fMRI data






§ INTRODUCTION

Network data analysis has received a lot of attention in recent literature, especially unsupervised analysis of a single network which is thought of as generated from an exchangeable random graph model, for example <cit.>  and many others.   This setting is applicable to a number of real life scenarios, such as social networks, but there are situations where the network nodes are labeled and therefore not exchangeable, and/or more than one network is available for analysis, which have received relatively less attention.  Here we focus on the setting motivated by brain connectomics studies, where a sample of networks is available from multiple populations of interest (for example, mentally ill patients and healthy controls).    In this setting, each unit in the population (e.g., a patient) is represented by their own network, and the nodes (brain regions of interest) are labeled and shared across all networks through a registration process that maps all individual brains onto a common atlas.   There are many classical statistical inference questions one can ask in this setting, for example, how to compare different populations <cit.>.   The question we focus on in this paper is a classification problem: given a training sample of networks with labeled nodes drawn from multiple classes, the goal is to learn the rules for predicting the class of a given network, and just as importantly, interpret these rules.


Network methods are a popular tool in the neuroscience literature <cit.>. A brain network represents connectivity between different locations of an individual's brain. How connectivity is defined varies with the type of imaging technology used and the conditions under which data were collected.   In this paper, we focus on functional connectivity, which is a measure of statistical association between each pair of locations in the brain, constructed from functional magnetic resonance imaging (fMRI) data, although the methods we develop are applicable to any sample of weighted networks with labeled nodes. In fMRI studies, BOLD (blood oxygen-level dependent) signal, a known correlate of underlying neural activity, is measured at a sequence of time points at many spatial locations in the brain, known as voxels,  resulting in a 4-dimensional data array, with three spatial dimensions and a time index. Brain networks constructed from fMRI data have been successfully used for various tasks, such as differentiating between certain illnesses, or between types of external stimuli  <cit.>, and contain enough information to identify individual subjects <cit.>.    Extensive statistical literature has focused on the analysis of raw fMRI data <cit.>, usually aiming to characterize brain activation patterns obtained from task-based fMRI experiments. In this paper, we focus on resting-state fMRI data, where no particular task is performed  and subjects are free to think about anything they want.  Thus registering the time dimension across different subjects is not possible.    The connectivity network approach, which averages over the time dimension in computing a measure of dependence between different voxels, is thus a natural choice, and has been widely used with multiple types of neuroimaging data.  






Two different datasets are analyzed in this paper, both of resting state fMRI studies on schizophrenic patients and healthy controls.  One dataset, COBRE (54 schizophrenics and 70 controls), is publicly available <cit.>;  another, which we will refer to as UMich data (39 schizophrenics and 40 controls), was collected internally in the last author's lab.     Having two datasets on the same disease allows us to cross-check models trained on one of them for classification on the other to check the robustness of our approach. The raw data arrays undergo pre-processing and registration steps, discussed in detail in the Appendix <ref>, along with additional details on data collection.     To construct a brain network from fMRI measurements, a set of nodes is chosen, typically corresponding to regions of interests (ROIs) from some predefined parcellation. In our analysis we use the parcellation of <cit.> (see Figure <ref>), which consists of 264 ROIs divided into 14 functional brain systems (in the COBRE data, node 75 is missing).  A connectivity measure is then computed for every pair of nodes, resulting in an adjacency matrix of size 264× 264. Many choices of connectivity measures are available  <cit.>;  perhaps the most commonly used one is the Pearson correlation coefficient between locations, computed by averaging over the time dimension.   It has been argued that partial correlations are a better measure of connectivity <cit.>, but the choice depends on the final goal of analysis. In this paper we follow the vast majority of the connectomics literature and measure connectivity on each individual by using marginal correlations between the corresponding time series (see Figure <ref>).  The correlations are then further rank-transformed and standardized;  see Appendix <ref> for details.    These transformations are intended to deal with subject-to-subject variability and the global signal regression issue <cit.>, and although they lose some information, we observed that on our datasets classification based on standardized ranks of marginal correlations outperformed classification based on other connectivity measures, such as marginal correlations. The methods we develop here are applicable to networks that encode any type of connectivity measure.


The problem of graph classification has been studied previously in other contexts, with a substantial literature motivated by the problem of classification of chemical compounds <cit.>, where graphs represent the compound's molecular structure. This setting is very different, with small networks of about 20 nodes on average, binary or categorical edges recorded with no noise, and different nodes corresponding to different networks, <cit.>.   
Classification methods for chemical compounds are usually based on finding certain discriminative patterns in the graphs, like subgraphs or paths  <cit.>, and using them as features for training a standard classification method  <cit.>.  Computationally, finding these patterns is only possible on small binary networks.   

Another type of methods are based on graph kernels <cit.>, which  define a similarity measure between two networks. These kernels combined with support vector machines (SVMs) have been successfully used on small networks <cit.>, but the curse of dimensionality makes local kernel methods unsuitable for large scale networks <cit.>.    On our datasets, graph kernel methods did not perform better than random guessing.




In the context of classifying  large-scale brain networks, two main approaches have been followed.  One approach is to reduce the network to its global summary measures such as the average degree, clustering coefficient, or average path length <cit.>,  and use those measures as features for training a classification method. Previous studies have reported significant differences on some of these network measures for groups of patients with certain brain diseases compared with healthy controls 
<cit.>, suggesting their usefulness as diagnostic biomarkers. However, global summary statistics collapse all local network information, which can harm the accuracy of classification and do not allow to identify local differences.   In our data analysis,  a method based on the network measures suggested in <cit.>  performed poorly for classification  (see Section <ref>).

An alternative approach to classification of large networks is to treat edge weights as a “bag of features”, vectorizing the unique elements of the adjacency matrix and ignoring the network nature of the data.    This approach can leverage many existing classification methods for vectors, and provides an interpretation at the edge level if variable selection is applied <cit.>.  Spatial correlation between edges connecting neighboring nodes can be incorporated <cit.>, although the effectiveness of this regularization will depend on the parcellation used to define nodes. Alternatively, an individual test can be used for each edge to find significant differences between two populations, with a multiple testing correction and without constructing a classifier at all <cit.>. While these methods can deliver good predictions, their interpretability is limited to individual edge selection, which is less scientifically interesting than identifying differentiating nodes or regions, and they cannot account for network structure. 


Taking the network structure into account can have benefits for both testing and classification settings. Some methods perform inference over groups of edges based on the community assignments of the nodes to which they are incident. For example, <cit.> introduced Network Contingency Analysis which begins with massive univariate testing at each edge, and then counts the number of superthreshold connections in each cell, a group of edges that connect nodes in two functional systems.  Nonparametric methods are then used to conduct inference on the count statistic for each cell, with multiple comparison correction for inference at the cell level.    Power can be improved  by applying a network-based  multiple testing dependence correction <cit.>.  For classification,  better interpretability and potentially accuracy can be obtained if we focus on understanding which brain regions or interactions between them are responsible for the differences.     In somewhat related work, <cit.> proposed to look for a minimal set of nodes which best explains the difference, though that requires solving a combinatorial problem.   Hypothesis testing on a type of graph average has also been proposed <cit.>. Bayesian nonparametrics approaches for modeling populations of networks allow to test for local edge differences between the groups <cit.>, but are computationally feasible only for small networks.  

Our goal in this paper is to develop a high-dimensional network classifier  that uses all the individual edge weights but also respects the network structure of the data and  produces more interpretable results.   To achieve this goal, we use structured sparsity penalties to incorporate the network information by penalizing both the number of edges and the number of nodes selected.  Although our main application here is classification of brain connectivity networks, our methods are applicable to any weighted graphs with labeled nodes, and to general prediction problems, not just classification.  

The rest of this paper is organized as follows.   In Section <ref>, we introduce our classifier and the structured penalties.  In Section <ref> we show how to efficiently solve the resulting convex optimization problem by a proximal algorithm, each step of which is a further optimization problem which we solve by the alternating direction method of multipliers (ADMM). 
The performance of our method is evaluated and compared with other methods using simulations in 
Section <ref>.
In Section <ref>, we analyze two brain connectivity datasets, each containing schizophrenic patients and healthy controls, and show that our regularization framework leads to state-of-the-art  accuracy while providing interpretable results, some of which are consistent with previous findings and some are new.  We conclude with a brief discussion in Section <ref>.







§ A FRAMEWORK FOR NODE SELECTION IN GRAPH CLASSIFICATION




 §.§ A penalized graph classification approach

We start from setting up notation.  All graphs we consider are defined on the same set of N labeled nodes. 
A graph can be represented with its adjacency matrix A∈^N× N.   We focus on graphs that are undirected (A_ij = A_ji) and contain no self-loops (A_ii = 0).   These assumptions are not required for the derivations below, but they match the neuroimaging setting and simplify notation. Our goal is  predicting a class label Y from the graph adjacency matrix A;  in this paper we focus on the binary classification problem where Y takes values {-1,1}, although extensions from binary to multi-class classification or real-valued responses are straightforward.   Throughout this paper, we use ·_p to denote the entry-wise ℓ_p norm, i.e., for a matrix A∈^m× n, A_p=(∑_i=1^m∑_j=1^n|A_ij|^p)^1/p.




A standard general approach  is to construct a linear classifier, which predicts the response Y from a linear combination of the elements of A, ⟨ A,B⟩ = B^TA, where we arrange the coefficients in a matrix B∈^N× N to emphasize the network nature of the predictors.    We focus on linear classifiers here because variable selection is at least as important as prediction itself in the neuroimaging application, and setting coefficients to 0 is a natural way to achieve this.    The coefficients are typically estimated from training data by minimizing an objective consisting of a loss function plus a penalty.  The penalty can be used to regularize the problem to  make the estimator well-defined in high-dimensional problems, to select important predictors, and to impose structure, and many such penalties have been proposed, starting from the lasso <cit.>.  Our focus is on designing a classifier in this framework that respects and utilizes the network nature of the predictors.    In brain networks in particular, neuroscientists believe that edges are organized in subnetworks, also called brain systems <cit.>, that carry out specific functions, and 
certain subnetworks are important for prediction <cit.>, although different studies tend to implicate different regions <cit.>.   Thus we aim to find nodes or subnetworks with good discriminative power, and hence select both the most informative nodes and edges.

Although the methods we develop here can be used on small networks,  our main focus is on the more challenging case of medium to large brain networks. In brain connectivity studies dealing with multiple subjects, while raw images may contain hundreds of thousands of voxels, they are commonly down-sampled according to a parcellation scheme with a coarser resolution, usually resulting in networks with hundreds or thousands of nodes representing ROIs (see for example <cit.>). This coarser resolution is essential for registration, as aligning different brains at a high resolution is much harder, but it still results in hundreds of thousands or millions of edges which serve as predictor variables.   Given the typical data sizes in this area of application, we focus on methods based on convex formulations, which allow for efficient and scalable implementations with convergence guarantees.  





Let {(A^(1),Y_1),…,(A^(n),Y_n)} be the training sample of undirected adjacency matrices with their class labels, 
and let  Y=(Y_1,…,Y_n).  
A  generic linear classifier described above is computed by finding the coefficients B defined by 

    B̂ =   _B∈{ℓ(B )  + Ω(B)},

where = {B∈^N× N:B=B^T, diag(B) = 0}, Ω is a penalty, and 

    ℓ(B) = 1/n∑_k=1^n ℓ̃(Y_k, A^(k); B )

is a loss function evaluated on the training data.  Our methodology can accommodate different choices of loss functions that can extend beyond classification problems (e.g., least squares or generalized linear models). The optimization algorithm presented in Section <ref> can work with any convex and continuously differentiable loss function, and further assumptions are required for consistency (see Section <ref>). In this paper, for the purpose of classification we use the logistic loss function in the simulations and data analysis, which is defined as 

    ℓ̃(Y_k, A^(k);B,b) = log{1+exp(-Y_k ⟨ A^(k), B⟩+b)} .

The threshold b is an additional parameter to be estimated.


To capture structural assumptions on important predictive edges, we focus on convex structured sparsity penalties <cit.> that  encourage a small number of active nodes, by which we mean nodes attached to at least one edge with a non-zero coefficient. One approach to finding a set of such nodes was proposed by <cit.>, who called it a signal-subgraph, and proposed finding the minimal set of nodes (called signal vertices) which together are incident to all selected edges (but not every node connected to a selected edge is a signal vertex).    Finding this set is a combinatorial optimization problem, and the set is not always uniquely defined.    Instead, we focus on convex formulations that allow for efficient computation and encourage small active node sets indirectly. 


Other convex penalties have been used for fMRI data as a way to enforce spatial smoothness in the solution  <cit.>. These methods assume that voxels are equally spaced in the brain, and neighboring voxels are highly correlated. In particular, <cit.> proposed penalties for brain network classification using these spatial assumptions.  Here, instead of enforcing a spatial regularization directly, we aim for a regularization that can be applied to any type of network data, and in particular to brain networks with coarse and/or uneven parcellations where enforcing spatial smoothness may not work as well.   In any case, the flexibility of convex optimization algorithms allows one to easily incorporate additional spatially-informed penalties if needed.





 §.§ Selecting nodes and edges through group lasso


To reflect the network structure of the predictors, we use a penalty that promotes a sparse classifier not only in the number of edges used, but also in the number of nodes.   The group lasso penalty  <cit.> is designed to eliminate a group of variables simultaneously.   Here we penalize the number of active nodes by treating all edges connected to one node as a group. Then eliminating this group (a row of coefficients in the matrix B) is equivalent to de-activating a node. The group penalty is defined as

    Ω_λ, ρ(B)= λ(∑_i=1^NB_(i)_2 + ρB_1),

where B_(i) denotes the vector of edge weights incident to the i-th node (or equivalently, the i-th row or column of B), and λ, ρ≥ 0 are tuning parameters.  Note that the constraint B=B^T makes the groups overlap, since a coefficient B_ij  belongs to groups associated with the nodes i and j, and therefore, the edge between nodes i and j would be selected only if both nodes are activated.    The second term in the penalty ρB_1 acts as the usual lasso penalty to promote sparsity inside the group <cit.>, allowing to select a subset of edges for an active node. Due to the overlap in the groups, this lasso penalty is usually necessary in order to produce sparse solution (see Proposition <ref>). The constraint diag(B) = 0 in (<ref>) is automatically enforced with this formulation. 


 
	An alternative to the constraint in the problem (<ref>) is to optimize over the set 
    = {B∈^N× N, diag(B) = 0}.
 Without the symmetry constraint and assuming undirected graphs, the penalty (<ref>) is equivalent to the overlapping group lasso formulation of <cit.>. This formulation has some advantages. Since it gives group lasso without overlaps, the lasso penalty ρB_1 is not required to obtain sparse solutions, and more efficient optimization algorithms exist for this case. This approach would loosely correspond to the idea of selecting signal nodes as in <cit.>, in the sense that an edge can be selected if at least one of its nodes is selected, and the second node could be inactive. The downside is that each edge now corresponds to two different coefficients B_ij and B_ji, the problem encountered by all variable selection methods that ignore symmetry, such as <cit.>. The standard solution for this problem, as suggested by <cit.>, is to take the average of the coefficients. 
	Intuitively, one would expect that the formulation using  would be better when the significant edges are incident to a small set of nodes, since both nodes have to be active for an edge to be selected, while using  may be better when for some nodes most of their edges are significant, creating “significant hubs”.  Since in our application we are primarily looking for discriminative brain subnetworks, we focus on the symmetrically constrained formulation for the rest of the paper.  We also found that in practice this second formulation results in less accurate classifiers for the neuroimaging data discussed in Section <ref>.


  The analogue to (<ref>) for directed graphs would assign coefficients B_ij and B_ij to the same group, resulting in the penalty 
	
    Ψ_λ, ρ(B)= λ(∑_i=1^N√(∑_j(B_ij^2+B_ji^2)) + ρB_1),

	where B∈ℬ̃. Alternatively, we can also use the formulation of Remark <ref>, by replicating the variables and estimating two matrices of coefficients, say B^(1) and B^(2), with the penalty 
	
    Ψ̃_λ, ρ(B^(1),B^(2)) = λ[∑_i=1^N√(∑_j{(B^(1)_ij)^2+(B^(2)_ji)^2}) + ρ(B^(1)_1 + B^(2)_1 )],

	with B^(1),B^(2)∈ℬ̃, and set the coefficients matrix to B=(B^(1) + B^(2))/2.   This formulation will again not directly select subnetworks as discussed in Remark <ref>. 




Finally, for numerical stability we add an extra ridge penalty term γ/2B_F^2=γ/2Tr(B^TB), with γ a small and fixed constant. There are several benefits of combining ridge and lasso penalties (see for example <cit.>).  The parameter γ can be potentially considered as an additional tuning parameter, but  here we only use a small fixed constant γ in order to avoid numerically degenerate solutions. In practice, the results are not sensitive to the exact value of γ.   

Putting everything together, to fit our graph classifier, we solve the problem 

    (B̂,b̂) =  _B∈ℬ,b∈   {1/n∑_k=1^nlog(1+exp(-Y_k ⟨ B,A^(k)⟩+b)) + γ/2B_F^2 +  . 
       .  λ( ∑_i=1^NB_i_2+ρB_1)}

for given values of λ and ρ, which will be chosen by cross-validation.





§ THE OPTIMIZATION ALGORITHM


Our optimization algorithm to solve the problem (<ref>) combines two common approaches to convex optimization: proximal algorithms and alternating direction method of multipliers (ADMM).  We use an accelerated version of the proximal algorithm <cit.> to solve the main problem (<ref>). In each step, we need to calculate a proximal operator, which is a further convex optimization problem solved with the ADMM algorithm.

The main optimization difficulty comes from the overlapping groups. Some algorithms have been proposed for this case, including a subgradient descent method <cit.>, which has a slow rate of convergence, or proximal algorithms based on smoothing the original problem <cit.>. Although smoothing yields fast algorithms, it is not clear that the sparsity pattern is preserved with those approximations. We follow an approach similar to <cit.> and <cit.>, but solve the proximal operator for the penalty (<ref>) directly using the ADMM method. This can potentially give a more accurate sparsity pattern, and the flexibility of the algorithm allows for additional penalties if desired, such as spatial smoothing similar to <cit.> (see Remark <ref>).


The main problem (<ref>) is solved with a proximal algorithm (see <cit.>). Recall that the proximal operator  for a function f is defined as _f(v) = _x{f(x)+1/2x-v_2^2}. Starting with an initial value B^(0)∈^N× N, a proximal algorithm solves the optimization problem (<ref>) 
by iteratively calculating the proximal operator of Ω=Ω_λ,ρ for a descent direction of the differentiable loss function ℓ.  We use an accelerated version of the algorithm <cit.>, which for each k=2,…, until convergence, performs the updates

    W^(k)    =  B^(k-1) + k-1/k+2(B^(k-1) - B^(k-2)) 
    
    	B^(k)    =  _t^(k)Ω{W^(k) - t^(k)∇ℓ(W^(k))}
        =  _B∈{1/2B-(W^(k)-t^(k)∇ℓ(W^(k)))_2^2 + t^(k)Ω(B) },

where ∇ℓ(W)∈^N× N is the gradient of the loss function ℓ at W and {t^(k)} is a sequence of positive values. If ∇ℓ is Lipschitz continuous, with L its Lipschitz constant, the sequence of values ℓ(B^(k))+Ω(B^(k)) converges to the optimal value at rate O(1/k^2) if t^(k)∈[0,1/L).    The value of t^(k) can be chosen using a backtracking search <cit.>, which decreases this value until the condition

    ℓ(B^(k)) ≤ℓ(W^(k)) + ⟨∇ℓ(W^(k)), B^(k)-W^(k)⟩  + 1/2t^(k)B^(k)-W^(k)_2^2

is satisfied. This procedure ensures that step sizes {t^(k)} become smaller as the algorithm progresses, until t^(k)<1/L. In practice, L might be large, which can make the algorithm slow to converge. It has been observed in other sparse high-dimensional problems that search strategies for t^(k) which allow for t^(k)  > 1/L when appropriate can actually speed up convergence  <cit.>. 
We use a strategy of this type, allowing t^(k) to increase by a factor of α≥ 1 if the relative improvement in the loss function on iteration k becomes small. We observed that this strategy significantly reduces the number of iterations until convergence. The entire procedure is summarized in Algorithm <ref> on the Appendix <ref>.



The logistic loss function of (<ref>) has an extra parameter b. Rather than including it as an unpenalized coefficient for a constant covariate, we use block coordinate descent and solve for b separately. This is convenient because the threshold b and the matrix of coefficients B may not be on the same scale. Thus, b can be updated by solving b^(k+1)  = _b∈ℓ(B^(k),b), which is easy to compute via Newton's method.

The proximal algorithm requires solving the proximal operator (<ref>), which has no closed form solution for the penalty (<ref>) under the symmetry constraint. Strategies based on smoothing this penalty have been proposed <cit.>. However, to allow for variable selection which results from non-differentiability of the penalty, we aim to solve the proximal operator directly using ADMM  (see <cit.> for a review). Note that if the symmetric constraint is relaxed as in Remark <ref>, the proximal operator has a closed form solution (see Remark <ref>).

The ADMM works by introducing additional constraints and performing coordinate descent in the corresponding augmented Lagrangian function. Setting Z=W^(k)-t^(k)∇ℓ(W^(k)) and t=t^(k), and introducing the variables Q,R∈^N× N, we can formulate (<ref>) as a convex optimization problem

    min_B̃,Q,R      1/2B̃-Z_2^2 + tλ(∑_i=1^NQ_(i)_2 + ρR_1)
       subject to      B̃ = Q,   B̃ = R,   B̃ = B̃^T,  diag(B̃) = 0.

The ADMM algorithm introduces the multipliers U,V∈^N× N and a penalty parameter μ>0 to perform gradient descent on the Lagrangian of (<ref>), given by ℒ_μ=ℒ_μ(B̃,Q,R,U,V) as

    ℒ_μ  =    1/2B̃-Z_2^2+ tλ(∑_i=1^NQ_(i)_2 + ρR_1) +  ⟨ U, B̃-Q⟩+
       ⟨ V, B̃-R⟩μ/2(B̃-Q_2^2+B̃-R_2^2 + B̃-B̃^T_2^2).

The value μ controls the gap between dual and primal feasibility. In practice, we observed that setting μ = 0.1 gives a good balance between primal and dual feasibility, although other self-tuning methods are available <cit.>. This function is optimized by coordinate descent, with each variable updated to minimize the value of ℒ_μ while all the other variables are fixed. This update has a closed form; the detailed steps of the ADMM are shown in Algorithm <ref> on the Appendix <ref>. These steps are performed until the algorithm converges within tolerance ϵ^>0. 
Note that ADMM will be performed in each iteration of the algorithm to solve (<ref>) and thus tolerance ϵ^ can be decreased as the algorithm progresses. On the other hand, performing only one iteration of algorithm (<ref>) gives a similar algorithm to the one of <cit.>.






	The ADMM  makes it very easy to incorporate additional penalties.  If Ψ is a new penalty, we can rewrite (<ref>) by introducing an additional parameter Q̃ so it becomes
	
    min_B̃,Q,Q̃,R      1/2B̃-Z^(k)_2^2 + tλ(∑_i=1^NQ_(i)_2 + ρR_1) + tΨ(Q̃)
       subject to      B̃ = Q,  B̃ = Q̃,   B̃ = R,   B̃ = B̃^T,  diag(B̃) = 0.

	We can obtain the Lagrangian formulation (<ref>) in a similar manner, and include new parameters in the ADMM updates, which can be performed efficiently as long as the proximal operator of Ψ has a closed form solution. This is in fact the case for some other penalties of interest, such as the GraphNet penalty <cit.>, which can incorporate spatial location information.  



	The alternative formulation for the graph penalty given in Remark <ref> corresponds to standard sparse group lasso <cit.>. In particular, we can still employ the proximal algorithms (<ref>) and (<ref>), but instead optimize over the set . Without the symmetric constraint on B, the overlap in the group lasso penalty disappears, and this vastly simplifies the problem. Using Theorem 1 of <cit.>, the update for B^(k) has a closed form solution given by
	
    Y^(k)    =  B^(k-1) + k-2/k(B^(k-1) - B^(k-2)) 
    
    		Z^(k)_ij    = (1- λρ/|Y^(k)_ij-t_k∇_ijℓ(Y^(k))_2)_+(Y_ij^(k)-t_k∇_ijℓ(Y^(k)))
    
    		B^(k)_(i)    =  (1-λ/Z^(k)_(i)_2)_+(Z^(k)_(i)),     i∈[N].









§ THEORY


In this section, we show that the solution of the penalized problem (<ref>) can recover the correct subgraph corresponding to the set of non-zero coefficients, and give its rates of convergence.  The theory is  a consequence of the results of <cit.> for establishing model selection consistency of regularized M-estimators under geometric  decomposability (see Appendix for details). We present explicit conditions for our penalty to work well, which depend on the data as well as the tuning parameters.

Let B^⋆⊂^N× N be the unknown parameter we seek to estimate, and we assume there is a set of active nodes 𝒢⊂[N] with |𝒢|=G, so that  B^⋆_ij=0 if i∈𝒢^C or j∈𝒢^C. We allow some edge weights inside the subgraph defined by 𝒢 to be zero, but we focus on whether the set 𝒢 is correctly estimated by the set 𝒢̂ of active nodes in B̂.   Denote by ℳ⊆^N× N the set of matrices where the only non-zero coefficients appear in the active subgraph, that is,

    ℳ={B∈^N× N| B_ij=0 for all i∈𝒢^C or j∈𝒢^C, B=B^T.}

There are two main assumptions on the loss function ℓ required for consistent  selection in high-dimensional models <cit.>. The first assumption is on the convexity of the loss function around B^⋆, while the second assumption bounds the size of the entries in the loss Hessian between the variables in the active subgraph and the rest. Let the loss Hessian ∇^2ℓ(B^⋆)∈^N× N⊗^N× N be defined by 

    ∇^2_(i,j),(k,l)ℓ(B) = ∂^2ℓ(B)/∂ B_ij∂ B_kl,

and define the matrix H_(i,j),𝒢∈^G× G with (i,j)∈(𝒢×𝒢)^C such that 




    (H_(i,j),𝒢)_k,l = Tr{(∇^2_(i,j),(𝒢,𝒢)ℓ(B^⋆))
    		Λ_(k,l),(·,·)},   1≤ k,l≤ G,,

where Λ∈^G× G⊗^G× G is a tensor such that Mat(Λ) is a pseudoinverse of Mat(∇^2_(𝒢,𝒢),(𝒢,𝒢)ℓ(B^⋆)), and Mat is the operation that unfolds the entries of a tensor Λ into a G^2× G^2 matrix. The matrix H_(i,j),𝒢 measures how well the variable corresponding to the edge (i,j) can be represented by the variables in the active subgraph.
assumpAssumption
[Restricted Strong Convexity] 
	There exists a set C⊂^N× N with B^⋆∈ C, and constants m>0, L̃<∞ such that
	
    ∑_i,jΔ_i,jTr{(∇^2_(i,j),(·,·)ℓ(B))Δ}
    		   ≥  m Δ_2^2,   ∀ B∈ C∩ℳ,Δ∈ C∩ℳ
    ∇^2ℓ(B) - ∇^2ℓ(B^⋆)_2    ≤L̃B-B^⋆_2,  ∀ B∈ C.


[Irrepresentability] There exists a constant 0<τ<1 such that
	
    max_i∈𝒢^C(∑_k=1^G(H_(i,j),𝒢)_k·_2)_j=1^N_2= 1-τ<0.

	
	
	
	
	
	


This version of the irrepresentability condition corresponds to the one usually employed in group lasso penalties <cit.>, but as we will see later, due to overlaps in the groups it further requires a lower bound on ρ to work for model selection.

The first two assumptions are stated directly as a function of the loss for a fixed design case, but they can be substituted with bounds in probability for the case of random designs. In order to obtain rates of convergence, we do require a distributional assumption on the first derivative of the loss. This assumption can be substituted with a bound on max_i∇ℓ(B^⋆)_(i)_2 (see Lemma <ref> in the Appendix).
[Sub-Gaussian score function]  Each pair in the sample (A,Y) is independent and comes from a distribution such that the entries of the matrix ∇ℓ̃(Y,A;B^⋆) are subgaussian.  That is, for all t>0 there is a constant σ^2>0 such that
	
    max_i,j(∇_ijℓ̃(Y, A; B^⋆)_∞>t) ≤ 2exp(-t^2/σ^2)  .


With these assumptions, we establish consistency and correct model selection. The proof is given in the Appendix..
propProposition 

	Suppose Assumptions <ref> and <ref> hold. 
	
		
  (a) Setting the penalty parameters as  λ = c_1σ√(log N/n)min{√(N)/1+ρ, 1/ρ} and ρ≥ 0 for some constant c_1>0, with probability at least 1-2/N the optimal solution of (<ref>) is unique and satisfies
		
    B̂-B^⋆_2 = O_P(σ N√(log N/n)).

		
  (b) Suppose Assumption <ref> also holds.  If  n> c_2G^2σ^2log N for a constant c_2>0, setting the penalty parameters as  λ = c_3σ√(log N/n)min{√(N)/1+ρ, 1/ρ} for some constant c_3>0, and 
		
    ρ > 1/τ - 1/√(G) ,

		then 
		
    B̂-B^⋆_2     = O_P(σ G√(log N/n)), 
    P(𝒢̂⊆𝒢)     = 1-2/N.

	


The part of the penalty associated with ρ causes the solution to be sparse. Due to the overlap in the groups, a small value of ρ will usually not result in zeros in the solution of the problem (<ref>).   The lower bound on ρ in (<ref>) ensures that the irrepresentability condition of <cit.> holds (see Lemma <ref> in the Appendix). 

Proposition <ref> ensures that, with high probability, all edges estimated to have non-zero weights are contained in the active subgraph, and quantifies the error in estimating the entries of B^∗. To ensure that all active nodes are recovered, at least one edge corresponding to each active node needs to have a non-zero weight. A similar result can be obtained to guarantee recovery of all active nodes under a stronger assumption that the magnitude of the non-zero entries of B^⋆ is bounded below by |B_ij^⋆|>c_4 G^2λ for a constant c_4. The condition in part (b) requires a sample size that grows faster than the size of the active subgraph, which in practice can be much smaller than the size of the network, making the method suitable for our applications in which the sample size is limited and the number of nodes is large. 





§ NUMERICAL RESULTS ON SIMULATED NETWORKS


In this section, we evaluate the performance of our method using synthetic networks. We are interested in assessing the ability of the method to correctly identify predictive edges, its classification accuracy, and comparisons to benchmarks.  
We compare the different methods' edge selection performance in simulations using area under the curve (AUC).   
Brain connectomic networks are characterized by organization of nodes into communities <cit.>, in which nodes within the same community tend to have stronger connections than nodes belonging to different communities. In order to generate synthetic networks that mimic this property, we  introduce community structure using the stochastic block model (SBM)  <cit.>.    Before generating edges, we assign each node a community label,  C_i, where C_i∈[K] for each i∈[N]. The node assignments are the same for all networks in the population. Given the community labels,  network edges are generated independently from a distribution that only depends on the community labels of the nodes associated with each edge. Since fMRI networks are real-valued networks, we generate edge weights from a Gaussian distribution, rather than the standard Bernoulli distribution normally used with the SBM.   Specifically, we draw each A_ij independently from  N(μ_C_iC_j, σ^2), with μ∈R^K× K defined by 

    μ_kl = {[      0.3,   if k=l,;       0.1 if k ≠ l, ].

and σ^2=0.18. These values were chosen to approximately match the distribution of edge weights in our datasets (see Section <ref>). We set the number of nodes N = 300, with K=12 communities of size 25 each. We work with undirected networks, so the adjacency matrices are symmetric, with 44,850 distinct edges. Although our method is able to scale to larger networks, this moderately sized setting is already highly computationally demanding for many of the comparison benchmarks.

To set up two different class distributions, we select a set of active nodes 𝒢 first, defined by the nodes corresponding to some communities selected at random.
Then, we alter a set of differentiating edges ℰ selected at random from 𝒢×𝒢 with probability p. For each edge (i,j)∈ℰ, the distribution in class Y=-1 is  N(μ_C_iC_j, σ^2), while the distribution in  class Y=1 is changed to  N(0.2, σ^2). Figure <ref> shows example expected adjacency matrices for each class. We then generate 50 networks from each class, resulting in a sample size of n=100. We vary  G=|𝒢| by changing the number of communities selected, and the value of p, to study the effect of the number of active nodes and the density of differentiating edges inside a subgraph.  The number of communities selected is set to 1 (|𝒢|=25), 2 (|𝒢|=50) and all communities (|𝒢|=300); note that in  the last scenario all nodes are active and hence the node structure is not informative at all. 




Since we are interested is identifying predictive edges and nodes, we use the AUC of the receiver operating characteristic (ROC) curve, for both edge and node selection. For each method, we calculate the ROC curve by changing its corresponding sparsity parameter to vary the number of edges selected . For a selection method ℳ and a sparsity parameter η let ℰ̂(ℳ, η) be the set of edges selected by ℳ, and 𝒢̂(ℳ, η) the set of active nodes corresponding to ℰ̂(ℳ, η). We calculate the edge false positive rate (EFPR) and edge true positive rate (ETPR) as

    EFPR(ℳ, η)  =  |ℰ̂(ℳ, η)∩ℰ^C |/|ℰ^C| ,    ETPR(ℳ, η)  =  |ℰ̂(ℳ, η)∩ℰ|/|ℰ| .

The node FPR and TPR are calculated similarly. 

We also evaluate the prediction accuracy of the methods. For each method, we use 5-fold cross-validation to select the best tuning parameter using the training data, and then compute the test error on a different dataset simulated under the same settings.   The AUC and test errors reported are averaged over 30 replications.  


Methods for benchmark comparisons on simulated networks were selected based on their good performance on real data  (see Section <ref>). For our method (GC), we vary the parameter ρ and compare results for two different values of λ, .05 (GC1) and 10^-4 (GC2). For unstructured regularized logistic regression, we use the elastic net <cit.>, with a fixed α=0.02 (ENet). The performance of elastic net is not very sensitive to different values of α, but the number of variables that the method is able to select with large values is limited (including the case of α=1 that corresponds to the lasso). A support vector machine with ℓ_1 penalty <cit.> is also included (SVML1) for comparison, and additionally we evaluate the classification error of the original support vector machines (SVM)  <cit.>.  For both SVMs, we use linear kernels, which performed better than nonlinear ones. We also consider an independent screening method for variable selection based on the two sample t-statistic (T-stat). Finally, we also compare with the signal-subgraph method (SS) <cit.> which is the only other method that takes into account the network structure of the predictor variables. Note that the signal subgraph is designed for binary networks, so in order to apply it we thresholded each edge at the population mean. For each method, we fit 10 different tuning parameters to change the sparsity of the solution.











Figure <ref> shows the values of the average AUC for selecting edges (top) and nodes (bottom).   For G=25 and 50, 
as the proportion of differentiating edges in the active subgraph increases, methods that take into account network structure (GC1, GC2 and SS) slightly improve their edge AUC, since enforcing node selection also results in better edge selection, while the edge AUC remains constant for unstructured methods (ENet, T-stat and SVML1).   On node selection, all methods improve the node AUC as the fraction of significant edges increases, but GC and SS have the largest gains. A similar trend is observed in classification error shown in Figure <ref>. All methods improve as the proportion of differentiating edges increases, but our method has the best performance overall.  Our method performed the best with the larger value of λ (GC1) on variable selection, particularly when the set of active nodes is smaller, but both values of λ give very good classification performance.  In the last scenario (G=300), all nodes are active so the node AUC is undefined, and the node structure is not informative at all. Although the performance of our method is no longer the best, it performs comparably to state of the art methods that do not use network structure.

 In terms of computing time, since there are many contributing factors including the software choice for implementation and the tuning parameters, a fair comparison is difficult.    We can roughly say that elastic net is the fastest, taking about a minute to run a cross-validation instance, while our method takes about 10 minutes on average, and the  signal-subgraph takes more than an hour.



§ APPLICATION TO SCHIZOPHRENIA DATA

We analyze the performance of the classifier on the two different brain fMRI datasets previously described in Section <ref>. 
The code of our classifier and the processed connectomics datasets are available at <https://github.com/jesusdaniel/graphclass>. 








 §.§ Classification results








First, we evaluate our method's classification accuracy. We use a nested 10-fold cross-validation to choose tuning parameters and estimate the test accuracy. The classifier is trained for a range of values of λ and ρ, with λ∈{10^-7, 10^-6.5, …,  10^-2} and ρ∈{10^-3, 10^-2.5, …,  10^2}. The value of γ in (<ref>) is set to 10^-5; we observed that setting γ to a small value  speeds up convergence without affecting the accuracy or sparsity of the solution.   Figure <ref> shows the average cross-validated accuracy, sparsity (fraction of zero coefficients) and node sparsity (fraction of inactive nodes), as a heat map over the grid of tuning parameter values.  We observe that  λ has little influence on sparsity, which is primarily controlled by ρ. Moreover, as Proposition <ref> suggests, values of ρ < 1 do not result in node selection.  As expected, accuracy generally decreases as the solution becomes sparser,  which is not uncommon in high-dimensional settings <cit.>. However, we can still achieve excellent accuracy with a substantially reduced set of features. In the COBRE dataset, the best accuracy is obtained with only 1886 edges (5.4%) but almost all nodes are active (260).   On the UMich data, 29733 edges (85.6%) achieve the best performance, and all nodes are active. Choosing parameters by cross-validation often tends to include too many noise variables <cit.>, as we also observed in simulations. A commonly used technique to report solutions that still achieve good accuracy with a substantially reduced set of features is the so-called “one-standard-error rule" <cit.>, in which one selects the most parsimonious classifier with cross-validation accuracy at most one standard error away from the best cross-validation accuracy.  Figure <ref> shows the solutions for each dataset obtained by this rule. Nodes are ordered by brain systems (see Figure <ref>). The fitted solution for COBRE has  549 non-zero coefficients (1.56%) and 217 active nodes (82.5%), while the UMich solution has 11748 non-zero entries (33.8%), and all nodes are active. Note that when many variables are selected, the magnitude of the coefficients becomes small due to the grouping effect of the penalty <cit.>. 






We also compared our method to benchmarks (Table <ref>), using the same methods as in the previous section and training and evaluating all methods using with the same nested 10-fold cross-validation. For SVM, we tested different kernels, including graph aware kernels <cit.>, but in most cases local kernel methods were no better than random guessing.  We additionally included random forests  and a method based on global and local network summaries previously proposed as features for classifying brain data <cit.>. For the latter, because our dataset is much larger, we only considered  global and node features proposed in <cit.>, which resulted in about 30,000 features per individual, and omitted edge features.     <cit.> evaluated their classifiers on a different parcellation of the COBRE data, and we do not include their methods since they are based on the assumption of equally spaced nodes and cannot be directly applied to our data. Their reported accuracy of 71.9% and 73.5% for the COBRE data is substantially lower than our method, although the results are not directly comparable.



Results in Table <ref> show that most methods performed better on the COBRE dataset than on the UMich dataset, which can be partially explained by the different sample sizes and possibly noise levels.  Besides differences in sample size and demographic characteristics (Table <ref>), the COBRE dataset is more homogeneous as it was collected using identical acquisition parameters, whereas the UMich dataset was pooled across five different experiments spanning seven years. 

Our method performs very well on both datasets, particularly among methods that  do variable selection. SVMs, which use the hinge loss, perform well too, and generally outperform methods using the logistic loss.   Our penalty can be combined with any loss, so we could also include our penalty combined with hinge loss which might potentially improve classification accuracy, but we do not pursue this direction, for two reasons:  one, our method is close to SVM + L1 as it is (better on COBRE, slightly worse on UMich but the difference is within noise levels), and because solutions based on logistic loss are generally considered more stable and preferable for variable selection <cit.>.   In Figure <ref>, we plot cross-validated classification accuracy of these methods as a function of the number of variables selected.  For the COBRE data, as we have observed before, good accuracy can be achieved with a fairly small number of edges, and the noisier UMich data requires more edges.   In all cases, our method uses fewer nodes than the others, as it is designed to do so.  



















Ultimately, assessing significance of the selected variables is necessary, which is in general a difficult task in high-dimensional settings and an active area of research (see for example <cit.>). In brain connectomics, it is particularly challenging to identify significant variables because of small sample sizes <cit.>.   Here we employ stability selection <cit.> which can be shown to control a type of false discovery rate by employing many rounds of random data splitting and calculating the probability of each variable being selected.  Some versions of this method have been theoretically studied, and upper bounds on the expected number of variables with a low selection probability that are included in the final solution (i.e., errors) have been derived under mild conditions <cit.>.    We implemented the version of stability selection proposed by <cit.>, with values of λ and ρ obtained by cross-validation on the COBRE data, and by the “one standard error rule" on the UMich dataset, since stability selection is most relevant to sparse solutions.   However, one of the advantages of stability selection is that it is not sensitive to the initial choice of tuning parameters, and changing tuning parameters  only slightly alters the ordering of variables with the largest selection probabilities.


The edges with the 15 largest selection probabilities are reported in Table <ref>. 
Using the results of <cit.> (equation 8), we estimated that the expected number of falsely selected variables (variables with a probability of selection smaller than the estimated) is bounded by 6.1 for the COBRE dataset and 9.7 for the UMich data,  which also suggests that results on the UMich data might be less reliable. While the two datasets yield somewhat different patterns of edge selection, it is notable that the default mode network (5)  was often selected in both.   This network has been consistently implicated in schizophrenia <cit.>, as well as other psychiatric disorders, possibly as a general marker of psychopathology <cit.>. In the COBRE dataset, edges were also selected from the fronto-parietal task control region (8), previously
linked to schizophrenia <cit.>. These results coincide with the findings of <cit.> on a different parcellation of the same data, which is an encouraging indication of robustness to the exact choice of node locations. Some of the variables with the highest estimated selection probabilities appear in the uncertain system (-1), in particular in the cell connecting it with salience system (9), which suggests that alternative parcellations that better characterize these regions may offer a better account of the schizophrenia-related changes.  
Additionally, sensory/somatomotor hand region (1) and salience system (9) also stand out in the UMich data, and these are networks that have also been implicated in schizophrenia <cit.>.

While results in Table <ref> do not fully coincide on the two datasets, there are clear commonalities.     Table <ref> compares classification accuracy when the classifier is trained on one dataset and tested on the other (with the exception of the intercept, since the datasets are not centered in the same way, which is fitted on a part of the test data, and the test error is then computed via 10-fold cross-validation).     While the accuracy is lower than when the same dataset is used for training and testing, as one would expect, it is still reasonably good and in fact better than some of the benchmark methods even when they train and test on the same data.  We again observe that the COBRE dataset is easier to classify.  




Figure <ref> shows the active nodes in the COBRE dataset (marked in green), corresponding to the endpoints of  the edges listed in Table <ref>.   We also identified a set of 25 nodes that are not selected in any of the sparse solutions with cross-validation accuracy within one standard error from the best solution (marked in purple).   These consistently inactive nodes are mostly clustered in two anatomically coherent regions.  







§ DISCUSSION

We have presented a method for classifying graphs with labeled nodes, motivated by brain connectomics but generally applicable to any setting with such graphs.   The distinct feature of our method is that it is graph-aware, aiming to select a sparse set of both edges and nodes, but it is general in the sense that it does not rely on the spatial structure of the brain.   The method is computationally efficient since the regularization we use is convex, and the solution is implemented with efficient optimization algorithms. These properties guarantee fast  convergence to the solution, making the methods scalable to networks with thousands of nodes, which is enough to deal with many of the brain atlases usually employed in neuroimaging (see for example <cit.>). Statistically, the rate of convergence depends on the number of active nodes only, not the total number of nodes, which allows for accurate results with even moderate sample sizes if the active node set is small.  


The results we obtained on the schizophrenia data are generally in agreement with previous studies.  In particular, the default mode network has been consistently implicated in schizophrenia and many other psychiatric disorders <cit.>.     While different subnetworks were implicated by the two different datasets, we are still able to predict the disease status fairly accurately by training on one dataset and testing on the other.   The differences between the two datasets may reflect real differences in samples collected at different sites and in different experiments, as significant pathophysiological heterogeneity occurs for all psychiatric diagnoses, or they may simply reflect type 2 errors.


Our methods work very generally with a sample of networks with labeled nodes and associated responses.   The many pre-processing steps inevitable when dealing with fMRI data always  add some uncertainty,  and pre-processing decisions can potentially affect downstream conclusions.    We aimed to somewhat mitigate this  by using ranks, which are more robust and showed a slightly better performance on our datasets.  Another option, when practical, is to compare multiple pre-processing pipelines, and/or multiple measures of connectivity, to further validate results.   Our method's independence of these particular choices and its computational efficiency make it an attractive option for such comparisons.  









§ ACKNOWLEDGMENTS


This research was supported in part by NSF grant DMS-1521551, ONR grant N000141612910, and a Dana Foundation grant to E. Levina, NSF training grant DMS-1646108 support for D. Kessler,  as well as by computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor.   S. F. Taylor's research  is supported by the National Institute of Mental Health (R01MH064148, R21MH086701, R21MH101676), the Boledovich Schizophrenia Research Fund and University of Michigan Clinical Translational Science Award (UL1RR024986).



apalike







§ OPTIMIZATION ALGORITHM DETAILS

The optimization procedure for solving the penalized prediction problem introduced above consists in a proximal algorithm, and the steps are detailed in Algorithm <ref>. Each step requires to solve a further convex optimization problem via ADMM. The exact solution of the steps of this method is shown in Algorithm <ref>.










§ PROOFS




Here we prove the bounds on Frobenius norm error and probability of support selection in Proposition <ref>, following the framework of <cit.> based on geometrical decomposability. A penalty Ω is geometrically decomposable if it can be written as 
    Ω(B)=h_A(B)+h_I(B) + h_E^⊥(B)

for all B, with A,I closed convex sets, E a subspace, and h_C the support function on C defined as h_C(B)=sup{⟨ Y,B⟩| Y∈ C.} .

The proof proceeds in the following steps.  
Lemma <ref> shows that  an equivalent form of our penalty (<ref>) is geometrically decomposable, allowing us to use the framework of <cit.>.
Lemma <ref>  shows the Assumption <ref> together with a lower bound on ρ imply that the irrepresentability assumption of <cit.> holds.  Assumption <ref> is directly on the entries of the loss Hessian, which simplifies the very general form of the assumption in <cit.>.  Lemma <ref> gives a bound on the entries of the loss gradient under the sub-gaussianity assumption <ref>.   Lemma <ref> gives explicit bounds for the compatibility constants that appear on Theorem 1 of <cit.>.   Finally, we combine these results to prove Proposition <ref>.

Without loss of generality, to simplify notation we assume that 𝒢={1,…,G}, that is, the active subgraph is in the first G rows of the matrix. 

	
	The penalty (<ref>) can be written as geometrically decomposable.  


   We use an equivalent formulation of the penalty in which every coefficient is penalized only once. Let B',B”∈^N× N be matrices such that the upper triangular part of B” and the diagonals of B' and B” are zero. Define
	
    Ω̃(B',B”) = ∑_i=1^NB'_(i)_2 + ρB”_1,

	and E = {(B',B”)∈^N× 2N :   B' = B'^T, B”_ij=B'_ij,  for i<j and B”_ij=0 for i≥ j}.
	Denote by R the transformation from ^N× N to ^N× 2N that replicates entries appropriately, 
	
    (RB)_ij = {[       B_ij if 1≤ j≤ N;   B_i(j-N)   if j> N.;            ].

	Therefore, for any B∈^N× N, we can uniquely define RB = (B',B”) such that Ω(B)=Ω̃(B',B”). 
	We then show that Ω̃ is geometrically decomposable. Moreover, for any (B',B”)∈ E we can define R^-1, so the penalties Ω and Ω̃ on E are equivalent. Define the sets A,I⊂^N× 2N such that
	
    A =   {(B',B”)  :    max_i∈𝒢B_(i)'_2 ≤ 1, max_i∈𝒢^CB_(i)'_2=0, . 
       .  max |B”_ij| ≤ρ, B_ij” = 0,(i,j)∈(𝒢×𝒢)^C },
    
    		I =    {(B',B”)  :  max_i∈𝒢^CB_(i)'_2 ≤ 1, max_i∈𝒢B_(i)'_2=0,.
       .  max |B”_ij| ≤ρ,  B_ij” = 0,(i,j)∈𝒢×𝒢} .
		
	Letting ⟨ Y,(B',B”)⟩=Tr(Y'B'^T) + Tr(Y”B”^T), combining the arguments of <cit.> for lasso and group lasso penalties, 
	
    h_A(B',B”) =     ∑_i∈𝒢B'_(i)_2 + ρ∑_(i,j)∈𝒢×𝒢|B”_ij|,
    
    		h_I(B',B”) =    ∑_i∈𝒢^CB'_(i)_2 + ρ∑_(i,j)∈(𝒢×𝒢)^C|B”_ij|,
    
    		h_E(B',B”) =    {[              0 if (B',B”)∈ E,;              ∞     otherwise. ].

	Hence,  Ω can be written as a geometrically decomposable penalty 
	
    Ω(B) = Ω̃(B',B”) = λ{h_A(B',B”)+h_I(B',B”)+h_E(B',B”)}.


We introduce some notation in order to state the irrepresentability condition of <cit.>. For a set F⊂^N× 2N and Y∈^N× 2N, denote by γ_F(Y) = inf{λ>0 : Y∈ F} the gauge function on C. Thus,

    γ_I(B',B”) = max{max_i∈𝒢^CB_(i)'_2, 1/ρmax_(i,j)∈(𝒢×𝒢)^C|B_ij”| } + 1_I(B',B”),

where 1_I(B)= 0 if B ∈ I and ∞ otherwise. Define 

    V(Z)=inf{γ_I(Y) :   Z-Y∈ E^⊥,Y∈^N× 2N}

for Z∈^N× 2N. Let ℳ̃=E∩span(I)^⊥ be the set of matrices with correct support in the extended space ^N× 2N, similarly to ℳ in (<ref>).
Denote by P_M and P_M^⊥ the projections onto ℳ̃ and ℳ̃^⊥. Define the function ℋ(Z):^N× N→^N× N as 

    ℋ(Z)_ij = {[ H_(i,j),𝒢 (P_M Z)_𝒢,𝒢               if j∈𝒢,;                     0            otherwise. ].

where H_(i,j),𝒢 is the matrix defined in (<ref>). The Irrepresentability Assumption 3.2 of <cit.> requires the existence of 0<τ̃<1 such that

    sup_Z∈ A V{P_M^⊥(Rℋ(Z) - Z)} < 1 - τ̃.

For a support function h, denote by ∂ h(M)=⋃_Y∈ M∂ h(Y) the set of subdifferentials of h in M. Note that ∂ h_A(M) = A, since 0∈ M and ∂ h_A(0) = A. 


	If Assumption <ref> holds and ρ>1/τ-1/√(G), then there exists 0<τ̃<1 such that (<ref>) holds.
	


	Since V is sublinear (Lemma 3.3 of <cit.>),
	
    sup_Z∈ A V{P_M^⊥(Rℋ(Z) - Z)}≤sup_Z∈ A V{P_M^⊥(Rℋ(Z))} + sup_Z∈ A V{P_M^⊥Z}.

	To bound the first term, note that E^⊥={(Z',Z”)|Z_ij'+Z_ji'+Z_ij”=0, j < i}.	Hence,
	
    V(Y',Y”) =    inf{γ(U',U”):  U_ij'-Y_ij^(1)+ U_ji'-Y_ji' + U_ij”-Y_ij”=0,  j < i }
    ≤   inf{γ(U',U”):   U_(i)' = Y_(i)', i∈𝒢^C;  U_𝒢^C,𝒢^C” = Y_𝒢^C,𝒢^C”;.
           . (U',U”)- (Y', Y”)∈ E^⊥}
    ≤   max{max_i∈𝒢^CY_(i)'_2, 1/ρ  Y_𝒢^C,𝒢^C”_∞}.

	Therefore,
	
    V(P_M^⊥(Rℋ(Z)))  ≤   max{max_i∈𝒢^C(P_M^⊥(Rℋ(Z)))_(i)^(1)_2, 1/ρ(P_M^⊥(Rℋ(Z)))^(2)_𝒢^C,𝒢^C_∞}
    
    		=    max_i∈𝒢^Cℋ(Z)_(i)_2,

	which implies that
	
    sup_Z∈ A V(P_M^⊥(Rℋ(Z)))  ≤   sup_Z∈ A{max_i∈𝒢^Cℋ(Z)_(i)_2}
    ≤   sup_B∈^G× G,B_(i)_2≤ 1{max_i∈𝒢^C(H_(i,j),𝒢B)_j=1^N_2}
    ≤   max_i∈𝒢^C(∑_k=1^G(H_(i,j),𝒢)_k·_2)_j=1^N_2  =  1-τ.

	
	
	Let Z=(Z',Z”)∈ A. Without loss of generality, assume that Z'_𝒢,𝒢=Z”_𝒢,𝒢=0 (note that these entries do not change V(P_M^⊥Z)). Therefore, P_M^⊥Z=Z. Hence, 
	
    V(Z) =    inf{γ(U',U”):   (U',U”)∈ I, (U',U”)-(Z',Z”)∈ E }
    
    		=    inf{γ(U',U”):  U'_ij + U”_ij=Z'_ji, 1≤ j ≤ G, G< i ≤ N }
    
    		=    inf{max{max_i∈𝒢^CU'_(i)_2,1/ρmax_j ∈𝒢
     i∈𝒢^C|U”_ij|} :  U'_ij + U”_ij=Z'_ji, j ∈𝒢, i∈𝒢^C }
    ≤   inf{max{max_i∈𝒢^CU'_(i)_2,1/ρmax_j ∈𝒢
     i∈𝒢^C|U”_ij|} :  U'_ij + U”_ij=1, j ∈𝒢, i∈𝒢^C}

	The last bound from  |Z'_ji|≤ 1 and no longer depends on Z.
	It is easy to see that the minimum is attained when, for each i>G, 
	
    U'_(i)_2 = 1/ρ|U”_ij|,     1≤ j≤ G,

	and therefore 
	
    V(Z) ≤√(G)/1+ρ√(G).
	
	Moreover, if Z^∗∈ A is defined such that (Z^∗)^(1)_G+1,i=1 for i=1,…,G and 0 elsewhere, then V(Z^∗) achieves this bound, which shows that ρ > 1- 1/√(G) is a necessary condition for the irrepresentability to hold, even in the case where the entries of the Hessian that denote the information between active and inactive edges is zero. Therefore, plugging the bounds (<ref>) and (<ref>) into equation (<ref>), we obtain (<ref>) holds as long as 1-τ +  √(G)/1+ρ√(G)<1, which implies that ρ>1/τ-1/√(G).	





The next lemma establishes a bound on the dual norm of Ω of the loss gradient function under a sub-Gaussian assumption. Let Ω^∗ denote the dual norm of Ω, so  Ω^∗(B) =sup{⟨ Y,B⟩| Y∈ C, Ω(Y)≤ 1.}.  
 Under Assumption <ref>,
	
    P(Ω^∗(∇ℓ(B^⋆))>t) ≤ 2N^2 min{exp( - n(1+ρ)^2t^2/N(σ^2)), exp( - nρ^2t^2/σ^2)}.

	


	By Hoeffding's inequality for sub-Gaussian variables, for all j,k and t>0,
	
    (|∇_jkℓ(B^⋆)|>t) ≤(|1/n∑_i=1^n∇_jkℓ_i(B^⋆)|>t) ≤  2exp(-nt^2/σ^2).

	Note that (1+ρ)∑_i=1^NB_(i)_2 ≤Ω(B). Let Φ(B)=1/1+ρmax_i=1,… NB_(i)_2. Thus,
	
    Ω^∗(B)≤sup_Y∈^N× N{Tr(YB):Ω_ρ=0(Y)≤1/1+ρ}=Φ(B).

	In a similar manner, ρB_1≤Ω(B). Setting Ξ(B) = 1/ρB_∞, we have
	
    Ω^∗(B)≤Ξ(B).

	Using (<ref>) and setting Λ=∇ℓ(B^⋆),
	
    P{Ω^∗(Λ)>t}   ≤   P{Φ(Λ) > t}
        =    P{max_1≤ i ≤ NΛ_(i)_2>(1+ρ)t}
       ≤   P{max_1≤ i ≤ Nmax_j≠ i|Λ_ij|>(1+ρ)t/√(N)}
       ≤   2N(N-1)exp{-n(1+ρ)^2t^2/2σ^2(N-1)},

	the last inequality obtained by arguments similar to Lemma 4.3 of <cit.>. In the same way, we can also bound the previous quantity using (<ref>) by
	
    P(Ω^∗(Λ)>t)    ≤   P(Ξ(Λ) > t)
        =    P(Λ_(i)_∞>ρ t)
       ≤   N(N-1)exp(-nρ^2t^2/2σ^2).

	Combining (<ref>) and (<ref>), we can obtain equation (<ref>). 





For a semi-norm Ψ:^N× N→, let κ_Ω be the compatibility constant between Ψ and the Frobenius norm, defined as

    κ_Ψ = sup{Ψ(B):B_2≤ 1, B∈ M},

and let κ_IC be the compatibility constant between the irrepresentable term and the dual norm Ω^∗ given by

    κ_IC = sup{V(P_M^⊥(RℋZ-Z)) : Ω^∗(Z)≤ 1 } .


 
	The following bounds on the compatibility constants hold: 
	
    κ_Ω    =    √(G) + ρ√(G(G-1)), 
    κ_Ω^∗   ≤   1/1+ρ,
    κ_IC   ≤    3 - τ.
	


 Note that Ω(Y) is maximized on {Y:Y_2≤ 1} when all entries of Y have magnitude equal to 1/√(G(G-1)).   Therefore 
	
    κ_Ω = G√(G-1/G(G-1)) + ρG(G-1)/√(G(G-1)) = √(G) + ρ√(G(G-1)).

	Similarly,  (<ref>) implies
	
    κ_Ω^∗≤sup{1/1+ρmax_i∈𝒢B_(i)_2:B_2≤ 1 }≤1/1+ρ.

	Finally, 
	
    κ_IC  =    sup{V(P_M^⊥(RℋZ-Z)) : Ω^∗(Z)≤ 1 }
    ≤   sup{V(P_M^⊥(RℋZ)) :  Ω^∗(Z)≤ 1 } + sup{V(P_M^⊥(Z)) :  Ω^∗(Z)≤ 1 }
    ≤    (1 - τ ) + 2 = 3-τ.

	






	Part (a). Since B̂ minimizes the problem (<ref>), 
	
    ℓ(B̂) + λΩ(B̂) ≤ℓ(B^⋆) + λΩ(B^⋆).

	Rearranging the terms, using Assumption <ref>, by the triangle inequality and the generalized Cauchy-Schwarz inequality,
	
    0  ≥   ℓ(B̂) - ℓ(B^⋆) + Ω(B̂) - Ω(B^⋆) 
    ≥   ⟨∇ℓ(B^⋆)^T, B̂-B^⋆⟩ + m/2B̂-B^⋆_2^2 - Ω(B̂ - B^⋆)
    ≥    -Ω(B̂ - B^⋆)Ω^∗(∇ℓ(B^⋆))- Ω(B̂ - B^⋆) + m/2B̂-B^⋆_2^2.

	By the argument for computing κ_Ω in (<ref>), Ω(Y)≤{√(N) + ρ√(N(N-1))}Y_2. Rearranging the terms in (<ref>),
	
    B̂-B^⋆_2≤2/m{√(N) + ρ√(N(N-1))}{λ + Ω^∗(B̂ - B^⋆)}.

	For any ρ, setting λ = 2 √(σ^2log N/n)min{√(N)/1+ρ, 1/ρ}, by Lemma <ref>, with probability at least 1-2/N, 
	
    B̂-B^⋆_2  ≤   4/m{√(N) + ρ√(N(N-1))}λ
    ≤   4/m√(σ^2log N/n){√(N) + ρ N}min{√(N)/1+ρ, 1/ρ}
    ≤   4/m√(σ^2log N/n) Nmin{1+ρ√(N), 1 +  1/ρ√(N)}
    ≤   4/m√(N^2σ^2log N/n).

	
	
	Part (b). Lemma <ref> gives a geometric decomposition of the penalty. Therefore, we can directly use Theorem 3.1 of <cit.>, since Lemma <ref> also ensures that their irrepresentability condition holds. Thus,
	
    B̂-B^⋆_2≤2/mκ_Ω(1 + τ/4κ_IC)λ,

	and 𝒢̂⊆𝒢 as long as 
	
    4κ_IC/τΩ^∗(∇ℓ(B^⋆)) < λ < m^2τ/2Lκ_Ω^2κ_Ω^∗κ_IC(1+τ/4κ_IC)^-2.

	Setting 
	
    λ = 8κ_IC/τ√(σ^2log N/n)min{√(N)/1+ρ, 1/ρ},
 
	using a similar argument than (<ref>), the left hand size of  (<ref>) holds with probability at least 1-2/N. The right hand side of (<ref>) holds as long as the sample size satisfies
	
    n > C(L, m, τ, κ_Ω, κ_Ω^∗, κ_IC) (√(G) + G)^2σ^2log N,

	with C(L, m, τ, κ_Ω, κ_Ω^∗, κ_IC)>0 a positive constant. Therefore, claims (<ref>) and (<ref>) follow. 
	





§ DATA AQUISITION AND PRE-PROCESSING





 §.§ Subjects and imaging




  §.§.§ The COBRE data


Raw anatomic and functional scans from 146 subjects (72 psychosis patients and 74 healthy control subjects) were downloaded from a public database (<http://fcon_1000.projects.nitrc.org/indi/retro/cobre.html>). Four subjects coded as ambidextrous (2 patients, 2 controls) were excluded to yield 70 psychosis patients and 72 controls for analysis. To enter the COBRE dataset, subjects had a diagnosis of either schizophrenia or schizoaffective disorder and were without histories of neurological disorder, mental retardation, severe head trauma with more than 5 minutes loss of consciousness and substance abuse/dependence within the last 12 months.

In the primary sample, two schizophrenic (SCZ) subjects and one healthy control (HC) subject had insufficient voxels in the cerebrospinal fluid (CSF) segmentation on the CSF, and they were dropped from additional analyses. Two additional SCZ subjects were excluded for scrub ratios (see discussion of scrubbing routine in fMRI Data Analysis) greater than 0.6, leaving 38 SCZ subjects and 42 HC subjects for the analysis. In the replication sample, 15 psychosis patients and two control subjects were excluded for scrub ratios greater than 0.6; one patient was excluded with incomplete data,  leaving 54 SCZ and 69 HC subjects for analysis (see Table <ref>). 

A full description of the imaging parameters for the COBRE dataset is available online at the link provided above and in  <cit.>.



  §.§.§ The UMich data


Subjects were selected from experiments conducted by Professor Stephan F. Taylor at the University of Michigan between 2004 and 2011 for task-based fMRI studies, which included resting state scans. Thirty-nine stable outpatients were selected with DSM-IV schizophrenia or schizoaffective disorder (SCZ) <cit.>. Forty healthy comparison (HC) subjects, without a lifetime history of Axis I psychiatric disorders <cit.>, were selected to approximate the age range, gender distribution and family education level of the patients. Prior to initial data collection, all subjects gave written, informed consent to participate in the protocol approved by the University of Michigan institutional review board (IRBMED).



MRI scanning occurred on a GE 3T Signa scanner (LX [8.3] release, General Electric Healthcare, Buckinghamshire, United Kingdom). Functional images were acquired with a T2*-weighted, reverse spiral acquisition sequence (gradient recalled echo, TE=30 msec, FA=90 degrees, field of view=22 cm, 40 slices, 3.0mm thick/0mm skip, equivalent to 64 x 64 voxel grid – yielding isotropic voxels 3 mm on edge). Because the data were acquired across different experiments, acquisition parameters differed slightly in the aggregate sample: 240 volumes at TR=1500 msec (11 SCZ, 10 HC), 180 volumes at TR=2000 msec (17 SCZ, 16 HC) and 240 volumes at 2000 msec (14 SCZ, 17 HC). Acquisitions were acquired in the resting state with eyes open and fixated on a large ‘plus’ sign projected on a monitor.








 §.§ Pre-processing


We first performed standard pre-processing steps. All scans were slice-time corrected and realigned to the 10th image acquired during a scanning session <cit.>. Subsequent processing was performed with the Statistical Parametric Mapping SPM8 package (Wellcome Institute of Cognitive Neurology, London). Anatomic normalization was done with the VBM8 toolbox in SPM8, using the high resolution structural scans obtained for both datasets. Normalizing warps were applied to the co-registered, functional volumes, which were re-sliced and smoothed with an 8 mm isotropic Gaussian smoothing kernel. To assess and manage movement, we calculated the frame-wise displacement (FD) <cit.>, for all 6 parameters of rotation and translation. We used a scrubbing routine to censor any frame with FD > 0.5 mm from the regression analysis described below, yielding a scrub ratio for each subject. Three-compartment segmentation of the high-resolution structural image from the VBM8 normalization was applied to the functional time series to extract cerebral spinal volume (CSF) and white matter (WM) compartments, which were then subjected to a principal component analysis to identify the top 5 components in each <cit.>, which should correspond to heart rate and respiratory effects on global signal <cit.>. Multiple regressions were applied to the time series to remove the following nuisance effects: linear trend, 6 motion parameters, their temporal derivatives, the quadratics of these 12 parameters, 5 components from the PCA of CSF, 5 components of PCA of WM, followed by band pass filtering from 0.01 – 0.1 Hz, and then motion scrubbing. For each 4D data set, time courses were then extracted from 10 mm diameter spheres based on the 264 sets of coordinates identified by <cit.>. From these time series, a cross-correlation matrix of Pearson r-values was obtained and Fisher's R-to-Z transformation was applied for each of the 264 nodes with every other node (for COBRE dataset, node 75 is missing). Finally, for each individual, edge weights were assigned to be ranks of these score, with edge scores ranked separately for each subject, and then these values were centered and standardized across the individuals. Ranks have been used previously in brain connectomic studies to reduce the effect of potential outliers <cit.>; although some information is lost with the rank transformation, we observed that while ranks do not increase the classification accuracy significantly, they tend to produce sparser solutions with a similar accuracy to Pearson correlations. 


