Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection

Sidhant Misra 1 Marc Vuffray 2 Andrey Y. Lokhov 1 Michael Chertkov 2 3

arXiv:1703.04886v2 [cs.LG] 26 Feb 2018

Abstract

We study the problem of reconstructing the graph
of a sparse Gaussian Graphical Model from independent observations, which is equivalent to
finding non-zero elements of an inverse covariance matrix. For a model of size p and maximum
degree d, the information theoretic lower bound
requires that the number of samples needed for recovering the graph perfectly is at least d log p/κ2 ,
where κ is the minimum normalized non-zero
entry of the inverse covariance matrix. Existing algorithms require additional assumptions to
guarantee perfect graph reconstruction, and consequently, their sample complexity is dependent
on parameters that are not present in the lower
bound. We show that even for very simple examples where the precision matrix is close to singularity, there is a huge gap between the samples
required by existing methods and the lower bound.
We propose an estimator, called S LICE, that consists of a cardinality constrained least-squares regression followed by a thresholding procedure.
We show that S LICE attains a sample complexity
of κ644 d log p, which differs from the lower bound
by only a factor proportional to 1/κ2 . To the best
of our knowledge, this is the first result in the
literature that provides a completely assumptionfree reconstruction guarantee and shows the existence of a polynomial-time algorithm whose sample complexity depends only on the parameters
present in the information theoretic lower bound.

1
Theoretical Division T-5, Los Alamos National Laboratory, Los Alamos, NM 87545, USA 2 Theoretical Division
T-4, Los Alamos National Laboratory, Los Alamos, NM
87545, USA 3 Skolkovo Institute of Science and Technology, 143026 Moscow, Russia. Correspondence to: Sidhant
Misra <sidhant@lanl.gov>, Marc Vuffray <vuffray@lanl.gov>,
Andrey Y. Lokhov <lokhov@lanl.gov>, Michael Chertkov
<chertkov@lanl.gov>.

1. Introduction
Graphical Models (GMs) are powerful and intuitive modelling tools to represent the dependency structure of a set
of random variables. In undirected graphical models, also
known as Markov Random Fields (MRFs), the sparsity pattern of the graph succinctly captures the conditional independence between the variables through the separation property.
In this paper, we study Gaussian Graphical Models (GGMs)
that are described on a graph G = (V, E) with |V| = p and
with maximum degree d. The probability density function
of the multivariate Gaussian distribution1 is given by


p
X
X
det(Θ)
1
P(x) =
exp −
θii x2i −
θij xi xj  ,
p
2
(2π) 2
i∈V

(i,j)∈E

(1)
where Θ is the inverse covariance matrix, also known as
the precision matrix. For the special case of GGMs, the
sparsity pattern of the precision matrix Θ and the graph G
is identical; (i, j) ∈ E iff θij 6= 0.
We consider the problem of reconstructing the underlying
graph G of an unknown sparse GGM from n i.i.d. samples
drawn from the distribution described by Eq. (1), where
samples {xki }i∈V are indexed by k ∈ {1, . . . , n}. We are
interested in finding tractable algorithms that provide an
estimate Ĝ of the graph G using as few samples as possible.
Moreover we require that with high-probability the graph
should be perfectly reconstructed with a given confidence,
meaning that for a prescribed δ > 0, we require that P(Ĝ =
G) > 1 − δ.
This problem has attracted great interest due to its application in numerous fields including protein interactions (Friedman, 2004; Jones et al., 2012), gene regulatory networks
(Basso et al., 2005; Menéndez et al., 2010) and neuroscience
(Huang et al., 2010; Varoquaux et al., 2010).
Numerous methods have been proposed in the literature
to reconstruct sparse GGMs. Meinhausen and Buhlman
in (Meinshausen & Bühlmann, 2006) adopted a regression
based approach and used the L ASSO (Tibshirani, 1996) to estimate the neighborhood of each vertex. Yuan and Lin (Yuan
1

Our analysis only relies on properties of the unbiased empirical covariance estimator and can then be easily extended to the
case of a distribution with non-zero mean.

Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection

& Lin, 2007) introduced the `1 regularized log-likelihood
estimator, commonly known as G RAPH L ASSO, which was
further explored in (d’Aspremont et al., 2008) and (Ravikumar et al., 2009). Anandkumar et al. in (Anandkumar et al.,
2012) used empirical conditional covariance thresholding
and proved structural consistency for walk-summable models. In (Cai et al., 2011), Cai, Liu and Luo introduced an
`1 regularized estimator called C LIME and later an adaptive version called ACLIME (Cai et al., 2016), that can be
implemented by linear programming. With suitable assumptions, all the above methods have been shown to successfully
exploit sparsity and reconstruct the underlying graph G perfectly with a sample complexity of O(log p).

Among other quantities, all the above algorithms exhibit
dependence on the condition number of Θ. The most prevalent assumption in the literature is the so-called Restricted
Eigenvalues (RE), which directly translates to the bound on
the condition number of Θ. As a result, it appears to be a
common belief that the RE condition is necessary to enable
model reconstruction. However, the condition number has
no impact on the IT lower bound in (2). To illustrate this,
consider the following very simple precision matrix


1
κ0
κ0
1
1 −  ,
Θ =  κ0
(4)
κ0 1 − 
1

An information theoretic lower bound for the minimum
number of samples required for perfect graph reconstruction
was derived in (Wang et al., 2010). The lower bound reads




p

 log p−d − 1
log
−
1
2
d
2


,
,
n∗ > max
dκ
dκ


