Computer Science
Computer Science-2021.gz:143208
Computer Science
aae7a9c2c8a7ac351295026cca85251a87330fa3
Learning Interpretable Models for Coupled Networks Under Domain Constraints

Modeling the behavior of coupled networks is challenging due to their intricate dynamics. For example in neuroscience, it is of critical importance to understand the relationship between the functional neural processes and anatomical connectivities. Modern neuroimaging techniques allow us to separately measure functional connectivity through fMRI imaging and the underlying white matter wiring through diffusion imaging. Previous studies have shown that structural edges in brain networks improve the inference of functional edges and vice versa. In this paper, we investigate the idea of coupled networks through an optimization framework by focusing on interactions between structural edges and functional edges of brain networks. We consider both types of edges as observed instances of random variables that represent different underlying network processes. The proposed framework does not depend on Gaussian assumptions and achieves a more robust performance on general data compared with existing approaches. To incorporate existing domain knowledge into such studies, we propose a novel formulation to place hard network constraints on the noise term while estimating interactions. This not only leads to a cleaner way of applying network constraints but also provides a more scalable solution when network connectivity is sparse. We validate our method on multishell diffusion and task-evoked fMRI datasets from the Human Connectome Project, leading to both important insights on structural backbones that support various types of task activities as well as general solutions to the study of coupled networks.


Introduction
To achieve the above goal, we consider an edge-by-edge formulation. We treat one set of edges in G 1 as predictors and the other set of edges in G 2 as response variables in a multivariate linear regression model. Past research for the above problem relies on the restrictive Gaussian assumption, which simplifies the problem but is difficult to justify, especially in the domain of brain architectures. Adopting Gaussian assumption on non-Gaussian data can significantly prevent the detection of conditional dependencies and may lead to incorrectly inferred relationships among variables.
The learning of relationships between two different modalities can be difficult without sufficient data. As a result, in sparser data settings, the ability to specify constraints based on domain knowledge can be beneficial. For example, in the case of brain data, functional edges have mainly local influences, and structural edges are more responsible for long-distance influences (Rubinov and Sporns 2010;Batista-García-Ramó and Fernández-Verdecia 2018). We want preferences encoded in domain knowledge to guide the selection of partial correlations of unexplained noise terms in the constructed model. Based on the above motivations, we propose a flexible and efficient framework CC-MRCE (Convex-set Constrained Multivariate Regression with Covariance Estimation) that simultaneously learns both regression coefficients between two coupled networks and the correlation structure of noise terms. In a departure from existing methods, our framework encodes domain knowledge as a set of convex constraints and adopts a pseudolikelihood-based neighborhoodselection objective in partial correlation estimation, which has been shown to be more robust to non-Gaussian data. Because of the CC-MRCE objective's bi-convex nature, we alternately solve a regression sub-problem and a constrained partial correlation sub-problem until convergence. The latter sub-problem requires feasible solutions under given domain constraints that we render tractable via a modified two-stage proximal gradient descent method.
We illustrate the use of our method in the context of the human brain. Brain data presents one of the greatest technical challenges in analysis and modeling due to a network-based characterization (Bassett and Sporns 2017;Bassett, Zurn, and Gold 2018), non-Gaussian nature of data, high dimensionality, a small number of samples, and the need to incorporate domain knowledge. We apply the pro-posed framework on the Human Connectome Project (HCP) dataset (Van Essen et al. 2013), where two coupled networks are constructed from fMRI scans (representing cortical activation) and diffusion scans (representing the anatomical substrate). We successfully predict a brain functional network from the given structural network; our method outperforms previous state-of-art methods, and our obtained models are easier to interpret. We investigate the structurefunction coupling for seven different tasks. Our findings agree with the nature of fMRI task and brain region functions in existing literature, thus validating our model's ability to discover meaningful couplings.
Our main contributions are as follows: • We propose a regularized multiple regression approach that adapts to non-Gaussian data. (Sections 2.1,2.2) • We incorporate prior domain knowledge to model estimation by formulating constraints into an optimization problem. (Section 2.3) • We develop a fast method based on nested FISTA for solving the proposed optimization problem. (Section 3) • We show the effectiveness of our model on HCP brain data using quantitative comparisons with existing approaches as well as a qualitative analysis. (Section 4.2)

