arXiv:1801.08512v1 [math.ST] 25 Jan 2018

Inference in high-dimensional graphical
models
Jana Janková

Sara van de Geer

Seminar for Statistics
ETH Zürich

Abstract
We provide a selected overview of methodology and theory for estimation and inference on the edge weights in high-dimensional directed and
undirected Gaussian graphical models. For undirected graphical models, two main explicit constructions are provided: one based on a global
method that maximizes the joint likelihood (the graphical Lasso) and one
based on a local (nodewise) method that sequentially applies the Lasso
to estimate the neighbourhood of each node. The estimators lead to confidence intervals for edge weights and recovery of the edge structure. We
evaluate their empirical performance in an extensive simulation study.
The theoretical guarantees for the methods are achieved under a sparsity condition relative to the sample size and regularity conditions. For
directed acyclic graphs, we apply similar ideas to construct confidence
intervals for edge weights, when the directed acyclic graph is identifiable.

1
1.1

Undirected graphical models
Introduction

Undirected graphical models, also known as Markov random fields, have become
a popular tool for representing network structure of high-dimensional data in a
large variety of areas including genetics, brain network analysis, social networks
and climate studies. Let G = (V, E) be an undirected graph with a vertex
set V = {1, 2, . . . , p} and an edge set E ⊂ V × V. Let X 0 = (X1 , X2 , ..., Xp )
be a random vector indexed by graph’s vertices. The joint distribution of X 0
belongs to the graphical model determined by G if Xj and Xk are conditionally
independent given all other variables whenever j and k are not adjacent in G.
The graph then encodes conditional independence structure among the entries
of X 0 .
If we assume that X 0 is normally-distributed with a covariance matrix Σ0 ,
one can show that the edge structure of the graph is encoded by the precision
0
matrix Θ0 := Σ−1
0 (assumed to exist). If Θij denotes the (i, j)-th entry of the
0
matrix Θ0 , it is well known that Θij = 0 ⇔ (i, j) 6∈ E. The non-zero entries in
1

the precision matrix correspond to edges in the associated graphical model and
the absolute values of these entries correspond to edge weights.
Therefore to estimate the edge structure of a Gaussian graphical model, we
consider the problem of estimating the precision matrix, based on a sample
of n independent instances X 1 , . . . , X n , distributed as X 0 . We are not only
interested in point estimation, but in quantifying the uncertainty of estimation
such as constructing confidence intervals and tests for edge weights. Confidence
intervals and tests can be used for identifying significant variables or testing
whether networks corresponding to different populations are identical.
The challenge arises due to the high-dimensional regime where the number of
unknown parameters can be much larger than the number of observations n. It
is instructive to first inspect the low-dimensional setting. In the regime when p is
fixed and the observations are normally distributed with EX i = 0, i = 1, . . . , n,
the sample covariance matrix Σ̂ := X T X/n (where X is the n × p matrix of
observations X 1 , . . . , X n ) is the maximum likelihood estimator of the covariance
matrix. The inverse of the sample covariance matrix Θ̂ := Σ̂−1 is the maximum
likelihood estimator of the precision matrix. Asymptotic linearity of Θ̂ follows
by the decomposition
Θ̂ − Θ0 = −Θ0 (Σ̂ − Σ0 )Θ0 + rem0 ,

(1)

where rem0 := −Θ0 (Σ̂ − Σ0 )(Θ̂ − Θ0 ) is the remainder term. The term rem0
can be bounded by Hölder’s inequality to obtain
krem0 k∞ ≤ kΘ0 (Σ̂ − Σ0 )k∞ Θ̂ − Θ0

,
1

where we used the notation kAk∞ = max
Pp1≤i,j≤p |Aij | for the supremum norm
of a matrix A and |||A|||1 := max1≤j≤p i=1 |Aij | for the `1 -operator norm. If
the fourth moments of X i ’s are bounded,
√ the decomposition (1) implies rates
of convergence kΘ̂ − Θ0 k∞ = OP (1/ n), where OP (1) denotes boundedness
√
in probability. The remainder term then satisfies krem0 k∞ = oP (1/ n), oP (1)
denoting convergence in probability to zero. Therefore, Θ̂ is indeed an asymptotically linear estimator of Θ0 and in this sense, we can say it is asymptotically
unbiased. Moreover, Θ̂ is asymptotically normal with a limiting normal distribution.
In high-dimensional settings, the sample covariance matrix does not perform well (see Johnstone (2001) and Johnstone and Lu (2009)) and if p > n,
it is singular with probability one. Various methods have been proposed that
try to reduce the variance of the sample covariance matrix at the price of introducing some bias. The idea of banding or thresholding the sample covariance
matrix was studied in Bickel and Levina (2008b), Bickel and Levina (2008a)
and El Karoui (2008). Methods inducing sparsity through Lasso regularization
were studied by another stream of works. These can be divided into two categories: global methods and local (nodewise) methods. Global methods estimate
the precision matrix typically via a regularized log-likelihood, while nodewise
methods split the problem into a series of linear regressions by estimating neighbourhood of each node in the underlying graph. A popular global method is the
2

`1 -penalized maximum likelihood estimator, known as the graphical Lasso. Its
theoretical properties were studied in a number of papers, including Yuan and
Lin (2007), Friedman, Hastie, and Tibshirani (2008), Rothman, Bickel, Levina,
and Zhu (2008) and Ravikumar, Raskutti, Wainwright, and Yu (2008). The
local approach on estimation of precision matrices in particular includes the
regression approach Meinshausen and Bühlmann (2006),Yuan (2010),Cai, Liu,
and Luo (2011),Sun and Zhang (2012) which uses a Lasso-type algorithm or
Dantzig selector (Candes and Tao, 2007) to estimate each column or a smaller
part of the precision matrix individually.
Inference for parameters in high-dimensional undirected graphical models
was studied in several papers. Multiple testing for conditional dependence in
Gaussian graphical models with asymptotic control of false discovery rates was
considered in Liu et al. (2013). The work Wasserman, Kolar, Rinaldo, et al.
(2014) proposes methodology for inference about edge weights based on BerryEsseen bounds and the bootstrap for certain special classes of high-dimensional
graphs. Another line of work (Ren, Sun, Zhang, Zhou, et al. (2015), Janková
and van de Geer (2014) and Janková and van de Geer (2016b)) proposes asymptotically normal estimators for edge weights in Gaussian graphical models based
on different modifications of initial Lasso-regularized estimators. In particular,
the paper Ren et al. (2015) proposes nodewise regression where each pair of variables, (Xi , Xj ), is regressed on all the remaining variables; this yields estimates
for the parameters of the joint conditional distribution of (Xi , Xj ) given all the
other variables. The papers Janková and van de Geer (2014) and Janková and
van de Geer (2016b) propose methodology inspired by the de-biasing approach in
high-dimensional linear regression that was studied in Zhang and Zhang (2014),
van de Geer, Bühlmann, Ritov, and Dezeure (2014) and Javanmard and Montanari (2014). This chapter discusses and unifies the ideas from the papers
Janková and van de Geer (2014) and Janková and van de Geer (2016b).
A different approach to structure learning in undirected graphical models is
the Hyvärinen score matching (see e.g. Drton and Maathuis for a discussion
of this approach). Methodology for asymptotically normal estimation of edge
parameters in pairwise (not necessarily Gaussian) graphical models based on
Hyvärinen scoring was proposed in Yu, Kolar, and Gupta (2016).

1.2

De-biasing regularized estimators

The idea of using regularized estimators for construction of asymptotically normal estimators is based on removing the bias associated
Pn with the penalty. Consider a real-valued loss function ρΘ and let Rn (Θ) := i=1 ρΘ (X i )/n denote the
average risk, given an independent sample X 1 , . . . , X n . Under differentiability
conditions, a regularized M-estimator Θ̂ based on the risk function Rn can often
be characterized by estimating equations
Ṙn (Θ̂) + ξ(Θ̂) = 0,

(2)

where Ṙn is the gradient of Rn and ξ(Θ̂) is a (sub-)gradient corresponding
to the regularization term, evaluated at Θ̂. The idea is to improve on the
3

initial estimator by finding a root Θ̂de-bias closer to the solution of estimating
equations without the bias term ξ(Θ̂), i.e. a new estimator Θ̂de-bias such that
Ṙn (Θ̂de-bias ) ≈ 0. A natural way is to define a corrected estimator Θ̂de-bias from
a linear approximation to Ṙn
Ṙn (Θ̂) + R̈n (Θ̂)(Θ̂de-bias − Θ̂) = 0.

(3)

In high-dimensional settings, the matrix R̈n (Θ̂) is typically rank deficient and
thus not invertible. Suppose that we have an approximate inverse denoted by
R̈n (Θ̂)inv . Then we can approximately solve (3) for Θ̂de-bias to obtain
Θ̂de-bias ≈ Θ̂ − R̈n (Θ̂)inv Ṙn (Θ̂),

(4)

