Differential Analysis of Directed Networks
Min Ren Dabao Zhang

Unified Perceptual Parsing for Scene Understanding
    Tete Xiao1*, Yingcheng Liu1*, Bolei Zhou2*, Yuning Jiang3, Jian Sun4
    December 30, 2023
========================================================================




We developed a novel statistical method to identify structural differences between networks characterized by structural equation models. We propose to reparameterize the model to separate the differential structures from common structures, and then design an algorithm with calibration and construction stages to identify these differential structures. The calibration stage serves to obtain consistent prediction by building the ℓ_2 regularized regression of each endogenous variables against pre-screened exogenous variables, correcting for potential endogeneity issue. The construction stage consistently selects and estimates both common and differential effects by undertaking ℓ_1 regularized regression of each endogenous variable against the predicts of other endogenous variables as well as its anchoring exogenous variables. Our method allows easy parallel computation at each stage. Theoretical results are obtained to establish non-asymptotic error bounds of predictions and estimates at both stages, as well as the consistency of identified common and differential effects. Our studies on synthetic data demonstrated that our proposed method performed much better than independently constructing the networks. A real data set is analyzed to illustrate the applicability of our method.




§  INTRODUCTION


It is of great importance and interest to detect sparse structural differences or differential structures between two cognate networks. For instance, the gene regulatory networks of diseased and healthy individuals may differ slightly from each other <cit.>, and identifying the subtle difference between them helps design specific drugs. Social networks evolve over times, and monitoring their abrupt changes may serve as surveillance to economic stability or disease epidemics <cit.>. However, addressing such practical problems demands differential analysis of large networks, calling for development of efficient statistical method to infer and compare complex structures from high dimensional data. In this paper, we focus on differential analysis of directed acyclic or even cyclic networks which can be described by structural equation models (SEMs).

Many efforts have been made towards construction of a single network via SEM. For example, both <cit.> and <cit.> employed genetic algorithms to search for the best SEM. Most recently, <cit.> employed a hierarchical Bayes approach to construct SEM-based networks. However, these approaches were designed for small or medium scale networks. For large-scale networks whose number of endogenous variables p exceeds the sample size n, <cit.> proposed a regularization approach to fit a sparse model. Because this method suffers from incapability of parallel computation, it may not be feasible for large networks. <cit.> proposed another penalization approach to fit the model in a node-wise fashion which alleviates the computational burden. Most recently, <cit.>, <cit.>, and <cit.> each proposed a two-stage approach to construct SEMs, with different algorithms designed at different stages. As shown by <cit.>, such a two-stage approach can have superior performance compared to other methods.

To the best of our knowledge, no algorithm has been proposed to conduct differential analysis of directed networks characterized by SEM. While a naive approach would separately construct each individual network and identify common and differential structures, this approach fails to take advantage of the commonality as well as sparse differential structures of the paired networks, leading to higher false positive rate or lower power. In this light, we introduce a novel statistical method, specially in the directed network regime, to conduct differential analysis of two networks via appropriate reparameterization of the corresponding models. There are two major features of our method. Firstly, we jointly model the commonality and difference between two networks explicitly. This helps us to gain dramatic performance improvements over the naive construction method. Secondly, benefiting from the flexible framework of SEMs, we are able to conduct differential analysis of directed networks. Most importantly, our method allow for both acyclic and cyclic networks. Compared to the other methods, directionality and allowing for cyclicity are crucial for many network studies, especially in constructing gene regulatory networks. As far as we know, our method is the first work on differential analysis of directed networks that enjoys the two promising features.

The rest of this paper is organized as follows. We first introduce the model and its identifiability condition in Section 2.1 and Section 2.2, respectively. Then, we present our proposed method of Reparameterization-based Differential analysis of directed Networks, termed as ReDNet, in Section <ref>. The theoretical justification of the proposed method is described in Section <ref>. Section <ref> includes our studies on synthetic data showing the superior performance of our method, as well as an analysis of the Genotype-Tissue Expression (GTEx) data sets. We conclude our paper with brief discussion in Section <ref>.





§  METHODS


Here we first introduce the model and its identification condition, and then describe our proposed 𝐑𝐞𝐃𝐍𝐞𝐭 method for identifying common and differential structures between two directed networks, followed with its theoretical justification.



 §.§ The Model


We consider two networks, each describing the dependencies among a common set of variables or nodes in a unique population. For each  node i∈{1, 2, …,p} in network k∈{1, 2}, its regulation structure can be represented by the following equation,

    𝐘^(k)_i_node  i  = 𝐘^(k)_-iγ^(k)_i_regulation by others + 𝐗^(k)ϕ^(k)_i_anchoring regulation+ ϵ^(k)_i_error,

where 𝐘^(k)_i is the i-th column of 𝐘^(k) and 𝐘^(k)_-i is the submatrix of 𝐘^(k) by excluding 𝐘^(k)_i, with 𝐘^(k) a n^(k)× p matrix. 𝐗^(k) is a n^(k)× q matrix with each column standardized to have ℓ_2 norm √(n^(k)). The vectors γ^(k)_i and ϕ^(k)_i encode the inter-nodes and anchoring regulatory effects, respectively. The index set of non-zeros of ϕ^(k)_i is known and denoted by 𝒜^(k)_i, in other words, 𝒜^(k)_i = supp(ϕ^(k)_i). The support set 𝒜^(k)_i indexes the direct causal effects for the i-th node, and can be prespecified based on the domain knowledge. However, the size of nonzero effect ϕ^(k)_i is unknown and can be estimated. Further property of 𝒜^(k)_i will be discussed in Section <ref>. All elements of the error term are independently distributed following a normal distribution with mean zero and standard deviation σ_i^(k). We assume that the matrix 𝐗^(k) and the error term ϵ^(k)_i are independent of each other. However 𝐘^(k)_-i and ϵ^(k)_i may correlate with each other. 𝐘^(k) and 𝐗^(k) include observed endogenous variables and exogenous variables, respectively.

By combining the p linear equations in (<ref>), we can rewrite the two sets of linear equations in a systematic fashion as two structural equation models below,

    𝐘^(1) = 𝐘^(1)Γ^(1) + 𝐗^(1)Φ^(1)+ ℰ^(1), 
    𝐘^(2) = 𝐘^(2)Γ^(2) + 𝐗^(2)Φ^(2)+ ℰ^(2),

where each matrix Γ^(k) is p × p with zero diagonal elements and represents the inter-nodes regulatory effects in the corresponding network. Specifically, excluding the i-th element (which is zero) from the i-th column of Γ^(k) leads to  γ^(k)_i. The q × p matrix Φ^(k) contains the anchoring regulatory effects and its i-th column is ϕ_i^(k). Each error term ℰ^(k) is n^(k)× p and has the error term ϵ^(k)_i as its i-th column.

Figure <ref> gives an illustrative example of networks with three nodes and one anchoring regulation per node for the structural equations in (<ref>). For example, with anchoring regulation on node Y_1, X_1 has a direct effect on node Y_1 but indirect effects on node Y_2 and Y_3 via Y_1.



For each network k, its full model in (<ref>) can be further transformed into the reduced form as follows,

    𝐘^(k) = 𝐗^(k)π^(k) +ξ^(k),

where the q × p matrix π^(k) = Φ^(k)(𝐈-Γ^(k))^-1 and the transformed error term ξ^(k) = ℰ^(k)(𝐈-Γ^(k))^-1. The reduced model (<ref>) reveals variables observed in ^(k) as instrumental variables which will be used later to correct for the endogeneity issue. Otherwise, directly applying any regularization based regression to equation (<ref>) will result in  non-consistent or suboptimal estimation of model parameters <cit.>.





 §.§ The Model Identifiability 
 

Here we introduce an identifiability assumption which helps to infer an identifiable system (<ref>) from available data. We assume that each endogenous variable is directly regulated by a unique set of exogenous variables as long as it regulates other endogenous variables. That is, any regulatory node needs at least one anchoring exogenous variable to distinguish the corresponding regulatory effects from association. Explicitly let ℳ^(k)_i0 denote the index set of endogenous variables which either directly or indirectly regulate the i-th endogenous variable in the k-th network. Thus, 𝒜_i^(k)⊆ℳ^(k)_i0. The model identification condition can be stated in the below.


Assumption 1. For any i=1,⋯, p, 𝒜_i^(k)∅ if there exists j such that i∈ℳ^(k)_j0. Furthermore, 𝒜_i^(k)∩𝒜_j^(k)=∅ as long as i j.


This assumption is slightly less restrictive than the one employed by <cit.>, and is a sufficient condition for model identifiability as it satisfies the rank condition in <cit.>. It can be further relaxed to allow nonempty 𝒜_i^(k)∩𝒜_j^(k) as long as each regulatory node has its own unique anchoring exogenous variables.

The above identifiability assumption not only identifies γ^(k)_i in model (<ref>) from π^(k) in model (<ref>) but also helps reveal regulatory directionality of the networks. As illustrated in Figure <ref>, we can not recover the directionality between nodes Y_1 and Y_2 without the extra information provided by the direct causal factors X_1 and X_2 because all four sub-networks consisting of Y_1 and Y_2 (without X_1 and X_2) will be Markov equivalent. The known set 𝒜_j^(k) serves as external prior knowledge which helps recover the directionality. In our two-stage construction of the differential network, the additional anchors X_1 and X_2 serve as instrumental variables in the calibration stage, since both X_1 and X_2 are independent of the error terms. The present direct causal effects from ^(k) together with Assumption 1 differentiates our approach from the classical graphical models <cit.> or the PC algorithm approaches <cit.>, since those methods either cannot recover edge directions or do not allow for cyclic structures due to lack of additional direct causal effects from ^(k).







 §.§ Two-Stage Differential Analysis of Networks 


Here we intend to develop a regularized version of the two-stage least squares. We first screen for exogenous variables and conduct ℓ_2 regularized regression of each endogenous variable against screened exogenous variables to obtain its good prediction which helps address the endogeneity issue in the following stage. At the second stage, we reparametrize the model to explicitly model the common and differential regulatory effects and identify them via the adaptive lasso method.





  §.§.§ The Calibration Stage