Problem Formulation
In this section, we first introduce existing works on multiple regression under Gaussian assumptions and then motivate our approach under non-Gaussian settings and domain constraints.

Multiple-output Regression Problem
Let D be an n-subject sample set in which all subjects share the same coupling (G 1 , G 2 ) but have different edge values. For subject i in D, let x (i) = (x (i) 1 , · · · , x (i) p ) be p-dimensional inputs that represent edge values in the first modality network G 1 , and y (i) = (y (i) 1 , · · · , y (i) p ) be pdimensional outputs 1 that stand for edge values in the second modality network G 2 . We assume that the inputs x i and outputs y i are correlated through a multivariate linear regression model: where B is the p × p regression coefficient matrix and its element β jk is the regression coefficient that measures the cross-modality impact of edge x j to edge y k , and (i) is the noise vector of subject i. The model can be expressed in the matrix form: where row i of X ∈ R n×p and Y ∈ R n×p are the structural and functional edge vectors x (i) and y (i) of subject i. A straightforward approach to estimating B is to solve p separate regression problems, assuming noise terms are independent and uncorrelated. Recently, advanced methods have been proposed to exploit the correlation in noise terms to improve the modeling. They accomplish the goal by introducing an assumption that noise terms (1) ,..., (n) are all i.i.d. Gaussian N (0, Ω −1 ) and then simultaneously estimating regression coefficients B and inverse covariance matrix Ω of the noise terms. Two popular methods along this direction are MRCE (Rothman, Levina, and Zhu 2010) and CGGM (Sohn and Kim 2012;Wytock and Kolter 2013;Yuan and Zhang 2014). The MRCE method considers the conditional distribution Y |X ∼ N (XB, Ω −1 ) and estimates both B and Ω by alternately minimizing the negative conditional Gaussian likelihood, with the 1 lasso penalty applied on the entries of B and Ω. The other method, CGGM, further assumes that X and Y are jointly Gaussian. Under such formulation, the conditional distribution of Y |X is given by N (−XΩ XY Ω −1 , Ω −1 ), which reparameterizes the regression coefficient B as −Ω XY Ω −1 . Compared with MRCE, the objective of CGGM is based on the negative conditional Gaussian likelihood as well, but is jointly convex for Ω XY and Ω, and therefore more friendly to computation.

Relaxing Gaussian Assumptions
Although MRCE and CGGM have received significant attention in solving multi-output regression problems, one drawback of these two approaches is the Gaussian assumption, especially for applications to brain data (Freyer et al. 2009;Hlinka et al. 2011;Eklund, Nichols, and Knutsson 2016). Recall that MRCE assumes the Gaussian noise, and CGGM further assumes joint Gaussian distribution over both inputs and outputs. We tested whether the HCP structural and functional data is Gaussian with a significance level of 0.05. The test rejects the Gaussian null hypothesis for 97.5% of structural edges and 36.3% of functional ones. Since our sample size is small, false negatives are more likely to occur (Columb and Atkinson 2016), namely failing to reject the Gaussian hypothesis when the underlying data is non-Gaussian. Therefore, the proportion of non-Gaussian data in brain networks is expected to be even higher. Thus, relying on Gaussian assumptions is likely to affect the constructed models negatively. To avoid a Gaussian assumption, we propose a pseudolikelihood approach for learning multi-output regression models by optimizing the following objective function: or in a neat matrix notion: where Ω = {ω jk } denotes the inverse covariance matrix, B = {β jk } denotes the coefficient matrix, and Ω D and Ω X denote the diagonal and off-diagonal parts of Ω. The proposed objective can be considered as a reparameterization of the Gaussian likelihood with Ω 2 and an approximation to the log-determinant term. It has been proven that under mild singularity conditions, such reparameterization can guarantee estimation consistency for distributions with sub-Gaussian tails (Peng et al. 2009;Khare, Oh, and Rajaratnam 2015). Next, we develop an optimization algorithm to minimize the objective. The objective function itself is not jointly convex for both variables B and X, but remains convex with respect to each of them while keeping the other fixed. Therefore, we adopt the alternating minimization idea. In the t-th iteration, we first fix B as the estimatedB (t−1) from the previous (t − 1)-th iteration, and calculate the empirical covariance matrix S of noise terms: Next, we estimate the inverse covariance matrix: Ω (t) with the given S (t−1) as a constant: Observe that the above subproblem follows CONCORD's original form, which is more robust to heavy-tailed data (Khare, Oh, and Rajaratnam 2015;Koanantakool et al. 2018) than the conventional Gaussian likelihood approach and can be efficiently solved using proximal gradient methods with a convergence rate of O(1/t 2 ) (Oh et al. 2014). Lastly, we keep Ω fixed at Ω (t) 2 and optimize the regression coefficients B: Note that subproblem (6) is convex when Ω (t) 2 is positive semi-definite. We present the above regression-based approach as CONCORD-MRCE in Algorithm 1 (pseudocode, see appendix).

