arXiv:1605.07244v1 [stat.ME] 24 May 2016

Optimal Estimation of Co-heritability in
High-dimensional Linear Models
Zijian Guo1 , Wanjie Wang1,2 , T. Tony Cai1 , and Hongzhe Li2
1

Department of Statistics, The Wharton School, University of Pennsylvania
2
Department of Biostatistics and Epidemiology, Perelman School of
Medicine, University of Pennsylvania

Abstract
Co-heritability is an important concept that characterizes the genetic associations within pairs of quantitative traits. There has been significant recent interest
in estimating the co-heritability based on data from the genome-wide association
studies (GWAS). This paper introduces two measures of co-heritability in the highdimensional linear model framework, including the inner product of the two regression vectors and a normalized inner product by their lengths. Functional de-biased
estimators (FDEs) are developed to estimate these two co-heritability measures. In
addition, estimators of quadratic functionals of the regression vectors are proposed.
Both theoretical and numerical properties of the estimators are investigated. In particular, minimax rates of convergence are established and the proposed estimators of
the inner product, the quadratic functionals and the normalized inner product are
shown to be rate-optimal. Simulation results show that the FDEs significantly outperform the naive plug-in estimates. The FDEs are also applied to analyze a yeast
segregant data set with multiple traits to estimate heritability and co-heritability
among the traits.
0

Zijian Guo is a Ph.D. student (Email: zijguo@wharon.upenn.edu). Wanjie Wang is a Postdoctoral
Research Fellow (Email: wanjiew@wharton.upenn.edu). T. Tony Cai is Dorothy Silberberg Professor of
Statistics (E-mail: tcai@wharton.upenn.edu). The research of T. Tony Cai was supported in part by NSF
Grants DMS-1208982 and DMS-1403708, and NIH Grant R01 CA127334. Hongzhe Li is Professor (E-mail:
hongzhe@upenn.edu). The research of Wanjie Wang and Hongzhe Li was supported in part by NIH grants
CA127334 and GM097505.

1

Keywords: Genetic correlations, genome-wide association studies, inner product, minimax
rate of convergence, quadratic functional.

1

Introduction

1.1

Motivation and Background

Genome-wide association studies (GWAS) have not only led to identification of thousands
of genetic variants or single nucleotide polymorphisms (SNPs) that are associated various
complex phenotypes (Manolio, 2010), they also provide important information to estimate
the genetic heritability, a key population parameter that can help understand the genetic
architecture of complex traits (Yang et al., 2010). Results from these GWAS have also shown
that many complex phenotype shared comment genetic variants, including various autoimmune diseases (Zhernakova et al., 2009) and psychiatric disorders (Lee et al., 2013). These
empirical evidence of shared genetic etiology for various phenotypes can inform nosology and
encourages the investigation of common pathophysiologies for related disorders that can be
explored for drug repositioning.
There has been a recent interest in estimating the genetic correlations between two traits
based on the genome-wide aggregate effects of causal variants affecting both traits. The
concept of co-heritability has been proposed to describe the genetic associations within pairs
of quantitative traits based on GWAS data. This is in contrast to the traditional approaches
estimating co-heritability based on twin or family studies. (Baradat, 1976) first generalized
the concept of heritability and defined the coefficient of genetic prediction between two
traits as the ratio of additive genetic covariance over the product of the phenotypic standard
deviation of each trait. This idea was recently extended to define co-heritability using GWAS
data by estimating the genetic risk scores (GRSs) for one phenotype and their association
with a second phenotype (Purcell et al., 2009; Wray et al., 2007). This definition is ad hoc
and does not account for the fact only a very small set of SNPs are expected to be associated
with either of the two traits.
Lee et al. (2012) and Yang et al. (2013) extended the mixed-effect model framework for

2

estimating the heritability to estimate genetic covariance and genetic correlation between
two traits. In their models, each individual’s phenotype is affected by a genetic random
effect, which is correlated across individuals by virtue of sharing some of the genetic variants affecting the traits, and an environmental random effect, which is uncorrelated across
individuals. They define co-heritability as the square-root of the ratio of the covariance of
the genetic random effects to the product of the total variances. The mixed-effect model
approach requires knowledge of the identity of the causal SNPs, and hence the covariance
matrix. This is however not available, Lee et al. (2012) and Yang et al. (2013) approximate
the genetic correlation between every pair of individuals across the set of causal SNPs by the
genetic correlation across the set of all genotyped SNPs. The concept of using an estimated
correlation matrix instead of the actual correlation matrix in estimating the heritability is
a questionable heuristic (Golan and Rosset, 2011). The very large number of SNPs used for
estimating the genetic correlations, most of them likely not causative, might lead to biased
estimate of correlations due the set of causal SNPs. This has been clearly demonstrated by
Golan and Rosset (2011).
The goal of this paper is to define two quantities that can be used for measuring the
co-heritability between a pair of traits based on GWAS data in the framework of highdimensional linear models. Instead of using the mixed-effects model framework, our definition of the co-heritability is similar in spirit to that of Baradat (1976). The proposed
estimation procedure involves selecting the relevant genetic variants, which overcomes the
problem of most of the SNPs in GWAS are not relevant to either of the two traits.

1.2

Problem Formulation

In this paper, a pair of trait values (y, w) are modeled as a linear combination of p genetic
variants and an error term that includes environmental and unmeasured genetic effects,
yn1 ×1 = Xn1 ×p βp×1 + n1 ×1

and wn2 ×1 = Zn2 ×p γp×1 + δn2 ×1 ,

(1)

where the rows Xi· are i.i.d. p-dimensional Sub-gaussian random vectors with mean µx
and covariance matrix Σ, the rows Zi· are i.i.d. p-dimensional Sub-gaussian random vectors
with mean µz and covariance matrix Γ and the error (, δ)| follows the multivariate normal
3

distribution with mean zero and covariance


σ12 In1 ×n1 0n1 ×n2


0n2 ×n1 σ22 In2 ×n2

and is assumed to be independent of X and Z. Throughout the paper there are no particular
assumptions on the relation between the random design matrices X and Z, where they can
be identical, correlated or independent.
In the study of genetic co-heritability, the pairs of traits y and w are assumed to have
mean zero, and the jth column of X, X·j , and the jth column of Z, Z·j , are the numerically
coded genetic markers at the jth genetic variant and are assumed to have variance 1. Under
this model, if the columns of X and Z are independent, for the i-th observation,
Var(yi ) =

X
j

βj2 + σ12 = kβk22 + σ12 ,

and Var(wi ) =

X
j

γj2 + σ22 = kγk22 + σ22 ,

therefore ||β||22 and ||γ||22 can then be interpreted as the narrow sense heritability.
Based on this model, one measure of genetic co-heritability is the inner product of the
regression coefficients
I (β, γ) = hβ, γi

(2)

which measures the shared genetic effects between these two traits. Alternatively, a normalized inner product, that is the ratio
R (β, γ) =

hβ, γi
1(kβk2 kγk2 > 0),
kβk2 kγk2

(3)

can also be used. The denominator kβk2 and kγk2 represent the total genetic variants
for two traits where the columns of the design matrix are independent. In the case where
one of kβk2 and kγk2 is vanishing, the ratio is defined to be zero, which means that there
is no correlation between two traits when one of the regression vector is zero. With this
normalization, R (β, γ) is always between -1 and 1 and can be used to compare co-heritability
among multiple pairs.
The focus of this paper is to develop estimators for I (β, γ) and R (β, γ) based on two
independent GWAS data with genotype data measured on the same set of genetic markers,
denoted by (yi , Xi· , i = 1, · · · , n1 ) and (wi , Zi· , i = 1, · · · , n2 ). A naive estimator is to
4

estimate β and γ first and then plug in the estimators of β and γ into the expressions (2)
and (3). For the problem of interest, usually there are many more genetic markers than the
sample size, that is p  max{n1 , n2 }. However, for any trait, one expect that only a few of
these markers have nonzero effects. One can apply any high-dimensional regression methods
such as Lasso (Tibshirani, 1996) and scaled Lasso (Sun and Zhang, 2012) and marginal
regression with screening (Fan et al., 2012; McCarthy et al., 2008) to estimate these sparse
regression coefficients. However, such plug-in estimators have several drawbacks in terms
of estimating heritability and co-heritability. The Lasso approach shrinks the estimation
towards 0, in particular, some weak effects might be shrunken to 0, yet the accumulation
of these weak effects may contribute significant to the trait variability. It is possible that
some genetic variants may have strong effects on one trait and weak effects on the other
trait. Due to shrinkage, the plug-in of Lasso type estimators fails to capture this part of
contribution to co-heritability from such genetic variants. Marginal regression calculates the
regression score between the trait and each single marker (i.e., y and X·j , 1 ≤ j ≤ p), and
screen for the large scores. This approach also suffers from the existence of weak effects, as
the marginal scores must be large enough to survive in the screening step.

1.3

Methods and Main Results

In this paper, a new estimator, functional de-biased estimator of inner product and normalized inner product (FDE) is proposed to estimate the co-heritability measures I (β, γ)
defined in (2) and R (β, γ) defined in (3). The proposed estimator for I (β, γ) involves two
steps, including estimating the inner product by the plug-in scaled Lasso estimator, and
then correcting the plug-in scaled Lasso estimator. Similar two-step procedures are proposed to estimate the quadratic functionals kβk22 and kγk22 . To estimate the normalized
inner product, we plug-in the estimators of the inner product and quadratic functionals into
the definition (3).
FDE is shown to achieve the minimax optimal convergence rates of estimating I (β, γ)
and R (β, γ). The correction step in FDE improves the estimation accuracy through correcting the shrinkage from Lasso and is crucial to achieve the optimal convergence rate. As

5

demonstrated in the simulation studies, FDE consistently outperforms the plug-in estimator
of the scaled Lasso estimator. In particular, FDE has significant improvement of estimating
the co-heritability measures in the existence of some genetic markers having weak effects
on one trait and strong effects on the other trait. In addition, FDE does not suffer from
dependency among genetic markers. FDE works for a broad class of dependency structure
of genetic markers, and it also allows dependency between two design matrices X and Z.
The theoretical analysis given in Section 3 shows that the optimal rate of convergence
for estimating I (β, γ) is
(kβk2 + kγk2 )



k log p
1
√ +
n
n



+

k log p
,
n

where p is the dimension, n is the sample size, and k is the maximum sparsity of β and γ.
The optimal rate depends not only on p, n and k, but also the signal strength kβk2 and
kγk2 . Restricting to the relatively strong signals with min{kβk2 , kγk2 } ≥ η0 ≥ Ck log p/n,
where C is some constant, the optimal convergence rate of estimating R (β, γ) is


k log p
1
1
1 k log p
√ +
+ 2
.
η0
n
η0 n
n
In contrast to estimating I (β, γ), the optimal rate scales to the inverse of the signal strength,
b (β, γ) proposed in Section 3 adaptively
represented by 1/η0 . The estimators bI (β, γ) and R
achieve the optimal rates for estimating I (β, γ) and R (β, γ), respectively. In addition,
the optimal rate of convergence for estimating the quadratic functionals Q (β) = kβk22 and

Q (γ) = kγk22 is also established.

1.4

Organization of the Paper

The rest of the paper is organized as follows. Section 2 presents the procedures for estimating
I (β, γ), Q (β), Q (γ), and R (β, γ) in detail. In Section 3, minimax rates of convergence for
the estimation problems are established and the proposed estimators are shown to attain
the optimal rates. In Section 4, simulation studies are conducted to evaluate the empirical
performance of the FDEs. A yeast cross data is used to illustrate the estimators in Section
5. The proofs of main theorems are present in Section 6. The remaining proofs and the
extended simulation studies are given in the supplementary materials.
6

1.5

Notation and Definitions

Basic notation and definitions that will be used in the rest of the paper are defined here. For
a matrix X ∈ Rn×p , Xi· , X·j , and Xi,j denote respectively the i-th row, j-th column, and
(i, j)-th entry of the matrix X, Xi,−j denotes the i-th row of X excluding the j-th coordinate,
and X−j denotes the sub-matrix of X excluding the j-th column. Let [p] = {1, 2, · · · , p}.
For a subset J ⊂ [p], XJ denotes the sub-matrix of X consisting of columns X·j with
j ∈ J and for a vector x ∈ Rp , xJ is the sub-vector of x with indices in J and x−J is
the sub-vector with indices in J c . For a vector x ∈ Rp , the `q norm of x is defined as
1
P
kxkq = ( qi=1 |xi |q ) q for q ≥ 0 with kxk0 denoting the cardinality of non-zero elements of

x and kxk∞ = max1≤j≤p |xj |. For a matrix A and 1 ≤ q ≤ ∞, kAkq = supkxkq =1 kAxkq

is the matrix `q operator norm. In particular, kAk2 is the spectral norm. For a symmetric
matrix A, λmin (A) and λmax (A) denote respectively the smallest and largest eigenvalue of
A. For a set S, |S| denotes the cardinality of S. For a ∈ R, a+ = max {a, 0} and sign(a)
is the sign of a, i.e., sign(a) = 1 if a > 0, sign (a) = −1 if a < 0 and sign (0) = 0. Define
1

the Sub-gaussian norm kxk2ψ2 of x ∈ Rp as kxk2ψ2 = supv∈S p−1 supq≥1 (E|v | x|q /q) q where

S p−1 is the unit sphere in Rp . The random vector x ∈ Rp is defined to be Sub-gaussian
if its corresponding Sub-gaussian norm is bounded; see Vershynin (2012) for more on Subgaussian random variables. For the design matriices X ∈ Rn1 ×p and Z ∈ Rn2 ×p , we define
b = X | X/n1 and Γ
b = Z | Z/n2 . Let
the corresponding sample covariance matrices as Σ
zα/2 denote the upper α/2 quantile of the standard normal distribution. For two positive

sequences an and bn , an . bn means an ≤ Cbn for all n and an  bn if an . bn and bn . an .
c and C are used to denote generic positive constants that may vary from place to place.

2

Estimation Methods

Details of the FDE estimation procedures are given in this Section. Estimation of I (β, γ) in
Section 2.1 is considered first and then estimation of the quadratic functionals Q (β) = kβk22

and Q (γ) = kγk22 is given in Section 2.2. An estimator of the normalized inner product
R (β, γ) is given in Section 2.3 by combining the estimators of I (β, γ), Q (β) and Q (γ)
developed in Sections 2.1 and 2.2.
7

2.1

Estimation of I (β, γ)