To address the endogeneity issue, we aim for good prediction of each endogenous variable following the reduced model in (<ref>). However, in the high-dimensional setting, the dimension q of 𝐗^(k) can be much larger than the sample size n^(k), and any direct prediction with all exogenous variables may not produce consistent prediction. Note that both <cit.> and <cit.> proposed to conduct variable selection with lasso or its variants and predict with selected exogenous variables. We here instead propose to first screen for exogenous variables with ISIS <cit.>, and then apply ridge regression to predict the endogenous variables with screened exogenous variables. While variable screening is more robust and provides higher coverage of true variables than variable selection, its combination with ridge regression puts less computational burden. Furthermore, as shown by <cit.>, ridge regression performs well in predicting the endogenous variables.

Let ℳ^(k)_i denotes the selected index set for i-th node in k-th network from the variable screening which reduces the dimension from q to d=|ℳ^(k)_i|. The Sure Independence Screening Property in <cit.> can be directly applied in our case to guarantee that ℳ^(k)_i covers the true set ℳ^(k)_i0 with a large probability.


Assumption 2. n^(1) and n^(2) are at the same order, i.e., n_min = min(n^(1), n^(2)) ≍ n^(1)≍ n^(2), and p ≍ q.


  Assuming Conditions 1-4 in the supplemental materials which restrict positive τ̃ and κ̃, under Assumption 2, there exists some θ∈ (0, 1-2κ̃-τ̃) such that, when d =|ℳ_i^(k)|= O((n_min)^1-θ), we have, for some constant C>0,

    ℙ( ℳ^(k)_i0⊆ℳ^(k)_i ) = 1 - 𝒪(exp{-C (n^(k))^1-2κ̃/log(n^(k))}).



Hereafter we assume that ℳ^(k)_i successfully covers the true set ℳ^(k)_i0 for convenience of stating the following assumptions and theorems. That is, the probability of successful screening is not incorporated into our assumptions or theorems in the below.

For node i in network k, let 𝐗^(k)_ℳ_i^(k) denotes the submatrix of 𝐗^(k) with prescreened columns which are indexed by ℳ_i^(k). With π_i^(k) denoting the i-th column of π^(k), the subvector of π_i^(k) indexed by ℳ_i^(k) will be simply denoted by π^(k)_ℳ_i^(k) without confusion. Such simplified notations will apply to other vectors and matrices in the rest of this paper.

With d pre-screened exogenous variables, we can apply ridge regression to the model

    𝐘^(k)_i = 𝐗_ℳ_i^(k) ^(k) π^(k)_ℳ_i^(k) +ξ^(k)_i,

to obtain the estimates π̂^(k)_ℳ_i^(k) of π^(k)_ℳ_i^(k), and predict 𝐘^(k)_i with 𝐘̂^(k)_i = 𝐗_ℳ_i^(k)^(k)π̂^(k)_ℳ_i^(k).





  §.§.§ The Construction Stage


With known 𝒜_i^(k), we can rewrite model (<ref>) as,

    𝐘^(k)_i= 𝐘^(k)_-iγ^(k)_i +𝐗^(k)_𝒜_i^(k)ϕ^(k)_𝒜_i^(k)+ϵ^(k)_i.

Before we use the predicted 𝐘^(k) to identify both common and differential regulatory effects across the two networks, we first reparametrize the model so as to define differential regulatory effects explicitly,

    β_i^- = γ_i^(1)-γ_i^(2)/2, β_i^+ =γ_i^(1)+γ_i^(2)/2.

Here β_i^- and β_i^+ represent the differential and average regulatory effects between the two networks, respectively. We need compare β_i^+ with β_i^- to identify the common regulatory effects, that is, effects of all regulations with nonzero values in β_i^+ but zero values in β_i^-.

Note that other differential analysis of networks may suggest a different reparametrization to identify common and differential regulatory effects. For example, in a typical case-control study, we may expect few structures in the case network mutated from the control network. While we are interested in identifying differential structures in the case network, we may be also interested in identifying baseline network structures in the control network. Therefore we may reparametrize the model with the regulatory effects in the control network, as well as the differential regulatory effects defined as the difference of regulatory effects between case and control networks. We want to point out that the method described here still applies and we can also derive similar theoretical results as follows.

Following the reparametrization in (<ref>), we can rewrite model (<ref>) as follows,

    [ 𝐘_i^(1); 𝐘_i^(2) ] =
    [  𝐘_-i^(1)  𝐘_-i^(1);  𝐘_-i^(2) -𝐘_-i^(2) ][ β_i^+; β_i^- ]+
          [ 𝐗_𝒜_i^(1)^(1)             0;             0 𝐗_𝒜_i^(2)^(2) ][ ϕ_𝒜_i^(1)^(1); ϕ_𝒜_i^(2)^(2) ]+
    [ ϵ_i^(1); ϵ_i^(2) ].


Denote

    𝐘_i =
    [ 𝐘_i^(1); 𝐘_i^(2) ],    𝐙_-i =
    [  𝐘_-i^(1)  𝐘_-i^(1);  𝐘_-i^(2) -𝐘_-i^(2) ], 
       β_i =
    [ β_i^+; β_i^- ],
       ϵ_i =
    [ ϵ_i^(1); ϵ_i^(2) ].

Further define the projection matrix for each network,

    𝐇_i^(k) = I_n^(k) - 𝐗_𝒜_i^(k)^(k)(𝐗_𝒜_i^(k)^(k)T𝐗_𝒜_i^(k)^(k))^-1𝐗_𝒜_i^(k)^(k)T.

Applying the projection matrix 𝐇_i = {𝐇_i^(1), 𝐇_i^(2)} to both sides of model (<ref>), we can remove the exogenous variables from the model and obtain,

    𝐇_i 𝐘_i = 𝐇_i 𝐙_-iβ_i + 𝐇_i ϵ_i.


To address the endogeneity issue, we predict 𝐙_-i by replacing its component 𝐘^(k)_-i with the predicted value Ŷ^(k)_-i from the previous stage, and then regressing 𝐇_i 𝐘_i against 𝐇_i 𝐙̂_-i with the adaptive lasso to consistently estimate β_i. That is, an optimal β_i can be obtained as,

    β̂_i  = β_iarg min {1/n ||𝐇_i 𝐘_i - 𝐇_i 𝐙̂_-iβ_i ||^2_2 +  λ_iω_i^T|β_i|_1 },

where |β_i|_1 is a vector taking elementwise absolute values of β_i, ω_i is the adaptive weights whose components are inversely proportional to the components of an initial estimator of β_i, and λ_i is the adaptive tuning parameter.

The two-stages algorithm is summarized in Algorithm <ref>. With the estimator β̂_i from the second stage, we can accordingly obtain estimators γ̂_i^(1) = β̂_i^+ + β̂_i^- and γ̂_i^(2) = β̂_i^+ - β̂_i^-.







 §.§ Theoretical Analysis
 

As shown in Theorem <ref>, a screening method like ISIS <cit.> can identify ℳ^(k)_i with size d = O(n_min^1-θ) which covers the true set ℳ^(k)_i0 with a sufficiently large probability. For the sake of simplicity and without loss of generality, in the following we assume ℳ^(k)_i0⊆ℳ^(k)_i.

We first investigate the consistency of predictions from the first stage. The consistency properties will be characterized by prespecified sequences f^(k)=o(n^(k)) but f^(k)→∞ as n^(k)→∞. We also denote f_max = f^(1) f^(2), i.e., max{f^(1),f^(2)}.

The following assumption is required for the consistency properties.


Assumption 3. For each network k, the singular values of 𝐈 - Γ^(k) are positively bounded from below, and there exist some positive constants c_1^(k) and c_2^(k) such that, for each node i, max_δ=1 (n^(k))^-1/2^(k)_ℳ_i^(k)δ≤ c_1^(k) and min_δ=1 (n^(k))^-1/2^(k)_ℳ_i^(k)δ≥ c_2^(k). Furthermore, the ridge parameter λ^(k)_i=o(n_min).


For the ease of exposition, we will omit the subscript ℳ_i^(k) from ^(k)_ℳ_i^(k) henceforth, and accordingly use π_i^(k) and π̂_i^(k) which include the zero components of excluded predictors.

Denote 𝐗 = {𝐗^(1),𝐗^(2)}, and

    𝐙 =
    [  𝐘^(1)  𝐘^(1);  𝐘^(2) -𝐘^(2) ],    Π =
    [  π^(1)  π^(1);  π^(2) -π^(2) ].

We use Π_j to denote the j-th column of the matrix Π and π^(k)_j to denote the j-th column of the matrix π^(k). We also use 𝐙̂ and Π̂ to denote the prediction of 𝐙 and estimate of Π, respectively. Note that, with the ridge parameter λ^(k)_i for the ridge regression taken on node i in network k, we have r_i^(k) = (λ^(k)_i)^2π^(k)_i^2/n^(k) and hence define r_max = 1≤ i ≤ pmax [r_i^(1) r_i^(2)]. Then the estimation and prediction losses at the first stage can be summarized in the following theorem.

 
Under Assumptions 1-3, for each j∈{1,2,…,2p}, there will exist some constant C_1 and C_2 such that, with probability at least 1- e^-f^(1) - e^-f^(2), 

1. Π̂_j - Π_j^2 ≤ C_1 (d r_max f_max)/n_min;

2. 𝐗(Π̂_j - Π_j)^2 ≤ C_2 (d r_max f_max).


The proof is detailed in the supplemental materials.

Note that these two sets of losses can be controlled by the same upper bounds across the two networks with probability at least 1 - e^-f^(1)+log (p) - e^-f^(2) + log (p). Therefore, f^(k) can be selected such that f^(k) - log(p)→∞, which will provide a probability approaching one to have the network-wide losses approaching zero.

Furthermore, the dimension p can be divergent up to an exponential order, say p=e^n_min^c for some c∈(0,1). We can set f^(1) = f^(2) = n_min^(1+c)/2 and, apparently, f^(k)=o(n_min) but f^(k) - log(p) = n_min^(1+c)/2-n_min^c →∞.

Since the ridge parameter λ^(k)_i=o(n_min), r_i^(k)=π^(k)_i^2 × o(n_min). Therefore, when all π^(k)_i are uniformly bounded, we have r_max = o(n_min). Otherwise, the ridge parameter λ^(k)_i should be adjusted accordingly to control both estimation and prediction losses.

Before we characterize the consistency of estimated regulatory effects on the second stage, we first introduce the following concept of restricted eigenvalue which is used to present an assumption.

The restricted eigenvalue of a matrix 𝐀 on an index set 𝒮 is defined as

    𝐀𝒮 = δ_𝒮^c≤ 3 δ_𝒮min𝐀δ/√(n)δ_𝒮.