Imposing Domain Constraints
Due to limited sample size of real-world datasets and their high dimensionality, incorporating accurate domain constraints can reduce the search space and avoid over-fitting. Under a linear mapping assumption, the partial correlation of response variables arises only from correlations in the noise terms. Therefore, Ω not only represents the inverse covariance of noise terms, but also equals the conditional inverse covariance of Y |X. The nonzero entries of Ω encode direct relationships among the target modality outputs Y that cannot be explained by weighted inputs XB of the source modality. It will be beneficial if the zero-vs-nonzero structure of Ω is partially given by domain experts and used as hard constraints in model estimation.
More formally, let M be a binary matrix that has the same dimensions as Ω. We can define a convex matrix set S M , containing all matrices that share the same set of zero entries with M . We can then improve the previous regressionbased approach to estimate Ω under the domain constraint that takes the form of Ω ∈ S M , written in an equivalent unconstrained convex form: where I{Ω ∈ S M } is an indicator function. This formulation can be extended to Ω ∈ C whenever C is a closed convex set of positive definite matrices.

Solution
In this section, we show how to adapt the previous solution when the inverse covariance Ω is constrained during estimation. Notice that the ideas of Algorithm 1 can be mostly used to solve Eq (7), except for the Ω-update step, which is now affected by the added constraints. The new Ω-update step needs to solve the following sub-problem: We follow the FISTA (Fast Iterative Soft-Thresholding Algorithm (Beck and Teboulle 2009)) approach that is used in CONCORD (Oh et al. 2014). This method utilizes an accelerated gradient algorithm using soft-thresholding as its proximal operator for the L 1 norm and achieves a fast O(1/t 2 ) convergence rate. Previous work (Oh et al. 2014) has also applied FISTA for partial correlation estimation and proved its efficiency. To adapt our constrained problem into the FISTA framework, we split our objective function (8) into a smooth part and a non-smooth part: For any symmetric matrix Ω, the gradient of the smooth function can be easily calculated as: ∇h 1 (Ω) = −2Ω −1 D + 2ΩS. With this formulation, we now adapt the FISTA iterative scheme to solve our network-constrained problem (8): where τ t is the step length and t denotes the iteration number. γ is a trade-off parameter that controls the extent to which the proximal operator maps points towards the minimum of h 2 (Ω), with larger values of γ associated with larger movement near the minimum. In these iterative steps, Θ t+1 is an expected position, updated purely by momentum. Within each loop, the algorithm first takes a gradient step of the estimated future position (Eq.10) and then applies the proximal mapping of a closed convex function h 2 (Ω).
In contrast to the standard FISTA approach, the composite function h 2 (Ω) consists of a sparsity penalty and a networkconstrained indicator function. More specifically, we can write down the explicit form of Eq. (11) according to the proximal operator definition: where A = Θ (t+1) − (nτ t /2)∇h(Θ (t+1) ). Instead of directly solving the original problem (12), we consider its dual problem as follows. Let matrix H be the dual variable of matrix Ω. We have: min where A D and A X denote the diagonal and off-diagonal part of A. Since the initial objective function above is convex in Ω X and concave in H X , we exchange the order of the minimum and maximum operator in which the inner minimization problem has an obvious solution through orthogonal projection theorem (Rockafellar 1970), written as where P C is defined as an projection operator: P C (Γ) = argmin R∈C R − Γ 2 F and its orthogonal projection operator P C ⊥ is defined as I − P C . In the special case that C = S M , projection P C (Γ) is equivalent to removing invalid nonzero entries of the input matrix Γ.
Inserting the optimal Ω X back into objective (13), we now obtain the final dual form of the problem (12): where any solutionĤ X to the dual problem corresponds to a primal solution through Eq.(14). Since the dual objective is continuously differentiable and constraints on l ∞norm are convex, we can again efficiently solve it with additional inner FISTA iterations, which is to minimize an equivalent composite objective min g 1 (H X ) + g 2 (H X ), where g 1 (H X ) is smooth and g 2 (H X ) is non-smooth: To adapt FISTA to the problem (15), the only thing left is to obtain the gradient of smooth function g 1 (H X ), for which we need the lemma below: Lemma 1. (Moreau 1965) If g is a closed proper convex function, and for any positive t, define a proximal opera- According to Lemma 1, we simply plug g(u) = δ(u ∈ C), and obtain that prox t (g)(x) = P C (x). Therefore the gradient ∇g 1 is calculated as: and the proximal mapping of function g 2 (Ω) becomes a projection of H X into the L ∞ -ball: This completes the modified Ω-update step in the constrained setting (referred to as Constrained-CONCORD). Its corresponding pseudocode is presented in Algorithm 3 (pseudocode, see appendix). The overall framework of simultaneously estimating Ω and B is presented in Algorithm 2 (pseudocode, see appendix), which we name as CC-MRCE.