Estimation of the inner product I (β, γ) = hβ, γi is considered first since I (β, γ) is of
significant interest in its own right. In addition, constructing a good estimator of I (β, γ) is
also a key step in developing an optimal estimator for the normalized inner product R (β, γ).
The scaled Lasso estimators for high-dimensional linear model (1) are defined through
the following optimization algorithm (Sun and Zhang, 2012),

and

b σ̂1 } = arg
{β,
{b
γ , σ̂2 } = arg

where λ0 =

√
2.01 log p.

min

β∈Rp ,σ1 ∈R+

min

γ∈Rp ,σ2 ∈R+

p
ky − Xβk22 σ1
λ0 X kX·j k2
+
+√
|βj |,
√
2n1 σ1
2
n1 j=1
n1
p
λ0 X kZ·j k2
kw − Zγk22 σ2
+
+√
|γj |,
√
2n2 σ2
2
n2 j=1
n2

(4)

(5)

A natural estimator of the inner product I (β, γ) can then be obtained by simply plugging
b given in (4) and (5) respectively. However, as we
in the scaled Lasso estimators βb and γ
b γ
b i is not a good estimate of I (β, γ) in general due to the bias in the Lasso or
will show, hβ,
scaled Lasso estimators and other reasons. An alternative is to use the de-biased estimator,

which was originally introduced for the construction of confidence intervals (Cai and Guo,

2016a; Javanmard and Montanari, 2014; van de Geer et al., 2014; Zhang and Zhang, 2014).
Plugging in the de-biased estimator does not lead to a good estimator either; it has large
variance. See Reamrk 2 for further discussions.
To construct an optimal estimator of I (β, γ), it is helpful to decompose the difference
b γ
b i and hβ, γi into three terms,
between hβ,
b γ
b γ
b i − hβ, γi = hb
b − γi − hβb − β, γ
b − γi.
hβ,
γ , βb − βi + hβ,

(6)

b − γi is “small”, but the first two terms
The last term on the right hand side, hβb − β, γ
b γ
b − γi can be large. Our approach is to first estimate these two terms
hb
γ , βb − βi and hβ,
b γ
b i to obtain the final estimate of I (β, γ).
and then subtract them from hβ,
b is given here. Since
Some intuition for estimating hb
γ , β − βi




1 |
b β − βb + 1 X | ,
X y − X βb = Σ
n1
n1
8

(7)

multiplying both sides of (7) by a vector u ∈ Rp yields



1
1 | |
|b
b
b
u X y − X β = u Σ β − β + u| X | ,
n1
n1

which can be written as


| 

1 | |
1
b = Σu
b −γ
b
u X y − X βb − hb
γ , β − βi
β − βb + u| X | .
n1
n1

(8)

(9)



If the vector u can be chosen such that the right hand side of (9) is “small”, then u| X | y − X βb /n1
b
is a good estimate of hb
γ , β − βi.

To construct u, note that the first term on the right hand side of (9) can be further

| 

b
b −γ
b −γ
b
b k∞ kβb − βk1 . To estimate hb
β − βb ≤ kΣu
γ , β − βi,
upper bounded as Σu

| 

b −γ
b
β − βb is controlled through
a projection vector u is constructed such that Σu
b −γ
b k∞ and the second term of (9) u| X | /n1 , which has mean 0,
constraining the term kΣu
b
is controlled through minimizing its variance σ12 u| Σu/n
1 . This leads to the following convex
b
optimization algorithm for identifying the projection vector u for estimating hb
γ , β − βi,


λ1
|b
b
b 1 = arg min u Σu : kΣu − γ
b k∞ ≤ kb
u
γ k2 √
,
(10)
n1
u∈Rp
√
where λ1 = 12λ2max (Σ) log p.
b 1 is
Remark 1. The solution of the above optimization problem might not be unique and u
defined as any minimizer of the optimization problem. The theory established in Section 3
still holds for any minimizer of (10). The optimization problem (10) is solved through its

equivalent Lagrange dual problem, which is computationally efficient and scales well to the
high-dimensional problem. See Step 2 in Table 1 for more details.


b is then estimated by u
b 1 is obtained, hb
b |1 X | y − X βb /n1 .
Once the projection vector u
γ , β−βi
b γ −γ
b i can be obtained via the convex alSimilarly, the projection vector for estimating hβ,
gorithm



λ2
|b
b
b
b
b 2 = arg min u Γu : kΓu − βk∞ ≤ kβk2 √
u
,
n2
u∈Rp
√
b γ−γ
b ) /n2 .
b i is estimated by u
b |2 Z | (w − Z γ
where λ2 = 12λ2max (Γ) log p. Then hβ,
The final estimator bI (β, γ) of I (β, γ) is given by


1
1
b γ
bI (β, γ) = hβ,
bi + u
b |1 X | y − X βb + u
b |2 Z | (w − Z γ
b) .
n1
n2
9

(11)

(12)

It is clear from the above discussion that the key idea for the construction of the final
b and
b 1 and u
b 2 such that hb
estimator bI (β, γ) is to identify the projection vectors u
γ , β − βi
b γ−γ
b i are well estimated. It will be shown in Section 3 that the estimator bI (β, γ) is
hβ,
adaptively minimax rate-optimal.

Remark 2. As mentioned before, simply plugging in the Lasso, scaled Lasso, or de-biased
estimator does not lead to a good estimator of I (β, γ). Another natural approach is to first
threshold the de-biased estimator to obtain a sparse estimator of the coefficient vectors (see
details in Zhang and Zhang (2014, Section 3.3), Guo et al. (2016, equation (10))) and then
plug in. This estimator is referred as the thresholded estimator. Simulations in Section
4 demonstrate that the proposed estimator defined in (12) outperforms the three plug-in
estimators using the scaled Lasso, de-biased, and thresholded estimators.

2.2

Estimation of Q (β) and Q (γ)

In order to estimate the normalized inner product R (β, γ), it is necessary to estimate the
quadratic functionals Q (β) = kβk22 and Q (γ) = kγk22 . To this end, we randomly split the


data (y, X) into two subsamples y (1) , X (1) with sample size n1 /2 and y (2) , X (2) with

sample size n1 /2 and the data (w, Z) into two subsamples w(1) , Z (1) with sample size n2 /2

and w(2) , Z (2) with sample size n2 /2.
b denote the optimizers of the scaled Lasso
With a slight abuse of notation, let βb and γ


algorithm (4) applied to y (1) , X (1) and (5) applied to w(1) , Z (1) , respectively. For the

scaled Lasso algorithms, the sample sizes n1 and n2 are replaced by n1 /2 and n2 /2, respecb 2 of Q (β) and kb
tively. Again, the simple plug-in estimator kβk
γ k2 of Q (γ) do not perform
2

2

well. Note that

b 2 − kβk2 = 2hβ,
b βb − βi − kβb − βk2 .
kβk
2
2
2

(13)

The second term on the right hand side of (13) is “small”, but the first can be large. The
b 2 does not estimate kβk2 well for this reason. Similar ideas as in
plug-in estimator kβk
2
2
b β − βi
b is
estimating hβ, γi are applied here to estimate kβk2 . Specifically, the term 2hβ,
2

b 2 to obtain the final estimate of kβk2 . To estimate
estimated first and then is added it to kβk
2
2
10

b β − βi,
b a projection vector uis identified such that the following difference is controlled,
hβ,




|  (2)

1
b β − βi
b = u| Σ
b (2) − βb β − βb + 1 u| X (2) | ,
y − X (2) βb − hβ,
u| X (2)
n1 /2
n1 /2
(14)

b (2) = X (2) | X (2) /(n1 /2). Let the projection vector u
b 3 be the solution to the followwith Σ
ing optimization algorithm

(

)

b ∞ ≤ kβk
b 2 pλ1
b (2) u : kΣ
b (2) u − βk
b 3 = arg min u| Σ
,
(15)
u
p
u∈R
n1 /2

  (2)
√
|
(2) |
(2) b
2
b
b
b3 X
y − X β /(n1 /2)
where λ1 = 12λmax (Σ) log p. We then estimate hβ, β−βi by u

and propose the final estimator of kβk22 as


|  (2)
1
|
2
(2)
(2)
b + 2b
b (β) = kβk
Q
u3
X
y − X βb
.
2
n1 /2
+

(16)

Similarly, the estimator of kγk22 is given by




| 1
(2) |
(2)
(2)
2
b
b
Z
w −Z γ
,
Q (γ) = kb
γ k2 + 2b
u4
n2 /2
+

where

(

b (2)
with Γ

λ2
b (2) u : kΓ
b (2) u − γ
b 4 = arg min u| Γ
b k∞ ≤ kb
u
γ k2 p
u
n2 /2


√
|
= Z (2)
Z (2) /(n2 /2) and λ2 = 12λ2max (Γ) log p.

)

(17)

,

(18)

Remark 3. Sample splitting is used here for the purpose of the theoretical analysis. In the
Simulation Section, the performance of the proposed estimator without sample splitting is
investigated; see Steps 5 - 7 in Table 1. The proposed estimator without sample splitting
performs even better numerically than with sample splitting since more observations are
b 2 and kb
b 3 and
used in constructing the initial estimators kβk
γ k2 and the projection vectors u
2

b4 .
u

2.3

2

Estimation of R (β, γ)

b (β) and Q
b (γ) constructed in Sections 2.1 and 2.2, a natural
Given the estimators bI (β, γ), Q

estimator for the normalized inner product R (β, γ) is given by


 |bI (β, γ) |


n
o 
b (β) Q
b (γ) > 0 , 1 ,
b (β, γ) = sign bI (β, γ) · min q
1 Q
R

 b
b
Q (β) Q (γ)
11

(19)

b (β) and Q
b (γ) are estimators of hβ, γi, kβk22 and kγk22 defined in (12), (16)
where bI (β, γ), Q
b (β) and Q
b (γ) is 0 if kβk2 and kγk2 are
and (17), respectively. It is possible that one of Q
2
2

close to zero. In this case, the normalized inner product R (β, γ) is estimated as 0. Since
b (β, γ) is truncated to ensure that it is
R (β, γ) is always between −1 and 1, the estimator R

within the range.

b (β),
The FDE algorithm without sample splitting for calculating the estimators bI (β, γ), Q
b (γ), and R
b (β, γ) is detailed in Table 1.
Q

3

Theoretical Analysis

The theoretical properties of the estimators introduced in the last section are investigated
in this secrtion. The upper bounds and lower bounds together show that the proposed
b (β), Q
b (γ), and R
b (β, γ) achieve the optimal rates of convergence.
estimators bI (β, γ), Q

3.1

Upper Bound Analysis

The samples sizes n1 and n2 are assumed to be of the same order, that is n1  n2 . Let
n = min{n1 , n2 } be the smallest of two sample sizes. The following assumptions are needed
to facilitate the theoretical analysis.
(A1) The covariance matrices Σ and Γ satisfy 1/M1 ≤ λmin (Σ) ≤ λmax (Σ) ≤ M1
and 1/M1 ≤ λmin (Γ) ≤ λmax (Γ) ≤ M1 , and the noise levels σ1 and σ2 satisfy
max{σ1 , σ2 } ≤ M2 , where M1 ≥ 1 and M2 > 0 are positive constants.
(A2) The `2 norms of the coefficient vectors β and γ are bounded away from zero in the
sense that
min {kβk2 , kγk2 } ≥ η0 ≥ C

r

k log p
,
n

where k = max{kβk0 , kγk0 }.

(20)

Assumption (A1) places a condition on the spectrum of the covariance matrices Σ and Γ
and an upper bound on the noise levels σ1 and σ2 . Assumption (A2) requires that the signal
has to be bounded away from zero by η0 , which is only used in the upper bound analysis of
the normalized inner product R (β, γ).
12

Table 1: FDE algorithm without sample splitting for estimating the inner product, lengths
and normalized inner product.

Input: design matrices: X, Z; response vectors: y, w; tuning parameters λ0 ,λ.
b (β), Q
b (γ) and R
b (β, γ).
Output: bI (β, γ), Q

Initial Lasso estimators:
1.

b from (4) and (5) with the tuning parameter λ0 .
Scaled Lasso: Calculate βb and γ

Inner product calculation:
2.

b 1 : Calculate u
b 1 = arg minu u| X | Xu/4n1 + u| γ
b + λt kuk1 , where
Projection vector u
√
b 1 can not be solved with λ| replaced
λt = λt−1 /1.5, and λ0 = λ/ n1 . Repeat until u
by λt+1 , or t ≥ 10.

3.

b 2 : Calculate u
b 2 = arg minu u| Z | Zu/4n2 + u| βb + λ| kuk1 , where
Projection vector u
√
b 2 can not be solved with λ| replaced
λ| = λt−1 /1.5, and λ0 = λ/ n1 . Repeat until u

by λt+1 , or t ≥ 10.
b γ
b 1+u
bi + u
b |1 X | (y − X β)/n
b |2 Z | (w − Z γ
b )/n2 .
4. Correction: bI (β, γ) = hβ,

Quadratic functional calculation:
5.

b 3 : Calculate u
b 3 = arg minu u| X | Xu/4n1 + u| βb + λ| kuk1 , where
Projection vector u
√
b 3 can not be solved with λ| replaced
λ| = λt−1 /1.5, and λ0 = λ/ n1 . Repeat until u
by λt+1 , or t ≥ 10.

6.

b 4 : Calculate u
b 4 = arg minu u| Z | Zu/4n2 + u| γ
b + λ| kuk1 , where
Projection vector u
√
b 1 can not be solved with λ| replaced
λ| = λt−1 /1.5, and λ0 = λ/ n2 . Repeat until u

by λt+1 , or t ≥ 10.

b 2 + 2b
b 1) ,
b (β) = kβk
7. Correction: Q
u|3 X | (y − X β/n
+

| |
2
b
b )/n2 .
Q (γ) = kb
γ k + 2b
u4 Z (w − Z γ
+

Ratio calculation:

b (β, γ) = sign(bI (β, γ)) · min
8. R



 n
q
o 
b
b
b
b
b
|I (β, γ) |/ Q (β) Q (γ) 1 Q (β) Q (γ) > 0 , 1 .

13

b (β)
The following theorem establishes the convergence rates of the estimators bI (β, γ), Q
b (γ), proposed in (12), (16) and (17), respectively.
and Q

Theorem 1. Suppose the assumption (A1) holds and k log p/n → 0 with k = max {kβk0 , kγk0 }.

Then for any fixed constant 0 < α < 0.5, with probability at least 1 − 2α − p−c , we have


bI (β, γ) − I (β, γ) ≤ C (kβk2 + kγk2 ) z√α/2 + k log p + C k log p ,
(21)
n
n
n


zα/2 k log p
k log p
b
Q (β) − Q (β) ≤ Ckβk2 √ +
+C
,
(22)
n
n
n
and

b (γ) − Q (γ) ≤ Ckγk2
Q

where c and C are positive constants.



zα/2 k log p
√ +
n
n



+C

k log p
,
n

(23)

The upper bound of estimating hβ, γi not only depends on k, n and p, but also scales to
the signal strengths kβk2 and kγk2 . For the estimation of the quadratic functional Q (β), the
convergence rate depends on kβk2 . The similar phenomenon holds for estimation of Q (γ).
b (β, γ) proposed
The following theorem establishes the convergence rate of the estimator R
in (19).

Theorem 2. Suppose the assumptions (A1) and (A2) hold and k log p/n → 0 with k =
max {kβk0 , kγk0 }. Then for any fixed constant 0 < α < 0.5, with probability at least 1 − α −
p−c , we have

b (β, γ) − R (β, γ) ≤ C 1
R
η0



zα/2 k log p
√ +
n
n



+C

1 k log p
,
η02 n

(24)

where c and C are positive constants.

In contrast to Theorem 1, Theorem 2 requires the extra assumption (A2) on the signal
strengths kβk2 and kγk2 . The convergence rate of estimating R (β, γ) is scaled to the inverse
of signal strengths, 1/η0 . This is different from the error bound in Theorem 1, where the
estimation accuracy is scaled to the `2 norm. The lower bound results established in Theorem
3 will demonstrate the necessity of Assumption (A2) for estimation of R (β, γ).

14

3.2

Minimax Lower Bound

This section establishes the minimax lower bounds that show the optimality of the proposed
estimators in Section 2. We first introduce the parameter spaces for θ = (β, Σ, σ1 , γ, Γ, σ2 ),
which is the product of parameter spaces for θβ = (β, Σ, σ1 ) and θγ = (γ, Γ, σ2 ). To
start with, we introduce the following parameter space for both θβ = (β, Σ, σ1 ) and θγ =
(γ, Γ, σ2 ),


1
G (k) = (β, Σ, σ1 ) : kβk0 ≤ k,
≤ λmin (Σ) ≤ λmax (Σ) ≤ M1 , σ1 ≤ M2 ,
M1

(25)

where M1 ≥ 1 and M2 > 0 are positive constants. The parameter space defined in (25)
requires the sparsity kβk0 ≤ k. The other conditions 1/M1 ≤ λmin (Σ) ≤ λmax (Σ) ≤ M1
and σ1 ≤ M2 are regularity conditions. Based on the definition (25), the parameter space
for (β, Σ, σ1 , γ, Γ, σ2 ) is defined as a product of two parameter spaces,
Θ(k) = G(k) × G(k) = {θ = (β, Σ, σ1 , γ, Γ, σ2 ) : (β, Σ, σ1 ) ∈ G(k), (γ, Γ, σ2 ) ∈ G(k)} .
(26)
Define a proper subspace G (k, η0 ) of G (k) as
G (k, η0 ) = G (k) ∩ {kβk2 ≥ η0 } ,

(27)

where η0 ≥ 0. Compared with G (k), the parameter space G (k, η0 ) is restricted to signals whose `2 norm is bounded away from 0 by η0 . A proper subspace of Θ(k) for θ =
(β, Σ, σ1 , γ, Γ, σ2 ) can be defined as
Θ(k, η0 ) = {θ = (β, Σ, σ1 , γ, Γ, σ2 ) : (β, Σ, σ1 ) ∈ G(k, η0 ), (γ, Γ, σ2 ) ∈ G(k, η0 )} .

(28)

The following theorem establishes the minimax lower bounds for the convergence rates of
estimating the inner product hβ, γi, the quadratic functionals Q (β) and Q (γ) and the
normalized inner product R (β, γ).
Theorem 3. Suppose k ≤ c min {n/log p, pν } for some constants c > 0 and 0 < ν < 1/2,
then
inf sup Pθ
bI θ∈Θ(k)







1
bI − I (β, γ) ≥ c0 (kβk2 + kγk2 ) √1 + k log p + k log p
≥ ,
n
n
4
n
15

(29)

inf sup Pθβ
b θβ ∈G(k)
Q





b − Q (β) ≥ c
Q

0



kβk2







1
k log p
√ +
n
n



k log p
+
n



1
≥ ,
4

(30)



1
k
log
p
k
log
p
1
b − Q (γ) ≥ c kγk2 √ +
inf sup Pθγ Q
+
≥ ,
(31)
b θγ ∈G(k)
n
n
4
n
Q



 
1
1 k log p
1
1
k log p
0
b
√ +
+ 2
,1
≥ , (32)
inf sup Pθ R − R (β, γ) ≥ c min
b θ∈Θ(k,η0 )
η0
n
η0 n
4
n
R
0

where c0 > 0 is a positive constant.

b (β),
Combined with Theorem 1, the above theorem implies that the estimators bI (β, γ), Q

b (γ) proposed in (12),(16) and (17) achieve the minimax lower bounds (29), (30) and (31),
Q

so they are minimax rate-optimal estimators of I (β, γ), Q (β) and Q (γ), respectively. For
p
estimation of R (β, γ), under the assumption that η0 ≥ C k log p/n, Theorem 2 shows
b (β, γ) proposed in (19) achieves the minimax lower bound 1/η0 ×
that the estimator R
√
b (β, γ) is the rate-optimal estimator
(1/ n + k log p/n) + 1/η02 × k log p/n in (32). Hence R
of R (β, γ) under the assumption that the signal strengths kβk2 and kγk2 are bounded away
p
from zero. However, if η0 ≤ c k log p/n, the estimation of R (β, γ) is not interesting since

the trivial estimator 0 will achieve the constant minimax lower bound in this case, which
demonstrates the necessity of Assumption (A2).

4

Simulation Evalulations and Comparisons

In this section, the performance of several estimators for I (β, γ) and R (β, γ) were compared
in two different settings. These estimators included plug-in scaled Lasso estimator (Sun and
Zhang, 2012), plug-in de-biased estimator (Cai and Guo, 2016a; Javanmard and Montanari,
2014; van de Geer et al., 2014; Zhang and Zhang, 2014), plug-in thresholded estimator
(Zhang and Zhang, 2014, Section 3.3) and the proposed estimator FDE. Details about these
estimators are listed as follows.
• Plug-in scaled Lasso estimator (Lasso): The inner product I (β, γ) is estimated by
b γ
b i and the normalized inner product R (β, γ) is estimated by
hβ,
h

b γ
b 2 kb
b i/ kβk
hβ,
γ k2
16

o
i n
b 2 kb
1 kβk
γ k2 > 0 .

• Plug-in de-biased estimator (De-biased): Denote the de-biased Lasso estimators as βe
e γ
e . The inner product I (β, γ) is estimated by hβ,
e i and the normalized inner
and γ
h
o
i n
e 2 ke
e 2 ke
e γ
e i/ kβk
product R (β, γ) is estimated by hβ,
γ k2 1 kβk
γ k2 > 0 .

• Plug-in thresholded estimator (Thresholded): Denote the thresholded estimators as
β̄ and γ̄. The inner product I (β, γ) is estimated by hβ̄, γ̄i and the normalized inner

 
product R (β, γ) is estimated by hβ̄, γ̄i/ kβ̄k2 kγ̄k2 1 kβ̄k2 kγ̄k2 > 0 .

• FDE: The inner product I (β, γ) is estimated by bI (β, γ) in (12) and the ratio R (β, γ)
b (β, γ) in (19). We consider FDE with sample splitting (FDE-S)
is estimated by R
b (β, γ).
and without sample splitting (FDE-NS) for R

Implementation of the de-biased, thresholded and FDE estimators requires the scaled
b in the initial step. The scaled Lasso estimator was implemented
Lasso estimators βb and γ

by the equivalent square-root Lasso algorithm (Belloni et al., 2011). The theoretical tuning
p
parameter is λ0 = 2.01 log p/n, which may be conservative in the numerical studies. In√
stead, the tuning parameter was chosen as λ0 = b 2.01 log p. However, the performances of
all estimators were evaluated across a grid of tuning parameter values b ∈ {.25, .5, .75, 1} (see
Supplementary Material, Section A.1). The results showed that b = .5 was a good choice for
p
all the estimators. Hence, λ0 = .5 2.01 log p/n was used for the numerical studies in this

section and Section 5. To implement the FDE algorithm, the other tuning parameter λ as
p
2.01 log p/n for the correction Steps 2,3,5 and 6 are given in Table 1 .

Comparisons of estimates of I (β, γ) and R (β, γ) are presented below. Results on es-

timating the quadratic functionals are presented in the Supplementary Material, Section
A.2. For each setting, with the parameters (p, n1 , n2 , s, s1 , s2 ), Σ, Γ, Fβ , Fγ specified, the
following steps were implemented:
1. Generate sets S1 ⊂ [p] and S2 ⊂ [p], with |S1 | = s1 , |S2 | = s2 and |S1 ∩ S2 | = s. For
β ∈ Rp and γ ∈ Rp , generate βj ∼Fβ and γk ∼Fγ , for j ∈ S1 , k ∈ S2 , and set βj = 0

and γk = 0, for j 6∈ S1 , k 6∈ S2 .
i.i.d

i.i.d

2. Generate Xi· ∼ N (0, Σ), 1 ≤ i ≤ n1 , and Zi· ∼ N (0, Γ), 1 ≤ i ≤ n2 .
17

i.i.d

i.i.d

3. Generate the noise i ∼ N (0, 1), 1 ≤ i ≤ n1 , and δi ∼ N (0, 1), 1 ≤ i ≤ n2 .
Generate the outcome as y = Xβ +  and w = Zγ + δ.
4. With X, y, Z, and w, estimate I (β, γ) and R (β, γ) through different estimators.
5. Repeat 2-4 for rep replications.
b
For a given quantity T, denote the estimator from l-th replication as T(X,
y, Z, w; l).

The mean squared error (MSE) is used to measure the performance of this estimator,
rep

MSE(T) =

1 X b
(T(X, y, Z, w; l) − T)2 .
rep l=1

(33)

Experiment 1. In this experiment, the parameters were set as follows: (p, n1 , n2 , rep) =
(600, 400, 400, 300), the sparsity parameters (s, s1 , s2 ) = (15, 30, 25), and the covariance matrices Σ and Γ satisfy Σij = Γij = (0.8)|i−j| . For given positive values τ1 and τ2 , the signals
of β satisfy that βji = (1 + i/s1 ) τ1 /2, for ji ∈ S1 , i = 1, 2, · · · , s1 , and the signals in γ satisfy
that γj = τ2 for j ∈ S2 . This simulation setting aimed to investigate the case where the coefficients for one regression are much larger than the other by varying the signal strength parameters as (τ1 , τ2 ) ∈ {(3.0, .1), (2.6, .2), (2.2, .3), (1.8, .4), (.1, 1.6), (.2, 1.4), (.3, 1.2), (.4, 1.0)}.
The results are summarized in Table 2. For all combinations of the signal strength
parameters, in terms of estimating the inner product I (β, γ), FDE consistently outperformed
the plugin estimates with Lasso and Thresholded Lasso. Moreover, with increasing difference
between τ1 and τ2 , the advantage of FDE over the plugin estimate using Lasso or thresholded
Lasso became larger. The same results were observed for estimation of the normalized inner
product R (β, γ), where FDE-NS had consistent better performance than other methods.
Although De-biased performed well in terms of estimating I (β, γ), it performed much worse
than FDE-NS for estimating R (β, γ).
As discussed in Section 2, the sample splitting of estimating the normalized inner product
is simply proposed to facilitate the theoretical analysis and might not be necessary for the
algorithm. Our simulation results indicated that the proposed estimator without sample
splitting (FDE-NS) performed quite well in all settings, even better than FDE-S, due to the
fact that more samples were used for estimation and correction steps. Such observations
18

led us to use the proposed estimator without sample splitting (FDE-NS) in the real data
analysis in Section 5.
Table 2: Estimation of the inner product I (β, γ) and the normalized inner product R (β, γ):
the truth and mean squared error (MSE) for the plug-in estimator with the scaled Lasso estimator (Lasso), plug-in estimator with the de-biased estimator (De-biased), plug-in estimator
with the thresholded estimator (Thresholded), and the proposed estimator bI (β, γ) (FDE),

proposed estimators with sample splitting (FDE-S) and without sample splitting (FDE-NS).
Strength parameters, (τ1 , τ2 )
Method

I (β, γ)

R (β, γ)

(1.8, .4)

(2.2, .3)

(2.6, .2)

(3, .1) (.1, 1.6)

(.2, 1.4)

(.3, 1.2)

(.4, 1)

Truth

8.088

7.414

5.841

3.370

1.797

3.145

4.044

4.493

Lasso

9.295

11.564

12.560

7.279

2.377

4.889

5.409

4.800

De-biased

1.733

2.191

2.324

1.386

.449

.838

.985

.886

Thresholded

2.029

3.377

6.463

5.789

1.877

3.024

2.432

1.546

FDE

1.847

2.471

2.662

2.118

.734

.995

1.028

.986

Truth

.5314

.5314

.5314

.5314

.5314

.5314

.5314

.5314

Lasso

.0023

.0075

.0332

.1260

.1574

.0624

.0227

.0087

De-biased

.0208

.0415

.0864

.1590

.1736

.1068

.0627

.0373

Thresholded

.0045

.0139

.0585

.1753

.0964

.0981

.0389

.0153

FDE-S

.0337

.0303

.0621

.0678

.2130

.1199

.0694

.0616

FDE-NS

.0036

.0064

.0163

.0580

.0892

.0237

.0116

.0061

Experiment 2. In this experiment, the parameters were set as follows: (p, n1 , n2 , rep) =
(800, 400, 400, 300), signal strength parameters (τ1 , τ2 ) = (.2, .1) and the covariance matrices
Σ and Γ satisfy Σij = Γij = (0.8)|i−j| . The signals of β follow that βji = (1 + i/s1 ) τ1 /2, for
ji ∈ S1 , i = 1, 2, · · · , s1 , and the signals in γ satisfy that γj = τ2 for j ∈ S2 .
This simulation setting aimed to investigate the relationship between the performance
of estimators and the signal sparsity level and vary the number of signals in β and γ as
(s1 , s2 ) ∈ {(40, 40), (50, 50), (60, 60), (70, 70), (80, 80), (90, 90), (100, 100), (110, 110)}, and fix
the number of common signals at s = 20. Since the number of the associated variants is
very large for both coefficient vectors, large values of τ1 and τ2 induce strong signals so that
19

all the methods perform well. Instead, the signal magnitude was fixed at τ1 = .2 and τ2 = .1
as weak signals.
The results are summarized in Table 3. Clearly, FDE outperformed the other methods.
When the signals became denser, the improvement of FDE over other methods was more
pronounced. For estimation of R (β, γ), the results showed that FDE-NS consistently outperformed other estimators. As the number of signals increased, the MSE corresponding to
FDE-NS decreased quickly.
Table 3: Estimation of the inner product I (β, γ) and the normalized inner product R (β, γ):
the truth and MSE for the plug-in estimator with the scaled Lasso estimator (Lasso), plug-in
estimator with the de-biased estimator (De-biased), plug-in estimator with the thresholded
estimator (Thresholded), and the proposed estimator bI (β, γ) (FDE), proposed estimators
with sample splitting (FDE-S) and without sample splitting (FDE-NS).

Sparsity parameter, s1 (s2 = s1 )

I (β, γ)

Method

40

50

60

70

80

90

100

110

Truth

.190

.170

.219

.212

.179

.221

.183

.221

Lasso

.032

.025

.039

.036

.024

.035

.023

.028

De-biased

.015

.015

.018

.017

.024

.027

.040

.066

Thresholded .027

.021

.031

.029

.020

.025

.018

.018

FDE

.020

.014

.021

.022

.011

.013

.008

.008

Truth

.4027

.2908

.3122

.2592

.1914

.2110

.1573

.1725

Lasso

.1157 .0517

.0539

.0370

.0166

.0180

.0097

.0079

.1267 .0601

.0659

.0411

.0160

.0173

.0063

.0059

.0732 .0504

.0262

.0277

.0155

.0142

R (β, γ) De-biased

Thresholded .1392

.0687

FDE-S

.1154 .1225

.0779

.0574

.0456

.0499

.0450

.0493

FDE-NS

.0847 .0340

.0368

.0294

.0115

.0091

.0055

.0047

20

5

Co-heritability Estimation of Yeast Cross Genome
Wide Association Data

Bloom et al. (2013) reported a large scale genome-wide association study of 46 quantitative
traits based on 1,008 Saccharomyces cerevisiae segregants between a laboratory strain and
a wine strain. The data set included 11,623 unique genotype markers. Since many of these
markers were highly correlated, Bloom et al. (2013) further selected a set of 4,410 markers
that are weakly dependent based in the linkage disequilibrium information. These markers
were selected by picking one marker closest to each centimorgan position on the genetic
map. The maker genotype was coded as 1 or −1, according to which strain it came from.
For each segregant, the traits of interest were the end-point colony size normalized by the
control growth under 46 different growth media, including Hydrogen Peroxide, Diamide,
Calcium, Yeast Nitrogen Base (YNB) and Yeast extract Peptone Dextrose (YPD), etc. To
demonstrate the co-heritability among these traits, 8 traits were considered, including the
normalized colony size under Calcium Chloride (Calcium), Diamide, Hydrogen Peroxide
(Hydrogen), Paraquat, Raffinose, 6 Azauracil (Azauracil), YNB, and YPD. Each trait was
normalized to have variance 1, so the quadratic norm represents the total genetic effects for
each trait.
FDE was applied to every pair of these 8 traits without sample splitting, and the results
are summarized in Table 4, including estimates of the inner product for all possible pairs,
the quadratic norm for each trait and the normalized inner product R (β, γ). The genetic
heritability of these traits ranged from 0.22 for Raffinose to 0.67 for YPD. About two thirds
of these pairs had an estimated co-heritability smaller than 0.1, indicating relatively weak
genetic correlations among these traits.
To further demonstrate the genetic relatedness among these pairs, for each trait, a tstatistic was calculated based on regressing the trait value y on genetic genetic marker X·j ,
for i ≤ j ≤ p. A larger absolute value of the t-statistic implied a stronger effect of the marker
on the trait. For any pair of traits, the scatter plot of the t-statistics provided a way of revealing the genetic relationship between them. Figure 1 shows the plots of several pairs of
the traits, including the pairs with a large positive R (β, γ), YPD v.s. YNB and Calcium v.s.
21

Table 4: FDE estimation for the inner product, the quadratic norm (top row of each block)
and the normalized inner product (bottom row of each block) among 8 traits of yeast cross.
Traits

Calcium

Diamide Hydrogen

Paraquat

Raffinose

Azauracil

YNB

Calcium

.3314

-.0189

-.1003

.0084

.0927

.0095

-.0656 .-0134

1

-.0286

-.1579

.0117

.1972

.0172

-.0968 -.0164

.4390

.0598

-.0039

.0500

.0446

-.0159 .0803

1

.0942

-.0053

.1065

.0809

-.0235 .0983

.4033

.0576

-.1040

.0601

.0672

.0637

1

.0799

-.2213

.1089

.0991

.0779

.5199

.0023

.0365

.1148

.1029

1

.0049

.0661

.1693

.1259

.2208

.0137

.0830

.0331

1

.0248

.1224

.0405

.3045

-.0259

.0703

1

-.0383

.0860

.4594

.4246

1

.5195

Diamide

Hydrogen

Paraquat

Raffinose

Azauracil

YNB

YPD

YPD

.6680
1

YNB, pairs with a large negative R (β, γ), Raffinose v.s. Hydrogen and Calcium v.s. Hydrogen, and pairs with R (β, γ) near to 0, including Paraquat v.s. Diamide and Paraquat v.s.
Raffinose. The plot clearly indicates a strong positive genetic correlation between YPD and
YNB. The correlation between Calcium and YNB is slightly smaller. Raffinose/Hydrogen
and Calcium/Hydrogen pairs clearly showed negative genetic correlation. There were several genetic variants with very large effects on Hydrogen, but they were not associated with
the other traits such as Raffinose and Calcium. The shared genetic variants were relatively
weak, leading to a correlation around -0.2. The two plots on the right show the pairs of trats
with weak correlations. These plots indicated that the proposed co-heritability measures can
22

15

10

10
10

5

0

Paraquat

Raffinose

YPD

5
5
0

−5

0

−5

−5
−10
−5

0
YNB

5

−10
−10

10

−5

0
5
Hydrogen

−10
−10

10

10

10

5

5

5

0

Paraquat

10

Calcium

Calcium

−10

0

−5

0
Diamide

5

10

0

−5

−5
−5
−10

−5

0
YNB

5

10

−10
−10

−10
−5

0
5
Hydrogen

10

−10

−5

0
5
Raffinose

10

15

Figure 1: Scatter plots of marginal regression t statistics for six pairs of traits, including the
pairs with positive correlation case (left panel), negative correlation (middle panel), and no
correlation (right panel).
indeed capture the genetic sharing among different related traits.

6
6.1

Proofs
Proof of Theorem 1

For simplicity of notation, we assume n1 = n2 and use n = n1 = n2 to represent the sample
size throughout the proof. The proofs can be easily generalized to the case n1  n2 . Without
loss of generality, we assume that the Sub-gaussian norm of random vectors Xi· and Zi· are

also upper bounded by M1 , that is, max kXi· k2ψ2 , kZi· k2ψ2 ≤ M1 .
Proof of (21)

23

Note the following important decomposition,


1
1
b γ
bI (β, γ) − I (β, γ) = hβ,
bi + u
b |1 X | y − X βb + u
b |2 Z | (w − Z γ
b ) − hβ, γi
n  
n




1
|1
|
|
|
b + u
b γ−γ
b ) − hβ,
b i − hβb − β, γ
b − γi
b 1 X y − X βb − hb
b 2 Z (w − Z γ
= u
γ , β − βi
n
n





1
1
b −γ
b − βb| (γ − γ
b |2 Z | δ + u
b ) − hβb − β, γ
b − γi.
b |1 Σ
b | β − βb + u
b |2 Γ
=b
u|1 X |  + u
n
n
(34)

The following lemmas are introduced to control the terms in (34) and similar results were

established in the analysis of Lasso, scaled Lasso and de-biased Lasso (Cai and Guo, 2016a;
Ren et al., 2013; Sun and Zhang, 2012; Ye and Zhang, 2010). The proofs of the following
lemmas can be found in the supplementary material, Section C.
Lemma 1. Suppose the assumption (A1) holds and k log p/n → 0 with k = max {kβk0 , kγk0 }.

Then with probability at least 1 − p−c , we have
r
r
log
p
log p
, kb
γ − γk1 ≤ Ck
,
kβb − βk1 ≤ Ck
n
n
and

kβb − βk2 ≤ C

r

k log p
,
n

kb
γ − γk2 ≤ C

r

k log p
,
n

(35)

(36)

where c and C are positive constants.
Lemma 2. Suppose the assumption (A1) holds and k log p/n → 0 with k = max {kβk0 , kγk0 }.

Then with probability at least 1 − p−c , we have





kβk0 log p
b 2 kγk0 log p .
b −γ
b − βb| (γ − γ
b |1 Σ
b | β − βb ≤ Ckb
b |2 Γ
b ) ≤ Ckβk
u
γ k2
and u
n
n
(37)
With probability at least 1 − 2α − p−c ,

and hence

zα/2
1
b |1 X |  ≤ Ckb
γ k2 √ ,
u
n
n

and


zα/2
k log p
b
b ≤ Ckb
y − X β − hb
γ , β − βi
γ k2 √ + Ckb
γ k2
,
n
n
n
1
b γ−γ
b 2 z√α/2 + Ckβk
b 2 k log p ,
b ) − hβ,
b i ≤ Ckβk
b |2 Z | (w − Z γ
u
n
n
n

1
b |1 X |
u



1
b 2 z√α/2 ,
b |2 Z | δ ≤ Ckβk
u
n
n

where c and C are positive constants.

24

(38)

(39)

By the decomposition (34) and the inequalities (36), (37) and (38), we obtain that


 z
k log p
k
log
p
α/2
b 2 + kb
bI (β, γ) − I (β, γ) ≤ C kβk
√ +
+C
.
γ k2
n
n
n

By (36), we establish (21).
Proof of (22) and (23)

The proof of (23) is similar to that of (22) and only the proof of (22) is present in the

  (2)
|
(2) b
(2) |
2
b
y − X β /(n/2)
u X
following. We introduce the estimator Q̄(β) = kβk + 2b
2

3

b (β) − Q (β) ≤ Q̄(β) − Q (β) .
and due to the fact that Q (β) is non-negative, we have Q

We decompose the difference between Q̄(β) and Q (β),


  (2)
| 1
2
2
(2) |
(2) b
b
Q̄(β) − Q (β) = kβk2 − kβk2 + 2b
u3
X
y −X β
n/2



|
1
b (2) − βb| β − βb + 2b
b |3 Σ
=2 u
u|3
X (2) (2) − kβb − βk22 .
n/2

Combined with the above argument, the upper bound (22) follows from (36) and the following
lemma. The proof of the following lemma can be found in the supplementary material Section
C.
Lemma 3. Suppose the assumption (A1) holds and k log p/n → 0 with k = max {kβk0 , kγk0 }.

Then with probability at least 1 − p−c − α,

and



b |3
u

|
1
b 2 z√α/2 ,
X (2) (2) ≤ Ckβk
n/2
n

b (2)
b |3 Σ
u

− βb|

where c and C are positive constants.

6.2



b
b 2 k log p ,
β − β ≤ Ckβk
n

(40)

(41)

Proof of Theorem 2



b ) /n. We introduce the new estimator
b |1 X | y − X βb /n and µ
b |2 Z | (w − Z γ
Let µ
e1 = u
e2 = u
b γ
bI (β, γ)
bi + µ
hβ,
e +µ
e2
q1
R̄(β, γ) = q
= q
.
b (β) Q
b (β)
b (β) Q
b (γ)
Q
Q
25

b (β, γ) − R (β, γ) ≤ R̄(β, γ) − R (β, γ) . We decomSince −1 ≤ R (β, γ) ≤ 1, we have R

pose the difference between R̄(β, γ) and R (β, γ),
*
+ 

b γ
bi + µ
b
hβ,
e1 + µ
e2
βb
β
γ
µ
e +µ
e
hβ, γi
γ
q
q
q
q
q 1 q2
=
−
,
+
−p
,
kβk2 kγk2
kβk22 kγk22
b (β) Q
b (γ)
b (β)
b (γ)
b (β) Q
b (γ)
Q
Q
Q
Q


*
+
*
+
b
b
βb
γ
β
γ
µ
e2
β
γ
γ

q
= q
,q
+
,q
+q
−
−
−
kβk2
kγk2
kβk2
kγk2
b
b
b
b
b
Q (β)
Q (γ)
Q (γ)
Q (β) Q (γ)

*
+
b
β
γ
µ
e1
β
.
q
−
,
+q
+ q
kβk
kγk
2
2
b (β)
b (β) Q
b (γ)
Q
Q
(42)
The following lemma controls the terms in the above equation and establishes the theorem.
The proof of the following lemma can be found in the supplementary material Section C.
Lemma 4. Suppose the assumptions (A1) and (A2) hold and k log p/n → 0 with k =
max {kβk0 , kγk0 }. Then with probability at least 1 − p−c − α, we have
*

*

βb

q
b (β)
Q

and
*

βb

b
β
γ
γ
q
−
−
,q
kβk
kγk
2
2
b (β)
b (γ)
Q
Q

γ
β
,
−
kβk2 kγk2

b
β
γ
γ
,q
−
kβk2
kγk
2
b (γ)
Q

+

+

+

µ
e1
1
q
+q
≤C
kβk2
b (β) Q
b (γ)
Q

µ
e3
1
q
+q
≤C
kγk
2
b (β) Q
b (γ)
Q

≤C


k log p
1
,
kβk2 kγk2 n

zα/2 k log p
√ +
n
n



+C

(43)


1
1
+
2
kβk2 kβk2 kγk2

k log p
,
n



k log p
,
n

(44)


zα/2 k log p
√ +
n
n



+C



1
1
+
kγk22 kβk2 kγk2

where c and C are positive constants.

By the decomposition (42), combining (43), (44) and (45), we have
bI (β, γ)
hβ, γi
q
−p
kβk22 kγk22
b (β) Q
b (β)
Q





zα/2 k log p
1
1
1
1
k log p
1
√ +
≤C
+
+C
+
+
.
2
2
kβk2 kγk2
n
kβk2 kγk2 kβk2 kγk2
n
n
Under the assumption (A2), the upper bound (24) follows from the above inequality.
26



(45)

6.3

Proof of (30) and (31) in Theorem 3

We first introduce the notations used in the proof of lower bound results. Let π denote
the prior distribution supported on the parameter space H. Let fπ (z) denote the density
function of the marginal distribution of the random variable Z with the prior π on H.
R
More specifically, fπ (z) = fθ (z) π (θ) dθ. We define the χ2 distance between two density
functions f1 and f0 by

Z

Z 2
(f1 (z) − f0 (z))2
f1 (z)
χ (f1 , f0 ) =
dz =
dz − 1
f0 (z)
f0 (z)
R
and the L1 distance by L1 (f1 , f0 ) = |f1 (z) − f0 (z)| dz. It is well known that
p
L1 (f1 , f0 ) ≤ χ2 (f1 , f0 ).
2

(46)

(47)

The proof of the lower bound is based on the following version of Le Cam’s Lemma

(LeCam (1973); Ren et al. (2013); Yu (1997)).
Lemma 5. Let T (θ) denote a functional on θ. Suppose that H0 = {θ0 }, H0 , H1 ⊂ Θ and
d = minθ∈H1 |T (θ) − T (θ0 )|. Let π denote a prior on the parameter space H1 . Then we
have
inf

sup

b θ∈H0 ∪H1
T

Pθ



b − T (θ) ≥ d
T
2



≥

1 − L1 (fπ , fθ0 )
.
2

(48)

The proofs of (30) and (31) are applications of Lemma 5. The key is to construct the
parameter spaces H0 = {θ0 }, H1 and the prior on H1 such that (i) H0 , H1 ⊂ Θ,(ii)the
L1 distance L1 (fπ , fθ0 ) is controlled and (iii)the distance d = minθ∈H1 |T (θ) − T (θ0 )| is
maximized. In the following, we provide the detailed proof of (30). The proof of (31) is very
similar to that of (30) and is omitted here. In the discussion of lower bound results, we will
assume µx = µz = 0 and the design Xi· and Zi· follow joint normal distribution.
We define the following two parameter spaces to facilitate the discussion,
(
)
r


k
log
p
M
2
and G W (k) = G (k) ∩ kβk2 ≤ C
.
G S (k) = G (k) ∩ kβk2 ≥
4
n
The lower bound (30) can be decomposed into the following three lower bounds,


k log p
1
0
b
inf sup Pθβ Q − Q (β) ≥ c kβk2
≥ ,
b θ ∈G S (k)
n
4
Q
β
27

(49)

(50)

inf

sup

b θ ∈G S (k)
Q
β

and

inf

sup

P θβ

b θ ∈G W (k)
Q
β



P θβ

b − Q (β) ≥ c kβk2 √1
Q
n



0

b − Q (β) ≥ c
Q

0k

log p
n





1
≥ ,
4

1
≥ ,
4

where c0 > 0 is a positive constant. Combining (50) and (51), we have



1
k log p
1
0
b
inf sup Pθβ Q − Q (β) ≥ c kβk2 √ +
≥ .
b θ ∈G S (k)
n
4
n
Q
β
On G S (k), we have
kβk2



1
k log p
√ +
n
n



(51)

(52)

(53)





1
k log p
k log p
= max kβk2 √ +
,
n
n
n

and on G W (k), we have





k log p
k log p
1
k log p
= max kβk2 √ +
,
.
n
n
n
n

Since both G S (k) and G W (k) are proper subspace of G (k), we can establish (30) by combining (53) and (52). In the following, we will establish the lower bounds (50), (51) and (52)
separately.
Proof of (50) Under the Gaussian random design model, Vi = (yi , Xi· ) ∈ Rp+1 fol-

lows a joint Gaussian distribution with mean 0. Let Σv denote the covariance matrix
For the indices of Σv , we use 0 as the index of 
yi and {1, · · · 
, p} as the in
v
v |
Σyy Σxy
 , where Σvyy ,
dices for (Xi1 , · · · , Xip ) ∈ Rp . Decompose Σv into blocks 
v
v
Σxy
Σxx
v
v
Σxx and Σxy denote the variance of yi , the variance of Xi· and the covariance of yi and

of Vi .

Xi· , respectively. Let Ω = Σ−1 denote the precision matrix. There exists a bijective
function h : Σv →
mapping h−1 : (β, Ω, σ1 ) → Σv , where
 (β, Ω, σ1 ) and the inverse

β | Ω−1 β + σ12 β | Ω−1
 and
h−1 ((β, Ω, σ1 )) = 
Ω−1 β
Ω−1
h(Σv ) = (Σvxx )−1 Σvxy , (Σvxx )−1 , Σvyy − Σvxy

|


(Σvxx )−1 Σvxy .

(54)

Based on the bijection, it is sufficient to control the χ2 distance between two multivariate
Gaussian distributions.
28

We introduce the null parameter space
(

G0 =

M2
θ0 = (β ∗ , I, σ0 ) : β ∗ = (β1∗ , η0 , 0, · · · , 0) with β1∗ >
and σ0 =
2

r

M2
8

)

⊂ G S (k) ,

where η0 ≥ 0. Define p1 = p − 2. Based on the mapping h, we have the corresponding null
parameter space for Σv , F0 = {Σv0 } where

(β ∗ )2 + η02 + σ02 β1∗ η0
0
 1


β1∗
1 0
0
v
Σ0 = 


η0
0 1
0

0
0 0 Ip1 ×p1






.




We then introduce the alternative parameter space for Σv , which will induce a parameter
space for (β, Ω, σ1 ) through the mapping h. Define F1 = {Σvα : α ∈ ` (p1 , k, ρ)}, where