For the i-th node, we use 𝒮_i to denote the non-zero indices of β_i, i.e., 𝒮_i = supp(β_i).
Further denote

    Π_-i =
    [  π_-i^(1)  π_-i^(1);  π_-i^(2) -π_-i^(2) ].


As in <cit.>, we impose the following restricted eigenvalue condition on the design matrix in (<ref>).


Assumption 4. There exists a constant ϕ_0 > 0 such that 𝐇_i𝐗Π_-i𝒮_i≥ϕ_0. Furthermore, ω_𝒮_i_∞≤ω_𝒮_i^c_-∞.


Let n=n^(1)+n^(2), c_max=c^(1)_1 c^(2)_1, and 𝐁 = [β_1,β_2,…, β_p]. The matrix norms · and ·_∞ are the maximum of column and row sums of absolute values of the matrix, respectively. For a vector, we define ·_∞ and ·_-∞ to be the maximum and minimum absolute values of its components. Then, we can derive the following loss bounds for the estimation and prediction at the second stage on the basis of Theorem <ref>.

 
Suppose that, for node i, the adaptive lasso at the second stage takes the tuning parameter λ_i≍ω_i_-∞^-1𝐁Π√( (d  r_max f_max)log(p) / n_min), and √((d  r_max f_max) / n) + c_maxΠ≤√(c_max^2Π^2+ϕ_0^2/(64C_2 |𝒮_i|)). Let h_n = (𝐁^2  1) ×( (nΠ^2/d)  (d r_max f_max))log(p). Under Assumptions 1-4, there exist positive constants C_3 and C_4 such that, with probability at least 1- 3e^-C_3 h_n + log(4pq) - e^-f^(1)+log(p) - e^-f^(2) + log(p),

1. Estimation Loss:

    β̂_i-β_i ≤ 8C_4 |𝒮_i| ×
          ω_𝒮_i_∞𝐁Π/ϕ_0^2 ω_i_-∞√((d  r_max f_max )log(p)/n_min);

2. Prediction Loss:

    1/n𝐇_i𝐙̂_-i (β̂_i - β_i  )^2 ≤ C_4^2 |𝒮_i|×
          ω_𝒮_i_∞^2 𝐁^2Π^2/ϕ_0^2 ω_i_-∞^2(d  r_max f_max)log(p)/n_min.




The main idea of the proof is to take advantage of the commonly used restricted eigenvalue condition and irrepresentable condition for lasso-type estimator. However, the design matrix in our case includes predicted values instead of the original one, which complicates the proof. We claim that the restricted eigenvalue and irrepresentable condition still hold for the predicted design matrix as long as the estimation and prediction losses are well controlled at the calibration stage. The proof is detailed in the supplemental materials.

The available anchoring regulators as required by Assumption 1 implies that both 𝐁>0 and Π>0, so h_n/log(p)→∞. That is, these loss bounds hold with a sufficient large probability with properly chosen f^(k).

The two sets of losses in Theorem <ref> can also be controlled across the whole system by the same upper bounds defined by replacing |𝒮_i| with s_max = max_i|𝒮_i|, with probability at least 1- 3e^-C_3 h_n + log(4 q) + 2log(p) - e^-f^(1)+2log(p) - e^-f^(2) + 2log(p). When both p and q are divergent up to an exponential order, say p≍ q≍ e^n_min^c for some c∈(0,1), we can set f^(1) = f^(2)=n_min^(1+c)/2 to guarantee the bounds at a sufficient large probability. However, the bounds are determined by (d  r_max f_max)log(p) which is o(n_min) only when c<min(1/3,θ). Therefore, if s_max also diverges up to n_min^c̃ with c̃<min(1/4,θ/2,1-θ), the losses can be well controlled for c<min((1-4c̃)/3, θ-2c̃).

Note that, with properly chosen f^(1) and f^(2), these losses are well controlled at o (n_min), revealing the fact that we need to have sufficient observations for each network for consistent differential analysis of the two networks.

Let W_i = diag{_i}. Denote ℐ_i = 1/nΠ_-i^T𝐗^T𝐇_i𝐗Π_-i and ℐ̂_i = 1/nΠ̂_-i^T𝐗^T𝐇_i𝐗Π̂_-i. Let  be a submatrix of ℐ_i with rows and columns both indexed by 𝒮_i, and  be a submatrix of ℐ_i with rows and columns indexed by 𝒮_i^c and 𝒮_i, respectively.  and  are similarly defined from ℐ̂_i. We further define the minimal signal strength b_i =j∈𝒮_imax |β_ij| and ψ_i = [-1]  W_𝒮_i.

The following assumption, reminiscent of the adaptive irrepresentable condition in <cit.>, helps investigate the selection consistency of regulatory effects.


Assumption 5. (Weighted Irrepresentable Condition) There exists a constant τ∈ (0,1) such that W_𝒮_i^c^-1[-1] W_𝒮_i < 1- τ.


 (Variable Selection Consistency) Denote Ŝ_i= { j: β̂_ij 0 }. Suppose that, for each node i,  is invertible, b_i > λ_iψ_i/(2-τ), and √((d  r_max f_max) / n) + c_maxΠ≤√(c_max^2Π^2+min(ϕ_0^2/64, τ (4-τ)^-1ω_i_-∞/ψ_i)/(C_2 |𝒮_i|)). Under Assumptions 1-5, there exists some constant C_5>0 such that 𝒮̂_i = 𝒮_i with probability at least 1- 3e^-C_5 h_n + log(4pq) - e^-f^(1)+log(p) - e^-f^(2) + log(p).


This theorem implies that our proposed method can identify both common and differential regulatory effects between the two networks with a sufficiently large probability. On the other hand, the assumed weighted irrepresentable condition means that the true signal should not correlate too much with irrelevant predictors so as to conduct a successful differential analysis. The corresponding proof is detailed in the supplemental materials.





§ EXPERIMENTS




 §.§ Synthetic Data Evaluation




Here we report on experiments with synthetic data to show the superior performance of our method. We compare the method  to a naive differential analysis which employs the 2SPLS method proposed by <cit.> to construct each network separately. Note that the 2SPLS method is modified here by applying ISIS to screen exogenous variables before conducting ridge regression to predict endogenous variables, making the naive differential analysis comparable to .

Synthetic data are generated from both acyclic and cyclic networks involving 1000 endogenous variables, with the sample size from 200 to 300. Each network includes a subnetwork of 50 endogenous variables, whose shared and differential structures will be investigated against its pair. On average, each endogenous variable has one regulatory effect in a sparse subnetwork, and three regulatory effects in a dense network. While each pair of subnetworks in comparison share many identical regulatory effects, they also share five regulatory effects but with opposite signs, and each network has five unique regulatory effects (so the total number of differential regulatory effects is 15). The nonzero regulatory effects were independently sampled from a uniform distribution over the range [-0.8,-0.3]∪[0.3,8]. While assuming each node is directly regulated by one exogenous variable, each exogenous variable was sampled from discrete values 0,1 and 2 with probabilities 0.25, 0.5 and 0.25, respectively. All of the noise terms were independently sampled from the normal distribution N(0,0.1^2). We also conducted differential analysis between two networks with both ^(1)^(2) and ^(1) = ^(2) as in practice the paired networks may or may not share identically valued exogenous variables.

We evaluate the the performance in terms of the false discovery rate (FDR), power and Matthews correlation coefficient (MCC) <cit.>. Let
TP, TN, FP and FN denote the numbers of true positives, true negatives, false positives, and false negatives, respectively.
MCC is defined as,

    MCC = TP×TN - FP×FN/√(  (TP + FP)(TP + FN)(TN + FP)(TN + FN) ).

-10pt
Here we refer nonzero effects as positives and zero effects as negatives. The MCC varies from 0 to 1 with larger values implying better variable selection.

In each differential analysis, the ridge regression employed the generalized cross validation <cit.> to select the ridge parameter, and the adaptive lasso used 10-fold cross-validation to choose its tuning parameter. Following the recommendation by <cit.>, (n^(k))^0.9 variables are screened by ISIS.

For each type of networks, 100 synthetic data sets were generated, and the differential analysis results are summarized in Figure <ref>. Overall, both  and the naive approach maintain high power in identifying differential regulatory effects. However, the naive approach fails to identify common regulatory effects and tends to report FDR over 80% on differential regulatory effects. Such a tendency to report false positives by the naive approach results in lower MCC, with dramatic decrease in identifying differential regulatory effects.

While both methods performed stably across networks with ^(1)^(2) and ^(1) = ^(2),  performed better in identifying differential regulatory effects from dense networks than sparse networks in terms of FDR and MCC. However, the naive approach tends to report even higher FDR and so much lower MCC when identifying differential regulatory effects from dense networks. Nonetheless, the naive approach fails to identify common regulatory effects for each type of networks so the corresponding FDR and MCC are undefined.

We also calculated the standard errors (SE) of the reported FDR, power, and MCC over 100 synthetic data sets (the results are not shown). They are all small with most at the scale of thousandth and others at the scale of hundredth. Therefore,  performed robustly in differential analysis of networks, and the 2SPLS approach by <cit.> performed also robustly in constructing single networks.





 §.§ The Genotype-Tissue Expression Data




We performed differential analysis of gene regulatory networks on two sets of genetic genomics data from the Genotype-Tissue Expression (GTEx) project <cit.>, with one collected from human whole blood (WB) and another one from human muscle skeletal (MS). The WB and MS data included genome-wide genetic and genotypic values from 350 and 367 healthy subjects, respectively. Both data sets were preprocessed following <cit.> and <cit.>, resulting in a total of 15,899 genes and 1,083,917 single nucleotide polymorphisms (SNPs) being shared by WB and MS.

Expression quantitative trait loci (eQTL) mapping <cit.> was conducted and identified 9875 genes with at least one marginally significant cis-eQTL (with p-value< 0.05). For each gene, we further filtered its set of cis-eQTL by controlling the pairwise correlation under 0.9 and keeping up to three cis-eQTL which have the strongest association with the corresponding gene expression. These cis-eQTL serve as anchoring exogenous variables for the genes, and expression levels of different genes are endogenous variables. At completion of preprocessing data, we have 9,875 endogenous variables and 23,920 exogenous variables.