Experiments
We present two sets of experiments. First, to understand our method's strengths and limitations, we utilize simulated datasets that allow us to inspect both reconstruction performance and model selection performance. We compare the proposed method CC-MRCE with other baselines when the underlying data distribution does not follow the Gaussian assumption. We show that both the non-Gaussian assumption and the network constraints contribute to the improvement of performance. Second, we conduct experiments on the Human Connectome Project (HCP) data (Van Essen et al. 2013). Our model offers a quantitative advantage over baseline methods for predicting functional networks from structural ones; at the same time, our results agree with existing neuroscience literature.

Application to Simulated Data
Data Generation. Using a similar approach to existing works (Peng et al. 2009;Khare, Oh, and Rajaratnam 2015;Rothman, Levina, and Zhu 2010), we generate our simulated dataset by first synthesizing two key model parameter matrices Ω 0 and B 0 , and then construct input, output, and noise terms (i.e. x (i) , y (i) and (i) 's).
Note that every positive-definite matrix has a Cholesky decomposition that takes the form of LL T , where L is a lower triangular matrix L, and if L is sparse enough then LL T is sparse as well. Therefore, we first sample a sparse lower triangular L with real and positive diagonal entries, and then generate our inverse covariance matrix with Ω 0 = LL T . The generated p×p positive definite matrix Ω 0 has 10% nonzeros entries and a condition number of 4.3. To demonstrate the robustness of proposed method CC-MRCE on non-Gaussian data, we sample the noise terms according to a multivariate t-distribution with zero mean and covariance matrix Σ = Ω −1 0 . Next, we generate a sparse coefficient matrix B 0 using the matrix element-wise product trick, In this construction approach, W has entries with independent draws from standard normal distribution N (0, 1). Each entry in K is drawn independently from a Bernoulli distribution that takes value 1 with probability s 1 . Q has rows that are either all one or all zero, which are determined by independent Bernoulli draws with success probability s 2 . Generating the sparse B 0 in this manner, we not only control its sparsity level, but also forcibly make (1 − s 2 )p predictors to be irrelevant for p responses, and guarantee that each relevant predictor is associated with s 1 p response variables.
In the following experiments, the probabilities s 1 and s 2 are chosen to be 0.15 and 0.8, the sample size n is fixed at 50, and the input and output dimensions p and q are both set to 20. The input x (i) 's are sampled from a multivariate normal distribution N (0, Σ X ) where (Σ X ) jk = 0.5 |j−k| , following previous works (Yuan and Lin 2007;Peng et al. 2010). The output y (i) 's are calculated as the linear model assumption in Eq.(1). We replicate the above process for independently generating a validation dataset of the same sample size. All penalty parameters are selected simultaneously and tuned according to the validation error.
Methods. In this heavy-tailed setting, we compare the performance of CC-MRCE to CGGM and MRCE, which are developed under Gaussian settings. The CGGM implementation that we used in this experiment is provided by (McCarter and Kim 2016) and has been optimized for largescale problems and limited memory. The MRCE implementation is provided by (Rothman, Levina, and Zhu 2010). Both CGGM and MRCE do not adapt to the network constraints we impose here.
Performance Evaluation. We evaluate the reconstruction performance of models using the conventional MSE error (in percentage). Correlation coefficients between predicted and ground-truth outputs are also provided along with their corresponding p-values. In addition, we use the relative area under the curve (AUC) of receiver operating characteristic (ROC) curves (Fawcett 2006;Friedman, Hastie, and Tibshirani 2010), with regards to Ω and B, as key measures to compare model selection performance of all these methods. The AUC of a perfect ROC curve, which would be 1, indicates an ideal recovery of gound truth zero-vs-nonzero structure in Ω (or B). However, models with large falsepositive rates (FPR) are barely meaningful in real scenarios. To focus on the initial portion of ROC curves, we control the FPRs simultaneously to be smaller than 0.2 for both Ω and B estimation. Thus, the maximum AUC value that a model can reach is just 0.2. For ease of comparison, we provide relative AUC values, divided by 0.2 to normalize to 1. For each method, we run the algorithm with at least 25 appropriate parameter pairs (λ 1 , λ 2 ) to get its ROC curve. Recall that all methods in this section are required to estimate 800 parameters given n = 50 samples. Test Results. Table 1 displays the test results of six different settings in all measures mentioned above. As can be seen, all four CC-MRCE variations (with different network constraints) obtain significantly smaller reconstruction MSE percentages, higher correlation coefficients, and smaller pvalues. Note that CGGM behaves the worst in all measures since it is deeply rooted in the Gaussian setting and is consequently misled by these assumptions. We also see that the performance of CC-MRCE gradually improves when the applied network constraints are more informative.
We also plot ROC curves for Ω and B estimation in Figure 1.(a) and Figure 1.(b), respectively. It is clear that CC-MRCE performs better than MRCE and CGGM, across different choices of network constraints. For Ω estimation, as expected, CC-MRCE with more strict constraints has steeper curves, suggesting that it recovers mostly correct partial relationships between variables with very few spurious connections, and therefore achieves higher AUC scores. CC-MRCE (perfect) behaves perfectly for Ω-ROC, by its definition. For the estimation of matrix B, a similar phe-nomena demonstrates that network constraints improve the learning of regression coefficients and lead to a better reconstruction performance. Without imposing network constraints, unconstrained formulation of CC-MRCE is likely to generate a biased estimate ofΩ on small datasets and can not recover ground truth features in B.

Application to Human Connectome Data
Problem Formulation. Many works in literature report on coupling between brain structural connectivities (SCs) and functional connectivities (FCs) for both resting state (Hagmann et al. 2008;Hermundstad et al. 2013;Honey et al. 2009Honey et al. , 2007 and task-evoked states (Hermundstad et al. 2013;Raichle 2015;Cole et al. 2014;Davison et al. 2015). One such recent work (Becker et al. 2018) maps SC to resting-state FC by aligning both the eigenvalues and eigenvectors of a subject's SC and FC matrices, and evaluates the mapping by reconstructing resting FCs from SCs.
Inspired by the domain observation (Batista-García-Ramó and Fernández-Verdecia 2018) that two functional connections sharing the same node are more likely to form a meaningful pathway or functional activation pattern, we constrain the partial correlations between non neighboring edges to be zero, i.e. Ω jk = 0 if e j and e k are not incident edges in the fMRI-constructed network.
We conducted experiments by reconstructing task FCs from SCs of 51 subjects with their fMRI data from HCP dataset under seven task states. We used the parcellation scheme in (Hagmann et al. 2008) with a spatial scale of 33, resulting in 83 brain regions and 3403 possible edges. Columns in predictor matrix, Equation (2), X ∈ R 51×3403 represent SCs, whose entries are numbers of white matter streamlines intersecting pairs of brain regions, and columns in response matrix Y ∈ R 51×3403 represent FCs, whose entries are functional correlations between cortical activities of brain regions. B denotes the mapping (or coupling) between SCs and FCs, and E denotes the part of FCs that cannot be explained by SCs. We also normalized SC values to (0, 1], a range comparable to FC. We ran a 10-fold crossvalidation of our model for each task, splitting data in a 9-1 train-validation ratio (46-5 split for the 51 subjects). We selected the optimal hyperparameters λ 1 and λ 2 by a 5 × 5 grid search in the log-scale between 10 −1.6 to 10 −0.4 , keeping the models with smallest Mean Squared Error (MSE) percentage on the validation sets, averaged across 10 folds. In our case, both λ 1 and λ 2 have optimal values around 0.1 across tasks. Aside from MSE, we also tested the Pearson correlation coefficient between predicted FCs and ground truth FCs, (referred to as Pearson's r-score, listed in Table 2) and minimum, maximum p-values. The reconstruction MSE percentage is below 1% for the training data and around 8% for the validation data. Strong and significant positive correlations are shown for both training (r-score around 0.6 to 0.8) and validation (r-score around 0.5) data. These results indicate our model's effectiveness in FC reconstruction by exploiting cross-modal coupling between SC-FC and domain prior knowledge on FC-FC relationships.
Performance Evaluation. Regarding the ability of finding accurate mappings between SCs and FCs, we com-   (Becker et al. 2018). For CGGM, the time required to run a 10-fold cross-validation with one set of hyperparameters on a single task ranges from one day to a week. So we ran three sets of hyperparameters on EMOTION and LAN-GUAGE, selected the best parameter pair, and fixed it for the rest of the five tasks. In particular, we chose penalty terms on B and E to be both 0.01. For VAE, both encoder and decoder consist of two fully connected layers, with latent variable dimensions being two. We use MSE as the training loss.
In our experiments, increasing the number of layers or latent dimensions of VAE did not improve the final performance. For the Spectral Mapping method, we follow the setup of the original paper, setting maximum path length k to seven. The 10-fold cross validation results of CGGM, VAE and Spectral Mapping method are reported in Table 2.
The assumption of Gaussian noise weakens CGGM's performance on the HCP dataset: it has unstable and large average MSE percentages across tasks, and the correlations between predictions and ground truth are very small and even negative, which are also not statistically significant with regards to p-values. On the other hand, VAE models have stable MSE percentage (around 80%) and average Pearson's r-scores with small standard deviations across all tasks. The correlations of VAE models are weak (all around 0.08), yet statistically significant with p-values constantly smaller than 0.0025 for all tasks. This shows VAE learns a slightly meaningful mapping, but with such a small sample size, deep learning models are unlikely to perform well. Lastly, the Spectral Mapping method is designed for maximizing the correlation between fMRI prediction and ground truth for brain data, so it performs well as for correlation, however the prediction values are off, resulting in high MSE percentages. In all, our regression-based model performs better in both correlation and value reconstruction, showing its superiority in prediction on non-Gaussian data with small sample sizes.
Result Interpretation. Apart from better reconstruction performance, our model also has the advantage of result interpretability. We can explore the SC-FC mapping through the resulting coefficient matrix Bs. Since our problem definition is FC = SC · B + E, the i th row in B corresponds to i th edge pair in the SC vector. In the following, we say an edge i exists if row i of B has nonzero entries. As the experiment is run under the 10-fold cross validation setting and each training partition may generate a different mapping B, we consider an entry in the common B to be nonzero if it is nonzero more than five times in these 10 trials. The results are shown in Figure 2. From the figure, we can see for every task, several rows in B have many more nonzero entries than the others. This indicates the existence of several significant structural edges being responsible for most of the functional activities. To test if this assumption is valid, we compared Bs from our model to randomly generated Bs with the same levels of sparsity. The resulting MSE percentages, although having large variances, often have smaller means than that of CGGM and VAE, implying the importance of sparsity level of the coupling. However, the resulting r-scores of predicted FCs using a random B is the lowest among all methods. Together with very large p-values, the results predicted by random coupling show no correlation between predictions and ground truth. This indicates that our models learned meaningful mapping information from SCs to FCs, and that structured sparsity of B is important for getting predictions besides the level of sparsity alone. We now analyze the B matrices for different cognitive tasks. During fMRI data acquisitions of all seven tasks, participants are presented with visual cues, either as images or videos, and they need to use motions such as pressing buttons to complete the tasks . Interestingly, apart from the LANGUAGE task, the mappings learned by our model predicts the strongest contribution of left precentral and left postcentral connection, which is on the motor cortex responsible for right-side body movement. From this, we assume most participants use their right hands to conduct the required finger movements for these tasks. All mappings also contain edges in the occipital lobe, complying with the visual nature of these tasks. The visualization of common edges that exist in all seven tasks is shown in Figure 3. This "backbone" roughly resembles the Default Mode Network (Greicius et al. 2009;Raichle 2015).
We then examined which structural edges contribute significantly to the functional activity under different tasks. For this, we plot the "high-contributing" edges in LANGUAGE, GAMBLING, RELATIONAL and WM tasks in Figure 4. An edge is considered as high-contributing if the number of nonzero entries of its corresponding row in B is more than half of the maximum number of nonzero entries of any row in B. From Figure 4, we notice although a common back- bone exists, structural connections in different brain regions are responsible for specific tasks (e.g. SCs in and around hippocampus area appear to be highly contributing to the WM FCs but not so for the LANGUAGE ones, which is consistent with the literature (Baddeley, Jarrold, and Vargha-Khadem 2011)). We also plot both entry-wise and edge-wise overlap ratio for the mapping of seven tasks in Figure 5. Another interesting phenomenon in Figure 2 is that the number of nonzero entries of B for LANGUAGE, GAM-BLING, RELATIONAL and WM are much larger than the other three tasks: EMOTION, MOTOR and SOCIAL, although the final model for each task have a similar level of prediction performance and similar hyperparameters. This is largely caused by the nature of non-overlapping Bs for EMOTION, MOTOR and SOCIAL tasks: their B overlap ratios are significantly smaller than the other four tasks as shown in Table 3. Here we define the overlap ratio as the number of nonzero entries in the common B (entry ij being nonzero if it's nonzero more than five times) over the number of nonzero entries in B ∪ = B 1 ∪ · · · ∪ B 10 with B k being the predicted B using the k th split in 10-fold cross validation. This is also why we omit these three tasks for Figure  4, as the predicted Bs are not stable across the population and only a few common connections show significant contributions. We assume this results from group heterogeneity when carrying out these tasks. Further studies with heterogeneous models for the population will be useful to verify this assumption.

Conclusions
In this paper, we proposed a regularized regression-based approach (CC-MRCE) for jointly learning linear models and partial correlations among variables under domain constraints. Motivated by the neuroscience application of pre- Figure 5: Entry-wise and edge-wise overlap ratio for the mapping of seven tasks. 5a considers entry-wise overlap of predicted B of different tasks. The value on position (task i, task j) is the entry-wise IoU (Intersection over Union) of task i's B and task j's B, i.e. number of nonzero entries in B i ∩ B j over number of nonzero entries in B i ∪ B j . 5b considers the the SC edges responsible for different tasks. An edge is considered to exist when its corresponding row in B has nonzero entries. The value in position (task i, task j) is the number of common SC edges of task i and task j over the number of SC edges of task j.
dicting functional brain activities from structural connections, the CC-MRCE method discards the Gaussian assumption and incorporates domain constraints into model estimation. We further developed a fast algorithm based on nested FISTA to solve the optimization problem. With synthetic data analysis, we demonstrated that both domain constraints and assumption of non-Gaussian data contribute to the performance improvement of CC-MRCE. Our experimental results on Human Connectome Project data show that CC-MRCE outperforms existing methods on prediction tasks and uncovers couplings that agree with existing neuroscience literature.  Table 3: Overlap ratios (%) of predicted B (SCs-FCs mapping) across 10 folds for seven tasks.
In every step of the while loop, CC-MRCE calls subfunction Constrained-CONCORD to estimate the partial correlation matrix Ω (t) under hard constraints, i.e. the solution Ω (t) shares the same zero and nonzero pattern as matrix E (shown in Algorithm 3). Such constraints can be replaced by any set of convex constraints on Ω (t) . Note that superscript t and t in Algorithm 3 are iteration counters of inner and outer stages in the Constrained-CONCORD algorithm, respectively. Lemma 2. Let L(g 1 ) be the Lipschitze constant of the gradient of objective function g 1 (H X ), then L(g 1 ) ≤ 2λ 2 2 γ 2 . Although we use line-search to pick a proper step length κ t in the inner-loop of Algorithm 3, it can be replaced with a constant step length κ t = 2λ 2 2 γ 2 according to the above lemma.