provided that the remainder term is small. We refer to the step (4) as the debiasing step since the correction term is proportional to the bias term. Generally
speaking, if the initial estimator Θ̂ and the approximate inverse of R̈n (Θ̂) are
consistent in a strong-enough sense, then the new estimator Θ̂de-bias will be a
consistent estimator of its population version Θ0 per entry at the parametric
rate. The de-biasing step (4) may be viewed as one step of the Newton-Raphson
scheme for numerical optimization.
In consecutive sections, we will look in detail at the bias of several particular
examples of regularized estimators, including the graphical Lasso (Yuan and
Lin (2007)) and nodewise Lasso (Meinshausen and Bühlmann (2006)). We now
provide a unified de-biasing scheme which covers both special cases treated
below (see also van de Geer (2016), Chapter 14). Suppose that a (possibly nonsymmetric) estimator Θ̂ is available which is an approximate inverse of Σ̂ in the
sense that the following condition is satisfied
Σ̂Θ̂ − I + η(Θ̂) = 0,

(5)

where η(Θ̂) is a bias term. This condition in some sense corresponds to the
estimating equations (2). We can express Θ̂ from (5) by straightforward algebra
which leads to the decomposition
Θ̂ + Θ̂T η(Θ̂) − Θ0 = −Θ0 (Σ̂ − Σ0 )Θ0 + rem0 + remreg ,

(6)

where
remreg = (Θ̂ − Θ0 )T η(Θ̂).
Compared to the regime with p fixed, there is an additional remainder remreg
corresponding to the bias term. Provided that the remainder terms rem0 and
remreg are small enough, we can take as a new, de-biased estimator, T̂ :=
Θ̂+Θ̂T η(Θ̂). The bias term η(Θ̂) can be expressed from (5) as η(Θ̂) = −(Σ̂Θ̂−I).
Hence we obtain
T̂ = Θ̂ + Θ̂T − Θ̂T Σ̂Θ̂.
(7)
Bounding the remainders rem0 and remreg in high-dimensional settings requires
more refined arguments than when p is fixed. Looking at the remainder remreg ,
we can again invoke Hölder’s inequality to obtain
kremreg k∞ = k(Θ̂ − Θ0 )T η(Θ̂)k∞ ≤
4

Θ̂ − Θ0

1

kη(Θ̂)k∞ .

Thus it suffices to control the rates of convergence of Θ̂ in |||·|||1 -norm and control
the absolute size of entries of the bias term.
√
Provided that the remainders are of small order 1/ n in probability, asymptotic normality per elements of T̂ is a consequence of asymptotic linearity and
can be established under tail conditions on X i ’s, by applying the Lindeberg’s
central limit theorem.

1.3

Graphical Lasso

If the observations are independent N (0, Σ0 )-distributed, the log-likelihood function is proportional to
`(Θ) := tr(Σ̂Θ) − log det(Θ).
The graphical Lasso (see Yuan and Lin (2007), d’Aspremont, Banerjee, and
El Ghaoui (2008), Friedman et al. (2008)) is based on the Gaussian log-likelihood function but regularizes it via an `1 -norm penalty on the off-diagonal
elements of the precision matrix. The diagonal elements of the precision matrix
correspond to certain partial variances and thus should not be penalized. The
graphical Lasso is defined by
Θ̂ = argminΘ=ΘT ,Θ0 tr(Σ̂Θ) − log det(Θ) + λkΘ− k1 ,

(8)

where λ is non-negative tuning parameter and we optimize over symmetric positive definite matrices, denoted by . Here Θ− represents the matrix obtained
by setting the diagonal elements of Θ to zero and kΘ− k1 is the `1 -norm of the
vectorized version of Θ− . The usefulness of the graphical Lasso is not limited
only to Gaussian settings; the theoretical results below show that it performs
well as an estimator of the precision matrix in a large class of non-Gaussian
settings.
A disadvantage of the graphical Lasso (8) is that the penalization does not
take into account that the variables have in general a different scaling. To take
these differences in the variances into account, we define a modified graphical
Lasso with a weighted penalty. To this end, let Ŵ 2 := diag(Σ̂) be the diagonal
matrix obtained from the diagonal of Σ̂. We let
X
Θ̂w = argminΘ=ΘT ,Θ0 tr(Σ̂Θ) − log det(Θ) +
Ŵii Ŵjj |Θij |.
(9)
i6=j

The weighted graphical Lasso Θ̂w is related to a graphical Lasso based on the
sample correlation matrix R̂ := Ŵ −1 Σ̂Ŵ −1 . To clarify the connection, we define
Θ̂norm = argminΘ=ΘT ,Θ0 tr(R̂Θ) − log det(Θ) + kΘ− k1 .

(10)

Then it holds that Θ̂w = Ŵ −1 Θ̂norm Ŵ −1 . The estimator Θ̂norm is of independent interest, if the parameter of interest is the inverse correlation matrix rather
than the precision matrix. The estimators Θ̂w and Θ̂norm are also useful from a
theoretical perspective as will be shown in the sequel.
5

We now apply the de-biasing ideas of Section 1.2 to the graphical Lasso
estimators defined above, demonstrating the procedure on Θ̂. By definition, the
graphical Lasso is invertible, and the Karush-Kuhn-Tucker (KKT) conditions
yield
Σ̂ − Θ̂−1 + λẐ = 0,
where
Ẑij = sign(Θ̂ij ) if Θ̂ij 6= 0,

and

kẐk∞ ≤ 1.

Multiplying by Θ̂, we obtain
Σ̂Θ̂ − I + λẐ Θ̂ = 0.
In line with Section 1.2 above, this implies the decomposition
Θ̂ + Θ̂T η(Θ̂) − Θ0 = −Θ0 (Σ̂ − Σ0 )Θ0 + rem0 + remreg ,
with the bias term η(Θ̂) = λẐ Θ̂. To control the remainder terms rem0 and
remreg , we need bounds for the `1 -error of Θ̂ and to control the bias term, it is
sufficient to control the upper bound kη(Θ̂)k∞ = kλẐ Θ̂k∞ ≤ λ Θ̂

.
1

Oracle bounds
Oracle results for the graphical Lasso were studied in Rothman et al. (2008)
under sparsity conditions and mild regularity conditions. In Ravikumar et al.
(2008), stronger results were derived under stronger regularity conditions (and
weaker sparsity conditions). Here we revisit these results and provide several
extensions.
We summarize the main theoretical conditions which require boundedness
of the eigenvalues of the true precision matrix and certain tail conditions.
Condition A1 (Bounded spectrum). The precision matrix Θ0 := Σ−1
exists
0
and there exists a universal constant L ≥ 1 such that
1/L ≤ Λmin (Θ0 ) ≤ Λmax (Θ0 ) ≤ L.
Condition A2 (Sub-Gaussianity). The observations X i , i = 1, . . . , n, are uniformly sub-Gaussian vectors, i.e. there exists a universal constant K > 0 such
that for every α ∈ Rp , kαk2 = 1 it holds

E exp |αT X i |2 /K 2 ≤ 2 (i = 1, . . . , n).
Under Condition A2, the Bernstein inequality implies concentration results
for Σ̂ as formulated in Lemma 1 below. The proof is omitted and may be found
in Bühlmann and van de Geer (2011) (Lemma 14.13). We denote the Euclidean
norm by k · k2 and the i-th column of a matrix A by Ai .

6

Lemma 1. Assume Condition A2 and that for non-random matrices A, B ∈
Rp×p it holds that kAi k2 ≤ M and kBi k2 ≤ M for all i = 1, . . . , p. Then for all
t > 0, with probability at least 1 − e−nt it holds that
r
√
2 log(2p2 ) log(2p2 )
T
2 2
kA (Σ̂ − Σ0 )Bk∞ /(2M K ) ≤ t + 2t +
+
.
n
n
To derive oracle bounds for the graphical Lasso, we rely on certain sparsity
conditions on the entries of the true precision matrix. To this end, we define for
j = 1, . . . , p,
Dj := {(i, j) : Θ0ij 6= 0, i 6= j},

dj := card(Dj ),

d := max |dj |.
j=1,...,p

The quantity dj is then the degree of a node Xj and d corresponds to the
maximum vertex degree in the associated graphical model (excluding vertex
self-loops). Furthermore, we define
S :=