We applied  to infer the differential gene regulation on a set of eighty target genes, which had largest changes on gene-gene correlation between the two tissues. We identified a total of 640 common and 572 differential regulations on the eighty target genes. To evaluate the significance of identified regulations, we bootstrapped 100 data sets, and conducted differential analysis on each bootstrap data set. As summarized in Table <ref>, 50, 43 and 34 differential regulatory effects were identified in over 70%, 80% and 90% of the bootstrap data sets, respectively.



The top five subnetworks bearing differential regulations on some of the eighty target genes were shown in Figure <ref>. We also constructed the differential networks using the naive approach (the results are not shown), and reported more regulations which cover the reported ones by . This concurs with our observation in the synthetic data evaluation that the naive approach tends to report higher false positives, especially for differential regulatory effects.





§ CONCLUSION
 

We have developed a novel two-stage differential analysis method named . The first stage, i.e., the calibration stage, aims for good prediction of the endogenous variables, and the second stage, i.e., the construction stage, identifies both common and differential network structures in a node-wise fashion. The key idea of  method is to appropriately reparametrize the independent models into a joint model so as to estimate differential and common effects directly. This approach can dramatically reduce the false discovery rate. In the experiments with synthetic data, we demonstrated the effectiveness of our method, which outperformed the naive approach with a large margin. Note that  allows independently conducting all ℓ_2 regularized regressions at the same time at the first stage, and all ℓ_1 regularized regressions at the same time at the second stage. Therefore,  not only permits parallel computation but also allows for fast subnetwork construction to avoid potential huge computational demands from differential analysis of large networks.

There are some interesting directions for future research. Firstly, it is worthwhile to explore other re-parametrization approaches such as baseline reparametrizaiton in a case-control study. Secondly, while we only consider differential analysis of two networks, it is possible to generalize our method to compare multiple networks, demanding more complex reparametrization. Finally, applying the proposed method for fully differential analysis of 53 tissues in the GTEx project still provides challenging computational and methodological issues.



§ ACKNOWLEDGMENTS


The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. The data used for the analysis described in this paper were obtained from dbGaP accession number phs000424.v7.p2 on 08/18/2017.

plainnat







Supplementary Materials

	12pt
	Differential Analysis of Directed Networks




There are five parts. Firstly, we collect in Section <ref> all notations used in our paper and here. We then describe the four conditions which help define the positive pair τ̃ and κ̃ for Theorem 1, and further prove Theorem 1 in Section <ref>. In Section <ref>, we prove Theorem 2 which provides bounds for both estimation and prediction losses at the calibration stage. In Section <ref>, we prove Theorem 3 which provides bounds for both estimation and prediction losses at the construction stage. In Section <ref>, we prove the variable selection consistency in Theorem 4.




§ NOTATIONS
 

Unless otherwise claimed, we will follow the notations defined here throughout the paper and supplementary materials.

For a vector, ||·||_2 and ||·||_1 denote the ℓ_2 and ℓ_1 norms, respectively; ·_∞ and ·_-∞ are defined to be the maximum and minimum absolute values of its components, respectively; |·|_1 implies taking element-wise absolute values of the vector so is itself a vector. For a matrix A=(a_ij)_m× n, A_1=max _1≤ j≤ n∑ _i=1^m|a_ij|, i.e., the maximum column sum of absolute values of its components, and A_∞=max _1≤ i≤ m∑ _j=1^n|a_ij|, i.e., the maximum row sum of absolute values of its components.

For a vector a and index set 𝒮, a_i , a_-i, and a_𝒮 denote the i-th entry, the subvector excluding the i-th entry in a, and the subvector of a indexed by 𝒮, respectively. For a matrix A, A_i and A_-i denote its i-th column and the submatrix of A excluding its i-th column, respectively. For a vector a_i and an index set 𝒮_i both sharing the same subscript, the subvector of a_i indexed by 𝒮_i is denoted by a_𝒮_i for simplicity. Similarly, the submatrix of a matrix A_i including columns indexed by the set 𝒮_i is denoted by A_𝒮_i for simplicity.

a ∨ b and a∧ b denote the maximum and minimum of a and b, respectively. λ_min(·) and λ_max(·) denote the minimum and maximum eigenvalues of the corresponding matrix, respectively. · denotes the expectation, and ℙ(·) denotes the probability of an event. Symbol ≍ denotes two terms at the same order. (·) denotes the trace of the corresponding matrix. For a set S, |S| denotes the number of its elements. For positive integers j and p, j|p denotes the remainder of j when divided by p.

Throughout the paper and here, C_1, C_2, …, c_1, c_2, …, c̃_1, c̃_2, ⋯, t_1, t_2, … are some positive constant numbers.





§ THE CONDITIONS AND PROOF OF THEOREM 1
 

For each k ∈{1,2}, the reduced model (3) includes p regression models, i.e., for i=1, 2, ⋯, p,

    𝐘_i^(k) = 𝐗^(k)π_i^(k) + ξ_i^(k).

Here we first state the four conditions in  Fan and Lv [2008] which restrict the positive pairs τ^(k) and κ^(k) so as to define τ̃=max{τ^(1),τ^(2)} and κ̃=max{κ^(1),κ^(2)} for Theorem 1, and then prove that we can successfully screen variables for each of the above linear regression model.

Denote Y_ji^(k), X_jl^(k), ξ_ji^(k), and π_ji^(k) as the j-th row of 𝐘_i^(k), 𝐗_l^(k), ξ_i^(k), and π_i^(k), respectively. Further denote Σ^(k) the variance-covariance matrix of the q random variables in observing 𝐗^(k). For any ℳ⊂{1, 2,⋯, q}, denote Σ^(k)_ℳ the variance-covariance matrix of the random variables in observing 𝐗^(k)_ℳ.

6ptCondition 1. Each ξ_ji^(k) is normally distributed with mean zero. (Σ^(k))^-1/2𝐗^(k)T is observed from a spherically symmetric distribution, and has the concentration property: there exist some constants c̃_1^(k),c̃_2^(k) > 1 and c̃_3^(k) > 0 such that, for any ℳ⊂{1, 2,⋯, q} with |ℳ|≥c̃_1^(k) n^(k), the eigenvalues of |ℳ|^-1𝐗^(k)_ℳ  (Σ^(k)_ℳ)^-1/2 (Σ^(k)T_ℳ)^-1/2𝐗^(k)T_ℳ are bounded either from above by c̃_2^(k) or from below by 1/c̃_2^(k) with probability at least 1-exp(-c̃_3^(k) n^(k)).

6ptCondition 2. var(Y_ji^(k)) = O(1). For some κ^(k)≥ 0, c̃_4^(k)>0, and c̃_5^(k) > 0,

    l ∈ℳ_i0^(k)min|π_li^(k)| ≥c̃_4^(k)/(n^(k))^κ^(k)  and  l ∈ℳ_i0^(k)min|cov( (π_li^(k))^-1 Y_ji^(k), X_jl^(k))| ≥c̃_5^(k).


6ptCondition 3. log(q) = O((n^(k))^c̃) for some c̃∈ (0, 1-2κ^(k)).

6ptCondition 4. There are some τ^(k)≥ 0 and c̃_6^(k) > 0 such that λ_max(Σ^(k) ) ≤c̃_6^(k) (n^(k))^τ^(k).
6pt


	Following the Sure Independence Screening Property by Fan and Lv [2008], there exists some θ^(k)∈ (0, 1-2κ^(k)-τ^(k)) such that, when d^(k) =|ℳ_i^(k)|= O((n^(k))^1-θ^(k)), we have, for some constant C>0,
	
    ℙ( ℳ^(k)_i0⊆ℳ^(k)_i ) = 1 - 𝒪(exp{-C (n^(k))^1-2κ^(k)/log(n^(k))}).

	Let θ = min(θ^(1),θ^(2)), then for d^(k)=|ℳ_i^(k)|≡  d=O(n_min^1-θ), we have
	
    ℙ( ℳ^(k)_i0⊆ℳ^(k)_i ) = 1 - 𝒪(exp{-C (n^(k))^1-2κ̃/log(n^(k))}).







§ PROOF OF THEOREM 2
 

