Honest confidence regions and optimality in
high-dimensional precision matrix estimation

arXiv:1507.02061v2 [] 20 Jul 2016

Jana Janková

Sara van de Geer

July 21, 2016

Abstract
We propose methodology for estimation of sparse precision matrices and statistical inference for their low-dimensional parameters in a high-dimensional setting
where the number of parameters p can be much larger than the sample size. We
show that the novel estimator achieves minimax rates in supremum norm and
the low-dimensional components of the estimator have a Gaussian limiting distribution. These results hold√uniformly over the class of precision matrices with
row sparsity of small order n/ log p and spectrum uniformly bounded, under
a sub-Gaussian tail assumption on the margins of the true underlying distribution. Consequently, our results lead to uniformly valid confidence regions for
low-dimensional parameters of the precision matrix. Thresholding the estimator
leads to variable selection without imposing irrepresentability conditions. The
performance of the method is demonstrated in a simulation study and on real
data.
Keywords: precision matrix sparsity inference asymptotic normality confidence
regions
Subject classification: 62J07 62F12

1

Introduction

We consider the problem of estimation of the inverse covariance matrix in a highdimensional setting, where the number of parameters p can significantly exceed the
sample size n. Suppose that we are given an n × p design matrix X, where the rows of
X are p-dimensional i.i.d. random vectors from an unknown distribution with mean
zero and covariance matrix Σ0 ∈ Rp×p . We denote the precision matrix by Θ0 := Σ−1
0 ,
assuming the inverse of Σ0 exists.
The problem of estimating the precision matrix arises in a wide range of applications.
Precision matrix estimation in particular plays an important role in graphical models
that have become a popular tool for representing dependencies within large sets of
variables. Suppose that we associate the variables X1 , . . . , Xp with the vertex set
V = {1, . . . , p} of an undirected graph G = (V, E) with an edge set E. A graphical
model G represents the conditional dependence relationships between the variables,
namely every pair of variables not contained in the edge set is conditionally independent given all remaining variables. If the vector (X1 , . . . , Xp ) is normally distributed,
1

each edge corresponds to a non-zero entry in the precision matrix (Lauritzen (1996)).
Practical examples of applications of graphical modeling include modeling of brain
connectivity based on FMRI brain analysis (Ng et al., 2013), genetic networks, financial data processing, social network analysis and climate data analysis.
A lot of work has been done on methodology for point estimation of precision matrices. We discuss some of the approaches below, but a selected list of papers includes
for instance Meinshausen and Bühlmann (2006); Friedman et al. (2008); Bickel and
Levina (2008); Yuan (2010); Cai et al. (2011); Sun and Zhang (2012). A common
approach assumes that the precision matrix is sufficiently sparse and employs the
`1 -penalty to induce a sparse estimator. The main goal of these works is to show
that, under some regularity conditions, the sparse estimator behaves almost as the
oracle estimator that has the knowledge of the true sparsity pattern.
Our primary interest in this paper lies not in point estimation, but we aim to quantify uncertainty of estimation by providing interval estimates for the entries of the
precision matrix. The challenge of this problem arises since asymptotics of regularized estimators which are the main tool in high-dimensional estimation is not easily
tractable (Knight and Fu, 2000), as opposed to the classical setting when the dimension of the unknown parameter is fixed.

1.1

Overview of related work

Methodology for inference in high-dimensional models has been mostly studied in
the context of linear and generalized linear regression models. From the work on
linear regression models, we mention the paper by Zhang and Zhang (2014) where a
semi-parametric projection approach using the Lasso methodology (Tibshirani, 1996)
was proposed, which was further developed and studied in van de Geer et al. (2013).
The approach leads to asymptotically normal estimation of the regression coefficients
and an extension of the method to generalized linear models is given in van de Geer
√
et al. (2013). The method requires sparsity of small order n/ log p in the highdimensional parameter vector and uses `1 -norm error bound of the Lasso. Further
alternative methods for inference in the linear model have been proposed and studied
in Javanmard and Montanari (2014), Belloni et al. (2014) and bootstrapping approach
was suggested in Chatterjee and Lahiri (2013), Chatterjee and Lahiri (2011).
Other lines of work on inference for high-dimensional models suggest post-model selection procedures, where in the first step a regularized estimator is used for model
selection and in the second step e.g. a maximum likelihood estimator is applied on
the selected model. In the linear model, simple post-model selection methods have
been proposed e.g. in Javanmard and Montanari (2013), Candes and Tao (2007).
These approaches are however only guaranteed to work under irrepresentability and
beta-min conditions (see Bühlmann and van de Geer (2011)). Especially in view of
inference, beta-min conditions which assume that the non-zero parameters are sufficiently large in absolute value, should be avoided.
In this paper, we consider estimation of precision matrices, which is a problem related to linear regression, however, it is a non-linear problem and thus it requires a
more involved treatment. One approach to precision matrix estimation is based on
regularization of the maximum likelihood in terms of the `1 -penalty. This approach is

2

typically referred to as the graphical Lasso, and has been studied in detail in several
papers, see e.g. Friedman et al. (2008), Rothman et al. (2008), Ravikumar et al.
(2008) and Yuan and Lin (2007). Another common approach to precision matrix
estimation is based on projections. This approach reduces the problem to a series
of regression problems and estimates each column of the precision matrix using a
Lasso estimator or Dantzig selector (Candes and Tao, 2007). The idea was first introduced in Meinshausen and Bühlmann (2006) as neighbourhood selection for Gaussian
graphical models and further studied in Yuan (2010), Cai et al. (2011) and Sun and
Zhang (2012).
Methodology leading to statistical inference for the precision matrix has been studied
only recently. The work Ren et al. (2015) proposes to use a more involved variation
of the regression approach to obtain an estimator which leads to statistical inference.
This approach leads to an estimator of the precision matrix which is elementwise
√
asymptotically normal, under row sparsity of order n/ log p, bounded spectrum
of the true precision matrix and Gaussian distribution of the sample. The paper
Janková and van de Geer (2015) proposes a method for statistical inference based
on the graphical Lasso. The work introduces a de-sparsified estimator based on the
graphical Lasso, which is also shown to be elementwise asymptotically normal.

1.2

Contributions and outline

We propose methodology leading to honest confidence intervals and testing for lowdimensional parameters of the precision matrix, without requiring irrepresentability
conditions or beta-min conditions to hold. Our work is motivated by the semiparametric approach in van de Geer et al. (2013) and is a follow-up of the work
Janková and van de Geer (2015). Compared to the previous work on statistical
inference for precision matrices, this methodology has several advantages. Firstly,
the estimator we propose is a simple modification of the nodewise Lasso estimator
proposed in Meinshausen and Bühlmann (2006). Hence the estimator is easy to implement and efficient solutions are available on the computational side. Secondly, the
novel estimator enjoys a range of optimality properties and leads to statistical inference under mild conditions. Firstly, the asymptotic distribution of low-dimensional
components of the estimator is shown to be Gaussian. This holds uniformly over
√
the class of precision matrices with row sparsity of order o( n/ log p), spectrum uniformly bounded in n and sub-Gaussian margins of the underlying distribution. This
results in honest confidence regions (Li, 1989) for low-dimensional parameters. The
proposed estimator achieves rate optimality as shown in Section 3.2. Moreover, the
de-sparsified estimator may be thresholded to guarantee variable selection without
imposing irrepresentable conditions. The computational cost of the method is order
O(p) Lasso regressions for estimation of all parameters and two Lasso regressions for
a single parameter.
The paper is organized as follows. Section 2 introduces the methodology. Section
3 contains the main theoretical results for estimation and inference and in Section
3.4 the suggested method is applied to variable selection. Section 4 provides a comparison with related work. Section 5 illustrates the theoretical results in a simulation

3

study. In Section 6, we analyze two real datasets and apply our method to variable
selection. Section 7 contains a brief summary of the results. Finally, the proofs were
deferred to the Supplementary material.
Notation. For a vector x = (x1 , . . . , xd ) ∈ Rd and p ∈ (0, ∞] we use kxkp to denote the p−norm of x in the classical sense. We denote kxk0 = |{i : xi 6= 0}|. For a
matrix A ∈ Rd×d we use the notations |||A|||∞ = maxi keTi Ak1 , |||A|||1 = AT ∞ and
kAk∞ = maxi,j |Aij |. The symbol vec(A) denotes the vectorized version of a matrix
A obtained by stacking the rows of A on each other. By ei we denote a p-dimensional
vector of zeros with one at position i. For real sequences fn , gn , we write fn = O(gn ) if
|fn | ≤ C|gn | for some C > 0 independent of n and all n > C. We write fn  gn if both
fn = O(gn ) and 1/fn = O(1/gn ) hold. Finally, fn = o(gn ) if limn→∞ fn /gn = 0. Furthermore, for a sequence of random variables xn we write xn = OP (1) if xn is bounded
in probability and we write xn = OP (rn ) if xn /rn = OP (1). We write xn = oP (1) if
xn converges in probability to zero.
P
Let
denote the convergence in distribution and → the convergence in probability.
Let Φ denote the cumulative distribution function of a standard normal random variable. By Λmin (A) and Λmax (A) we denote the minimum and maximum eigenvalue
of A, respectively. Let a ∨ b, a ∧ b denote max(a, b), min(a, b), respectively. We use
letters C, c to denotes universal constants. These are used in the proofs repeatedly
to denote possibly different constants.

2

De-sparsified nodewise Lasso

Our methodology is a simple modification of the nodewise Lasso estimator proposed in
Meinshausen and Bühlmann (2006). The idea is to remove the bias term which arises
in the nodewise Lasso estimator due to `1 -penalty regularization. This approach in
inspired by literature on semiparametric statistics Bickel et al. (1993); van der Vaart
(2000). We note several papers have used this idea in the context of high-dimensional
sparse estimation, see Zhang and Zhang (2014); van de Geer et al. (2013); van de Geer
(2016); Javanmard and Montanari (2014); Janková and van de Geer (2015).
We first summarize the nodewise Lasso method introduced in Meinshausen and
Bühlmann (2006) and discuss some of its properties. This method estimates an
unknown precision matrix using the idea of projections to approximately invert the
sample covariance matrix. For each j = 1, . . . , p we define the vector γj = {γj,k , k 6= j}
as follows
γj := arg min EkXj − X−j γk22 /n
(1)
γ∈Rp−1