p
[

Dj ,

s :=

j=1

p
X

dj ,

j=1

thus S denotes the overall off-diagonal sparsity pattern and s is the overall
number of edges (excluding self-loops).
The following theorem is an extension of the result for the graphical Lasso
in Rothman et al. (2008) to the `1 -norm. The paper Rothman et al.
P (2008)
derives rates in Frobenius norm k · kF , which is defined as kAk2F := i,j |A2ij |
for a matrix A.
Theorem 1 (Regime p  n). Let Θ̂ be the minimizer defined by (8). Assume
Conditions A1 and A2. Then for λ satisfying 2λ0 ≤ λ ≤ 1/(8LcL ) and 8c2L sλ2 +
8cL pλ20 ≤ λ0 /(2L), on the set kΣ̂ − Σ0 k∞ ≤ λ0 , it holds that
2
2
2
kΘ̂ − Θ0 k2F /cL + λkΘ̂− − Θ−
0 k1 ≤ 8cL sλ + 8cL pλ0 ,

and
Θ̂ − Θ0

1

≤ 16c2L (p + s)λ,

where cL = 8L2 .
The slow rate in the result above arises from the part of the estimation error
tr[(Σ̂ − Σ0 )(Θ̂ − Θ0 )] which is related to the diagonal elements of the precision
matrix. However, proper normalizing removes this part of the estimation error.
The following theorem derives an extension of Rothman et al. (2008) for the
normalized graphical Lasso Θ̂norm . Denote the true correlation matrix by R0
and let K0 := R0−1 denote the inverse correlation matrix.
Theorem 2 (Regime p  n). Assume that Conditions A1 and A2 hold. Then
for λ satisfying 2λ0 ≤ λ ≤ 1/(8L2 ) and 8c2L sλ2 ≤ λ0 /(2L), on the set kR̂ −
R0 k∞ ≤ λ0 it holds for some constant CL > 0 that
−
2
2
kΘ̂norm − K0 k2F + λkΘ̂−
norm − K0 k1 ≤ 8cL sλ ,

7

Θ̂norm − K0

1

≤ 8cL sλ2 + 8c2L sλ.

Θ̂w − Θ0

1

≤ CL sλ,

where cL = 8L2 .
Using the normalized graphical Lasso leads to faster rates in Frobenius norm
and `1 -norm as shown above. The rates for Θ̂w in |||·|||1 -norm can be then
established immediately. To derive a high-probability bound for kR̂ − R0 k∞ , we
may apply Lemma 1 together with Hölder’s inequality to obtain kR̂ − R0 k∞ =
p
p
OP ( log p/n). Hence, Θ̂w − Θ0
= OP (s log p/n).
1

Remark 1. The above result requires a strong condition on the sparsity in
Θ0 , i.e. there can only be very few non-zero coefficients due to the restriction
8c2L sλ2 ≤ λ0 /(2L). This condition guarantees that a margin condition is satisfied.
√ An example of a graph satisfying the condition is a star graph with order
n edges.
Asymptotic normality
Once oracle results in `1 -norm are available, we can easily obtain results on
asymptotic normality of the de-biased estimator 2Θ̂ − Θ̂Σ̂Θ̂ for the graphical
Lasso and its weighted version. We denote the asymptotic variance of the de2
denotes convergence
biased estimator by σij
:= nvar((Θ0i )T Σ̂Θ0j ). The arrow
in distribution and for a matrix A we denote (A)ij its (i, j)-entry.
p
Theorem 3√(Regime p  n). Assume Conditions A1, A2, λ  log p/n and
√
that (p + s) d = o( n/ log p). Then it holds that
2Θ̂ − Θ̂Σ̂Θ̂ − Θ0 = −Θ0 (Σ̂ − Σ0 )Θ0 + rem,
where

(11)



√
√
kremk∞ = OP 24(8L2 )2 L(p + s) d + 1λ2 = oP (1/ n).

Moreover, for i, j = 1, . . . , p,
√
n(2Θ̂ − Θ̂Σ̂Θ̂ − Θ0 )ij /σij

N (0, 1).

The result of Theorem 3 gives us tools to construct approximate confidence
intervals and tests for individual entries of Θ0 . However, we need a consistent
estimator of the asymptotic variance σij . For the Gaussian case, one may take
2
σ̂ij
:= Θ̂ii Θ̂jj + Θ̂2ij . We omit the proof of consistency of σ̂ij and point the reader
to Janková and van de Geer (2016b), where other distributions are treated as
well. Moreover, Theorem 3 implies parametric
p rates of convergence for estimation of individual entries and a rate of order log p/n for the error in supremum
norm. Theorem 3 requires a stronger sparsity condition that the corresponding
oracle-type inequality in Theorem 1. This is to be expected as will be argued
in Section 1.7.
Using the weighted graphical Lasso, the results of Theorem 3 can be established under weaker conditions as shown in the following theorem.
8

√
√
Theorem 4 (Regime
p  n). Assume Conditions A1, A2 and s d = o( n/ log p).
p
Then for λ  log p/n, the asymptotic linearity (11) holds with Θ̂w , where


√
√
kremk∞ = OP 12(8L2 )2 s d + 1λ2 = oP (1/ n).
Moreover, for i, j = 1, . . . , p,

√

n(2Θ̂w + Θ̂w Σ̂Θ̂w − Θ0 )ij /σij

N (0, 1).

If the parameter of interest is instead the inverse correlation matrix, we
formulate a partial result below.
p
Proposition
1
(Regime
p

n).
Assume
Conditions
A1,
A2,
λ

log p/n
√
√
and that s d = o( n/ log p). Then it holds
2Θ̂norm − Θ̂norm R̂Θ̂norm − K0 = −K0 (R̂ − R0 )K0 + rem,


√
√
kremk∞ = OP 12(8L2 )2 Ls d + 1λ2 = oP (1/ n).
To claim asymptotic normality of K0 (R̂ − R0 )K0 per entry would require
extensions of central limit theorems to high-dimensional settings (see Chernozhukov, Chetverikov, and Kato (2014)) and an extension of the δ-method.
We do not study these extensions in the present work. To give a glimpse,
in
√ the regime when p is fixed, by the central limit theorem it follows that
n(Σ̂ − Σ0 )
Np2 (0, C), where C is the asymptotic covariance matrix. Then
we may apply the δ-method with hij (Σ) := (Ki0 )T diag(Σ)−1/2 Σdiag(Σ)−1/2 Kj0
√
to obtain n(hij (Σ̂) − hij (Σ0 ))
N (0, ḣ(Σ0 )T C ḣ(Σ0 )).
Finally, we show that the sparsity conditions in the above results may be
further relaxed under a stronger regularity condition on the true precision matrix. We provide here a simplified version of the result in Janková and van de
Geer (2014) which assumes an irrepresentability condition on the true precision
matrix. Let κΣ0 be the `∞ -operator norm of the true covariance matrix Σ0 , i.e.
κΣ0 = |||Σ0 |||1 . Let H0 be the Hessian of the expected Gaussian log-likelihood
evaluated at Θ0 , i.e. H0 = Σ0 ⊗ Σ0 . When calculating the Hessian matrix, we
treat the precision matrix as non-symmetric; this will allow us to accommodate
non-symmetric estimators as well. For any two subsets T and T 0 of V × V, we
use HT0 T 0 to denote the |T | × |T 0 | matrix with rows and columns of H0 indexed
0 −1
by T and T 0 respectively. Define κH0 = (HSS
)
.
1
Condition A3. (Irrepresentability condition) There exists α ∈ (0, 1] such that
0
0 −1
maxe∈S c kHeS
(HSS
) k1 ≤ 1 − α, where S c is the complement of S.
Condition A3 is an analogy of the irrepresentable condition for variable selection in linear regression (see van de Geer and Bühlmann (2009)). If we
define the zero-mean edge random variables as Y(i,j) := Xi Xj − E(Xi Xj ), then
the matrix H0 corresponds to covariances of the edge variables, in particular
0
0
H(i,j),(k,l)
+ H(j,i),(k,l)
= cov(Y(i,j) , Y(k,l) ). Condition A3 means that we require
that no edge variable Y(j,k) , which is not included in the edge set S, is highly
correlated with variables in the edge set (see Ravikumar et al. (2008)). The
9

parameter α then is a measure of this correlation with the correlation growing
when α → 0. Some examples of matrices satisfying the irrepresentable condition
may be found in Janková and van de Geer (2016b).
Theorem 5 √
(Regime p  n). Assume that Conditions A1, A2 and A3
p are satisfied, d = o( n/ log p), κΣ0 = O(1) and κH0 = O(1). Then for λ  log p/n,
the asymptotic
linearity (11) holds with Θ̂, where kremk∞ = OP (d log p/n) =
√
oP (1/ n). Moreover,
√
n(2Θ̂ − Θ̂Σ̂Θ̂ − Θ0 )ij /σij
N (0, 1).
The proof of Theorem 5 may be found in Janková and van de Geer (2014).
We remark that
√ under the irrepresentability condition, one can show that |Θ̂ij −
Θ0ij | = OP (1/ n) (see Ravikumar et al. (2008)). This means that one could use
the graphical Lasso itself to construct confidence intervals of asymptotically optimal (parametric) size. However, this holds under the strong irrepresentability
condition which is often violated in practice.
Comparing the results obtained for the (weighted) graphical Lasso, we see
that the strongest result was attained√under the irrepresentable condition and
under the sparsity condition d = o( n/ log p). An analogous result has not
yet been obtained for the graphical Lasso without the irrepresentable condition
(under the same sparsity condition). However, without the irrepresentable condition, we showed the
results for the weighted graphical Lasso under the
√ same √
sparsity condition s d = o( n/ log p). In the next section, we will consider a
procedure based on pseudo-likelihood, for which we can derive identical
√ results
under weaker conditions, namely under the sparsity condition d = o( n/ log p)
and under the Conditions A1 and A2.

1.4

Nodewise square-root Lasso

An alternative approach to estimate the precision matrix is based on linear regression. The idea of nodewise Lasso is to estimate each column of the precision
matrix by doing a projection of every column of the design matrix on all the remaining columns. While this is a pseudo-likelihood method, the decoupling into
linear regressions gains more flexibility in estimating the individual scaling levels compared to the graphical Lasso which aims to estimate all the parameters
simultaneously. Moreover, by splitting the problem up into a series of linear
regressions, the computational burden is reduced compared to the graphical
Lasso.
In low-dimensional settings, regressing each column of the design matrix on
all the other columns would simply recover the inverse of the sample covariance
matrix Σ̂. However, due to the high-dimensionality of our setting, the matrix Σ̂ 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 square-root Lasso (Belloni, Chernozhukov,
and Wang (2011)).

10

The theoretical motivation can be understood in greater detail from the
population version of the method. For each j = 1, . . . , p, we define the vector
0
of partial correlations γj0 = {γj,k
, k 6= j} as follows
γj0 := argminγ∈Rp−1 EkXj − X−j γk22 /n,

(12)

and we denote the noise level by τj2 = EkXj − X−j γj0 k22 /n. Then one may show
that the j-th column of Θ0 can be recovered from the partial correlations γj0 and
the noise level τj using the following identity: Θ0j = (−γj,1 , . . . , −γj,j−1 , 1, −γj,j+1 , . . . ,
−γj,p )T /τj2 .
Hence, the idea is to define for each j = 1, . . . , p the estimators of the partial
correlations, γ̂j = {γ̂j,k , k = 1, . . . , p, j 6= k} ∈ Rp−1 using, for instance, the
square-root Lasso with weighted penalty,
γ̂j := argminγ∈Rp−1 kXj − X−j γk2 /n + 2λkŴ−j γk1 ,

(13)

where by A−j we denote the matrix A without its j-th column. We further
define estimators of the noise level
τ̂j2 := kXj − X−j γ̂j k22 /n,

τ̃j2 := τ̂j2 + λτ̂j kγ̂j k1 ,

for j = 1, . . . , p. Finally, we define the nodewise square-root Lasso estimator


1/τ̃12 −γ̃1,2 /τ̃12
...
−γ̃1,p /τ̃12
 −γ̃2,1 /τ̃22
1/τ̃22
...
−γ̃2,p /τ̃22 


Θ̂ := 
(14)

..
..
..
..


.
.
.
.
−γ̃p,1 /τ̃p2

−γ̃p,p−1 /τ̃p2

...

1/τ̃p2

An equivalent way of formulating the definitions above is
(γ̂j , τ̂j ) = argminγ∈Rp−1 ,τ ∈R kXj − X−j γk22 /n/(2τ ) + τ /2 + 2λkŴ−j γk1 . (15)
Alternative versions of the above estimator were considered in the literature.
One may use the Lasso (Tibshirani (1996)) instead of the square-root Lasso (as
in Janková and van de Geer (2016b)) or the Dantzig selector (see van de Geer
(2016)). Furthermore, one may define the nodewise square-root Lasso with τ̂j
in place of τ̃j .
The properties of the column estimator Θ̂j were studied in several papers
(following Meinshausen and Bühlmann (2006)) and it has been shown to enjoy
oracle properties under the Conditions A1, A2 and under the sparsity condition
d = O(n/ log p) (see van de Geer et al. (2014), where a similar version was
considered).
In line with Section 1.2, we consider a de-biased version of the nodewise
square-root Lasso estimator. The KKT conditions for the optimization problem
(13) give
T
−τ̂j X−j
(Xj − X−j γ̂j )/n + λκ̂j = 0,
(16)

11

for j = 1, . . . , p, where κ̂j is an element of the sub-differential of the function
γj 7→ kγj k1 at γ̂j , i.e. for k ∈ {1, . . . , p} \ {j},
κ̂j,k = sign(γ̂j,k ) if γ̂j,k 6= 0,

and

kκ̂j k∞ ≤ 1.

If we define Ẑj to be a p × 1 vector
Ẑj := (κ̂j,1 , . . . , κ̂j,j−1 , 0, κ̂j,j+1 , . . . , κ̂j,p ),
then the KKT conditions may be equivalently stated as follows
Σ̂Θ̂j − ej − λ

τ̂j
Ẑj = 0.
τ̃j2

Let Ẑ be a matrix with columns Ẑj for j = 1, . . . , p, τ̂ be a diagonal matrix
with elements (τ̂1 , . . . , τ̂p ) and similarly τ̃ := diag(τ̃1 , . . . , τ̃p ). As in Section 1.2,
this yields the decomposition (6) with a bias term η(Θ̂) := ẐΛτ̂ τ̃ −2 . The bias
term can then be controlled as
kη(Θ̂)k∞ ≤ λkτ̂ k∞ kτ̃ −2 k∞ kẐk∞ ≤ λ max τ̂j /τ̃j2 .
1≤j≤p

Theorem √
6 (Regime p  n). Suppose that Conditions A1, A2 are satisfied
and
d
=
o(
n/ log p). Let Θ̂node be the estimator defined in (14) and let λ 
p
log p/n. Then it holds
Θ̂ + Θ̂T − Θ̂T Σ̂Θ̂ − Θ0 = −Θ0 (Σ̂ − Σ0 )Θ0 + rem,
√
where kremk∞ = OP (dλ2 ) = oP (1/ n). Moreover,
√
n(Θ̂ + Θ̂T − Θ̂T Σ̂Θ̂ − Θ0 )ij /σij
N (0, 1).
When the parameter of interest is the inverse correlation matrix, we can
use the normalized version of the nodewise square-root Lasso and we obtain an
analogous result.
Proposition
2 (Regime p √n). Suppose that Conditions A1, A2 are satisfied,
p
let λ  log p/n and d = o( n/ log p). Then
Θ̂norm + Θ̂Tnorm − Θ̂Tnorm R̂Θ̂norm − Θ0 = −K0 (R̂ − R0 )K0 + rem,
√
where kremk∞ = oP (1/ n).

1.5

Computational view

For the nodewise square-root Lasso, we need to solve p square-root Lasso regressions, which can be efficiently handled using interior-point methods with polynomial computational time or first-order methods (see Belloni et al. (2011)). Alternatively to nodewise square-root Lasso, the nodewise Lasso studied in Janková
12

and van de Geer (2016b) may be used, which requires selection of a tuning
parameter for each of the p regressions. This can be achieved e.g. by crossvalidation and can be implemented efficiently using parallel methods (Efron,
Hastie, Johnstone, and Tibshirani (2004)). The graphical Lasso presents a more
computationally challenging problem; we refer the reader to e.g. Mazumder and
Hastie (2012). The computation of the de-biased estimator itself only involves
simple matrix addition and multiplication.

1.6

Simulation results

We consider a setting with n independent observations generated from Np (0, Θ−1
0 ),
where the precision matrix Θ0 follows one of the three models:
1. Model 1: Θ0 has two blocks of equal size and each block is a five-diagonal
matrix with elements (1, 0.5, 0.4) and (2, 1, 0.6), respectively.
2. Model 2: Θ0 is a sparse precision matrix generated using the R package GGMselect (using the function simulateGraph() with parameter 0.07).
The matrix was converted to a correlation matrix with the function cov2cor().
3. Model 3: Θ0ij = 0.5|i−j| , i, j = 1, . . . , p.
We consider 6 different methods: the de-biased estimator based on the
(1)
(2)
(3)
(4)

graphical Lasso (glasso)
weighted graphical Lasso (glasso-weigh),
nodewise square-root Lasso (node-sqrt) as defined above,
nodewise square-root Lasso with alternative τ̃ as in Sun and Zhang (2012)
(node-sqrt-tau)
(5) nodewise Lasso as in Janková and van de Geer (2016b) (node)
and we also consider
(6) the maximum likelihood estimator (MLE).
Furthermore, as a benchmark we report the oracle estimator (oracle) which
applies maximum likelihood using the knowledge of true zeros in the precision
matrix. We also report the target coverage and the efficient asymptotic variance
of asymptotically regular estimators (see Janková and van de Geer (2016a))
(perfect).
For the graphical Lasso (1) and the weighted graphical Lasso (2) we choose
the tuning parameter by maximizing the likelihood on a validation data set (a
new data of the size n). For methods (3), (4) and (5), the universal choice
p
log p/n is used.
We display results on confidence intervals for nominal coverage 95% and for
normally distributed observations in Tables 1, ??, 2 and 3. For other nominal
coverages, we obtain similar performance (these results are not reported). For
other than Gaussian distributions, we refer the reader to the simulation results in
13

Figure 1: Visualization of graphical models used in simulations. Models 1,2 and
3 from left to right. For Model 3, we only plot edges with a weight greater than
0.1.

Janková and van de Geer (2016b). Firstly, the results of the simulations suggest
that the de-biased estimators perform significantly better than the maximum
likelihood estimator even though p < n and secondly, the nodewise methods
seem to outperform the graphical lasso methods in our settings.

1.7

Discussion

We have shown several constructions of asymptotically linear estimators of the
precision matrix (and inverse correlation matrix) based on regularized estimators, which immediately lead to inference in Gaussian graphical models. Efficient algorithms are available for both methods as discussed in Section 1.5. The
constructed estimators achievep
entrywise estimation at the parametric rate and
a rate of convergence of order log p/n in supremum norm.
To provide a brief comparison of the two methods analyzed above, both theoretical and computational results seem in favor of the de-sparsified nodewise
Lasso. Theoretical results for nodewise Lasso in the
√ regime p  n only need
the mild conditions Conditions A1, A2 and d = o( n/ log p) and are uniform
over the considered model. Moreover, the de-sparsified nodewise Lasso may be
thresholded again to yield recovery of the set S with no false positives, and under a beta-min type condition, exact recovery of the set S, asymptotically, with
high probability. The graphical Lasso requires that we impose the strong irrepresentability condition in the high-dimensional regime. However, the graphical
Lasso might be preferred on the grounds that it does not decouple the likelihood.
Moreover, the graphical Lasso estimator is always strictly positive definite and
thus yields an estimator of the covariance matrix as well. The invertibility of
the nodewise Lasso has not yet been explored. √
We remark that the sparsity condition d = o( n/ log p) implied by our analysis is stronger than the condition needed for oracle inequalities and recovery,
namely d = o(n/ log p). However, one can show that this sparsity condition
is essentially necessary for asymptotically normal estimation. This follows by
inspection of the minimax rates (see Ren et al. (2015)).

14

Table 1: Average coverages and lengths of confidence intervals over the active
set S0 := S ∪ {1, . . . , p} and its complement S0c , over 100 realizations. The
average value of the tuning parameters is reported in the last column. The
benchmark “estimators” are labeled by a star. The significance level is 0.05.
Model 1: Block 1: (1, 0.5, 0.4), Block 2: (2, 1, 0.6)
p = 100, n = 200

1
2
3
4
5
6
7
8

Method
glasso
glasso-weigh
node-sqrt
node-sqrt-tau
node
MLE
oracle*
perfect*

Coverage
S0
S0c
77.19 98.07
35.02 98.65
89.92 94.02
83.48 97.40
90.58 96.77
20.92 84.27
94.95
95.00 95.00

Length
S0
S0c
0.36 0.32
0.31 0.27
0.48 0.42
0.38 0.33
0.41 0.35
0.97 0.81
0.49 0.40
0.48 0.40

Average λ
0.088
0.088
0.152
0.152
0.152
-

p = 100, n = 400

1
2
3
4
5
6
7
8

Method
glasso
glasso-weigh
node-sqrt
node-sqrt-tau
node
MLE
oracle*
perfect*

Coverage
S0
S0c
84.28 97.53
46.22 98.41
91.57 94.40
87.11 97.13
91.48 96.40
41.42 91.29
94.87
95.00 95.00

15

Length
S0
S0c
0.27 0.23
0.24 0.20
0.34 0.29
0.28 0.24
0.30 0.25
0.41 0.38
0.34 0.29
0.34 0.29

Average λ
0.049
0.049
0.107
0.107
0.107
-

Table 2: Average coverages and lengths of confidence intervals over the active
set S0 := S ∪ {1, . . . , p} and its complement S0c , over 100 realizations. The
average value of the tuning parameters is reported in the last column. The
benchmark “estimators” are labeled by a star. The significance level is 0.05.
Model 2
p = 100, n = 400

1
2
3
4
5
6
7
8

Method
glasso
glasso-weigh
node-sqrt
node-sqrt-tau
node
MLE
oracle*
perfect*

Coverage
S0
S0c
64.17 98.65
16.80 98.56
87.23 94.43
89.81 97.23
38.19 99.07
50.98 91.22
98.51
95.00 95.00

Length
S0
S0c
0.16 0.15
0.05 0.05
0.24 0.21
0.20 0.18
0.10 0.10
0.30 0.26
0.23 0.20
0.22 0.20

Average λ
0.067
0.067
0.107
0.107
0.107
-

Table 3: Average coverage and length of confidence intervals over all the entries
and an average value of the tuning parameter λ. The significance level is 0.05.
Model 3: Θij = 0.5i−j
p = 100, n = 200
1
2
3
4
5
6
7

Method
glasso
glasso-weigh
node-sqrt
node-sqrt-tau
node
MLE
perfect*

Coverage
90.43
75.81
93.36
92.91
89.88
80.41
95.00

16

Length
0.19
0.33
0.28
0.22
0.20
0.56
0.28

Average λ
0.138
0.138
0.152
0.152
0.152
-

2

Directed acyclic graphs

In this section, we use the de-biasing ideas to construct confidence intervals
for edge weights in directed acyclic graphs (abbreviated as DAGs). A directed
acyclic graph is a directed graph (we distinguish between edges (j, k) and (k, j))
without directed cycles. We consider the Gaussian DAG model, where the DAG
represents the probability distribution of a random vector (X1 , . . . , Xp ) with a
Gaussian distribution N (0, Σ0 ), where Σ0 ∈ Rp×p is an unknown covariance
matrix. A Gaussian DAG may be represented by the linear structural equations
model
X
0
Xj =
βk,j
Xk + j ,
j = 1, . . . , p,
k ∈ pa(j)

where 1 , . . . , p are independent and j ∼ N (0, (ωj0 )2 ). The set pa(j) is called
the set of parents of a node j and it contains all nodes k ∈ {1, . . . , p} such that
there exists a directed edge k → j.
0
Our aim is to construct confidence intervals for edge weights βk,j
. However,
0
without further conditions, the DAG and the βk,j ’s may not be identifiable
from the structural equations model. To ensure identifiability, we assume that
the error variances are equal: ωj0 = ω0 for all j = 1, . . . , p. We remark that
one might equivalently assume that the error variances are all known up to a
multiplicative constant. In this setting, the DAG is identifiable as shown in
Peters and Bühlmann (2014). Our strategy is to use a two-step procedure: in
the first step we use the estimator proposed in van de Geer and Bühlmann
(2013) to estimate the ordering of the variables and in the second step, we use
a de-biased version of nodewise regression to construct the confidence intervals.
Given an n × p matrix X = [X1 , . . . , Xp ], with rows being n independent
observations from the structural equations model, one may rewrite the above
model in a matrix form
X = XB0 + E,
0
0
where B0 := (βk,j
) is a p × p matrix with βj,j
= 0 for all j, and E is an
n × p matrix of noise vectors E := (1 , . . . , p ) with columns j independent of
0
Xk whenever βk,j
6= 0. The rows of E are independent N (0, ω02 I)-distributed
random vectors. The model then implies that X has covariance matrix

Σ0 = ω02 ((I − B0 )−1 )T (I − B0 )−1 .
We define the precision matrix (assumed to exist) by Θ0 := Σ−1
0 . Notice that
Θ0 =

1
(I − B0 )(I − B0 )T .
ω02

We further consider the class of precision matrices corresponding to DAGs. That
is, we let
1
Θ := Θ(B, ω) = 2 (I − B)(I − B)T ,
ω
where (B, ω) is such that there exists a DAG representing the distribution
N (0, Σ) with Σ = ω 2 ((I − B)−1 )T (I − B)−1 . This means that ω > 0 and
17

B can be written as a lower-diagonal matrix, up to permutation of rows. Further we let sB denote the number of nonzero entries in B, which corresponds
to the number of edges in the DAG. Moreover, we denote by B the set of all
edge weights B of DAGs with parameters (B, ω) which have at most αn/ log p
incoming edges (parents) at each node, where α > 0 is given.

2.1

Maximum likelihood estimator with `0 -penalization

In the first step, we use an `0 -penalized maximum likelihood estimator to estimate the DAG. Let Σ̂ = X T X/n be the Gram matrix based on the design matrix
X. The minus log-likelihood is proportional to `(Θ) = trace(ΘΣ̂) − log det(Θ).
Consider the penalized maximum likelihood estimator proposed in van de Geer
and Bühlmann (2013),
(B̂, ω̂)

:=

argminB,ω { `(B, ω) + λ2 sB : Θ = Θ(B, ω), for some DAG
with parameters (B, ω) where B ∈ B},