(β1∗ )2 + η02 + σ02 β1∗ η0 ρ0 α|




∗
|

β1
1 0
α 
v

,
Σα = 
(55)


η0
0 1
0 


ρ0 α
α 0 Ip1 ×p1

with ρ0 = β1∗ + σ0 and

` (p1 , k, ρ) = {α : α ∈ Rp1 , kαk0 = k − 2, αi ∈ {0, ρ} for 1 ≤ i ≤ p1 } .

(56)

Then we construct the corresponding parameter space G1 for (β, Ω, σ1 ), which is induced by
the mapping h and the parameter space F1 ,
G1 = {(β, Ω, σ1 ) : (β, Ω, σ1 ) = h (Σv )

for Σv ∈ F1 } .

(57)

For Σvα , the corresponding (β, Ω, σ1 ) defined as h (Σvα ) satisfies
β1 =

−kαk22 ρ0 + β1∗
σ0
, β2 = η0 , β−{1,2} = (ρ0 − β1 ) α =
α,
2
1 − kαk2
1 − kαk22
σ12 = σ02 −

kαk22 (ρ0 − β1∗ )2
≤ σ02 ≤ M2 ,
1 − kαk22

max {|λmin (Ω) − 1| , |λmax (Ω) − 1|} ≤ kαk2 .
29

(58)

(59)
(60)

When kαk2 is chosen such that kαk2 ≤ min {1 − 1/M1 , M1 − 1} , then we have
1
≤ λmin (Ω) ≤ λmax (Ω) ≤ M1 .
M1

(61)

The difference between β1 and β1∗ is
β1 − β1∗ =
By taking ρ =

p

kαk22 (β1∗ − ρ0 )
−σ0 kαk22
=
.
1 − kαk22
1 − kαk22

log(4p1 /k 2 )/2n, we have |β1 − β1∗ | ≤ Ck log p/n and hence β1 ≥ β1∗ −

Ck log p/n ≥ M2 /4. Combined with (59), (58) and (61), we show that G1 ⊂ G S (k). Let
π denote the uniform prior over the parameter space G1 induced by the uniform prior of α
p
over ` (p1 , k, ρ) where ρ = log(4p1 /k 2 )/2n. The control of L1 (fπ , fθ0 ) is established in the
following lemma, which follows from Lemma 2 of Cai and Guo (2016a) and is established in
(7.21) of Cai and Guo (2016a).
Lemma 6. Suppose that k ≤ c {n/log p, pγ }, where 0 ≤ γ < 12 and c is a sufficient small posp
itive constant. For ρ = log(4p1 /k 2 )/2n, we establish that kαk2 ≤ min {1 − 1/M1 , M1 − 1}
and