and denote ηj := Xj − X−j γj and the noise level by τj2 = EηjT ηj /n. We define the
column vector Γj := (−γj,1 , . . . , −γj,j−1 , 1, −γj,j+1 , . . . , −γj,p )T . Then one may show
Θ0 = (Θ01 , . . . , Θ0p ) = (Γ1 /τ12 , . . . , Γp /τp2 ),

(2)

where Θ0j is the j-th column of Θ0 . Hence the precision matrix Θ0 may be recovered
from the partial correlations γj,k and from the noise level τj2 . In our problem, we are
4

only given the design matrix X. The idea of nodewise Lasso is to estimate the partial
correlations and the noise levels by doing a projection of every column of the design
matrix on all the remaining columns. In low-dimensional settings, this procedure
would simply recover the sample covariance matrix XT X/n. However, due to the
high-dimensionality of our setting, the matrix XT X/n is not invertible and we can
only do approximate projections. If we assume sparsity in the precision matrix (and
thus also in the partial correlations), this idea can be effectively carried out using the
Lasso. Hence, for each j = 1, . . . , p define the estimators of the regression coefficients,
γ̂j = {γ̂j,k , k = 1, . . . , p, j 6= k} ∈ Rp−1 , as follows
γ̂j := arg min kXj − X−j γk22 /n + 2λj kγk1 .
γ∈Rp−1

(3)

We further define the column vectors
Γ̂j := (−γ̂j,1 , . . . , −γ̂j,j−1 , 1, −γ̂j,j+1 , . . . , −γ̂j,p )T ,
and estimators of the noise level
τ̂j2 := kXj − X−j γ̂j k22 /n + λj kγ̂j k1 ,
for j = 1, . . . , p. Finally, we define the j-th column of the nodewise Lasso estimator
Θ̂ as
Θ̂j := Γ̂j /τ̂j2 .
(4)
The estimator Θ̂j of the precision matrix was studied in several papers (following
Meinshausen and Bühlmann (2006)) and has been shown to enjoy oracle properties
under mild conditions on the model. These conditions include bounded spectrum
of the precision matrix, row sparsity of small order n/ log p and a sub-Gaussian distribution of the rows of X (alternatively to sub-Gaussianity, one may assume that
the covariates are bounded as in van de Geer et al. (2013)). Our approach uses the
nodewise Lasso estimator as an initial estimator. The next step involves de-biasing or
de-sparsifying, which may be viewed as one step using the Newton-Raphson scheme
for numerical optimization. This is equivalent to “inverting” the Karush-Kuhn-Tucker
(KKT) conditions by the inverse of the Fisher information as in van de Geer et al.
(2013). The challenge then also comes from the need to estimate the Fisher information matrix which is a p2 × p2 matrix. We show that the estimator Θ̂ can be used in a
certain way to create a surrogate of the inverse Fisher information matrix. Since the
estimator Θ̂j can be characterized by its KKT conditions, it is convenient to work
with these conditions to derive the new de-sparsified estimator. Consider hence the
KKT conditions for the optimization problem (3)
− XT−j (Xj − X−j γ̂j )/n + λj κ̂j = 0,

(5)