(17)

where λ ≥ 0 is a tuning parameter. The estimator is denoted by Θ̂ = Θ(B̂, ω̂)
and it has ŝ := sB̂ edges. Calculating the `0 -penalized maximum likelihood
estimator over the class of DAGs is a computationally intensive task, especially
because it involves a search through a class of DAGs under a non-convex constraint of acyclicity of the graph and due to the `0 -penalty. For large scale problems, greedy algorithms may be used, see e.g. Chickering (2002); Hauser and
Bühlmann (2012). The reason for using the `0 -penalty instead of `1 -penalization
in the definition of (17) was discussed in van de Geer and Bühlmann (2013). The
`1 -penalty leads to an objective function which is not constant over equivalent
DAGs encoding the same distribution. The `0 -penalization leads to invariant
scores over equivalent DAGs. The theoretical properties of Θ̂ were studied in
van de Geer and Bühlmann (2013) under the conditions summarized below. We
remark that the paper van de Geer and Bühlmann (2013) primarily studies the
estimator (17) with unequal variances and shows that the estimator converges
to some member of the Markov equivalence class (cf. Pearl (2016)) of a DAG
with a minimal number of edges, under certain conditions.
To make their result precise, we define some further notions. For any vector
β ∈ Rp , let kXβk := (β T Σ0 β)1/2 . By an ordering of variables we mean any
permutation of the set {1, . . . , p}. For any ordering of the variables, π, we let
B̃(π) be the matrix obtained by doing a Gram-Schmidt orthogonalization of the
columns of X in the ordering given by π, with respect to the norm k·k. Moreover,
let Ω̃0 (π) = (I − B̃0 (π))T Σ0 (I − B̃0 (π)) = diag((ω̃10 (π))2 , . . . , (ω̃p0 (π))2 ). We
restate the conditions assumed in van de Geer and Bühlmann (2013).
Condition B1. There exists a universal constant L ≥ 1 such that
1/L ≤ Λmin (Σ0 ) ≤ Λmax (Σ0 ) ≤ L.