4κ2
log 1 + 1−κ
− 1+(d−1)κ

where 1 −  > κ0 > 0. The minimum normalized edge
strength for the above matrix is given by κ = κ0 . On the
(Θ)
other hand the condition number is given by λλmax
≥ −1 .
min (Θ)
This means that if we keep κ0 fixed and let  → 0, the IT
bound in (2) remains fixed whereas the condition number,
and hence the sample complexity of existing algorithms
diverges to ∞.

(2)
where the parameter κ denotes the minimum normalized
edge strength and is given by
|θij |
κ = min p
.
(i,j)∈E
θii θjj

(3)

The inverse square dependence on κ in (2) appears due to
the difficulty in detecting edges of low strength from their
absence. Noticeably, the IT lower bound only depends on
the dimension p, sparsity d, and minimum edge strength κ.
In contrast, existing estimators in the literature always require additional assumptions to prove perfect graph reconstruction. The analysis of L ASSO in (Meinshausen &
Bühlmann, 2006) assumes certain incoherence properties
of the matrix Θ, reminiscent of the compressed sensing
problem. In (Ravikumar et al., 2009), a variant of the incoherence condition is assumed, which by itself can be hard to
verify. The proof of C LIME (Cai et al., 2011) and ACLIME
(Cai et al., 2016) require that the eigenvalues of the precision
matrix are bounded. The analysis of the conditional covariance thresholding algorithm (Anandkumar et al., 2012) was
conducted only for the class of so-called walk-summable
models. More recently, non-convex optimization based
methods (Johnson et al., 2012; Wang et al., 2016) have been
proposed to improve sample complexity but they require
that the eigenvalues of Θ are bounded. The very popular PC
algorithm for GM selection (Kalisch & Bühlmann, 2007;
Spirtes et al., 2000) requires the so-called strong faithfulness
assumption, which has been shown to eliminate a non-trivial
volume in the space of all precision matrices (Uhler et al.,
2013). Due to these additional assumptions, the sample
complexity of the algorithms depends on extra parameters
that are not present in the IT lower bound (2).

With the above illustration, we see that there is a gap between the information theoretic lower bound and the sample
complexity of known algorithms. It remains an open question whether the IT lower bound in (2) is loose and fails
to capture dependence on additional parameters, or if the
existing algorithms fail to achieve the same parametric dependence as the lower bound in their sample complexity.
In this manuscript, we propose a two stage algorithm consisting of a sparse least-squares regression followed by a
thresholding procedure. We call this algorithm Sparse Leastsquares Inverse Covariance Estimator (S LICE). In a completely assumption-free way we prove that S LICE successfully reconstructs the graph G perfectly
 d+1with
 probability
greater than 1 − δ using d + κ324 log 4p δ
samples. The
sample complexity of our algorithm is off from the lower
bound only by a O(1/κ2 ) factor. More importantly, it depends only on the same set of parameters as the lower bound,
thus closing the information theoretic gap with respect to
parametric dependence. To the best of our knowledge, this
is the first result in the literature that shows the existence of
an algorithm that achieves the same parametric dependence
as the IT lower bound with a completely assumption-free
proof.
The first step in the algorithm adopts a regression perspective and similar to the L ASSO, minimizes an empirical least
squares objective for each vertex. However, instead of using
an `1 penalty we use a cardinality constraint, often referred
to as an `0 constraint, which enforces that the number of
non-zero elements of the regression vector is at most d. The
second step, based on a product and threshold procedure is
unique to our algorithm and is critical to keep the sample
complexity low.

Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection

The computational complexity of algorithm is dictated primarily by the heavier non-convex regression step. In its
bare form, this step is equivalent to an exhaustive search
over all possible size d neighborhoods of a vertex, resulting
in a computational complexity of O(pd+1 ). However, as
we show in our numerical experiments, the problem can be
solved much faster in practice by recasting it as a Mixed Integer Quadratic Program (MIQP), an approach that has been
explored previously in the compressed sensing literature
(Bertsimas et al., 2016). The MIQP formulation leverages
significant recent advances in mixed integer programming
technology and availability of powerful commercial solvers.
The paper is organized as follows: In Section 2 we introduce
our algorithm S LICE and its sub-routines. In Section 3 we
provide rigorous mathematical guarantees on the sample
complexity required to achieve perfect graph reconstruction
with our algorithm. Proofs and analysis of the optimization step and the thresholding procedure can be found in
Section 4 and in Section 5 respectively. Section 6 contains
numerical simulations on synthetic and real data. Section 7
contains conclusions and discussions on future work.

2. Algorithm

s.t.

kβi k0 ≤ d,

(6)