Note that ξ^(k) = ℰ^(k) (I-Γ^(k))^-1 for k∈{1, 2}. Suppose that the singular values of both (I-Γ^(k)) are positively bounded from below by a constant c. Denote σ_i^(k)2=(ϵ_ji^(k)) and σ̃_i^(k)2=(ξ_ji^(k)). Then σ̃_i^(k)≤σ_pmax/c  = 1≤ i ≤ pmax  (σ_i^(1)σ_i^(2))/c.


	Under Assumptions 1-3, for each network k ∈{1,2} in the calibration step, there exist positive constants C_1^(k) and C_2^(k) such that, with probability at least 1-e^-f^(k),
	
		
  * (Estimation Loss) π̂^(k)_i - π^(k)_i^2 ≤ C_1^(k)(r_i^(k) d  f^(k))/n^(k);
		
  * (Prediction Loss) ^(k)(π̂^(k)_i  - π^(k)_i )^2 / n^(k)≤ C_2^(k)(r_i^(k) d  f^(k))/n^(k).
	



	We have the closed form ridge estimator π̂_ℳ_i^(k)^(k) for the linear model 𝐘_i^(k) = 𝐗_ℳ_i^(k)^(k)π_ℳ_i^(k)^(k) + ξ_i^(k).
	
    π̂_ℳ_i^(k)^(k)= (𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k) +λ^(k)_iI_d)^-1𝐗_ℳ_i^(k)^(k)T𝐘_i^(k),

	where λ^(k)_i is the ridge tuning parameter. Plugging in the equation 𝐘^(k)_i = 𝐗_ℳ_i^(k)^(k)π_ℳ_i^(k)^(k) +ξ_i^(k), we have
	
    π̂_ℳ_i^(k)^(k)=    {(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k) +λ^(k)_iI_d)^-1𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)π_ℳ_i^(k)^(k)}
       + {(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k) +λ^(k)_iI_d)^-1𝐗_ℳ_i^(k)^(k)Tξ_i^(k)}.

	The difference between the ridge estimator π̂_ℳ_i^(k)^(k) and the true π_ℳ_i^(k)^(k) can be written as
	
    π̂_ℳ_i^(k)^(k)- π_ℳ_i^(k)^(k)   = -λ^(k)_i(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k) +λ^(k)_iI_d)^-1π_ℳ_i^(k)^(k) + (𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k) +λ^(k)_iI_d)^-1𝐗_ℳ_i^(k)^(k)Tξ_i^(k).

	
	For simplicity, we denote the composite forms of  π_ℳ_i^(k)^(k) and 𝐗_ℳ_i^(k)^(k) as follows,
	
    π̃_ℳ_i^(k)^(k) = -λ^(k)_i(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k) +λ^(k)_iI_d)^-1π_ℳ_i^(k)^(k); 
          𝐗̃_ℳ_i^(k)^(k) = 𝐗_ℳ_i^(k)^(k)(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)+λ^(k)_iI_d)^-1.

	Then we have the following simplified form of the difference,
	
    π̂_ℳ_i^(k)^(k)- π_ℳ_i^(k)^(k) = π̃_ℳ_i^(k)^(k) + 𝐗̃_ℳ_i^(k)^(k)Tξ_i^(k).

	
	To obtain the ℓ_2 norm losses of estimation and prediction, we write
	
    π̂_ℳ_i^(k)^(k)- π_ℳ_i^(k)^(k)^2
       =   π̃_ℳ_i^(k)^(k)Tπ̃_ℳ_i^(k)^(k)_T_21  + 2π̃_ℳ_i^(k)^(k)T𝐗̃_ℳ_i^(k)^(k)Tξ_i^(k)_T_22 + ξ_i^(k)T𝐗̃_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)Tξ_i^(k)_T_23,
    𝐗_ℳ_i^(k)^(k)(π̂_ℳ_i^(k)^(k)- π_ℳ_i^(k)^(k))^2
       =   π_ℳ_i^(k)^(k)(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)) π̃_ℳ_i^(k)^(k)_T_24 + 2π̃_ℳ_i^(k)^(k)T(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)) 𝐗̃_ℳ_i^(k)^(k)Tξ_i^(k)_T_25
           + ξ_i^(k)T𝐗̃_ℳ_i^(k)^(k)(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)) 𝐗̃_ℳ_i^(k)^(k)Tξ_i^(k)_T_26.

	
	Firstly, we will derive the bound for T_24, T_25 and T_26 terms, then we can obtain similar results for term T_21, T_22 and T_23 by simply removing the matrix 𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k). Denote the singular value decomposition 𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k) = U^(k)T_i V^(k)_i U^(k)_i, where U^(k)_i is a unitary matrix, V^(k)_i is a diagonal matrix with eigenvalues v_i. Therefore, the shared component of π̃_ℳ_i^(k)^(k) and 𝐗̃_ℳ_i^(k)^(k) can be rewritten as
	
    (𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k) +λ^(k)_iI_d)^-1 = U^(k)T_i(V^(k)_i + λ^(k)_i I_d)^-1U^(k)_i.

	
	By Assumption 3, there are some constants c_1,c_2 such that max_δ=1 (n^(k))^-1/2𝐗_ℳ_i^(k)^(k)δ≤ c_1 and min_δ=1 (n^(k))^-1/2𝐗_ℳ_i^(k)^(k)δ≥ c_2. Thus, λ_max(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)) < c_1^2 n^(k) and λ_min(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)) > c_2^2 n^(k). That is, v_j ≍ n^(k) for each eigenvalue. Let b = U^(k)_iπ_ℳ_i^(k)^(k), then b = π_ℳ_i^(k)^(k). Noting that λ^(k)_i = o(n^(k)) in Assumption 3, we can bound the term T_24 as follows,
	
    T_24   = π̃_ℳ_i^(k)^(k)T(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)) π̃_ℳ_i^(k)^(k) = λ^(k)2_i b^T V^(k)_i(V^(k)_i+λ^(k)_iI_d)^-1V^(k)_i(V^(k)_i+ λ^(k)_iI_d)^-1b
       = λ^(k)2_i∑_j=1^dv_j b_ij^2/(v_j+λ^(k)_i)^2 = λ^(k)2_iπ^(k)_ℳ_i^(k)^2/n^(k) = r^(k)_i.

	Similarly, removing the term 𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k), we have
	
    T_21 = λ^(k)2_iπ_ℳ_i^(k)^(k)^2/n^(k) = r_i^(k)/n^(k).

	
	Noting that T_25 follows a Gaussian distribution, we can write the probability of deviation of T_25 with the classical Gaussian tail inequality, for any positive number t,
	
    ℙ( T_25≤ t ) ≥ 1- exp(-1/2 t^2/(T_25)).

	Furthermore,
	
    (T_25)    =    4σ̃^(k)2_iπ̃^(k)T_(i)(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)) 𝐗̃_ℳ_i^(k)^(k)T𝐗̃_ℳ_i^(k)^(k)(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)) π̃_ℳ_i^(k)^(k)
       =    4σ̃^(k)2_iλ^(k)2_i b^T (V +λ^(k)_iI_d)^-1V^(k)_i(V^(k)_i+λ^(k)_iI_d)^-1
          × V^(k)_i(V^(k)_i+λ^(k)_iI_d)^-1V^(k)_i(V^(k)_i+λ^(k)_iI_d)^-1b 
       =    4σ̃^(k)2_iλ^(k)2_i∑_j=1^dv_j^3 b_ij^2/(v_j+λ^(k)_i)^4=σ̃^(k)2_iλ^(k)2_iπ_ℳ_i^(k)^(k)^2/n^(k) = σ̃^(k)2_i r_i^(k).

	Letting t = √(2 (T_25) (f^(k) +log 2)), we obtain that, with probability at least 1- e^-f^(k)/2,
	
    T_25 = √(r_i^(k) f^(k)).
	
	Similarly, removing 𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k), we can obtain that, concurring with (<ref>),
	
    T_22 = √(r_i^(k) f^(k))/n^(k).

	
	The term T_26 follows a non-central χ^2 distribution. We can invoke the Hanson-Wright inequality [Rudelson et al., 2013] to bound the probability of its extreme deviation, for some constant t_2>0, 
	
    ℙ(T_26≤T_26 + t)
       ≥    1- exp{-t^2 t_2/σ̃^(k)4_i𝐗̃_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)T^2∧-t t_2/σ̃^(k)2_i𝐗̃_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)T}.

	To understand this probabilistic bound, we need to calculate T_26 and the two involved norms. Firstly,
	
    T_26   =   σ̃^(k)2_i(𝐗̃_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)T) 
       =   σ̃^(k)2_i(V^(k)_i(V^(k)_i+λ^(k)_iI_d)^-1V^(k)_i(V^(k)_i+λ^(k)_iI_d)^-1) 
       =   σ̃^(k)2_i∑_j=1^dv_j^2/(v_j+λ^(k)_i)^2 = d σ̃^(k)2_i.

	The Frobenius norm can be simplified as follows,
	
    𝐗̃_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)T^2
       =   (𝐗̃_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)T𝐗̃_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)T)  
       =   (((𝐗_ℳ_i^(k)^(k))^T𝐗_ℳ_i^(k)^(k)) (𝐗̃_ℳ_i^(k)^(k))^T 𝐗̃_ℳ_i^(k)^(k) ((𝐗_ℳ_i^(k)^(k))^T𝐗_ℳ_i^(k)^(k)) (𝐗̃_ℳ_i^(k)^(k))^T𝐗̃_ℳ_i^(k)^(k))  
       =   (V^(k)_i(V^(k)_i+λ^(k)_iI_d)^-1 V^(k)_i(V^(k)_i+λ^(k)_iI_d)^-1V^(k)_i (V^(k)_i+λ^(k)_iI_d)^-1V^(k)_i(V^(k)_i+λ^(k)_iI_d)^-1) 
       =   ∑_j=1^dv_j^4/(v_j+λ^(k)_i)^4 = d.

	Note that λ_max(𝐗_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T) ≍ n^(k), then, the operator norm can be simplified as follows,
	
    𝐗̃_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)T
       =   𝐗_ℳ_i^(k)^(k)(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)+λ^(k)_iI_d)^-1𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)(𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)+λ^(k)_iI_d)^-1𝐗_ℳ_i^(k)^(k)T
       =   λ_max(𝐗_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T)/n^(k)2 = 1.

	Letting
	
    t =    √(σ̃^(k)4_i𝐗̃_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)T^2× (f^(k)+log 2 )/t_2)
       ( σ̃^(k)2_i𝐗̃_ℳ_i^(k)^(k)𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k)𝐗̃_ℳ_i^(k)^(k)T× (f^(k) + log 2)/t_2 ),

	and combining (<ref>), (<ref>), (<ref>), and (<ref>), we obtain that, with probability at least 1-e^-f^(k)/2,
	
    T_26 = d √(d f^(k)) f^(k).

	Similarly, removing 𝐗_ℳ_i^(k)^(k)T𝐗_ℳ_i^(k)^(k), we can obtain that, concurring with (<ref>),
	
    T_23 = (d √( d  f^(k)) f^(k)) /n^(k).

	Collecting the bounds (<ref>), (<ref>), (<ref>) and noting the definition of 𝐗_ℳ_i^(k)^(k) and π_ℳ_i^(k)^(k),  we conclude there exists some constant C_2^(k)>0 such that, with probability at least 1-e^-f^(k),
	
    1/n^(k)𝐗^(k)( π̂^(k)_i- π^(k)_i  ) ^2=1/n^(k)𝐗_ℳ_i^(k)^(k)( π̂_ℳ_i^(k)^(k) - π_ℳ_i^(k)^(k) ) ^2 ≤ C_2^(k)r^(k)_i d  f^(k)/n^(k).

	Similarly, collecting the bound (<ref>), (<ref>) and (<ref>), we conclude there exists some constant C_1^(k)>0 such that, with probability at least 1-e^-f^(k),
	
    π̂^(k)_i - π^(k)_i ^2=
    		π̂_ℳ_i^(k)^(k) - π_ℳ_i^(k)^(k)^2 ≤ C_1^(k)r_i^(k) d  f^(k)/n^(k).

	This concludes the proof of Lemma <ref>.



To bound the estimation loss, we write

    Π̂_j - Π_j^2  = π̂^(1)_j|p - π^(1)_j|p^2 + π̂^(2)_j|p - π^(2)_j|p^2,