18

Condition B2. There exists a constant ηω > 0 such that for all π such that
Ω̃0 (π) 6= ω02 I it holds
p

1X 0
(|ω̃ (π)|2 − ω02 )2 > 1/ηω .
p i=1 j
Condition B3. There exists a sufficiently small constant α∗ such that p ≤
α∗ n/ log p.
Condition B2 is an “omega-min” condition: it imposes that if one uses the
wrong permutation then the error variances are far enough from being equal.
Under the above conditions, the `0 -penalized maximum likelihood estimator
with high-probability correctly estimates the ordering of the variables as shown
in van de Geer and Bühlmann (2013). Let π0 be an ordering of the variables such
that a Gram-Schmidt orthogonalization of the columns of X in the order given
by π0 with respect to the norm k · k yields B0 . Denote the ordering of variables
estimated by the `0 -penalized maximum likelihood estimator by π̂. Then the
result in van de Geer and Bühlmann (2013) states that under Conditions B1,
B2 and B3, with high-probability it holds that π̂ = π0 .

2.2

Inference for edge weights

Given that we have recovered the true ordering π0 , the problem reduces to estimation of regression coefficients in a nodewise regression model, where each
variable is a function of a known set of its “predecessors”. Therefore to con0
struct asymptotically normal estimators for the βk,j
’s, we may use a nodewise
regression approach.
An estimated ordering π̂ yields estimates b
p(j) of the predecessor sets for
each node j = 1, . . . , p. If we have recovered the true ordering, that is π̂ = π0 ,
the estimated predecessor sets b
p(j) are equal to the true predecessor sets, which
are supersets of the parent sets pa(j) for each j = 1, . . . , p. Consequently, given
the predecessor sets, we may obtain a new estimator for the edge weights by
regressing the j-th variable Xj on all its predecessors. The predecessor sets
p(j) might be as large as p − 1, therefore it is necessary to use regularization.
b
Then we use the de-biasing technique in a similar spirit as in Section 1.2. We
remark that in the initial step, one may use any estimator which guarantees
exact recovery of the ordering π0 .
For any non-empty subset T ⊆ {1, . . . , p}, we denote by XT the n × |T |
matrix formed by taking the columns Xk of X such that k ∈ T. We define the
nodewise regression estimator as proposed in Janková and van de Geer (2016b)
(altenatively, one may use the nodewise square-root Lasso studied in Section
1.4)
β̂j = argmin kXj − Xbp(j) βk22 /n + 2λj kβk1 .
β∈R|bp(j)|