1
L1 (fπ , fθ0 ) ≤ .
4

(62)

To apply Lemma 5, we consider the functional T(θ) = kβk22 and calculate the distance
d = β12 + (

σ0
σ0
σ0
)2 kαk22 − (β1∗ )2 =
kαk22 β1 + β1∗ −
.
2
2
1 − kαk2
1 − kαk2
1 − kαk22

(63)

By taking β1∗ = M2 /2 > 2σ0 and η0 = 0, we have β1 ≥ M2 /2 − Ck log p/n and

σ
σ0
log p
log p
kαk22 β1 + β1∗ −
≥ ck
σ0 max {|β1 | , |β1∗ |} ≥ ck
σ0 max {kβk2 , kβ ∗ k2 } ,
2
2
1 − kαk2
1 − kαk2
n
n

∗
where the last inequality follows from the fact that kβ−1 k2 ≤ Ck log pσ0 /n and β−1
= 0.

Combined with (62), an application of Lemma 5 leads to (50).
Proof of (51) We construct the following parameter spaces,


M2
∗
∗
∗
∗
G0 = θ0 = (β , I, σ0 ) : β = (β1 , η0 , 0 · · · , 0) , and β1 ≥
⊂ G S (k)
2



¯
∗
G1 = θ1 = (β, I, σ0 ) : β = β1 + σ0 √ , η0 , 0 · · · , 0
⊂ G S (k) ,
n