where π^(k)_j|p and π̂^(k)_j|p are the j|p columns of π^(k) and π̂^(k), respectively. Following the bounds in Lemma <ref> for both networks, we obtain the overall estimation bound as, with probability at least 1 - e^-f^(1) - e^-f^(2),

    Π̂_j - Π_j^2    ≤ C_1^(1)r_j|p^(1) d  f^(1)/n^(1) + C_1^(2)r_j|p^(2) d  f^(2)/n^(2)
       ≤(C_1^(1) + C_1^(2)) (r_j|p^(2) d  f^(2)) (r_j|p^(2) d  f^(2)) /n^(1) n^(2)
        = C_1 d(r_j|p^(1) r_j|p^(2)) (f^(1) f^(2))/n^(1) n^(2)≤ C_1 d r_max f_max/n^(1) n^(2),

where C_1 = C_1^(1) + C_1^(2). Similarly, we write the prediction bound as, with probability at least 1 - e^-f^(1) - e^-f^(2),

    𝐗(Π̂_j - Π_j)^2     = X^(1)(π̂^(1)_j|p - π^(1)_j|p  ) ^2 + X^(2)(π̂^(2)_j|p - π^(2)_j|p )^2
       ≤  C_2^(1){ r_j|p^(1) d  f^(1) + C_2^(2) r_j|p^(2) d  f^(2)}
       ≤ C_2 { d(r_j|p^(1) r_j|p^(2)) (f^(1) f^(2)) }≤ C_2 {d  r_max f_max},

where C_2 = C_2^(1) + C_2^(2) and r_max = 1≤ i ≤ pmax (r_i^(1) r_i^(2)). This concludes the proof of Theorem 2.





§ PROOF OF THEOREM 3
 

Let c_max = c_1^(1) c_1^(2), and further denote

    g_n= C_2 d r_max f_max/n+ 2c_max C_2 Π√( d r_max f_max/n).


 
	Suppose that, for node i,
	
    √((d  r_max f_max) / n) + c_maxΠ≤√(c_max^2Π^2+ϕ_0^2/(64C_2 |𝒮_i|)).

	Under Assumptions 1-3, we have 𝐇_i 𝐗Π̂_-i𝒮_i≥ϕ_0/2 with probability at least 1- e^-f^(1)+log p - e^-f^(2)+log p.



	The inequality (<ref>) implies that g_n≤ϕ_0^2/(64|𝒮_i|).
	
	For any index set 𝒮_i and vector δ, note the definition of ϕ_re(·), then, we have that δ^2 ≤ (δ_𝒮_i^c + δ_𝒮_i)^2 ≤ (3√(|𝒮_i|)δ_𝒮_i + √(|𝒮_i|)δ_𝒮_i)^2 = 16 |𝒮_i| δ_𝒮_i^2. we also have
	
    δ^T((𝐇_i𝐗Π̂_-i)^T  (𝐇_i𝐗Π̂_-i) -  (𝐇_i𝐗Π_-i)^T (𝐇_i𝐗Π_-i) ) δ/nδ_𝒮_i^2
       ≤   δ^2/nδ_𝒮_i ^2  j_1,j_2max| (𝐇_i𝐗Π̂_j_1)^T  (𝐇_i𝐗Π̂_j_2) -  (𝐇_i𝐗Π_j_1 )^T (𝐇_i𝐗Π_j_2 )| 
       ≤   16 |𝒮_i|/nj_1,j_2max| (𝐇_i𝐗Π̂_j_1)^T  (𝐇_i𝐗Π̂_j_2) -  (𝐇_i𝐗Π_j_1 )^T (𝐇_i𝐗Π_j_2)|.

	Note that,
	
    (𝐇_i𝐗Π̂_j_1)^T  (𝐇_i𝐗Π̂_j_2) -  (𝐇_i𝐗Π_j_1 )^T (𝐇_i𝐗Π_j_2)
       =   (Π̂_j_1 -  Π_j_1 )^T𝐗^T𝐇_i𝐗 (Π̂_j_2 -  Π_j_2)_T_31 + (Π̂_j_1 -  Π_j_1 )^T𝐗^T𝐇_i𝐗Π_j_2_T_32 + (𝐗Π_j_1 )^T𝐇_i𝐗 (Π̂_j_2 -  Π_j_2)_T_33.

	We will derive the bounds for each of these three terms separately. With 𝐇_i a projection matrix, we have λ_max(𝐇_i) = 1. We can obtain that
	
    |T_31|    ≤   𝐇_i𝐗 (Π̂_j_1 - Π_j_1 )×𝐇_i𝐗 (Π̂_j_2 -  Π_j_2 )
       ≤   λ_max(𝐇_i)  𝐗 (Π̂_j_1 - Π_j_1)×𝐗 (Π̂_j_2 - Π_j_2 )
       =   𝐗 (Π̂_j_1 - Π_j_1 )×𝐗 (Π̂_j_2 - Π_j_2 ).

	Note that
	|T_32| ≤𝐗 Π_j_2𝐇_i𝐗 (Π̂_j_1 -  Π_j_1), and following Assumption 3, we have
	
    𝐗 Π_j_2^2    =   X^(1)π^(1)_j|p^2 + X^(2)π^(2)_j|p^2
       ≤    (c_1^(1))^2 n^(1)π^(1)_j|p^2 + (c_1^(2))^2 n^(2)π^(2)_j|p^2 
       ≤    c_max^2 n (π^(1)_j|p^2 +π^(2)_j|p^2)
       ≤    c_max^2 n (π^(1)_j|p +π^(2)_j|p)^2
       ≤    c_max^2 nΠ^2.

	Therefore,
	
    |T_32| ≤𝐗Π_j_2𝐇_i𝐗 (Π̂_j_1 -  Π_j_1 )≤ c_max√(n)Π𝐗 (Π̂_j_1 - Π_j_1 ).

	Similarly, we can have
	
    |T_33| ≤ c_max√(n)Π𝐗 (Π̂_j_2 - Π_j_2).

	Theorem 2 leads to the following, with probability at least 1- e^-f^(1)+log(p) - e^-f^(2)+log(p),
	
    |T_31|/n   ≤ C_2  d r_max f_max/n,
    |T_32|/n   ≤ c_maxC_2 Π√( d r_max f_max/n),
    |T_33|/n   ≤ c_maxC_2 Π√( d r_max f_max/n).

	Putting the above three inequalities together, we have,
	
    δ^T( (𝐇_i𝐗Π̂_-i)^T  (𝐇_i𝐗Π̂_-i) -  (𝐇_i𝐗Π_-i)^T (𝐇_i𝐗Π_-i) ) δ/nδ_𝒮_i^2 
       ≤     16 |𝒮_i| ×|T_31| + |T_32| + |T_33|/n = 16 |𝒮_i| g_n≤ 16 |𝒮_i| ϕ_0^2/64 |𝒮_i| = ϕ_0^2/4.

	Together with Assumption 4, we have 𝐇_i𝐗Π̂_-k𝒮_k≥ϕ_0/2. This concludes the proof of Lemma <ref>.