19

(18)

The Karush-Kuhn-Tucker conditions for the above optimization problem give
−XbpT(j) (Xj − Xbp(j) β̂j )/n + λj Ẑj = 0,
where the entries of Ẑj satisfy Ẑk,j = sign(β̂k,j ) if β̂k,j 6= 0, and kẐj k∞ ≤ 1 (β̂k,j
denotes the k-th entry of β̂j ). Similarly as in the case of undirected graphical
models, we can define a de-biased estimator. The Hessian matrix of the risk
function in (18) is given by
Σ̂bp(j) := XbpT(j) Xbp(j) /n.
To find a surrogate inverse for Σ̂bp(j) , we construct Θ̂bp(j) using the nodewise
Lasso with tuning parameters λk,j for k ∈ b
p(j). Using Θ̂bp(j) , we define the
de-biased estimator
b̂j := β̂j + Θ̂bTp(j) XbpT(j) (Xj − Xbp(j) β̂j )/n.

(19)

Theorem 7 below shows that the entries of the de-biased estimator are
asymptotically normal. To formulate the result, we define Θ0p(j) to be the matrix
obtained by taking the rows and columns of Θ0 contained in the true predecessor
set p(j). Denote the k-th column of Θ0p(j) by Θ0p(j),k . To provide asymptotically
0
0
: k ∈ p(j)), we need to imnormal estimators for the parameters βp(j)
= (βk,j
pose a sparsity condition on the sizes of the parent sets, which will be denoted
by dj = |pa(j)|.
Theorem
7 (Regime p ≤ n). Let B̂ be the estimator defined by (17) with
p
λ  log p/n and denote the predecessor sets estimated based on B̂ by p̂(j) for
j = 1, . . . , p. p
Let b̂j be defined in (19) with sufficiently large tuning parameters
p(j). Assume Conditions
λj  λk,j  log p/n, uniformly in j, k, where k ∈ b
B1,
B2
and
B3
are
satisfied
with
1/(|α
|
+
|η
|)
=
O(1)
and assume that dj =
∗
ω
√
o( n/ log p). Then it holds
0
T
b̂j − βp(j)
= (Θ0p(j) )T Xp(j)
j /n + rem,

where

(20)

√
kremk∞ = oP (1/ n).

Furthermore, for every k ∈ p(j),
√
0
n(b̂k,j − βk,j
)/σk,j

N (0, 1),

where the asymptotic variance of the de-sparsified estimator is given by
2
T
σk,j
:= nvar((Θ0p(j),k )T Xp(j)
j ) = ω02 (Θ0p(j) )kk .

The result of Theorem 7 can be used to construct confidence intervals for the
0
edge weights βk,j
. To estimate the asymptotic variance, we may define ω̂j2 :=
2
kXj − Xbp(j) β̂j k22 /n and σ̂k,j
:= ω̂j2 (Θ̂bp(j) )kk . The consistency of this estimator
may be easily checked.
20

3

Conclusion

We have provided a unified approach to construct asymptotically linear and normal estimators of low-dimensional parameters of the precision matrix based on
regularized estimators. These estimators allow us to construct confidence intervals for edge weights in high-dimensional Gaussian graphical models and, under
an identifiability condition, for edge weights in the high-dimensional Gaussian
DAG model.
For Gaussian graphical models, we provided two explicit simple constructions: one based on a global method using the graphical Lasso and the second
based on a local method using nodewise Lasso regressions. Efficient computational methods are available for both methods as discussed in Section 1.5. The
constructed estimators are asymptotically normal per entry, achieving the efficient asymptotic variance from the parametric setting. For a detailed analysis
of semi-parametric efficiency bounds in Gaussian graphical models, we refer to
Janková and van de Geer (2016a). For testing hypothesis about a set of edges,
the usual multiple testing corrections may be used although in practical applications, these might turn out to be too conservative. More efficient methods
for multiple testing in this setting are yet to be developed. While throughout
the presented results we have imposed “exact” sparsity constraints on the underlying parameters, we remark that the results might as well be extended to
models which are only approximately sparse (see e.g. Bühlmann and van de
Geer (2011)).
Our main interest lied in developing methodology for graphical models representing continuous random vectors. However, many applications involve discrete
graphical models, where random variables Xj at each vertex j ∈ V take values
in a discrete space. A popular family of distributions for the binary case where
Xj ∈ {−1, 1} is the Ising model. This model finds applications in statistical
physics, neuroscience or modeling of social networks. The Ising model can be
efficiently estimated via a nodewise method: the individual neighbourhoods can
be estimated with `1 -penalized logistic regression as proposed in Ravikumar,
Wainwright, Lafferty, et al. (2010). Logistic regression falls into the framework
of generalized linear models for which the de-biasing methodology was proposed
in van de Geer et al. (2014). Consequently, one may compute the neighbourhood
estimator via `1 -penalized logistic regression and then compute the de-biased
estimator along the lines of van de Geer et al. (2014).
For directed acyclic graphs, we showed that confidence intervals for edge
weights may be constructed for the Gaussian DAG when it is identifiable and
p ≤ α∗ n/ log p. To this end, we require that the error variances in the structural
equations model are equal, or known up to a multiplicative constant. If the
variance are not equal, the model may not be identifiable and work on inference
in this setting is yet to be developed.

21

4
4.1

Proofs
Proofs for undirected graphical models

Lemma 2. Assume that 1/L ≤ Λmin (Θ0 ) ≤ Λmax (Θ0 ) ≤ L for some constant
L ≥ 1. Let E(∆) := tr[∆Σ0 ] − [log det(∆ + Θ0 ) − log det(Θ0 )]. Then for all ∆
such that k∆kF ≤ 1/(2L), E(∆) is well defined and
E(∆) ≥

1
k∆k2F .
2(L + 1/(2L))2

(21)

Proof of Lemma 2. First we show that E(∆) is well defined for all ∆ such that
k∆kF ≤ 1/(2L). To this end, we need to check that Λmin (Θ0 + ∆)
p ≥ c1 for some
c1 > 0. Denote the spectral norm of a matrix M by kM k := Λmax (M M T ).
We have
Λmin (Θ0 + ∆) = min xT (Θ0 + ∆)x ≥ Λmin (Θ0 ) − k∆kF ≥ 1/(2L),
kxk2 =1

T

where we used that |x ∆x| ≤ k∆kxT x and that k∆k ≤ k∆kF .
A second order Taylor expansion with remainder in integral form yields
log det(∆ + Θ0 ) − log det(Θ0 )
Z 1

= tr(∆Σ0 ) − vec(∆)T
(1 − v)(Θ0 + v∆)−1 ⊗ (Θ0 + v∆)−1 dv vec(∆).
0

Then for all ∆ such that k∆kF ≤ 1/(2L), it holds
Z 1

E(∆) = vec(∆)T
(1 − v)(Θ0 + v∆)−1 ⊗ (Θ0 + v∆)−1 dv vec(∆),
0

where ⊗ denotes the Kronecker product and the remainder in the Taylor expansion is in the integral form. Using the fact that the eigenvalues of Kronecker
product of symmetric matrices is the product of eigenvalues of the factors, it
follows for all ∆ such that k∆kF ≤ 1/(2L)
Z 1

Λmin
(1 − v)(Θ0 + v∆)−1 ⊗ (Θ0 + v∆)−1 dv
0

Z
≥

1

(1 − v)Λ2min ((Θ0 + v∆)−1 )dv

0

≥
≥

1
min Λ2 ((Θ0 + v∆)−1 )
2 0≤v≤1 min
1
min
Λ2 ((Θ0 + ∆)−1 ).
2 ∆:k∆kF ≤1/(2L) min

Next we obtain
−2
Λ2min ((Θ0 + ∆)−1 ) = Λ−2
≥
max (Θ0 + ∆) ≥ (kΘ0 k + k∆k)

1
> 0,
(L + 1/(2L))2

where we used k∆k ≤ k∆kF ≤ 1/(2L). Finally this yields that E(∆) ≥
1
2
2(L+1/(2L))2 k∆kF for all k∆kF ≤ 1/(2L), as required.
22

Proof of Theorems 1 and 2. We will prove both Theorem 1 and Theorem 2 at
the same time, since the proofs only differ slightly. For the proof of Theorem
2, one has to replace Σ̂, Σ0 , Θ̂, Θ0 in the proof below by Γ̂, Γ0 , Θ̂norm , K0 ,
respectively.
Let Θ̃ := αΘ̂ + (1 − α)Θ0 , where α :=

M
M +kΘ̂−Θ0 kF

, for some M > 0 to

be specified later. The definition of Θ̃ implies that kΘ̃ − Θ0 kF ≤ M. By the
convexity of the loss function and by the definition of Θ̂, we have
tr(Θ̃Σ̂) − log det(Θ̃) + λkΘ̃− k1 ≤ tr(Θ0 Σ̂) − log det(Θ0 ) + λkΘ−
0 k1 .

(22)

Denote ∆ = Θ̃ − Θ0 and let
E(∆) := tr(∆Σ0 ) − [log det(∆ + Θ0 ) − log det(Θ0 )].
The inequality (22) implies the basic inequality
E(∆) + λkΘ̃− k1 ≤ −tr[∆(Σ̂ − Σ0 )] + λkΘ−
0 k1 .
On the set {kΣ̂ − Σ0 k∞ ≤ λ0 }, we can bound the empirical process term by
|tr[∆(Σ̂ − Σ0 )]|

≤

+
−
+
+
kΣ̂− − Σ−
0 k∞ k∆ k1 + kΣ̂ − Σ0 k2 k∆ kF

+
≤ λ0 k∆− k1 + kΣ̂+ − Σ+
0 k2 k∆ kF .

(23)

In what follows, we work on the set {kΣ̂ − Σ0 k∞ ≤ λ0 }.
We now choose M such that M ≤ 1/(2L), this then implies kΘ̃−Θ0 kF ≤ 1/(2L).
But then Lemma 2 implies that E(Θ̃ − Θ0 ) is well defined and
E(Θ̃ − Θ0 ) ≥ ckΘ̃ − Θ0 k2F ,

(24)

where one can take c := 1/(8L2 ). Using bounds (23) and (24), we obtain from
the basic inequality
+
+
+
ck∆k2F + λkΘ̃− k1 ≤ λ0 k∆− k1 + λkΘ−
0 k1 + kΣ̂ − Σ0 k2 k∆ kF

By the triangle inequality and taking λ ≥ 2λ0 , we obtain
−
+
+
+
2ck∆k2F + λkΘ̃−
S c k1 ≤ 3λk∆S k1 + 2kΣ̂ − Σ0 k2 k∆ kF

Consequently,
2ck∆k2F + λk∆− k1

≤
≤

+
+
+
4λk∆−
S k1 + 2kΣ̂ − Σ0 k2 k∆ kF
√
+
+
+
4λ sk∆−
S kF + 2kΣ̂ − Σ0 k2 k∆ kF

+ 2
2
+
+ 2
≤ 8sλ2 /c2 + ck∆−
S kF /2 + 8kΣ̂ − Σ0 k2 /c + ck∆ kF /2.
2
Taking M such that λ0 M ≥ 8sλ2 /c2 + 8kΣ̂+ − Σ+
0 k2 /c,
2
ck∆k2F + λk∆− k1 ≤ 8sλ2 /c2 + 8kΣ̂+ − Σ+
0 k2 /c ≤ λ0 M.

23

Taking M ≥ 4λ0 /c,
k∆k2F ≤ λ0 M/c ≤ M 2 /4.
But then k∆kF ≤ M/2. The definition of Θ̃ in turn implies that kΘ̂−Θ0 kF ≤ M,
and we can repeat all the arguments with Θ̂ in place of Θ̃. Repetition of the
arguments leads to the oracle inequality
+ 2
2 2
+
ckΘ̂ − Θ0 k2F + λkΘ̂− − Θ−
0 k1 ≤ 8sλ /c + 8kΣ̂ − Σ0 k2 /c.

Finally we distinguish the case of non-normalized graphical Lasso (based on the
covariance matrix) and the normalized graphical Lasso (based on the correlation
matrix). We have for the case of
a) normalized graphical Lasso: Σ̂+ − Σ+
0 = 0 (recall here that Σ̂ ≡ R̂, Σ0 ≡ R0 )
and the oracle inequality gives
2 2
ckΘ̂ − Θ0 k2F + λkΘ̂− − Θ−
0 k1 ≤ 8sλ /c .