where the `0 -norm counts the number of non-zero components and Σ̂ denotes the empirical covariance matrix whose
components are given by
n

1X k k
xi xj .
n

Σ̂ij =

(7)

k=1

As the notation suggests, the objective in (6) can be rewritten
only in terms of Σ̂ as
Li (βi , Σ̂) = βiT Σ̂īī βi + 2Σ̂iī βi + Σ̂ii ,

(8)

where ī = [p] \ i.
Step 2: Product and Threshold
Once the estimates β̂i have been obtained for all i ∈ V, we
estimate the edge-set Ê through the following thresholding
procedure


q
Ê = (i, j) ∈ V × V : |β̂ij × β̂ji | > κ/2 .
(9)
The estimated graph is then declared as Ĝ = (V, Ê).

In this section we describe the details of the S LICE algorithm, and the interpretation and implementation of its subroutines.
2.1. Details of the S LICE Algorithm
In this section we specify the details of the algorithm,
which consists of a minimization routine followed by a
post-processing routine.
Step 1: Least Squares with `0 -constraint
Our first step consists of minimizing the well-known
quadratic least squares objective with an added `0 constraint.
This step is based on the interpretation of a multivariate
Gaussian distribution where each variable can be expressed
as a noisy linear combination of its neighbors,
X θij
Xi =
− Xj + ηi ,
(5)
θii
j6=i

where ηi is independent from X[p]\i , where [p] =
{1, 2, . . . , p}. The above representation suggests the use
of sparse linear regression, e.g., see (Meinshausen &
Bühlmann, 2006), where the regression coefficients corθ
respond to βij = θij
, and the non-zero values of the coeffiii
cients determine the graph structure.
Define βi = (βij )j6=i . For each i ∈ V, let
2

n
1 X k X
βij xj  ,
β̂i = argmin Li (βi , Σ̂) =
xi +
n
βi ∈Rp−1
k=1

j6=i

2.2. Interpretation as Minimum Conditional Variance
For a given i ∈ V, the `0 -constraint in (6) is tantamount
to an exhaustive search over all candidate neighborhoods
of vertex i, denoted by B̂i , of size exactly d. Note that
neighborhoods of size strictly less than d can be eliminated
from the search as they are strictly contained in bigger candidate neighborhoods. Thus, the optimization in (6) can be
rewritten as
min
B̂i ⊆[p]\i : |B̂i |=d

βi ∈

min

R

p−1 :Supp(β )⊂B̂
i
i

Li (βi , Σ̂).

(10)

Since we will be restricting ourselves to the case when
d < n, each d × d submatrix of Σ̂ has full rank, and the
inner minimization in (10) can be explicitly resolved to get
β̂i (B̂i ) = −Σ̂−1
Σ̂B̂i i .
B̂ B̂
i

(11)

i

The corresponding optimal value is given by
L∗i (B̂i , Σ̂) = Li (β̂i (B̂i ), Σ̂) = Σ̂ii − Σ̂iB̂i Σ̂−1
Σ̂
B̂i B̂i B̂i i
h
i
= Σ̂−1
,
(12)
(iB̂ )(iB̂ )
i

i

11