(64)

where ¯ is a small positive constant to be chosen later.
The proof of the following lemma can be found in the supplementary material Section C.
30

Lemma 7. If ¯ =

p

log (17/16) /2, then we have
1
L1 (fπ , fθ0 ) ≤ .
4

(65)

To apply Lemma 5, we take η0 = 0 and calculate the distance
¯
1
¯2
d = β12 − (β1∗ )2 = 2β1∗ σ √ + σ 2
≥ c max{kβk2 , kβ ∗ k2 } √ .
n
n
n
Applying Lemma 5, we establish (51).
Proof of (52) We introduce the following null and alternative parameter spaces,
G0 = {(β, I, σ1 ) : β = (η0 , 0, 0, · · · , 0)}

(66)

G1 = {(β, I, σ1 ) : β = (η0 , 0, α) with α ∈ `1 (p, k, ρ)} ,
where


`1 (p, k, ρ) = α : α ∈ Rp−2 , kαk0 = k, αi ∈ {0, ρ} .

(67)

Let π denote the prior over the parameter space G1 induced by the uniform prior of α over
p
log(4p1 /k 2 )/2n. The control of L1 (fπ , fθ0 ) is established in the
`1 (p, k, ρ) where ρ =
following lemma, which follows from Lemma 7 of Cai and Guo (2016b) and is established in
(1.6) of Cai and Guo (2016b).
Lemma 8. Suppose that k ≤ c {n/log p, pγ }, where 0 ≤ γ <
p
positive constant. For ρ = log(4p1 /k 2 )/2n, we establish

1
2

and c is a sufficient small

1
L1 (fπ , fθ0 ) ≤ .
8

(68)

By specifying η0 = 0, we have G0 and G1 defined in (66) are proper subspaces of the

parameter space G W (k). To apply Lemma 5, we calculate the distance d = kαk22 ≥ ck log p/n.
Applying Lemma 5, we establish (52).

SUPPLEMENTARY MATERIAL
Title: Supplement to “Optimal Estimation of Co-heritability in High-dimensional Linear
Regressions”. (.pdf file)
31

References
Baradat, P. (1976), “Use of juvenile-mature relationships in individual selection including information from relatives,” Proceeding of IUFRO Meeting on advanced generation breeding,
Bordeax, pp. 121–138.
Belloni, A., Chernozhukov, V., and Wang, L. (2011), “Square-root lasso: pivotal recovery of
sparse signals via conic programming,” Biometrika, 98(4), 791–806.
Bloom, J. S., Ehrenreich, I. M., Loo, W. T., Lite, T.-L. V., and Kruglyak, L. (2013), “Finding
the sources of missing heritability in a yeast cross,” Nature, 494(7436), 234–237.
Cai, T. T., and Guo, Z. (2016a), “Confidence intervals for high-dimensional linear regression:
minimax rates and adaptivity,” The Annals of Statistics, To appear.
Cai, T. T., and Guo, Z. (2016b), “Supplement to “Confidence intervals for high-dimensional
linear regression: minimax rates and adaptivity”,” The Annals of Statistics, To appear.
Fan, J., Han, X., and Gu, W. (2012), “Estimating false discovery proportion under arbitrary
covariance dependence,” Journal of the American Statistical Association, 107(499), 1019–
1035.
Golan, D., and Rosset, S. (2011), “Accurate estimation of heritability in genome wide studies
using random effects models,” Bioinformatics, 27, i317–i323.
Guo, Z., Kang, H., Cai, T. T., and Small, D. S. (2016), “Confidence Interval for Causal
Effects with Possibly Invalid Instruments Even After Controlling for Many Confounders,”
arXiv preprint arXiv:1603.05224, .
Javanmard, A., and Montanari, A. (2014), “Confidence intervals and hypothesis testing for
high-dimensional regression,” The Journal of Machine Learning Research, 15(1), 2869–
2909.
LeCam, L. (1973), “Convergence of estimates under dimensionality restrictions,” The Annals
of Statistics, 1(1), 38–53.
32

Lee, S. H., Ripke, S., Neale, B. M., Faraone, S. V., Purcell, S. M., Perlis, R. H., Mowry,
B. J., Thapar, A., Goddard, M. E., and Witte, J. S. (2013), “Genetic relationship between
five psychiatric disorders estimated from genome-wide SNPs,” Nature Genetics, 45.
Lee, S. H., Yang, J., Goddard, M. E., Visscher, P. M., and Wray, N. R. (2012), “Estimation of pleiotropy between complex diseases using single-nucleotide polymorphism-derived
genomic relationships and restricted maximum likelihood,” Bioinformatics, 28(19), 2540–
2542.
Manolio, T. (2010), “Genomewide association studies and assessment of the risk of disease,”
New England Journal of Medicine, 363, 166–176.
McCarthy, M. I., Abecasis, G. R., Cardon, L. R., Goldstein, D. B., Little, J., Ioannidis,
J. P., and Hirschhorn, J. N. (2008), “Genome-wide association studies for complex traits:
consensus, uncertainty and challenges,” Nature Reviews Genetics, 9(5), 356–369.
Purcell, S. M., Wray, N. R., Stone, J. L., Visscher, P. M., O’Donovan, M. C., Sullivan,
P. F., Sklar, P., Ruderfer, D. M., McQuillin, A., Morris, D. W. et al. (2009), “Common
polygenic variation contributes to risk of schizophrenia and bipolar disorder,” Nature,
460(7256), 748–752.
Ren, Z., Sun, T., Zhang, C.-H., and Zhou, H. H. (2013), “Asymptotic normality and optimalities in estimation of large gaussian graphical model,” arXiv preprint arXiv:1309.6024,
.
Sun, T., and Zhang, C.-H. (2012), “Scaled sparse linear regression,” Biometrika, 101(2), 269–
284.
Tibshirani, R. (1996), “Regression shrinkage and selection via the lasso,” Journal of the
Royal Statistical Society. Series B (Methodological), 58(1), 267–288.
van de Geer, S., Bühlmann, P., Ritov, Y., and Dezeure, R. (2014), “On asymptotically optimal confidence regions and tests for high-dimensional models,” The Annals of Statistics,
42(3), 1166–1202.
33

Vershynin, R. (2012), “Introduction to the non-asymptotic analysis of random matrices,” in
Compressed Sensing: Theory and Applications, eds. Y. Eldar, and G. Kutyniok Cambridge
University Press, pp. 210–268.
Wray, N. R., Goddard, M. E., and Visscher, P. M. (2007), “Prediction of individual genetic
risk to disease from genome-wide association studies,” Genome Research, 17(10), 1520–
1528.
Yang, J., Benyamin, B., McEvoy, B. P., Gordon, S., Henders, A. K., Nyholt, D. R., Madden,
P. A., Heath, A. C., Martin, N. G., and Montgomery, G. W. (2010), “Common SNPs
explain a large proportion of the heritability for human height,” Nature Genetics, 42, 565–
569.
Yang, L., Neale, B. M., Liu, L., Lee, S. H., Wray, N. R., Ji, N., Li, H., Qian, Q., Wang, D.,
Li, J. et al. (2013), “Polygenic transmission and complex neuro developmental network
for attention deficit hyperactivity disorder: Genome-wide association study of both common and rare variants,” American Journal of Medical Genetics Part B: Neuropsychiatric
Genetics, 162(5), 419–430.
Ye, F., and Zhang, C.-H. (2010), “Rate Minimaxity of the Lasso and Dantzig Selector for
the `q Loss in `r Balls,” The Journal of Machine Learning Research, 11, 3519–3540.
Yu, B. (1997), “Assouad, Fano, and Le Cam,” in Festschrift for Lucien Le Cam Springer,
pp. 423–435.
Zhang, C.-H., and Zhang, S. S. (2014), “Confidence intervals for low dimensional parameters
in high dimensional linear models,” Journal of the Royal Statistical Society: Series B
(Statistical Methodology), 76(1), 217–242.
Zhernakova, A., van Diemen, C., and Wijmenga, C. (2009), “Detecting shared pathogenesis
from the shared genetics of immune-related diseases,” Nature Reviews Genetics, 10, 43–45.

34

Supplement to “Optimal Estimation of Co-heritability
in High-dimensional Linear Regressions”∗
Zijian Guo1 and Wanjie Wang1,2 T. Tony Cai1 Hongzhe Li2
1

Department of Statistics, The Wharton School, University of Pennsylvania
2
Department of Biostatistics and Epidemiology, Perelman School of
Medicine, University of Pennsylvania

Abstract
The supplementary materials to the paper “Optimal Estimation of Co-heritability
in High-dimensional Linear Regressions” present extended simulations in Section A.
In Section B, we prove (29) and (32) in Theorem 3. We also provide detailed proofs
of extra lemmas in Section C.

A

Simulation Results

A.1

Choice of λ0

In this section, we explore the choice of the tuning parameter λ0 for the scaled Lasso esb , which is required in the Lasso approach, and also in the initial step of
timators βb and γ
De-biased, Thresholded and FDE estimators. For these approaches, the theoretical tuning
√
parameter λ0 = 2.01 log p could be conservative in numerical studies. We would like to
choose a tuning parameter such that all the proposed estimators have a reasonable performance.
1

In this experiment, we implemented the same steps as Section 4, with estimators Lasso,
De-biased, Thresholded and FDE (FDE-S and FDE-NS for R (β, γ)) with tuning parameters
λ0 . Parameters were set the same as Experiment 1, that (p, n1 , n2 , rep) = (600, 400, 400, 300),
(s, s1 , s2 ) = (15, 30, 25), and the covariance matrices Σ and Γ satisfy Σij = Γij = (0.8)|i−j| .
For given values τ1 and τ2 , the signals of β satisfy that βji = (1 + i/s1 ) τ1 /2, for ji ∈ S1 ,
i = 1, 2, · · · , s1 , and the signals in γ satisfy that γj = τ2 for j ∈ S2 . We consider the
signal strength parameters as (τ1 , τ2 ) ∈ {(3.0, .1), (1.8, .4)}, including the settings that one
coefficient vector has much stronger signals than the other, and that the signals in two
coefficient vectors are close. To investigate the effect of tuning paramters, we take λ0 =
√
b 2.01 log p, and vary b ∈ {.25, .5, .75, 1}.
The results are summarized in Figure A1. For the estimation of I (β, γ), the results show
that as b increased from .25 to .5, the MSE of most approaches decreased except Lasso.
When b increased from .5 to 1, the MSE of Lasso increased significantly, and the MSE of
other estimators are relatively stable. For the estimation of R (β, γ), all the approaches
perform worse as b increases from .5 to 1. When b increases from .25 to .5, some approaches
work better, and the other approaches do not change much. Combining the results for both
I (β, γ) and R (β, γ), b = .5 is a generally good choice of tuning parameter such that most
√
approaches have almost the best performance. Therefore, we take λ0 = .5 2.01 log p in our
simulation studies.

A.2

Simulation results for quadratic functionals

In this section, we examine the performance of several estimators of the quadratic functional
Q(β), including the plug-in scaled Lasso estimator (Lasso), the plug-in de-biased estimator
(De-biased), the plug-in thresholded Lasso estimator (Thresholded), and the proposed estimator FDE with sample splitting (FDE-S) and without sample splitting (FDE-NS). Details
about these estimator are listed as follows.
b 2.
• Lasso: The quadratic functional kβk22 is estimated by kβk
2

e The quadratic functional
• De-biased: Denote the de-biased Lasso estimators as β.
e 2.
kβk2 is estimated by kβk
2

2

2

50
45
40

0.25

Lasso
De−biased
Thresholded
FDE

0.2

35

Lasso
De−biased
Thresholded
FDE−S
FDE−NS

0.15

30
25

0.1

20
15

0.05

10
5

0
0
0.25

0.5

0.75

1

0.25

Tuning Parameter

0.5

0.75

1

Tuning Parameter

12
11
10
9

Lasso
De−biased
Thresholded
FDE

0.35
0.3

8

Lasso
De−biased
Thresholded
FDE−S
FDE−NS

0.25

7
0.2

6
5

0.15

4
0.1
3
2

0.05

1
0.25

0.5

0.75

0
0.25

1

Tuning Parameter

0.5

0.75

1

Tuning Parameter

Figure A1: MSE for estimators Lasso (star), De-biased (square), Thresholded (circle), FDE
(upward-pointing triangle), FDE-S (upward-pointing triangle) and FDE-NS (downwardpointing triangle) for tuning parameters with b ∈ {0.25, 0.5, 0.75, 1}. Top left: Estimation of
I (β, γ) for (τ1 , τ2 ) ∈ (1.8, 0.4); Top right: Estimation of R (β, γ) for (τ1 , τ2 ) = (1.8, 0.4); Bottom left: Estimation of I (β, γ) for (τ1 , τ2 ) = (3, 0.1); Bottom right: Estimation of R (β, γ)
for (τ1 , τ2 ) = (3, 0.1).

3

• Thresholded: Denote the thresholded estimator as β̄. The quadratic functional kβk22
is estimated by kβ̄k22 .

• FDE: Estimate kβk22 as (16) in the main paper. We consider the method with data
splitting (FDE-S) and without data splitting (FDE-NS).
For each setting, with parameters (p, n, s), Σ and Fβ specified, the following steps were
implemented:
1. Generate sets S ⊂ [p], with |S| = s. For β ∈ Rp , generate βj ∼Fβ for j ∈ S1 , and set
βj = 0 for j 6∈ S1 .
i.i.d

i.i.d

2. Generate Xi· ∼ N (0, Σ), 1 ≤ i ≤ n. Generate the noise i ∼ N (0, 1), 1 ≤ i ≤ n,
and generate the outcome as y = Xβ + .
3. With X and y, estimate kβk22 through different estimators.
4. Repeat 2-3 for rep repetitions.
We use MSE in (33) in the main paper to measure the performance of each estimator. Similar
to Section 4, we implement Experiments 1 and 2.
Experiment 1. This experiment aimed to explore the effect of signal strength for two
types of Fβ . It contains two sub-experiments.
In Experiment 1a, the parameters were set as follows: (p, n, s, rep) = (600, 400, 30, 300)
and the covariance matrix Σ satisfies Σ(i, j) = (0.8)|i−j| . The signals of β follows that
βji = (1 + i/s)τ /2, i = 1, 2, · · · , s. To investigate how the signal strength affects the
estimation, vary the signal strength parameter τ ∈ {.1, .2, .3, .4, 1.8, 2.2, 2.6, 3}.
In Experiment 1b, we use the same setting, except that s = 25, Fβ is a point mass at τ ,
and τ ∈ {.1, .2, .3, .4, 1.0, 1.2, 1.4, 1.6}.
The results are summarized in Tables 5 and 6. FDE-NS has consistently outperformed
Lasso, De-biased and Thresholded when the signals are weak. The performances of FDE-NS,
Thresholded and De-biased were very similar when the signals are strong.
Experiment 2. This experiment aimed to examine the effect of signal sparsity on
the estimators. The parameters were set as follows: (p, n, rep) = (800, 400, 300), and the
4

Table 5: Estimation of the quadratic functional Q (β) in Experiment 1a: the truth and
MSE for the plug-in estimator with the scaled Lasso estimator (Lasso), plug-in estimator
with the de-biased estimator (De-biased), plug-in estimator with the thresholded estimator
(Thresholded), and the proposed estimators with sample splitting (FDE-S) and without
sample splitting (FDE-NS).
Strength parameter, τ
Method

.1

.2

.3

.4

1.8

2.2

2.6

3.0

Truth

.179

.715

1.609

2.860

57.920

86.522

120.844

160.888

Lasso

.014

.210

.882

2.224

95.226

147.032

210.047

273.340

4.430 4.333 3.773

3.044

4.953

10.301

17.790

26.021

De-biased
Thresholded

.096

.062

.055

.101

14.025

22.824

33.764

44.588

FDE-S

.024

.287 1.079

2.702

152.593

192.361

337.549

351.409

FDE-NS

.012

.015

.277

9.725

21.173

21.610

40.177

.100

Table 6: Estimation of the quadratic functional Q (β) in Experiment 1b: the truth and
MSE for the plug-in estimator with the scaled Lasso estimator (Lasso), plug-in estimator
with the de-biased estimator (De-biased), plug-in estimator with the thresholded estimator
(Thresholded), and the proposed estimators with sample splitting (FDE-S) and without
sample splitting (FDE-NS).
Strength parameter, τ
Method

.2

.3

.4

1.0

1.2

1.4

1.6

Truth

.250 1.000

2.000

4.000

25.000

36.000

49.000

64.000

Lasso

.025

.370

1.423

3.290

31.125

46.710

65.449

87.341

4.400

4.246

3.463

2.679

.822

1.324

2.387

4.011

Thresholded

.091

.057

.063

.146

3.553

5.812

8.632

12.012

FDE-S

.025

.370

1.423

3.290

31.125

46.710

65.449

87.341

FDE-NS

.009

.031

.132

.298

2.913

5.576

6.407

8.067

De-biased

.1

5

covariance matrix Σ has elements Σij = (0.8)|i−j| . Two settings of Fβ were considered: I.
the signals of β follows that βji = (1 + i/s)τ /2, i = 1, 2, · · · , s, τ = .2; and II. βj = τ ,
j ∈ S, τ = .1. To investigate the effect of signal sparsity, vary the number of signals
s ∈ {40, 50, 60, 70, 80, 90, 100, 110}.
The results are summarized in Table 7. FDE-NS has consistently outperformed other
methods for Settings I and II. In Setting II, where the signals were weaker than Setting I,
the improvement of FDE-NS over other methods was larger. Compared to other methods,
De-biased works much worse, which results from the large noise of the de-biased estimator.
Table 7: Estimation of the quadratic functional Q (β): the truth and MSE for the plug-in
estimator with the scaled Lasso estimator (Lasso), plug-in estimator with the thresholded
estimator (Thresholded), and the proposed estimators with sample splitting (FDE-S) and
without sample splitting (FDE-NS).
Sparsity parameter, s
Method

I

II

40

50

60

70

80

90

100

110

Truth

.5535

.6868

.8201

.9534

1.0868

1.2201

1.3534

1.4867

Lasso

.1347

.2179

.2688

.3609

.4634

.5285

.6236

.7331

7.1867

7.2886

7.4523

7.4733

7.6875

8.0140

8.1476

8.3544

Thresholded

.2012

.1717

.1997

.1827

.1779

.2052

.2072

.2331

FDE-S

.1788

.2843

.3573

.4949

.6078

.7237

.8509

1.0035

FDE-NS

.0134

.0264

.0296

.0527

.0790

.0919

.0954

.1074

Truth

.4000

.5000

.6000

.7000

.8000

.9000

1.0000

1.1000

Lasso

.0683

.1150

.1401

.1747

.2326

.2622

.2893

.2988

7.2395

7.5467

7.6953

7.9302

8.1280

8.5039

8.6979

8.9836

Thresholded

.2078

.1973

.2143

.2212

.2370

.2589

.2999

.3806

FDE-S

.0683

.1150

.1401

.1747

.2326

.2622

.2893

.2988

FDE-NS

.0057

.0105

.0113

.0113

.0148

.0189

.0264

.0234

De-biased

De-biased

6

B

Proof of (29) and (32) in Theorem 3

We will apply Lemma 5 and establish the lower bounds (29) and (32) in this section. Let
f (X, y, Z, w) denote the joint density function of X, y, Z, w, f (X, y) denote the joint
density function of X, y and f (Z, w) denote the joint density function of Z, w. The following lemma reduces the two sample L1 distance to one sample L1 distance. The proof of
the following lemma can be found in the supplementary material Section C.
Lemma 9. Assuming independence between (X, y) and (Z, w), then we have
L1 (fπ (X, y, Z, w) , fθ0 (X, y, Z, w)) ≤ L1 (fπ (X, y) , fθ0 (X, y))+L1 (fπ (Z, w) , fθ0 (Z, w)) .
(69)

Hence, by restricting to the cases where (X, y) and (Z, w) are independent, it is sufficient
to control the L1 distances L1 (fπ (X, y) , fθ0 (X, y)) and L1 (fπ (Z, w) , fθ0 (Z, w)). In the
following, we construct the parameter space H0 and H1 and apply Lemma 5 and Lemma 9.
The lower bound (32) in the main paper can be decomposed into the following three lower
bounds,
inf

sup

b θ∈Θ(k,η0 )
R

Pθ



b − R (β, γ) ≥ c0 1 k log p
R
η0 n



1
≥ ,
4

for η0 ≥ C

p
k log p/n,

(70)


p
1
0 1 1
b
inf sup Pθ R − R (β, γ) ≥ c √
(71)
≥ , for η0 ≥ C k log p/n,
b θ∈Θ(k,η0 )
η0 n
4
R



1
k
log
p
1
0
b − R (β, γ) ≥ c min
inf sup Pθ R
,1
≥ ,
(72)
2
b θ∈Θ(k,η0 )
η0 n
4
R
p
where c0 > 0 is a positive constant. Assuming η0 ≥ C k log p/n, combining (70), (71) and


(72), we can establish the following lower bound




1
1
k
log
p
1
k
log
p
1
0
0
b − R (β, γ) ≥ c
√ +
inf sup Pθ R
+c 2
≥ .
b θ∈Θ(k,η0 )
η0
n
η0 n
4
n
R

In the case η0 ≤ C

p
k log p/n, by (72), we establish the lower bound
inf

sup

b θ∈Θ(k,η0 )
R

Pθ




b − R (β, γ) ≥ c0 ≥ 1 .
R
4
7

(73)

(74)

Combining (73) and (74), we establish (32) in the main paper.
To facilitate the discussion of the lower bound (29) in the main paper, we introduce the
following three disjoint parameter spaces,
(
r

) 

k log p
M2
Θ1 (k) = Θ (k) ∩ kβk2 ≤ C
∩ kγk2 ≥
,
n
4
)
r

 (
M2
k log p
Θ2 (k) = Θ (k) ∩ kβk2 ≥
∩ kγk2 ≤ C
,
4
n
) (
)
(
r
r
k log p
k log p
∩ kγk2 ≤ C
.
Θ3 (k) = Θ (k) ∩ kβk2 ≤ C
n
n

(75)

The lower bound (29) in the main paper can be decomposed into the following three
lower bounds,



k log p
1
0
b
inf
sup
Pθ I − I (β, γ) ≥ c (kβk2 + kγk2 )
≥ ,
bI θ∈Θ1 (k)∪Θ2 (k)
n
4


1
1
0
inf
sup
Pθ bI − I (β, γ) ≥ c (kβk2 + kγk2 ) √
≥ ,
bI θ∈Θ1 (k)∪Θ2 (k)
4
n


k log p
1
inf sup Pθ bI − I (β, γ) ≥ c0
≥ ,
bI θ∈Θ3 (k)
n
4

where c0 > 0 is a positive constant. Combining (76) and (77), we can establish



1
k
log
p
1
0
inf
sup
Pθ bI − I (β, γ) ≥ c (kβk2 + kγk2 ) √ +
≥ .
bI θ∈Θ1 (k)∪Θ2 (k)
n
4
n

(76)

(77)
(78)

(79)

On Θ1 (k) ∪ Θ2 (k), we have






1
k log p
1
k log p
k log p
= max (kβk2 + kγk2 ) √ +
,
(kβk2 + kγk2 ) √ +
n
n
n
n
n
and on Θ3 (k), we have




k log p
1
k log p
k log p
= max (kβk2 + kγk2 ) √ +
,
.
n
n
n
n
Since both Θ1 (k) ∪ Θ2 (k) and Θ3 (k) are proper subspace of Θ (k), we can establish (29) by
combining (79) and (78).
We will first present the proofs of (70) and (76), which are similar in terms of construction
of the null and alternative parameter spaces. Then we will present proofs for (71) and (77)
8

and for (72) and (78). Define γ ∗ = (γ1∗ , 0, 0 · · · , 0) for some γ1∗ ≥ 0.
Proof of (70) and (76) We define the null and alternative spaces as
H0 = {(β ∗ , I, σ0 , γ ∗ , I, σ2 ) : β ∗ = (β1∗ , η0 , 0 · · · , 0)} ,
H1 = {(β, Ω, σ1 , γ ∗ , I, σ2 ) : (β, Ω, σ1 ) ∈ G1 } ,

(80)

where G1 is defined in (57) in the main paper. Let π denote the uniform prior over the
parameter space H1 , which is induced by the uniform prior of α over `1 (p, k, ρ) defined in
p
(56) in the main paper. Taking ρ = log(4p1 /k 2 )/2n, by Lemma 9 and Lemma 6, we

establish that

1
L1 (fπ , fθ0 ) ≤ .
4

(81)

1. To establish the lower bound for the ratio, we take β1∗ = 0 and then we have β1 =
−kαk22 σ0 /(1 − kαk22 ) and specify γ1∗ = η0 and hence H0 , H1 ⊂ Θ(k, η0 ). We calculate
the distance as
d=

k logn p
β1 η 0
− 0 ≥ cq
.
kβk2 |η0 |
log p
2
η +k
0

p
For η0 ≥ C k log p/n, we establish (70).

n

2. To establish the lower bound for the inner product, we take β1∗ = 0 and η0 = 0 and
p
then we have β1 = −kαk22 σ0 /(1 − kαk22 ), kβ ∗ k2 = 0 and kβk2 ≤ klog p/n. We
take γ1∗ = M2 /2 and hence H0 , H1 ⊂ Θ1 (k). We calculate
d=

|β1 γ1∗

kαk22 σ0
k log p
− 0| = kγ k2
≥ ckγ ∗ k2
.
2
1 − kαk2
n
∗

Applying Lemma 5, we establish


k log p
1
b
inf sup Pθ I − I (β, γ) ≥ ckγk2
≥ .
bI θ∈Θ1 (k)
n
4

By symmetry, we can also establish


k log p
1
b
inf sup Pθ I − I (β, γ) ≥ ckβk2
≥ .
bI θ∈Θ2 (k)
n
4

Combining (82) and (83), we establish (76).

9

(82)

(83)

Proof of (71) and (77) We define the null and alternative spaces as
H0 = {(β ∗ , I, σ0 , γ ∗ , I, σ2 ) : β ∗ = (β1∗ , η0 , 0 · · · , 0)} ,



¯
∗
∗
H1 = (β, I, σ0 , γ , I, σ2 ) : β = β1 + σ0 √ , η0 , 0 · · · , 0
.
n

(84)

By Lemma 9 and Lemma 7, we establish that
1
L1 (fπ , fθ0 ) ≤ .
4

(85)

1. To establish the lower bound for the ratio, we take β1∗ = 0 and γ1∗ = η0 and hence
H0 , H1 ⊂ Θ(k, η0 ). We calculate the distance as
1

√
β1 η 0
n
d=
− 0 ≥ cq
.
kβk2 |η0 |
1
2
η +c
0

n

p
√
For η0 ≥ C k log p/n  1/ n, we establish (71).

2. To establish the lower bound for the inner product, we take β1∗ = 0 and η0 = 0
√
and then we have kβ ∗ k2 = 0 and kβk2 ≤ 1/ n. We take γ1∗ = M2 /2 and hence
H0 , H1 ⊂ Θ1 (k). We calculate the distance as
σ0
d = |β1 γ1∗ − 0| ≥ ckγ ∗ k2 √ .
n
Applying Lemma 5, we establish


1
1
inf sup Pθ bI − I (β, γ) ≥ ckγk2 √
≥ .
bI θ∈Θ1 (k)
4
n

By symmetry, we can also establish


1
1
≥ .
inf sup Pθ bI − I (β, γ) ≥ ckβk2 √
bI θ∈Θ2 (k)
4
n

(86)

(87)

Combining (86) and (87), we establish (77).

Proof of (72) and (78) We introduce the following null and alternative,
H0 = {(β, I, σ1 , γ, I, σ2 ) : β = (η0 , 0, 0, · · · , 0) , γ = (0, η0 , 0, · · · , 0)} ⊂ Θ(k, η0 ) ⊂ Θ(k)
H1 = {(β, I, σ1 , γ, I, σ2 ) : β = (η0 , 0, α) , γ = (0, η0 , α) , with α ∈ `1 (p, k, ρ)} ⊂ Θ(k, η0 ) ⊂ Θ(k),
(88)
10

where `1 (p, k, ρ) is defined in (67) in the main paper. Let π denote the prior over the
parameter space H1 induced by the uniform prior of α over `1 (p, k, ρ) defined in (67) in the
p
main paper. Taking ρ = log(4p1 /k 2 )/2n, by Lemma 8, we have
L1 (fπ (y, X), fθ0 (y, X)) ≤

1
8

1
and L1 (fπ (w, Z), fθ0 (w, Z)) ≤ .
8

By Lemma 9, we have
1
L1 (fπ (y, X, w, Z), fθ0 (y, X, w, Z)) ≤ .
4

(89)

1. To establish the lower bound for the normalized inner product, we calculate the
distance d as
kαk22

d = p
2 − 0 ≥ c min
η02 + kαk22




1 k log p
,1 .
η02 n

(90)

Applying Lemma 5, we establish (71).

2. To establish the lower bound for the inner product, we take η0 = 0 and hence H0 , H1 ⊂
Θ3 (k). We calculate the distance
d = kαk22 − 0 ≥ c

k log p
.
n

Applying Lemma 5, we establish (78).

C

Supplementary Proof

In this section, we prove the extra lemmas. The proofs of Lemma 1 and 2 follow from the
results and analysis in Cai and Guo (2016a); Ren et al. (2013); Sun and Zhang (2012); Ye
and Zhang (2010).
Proof of Lemma 1
In the following, we establish upper bounds for kβb − βk1 , kβb − βk2 and the corresponding

proofs for kb
γ −γk1 and kb
γ −γk2 are similar and omitted here. The proof has been established
in Sun and Zhang (2012); Ye and Zhang (2010) for fixed designs under certain assumptions
for the design. In the following, we establish that the conditions are satisfied with high
11

probability for the Sub-gaussian random designs considered in this paper.
By normalizing the columns of X and the true sparse vector β, the linear regression model
can be expressed as
y = Hd + ,

with H = XD, d = D −1 β and  ∼ N (0, σ12 I),

(91)

√
where D = diag ( n/kX·j k2 )j∈[p] denotes the p × p diagonal matrix with (j, j) entry to be
√
√
n/kX·j k2 . Take δ0 = 1.0048, we have λ0 = (1.01) 2δ0 log p. Take 0 = 202, C1 = 2.25,
c0 = 1/6 and C0 = 3. Rather than use the constants directly in the following discussion, we
use δ0 , 0 , C1 , C0 and c0 to represent the above fixed constants in the following discussion.
We also assume that log p/n ≤ 1/25 and δ0 log p > 2. We review the following definitions of
restricted eigenvalue introduced in Bickel et al. (2009),
κ(X, k, α0 ) =

min

min

kXηk2
√
.
nkηJ0 k2

(92)

min

min

kXηk2
√
,
nkηJ01 k2

(93)

η6=0,
J0 ⊂{1,··· ,p},
kηJ c k1 ≤α0 kηJ0 k1
|J0 |≤k
0

and
κ(X, k, s, α0 ) =

η6=0,
J0 ⊂{1,··· ,p},
kηJ c k1 ≤α0 kηJ0 k1
|J0 |≤k
0

where J1 denotes the subset corresponding to the s largest in absolute value coordinates of
√
η outside of J0 and J01 = J0 ∪ J1 . Define σ1ora = ky − Xβk2 / n, and
√
√
kλ0
.
(94)
τ = 1 + 0 √
nκ(H, k, 1 + 20 )
Let Ω = Σ−1 denote the precision matrix of the design matrix. To facilitate the proof, we
define the following events for the random design X and the error ,


2 1
kX·j k2
7p
√
M1 for 1 ≤ j ≤ p ,
G1 =
< √
<
5 M1
5
n
(
)
r
(σ1ora )2
log p
log p
G2 =
−1 ≤2
+2
,
σ12
n
n
(
)
1
kXvk2
3
G3 =
≤√
≤
, where v = Ωb
γ,
1
2
2
b k2
nkΩ 2 γ
(
)
c
G4 = κ(X, k, k, α0 ) ≥ p
,
λmax (Ω)
12

G5 =
S1 =

(


|

kH k∞
≤ σ1
n

r

2δ0 log p
n

)

,


kH | k∞
0 − 1
ora
≤ σ λ0
(1 − τ ) ,
n
0 + 1

S2 = {(1 − ν0 ) σ̂ ≤ σ ≤ (1 + ν0 )σ̂1 } ,
n
p
√ o
b −γ
b | k∞ ≤ kb
B1 = kb
γ | ΩΣ
γ k2 λ1 / n , where λ1 = 4C0 M12 log p.
b = X | X/n and γ
b is the scaled Lasso estimator defined in (5) in the main paper.
where Σ

Define G = ∩5i=1 Gi and S = ∩2i=1 Si . We will check the assumptions in Corollary 1 in Sun
and Zhang (2012) and apply equation (23) in Sun and Zhang (2012). By the definition of
τ ∗ in Sun and Zhang (2012), we have τ ∗ ≤ τ where τ is defined in (94). Hence, on the event

S1 , equation (23) in Sun and Zhang (2012) holds. By equation (27) and the discussion after
(27) in Ye and Zhang (2010), we obtain
√
kλ0 σ1
b
.
kβ − βk2 ≤ C √ 2
nκ (H, k, k, 1 + 20 )

Similar to the proof of Lemma 13 in Cai and Guo (2016b), we establish



n
max kX·j k2
2
2
κ (H, k, k, 1 + 20 ) ≥
κ X, k, k, (1 + 20 )
.
max kX·j k22
min kX·j k2

(95)

(96)

Hence, on the event G ∩ S, we establish the upper bound for kβb − βk2 in (36) in the main

paper. By the discussion in Sun and Zhang (2012); Ye and Zhang (2010) or Lemma 5 in
Cai and Guo (2016a), we can establish the upper bound for kβb − βk1 in (35) in the main
paper. The remaining proof of Lemma 1 relies on the following lemma, which controls the
probability of events G, S and B1 .
Lemma 10. Suppose the assumption (A1) in the main paper holds and k log p/n → 0 with
k = max {kβk0 , kγk0 }. Then
P (G) ≥ 1 −

6
1
− 2p1−C1 − √
p1−δ0 − c0 exp (−cn) ,
p
2 πδ0 log p

(97)

and
1
P (G ∩ S) ≥ P (G) − 2 exp (−cn) − c00 √
p1−δ0 ,
log p
13

(98)

and
2

P (B1 ) ≥ 1 − 2p1−c0 C0 ,

(99)

where c, c0 and c00 are universal positive constants.
Proof of Lemma 2




b −γ
b |1 Σ
b | β − βb ,|b
In the following, we establish the upper bound for u
u|1 X | /n| and




b and the corresponding proofs for u
b − βb| (γ − γ
b) ,
b |2 Γ
b |1 X | y − X βb /n − hb
γ , β − βi
u

b γ
b |2 Z | (w − Z γ
b ) /n − hβ,
b − γi are similar and omitted here.
|b
u|2 Z | δ/n| and u
Note that






1
1
b =u
b −γ
b |1 X | y − X βb − hb
b |1 X |  + u
b |1 Σ
b | β − βb .
u
γ , β − βi
n
n

It is sufficient to establish (38) and (37) in the main paper and then (39) in the main paper
p
b −γ
b | k∞ ≤
follows. On the event S ∩ G ∩ B1 , we have kβb − βk1 ≤ Ck log p/n and kb
u|1 Σ




√
b −γ
b |1 Σ
b | β − βb in (37) in the main paper follows from the
kb
γ k2 λ1 / n. The control of u



b 1 and Lemma 1. Since u
b −γ
b −γ
b |1 Σ
b | k∞ kβ − βk
b 1 only
b | β − βb ≤ kb
inequality u
u|1 Σ

depends on X, Z and δ, which is independent of , we obtain


1
|1
|
2 2
b 1 X  | Z, X, δ ∼ N 0, 2 kX u
b 1 k2 σ 1
u
n
n
and hence


1
b 1 k2 σ1
≥ zα/2 kX u
P
n
n

Z 
1
|1
|
b 1 X  ≥ zα/2 kX u
b 1 k2 σ1 | Z, X, δ dZdXdδ = 1 − α.
= P u
n
n


1
b |1 X | 
u

(100)

b 1 defined
On the event B1 , Ωb
γ belongs to the feasible set of (10). By the definition of u
b 1 k22 /n ≤ kXΩb
in (10), we have kX u
γ k22 /n. On the event G3 , the term kXΩb
γ k22 /n can be

further upper bounded by CkΩb
γ k22 ≤ Ckb
γ k22 . Combined with (100), we establish (38) in

the main paper.

Proof of Lemma 3
The proof of this lemma is similar to that of Lemma 2. We only provide the sketches of the
proofs in the following. The proof of (40) in the main paper is similar to (38) in the main
14

paper and follows from the fact
2
b |3
u

n

X


(2) |



(2)

| X, Z, δ, 

(1)

∼N





4 |
(2)
(2) |
b3 .
b X
0, 2 u
X u
n 3

The proof of (41) in the main paper is similar to (37) in the main paper and follows from the



b 1 and the fact that kb
b (2) − βb| β − βb ≤ kb
b (2) − βb| k∞ kβ − βk
b (2) −
b |3 Σ
inequality u
u|3 Σ
u|3 Σ
p
b 2 log p/n.
βb| k∞ ≤ Ckβk

Proof of Lemma 4

By (22) in the main paper, on the event S ∩ G ∩ B1 , we have


p
p


q
Ckβk2 √1n + k log
+ C k log
n
n
1
k log p
1 k log p
b
q
Q (β) − kβk2 ≤
≤C √ +
+C
.
n
kβk2 n
n
b
Q (β) + kβk2
(101)
Similarly, on the event S ∩ G ∩ B1 , we have


q
1
k log p
1 k log p
b
.
Q (γ) − kγk2 ≤ C √ +
+C
n
kγk2 n
n

(102)

Combined with Assumption (A2) in the main paper, we establish that, on the event S∩G∩B1 ,
q
q
b
b (γ)
Q (β)
Q
1
1
−1 ≤
−1 ≤ .
and
(103)
kβk2
2
kγk2
2
By Lemma 1 and Assumption (A2) in the main paper, we can also establish that, on the
event S ∩ G ∩ B1 ,

b 2
1
kβk
−1 ≤
kβk2
2

and

kb
γ k2
1
−1 ≤ .
kγk2
2

(104)

By the triangle inequality, we have
βb
β
q
−
b (β) kβk2
Q

2

βb
βb
βb
β
= q
−
+
−
b (β) kβk2 kβk2 kβk2
Q

2

q
b (β) − kβk2
Q
1
b 2 q
≤ kβk
+
kβb − βk2 .
kβk2
b
Q (β)kβk2

By Assumption (A2) in the main paper, (101), (103) and (104), on the event S ∩ G ∩ B1 ,
βb

q
b (β)
Q

β
−
kβk2

≤C
2

1
kβk2



1
k log p
√ +
n
n



1 k log p
1
+
+
2
kβk2
kβk2 n

r

k log p
n

!

1
≤C
kβk2

r

k log p
.
n

(105)

15

Similarly, on the event S ∩ G ∩ B1 , we have
b
γ

q
b (γ)
Q

γ
−
kγk2

2

1
≤C
kγk2

r

k log p
.
n

(106)

Hence, (43) in the main paper follows from (105) and (106). Note that

*
+
*
+
b
b
b
β
β
µ
e1
β−β 
γ
µ
e1
γ
q
q
q
−
= q
− q
,q
,
+q
b (β) kβk2 kγk2
b (β) Q
b (γ)
b (β) Q
b (γ)
b (γ)
b (β)
Q
Q
Q
Q
Q
*
+ *
+
b−β
b
b
β
β
γ
β
γ
.
− q
,q
−
,q
−
kγk2
kβk2
b
b
b
Q (γ)
Q (β)
Q (β)
(107)

By (39) in the main paper, we have

and hence
µ
e1
q
q
−
b
b
Q (β) Q (γ)

b ≤ Ckb
γ k2
µ
e1 − hb
γ , β − βi
*

β − βb
q
,q
b (γ)
b (β)
Q
Q
b
γ

+



1
k log p
√ +
n
n



,



kb
γ k2
1
k log p
q
√ +
≤C q
n
n
b
b
Q (β) Q (γ)


k log p
1
1
√ +
≤C
,
kβk2
n
n

(108)

(109)

where the last inequality follows from (103) and (104). Note that
*
+ *
+
b
βb − β
βb
γ
γ
β
q
,q
−
,q
−
kγk2
kβk2
b
b
b
Q (γ)
Q (β)
Q (β)

 
+ (110)
*
+ *
b
b
b
b
γ
γ
β−β
γ  β−β   β
β 
= q
−
,q
+
, q
− q
−
.
kγk2
b (γ) kγk2
b (β)
b (β)
b (β) kβk2
Q
Q
Q
Q
We further have
+
*
b
b
γ
γ
βb − β
γ
γ
q
−
,q
≤ q
−
b (γ) kγk2
b (β)
b (γ) kγk2
Q
Q
Q
16

2

kβb − βk2
1
k log p
q
≤C
,
kβk2 kγk2 n
b
Q (β)
(111)

where the last inequality follows from (103) and (106). Note that
 

+
*
b
b
γ  β−β   β
1
β 
1
− q
, q
−
≤ kβk2 q
−
kγk2
b (β)
b (β) kβk2
b (β) kβk2
Q
Q
Q


1
1
k log p
1 k log p
√ +
≤C
+C
,
kβk2
n
kβk22 n
n

(112)

where the first inequality follows from Cauchy inequality and the last inequality follows from
(101) and (103). Combining (109), (110), (111) and (112), we establish (44) in the main
paper. Similarly, we establish (45) in the main paper.
Proof of Lemma 7 By (7.24) and (7.25) in Cai and Guo (2016a), if 2kβ ∗ − βk22 /σ02 <
log 2/2, then we have


− n2

2kβ ∗ − βk22
2nkβ ∗ − βk22
χ (fθ1 (y, X) , fθ0 (y, X)) = 1 −
− 1 ≤ exp
− 1, (113)
σ02
σ02
2

where the first equality follows from the moment generating function of χ2 distribution and
the second inequality follows from the inequality 1/(1 − x) ≤ exp(2x) for x ∈ [0, log 2/2] .
p
Taking ¯ = log (17/16) /2, we establish (65) in the main paper.

Proof of Lemma 9

By the definition of L1 distances and the independence between (X, y) and (Z, w), we have
Z
L1 (fπ (X, y, Z, w) , fθ0 (X, y, Z, w)) = |fπ (X, y, Z, w) − fθ0 (X, y, Z, w)| dXdydZdw
Z
= |fπ (X, y) fπ (Z, w) − fθ0 (X, y) fθ0 (Z, w)| dXdydZdw.
By the triangle inequality, the last term in the above inequality can be further upper bounded
by
Z

|fπ (X, y) (fπ (Z, w) − fθ0 (Z, w))| dXdydZdw +

Z

(|fπ (X, y) − fθ0 (X, y)) fθ0 (Z, w)| dXdydZdw

= L1 (fπ (X, y) , fθ0 (Z, w)) + L1 (fπ (Z, w) , fθ0 (Z, w)) ,

where the equality follows from
Proof of Lemma 10

R

|fπ (X, y)| dXdy = 1 and

R

(114)

|fθ0 (Z, w)| dZdw = 1.

The proof of Lemma 10 is a generalization of Lemma 4 in Cai and Guo (2016a). In the
following, we extend the Gaussian design in Cai and Guo (2016a) to Sub-gaussian design
17

considered in this paper. Since the error of the regression is still assumed to be Gaussian, it
is sufficient to establish the probability bound of G1 , G3 and G4 for the Sub-gaussian design
matrix. We first control the probability of G1 . By Corollary 5.17 in Vershynin (2012) and
the union bound, we have

 2
 



1
1


2
2
P max
kX·j k2 − EkX·j k2 ≥  ≤ 2p exp − min
n ,
,
1≤j≤p n
6
K2 K
p
where K = 4M1 . Taking  = 12M1 log p/n, we have
!
r

1
1
log
p
P max
kX·j k22 − EkX·j k22 ≥ 12M1
≤ 2p− 2 ,
1≤j≤p n
n

(115)

1

and hence P (G1 ) ≥ 1 − 2p− 2 .

By Theorem 1.6 in Zhou (2009), if n ≥ c0 M14 /θ2 ×max {9(2 + α0 )2 M1 k log (5ep/k) , 9 log p} ,

then with probability at least 1 − 2 exp (−cθ2 n/M14 ) , for all η such that there exists |J0 | ≤ k
and kηJ0c k1 ≤ α0 kηJ0 k1 ,
kXΩηk2
1−θ ≤ √
≤ 1 + θ,
1
nkΩ 2 ηk2

and 1 − θ ≤ √

kXηk2
1

nkΣ 2 ηk2

≤ 1 + θ.

(116)

By the discussion after equation (7) in Sun and Zhang (2012), the scaled Lasso estimator
p
b can be viewed as the solution of Lasso algorithm with tuning parameter 2.01 log p/nb
γ
σ2 .

Combined with Corollary B.2 in the paper Bickel et al. (2009), we can show that the Lasso
estimator satisfies the property that kb
γJ0c k1 ≤ α0 kb
γJ0 k1 where J0 = supp(γ) and hence
kXΩb
γ k2
−1 ≤
sup
√
1
b k2
nkΩ 2 γ
J0 ⊂{1,··· ,p},
|J0 |≤k

sup
η6=0,
kηJ c k≤α0 kηJ0 k1

kXΩηk2
−1 .
√
1
nkΩ 2 ηk2

(117)

0

b , if n ≥ Ck log p, then we establish P (G3 ) ≥ 1 − c exp(−c0 n)
By taking θ = 1/2 and η = γ
by the first half of (116) and (117). By the definition of (93) and the second half of (116),
we establish P (G4 ) ≥ 1 − 2 exp (−cn) .

References
Bickel, P. J., Ritov, Y., and Tsybakov, A. B. (2009), “Simultaneous analysis of Lasso and
Dantzig selector,” The Annals of Statistics, 37(4), 1705–1732.
18

Cai, T. T., and Guo, Z. (2016a), “Confidence intervals for high-dimensional linear regression:
minimax rates and adaptivity,” The Annals of Statistics, To appear.
Cai, T. T., and Guo, Z. (2016b), “Supplement to “Confidence intervals for high-dimensional
linear regression: minimax rates and adaptivity”,” The Annals of Statistics, To appear.
Ren, Z., Sun, T., Zhang, C.-H., and Zhou, H. H. (2013), “Asymptotic normality and optimalities in estimation of large gaussian graphical model,” arXiv preprint arXiv:1309.6024,
.
Sun, T., and Zhang, C.-H. (2012), “Scaled sparse linear regression,” Biometrika, 101(2), 269–
284.
Vershynin, R. (2012), “Introduction to the non-asymptotic analysis of random matrices,” in
Compressed Sensing: Theory and Applications, eds. Y. Eldar, and G. Kutyniok Cambridge
University Press, pp. 210–268.
Ye, F., and Zhang, C.-H. (2010), “Rate Minimaxity of the Lasso and Dantzig Selector for
the `q Loss in `r Balls,” The Journal of Machine Learning Research, 11, 3519–3540.
Zhou, S. (2009), “Restricted eigenvalue conditions on Subgaussian random matrices,” arXiv
preprint arXiv:0912.4045, .

19