(Basic Inequality) 
	Let η_i =  2n^-1𝐙̂_-i^T 𝐇_i ϵ_i - 2 n^-1𝐙̂_-i^T 𝐇_i (𝐙̂_-i - 𝐙_-i)  β_i and
	
    ℰ(λ_i) = {W_i^-1η_i≤λ_i/2 },

	for λ_i specified in Theorem 3. Under Assumptions 1-2, with h_n defined in Theorem 3, there exit a positive constant C_3 > 0 such that
	
    ℙ(ℰ(λ_i)) ≥ 1-e^-C_3 h_n +log(4q) - e^-f^(1)+log(p) - e^-f^(2) + log(p).

	Concurring with event ℰ(λ_i), we have the following basic inequality,
	
    n^-1𝐇_i𝐙̂_-i (β̂_i - β_i ) ^2 + λ_i_i^T |β̂_i|_1 ≤λ_i_i^T |β_ i |_1 +η_i^T (β̂_i - β_i).




	Letting
	
    ξ_-i =
    		[  ξ_-i^(1)  ξ_-i^(1);  ξ_-i^(2) -ξ_-i^(2) ],

	we have 𝐙_-i = 𝐗Π_-i +ξ_-i. With 𝐙̂_-i = 𝐗Π̂_-i, we get
	
    η_i   =   2/nΠ̂_-i^T 𝐗^T 𝐇_iϵ_i - 2/nΠ̂_-i^T 𝐗^T 𝐇_i (𝐗Π̂_-i - 𝐗Π_-i -ξ_-i) β_i
       =   2/n (Π̂_-i - Π_-i)^T 𝐗^T𝐇_iϵ_i_T_34 + 2/nΠ_-i^T𝐗^T𝐇_iϵ_i_T_35
          + 2/n (Π̂_-i - Π_-i)^T𝐗^T 𝐇_i ξ_-iβ_i _T_36+ 2/nΠ_-i^T 𝐗^T 𝐇_i ξ_-iβ_i _T_37
          - 2/n(Π̂_-i - Π_-i)^T𝐗^T𝐇_i𝐗(Π̂_-i -Π_-i) β_i _T_38 - 2/nΠ_-i^T𝐗^T𝐇_i𝐗(Π̂_-i-Π_-i) β_i _T_39.

	We aim to bound each of these six terms by λ_i/12 either probabilistically or deterministically.
	
	Firstly, for some constant t_λ> 0, we choose the adaptive lasso tuning parameter as below,
	
    λ_i = t_λω_i_-∞^-1𝐁Π√((d  r_max f_max)log(p)/n_min).

	Denoting the j-th column of 𝐗 by X_· j, we have X_· j^T X_· j = n^(k) for k∈{1,2} due to standardization. Furthermore,
	
    (1/n X_· j^T 𝐇_i ϵ_i) ≤1/n^2 X_· j^T 𝐇_i X_· jσ_pmax^2≤n^(k)/n^2σ_pmax^2 ≤1/nσ_pmax^2.

	
	For T_34, via the classical Gaussian tail inequality, we have
	
    ℙ(W_i^-1 T_34≥λ_i/12)    ≤   ℙ(  2/n(Π̂_-i  - Π_-i)^T  𝐗^T𝐇_iϵ_i≥λ_iω_i_-∞/12) 
       ≤   ℙ( (Π̂_-i  - Π_-i)^T 2/n𝐗^T𝐇_iϵ_i≥λ_iω_i_-∞/12) 
       ≤   ℙ(2/n𝐗^T𝐇_iϵ_i≥λ_iω_i_-∞/12 δ_Π) ≤ 2q exp{- n λ_i^2ω_i_-∞^2/1152σ^2_pmaxδ_Π^2}
       ≤    2q· p^- n/d t_1 𝐁^2Π^2 ≤ 2q· p · p^- t_1 𝐁^2 n/dΠ^2,

	where t_1 = t^2_λ/(2304C_1 σ^2_pmax), and δ_Π is the maximum estimation loss of the first stage. The last inequality is obtained based on the following bound of δ_Π. Following Theorem 2, δ_Π satisfies the following inequality with probability at least 1 - e^-f^(1)+log(p) - e^-f^(2) + log(p),
	
    δ_Π^2 = 1≤ j ≤ 2pmaxΠ̂_j - Π_j^2 ≤1≤ j ≤ 2pmax(2dΠ̂_j - Π_j^2 )
    			≤ 2C_1 d {d r_max f_max/n_min}.

	Note that the first inequality of (<ref>) holds, since Π̂ and Π have at most 2d non-zeros based on our assumptions and the screening in the calibration step.
	
	Similarly, for the second term T_35, we have that, with t_2 = (t_λ)^2/1152σ_pmax^2,
	
    ℙ(W_i^-1 T_35≥λ_i/12)    ≤   ℙ(2/nΠ_-i^T  𝐗^T𝐇_iϵ_i≥λ_iω_i_-∞/12) 
       ≤   ℙ( Π_-i^T2/n𝐗^T𝐇_iϵ_i≥λ_iω_i_-∞/12) 
       ≤   ℙ( 2/n𝐗^T𝐇_iϵ_i≥λ_iω_i_-∞/12 Π_-i^T ) ≤ 2q exp{ - nλ_i^2ω_i_-∞^2/1152σ^2_pmaxΠ_-i^T^2}
       =    2q · p^- t_2 𝐁^2 (d  r_max f_max) n/n_min≤ 2q · p · p^- t_2 𝐁^2 (d  r_max f_max) n/n_min.

	For the third term T_36, we write
	
    ℙ(W_i^-1 T_36≥λ_i/12)    ≤   ℙ( (Π̂_-i -Π_-i)^T2/n𝐗^T𝐇_i ξ_-iβ_i ≥λ_iω_i_-∞/12) 
       ≤   ℙ(δ_Π×j_1,j_2max |2/n X_· j_1^T𝐇_iξ_j_2| ×β_i≥λ_iω_i_-∞/12) 
       ≤   ℙ(j_1,j_2max |2/n X_· j_1^T𝐇_iξ_j_2 | ≥λ_iω_i_-∞/12δ_Πβ_i) ≤ 2q · 2p exp{-nλ_i^2ω_i_-∞^2/1152 σ̃_pmax^2 δ_Π^2β_i^2}
       =    4q· p · p^-t_3 Π^2 n/d,

	where σ̃_pmax^2 = imax(σ̃_i^(1)σ̃_i^(2)), (1/n X_· j_1^T𝐇_iξ_j_2 ) ≤σ̃_pmax^2/n and  t_3 = t_λ^2/2304 C_1σ̃_pmax^2. Similarly, with t_4 = t_λ^2/1152σ̃_pmax^2, we write T_37 term as
	
    ℙ(  W_i^-1 T_37≥λ_i/12)    ≤    2q · 2p ·exp{-nλ_i^2 ω_i_-∞^2/1152σ̃_pmax^2Π_-i^T^2β_i^2}
       =    4q· p · p ^- t_4 (d  r_max f_max) n/n_min.

	For the deterministic term T_38, choosing t_λ≥ 12C_2Π^-1√((d r_max f_max)/(nlog(p))), along with Cauchy-Schwarz Inequality, we have
	
    W_i^-1 T_38   ≤   β_iω_i_-∞^-1/nj_1,j_2max | (Π̂_j_1 - Π_j_1)^T𝐗^T𝐇_i𝐗(Π̂_j_2 -Π_j_2)|
       ≤   β_iω_i_-∞^-1/nj_1,j_2max{𝐇_i𝐗(Π̂_j_1 -Π_j_1)𝐇_i𝐗(Π̂_j_2 -Π_j_2)}
       ≤   β_iω_i_-∞^-1/nj_1,j_2max{λ_max(𝐇_i) 𝐗(Π̂_j_1 -Π_i_1)𝐗(Π̂_j_2 -Π_j_2)}
       ≤   β_iω_i_-∞^-1/nj_1,j_2max{𝐗(Π̂_j_1 -Π_j_1)𝐗(Π̂_j_2 -Π_j_2)}
       ≤   β_iω_i_-∞^-1 C_2 d r_max f_max/n≤λ_i/12×( 12C_2/ t_λΠ√( d  r_max f_max/n log(p)))  ≤λ_i/12.

	Similarly, we choose t_λ≥ 24 √(C_2 n_min/(nlog(p))), and take Theorem 2 to obtain
	
    W_i^-1 T_39   ≤    2β_iΠ_-i^Tω_i_-∞^-1/nj_1,j_2max | X_· j_1^T𝐇_i𝐗(Π̂_j_2 -Π_j_2)|  
       ≤    2β_iΠ_-i^Tω_i_-∞^-1/√(n)j_2max𝐇_i𝐗(Π̂_j_2 -Π_j_2)
       ≤    2β_iΠ_-i^Tω_i_-∞^-1/√(n)j_2max𝐗(Π̂_j_2 -Π_j_2)≤λ_i/12×(  24/ t_λ√(C_2 n_min/nlog(p))) ≤λ_i/12.

	Note that n≥ n_min. Putting together the probabilistic bounds (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>), along with union bound, there exist a constant C_3 > 0 such that
	
    ℙ(ℰ(λ_i)) ≥ 1-3e^-C_3 h_n +log(4pq) - e^-f^(1)+log(p) - e^-f^(2) + log(p).

	
	Next we will establish the basic inequality, concurring with the event ℰ(λ_i).
	
	Since the estimator β̂_i from the adaptive lasso minimizes the corresponding objective function, we have
	
    1/n𝐇_i𝐘_i -𝐇_i 𝐙̂_-iβ̂_i  + λ_i_i^T |β̂_i|_1 
    ≤1/n𝐇_i𝐘_i -𝐇_i 𝐙̂_-iβ_i    +  λ_i_i^T  |β_i |_1.

	Because 𝐇_i𝐘_i = 𝐇_i 𝐙_-iβ_i   +𝐇_iϵ_i, we can rewrite
	
    𝐇_i𝐘_i -𝐇_i 𝐙̂_-iβ̂_i ^2
       =   𝐇_i 𝐙_-iβ_i   +𝐇_iϵ_i -𝐇_i 𝐙̂_-iβ̂_i ^2  
       =   𝐇_iϵ_i^2 - 2 ϵ_i^T𝐇_i(  𝐙̂_-iβ̂_i -𝐙_-iβ_i   )+ 𝐇_i 𝐙̂_-iβ̂_i -𝐇_i𝐙̂_-iβ_i + 𝐇_i𝐙̂_-iβ_i   -𝐇_i 𝐙_-iβ_i ^2  
       =   𝐇_iϵ_i^2 - 2 ϵ_i^T𝐇_i(  𝐙̂_-iβ̂_i -𝐙_-iβ_i   )+ 𝐇_i𝐙̂_-i (β̂_i -  β_i  )^2 +𝐇_i(𝐙̂_-i -𝐙_-i) β_i ^2  
          +  2 β_i  ^T (𝐙̂_-i -𝐙_-i)^T𝐇_i𝐙̂_-i(β̂_i - β_i ).

	Similarly we can rewrite
	
    𝐇_i𝐘_i - 𝐇_i𝐙̂_-iβ_i ^2    =   𝐇_i𝐙_-iβ_i   +𝐇_iϵ_i -𝐇_i𝐙̂_-iβ_i ^2  
       =   𝐇_iϵ_i^2 + 𝐇_i(𝐙̂_-i-𝐙_-i) β_i ^2 - 2ϵ_i^T𝐇_i(𝐙̂_-i -𝐙_-i) β_i.

	Plugging equations (<ref>) and (<ref>) into (<ref>), we then have
	
    1/n𝐇_i𝐙̂_-i (β̂_i - β_i  ) ^2 +λ_i_i^T |β̂_i|_1
       ≤   λ_i_i^T |β_i|_1 +(2/n𝐙̂_-i^T 𝐇_iϵ_i - 2/n𝐙̂_-i^T 𝐇_i (𝐙̂_-i - 𝐙_-i)  β_i )^T(β̂_i - β_i  )  
       =   λ_i_i^T |β_i|_1 + η_i ^T(β̂_i - β_i  ).

	Thus, the basic inequality is established. This concludes the proof of Lemma <ref>.


Conditioning on the event ℰ(λ_i), we remove the random term η_i from the basic inequality as

    1/n𝐇_i𝐙̂_-i (β̂_i - β_i)^2
       ≤   λ_i_i^T |β_ i|_1 - λ_i_i^T |β̂_i|_1  + η_i^T (β̂_i - β_i) 
       ≤   λ_iω^T_𝒮_i |β_𝒮_i|_1 -
    	λ_iω^T_𝒮_i |β̂_𝒮_i|_1 -λ_iω^T_𝒮_i^c |β̂_𝒮_i^c|_1 + η_𝒮_i^c^T(β̂_𝒮_i^c) +  η_𝒮_i^T(β̂_𝒮_i -β_𝒮_i) 
       ≤   λ_iω^T_𝒮_i |β̂_𝒮_i -
    	β_𝒮_i|_1 -λ_iω^T_𝒮_i^c |β̂_𝒮_i^c|_1 + λ_i/2ω^T_𝒮_i^c|β̂_𝒮_i^c|_1 +  λ_i/2ω^T_𝒮_i |β̂_𝒮_i -β_𝒮_i|_1  
       ≤   3/2λ_iω^T_𝒮_i |β̂_𝒮_i -
    	β_𝒮_i|_1- 1/2λ_iω^T_𝒮_i^c |β̂_𝒮_i^c|_1 
       ≤   3/2λ_iω_𝒮_i_∞β̂_𝒮_i - β_𝒮_i - 1/2λ_iω_𝒮_i^c_-∞β̂_𝒮_i^c.

The fact that 𝐇_i𝐙̂_-i (β̂_i - β_i)^2 is always positive leads to

    ω_𝒮_i^c_-∞β̂_𝒮_i^c≤ 3ω_𝒮_i_∞β̂_𝒮_i -
    		β_𝒮_i,