where in the last line we use the matrix inversion lemma.
The expression in (12) is the empirical conditional variance
of Xi conditioned on XB̂i ,
d i |X ).
L∗i (B̂i , Σ̂) = Var(X
B̂i

(13)

Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection

Thus, the first step of the S LICE algorithm is equivalent to
finding a candidate neighborhood B̂i ∈ [p] for vertex i that
minimizes the empirical conditional variance. To see why
this is a suitable objective, it is helpful to consider the limit
of large number of samples, when the empirical conditional
variance is replaced by the true one. In this case, any set B̂i
that minimizes the conditional variance must be a superset
of the neighborhood of i. This follows from the observation
that conditioning reduces variance in the case of multivariate
Gaussians, and the separation property of GMs,

Var (Xi | XA ) ≥ Var Xi | X[p]\{i} , ∀A ⊆ [p] \ {i},
(14)

Var Xi | X[p]\{i} = Var (Xi | XBi ) ,
(15)
where Bi ⊂ [p] denots the set of neighbors of i.
We remark here that (14) is a property specific to Gaussians, but can also be thought of as a special case of
the more general information theoretic inequality which
states that conditioning reduces entropy. Indeed, (14)
and (15) holds for any graphical model if we replace
Var(.) by the Shannon entropy function H(.), which are
equivalent in the case of Gaussians since H (Xi | XA ) =
1
2 log [2πeVar (Xi | XA )]. Incidentally, minimizing conditional entropy, more popularly referred to as the pseudolikelihood loss function has been proved to reconstruct graph
structure in (Ravikumar et al., 2010) under incoherence assumption, and recently in (Lokhov et al., 2017) without any
assumptions.

3. Sufficient Conditions for Structure
Recovery
In this section we present our main result which states that
the S LICE estimator successfully recovers the graph structure G with a number of samples that depends only on κ, d
and p.
Theorem 1. Given δ > 0, the probability of perfect graph
reconstruction using S LICE is lower bounded as P(Ĝ =
G) > 1 − δ, provided that the number of samples satisfies
 d+1 
32
4p
n − d > 4 log
.
(16)
κ
δ
We prove Theorem 1 through the results below that provide
guarantees for each step of the S LICE estimator.
Proposition 1 (Optimal support contains the true support).
For each i ∈ V, let B̂i ⊂ [p] be the support of the optimal
solution in (6) and let Bi ⊂ [p] be the neighbors of i. Then
for any δ > 0, the support B̂i satisfies Bi ⊆ B̂i with
probability greater than 1 − δ/2, provided that the number
of samples satisfies
 d+1 
32
4p
n − d > 4 log
.
(17)
κ
δ

Proposition 2 (Post-processing Proposition). Assume that


8dp
64
n − d > 2 log
.
(18)
κ
δ
Then with probability greater than 1 − δ/2, the post processing procedure consisting of Product and Threshold terminates with exactly the correct support.
Theorem 1 follows by combining Proposition 1 and Proposition 2 and applying the union bound.

4. Proof of Proposition 1
To prove Proposition 1, it is sufficient to show that if Bi is
not a subset of B̂i then the corresponding objective value
L∗i (B̂i , Σ̂) in (12) is too large. We establish this in two steps.
We first show that when Σ̂ is replaced by the true underlying covariance matrix Σ in (12), then the corresponding
objective value L∗i (Bi , Σ) is smaller than L∗i (B̂, Σ) by a
factor that depends on κ, for any |B̂| = d that does not
contain Bi (see Lemma 1). We then show that the ordering
is preserved even with the estimation error that results when
we replace Σ by Σ̂ provided that we have sufficient number
of samples as specified by Proposition 1. These statements
are formalized in the following lemmas.
Lemma 1 (Multiplicative gap in noiseless optimal solutions). Fix i ∈ V and let Bi ⊂ [p] be the neighbors of i.
Let B̂ ⊂ [p] be any subset such that |B̂| = d and Bi 6⊆ B̂.
Then
L∗i (B̂, Σ) ≥ L∗i (Bi , Σ)(1 − κ2 )−1 .

(19)

Lemma 2 (Large deviations on L∗i (., Σ̂)). Fix i ∈ V and
let 0 <  < 1. Then for every subset A ⊂ [p] \ {i} with
|A| = d, we have
(1 − )L∗i (A, Σ) ≤ L∗i (A, Σ̂) ≤ (1 + )L∗i (A, Σ), (20)
 −(n−d)2 /8
with probability at least 1 − 2 p−1
.
d e
We show that Proposition 1 follows from the above two
lemmas.
Proof of Proposition 1. Combining Lemma 1 and Lemma 2
and using  = κ2 /2 we have for any i ∈ V that the sequence
of inequalities
L∗i (Bi , Σ̂) < (1 + )L∗i (Bi , Σ) < (1 + )(1 − κ2 )L∗i (B̂, Σ)
<

(1 − κ2 )(1 + ) ∗
Li (B̂i , Σ̂) < L∗i (B̂, Σ̂),
1−

is satisfied for all {B̂ ⊂ [p] \ {i} : |B̂| = d, Bi 6⊆ B̂} with
 −(n−d)κ4 /32
probability at least 1 − 2 p−1
. To complete
d e
the proof we take the union bound for all i ∈ [p] to get




p − 1 −(n−d)κ4 /32
P ∃i ∈ [p] : Bi 6⊆ B̂i < 2p
e
< δ/2,
d

Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection

≥ 1 − 4dpe−(n−d+1)κ
where the last
 inequality follows from n − d
4pd+1
32
log
.
κ4
δ

>

5. Proof of Proposition 2
In this section, we prove the post-processing proposition.
From Proposition 1 we have that for all i ∈ V the estimated
neighborhoods B̂i in Step 1 of S LICE satisfies Bi ⊂ B̂i . The
post processing aims at eliminating the vertices u ∈ B̂i \ Bi
so that we obtain the exact neighborhood. For the rest of
the section we assume that the statement in Proposition 1
holds. Our proof of Proposition 2 is based on showing
that the estimates β̂ij are sufficiently accurate. First we
show that the random quantities
h i β̂ij estimated in step one
θ

have the correct mean E β̂ij = θij
, and are conditionally
ii
distributed as normal random variables. Then we show that
although the variance of the normal random variable β̂ij
is dependent on the conditioning, it can be bounded from
above with high probability. These steps are contained in
the following lemmas.
Lemma 3. Fix i ∈ V. The conditional distribution of β̂ij
for any j ∈ B̂i is given by


Θij −1  −1 
,
(21)
,Θ
Σ̂B̂ B̂
β̂ij | Σ̂B̂i B̂i ∼ N
i i jj
Θii ii
where N (., .) denotes the normal distribution.


Lemma 4. For any  > 0 the random variable Σ̂−1
B̂ B̂
i

i

jj

satisfies the following inequality
h

i
(n−d+1)2
−1
−1
8
P Σ̂B̂ B̂
> (1 − ) Θjj ≤ e−
. (22)
i

i

jj

We combine the above two lemmas to prove that the estiθ
mated quantities β̂ij are close to θij
.
ii
Lemma 5. Fix any i ∈ V and j ∈ B̂i . The estimates β̂ij
satisfy the following inequality
!
r
2
θij
κ θjj
P |β̂ij −
|≥
≤ 4e−(n−d+1)κ /64 , (23)
θii
4 θii
where κ denotes the minimum normalized edge strength
defined in (3).
We show that Proposition 2 follows from Lemma 5.
Proof of Proposition 2. Using Lemma 5 and using the
union bound for all i ∈ V and j ∈ B̂i , we get
!
r
θij
κ θjj
P |β̂ij −
|≤
∀j ∈ B̂i , i ∈ V
(24)
θii
4 θii

2

/64

(a)

δ
≥ 1− ,
2

(25)

where 
the implication
(a) is obtained by using n − d >

8dp
64
in the premise of Proposition 2.
κ2 log
δ
Using (25) for both i and j we get with probability greater
than least 1 − 2δ ,
s
!
!
r
|θij | κ θjj
|θij | κ θii
|β̂ij β̂ji | ≥
(26)
−
−
θii
4 θii
θjj
4 θjj
!2
|θij |
κ
= p
−
.
(27)
4
θii θjj
From
q (27), we get that for (i, j) ∈ E the estimates satisfy
|β̂ij ||β̂ji | ≥ 3κ/4 > κ/2. An identical argument can be
q
used to show that |β̂ij ||β̂ji | ≤ κ/4 < κ/2. This proves
that the post-processing step recovers the exact support.

6. Numerical Simulations
6.1. Illustration of the condition number problem
In this section, we construct a very simple counterexample,
consisting of a sequence of matrices with growing condition
number λλmax
but fixed minimum normalized edge strength
min
κ. The primary purpose of this experiment is to demonstrate that the sample complexity of existing reconstruction
algorithms are indeed sensitive to the condition number of
the precision matrix Θ, whereas the sample complexity dictated by the IT lower bound in (2) as well as our proposed
algorithm S LICE shows no such dependence.
The counter example sequence inspired by (4) consists of
a triangle with two weak links and one stronger link and a
collection of independent nodes. This family of GGMs is
parametrized by the following inverse covariance matrix,


1
κ
κ
0

 κ
1
1−
0
,
Θκ,,σ = 
 κ 1−

1
0
1
0
0
0 σ2 I(p−3)×(p−3)
(28)
where 1 −  is the strength of the strong link, κ < 1 − 
is the strength of the weak links and σ 2 is the variance of
the independent nodes. This family of graphs are chosen
such that κ in (28) corresponds to the minimum normalized
edge strength in (3). Note that the maximum degree is
d = 2. This problem can be interpreted as detecting a
triangle within a cloud of independent nodes, a situation
that is very plausible in practice.
The simulations are performed for matrix dimension p =
200 and n = 175 samples which satisfies n < p. We repeat

Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection

the reconstruction procedure 50 times with independent
samples for different values of σ 2 ∈ {1, . . . , 104 } while
κ = 0.4 and  = 0.01 are fixed. The regularizer parameters in ACLIME, L ASSO and G RAPH L ASSO have been
optimized to yield the best possible results for each value of
σ 2 , an advantage that cannot be availed in practice. S LICE
inherently does not have this issue.
For each algorithm we compute its estimate of the normalized link values (1,2) and (1,4),
r
r
θ12 θ21
θ14 θ41
, γ14 =
.
(29)
γ12 =
θ11 θ22
θ11 θ44

succeed with high-probability and are able to correctly identify that there is a link (1, 2) and no link between (1, 4).
However for larger value of σ 2 , the probability of failure
of ACLIME, L ASSO and G RAPH L ASSO is close to one
while S LICE remains insensitive to changes in σ 2 . This
simple example highlights that when the sample complexity of algorithms depends on parameters not present in the
information theoretic bound, the graph reconstruction can
be adversarially affected even by the presence of additional
independent nodes.
1

Pfailure

0.8

We declare that an algorithm fails to reconstruct the graph
whenever γ12 ≤ γ14 : if this condition is satisfied, then links
(1,2) or/and (1,4) are incorrectly reconstructed regardless
of the thresholding procedure. Note that this choice of reconstruction failure criterion is quite generous. It is very
unlikely that one can devise a successful thresholding procedure solely based on the criterion γ12 ≤ γ14 when γ12
and γ14 are close to each other. This is particularly true
for several reconstructions provided by G RAPH L ASSO and
ACLIME as illustrated in Figure 1, whereas the procedure
appears to provide no advantage to SLICE. Note that we
compare normalized link strengths γij , which are invariant
to rescaling of the Θ matrix, instead of matrix element ratios
βij or matrix elements θij . Thus reconstruction based on
the latter quantities would fail for some rescaling of Θ.
γ12

0.5

SLICE

0.4

0.3

0.3

0.2

0.2

0.1

γ14

0
0

0.1

γ12

1

0.2

0.3

γ12

0.5

0.4

0.4

0.1

γ14

0

0.5

0

0.1

γ12

GRAPH LASSO

0.4

0.8

LASSO

0.2

0.3

0.4

0.5

ACLIME

0.3

0.6

0.2

0.4
0.2

γ14

0
0

0.2

0.4

0.6

0.8

1

0.1

γ14

0
0

0.1

0.2

0.3

2

0.4

√

Figure 1. Illustration of reconstruction failure for σ = 1000
We show the scatter plot of reconstructed values γ12 and γ14 obtained through 50 trial reconstructions. We declare that reconstruction fails when γ12 < γ14 which corresponds to points in the lower
right part of the graphs (highlighted in red). S LICE demonstrates an
almost ideal behavior with γ14 = 0 and γ12 > κ/2 = 0.2. Note
that G RAPH L ASSO and ACLIME systematically yield γ14 > γ12
although (1,4) is not an existing edge.

Simulation results are summarized in Figure 2 where the
probability of failure is plotted against the variance σ 2 of the
independent nodes. For σ 2 close to one, all four algorithms

ACLIME

0.6

GRAPH LASSO

0.4

LASSO
SLICE

0.2

σ2

0
100

101

102

103

104

Figure 2. Detecting a triangle in a cloud of independent nodes.
Empirical probability of failure averaged over 50 trial reconstructions with a fixed number of samples n = 175. The values for all
algorithms except S LICE have been optimized over regularization
coefficients. All algorithms except SLICE fail for large σ 2 values.

6.2. Tests on synthetic and real data
In this section, we present several tests on synthetic and real
data. Our goal is merely to illustrate that the use of modern Mixed-Integer Quadratic Programming (MIQP) solvers
such as Gurobi (Gurobi Optimization, 2016) allows one to
run SLICE in a reasonable time even on relatively large
realistic problems.
As a first test, we run SLICE on synthetic random graph
instances of different degrees (d = 3 and d = 4) and sizes
(p = 10, p = 100 and p = 1000). The link strengths γij
have been randomly generated in the ranges [0.2, 0.4] for
d = 3 and [0.2, 0.3] for d = 4 instances. The family of regular random graphs has been chosen to eliminate potential
dependencies on the heterogeneity in the degree distributions. For implementation, we used one possible MIQP
formulation presented in the Supplementary Material, and
the JuMP framework (Dunning et al., 2017) in julia for running the Gurobi solver. The running times for SLICE with
n = 104 samples for each problem instance are presented
in the Table 1. Notice that the practical scaling of running
times is significantly better than what one would expect
from the worst-case complexity O(pd+1 ) for the full graph
reconstruction.
For the illustration on real data, we use the biological data
set related to the Riboflavin production with B. subtilis.
This data set contains the logarithm of the Riboflavin production rate alongside the logarithms of normalized expres-

Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection
Table 1. Comparison of running times for SLICE on various
regular random graphs with n = 104 samples: MIQP Gurobi
solver time and total running time for reconstruction of the neighborhood of one node, and total time for learning the entire graph.
G RAPH
(p, d)

1 NODE
(G UROBI )

1 NODE
( TOTAL )

F ULL PROBLEM
( TOTAL )

(10, 3)
(10, 4)
(100, 3)
(100, 4)
(1000, 3)
(1000, 4)

0.01 SEC
0.03 SEC
0.03 SEC
0.04 SEC
15.7 SEC
92.3 SEC

2.7 SEC
2.8 SEC
2.7 SEC
2.8 SEC
19.3 SEC
96 SEC

7.2 SEC
7.6 SEC
19.8 SEC
21.7 SEC
18 HOURS
29.3 HOURS

tically impact the reconstruction even for simple models as
illustrated in this paper.
In future work, we would like to perform a more exhaustive
numerical study of S LICE with comparison against state
of art methods including the recently developed “Interaction Screening” algorithm (Vuffray et al., 2016) that was
shown to be optimal for reconstructing discrete graphical
models. On the theoretical side, it would be interesting
to know if the information theoretic lower bound on the
sample complexity is perfectly tight with respect to scaling
in all involved parameters, and if this scaling is efficiently
achievable. This requires either improving the IT bound or
to devise a polynomial time reconstruction algorithm that
closes the gap.

A. Proofs
In this Appendix we prove Lemma 3, 4 and 5. In the process we make heavy use of the properties of the inverse
Wishart distribution based on the properties of the Schur
complement stated in the following lemma.
Lemma 6 ((Ouellette, 1981)). Let X ∈ k×k ∼ W (V, l)
be a random matrix distributed according to the Wishart
distribution with parameter V  0 and order l > k − 1. Let
Y = X −1 with Y ∼ W −1 (U, l) where U = V −1 . Let




X11 X12
V11 V12
X=
, V =
X21 X22
V21 V22

R

Figure 3. Graph learned with SLICE from Riboflavin data set.
This real-world data set (Bühlmann et al., 2014) contains p = 101
variables and n = 71 samples. In the reconstruction procedure,
the maximum degree has been set to d = 6.

sion levels of 100 genes that are most responsive to the
Riboflavin production. Hybridization under different fermentation conditions lead to the acquisition of n = 71
samples, see (Bühlmann et al., 2014) for more details and
raw data. The graph reconstructed with SLICE and constraint d = 6 is depicted in the Figure 3. It took about 2.5
days for the algorithm to learn this graph (with the proof of
optimality of the obtained solution) in this high-dimensional
regime. Notice that again the practical running time for
SLICE using MIQP technology is much lower than the one
required to search over the 1014 candidate neighborhoods
of size d = 6. This example is a perfect illustration of a
trade-off between sample and algorithmic complexity in
real-world problems where the collection of samples might
be very costly.

7. Conclusions and Path Forward
In this paper we propose a polynomial time algorithm, called
S LICE, that perfectly recovers the graph of arbitrary Gaussian graphical models. Contrary to existing algorithms, we
prove that the number of samples needed by S LICE neither
depends on the incoherence properties of the precision matrix nor on its conditioning number, but only on parameters
present in the information theoretic lower-bound. This is
critical when samples are scarce or expensive and can dras-

be any compatible block matrix representation of X and
V . Consider block representations for Y, U with the same
dimensions k1 , k2 that satisfy k1 + k2 = k. Then,
(a) The Schur complement of Y11 is distributed as
−1
−1
−1
Y11 − Y12 Y22
Y21 = X11
∼ W −1 (V11
, l)
−1
= W −1 (U11 − U12 U22
U21 , l),
−1
−1
(b) The random matrix Y22
Y21 conditioned on X11
is
distributed as a matrix normal distribution

−1
−1
−1
−1
−1
Y22
Y21 | X11
∼ N U22
U21 , X11
⊗ U22
.

Proof of Lemma 2. Fix i ∈ V and A ⊂ [p] \ {i} with |A| =
d. Using properties of the Wishart distribution we get that,
−1
2
Σ̂ii − Σ̂iA Σ̂−1
AA Σ̂Ai ∼ (Σii − Σi,A ΣAA ΣAi )χn−d (30)

= L∗i (A, Σ)χ2n−d ,

(31)

where χ2t denotes the standard Chi-squared distribution with
t degrees of freedom. Using the Chernoff bound,
2

P(χ2n−d ≶ 1 ∓ ) < e−(n−d)( 2 − 2 log(1+)) < e−(n−d) /8 .
(32)


1

The proof is completed by using the union bound for all
A ⊂ [p] \ i with |A| = d.

Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection

Proof of Lemma 1. From (13) we have that
(a)

(b)

−1
L∗i (Bi , Σ) = Var(Xi |XBi ) = Var(Xi |X[p]\i ) = θii
,
(33)

Proof of Lemma 4. From Lemma 6, the random matrix
Σ̂B̂i B̂i is distributed according to the Wishart distribution


Σ̂B̂i B̂i ∼ W ΣB̂i B̂i , n . Hence

where (a) follows from the separation property of graphical
models, and (b) follows from (12). Similarly,

h

L∗i (B̂, Σ) = Var(Xi |XB̂ ).

(a)

Using the law of total variance we get that
h
i
Var(Xi |XB̂ ) = E Var(Xi |XBi ∪B̂ )|XB̂
i

 h
+ Var E Xi |XBi ∪B̂ |XB̂


X
1
1
θij Xj | XB̂  .
=
+ 2 Var 
θii
θii

(34)

Σ̂−1
B̂ B̂
i

i −1
jj

i

= Σ̂jj − Σ̂j(B̂i \j) Σ̂−1
(B̂ \j)(B̂

∼ (Σjj − Σj(B̂i \j) Σ−1
(B̂ \j)(B̂
i

(42)

(35)


αj−1 = Σjj − Σj(B̂i \j) Σ−1
(B̂ \j)(B̂
i

The probability P

h
i
Σ̂−1
B̂ B̂

(a)

−1
= θii

i

≤ e−(n−d+1)

jj



2

/8



(44)

,

(45)

where (a) follows from (43) and (32).

θij Xj | X[p]\{i,u} 

j∈Bi ∪B̂



−1
(43)

> (1 − )−1 Θjj

= P χ2n−d+1 < (1 − )αj−1 Θ−1
jj


X

Σ(B̂i \j)j

(36)

j∈Bi ∪B̂
−2
≥ θii
Var 

i \j)

= Ψjj − Ψ2ji Ψ−1
ii ≤ Ψjj ≤ Θjj .

i

(a)

Σ(B̂i \j)j )χ2n−d+1

where (a) follows from Lemma 6 and the constant αj can
be bounded from above for every j ∈ B̂i as

Let u ∈ Bi \ B̂. From above, we get



i \j)

Σ̂(B̂i \j)j

, αj χ2n−d+1 ,

j∈Bi ∪B̂

Var(Xi |XB̂ ) − Var(Xi |XBi )


X
−2
θij Xj | XB̂ 
= θii
Var 

i \j)

i

−1 (b)
θii θuu
κ2
. (37)
−
1
≥
Var(X
|X
)
i
B
i
2
θiu
1 − κ2

The inequality (a) follows from the fact that conditioning
reduces variance in Gaussian and observing that B̂ ⊆ [p] \
{i, u}. The inequality (b) follows from (3).
Proof of Lemma 3. For any i ∈ V let Ψ̂(iB̂i )(iB̂i ) =

−1

−1
Σ̂(iB̂i )(iB̂i )
and Ψ(iB̂i )(iB̂i ) = Σ(iB̂i )(iB̂i )
. Using the block matrix decomposition for matrix inverse
Ψ(iB̂i )(iB̂i ) = Θ(iB̂i )(iB̂i ) − Θ(iB̂i )D Θ−1
DD ΘD(iB̂i ) , (38)
where D = [p] \ {i ∪ B̂}. Since Bi ⊆ B̂i , we must have
ΘiD = 0. Hence the matrix Ψ(iB̂i )(iB̂i ) satisfies the following properties for all i and j ∈ Bi :
Ψ(i)(iBi ) = Θ(i)(iBi ) ,

Ψjj ≤ Θjj .

(39)

i
h
Proof of Lemma 5. Define the event E = Σ̂−1
B̂ B̂

≤

(1 − 1 )−1 Θjj . We bound the deviation of β̂ij from

as

i

where (a) follows from (39).

(41)

r

√
√ 
2
≤ 2Φc  1 − 1 n + e−(n−d+1)1 /8 ,

(a)

(46)

where (a) follows by bounding the first term using Lemma 3
and the definition of the event E, and bounding the second
term by the probability of the event E using Lemma 4.
Setting  = κ/4 and 1 = κ/2, and a tight upper bound on
the Gaussian CDF, we get
√
√ 
2
Φc  1 − 1 n + e−(n−d+1)1 /8
≤ e−κ

(40)

jj
θij
θii

!
!
r
θjj
θjj
θij
= P |β̂ij −
|≥
|E
θii
θii
θii
!
r
θjj
θij
c
|≥
|E
+ P |β̂ij −
θii
θii

θij
P |β̂ij −
|≥
θii

2

From Lemma 6, part (b) we get that


Ψij −1  −1 
β̂ij | Σ̂B̂i B̂i ∼ N
,Ψ
Σ̂B̂ B̂
i i jj
Ψii ii


 
Θij −1
(a)
−1
= N
,Θ
Σ̂B̂ B̂
,
i i jj
Θii ii

i

n/64

2

+ e−(n−d+1)κ

/32

2

≤ 2e−(n−d+1)κ

Combining (46) and (47) we complete the proof.

/64

.
(47)

Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection

References
Anandkumar, Animashree, Tan, Vincent YF, Huang, Furong,
and Willsky, Alan S. High-dimensional gaussian graphical model selection: Walk summability and local separation criterion. Journal of Machine Learning Research, 13
(Aug):2293–2337, 2012.
Basso, Katia, Margolin, Adam A, Stolovitzky, Gustavo,
Klein, Ulf, Dalla-Favera, Riccardo, and Califano, Andrea.
Reverse engineering of regulatory networks in human b
cells. Nature genetics, 37(4):382–390, 2005.
Bertsimas, Dimitris, King, Angela, Mazumder, Rahul, et al.
Best subset selection via a modern optimization lens. The
Annals of Statistics, 44(2):813–852, 2016.
Bühlmann, Peter, Kalisch, Markus, and Meier, Lukas. Highdimensional statistics with a view toward applications in
biology. Annual Review of Statistics and Its Application,
1:255–278, 2014.
Cai, T Tony, Liu, Weidong, and Zhou, Harrison H. Estimating sparse precision matrix: Optimal rates of convergence
and adaptive estimation. The Annals of Statistics, 44(2):
455–488, 2016.
Cai, Tony, Liu, Weidong, and Luo, Xi. A constrained `1
minimization approach to sparse precision matrix estimation. Journal of the American Statistical Association, 106
(494):594–607, 2011.
d’Aspremont, Alexandre, Banerjee, Onureena, and
El Ghaoui, Laurent. First-order methods for sparse covariance selection. SIAM Journal on Matrix Analysis and
Applications, 30(1):56–66, 2008.
Dunning, Iain, Huchette, Joey, and Lubin, Miles. Jump: A
modeling language for mathematical optimization. SIAM
Review, 59(2):295–320, 2017.
Friedman, Nir. Inferring cellular networks using probabilistic graphical models. Science, 303(5659):799–805,
2004.
Gurobi Optimization, Inc. Gurobi optimizer reference manual, 2016. URL http://www.gurobi.com.
Huang, Shuai et al. Learning brain connectivity of
Alzheimer’s disease by sparse inverse covariance estimation. NeuroImage, 50(3):935–949, 2010.
Johnson, Christopher, Jalali, Ali, and Ravikumar, Pradeep.
High-dimensional sparse inverse covariance estimation
using greedy methods. In Artificial Intelligence and Statistics, pp. 574–582, 2012.

Jones, David T, Buchan, Daniel WA, Cozzetto, Domenico,
and Pontil, Massimiliano. Psicov: precise structural contact prediction using sparse inverse covariance estimation
on large multiple sequence alignments. Bioinformatics,
28(2):184–190, 2012.
Kalisch, Markus and Bühlmann, Peter.
Estimating
high-dimensional directed acyclic graphs with the pcalgorithm. Journal of Machine Learning Research, 8
(Mar):613–636, 2007.
Lokhov, Andrey Y, Vuffray, Marc, Misra, Sidhant,
and Chertkov, Michael. Optimal structure and parameter learning of Ising models. arXiv preprint
arXiv:1612.05024, 2017.
Meinshausen, Nicolai and Bühlmann, Peter.
Highdimensional graphs and variable selection with the lasso.
The annals of statistics, pp. 1436–1462, 2006.
Menéndez, Patricia, Kourmpetis, Yiannis AI, ter Braak,
Cajo JF, and van Eeuwijk, Fred A. Gene regulatory networks from multifactorial perturbations using graphical
lasso: application to the dream4 challenge. PloS one, 5
(12):e14147, 2010.
Ouellette, Diane Valerie. Schur complements and statistics.
Linear Algebra and its Applications, 36:187–295, 1981.
Ravikumar, Pradeep, Raskutti, Garvesh, Wainwright, Martin J, and Yu, Bin. Model selection in gaussian graphical
models: High-dimensional consistency of `1 -regularized
MLE. In Advances in Neural Information Processing
Systems 21, pp. 1329–1336, 2009.
Ravikumar, Pradeep, Wainwright, Martin J, Lafferty,
John D, et al. High-dimensional Ising model selection
using `1 -regularized logistic regression. The Annals of
Statistics, 38(3):1287–1319, 2010.
Spirtes, Peter, Glymour, Clark N, and Scheines, Richard.
Causation, prediction, and search. MIT press, 2000.
Tibshirani, Robert. Regression shrinkage and selection via
the lasso. Journal of the Royal Statistical Society. Series
B (Methodological), pp. 267–288, 1996.
Uhler, Caroline, Raskutti, Garvesh, Bühlmann, Peter, Yu,
Bin, et al. Geometry of the faithfulness assumption in
causal inference. The Annals of Statistics, 41(2):436–463,
2013.
Varoquaux, Gael, Gramfort, Alexandre, Poline, JeanBaptiste, and Thirion, Bertrand. Brain covariance selection: better individual functional connectivity models
using population prior. In Advances in Neural Information Processing Systems 23, pp. 2334–2342, 2010.

Assumption-Free and Near-Optimal Sparse Inverse Covariance Selection

Vuffray, Marc, Misra, Sidhant, Lokhov, Andrey, and
Chertkov, Michael. Interaction screening: Efficient and
sample-optimal learning of Ising models. In Advances
in Neural Information Processing Systems 29, pp. 2595–
2603, 2016.
Wang, Lingxiao, Ren, Xiang, and Gu, Quanquan. Precision
matrix estimation in high dimensional gaussian graphical
models with faster rates. In Proceedings of the 19th
International Conference on Artificial Intelligence and
Statistics, pp. 177–185, 2016.
Wang, Wei, Wainwright, Martin J, and Ramchandran, Kannan. Information-theoretic bounds on model selection for
gaussian markov random fields. In IEEE International
Symposium on Information Theory Proceedings (ISIT),
pp. 1373–1377, 2010.
Yuan, Ming and Lin, Yi. Model selection and estimation
in the gaussian graphical model. Biometrika, pp. 19–35,
2007.