b) non-normalized graphical Lasso: we can bound
√
√
kΣ̂+ − Σ+
pkΣ̂+ − Σ+
pλ0 .
0 k2 ≤
0 k∞ ≤
Hence the oracle inequality gives
2 2
2
ckΘ̂ − Θ0 k2F + λkΘ̂− − Θ−
0 k1 ≤ 8sλ /c + 8pλ0 /c.

To show the second statement of the theorems, we use the above oracle inequalities and the following upper bound
Θ̂ − Θ0

≤
1

≤

0 −
max |Θ̂jj − Θ0jj | + kΘ̂−
j − (Θj ) k1

j=1,...,p

kΘ̂ − Θ0 kF + kΘ̂− − Θ−
0 k1 .

To show the third statement of Theorem 2, we use the upper bound
Θ̂w − Θ0

Ŵ −1 Θ̂norm Ŵ −1 − W0−1 K0 W0−1

=
1

≤

kŴ k2∞

Θ̂norm − K0

1

1

+ kŴ − W0 k∞ |||K0 |||1 kŴ k∞

+kW0 k∞ |||K0 |||1 kŴ − W0 k∞

Proof of Theorem 3. Denote CL := 16(8L2 )2 . Using the results of Theorem 1,
we obtain
kremk∞

≤
≤

k(Θ̂ − Θ0 )T (Σ̂Θ0 − I)k∞ + k(Θ̂ − Θ0 )T λẐ Θ̂k∞
Θ̂ − Θ0

kΣ̂ − Σ0 k∞ |||Θ0 |||1 + Θ̂ − Θ0 λkẐk∞ Θ̂
1
1
√
√
≤ CL (p + s)λ d + 1Λmax (Θ0 )λ0 + 2CL (p + s)λ d + 1Λmax (Θ0 )λ
√
3
≤
CL L(p + s) d + 1λ2 .
2
1

24

p
p
Taking λ  log p/n, Lemma 1 implies kΣ̂ − Σ0 k∞√= OP ( log p/n). Then by
the sparsity condition, we obtain kremk∞ = OP (1/ n). By Conditions A1 and
A2, the random variable (Θ0 (Σ̂ − Σ0 )Θ0 )ij has bounded fourth moments and
asymptotic normality per entry follows by application of Lindeberg’s central
limit theorem for triangular arrays (see Janková and van de Geer (2016b) for
more details).
Proof of Theorem 4. For the remainder,
we obtain similarly as in the proof of
√
Theorem 1, kremk∞ ≤ 32 CL Ls d + 1λ2 . Asymptotic normality follows by analogous arguments.
Proof of Proposition 1. Denote CL := 16(8L2 )2 . Using the results of Theorem
2, we obtain
kremk∞

≤
≤
≤

Θ̂norm − K0 kΓ̂K0 − Ik∞ + Θ̂norm − K0 λkẐk∞ Θ̂norm
1
1
√
√
CL sλ d + 1Λmax (Θ0 )λ0 + CL sλ d + 1Λmax (Θ0 )λ
√
3
CL Ls d + 1λ2 .
2

The sparsity condition implies the result.
Proof of Theorem 6. The proof follows along the same lines as the proof of Theorem 1 in Janková and van de Geer (2016b). The only difference is that here we
consider a weighted Lasso to estimate the partial correlations and the estimator τ̂j is defined slightly differently. But for the weighted Lasso (with weights
bounded away from zero and bounded from above with high probability), oracle
inequalities of the same order can be obtained, see Section 6.9 in Bühlmann and
van de Geer (2011), i.e.
kX−j (γ̂j − γj0 )k22 /n + λkγ̂j − γj0 k1 = OP (dj log p/n).
For the estimator of variance we have
|τ̂j2 − τj2 |

≤

kX−j (γ̂j − γj0 )k22 /n + 2|(Xj − X−j γj0 )T X−j (γ̂j − γj0 )/n|

= kX−j (γ̂j − γj0 )k22 /n + 2k(Xj − X−j γj0 )T X−j /nk∞ kγ̂j − γj0 k1
√
= OP (1/ n).
The rest of the proof follows as in Janková and van de Geer (2016b).

4.2

Proofs for directed acyclic graphs