which, following Assumption 4, further implies that

    β̂_𝒮_i^c -β_𝒮_i^c≤ 3 β̂_𝒮_i -
    		β_𝒮_i.


The above inequality, as well as the last inequality in (<ref>), implies that

    1/n𝐇_i𝐙̂_-i (β̂_i - β_i)^2
       ≤   3/2λ_iω_𝒮_i_∞β̂_𝒮_i -
    		β_𝒮_i≤3/2λ_iω_𝒮_i_∞√(|𝒮_i|)β̂_𝒮_i -
    		β_𝒮_i
       ≤   3/2λ_iω_𝒮_i_∞√(|𝒮_i|)2 𝐇_i𝐙̂_-i (β̂_i - β_i)/√(n)ϕ_0,

where the last inequality follows Assumption 4 and Lemma <ref>. The above inequality leads to that,

    1/n𝐇_i𝐙̂_-i (β̂_i - β_i)^2 ≤9 ( ω_𝒮_i_∞ )^2/ϕ_0^2|𝒮_i| λ_i^2.

Plugging in (<ref>), and letting C_4 = 3t_λ, we obtain that

    1/n𝐇_i𝐙̂_-i (β̂_i - β_i)^2 ≤ C_4^2 ω_𝒮_i_∞^2 𝐁^2Π^2/ϕ_0^2ω_i_-∞^2   |𝒮_i|  (d  r_max f_max)log(p)/n_min.


Taking this inequality, we can follow Assumption 4 and Lemma <ref> to derive that

    β̂_i -β_i   =   β̂_𝒮_i^c+β̂_𝒮_i - β_𝒮_i≤(3ω_𝒮_i_∞/ω_𝒮_i^c_-∞ +1) β̂_𝒮_i - β_𝒮_i
       ≤   (3ω_𝒮_i_∞/ω_𝒮_i^c_-∞  +1)√(|𝒮_i|)2 𝐇_i𝐙̂_-i(β̂_i - β_i)/√(n)ϕ_0
       ≤   (3ω_𝒮_i_∞/ω_𝒮_i^c_-∞  +1)√(|𝒮_i|)2C_4 ω_𝒮_i_∞𝐁Π/ϕ_0^2 ω_i_-∞√(|𝒮_i|)√((d  r_max f_max)log(p)/n_min)
       ≤    8C_4 ω_𝒮_i_∞𝐁Π/ϕ_0^2 ω_i_-∞ |𝒮_i| √((d  r_max f_max)log(p)/n_min),

where the last inequality follows Assumption 4. Since the inequality (<ref>) concurs with the event ℰ(λ_i), the above prediction and estimation bounds hold with probability at least 1- 3e^-C_3 h_n + log(4pq) - e^-f^(1)+log(p) - e^-f^(2) + log(p). This completes the proof of Theorem 3.





§ PROOF OF THEOREM 4
 

 
	Suppose that, for node i,
	
    √((d  r_max f_max) / n) + c_maxΠ≤√(c_max^2Π^2+min(ϕ_0^2/64, τ(4-τ)^-1ω_i_-∞/ψ_i)/(C_2 |𝒮_i|)).

	Under Assumptions 1-5, we have W^-1_𝒮^c_i([-1])W_𝒮_i≤ 1- τ/2 with the probability at least 1- e^-f^(1)+log(p)- e^-f^(2)+log(p).



	The inequality (<ref>) implies that ψ_iω_i_-∞^-1 |𝒮_i| g_n ≤τ/4 - τ.
	
	By the inequalities (<ref>) and (<ref>) in the proof of Lemma <ref> and union bound, we have that, with probability at least 1- e^-f^(1)+log(p)- e^-f^(2)+log(p),
	
    j_1,j_2max{1/n|(𝐇_i𝐗Π̂_j_1)^T  (𝐇_i𝐗Π̂_j_2) -  (𝐇_i𝐗Π_j_1 )^T (𝐇_i𝐗Π_j_2 )| }≤ g_n.

	With the definitions of infinity norm ·, , and , we can obtain the following inequality indexed by set 𝒮_i,
	
    ψ_iW_𝒮_i^-1( - )≤ψ_iω_𝒮_i_-∞^-1 - ≤ψ_iω_𝒮_i_-∞^-1 |𝒮_i| g_n ≤τ/4 - τ.

	Similarly we can obtain the following bound indexed by the complement set 𝒮^c_i,
	
    ψ_iW_𝒮^c_i^-1([] - )≤ψ_iω_𝒮_i^c_-∞^-1 |𝒮_i| g_n ≤τ/4 - τ.
	
	Applying the matrix inversion error bound in Horn and Johnson [2012] and the triangular inequality, we have that 
	
    [-1] W_𝒮_i   ≤   [-1] W_𝒮_i + [-1] W_𝒮_i -[-1]W_𝒮_i
       ≤   ψ_i  + ψ_iW_𝒮_i^-1([] - )/1- ψ_iW_𝒮_i^-1([] - )ψ_i≤ψ_i + τ/4-2τψ_i≤4-τ/4-2τψ_i.

	
	Also note that we can rewrite
	
    W^-1_𝒮_i^c([-1]  - [-1])W_𝒮_i
       =    W^-1_𝒮_i^c( - ) [-1] W_𝒮_i + W^-1_𝒮_i^c[-1] W_𝒮_i W_𝒮_i^-1( - ) [-1] W_𝒮_i.

	Then, it follows from (<ref>), (<ref>), (<ref>) and Assumption 5 that
	
    W^-1_𝒮_i^c(  [-1]  - [-1]  )W_𝒮_i
       ≤   W^-1_𝒮_i^c(  - )[-1] W_𝒮_i
           + W^-1_𝒮_i^c[-1]W_𝒮_iW_𝒮_i^-1(-)[-1] W_𝒮_i≤τ/2.

	Therefore, together with Assumption 5 again, we can conclude that W^-1_𝒮_i^c([-1]) W_𝒮_i≤ 1- τ/2.
	
	This concludes the proof of Lemma <ref>.


The optimality of β̂_i in the adaptive lasso step and KKT condition lead to

    -2/n (𝐇_i𝐙̂_-i)^T(𝐇_i𝐘_i - 𝐇_i𝐙̂_-iβ̂_i) + λ_i W_iα_i =0,

where α_i∈ℝ^2p-2, satisfying that α_i≤ 1 and α_ijI(β̂_ij 0) = sign(β̂_ij).

Plug in the equation 𝐇_i_i = 𝐇_i𝐙_-iβ_i + 𝐇_iϵ_i, we can have that

    𝐇_i𝐘_i - 𝐇_i𝐙̂_-iβ̂_i   =   𝐇𝐙_-iβ_i +𝐇_iϵ_i - 𝐇_i𝐙̂_-iβ̂_i
       =   𝐇_iϵ_i + 𝐇_i𝐙_-iβ_i - 𝐇_i𝐙̂_-iβ_i + 𝐇_i𝐙̂_-iβ_i - 𝐇_i𝐙̂_-iβ̂_i
       =   𝐇_iϵ_i - 𝐇_i(𝐙̂_-i - 𝐙_-i)β_i -𝐇_i𝐙̂_-i(β̂_i-β_i).

This, along with KKT condition (<ref>), leads to

    2ℐ̂_i(β̂_i - β_i) -η_i = -λ_i W_iα_i,

where η_i is defined in Lemma <ref>.

Letting β̂_𝒮_i^c = β_𝒮_i^c = 0, equation (<ref>) can be decomposed as

    2ℐ̂_i,11(β̂_𝒮_i - β_𝒮_i) -η_𝒮_i   = -λ_iW_𝒮_iα_𝒮_i, 
    
    			2ℐ̂_i,21(β̂_𝒮_i - β_𝒮_i) -η_𝒮_i^c    = -λ_iW_𝒮_i^cα_𝒮_i^c.

We can solve for β̂_𝒮_i from the first equation of (<ref>) as

    β̂_𝒮_i - β_𝒮_i   =    2^-1ℐ̂_i,11^-1(η_𝒮_i -λ_iW_𝒮_i^T α_𝒮_i) = 2^-1ℐ̂_i,11^-1 W_𝒮_i   (W_𝒮_i^-1η_𝒮_i -λ_iα_𝒮_i).


Following the similar strategy in the proof of Lemma <ref>, we can prove that there exists a constant C_5 > 0 such that W_i^-1η_i≤τ/4-τλ_i with probability at least 1- 3e^-C_5 h_n + log (4q) + log (p)- e^-f^(1)+log (p) - e^-f^(2) + log (p). Thus, together with α_𝒮_i≤ 1, we obtain the infinity norm estimation loss on the true support set 𝒮_i

    β̂_𝒮_i - β_𝒮_i   ≤    2^-1ℐ̂_i,11^-1 W_𝒮_i ( W_𝒮_i^-1η_𝒮_i + λ_i  ) 
       ≤    2^-14-τ/4-2τψ_i4/4-τλ_i = λ_iψ_i/2-τ≤j∈𝒮_imin |β_ij|=b_i,

where the last inequality comes from the condition on the minimal signal strength b_i. The above inequality implies sign( β̂_𝒮_i) = sign( β_𝒮_i).

Plugging (<ref>) into the left hand side of the second equation in (<ref>), we can verify that

    W_𝒮_i^c^-1ℐ̂_i,21 (ℐ̂_i,11)^-1(η_𝒮_i-λ_iW_𝒮_iα_𝒮_i) -  W_𝒮_i^c^-1η_𝒮_i^c
       ≤    W_𝒮_i^c^-1ℐ̂_i,21ℐ̂_i,11^-1 W_𝒮_i ( W^-1_𝒮_iη_𝒮_i +λ_i) +  W_𝒮_i^c^-1η_𝒮_i^c
       ≤    (1-τ/2)(4/(4-τ))λ_i + τ/(4-τ) λ_i = λ_i.

Therefore, we have constructed a solution β̂_i which satisfies the KKT condition (<ref>) and sign(β̂_i) = sign(β_i), that is, 𝒮̂_i = 𝒮_i. This completes the proof of Theorem 4.



§ REFERENCES


	
   Jianqing Fan and Jinchi Lv. Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(5): 849–911, 2008.
	
   Roger A Horn and Charles R Johnson. Matrix Analysis. Cambridge University Press, 2012.
	
    Mark Rudelson, Roman Vershynin, et al. Hanson-wright inequality and sub-gaussian concentration. Electronic Communications
	in Probability, 18, 2013.