for j = 1, . . . , p, where κ̂j is the sub-differential of the function γj 7→ kγj k1 at γ̂j , i.e.
(
sign(γ̂j,k )
if γ̂j,k 6= 0
κ̂j,k =
aj,k ∈ [−1, 1] otherwise,

5

where k ∈ {1, . . . , p} \ {j}. If we define Ẑj to be a p × 1 vector
Ẑj := (κ̂j,1 , . . . , κ̂j,j−1 , 0, κ̂j,j+1 , . . . , κ̂j,p )/τ̂j2 ,
then the KKT conditions may be equivalently stated as follows
Σ̂Θ̂j − ej − λj Ẑj = 0, for j = 1, . . . , p,

(6)

where Σ̂ = XT X/n is the sample covariance matrix. This is shown in Lemma 11 in
the Supplementary material. Consequently, the KKT conditions (6) imply a bound
kΣ̂Θ̂j − ej k∞ ≤ λj /τ̂j2 for each j = 1, . . . , p, which will be useful later. Note that the
KKT conditions may be equivalently summarized in a matrix form as Σ̂Θ̂−I−ẐΛ = 0,
where the columns of Ẑ are given by Ẑj for j = 1, . . . , p and Λ is a diagonal matrix
with elements (λ1 , . . . , λp ).
Multiplying the KKT conditions (6) by Θ̂i , we obtain
Θ̂Ti (Σ̂Θ̂j − ej ) − Θ̂Ti λj Ẑj = 0.
Then we note that adding Θ̂ij − Θ0ij to both sides and rearranging we get
Θ̂ij − Θ̂Ti λj Ẑj − Θ0ij

= Θ̂ij − Θ̂Ti (Σ̂Θ̂j − ej ) − Θ0ij
˜ ij ,
= −(Θ0 )T (Σ̂ − Σ0 )Θ0 + ∆
i

(7)

j

˜ ij = −(Θ̂i − Θ0 )T (Σ̂Θ̂j − ej ) − (Θ̂j − Θ0 )T (Σ̂Θ0 − ei ) is a term which can
where ∆
i
j
i
be shown to be oP (n−1/2 ) under certain conditions (Lemma 1). Hence we define the
de-sparsified nodewise Lasso estimator
T̂ := Θ̂ − Θ̂T (Σ̂Θ̂ − I) = Θ̂ + Θ̂T − Θ̂T Σ̂Θ̂.

3

(8)

Theoretical results

In this part, we inspect the asymptotic behaviour of the de-sparsified nodewise Lasso
estimator (8). In particular we consider the limiting distribution of individual entries
of T̂ and show the convergence to the Gaussian distribution is uniform over the considered model. For construction of confidence intervals, we consider estimators of the
asymptotic variance of the proposed estimator, both for Gaussian and sub-Gaussian
design. We derive convergence rates of the method in supremum norm and consider
application to variable selection.
For completeness, in Lemma 6 in the Supplementary material, we summarize convergence rates of the nodewise Lasso estimator. The result is essentially the same
as Theorem 2.4 in van de Geer et al. (2013) so the proof of the common parts is
omitted. Recall that V = {1, . . . , p} and we define the row sparsity by sj := kΘ0j k0 ,
maximum row sparsity by s := max1≤j≤p sj and the coordinates of non-zero entries
of the precision matrix by S0 := {(i, j) ∈ V × V : Θ0ij 6= 0}. For the analysis below,
we will need the following conditions.

6

A1 (Bounded spectrum) The inverse covariance matrix Θ0 := Σ−1
0 exists and there
exists a universal constant L ≥ 1 such that
1/L ≤ Λmin (Θ0 ) ≤ Λmax (Θ0 ) ≤ L.
A2 (Sparsity)

s log p
n

= o(1).

A3 (Sub-Gaussianity condition) Suppose that the design matrix X has uniformly
sub-Gaussian rows Xi , i.e. there exists a universal constant K such that

sup
E exp |αT Xi |2 /K 2 ≤ 2 (i = 1, . . . , n).
α∈Rp :kαk2 ≤1

The lower bound in A1 guarantees that the noise level τj2 = 1/Θ0jj does not diverge.
The upper bound (equivalently lower bound on eigenvalues of Σ0 ) guarantees that
the compatibility condition (see Bühlmann and van de Geer (2011)) is satisfied for
the matrix Σ0−j,−j , which is the true covariance matrix Σ0 without the j-th row
and j-th column. The sub-Gaussianity condition A3 is used to obtain concentration
results which are crucial to our analysis. Condition A3 is also used to ensure that the
compatibility condition is satisfied for Σ̂ with high probability (see Bühlmann and
van de Geer (2011)). Conditions A1, A2 and A3 are the same conditions as used in
van de Geer et al. (2013) to obtain rates of convergence for the nodewise regression
estimator. Define the parameter set
G(s) := {Θ ∈ Rp×p : max kΘi k0 ≤ s, A1 is satisfied}.
1≤i≤p

The following lemma shows that the proposed estimator T̂ can be decomposed into
√
a pivot term and a term which is of small order 1/ n with high probability.
Lemma 1. Suppose
q that Θ̂ is the nodewise Lasso estimator with regularization parameters λj ≥ c logn p , uniformly in j, for some sufficiently large constant c > 0.
Suppose that A2 and A3 are satisfied. Then for each (i, j) ∈ V × V it holds
√
√
n(T̂ij − Θ0ij ) = − n(Θ0i )T (Σ̂ − Σ0 )Θ0j + ∆ij ,
(9)
where there exists a constant C > 0 such that


slog p
lim sup P
max |∆ij | ≥ C √
= 0.
n→∞ Θ ∈G(s)
i,j=1,...,p
n
0
From Lemma 1 it follows that we need to assume stronger sparsity condition than
√
A2 for the remainder term ∆ij to be negligible after normalization by n. This is
accordance with other literature on the topic, see van de Geer et al. (2013), Ren et al.
(2015). Hence we introduce the following strengthened sparsity condition.
A2∗

s√
log p
n

= o(1).

7

The next result shows that the elements of T̂ are indeed asymptotically normal. To
this end, we further define the asymptotic variance
T

2
σij
:= var((Θ0i ) X1 X1T Θ0j ).

In some of the results to follow, we shall assume a universal lower bound on σij as
follows.
A4 There exists a universal constant ω > 0 such that σij ≥ ω.
Assumption A4 is satisfied e.g. under Gaussian design and A1. Denote a parameter
set
G̃(s) := {Θ ∈ Rp×p : max kΘi k0 ≤ s, A1, A4 are satisfied}.
1≤i≤p

Theorem 1 (Asymptotic normality).qSuppose that Θ̂ is the nodewise Lasso estimator

with regularization parameters λj  logn p uniformly in j. Suppose that A2∗ and A3
are satisfied. Then for every (i, j) ∈ V × V and z ∈ R it holds

√
n(T̂ij − Θ0ij )/σij ≤ z − Φ(z)| = 0.
lim sup |PΘ0
n→∞

Θ0 ∈G̃(s)

To construct confidence intervals, a consistent estimator of the asymptotic variance σij is required. Consistent estimators of σij are discussed in Section 3.1 (see
Lemmas 2 and 3). Hence Theorem 1 implies uniformly valid asymptotic confidence
√
intervals Iα := [T̂ij ± Φ−1 (1 − α/2)σ̂ij / n], i.e.
lim

n→∞

sup |PΘ0 (Θ0ij ∈ Iα ) − (1 − α)| = 0.
Θ0 ∈G̃(s)

The result also enables testing hypotheses about individual elements of the precision
matrix. For testing multiple hypothesis simultaneously, we may use the standard
procedures such as Bonferroni-Holm procedure (see van de Geer et al. (2013)).

3.1

Variance estimation

For the case of Gaussian observations, we may easily calculate the theoretical variance
and plug in the estimate Θ̂ in place of the unknown Θ0 as is displayed in Lemma 2
below.
Lemma 2. Suppose that assumptions A2 and A3 are satisfied and assume that the
rows of the design matrix X are independent
N (0, Σ0 )-distributed. Let Θ̂ be the
p
nodewise Lasso estimator and let λj ≥ cτ log p/n uniformly in j for some τ, c > 0.
2 := Θ̂ Θ̂ + Θ̂2 we have
Then for σ̂ij
ii jj
ij

sup P
Θ0 ∈G(s)

max

i,j=1,...,p

2
|σ̂ij

−

2
σij
|

p
≥ Cτ s log p/n

for some constants Cτ , c1 , c2 > 0.
8



≤ c1 p1−τ c2 ,

√
2 − σ 2 | = o (1/n1/4 ).
Lemma 2 implies that under s = o( n/ log p), we have a rate |σ̂ij
P
ij
If Gaussianity is not assumed, we may replace the estimator of the variance with
the empirical version, and plug in Θ̂ in place of the unknown Θ0 . Thus we take the
2 , where Θ̂ is the nodewise regression estimator
following estimator of σij
n

2
σ̂ij

1X T
:=
(Θ̂i Xk XkT Θ̂j )2 − Θ̂2ij .
n

(10)

k=1

The following lemma justifies this procedure under A1, A2∗ , A3.
Lemma 3. Suppose that the assumptions A2∗ and A3 are satisfied and for some
4
1− = 0. Let Θ̂ be the nodewise Lasso
 > 0, it holds that limn→∞
p log (p ∨ n)/n
estimator and let λj ≥ cτ log p/n uniformly in j for some τ, c > 0. Let σ̂ij be the
estimator defined in (10). Then for all η > 0


2
2
lim sup P
max |σ̂ij − σij | ≥ η = 0.
n→∞ Θ ∈G(s)
0

3.2

i,j=1,...,p

Rates of convergence

The de-sparsified estimator achieves optimal rates of convergence in supremum norm.
Observe first that for the nodewise regression estimator it holds by (7), Lemma 1 and
Lemma 10 in the Supplementary material that



log p
1
0
T
.
Θ̂ij − Θij = Θ̂i (Σ̂Θ̂j − ej ) + OP max √ , s
n
n
By Hölder’s inequality and the KKT conditions it follows
|Θ̂Ti (Σ̂Θ̂j − ej )| ≤ λj kΘ̂i k1 /τ̂j2 .
Consequently, for the rates of convergence of the nodewise Lasso in supremum norm
we find



log p
2 1
.
kΘ̂ − Θ0 k∞ = OP max
max λj kΘ̂i k1 /τ̂j , √ , s
i,j=1,...,p
n
n
De-sparsifying the estimator Θ̂ as in (8) removes the term involving λj kΘ̂i k1 /τ̂j2 in
the above rates.
Theorem 2 (Rates of convergence). Assume that A2 and A3 are satisfied. Let τ > 0
and let T̂pbe the de-sparsified nodewise Lasso estimator with regularization parameters
λj ≥ cτ log p/n for some sufficiently large constant c > 0, uniformly in j. Then
there exist constants Cτ , c1 , c2 > 0 such that



1
log p
sup P |T̂ij − Θ0ij | ≥ Cτ max √ , s
≤ c1 e−c2 τ .
n
n
Θ0 ∈G(s)
and

(r
sup P kT̂ − Θ0 k∞ ≥ Cτ max
Θ0 ∈G(s)

9

log p log p
,s
n
n

)!
≤ c1 p1−c2 τ .

We compare the results of Theorem 2 with results on optimal rates of convergence
derived for Gaussian graphical models in Ren et al. (2015). Suppose that the observations are Gaussian, i.e. X1 , . . . , Xn ∼ N (0, Σ0 ). For s ≤ C0 n/ log p for some
C0 > 0 and p ≥ sν for some ν > 2 it holds (see Ren et al. (2015))



1
log p
0
√
> c1 > 0,
inf sup P |Θ̂ij − Θij | > max C1
, C2 s
n
n
Θ̂ij Θ0 ∈G(s)
and
(
inf sup P kΘ̂ − Θ0 k∞ > max
Θ̂ Θ0 ∈G(s)

C10

r

log p 0 log p
, C2 s
n
n

)!
> c2 > 0,

where C1 , C2 , C10 , C20 are positive constants depending on ν and C0 only. As follows
from Theorem 2, the de-sparsified nodewise Lasso attains the lower bound on rates
and thus is in this sense optimal (considering the Gaussian setting).

3.3

Other de-sparsified estimators

The de-sparsification may work for other estimators of the precision matrix, provided that certain conditions are satisfied. This is formulated in Lemma 4 below. A
particular example of interest is the square-root nodewise Lasso estimator, that will
be discussed below. This estimator has the advantage that it is self-scaling in the
variance, similarly as the square-root Lasso (Belloni et al., 2011) on which it is based.
Lemma 4. Assume that for some estimator Ω̂ = (Ω̂1 , . . . , Ω̂p ) it holds
p
p
max kΩ̂j − Θ0j k1 = OP (s log p/n), kΣ̂Ω̂ − Ik∞ = OP ( log p/n).
j=1,...,p

(11)

Then for T̂ := Ω̂ + Ω̂T − Ω̂T Σ̂Ω̂ it holds under A1, A2∗ , A3
p
kT̂ − Θ0 k∞ = OP (max{s log p/n, log p/n}).
Moreover,

√

n(T̂ij − Θ0ij )/σij

N (0, 1).

We briefly consider nodewise regression with the square-root Lasso as an example.
The square-root Lasso estimators may be defined via
γ̂j := arg min kXj − X−j γk2 /n + 2λ0 kγk1 ,
γ∈Rp−1

for j = 1, . . . , p. Define τ̂j2 := kXj − X−j γ̂j k22 /n and τ̃j2 := τ̂j2 + λ0 τ̂j kγ̂j k1 . The
nodewise square-root Lasso is then given by Θ̂j,sqrt := Γ̂j /τ̃j2 , where
Γ̂j := (−γ̂j,1 , . . . , −γ̂j,j−1 , 1, −γ̂j,j+1 , . . . , −γ̂j,p )T .
Note that compared to the nodewise Lasso, the difference lies in estimation of the
partial correlations, where we used the square-root Lasso that “removes the square”

10

from the squared loss. The Karush-Kuhn-Tucker conditions similarly as in Lemma
11 in Supplementary material give
Σ̂Θ̂j,sqrt − ej −

τ̂j
λ0 Z̃j = 0,
τ̃j2

where Z̃j := (κ̃j,1 , . . . , κ̃j,j−1 , 1, κ̃j,j+1 , . . . , κ̃j,p )T and κ̃j,i is the sub-differential of the
function β 7→ kβk1 with respect to βi , evaluated at γ̂j . The paper Belloni et al. (2011)
further shows that the `1 -rates for the square-root Lasso satisfy condition (11). The
de-sparsified estimator may then be defined in the same way as in (8). Then the
conditions of Lemma 4 are satisfied and this implies that a de-sparsified nodewise
square-root Lasso achieves the same rates as the de-sparsified nodewise Lasso and
thus is also rate-optimal.

3.4

The thresholded estimator and variable selection

The de-sparsified estimator can be used for variable selection without imposing irrepresentable conditions. Under mild conditions, the procedure leads to exact recovery
of the coefficients that are sufficiently larger in absolute value than the noise level.
The following corollary is implied by Theorem 2 and Lemmas 2 and 3.
Corollary 1. Let Θ̂ be obtainedpusing the nodewise Lasso and T̂ be defined as in (8)
with tuning parameters λj ≥ cτ log p/n uniformly in j, for some c, τ > 0. Assume
that conditions A1, A2∗ , A3 and A4 are satisfied. Let σ̂ij , i, j = 1, . . . , p be the
estimator from Lemma 3 and assume that log4 (p ∨ n)/n1− = o(1) for some  > 0.
Then there exists some constant Cτ > 0 such that
p
lim P( max |T̂ij − Θ0ij |/σ̂ij ≥ Cτ log p/n) = 0.
n→∞

i,j=1,...,p

If, in addition, the rows of X are N (0, Σ0 )-dsitributed and σ̂ij , i, j = 1, . . . , p is instead
the estimator from Lemma 2, then there exist constants c1 , c2 , Cτ such that
p
P( max |T̂ij − Θ0ij |/σ̂ij ≥ Cτ log p/n) ≤ c1 p1−c2 τ .
i,j=1,...,p

Corollary 1 implies that we may define the re-sparsified estimator
T̂ijthresh := T̂ij 1|T̂

ij |>Cτ σ̂ij

√

log p/n

,

where σ̂ij is defined as in Corollary 1. Denote Ŝ thresh := {(i, j) ∈ V × V : T̂ijthresh 6= 0}.
p
Denote S0act := {(i, j) ∈ V × V : |Θ0ij | ≥ 2Cτ σij log p/n}. Then it follows directly
from Corollary 1 that with high probability
S0act ⊂ Ŝ thresh ⊂ S0 .
The inclusion S0act ⊂ Ŝ thresh represents that T̂ thresh correctly identifies all the nonzero parameters which are above the noise level. The inclusion Ŝ thresh ⊂ S0 means
that there are no false positives. If for all (i, j) ∈ S0 it holds
p
(12)
|Θ0ij | ≥ 2Cτ σij log p/n,
then we have exact recovery, i.e. with high probability: Ŝ thresh = S0 .
11

4

Further comparison to previous work

Closely related is the paper Janková and van de Geer (2015), where asymptotically
normal estimation of elements of the concentration matrix is considered based on
the graphical Lasso. While the analysis follows the same principles, the estimation
method used here does not require the irrepresentability condition (Ravikumar et al.,
2008) that is assumed in Janková and van de Geer (2015). Hence we are able to show
that our results hold uniformly over the considered class. Furthermore, regarding the
computational cost, our method uses Lasso regressions, which can be implemented
using fast algorithms as in Efron et al. (2004). In comparison, the graphical Lasso
method presents a more challenging computational problem, for more details see e.g.
Mazumder and Hastie (2012).
Another related work is the paper Ren et al. (2015). This paper suggests an estimator for the precision matrix which is shown to have one dimensional asymptotically
normal components with asymptotic variance to Θ0ii Θ0jj + (Θ0ij )2 . The assumptions
and results used in the paper are essentially identical with our assumptions and theoretical results in the present paper. However, there are some differences. The paper
Ren et al. (2015) assumes Gaussianity of the underlying distribution, while we only
require sub-Gaussianity of the margins. Another difference is in the construction of
the estimators. Both approaches use regression to estimate the elements of the precision matrix, but the paper Ren et al. (2015) concentrates on estimation of the joint
distribution of each pair of variables (X i , X j ) for i, j = 1, . . . , p. Thus it is computationally more intensive as it requires O(ps) high-dimensional regressions (see Ren
et al. (2015)), while our methodology only requires O(p).

5

Simulation results

In this section we report on the performance of our method on simulated data and
provide a comparison to another methodology. The random sample X1 , . . . , Xn satisfies EXi = 0, var(Xi ) = Θ−1
0 , where the precision matrix Θ0 = five-diag(ρ0 , ρ1 , ρ2 )
is defined by


ρ0 if i = j,



ρ
if |i − j| = 1,
1
Θ0ij =

ρ2 if |i − j| = 2,



0
otherwise.
We consider the settings S1 = (ρ0 , ρ1 , ρ2 ) = (1, 0.3, 0) and S2 = (ρ0 , ρ1 , ρ2 ) =
(1, 0.5, 0.3). The second setting (1, 0.5, 0.3) is further adjusted by randomly perturbing each non-zero off-diagonal element of Θ0 by adding a realization from the uniform
distribution on the interval [−0.05, 0.05]. We denote this new perturbed model by
(1, 0.5, 0.3)U . Hence the second precision matrix was chosen randomly. The sparsity
√
assumption requires s = o( n/ log p). We have chosen the sample sizes for numerical experiments according to the sparsity assumption (for this purpose, we ignored
possible constants in the sparsity restriction), i.e. n ≥ s2 log2 p.

12

5.1
5.1.1

Asymptotic normality and confidence intervals for individual
parameters
The Gaussian setting

In this section, we consider normally-distributed observations, X i ∼ N (0, Θ−1
0 ), for
√
i = 1, . . . , n. In Figure 1 we display histograms of n(T̂ij − Θ0ij )/σ̂ij for (i, j) ∈
{(1, 1), (1, 2), (1, 3)}, where T̂ is defined in (8) and the empirical variance σ̂ij is estimated as suggested by Lemma 2. Superimposed is the density of N (0, 1).
Secondly, we investigate the properties of confidence intervals constructed using the
de-sparsified nodewise Lasso. For comparison, we also provide results using confidence intervals based on the de-sparsified graphical Lasso introduced in Janková and
van de Geer (2015). The coverage and length of the confidence interval were estimated
by their empirical versions,
√
α̂ij := PN 1{Θ0,ij ∈Iij,α } and `ˆij := PN 2Φ−1 (1 − α/2)σ̂ij / n,
respectively, using N = 300 random samples. For a set A ⊂ V × V , we define the
average coverage over the set A (and analogously average length avglengthA ) as
X
avgcovA :=
α̂ij /|A|.
(i,j)∈A

We report average coverages over the sets S0 and S0c . These are denoted by avgcovS0
and avgcovS0c , respectively. Similarly, we calculate average lengths of confidence intervals for each parameter Θ0ij from N = 300 iterations and report avglengthS0 and
avglengthS0c .
The results of the simulations are shown in Tables 1 and 2. The target coverage level
is 95%. The methodology for the choice of the tuning parameters was used as follows
(see Ren et al. (2015)), for both methods,
p
√
ŝ = n/ log p, B = qt(1 − ŝ/(2p), n − 1), λ = B/ n − 1 + B 2 ,
(13)
where qt(β, n − 1) denotes the β-quantile of a t-distribution with n − 1 degrees of
freedom.
5.1.2

A sub-Gaussian setting

In this section, we consider a design matrix with rows having a sub-Gaussian distribution other than the Gaussian distribution. Let U := (U1 , . . . , Un ) be an n × p matrix
with jointly independent
entries generated from a continuous uniform distribution on
√ √
the interval [− 3, 3]. Further consider a matrix Θ0 := five-diag(1, 0.3, 0) and let
Σ0 = Θ−1
0 . Then we define
1/2
Xi := Σ0 Ui
for i = 1, . . . , n. Then the expectation of Xi is zero and the covariance matrix of Xi
is exactly Σ0 and the precision matrix is Θ0 . It follows by Hoeffding’s inequality that
Xi defined as above is sub-Gaussian with a universal constant K > 0.
A further difference compared to the simulations in Section 5.1.1 is that we now
13

0.3

Density

0.2

0.3
Density

0.2

0.4
0.3
−3

−2

−1
0
1
^
^ 11
n (T11 − θ11) σ

2

3

0.1
0.0

0.0

0.0

0.1

0.1

0.2

Density

0.5

0.4

0.4

0.6

0.5

Asymptotic normality in the Gaussian setting

−3

−2

−1
0
1
^
^ 12
n (T12 − θ12) σ

2

3

−3

−2

−1
0
1
^
^ 13
n (T13 − θ13) σ

2

3

√

n(T̂ij − Θ0ij )/σ̂ij , (i, j) ∈ {(1, 1), (1, 2), (1, 3)}. The sample size was
n = 500 and the number of parameters p = 100. The nodewise regression estimator was calculated
300 times. The setting is S1 = (1, 0.3, 0).

Figure 1: Histograms for

Gaussian setting: Estimated coverage probabilities and lengths
Setting S1 = (1, 0.3, 0)

S0
avgcov

S0
avglength

S0c
avgcov

S0c
avglength

p

n

100

191

D-S NW
D-S GL

0.945
0.931

0.302
0.293

0.963
0.974

0.262
0.254

200

253

D-S NW
D-S GL

0.947
0.928

0.267
0.254

0.963
0.976

0.232
0.220

300

293

D-S NW
D-S GL

0.949
0.928

0.238
0.236

0.965
0.977

0.220
0.205

400

324

D-S NW
D-S GL

0.948
0.925

0.246
0.228

0.965
0.981

0.230
0.223

Table 1: A table showing a comparison of de-sparsified nodewise Lasso (D-S NW) and de-sparsified
graphical Lasso (D-S GL). Parameter p takes values 100, 200, 300, 400 and the corresponding values
n are given by n = ds2 log2 pe, where s = 3. The regularization parameter was chosen as described
in (13). The number of generated random samples was N = 300.

estimate the variance of the de-sparsified estimator using the formula proposed in
(10) for sub-Gaussian settings:
n

2
σ̂ij

1X T
:=
(Θ̂i Xk XkT Θ̂j )2 − Θ̂2ij ,
n

(14)

k=1

where Θ̂ is the nodewise Lasso. The regularization parameters for the nodewise Lasso
are used in accordance with (13). Figure 2 again displays the histograms related to
several entries of the de-sparsified nodewise Lasso. Results related to the constructed
confidence intervals are summarized in Table 3. The results demonstrate that the
14

Gaussian setting: Estimated coverage probabilities and lengths
Setting S2 = (1, 0.5, 0.3)U

S0
avgcov

S0
avglength

S0c
avgcov

S0c
avglength

D-S NW
D-S GL

0.896
0.781

0.164
0.153

0.975
0.980

0.146
0.137

702

D-S NW
D-S GL

0.868
0.729

0.142
0.133

0.976
0.982

0.126
0.119

300

814

D-S NW
D-S GL

0.863
0.712

0.131
0.124

0.976
0.984

0.117
0.110

400

898

D-S NW
D-S GL

0.859
0.709

0.125
0.118

0.976
0.984

0.111
0.105

p

n

100

531

200

Table 2: A table showing a comparison of de-sparsified nodewise Lasso (D-S NW) and the desparsified graphical Lasso (D-S GL). Parameter p takes values 100, 200, 300, 400 and the corresponding
values n are given by n = ds2 log2 pe, where s = 5. The regularization parameter was chosen as
described in (13). The number of generated random samples was N = 300.

de-sparsified nodewise Lasso performs relatively well even under this non-Gaussian
setting.

−3

−2

−1 0
1
2
^
^ 11
n (T11 − θ11) σ

3

0.5
Density

0.0

0.1

0.2

0.3

0.4

0.5
0.3

Density

0.0

0.0

0.1

0.2

0.2
0.1

Density

0.3

0.4

0.4

Asymptotic normality in the sub-Gaussian setting

−3

−2

−1 0
1
2
^
^ 12
n (T12 − θ12) σ

3

−3

−2

−1 0
1
2
^
^ 13
n (T13 − θ13) σ

3

√

n(T̂ij − Θ0ij )/σ̂ij , (i, j) ∈ {(1, 1), (1, 2), (1, 3)}. The sample size was
n = 500 and the number of parameters p = 100. The nodewise regression estimator was calculated
300 times. The setting is S1 = (1, 0.3, 0).

Figure 2: Histograms for

5.2

Variable selection

For variable selection as suggested in Corollary 1, we compare the de-sparsified nodewise Lasso and the de-sparsified graphical Lasso. The setting is again as in Section
5.1.1. Average true positives and false positives over 100 repetitions are reported.
Choice of the tuning parameters is according to (13) and the thresholding level is

15

Sub-Gaussian setting: Estimated coverage probabilities and lengths
Setting S1 = (1, 0.3, 0)

S0
avgcov

S0
avglength

S0c
avgcov

S0c
avglength

D-S NW
D-S GL

0.906
0.811

0.234
0.190

0.949
0.944

0.249
0.216

253

D-S NW
D-S GL

0.909
0.791

0.203
0.165

0.950
0.946

0.217
0.187

300

293

D-S NW
D-S GL

0.911
0.765

0.189
0.152

0.950
0.947

0.202
0.173

400

324

D-S NW
D-S GL

0.911
0.740

0.180
0.143

0.951
0.947

0.192
0.164

p

n

100

191

200

Table 3: A table showing a comparison of de-sparsified nodewise Lasso (D-S NW) and the desparsified graphical Lasso (D-S GL). Parameter p takes values 100, 200, 300, 400 and the corresponding
values n are given by n = ds2 log2 pe, where s = 3. The regularization parameter was chosen as
described in (13). The number of generated random samples was N = 300.

given by
r

log p
,
(15)
n
taking ν = 1 for the de-sparsified nodewise regression, ν = 0.5 for the de-sparsified
graphical Lasso. We take σ̂ij = Θ̂ii Θ̂jj + Θ̂2ij as in Lemma 2. The results of this
simulation experiment are summarized in Table 4.
λthresh = σ̂ij

2ν

Estimated true positives (TP) and false positives (FP)
Setting S1 = (1, 0.5, 0.4)

TP

TP rate %

FP

FP rate %

p = 100
|S0 | = 494

D-S NW
D-S GL

494
493.98

100.0
99.999

0
0

0
0

p = 200
|S0 | = 994

D-S NW
D-S GL

994
993.62

100.0
99.961

0
0

0
0

p = 300
|S0 | = 1494

D-S NW
D-S GL

1494
1492.42

100.0
99.894

0
0

0
0

p = 400
|S0 | = 1994

D-S NW
D-S GL

1994.00
1989.08

100.0
99.753

0
0

0
0

Table 4: Estimated true positives (TP) and false positives (FP) for the de-sparsified nodewise
regression estimator (D-S NW) and for the D-S GL estimator. The sample size n = 400 was held
constant for all the values of p; the number of repetitions was N = 100. The thresholding levels was
chosen as in (15).

6

Real data experiments

We consider two real datasets, where we model the conditional independence structure
of the covariates using a graphical model. In particular, we aim to do edge selection
16

and we estimate the edge structure of the graphical model using the de-sparsified
nodewise Lasso. The first dataset is the Prostate Tumor Gene Expression dataset,
which is available in the R package spls. The second dataset is about riboflavin
(vitamin B2 ) production by bacillus subtilis. The dataset is available from the R
package hdi.
For both datasets, the procedure is essentially identical. We only consider the first 500
covariates which have the highest variances. In the first step, we split the sample and
use 10 randomly chosen observations to estimate the variances of the 500 variables.
With the estimated variances, we scale the design matrix containing the remaining
observations. We calculate the nodewise Lasso using the tuning parameter as in the
simulation study, and then calculate the de-sparsified nodewise Lasso. We threshold
√
the de-sparsified nodewise Lasso at the level Φ−1 (1−α/(2p2 ))σ̂ij / n, where α = 0.05
2 = Θ̂ Θ̂ + Θ̂2 is an estimate of the asymptotic variance calculated under the
and σ̂ij
ii ii
ij
assumption of normality and using the nodewise Lasso estimator Θ̂.
The first dataset contained observations of p = 4088 logarithms of genes expression
levels from n = 71 genetically engineered mutants of bacillus subtilis. We
 considered
500
500 variables with the highest variances, hence a full graph contains 2 edges. The
de-sparsified nodewise Lasso identified 20 edges as significant. For comparison, the
de-sparsified graphical Lasso introduced in Janková and van de Geer (2015) identified
5 edges as significant. It is worth pointing out that the set of edges selected by the desparsified graphical Lasso is a subset of the edges selected by de-sparsified nodewise
Lasso.
The second dataset contained n = 102 observations on p = 6033 variables. We used
the procedure above to do edge selection using the de-sparsified nodewise Lasso. Our
analysis identified 108 edges as significant using the de-sparsified nodewise Lasso.
For comparison, the de-sparsified graphical Lasso identified 28 edges as significant.
Again, the set of edges selected by the de-sparsified graphical Lasso is a subset of the
edges selected by de-sparsified nodewise Lasso.

7

Conclusions

We proposed a methodology for low-dimensional inference in high-dimensional graphical models. The method, called the de-sparsified nodewise Lasso, is easy to implement
and computationally competitive with the state-of-art methods. We studied asymptotic properties of the de-sparsified nodewise Lasso under mild conditions on the
model. The de-sparsified nodewise Lasso enjoys rate optimality in supremum norm
and leads to exact variable selection under beta-min conditions and mild conditions
on the model. We demonstrated its performance on several models in a simulation
study and on two real datasets. These numerical studies showed that it performs
well in a variety of settings, including non-Gaussian settings. Further open questions
concern for instance the asymptotic efficiency of the proposed estimator, similarly as
in the low dimensional settings.

17

References
Belloni, A., Chernozhukov, V., and Hansen, C. (2014). Inference on treatment effects
after selection amongst high-dimensional controls. Rev. Econ. Stud., 81(2):608–650.
Belloni, A., Chernozhukov, V., and Wang, L. (2011). Square-root Lasso: Pivotal
recovery of sparse signals via conic programming. Biometrika, 98(4):791–806.
Bickel, P. J., Klaassen, C. A., Ritov, Y., and Wellner, J. A. (1993). Efficient and
adaptive estimation for semiparametric models. Springer.
Bickel, P. J. and Levina, E. (2008). Covariance regularization by thresholding. Ann.
Statist., 36(6):2577–2604.
Bühlmann, P. and van de Geer, S. (2011). Statistics for high-dimensional data.
Springer.
Cai, T., Liu, W., and Luo, X. (2011). A constrained l1 minimization approach to
sparse precision matrix estimation. J. Amer. Statist. Assoc., 106:594–607.
Candes, E. and Tao, T. (2007). The dantzig selector: Statistical estimation when p
is much larger than n. Ann. Statist., 35(6):2313–2351.
Chatterjee, A. and Lahiri, S. N. (2011). Bootstrapping lasso estimators. J. Amer.
Statist. Assoc., 106(494):608–625.
Chatterjee, A. and Lahiri, S. N. (2013). Rates of convergence of the adaptive LASSO
estimators to the oracle distribution and higher order refinements by the bootstrap.
Ann. Statist., 41(3).
Durrett, R. (2010). Probability: theory and examples. Cambridge University Press.
Efron, B., Hastie, T., Johnstone, I., and Tibshirani, R. (2004). Least Angle Regression. Ann. Statist., 32(2):407–451.
Friedman, J., Hastie, T., and Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9:432–441.
Janková, J. and van de Geer, S. (2015). Confidence intervals for high-dimensional
inverse covariance estimation. Electron. J. Statist., 9:1205–1229.
Javanmard, A. and Montanari, A. (2013). Model selection for high-dimensional regression under the generalized irrepresentability condition. In Burges, C., Bottou,
L., Welling, M., Ghahramani, Z., and Weinberger, K., editors, Advances in Neural
Information Processing Systems 26, pages 3012–3020.
Javanmard, A. and Montanari, A. (2014). Confidence intervals and hypothesis testing
for high-dimensional regression. J. Mach. Learn. Res., 15(1):2869–2909.
Knight, K. and Fu, W. (2000). Asymptotics for lasso-type estimators. Ann. Statist.,
28(5):1356–1378.

18

Lauritzen, S. L. (1996). Graphical Models. Clarendon Press, Oxford.
Li, K. C. (1989). Honest confidence regions for nonparametric regression. Ann.
Statist., 17(3):1001–1008.
Mazumder, R. and Hastie, T. (2012). The Graphical Lasso: New Insights and Alternatives. Electron. J. Statist., page 21252149.
Meinshausen, N. and Bühlmann, P. (2006). High-dimensional graphs and variable
selection with the lasso. Ann. Statist., 34(3):1436–1462.
Ng, B., G. Varoquaux, J.-B. P., and Thirion, B. (2013). A novel sparse group gaussian
graphical model for functional connectivity estimation. Information Processing in
Medical Imaging.
Ravikumar, P., Raskutti, G., Wainwright, M. J., and Yu, B. (2008). High-dimensional
covariance estimation by minimizing l1-penalized log-determinant divergence. Electron. J. Statist., 5:935–980.
Ren, Z., Sun, T., Zhang, C.-H., and Zhou, H. H. (2015). Asymptotic normality
and optimalities in estimation of large gaussian graphical models. Ann. Statist.,
43(3):991–1026.
Rothman, A. J., Bickel, P. J., Levina, E., and Zhu, J. (2008). Sparse permutation
invariant covariance estimation. Electron. J. Statist., 2:494–515.
Sun, T. and Zhang, C.-H. (2012). Sparse matrix inversion with scaled Lasso. J. Mach.
Learn. Res., 14:3385–3418.
Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. R. Stat.
Soc. Ser. B. Stat. Methodol., 58:267–288.
van de Geer, S. (2016). Worst possible sub-directions in high-dimensional models.
Journal of Multivariate Analysis, 146:248–260.
van de Geer, S., Bühlmann, P., Ritov, Y., and Dezeure, R. (2013). On asymptotically
optimal confidence regions and tests for high-dimensional models. Ann. Statist.,
42(3):1166–1202.
van der Vaart, A. (2000). Asymptotic Statistics. Cambridge University Press.
Yuan, M. (2010). High dimensional inverse covariance matrix estimation via linear
programming. J. Mach. Learn. Res., 11:2261–2286.
Yuan, M. and Lin, Y. (2007). Model selection and estimation in the gaussian graphical
model. Biometrika, page 117.
Zhang, C.-H. and Zhang, S. S. (2014). Confidence intervals for low-dimensional parameters in high-dimensional linear models. J. R. Stat. Soc. Ser. B. Stat. Methodol.,
76:217–242.

19

Supplementary material
We summarize some preliminary material in Section 8 and proofs of all the results in
Section 9. In appendices A and B, we summarize some known results.

8

Preliminary results: Rates of convergence of the nodewise Lasso

We provide a brief overview of the results for the nodewise Lasso (derived in van
de Geer et al. (2013)) in the following two Lemmas. Lemma 5 below follows from
classical concentration results for sub-Gaussian random variables (see Bühlmann and
van de Geer (2011)). The proof of Lemma 6 below can be found in van de Geer et al.
(2013).
Lemma 5. Suppose that Xj , j = 1, . . . , p, with EXj = 0 and var(Xj ) = Σ0 satisfies
the eigenvalue condition A1 and the sub-Gaussianity condition A3. Then there exist
constants c1 , c2 , c > 0 such that for any τ sufficiently large
P(Tjc ) ≤ c1 p−τ c2 ,
where
Tj

p
:= {kΣ̂Θ0j − ej k∞ /n ≤ cτ log p/n,
p
kΣ̂ − Σ0 k∞ ≤ cτ log p/n,
p
|ηjT ηj /n − τj2 | ≤ cτ log p/n}.

Moreover, for the set ∩pj=1 Tj we get by the union bound
P((∩pj=1 Tj )c ) = P(∪pj=1 Tjc ) ≤ p max P(Tjc ) ≤ c1 p1−τ c2 .
j=1,...,p

Lemma 6 (a version of Theorem 2.4 in van de Geer et al. (2013)). Let τ > 0 and
suppose q
that Θ̂ is the nodewise Lasso estimator (4) with regularization parameters

λj ≥ cτ logn p uniformly in j, for some sufficiently large constant c > 0. Suppose
that A1, A2, A3 are satisfied. Then there exists a constant Cτ > 0 such that on the
set Tj defined in Lemma 5 it holds
p
p
kΘ̂j − Θ0j k1 ≤ Cτ s log p/n, |τ̂j2 − τj2 | ≤ Cτ s log p/n,
kΘ̂j − Θ0j k22 ≤ Cτ slog p/n.

Furthermore, on the set ∩pj=1 Tj we have
p
max kΘ̂j − Θ0j k1 ≤ Cτ s log p/n,

max |τ̂j2 − τj2 | ≤ Cτ

j=1,...,p

j=1,...,p

max kΘ̂j − Θ0j k22 ≤ Cτ slog p/n.

j=1,...,p

20

p
s log p/n,

9

Proofs for Section 3

Lemma 1. We first derive a bound for |∆ij | (which will be useful later) and then a
bound for k∆k∞ .
Let τ > 0 and consider the set Tj defined in Lemma 5. Note that by Lemma 5
(and since kΘ0i k2 ≤ L, kΘ0j k2 ≤ L) we have P(Ti ∩ Tj ) ≤ c1 p−c2 τ for some constants
c1 , c2 > 0. We now condition on the set Ti ∩ Tj .
By the definition of T̂ , we have
T̂ij − Θ0ij

= Θ̂ij − Θ0ij − Θ̂Ti λj Ẑj
= −(Θ0i )T (Σ̂ − Σ0 )Θ0j
+ (Σ̂Θ0i − ei )T (Θ0j − Θ̂j ) + (Θ0i − Θ̂i )T (Σ̂Θ̂j − ej ) .
{z
}
|
{z
} |
rem2

rem1

For the first remainder, we obtain
|rem1 | = |(Σ̂Θ0i − ei )T (Θ0j − Θ̂j )| ≤ kΣ̂Θ0i − ei k∞ kΘ0j − Θ̂j k1 .
Under A1 and A3, by Lemma 10 in Appendix A, we have
p
max kΣ̂Θ0j − ej k∞ = OP ( log p/n).
j=1,...,p

Consequently, and using `1 -rates of convergence from Lemma 6, we obtain for some
constant Cτ > 0
|rem1 | ≤ Cτ slog p/n.
q
p
2
2
and since 1/τj2 = O(1), we have 1/τ̂j2 =
By Lemma 6, we have |τ̂j − τj | ≤ Cτ s log
n
O(Cτ ). Hence for the second remainder, by the KKT conditions and Lemma 6, we
obtain
|rem2 | = |(Θ0i − Θ̂i )T (Σ̂Θ̂j − ej )|
≤ kΣ̂Θ̂j − ej k∞ kΘ0i − Θ̂i k1
p
≤ λj /τ̂j2 (1 + kκ̂j k1 )s log p/n = O(Cτ s log p/n).
˜ ij | = O(slog p/n) on the set Ti ∩ Tj .
˜ ij := rem1 + rem2 , we have |∆
Therefore, for ∆
Now condition on the set ∩pj=1 Tj . By Lemma 5, we have for some constants c1 , c2 > 0
it holds
P((∩pj=1 Tj )c ) ≤ c1 p1−c2 τ .
We again have the decomposition
T̂ − Θ0 = −Θ0 (Σ̂ − Σ0 )Θ0 + (Θ0 Σ̂ − I)(Θ0 − Θ̂) + (Θ0 − Θ̂)T (Σ̂Θ̂ − I) .
|
{z
}
˜
∆

21

(16)

But then we obtain
˜ ∞ ≤ k(Θ0 Σ̂ − I)(Θ0 − Θ̂)k∞ + k(Θ0 − Θ̂)T (Σ̂Θ̂ − I)k∞
k∆k
≤ kΘ0 Σ̂ − Ik∞ Θ0 − Θ̂
≤

1

+

Θ0 − Θ̂

max kΣ̂Θ0j − ej k∞ Θ0 − Θ̂

j=1,...,p

+ Θ0 − Θ̂

1

kΣ̂Θ̂ − Ik∞

1

max λj (1 + kκ̂j k1 )/τ̂j2 .

1 j=1,...,p

p
By Lemma 6, we have Θ0 − Θ̂
= O(Cτ s log p/n) and maxj=1,...,p |τ̂j2 − τj2 | =
1
√
O(Cτ maxj=1,...,p sλj ). Hence, and by the last display we get
√
√
˜ ∞ = O(Cτ s log p/ n).
k∆k∞ = nk∆k
Therefore for some constant Cτ > 0 we obtain
√
P(k∆k∞ ≥ Cτ s log p/ n) ≤ P((∩pj=1 Tj )c ) ≤ c1 p1−c2 τ .
Now note that the constants c1 , c2 are universal, therefore we can take the supremum
over G(s) to obtain
√
sup P(k∆k∞ ≥ Cτ s log p/ n) ≤ c1 p1−c2 τ .
Θ0 ∈G(s)

If we choose τ > 0 sufficiently large so that 1 − c2 τ < 0 then
√
lim sup P(k∆k∞ ≥ Cτ s log p/ n) ≤ lim c1 p1−c2 τ = 0.
n→∞ Θ ∈G(s)
0

n→∞

Theorem 1. By (9), it holds
√
n(T̂ij − Θ0ij )/σij = −Z̃/σij + rem/σij ,
where
Z̃ :=

√

n(Θ0i )T (Σ̂ − Σ0 )Θ0j =:

n
√ X
1
((Θ0i )T Xk XkT Θ0j − Θ0ij ),
n
n
k=1

√

−

ei )T (Θ0j

−

Θ̂0j )

√

n(Θ0i − Θ̂0i )T (Σ̂Θ̂0j − ej ).


√ p = oP (1). Note that under A3,
First observe that by Lemma 1, |rem| = OP s log
n
the results of Lemma 1 are uniform in Θ0 .
and rem :=

n(Σ̂Θ0i

+

Denote Zk := (Θ0i )T Xk XkT Θ0j − Θ0ij . To show Z̃/σij
N (0, 1), we apply BerryEsseen theorem (see e.g. Durrett (2010)). By Lemma 8 in Appendix A (note that
kΘ0i k2 ≤ L, kΘ0j k2 ≤ L), we have a moment bound for m ≥ 2,
E|Zk |m /(2σ0 L2 K)m ≤
22

m!
(K/σ0 )m−2 ,
2

for some constant σ0 > 0. Under A1, for every m fixed we have supG̃ E|Zk |m = O(1).
2 = E(Z 2 ) ≥ ω 2 for some universal constant ω > 0. Therefore,
We have σij
1
|PΘ0 (Z̃/σij < z) − Φ(z)| ≤

3E|Z1 |3
C
√ ≤√ ,
n
(E|Z1 |2 )3/2 n

where C does not depend on Θ0 and n, which concludes the proof.

9.1

Proofs for Section 3.1

Lemma 2. Since X1 ∼ N (0, Σ0 ), then Z := Θ0 X ∼ N (0, Θ0 ). It is well known that
E(Zi2 Zj2 ) = Θ0ii Θ0jj + 2(Θ0ij )2 . Consequently, we get
T

2
σij
= var((Θ0i )T X1 (Θ0j ) X1 ) = var((Θ0 X1 )i (Θ0 X1 )j ) = Θ0ii Θ0jj + (Θ0ij )2 .

We next use the results of Lemmas 5 and 6. The set ∩pi=1 Tj from Lemma 5 holds
ˆ ij :=
with probability at least 1 − c1 p1−τ c2 , for some constants c1 , c2 > 0. We define ∆
p
0
Θ̂ij − Θij . Conditioning on the set ∩i=1 Tj we thus get (using Lemma 6 (Euclidean
norm bound) and condition A1)
2
2
max |σ̂ij
− σij
| ≤

i,j=1,...,p

≤

max |Θ̂ii Θ̂jj − Θ0ii Θ0jj | + max |Θ̂2ij − (Θ0ij )2 |

i,j=1,...,p

i,j=1,...,p

ˆ ii ∆
ˆ jj +
max |∆

i,j=1,...,p

ˆ jj
Θ0ii ∆

ˆ ii |
+ Θ0jj ∆

ˆ ij (∆
ˆ ij + 2Θ0ij )|
max |∆
p
≤ Cτ slog p/n,
+

i,j=1,...,p

(17)

for some constant Cτ > 0.
Lemma 3. To simplify notation, we write Θ := Θ0 . We have
)
( n
1X T
2
2
T
2
2
max |σ̂ij − σij | =
max |
(Θ̂i Xk Xk Θ̂j ) − Θ̂ij
i,j=1,...,p
i,j=1,...,p
n
k=1

− E(Θi X1 X1T Θj )2 − Θ2ij |
n
1X T
≤
max |
(Θ̂i Xk XkT Θ̂j )2 − E(Θi X1 X1T Θj )2 |
i,j=1,...,p n
| k=1
{z
}
I

+ max | Θ̂2ij − Θ2ij |.
i,j=1,...,p |
{z }
II

By A1, |Θij | = O(1) and by Lemma 6, we have
max kΘ̂i − Θi k2 = OP (

i,j=1,...,p

23

p
slog p/n).

Hence
p
max |II| ≤ 2 max |Θij ||Θ̂ij − Θij | + |Θ̂ij − Θij |2 = OP ( slog p/n).

i,j=1,...,p

i,j=1,...,p

For the first term, we have
n

1X T
max |I| =
max |
(Θ̂i Xk XkT Θ̂j )2 − E(Θi X1 X1T Θj )2 |
i,j=1,...,p
i,j=1,...,p n
k=1

n
1X T
≤
max |
(Θ̂i Xk XkT Θ̂j )2 − (ΘTi Xk XkT Θj )2 |
i,j=1,...,p n
| k=1
{z
}
rem1

n
1X T
+ max |
(Θi Xk XkT Θj )2 − E(Θi Xk XkT Θj )2 | .
i,j=1,...,p n
k=1
{z
}
|
rem2

By symmetrization,
n

E|rem2 | ≤ 2E max |
i,j=1,...,p

1X T
(Θi Xk XkT Θj )2 k |,
n
k=1

where k , k = 1, . . . , n is a sequence of independent Radechamer random variables,
independent of X. Let K > 0 and consider truncation by K as follows
n

1X T
E|rem2 | ≤ 2E max |
(Θi Xk XkT Θj )2 k 1(ΘT Xk X T Θj )2 ≤K |
i
k
i,j=1,...,p n
k=1
{z
}
|
r1

n
1X T
+2E max |
(Θi Xk XkT Θj )2 k 1(ΘT Xk X T Θj )2 >K |,
i
k
i,j=1,...,p n
k=1
|
{z
}
r2

The term r1 can be bounded using Hoeffding’s inequality since for i, j = 1, . . . , p and
k = 1, . . . , n it holds
|(ΘTi Xk XkT Θj )2 k 1(ΘT Xk X T Θj )2 ≤K | ≤ K.
i

k

Thus for any a > 0
n

P( max |
i,j=1,...,p

≤e

1X T
(Θi Xk XkT Θj )2 k 1(ΘT Xk X T Θj )2 ≤K | ≥ a)
i
k
n

k=1
−na2 /(2K 2 )

.

Hence integrating over a, we get for the expectation
√
Z ∞
2π K
√ .
Er1 =
P(r1 > a)da ≤
2
n
0
24

Now we bound Er2 . Rewriting the expectation and using the union bound, we obtain
Z ∞
P(r2 > x)dx
Er2 =
0
Z ∞
≤ p2 n
max max P(|(ΘTi Xk XkT Θj )2 k 1(ΘT Xk X T Θj )2 >K | > x)dx
i
k
i,j=1,...,p
k=1,...,n
Z0 ∞
≤ p2 n
max max P((ΘTi Xk XkT Θj )2 > max(x, K))dx
0

i,j=1,...,p k=1,...,n

Next by the sub-Gaussianity condition A3, we have

1/2
P (ΘTi Xk XkT Θj )2 > max(x, K) ≤ c1 e−c2 max(x,K) ,
for some constants c1 , c2 > 0. Hence
Z ∞
1/2
2
c1 e−c2 max(x,K) dx
Er2 ≤ p n
0


2
1
1/2
2
−c2 K 1/2
2
1/2
= p nc1 Ke
+ p nc1
K
+
e−c2 K .
c2
c2
We choose K = n(1−)/2 where the  is taken from the assumptions of this lemma.
Then we get
E|rem2 | ≤ Er1 + Er2
√
2π −/2
(1−)/4
n
+ p2 nc1 n(1−)/2 e−c2 n
≤
2


2
1
(1−)/4
2
(1−)/4
+ p nc1
n
+
e−c2 n
.
c2
c2
Then by the assumption log(p ∨ n)/n(1−)/4 = o(1), we obtain that the above bound
converges to zero for n → ∞. This implies that |rem2 | = oP (1).
To bound maxi,j=1,...,p |rem1 |, we denote Wi,k := ΘTi Xk and Ŵi,k := Θ̂Ti Xk for i =
1, . . . , p and k = 1, . . . , n. Then we can rewrite
n

rem1 = |

1X
(Ŵi,k Ŵj,k )2 − (Wi,k Wj,k )2 |.
n
k=1

Observe that
n

rem1 = |

1X
(Ŵi,k Ŵj,k )2 − (Wi,k Wj,k )2 |
n
k=1

n
1X
= |
{(Ŵi,k − Wi,k )(Ŵj,k − Wj,k )
n
k=1

+ (Ŵi,k − Wi,k )Wj,k
+ Wi,k (Ŵj,k − Wj,k ) + Wi,k Wj,k }2
−(Wi,k Wj,k )2 |.
25

To bound |rem1 |, we make use of the last expression where we expand the curly
bracket {·}2 . We have
n

1X
{(Ŵi,k − Wi,k )(Ŵj,k − Wj,k ) + (Ŵi,k − Wi,k )Wj,k
n

rem1 =

k=1

+ Wi,k (Ŵj,k − Wj,k ) + Wi,k Wj,k }2 − (Wi,k Wj,k )2
n
1X
(Ŵi,k − Wi,k )2 (Ŵj,k − Wj,k )2
n

=

k=1

+ 2(Ŵi,k − Wi,k )2 (Ŵj,k − Wj,k )Wj,k
+2(Ŵi,k − Wi,k )(Ŵj,k − Wj,k )2 Wi,k
+ 4(Ŵi,k − Wi,k )(Ŵj,k − Wj,k )Wi,k Wj,k
2
2
+(Ŵi,k − Wi,k )2 Wj,k
+ 2(Ŵi,k − Wi,k )Wj,k
Wi,k
2
2
(Ŵj,k − Wj,k )Wj,k .
(Ŵj,k − Wj,k )2 + 2Wi,k
+Wi,k

By (iteratively) applying the Cauchy-Schwarz (C-S) inequality to eachPsummation
term in the last
display, we can bound each term by a term involving n1 nk=1 (Ŵl,k −
P
4 , l = i, j. Consequently, it suffices to find a rate for
Wl,k )4 and n1 nk=1 Wl,k
n

1X
(Ŵl,k − Wl,k )4
l=1,...,p n
max

k=1

and to show that maxl=1,...,p

Pn

1
n

4
k=1 Wl,k

= OP (1). Then
n

max |rem1 | = OP ( max (

i,j=1,...,p

l=1,...,p

1X
(Ŵl,k − Wl,k )4 )1/4 )).
n
k=1

We now show that
n

√
1X
(Ŵl,k − Wl,k )4 = OP ((s log p/ n)2 ),
l=1,...,p n
max

(18)

k=1

and that

n

1X 4
Wl,k = OP (1).
l=1,...,p n
max

(19)

k=1

The claim (18) follows by Lemma 7 below. We now show (19). When bounding
maxi,j=1,...,p |rem2 | above, we have shown that for all a > 0
n

1X 2
2
2
2
lim P( max |
Wi,k Wj,k
− EWi,k
Wj,k
| ≥ a) = 0.
n→∞
i,j=1,...,p n
k=1

2 W 2 = O(1), it follows that also (note here that the case when
But then since EWi,k
j,k
i = j is covered)
n
1X 2
2
lim P( max
Wi,k Wj,k
≥ C2 ) = 0,
n→∞
i,j=1,...,p n
k=1

26

for some constant C2 > 0.
Collecting all the results above, for any η > 0 we have
2
2
lim P( max |σ̂ij
− σij
| ≥ η) = 0.

n→∞

i,j=1,...,p

Lemma 7. Under the assumptions of Lemma 3 (and using notation of Lemma 3),
on the set ∩pj=1 Tj it holds
n

√
1X
max
(Ŵl,k − Wl,k )4 = O((s log p/ n)2 ).
l=1,...,p n
k=1

Proof. We will first show that
n

1X
(Ŵl,k − Wl,k )2 = O(s log p/n).
l=1,...,p n
max

k=1

To see this, observe that for l ∈ {i, j},
n

1X
max
(Ŵl,k − Wl,k )2 =
l=1,...,p n
k=1

≤

max kX(Θ̂l − Θl )k22 /n,

l=1,...,p

max (Θ̂l − Θl )T Σ0 (Θ̂l − Θl )

l=1,...,p

+ max |(Θ̂l − Θl )T (Σ̂ − Σ0 )(Θ̂l − Θl )|
l=1,...,p

≤ Λmax (Σ0 ) max kΘ̂l − Θl k22
l=1,...,p

+ kΣ̂ − Σ0 k∞ max kΘ̂l − Θl k21 .
l=1,...,p

By the above bound and by Lemma 6, we observe
n

1X
(Ŵl,k − Wl,k )2 = OP (s log p/n),
l=1,...,p n
max

k=1

as required. Rewriting this as
n

√
1 X
max √
(Ŵl,k − Wl,k )2 = O(s log p/ n),
l=1,...,p
n
k=1

and taking square of both sides, we obtain
n

n

1X
max
(Ŵl,k − Wl,k )4 ≤
l=1,...,p n
k=1

1 X
max √
(Ŵl,k − Wl,k )2
l=1,...,p
n
k=1
√
= O((s log p/ n)2 ).

27

!2

9.2

Proofs for Section 3.2

Theorem 2. By the decomposition (16), we have
T̂ − Θ0 = −Θ0 (Σ̂ − Σ0 )Θ0 + (Θ0 Σ̂ − I)(Θ0 − Θ̂) + (Θ0 − Θ̂)T (Σ̂Θ̂ − I) .
{z
}
|
˜
∆

Under A1, A2, A3 by the proof of Lemma 1, there exist c1 , c2 such that with probability at least c1 p1−c2 τ it holds that
˜ ∞ = O(Cτ s log p/n).
k∆k
By Lemma 9 in Appendix A, with probability at least c1 e−c2 τ it holds
√
|(Θ0i )T (Σ̂ − Σ0 )Θ0j | = O(Cτ / n).
Hence there exists a constant Cτ > 0 such that



1
log p
0
P |T̂ij − Θij | ≥ Cτ max √ , s
≤ c3 max{p1−c2 τ , e−c2 τ }.
n
n
But then since the constants c2 , c3 are universal, we can take the supremum over the
model. By taking τ sufficiently large it follows that



1
log p
0
sup P |T̂ij − Θij | ≥ Cτ max √ , s
≤ c4 e−τ .
n
n
Θ0 ∈G(s)
We proceed to show the second result of the Theorem. We have by Lemma 10 in
Appendix A that there exist constants c1 , c2 such that with probability at least c1 p−c2 τ
it holds
p
kΘ0 (Σ̂ − Σ0 )Θ0 k∞ = O(Cτ log p/n).
˜ ∞ + kΘ0 (Σ̂ − Σ0 )Θ0 k∞ . But then there exist constants
We have kT̂ − Θ0 k∞ ≤ k∆k
Cτ , c1 , c2 > 0 such that
(r
)!
log p log p
P kT̂ − Θ0 k∞ ≥ Cτ max
,s
≤ c1 p1−c2 τ .
n
n

9.3

Proofs for Section 3.3

Lemma 4. First note that
T̂ − Θ0 = Ω̂ − Θ0 − Ω̂T (Σ̂Ω̂ − I)
= Ω̂ − Θ0 − ΘT0 (Σ̂Ω̂ − I) − (Ω̂ − Θ0 )T (Σ̂Ω̂ − I)
= −Θ0 (Σ̂ − Σ0 )Θ0 − (Θ0 Σ̂ − I)(Ω̂ − Θ0 ) − (Ω̂ − Θ0 )T (Σ̂Ω̂ − I) .
|
{z
} |
{z
}
rem1

Under A1 and A3 we have
√
max kΣ̂Θ0j − ej k∞ = OP (1/ n).

j=1,...,p

28

rem2

For the first remainder we then have
krem1 k∞ = k(Θ0 Σ̂ − I)(Ω̂ − Θ0 )k∞ ≤ kΘ0 Σ̂ − Ik∞ Θ0 − Ω̂

1

= OP (s log p/n).
For the second remainder we have
krem2 k∞ = k(Ω̂ − Θ0 )T (Σ̂Ω̂ − I)k∞ ≤ kΣ̂Ω̂ − Ik∞ Ω̂ − Θ0

1

= OP (s log p/n).
Collecting the above statements yields for i, j = 1, . . . , p
T̂ij − Θ0ij

and hence we have

A

= (Θ0i )T (Σ̂ − Σ0 )Θ0j + OP (s log p/n)
√
= (Θ0i )T (Σ̂ − Σ0 )Θ0j + oP (1/ n).
√

n(T̂ij − Θ0ij )/σij

N (0, 1).

Concentration results for sub-Gaussian design

Lemma 8. Let α, β ∈ Rp such that kαk2 ≤ M, kβk2 ≤ M. Let Xk ∈ Rp satisfy the
sub-Gaussianity assumption A3 with a constant K > 0. Then for m ≥ 2,
E|αT Xk XkT β − EαT Xk XkT β|m /(2M 2 K 2 )m ≤

m!
.
2

Lemma 9. Let α, β ∈ Rp such that kαk2 ≤ M, kβk2 ≤ M. Let Xk ∈ Rp satisfy the
sub-Gaussianity assumption A3 with a constant K > 0. For all t > 0

√ 
P |αT Σ̂β − αT Σ0 β|/(2M 2 K 2 ) > t + 2t ≤ 2e−nt .
Lemma 10. Assume kαi k2 ≤ M, kβk2 ≤ M for all i = 1, . . . , p and A3 with K. For
all t > 0 it holds
!
r
√
2 log(2p) log(2p)
T
2 2
P max |αi (Σ̂ − Σ0 )β|/(2M K ) > t + 2t +
+
≤ e−nt .
i=1,...,p
n
n
For the proofs of Lemmas 8, 9, 10, see Janková and van de Geer (2015).

29

B

Rates of convergence of the nodewise Lasso

Lemma 11. Let Θ̂j be obtained as in (4). Then it holds Σ̂Θ̂j − ej − λj Ẑj = 0.
Lemma 11. The KKT conditions (5) give
− XT−j (Xj − X−j γ̂j )/n + λj κ̂j = 0.

(20)

Multiplying (20) by γ̂jT and since γ̂jT κ̂j = kγ̂j k1 , we obtain
τ̂j2 = XjT (Xj − X−j γ̂j )/n.

(21)

Rewriting Xj − X−j γ̂j = XΓ̂j = XΘ̂j τ̂j2 , from (21) we obtain
XjT XΘ̂j /n = 1.

(22)

Similarly rewriting (20) gives
XT−j XΘ̂j /n =

λj
κj .
τ̂j2

(23)

Combining (22) and (23) we obtain XT XΘ̂j /n − ej = λj Ẑj , and rearranging, we get
Σ̂Θ̂j − ej − λj Ẑj = 0 as required.

30