Proof of Theorem 7. By Theorem 5.1 in van de Geer and Bühlmann (2013), we
have under the conditions of the theorem that π̂ = π0 with high probability.
Then also b
p(j) = p(j) for all j, with high probability. Therefore, the estimated
p(j) in the definitions of β̂j and Θ̂bp(j),k , k ∈ b
b
p(j) (and elsewhere) can be replaced
25

1

by p(j). The nodewise Lasso then yields oracle estimators
p β̂j and Θ̂k , k ∈ p(j)
under the Condition B1 and under the sparsity dj = o( n/ log p) (see Janková
and van de Geer (2016b)). This gives in particular that for all j = 1, . . . , p
max kΘ̂p(j),k − Θ0p(j),k k1 = OP (max dj λj ),
k

k

max kΣ̂p(j) Θ̂p(j),k − ek k∞ = OP (max λj ),
k

k

kβ̂j −

βj0 k1

= OP ( max dj λj ).
j=1,...,p

We can write the decomposition
0
T
b̂k,j − βk,j
= (Θ0p(j),k )T Xp(j)
j /n + remk,j ,

(25)

T
where remk,j = (Θ̂p(j),k −Θ0p(j),k )T Xp(j)
j /n−(Σ̂p(j) Θ̂p(j),k −ek )T (β̂j −βj0 ). First
note that by normality and by the independence of Xp(j) and j (which follows
T
by the independence of j ’s and acyclicity of the graph), it holds kXp(j)
j /nk∞ =
p
OP ( log p/n). By Hölder’s inequality

max |remk,j |
k

≤

T
max kΘ̂p(j),k − Θ0p(j),k k1 kXp(j)
j /nk∞
k

+ max kΘ̂Tp(j),k Σ̂p(j) − ek k∞ kβ̂j − βj0 k1
k
√
= OP (dj log p/n) = oP (1/ n),
√
where we used the sparsity assumption dj = o( √n/ log p). Thus we have shown
that the remainder in (25) is of small order 1/ n. Then applying Lindeberg’s
central limit theorem for triangular arrays and by Conditions A1 and A2,
√
T
(Θ0p(j),k )T Xp(j)
j /(σk,j n)
N (0, 1),
which shows the claim.

References
A. Belloni, V. Chernozhukov, and L. Wang. Square-Root Lasso: Pivotal Recovery of Sparse Signals via Conic Programming. Biometrika, 98(4):791–806,
2011.
P. J. Bickel and E. Levina. Regularized estimation of large covariance matrices.
Annals of Statistics, 36(1):199–227, 2008a.
P. J. Bickel and E. Levina. Covariance regularization by thresholding. Annals
of Statistics, 36(6):2577–2604, 2008b.
P. Bühlmann and S. van de Geer. Statistics for high-dimensional data. Springer,
2011.
26

T. Cai, W. Liu, and X. Luo. A constrained l1 minimization approach to sparse
precision matrix estimation. Journal of the American Statistical Association,
106, 2011.
E. Candes and T. Tao. The dantzig selector: Statistical estimation when
p is much larger than n. Annals of Statistics, 35(6):2313–2351, 12 2007.
doi: 10.1214/009053606000001523. URL http://dx.doi.org/10.1214/
009053606000001523.
V. Chernozhukov, D. Chetverikov, and K. Kato. Central Limit Theorems and
Bootstrap in High Dimensions. ArXiv: 1412.3661, 2014.
D.M. Chickering. Optimal structure identification with greedy search. Journal
of Machine Learning Research, 3:507–554, 2002.
Alexandre d’Aspremont, Onureena Banerjee, and Laurent El Ghaoui. Firstorder methods for sparse covariance selection. SIAM Journal on Matrix Analysis and Applications, 30(1):56–66, 2008. ISSN 0895-4798. doi:
10.1137/060670985. URL http://dx.doi.org/10.1137/060670985.
Mathias Drton and Marloes H Maathuis. Structure learning in graphical modeling. Annual Review of Statistics and Its Application, 4:365–393.
B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least Angle Regression.
Ann. Statist., 32(2):407–451, June 2004.
N. El Karoui. Operator norm consistent estimation of large dimensional sparse
covariance matrices. Annals of Statistics, 36(6):2717–2756, 2008.
J. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covariance estimation
with the graphical lasso. Biostatistics, 9:432–441, 2008.
A. Hauser and P. Bühlmann. Characterization and greedy learning of interventional Markov equivalence classes of directed acyclic graphs. Journal of
Machine Learning Research, 13:2409–2464, 2012.
J. Janková and S. van de Geer. Confidence intervals for high-dimensional inverse
covariance estimation. Electronic Journal of Statistics, 9(1):1205 –1229, 2014.
J. Janková and S. van de Geer. Semi-parametric efficiency bounds for highdimensional models. ArXiv:1601.00815, 2016a.
J. Janková and S. van de Geer. Honest confidence regions and optimality for
high-dimensional precision matrix estimation. TEST, 26(1):143–162, 2016b.
A. Javanmard and A. Montanari. Confidence intervals and hypothesis testing
for high-dimensional regression. Journal of Machine Learning Research, 15
(1):2869–2909, 2014.

27

I. M. Johnstone and A. Y. Lu. On consistency and sparsity for principal components analysis in high dimensions. Journal of the American Statistical
Association, 104(486):682–693, 2009.
Iain M. Johnstone. On the distribution of the largest eigenvalue in principal
components analysis. Annals of Statistics, 29(2):295–327, 04 2001. doi: 10.
1214/aos/1009210544.
Weidong Liu et al. Gaussian graphical model estimation with false discovery
rate control. Annals of Statistics, 41(6):2948–2978, 2013.
R. Mazumder and T. Hastie. The Graphical Lasso: New Insights and Alternatives. Electronic Journal of Statistics, 6:2125–2149, 2012.
N. Meinshausen and P. Bühlmann. High-dimensional graphs and variable
selection with the lasso. Annals of Statistics, 34(3):1436–1462, 06 2006.
doi: 10.1214/009053606000000281. URL http://dx.doi.org/10.1214/
009053606000000281.
J. Pearl. Causality: Models, Reasoning and Inference. Cambridge University
Press, 2016.
J. Peters and P. Bühlmann. Identifiability of Gaussian structural equation models with equal error variances. Biometrika, 101:219–228, 2014.
P. Ravikumar, G. Raskutti, M. J. Wainwright, and B. Yu. High-dimensional
covariance estimation by minimizing l1-penalized log-determinant divergence.
Electronic Journal of Statistics, 5:935–980, 2008.
Pradeep Ravikumar, Martin J Wainwright, John D Lafferty, et al. Highdimensional Ising model selection using 1-regularized logistic regression. Annals of Statistics, 38(3):1287–1319, 2010.
Zhao Ren, Tingni Sun, Cun-Hui Zhang, Harrison H Zhou, et al. Asymptotic
normality and optimalities in estimation of large gaussian graphical models.
Annals of Statistics, 43(3):991–1026, 2015.
A. J. Rothman, P. J. Bickel, E. Levina, and J. Zhu. Sparse permutation invariant
covariance estimation. Electronic Journal of Statistics, 2:494–515, 2008. doi:
10.1214/08-EJS176. URL http://dx.doi.org/10.1214/08-EJS176.
T. Sun and C.-H. Zhang. Sparse matrix inversion with scaled Lasso. Journal of
Machine Learning Research, 14:3385–3418, 2012.
R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of The
Royal Statistical Society: Series B, 58:267–288, 1996.
S. van de Geer. Estimation and Testing under Sparsity: École d’Été de SaintFlour XLV. Springer, 2016.

28

S. van de Geer and P. Bühlmann. On the conditions used to prove oracle
results for the lasso. Electronic Journal of Statistics, 3:1360–1392, 2009. doi:
10.1214/09-EJS506. URL http://dx.doi.org/10.1214/09-EJS506.
S. van de Geer and P. Bühlmann. `0 -penalized maximum likelihood for sparse
directed acyclic graphs. Annals of Statistics, 41(2):536–567, 2013.
S. van de Geer, P. Bühlmann, Y. Ritov, and R. Dezeure. On asymptotically
optimal confidence regions and tests for high-dimensional models. Annals of
Statistics, 42(3):1166–1202, 2014.
Sara van de Geer. Worst possible sub-directions in high-dimensional models.
Journal of Multivariate Analysis, 146:248–260, 2016.
Larry Wasserman, Mladen Kolar, Alessandro Rinaldo, et al. Berry-Esseen
bounds for estimating undirected graphs. Electronic Journal of Statistics,
8(1):1188–1224, 2014.
Ming Yu, Mladen Kolar, and Varun Gupta. Statistical inference for pairwise
graphical models using score matching. In D. D. Lee, M. Sugiyama, U. V.
Luxburg, I. Guyon, and R. Garnett, editors, Advances in Neural Information
Processing Systems 29, pages 2829–2837. Curran Associates, Inc., 2016.
Ming Yuan. High dimensional inverse covariance matrix estimation via linear
programming. The Journal of Machine Learning Research, 11:2261–2286,
2010.
Ming Yuan and Yi Lin. Model selection and estimation in the Gaussian graphical
model. Biometrika, 94(1):19–35, 2007.
C.-H. Zhang and S. S. Zhang. Confidence intervals for low-dimensional parameters in high-dimensional linear models. Journal of the Royal Statistical
Society: Series B, 76:217–242, 2014.

29

